简介：这份资源面向地球物理勘探、信号处理与计算物理方向的学习者和研究人员，聚焦声波数值模拟中的PML边界处理与高阶有限差分算法，用于缓解有限差分求解波动方程时的数值频散问题，提升地下介质中声波传播模拟的精度。压缩包内共1个文件，为m格式的MATLAB脚本，整体约2KB，代码结构紧凑，便于直接阅读与二次修改。内容涵盖网格设置、PML吸收边界条件定义、高阶差分算子实现以及模拟结果可视化等关键环节，读者可据此理解开放边界下反射抑制与频散控制的实现思路，并对照自身模型调试参数、验证算法效果。目前已有153人学习下载，适合具备MATLAB基础、希望深入掌握声波正演模拟技术的工程师与研究人员参考。
1. 从 shengbo.rar 说起：PML 边界、高阶差分与频散到底卡在哪
很多人第一次拿到shengbo.rar这类声波有限差分模拟包，跑通一炮记录就以为完事了，结果把震源挪到模型边缘，或者把模拟时长拉长，波场里立刻冒出莫名其妙的假反射和尾巴——这就是 PML 边界和频散在收“学费”。声波模拟的核心链路其实就三件事：用高阶差分把波动方程离散化、用 PML 吸收边界把外向波吃掉、用足够密的网格和合适的差分阶数把数值频散压到可接受范围。这三者互相牵制：网格越粗越省算力，但频散越重；差分阶数越高越精确，但边界和稳定性条件更苛刻。这篇笔记面向做地震波正演、声学成像、超声仿真的一线同行，把shengbo.rar背后这套方案拆成能直接复现的步骤，参数怎么设、坑在哪，一次讲清。
2. 声波方程高阶差分：从二阶到八阶的离散化选择
2.1 为什么高阶差分能压频散
声波方程里空间二阶导数的差分近似，本质是用有限个网格点的加权和去逼近连续导数。二阶中心差分只用到相邻两个点，它的截断误差是 $O(\Delta x^2)$，在每波长网格点数偏少时，相速度会明显偏离真实速度，表现为短波长成分跑得快或慢，这就是数值频散。高阶差分把模板拉长到 2N 个点，截断误差降到 $O(\Delta x^{2N})$，同样的网格密度下频散误差小一个量级。
我一般用 Taylor 展开系数表来生成差分系数，而不是手抄。下面这段 Python 直接算出 2/4/6/8 阶的系数，方便你换阶数时不用查表：
import numpy as np from math import factorial def fd_coeffs(order): """order: 2,4,6,8 空间差分阶数，返回半侧系数""" N = order // 2 A = np.zeros((N, N)) b = np.zeros(N) for i in range(N): for j in range(N): A[i, j] = (2 * (j + 1) - 1) ** (2 * i + 1) b[i] = 1.0 if i == 0 else 0.0 c = np.linalg.solve(A, b) return c # c[0] 对应最近邻点系数 for o in [2, 4, 6, 8]: print(o, fd_coeffs(o))逻辑说明：方程来自偶数阶导数的 Taylor 展开，只保留奇数阶导数项，解出的c是半侧系数，实际模板是[-c[::-1], 0, c]。参数说明：order越大，模板越长，单点计算量线性增加，但允许的每波长网格点数下降。经验值是每波长至少 5 个点用二阶差分就够，但如果你要模拟长距离传播，建议 6 阶以上，把每波长点数压到 3～4 个还能接受。
2.2 时间推进与稳定性条件
空间离散完，时间上用二阶显式差分，得到经典的蛙跳格式。稳定性由 CFL 条件约束：
$$ \Delta t \le \frac{2}{\pi v_{max}} \cdot \frac{1}{\sqrt{\sum |c_i|}} \cdot \Delta x $$
实际写代码时我习惯留 0.8 的安全系数。下面是最小可跑的时间步进核心：
# 假设 p_prev, p_curr 为前两时刻波场，c 为差分系数 def step(p_prev, p_curr, v, c, dx, dt): lap = np.zeros_like(p_curr) for i, ci in enumerate(c, start=1): lap += ci * (np.roll(p_curr, i, axis=0) + np.roll(p_curr, -i, axis=0) + np.roll(p_curr, i, axis=1) + np.roll(p_curr, -i, axis=1)) lap /= dx * dx p_next = 2 * p_curr - p_prev + (v * dt) ** 2 * lap return p_next逻辑说明：np.roll实现周期边界，正式跑之前要换成 PML 区域。参数说明：dx是网格间距，dt由 CFL 决定，v是速度模型。注意v*dt/dx这个比值直接决定稳定性，超过 0.5 左右就会发散，别问我怎么知道的。
3. PML 吸收边界：把假反射压到 -60 dB 以下
3.1 PML 的拉伸坐标思路
PML 不是简单地把边界值置零，而是在边界区域引入复数坐标拉伸，让外向波在进入 PML 后迅速衰减。实现上通常用分裂场法或辅助微分方程法。我常用的是把空间导数除以一个拉伸因子：
$$ \frac{\partial}{\partial x} \rightarrow \frac{1}{s_x}\frac{\partial}{\partial x}, \quad s_x = 1 + \frac{\sigma_x}{i\omega} $$
在时域里这等价于引入记忆变量。下面是一个简化的 PML 衰减因子设置：
def build_pml(nx, nz, npml=20, alpha=0.001, vmax=3000): """返回 x/z 方向的衰减剖面 sigma""" sigma_x = np.zeros(nx) sigma_z = np.zeros(nz) for i in range(npml): d = (npml - i) / npml sigma_x[i] = alpha * vmax * d ** 2 / npml sigma_x[-i-1] = sigma_x[i] sigma_z[i] = alpha * vmax * d ** 2 / npml sigma_z[-i-1] = sigma_z[i] return sigma_x, sigma_z逻辑说明：d是归一化距离，平方衰减是多项式吸收的常见形式，alpha控制最大衰减强度。参数说明：npml一般取 10～30 层，太薄吸收不干净，太厚浪费算力；alpha在 0.001～0.01 之间调，太大反而会因为阻抗突变产生反射。
3.2 边界反射的验证方法
PML 调没调好，不能靠眼睛看。我一般做一个均匀介质点源实验，把接收器放在边界内侧，看尾波能量。如果 PML 正常，边界反射应该低于主波 -60 dB。下面这段代码算反射能量比：
def pml_reflection_ratio(trace, t0_idx, t1_idx): """t0_idx 主波到达时刻，t1_idx 之后为尾波窗口""" main = np.max(np.abs(trace[t0_idx:t0_idx+50])) tail = np.max(np.abs(trace[t1_idx:t1_idx+200])) return 20 * np.log10(tail / main)逻辑说明：主波窗口取峰值附近，尾波窗口取主波之后足够远的位置。参数说明：如果比值高于 -40 dB，说明 PML 层数不够或alpha偏小；如果出现高频振荡，可能是alpha太大导致数值不稳定。这个验证步骤别省，我见过太多人直接上模型，结果边界反射和真实反射混在一起，解释的时候全乱套。
4. 频散控制：网格、阶数、震源频率的三角平衡
4.1 每波长网格点数的经验下限
频散的本质是不同频率成分的相速度不一致。控制频散最直接的办法是保证每波长有足够多的网格点。经验公式：
$$ G = \frac{v_{min}}{f_{max} \cdot \Delta x} $$
$G$ 是每波长网格点数。二阶差分要求 $G \ge 10$，四阶 $G \ge 6$，六阶 $G \ge 4$，八阶 $G \ge 3$。低于这个值，频散尾巴会明显到影响走时拾取。下面是一个快速检查脚本：
def check_dispersion(vmin, fmax, dx, order): G = vmin / (fmax * dx) limits = {2: 10, 4: 6, 6: 4, 8: 3} ok = G >= limits[order] print(f"G={G:.2f}, 要求>={limits[order]}, {'通过' if ok else '频散风险高'}") return ok逻辑说明：vmin取模型最低速度，fmax取震源最高有效频率。参数说明：如果检查不通过，优先减小dx或提高差分阶数，而不是降fmax，因为降频率会损失分辨率。
4.2 震源函数与频散的关系
震源主频越高，频散越容易暴露。我一般用 Ricker 子波，主频f0和最高有效频率的关系约为fmax ≈ 2.5 * f0。设置震源时把f0和网格参数一起算：
def ricker(f0, dt, nt, delay=1.5/f0): t = np.arange(nt) * dt - delay return (1 - 2 * (np.pi * f0 * t) ** 2) * np.exp(-(np.pi * f0 * t) ** 2)逻辑说明：delay保证子波从零平稳起步。参数说明：f0选完后用fmax = 2.5*f0去跑上面的频散检查。如果模型速度变化大，用vmin检查，别用平均速度，否则低速层里频散会翻车。
5. 避坑与排查：PML 边界和频散最常见的 5 个翻车现场
5.1 边界处波场突然归零，出现强反射
现象：波传到模型边缘时，波场被硬截断，产生和真实反射同量级的假反射。原因：PML 层没生效，或者 PML 区域内的速度模型和内部不一致。解决：检查 PML 区域的更新方程是否用了拉伸因子，确认sigma剖面在边界处平滑过渡，别在 PML 内放速度异常体。
5.2 模拟跑几百步后能量指数增长
现象：波场幅值越来越大，最后溢出。原因：CFL 条件不满足，dt太大。解决：按dt = 0.8 * dx / (vmax * sqrt(sum(|c|)))重算，vmax取模型最大速度，别用平均速度。如果速度对比度大，vmax可能比你想象的高很多。
5.3 频散尾巴拖得很长，走时拾取偏晚
现象：子波后面跟着一串振荡，越远越明显。原因：每波长网格点数不够，或者差分阶数偏低。解决：用check_dispersion检查G值，不达标就减小dx或升阶。注意升阶会增加计算量，先试减小dx，因为dx减半计算量涨 4 倍，但升阶只涨模板长度。
5.4 PML 吸收效果随入射角变差
现象：大角度入射的波在 PML 里反射明显。原因：多项式吸收对掠射角效果差。解决：把alpha调大一点，或者改用卷积 PML（CPML）。CPML 实现复杂一些，但大角度吸收能再压 10～20 dB。如果模拟区域不大，直接加宽 PML 层数也能缓解。
5.5 换模型后结果对不上，但参数没改
现象：同样的差分阶数和 PML 参数，换个速度模型就出问题。原因：vmax变了，CFL 和 PML 的alpha都依赖vmax。解决：每次换模型重新算dt和sigma，别把参数写死在代码里。我习惯把vmax、f0、dx作为输入参数，自动推导dt和 PML 剖面。
6. 进阶技巧：用波数域检验频散曲线，把参数调优变成可量化的事
前面讲的都是经验规则，但真要说服自己或审稿人，得拿出频散曲线。做法很简单：在均匀介质里放一个点源，沿一条测线记录波场，做二维 FFT 得到频率-波数谱，再和理论频散关系对比。理论相速度是常数，数值相速度偏离多少，一眼就能看出来。
import numpy as np def dispersion_curve(traces, dx, dt): """traces: (nt, nx) 波场记录""" nt, nx = traces.shape F = np.fft.fft2(traces) F = np.fft.fftshift(F) k = np.fft.fftshift(np.fft.fftfreq(nx, d=dx)) * 2 * np.pi f = np.fft.fftshift(np.fft.fftfreq(nt, d=dt)) return f, k, np.abs(F) def phase_velocity_error(f, k, v_true): """在能量脊上取点算相速度误差""" # 简化：取每个频率最大能量对应的 k idx = np.argmax(np.abs(F), axis=1) k_peak = k[idx] v_num = 2 * np.pi * f / (k_peak + 1e-12) return np.abs(v_num - v_true) / v_true逻辑说明：fft2把时空域转到频率-波数域，能量脊对应频散曲线。参数说明：dx和dt要和模拟一致，v_true是均匀介质速度。如果误差在有效频带内超过 1%，说明频散控制不够，回去调dx或阶数。
我自己的习惯是：每换一个速度模型，先跑一个均匀介质测试，把频散曲线和 PML 反射比都过一遍，再上真实模型。这两个检查加起来不到十分钟，但能省掉后面几小时的返工。shengbo.rar这类包给的是骨架，参数和验证得自己补上，别指望默认配置能适配所有场景。希望帮到你。
本文还有配套的精品资源，点击获取
