基于UPML邊界條件的LOD-FDTD　方法GPR波的數(shù)值模擬*

韓曉冰，魏海亮，李　奇

(西安科技大學(xué) 通信與信息工程學(xué)院，陜西 西安 710054)

?

基于UPML邊界條件的LOD-FDTD方法GPR波的數(shù)值模擬*

韓曉冰，魏海亮，李奇

(西安科技大學(xué) 通信與信息工程學(xué)院，陜西 西安 710054)

提出局部一維時域有限差分方法(LOD-FDTD)和單軸各向異性理想匹配層(UPML)邊界條件相結(jié)合的LOD-FDTD探地雷達(GPR)電磁波傳播數(shù)值模擬算法。通過將Maxwell方程離散化，導(dǎo)出三維GPR波LOD-FDTD算法迭代差分形式，引入UPML邊界條件、給出公式推導(dǎo)。在此基礎(chǔ)上，建立UPML邊界條件應(yīng)用模型、LOD-FDTD算法性能分析模型，驗證UPML邊界條件在數(shù)值模擬中對GPR波具有良好的吸收特性和算法在時間步不受CFL條件的限制下計算的高效性。建立LOD-FDTD方法的GPR波的正演模型，得到波場快照及其雷達剖面圖，為了解電磁波在復(fù)雜介質(zhì)中的傳播特性和雷達數(shù)據(jù)的解釋提供了依據(jù)。

GPR；UPML；LOD-FDTD方法；波場快照；雷達剖面圖

0　引　言

探地雷達(GroundPenetratingRadar)是用高頻無線電波來確定介質(zhì)內(nèi)部物質(zhì)分布規(guī)律的一種地球物理方法[1]。該方法具有高分辨率、高效率、無損探測、結(jié)果直觀等優(yōu)點，在工程建設(shè)、水文地質(zhì)調(diào)查、考古研究、礦產(chǎn)勘查、軍事等諸多領(lǐng)域已得到廣泛的應(yīng)用。

時域有限差分(FDTD)方法是研究電磁場在介質(zhì)中傳播的主要方法之一，F(xiàn)DTD法以其簡單靈活的特點得了廣泛的應(yīng)用。秦臻等[2]利用FDTD算法結(jié)合改進Mur吸收邊界條件進行了有耗介質(zhì)中的雷達正演模擬；劉四新等[3]利用FDTD算法結(jié)合改進的BerengerPML邊界條件實現(xiàn)了在有損耗頻散介質(zhì)中的正演模擬；馮晅等[4]在單軸各向異性理想匹配層(UPML)邊界條件下模擬了全極化探地雷達的正演模型；馮德山等[5]在UPML媒質(zhì)中利用交替方向隱式有限差分方法(ADI-FDTD)實現(xiàn)了GPR全波場的數(shù)值模擬。

FDTD方法進行探地雷達的數(shù)值模擬，擁有便于編程，計算效率高，計算精度高的優(yōu)勢。但是Courant-Friedrich-Levy(CFL)穩(wěn)定條件卻是限制該方法的最大問題。CFL穩(wěn)定條件限定了時間步長與空間步長的關(guān)系導(dǎo)致傳統(tǒng)的FDTD算法在解決復(fù)雜模型時，計算時間過長，影響計算效率。尋找更高效率與精度的算法是今后研究的一個重點。

