基于UPML边界条件的LOD-FDTD 方法GPR波的数值模拟*

韩晓冰，魏海亮，李　奇

(西安科技大学 通信与信息工程学院，陕西 西安 710054)



基于UPML边界条件的LOD-FDTD方法GPR波的数值模拟*

韩晓冰，魏海亮，李奇

(西安科技大学 通信与信息工程学院，陕西 西安 710054)

提出局部一维时域有限差分方法(LOD-FDTD)和单轴各向异性理想匹配层(UPML)边界条件相结合的LOD-FDTD探地雷达(GPR)电磁波传播数值模拟算法。通过将Maxwell方程离散化，导出三维GPR波LOD-FDTD算法迭代差分形式，引入UPML边界条件、给出公式推导。在此基础上，建立UPML边界条件应用模型、LOD-FDTD算法性能分析模型，验证UPML边界条件在数值模拟中对GPR波具有良好的吸收特性和算法在时间步不受CFL条件的限制下计算的高效性。建立LOD-FDTD方法的GPR波的正演模型，得到波场快照及其雷达剖面图，为了解电磁波在复杂介质中的传播特性和雷达数据的解释提供了依据。

GPR；UPML；LOD-FDTD方法；波场快照；雷达剖面图

0　引　言

探地雷达(GroundPenetratingRadar)是用高频无线电波来确定介质内部物质分布规律的一种地球物理方法[1]。该方法具有高分辨率、高效率、无损探测、结果直观等优点，在工程建设、水文地质调查、考古研究、矿产勘查、军事等诸多领域已得到广泛的应用。

时域有限差分(FDTD)方法是研究电磁场在介质中传播的主要方法之一，FDTD法以其简单灵活的特点得了广泛的应用。秦臻等[2]利用FDTD算法结合改进Mur吸收边界条件进行了有耗介质中的雷达正演模拟；刘四新等[3]利用FDTD算法结合改进的BerengerPML边界条件实现了在有损耗频散介质中的正演模拟；冯晅等[4]在单轴各向异性理想匹配层(UPML)边界条件下模拟了全极化探地雷达的正演模型；冯德山等[5]在UPML媒质中利用交替方向隐式有限差分方法(ADI-FDTD)实现了GPR全波场的数值模拟。

FDTD方法进行探地雷达的数值模拟，拥有便于编程，计算效率高，计算精度高的优势。但是Courant-Friedrich-Levy(CFL)稳定条件却是限制该方法的最大问题。CFL稳定条件限定了时间步长与空间步长的关系导致传统的FDTD算法在解决复杂模型时，计算时间过长，影响计算效率。寻找更高效率与精度的算法是今后研究的一个重点。

无条件稳定的隐式计算思想[6]的ADI-FDTD[7]算法和LOD-FDTD[8-9]算法的时间步不受CFL条件的限制，把n→n+1的过程分成n→n+1/2和n+1/2→n+1 2个子时间过程。ADI-FDTD算法2个时间步分别采用向前和向后差分格式，误差互补；LOD-FDTD算法被证明与 ADI-FDTD方法有相同的计算精度，且计算时间比 ADI-FDTD方法减少55%[10]。

文中组织如下：第1节推导出Maxwell方程组的LOD-FDTD算法的差分公式；第2节推导出UPML匹配层的LOD-FDTD算法的差分公式；第3节给出UPML吸收边界模型、LOD-FDTD算法性能分析模型、复杂介质电磁波正演模型；第4节给出结论。

1　LOD-FDTD差分格式

三维空间的Maxwell方程组可表示矩阵形式[11]

(1)

其中向量φ=[Ex,Ey,Ez,Hx,Hy,Hz)，矩阵A，B如文献[12]中式(33)和式(34)所示。E为电场强度(V/m)，为磁场强度(A/m)。

(2)

离散化过程中对空间偏导数采用中心差分格式。n→n+1/2子时间步中Ex和Hz分量的计算公式如下所示，其它分量可做类似处理

(3)

(3)式都包含n+1/2时刻的电场或磁场值，首先隐式更新电(磁)场分量，然后显示更新磁(电)分量。隐式更新过程如下所示。

(4)

转化后等价于一个矩阵方程的求解，其中系数矩阵是一个带状对角矩阵，带宽为2N-1随阶数升高而增大。

