弹性波方程正演模拟：从交错网格到PML吸收边界的完整实践指南

简介弹性波方程是地球物理学中描述地震波传播的核心方程正演模拟则借助数值求解方法预测地震波在介质中的传播路径、波形振幅与到达时刻。这份面向地震波数值模拟学习者与科研人员的MATLAB算法资源利用MATLAB的矩阵运算特性采用10阶差分精度对弹性波方程进行离散求解相比低阶格式能更细腻地刻画高频波场细节适合用来学习高阶有限差分方法的实现思路。压缩包内仅含1个MATLAB脚本文件包体约2KB代码短小精悍便于逐行研读和修改测试但由于未加入吸收边界条件初次运行时更偏重展示无边界情形下的波场基本行为。目前已有133人学习资源可作为高阶地震波模拟的入门样例在此基础上自行添加PML或人工边界条件即可进一步优化模拟结果在实际地质构造中的物理合理性。1. 弹性波方程正演模拟为什么值得做声波近似掩盖了太多真实信息接手过实际工区正演任务的人大概都经历过这种别扭声波方程正演跑得飞快模拟出来单分量地震记录也像模像样可一旦把振幅、相位谱拿去和野外观测数据对比立刻就能看出不对劲——横波信息完全空白AVO 响应不对转换波更是无从谈起。问题不在代码写错了而在物理模型本身。声波方程假设介质只传播纵波地层内没有剪切应力这在储层描述、裂缝预测和流体识别这些真正需要弹性参数约束的场景里等于主动丢掉了一大半有效信息。弹性波方程的地震波场正演模拟简单说就是在给定速度模型、密度模型和衰减参数的条件下用数值方法同时求解纵波和横波在介质中的传播过程输出多分量波场快照和合成地震记录。它解决的问题很具体帮你在没有野外采集数据之前就能验证观测系统能不能收到目标层的转换波帮全波形反演准备敏感度核函数帮微地震事件定位判断震源机制产生的波场特征。适合谁做地震采集设计、VSP 资料处理、FWI 和微地震解释的工程师和研究生只要你的工作对象是弹性波而不是纯声波这篇文章里网格参数怎么定、边界怎么吸收、震源怎么加载这些事都会直接落在你的日活里。2. 弹性波方程的理论基础从双曲型方程组到交错网格有限差分2.1 一阶速度-应力方程为什么数值实现几乎都选它弹性波正演可选的方程形式有两种二阶位移方程和一阶速度-应力方程。教科书里推导方便多用位移表示但真正落到数值计算主流实现几乎都选择一阶速度-应力方程。原因在离散误差上。二阶位移方程直接对空间求二阶导需要至少三个网格点才能构造中心差分而一阶速度-应力方程组里空间导数和时间导数都只到一阶配合交错网格可以让空间差分模板只跨半个网格步长相同网格密度下相位误差更小各向异性伪影更少。对于含自由表面的模型一阶形式处理应力自由边界条件也更直接。在二维各向同性介质中一阶速度-应力方程的标准形式是∂vx/∂t (1/ρ) * (∂σxx/∂x ∂σxz/∂z) ∂vz/∂t (1/ρ) * (∂σxz/∂x ∂σzz/∂z) ∂σxx/∂t (λ2μ) * ∂vx/∂x λ * ∂vz/∂z ∂σzz/∂t λ * ∂vx/∂x (λ2μ) * ∂vz/∂z ∂σxz/∂t μ * (∂vx/∂z ∂vz/∂x)这里 vx、vz 是质点振动速度分量σxx、σzz 是正应力σxz 是剪切应力ρ 是密度λ 和 μ 是拉梅常数。实际程序里一般不直接存 λ 和 μ而是换算成纵波速度 vp、横波速度 vs 和密度 ρ空间任意点的 λ ρ(vp² - 2vs²)μ ρvs²。给参数做换算的时候注意一个常见误用有人为了省事直接在空间上线性插值 vp 和 vs这种做法会在介质分界面处产生虚假的薄层反射因为 vp² 的线性平均并不等于平均介质中的等效速度。我一般先对 vp、vs、ρ 分别做插值再算 λ 和 μ如果模型本身是逐点给定的块状模型插值顺序带来的误差在低频段不明显但涉及高角度界面时还是严格一点好。2.2 交错网格上的空间差分离散精度与实现的平衡交错网格的核心思想是把不同物理量放在半个网格步长的错位位置上。以二维笛卡尔网格为例法向应力 σxx、σzz 定义在整网格点上速度分量 vx 定义在 x 方向半网格偏移点vz 定义在 z 方向半网格偏移点剪切应力 σxz 定义在 x 和 z 同时偏移半网格的位置。这样每个空间导数在相邻半个网格点之间用中心差分构造只需要两点模板就能达到二阶精度。空间差分阶数是精度和计算量的博弈。四阶空间差分每阶导数需要 4 个网格点是性价比最高的选择能在每波长 56 个网格点的情况下把数值频散压到肉眼不可见要模拟长传播距离或做高精度全波形反演可以升到八阶甚至十阶但需要额外处理边界附近的单向差分。时间推进上多数实现用二阶龙格-库塔或蛙跳格式时间步长受 CFL 条件限制二维情况下大约要满足 vmax * dt / dx 0.60.7四阶空间差分实际经验值更紧通常取 0.40.5 更稳妥。一个常见折衷空间上做高阶时间上仍保持二阶。因为时间阶数不够会在波形尾部造成轻微的非对称误差而空间阶数不够会造成波长尺度内的振荡尾巴两种误差的特征不同前者更容易靠吸收边界压住后者会在全剖面到处漏。这一节落到代码里核心就是构造差分系数。四阶中心差分在均匀网格上的系数是 [1/12, -2/3, 0, 2/3, -1/12]/dx交错网格上的半网格差分系数会不一样要按物理量的空间位置分两套系数处理。你要是直接拿整网格系数套到交错网格上波场会立刻出现棋盘状高频振荡。3. 跑通第一个弹性波正演模拟模型、参数与最小代码骨架3.1 设计一个两层介质模型密度、纵波速度、横波速度怎么给第一步先不追求复杂构造用最简单的水平两层模型验证程序正确性。假定模型尺寸 2000 m × 2000 m网格间距 dx dz 10 m即 200 × 200 个网格点。上层介质 vp 2500 m/svs 1400 m/sρ 2200 kg/m³代表风化层下层介质 vp 4000 m/svs 2300 m/sρ 2400 kg/m³代表致密砂岩。分界面放在深度 1000 m 处。参数给定时注意泊松比的约束。在各向同性弹性介质中vp / vs 的下限是 √2低于这个值意味着泊松比为负能量上不稳定上限没有硬约束但 vp / vs 超过 3.5 在沉积岩里是极端的数值上会让 λ 远大于 μ导致正应力更新方程中两项大数相消浮点误差被放大。模型用 NumPy 数组初始化每个物理量一个数组不压缩存储import numpy as np nx, nz 200, 200 dx dz 10.0 vp np.ones((nz, nx)) * 2500.0 vs np.ones((nz, nx)) * 1400.0 rho np.ones((nz, nx)) * 2200.0 vp[nz//2:, :] 4000.0 vs[nz//2:, :] 2300.0 rho[nz//2:, :] 2400.0 mu rho * vs**2 lam rho * vp**2 - 2.0 * mu这段代码里nz//2是数组行方向的中线对应深度 1000 m。λ 和 μ 用vp**2而不是vp直接参与计算目的是保持拉梅常数与速度之间的平方关系后续更新应力时不要写成(lam 2*mu) * dvx_dx lam * dvz_dz以外的东西拉梅常数算错最常见就是把2*mu的系数漏掉或者把lam乘到剪切项上。3.2 稳定性条件与网格尺寸为什么 CFL 数不能随便取网格尺寸和时步的决策直接决定一个正演算例从「能跑」变成「跑得对」。网格间距的约束来自频散控制。通常要求最短波长内至少有 56 个网格点而最短波长取决于横波速度和震源最高频率λ_min vs_min / f_max。如果 vs_min 是 1400 m/s震源最高频率 30 Hzλ_min ≈ 47 m网格间距取 10 m 大约 4.7 个点勉强达到四阶差分的及格线想要更好的波形保真度要么网格加密到 78 m要么把震源主频降到 20 Hz。时间步长的约束来自 CFL 条件。对二维四阶空间差分经验上限在 0.45 附近这里 CFL 数的定义是 vmax * dt / dx。vmax 取整个模型的最大纵波速度 4000 m/sdx 取 10 m 时vmax np.max(vp) cfl 0.45 dt cfl * dx / vmax print(dt)计算得到 dt ≈ 0.001125 s取整为 0.001 s 安全余量更大。注意这里必须用最大纵波速度而不是平均速度或最大横波速度因为弹性波正演里纵波总是跑得比横波快如果按 vs_max 去定 dt模拟会在纵波经过高速度区时直接数值爆炸。时间步数等于模拟时长除以 dt模拟 1 s 的波场记录需要 1000 步这个数量级适合在笔记本上调试。有些实现里会建议用二阶时间精度的 CFL 上限 0.7 甚至 0.9那是针对声波标量方程和低阶差分格式的结论弹性波方程里横波与纵波耦合、自由表面边界条件又会让局部不稳定性提前出现我自己的习惯是把 CFL 控制在 0.350.45慢是慢一点但不用半夜爬起来看程序有没有发散。3.3 最小可运行代码骨架时间循环里到底在算什么用 Python 写一个纯 NumPy 实现不考究性能目的是把每个计算环节摊开看清楚。主循环里只做五件事更新切应力对速度的贡献、更新正应力对速度的贡献、更新正应力、更新剪切应力、加载震源项。nt 1000 vx np.zeros((nz, nx)) vz np.zeros((nz, nx)) sxx np.zeros((nz, nx)) szz np.zeros((nz, nx)) sxz np.zeros((nz, nx)) # 四阶交错网格差分系数dx 方向 c1 9.0 / 8.0 / dx c2 -1.0 / 24.0 / dx for it in range(nt): # 由应力更新速度分量 dvx_dx c1 * (sxx[2:-2, 3:-1] - sxx[2:-2, 1:-3]) \ c2 * (sxx[2:-2, 5:-1] - sxx[2:-2, 0:-4]) # 注意这里垂直方向用 dz 对应系数 vx[2:-2, 2:-2] dt / rho[2:-2, 2:-2] * dvx_dx # 同理更新 vz由 szz 和 sxz 的组合 # ... 此处省略 vz、sxx、szz、sxz 的同类更新 # 震源项以集中力形式加载到 vz 分量 vz[src_z, src_x] dt * source_wavelet[it] / rho[src_z, src_x]上面代码里最值得解释的是索引偏移。交错网格差分时速度点和应力点不在同一套网格索引上卷积模板要在前后各错开两个点四阶精度用 4 个点所以边界上各留两个网格点不做差分更新这些点的值需要另外通过吸收边界或单向差分补充。很多初版实现跑出来的图边界上有一圈高振幅亮斑就是没处理这层索引内外之差导致的。实际跑这种规模的模型纯 NumPy 双层循环慢得离谱正确做法是把所有运算写成数组切片的广播形式利用 NumPy 在 C 层面做循环。上面代码里刻意写了半截因为完整实现需要把五个分量方程全部展开建议自己动手补齐 vx 在 z 方向的应力贡献、vz 的 x/z 两个方向贡献以及三个应力分量的速度散度组合补完这一步你对弹性波正演的理解能超过大半只跑过现成程序的人。4. 吸收边界与震源加载让模拟结果真正可用的两个环节4.1 完全匹配层PML怎么加才能压住人为边界反射模拟区域的四条边如果不做处理波传到边界会被反射回来叠加在有效波场上形成假同相轴。简单吸收带的思路是把边界区域内的波场振幅乘以衰减系数但这会带来边界反射和频散工程上长期使用且效果可靠的是完全匹配层。PML 的核心思想是把波动方程在吸收层内的解替换为一种在数学上与内部解匹配、但在吸收层内指数衰减的变换形式。常见的实现方法是复坐标伸缩法对吸收层内的空间导数算子做一个拉普拉斯域变换引入阻尼剖面 σ(x)。实现时通常需要把速度和应力在 PML 区域内拆分成两个方向的分量比如 vx 的 x 方向贡献和 z 方向贡献分别用不同的衰减系数更新这就让代码量一下翻倍。实际建模时长宽各加 20 个网格点的 PML 层总厚度大约 200 m。衰减剖面从内到外按二次或三次函数从 0 增加到最大值最大值通常取 σ_max -3 * vp_max * log(0.001) / (2 * L_pml)其中 L_pml 是 PML 层物理厚度。注意如果 σ_max 取大了PML 层本身会变成一个引射反射界面取小了又压不住低频反射。判断 PML 是否调好的标准很简单看波场快照在边界附近有没有圆弧状弱能量——有则说明衰减还不够或剖面过渡不光滑。一个更隐蔽的坑在介质参数进入 PML 时发生突变。如果 PML 区域的 vp、vs、ρ 直接沿用到内部边缘的赋值而内部模型在边缘处本身就有横向不均匀性PML 和内部之间的阻抗差会生成虚假反射。稳妥做法是把 PML 区域内的材料参数设成内部边缘值的平滑延拓。4.2 震源子波与加载方式爆炸源、集中力源的区别震源类型决定了波场中 P 波和 S 波能量的比例关系这是弹性波正演和声波正演一个明显的操作差异所在。常见两种源爆炸源只在正应力分量中加载不加载剪切应力理论上只辐射 P 波在均匀各向同性介质中不应有 S 波。垂直集中力源同时扰动正应力和剪切应力天然产生 P 波和 SV 波更接近可控震源或微地震事件的辐射特性。加载时建议用主频 25 Hz 的雷克子波def ricker(freq, dt, nt): t np.arange(nt) * dt t0 1.0 / freq factor (np.pi * freq)**2 return (1.0 - 2.0 * factor * (t - t0)**2) * \ np.exp(-factor * (t - t0)**2)子波零延迟时间 t0 取一个主频周期1/freq确保 t0 时振幅接近零避免零时刻直接加载一个巨大冲量否则波场中会混入直流分量后面的走时分析全是错的。加载位置与方式震源网格点下标 (iz_src, ix_src) 通常放在自由表面以下 510 个网格点避免自由表面反射和直达波在起跳点纠缠。把子波按时间步逐点加到对应分量的网格上使用如下形式source_wavelet ricker(25.0, dt, nt) for it in range(nt): # 在震源点添加垂直集中力源 vz[src_z, src_x] dt * source_wavelet[it] / rho[src_z, src_x] # 如果采用爆炸源则改为 # sxx[src_z, src_x] dt * source_wavelet[it] # szz[src_z, src_x] dt * source_wavelet[it]注意分母上的密度。集中力源加载的是力密度单位体积力物理量纲上是加速度所以加速度 力 / 密度而爆炸源加载的是应力直接叠加不除密度。这个细节搞反的后果是振幅尺度差一个密度量级且 P/S 波相对振幅全乱。还有一类源是纯剪应力源用于模拟横波辐射但实际数据处理中用得少。如果你做微地震正演震源要换成矩张量加载规则是给九个应力分量同时加不同权重的子波。初次练习不必碰这个先用垂直集中力源跑通主流程最重要。5. 弹性波正演模拟避坑5 个高频翻车现场与排查思路5.1 模拟爆炸波场在几百步内出现高频振荡现象程序运行到 200400 步时整个波场快照布满细密的棋盘状纹路振幅指数级增长边界反射还没走到就变成噪声场了。原因时间步长超过 CFL 限值或者局部区域 vp 异常大比如模型里混入了 6000 m/s 的高速体但 dt 按平均值算。另一种可能是不小心把整网格差分系数用到了交错网格上造成高频截断不满足。解决重新按 vmax max(vp) 计算 dtCFL 取 0.4 以下重跑。如果高频振荡只出现在震源附近并随着距离增长而扩散检查空间差分系数 c1、c2 的正负号——四阶交错网格系数的符号顺序写错一个波场会在一个波长的距离内彻底散掉。5.2 近震源处有拖尾振荡像梳子齿一样现象波场快照中震源周围有一圈周期性的细密波纹不是正常环形波前。原因空间网格太粗最短波长内只有两三个网格点四阶差分无法分辨该波长成分数值频散把能量在网格尺度上散开。解决加密网格到最短波长 6 个点以上或降低震源主频。如果不想同时增加计算量可以只把震源附近区域加密但交界面处要插值衔接麻烦不如先整体加密看结果趋势再决定。另一种有效手段是把空间差分阶数从四阶提到八阶但对应的 PML 层差分公式要同步改不能只换系数。5.3 模拟记录里背景有一条水平直线状强振幅现象整个模拟时间窗内固定某一深度处的所有检波器都收到一个水平延伸的强能量条带走时与震源距离无关。原因这是经典的「零频漂移」常见于采用爆炸源但未让子波均值归零。雷克子波如果截断时间窗太短比如只取主频半个周期子波本身带有直流分量加载到应力分量后形成静态偏移波场快照上就是一条水平的直流亮线。检查子波均值是否为零可加一步source_wavelet - np.mean(source_wavelet)我见过不少程序在子波上的这个减法比吸收边界调参更能改善记录质量先做这个再碰 PML。5.4 转换波看得到但极性反了现象合成三分量记录中PS 波转换横波同相轴的极性符号与理论相反或者与纵波记录的相对极性关系不对。原因应力与速度分量之间的符号约定不一致。速度-应力方程写成∂vx/∂t (1/ρ) * ∂σxx/∂x还是 (1/ρ) * (∂σxx/∂x ...)以及剪切应力更新方程中∂vx/∂z与∂vz/∂x的相对正负号不同教材有不同约定。换个符号体系波场极性就整体翻转。解决方式是锚定一个参考解来校验在均匀半空间模型中计算垂直分量上 P 波初动方向P 波到达时质点运动方向应与传播方向同向即向下传播时向下运动按这个基准去校正符号别相信代码注释里写的极性说明。5.5 边界反射压不住PML 失效现象PML 区域内部波场仍有弧形反射回传到内部区域反射波振幅比理论值高一个量级而且反射看起来来自 PML 层内某个位置而不是外部边界。原因最常见是 PML 参数剖面设计错了。衰减系数 σ 从内边界处直接跳到最大值产生了阻抗突变面或 PML 区域的空间差分阶数低于内部区域导致入射波在 PML 入口处发生频散还有一种是把 PML 内部的速度参数设置与内部不连续形成人工反射界面。解决先检查 σ 从 0 平滑过渡到 σ_max 的形态用二次或三次函数再确认 PML 区域差分方式与内部一致最后将 PML 厚度加大到 30 个网格点。调好后做一个空模型测试均匀半空间震源放在中心看边界反射振幅是否低于主波振幅的 1%。6. 用波场快照与合成记录验证模拟走时、振幅与偏振的综合检验跑通一个正演程序只是第一步真正让模拟结果具备说服力的是它能不能通过三层验证解析解对比、走时一致性、偏振特性。我的习惯是每换一种新介质模型先做这三件事再谈后续应用。第一层解析解对比。在各向同性均匀半空间中兰姆问题Lambs problem的解析解给出了垂向集中力源产生的波场精确表达式这是「金标准」。不要拿自己离散模型的数值解直接跟解析解逐点对比因为网格频散和震源离散化总会带来细微差异。我一般取相距震源数十个波长的检波点对比 P 波和 S 波到达时刻的均方根误差误差小于一个时间步长即可认为程序正确。第二层走时检验。在两层模型中由射线理论手工计算 P 波和 PS 转换波的理论走时与合成记录中拾取的走时对比。P 波反射走时公式 t_P sqrt((2h)² x²) / vpPS 转换波走时要分段计算下行 P 波和上行 S 波路径没有解析闭式解可以用简单的二分法求 Snell 定律条件下的最小走时。偏差超过两个时间步长时优先怀疑震源子波零时刻定义和检波器基准面高程。第三层偏振检验。弹性波的三分量合成记录里P 波质点振动方向平行于传播方向S 波垂直于传播方向。取直达 P 波的时窗对水平分量和垂直分量作质点运动图hodogramP 波应该是一条过原点的直线斜率等于入射角的正切如果看到椭圆形状说明三分量增益不匹配或者模拟输出坐标系定义错了。这三层验证做完才能放心把模拟结果用于 AVA 响应分析或者观测系统参数论证。过去我在一个 VSP 正演项目里吃过亏当时跨过了解析解对比直接做观测系统设计结果合成波形里 P 波初动极性没问题但后续的转换波组参与野外数据的对比始终差半拍最后排查发现是 Vz、Vx 分量的坐标轴方向定义与采集记录约定不一致白做了一个月的观测方案。后来就把偏振检验作为每次正演任务结束前必须跑一遍的例行程序再没出过这类错。弹性波正演到这里骨架就齐了方程、网格、CFL、PML、震源、验证一层一层搭起来每一步都有对应的检查办法。调参遇到玄学问题别急着改代码先回去看子波均值和 CFL 数这两个地方占八成翻车原因剩下的二成用参考解对比定位比反复试参数更接近问题本质。希望这篇笔记能帮你少走几段弯路。本文还有配套的精品资源点击获取