吸收边界条件 (Absorbing Boundary Condition) 是探地雷达 (Ground Penetrating Radar) 数值模拟及其他计算电磁学领域的研究热点.以有限的网格空间来模拟开放的空间边界，不完善的问题空间数值截断会产生严重的数值反射，甚至恶化计算结果 (Taflove and Hagness, 2000；葛德彪和闫玉波，2005；王长清和祝西里，2014)，严重影响探地雷达正演计算的精度.常见的吸收边界条件主要有：利用模零化微分算子建立的辐射边界条件 (Kriegsmann and Morawetz, 1978；Bayliss and Turkel, 1980)；基于单行波方程的各阶近似Mur吸收边界条件 (Mur, 1981a, 1981b)；将磁场分量也参与计算的超吸收边界条件 (Fang and Mei, 1988；Mei and Fang, 1992)；建立衰减层的Sarma吸收边界条件 (Sarma et al., 1998；王月英和宋建国，2007).自Berenger (1994, 1996) 提出了基于分裂场的完全匹配层 (Perfect Match Layer) 吸收边界条件以来，国内外学者对吸收边界条件开展了大量研究：Kuzuoglu和Mittra (1996)提出了复频移完全匹配层 (CFS-PML)，实现了对低频凋落波的有效吸收；Wang和Tang (2003)将非分裂的PML吸收边界条件应用于地震波的数值模拟中，提高了边界模拟精度；Roden和Gedney (2000)提出了卷积完全匹配层 (CPML) 边界条件，并开展了将其应用于任意介质的电磁波数值模拟研究；Berenger (2002)应用CFS-PML实现了对导波管中隐失波的高效吸收；刘四新和曾昭发 (2006)开展了广义完全匹配层 (GPML) 与常规PML边界条件应用于频散介质的GPR正演研究；Komatitsch和Martin (2007)将非分裂卷积的PML吸收边界条件应用于地震波波动方程中，实现了对临界入射波的有效吸收；Giannopoulos (2008)提出了递归积分完全匹配层 (RIPML) 边界条件，在CPML边界条件的基础上节约了计算内存；冯德山等 (2011)将单轴各向异性完全匹配层 (UPML) 边界条件和交替方向隐式时域有限差分法结合，对GPR进行了正演研究，有效地减少了非物理的雷达波反射；Li等 (2012)将基于复频移坐标变换的RIPML扩展到三维高阶GPR有限差分正演模拟中；冯德山等 (2012)提出了结合透射边界和Sarma边界的混合边界条件，改善了有限元法GPR正演的边界反射；薛昭等 (2014)采用间断Galerkin有限元法对弹性波进行数值模拟，推导了基于复频移拉伸坐标变换的新PML吸收边界条件 (CFS-NPML)，实现了对面波等震相的人为边界反射的有效吸收.

由此可见，吸收边界条件的算法种类诸多，且各有优势，针对不同的领域和不同的数值计算方法应选择相适应的吸收边界条件，且该研究处在不断完善和优化的过程.在这些数值计算方法中，标准交错网格的有限差分 (FDTD) 具有高效迭代、编程简单、节约内存等特点，是计算电磁学数值模拟中最常用的方法.但是该算法必须满足非常严格的CFL稳定性条件，若在单位时间步长内，电磁波传播的距离超过了单元网格的尺寸，在计算后期将会引发严重的数值色散 (Taflove and Hagness, 2000).为了解决严格的CFL稳定性条件限制，Namiki (1999, 2000) 提出了交替方向隐式FDTD算法，在不满足CFL条件下获得了稳定精确的解；Saenger和Shapiro (2002)，Saenger和Bohlen (2004)等采用旋转交错网格 (RSG) 有限差分法模拟了地震波在黏弹性介质和各向异性介质中的传播；严红勇和刘洋 (2012)提出了隐式交错的FDTD算法，以较小的计算成本精确地反映地下介质的真实情况；黄超和董良国 (2009)等提出了一种空间网格可任意奇数倍变化与时间步长任意变化的标准交错网格FDTD算法，有效地降低了空间与时间上的过采样问题，提高了模拟精度；张慧和李振春 (2011)则在时间上进行了局部变化，提出了双变网格FDTD算法，并给出了该算法的误差分析；方宏远和林皋 (2013)推导了一阶显式辛分块龙格库塔的辛算法，在获得与传统FDTD算法精度相同的情况下节约了计算时间；冯德山等 (2014)采用了粗细网格交换的亚网格技术，在计算精度和效率上取得较好平衡.

