弹性波正演模拟：交错网格有限差分、参数定标与工程实践

简介：弹性波方程正演模拟是地震学与地球物理勘探中的基础环节，这份压缩包面向相关专业学生与科研人员，提供基于10阶精度差分算法的MATLAB实现脚本。包内仅含1个m文件，压缩包仅2KB，文件虽小却展示了高阶精度离散化弹性波方程的核心思路，可作为学习高阶有限差分正演模拟的入门参考。已有133人学习下载。脚本未设置边界条件，直接反映了理想化介质中的波场传播特征，读者可结合该代码体会波动方程离散化、递推求解与波场快照输出的实现细节，并在此基础上自行扩展PML或人工边界条件，进一步研究实际地质模型的反射与衰减特性。整体上，这份资源适合对地震波数值模拟感兴趣、具备一定MATLAB基础的初学者，代码结构简洁，便于阅读与改进练习。

1. 弹性波方程正演模拟：先搞清楚成本从哪来

弹性波方程的地震波场正演模拟，很多团队一上来就盯着方程本身，真正让项目卡住的往往是网格、步长和内存的三角关系。同样一个二维模型，震源主频从 30Hz 提到 60Hz，空间网格间距要缩一半，时间步长跟着缩，内存压力变成原来的 8 倍以上，迭代步数还会翻倍。正演模拟的本质不是把波场算出来，而是「在可接受的代价内把波场算准」。

这里的输入是速度模型、密度模型和震源子波，输出是各时刻的波场快照和检波点处的合成地震记录。它服务于地震成像的反演算子、观测系统设计、复杂介质响应分析，也常被拿来验证新反演算法。适合做地震资料处理、主动源勘探、城市噪声模拟的工程师，以及想转入数值计算领域的平台研发。下面从方程选择、离散实现到参数定标，按一套可复现的标准流程走。

2. 从位移二阶到速度-应力一阶：弹性波方程的两种主流离散形式

2.1 位移方程直观但难离散，速度-应力方程才是工程主力

弹性波动最直观的写法是位移形式。各向同性介质中，位移向量 u 满足：

ρ ∂²u/∂t² = ∂/∂x_j [λ δ_ij ∂u_k/∂x_k + μ (∂u_i/∂x_j + ∂u_j/∂x_i)] + f_i

未知量少，二维只有两个分量，这个优点在做解析推导时很诱人。但离散时问题就来了：方程里出现二阶空间导数，在标准规则网格上需要 9 点甚至更多模板；在介质分界面处，二阶导数要求位移场具备更高的连续性，而实际界面两侧位移连续、应力连续，速度却可能间断。用二阶导数离散，界面附近的误差会被放大，直接体现在合成记录的高频尾巴上。

工程上更常用的是一阶速度-应力方程，把速度 v 和应力张量 τ 一并作为未知量：

ρ ∂v_i/∂t = ∂τ_ij/∂x_j + f_i∂τ_ij/∂t = λ δ_ij (∂v_k/∂x_k) + μ (∂v_i/∂x_j + ∂v_j/∂x_i) + M_ij

其中 f 是力源项，M_ij 是矩张量源项，爆炸源时取 M_xx = M_zz = 常数。这样处理有三个实打实的好处：

• 只需要一阶空间导数，交错网格上可以用相邻两个网格点对半格点做中心差分，数值各向异性显著降低
• 应力分量的边界条件天然对应自由表面的零应力条件，不需要对位移场做特殊约束
• 扩展到各向异性介质时，只需改应力更新式中的弹性系数张量，整体结构不动

代价是未知量变多：二维需要 vx、vz、τxx、τzz、τxz 五个数组，三维需要九个数组。内存占用比位移形式大约多 1.5 倍，但换来的稳定性和精度收益，在大规模三维计算中仍是最主流的选择。

2.2 时间递推用蛙跳格式，内存最少

