复频率参数完全匹配层吸收边界在瞬变电磁法正演中的应用

李展辉，黄清华

北京大学地球与空间科学学院地球物理学系，北京 100871

1 引言

时域有限差分（Finite difference time domain，FDTD）是瞬变电磁法（Transient electromagnetic，TEM）正演常用方法之一.Goldman和Stoyer（1983）首次使用隐式FDTD模拟三维轴对称模型，随后Oristaglio和Hohmann（1984）使用修改的Du Fort-Frankel格式模拟了TEM二维模型，Wang和Hohmann（1993）在Du Fort-Frankel格式的基础上发展了一种TEM三维FDTD正演方法，为后续部分学者借用（孙怀凤等，2013；许洋铖等，2012；Commer and Newman，2004；Endo and Noguchi，2002）.FDTD对TEM的三维正演是在有限模型空间内模拟能扩散到更大区域的瞬变场，因而需要引入截断边界条件.正演中通常采用狄利克雷边界条件（Dirichlet boundary condition，DBC）（Wang and Hohmann，1993），并设置边界足够远.DBC在扩散场未抵达或刚抵达边界时效果较好，但随着延迟时间的推移，反射场将对原扩散场产生较大的干扰.要消除这种干扰，则需要更大的模型，这将增加内存和计算量.因而一些学者开始研究其他的吸收边界来代替DBC，比如Mur吸收边界（秦臻和胡文宝，2004；姬金祖，2008；岳建华，2007）以及在其基础上修改而来的新吸收边界条件（杨海燕和岳建华，2009）.但是他们并没有对吸收边界的效果作详细的讨论，也没有给出与理论解的对比.李墩柱（2010）将电磁波正演中应用最普遍的完全匹配层（Perfectly Matched Layer，PML）吸收边界（Berenger，1994，1996；李展辉等，2009；Chew and Weedon，1994）应用到TEM半空间模型正演中，然而文中模型的水平范围为［-3000m，3000m］，半空间电导率为0.01S·m-1，但是最长延迟时间仅为8ms，涡流电场形成的“烟圈”还远未到达边界，晚期吸收效果并不明确，在参数设置上也仅仅说明需要多次尝试.本文深入研究了PML在TEM正演中的吸收效果，并最终采用了其中的一种——复频率参数完全匹配层（Complex Frequency Shifted Perfectly Matched Layer，CFS-PML）（Kuzuoglu and Mittra，1996；Roden and Gedney，2000；姜永金等，2004；Berenger，2012）.首先通过推导扩散场在常规PML内部的平面波解的形式，指出了常规PML在TEM正演中失效的原因，然后给出了扩散场在CFS-PML内部的平面波解，研究了CFS-PML各项参数对吸收效果的影响，并给出了CFS-PML的参数选择方案.最后将CFSPML应用到全空间和半空间模型中以验证CFS-PML的有效性.同时为了提高半空间模型的正演效率，本文采用了李展辉和黄清华（2011）提出的新方法，而全空间模型正演则采用普遍使用的Du Fort-Frankel格式.

2 理论推导

2.1 含PML的Maxwell频域方程组

频率域带PML的旋度方程可以表示为（Chew and Weedon，1994）

其中

常规PML中（Chew and Weedon，1994），

其中κu为网格延拓因子，由于本文采用非等距网格，κu可直接包含于网格中.因此在下文中所有的κu值都设置为1.令ε′＝εΔ-jσ／ω，可得

2.2 常规PML和CFS-PML的吸收性质

设式（5）的平面波解为

定义复波速

可得到频散方程

对于常规PML，式（6）可写为

其中

而TEM的低频近似满足

那么复波速可以近似为

将式（13）代入到式（10）中，并且考虑在x方向吸收，设置 （σx，σy，σz）＝ （σx，0，0），可得

设

为PML的吸收系数，APML越小，代表吸收越强烈.从式（15）中可以看出，随着角频率ω的减小，APML也会随着减小，而当ω→0时，亦有APML→0.这意味着PML对低频场的吸收随着频率的降低而更为强烈，过强的吸收在数值计算中会引起大的异常（Berenger，2000，1996）.

对于CFS-PML，经过类似的推导之后可以得到

从式（16）可以看出，αx的出现使得CFS-PML避免了在ω接近0时出现吸收非常强烈的现象.

设置CFS-PML的一个参考频率：

把式（17）代入式（16），得

当f≫fαx时，式（18）近似为式（14），CFS-PML退化成常规PML；当f≪fαx时，式（18）近似为

这仅仅是x方向的网格延拓.