本文在Roden和Gedney，Giannopoulos，Saenger和Shapiro等的研究基础上，提出了基于非分裂卷积完全匹配层 (UCPML) 的RSG-FDTD方法.数值稳定性分析表明，该方法弱化了CFL稳定性条件的严格限制，提高了全波场数值模拟的计算效率.并以广泛存在于浅地表的软弱夹层为例，建立了衰减夹层数值模型，重点研究了持续时间较长的低频隐失波的传播特征和衰减规律，并分别从波场快照、单道波信号、时间域反射误差和频率域反射误差等方面分析了该方法的全局反射误差，UCPML有效地吸收了主区域中衰减较慢的低频隐失波，减少了有限网格空间的边界反射，更真实地模拟了开放空间.

与标准交错网格或传统FDTD相比，旋转交错网格FDTD采用不同的空间交错策略，即先计算单元网格对角线方向上的场分量及相关介电参数的差分值，然后用这些对角线方向上的差分值的线性组合来计算沿坐标轴方向上的差分值 (Saenger and Shapiro, 2002；陈浩等，2006；严红勇和刘洋，2012)，RSG差分格式在旋转坐标系 (x′, z′) 中写成如下形式：

(1)

式中z′和x′分别代表对角线方向，Δr表示单元网格中对角线的长度.根据 (1) 式可将原标准交错网格方向上的空间偏导表示成对角线方向上偏导数的线性组合，即为：

(2)

式中为RSG中的差分算子.不同的空间交错策略详见图 1、图 2所示.

图 1

Fig. 1

图 1 标准交错网格单元 Fig. 1 Cells of standard staggered grid

图 2

Fig. 2

图 2 旋转交错网格单元 Fig. 2 Cells of rotated staggered grid

以H (x, z, t)的TM_y极化模式为例，标准网格中H (x, z, t)的中心差分算子为：

(3)

联合 (2)、(3) 式可得RSG网格中的H (x, z, t)差分算子：

(4)

在RSG坐标系统中，(4) 式可以相似地写为：

(5)

在TM_y模式下，RSG网格中场分量E_y、H_x、H_z的布局如图 3所示.

图 3

Fig. 3

图 3 TMy模式的旋转交错网格场分量 Fig. 3 Field component of RSG in TMy mode

为改善对低频隐失波的吸收效果，在复系数拉伸坐标参数S_i(Kuzuoglu and Mittra, 1996) 中，除了定义衰减吸收系数σ_i外，再引入 (k_i, α_i) 两个自由可变因子，重新定义复系数拉伸坐标参数S_i：

(6)

Berenger (1994)定义k_i=1，α_i=0，而在UCPML中，仅假定σ_i和α_i为正实数，k_i≥1，定义S_i^-1的逆拉普拉斯变换为S_i，则有：

(7)

式中ζ_i(t)的离散信号定义为：

(8)

其中：

(9)

同样以H (x, z, t)的TM_y极化模式为例，Maxwell方程组中TM_y波的独立方程组为：

(10)

式中：σ^m为磁导率 (H·m^-1)，σ^e为电导率 (S·m^-1).

将其引入到上述的复系数拉伸的坐标系中，则其频率域中的表达式为：

(11)

利用逆拉普拉斯变换，方程组 (10) 中计算电场分量E_y的时间域形式为 (H_x, H_y类似)：

(12)

联合 (5)、(12) 式，重写 (12) 式，得到RSG网格中的时间域表达式：

(13)

结合 (8)、(9) 式，当取横纵向相同的均匀网格时，得到 (13) 式在时间和空间上的RSG离散方程：

(14)

为改善计算效率，引入求和辅助变量表达式ψ(i, j)，并定义：

(15)

(16)

其中：(i=x或z).

则 (14) 式的离散卷积可以递归形式进行差分离散：

(17)

在RSG算法的每一个时间步迭代前，先计算参数σ_i，α_i，k_i，a_i和b_i并储存，α_i是针对吸收低频隐失波的自由可变因子，设置为从外边界向内边界逐渐减小至零，为简化程序，引入参数符号，令系数：

(18)

对于电场分量E_y的RSG离散差分的UCPML更新方程为：

(19)

由离散方程 (17) 和 (19) 的显式差分格式可知，其数值计算必须满足相应的数值色散稳定性条件 (葛德彪和闫玉波，2005；王长清和祝西里，2014)，如图 3的场分量示意图所示，与标准交错网格相比较，在获取同等精度的情况下，RSG网格可以选择较大的时间步长，即可以选取的时间步长为标准交错网格的倍，较大的时间步长将提高主区域所有网格节点的计算效率，尤其是在精细刻画模型特征的细网格、高阶RSG或三维的数值计算中优势更为明显，其相应的数值色散稳定性条件为：

(20)

其中v为波速，dt为时间步长，Δx′和Δz′为RSG网格中的空间步长，以均等空间步长Δx′=Δz′=d为例，则在RSG方法中，时间步长的稳定性条件为：

(21)

而在标准交错网格中，对应的数值稳定性条件为：

