企业官网建设流程全解析

简介：面向地球物理勘探、声学建模与地震模拟等方向的研究者，这套基于Matlab的二维波动方程数值模拟代码提供了一套完整的PML（完美匹配层）边界处理方案，能够有效吸收边界入射波、避免虚假反射对波场结果的污染。压缩包共7个文件，全部为M脚本，代码按模块拆分配置参数、波场更新、震源定义等子函数，主程序负责网格初始化与时间步进迭代，整体包体仅7KB，结构紧凑便于阅读。已有245人学习浏览，适合掌握波动方程与有限差分基础的科研人员或高年级本科生，用于快速搭建二维波场模拟框架。通过学习可掌握PML层厚选取、衰减因子配置、稳定性调节等关键细节，并获得一套可扩展的Matlab波场脚本，便于后续更换震源类型或吸收边界方案开展对比实验。

1. 二维波场模拟从固定边界换到 PML：反射从 -20 dB 到 -80 dB 的差距

把二维波动方程放到有限差分网格里跑，第一件让人头疼的事不是震源怎么加，而是网格边界会反射。早期我做地震波场模拟时用固定边界，点源激发后不到几百个时间步，边界反射就跟着直达波混在一起，波场快照上全是同心圆叠加的假象，振幅看着像“二次激发”。后来换成一阶吸收边界，好一些，但斜入射波和低频成分依然压不掉。真正解决问题的还是 PML（Perfectly 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 在离散化后对掠射波和低频分量并不稳定，尤其是模拟时间拉长时，边界层内会出现数值发散。CPML（Convolutional 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和power决定吸收包络，kappa和alpha修正离散误差。实际中我的调参顺序是：先固定power=2，把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} / (n+1)})，要让反射低于 -60 dB，sigma_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不是越大越好。层内网格点增加会带来额外的内存和计算量，也会让衰减曲线拉长，如果模拟时间不够长，波还没走出条带，边界效果反而看不出来。用一张表概括我常用的配置：

推荐的最小层数是 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，一维模拟都能稳定，但放到二维可能因为角点处双方向拉伸叠加而不稳定。所以在二维正式跑之前，先在一维确认参数距离稳定性边界足够远，再回到二维做最后的微调。

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