一阶速度-应力方程在时间方向上的标准做法是蛙跳格式：速度更新到 n+1/2 时刻，应力更新到 n+1 时刻，交替推进。更新式写作：

v^(n+1/2) = v^(n-1/2) + Δt/ρ · D τ^n + Δt/ρ · f^nτ^(n+1) = τ^n + Δt · C : D v^(n+1/2) + Δt · M^(n+1/2)

D 是空间差分算子，C 是弹性刚度张量。蛙跳格式需要保存的速度和应力都只有一份历史数据，内存开销最小，计算量也最小。它在一阶双曲型问题上的稳定性由 CFL 条件约束，第 4 章会具体计算。

有人倾向用四阶 Runge-Kutta 做时间积分，认为精度更高。实际上在同样的网格配置下，RK4 能把时间步长放宽 2 到 3 倍，但每步要做 4 轮空间导数，临时数组至少翻一倍，综合算力往往比蛙跳多 50% 以上。对绝大多数勘探频率范围的弹性波正演模拟，二阶蛙跳已经足够，只有做长时间波场传播、需要严格控制时间频散的研究型代码，才会换高阶时间格式。

2.3 交错网格有限差分、有限元、谱元怎么选

空间离散方法的选择直接影响网格规模和生产效率。把四类方法的定位放在一张表里看更清楚：

从表里能看出，交错网格有限差分在每波长网格数上不是最省的，但胜在模型参数直接落在格点上、无需网格剖分，复杂程度可控。后面的实现都围绕它展开。

3. 用交错网格有限差分把弹性波方程跑起来：时间递推、震源与 PML

3.1 二维最小实现：五个数组交替更新

先固定二维各向同性介质的一阶速度-应力方程组分量形式：

ρ ∂vx/∂t = ∂τxx/∂x + ∂τxz/∂zρ ∂vz/∂t = ∂τxz/∂x + ∂τzz/∂z∂τxx/∂t = (λ+2μ) ∂vx/∂x + λ ∂vz/∂z∂τzz/∂t = (λ+2μ) ∂vz/∂z + λ ∂vx/∂x∂τxz/∂t = μ (∂vx/∂z + ∂vz/∂x)

交错网格的排布规则是：vx 放在半网格点 (i+1/2, j)，vz 放在 (i, j+1/2)，τxx、τzz 放在整格点 (i, j)，τxz 放在 (i+1/2, j+1/2)。这样每个应力对速度的导数都在对应速度点正中央取值，空间精度比普通网格高半阶，也不需要额外插值。

下面是一段可直接运行的 NumPy 核心更新代码。为便于展示，空间导数用前后向差分组合来模拟半格偏移，生产级代码建议严格按交错网格定义做偏移：

import numpy as np # 模型参数（每个网格点一个值），单位统一为 Pa、kg/m^3 lam = np.full((nz, nx), 2.0e9) # 拉梅第一参数 mu = np.full((nz, nx), 2.0e9) # 剪切模量 rho = np.full((nz, nx), 2000.0) # 密度 vx = np.zeros((nz, nx)) # 水平速度 vz = np.zeros((nz, nx)) # 垂直速度 txx = np.zeros((nz, nx)) # 正应力 xx tzz = np.zeros((nz, nx)) # 正应力 zz txz = np.zeros((nz, nx)) # 剪应力 xz dt, dx, dz = config['dt'], config['dx'], config['dz'] for it in range(config['nt']): # 应力 -> 速度，2 阶交错差分近似 vx[1:-1, 1:-1] += (dt / rho[1:-1, 1:-1]) * ( (txx[1:-1, 1:] - txx[1:-1, :-1]) / dx + (txz[1:, 1:-1] - txz[:-1, 1:-1]) / dz ) vz[1:-1, 1:-1] += (dt / rho[1:-1, 1:-1]) * ( (txz[1:-1, 1:] - txz[1:-1, :-1]) / dx + (tzz[1:, 1:-1] - tzz[:-1, 1:-1]) / dz ) # 速度 -> 应力 txx[1:-1, 1:-1] += dt * ( (lam[1:-1, 1:-1] + 2 * mu[1:-1, 1:-1]) * (vx[1:-1, 1:] - vx[1:-1, :-1]) / dx + lam[1:-1, 1:-1] * (vz[1:, 1:-1] - vz[:-1, 1:-1]) / dz ) tzz[1:-1, 1:-1] += dt * ( (lam[1:-1, 1:-1] + 2 * mu[1:-1, 1:-1]) * (vz[1:, 1:-1] - vz[:-1, 1:-1]) / dz + lam[1:-1, 1:-1] * (vx[1:-1, 1:] - vx[1:-1, :-1]) / dx ) txz[1:-1, 1:-1] += dt * mu[1:-1, 1:-1] * ( (vx[1:-1, 1:] - vx[1:-1, :-1]) / dz + (vz[1:, 1:-1] - vz[:-1, 1:-1]) / dx )

