简介针对地球物理勘探、地震波模拟与计算声学中的声波正演需求这份资源为研究生、工程师或科研人员提供了一个基于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(i1/2) - (1/24)f(i3/2) ]资源包里通常提供二阶、四阶、八阶三档可选。我一般建议常规正演用四阶PML 区域内外保持一致做高精度基准解时用八阶。高阶差分系数可以通过泰勒展开统一求解下面给出一个通用生成函数function coef diff_coef(N) % N 为半模板长度N1 对应二阶N2 对应四阶N4 对应八阶 % 求一阶导数中心差分系数在 x0 处恰等于导数 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, dx, dz, dt, nt差分更新函数计算空间导数更新 vx/vz/p差分阶数PML 初始化函数生成吸收系数在边界区域做衰减pml_len, R震源加载函数注入雷克子波或爆炸源震源位置、主频 f0画图脚本输出波场快照、接收点记录输出间隔、接收点位置拿到包之后先不要急着跑最大网格。我习惯先把 nx,nz 改成 200×200pml_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 网格步长频散的硬约束网格步长是压制频散的第一道闸门。经验准则是每个最小波长至少有 58 个网格点。最小波长由介质最小速度和震源最高有效频率决定dx ≤ c_min / (8 * f_max)其中 f_max 通常取主频的 23 倍。比如 c_min2000 m/sf030 Hz按 f_max90 Hz 算dx 应该不超过 2.78 m。资源包默认 dx5 m配合 30 Hz 主频其实偏粗所以必须用高阶差分来补偿。如果你用的是二阶差分那这个 dx 下频散会非常明显波形尾部会有长一的振荡尾巴。4.2 时间步长与 CFL 稳定性时间步长受 CFL 条件约束。二维速度-应力格式的稳定条件是c * dt * sqrt(1/dx² 1/dz²) ≤ 1若 dxdzh则 cdt/h ≤ 1/√2。实际工程中我会再留余量取 cdt/h 在 0.40.6 之间。资源包默认 dt0.001dx5c2000算出来 c*dt/h0.4是一个保守但安全的组合。如果你缩小网格步长一定要同步缩小时间步长否则高频振荡可能当场发散成 NaN。4.3 差分阶数用计算量换频散精度不同差分阶数的空间模板长度和频散表现有明显差异差分阶数空间模板点同等精度可放宽 dx计算成本频散表现二阶3基准低明显四阶5约 1.5 倍中较微弱八阶9约 2 倍高很微弱资源包里一般有一个全局参数 order设成 2 或 4。我建议如果只是看波场形态、验证网格边界四阶足够如果要做定量频散误差分析建议跑一版八阶作为参考解。切换阶数时要注意PML 区域内的差分系数也必须同步切换否则 PML 内部和外部网格阻抗不匹配会产生新的数值反射。4.4 PML层数与吸收系数PML 层数太少吸收曲线还没衰减到足够低就碰到外边界能量会反弹回来层数太多计算区域被压缩。我习惯取 1220 层。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 降低到 1520 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 扫描差分阶数比较频散差异搞清楚了基础流程这套资源最值得做的事就是用它定量对比二阶、四阶、八阶的频散误差。固定 dx5、dt0.001、f030 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 层数打印出来顺手画一条能量衰减曲线确认边界不是默认无效参数。搞明白这两个指标声波模拟的波场图才算真正可信希望帮到你。本文还有配套的精品资源点击获取