企业官网建设流程全解析

简介：这组MATLAB代码聚焦五种完美匹配层（PML）实现与效果对比，适用于FDTD等数值仿真中需要抑制边界反射的电磁波、声波或弹性波问题，也适合刚接触吸收边界条件的研究者或高年级本科生对照学习。压缩包共23个文件，包含12个.m脚本、10个.mat数据文件和1份PDF文档，体积仅301KB，代码中分别实现了Smith-PML、Steklov-Poincaré PML、多层PML、自适应PML与改进的Bermúdez PML等不同思路，并配有扫参脚本和反射计算函数，可直观比较反射系数、衰减效果与计算效率。目前已有631人学习，包内附有PML Comparison.pdf对比说明和对应的TMz仿真数据，便于快速复现五类PML在12GHz中心频率、24GHz带宽下的表现，可为实际建模中选取合适PML提供直接参考。

1. 为什么要在MATLAB里对比5种PML

做二维FDTD电磁波仿真的工程师，多半已经学会了“在区域四周加PML”来压缩反射，既然PML在理论上接近完美吸收，为什么还要对比？因为“接近完美”只在离散网格、垂直入射、有限带宽这三个条件下成立。换一种PML实现，改一行系数，可能在斜入射、低频或长时间步进时反射误差相差两个数量级，而这类误差会直接进入后续的近场、远场提取结果。这里用同一个2D TE波FDTD框架，把Berenger分裂场PML、UPML、CPML、NPML以及一个经常被误认成PML的指数衰减层同时实现，给出可运行的MATLAB代码、参数设置和量化对比方法。这组内容不是要告诉你“哪个最好”，而是给出一个能自己复现、能在论文里画对比曲线的基准。

2. 五种PML的原理选型和MATLAB关键递归

2.1 Berenger分裂场PML：作为对比基准的经典实现

1994年Berenger提出的分裂场PML，核心思路是把边界区域内的场分量按坐标方向劈裂，让波在不同方向拥有不同的虚衰减率，从而在没有真实反射界面的情况下把外向波“吞掉”。在二维TE波（Ez、Hx、Hy）中，常见做法是对Hx和Hy各保存两个子分量，Ez保持单一场，靠不同方向上的电导率σ_x、σ_y驱动衰减。

% 以x方向PML区域内的Hy为例：Hy拆成Hy_x与Hy_z Hy_x(i, j) = b_x(i) .* Hy_x(i, j) - c_x(i) .* (Ez(i+1, j) - Ez(i, j)) / dx; Hy_z(i, j) = Hy_z(i, j) + (Ez(i, j+1) - Ez(i, j)) / dy; Hy(i, j) = Hy_x(i, j) + Hy_z(i, j);

这里的b_x、c_x是预计算好的衰减系数数组，写法取自分裂场PML的常见离散格式。逻辑上，Hy_x只由Ez沿x方向的差分驱动并带上σ_x衰减，Hy_z由y方向差分驱动、不受σ_x影响，再把两者相加还原出Hy。参数要点：σ_x在同一层内用一个标量即可，但必须沿层号渐变，不能直接从0跳到σ_max；c_x通常取 (1 - b_x) 的某种比例，不必每一层独立推导系数，只要保证b_x、c_x由同一σ_x生成。这一实现是后续四种的对照基准，它的弱点是角点交界处两个方向的σ同时作用容易写错，且长时间仿真中会出现通常所说的late-time reflection，也就是低频尾反射。

2.2 CPML：卷积PML，工程上默认选项

CPML（Convolutional PML）由Roden与Gedney在2000年提出，是目前电磁仿真里默认使用的PML形式。它把坐标伸缩从频域算子改写成时域卷积，又用指数递归近似卷积结果，于是每个方向只需维护一两个辅助变量，不分裂场分量，修改起来比分裂场干净得多。

% CPML系数与辅助变量psi的递归（以x方向、Ez的x偏导为例） b_x = exp( -(sig_x ./ kappa_x + alpha_x) * dt / eps0 ); c_x = sig_x ./ (sig_x .* kappa_x + kappa_x.^2 .* alpha_x) .* (b_x - 1); psi_ez_x = b_x .* psi_ez_x + c_x .* (Ez(:, 2:nx) - Ez(:, 1:nx-1)) / dx;