代码前半段用应力差分更新速度，后半段用速度差分更新应力，dt 是共同因子。交错网格的关键就在于每个差分方向都落在对应半格点上，(txx[i, j+1] - txx[i, j]) / dx近似的是 txx 在 vx 点 (i, j+1/2) 处的 x 方向导数，正好不需要插值。

有三点需要注意。第一，lam、mu、rho 放在格点上，但 τxz 相位处的 μ 要做谐波平均，直接取算术平均会在高速对比界面上产生寄生振荡。第二，这个版本没有吸收边界，波到边界会被反射，实际计算不能直接裸跑。第三，代码空间精度是 2 阶，生产规模建议升到 4 阶或 8 阶，第 4 章会给出模板。

3.2 Ricker 子波震源：主频、延时和加载方式

震源函数最常用 Ricker 子波，时间域形式：

s(t) = (1 - 2 (π f0 (t - t0))²) · exp(-(π f0 (t - t0))²)

频谱峰值在 f0 处，高频端约在 2.5 f0 处截止。实现只需几行：

def ricker(f0, t, delay=None): if delay is None: delay = 1.0 / f0 p = np.pi * f0 * (np.asarray(t) - delay) return (1.0 - 2.0 * p * p) * np.exp(-p * p)

delay 一般取 1 到 1.5 个主频周期，保证子波在 t=0 时幅值接近零，避免加载瞬间产生硬激励。震源加载分两类：力源和爆炸源。

力源直接把子波乘以方向向量加到速度数组的震源格点上，方向由源的方向决定。爆炸源则是加到正应力上，而且不能把同一个子波同时加到 τxx 和 τzz，否则会破坏各向同性压力平衡，产生不真实的纯膨胀波源：

# 压力源：按 (λ+2μ) 与 λ 的比例分配 src_xx = -s * (lam[sx, sz] + 2.0 * mu[sx, sz]) * scale src_zz = -s * lam[sx, sz] * scale txx[sx, sz] += dt * src_xx tzz[sx, sz] += dt * src_zz

scale 是幅度标定因子，把应力单位换算到目标量级。这个分配保证源处只有膨胀分量、无剪切分量，激发的是纯纵波；如果分配比例写反，波场里会混入明显的人为横波，这是新手最容易看错的地方。

3.3 PML 吸收边界：为什么海绵边界不划算

波场传到模型边界必须被吸收，否则反射波会把合成记录搅乱。最简单的海绵衰减边界在边界外侧加指数衰减带，但想达到 -30dB 衰减通常要 100 层以上，计算冗余太大。成熟方案是 PML（完全匹配层），通过坐标拉伸把波场映射到衰减空间，配合递归卷积实现，一般 10 到 30 层就能衰减 40dB。

实现 PML 要给边界区域配置阻尼剖面。阻尼系数按多项式增长：

α(x) = α0 · (x / L)^p

其中 L 是 PML 层数，p 是增长阶数。α0 由目标反射系数决定，经验公式取：

