方宏遠，李 健，鐘燕輝，王復(fù)明

(1．鄭州大學(xué) 水利與環(huán)境學(xué)院，河南 鄭州450001;2．水利與交通基礎(chǔ)設(shè)施安全防護河南省協(xié)同創(chuàng)新中心，河南 鄭州450001)

0 引言

近年來，探地雷達憑借快速、連續(xù)、無破損等優(yōu)點已經(jīng)廣泛應(yīng)用于道路工程質(zhì)量檢測中． 精確測得路面結(jié)構(gòu)各層材料的介電參數(shù)，是保證檢測質(zhì)量的關(guān)鍵之一． 獲取材料介電參數(shù)的方法一般有兩種:①檢測法，即選取代表點鉆芯取樣，測試其介電參數(shù)，以檢測結(jié)果作為一段距離內(nèi)道路材料的介電參數(shù)，這種以點代面的方法不能反映整體路面材料介電特性;②反演方法，即基于探地雷達實測回波信號，利用優(yōu)化算法反算路面結(jié)構(gòu)層參數(shù)，這種方法檢測精度高，廣泛應(yīng)用于探地雷達檢測中．基于實測信號反演路面結(jié)構(gòu)層參數(shù)的關(guān)鍵步驟之一就是建立探地雷達電磁波在路面結(jié)構(gòu)層中傳播的精確高效數(shù)值模型． 傳統(tǒng)模擬探地雷達電磁波在地下結(jié)構(gòu)中傳播的方法有很多，主要包括傳遞矩陣方法［1］、射線追蹤方法［2］、時域偽譜方法［3］、無網(wǎng)格方法［4］、時域有限差分方法(FDTD)［5］、ADI-FDTD 方法［6］等． 這些方法雖然都可以模擬探地雷達電磁波在層狀結(jié)構(gòu)中的傳播，但也存在一些問題．比如傳遞矩陣方法對于有耗介質(zhì)容易出現(xiàn)指數(shù)溢出問題，時域有限差分方法受到CFL 穩(wěn)定條件的限制，計算效率受到制約，ADI-FDTD 方法雖然打破了CFL 條件的制約，但過大的時間步長會導(dǎo)致數(shù)值色散性顯著增加．基于邊值問題的精細積分方法是一種無條件穩(wěn)定的高精度數(shù)值算法，其數(shù)值精度可以達到計算機意義上的解析解，且非常適合求解層狀結(jié)構(gòu)中波動問題［7-8］．筆者采用精細積分方法建立探地雷達電磁波在路面結(jié)構(gòu)層中傳播的數(shù)值模型． 并通過與FDTD 方法以及實測回波信號進行對比，驗證筆者所建數(shù)值模型的精度和效率．

1 控制方程

探地雷達電磁波在層狀結(jié)構(gòu)中傳播過程如圖1 所示，發(fā)射天線向路面結(jié)構(gòu)層發(fā)射電磁脈沖，電磁波在路面結(jié)構(gòu)中遇到介質(zhì)分界面發(fā)生反射和透射，反射波到達地表后由接收天線接收．探地雷達電磁波在層狀結(jié)構(gòu)中的傳播過程滿足Maxwell 方程組，由電磁場理論可知，二維各向同性介質(zhì)中，頻域TE 波方程可以表示為

式中:H 表示磁場;E 表示電場;ε 表示介電常數(shù);σ 表示電導(dǎo)率;σm表示導(dǎo)磁率;μ 表示磁導(dǎo)率;ω表示角頻率;i 表示虛數(shù)單位．

平面問題中電場和磁場的頻域列式可寫為E(x，z)=E(z)·exp(ikx)． (3)H(x，z)=H(z)·exp(ikx)． (4)式中:k 表示波數(shù)．

將式(3)、(4)代入方程(1)、(2)中可得:

方程(5)和(6)可表示為

