劉世奇,李 鈺

（1.湖南省水利廳，湖南 長沙 410001；2.華東理工大學(xué)信息學(xué)院，上海 200237）

0 前 言

堤壩隱患一直是水利建設(shè)中的大事，探地雷達(dá)是堤壩隱患檢測中一種十分重要的方法。對探地雷達(dá)電磁波的傳播過程進(jìn)行計算機(jī)仿真有利于水利工程師更好地使用探地雷達(dá)進(jìn)行堤壩的隱患排查。目前電磁波探測的技術(shù)原理有很多種，如Goodman和Mcchan等人于1995年提出的基于射線的方法,Zeng等人[1]提出的基于頻域的方法以及Casper等人提出的偽譜方法。但是，這些方法都有一定的局限性，只有在各自特定的條件下才能取得比較好的性能。目前為止，有限差分法（FDTD）在電磁波媒質(zhì)傳播的模擬中應(yīng)用得最為廣泛[2]，效果也比較好，并且適合任意復(fù)雜媒質(zhì)中波的傳播。本文在Matlab軟件中通過建立隱患的理論模型，利用FDTD方法來進(jìn)行仿真，得出了電磁波在到達(dá)隱患表面時的一些傳輸特征，對實(shí)際堤壩檢測中及時發(fā)現(xiàn)隱患具有重要的理論指導(dǎo)意義。

1 探地雷達(dá)計算機(jī)仿真的理論基礎(chǔ)

探地雷達(dá)的數(shù)值模擬可以從兩個麥克斯韋旋度方程開始[3，4]：

由于仿真空間不可能是無限的，所以在本程序中采用了卷積PML吸收邊界條件來抑制邊界區(qū)域波的反射。在這一邊界條件下，計算過程中就可以將有效邊界區(qū)域和吸收邊界區(qū)域分離開來考慮，對程序的編制比較方便。同時，諸多研究也表明基于卷積PML的吸收邊界方法具有很好的性能，吸收邊界層的網(wǎng)格可以比較少，對計算復(fù)雜性也帶來很多好處。為應(yīng)用PML吸收邊界條件[5]，采用了帶伸縮因子的復(fù)坐標(biāo)空間。因此，可以將微分操作符?表示成：

其中：

上式中l(wèi)=x,y,z。

在探地雷達(dá)數(shù)值模擬過程中，考慮波的傳播方向?yàn)?y方向，與x-z平面垂直，這時，探地雷達(dá)工作在TM模式下，涉及的場量為、、三個量。所以麥克斯韋方程可以進(jìn)一步簡化為：

在有效計算區(qū)域內(nèi)部，假設(shè)sx和sz都為1，在吸收邊界網(wǎng)格內(nèi)這兩個參數(shù)可以通過式（4）進(jìn)行計算。為用時域有限差分法來逼近式（5）～（7），采用了經(jīng)典的Yee網(wǎng)格。計算過程中，所有空間導(dǎo)數(shù)利用四階差分進(jìn)行逼近，以達(dá)到盡可能高的精確度。時間導(dǎo)數(shù)采用二階差分進(jìn)行逼近。

時域有限差分法應(yīng)用的一個重要步驟是進(jìn)行時間離散參數(shù)Δt和空間離散參數(shù)的選擇。假如空間離散參數(shù)Δx、Δz和Δt比較大的話，那么網(wǎng)格數(shù)目就會比較少，計算速度就會比較快。但是，參數(shù)Δt過大將導(dǎo)致數(shù)值不穩(wěn)定，參數(shù)Δx和Δz過大將導(dǎo)致空間采樣率不夠，發(fā)生數(shù)值色散。為了使整個正演模擬計算過程數(shù)值穩(wěn)定，采用以下原則進(jìn)行時間離散參數(shù)的選擇[2]：

其中，μmin和εmin分別為計算模型中媒質(zhì)的最小磁導(dǎo)率和最小介電常數(shù)。

對空間離散參數(shù)Δx和Δz的選擇，在程序中按照保證電磁波在媒質(zhì)中傳播一個波長上至少被空間采樣五次這樣的原則來確定。