瞬变场的优势频率会随着时间逐渐变低.当优势频率足够低并进入PML内部时，PML会对这些优势频率产生过强的吸收进而使整个瞬变场产生异常.而CFS-PML能避免这种情况，但是σx／αx这个比例非常重要，从式（19）可以看出，对于f≪fαx频率部分只是一个网格延拓因子，如果过大同样会引起数值异常.为此我们可以设置另一个参考频率：

2.3 含CFS-PML的Maxwell时域方程组

时域方程组可以用卷积的形式表达出来（Roden and Gedney，2000），以Ex为例：

其中¯sy（t）和¯sz（t）分别是1／sy（ω）和1／sz（ω）的拉普拉斯逆变换：

其中δ（t）是狄拉克函数，ζu（t）可以表述如下：

其中u（t）为单位阶跃函数.那么最终式可以表示为

其中

其他分量可以类似地表达出来.

2.4 CFS-PML参数选取研究

其中L为PML的厚度.假设场沿x方向传播，X＝x.定义CFS-PML吸收系数为

式中

式（28）中存在两个可控参数：fαx和L.吸收系数A越小代表额外吸收越强烈.很明显，吸收系数A是CFS-PML厚度L的递减函数，L越大，吸收越强烈.下面将对参数fαx选取作具体的分析.

假设介质是均匀的，且电导率为σ＝1×10-3S·m-1，场源到PML内边界的距离为L1＝1400m，PML厚度为L2＝600m，关注的最大时间为tmax＝1ms，那么抵达边界并穿过PML返回场内的最低频率为

图1 不同参考频率下CFS-PML对扩散场的额外吸收系数比较Fig.1 The additional absorbing coefficient of CFS-PML under different reference frequencies

3 数值模拟结果及讨论

本文除巷道带异常体模型外，所有正演均采用相同的网格：网格数量为128×128×128，中央网格大小为Δxmin＝Δymin＝Δzmin＝10m，从中央到四周逐渐增加网格步长，最大处为Δxmax＝10Δxmin，Δymax＝10Δymin，Δzmax＝7Δzmin.图2描述了网格的主要特征.CFS-PML的层数为8层，已经包含在上述网格中.

图2 瞬变电磁三维正演网格示意图Fig.2 A sketch map of the mesh for TEM 3Dforward modeling

3.1 均匀全空间

为了方便与理论解对比，我们采用了均匀全空间模型.同时为了探讨CFS-PML对不同介质的适应性，我们选取了两种不同的介质：σ＝1×10-2S·m-1和σ＝1×10-3S·m-1.以水平正方形线框作为发射线框，线框边长L＝70m，垂直方向位于z＝0处，水平方向位于x-y平面的正中央.采用脉冲形式的发射电流，最终比对垂直磁场Hz.

选择的测线为z＝0和y＝0平面的交线.我们定义在CFS-PML与正常区域交接处Hz易号的时间点为早期、晚期的分界点.首先对σ＝1×10-2S·m-1的三个早期时间点进行对比，根据牛之链（2007）中式（2-33）计算得到早晚期分界延迟时间约为4.5ms，那么首先选取t＝0.35ms，1.23ms和2.64ms三个时间点.图3为在上述三个时间点测线上分别使用CFS-PML、DBC得到的Hz值与理论值的对比.可以看出，在早期三者几乎一致.这是因为扩散场尚未或者刚刚抵达边界，边界条件的区别未能体现出来.因此在下文中，本文将着重扩散晚期的对比.

图4为t＝6.80ms，12.50ms，和20.00ms时分别使用CFS-PML，DBC得到的Hz值和理论解的对比.可以看出t＝6.8ms时使用DBC的结果在边界处已经出现了较大的异常，随后这种异常开始扩散到场的中央，结果越来越偏离理论解.而使用CFS-PML吸收边界时，不管在中央还是在边界，Hz值依然和理论解符合得很好.CFS-PML的吸收作用在晚期非常显著.

图5为σ＝1×10-3S·m-1情况下不同时刻使用CFS-PML、DBC得到的Hz值和理论解之间的对比，同时加入了CFS-PML中设置fαx＝∞所得的结果.σ＝1×10-3S·m-1模型中场的扩散速度要明显快于σ＝1×10-2S·m-1的情况，但是两种情况下的场在一定的条件下具有相似性（王华军，2008）.根据这种相似性，t＝2ms时场的形态类似于σ＝1×10-2S·m-1时t＝20ms情况下场的形态，因此本文将σ＝1×10-3S·m-1情况下最晚的一个时间点设定为t＝2ms.从图5可以看出，在晚期CFS-PML的优势依然非常明显，而且经过适当设置后的CFS-PML也明显优于设置fαx＝∞的CFS-PML.

