探地雷達的時域仿真方法研究

王道平， 況曉靜， 范程華

（1.合肥師范學院 教務處，安徽 合肥230061；2.合肥師范學院 電子信息工程學院，安徽 合肥230601）

1 引言

探地雷達（GPR）應用非常廣泛，現已解決了很多工程實際問題，成為淺層勘探地下的高分辨率成像的有力工具。目前國外學者對模擬GPR數值模型研究較多，1995年X.Zeng采用了頻域的方法［1］；1998年H.W.Chen做了二維時域有限差分法的分析［2］；1999年 K.J.Ellefsen采用了積分變換的方法［3］。然而在國內，這一方面的研究工作開展的相對較少。本文采用二維高階時域有限差分結合基于復坐標延伸的遞歸卷積吸收邊界條件（CPML）技術［4］來模擬GPR，給出并通過具體計算實例驗證了所提算法的穩定性條件及算法的有效性。

2 高階時域有限差分近似

GPR技術可以大致分為兩種主要的操作模式：（1）表面反射探測。發射和接收天線分別被放置在表面和地下，通過改變電特性成像。（2）鉆孔式測量。指一個或多個天線被放置在鉆孔中，則地下的特性將被斷層X光攝影裝置測量出來。建立GPR數值模型為探索地下內部結構和GPR數據之間的聯系提供了一種手段。采集到的數據信息來表示地下內部的結構特性。

對于表面反射GPR勘測，采用TM波模式，天線被放置于X－Z測量平面；對于交叉孔和垂直雷達測量，采用TE波模式，天線包含測量平面。在此，我們運用高階時域有限差分方法建立GPR的離散模型。首先引入頻域麥克斯韋的兩個旋度方程

對上述Maxwell方程進行差分離散，傳統的FDTD方法，時間和空間上均采用具有二階精度的中心差分法近似，雖然場分量迭代公式簡單，但數值色散誤差較大。這里采用高階FDTD方法具有良好的數值色散特性，該方法的思想是：在時間上采用二階差分近似，而在空間上則采用四階精度的離散格式。該算法的具體離散過程可如下簡單的表示為時間上：

空間上：

3 數值穩定性分析

在FDTD仿真中，為了確保結果正確，最重要的是選擇合適的時間和空間步長。選擇Δx，Δz，Δt盡可能的大，但Δt過大，則FDTD仿真的數值結果發散。通過穩定性分析，可得高階FDTD方法的穩定性條件

4 CPML吸收邊界

考慮PML邊界條件時，采用基于復坐標延伸空間變量［4］的完全匹配層吸收邊界條件（CPML），（1）式中的▽算子改為以下形式：

其中ck＝κk＋σk／（αk＋iωε0），k＝x，y，z；稱為復坐標延伸空間變量，且僅在k方向上變化［5］。強調說明這里的σk，κk和αk不是真實的電特性參數，它們通過復坐標延伸增加了額外自由度，便于PML邊界的設定。如果選擇延伸參數ck＝1，即σx＝σz＝0，κx＝κz＝1，則TM與TE波則為標準模式。所以在模擬網格的內部設置cx，cz＝1，在PML邊界區，cx，cz設置為復值，即σx和σz大于零，則波會被吸收；設置κx和κz大于1，或者設置αx和αz大于0，會吸收凋零波。在這里采用卷積PML（CPML）方法［4］，利用1／ck在時域中的表達，避免了在其它PML方法中常見的電場和磁場組成部分的分裂，在邊界區域只需要改變很少的坐標延伸變量，而不需要修正FDTD方程等式。選擇CPML方法的優點是與媒質無關，即這個方法可以用在任何媒質中。

在理論上，為了更好的吸收電磁波，PML區域中的σx，σz，κx，κz值應盡可能的大。然而在實際離散的FDTD空間，若電特性改變大，會發生反射，因此PML參數必須被設置成漸變的，從網格內部到邊緣逐漸增加。即

其中κk，σk只在k方向上發生變化，因此這個坐標變量是一維函數。本文計算實例中設置m＝4，κkmax＝5。加入CPML吸收邊界條件后，通過遞歸卷積技術，場分量迭代格式改寫為：