逻辑说明：第一行是卷积指数的衰减因子，里面sig_x是PML电导率，kappa_x是坐标拉伸系数，alpha_x是频移因子；第二行把衰减因子转换成卷积幅度系数c_x；第三行把当前差分依次卷积进去。注意b_x、c_x是一维数组，按层预计算，psi_ez_x是二维辅助数组，维度与PML内网格一致。参数要点：起始配置建议kappa_x=1、alpha_x=0，先用最简形式跑通，再按后续章节的调参方式优化；sig_x则与分裂场共用同一多项式剖面，保证对比公平性。CPML对斜入射和掠射角的容忍度明显高于分裂场，是网格不细时最稳的选择。

2.3 UPML：用各向异性介质张量吸收边界

UPML的出发点是“各向异性介质”视角：把PML区看作一层拥有复相对介电常数张量和复磁导率张量的介质，沿边界的切向与法向分量使用不同张量分量，从而在介质层内部实现阻抗连续。MATLAB里常见的做法是预先计算一维系数数组ca_x、ca_y、cb_x、cb_y，在更新电场和磁场时根据当前格点位置选取对应系数。

% 左边界PML内，Ez更新时先取材料系数数组 Ez(i, j) = ca_x(i) .* Ez(i, j) + cb_x(i) .* ( (Hy(i, j) - Hy(i, j-1)) / dx ... - (Hx(i+1, j) - Hx(i, j)) / dy );

这段代码与常规FDTD更新式几乎一样，区别只在系数ca_x、cb_x沿PML层数渐变，而不是常数1。逻辑说明：ca_x对应时间项，cb_x对应旋度项，二者组合保证了PML介质内的波阻抗与主域一致。参数要点：实现时最容易出错的是“方向对应关系”——x方向PML内，切向是Ey、Hz，法向是Ex，而本处2D TE的坐标写得再简化也不能把ca_x、cb_x互相颠倒。UPML对垂直入射的反射性能与CPML相当，但推导形式对并行计算和色散介质耦合更友好；作为对比组，它的行为能帮判断“吸收好坏是否来自卷积近似”。

2.4 NPML：近PML，代码最简的一类

NPML由Cummer提出，思路是把PML看成对场分量做复坐标映射，但在每个格点只需要一组预计算好的插值权重，就可以把标准FDTD的导数换成复坐标下的组合导数。它和CPML解决的是同一个问题，但实现上更接近“改差分模板”而不是“维护辅助变量”。

% NPML: 用预计算的interp_fx、interp_gx替换标准x方向差分 dEz_dx = interp_fx(i) .* (Ez(i+1, j) - Ez(i, j)) / dx ... + interp_gx(i) .* (Ez(i, j) - Ez(i-1, j)) / dx;

这里interp_fx、interp_gx是初始化时根据σ_x、σ_y计算出的两层权重，二者之和应保持为1，否则模板会引入常数偏移。逻辑说明：权重随PML层位置变化，在靠近主域的一侧与标准差分几乎一致，在外边界一侧则加权混合，从而在数学上等价于坐标拉伸。参数要点：NPML不需要psi数组，内存占用比CPML少一截，但对掠射角的反射通常比CPML高几个分贝；在二维均匀网格中它是最容易改写的PML，适合先写出来做正确性验证。权重计算在初始化时完成一次，主循环里没有任何指数或除法运算，因此单步耗时最低。

2.5 指数衰减层：不是PML，但必须留一个对照组

还有一类代码在MATLAB问答社区里流传很广：把边界区场值每步乘一个指数衰减因子exp(-sigma*dt)，再继续正常更新。它在形式上“有PML三个字母”，实际只是外加阻尼，遇到波阻抗变化照样会反射。这里保留它不是为了宣传，而是让对比图里有一条“非PML”底线，用来确认其它四种确实在吸收机制上比单纯衰减高一个量级。

% DECAY层：在PML区直接衰减场（非PML） Ez(pml_x, :) = Ez(pml_x, :) .* exp(-sigma_x .* dt);

逻辑说明：这就是一阶衰减，加入后低频分量衰减慢、高频分量衰减快，频谱响应不平坦；与真正的PML相比，它缺少阻抗匹配，因此在边界处会有一次可见反射。参数要点：sigma_x可以沿用同样的剖面，但即使剖面对，也无法获得PML级别的宽频吸收。把它作为第5种实现放进统一框架，正是为了回答“如果不写PML、只加衰减，损失到底多大”这个常见问题。

2.6 五种实现的结构对比速查表

