简介：针对地球物理勘探、地震波模拟与计算声学中的声波正演需求，这份资源为研究生、工程师或科研人员提供了一个基于MATLAB的简洁示例，核心聚焦PML完美匹配层边界、高阶有限差分与频散抑制问题。包内仅有1个m脚本文件，rar压缩包整体约2KB，极其精简；代码框架可演示如何设置网格、施加PML吸收边界并执行高阶差分运算，适合已有有限差分基础、正在调试声波模拟程序的读者参考。目前该资源已有153人学习。借助示例，读者可直观对比低阶与高阶差分在数值频散控制上的差异，理解PML在开放边界截断中的作用，并可将代码骨架迁移至自己的正演模拟或反演成像流程中，节省从零搭建算法的时间。由于代码紧凑，也便于逐行阅读和修改，用于教学演示或快速原型验证。

1. 声波有限差分模拟：频散和PML边界是两个绕不过去的坎

做声波有限差分模拟时，最容易翻车的不是震源设置，而是波场图画出来之后，波形像被风吹散一样，高频在前面跑、低频在后面拖，同一个波越传越宽；如果边界再不做吸收处理，网格四周的硬反射还会把有效波彻底搅乱。这套 shengbo.rar 资源给了一条完整的声波模拟链路：PML 边界吸收、高阶差分压制频散、从震源加载到波场快照都有现成代码。它适合正在做声波正演、准备波场特征分析，或者需要拿正演内核去喂反演算法的人。新手能直接跑通，熟手可以拿它做参数对照实验的底座。

2. 波动方程离散化：二阶差分为何产生频散，PML层如何吸收边界反射

2.1 二阶声波方程与速度-应力一阶格式

声波模拟的核心是求解二阶波动方程：

∂²p/∂t² = c² ∇²p

其中 p 是压力，c 是介质声速。如果直接对时间和空间做二阶中心差分，很容易得到一个有明显数值频散的格式。原因在于二阶中心差分算子的相速度与波数有关，波数越大（对应波长越短），相速度偏离真实值越多，于是短波长成分越跑越慢，在波场图上表现为波形边部出现振荡式拖尾。

工程上我一般先把二阶方程改写成速度-应力一阶方程组：

∂v/∂t = -(1/ρ) ∂p/∂x
∂p/∂t = -ρc² ∂v/∂x

这里 v 是质点振速。一阶格式的好处是边界条件处理更直观，特别是做 PML 吸收层时，衰减项可以直接挂在速度更新方程上。资源包里采用的正是这种速度-应力一阶格式，和很多商业正演软件的内部思路一致。下面是网格与场量初始化的常见写法：

% 声波模拟基础参数定义 nx = 600; nz = 600; % 网格数 dx = 5; dz = 5; % 网格步长，单位 m c = 2000; % 声速，m/s rho = 1000; % 密度，kg/m^3 nt = 1200; % 时间步数 dt = 0.001; % 时间步长，s % 场量初始化 vx = zeros(nx, nz); % x 方向质点速度 vz = zeros(nx, nz); % z 方向质点速度 p = zeros(nx, nz); % 压力场

为什么先定义速度再定义压力？因为时间推进时，速度更新依赖压力的空间导数，压力更新又依赖速度的空间导数，两者交替迭代。记住 p 和 v 的时间步差半个步长，这在交错网格实现中尤其重要，资源包内已经处理好了这一步，自己重写时容易忽略。

2.2 高阶空间差分：系数怎么算，频散怎么降

空间导数的差分阶数直接决定频散程度。二阶中心差分相当于用相邻三点拟合一阶导，误差项含 Δx²。四阶差分用五点模板，能显著压低相速度误差。交错网格下，一阶导数的四阶近似为：

∂f/∂x ≈ (1/Δx)[ (1/24)f(i-3/2) - (9/8)f(i-1/2) + (9/8)f(i+1/2) - (1/24)f(i+3/2) ]

资源包里通常提供二阶、四阶、八阶三档可选。我一般建议：常规正演用四阶，PML 区域内外保持一致；做高精度基准解时用八阶。高阶差分系数可以通过泰勒展开统一求解，下面给出一个通用生成函数：

function coef = diff_coef(N) % N 为半模板长度：N=1 对应二阶，N=2 对应四阶，N=4 对应八阶 % 求一阶导数中心差分系数，在 x=0 处恰等于导数 L = 2*N - 1; % 模板长度 A = zeros(L, L); for i = 1:L x = (i - N); for j = 1:L A(i, j) = x^(j-1); end end b = zeros(L, 1); b(2) = 1; % 目标：p'(0) = 1 coef = A\b; end