(22)

为获取更精确的FD数值模拟结果，通常可以采用更高阶的FD算法或增加最小波长所占据的空间网格数量来实现，显然，在RSG算法中，较大的时间步长能获取与标准交错网格精度相同的数值解，且该数值解的求解过程更稳定.

为深入研究UCPML对波前传播后期的低频隐失波的吸收效果，设置了一个衰减夹层模型，如图 4所示，其中计算区域为0.8 m×0.8 m的正方形，衰减夹层为0.8 m×0.2 m的矩形，计算区域边界的UCPML及常规PML吸收边界条件的厚度均为5个网格，为了更细致地比较边界区域附近的吸收性能，将接收点和激励源放置在距离边界均为2个网格的位置，即接收点R的坐标为 (0.002，0.002)，激励源T的坐标为 (0.798，0.798).衰减介质的介电参数分别为相对介电常数ε_a=20.0，磁导率μ_a=5.0 H·m^-1，电导率σ_a=10.0 S·m^-1，衰减系数为20 dB·m^-1，激励源函数采用高斯一阶偏导的Ricker子波.为了细致地研究Ricker子波各频率成分在传播后期的绕射特征和衰减规律，设置空间步长为Δx′=Δz′=0.001 m，中心频率为f₀=500 MHz，其函数形式如下：

图 4

Fig. 4

图 4 含衰减夹层的正演模型 Fig. 4 Forward model with an attenuation interlayer embedded

(23)

式中：T₀为时间移动量，τ为决定脉冲在时域宽度的参数.

为对比两种吸收边界条件的吸收效果，计算接收点R处接收到的信号的能量振幅和反射误差，建立针对衰减夹层正演模型的参考模型，如图 5所示.与衰减夹层正演模型相比，参考模型的尺寸、介电参数、激励源与接收点的坐标、激励源信号函数均保持不变，仅将原模型的计算区域扩大4倍，且参考模型四周的边界均采用全反射的理想电导体 (PEC) 边界，以确保当激励起来的场传播至PEC边界上时能全反射回主区域.由于参考模型的问题空间为原模型的4倍，反射波达到接收点R处需要一定时间，因此，接收点R在反射波到达之前所记录的场值，与边界是开放时所记录的场值是相同的.

图 5

Fig. 5

图 5 含衰减夹层的参考模型 Fig. 5 Reference model with an attenuation interlayer embedded

针对衰减夹层正演模型，分别采用UCMPL和常规PML吸收边界条件进行了计算，并记录了不同迭代时刻E_y分量的波场快照，场值已归一化，如图 6所示.比较图 6a与图 6b，图 6c与图 6d可知，在波前面未传播至衰减夹层下界面的传播时间段，两种吸收边界条件对反射波的吸收效果表现相当，从边界上记录的场值能量幅值上看，UCPML的吸收效果略优于常规PML，而当雷达波前面经过衰减夹层后，计算区域中产生了大量的反射波和绕射波，尤其是传播后期产生的衰减较慢的低频隐失波，见图 6f中箭头所示区域，在该传播时间段内，UCPML吸收边界条件的高效吸收特性开始凸显，如图 6e三个边界上记录的场值能量所示.在图 6g与图 6h的10 ns时刻的波场快照中，对比计算区域中四个边角及四个边界平面上的绕射波能量可知，UCPML边界条件已经基本完成了主计算区域中的绕射波、反射波及低频隐失波的吸收，仅在低倾角入射的边角位置未完全吸收，如图 6g中5处箭头所示，而常规PML边界条件在边角和边界上仍存在幅值较大的非物理反射，且在主计算区域中，大量的低频隐失波仍在衰减夹层与边界之间进行多次反射和传播，见图 6h中5处箭头所示，表明在全时波场数值模拟的后期，UCPML边界条件对低频隐失波的吸收效果优于常规PML.

图 6

Fig. 6

图 6 UCPML (a, c, e, g) 与常规PML吸收边界条件 (b, d, f, h) 下Ey分量不同时刻的波场快照 (a, b) 4 ns时刻；(c, d) 6 ns时刻；(e, f) 8 ns时刻；(g, h) 10 ns时刻. Fig. 6 Snapshots of field component Ey at different times with UCPML (a, c, e, g) and regular PML (b, d, f, h)

分别采用UCPML和常规PML吸收边界条件，在图 4衰减夹层正演模型中的接收点R处 (0.002，0.002) 计算了所有迭代时间步的单道波信号，见图 7所示.在第2200步至第4300步的时间步区间，两曲线的信号重合较好，表明在波前面穿透衰减夹层后到达接收点R处的时间段，两种吸收边界条件具有相当的吸收效果，随着绕射波、层间多次波和隐失波的出现，在波前传播后期的某时间段内，两种吸收边界条件均出现了较大振幅的接收信号，见图中第4350步至第4500步时间段内记录的波场值曲线，其中矩形框内的曲线为该时间步范围段内信号放大的曲线，对比矩形框中放大的两个信号可知，UCPML吸收边界条件在吸收传播后期的隐失波方面具有明显优势.