無條件穩(wěn)定的隱式計算思想[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節(jié)推導(dǎo)出Maxwell方程組的LOD-FDTD算法的差分公式；第2節(jié)推導(dǎo)出UPML匹配層的LOD-FDTD算法的差分公式；第3節(jié)給出UPML吸收邊界模型、LOD-FDTD算法性能分析模型、復(fù)雜介質(zhì)電磁波正演模型；第4節(jié)給出結(jié)論。

1　LOD-FDTD差分格式

三維空間的Maxwell方程組可表示矩陣形式[11]

(1)

其中向量φ=[Ex,Ey,Ez,Hx,Hy,Hz)，矩陣A，B如文獻[12]中式(33)和式(34)所示。E為電場強度(V/m)，為磁場強度(A/m)。

(2)

離散化過程中對空間偏導(dǎo)數(shù)采用中心差分格式。n→n+1/2子時間步中Ex和Hz分量的計算公式如下所示，其它分量可做類似處理

(3)

(3)式都包含n+1/2時刻的電場或磁場值，首先隱式更新電(磁)場分量，然后顯示更新磁(電)分量。隱式更新過程如下所示。

(4)

轉(zhuǎn)化后等價于一個矩陣方程的求解，其中系數(shù)矩陣是一個帶狀對角矩陣，帶寬為2N-1隨階數(shù)升高而增大。

2　PML邊界條件

FDTD方法進行電磁波傳播的數(shù)值模擬時，吸收邊界條件是一個重要的環(huán)節(jié)，PML邊界條件包括BPML，CPML，UPML等邊界條件。BPML吸收層需要將場分量進行分裂處理，增加了計算的復(fù)雜度；CPML吸收層則對低頻分析有一定的優(yōu)越性；UPML吸收層不需要對場分量進行分裂處理，具有寬頻帶吸收特性，計算簡單，易于實現(xiàn)。針對文中的無條件穩(wěn)定的LOD-FDTD算法，選取與之相匹配的UPML邊界條件。

由電磁場與電磁波的理論可知，在單軸各向異性PML媒質(zhì)中Maxwell 2個旋度方程可以表示成下面的形式[13]

(5)

式中si=κi+σi/jwε0，i=x,y,z；si是有損耗單軸各向異性介質(zhì)在i方向的相對介電常數(shù)張量和相對導(dǎo)磁率張量；參數(shù)κi是用來吸收到達UPML層的凋落波；參數(shù)σi是PML區(qū)域的衰減因子。

LOD-FDTD差分方法應(yīng)用到UPML吸收層，2個子過程要求首先隱式更新3個電(磁)場分量，然后顯示更新3個磁(電)分量。以n→n+1/2子時間步加以說明。

由式(5)可以得到n→n+1/2子時間步中Ex和Hz相關(guān)的方程，為方便計算引入Dx=ε1szEx/sx,Bz=μ1syHz/sz，中間參量。方程如下所示

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

(10)

(11)

3　數(shù)值分析

模型一：UPML邊界條件應(yīng)用實例

為驗證UPML邊界條件的吸收效果，在三維真空空間中進行點偶極子輻射。FDTD計算區(qū)域為40×40×40個元胞，設(shè)垂直點偶極子位于計算區(qū)域的中心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條邊上產(chǎn)生了強的反射波。由圖2波場快照可知，雷達波等值線是非常規(guī)則的同心圓，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圖中也說明了均勻介質(zhì)中雷達波前以球面波輻射傳播的特點。

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

為進一步研究LOD-FDTD算法的計算精度和計算速度，在三維真空空間中進行點偶極子輻射。FDTD模擬區(qū)域為2 m×2 m×2 m，使用模型一的網(wǎng)格剖分方式，模擬區(qū)域四周被UPML吸收層包圍，UPML邊界設(shè)為10層網(wǎng)格?？疾煊^察點(20，30，20.5)的電場值Ez.

引入誤差模量

其中Ez,i是在iΔt時刻的電場分量Ez不同F(xiàn)DTD算法的仿真解，Ez,ref對應(yīng)時刻電場分量的解析式解。定義CFLN為CFL穩(wěn)定條件中所用的最大時間步長的倍數(shù)，圖3為傳統(tǒng)FDTD算法和不同CFLN條件下LOD-FDTD算法的誤差模量時間變化圖。圖3所示，LOD-FDTD算法的計算精度與傳統(tǒng)FDTD算法相比略有不足，且計算精度隨CFLN的增大而變小。這可理解為LOD-FDTD算法迭代時間增倍、迭代次數(shù)縮倍造成計算精度細微變差，但當(dāng)CFLN小于等于4時，LOD-FDTD算法仍可認為是一種高精度算法。

LOD-FDTD算法采用隱式計算思想具有無條件穩(wěn)定特性，突破了一般FDTD算法的CFL穩(wěn)定條件的限制，在基本保持計算精度的同時提高了計算的速度。

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

模型三：探地雷達正演模擬

目前探地雷達在地表發(fā)射高頻無線電波，并用接收天線在地面接收反射回來的的信號。采用發(fā)射、接收天線以零間距沿側(cè)線同步移動的剖面法探測方式來模擬探地雷達的的時間剖面掃描圖像。為了方便顯示，這里只考慮地下半空間的情況。

在三維空間中進行點偶極子輻射。FDTD模擬區(qū)域為16m×16m×16m，仍采用模型一的網(wǎng)格剖分方式。模擬區(qū)域四周被UPML吸收層包圍，UPML邊界設(shè)為10層網(wǎng)格。模擬區(qū)域是上層是素土層，下層是煤層，煤層成階梯狀分布，素土層相對介電常數(shù)為2，電導(dǎo)率為σ=2×10-8S/m，煤層相對介電常數(shù)為6，電導(dǎo)率為σ=2×10-8S/m.左方(3.5，8，6)m素土層和煤層交界處存在高2.5m寬0.5m的垂直真空裂縫，素土層中間方位(8.5，8，6)m存在的圓柱形空洞和位于右方(11.5，8，7)m金屬圓柱體，金屬圓柱體相對介電常數(shù)為8，電導(dǎo)率為σ=3.72×107S/m.模擬區(qū)域z方向的二維截面圖如圖3所示。輻射源采用模型一中的高斯脈沖，設(shè)測線位于深度3m處。

圖4　二維截面圖Fig.4　2-d section

圖5　500Δt時刻電場Ez快照Fig.5　Ez snapshots at 500Δt

圖5是500Δt時刻電場Ez二維快拍圖。由圖可見雷達波前以球面波輻射傳播，在依次遇到真空圓柱體、素土層煤層分界面、金屬圓柱體、真空裂縫后產(chǎn)較強的反射波后，雷達波前仍繼續(xù)向下輻射。

圖6　雷達波正演掃描圖Fig.6　Radar wave forward scan

圖6是模型的正演掃描圖。由圖顯示煤層的梯狀分布清晰可辯，垂直真空裂縫的上端和與分界面處都產(chǎn)生了較強的繞射波，2個圓柱體產(chǎn)生了較強的反射波，并觀察到在圓柱體上表面和素土層與煤層的分界面之間存在一強反射點。兩圓柱體的反射波是有區(qū)別的，金屬圓柱體產(chǎn)生的繞射波能量更強，幅度更大。

通過模擬雷達波在復(fù)雜介質(zhì)傳播中電場快照正演掃描圖，準(zhǔn)確地觀察了電磁波在介質(zhì)中的傳播規(guī)律，為今后的GPR資料解讀和圖像處理提供了依據(jù)。

4　結(jié)　論

1)為研究雷達電磁波的三維傳播特性，文中引入了無條件穩(wěn)定的LOD-FDTD算法，給出雷達電磁波的三維LOD-FDTD算法的迭代計算公式，指明其無條件穩(wěn)定的隱式計算思想；