这段代码的思路很直白：用一个长度 L 的差分模板近似导数，要求它对常数项、一次项、二次项直到 L-1 次项都有精确的加权效果，于是每个目标导数对应一个线性方程。这里的参数 N 是半模板长度，真正使用时要把系数按交错网格错半个点处理。得到的 coef 就是各节点的权重，更新波场时把这些系数点乘到压力或速度的邻域上。

2.3 PML吸收边界：复坐标拉伸与辅助方程

PML 的基本原理是在计算区域外增加数层吸收层，通过复坐标拉伸把传播波变成衰减波。频域里空间导数变为：

∂/∂x → (1/(1 + σ_x(x)/iω)) ∂/∂x

其中 σ_x 在吸收层内从零逐渐增大，吸收层外为零。时域实现通常有两种做法：一种是在速度更新方程中加入与 σ 相关的记忆变量；另一种是直接在每一步把全场速度乘以 exp(-σdt)。资源包采用的是辅助微分方程方式，比简单乘指数更稳定，尤其适合长时间模拟。

σ_max 的经验公式为：

σ_max = -(3c)/(2Δx) ln(R)

其中 R 是目标反射系数，通常取 0.001。PML 层内 σ 按距离平方渐变，这样可以减少层内数值反射。资源包里的 make_pml 函数大致如下：

function sig = make_pml(ncells, pml_len, c, dx, R) % 生成一维 PML 的 sigma 系数 sig = zeros(1, ncells); sigma_max = -(3*c/(2*dx)) * log(R); for i = 1:pml_len d = (pml_len - i) / pml_len; % 从内边界到外边界归一化距离 sig(i) = sigma_max * d^2; % 平方渐变 end % 注意这里只给了左边 PML，右边和两端的 PML 需要按镜像方式补全 end

使用这个系数时，注意它只作用在 PML 区域内的速度更新项里，不能对内部区域做任何改动。资源包里把 P 和 V 的更新都封装成了独立函数，PML 区域和内部区域共用一个差分算子，只额外乘上衰减项。这一步做对，边界反射基本肉眼不可见。

3. 解压 shengbo.rar 到跑出第一张波场图：文件结构与复现步骤

3.1 包内文件怎么分工

解压后一般是几类文件。常见布局如下表，实际文件名可能略有差异，但职责分工是一致的：

拿到包之后，先不要急着跑最大网格。我习惯先把 nx,nz 改成 200×200，pml_len 改成 12，震源主频调到 20 Hz，跑一个几十步看看波形形态是否正常，确认整个链路没有断掉，再放大网格。

3.2 最小可运行示例

资源包里主程序封装得比较好，核心循环一般是这样：

