多源并发下Mur二阶吸收边界和非分裂递归卷积完全匹配层对比研究

崔凡，陈毅，薛晗鹏，彭苏萍，杜云飞

(1.中国矿业大学(北京) 地球科学与测绘工程学院，北京 100083;2.中国矿业大学(北京) 煤炭资源与安全开采国家重点实验室，北京 100083)

0 引言

目前，在时域有限差分方法中运用边界条件来截断网格存在很多方法。Bayliss和Turkel[1]采用外行波的模拟法作为吸收边界条件，实现了对电磁波的简单吸收；Engquist和Majda[2]提出了基于单向波动方程的Engquist-Majda吸收边界条件；Mur[3]给出了波动方程的各阶近似及差分形式，表明了在高阶近似时，对电磁波的吸收效果较好，然而对于各向异性介质和色散介质，以二阶近似Mur作为吸收边界条件会在网格截断处发生强反射。之后，Berenger[4]提出了完全匹配层吸收边界条件，他将电磁波在吸收边界区域进行波场分裂，并对各个分裂场赋予不同的损耗值，这相当于在FDTD网格外加入一层吸收媒介，它具有不依赖于外向传播的电磁波入射角及频率的波阻抗，因此，进入完全匹配层内的电磁波快速地衰减，实现了对不同频率、不同角度的入射波较好的吸收效果。有很多研究发现[5-7]，Berenger完全匹配层作为吸收边界条件比旁轴近似吸收边界条件、Higdon吸收边界条件[8]、廖氏吸收边界条件[9]和指数衰减吸收边界条件[10]具有更好地吸收效果。但是，Berenger的PML理论违背了波的折射原理，在Maxwell方程中是不能实现的，同时，电磁场的分裂增加了数值计算的难度，计算效率低，而且只对行波具有吸收效果，对空域中衰减的隐失波、低频波以及入射角度较小的掠射波吸收效果差[11]。为了改善边界条件的吸收效果，Sacks等[12]，Gedney等[13]提出了单轴各向异性完全匹配层作为吸收边界条件，该方法在数学上与Berenger提出的经典完全匹配层等价，它是基于Maxwell方程的，不分裂波场，而是在吸收系数中引入线性吸收因子，实现了对隐失波、低频波以及在PML层内凋落波等干扰波的吸收，该方法作为良好的吸收边界条件在差分算法中被广泛使用。如肖明顺等[14]、詹应林等[15]将各向异性完全匹配层应用于二维探地雷达的正演中；中南大学冯德山等[16-18]实现了在二维和三维空间中将各向异性完全匹配层引入交替方向隐式时域有限差分算法中(ADI-FDTD)；吉林大学的李静等[19]人开展了三维各向异性完全匹配层作为吸收边界条件的高阶时域有限差分正演。同一时间段内，Chew等[20]提出了拉伸坐标完全匹配层理论；Kuzuoglu等[21]在复频域内的完全匹配层变量Sk中引进低频分量吸收参数，将复平面的极点从实轴移动到虚轴，加强了对低频波和隐失波的吸收。然而到现研究阶段，所有的吸收边界条件都是在单个激励源的条件下对电磁波进行数值模拟。

本文基于Maxwell方程推导了一种在时域有限差分算法中运用递推卷积求解电磁场的方法，该方法在离散网格条件下通过递推形式计算卷积，不分裂变量，直接计算卷积，避免了对卷积求解的复杂计算。并将该方法运用到多源并发下对电磁波的吸收，为阵列探地雷达数值仿真做了铺垫。

1 基本原理

1.1 坐标伸缩变化下Maxwell卷积方程

非递归卷积完全匹配层主要求解是在麦克斯韦方程的频率域进行，对连续和时谐场，麦克斯韦旋度支配方程可写为：

(2)

在PML中电导率不是由物理空间模型的基本参数决定的，而是设置电导率参数使网格截断的反射最小。在此意义上，电导率σ具有任意性，将一个从属于特定位置的相对介电常数来规范电导率，在分裂场PML条件下定义：

(3)

当频率趋于0时，Sw没有意义。这将导致低频引发数值异常，为修正这个问题，增加额外的参数来保证频率趋于0时Sw值的有限性。引入更一般化的坐标伸缩因子Si：

(4)

式中：Ki是改善PML对表面波吸收的参数；σi为PML层内k方向电导率参数；αi为改善PML对低频分量的吸收参数；σi和αi都为大于零的正实数；Ki≥1。在三维空间下，i∈{x,y,z}，每一个Si项总是与i∈{x,y,z}方向的道数是成对的。有耗介质中，在TM极化模式下(Hx,Hy,Ez)，三维空间下电场频域方程表达式为：