下表从实现结构角度给出一个快速判断，具体数值会随源、网格与厚度变化，后面章节再给更严谨的测量方法。

选型结论可以很直白：不在色散介质里做研究时，CPML是最省心的默认；要写教学代码或快速验证算法正确性，NPML能最快跑通；想在论文里展示不同吸收边界的差异，分裂场和指数衰减层作为对照组价值最高。

3. 用同一个FDTD框架跑通5种PML

3.1 最小可运行的MATLAB主循环骨架

五种PML放在一起时，最容易出的问题是“变量名一样但含义不同”，因此先把统一接口定下来。下面是一个最小但结构完整的MATLAB主循环，模式是主域更新、边界更新分离，五种PML各自实现同一个函数签名。

function [probe, param] = run_pml_cmp(pml_type, nx, nt) % pml_type可选 'SF' 'CPML' 'UPML' 'NPML' 'DECAY' npml = 10; dx = 1e-3; dt = dx / (phys_c * sqrt(2)); % CFL=1 [Hx, Hy, Ez] = deal(zeros(nx, nx)); probe = zeros(nt, 1); px = round(nx/2); py = round(nx/2); for n = 1:nt Ez(px, py) = Ez(px, py) + sin(2*pi*1e9*n*dt)^2; % 点源 Hx(1:end-2,:) = Hx(1:end-2,:) - dt/(mu0*dx) * diff(Ez, 1, 1); Hy(:,1:end-2) = Hy(:,1:end-2) + dt/(mu0*dx) * diff(Ez, 1, 2); Ez(2:end-1,2:end-1) = Ez(2:end-1,2:end-1) + ... dt/eps0/dx * (diff(Hy, 1, 2) - diff(Hx, 1, 1)); switch lower(pml_type) case 'cpml' [Ez, Hx, Hy] = cpml_step(Ez, Hx, Hy, pml_c, dt, npml); case 'sf' [Ez, Hx, Hy] = sf_step(Ez, Hx, Hy, pml_c, dt, npml); case 'upml' [Ez, Hx, Hy] = upml_step(Ez, Hx, Hy, pml_c, dt, npml); case 'npml' [Ez, Hx, Hy] = npml_step(Ez, Hx, Hy, pml_c, dt, npml); case 'decay' [Ez, Hx, Hy] = decay_step(Ez, Hx, Hy, pml_c, dt, npml); end probe(n) = Ez(px, py); end

代码说明：主循环沿用经典2D FDTD的差分顺序，即先更新H、再更新E，PML步骤在每次E更新后执行。pml_c是调用方预先构造的系数结构体，包含sig_x、sig_y、kappa_x、kappa_y、alpha_x、alpha_y，五种实现都从同一剖面给出，这样单步耗时和反射误差的差异才能全部归因到算法本身而不是参数不一致。注意点：源放在中心格点且只在初期激活，避免持续激励把反射波淹没在直达波里；探针记录点要在物理区内而非PML内部。启动一个仿真的命令是：

matlab -batch "run_pml_cmp('CPML', 120, 1000)"

matlab -batch在R2019a及之后都可用，适合无GUI的批处理对比。这里要求每种PML的边界函数共用同一套pml_c结构，各函数内只用索引区分“下边界、上边界、左右边界和角点”，这是整个对比脚本中最容易返工的地方。

3.2 五种PML必须共享同一份参数剖面

也许有人会在对比时为了“让每种PML表现好一点”而单独调参，这会直接破坏对比公平性。正确的做法是：先确定源类型、网格步长、PML厚度和σ剖面，让五种实现全部使用同一组预计算参数；对实现本身预留的额外自由度（如CPML的kappa与alpha、UPML的系数表）保持默认一致，除非单独标注“这是最优参数下的性能”。这样才能画出有意义的对比曲线。

3.3 从主循环里看各实现的工作量差异

从上述框架可以直观看到，五种PML在单步内的额外操作分别是：分裂场需要对Hx、Hy各自做两次带系数递归；CPML需要更新psi并对Ez加上卷积修正；UPML需要把标准更新后的场再乘一次材料系数并更新边界区；NPML只需用权重模板替换差分；DECAY层只做一次乘法。单从MATLAB向量化角度看，NPML与DECAY的边界代码最短，CPML其次，分裂场和UPML最长。我一般会根据步数规模来选择：单次仿真在万步以内，优先用CPML并接受它的psi数组；要扫几千组参数时，换NPML把单位步长开销压缩下来，能省出几小时等待。