卷積PML吸收邊界條件[6,7]能使電磁波在經(jīng)過少量吸收邊界網(wǎng)格后就得到衰減，另外在實(shí)現(xiàn)上也比較簡單，只要對邊界區(qū)域的坐標(biāo)變量參數(shù)進(jìn)行設(shè)置就可以了，在實(shí)際的波傳播區(qū)域完全可以不受影響，另外，使計算過程中采用并行運(yùn)算成為可能。在實(shí)際的波傳播區(qū)域，kx和kz等于1，同時σx和σz等于0。而在邊界區(qū)域σx和σz大于0，這樣可使電磁波在進(jìn)入吸收邊界區(qū)域后很快衰減。程序中，波的傳播參數(shù)是按照以下原則設(shè)置的：

上面一行對應(yīng)于波的傳播區(qū)域，下面一行對應(yīng)于波的吸收邊界區(qū)域。

同理，將電導(dǎo)率σl設(shè)置為

從式（9）和（10）還可以看出，電磁波越深入邊界層就會衰減得越快，這也是希望得到的結(jié)果。

2 探地雷達(dá)計算機(jī)仿真流程

探地雷達(dá)的計算機(jī)仿真過程其實(shí)是一個解一定邊界條件下的麥克斯韋方程的過程。本文設(shè)計的整個仿真方法可以以圖1所示的流程來表示。

圖1 探地雷達(dá)模擬計算過程Fig.1 Calculation process of the numerical simulation

3 計算模型與結(jié)果

3.1 計算模型的建立

圖2中三個大小不同、埋于地下不同深度的黃色球狀物是本研究中需要探測的三個目標(biāo)對象。整個模型在縱向基本可以分成三層，Z＞0部分為空氣層，本文對空氣層設(shè)置的介電常數(shù)為ε0=1，磁導(dǎo)率為μ0=1以及電導(dǎo)率σ0=0。與空氣接觸的地下一層為介電常數(shù)為10的一種介質(zhì)，地下第二層為介電常數(shù)為25的媒質(zhì)。圖中三個圓形物質(zhì)從左到右分別為A、B、C，介電常數(shù)均為16，是本文需要探測的三個目標(biāo)。為研究問題的簡化，假定兩種地下媒質(zhì)的磁導(dǎo)率均為1，電導(dǎo)率為最上層的空氣最小，而地下兩層媒質(zhì)與介電常數(shù)對應(yīng)具有不同的電導(dǎo)率。這三層媒質(zhì)用以表示探地雷達(dá)探測過程中遇到的基本地下結(jié)構(gòu)。雖然這種結(jié)構(gòu)有一定的理想化，但還是對實(shí)際地質(zhì)結(jié)構(gòu)的一種抽象，對進(jìn)行探地雷達(dá)的正演模擬研究很有代表性，對研究探地雷達(dá)的探測機(jī)理具有十分重要的作用。圖2描述了FDTD仿真的基本模型，以不同的顏色表示不同的介質(zhì)特征參數(shù)。

圖2 探地雷達(dá)數(shù)值模擬的數(shù)學(xué)模型Fig.2 Numerical model for simulation of GPR

3.2 計算結(jié)果和分析

在探地雷達(dá)的發(fā)射源選擇上，選用了Black-man-Harris脈沖源，在程序中通過將該一定頻率的脈沖源加到電場E→y上來實(shí)現(xiàn)。在整個模擬計算過程中，所采用的電磁波頻率為100 MHz。發(fā)射頻率和波的形式對探測結(jié)果的影響將在下一階段的研究中展開。

從圖3可以看到電磁波的傳播速度與介質(zhì)特征參數(shù)之間的關(guān)系。在空氣中電磁波的傳播速度明顯快于地下一層介質(zhì)。在Z＞0部分的空氣介質(zhì)中，電磁波傳播得很快，而在Z＜0的地下介質(zhì)部分，電磁波的傳播速度比較慢，這是由于地下介質(zhì)的介電常數(shù)εr1大于空氣中的介電常數(shù)ε0，而波的傳播速度與媒質(zhì)的介電常數(shù)的平方根成反比。