% shengbo_mini.m —— 最小可运行声波模拟 clear; clc; % 模型参数 nx = 400; nz = 400; dx = 5; dz = 5; dt = 0.001; nt = 800; c = 2000; rho = 1000; % PML 参数 pml_len = 20; R = 0.001; sigx = make_pml(nx, pml_len, c, dx, R); % 需要左右拼接 sigz = make_pml(nz, pml_len, c, dz, R); % 场量 p = zeros(nx, nz); vx = zeros(nx, nz); vz = zeros(nx, nz); % 震源位置与主频 sx = round(nx/2); sz = round(nz/2); f0 = 30; % Hz % 接收点位置 rec = round(linspace(10, nx-10, 20)); seismo = zeros(nt, length(rec)); % 主循环 for it = 1:nt % 雷克子波加载到压力场 tau = (it-1)*dt - 1/f0; src = (1 - 2*(pi*f0*tau)^2) * exp(-(pi*f0*tau)^2); p(sx, sz) = p(sx, sz) + src; % 压力空间梯度（简化中心差分，完整版在包内） [dpdx, dpdz] = gradient(p, dx, dz); % 速度更新 vx = vx - (dt/rho) * dpdx; vz = vz - (dt/rho) * dpdz; % 压力更新（利用速度散度） dvxdx = gradient(vx, dx, dz); dvzdz = gradient(vz, dx, dz); p = p - rho*c^2*dt*(dvxdx + dvzdz); % 记录接收点 seismo(it, :) = p(rec, 1); % 每 200 步画一张快照 if mod(it, 200) == 0 imagesc(p'); axis equal; colorbar; title(sprintf('time %.3f s', it*dt)); drawnow; end end

这段代码里，gradient 是 MATLAB 内置的二阶精度梯度算子，所以只能用来演示流程。资源包里的正式版本针对速度-应力交错网格做了索引对齐，差分精度更高，直接运行不会出现这里可能有的半个网格偏移。把 pml_len、f0、dx 这三个参数理解成入口，其他都是内部细节。

3.3 波场快照与地震记录输出

运行完成后，你手里会得到 seismo 矩阵，每一列就是一个接收点的压力时间序列。保存结果时，我一般用：

save('sim_result.mat', 'seismo', 'dx', 'dt', 'c', 'rho', 'f0');

这样后续做频散分析、PML 吸收效果验证时，可以直接读这个 mat 文件，不需要重新跑一次长计算。波场快照如果不想实时画，可以每个 it 保存一个优化后的 uint8 崔值，最后统一合成视频，好处是快速回放时不会因为重算拖慢速度。

4. 参数选择的四个关键旋钮：网格步长、时间步长、差分阶数与PML层数

4.1 网格步长：频散的硬约束

网格步长是压制频散的第一道闸门。经验准则是：每个最小波长至少有 5～8 个网格点。最小波长由介质最小速度和震源最高有效频率决定：

dx ≤ c_min / (8 * f_max)

其中 f_max 通常取主频的 2～3 倍。比如 c_min=2000 m/s，f0=30 Hz，按 f_max=90 Hz 算，dx 应该不超过 2.78 m。资源包默认 dx=5 m，配合 30 Hz 主频其实偏粗，所以必须用高阶差分来补偿。如果你用的是二阶差分，那这个 dx 下频散会非常明显，波形尾部会有长一的振荡尾巴。

4.2 时间步长与 CFL 稳定性

时间步长受 CFL 条件约束。二维速度-应力格式的稳定条件是：

c * dt * sqrt(1/dx² + 1/dz²) ≤ 1

若 dx=dz=h，则 cdt/h ≤ 1/√2。实际工程中我会再留余量，取 cdt/h 在 0.4～0.6 之间。资源包默认 dt=0.001，dx=5，c=2000，算出来 c*dt/h=0.4，是一个保守但安全的组合。如果你缩小网格步长，一定要同步缩小时间步长，否则高频振荡可能当场发散成 NaN。

4.3 差分阶数：用计算量换频散精度

不同差分阶数的空间模板长度和频散表现有明显差异：

资源包里一般有一个全局参数 order，设成 2 或 4。我建议：如果只是看波场形态、验证网格边界，四阶足够；如果要做定量频散误差分析，建议跑一版八阶作为参考解。切换阶数时要注意，PML 区域内的差分系数也必须同步切换，否则 PML 内部和外部网格阻抗不匹配，会产生新的数值反射。

4.4 PML层数与吸收系数

PML 层数太少，吸收曲线还没衰减到足够低就碰到外边界，能量会反弹回来；层数太多，计算区域被压缩。我习惯取 12～20 层。R 取 0.001 到 0.0001。R 越小，σ_max 越大，理论上吸收越强，但 σ_max 过大会在 PML 层内产生离散反射，所以不建议 R 小于 0.0001。修改 PML 参数后，记得检查 make_pml 返回的 σ 曲线是否单调递增，如果不是单调的，说明 pml_len 与 sigma_max 的匹配出问题了。

5. 避坑与排查：我在这套模拟上翻过的五个车

5.1 现象：波场在计算区域边缘出现明显“硬反射”

原因：PML 参数没生效，或者 PML 层内场量更新只改了响应，却没有在速度更新里加入衰减项。最常见的是 make_pml 函数只生成了左边界系数，右边界和上下边界都是零。

解决：print 出 sigx 和 sigz，检查四个方向 PML 区域的系数是否非零；确认速度更新循环里，对 PML 内部的 vx/vz 额外乘了 exp(-σ*dt) 或者加入了记忆变量项。我通常会写一个自检函数，统计边界前 5 层的 σ 总和，非零才继续运行。

5.2 现象：震源附近高频振荡发散，几步后变成 NaN

原因：时间步长不满足 CFL 条件，或者雷克子波主频太高导致最小波长只占两个网格点。

解决：先算 c*dt/h，如果大于 0.7 就减小 dt；再把 f0 降到 20 Hz 试一次。资源包自带了一个 check_stability 函数，会在运行时打印当前 CFL 数。遇到 NaN 还有一招：把 dt 改成原值的一半，看是否立刻正常，这能快速定位是稳定性问题而不是算法 bug。

5.3 现象：波形随时间逐渐“糊开”，单频波尾部振荡明显

原因：空间差分阶数太低，加上网格步长相对波长太粗。二阶差分对 30 Hz 主频、5 m 步长的模拟，频散误差累积很快。

解决：把差分阶数从 2 换到 4，同时检查 dx 是否满足第 4.1 节的波长准则。如果还想保留原有 dx，就把 f0 降低到 15～20 Hz，用低频源压制频散。我也遇到过用户直接上八阶差分但 dx 依然 20 m 的情况，这时频散确实小了，但计算量也上去了，效果未必比四阶加小步长更好。

5.4 现象：PML 层内出现条状强振幅，像“第二震源”

原因：PML 区域的空间差分使用了单边差分或与内部不一致的差分算子，导致阻抗突变。另外，如果直接把压力场 p 在边界处乘衰减因子，而不是在速度更新中加入衰减，也会产生这种假源。

解决：确认 PML 内外使用同一个 diff_coef 生成的差分系数；衰减项只作用在速度场，不直接乘 p。检查 make_pml 的渐变函数是不是 d²，用 d 的线性渐变虽然简单，但在层内会产生更多反射。

5.5 现象：网格稍大（800×800）就跑得很慢，几分钟才一步

原因：代码里用了双层 for 循环，或者每步都在画高分辨率 imagesc 图。MATLAB 在循环里逐点更新压力场，速度极慢。

解决：把内层循环向量化成矩阵运算。资源包一般已经向量化，但如果你在复现时是自己写的循环，可以对照以下写法：

% 向量化速度更新示例（去掉 PML 项后） vx(2:end-1,:) = vx(2:end-1,:) - (dt/rho) * (p(3:end,:) - p(1:end-2,:)) / (2*dx); vz(:,2:end-1) = vz(:,2:end-1) - (dt/rho) * (p(:,3:end) - p(:,1:end-2)) / (2*dz);

这比 for 循环快两个数量级。另外把画图间隔从每一幅改为每 200 步，或者直接不开 figure 窗口，只保存快照文件，速度立刻改善显著。

6. 进阶用法：用同一份代码做定量频散和 PML 吸收效果验证

6.1 扫描差分阶数，比较频散差异

搞清楚了基础流程，这套资源最值得做的事就是用它定量对比二阶、四阶、八阶的频散误差。固定 dx=5、dt=0.001、f0=30 Hz、模型大小不变，分别跑三个阶数，在同一个接收点记录波形：

orders = [2 4 8]; figure; hold on; for k = 1:3 config.order = orders(k); % 调用资源包内封装好的传播函数，返回接收点压力记录 [seismo, t] = run_simulation(config); % 归一化并做轻微纵移，方便比较 seismo = seismo / max(abs(seismo)); plot(t, seismo + (k-1)*1.5, 'LineWidth', 1.2); end legend('2阶','4阶','8阶'); xlabel('时间 s'); ylabel('归一化压力');

运行后你会看到：二阶结果的波前面明显更“软”，到达时刻偏晚，尾部振荡多；八阶的波前面陡峭，到达时刻更接近理论值。这个对比图贴在报告里很有说服力，也能帮自己判断当前模型参数下差分阶数取舍。

6.2 用能量衰减曲线验证 PML 质量

验证 PML 吸收效果，不用靠肉眼，我常用两个探针：一个放在内部区域中心，一个放在 PML 内边界内侧 5 格处。记录两个位置的压力均方根随时间变化。如果 PML 工作正常，边缘探针的能量应随传播快速衰减；若出现明显台阶或回升，说明边界反射折返进来了。

% 记录 PML 内侧一点的能量包络 probe_x = pml_len + 5; probe_z = round(nz/2); p_energy = sqrt(mean(p(probe_x, :).^2, 2)); % 每次时间步取均方根 % 理想情况下 p_energy 应随时间单调下降，若出现上扬尖峰，则检查 PML 层数

从那以后，我每次跑这套正演，都会先花两分钟把差分阶数和 PML 层数打印出来，顺手画一条能量衰减曲线，确认边界不是默认无效参数。搞明白这两个指标，声波模拟的波场图才算真正可信，希望帮到你。

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