2　PML边界条件

FDTD方法进行电磁波传播的数值模拟时，吸收边界条件是一个重要的环节，PML边界条件包括BPML，CPML，UPML等边界条件。BPML吸收层需要将场分量进行分裂处理，增加了计算的复杂度；CPML吸收层则对低频分析有一定的优越性；UPML吸收层不需要对场分量进行分裂处理，具有宽频带吸收特性，计算简单，易于实现。针对文中的无条件稳定的LOD-FDTD算法，选取与之相匹配的UPML边界条件。

由电磁场与电磁波的理论可知，在单轴各向异性PML媒质中Maxwell 2个旋度方程可以表示成下面的形式[13]

(5)

式中si=κi+σi/jwε0，i=x,y,z；si是有损耗单轴各向异性介质在i方向的相对介电常数张量和相对导磁率张量；参数κi是用来吸收到达UPML层的凋落波；参数σi是PML区域的衰减因子。

LOD-FDTD差分方法应用到UPML吸收层，2个子过程要求首先隐式更新3个电(磁)场分量，然后显示更新3个磁(电)分量。以n→n+1/2子时间步加以说明。

由式(5)可以得到n→n+1/2子时间步中Ex和Hz相关的方程，为方便计算引入Dx=ε1szEx/sx,Bz=μ1syHz/sz，中间参量。方程如下所示

离散以上公式，得到Ex和Hz的隐式差分方程

(10)

(11)

3　数值分析

模型一：UPML边界条件应用实例

为验证UPML边界条件的吸收效果，在三维真空空间中进行点偶极子辐射。FDTD计算区域为40×40×40个元胞，设垂直点偶极子位于计算区域的中心Ez(20，20，20.5)处，四周被PML吸收层包围。元胞尺寸为5 cm×5 cm×5 cm，即δ=5 cm，按照δ=2cΔt，时间间隔为Δt=83.333 ps.考察z=20δ平面的电场值Ez.计算中辐射源采用高斯脉冲

(12)

图1为未加入UPML边界条件Ez在z=20δ平面80Δt时时刻的波场快照，图2为加入UPML边界条件Ez在z=20δ平面80Δt时时刻的波场快照，图1可以看出无边界条件下波场分量在截断边界的4条边上产生了强的反射波。由图2波场快照可知，雷达波等值线是非常规则的同心圆，4个边、角没有明显的反射，表明UPML层的吸收效果很好雷达波无反射地进入完全匹配层后被迅速衰减掉。

图1　未加入边界条件波场快照Fig.1　Wave field snapshots without the boundary conditions

图2　加入边界条件的波场快照Fig.2　Wave field snapshots without the boundary conditions

对比图1，2表明：UPML层吸收效果良好，适合作为LOD-FDTD算法的边界条件。同时2图中也说明了均匀介质中雷达波前以球面波辐射传播的特点。

模型二：LOD-FDTD算法性能分析

为进一步研究LOD-FDTD算法的计算精度和计算速度，在三维真空空间中进行点偶极子辐射。FDTD模拟区域为2 m×2 m×2 m，使用模型一的网格剖分方式，模拟区域四周被UPML吸收层包围，UPML边界设为10层网格。考察观察点(20，30，20.5)的电场值Ez.

引入误差模量

其中Ez,i是在iΔt时刻的电场分量Ez不同FDTD算法的仿真解，Ez,ref对应时刻电场分量的解析式解。定义CFLN为CFL稳定条件中所用的最大时间步长的倍数，图3为传统FDTD算法和不同CFLN条件下LOD-FDTD算法的误差模量时间变化图。图3所示，LOD-FDTD算法的计算精度与传统FDTD算法相比略有不足，且计算精度随CFLN的增大而变小。这可理解为LOD-FDTD算法迭代时间增倍、迭代次数缩倍造成计算精度细微变差，但当CFLN小于等于4时，LOD-FDTD算法仍可认为是一种高精度算法。

LOD-FDTD算法采用隐式计算思想具有无条件稳定特性，突破了一般FDTD算法的CFL稳定条件的限制，在基本保持计算精度的同时提高了计算的速度。

图3　算法精度分析图Fig.3　Algorithm precision analysis

模型三：探地雷达正演模拟