(5)

将方程(5)通过傅里叶逆变换转化到时域，由于坐标伸缩因子与频率无关，所以在时域内方程右边存在卷积，即：

(7)

对坐标伸缩因子的倒数求傅里叶逆变换有：

其中δ(t)是Dirac冲击函数，u(t)是单位阶跃函数。定义：

i=(x,y,z) 。

(9)

因为Dirac冲击函数和其他函数的卷积仍为原函数，即δ(t)*f(t)=f(t)，将式(10)代入式(6)有：

(11)

式(11)为电场z分量的时域递推卷积方程，直接在递推方程(14)中计算时域卷积效率是很低的，为有效计算卷积，下面给出递归卷积求解原理和在时域有限差分方法中的实现步骤。

1.2 基于递归卷积的非分裂场PML的FDTD实现

(12)

式中：Eui表明这个函数将会在电场Eu分量上更新，并与i(i=x,y,z)方向的空间导数有关。假设积分变量τ在式(12)中是连续变化的，由于时域有限差分方法的离散性，所以Hv场是离散变化的。用连续变量t表示Hv时，只取离散值，即

(13)

卷积包含f(qΔt-τ)函数。在时间步为q=0时，函数值为f(-τ)，如图1b所示。在所有采样点fn都关于原点反转对称，在-Δt≤τ<0上函数值为f1，在-2Δt≤τ<-Δt上函数值为f2，在-3Δt≤τ<-2Δt上函数值为f3，其余各时刻的函数值以此类推。

图1 函数f(t)阶梯表示法Fig.1 Stepped representation of function f(t)

图1c展示了在特定时刻q=4时f(qΔt-τ)的函数值。如图所示，扩展到τ=0右侧右边的第一个脉冲具有fq的值，扩展到τ=Δt右边的值为fq-1，扩展到τ=2Δt右边的值为fq-2,以此类推。此时移动函数可表示为：

(14)

对磁场求导数有：

(15)

(16)

式(16)中，脉冲函数pk(t)用来确定积分上下限。考虑如下特殊积分：

(17)

把式(16)中的ζi(t)离散冲击响应定义为：

(18)

式中：

i=(x,y,z) ，

(19)

将式(17)、(18)、(19)、(20)代入式(16)中有：

(21)

将k=0项和其他项分离有：

(22)

令n=k-1代替指数项，即k=n+1：

(23)

用k来表示指数项n有：

(24)

对比式(24)和式(21)有：

(25)

空间离散下的电场递推式，按照Yee氏网格将式(25)进行时间和空间上的离散，可得：

(26)

上式中Z0i(k)项含有卷积的计算，离散卷积计算十分复杂，但是同时Z0i(k)项是简单的指数形式，从而它们的和可以通过递归卷积来得到，引入一组新的辅助表达式φi：

(i=x,y)

(27)

将φi代入式(26)经整理后得到电场Ez分量的递推表达式：

(28)

式(28)中：σ为电导率，ε为介质的介电常数，m=(i,j,k)，在仿真时令模拟区域的K值为1，在PML层内使K值渐进变化，PML层为有限的厚度单元，σ、K和α在PML层中单调变化：

(29)

(30)

(31)

1.3 Mur二阶吸收边界条件的FDTD实现

在二维直角坐标系下，行波的平面波解为：

f(x,y,t)=Aexp[j(ωt-αxx-αyy)] ，

(32)

假设在x=0的位置设置截断边界，则在x≥0的区域会同时存在入射波和反射波。因此式(32)可写为：

(33)

可将上式分解为左行波f-(入射波)和右行波f+(反射波)之和：

(34)

令微分算子L使得L±f±=0，其中：

(36)

对左行波使得L-f-|x=0中并令f=Ez，化简后得到：

(37)

(38)

在(i+1/2,j)点和t=(n+1/2)Δt时刻对式(62)进行差分离散并进行线性插值可得到左截断边界处Mur的二阶吸收边界条件递推式：

(39)

式(39)中，v为电磁波传播速度，μ为磁导率，Δx、Δy、Δt分别为x、y方向的空间步长和时间步长。式(39)给出的Mur二阶吸收边界条件递推公式没有在边角点处对电场值进行修正，因此，在4个边角位置处发生明显的波反射，如图2所示，给出了激励源位于模型区域中心位置处，Ez波场在600个时间步的变化，在第400个时间步时，由于没有对边角进行修正，此时4个角点位置处都发生了明显的反射，在500和600时间步时刻的研究区域内存在4个角点反射的波场。