2)在LOD-FDTD算法中成功引入UPML吸收邊界條件，且不破壞LOD-FDTD算法的無條件穩(wěn)定性，使模型免受截取邊界強反射的影響，使計算更加便捷；

3)設(shè)計了UPML吸收邊界模型、LOD-FDTD算法性能分析模型，驗證了UPML吸收邊界具有良好的吸收特性及LOD-FDTD算法克服常規(guī)FDTD算法受CFL條件的限制，可選取較大的時間步長，雖然需要迭代求解大型線性方程，仍然提高了計算的效率。設(shè)計了復(fù)雜介質(zhì)電磁波正演模型，進行模型的數(shù)值模擬，得到了復(fù)雜介質(zhì)模型電場快照和掃面圖，了解了電磁波在復(fù)雜介質(zhì)中的傳播特性，有助于地質(zhì)雷達資料的解釋。

References

[1]李大心.探地雷達方法與應(yīng)用[M].北京：地質(zhì)出版社，1994.

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

[2]秦臻，胡文寶.有耗介質(zhì)FDTD計算的改進的Mur吸收邊界條件[J].漢江石油學(xué)院學(xué)報，2004，26(2)：76-77.

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

[3]劉四新，曾昭發(fā).頻散介質(zhì)地質(zhì)雷達波傳播的數(shù)值模擬[J].地球物理學(xué)報，2007，50(1)：320-326.

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

[4]馮晅，鄒立龍，劉財.全極化探地雷達正演模擬[J].地球物理學(xué)報，2011，20(4)：49-357.

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

[5]馮德山，陳承申，戴前偉.基于UPML邊界條件的交替方向隱式有限差分法GPR全波場數(shù)值模擬[J].地球物理學(xué)報，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責(zé)任編輯：高佳

陜西教育科技研究計劃(145K1470)

韓曉冰(1968-)，男，陜西西安人，教授，E-mail:hanxiaobing@xust.edu.cn

P 631

A