VTI介质弹性波正演：MATLAB高阶交错网格与PML吸收边界实现

简介本资源是一套基于MATLAB实现VTI介质弹性波方程高阶交错网格有限差分正演模拟的完整代码方案面向电子信息工程、地球物理、计算数学等专业的本科生与研究生适用于课程设计、期末大作业及毕业设计中对各向异性介质波动数值模拟的学习与实践。压缩包共2个文件1个核心MATLAB脚本.m文件用于算法实现与参数调控1个PNG图像文件展示典型模拟结果总大小仅16KB轻量易部署支持MATLAB 2014a至2021a多版本直接运行。已有154人学习下载。代码采用参数化编程设计关键物理参数如各向异性系数、网格步长、PML厚度等均集中可调注释详尽、逻辑清晰内置PML吸收边界条件以有效压制边界反射显著提升模拟精度附带可直接运行的案例数据无需额外配置即可复现弹性波在VTI介质中的传播过程为理解地震波传播机理与正演建模提供可靠仿真支撑。 做地球物理正演的人应该都清楚一套能稳定跑出VTI介质弹性波场的有限差分程序不只是课程作业里的一个环节更是后面做波形反演、偏移成像、采集观测系统设计的基础支撑。这个项目的标题看起来很长其实拆开就是三件事MATLAB写代码、VTI介质弹性波方程、高阶交错网格有限差分加PML吸收边界。我刚开始接触的时候也踩了不少坑从网格排布搞到半夜、PML边界反射压不下去、到后面波场一算就炸都是一点点调试过来的。这篇就把整个过程完整复盘一遍从方程推导到代码结构再到调试经验给正要上手或者正在挣扎的同学一份可直接照做的参考。1. 项目思路拆解VTI正演到底在解决什么问题1.1 为什么选VTI介质而不是各向同性很多初学者会问为什么要折腾各向异性真实地层里页岩层、薄互层、裂缝发育区地震波速度随传播方向其实是有明显变化的。VTIVertical Transverse Isotropy垂直横向各向同性是描述这种各向异性最常用、也最基础的一种模型它假设介质在水平方向是无限对称的但垂直方向性质不同。理论上只需要5个独立弹性常数就能描述完整刚度张量在二维x-z平面内实际真正参与波场计算的只有4个c11、c13、c33、c44。相比完全各向异性需要21个参数的复杂程度VTI是在物理贴近度和计算成本之间取得很好平衡的选择。这也正是实际勘探中用的最多的各向异性模型之一。很多反演和成像流程的起点都是用VTI正演生成合成记录检验算法的鲁棒性。如果你的正演只支持各向同性那后面做各向异性偏移、参数反演的时候根本没法衔接。1.2 MATLAB在这个项目里的定位与优势我一直不主张拿MATLAB去做超大网格、超长时间的正演——那是Fortran或C的领地。但作为研究验证和教学验证MATLAB有不可替代的优势矩阵操作天然友好画波场快照一个imagesc就搞定调试时能实时看到变量变化而且代码可读性高方便跟论文公式逐行对照。这个项目选MATLAB正适合用来把算法核心搞清楚。当然MATLAB写有限差分有个大坑如果你按着C语言的习惯用嵌套循环逐个网格点更新那3D模型基本跑不动2D模型也会慢得让人抓狂。正确的做法是尽量向量化用数组切片和矩阵运算替换底层循环。后面我会单独讲性能优化这里先提个醒。网格规模在1000×1000以内、时间步在2000步以内的模拟向量化的MATLAB代码完全可以做到运行时间内完成这才体现了MATLAB作为设计与验证语言的价值。1.3 项目的完整技术栈一览这个项目涉及的知识点其实可以列出很长一串弹性波一阶速度-应力方程、VTI刚度矩阵、交错网格空间差分、时间二阶差分、雷克子波震源加载、PML边界条件、CFL稳定性条件、数值频散抑制、MATLAB向量化编程、波场快照与地震记录可视化。每一个点单独拉出来都能写一篇详细文章但它们在正演程序里是环环相扣的。我建议读者不要只盯着某一个环节而是在动手写代码前先把整套逻辑串起来否则很可能出现“方程写对了但网格排布错位”“PML公式对但参数没调好”这类核心bug。2. VTI介质弹性波方程从数学到物理2.1 速度-应力一阶方程组的推导做有限差分正演推荐使用一阶速度-应力方程而不是二阶位移方程。原因是交错网格天然适合一阶方程组空间差分只需要计算一阶导数可以把精度做得很高同时变量物理意义明确速度分量和应力分量直接对应网格点上需要更新的状态。二维情况下x为水平方向z为垂直方向z向下为正VTI弹性介质的速度-应力方程可以写成这样ρ ∂vx/∂t ∂σxx/∂x ∂σxz/∂zρ ∂vz/∂t ∂σxz/∂x ∂σzz/∂z∂σxx/∂t c11 ∂vx/∂x c13 ∂vz/∂z∂σzz/∂t c13 ∂vx/∂x c33 ∂vz/∂z∂σxz/∂t c44 (∂vx/∂z ∂vz/∂x)四个弹性常数加一个密度ρ就是VTI二维正演的全部材料参数。对比各向同性的情况——那里是λ和μ两个参数现在变成4个多出来的部分正是描述各向异性强度所必需的。这组方程的物理含义是速度场的变化取决于当前应力场的空间梯度应力场的变化又取决于速度场的空间梯度。二者互相耦合波场就在这种“速度-应力交替更新”中一步步向前传播。这也决定了时间步进的方式速度在n1/2时刻更新应力在n1时刻更新交替推进这就是经典的蛙跳格式。2.2 从Thomsen参数到刚度系数实际建模时大家很少直接给c11、c13这种量级的刚度值而是给出Vp0、Vs0、Thomsen各向异性参数ε和δ再由公式换算成刚度系数。这个换算关系值得写清楚因为如果你用的是别人给的各向异性模型文件里面大概率存的是Thomsen参数而不是刚度张量。基本换算关系如下c33 ρ Vp0²c44 ρ Vs0²c11 c33 (1 2ε)c13 sqrt(2δ c33 (c33 - c44) (c33 - c44)²) - c44举一个具体数字设密度ρ2200 kg/m³Vp03500 m/sVs02000 m/sε0.2δ0.1。那么c332.695e10 Pac448.8e9 Pac113.773e10 Pac13按公式算出来大约为1.187e10 Pa。把这些参数作为初始文件输入程序运行出来的波场就能体现VTI介质的典型特征qP波和qSV波不再像各向同性介质那样呈圆对称传播而是呈现随角度变化的方向速度差异。注意实际代码里我建议直接用刚度系数做内部计算只在数据输入输出时保留Thomsen参数格式。这样既方便别人理解你的模型又避免了程序内部频繁换算带来的混乱。2.3 CFL稳定性条件与网格参数设计有限差分不是随便给个时间步长就能跑。显式时间推进要求时间步长满足CFL条件否则波场会数值发散。对于交错网格二维弹性波方程稳定性条件可以写成Δt ≤ 1 / (Vmax × sqrt(1/dx² 1/dz²))Vmax是介质中最大相速度dx和dz是空间网格间距。如果dxdz10mVmax4000m/s那么Δt上限约为1/(4000×sqrt(2)/10) 1.77 ms。实际取值一般取上限的60%到80%因为理论推导假设的是均匀介质且忽略了数值色散对高频分量的影响留一点安全余量是必须的。空间网格密度则要满足每个最小波长至少覆盖6到10个网格点。最小波长等于Vmin除以最高有效频率。如果Vmin2500 m/s震源主频30Hz有效频率上限约60Hz那么最小波长约42m对应网格间距应该控制在5到7m以内。网格太粗跑出来会看到明显的“频率尾随”现象——波前面拖着大尾巴这就是数值频散。这块我在后面常见问题里还会再提。3. 高阶交错网格为什么非“交错”不可3.1 普通网格与交错网格的差别如果直接用一个普通网格把所有变量定义在同一位置那么你要计算应力对x的偏导时就要跨过一个网格点空间差分实际上的“有效步长”是2dx精度损失明显。交错网格的巧妙之处在于把不同的变量放在错开半个网格点的位置上使得每个导数都能在相邻两个采样点之间近似相当于把差分步长缩短了一半精度和稳定性都更好。具体到弹性波二阶方程常见的排布方式是σxx和σzz定义在整数网格点(i, j)vx定义在水平方向半格点(i1/2, j)vz定义在垂直方向半格点(i, j1/2)σxz定义在两个方向都是半格点的位置(i1/2, j1/2)这种错开半个格点的布局正好让每个偏导数在求值时都处于“对称采样”的位置自然比普通网格更加精确。这也是为什么绝大多数的地震波有限差分程序都采用交错网格这个设计从20世纪90年代Virieux等人的工作就奠定了到现在仍旧是主流。3.2 高阶空间差分系数与精度对比交错网格上的一阶导数可以用多个相邻点的加权组合来近似阶数越高需要的点数越多但对波数的分辨率越好相同网格间距下数值频散越小。一阶导数的高阶差分公式可以写成∂f/∂x ≈ (1/dx) Σ an [f(x (n-1/2)dx) - f(x - (n-1/2)dx)]n1,2,...,N这里的系数an是精确确定的。2阶、4阶、6阶严格说是2N阶的前几个系数如下阶数a1a2a3a42阶1---4阶9/8-1/24--6阶75/64-25/3843/640-8阶1225/1024-245/307249/5120-5/7168我自己实际用的最多的是4阶和8阶。4阶存储量小、写起来简单适合模型规模大、频率不高的工业测试8阶数值频散更小可以容忍更大的网格间距在小范围精细模型上效果很好。你完全可以根据网格密度切换阶数代码里把差分系数放在数组里改起来特别方便。3.3 离散格式的时间推进空间差分确定后时间推进采用二阶中心差分。完整的一次更新流程是用当前应力场计算vx和vz在新半时刻的值用更新后的速度场计算各个应力分量在新整时刻的值交替进行直到达到总模拟时长用伪代码表示就是for it 1 : nt % 速度更新应力场已知 vx vx dt/rho .* (Dxx(sxx) Dxz(sxz)); vz vz dt/rho .* (Dxz(sxz) Dzz(szz)); % 应力更新速度场已知 sxx sxx dt .* (c11 .* Dxx(vx) c13 .* Dzz(vz)); szz szz dt .* (c13 .* Dxx(vx) c33 .* Dzz(vz)); sxz sxz dt .* c44 .* (Dxz(vx) Dxz(vz)); % 震源加载、边界吸收、波场记录 end这里的Dxx、Dzz、Dxz都是空间差分算子实现方式就是按照上一节说的交错网格位置对数组做切片加减运算。注意时间步进只做到二阶精度很少有人在时间方向上也用高阶格式。因为时间上的高阶格式需要保存多个历史时刻的波场内存消耗成倍增加收益却没有空间阶数提高来得明显。工程实践中空间4阶时间2阶是长期验证过的黄金组合。4. PML吸收边界把人工反射压到看不见4.1 人工边界反射为什么必须处理数值模拟只能用有限大小的区域代表无限大地层人为截断的边界如果不加处理波场一传到边界就会反射回来和你关心的有效信号混在一起观测记录直接变成一锅粥。早期有人用简单的衰减边界也有人用旁轴近似但效果都不稳定。PMLPerfectly Matched Layer完美匹配层出现之后边界吸收问题有了一个理论上几乎无反射的优雅方案现在基本成了标配。PML的思路是在计算区域外围加一圈损耗材料层波进入这个区域后被按指数方式衰减吸收理论上在截断界面处可以做到零反射。要做到效果好PML层的厚度、衰减剖面的形状、参数取值都有讲究。我见过不少同学花了很多功夫把内部区域代码写得很漂亮结果PML参数没调好波从边界反射回来整个剖面密密麻麻全是假轴还以为是自己方程写错了。4.2 CPML的实现原理传统的分裂式PML需要把每个波场分量按方向分裂开方程数量瞬间翻倍写起来繁琐。现在更推荐使用CPMLConvolutional PML卷积完美匹配层。CPML不需要分裂变量只需要引入一组辅助变量记录边界的“记忆效应”实现起来清晰很多。CPML的核心思想是对空间导数做复坐标拉伸替换∂/∂x → (1/κx) ∂/∂x ψx其中ψx是引入的辅助变量在时间域通过一个递归公式更新ψx^(n1) bx ψx^n ax (∂u/∂x) / dx这里的系数bx和ax由PML的衰减剖面决定bx exp(-(dx αx) Δt)ax dx (bx - 1) / (dx αx)其中dx是PML区域的吸收强度剖面κx是坐标伸缩系数αx是一个帮助吸收低频信号的小参数。衰减剖面通常取为d(x) d_max × (x/L)^NN一般取2或3d_max的经验公式为d_max -(N1) Vmax ln(R) / (2L)R是期望反射系数取0.001或0.0001L是PML厚度单位米Vmax是介质最大波速。举个例子若Vmax4000m/sPML厚度L100mN2R0.001则d_max约等于-3×4000×(-6.9)/(200) ≈ 414 s⁻¹。这个量级的d_max配合内部平滑衰减边界反射可以被压制到远低于有限差分本身的数值噪声肉眼几乎无法察觉。4.3 在交错网格里实现CPML需要注意的细节CPML实现中最容易出错的点是半格点位置上的辅助变量也要跟着交错。例如vx在(i1/2, j)上那它对x求导时使用的ψx就应该也定义在(i1/2, j)附近而σxz在(i1/2, j1/2)上它的x方向辅助变量和z方向辅助变量需要分别定义内存和更新逻辑都要对应上。我建议用三个三维数组分别保存ψxx、ψzz、ψxz每个方向一个记忆变量来统一管理。具体来说速度更新时对σxx的x导数、σxz的z导数要加ψ应力更新时对vx的x导数、vz的z导数、vx的z导数、vz的x导数都要加ψ一句话PML区域内的每个空间导数都要“过一遍”记忆变量更新而不是只在边界那一层做处理。PML厚度我一般取20到30个网格点太少吸收效果差太多浪费计算量。对于验证性测试20个点配合二阶剖面已经足够干净。5. MATLAB程序架构与关键实现细节5.1 主程序结构与模块划分写正演程序不是把公式堆在一起就行好的结构能让你后期调试和扩展事半功倍。我这个项目的MATLAB代码按照功能拆成几个模块主脚本只负责串联整个流程%% 主函数 % 1. 参数设置网格、速度模型、震源、PML % 2. 初始化波场数组、差分系数、PML系数 % 3. 时间循环速度更新、应力更新、震源加载、PML更新、波场快照 % 4. 地震记录保存与绘图一个清晰的主函数大概长下面这个样子clear; clc; close all; % 模型参数 nx 600; nz 600; dx 10; dz 10; dt 1e-3; nt 1500; % 模型介质参数均匀VTI模型 rho 2200 * ones(nz, nx); c11 3.773e10 * ones(nz, nx); c13 1.187e10 * ones(nz, nx); c33 2.695e10 * ones(nz, nx); c44 8.800e9 * ones(nz, nx); % 震源参数 f0 25; t0 1 / f0; xs nx/2; zs nz/2; % 边界与PML参数 npml 20; ...虽然某些参数在后面会被替换成更精确的空间变化模型但这个骨架已经足够支撑完整的正演流程。后续要做分层模型只需要把c11、c33这些参数从常数改成按深度区分的数组即可。5.2 差分算子与交错网格的向量化实现交错网格差分在MATLAB里用数组切片来实现效率远比循环高。以4阶精度的x方向一阶导数作用于定义在整数网格上的变量为例核心代码是function dfdx diff_x_4(f, dx) % f定义在(i, j)上输出df/dx也定义在(i, j)上 dfdx zeros(size(f)); dfdx(:, :, 1:end) ... dfdx(:, 3:end-2) (9/8)*(f(:, 4:end-1) - f(:, 2:end-3)) / dx ... - (1/24)*(f(:, 5:end) - f(:, 1:end-4)) / dx; end实际代码还要处理边界上的单边差分但一般情况下PML区域内的导数计算会单独处理内部区域用这个格式已经足够。这里的关键是搞清楚每个变量在网格上的精确位置把切片的偏移量和交错位置对应起来。我在调试时吃过不少亏最后的方法是画一张网格位置示意图贴在屏幕旁边每次写切片都对照着看。5.3 震源加载与地震记录震源采用雷克子波Ricker wavelet这是地震正演里最常用的震源时间函数for it 1:nt t (it-1) * dt; source (1 - 2*pi^2*f0^2*(t-t0)^2) * exp(-pi^2*f0^2*(t-t0)^2); % 激振力源加载到vx或szz上 sxx(zs, xs) sxx(zs, xs) source; szz(zs, xs) szz(zs, xs) source; ... end实际震源可以设置成加在应力分量上的力源或者加在速度分量上的位移源。最直接的测试方式是加爆炸源同时激发σxx和σzz这样在均匀VTI介质中激发的波场以qP波为主比较容易判断波场形态是否合理。记录地震记录的方式是在地表或特定深度放一排接收点每个时刻把相应位置的某个波场分量记录下来seis(it, :) szz(rec_z, rec_x); % 记录垂直分量一个常见的做法是把接收点放行在z10m的浅层位置每5个网格点一个检波器这样能得到类似野外勘探的单炮记录。5.4 性能优化让MATLAB跑得更快一点我前面反复强调向量化这里说几个实测有效的经验循环内不要做动态数组扩容所有波场数组在进入时间循环前就分配好内存不要在时间循环里直接调用plot或imagesc画图每20-50步截一张快照存数组等模拟结束后统一绘图差分系数的计算放在循环外避免重复计算如果模型是均匀的尽量利用“标量系数与数组乘法的分离”把系数乘到差分结果上而不是每个循环重新读取参数数组在MATLAB里用单精度而不是双精度能显著减少内存带宽压力对于大部分正演精度来说完全足够我实际跑600×600网格、1500步的均匀VTI模型向量化的代码在一般的笔记本上大概需要两三分钟。如果换成循环写法估计跑完需要几小时。差别就是这么明显。6. 正演结果分析与精度验证6.1 均匀VTI模型波场形态验证先用均匀模型做测试好处是理论预期非常明确。VTI介质中的点源激发的波场包含两个传播模态qP波准纵波和qSV波准横波。qP波波前面是椭圆形长轴方向沿水平方向——因为c11大于c33水平方向传播速度更快。qSV波波前形态更加复杂在内侧还会出现三叉区triplication特征这是各向异性介质中群速度与相速度方向差异导致的经典现象。我用程序跑出来的波场快照在250ms时刻能看到明显的椭圆形外圈那是对应qP波的波前内圈有一团形态比较复杂的波场对应qSV波。如果换成各向同性介质波前就是标准圆弧所以这个差异是验证各向异性实现是否正确的最直接依据。6.2 PML效果检查PML效果怎么看最直观的方法是跑一个同参数的大模型让边界足够远使得边界反射到达记录组的时间晚于整个模拟时长把它的结果作为参考然后对比小模型加PML的结果。两者一致说明PML吸收OK。另一个快速办法是直接观察边界附近的波场快照正常PML工作的情况下波场到达PML内部后幅度迅速衰减内部边界处不应该出现明显的二次弧形反射。如果你发现边界处有明显的反射圆弧回来先别急着调公式。常见情况是PML区的衰减剖面d_max给得太小或者PML层数太少。简单调整办法是提高R的阶数从0.001改到0.0001、增大PML厚度或者把剖面指数N从2改成3。反向也可能出问题d_max太大造成数值离散化误差增大这时反而会产生数值反射。我遇到过一个非常典型的案例就是d_max取值过大导致边界区域波速虚部过大出现比正常反射还严重的伪信号。这个坑值得大家留意。6.3 与解析解或参考解的对比要做定量验证可以找一个简单的VTI模型与半解析方法比如广义射线法、反射透射系数矩阵方法对比。但对初学者来说建议用一种更简单的自洽性验证把ε和δ都设成0VTI退化为各向同性介质这时用传统的各向同性解析解比如声波或弹性波的格林函数来对照你的结果是否合理。如果VTI代码在这组参数下能还原各向同性的结果说明主要代码框架是可靠的剩下的误差就出在各向异性参数处理上排查范围就缩小了很多。我这里记录一个更直接的定量验证在2D均匀弹性介质中点源激发的P波到达某个观测点的走时可以用已知波速计算。比如距离800mVp3500m/s理论走时约0.229s。从合成地震记录上读到同相轴的峰值位置误差在几个采样点以内就算正常。我的程序实测误差基本在2个时间步长以内说明时间推进和空间离散都处理得比较准确。7. 调试经验与常见问题速查7.1 波场发散、数值爆炸怎么排查波场跑到一半开始出现NaN或者振幅指数增长这是所有有限差分玩家都会遇到的头号问题。我的排查顺序是先检查CFL条件是否满足很多爆炸就是dt取值太大引起的检查差分系数符号是否搞反交错网格的系数次序错一位就会导致差分格式不稳定检查震源加载处是否合理加载幅度过大也会在局部引发数值自激检查PML区域的d_max是否异常内部参数突变会造成波场扭曲如果程序在初始几步就爆炸多半是差分算子有bug如果是在一两百步后才开始发散重点检查PML区域和时间步长这个区分能帮你快速缩小排查范围。7.2 边界反射压不下去或PML区域出现毛刺边界反射问题是最让人头痛的。这里给出一个我自己屡试不爽的排错路径第一步把PML区域的所有衰减项直接关闭设d_max0看看会发生什么。如果这时候边界反射非常强烈且规律性极强说明PML开关逻辑本身没问题第二步逐步增大PML厚度从10层开始测直到反射消失第三步检查PML内部是否有“硬边界”式的网格编号错误——比如把PML区域的第1个网格和内部区域的计算公式搞混实际经验里80%的PML问题是出在“变量交错位置与记忆变量错位”上。我的建议是先在1D问题里验证PML代码逻辑再扩展到2D排查成本能降低一个量级。7.3 MATLAB常用调试技巧调试有限差分程序最实用的就是加断点分步检查波场快照的对称性。均匀模型下波场必须关于震源位置对称如果不对称十有八九是差分算子的切片索引写错了。用imagesc画一下波场肉眼看是否对称比打印数值快得多。另外把核心差分函数单独拎出来做单元测试。比如构造一个简单函数f(x)sin(kx)解析导数是k cos(kx)把差分结果和解析结果对比。这个测试一旦通过你就能确定差分算子本身没问题后面出错就是网格布局或者PML逻辑的问题。7.4 内存溢出与计算时间过长2D模型的内存压力虽然比3D小但如果你把每个波场分量都存成double数组600×600×5个变量就已经有14.4MB如果还要存每个时刻的快照很快内存就吃紧。我的建议是波场本身用单精度读取显示时再转double波场快照只存你需要的那几个时刻不要每步都存接收记录用double存没问题但不要把所有时刻的完整波场都留在内存里如果真的需要长时间大模型考虑用MATLAB Coder把核心循环转成mex文件能获得接近C语言的性能我个人在实际操作中的体会是很少有什么bug是孤立的技术难题多数都是“网格位置没搞清、差分系数配错、PML参数不合理、循环里不小心改了不该改的数组”这类低级错误的排列组合。把排查顺序固定下来一条一条过反而比苦思冥想更有效。这个程序后续还可以往很多方向扩展比如加自由地表边界做弹性波全波形模拟、改成声波方程做快速正演、加进GPU加速用大型模型去试真实测线。MATLAB本身在这类研究的原型验证阶段效率极高把这个VTI正演基础打好后面无论往哪个方向深入心里都有底。本文还有配套的精品资源点击获取