图3 全空间σ＝1×10-2S·m-1模型中，t＝0.35ms，1.23ms和2.64ms时使用CFS-PML吸收边界、DBC的Hz分量和理论解的对比Fig.3 The comparison among the Hzcomponents obtained from CFS-PML，DBC，and the analytical solutions at the times of 0.35ms，1.23ms，and 2.64ms of a whole space model withσ＝1×10-2S·m-1

图4 全空间σ＝1×10-2S·m-1模型中，t＝6.80ms，12.50ms和20.00ms时使用CFS-PML吸收边界条件、DBC的Hz分量和理论解的对比Fig.4 The comparison among the Hzcomponents obtained from CFS-PML，DBC，and the analytical solutions at the times of 6.80ms，12.50ms，and 20.00msof a whole space model withσ＝1×1 0-2 S·m-1

图5 全空间σ＝1×10-3S·m-1模型中，t＝0.680ms和2.00ms时使用CFS-PML吸收边界条件、DBC的Hz分量和理论解的对比，对比中还加入了在CFS-PML中设置fαx＝ ∞ 所得的Hz分量Fig.5 The comparison among the Hzcomponents obtained from CFS-PML，DBC，and the analytical solutions at the times of 0.68ms and 2.00ms of a whole space model withσ＝1×10-3S·m-1.The Hz component obtained from CFS-PML with fαx＝ ∞is also included for reference

为显示CFS-PML内部场的情况，我们将σ＝1×10-3S·m-1模型中2ms时CFS-PML内外的场都描绘出来并与理论解进行对比，如图6所示.可以看出，CFS-PML内部的场与理论解截然不同.CFS-PML内部出现了Hz的正负变化，意味着涡流的存在，而理论解并没有Hz符号的变化.CFS-PML虽然不能将扩散场完全吸收，但是其在吸收的同时以某种方式将涡流限制在了CFS-PML内部，使得正常区域的扩散场不受影响.

3.2 均匀半空间

本文采用了李展辉和黄清华（2011）介绍的方法在半空间模型中引入空气层.由于空气层中波速远大于地下场的扩散速度，由式（10）可知，针对地下介质设计的CFS-PML在空气层中对波的吸收基本为零，因此在空气层中需要在CFS-PML介质内设置一定的电导率进行吸收.半空间模型中，地下介质的电导率设置为σ＝1×10-2S·m-1，空气层中电导率设置为1×10-6S·m-1，并在靠近边界处逐渐增加.地空边界设置在z＝0处，激励源的位置与形式均与全空间模型一致.

图7显示了使用CFS-PML吸收边界、DBC所得到的Hz分量和理论解之间的对比.从图中可以看出，t＝6.80ms时使用DBC的Hz分量在两端已经偏离了理论解.随着时间的推移，误差越来越大，在20ms时，误差已经在一个量级以上.而使用CFS-PML的Hz与理论解的吻合度大幅优于使用DBC的情况，但是也有一定的误差，而且这种误差也随时间的推移缓慢地增加.这可能是CFS-PML不能对空气层中的场产生有效的吸收造成的，越到晚期越明显.图8详细展示了线圈中心Hz相对误差随时间变化.从图中可以看出，线圈中心相对误差随着时间逐渐增加，在20ms时误差达到了13%.实际应用中需要根据精度要求调整模型的大小以控制相对误差.

图6 全空间σ＝1×10-3S·m-1模型CFS-PML内外Hz分量展示，并与理论解比较以体现CFS-PML内部场特征Fig.6 The Hzcomponent in and out the CFS-PML of a whole space model withσ＝1×10-3S·m-1 at 2ms

3.3 巷道模型

全空间瞬变电磁法典型的应用之一为煤矿巷道内地下水的探测（李宇等，2012；于景 等，2011；郭纯等，2006）.模型如图9所示.地下水所在的位置电导率较高，因此模型中设定围岩电导率为0.05S·m-1，异常体电导率为0.5S·m-1.巷道内通常很狭小，因此发射线圈设定为3m×3m，并采用重叠回线的方式记录感应电压.由于该模型不存在理论解，为了获取类似理论解的参考解，我们将原网格x，y方向各增加60个网格，在z方向增加40个网格.总网格数量达到188×188×168，以保证在我们关心的时间窗内，边界反射场没有进入观测区域.考虑到模型尺寸问题，最大的延迟时间设置为0.6ms，因为根据扩散深度公式（牛之琏，2007）