图2 未修正角点不同时间步的Ez波场Fig.2 Ez wave field at different time steps of uncorrected corners

1.4 Mur二阶吸收边界条件的角点修正

(40)

创建新坐标系ηoξ，新坐标系于原始坐标系逆时针旋转45°，在新坐标系下Mur一阶近似公式可表示为：

(41)

图3 TM波左下角点修正示意Fig.3 Schematic diagram of TM wave lower left corner correction

由于P点与原点在同一离散网格内，忽略波振幅衰减，利用线性插值可得到：

。(42)

根据式(40)结合式(42),当Δx≠Δy时，有：

(43)

根据式(43)，修正角点位置处的电场值取决于该角点相邻一个时间步位置处前一时刻的电场值和该点的速度值。

2 模型验证

2.1 均匀介质

为验证在均匀介质下Mur二阶近似吸收边界条件和基于递归卷积完全匹配层作为吸收边界条件在多源并发下对电磁波的吸收效果，设置模拟区域网格大小为200×200，四周PML层厚度为10个网格。模拟区域的物性参数为：εr=3.5,ur=1,σ=2.5 mS/m。两个激励源采用主频为900 MHz的布莱克曼—哈里斯脉冲，两个激励源的位置沿x轴相距80个空间步长，空间步长为0.006 m，时间步长为 0.015 ns。不同时间步的Ez波场如图4所示。图4为应用递归卷积完全匹配层作为吸收边界条件的Ez波场快照，仿真时间持续400个时间步长。在第100时间步时(图4a)，两个并发激励源电磁波相遇，经相互干涉作用后继续扩散；在时间步200时(图4b)，部分电磁波进入PML层被吸收掉； 时间步300时(图4c)两个激励源产生的扩散电磁波到达彼此对面PML层，此时部分波到达边界进入PML层无任何反射发生；时间步400时(图4d)，波扩散离开研究区域，在整个仿真持续时间段内，无任何明显反射发生，电磁场在网格截断位置处被PML层有效吸收掉了。

图5给出了Mur二阶吸收边界条件不同时间步Ez波场快照，二阶Mur吸收边界条件是运用微分二阶近似来求解单向行波方程，它是将微分方程进行泰勒级数展开后去掉高阶项的近似算法，此时反射系数不为零，因此不能实现对波的良好吸收。在两个激励源并发的条件下，电磁波相互干涉，重构后的波场形成小方形包络，如图5c箭头所示位置，随着时间推移，小方形包络在水平方向上扩展开来，以平面波束形式达到网格边界位置。这时，Mur二阶吸收边界对电磁波吸收效果较差，会在边界处发生反射。在前300个时间步内(图5a,5b,5c)，由单个激励源产生的电磁波在均匀介质中扩散，以球面的形式到达网格边界处， 在第400个时间步时(图5d)，经过干涉后的波到达边界，此时部分波返回研究区域；在时间步500时，波已经扩散离开研究区，然而因吸收不完全导致部分反射波返回研究区域，如图5e、5f所示，此时波在边界位置发生多次反射，形成谐振，对下一时间段内的有效波场造成二次干扰。

图4 均匀介质递归卷积完全匹配层不同时间步Ez波场快照Fig.4 Snapshots of Ez wave field at different time steps of the homogeneous medium recursive convolution perfectly matched layer

图5 均匀介质Mur二阶吸收边界不同时间步Ez波场快照Fig.5 Snapshots of Ez wave field at different time steps of Mur second-order absorbing boundary in homogeneous medium

2.2 低洼模型

在复杂地质构造条件下，在地下介质空间中扩散的电磁波会产生反射、绕射现象。绕射波会以不同的角度入射到匹配层中，为体现两种不同的吸收边界对波的吸收效果，构造了复杂低洼地质模型，模拟现实地质环境中的正、逆断层。设置模型区域网格大小为300×400，模拟区域的物性参数为：上层介质相对介电常数ε1=8，电导率σ1=0.003 S/m，下层介质相对介电常数ε2=5，电导率σ2=0.001 S/m。所有激励源采用主频为900 MHz的布莱克曼—哈里斯脉冲，空间步长为0.006 m，时间步长为 0.015 ns，26个激励源的位置沿x轴分布，每个激励源相距5个空间步长，其模型的空间分布如图6所示。对于递归卷积完全匹配层在四周设置10个网格大小的PML层厚度。为方便更加清晰显示波场信息，下面将使用灰度图显示不同时间步Ez波场的信息。

图6 低洼模型结构示意Fig.6 Schematic diagram of low-lying model structure

