反射式GPR中地下異構介質的FDTD二維建模與仿真

摘 要:在此利用地下表層介質物理特性具有高斯統計分布這一分布特性在傳統的均勻空間加入高斯隨機信號，在進行探測時這些背景信號就會影響GPR探測響應，從而得到相應的地質圖像，建立了各向異性地質模型，理論上通過這種建模方法可以創建任何復雜的地質模型，這將為研究地下介質提供更豐富的地質資料。 在最后用反射式GPR對所建地質模型進行了 FDTD的二維Matlab仿真，取得了預想的效果，仿真圖像形象地證明了這一各向異性地質模型的可行性。

關鍵詞:各向異性地質模型; 反射式GPR; Matlab仿真; FDTD

中圖分類號:TP3919文獻標識碼:A

文章編號:1004-373X(2010)08-0151-03

FDTD 2-D Modeling and Simulation of Undergound Heterogeneity Media in Reflective GPR

ZHOU Mei-li， BAI Zong-wen， LI Hai-jun

(College of Physics and Electronic Information， Yan’an University， Yan’an 716000， China)

Abstract: The physical properties in the near surface of underground has Gauss statistical distribution. Therefore， when probing the underground with GPR technique， the Gauss random variables may add into the homogeneous space by which the heterogeneity space is created. Then those background signals may influence the GPR detection response， by which the responding geologic images are achieved. Theoretically， a complex geologic model can be established with the method， which can enrich the geologic explanation data. The geologic model using the reflection-model GPR technical， and the simulating images are provided by FDTD simulation.

Keywords: heterogeneity geologic model; reflective GPR; Matlab simulation; FDTD

0 引 言

以往的論文研究主要集中在GPR均勻介質的建模仿真，然而近地表介質的性質經常發生變化，隨著采樣點的增加，這種異質性隨之增加[1]。本文采用FDTD方法來模擬GPR地質數據，沿用了原有的PML空間邊界模型，但是每種邊界內的物理模型數據由高斯統計分布來決定，這種統計分布比傳統的單值方法更能準確地表示近地表層物理分布特性[2]。本文提出的方法(正演模擬結果)是利用地下介質的分布特點對傳統的均勻空間加入高斯隨機噪聲近似各向異性介質那樣影響GPR響應，從而創建更加真實的地質模型。

文中改進后的各向異性介質模型與原有模型相比，能更真實地反映地下介質情況。為了描述簡單，本文只討論了二維空間的地質模型仿真圖像，但這一方法完全可以應用到三維空間中，也可根據計算機視覺的研究方法實現二位圖像的三維重建[3]，從而得到其三維情況的信息。在文中引用了原有的模型和改進模型的仿真圖像，根據仿真圖像的比較就可判斷相應的地質分布情況不再均勻。那么這樣的模型如果應用到GPR探測仿真，就可獲得更多的地質解釋資料[4-5]，理論上可以解釋任何一種地質分布情況。這將對探測地下介質分布有很大的幫助。

1 原有的地質模型[6]

原有模型是一個用于TM模式的反射GPR地質模型，如圖1所示。地下介質分為兩層，寬度20 m，深度10 m，其中上層表示泥沙滲流帶，其εr=9，σ=1 ms/m，下層表示飽和帶的介質εr=25，σ=5 ms/m，上層中有三個大小互不相同的異質區域，其中εr=16，σ=1 ms/m。對于所有介質，μ的值設定為大小等于自由空間的值μ0，空氣和地面的交界面包含在模型z=0的平面中，通過的向網格內部增加薄層來簡單地實現，令其εr=1，σ=0 ms/m。沿空氣和地面的交界面處可放置反射式探地雷達系統裝置，即可進行GPR仿真。

這一模型與傳統的均勻半空間模型[7]比較已考慮到地下介質中的快狀障礙物或硬塊的存在，可是實際情況往往不僅如此，近地表介質物理特性由于介質含水率差異，植被均異、材料類型差異等因素形成介質各向異性[8]，即在同一點處沿各個方向上其電特性參數不同。鑒于這種情況有必要對地下介質模型做徹底性的改進，使其更加接近實際。