图7 半空间模型中，t＝6.80ms，12.50ms和20.00ms时使用CFS-PML吸收边界条件、DBC的Hz分量和理论解的对比Fig.7 The comparisons among the Hzcomponents obtained from CFS-PML，DBC，and the analytical solutions at the times of t＝6.80ms，12.50ms，and20.00ms of a half space model

图8 回线中心Hz的相对误差随时间分布图Fig.8 The relative errors of Hzvarying over time at the center of the rectangular loop

计算可得此时的扩散深度为124m，已经超出了本模型范围.式（31）中ρ为电阻率，t为延迟时间.结果如图10所示.从图中可以看出CFS-PML的结果与参考解一直符合得很好，最终的相对误差不超过6%.

图9 巷道内含异常体模型，模型关于y＝0对称Fig.9 Underground roadway model with an anomalous body，which is symmetric about the y＝0plane

图10 巷道模型的重叠回线测量垂直感应电动势结果图，左下角的小图为相对于参考解的相对误差Fig.10 Vertical emf of the underground roadway model.The relative errors to the reference solution are shown in the embedded figure

4 结论

本文将CFS-PML应用到TEM三维FDTD正演中，通过适当地设置CFS-PML参数保证了数值计算的稳定性和良好的吸收效果.全空间和半空间数值结果表明CFS-PML在瞬变场早期与传统的狄利克雷边界条件的结果一致，晚期则明显优于狄利克雷边界条件的结果，并且能有效地节省计算时间和存储空间.但是针对地下扩散场设计的CFSPML并不能有效地吸收空气中的电磁场，导致半空间模型结果到后期也产生一定的误差，并且随着时间缓慢的增长.实际应用中需要根据TEM正演的最大延迟时间以及对误差的容忍度适当调整模型的大小，以确保误差在要求之内.本文最后将CFSPML应用到与实际情况贴近的巷道模型中，采用重叠回线测量感应电动势，并与参考解进行对比和误差分析，再一次证明了CFS-PML的有效性.

Bérenger J P.1994.A perfectly matched layer for the absorption of electromagnetic waves.Journal of Computational Physics，114（2）：185-200.

Bérenger J P.1996.Three-dimensional perfectly matched layer for the absorption of electromagnetic waves.Journal of Computational Physics，127（2）：363-379.

Bérenger J P.1996.Perfectly matched layer for the FDTD solution of wave-structure interaction problems.IEEE Transactions on Antennas and Propagation，44（1）：110-117.

Bérenger J P.2000.Numerical reflection of evanescent waves by PMLs：origin and interpretation in the FDTD case.Expected consequences to other finite methods.Int.J.Numer.Model，13（2-3）：103-114.

Bérenger J P.2012.An optimized CFS-PML for wave-structure interaction problems.IEEE Transactions on Electromagnetic Compatibility，54（2）：351-358.

Chew W C，Weedon W H.1994.A 3Dperfectly matched medium from modified Maxwell′s equations with stretched coordinates.Microwave and Optical Technology Letters，7（13）：599-604.

Commer M，Newman G.2004.A parallel finite-difference approach for 3Dtransient electromagnetic modeling with galvanic sources.Geophysics，69（5）：1192-1202.

Endo M，Noguchi K.2002.Three-dimensional modeling considering the topography for the case of the time-domain electromagnetic method.Amsterdam：Elsevier Science BV，35：85-107.

Goldman M M，Stoyer C H.1983.Finite-difference calculations of the transient field of an axially-symmetric earth for vertical magnetic dipole excitation.Geophysics，48（7）：953-963.

Guo C，Liu B Z，Bai D H.2006.Prediction of water disasters ahead of tunneling in coal mine using continuous detection by UWTEM.Seismology and Geology（in Chinese），28（3）：456-462.

Ji J Z，Gao X，Liu Z H.2008.Simplified format application of 2D two-rank Mur absorbing boundary condition.Journal of Shenyang Institute of Aeronautical Engineering（in Chinese），25（2）：10-13.

Jiang Y J，Mao J J，Cai S L.2004.Implementation and analysis of the perfectly matched layer media with CFS for the PSTD method.Journal of Microwaves（in Chinese），20（4）：36-39.

Kuzuoglu M，Mittra R.1996.Frequency dependence of the constitutive parameters of causal perfectly matched anisotropic absorbers.IEEE Microwave and Guided Wave Letters，6（12）：447-449.

Li D Z.2010.Numerical Simulation of Transient Electromagnetic Fields（in Chinese）［Master′s thesis］.Beijing：Peking University.

Li Y，Tong B，Guo C F，etal.2012.Study on application and denoising of mine full space transient electromagnetic method.Coal Engineering（in Chinese），（11）：29-31.

