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

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

0 引言

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

1 控制方程

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

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

平面問(wèn)題中電場(chǎng)和磁場(chǎng)的頻域列式可寫(xiě)為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 對(duì)z 的偏導(dǎo)數(shù);T 為2 ×2 矩陣，矩陣元素為T(mén)11= T22=0，T12=σm+iωμ+k2/(iωε+σ)，T21=σ+iωε．

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

2 基于邊值問(wèn)題的精細(xì)積分方法

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

Eb=FEa-GHb． (8)

Ha=QEa+PHb． (9)

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

假設(shè)Ea和Ha已知，將方程(8)和(9)對(duì)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是相互獨(dú)立不相關(guān)的任意變量，方程(14)與(15)如果恒成立，那么Ea和Hb的系數(shù)必須為0，由此可得:

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

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

Eb=F1Ea-G1Hb; (17)

Ha=Q1Ea+P1Hb; (18)

Ec=F2Eb-G2Hc; (19)

Hb=Q2Eb+P2Hc． (20)

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

Ec=FcEa-GcHc; (21)

Ha=QcEa+PcHc． (22)

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

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

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

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

方程(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，然后再將每個(gè)子層分為2N(N一般取20)等厚的微層，微層層厚t=dsubi/2N，由于此時(shí)層厚t 非常小，微層的區(qū)段變量F，G，P 和Q 可由泰勒級(jí)數(shù)展開(kāi)求得

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

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

這里應(yīng)注意，由于t 非常小，導(dǎo)致P'(t)和F'(t)非常小，由于計(jì)算機(jī)存在截?cái)嗾`差，直接將P'(t)和F'(t)與1 相加會(huì)導(dǎo)致精度損失殆盡，因此需對(duì)P'(t)和F'(t)進(jìn)行單獨(dú)儲(chǔ)存計(jì)算，方程(30)為:

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

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

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

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

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

3．2 上部邊界條件

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

式中:下標(biāo)u 代表上半無(wú)限空間介質(zhì);Λu=diag(λu1，λu2)代表矩陣Tu的特征值矩陣;Mu=(αu1，αu2)代表相應(yīng)特征向量矩陣;αui(i =1，2)為對(duì)應(yīng)于特征值λui的特征向量．這里需注意特征值的排序，應(yīng)保持第一行特征值λu1為正，對(duì)應(yīng)于沿z 軸向上傳播的波，第二行的特征值λu2為負(fù)，對(duì)應(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 向量，第一行變量表示向上傳播的波，第二行變量表示向下傳播的波．對(duì)于上半無(wú)限空間來(lái)講，不存在向下傳播的波，即bu(z0)第二行為0，因此，

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

3．3 下部邊界條件

類(lèi)似的下半無(wú)限空間的解可表示為

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

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

精細(xì)積分算法計(jì)算流程圖如圖2 所示．

圖2 精細(xì)積分算法計(jì)算流程圖Fig．2 Flow chat of PIM

4 數(shù)值算例

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

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

由表1 可知，精細(xì)積分算法的計(jì)算精度與可以達(dá)到計(jì)算機(jī)意義上的解析解．

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

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

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

圖3 入射波形Fig．3 Incident wave

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

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

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

5 結(jié)論

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

［1］ 鄭宏興，葛德彪． 廣義傳播矩陣法分析分層各向異性材料對(duì)電磁波的反射與透射［J］． 物理學(xué)報(bào)，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］ 馮德山，王洪華，戴前偉． 基于無(wú)單元Galerkin 法探地雷達(dá)正演模擬［J］． 地球物理學(xué)報(bào)，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．