其中參數Νa，Νbk，Νc，Μbk，Μc的定義如下：

參數（i，j）＝（iΔx，jΔz）。這里 Δx 表示x 方向的離散間隔，Δz表示z方向的離散間隔，且（8）式中的ΦHxz，ΦHzx，ΦEyx的定義為：

其中Ak，Bk是CPML的修正系數，且在模擬網格中隨位置而發生變化。

5 實例分析

5.1 反射GPR

下面是TM模式下的電特性模型反射GPR的例子。如圖1地下分層介質，第一層代表滲流沙區ε＝8，δ＝1mS／m，第二層代表材料飽和區ε＝16，δ＝1mS／m第三層代表金屬層ε＝1，δ＝3000mS／m。在z＝0處是代表空氣和地表交界面，設置ε＝1，δ＝0mS／m。源和接收被每隔0.2mm放置在地面上。采用Blackman－Harris源脈沖源，主頻為100MHz。取Δx＝0.04m，Δt＝0.075ns。圖2－圖5顯示了源在10m處取不同時刻的Ey場的傳播圖。如圖2所示，在t＝30ns時，我們捕捉到場從源向外傳播還沒有遇到分層媒質。如圖3所示，在t＝60ns，波到達了第二層，并且開始有部分波反射。如圖4所示，在t＝90ns，波到達了第三層金屬層，波被全部散射回去。如圖5所示，在t＝105ns，波的傳播非常復雜。

圖1 TM模式下的電特性模型

圖2 t＝30ns Ey 場分布

圖3 t＝60ns Ey 場分布

圖4 圖4 t＝90ns Ey 場分布

圖5 t＝105ns Ey 場分布

圖6為時間步長取0.085ns的Ey場分布。在t＝12.75ns時，Ey場分量的數量級達到1030，可見超出最大時間步長時，隨著迭代步數的增加，Ey場分量將會嚴重發散。通過此計算實例驗證了所提算法的穩定性條件的有效性。

圖6 取0.085ns，t＝12.75ns時場分布

5.2 垂直探測GPR

考慮相對介電常數被設置在20～32之間，σ＝5mS／m，μ等于自由空間中的磁導率。源被放置地下1m到10m深，間隔為0.5m，位于x＝0.5m處。接收分別放置在與源相同的深度，位于x＝6m處。主頻為100MHz的Blackman－Harris脈沖作為源函數。取Δx＝0.025m，Δt＝0.02ns。圖7－圖8顯示了在不同時刻的Ez場分量，當Ez源在z＝6m處可以看到波形較圓，這是因為波速在不同的媒質中改變較小。同時可以看到由媒質特性的改變造成了小部分的反射。

圖7 t＝60nsEy場分布

圖8 t＝100nsEy場分布

6 結論

本文提供了一個相對容易理解的二維GPR模型，采用二維高階FDTD數值模擬近似方法來模擬GPR，表現了較高的數值特性。并且在邊界處采用遞歸卷積技術CPML吸收邊界條件，提高了運算效率。同時分析數值穩定性，得到了最大空間和時間離散間隔。最后進行實例仿真，結果與理論分析很好的吻合，證明了這種方法的可行性。

［1］Zeng X，McMechan G A，Cai J，Chen H W.Comparison of ray and Fourier methods for modeling monostatic groundpenetrating radar profiles ［J］. Geophysics.1995，60：1727－1734.

［2］Swick D，Chen H W，Huang T M.Finite difference time domain simulation of GPR data［J］.Journal of Applied Geophysics.1998，40：139－163.

［3］Ellefsen K J.Effects of layered sediments on the guided wave in crosswell radar data［J］.Geophysics.1999，64：1698－1707.

［4］Roden J A，Gedney S D.Convolution PML（CPML）：an efficient FDTD implementation of the CFS－PML for arbitrary media［J］.Microwave and Optical Technology Letters.2000，27：334－339.

［5］Chew W C，Weedon W H.A 3－D perfectly matched medium from modified Maxwell’s equations with stretched coordinates［J］.Microwave and Optical Technology Letters.1994，7：599－604.