目前探地雷达在地表发射高频无线电波，并用接收天线在地面接收反射回来的的信号。采用发射、接收天线以零间距沿侧线同步移动的剖面法探测方式来模拟探地雷达的的时间剖面扫描图像。为了方便显示，这里只考虑地下半空间的情况。

在三维空间中进行点偶极子辐射。FDTD模拟区域为16m×16m×16m，仍采用模型一的网格剖分方式。模拟区域四周被UPML吸收层包围，UPML边界设为10层网格。模拟区域是上层是素土层，下层是煤层，煤层成阶梯状分布，素土层相对介电常数为2，电导率为σ=2×10-8S/m，煤层相对介电常数为6，电导率为σ=2×10-8S/m.左方(3.5，8，6)m素土层和煤层交界处存在高2.5m宽0.5m的垂直真空裂缝，素土层中间方位(8.5，8，6)m存在的圆柱形空洞和位于右方(11.5，8，7)m金属圆柱体，金属圆柱体相对介电常数为8，电导率为σ=3.72×107S/m.模拟区域z方向的二维截面图如图3所示。辐射源采用模型一中的高斯脉冲，设测线位于深度3m处。

图4　二维截面图Fig.4　2-d section

图5　500Δt时刻电场Ez快照Fig.5　Ez snapshots at 500Δt

图5是500Δt时刻电场Ez二维快拍图。由图可见雷达波前以球面波辐射传播，在依次遇到真空圆柱体、素土层煤层分界面、金属圆柱体、真空裂缝后产较强的反射波后，雷达波前仍继续向下辐射。

图6　雷达波正演扫描图Fig.6　Radar wave forward scan

图6是模型的正演扫描图。由图显示煤层的梯状分布清晰可辩，垂直真空裂缝的上端和与分界面处都产生了较强的绕射波，2个圆柱体产生了较强的反射波，并观察到在圆柱体上表面和素土层与煤层的分界面之间存在一强反射点。两圆柱体的反射波是有区别的，金属圆柱体产生的绕射波能量更强，幅度更大。

通过模拟雷达波在复杂介质传播中电场快照正演扫描图，准确地观察了电磁波在介质中的传播规律，为今后的GPR资料解读和图像处理提供了依据。

4　结　论

1)为研究雷达电磁波的三维传播特性，文中引入了无条件稳定的LOD-FDTD算法，给出雷达电磁波的三维LOD-FDTD算法的迭代计算公式，指明其无条件稳定的隐式计算思想；

2)在LOD-FDTD算法中成功引入UPML吸收边界条件，且不破坏LOD-FDTD算法的无条件稳定性，使模型免受截取边界强反射的影响，使计算更加便捷；

3)设计了UPML吸收边界模型、LOD-FDTD算法性能分析模型，验证了UPML吸收边界具有良好的吸收特性及LOD-FDTD算法克服常规FDTD算法受CFL条件的限制，可选取较大的时间步长，虽然需要迭代求解大型线性方程，仍然提高了计算的效率。设计了复杂介质电磁波正演模型，进行模型的数值模拟，得到了复杂介质模型电场快照和扫面图，了解了电磁波在复杂介质中的传播特性，有助于地质雷达资料的解释。

References

[1]李大心.探地雷达方法与应用[M].北京：地质出版社，1994.

LIDa-xin.Groundpenetratingradar(GPR)methodandapplication[M].Beijing：GeologicalPress，1994.

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

QINZhen，HUWen-bao.ModifiedMurabsorbingboundaryconditionsforlossymediumFDTDcalculation[J].JournalofJianghanPetroleumInstitute，2004，26(2)：76-77.

[3]刘四新，曾昭发.频散介质地质雷达波传播的数值模拟[J].地球物理学报，2007，50(1)：320-326.

LIUSi-xin，ZENGZhao-fa.Numericalsimulationofdispersionmediumgeologicalradarwavetransmission[J].JournalofGeophysics，2007，50(1)：320-326.

[4]冯晅，邹立龙，刘财.全极化探地雷达正演模拟[J].地球物理学报，2011，20(4)：49-357.

FENGXuan，ZOULi-long，LIUCai.Allgeochemicalpenetratingradarforwardmodeling[J].JournalofGeophysics，2011，20(4)：49-357.

