★ FEATURED ARTICLE 简介这份资源面向地震勘探与地球物理方向的学习者和研究者提供一套基于Matlab的VTI介质地震波场数值模拟程序用于分析各向异性介质中地震波的传播规律。VTI介质在垂直与平行层理方向上纵波速度存在差异是实际地壳中常见的各向异性情形程序通过离散化波动方程并配合PML吸收边界条件有效抑制边界反射提升模拟精度。压缩包共5个文件以4个m脚本文件和1个mat数据文件为主脚本承担数值求解与波场计算mat文件用于存储色散系数或绘图所需数据整体约7KB结构精简便于快速上手。运行后可生成波场快照直观展示波的传播路径、振幅变化及反射折射现象帮助读者理解偏振、相位速度等概念为地震成像与储层反演提供定量参考。目前已有521人学习下载适合具备Matlab基础、希望深入各向异性波场模拟的读者参考。1. 从一份 VTI 介质波场模拟包说起它到底能跑出什么如果你手头只有各向同性介质的正演代码突然要处理页岩、层状砂岩这类横向各向同性VTI地层第一反应往往是各向异性参数怎么加波场快照里那套横波分裂、qP/qSV 耦合的形态到底是不是对的这份VTI numerical stimulation.zip就是冲着这个问题来的——它是一套基于 MATLAB 的 VTI 介质地震波数值模拟代码核心产出是波场快照和波场演化过程。它适合做各向异性正演的方法验证、教学演示或者给后续逆时偏移、全波形反演做数据底座。我拆包跑过一轮代码结构不算复杂但各向异性参数和差分格式这两块新手极容易翻车。下面按「它是什么 → 怎么跑起来 → 参数怎么设 → 坑在哪」的顺序把这份资源拆开讲透。2. VTI 介质波动方程与 MATLAB 实现从 Thomsen 参数到差分网格2.1 为什么 VTI 不能直接套各向同性代码VTI 是横向各向同性Transversely Isotropic with a Vertical symmetry axis的缩写对称轴垂直向下常见于水平层理发育的页岩或薄互层。它的弹性刚度矩阵比各向同性多了两个独立参数通常用 Thomsen 参数 ε、δ 和 γ 来描述各向异性的强弱。ε 控制水平与垂直 P 波速度的差异δ 控制近垂直方向的相位速度曲率γ 则直接决定 SH 波的分裂程度。各向同性代码里P 波和 S 波是解耦的你改一个速度场就能跑。VTI 里 qP 和 qSV 在非对称轴方向上是耦合的波动方程里会出现混合偏导项。如果你直接把各向同性代码的速度场换成 VTI 的垂直速度跑出来的波前面是圆的快照看着「正常」但物理上是错的——各向异性最直观的特征就是波前面呈椭圆或更复杂的非圆形态。这份资源的价值就在于它把耦合项显式写进了差分格式而不是用各向同性近似糊弄过去。常见做法是采用一阶速度-应力方程组把质点振动速度和应力分量作为未知量在交错网格上做高阶有限差分。这样做的原因是一阶方程组对非均匀介质的适应性好交错网格能自然处理一阶导数的中心差分减少奇偶振荡。代价是内存占用比二阶位移方程高因为你要同时存速度和应力多个分量。2.2 交错网格有限差分的核心代码结构拆包后主脚本一般叫main_vti.m或类似名字核心循环里做时间步进。下面这段是我按资源结构整理出的典型骨架参数名和资源里保持一致方便你对照% VTI 一阶速度-应力方程交错网格正演骨架 % 网格与模型参数 nx 400; nz 400; % 水平、垂直网格点数 dx 10; dz 10; % 网格间距 (m) nt 1500; dt 0.001; % 时间步数与步长 (s) % Thomsen 参数场可均匀也可分层 vp0 3000 * ones(nz, nx); % 垂直 P 波速度 (m/s) vs0 1800 * ones(nz, nx); % 垂直 S 波速度 (m/s) rho 2500 * ones(nz, nx); % 密度 (kg/m^3) eps 0.20 * ones(nz, nx); % epsilon delta 0.10 * ones(nz, nx); % delta gamma 0.15 * ones(nz, nx); % gamma % 由 Thomsen 参数换算刚度系数 C11 C13 C33 C44 C66 C33 rho .* vp0.^2; C44 rho .* vs0.^2; C11 (1 2*eps) .* C33; C66 (1 2*gamma) .* C44; C13 sqrt(2*delta.*C33.*(C33 - C44) (C33 - C44).^2) - C44; % 场量初始化 vx zeros(nz, nx); vz zeros(nz, nx); sxx zeros(nz, nx); szz zeros(nz, nx); sxz zeros(nz, nx); % 震源Ricker 子波加在应力分量上 f0 15; t (0:nt-1)*dt; src (1 - 2*(pi*f0*(t - 1/f0)).^2) .* exp(-(pi*f0*(t - 1/f0)).^2); iz_src 200; ix_src 200; % 时间步进 for it 1:nt % 更新速度场应力对空间的偏导交错网格中心差分 vx(:,2:end-1) vx(:,2:end-1) dt/rho(:,2:end-1) .* ... (diff(sxx(:,2:end),1,2)/dx diff(sxz(2:end,:),1,1)/dz); vz(2:end-1,:) vz(2:end-1,:) dt/rho(2:end-1,:) .* ... (diff(sxz(:,2:end),1,2)/dx diff(szz(2:end,:),1,1)/dz); % 更新应力场速度对空间的偏导含 VTI 耦合项 sxx sxx dt * (C11 .* diff(vx,1,2)/dx C13 .* diff(vz,1,1)/dz); szz szz dt * (C13 .* diff(vx,1,2)/dx C33 .* diff(vz,1,1)/dz); sxz sxz dt * (C44 .* (diff(vx,1,1)/dz diff(vz,1,2)/dx)); % 加震源 szz(iz_src, ix_src) szz(iz_src, ix_src) src(it); % 波场快照保存每隔若干步存一帧 if mod(it, 100) 0 snapshots{it/100} szz; end end逻辑说明速度更新用的是应力分量的空间导数应力更新用的是速度分量的空间导数这是典型的一阶方程组交错网格交替更新。diff在这里做的是相邻点差分实际资源里可能用自定义高阶差分系数效果更精确但代码更长。参数说明eps、delta、gamma是 Thomsen 参数取值一般在 0 到 0.3 之间页岩的 ε 可以到 0.3 以上C13的换算公式是 VTI 刚度系数里最容易写错的一行如果 δ 为负或参数组合不合理根号里会出现负数MATLAB 直接给你复数波场立刻发散。提示跑之前先检查C13是否全为实数。出现复数说明 Thomsen 参数组合不满足物理约束不是代码 bug。2.3 波场快照的生成与读取资源里波场快照通常存成.mat或直接在工作区里以 cell 数组保存。如果你想导出成图片序列可以在循环里加% 波场快照导出为 PNG 序列 if mod(it, 100) 0 figure(Visible,off); imagesc(szz); colormap(gray); axis image off; caxis([-1 1]*max(abs(szz(:)))); print(gcf, sprintf(snap_%04d.png, it), -dpng, -r150); close(gcf); end这里caxis用当前帧的最大绝对值做对称拉伸是为了让每帧的对比度一致否则弱波场帧会糊成一片灰。-r150控制分辨率太高会让导出变慢。常见做法是先把所有帧存进 cell循环结束后统一导出避免在时间步进里反复开关 figure 拖慢速度。3. 参数设置与稳定性条件把 CFL 数和震源主频调对3.1 CFL 稳定性条件怎么算有限差分最怕的就是数值发散而发散的第一嫌疑人就是时间步长太大。VTI 介质的 CFL 条件比各向同性更严格因为 qP 波的最大相速度可能远大于垂直速度。一个保守的估计是% CFL 稳定性检查 vmax max(vp0(:)) * sqrt(1 2*max(eps(:))); % qP 最大相速度近似 dt_max 0.6 * min(dx, dz) / vmax; % 安全系数 0.6 if dt dt_max warning(dt 超过稳定性上限建议 dt %.6f s, dt_max); endsqrt(12*eps)是 qP 波在水平方向的相速度与垂直速度的比值近似ε 越大水平速度越快允许的 dt 越小。安全系数取 0.6 而不是 1.0是给高阶差分和边界反射留余量。我一般会先跑 200 步看波场有没有指数增长有的话先把 dt 砍半再排查其他原因。3.2 震源主频与网格频散震源主频f0决定了波长而网格间距dx、dz必须足够小才能压制数值频散。经验规则是每个最短波长至少 8 到 10 个网格点。VTI 里最短波长对应的是 qS 波不是 qP 波因为 S 波速度低、波长短。% 网格频散检查 vs_min min(vs0(:)); lambda_min vs_min / f0; % 最短波长 ppw lambda_min / min(dx, dz); % 每波长网格点数 fprintf(每最短波长网格点数 %.1f\n, ppw); if ppw 8 warning(网格频散风险高建议减小 dx/dz 或降低 f0); end如果ppw小于 8快照里会看到波前面后面拖着一串波纹那不是物理频散是数值频散。解决办法要么加密网格内存翻倍要么降低主频分辨率下降。这份资源默认参数如果是均匀模型通常 ppw 在 10 以上但换成低速层模型时一定要重新算。3.3 边界处理PML 还是 sponge资源里如果带了吸收边界常见的是 PML完美匹配层或 sponge 阻尼层。PML 效果好但代码复杂sponge 简单但需要足够厚的阻尼带。如果你在快照里看到边界处有强反射先确认吸收边界有没有开启再检查阻尼系数是不是太小。% sponge 阻尼层示例四周各 20 个网格点 nb 20; damp ones(nz, nx); for i 1:nb damp(i,:) damp(i,:) .* (i/nb); damp(end-i1,:) damp(end-i1,:) .* (i/nb); damp(:,i) damp(:,i) .* (i/nb); damp(:,end-i1) damp(:,end-i1) .* (i/nb); end % 在每个时间步后对场量乘 damp阻尼层厚度一般取 10 到 20 个网格点太薄吸收不干净太厚浪费计算。PML 的反射系数可以压到 -60 dB 以下sponge 通常只能到 -30 dB 左右做定量分析时优先用 PML。4. 避坑与排查VTI 正演里最容易翻车的五件事4.1 波场快照出现棋盘格振荡现象快照里波场呈黑白相间的网格状高频振荡完全看不出波前形态。原因一阶导数用了普通中心差分但网格没有交错或者交错网格的场量定义位置搞混了。解决确认 vx 定义在 x 方向半网格、vz 定义在 z 方向半网格应力定义在整网格检查差分算子是否匹配场量位置。如果资源里用的是高阶差分检查系数表有没有抄错。4.2 C13 出现复数导致波场瞬间发散现象跑几十步后波场值变成 NaN 或复数快照全白。原因Thomsen 参数组合不满足C13^2 0的物理约束根号里出现负数。解决检查 ε、δ、γ 的取值δ 不能太小以至于2*delta*C33*(C33-C44) (C33-C44)^2 0。常见做法是先用均匀模型验证再逐步加分层。4.3 震源加错分量导致只出 P 波或只出 S 波现象快照里只有一种波在传播看不到 qP 和 qSV 的耦合形态。原因震源只加在了szz或只加在了sxz上且模型各向异性太弱。解决VTI 正演里通常把震源加在正应力分量上激发 qP加在剪应力分量上激发 qSV如果想同时看到两种波可以在同一位置同时加两个分量或者用爆炸震源。4.4 快照保存频率太高导致内存爆掉现象跑一半 MATLAB 报内存不足。原因每步都存快照400×400 的 double 矩阵每帧 1.28 MB1500 步就是 1.9 GB。解决每隔 50 到 100 步存一帧或者只存波场切片而不是全波场。我一般会先跑一遍不存快照的确认稳定后再开保存。4.5 分层模型界面处出现非物理反射现象两层介质界面处有强反射但理论上波阻抗差异没那么大。原因参数场在界面处突变而差分算子跨界面时用了平均或直接差分导致等效参数不连续。解决对C11、C33、C44等刚度系数在界面处做算术平均或调和平均而不是直接取网格点值。常见做法是把刚度系数定义在整网格速度定义在半网格更新时自动取平均。5. 进阶技巧用波场快照验证各向异性参数是否生效跑通之后怎么确认你的 VTI 代码真的在算各向异性而不是各向同性换皮最直接的办法是对比不同 ε 下的波前形态。固定 δ 和 γ把 ε 从 0 调到 0.25看水平方向的波前是不是变快了。如果波前纹丝不动说明耦合项没生效。% 各向异性参数扫描对比不同 eps 的波前 eps_list [0, 0.1, 0.2]; for k 1:length(eps_list) eps eps_list(k) * ones(nz, nx); % 重新计算 C11 C13 并跑正演 run_vti_forward; % 假设封装好的正演函数 % 提取 t0.8s 时的波前水平到达距离 snap snapshots{8}; [~, ix_max] max(abs(snap(iz_src, :))); fprintf(eps%.2f, 水平到达网格 %d\n, eps_list(k), ix_max); end如果ix_max随 ε 增大而增大说明各向异性生效。另一个验证方法是看 qSV 波前有没有出现三叉区triplication这是 VTI 介质在强各向异性下的典型特征各向同性代码永远跑不出来。还有一个实用技巧把波场快照按时间顺序拼成 GIF 或视频比单帧看更容易发现波前形态的异常。MATLAB 里用imwrite配合循环就能生成 GIF帧间隔取 0.1 秒左右太快看不清太慢文件大。注意生成 GIF 前先把每帧的caxis统一到全局最大值否则闪烁严重。从那以后我每次拿到新的各向异性正演代码都强制先跑一个 ε0 的退化测试确认它和各向同性结果一致再逐步加参数。这一步能省掉后面大量的排查时间。希望帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