式中:v=［Ex，Hy］T;v'表示向量v 對z 的偏導(dǎo)數(shù);T 為2 ×2 矩陣，矩陣元素為T11= T22=0，T12=σm+iωμ+k2/(iωε+σ)，T21=σ+iωε．

圖1 層狀結(jié)構(gòu)中探地雷達電磁波傳播示意圖Fig．1 Diagram of GPR wave propagation in layered structure

2 基于邊值問題的精細積分方法

對于線性系統(tǒng)，任意區(qū)間［za，zb］滿足:

Eb=FEa-GHb． (8)

Ha=QEa+PHb． (9)

式中:Ea和Ha表示z =za處的電場和磁場;Eb和Hb表示z=zb處的電場和磁場，為方便描述，省略電場下標(biāo)x 和磁場下標(biāo)y;F，G，P 和Q 為待求的復(fù)數(shù)變量，其值僅與za和zb相關(guān)．

假設(shè)Ea和Ha已知，將方程(8)和(9)對zb求導(dǎo)可得:

E'b=F'Ea-G'Hb-GH'b; (10)

0 =Q'Ea +P'Hb +PH'b． (11)

zb端的控制方程可以表示為

E'b=T11Eb+T12Hb． (12)

H'b=T21Eb+T22Hb． (13)

聯(lián)立方程(8)～(13)可得:

(F' -T11F-GT21F)Ea+(-G' -

T12-GT22+T11G+GT21G)Hb=0． (14)