α0 = -( (p+1) · vmax / (2 · L) ) · ln(R)

R 取 0.001 时，典型参数配置如下：

纯 PML 在大角度掠入射时仍有反射，工程上更多用卷积 PML（CPML）。CPML 的实现要点是在每个空间导数后增加一个记忆变量递归项：

# CPML 记忆变量，以 τxx 对 x 的导数为例 grad = (txx[1:-1, 1:] - txx[1:-1, :-1]) / dx psi_vx = b_x * psi_vx + a_x * grad vx[1:-1, 1:-1] += (dt / rho[1:-1, 1:-1]) * (grad + psi_vx)

更新系数 a、b 随网格位置变化：

b = exp(-(α + α_max) · dt)a = α_max · (b - 1) / (α + α_max)（α + α_max 不为 0 时）

CPML 比经典 PML 多几行代码，但对入射角不敏感，稳定性也更好。检查 PML 是否生效的办法是在模型角落放一个虚检波器，记录波场尾部能量，若尾部能量比主波振幅低 40dB 以上即为合格。

4. 主频、网格间距与 CFL：弹性波正演模拟的三个必调参数

4.1 网格间距由最小横波速度决定：每最小波长 8 到 12 个点

网格间距不能拍脑袋。限制条件来自最小速度：弹性波场里横波速度最低，相同频率下横波波长最短，需要最高分辨率。工程标准是每最小波长至少 8 个网格点，更严格的项目按 10 到 12 个点控制。

设震源主频 f0，有效最高频 fmax ≈ 2.5 f0，则：

h ≤ Vs_min / (fmax · ppw)

ppw 是每最小波长的网格点数，取 8 到 12。看一个典型例子的开销增长：

格数从 250 涨到 2000，是 8 倍关系；时间步长还要再缩 8 倍，迭代步数涨 8 倍，总计算量是 64 倍的关系。这就是主频翻一倍、算力需求涨两个数量级的原因。做三维时这个倍率还要乘一次稀疏度相关项，所以很多三维生产项目把主频压在 10 到 15Hz，是为了让正演模拟在划算的区间内跑完。

4.2 CFL 条件和时间步长：dt 被稳定性约束，不是越小越好

二阶蛙跳格式的稳定性条件是 CFL 条件。二维交错网格下：

dt ≤ h / (Vp_max · sqrt(2))

三维是 sqrt(3) 而不是 sqrt(2)。Vp_max 取模型最大纵波速度，因为纵波速度最快，最先触碰稳定性边界。实际加 0.6 到 0.7 的安全系数：

dt = 0.6 · h / (Vp_max · sqrt(2))

举例：Vp_max = 3000 m/s，h = 5 m，理论上限 dt = 5/(3000·1.414) = 1.18 ms，取 0.6 倍后是 0.71 ms。若模拟时长为 1 秒，需要约 1400 步。

取 0.6 而不是 0.95 的原因有三个：模型速度不是常数，插值后局部波速略高于给定值；PML 区的阻尼项改变局部有效波速；非均匀介质的实际频散关系偏离均匀网格假设。留这个裕度几乎不增加成本，不建议压线运行。

不稳定的典型现象是波场中出现颗粒状高频噪声，振幅指数增长，很快覆盖全模型。排查顺序固定为：

1. 检查 dt 是否满足 CFL，Vp_max 是否取到最大值，包括 PML 内部
2. 检查 PML 层数是否足够，阻尼系数是否算错
3. 检查自由表面处理：真空法把密度设成极小值而不是 0
4. 检查震源加载是否在源点产生超出模板分辨能力的陡峭突变

4.3 数值频散：从波形尾部判断网格是否够细

弹性波有限差分的数值频散表现很有特征：波前后面拖一段高频振荡，振幅不大但很扎眼，像给主波梳了个锯齿头。P 波比 S 波明显，因为 P 波速度高，同样格距下每波长格数更少。