Li Z H，Huang Q H，Wang Y B.2009.A 3-D staggered grid PSTD method for borehole radar simulations in dispersive media.Chinese J.Geophys.（in Chinese），52（7）：1915-1922.

Li Z H，Huang Q H.2011.A 3Dcurvilinear staggered-grid finite difference time-domain method with an impulse source for transient electromagnetic method.Chinese Geophysics（in Chinese），240.

Niu Z L.2007.Principle of Time Domain Electromagnetic Method（in Chinese）.Changsha：Central South University Press.

Oristaglio M L，Hohmann G W.1984.Diffusion of electromagnetic fields into a two-dimensional earth-a finite-difference approach.Geophysics，49（7）：870-894.

Qin Z，Hu W B.2004.Modified Mur absorbing boundary conditions for lossy medium FDTD calculation.Journal of Jianghan Petroleum Institute（in Chinese），26（2）：76-77.

Roden J A，Gedney S D.2000.Convolution PML（CPML）：An efficient FDTD implementation of the CFS-PML for arbitrary media.Microwave and Optical Technology Letters，27（5）：334-339.

Sun H F，Li X，Li S C，etal.2013.Three-dimensional FDTD modeling of TEM excited by a loop source considering ramp time.Chinese J.Geophys.（in Chinese），56（3）：1049-1064.

Wang H J.2008.Time domain transient electromagnetism all time apparent resistivity translation algorithm.Chinese J.Geophys.（in Chinese），51（6）：1936-1942.

Wang T，Hohmann G W.1993.A finite-difference，time-domain solution for 3-dimensional electromagnetic modeling.Geophysics，58（6）：797-809.

Xu Y C，Lin J，Li S Y，etal.2012.Calculation of full-waveform airborne electromagnetic response with three-dimension finitedifference solution in time-domain.Chinese J.Geophys.（in Chinese），55（6）：2105-2114.

Yang H Y，Yue J H.2009.Application of absorbing boundary condition in whole-space computation of transient electromagnetic response.Journal of China University of Mining ＆Technology（in Chinese），38（2）：263-268.

Yu J C，Liu Z Q，Liao J J，etal.2011.Application of full space transient electromagnetic method to mine water prevention and control.Coal Science and Technology（in Chinese），39（9）：110-113.

Yue J H，Yang H Y，Hu B.2007.3Dfinite difference time domain numerical simulation for TEM in-mine.Progress in Geophysics（in Chinese），22（6）：1904-1909.

附中文参考文献

郭纯，刘白宙，白登海.2006.地下全空间瞬变电磁技术在煤矿巷道掘进头的连续跟踪超前探测.地震地质，28（3）：456-462.

姬金祖，高旭，刘战合.2008.二维Mur二阶吸收边界条件简化格式的应用.沈阳航空工业学院学报，25（2）：10-13.

姜永金，毛钧杰，柴舜连.2004.CFS-PML边界条件在PSTD算法中的实现与性能分析.微波学报，20（4）：36-39.

李墩柱.2010.地形起伏下的瞬变电磁场数值模拟［硕士论文］.北京：北京大学.

李宇，童兵，郭昌放等.2012.矿井全空间瞬变电磁法应用及其消噪研究.煤炭工程，（11）：29-31.

李展辉，黄清华，王彦宾.2009.三维错格时域伪谱法在频散介质井中雷达模拟中的应用.地球物理学报，52（7）：1915-1922.

李展辉，黄清华.2011.瞬变电磁正演之带激励源的三维曲线交错网格时域有限差分方法.中国地球物理，240.

牛之琏.2007.时间域电磁法原理.长沙：中南大学出版社.

秦臻，胡文宝.2004.用于有耗介质FDTD计算的改进的Mur吸收边界条件.江汉石油学院学报，26（2）：76-77.

孙怀凤，李貅，李术才等.2013.考虑关断时间的回线源激发TEM三维时域有限差分正演.地球物理学报，56（3）：1049-1064.

王华军.2008.时间域瞬变电磁法全区视电阻率的平移算法.地球物理学报，51（6）：1936-1942.

许洋铖，林君，李肃义等.2012.全波形时间域航空电磁响应三维有限差分数值计算.地球物理学报，56（6）：2105-2114.

杨海燕，岳建华.2009.吸收边界条件在全空间瞬变电磁计算中的应用.中国矿业大学学报，38（2）：263-268.

于景邨，刘振庆，廖俊杰等.2011.全空间瞬变电磁法在煤矿防治水中的应用.煤炭科学技术，39（9）：110-113.

岳建华，杨海燕，胡搏.2007.矿井瞬变电磁法三维时域有限差分数值模拟.地球物理学进展，22（6）：1904-1909.