圖1 原有地質模型

2 改進的地質模型

如上所述在局部范圍內由于濕度、植被等導致介質各向異性[9-10]，對于介質的這種電特性可用高斯隨機分布來描述。即在垂直z方向上加上高斯隨機信號，使得沿z方向介質各向異性。

在加高斯隨機信號時應該遵從的原則就是標準偏差應小于原有介質模型中一個元胞的長度以免發生空間色散，也就是對這種介質的改進從標準偏差0.1~1之間進行選擇來確保原有介質參數的主導性，這也符合實際情況，加入高斯信號后各介質的電性參數為:



εrm=εr+st×randon(1)



式中:εrm是改進后介質的相對介電常數;εr是原來介質的相對介電常數;st為標準偏差，random為高斯隨機變量。改進之后的介質模型Matlab仿真下列圖形。圖2，圖3分別是當st值為0.1，0.75時的各向異性介質的仿真圖。

圖2 st為0.1時的介質模型

3 反射式GPR仿真

模型中輸入信號為Ricker子波，中心頻率為100 MHz。位于空氣和介質的水平交界面x=10 m處。本文利用FINDDX.M程序產生的最大空間離散步長為0.043 2 m，選Δx=Δy=0.04 m，最大的時間步長由findtt.m來確定，為0.080 1 ns，采用Δt=0.08 ns，采用PML吸收邊界，用FDTDo(2，4)方案進行數值模擬Ricker子波，即所有的場分量的空間展開都保留四階精度的有限差分近似表達式，時間展開保留二階精度有限差分近似。模擬原理如下:



Ey|n+1i，j=Ca|i，j\\+Cbx|i，j\\+

Cbz|i，j\\+Ccx|i，j\\(2)

Hx|n+1/2i，j+1/2=Hx|n+1/2i，j+1/2-Dbz|i，j+1/2\\-Dbz|i，j+1/2\\(3)

Hz|n+1/2i+1/2，j=Hz|n-1/2i+1/2，j-Dbz|i+1/2，j\\-Dcz|i+1/2，j\\(4)



圖3 st為0.75時的介質模型

式中:下標表示空間位置，上標表示時間，采用FDTD模型模擬這些方程就是通過電場和磁場在時間順序上交替抽樣向前推移。式(2)~式(4)中的更新系數Ca，Cbx，Cbz，Cc，Dbx，Dbz，Dc分別表示電特性和各網格參數，其形式給定如下:



Ca=\\(1+σΔt/2ε)-1

Cbk=(Δt/ε)\\(24kk+Δt)-1

Cc=(Δt/ε)(1+σΔt/2ε)-1

Dbk=(Δt/μ)(24kkΔk)-1

Dc=Δt/μ(5)

上述這些系數都是空間位置的函數，與介質的ε，μ，δ，kx，kz有關，所以它們都是隨時間推移而變化的。麥克斯韋方程組仿真模擬采用遞歸卷積技術轉化為\\式(2)~式(4)，式中:



ΨHzx|ni+1/2，j=Bx|i+1/2，j\\+Ax|i+1/2\\ΨEyz|n+1/2i，j

=Bz|i，j-ΨEyz|n-1/2i，j+Ax|i，j

ΨHxxz|ni，j+1/2=Bz|i，j+1/2\\+Az|i，j+1/2

ΨEyx|n+1/2i，j=Bx|i，j-ΨEyxz|n-1/2i，j+Ax|i，j\\(6)

式中:



Ak=σkσkkk+αkk2k(Bk-1)，Bk=eΔtε0(σk/kk+αk)(7)



是PML更新系數，其隨網格節點位置變化而變化。從方程(6)可以看出，在當前時刻的卷積因子的值是由前一時刻值計算得出的。因此，在FDTD仿真期間，必須將ΨHxA，ΨHzx，ΨEyx和ΨEyzA以及各場分量的前一時刻值存儲起來以便計算。

在對式(2)~式(4)進行FDTD計算機模擬仿真期間，為了引進電場源，在空間網格節點位置加入脈沖源來更新Ey分量。