出现频散时，优先提升空间精度而不是盲目加密网格。空间差分从 2 阶提到 4 阶，在相同网格下频散显著减少，计算量只增加约一倍；加密网格到两倍需要 8 倍算力（三维）。一维 4 阶差分模板是：

∂f/∂x ≈ (1/12 · f_{i-2} - 2/3 · f_{i-1} + 2/3 · f_{i+1} - 1/12 · f_{i+2}) / h

边界附近退化为 2 阶或单边格式。确认网格够不够细，最稳妥的验证是收敛性检查：把 h 减半跑同一模型，比较同一接收点波形。两条波形几乎重合说明已收敛；若尾部差异大，说明原网格偏粗，加密是唯一出路。

4.4 参数自检脚本：跑大模型前先算一遍预算

把前三节写成一段快速预算脚本，任何模型换参数之前先跑一遍，能避免把几天算力浪费在错误配置上：

def wave_parameter_report(vs_min, vp_max, f0, nx, nz, nt_sim, ppw=10, cfl=0.6, pml=20, fp32=True): fmax = 2.5 * f0 h = vs_min / (fmax * ppw) dt = cfl * h / (vp_max * (2 ** 0.5)) nt = int(nt_sim / dt) bytes_per_point = 4 if fp32 else 8 core_mem = 5 * (nx + 2 * pml) * (nz + 2 * pml) * bytes_per_point total_mem_gb = core_mem / 1024 ** 3 print(f"建议网格间距 h = {h:.3f} m") print(f"建议时间步长 dt = {dt * 1000:.3f} ms") print(f"模拟 {nt_sim:.3f} s 需要 {nt} 步") print(f"核心数组内存 ≈ {total_mem_gb:.2f} GB (不含 PML 辅助变量)") return h, dt, nt, total_mem_gb wave_parameter_report(vs_min=800, vp_max=3200, f0=10, nx=800, nz=800, nt_sim=1.5)

输出四个最关键的数字：网格间距、时间步长、迭代步数和核心内存。PML 辅助变量、临时数组和输出快照要另算。看到 h 或 dt 不符合预期，回到 4.1 和 4.2 检查模型参数是否写错；内存超出机器上限，就该考虑降主频、换单精度或换并行方案。

5. 把内存和算力花在刀刃上：弹性波正演模拟的实测加速与验证

5.1 数组布局和精度：float32 + 连续内存是底线

用 float32 存储波场数组，误差在 1e-7 量级，对勘探频段的波场计算完全够用，内存和访存带宽直接减半。NumPy 里按 Fortran 连续（order='F'）创建数组，让列方向连续，配合按列更新，缓存命中率比默认 C 顺序高不少。性能差异在三维大网格上能拉到 30% 以上。追求更高吞吐时，把 5 个核心数组合并成一张连续的大数组，避免多次查表。

5.2 快照输出按需落盘，不要每步都写

IO 是波场正演被忽略的性能瓶颈。先算好每秒需要的快照张数，比如 50 张，然后设置输出间隔snap_interval = nt // 50，在时间循环里用取模判断触发。检波器记录则在每一步做插值叠加到合成记录数组里，最后一次性写盘。这样快照只写需要的帧数，合成地震记录不经过中间文件，能省下大量磁盘和写入时间。

5.3 均匀介质解析解对比是最后的守门员

参数配置齐了不代表算对了。最直接的验证是在均匀介质中放一个爆炸源，把数值解和解析解对比。均匀全空间里 P 波到达时间为距离除以 Vp，S 波到达时间为距离除以 Vs。取波形第一个峰的到达时刻对比，偏差应小于一个采样间隔；振幅 RMS 偏差控制在 5% 以内。如果 P 波对的、S 波慢半拍，问题大概率出在拉梅参数到纵横波速度的换算；如果两者都偏早或偏晚，检查 dt 和震源延时设置。

验证通过后再进入复杂模型，才算把弹性波正演模拟的参数体系跑完整。

本文还有配套的精品资源，点击获取