图 7

Fig. 7

图 7 UCPML (实线) 与常规PML (虚线) 下接收点R处记录的单道波信号 Fig. 7 Single scan signal captured at receivers R with UCPML (solid line) and regular PML (dotted line)

针对实际的正演模型和扩大的参考模型，在接收点R处分别计算并储存了正演过程中不同时刻的E_y分量波场值，如图 8所示，其中虚线为参考模型中储存的E_y分量波场值，实线为实际正演模型中储存的E_y分量波场值.如图所示，两曲线在R处接收到波场的时间段重合一致，表明参考模型中数值解的正确性，且在11.4~15.7 ns之间雷达波的传播时段，参考模型曲线波场值大于实际模型曲线波场值，即在雷达波波前传播的后期，显示了吸收边界条件对接收点R处的吸收效果.同时，在16 ns之后的任意传播时间段，并未产生任何波场信号的扰动，表明参考模型的数值可以作为不同吸收边界条件的对比参考值，证明了该参考模型的有效性.

图 8

Fig. 8

图 8 实际模型 (实线) 与参考模型 (虚线) 中接收点R处的Ey波场分量 Fig. 8 Electric field component Ey at receivers R in real model (solid line) and reference model (dotted line)

为了对比不同时刻UCPML与常规PML吸收边界条件在反射边界处 (距离左边界和下边界均为2个网格) 的吸收效果差异，分别计算了接收点R处UCPML与常规PML吸收边界条件在时间域的反射误差 (TDRE)δ_t(Taflove and Hagness, 2000)：

(24)

式中δ_t为时间域反射误差，E_yⁿ为不同吸收边界条件下在R处记录的E_y分量波场值，E_y^ref为参考模型在R处记录的E_y分量波场值，为参考模型中计算的E_y波场分量中绝对值的极大值.

如图 9所示，虚线为常规PML的反射误差，实线为UCPML的反射误差，在波前传播的前期，两曲线重合一致，表明在还未产生衰减夹层绕射的时间段，两种吸收边界条件具有相当的吸收效果.随着波前传播经过衰减夹层及绕射回波的产生，UCPML的高效吸收效果优势开始体现，其反射误差明显小于常规PML，在波前传播后期及持续时间较长的隐失波出现，UCPML反射误差的收敛速度也优于常规的PML，体现了UCPML吸收边界条件对持续时间较长的隐失波高效的吸收效果.

图 9

Fig. 9

图 9 UCPML (实线) 与常规PML (虚线) 的时间域反射误差 Fig. 9 Time domain reflection errors (TDRE) with UCPML (solid line) and regular PML (dotted line)

为对比两种吸收边界条件在频率范围内的吸收效果，计算了接收点R处在频率域内的反射误差 (FDRE)δ_f(Taflove and Hagness, 2000)：

(25)

其中F为傅里叶变换.如图 10所示，对于1.5 GHz至2.0 GHz的高频率成分电磁波，两种吸收边界条件的反射误差相差并不大，且均出现了反射误差的数值振荡现象，UCPML的反射误差极大值和平均值略小于常规PML，表明UCPML对高频率成分波的吸收效果略优于常规PML；而对于20 MHz至250 MHz的低频率成分电磁波，UCPML曲线的反射误差则明显低于常规PML，表明UCPML吸收边界条件对衰减夹层模型中的低频率成分波吸收更好.

图 10

Fig. 10

图 10 UCPML (实线) 与常规PML (虚线) 的频率域反射误差 Fig. 10 Frequency domain reflection errors (FDRE) with UCPML (solid line) and regular PML (dotted line)

本文结合RSG差分算法和UCPML吸收边界条件，给出了RSG差分格式，推导了探地雷达TM_y激化模式下基于UCPML的RSG差分更新方程组，讨论了相应的数值色散稳定性条件，即在选取较大时间步长的情况下，获得了与标准交错网格相同精度的数值解，提高了探地雷达正演计算效率.同时对衰减夹层模型进行了数值模拟，重点分析了UCPML与常规PML边界条件的吸收性能、时间域全局反射误差和频率域全局反射误差.结果表明，两种吸收边界条件对早期的高频绕射波成分的吸收性能差异不大，而在衰减夹层模型的全时波场数值模拟的后期，使用常规PML边界条件时，主计算区域中的绕射效应较显著，且该绕射效应持续时间较长；全局反射误差的结果表明了UCPML的吸收效果优于常规PML，且更能有效地吸收该时间段出现的低频隐失波.