(Q' +PT21F)Ea+(-PT21G+

P' +PT22)Hb=0． (15)

由于變量Ea和Hb是相互獨立不相關(guān)的任意變量，方程(14)與(15)如果恒成立，那么Ea和Hb的系數(shù)必須為0，由此可得:

2．1 相鄰區(qū)段方程的組裝

對于區(qū)段［za，zb］和［zb，zc］，滿足方程(8)和(9):

Eb=F1Ea-G1Hb; (17)

Ha=Q1Ea+P1Hb; (18)

Ec=F2Eb-G2Hc; (19)

Hb=Q2Eb+P2Hc． (20)

對于合并后的區(qū)段［za，zc］，也應(yīng)滿足方程(8)和(9):

Ec=FcEa-GcHc; (21)

Ha=QcEa+PcHc． (22)

聯(lián)立方程(17)和(20)解得:

將方程(23)和(24)代入方程(18)和(21)整理可得:

對比方程(19)、(20)和方程(25)、(26)左右端項可得區(qū)段變量合并方程組(27):

2．2 區(qū)段變量計算

方程(27)給出了由相鄰區(qū)段［za，zb］和［zb，zc］區(qū)段變量F1，G1，P1，Q1和F2，G2，P2，Q2求解合并后區(qū)段［za，zc］區(qū)段變量Fc，Gc，Pc，Qc的方法，但是如何求解初始區(qū)段變量F1，G1，P1，Q1和F2，G2，P2，Q2仍然未知．

假設(shè)任意一層介質(zhì)i，其厚度di=zi-zi-1，首先將其劃分為2M(M 為正整數(shù))等厚子層，厚度=di/2M，然后再將每個子層分為2N(N一般取20)等厚的微層，微層層厚t=dsubi/2N，由于此時層厚t 非常小，微層的區(qū)段變量F，G，P 和Q 可由泰勒級數(shù)展開求得

式中:qi，gi，ji，yi(i=1，2，3，4)為復(fù)數(shù)變量．

將方程(28)代入方程(16)對比各階系數(shù)可得:

這里應(yīng)注意，由于t 非常小，導(dǎo)致P'(t)和F'(t)非常小，由于計算機存在截斷誤差，直接將P'(t)和F'(t)與1 相加會導(dǎo)致精度損失殆盡，因此需對P'(t)和F'(t)進行單獨儲存計算，方程(30)為:

得到微層區(qū)段變量Ft，Gt，Pt和Qt后，利用方程(30)即可求出合并區(qū)段矩陣Fc，Gc，Pc和Qc，由于同一介質(zhì)層內(nèi)參數(shù)相同，所以每次合并后，層數(shù)會減少一半，經(jīng)過N 次合并后即可得到子層的區(qū)段變量Fsub，Gsub，Psub和Qsub，然后利用子層區(qū)段變量合并求得層i 的區(qū)段變量Fi，Gi，Pi和Qi，這里需要說明的是，此時層厚已經(jīng)不是一個很小的數(shù)，可以采用式(27)進行合并計算． 以此類推求得各層的區(qū)段矩陣，然后再利用式(27)合并得到整個n 層結(jié)構(gòu)的區(qū)段變量Ftotal，Gtotal，Ptotal和Qtotal．

3 激勵源設(shè)置與邊界條件

3．1 激勵源設(shè)置

設(shè)入射波源為第1 層介質(zhì)和上半無限空間分界面處(z=z0)處電場和磁場的間斷，即

式中:Em(k，ω)和Hm(k，ω)分別為電場和磁場的間斷值．z+0表示z 軸正方向一側(cè)z=z0界面，z-0表示z 軸負(fù)方向一側(cè)z=z0界面．

3．2 上部邊界條件

如圖1 所示，n 層結(jié)構(gòu)上部為半無限空間，其狀態(tài)方程可以表示為:

式中:下標(biāo)u 代表上半無限空間介質(zhì);Λu=diag(λu1，λu2)代表矩陣Tu的特征值矩陣;Mu=(αu1，αu2)代表相應(yīng)特征向量矩陣;αui(i =1，2)為對應(yīng)于特征值λui的特征向量．這里需注意特征值的排序，應(yīng)保持第一行特征值λu1為正，對應(yīng)于沿z 軸向上傳播的波，第二行的特征值λu2為負(fù)，對應(yīng)于沿z 軸向下傳播的波．

設(shè)bu(z)=M-1u vu(z)，方程(32)可化為:

由一階常微分方程理論可知，方程(33)的解可以表示為

bu(z)=exp［Λu(z-z0)］·bu(z0)，z ＞z0． (34)

其中bu(z0)是2 ×1 向量，第一行變量表示向上傳播的波，第二行變量表示向下傳播的波．對于上半無限空間來講，不存在向下傳播的波，即bu(z0)第二行為0，因此，

式中:下標(biāo)0 代表向量v 在z =z0處的值，bu1代表bu(z0)中第一行元素． 聯(lián)立求解方程(35)可得上半無限空間輻射邊界條件:

3．3 下部邊界條件

類似的下半無限空間的解可表示為

bd(z)=exp［Λd(zd-z)］·bd(zn)，z ＜zn． (37)式中:下標(biāo)d 表示下半無限空間介質(zhì)，由于下半無限空間中沒有向上傳播的波，即bd(zn)的第一行元素為0，即:

式中:下標(biāo)n 代表向量v 在z =zn處的值，bd2代表bd(zn)中第二行元素．求解方程(38)，可得下半空間輻射邊界條件:

精細積分算法計算流程圖如圖2 所示．

圖2 精細積分算法計算流程圖Fig．2 Flow chat of PIM

4 數(shù)值算例

(1)單層結(jié)構(gòu)反射系數(shù)． 通過一個單層結(jié)構(gòu)來驗證本算法的精度，設(shè)入射波方向與介質(zhì)分界面垂直．即電場與磁場強度與波數(shù)k 不相關(guān)，那么方程(7)系統(tǒng)矩陣可簡化為T11= T22=0，T12=σm+iωμ，T21=σ +iωε． 式(31)中電場和磁場的間斷值可簡化為Em(ω)和Hm(ω)．

設(shè)單層結(jié)構(gòu)介質(zhì)相對介電常數(shù)取9，電導(dǎo)率取0．02 s/m，上部半無限空間介質(zhì)相對介電常數(shù)取1，電導(dǎo)率為0，下半無限空間介質(zhì)相對介電常數(shù)取30，電導(dǎo)率取0．05 s/m． 所有介質(zhì)磁導(dǎo)率都等于真空磁導(dǎo)率． 表1 給出了單層結(jié)構(gòu)上表面處的反射系數(shù)，

由表1 可知，精細積分算法的計算精度與可以達到計算機意義上的解析解．

表1 單層結(jié)構(gòu)反射系數(shù)計算結(jié)果Tab．1 Reflection coefficients of single layered structure

(2)四層半剛性路面結(jié)構(gòu)． 表2 給出了四層半剛性路面結(jié)構(gòu)介電參數(shù)取值(基于系統(tǒng)識別方法反演獲得)． 圖3 為時域入射波信號，由金屬板反射試驗獲得．圖4 和圖5 分別給出了基于精細積分方法和FDTD 方法得到的模擬波形和實測波形的對比圖．

表2 四層半剛性路面結(jié)構(gòu)參數(shù)設(shè)置統(tǒng)計表Tab．2 Parameters of four-layer semi-rigid pavement

圖3 入射波形Fig．3 Incident wave

由圖4 和圖5 可知，兩種方法計算結(jié)果與實測波形基本吻合，但計算時間精細積分方法需要0．417 3 s，F(xiàn)DTD 方法需要0．589 5 s，精細積分方法計算效率比傳統(tǒng)FDTD 方法高40%左右．

圖4 精細積分方法模擬結(jié)果與實測波形對比圖Fig．4 Comparison of the PIM and measured waveform

圖5 時域有限差分方法模擬結(jié)果與實測波形對比圖Fig．5 Comparison of FDTD and measured waveform

5 結(jié)論

筆者基于兩端邊值問題精細積分方法建立層狀結(jié)構(gòu)中探地雷達電磁波傳播數(shù)值模型，該模型可以精確模擬探地雷達電磁波在層狀結(jié)構(gòu)中的傳播過程．通過與路面實測探地雷達回波信號對比可知，模擬波形在波幅、時延等方面與實測波形吻合良好．此外，相比于傳統(tǒng)時域有限差分方法，同等精度下，本算法可節(jié)省約40%的計算時間．

［1］ 鄭宏興，葛德彪． 廣義傳播矩陣法分析分層各向異性材料對電磁波的反射與透射［J］． 物理學(xué)報，2000，49(9):1702 -1706．

［2］ HUANG Yue-qin，ZHANG Jian-zhong，LIU Qinghuo． Three-dimensional GPR ray tracing based on wavefront expansion with irregular cells ［J］． IEEE Transaction on Geoscience and Remote Sensing，2011，49(2):679 -687．

［3］ FAN Guo-xing，LIU Qing-huo，HESTHAVEN J S．Multidomain pseudospectral time-domain simulations of scattering by objects buried in lossy media［J］． IEEE Transaction on Geoscience and Remote Sensing，2002，40(6):1366 -1373．

［4］ 馮德山，王洪華，戴前偉． 基于無單元Galerkin 法探地雷達正演模擬［J］． 地球物理學(xué)報，2013，56(1):298 -308．

［5］ IRVING J，KNIGHT R． Numerical modeling of ground penetrating radar in 2-D using MATLAB ［J］． Computers ＆ Geosicences，2005，32(9):1247 -1258．

［6］ FENG De-shan，DAI Qian-wei． GPR numerical simulation of full wave field based on UPML boundary condition of ADI - FDTD［J］． NDT ＆ E International，2011，44(6):495 -504．

［7］ ZHONG Wan-xie． Combined method for the solution of asymmetric Riccati differential equations ［J］． Computer Methods in Applied Mechanics and Engineering，2001，191(1/2):93 -102．

［8］ ZHONG Wan-xie，LIN Jia-hao，GAO Qiang． The precise computation for wave propagation in a stratified materials［J］． Journal for Numerical Methods in Engineering，2004，60:11 -25．