王道平， 況曉靜， 范程華

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

1 引言

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

2 高階時域有限差分近似

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

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

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

空間上：

3 數(shù)值穩(wěn)定性分析

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

4 CPML吸收邊界

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

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

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

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

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

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

其中Ak，Bk是CPML的修正系數(shù)，且在模擬網(wǎng)格中隨位置而發(fā)生變化。

5 實例分析

5.1 反射GPR

下面是TM模式下的電特性模型反射GPR的例子。如圖1地下分層介質(zhì)，第一層代表滲流沙區(qū)ε＝8，δ＝1mS／m，第二層代表材料飽和區(qū)ε＝16，δ＝1mS／m第三層代表金屬層ε＝1，δ＝3000mS／m。在z＝0處是代表空氣和地表交界面，設(shè)置ε＝1，δ＝0mS／m。源和接收被每隔0.2mm放置在地面上。采用Blackman－Harris源脈沖源，主頻為100MHz。取Δx＝0.04m，Δt＝0.075ns。圖2－圖5顯示了源在10m處取不同時刻的Ey場的傳播圖。如圖2所示，在t＝30ns時，我們捕捉到場從源向外傳播還沒有遇到分層媒質(zhì)。如圖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場分量的數(shù)量級達到1030，可見超出最大時間步長時，隨著迭代步數(shù)的增加，Ey場分量將會嚴重發(fā)散。通過此計算實例驗證了所提算法的穩(wěn)定性條件的有效性。

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

5.2 垂直探測GPR

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

圖7 t＝60nsEy場分布

圖8 t＝100nsEy場分布

6 結(jié)論

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

［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.