[5]冯德山，陈承申，戴前伟.基于UPML边界条件的交替方向隐式有限差分法GPR全波场数值模拟[J].地球物理学报，2010，53(10)：2 484-2 496.

FENGDe-shan，CHENCheng-shen，DAIQian-wei.BasedonUPMLboundaryconditionsofalternatingdirectionimplicitfinitedifferencemethod，GPRfullwavefieldnumericalsimulation[J].JournalofGeophysics，2010，53(10)：2 484-2 496.

[6]ShibayamaJ，MurakiM，YamauchiJ，etal．EfficientimplicitFDTDalgorithmbasedonlocallyone-dimensionalscheme[J].ElectronicsLetters，2005，41(19)： 1 046-1 047．

[7]KONGYong-dan，CHUQing-xin，LIRong-lin.Efficientunconditionallystableone-stepleapfrogADI-FDTDmethodwithlownumericaldispersion[J].IETMicrowavesAntennasandPropagation，2014，8(5)：337-345.

[8]RanaMM，MohanAS.CentreforHealthTechnol.ConvolutionalperfectlymatchedlayerABCfor3-DLOD-FDTDusingfundamental[J].IEEEMicrowaveandWirelessComponentsLetters，2013，23(8)：388-399.

[9]WANGJian-bao，ZHOUBi-hua.Efficiency-improvedLOD-FDTDmethodforsolvingperiodicstructuresatobliqueincidence[J].IEEEMicrowaveandWirelessComponentsLetters，2013，23(8)：521-523.

[10]LIUQi-feng，CHENZhi-zhang，YINWen-yan.Anefficiencyunconditionallystablethree-simensionalLOD-FDTDmethod[J].IEEEMTT-SInternational，MicrowaveSymposiumDigest，2008：45-48.

[11]LiuQF，ChenZZ，YinWY．Anefficientunconditionallystablethree-dimensionalLOD-FDTDmethod[C]//Proceedingsof2008IEEEMTT-SInternationalMicrowaveSymposiumDigest.Atlanta，GA：IEEE，2008：45-48．

[12]TanEL.Fundamentalschemesforefficientunconditionallystableimplicitfinite-differenceTime-domainmethods[J].IEEETransactionsonAntennasandPropagation，2008，56(1)：170-177.

[13]葛德彪，闰玉波．电磁波时域有限差分方法[M].西安：西安电子科技出版社,2011.

GEDe-biao，RUNYu-bo.Electromagneticfinitedifferencetimedomainmethod[M].Xi’an：Xi’anElectronicScienceandTechnologyPress，2011.

GPR wave numerical simulation based on the UPML boundary conditions of LOD-FDTD method

HAN Xiao-bing，WEI Hai-liang，LI Qi

(CollegeofCommunicationandInformationEngineering,Xi’anUniversityofScienceandTechnology,Xi’an710054,China)

ForsimulatingthetransmissionofGroundPenetratingRadar’selectromagneticwaveincomplexmedium，aLOD-FDTDalgorithmofthecombinationofLOD-FDTDandUPMLboundaryconditionsisputforward.Throughthediscretizationofthethree-dimentionalMaxwell’sequation，sub-timestepiterativefinitedifferentialformofLOD-FDTDalgorithmisderived.IntroducingUPMLboundaryconditionsanditsformuladerivation，thentheUPMLboundaryconditionapplicationmodelandLOD-FDTDalgorithmprecisionanalysismodelareestablishedtoprovetheabsorbingvalidationofUPMLboundaryconditioninnumericalsimulationofGPRwaveandtheefficiencyofthealgorithmunderthetimestepbeingnotrestrictedbyCFLcondition.BasedonestablishingtheforwardmodelofGPRwavebyLOD-FDTDmethod，wavefieldsnapandshotradarprofilearecaptured，whichcanbeusedtobetterunderstandthetransmissioncharacteristicsofelectromagneticwaveincomplexmediumandtheradardatainterpretation.

GPR；UPML；LOD-FDTDalgorithm；wavefieldsnap；radarprofile

10.13800/j.cnki.xakjdxxb.2016.0516

1672-9315(2016)05-0703-06

2015-09-24责任编辑：高佳

陕西教育科技研究计划(145K1470)

韩晓冰(1968-)，男，陕西西安人，教授，E-mail:hanxiaobing@xust.edu.cn

P 631

A