简介面向地球物理勘探、声学建模与地震模拟等方向的研究者这套基于Matlab的二维波动方程数值模拟代码提供了一套完整的PML完美匹配层边界处理方案能够有效吸收边界入射波、避免虚假反射对波场结果的污染。压缩包共7个文件全部为M脚本代码按模块拆分配置参数、波场更新、震源定义等子函数主程序负责网格初始化与时间步进迭代整体包体仅7KB结构紧凑便于阅读。已有245人学习浏览适合掌握波动方程与有限差分基础的科研人员或高年级本科生用于快速搭建二维波场模拟框架。通过学习可掌握PML层厚选取、衰减因子配置、稳定性调节等关键细节并获得一套可扩展的Matlab波场脚本便于后续更换震源类型或吸收边界方案开展对比实验。1. 二维波场模拟从固定边界换到 PML反射从 -20 dB 到 -80 dB 的差距把二维波动方程放到有限差分网格里跑第一件让人头疼的事不是震源怎么加而是网格边界会反射。早期我做地震波场模拟时用固定边界点源激发后不到几百个时间步边界反射就跟着直达波混在一起波场快照上全是同心圆叠加的假象振幅看着像“二次激发”。后来换成一阶吸收边界好一些但斜入射波和低频成分依然压不掉。真正解决问题的还是 PMLPerfectly Matched Layer完美匹配层尤其是加了复频率移位的 CPML 形式理论上可以让反射系数降到 -80 dB 以下工程上至少比固定边界干净一个量级。wavefield.zip里那套mat2d_pml代码就是给二维模拟配上 CPML 边界的一整套 Matlab 实现。main_fd.m负责时间推进defparamcpml.m和defparamcpml1d.m分别定义二维、一维边界参数defwave2d_pml.m做波场更新defsrc.m生成震源。这套文件适合做地震波传播、声学建模、探地雷达这类场景也适合第一次把边界条件换成 PML 的人拿来当骨架理论怎么落进参数、参数怎么落进差分格式、最后怎么验证边界到底有没有吸干净。2. 波动方程的时间推进与复坐标拉伸CPML 如何把边界吸收变成解析问题2.1 一阶速度-应力格式比二阶位移格式更适合接 PML二维弹性波模拟最常见的是二阶位移格式[ \rho \frac{\partial^2 u}{\partial t^2} \nabla \cdot \sigma f ]这个格式写起来简单但边界条件处理起来不够直接。PML 通常要对各方向导数做复坐标拉伸二阶位移格式要做多次空间导数拉伸因子会叠进去推导和编程都容易乱。所以我一般建议看wavefield.zip里的defwave2d_pml.m时先确认它走的是不是一阶速度-应力格式[ \rho \frac{\partial v_x}{\partial t} \frac{\partial \sigma_{xx}}{\partial x} \frac{\partial \sigma_{xz}}{\partial z} ][ \frac{\partial \sigma_{xx}}{\partial t} (\lambda 2\mu) \frac{\partial v_x}{\partial x} \lambda \frac{\partial v_z}{\partial z} ]速度-应力格式把二阶方程拆成两个一阶方程组每次只对单个方向求导PML 的拉伸因子可以独立作用在 (x)、(z) 方向边界层的各向异性天然就能表达。main_fd.m里如果按波场分量逐个更新基本就是这个结构。2.2 复坐标拉伸把反射问题转换成吸收问题PML 的核心不是直接在边界加吸收系数而是把波动方程里的空间坐标换成复坐标。以声波方程为例原本的平面波解是 (e^{i(kx - \omega t)})如果令[ \tilde{x} x - \frac{i}{\omega} \int_0^x \sigma(s) ds ]代回波动方程后解会多出一个衰减因子 (e^{-\frac{1}{\omega} \int \sigma ds})。也就是说只要在边界区域让 (\sigma) 从内到外逐渐增大进入 PML 层的波就会被解析地衰减掉而不会产生反射。(\sigma) 的渐变过程越光滑内层和边界层的阻抗匹配就越好。这里有一个关键点(\sigma) 必须平滑过渡不能从 0 直接跳到一个大数。常见的做法是沿 PML 层内距离 (d) 做幂函数渐变[ \sigma(d) \sigma_{max} \left( \frac{d}{L_{pml}} \right)^n ]其中 (n) 一般取 2 或 3。defparamcpml.m里通常会有power参数就是控制这个幂次。幂次太低过渡太陡幂次太高内部薄层衰减不足。2.3 CPML 是用辅助变量替代频率相关项原始的 PML 在离散化后对掠射波和低频分量并不稳定尤其是模拟时间拉长时边界层内会出现数值发散。CPMLConvolutional PML的做法是把复坐标拉伸改写为[ \frac{\partial}{\partial \tilde{x}} \frac{1}{\kappa_x} \frac{\partial}{\partial x} \zeta_x(t) * \left( \frac{\partial}{\partial x} \right) ]其中 (\zeta_x(t)) 是一个时间域衰减函数可以通过辅助变量 (\psi) 递推实现不需要在频域计算。递推式类似psi_x b_x * psi_x a_x * (grad_x);这里的psi_x就是 CPML 的记忆变量b_x和a_x由 (\sigma_x)、(\kappa_x)、(\alpha_x) 决定。采用 CPML 后边界层的参数从一组变成三组(\sigma) 控制衰减、(\kappa) 改善掠射波吸收、(\alpha) 抑制低频虚假反射。defparamcpml.m里如果同时出现了sigma_max、kappa_max、alpha_max那基本可以确定是 CPML 而不是经典 PML。2.4 PML 参数之间的耦合关系参数作用取值参考调参影响sigma_max最大衰减系数约 0.7 ~ 1.0 倍 ((n1) / (L_{pml} \cdot \Delta x))太小吸收不足太大反射增大kappa_max复坐标拉伸强度1 ~ 20提高对掠射波的吸收过大会导致差分不稳定alpha_max复频率平移系数约 (\pi f_0)压制低频虚假反射过大会减弱主频吸收power渐变幂次2 ~ 3决定衰减曲线形状L_pmlPML 层网格数10 ~ 20层越厚参数过渡越平缓这些参数不是独立调节的。sigma_max和power决定吸收包络kappa和alpha修正离散误差。实际中我的调参顺序是先固定power2把L_pml加厚到 15 个网格再调sigma_max到边界反射能量最低最后再动alpha_max。3. main_fd.m 主流程与 defparamcpml.m 的参数拼装3.1 压缩包里的文件各自在流程中的位置wavefield.zip解压后能看到一组文件名命名并不是随意起的。mat2d_pml是主目录里面main_fd.m是入口defparamcpml.m定义二维模型参数和 CPML 边界参数defparamcpml1d.m是一维版本用于快速试边界defsrc.m生成震源defwave.m和defwave2.m是无边界或不同边界版本的波场更新defwave2d_pml.m是带 PML 的波场更新核心。我把这套流程理解成三层。第一层是模型定义defparamcpml.m把网格数、网格间距、时间步长、PML 层数打包成一个参数结构体。第二层是波场初始化main_fd.m调用defsrc.m得到震源时间序列并把波场数组置零。第三层是迭代更新defwave2d_pml.m在每个时间步里先更新速度分量再更新应力分量同时更新 PML 层内的辅助变量。3.2 main_fd.m 的关键循环main_fd.m的常见主体循环可以简化为下面的伪代码实际文件里可能多了存储和可视化逻辑% main_fd.m 主体循环简化结构 cfl 0.25; % CFL 数受稳定性条件限制 dx param.dx; % 空间网格间距 dt cfl * dx / param.cmax; % 时间步长满足 Courant 条件 nt param.nt; % 总时间步 for it 1:nt t (it - 1) * dt; src defsrc(t, param.t0, param.f0); % 当前时刻震源值 % 用速度-应力格式更新速度分量 vx(2:end-1, :) vx(2:end-1, :) dt / rho ... * (sxx(2:end, :) - sxx(1:end-1, :)) / dx; vz(:, 2:end-1) vz(:, 2:end-1) dt / rho ... * (szz(:, 2:end) - szz(:, 1:end-1)) / dx; % 在震源位置加入点源作用力 vx(isrc, jsrc) vx(isrc, jsrc) dt * src / rho; % 更新应力分量 sxx(2:end-1, :) sxx(2:end-1, :) dt ... * (lambda 2 * mu) .* (vx(2:end-1, :) - vx(1:end-2, :)) / dx; % 调用 defwave2d_pml.m 更新 PML 层内的辅助变量 [vx, vz, sxx, szz, psi] defwave2d_pml(... vx, vz, sxx, szz, psi, param, dt); if mod(it, param.snap_interval) 0 % 保存波场快照用于后续分析 snapshot(:, :, it / param.snap_interval) sxx; end end这里代码的意图很明确先把内部网格用普通有限差分更新再在震源位置注入力最后对 PML 层做修正。PML 修正放在内部网格更新之后是为了让辅助变量能拿到当前时刻已经更新完的导数项。代码里几个关键参数param.cmax是最大波速用来计算时间步长param.t0和param.f0是震源时间函数的主频和延迟时间param.snap_interval控制每隔多少步保存一帧波场。给点源加力时src / rho的写法是把力等效成加速度源如果代码里用的是位移源单位会不同后处理时要看清输出的是应力分量sxx还是位移分量。3.3 defparamcpml.m 中每条参数为什么是这么一组以我读过的类似参数文件为例defparamcpml.m的核心输出是一个param结构体% defparamcpml.m 中的典型参数定义 param.nx 301; % 横向网格数 param.nz 301; % 纵向网格数 param.dx 2.0; % 网格间距单位 m param.dz 2.0; param.cmax 2500; % 最大波速用于计算 dt param.f0 15; % 震源主频单位 Hz param.t0 0.15; % 震源延迟时间 % CPML 边界参数 param.npml 15; % PML 层厚度网格数 param.power 2; % sigma 渐变幂次 param.sigma_max 0.85 * (param.power 1) ... / (param.npml * param.dx); % 最大衰减系数 param.kappa_max 7; % 复坐标拉伸最大值 param.alpha_max pi * param.f0;% 复频率平移最大值sigma_max的计算不是经验拍脑袋而是来自平面波在 PML 层内的反射系数公式。理论反射系数约等于 (e^{-2 \sigma_{max} L_{pml} / (n1)})要让反射低于 -60 dBsigma_max就需要满足上面的量级关系。alpha_max取pi * f0是一个通用起点它的物理意义是让低频截止区域落在主频之外避免主频附近的波被错误衰减。npml 15听起来不多但相当于给边界区加了 15 个网格的缓冲。对于 2 米网格间距PML 物理厚度是 30 米等于主波长的 0.3 倍。如果模型里包含低速层波长短PML 层需要更厚。我会根据最小波长调整% 按最小波长估算 PML 厚度 lambda_min param.cmin / param.f0; param.npml max(10, ceil(lambda_min / param.dx));这样能保证 PML 层内至少容纳一个最短波长。3.4 从 defparamcpml1d.m 到二维 pml 的维数对齐defparamcpml1d.m和defparamcpml.m看起来参数格式相似但有一个容易忽略的坑一维模型里kappa和alpha是一维数组二维模型里需要扩展成二维 mask。defwave2d_pml.m内部一般会生成两个方向上的衰减剖面% 在 x 和 z 方向构建 PML 衰减系数 for ix 1:param.nx d min(ix - 1, param.nx - ix); % 到边界的距离 d min(d, param.npml); % 超过 PML 层则封顶 sigma_x(ix) param.sigma_max * (d / param.npml)^param.power; end注意这里d的计算方式如果网格横跨1:nx左侧边界距离是ix-1右侧是nx-ix取两者较小值。这个对称渐变能保证两个方向的边界衰减率一致。二维sigma再通过外积扩展sigma_2d sigma_x * ones(1, param.nz); % 横向衰减 sigma_3d sigma_2d ones(param.nx, 1) * sigma_z; % 叠加纵向衰减叠加而不是相乘是因为 CPML 辅助变量的递推公式中对 x 和 z 方向的衰减是独立相加的。如果错误地用了相乘边界区域的等效衰减会指数增大导致边界处提前出现数值反射。4. 用 defsrc.m 和背景波场判断 PML 层是否真的“无反射”4.1 通过波场快照区分边界伪反射与残留多次波PML 调好后第一步是看波场快照。通常在模拟时间接近nt时如果 PML 正常波场进入边界区后振幅会平滑递减内部网格区域不应出现与直达波走时规律不符的新弧形波前。伪反射的特征是从边界位置产生一个与入射波极性相反的波前且它的曲率中心在边界外侧。我习惯在运行时同时输出两个量一个是内部区域的最大振幅随时间变化曲线一个是边界层内 3 个网格点处的振幅曲线。如果内部振幅曲线在波离开震源后存在异常回升哪怕幅度只有主瓣的 1%也要怀疑 PML。下面是一段用于量化反射的代码思路% 计算 PML 反射系数边界附近接收点与内部参考点能量比 function refl calc_pml_reflection(record, ref_idx, rec_idx) % record: 单个波场分量随时间演化每一列是一个网格点 e_ref sum(record(:, ref_idx).^2); % 参考点总能量 e_rec sum(record(:, rec_idx).^2); % 边界附近接收点总能量 refl_db 10 * log10(e_rec / e_ref); end把接收点放在距离 PML 内边界 5 个网格处参考点放在震源同侧内部区域。比较两个能量值时要保证计算窗口内不含震源直接激发时段否则直达波主导能量反射信号被淹没。一般从首次经过开始取时间窗比如直达波到参考点之后再取后续 2 倍时窗。4.2 检查 PML 层内的衰减曲线PML 层内的波场振幅理论上应该沿衰减方向单调下降。在defwave2d_pml.m返回波场后取一条穿过 PML 层的剖面线绘制振幅与到边界距离的关系% 沿 x 方向取一条 z 固定的剖面 profile squeeze(sxx(:, 150, 200)); % 取某时刻波场快照 x_axis (0:param.nx-1) * param.dx; % 截取进入 PML 层的部分 pml_start param.nx - param.npml; y profile(pml_start:end); x x_axis(pml_start:end);理想曲线上振幅至少下降 1 到 2 个量级。如果曲线出现先降后升通常是alpha参数过大或sigma过小如果曲线下降过早说明边界区域的有效介质属性与内部分离会产生反射。还有一个常见的定式PML 层内波场出现“拖尾振荡”多半是因为kappa_max设置过大导致更新方程的最不利特征值超过 1。这时候可以先把kappa_max降到 1确认稳定性再逐步增大。4.3 PML 层至少留多少个网格点npml不是越大越好。层内网格点增加会带来额外的内存和计算量也会让衰减曲线拉长如果模拟时间不够长波还没走出条带边界效果反而看不出来。用一张表概括我常用的配置网格间距 / 最小波长比推荐 PML 层数说明波长覆盖 5 个网格以下20 层以上低频模型低速层需加厚波长覆盖 5 ~ 10 个网格15 层大多数简单模型够用波长覆盖 10 个网格以上10 ~ 15 层可适当减少但要保证衰减渐变推荐的最小层数是 10。低于这个数sigma要在很短距离内从 0 升到很大离散误差会以虚假反射形式漏出来。5. 一维 CPML 参数文件值多少钱先跑通 defparamcpml1d.m 再回二维5.1 用 1D 测试把二维参数扫描时间缩短一个量级defparamcpml1d.m是这套压缩包里最容易被忽略的文件。二维模拟调参数每跑一次主程序要计算nx * nz个网格点快照数量又多一次全流程可能要几分钟甚至更久。而一维 CPML 只要几十个网格同一个模型几秒钟就能跑完。我通常会把二维模型的纵向切面抽出来在一维模型里先测sigma_max、alpha_max和npml得到一组可行参数后再放回二维。具体操作是先改defparamcpml1d.m里的cmin、cmax、f0让一维模型波速和二维模型一致然后跑一段模拟输出一维波场末端边界区域的最大振幅。一维里最优的alpha_max往往和二维最优值相差不大可以直接复用。% 从一维参数文件读取最优 alpha用于二维 alpha_1d_opt param_1d.alpha_max; param_2d.alpha_max alpha_1d_opt; % 直接搬过来如果一维里边界反射已经压到 -70 dB 以下二维通常不会差太多如果一维就很差二维调半天也是浪费时间。5.2 把一维衰减剖面画出来二维 PML 参数按比例放大一维测试结束后我会把一维波场在 PML 层内的衰减曲线导出来当作二维调参的参考。比如在defparamcpml1d.m里加入记录每个网格点峰值振幅的代码% 记录一维波场中每个网格点的最大绝对振幅 peak_amp max(abs(wavefield(1:end)), [], 2); peak_amp peak_amp / max(peak_amp); % 归一化绘制这条曲线时重点看衰减是否在 PML 层内均匀发生。如果衰减曲线在靠近内部区域时就开始明显下降说明sigma_max偏大波在真正的 PML 支撑带之前就被“挡”了一下等效反射源提前出现。正确的曲线是内部区域平坦PML 层内平滑下降末端接近数值精度。二维模型里由于几何扩散和角点效应边界衰减会和一维略有差异但参数基本不用大改只需要把npml按网格间距做比例缩放。一维模型比二维多一个优势可以快速找到稳定临界点。比如kappa_max从 1 一直加到 30一维模拟都能稳定但放到二维可能因为角点处双方向拉伸叠加而不稳定。所以在二维正式跑之前先在一维确认参数距离稳定性边界足够远再回到二维做最后的微调。本文还有配套的精品资源点击获取 波动方程模拟中的CPML完美匹配层：从-20dB到-80dB的边界吸收技术
简介面向地球物理勘探、声学建模与地震模拟等方向的研究者这套基于Matlab的二维波动方程数值模拟代码提供了一套完整的PML完美匹配层边界处理方案能够有效吸收边界入射波、避免虚假反射对波场结果的污染。压缩包共7个文件全部为M脚本代码按模块拆分配置参数、波场更新、震源定义等子函数主程序负责网格初始化与时间步进迭代整体包体仅7KB结构紧凑便于阅读。已有245人学习浏览适合掌握波动方程与有限差分基础的科研人员或高年级本科生用于快速搭建二维波场模拟框架。通过学习可掌握PML层厚选取、衰减因子配置、稳定性调节等关键细节并获得一套可扩展的Matlab波场脚本便于后续更换震源类型或吸收边界方案开展对比实验。1. 二维波场模拟从固定边界换到 PML反射从 -20 dB 到 -80 dB 的差距把二维波动方程放到有限差分网格里跑第一件让人头疼的事不是震源怎么加而是网格边界会反射。早期我做地震波场模拟时用固定边界点源激发后不到几百个时间步边界反射就跟着直达波混在一起波场快照上全是同心圆叠加的假象振幅看着像“二次激发”。后来换成一阶吸收边界好一些但斜入射波和低频成分依然压不掉。真正解决问题的还是 PMLPerfectly Matched Layer完美匹配层尤其是加了复频率移位的 CPML 形式理论上可以让反射系数降到 -80 dB 以下工程上至少比固定边界干净一个量级。wavefield.zip里那套mat2d_pml代码就是给二维模拟配上 CPML 边界的一整套 Matlab 实现。main_fd.m负责时间推进defparamcpml.m和defparamcpml1d.m分别定义二维、一维边界参数defwave2d_pml.m做波场更新defsrc.m生成震源。这套文件适合做地震波传播、声学建模、探地雷达这类场景也适合第一次把边界条件换成 PML 的人拿来当骨架理论怎么落进参数、参数怎么落进差分格式、最后怎么验证边界到底有没有吸干净。2. 波动方程的时间推进与复坐标拉伸CPML 如何把边界吸收变成解析问题2.1 一阶速度-应力格式比二阶位移格式更适合接 PML二维弹性波模拟最常见的是二阶位移格式[ \rho \frac{\partial^2 u}{\partial t^2} \nabla \cdot \sigma f ]这个格式写起来简单但边界条件处理起来不够直接。PML 通常要对各方向导数做复坐标拉伸二阶位移格式要做多次空间导数拉伸因子会叠进去推导和编程都容易乱。所以我一般建议看wavefield.zip里的defwave2d_pml.m时先确认它走的是不是一阶速度-应力格式[ \rho \frac{\partial v_x}{\partial t} \frac{\partial \sigma_{xx}}{\partial x} \frac{\partial \sigma_{xz}}{\partial z} ][ \frac{\partial \sigma_{xx}}{\partial t} (\lambda 2\mu) \frac{\partial v_x}{\partial x} \lambda \frac{\partial v_z}{\partial z} ]速度-应力格式把二阶方程拆成两个一阶方程组每次只对单个方向求导PML 的拉伸因子可以独立作用在 (x)、(z) 方向边界层的各向异性天然就能表达。main_fd.m里如果按波场分量逐个更新基本就是这个结构。2.2 复坐标拉伸把反射问题转换成吸收问题PML 的核心不是直接在边界加吸收系数而是把波动方程里的空间坐标换成复坐标。以声波方程为例原本的平面波解是 (e^{i(kx - \omega t)})如果令[ \tilde{x} x - \frac{i}{\omega} \int_0^x \sigma(s) ds ]代回波动方程后解会多出一个衰减因子 (e^{-\frac{1}{\omega} \int \sigma ds})。也就是说只要在边界区域让 (\sigma) 从内到外逐渐增大进入 PML 层的波就会被解析地衰减掉而不会产生反射。(\sigma) 的渐变过程越光滑内层和边界层的阻抗匹配就越好。这里有一个关键点(\sigma) 必须平滑过渡不能从 0 直接跳到一个大数。常见的做法是沿 PML 层内距离 (d) 做幂函数渐变[ \sigma(d) \sigma_{max} \left( \frac{d}{L_{pml}} \right)^n ]其中 (n) 一般取 2 或 3。defparamcpml.m里通常会有power参数就是控制这个幂次。幂次太低过渡太陡幂次太高内部薄层衰减不足。2.3 CPML 是用辅助变量替代频率相关项原始的 PML 在离散化后对掠射波和低频分量并不稳定尤其是模拟时间拉长时边界层内会出现数值发散。CPMLConvolutional PML的做法是把复坐标拉伸改写为[ \frac{\partial}{\partial \tilde{x}} \frac{1}{\kappa_x} \frac{\partial}{\partial x} \zeta_x(t) * \left( \frac{\partial}{\partial x} \right) ]其中 (\zeta_x(t)) 是一个时间域衰减函数可以通过辅助变量 (\psi) 递推实现不需要在频域计算。递推式类似psi_x b_x * psi_x a_x * (grad_x);这里的psi_x就是 CPML 的记忆变量b_x和a_x由 (\sigma_x)、(\kappa_x)、(\alpha_x) 决定。采用 CPML 后边界层的参数从一组变成三组(\sigma) 控制衰减、(\kappa) 改善掠射波吸收、(\alpha) 抑制低频虚假反射。defparamcpml.m里如果同时出现了sigma_max、kappa_max、alpha_max那基本可以确定是 CPML 而不是经典 PML。2.4 PML 参数之间的耦合关系参数作用取值参考调参影响sigma_max最大衰减系数约 0.7 ~ 1.0 倍 ((n1) / (L_{pml} \cdot \Delta x))太小吸收不足太大反射增大kappa_max复坐标拉伸强度1 ~ 20提高对掠射波的吸收过大会导致差分不稳定alpha_max复频率平移系数约 (\pi f_0)压制低频虚假反射过大会减弱主频吸收power渐变幂次2 ~ 3决定衰减曲线形状L_pmlPML 层网格数10 ~ 20层越厚参数过渡越平缓这些参数不是独立调节的。sigma_max和power决定吸收包络kappa和alpha修正离散误差。实际中我的调参顺序是先固定power2把L_pml加厚到 15 个网格再调sigma_max到边界反射能量最低最后再动alpha_max。3. main_fd.m 主流程与 defparamcpml.m 的参数拼装3.1 压缩包里的文件各自在流程中的位置wavefield.zip解压后能看到一组文件名命名并不是随意起的。mat2d_pml是主目录里面main_fd.m是入口defparamcpml.m定义二维模型参数和 CPML 边界参数defparamcpml1d.m是一维版本用于快速试边界defsrc.m生成震源defwave.m和defwave2.m是无边界或不同边界版本的波场更新defwave2d_pml.m是带 PML 的波场更新核心。我把这套流程理解成三层。第一层是模型定义defparamcpml.m把网格数、网格间距、时间步长、PML 层数打包成一个参数结构体。第二层是波场初始化main_fd.m调用defsrc.m得到震源时间序列并把波场数组置零。第三层是迭代更新defwave2d_pml.m在每个时间步里先更新速度分量再更新应力分量同时更新 PML 层内的辅助变量。3.2 main_fd.m 的关键循环main_fd.m的常见主体循环可以简化为下面的伪代码实际文件里可能多了存储和可视化逻辑% main_fd.m 主体循环简化结构 cfl 0.25; % CFL 数受稳定性条件限制 dx param.dx; % 空间网格间距 dt cfl * dx / param.cmax; % 时间步长满足 Courant 条件 nt param.nt; % 总时间步 for it 1:nt t (it - 1) * dt; src defsrc(t, param.t0, param.f0); % 当前时刻震源值 % 用速度-应力格式更新速度分量 vx(2:end-1, :) vx(2:end-1, :) dt / rho ... * (sxx(2:end, :) - sxx(1:end-1, :)) / dx; vz(:, 2:end-1) vz(:, 2:end-1) dt / rho ... * (szz(:, 2:end) - szz(:, 1:end-1)) / dx; % 在震源位置加入点源作用力 vx(isrc, jsrc) vx(isrc, jsrc) dt * src / rho; % 更新应力分量 sxx(2:end-1, :) sxx(2:end-1, :) dt ... * (lambda 2 * mu) .* (vx(2:end-1, :) - vx(1:end-2, :)) / dx; % 调用 defwave2d_pml.m 更新 PML 层内的辅助变量 [vx, vz, sxx, szz, psi] defwave2d_pml(... vx, vz, sxx, szz, psi, param, dt); if mod(it, param.snap_interval) 0 % 保存波场快照用于后续分析 snapshot(:, :, it / param.snap_interval) sxx; end end这里代码的意图很明确先把内部网格用普通有限差分更新再在震源位置注入力最后对 PML 层做修正。PML 修正放在内部网格更新之后是为了让辅助变量能拿到当前时刻已经更新完的导数项。代码里几个关键参数param.cmax是最大波速用来计算时间步长param.t0和param.f0是震源时间函数的主频和延迟时间param.snap_interval控制每隔多少步保存一帧波场。给点源加力时src / rho的写法是把力等效成加速度源如果代码里用的是位移源单位会不同后处理时要看清输出的是应力分量sxx还是位移分量。3.3 defparamcpml.m 中每条参数为什么是这么一组以我读过的类似参数文件为例defparamcpml.m的核心输出是一个param结构体% defparamcpml.m 中的典型参数定义 param.nx 301; % 横向网格数 param.nz 301; % 纵向网格数 param.dx 2.0; % 网格间距单位 m param.dz 2.0; param.cmax 2500; % 最大波速用于计算 dt param.f0 15; % 震源主频单位 Hz param.t0 0.15; % 震源延迟时间 % CPML 边界参数 param.npml 15; % PML 层厚度网格数 param.power 2; % sigma 渐变幂次 param.sigma_max 0.85 * (param.power 1) ... / (param.npml * param.dx); % 最大衰减系数 param.kappa_max 7; % 复坐标拉伸最大值 param.alpha_max pi * param.f0;% 复频率平移最大值sigma_max的计算不是经验拍脑袋而是来自平面波在 PML 层内的反射系数公式。理论反射系数约等于 (e^{-2 \sigma_{max} L_{pml} / (n1)})要让反射低于 -60 dBsigma_max就需要满足上面的量级关系。alpha_max取pi * f0是一个通用起点它的物理意义是让低频截止区域落在主频之外避免主频附近的波被错误衰减。npml 15听起来不多但相当于给边界区加了 15 个网格的缓冲。对于 2 米网格间距PML 物理厚度是 30 米等于主波长的 0.3 倍。如果模型里包含低速层波长短PML 层需要更厚。我会根据最小波长调整% 按最小波长估算 PML 厚度 lambda_min param.cmin / param.f0; param.npml max(10, ceil(lambda_min / param.dx));这样能保证 PML 层内至少容纳一个最短波长。3.4 从 defparamcpml1d.m 到二维 pml 的维数对齐defparamcpml1d.m和defparamcpml.m看起来参数格式相似但有一个容易忽略的坑一维模型里kappa和alpha是一维数组二维模型里需要扩展成二维 mask。defwave2d_pml.m内部一般会生成两个方向上的衰减剖面% 在 x 和 z 方向构建 PML 衰减系数 for ix 1:param.nx d min(ix - 1, param.nx - ix); % 到边界的距离 d min(d, param.npml); % 超过 PML 层则封顶 sigma_x(ix) param.sigma_max * (d / param.npml)^param.power; end注意这里d的计算方式如果网格横跨1:nx左侧边界距离是ix-1右侧是nx-ix取两者较小值。这个对称渐变能保证两个方向的边界衰减率一致。二维sigma再通过外积扩展sigma_2d sigma_x * ones(1, param.nz); % 横向衰减 sigma_3d sigma_2d ones(param.nx, 1) * sigma_z; % 叠加纵向衰减叠加而不是相乘是因为 CPML 辅助变量的递推公式中对 x 和 z 方向的衰减是独立相加的。如果错误地用了相乘边界区域的等效衰减会指数增大导致边界处提前出现数值反射。4. 用 defsrc.m 和背景波场判断 PML 层是否真的“无反射”4.1 通过波场快照区分边界伪反射与残留多次波PML 调好后第一步是看波场快照。通常在模拟时间接近nt时如果 PML 正常波场进入边界区后振幅会平滑递减内部网格区域不应出现与直达波走时规律不符的新弧形波前。伪反射的特征是从边界位置产生一个与入射波极性相反的波前且它的曲率中心在边界外侧。我习惯在运行时同时输出两个量一个是内部区域的最大振幅随时间变化曲线一个是边界层内 3 个网格点处的振幅曲线。如果内部振幅曲线在波离开震源后存在异常回升哪怕幅度只有主瓣的 1%也要怀疑 PML。下面是一段用于量化反射的代码思路% 计算 PML 反射系数边界附近接收点与内部参考点能量比 function refl calc_pml_reflection(record, ref_idx, rec_idx) % record: 单个波场分量随时间演化每一列是一个网格点 e_ref sum(record(:, ref_idx).^2); % 参考点总能量 e_rec sum(record(:, rec_idx).^2); % 边界附近接收点总能量 refl_db 10 * log10(e_rec / e_ref); end把接收点放在距离 PML 内边界 5 个网格处参考点放在震源同侧内部区域。比较两个能量值时要保证计算窗口内不含震源直接激发时段否则直达波主导能量反射信号被淹没。一般从首次经过开始取时间窗比如直达波到参考点之后再取后续 2 倍时窗。4.2 检查 PML 层内的衰减曲线PML 层内的波场振幅理论上应该沿衰减方向单调下降。在defwave2d_pml.m返回波场后取一条穿过 PML 层的剖面线绘制振幅与到边界距离的关系% 沿 x 方向取一条 z 固定的剖面 profile squeeze(sxx(:, 150, 200)); % 取某时刻波场快照 x_axis (0:param.nx-1) * param.dx; % 截取进入 PML 层的部分 pml_start param.nx - param.npml; y profile(pml_start:end); x x_axis(pml_start:end);理想曲线上振幅至少下降 1 到 2 个量级。如果曲线出现先降后升通常是alpha参数过大或sigma过小如果曲线下降过早说明边界区域的有效介质属性与内部分离会产生反射。还有一个常见的定式PML 层内波场出现“拖尾振荡”多半是因为kappa_max设置过大导致更新方程的最不利特征值超过 1。这时候可以先把kappa_max降到 1确认稳定性再逐步增大。4.3 PML 层至少留多少个网格点npml不是越大越好。层内网格点增加会带来额外的内存和计算量也会让衰减曲线拉长如果模拟时间不够长波还没走出条带边界效果反而看不出来。用一张表概括我常用的配置网格间距 / 最小波长比推荐 PML 层数说明波长覆盖 5 个网格以下20 层以上低频模型低速层需加厚波长覆盖 5 ~ 10 个网格15 层大多数简单模型够用波长覆盖 10 个网格以上10 ~ 15 层可适当减少但要保证衰减渐变推荐的最小层数是 10。低于这个数sigma要在很短距离内从 0 升到很大离散误差会以虚假反射形式漏出来。5. 一维 CPML 参数文件值多少钱先跑通 defparamcpml1d.m 再回二维5.1 用 1D 测试把二维参数扫描时间缩短一个量级defparamcpml1d.m是这套压缩包里最容易被忽略的文件。二维模拟调参数每跑一次主程序要计算nx * nz个网格点快照数量又多一次全流程可能要几分钟甚至更久。而一维 CPML 只要几十个网格同一个模型几秒钟就能跑完。我通常会把二维模型的纵向切面抽出来在一维模型里先测sigma_max、alpha_max和npml得到一组可行参数后再放回二维。具体操作是先改defparamcpml1d.m里的cmin、cmax、f0让一维模型波速和二维模型一致然后跑一段模拟输出一维波场末端边界区域的最大振幅。一维里最优的alpha_max往往和二维最优值相差不大可以直接复用。% 从一维参数文件读取最优 alpha用于二维 alpha_1d_opt param_1d.alpha_max; param_2d.alpha_max alpha_1d_opt; % 直接搬过来如果一维里边界反射已经压到 -70 dB 以下二维通常不会差太多如果一维就很差二维调半天也是浪费时间。5.2 把一维衰减剖面画出来二维 PML 参数按比例放大一维测试结束后我会把一维波场在 PML 层内的衰减曲线导出来当作二维调参的参考。比如在defparamcpml1d.m里加入记录每个网格点峰值振幅的代码% 记录一维波场中每个网格点的最大绝对振幅 peak_amp max(abs(wavefield(1:end)), [], 2); peak_amp peak_amp / max(peak_amp); % 归一化绘制这条曲线时重点看衰减是否在 PML 层内均匀发生。如果衰减曲线在靠近内部区域时就开始明显下降说明sigma_max偏大波在真正的 PML 支撑带之前就被“挡”了一下等效反射源提前出现。正确的曲线是内部区域平坦PML 层内平滑下降末端接近数值精度。二维模型里由于几何扩散和角点效应边界衰减会和一维略有差异但参数基本不用大改只需要把npml按网格间距做比例缩放。一维模型比二维多一个优势可以快速找到稳定临界点。比如kappa_max从 1 一直加到 30一维模拟都能稳定但放到二维可能因为角点处双方向拉伸叠加而不稳定。所以在二维正式跑之前先在一维确认参数距离稳定性边界足够远再回到二维做最后的微调。本文还有配套的精品资源点击获取