圖3 不同媒質(zhì)中波的傳播速度Fig.3 Wave velocity in different medium

圖4 目標(biāo)物A的探測Fig.4 Detection of objective A

圖4描述了當(dāng)探地雷達(dá)的發(fā)射源位于坐標(biāo)（1/96，0，0）位置時對目標(biāo)對象A的探測過程。從圖4（a）可以看到，當(dāng)電磁波經(jīng)過36 ns的傳播后碰到了探測目標(biāo)A，這時侯一部分波被目標(biāo)對象A反射，另一部分波繼續(xù)向前傳播。發(fā)生反射的位置與圖2探測模型中設(shè)置的目標(biāo)對象A的位置基本一致。圖4（b）描述了在圖4（a）的基礎(chǔ)上，電磁波繼續(xù)向前傳播12 ns后的電場Ey圖。從圖中可以看到，反射波經(jīng)過一段時間的傳播后已經(jīng)衰減了很多，但前向波仍然繼續(xù)傳播。這說明反射波相對于透射波比較弱，在媒質(zhì)中衰減得更快。從這三個圖可以看到，正演計算模擬基本能夠比較準(zhǔn)確地探測到目標(biāo)對象A，電磁波的波峰和波谷交替變化向前傳播，當(dāng)媒質(zhì)特征（介電常數(shù)、磁導(dǎo)率、電導(dǎo)率）變化時會發(fā)生反射，反射波的產(chǎn)生與否可以作為判斷是否有探測對象存在的依據(jù)。

圖5描述了當(dāng)探地雷達(dá)的發(fā)射源位于（41/96,0,0）時對目標(biāo)對象B的探測過程。從圖5（a）可以看到，當(dāng)電磁波從波源開始傳播過48 ns后，開始探測到目標(biāo)對象B，因?yàn)檫@個時候開始，與圖4一樣有反射波的出現(xiàn)。而圖4中反射波開始出現(xiàn)的時刻是在36 ns，這是由于目標(biāo)對象B比目標(biāo)對象A埋得更深，電磁波需要花費(fèi)更多的時間才能傳播到。從這個圖中還可以看到，反射波出現(xiàn)的位置與圖2模型中設(shè)置的目標(biāo)對象B的位置基本一致，這進(jìn)一步說明了模擬計算是有效的。圖5（b）為電磁波繼續(xù)傳播4 ns后的場強(qiáng)圖，可以明顯看到波的傳播，電磁波的傳播速度很快。

圖5 目標(biāo)物B的探測Fig.5 Detection of objective B

4 結(jié) 語

本文利用MATLAB對探地雷達(dá)中電磁波的傳播規(guī)律進(jìn)行了計算機(jī)仿真并進(jìn)行了詳細(xì)的分析。從仿真結(jié)果可以看到，建立隱患模型的仿真方法是很有效的，水利工作者利用探地雷達(dá)進(jìn)行大壩隱患探測工作的預(yù)測時可參考探討的方法。

[1]Annan AP.GPR-history，trends and future developments[J].Subsurface Sensing Technologies and Applications,2002,3(4):253-270.

[2]葛德彪.電磁波時域有限差分方法[M].西安：西安電子科技大學(xué)出版社，2002：25-26.

[3]曾昭發(fā)，劉四新.探地雷達(dá)方法原理及應(yīng)用[M].科學(xué)出版社，2006.

[4]李想堂，王端宜，張肖寧.探地雷達(dá)技術(shù)的回顧與展望[J].無損檢測，2006，28（9）:479-483.

[5]Yee K S.Numerical solution of initial boundary value prob-lem involving Maxwell equations isotropicmedia[J].IEEE Trans.Antennas Propagat，May 1996,AP-14(3):302-307.

[6]Engquist B,Majda A.Absorbing boundary conditions for the numerical simulation of waves[J].Math.Comp.,1977,31:629-651.

[7]Daniels D J.A study of surface penetrating radar and its ap-plication[J].Microwave Jr,1994,37(2):68-82.