简介弹性波方程正演模拟是地震学与地球物理勘探中的基础环节这份压缩包面向相关专业学生与科研人员提供基于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 时间递推用蛙跳格式内存最少一阶速度-应力方程在时间方向上的标准做法是蛙跳格式速度更新到 n1/2 时刻应力更新到 n1 时刻交替推进。更新式写作v^(n1/2) v^(n-1/2) Δt/ρ · D τ^n Δt/ρ · f^nτ^(n1) τ^n Δt · C : D v^(n1/2) Δt · M^(n1/2)D 是空间差分算子C 是弹性刚度张量。蛙跳格式需要保存的速度和应力都只有一份历史数据内存开销最小计算量也最小。它在一阶双曲型问题上的稳定性由 CFL 条件约束第 4 章会具体计算。有人倾向用四阶 Runge-Kutta 做时间积分认为精度更高。实际上在同样的网格配置下RK4 能把时间步长放宽 2 到 3 倍但每步要做 4 轮空间导数临时数组至少翻一倍综合算力往往比蛙跳多 50% 以上。对绝大多数勘探频率范围的弹性波正演模拟二阶蛙跳已经足够只有做长时间波场传播、需要严格控制时间频散的研究型代码才会换高阶时间格式。2.3 交错网格有限差分、有限元、谱元怎么选空间离散方法的选择直接影响网格规模和生产效率。把四类方法的定位放在一张表里看更清楚方法空间精度典型配置网格需求适合场景主要缺点交错网格有限差分2 到 8 阶每最小波长 8 到 12 个点大规模三维、快速迭代、各向异性程度不极端自由表面需真空法或镜像法处理边界参数要调有限元低阶到高阶网格尺寸与波长比约 1/8复杂地形、强不规则界面需要组装和求解大规模线性系统谱元法高斯-洛巴托点 N4 到 8每最小波长 3 到 5 个点高保真强散射模拟、谱元-有限元耦合实现复杂度高曲线网格生成难伪谱法FFT 全局算子每波长 2 到 3 个点均匀背景大尺度波场介质剧烈变化时产生吉布斯振荡从表里能看出交错网格有限差分在每波长网格数上不是最省的但胜在模型参数直接落在格点上、无需网格剖分复杂程度可控。后面的实现都围绕它展开。3. 用交错网格有限差分把弹性波方程跑起来时间递推、震源与 PML3.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 放在半网格点 (i1/2, j)vz 放在 (i, j1/2)τxx、τzz 放在整格点 (i, j)τxz 放在 (i1/2, j1/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, j1] - txx[i, j]) / dx近似的是 txx 在 vx 点 (i, j1/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, delayNone): 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 个主频周期保证子波在 t0 时幅值接近零避免加载瞬间产生硬激励。震源加载分两类力源和爆炸源。力源直接把子波乘以方向向量加到速度数组的震源格点上方向由源的方向决定。爆炸源则是加到正应力上而且不能把同一个子波同时加到 τ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_zzscale 是幅度标定因子把应力单位换算到目标量级。这个分配保证源处只有膨胀分量、无剪切分量激发的是纯纵波如果分配比例写反波场里会混入明显的人为横波这是新手最容易看错的地方。3.3 PML 吸收边界为什么海绵边界不划算波场传到模型边界必须被吸收否则反射波会把合成记录搅乱。最简单的海绵衰减边界在边界外侧加指数衰减带但想达到 -30dB 衰减通常要 100 层以上计算冗余太大。成熟方案是 PML完全匹配层通过坐标拉伸把波场映射到衰减空间配合递归卷积实现一般 10 到 30 层就能衰减 40dB。实现 PML 要给边界区域配置阻尼剖面。阻尼系数按多项式增长α(x) α0 · (x / L)^p其中 L 是 PML 层数p 是增长阶数。α0 由目标反射系数决定经验公式取α0 -( (p1) · vmax / (2 · L) ) · ln(R)R 取 0.001 时典型参数配置如下参数建议范围对波场的影响PML 层数 L15 到 30太小会穿透反射太大内侧阻尼陡增大角度入射时反射增长阶数 p2 到 3p 越大衰减越快但过陡会导致阻抗渐变不足目标反射系数 R1e-3 左右越小越强吸收但要更多层纯 PML 在大角度掠入射时仍有反射工程上更多用卷积 PMLCPML。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。看一个典型例子的开销增长f0fmax2.5f0Vs_min1000 m/s 时的 h2000 m 模型格数5 Hz12.5 Hz8 m按 10 点25010 Hz25 Hz4 m50020 Hz50 Hz2 m100040 Hz100 Hz1 m2000格数从 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/sh 5 m理论上限 dt 5/(3000·1.414) 1.18 ms取 0.6 倍后是 0.71 ms。若模拟时长为 1 秒需要约 1400 步。取 0.6 而不是 0.95 的原因有三个模型速度不是常数插值后局部波速略高于给定值PML 区的阻尼项改变局部有效波速非均匀介质的实际频散关系偏离均匀网格假设。留这个裕度几乎不增加成本不建议压线运行。不稳定的典型现象是波场中出现颗粒状高频噪声振幅指数增长很快覆盖全模型。排查顺序固定为检查 dt 是否满足 CFLVp_max 是否取到最大值包括 PML 内部检查 PML 层数是否足够阻尼系数是否算错检查自由表面处理真空法把密度设成极小值而不是 0检查震源加载是否在源点产生超出模板分辨能力的陡峭突变4.3 数值频散从波形尾部判断网格是否够细弹性波有限差分的数值频散表现很有特征波前后面拖一段高频振荡振幅不大但很扎眼像给主波梳了个锯齿头。P 波比 S 波明显因为 P 波速度高同样格距下每波长格数更少。出现频散时优先提升空间精度而不是盲目加密网格。空间差分从 2 阶提到 4 阶在相同网格下频散显著减少计算量只增加约一倍加密网格到两倍需要 8 倍算力三维。一维 4 阶差分模板是∂f/∂x ≈ (1/12 · f_{i-2} - 2/3 · f_{i-1} 2/3 · f_{i1} - 1/12 · f_{i2}) / h边界附近退化为 2 阶或单边格式。确认网格够不够细最稳妥的验证是收敛性检查把 h 减半跑同一模型比较同一接收点波形。两条波形几乎重合说明已收敛若尾部差异大说明原网格偏粗加密是唯一出路。4.4 参数自检脚本跑大模型前先算一遍预算把前三节写成一段快速预算脚本任何模型换参数之前先跑一遍能避免把几天算力浪费在错误配置上def wave_parameter_report(vs_min, vp_max, f0, nx, nz, nt_sim, ppw10, cfl0.6, pml20, fp32True): 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_min800, vp_max3200, f010, nx800, nz800, nt_sim1.5)输出四个最关键的数字网格间距、时间步长、迭代步数和核心内存。PML 辅助变量、临时数组和输出快照要另算。看到 h 或 dt 不符合预期回到 4.1 和 4.2 检查模型参数是否写错内存超出机器上限就该考虑降主频、换单精度或换并行方案。5. 把内存和算力花在刀刃上弹性波正演模拟的实测加速与验证5.1 数组布局和精度float32 连续内存是底线用 float32 存储波场数组误差在 1e-7 量级对勘探频段的波场计算完全够用内存和访存带宽直接减半。NumPy 里按 Fortran 连续orderF创建数组让列方向连续配合按列更新缓存命中率比默认 C 顺序高不少。性能差异在三维大网格上能拉到 30% 以上。追求更高吞吐时把 5 个核心数组合并成一张连续的大数组避免多次查表。5.2 快照输出按需落盘不要每步都写IO 是波场正演被忽略的性能瓶颈。先算好每秒需要的快照张数比如 50 张然后设置输出间隔snap_interval nt // 50在时间循环里用取模判断触发。检波器记录则在每一步做插值叠加到合成记录数组里最后一次性写盘。这样快照只写需要的帧数合成地震记录不经过中间文件能省下大量磁盘和写入时间。5.3 均匀介质解析解对比是最后的守门员参数配置齐了不代表算对了。最直接的验证是在均匀介质中放一个爆炸源把数值解和解析解对比。均匀全空间里 P 波到达时间为距离除以 VpS 波到达时间为距离除以 Vs。取波形第一个峰的到达时刻对比偏差应小于一个采样间隔振幅 RMS 偏差控制在 5% 以内。如果 P 波对的、S 波慢半拍问题大概率出在拉梅参数到纵横波速度的换算如果两者都偏早或偏晚检查 dt 和震源延时设置。验证通过后再进入复杂模型才算把弹性波正演模拟的参数体系跑完整。本文还有配套的精品资源点击获取