4. 对比实验：把吸收好坏变成可写进报告的数字

4.1 用大小域差值算归一化反射误差

吸收边界性能的标准测量方法是“同源同探针双仿真”：第一次用足够大的计算域并加PML，认为探针处不会有边界反射参与，作为参考；第二次把域缩小到目标尺寸且同样加PML，探针位置不变，两次结果的差就是边界反射波。归一化表达式写成峰值误差和时域平均误差两种更合理。

% 假设ref_probe来自大域，small_probe来自目标域，二者长度一致 err = small_probe - ref_probe; norm_peak = 20 * log10(max(abs(err)) / max(abs(ref_probe))); norm_avg = 20 * log10(sqrt(sum(err.^2)) / sqrt(sum(ref_probe.^2)));

代码逻辑很清楚：峰值误差反映单次反射的最大瞬时贡献，适合观察瞬时脉冲；平均误差反映长时间仿真的整体污染水平，适合判断晚时反射。参数说明：大域的边长至少比小域多出2倍PML厚度再加20格缓冲，避免大域自身PML反射在观察时段内先到达探针；探针离边界距离要固定，否则不同PML对近场的感应差异会混入结果。量化时建议把误差信号在时域持续累积到源完全熄灭后再延长至少1000步，以便让晚时反射充分进入统计。

4.2 一个可复现的典型量级参考表

下面的量级取自均匀网格2D TE、点源、PML厚度10层、源频带与网格满足dx不超过lambda/10时的常见结果，具体项目里的绝对值会有偏移，但彼此之间的相对次序一般是稳定的。

解释下这张表怎么读：DECAY的峰值在-20dB附近，说明单纯衰减会造成显著边界回波；SF-PML如果长时间运行，峰值和平均误差之间的差距会变大，NPML的差距相对小；CPML在垂直入射且σ剖面给当时，能稳定压到-70dB以下。参数说明：如果换成40度斜入射，DECAY基本不变，NPML和分裂场可能各上升6到10dB，CPML通常仍低于-50dB，所以表格只适用于默认垂直入射。

4.3 三个可以直接落地的规律

第一，峰值反射误差主要由靠近主域的PML区域决定，因此σ剖面在前两层里必须缓慢上升，而不是把σ_max直接顶到第一个格点。第二，晚时误差主要由低频分量决定，CPML导入alpha_x之后能显著压低，分裂场即使增加厚度也很难改善，这是机制上的差异。第三，PML内部有什么物理场不重要，探针和结果提取区域必须留在物理域内，否则会把吸收率误判成反射误差。这三点在后续章节是反复出现的检查项。

5. 参数该怎么设：厚度、σ剖面、κ和α

5.1 一组可靠起点参数

PML参数不需要每次从零开始猜。常见做法是从一套经验起点出发，先跑通、再根据频谱和误差曲线微调。下面的起点适用于均匀网格、无耗背景介质。

这个经验式的直观含义：真空中eta0约377欧姆，介质中要换成介质波阻抗，也就是eta0除以sqrt(epsilon_r)（磁介质按对应修正），这样sigma_max会随介质波阻抗自动缩放。参数要点：厚度不是线性决定性能的，从10层加到16层收益可观，从16层加到22层收益会明显变缓，看曲线时可以早点停止。

5.2 σ、κ、α剖面生成代码与绘图检查

参数剖面是PML实现的入口，也是所有代码中最值得先画出来看的部分。下面是生成一维剖面并立即用MATLAB画图的代码：

% 生成PML剖面，npml为层数，m为多项式阶数 profile = (0.5 : npml-0.5) / npml; % 每层的中心位置 sigma_x = sigma_max * profile .^ m; kappa_x = 1 + (kappa_max - 1) * profile .^ m; alpha_x = alpha_max * (1 - profile) .^ m; figure; subplot(3,1,1); plot(1:npml, sigma_x); title('sigma'); subplot(3,1,2); plot(1:npml, kappa_x); title('kappa'); subplot(3,1,3); plot(1:npml, alpha_x); title('alpha');

逻辑说明：profile用每层中心位置而不是边界位置，避免在PML起始处出现半格的零厚度；sigma_x从最内层向外递增，kappa_x同向递增，alpha_x则反向递减，这样在物理域与PML交界处alpha最大、sigma最小，交界阻抗连续。参数说明：如果画出的sigma曲线在最外层还继续增大而不饱和，说明npml不够或m偏大，此时反射来自截断边界而不是PML内部；如果sigma前两格增长过快，则主域边界会有明显的一阶反射峰值。

