反射式GPR中地下異構(gòu)介質(zhì)的FDTD二維建模與仿真

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

關(guān)鍵詞:各向異性地質(zhì)模型; 反射式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均勻介質(zhì)的建模仿真，然而近地表介質(zhì)的性質(zhì)經(jīng)常發(fā)生變化，隨著采樣點的增加，這種異質(zhì)性隨之增加[1]。本文采用FDTD方法來模擬GPR地質(zhì)數(shù)據(jù)，沿用了原有的PML空間邊界模型，但是每種邊界內(nèi)的物理模型數(shù)據(jù)由高斯統(tǒng)計分布來決定，這種統(tǒng)計分布比傳統(tǒng)的單值方法更能準確地表示近地表層物理分布特性[2]。本文提出的方法(正演模擬結(jié)果)是利用地下介質(zhì)的分布特點對傳統(tǒng)的均勻空間加入高斯隨機噪聲近似各向異性介質(zhì)那樣影響GPR響應(yīng)，從而創(chuàng)建更加真實的地質(zhì)模型。

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

1 原有的地質(zhì)模型[6]

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

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

圖1 原有地質(zhì)模型

2 改進的地質(zhì)模型

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

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



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



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

圖2 st為0.1時的介質(zhì)模型

3 反射式GPR仿真

模型中輸入信號為Ricker子波，中心頻率為100 MHz。位于空氣和介質(zhì)的水平交界面x=10 m處。本文利用FINDDX.M程序產(chǎn)生的最大空間離散步長為0.043 2 m，選Δx=Δy=0.04 m，最大的時間步長由findtt.m來確定，為0.080 1 ns，采用Δt=0.08 ns，采用PML吸收邊界，用FDTDo(2，4)方案進行數(shù)值模擬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時的介質(zhì)模型

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



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)

上述這些系數(shù)都是空間位置的函數(shù)，與介質(zhì)的ε，μ，δ，kx，kz有關(guān)，所以它們都是隨時間推移而變化的。麥克斯韋方程組仿真模擬采用遞歸卷積技術(shù)轉(zhuǎn)化為\\式(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更新系數(shù)，其隨網(wǎng)格節(jié)點位置變化而變化。從方程(6)可以看出，在當前時刻的卷積因子的值是由前一時刻值計算得出的。因此，在FDTD仿真期間，必須將ΨHxA，ΨHzx，ΨEyx和ΨEyzA以及各場分量的前一時刻值存儲起來以便計算。

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

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

圖4 GPR的FDTD仿真流程圖

4 仿真圖像

由于篇幅限制本文僅對簡單異構(gòu)介質(zhì)模型當st=0.1時進行反射式GPR的FDTD仿真，模型中輸入信號源為Ricker子波，中心頻率為100 MHz。在空氣和介質(zhì)交界面沿x方向每0.2鉆孔測量。利用finddx.m程序產(chǎn)生的最大空間離散步長為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邊界吸收條件，從圖中可以看到任意平面內(nèi)都沒有發(fā)現(xiàn)源自模型邊界區(qū)域的反射，可見吸收效果非常好。

5 結(jié) 論

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

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

6 結(jié) 語

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

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

參考文獻

[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混合編程實現(xiàn)圖像三維重建\\. 微計算機信息， 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]曾昭發(fā). 探地雷達方法原理及應(yīng)用\\. 北京: 科學(xué)出版社， 2006.

[6]周輝. 地質(zhì)雷達數(shù)據(jù)處理現(xiàn)狀和展望\\. 地學(xué)前緣， 2001， 8(2): 230-231.

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

[8]葛德彪， 閆玉波. 電磁波時域有限差分方法\\. 西安: 西安電子科技大學(xué)出版社， 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]譚廷棟. 測井學(xué)\\. 北京: 石油工業(yè)出版社， 1998.