在多个激励源并发的情况下，采用相同主频的脉冲电磁波，此时在近地表区域，由单个激励源产生的球面波在空间上进行波场重构，子波波前形成平面电磁波束，如图7a,7b所示。平面电磁波束在地下空间扩散，在第400个时间步时遇到第二层介质产生绕射波，如图7c箭头指示位置；第500个时间步时，低洼模型底部产生的反射波如图7d标记所示，此时反射波的能量明显高于绕射波能量；两种吸收边界条件下，电磁波在介质区域内传播规律相同，有相同的波场快照；在第700个时间步时，平面电磁波束到达底部边界，基于递归卷积完全匹配层作为吸收边界条件的波场快照如图8a所示，此时平面电磁波束进入完全匹配层被有效地吸收掉，而以二阶Mur作为吸收边界条件在底部网格截断处发生了反射，部分波返回介质区域，如图8b圆圈位置，这将对能量相对较弱的绕射波造成了二次干扰，影响有效波场。

为了对比在不同偏移距下，两种吸收边界条件对多源并发电磁波的吸收效果，设置了宽角法观测低洼模型，将26个发射天线(激励源)置于低洼位置正上方，每个发射天线间隔两个空间步长，接收天线沿x轴移动，采集200道雷达数据。天线参数依然采用900 MHz主频的布莱克曼—哈里斯脉冲，多个发射天线无延时同时发射相同信号脉冲电磁波，得到在多源并发下，以递归卷积完全匹配层作为吸收边界条件的数据记录如图9a所示。此时，由于采用的是相同极化方向的激励信号，直达波、反射波与单个激励源雷达数据记录(如图10a所示)相比加大了子波时宽，增强了来自地下界面的反射信号，此时在直达波左、右两翼无明显畸变，无虚假反射。而Mur二阶吸收边界条件在大偏移距处有强烈的虚假反射记录；这是由于偏移距越大，波入射到匹配层的入射角度越大，掠射情况越严重。而在多个激励源的情况下，在大偏移距下，虚假反射会更加明显；这是由于加大了子波的时宽和介质的不均一性导致的色散造成的。如图9b箭头所指示的位置，在直达波的左右翼都发生了明显的波形畸变和虚假反射记录。虽然在单个天线时的雷达数据记录上，直达波左右两翼畸变不明显，但是也产生了明显的虚假反射信号，如图10b箭头所指示的位置。

图7 低洼模型不同时间步Ez波场快照Fig.7 Snapshots of the Ez wave field of the low-lying model at different time steps

a—递归卷积完全匹配层波场快照;b—Mur二阶吸收边界波场快照a—recursive convolution complete matching layer wave field snapshot;b—snapshot of Mur second-order absorbing boundary wave field图8 低洼模型700时间步Ez波场快照Fig.8 Snapshot of the low-lying model at 700 time step Ez wave field

a—递归卷积完全匹配层数据记录;b—Mur二阶吸收边界数据记录a—recursive convolution complete matching layer data record;b—Mur second-order absorbing boundary data record图9 多源并发雷达数据记录Fig.9 Multi-source concurrent radar data recording

3 结论与展望

1)在均匀介质下，Mur二阶吸收边界条件对多源并发下电磁波的吸收不完全，会在网格截断处产生反射回波返回波场模拟区域。而基于非分裂的递推卷积完全匹配层能很好地吸收多源并发下的电磁波，有效改善边界位置处对干涉后的电磁波吸收效果。

2)在复杂地质构造下，电磁波在地下介质中发生反射、绕射后以不同的入射角度到达网格截断处，传统的Mur二阶吸收边界条件在网格截断位置直接对微分方程进行近似求解，不能很好地消除强电磁反射，在边界网格截断位置处会有部分反射波返回研究区域对有效波造成二次干扰，特别是在大偏移距下，会造成直达波波形畸变，形成虚假反射记录；而基于非分裂递归卷积的完全匹配层作为吸收边界条件能很好地抑制直达波波形畸变现象，消除虚假反射记录。

3)基于非分裂递归卷积的完全匹配层作为吸收边界条件在不分裂波场的情况下，又避免了直接对卷积的求解的复杂运算，提升了计算的效率，改善了吸收的效果。

本文详细介绍了非分裂递归卷积完全匹配层和Mur二阶吸收边界结合时域有限差分的实现方法，并将它们运用到多个激励源并发下探地雷达正演模拟，基于C语言和matlab平台开发了相应的数值计算程序，为阵列探地雷达的数值模拟研究提供了帮助。

致谢:感谢煤炭资源与安全开采国家重点实验室提供的计算设备，感谢评审的专家提出的宝贵意见。