5.3 kappa和alpha的参数踩坑

这两个参数容易调反。kappa_x增大会让坐标拉伸更强，但同时也放大该层内差分误差；kappa_max从10往上调，对垂直入射几乎无感，对45度以上掠射角才有效，因此不要一上来就设成大数。alpha_x的目的是把坐标伸缩的极点从0频率挪到非零频率，减少直流和低频分量的晚时反射，但它会削弱极低频的吸收能力——alpha_max取到经验式上限附近时，近DC分量可能反而反射回来。调参顺序推荐：先固定alpha=0，把npml和m调到垂直入射峰值达到-70dB；加入掠射角测试后需要提一档kappa，最后再看长时间误差频谱决定是否注入alpha，每一步都复用4.1的测量脚本，而不是同时调整多个参数。

5.4 网格离散度决定PML性能上限

PML性能不是孤立指标。网格越粗，差分对斜入射和近掠射角的数值相速误差越大，PML的匹配也就越不准。经验上至少保持dx不超过最短目标波长的十分之一，做掠射角对比时建议到二十分之一。时间步长按CFL条件取，2D均匀网格约dt = dx除以（c乘以sqrt(2)），不满足时PML层内部的指数系数会出现虚部误差，表现为仿真后期边界区域数值发散。判断网格是否够细的笨办法：把同一PML配置放到两倍细网格上，如果峰值反射改善量小于3dB，说明当前网格已经不是主要瓶颈。

6. 验证吸收质量与三个高频坑

6.1 用频谱而不是只盯时域判PML好坏

时域峰值误差只能反映“某一瞬间最大反射”，低频晚时反射往往在波形上只是缓慢漂移，肉眼很难察觉。更好的验证方式是做一次FFT，看反射误差的频谱是否在全频带上都低于目标值。

% 探针误差信号err做FFT，归一化到峰值并换算成dB频谱 Nf = numel(err); f = (0:Nf-1)/Nf/dt; spec = abs(fft(err)); spec_dB = 20*log10(spec / max(spec(2:end))); % 忽略直流分量 semilogx(f(2:Nf/2), spec_dB(2:Nf/2)); ylim([-100 0]);

代码逻辑：err来自4.1的大域与小域差值，它本身就是反射波序列；FFT后能看出反射能量集中在低频还是高频，从而判断是晚时反射还是离散误差；用log横轴更适合观察低频端。参数说明：采样时长要覆盖源熄灭后的多轮反射，否则FFT频率分辨率不够，低频段会出现虚假起伏。验证标准可以这样定：目标频带内反射谱低于-60dB的PML配置用于定量仿真才算合格，-40dB以上的配置只能用于波形演示。

6.2 坑一：sigma剖面第一格就顶满

许多初版代码会把sigma_max直接赋给PML每一层，导致主域与PML交界处出现阻抗阶跃，波还没进入PML就先反射一部分。判断方法是看5.2那三张曲线图：sigma应从0缓慢爬升，不是从最高值起始。修正方法是把sigma_x = sigma_max * profile.^m里的profile从0到1渐变，并保证最内层sigma相对sigma_max低于1%量级。

6.3 坑二：角点区域沿用一个方向公式

边界段与角点段必须分开处理。边段只需要一个方向的sigma和一个方向的kappa；角点段位于两个方向的PML交叠里，必须同时使用两套方向系数。常见做法是用mask矩阵分别标记上、下、左、右与四个角点，各自套不同更新式；如果只写边段不写角点，角点会成为二次辐射源，表现为延迟一段时间后从角落返回一圈明显的弧状反射。检查方法是把误差场画成空间云图，看是否有一圈圆弧从角点方向开始传播。

6.4 坑三：把探针放进PML里读数

想“直接看到吸收效果”而把探针放在PML区，读出来的曲线更像衰减包络，不是反射误差，会把PML的缺点掩盖。探针必须位于物理域内，离PML内边界至少留3到5格；物理域余量越小，近场差异越容易被误判为吸收差。最后给一个固定调参路径：用CPML先跑垂直入射，固定npml=10、m=3、alpha_max=0，把sigma_max按5.1公式初始化；确认峰值反射低于-70dB后，加45度斜入射复核，若恶化超过5dB再逐格上调kappa_max，直到频谱反射在目标频带稳定低于-60dB。

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