劉世奇,李 鈺

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

0 前 言

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

1 探地雷達(dá)計算機仿真的理論基礎(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方向，與x-z平面垂直，這時，探地雷達(dá)工作在TM模式下，涉及的場量為、、三個量。所以麥克斯韋方程可以進一步簡化為：

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

時域有限差分法應(yīng)用的一個重要步驟是進行時間離散參數(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)定，采用以下原則進行時間離散參數(shù)的選擇[2]：

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

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

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

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

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

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

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

探地雷達(dá)的計算機仿真過程其實是一個解一定邊界條件下的麥克斯韋方程的過程。本文設(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)有一定的理想化，但還是對實際地質(zhì)結(jié)構(gòu)的一種抽象，對進行探地雷達(dá)的正演模擬研究很有代表性，對研究探地雷達(dá)的探測機理具有十分重要的作用。圖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上來實現(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，因為這個時候開始，與圖4一樣有反射波的出現(xiàn)。而圖4中反射波開始出現(xiàn)的時刻是在36 ns，這是由于目標(biāo)對象B比目標(biāo)對象A埋得更深，電磁波需要花費更多的時間才能傳播到。從這個圖中還可以看到，反射波出現(xiàn)的位置與圖2模型中設(shè)置的目標(biāo)對象B的位置基本一致，這進一步說明了模擬計算是有效的。圖5（b）為電磁波繼續(xù)傳播4 ns后的場強圖，可以明顯看到波的傳播，電磁波的傳播速度很快。

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

4 結(jié) 語

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