綜上是FDTD的二維建模與仿真原理，以下將對二維反射式GPR模型探測進行模擬仿真，軟件仿真流程如圖4所示。

圖4 GPR的FDTD仿真流程圖

4 仿真圖像

由于篇幅限制本文僅對簡單異構介質模型當st=0.1時進行反射式GPR的FDTD仿真，模型中輸入信號源為Ricker子波，中心頻率為100 MHz。在空氣和介質交界面沿x方向每0.2鉆孔測量。利用finddx.m程序產生的最大空間離散步長為0.043 2 m，選Δx=Δy=0.04 m，最大的時間步長由findtt.m來確定，為0.080 1 ns，采用Δt=0.08 ns，得仿真圖像如圖5所示。

圖5表示輸入信號源位于x=10 m處進行的FDTD仿真效果。圖中描述的是Ey場分量截面圖。這里采用的是PML邊界吸收條件，從圖中可以看到任意平面內都沒有發現源自模型邊界區域的反射，可見吸收效果非常好。

5 結 論

由上述仿真圖像分析可得:在t=30 ns時，沒有反射波，表明傳播波還沒有遇到地下異構介質;在t=50 ns時，很明顯波在上層介質中部遇到了較大的異常塊，部分波被反射回來。由于直達波在空氣中傳播較介質中快，所以此時源自空氣的散射波也在這些點留下模型域的邊界。在t=70 ns時，反射波已回到空氣和地面的交界面。在t=110 ns時，發現源自兩層介質交界面處的反射。圖5所描述波的傳播情況正好驗證了上述異構介質模型結構，可見利用這種方法理論上可以創建任意的地質資料以供參考研究。由于篇幅的限制，在這里僅對圖2這一地質模型進行仿真，即使對于如同圖2這樣一個簡單的地球模型，波的傳播迅速變得非常復雜，可見在仿真其他參雜更大噪聲信號的模型時地質圖像會更復雜，下面的工作就可以從這一方向考慮。

圖5 不同時刻Ey分量傳播情況的記錄

6 結 語

利用地下表層介質物理特性具有高斯統計分布這一分布特性在傳統的均勻空間加入高斯隨機信號，在進行探測時這些背景信號就會影響GPR探測響應，從而得到相應的地質圖像， 建立了各向異性地質模型，理

論上通過這種建模方法可以創建任何復雜的地質模型，在最后用反射式GPR對所建地質模型進行了 FDTD的二維Matlab仿真，取得了預想中的效果，仿真圖像形象地證明了這一各向異性地質模型的可行性。

參考文獻

[1]JENNIFER Jane Holt. Finite difference time domain modeling of dispersion from heterogeneous ground properties in ground penetrating radar\\. USA: Ohio State University， 2004.

[2]DANIELS D J. Surface penetrating radar\\. London: IEE Press， 1996.

[3]孫惠， 孫莉. Visual C++與Matlab混合編程實現圖像三維重建\\. 微計算機信息， 2007， 23(3): 294296.

[4]PETERS J R， DANIELS D J， YOUNG J D. Ground penetrating radar as a subsurface environmental sensing tool\\. Proc. of IEEE， 1994(12): 1802-1821.

[5]曾昭發. 探地雷達方法原理及應用\\. 北京: 科學出版社， 2006.

[6]周輝. 地質雷達數據處理現狀和展望\\. 地學前緣， 2001， 8(2): 230-231.

[7]KNIGHT J I R. Numerical modeling of ground-penetrating radar in 2-D using Matlab\\. ComputersGeosciences， 2005.

[8]葛德彪， 閆玉波. 電磁波時域有限差分方法\\. 西安: 西安電子科技大學出版社， 2002.

[9]JENNIFER Jane Holt. Finite difference time domain modeling of dispersion from heterogeneous ground properties in ground penetrating radar\\. USA: Ohio State University，2004:19-20.

[10]LUNT I A， HUBBARD S S， Soil moisture content estion using gound-penetrating radar reflection data\\. Journal of Hydrology， 2005， 307: 254-269.

[11]譚廷棟. 測井學\\. 北京: 石油工業出版社， 1998.