基于GPR數(shù)值模擬正演的地下介質(zhì)分析方法

宋正淳，邱哲理，張雪冰，孫 超，胡 昕，覃育祥

吉林建筑大學(xué) 測(cè)繪與勘查工程學(xué)院,長(zhǎng)春 130118

0 引言

探地雷達(dá)[1]探測(cè)地下介質(zhì)多為0.5 m～5 m區(qū)間內(nèi),具有無(wú)損探測(cè)、探測(cè)速度快且探測(cè)過(guò)程連續(xù)、分辨率高、操作方便靈活、探測(cè)費(fèi)用低等優(yōu)點(diǎn).因其優(yōu)點(diǎn)故目前廣泛用于考古、礦產(chǎn)勘查、巖土工程勘探、工程質(zhì)檢、建筑結(jié)構(gòu)檢測(cè)以及軍事勘查等眾多領(lǐng)域中.探地雷達(dá)其工作原理是利用天線激發(fā)高頻脈沖電磁波傳至地下有耗介質(zhì)中,地下各層有耗介質(zhì)各不相同且介質(zhì)具有不均勻性、各向異性、強(qiáng)衰減性，故接收的電磁波信號(hào)波動(dòng)也不同.利用其電磁信號(hào)的差異性,可將一道道回波信號(hào)接收的振幅與時(shí)間、頻率等參數(shù)來(lái)分析和推斷被探測(cè)土體的地下構(gòu)造并數(shù)據(jù)處理繪制成圖像.而本文研究的gprMax[2]軟件則是將電磁波信號(hào)模擬仿真,仿真信號(hào)波動(dòng)數(shù)據(jù)處理組成B-scan圖[3],屆時(shí)可通過(guò)圖形研究地下不同介質(zhì)的位置及形態(tài).其方法是基于時(shí)域有限差分(Finite difference time domain,縮寫為FDTD)[4]的麥克斯韋[5]方程求解電磁信號(hào)強(qiáng)度.FDTD方法是將仿真區(qū)域拆分成一個(gè)個(gè)空間網(wǎng)格,用數(shù)值有限差分方法逼近麥克斯韋方程,并將時(shí)間離散化,再根據(jù)已知的初始邊界條件和介質(zhì),通過(guò)數(shù)據(jù)計(jì)算處理,將一段時(shí)間內(nèi)的反射信號(hào)數(shù)據(jù)進(jìn)行整合,即可將仿真模型的二維或三維的數(shù)據(jù)可視化.介電常數(shù)存在差異的兩個(gè)界面會(huì)反射電磁波,觀察生成的B-scan圖的電磁波現(xiàn)象,通過(guò)兩組實(shí)驗(yàn)對(duì)比觀察其極化[6]方向和相位[7]變化的規(guī)律.

1 探地雷達(dá)數(shù)據(jù)模擬

1.1 探地雷達(dá)數(shù)據(jù)模擬基本原理

從宏觀上來(lái)說(shuō),所有的電磁現(xiàn)象均可由麥克斯韋方程解釋.以下均是一階偏微分方程,用來(lái)表示基本電磁場(chǎng)量與其源的關(guān)系:

×E=-?B?t

(1)

(2)

(3)

(4)

式中，E是電場(chǎng)強(qiáng)度,V/m;H是磁場(chǎng)強(qiáng)度,V/m;B是磁通密度,Wb/m2;D是電位移,C/m;J是電流密度,A/m2;q是電荷密度,C/m3;t是時(shí)間,s.除上述式子外,其媒質(zhì)本構(gòu)關(guān)系式是:

D=εE，B=μH，J=σE

(5)

式中，ε是相對(duì)介電常數(shù);μ是磁導(dǎo)率,H/m;σ是電導(dǎo)率,S/m.基于麥克斯韋方程組下,根據(jù)模型的位置及時(shí)間的條件下求解上述方程式,得到探地雷達(dá)模擬數(shù)據(jù).將模擬的電磁信號(hào)數(shù)據(jù)分析處理即可得到A-scan,B-scan,C-scan圖像,B-scan圖像即是本文主要分析研究的重點(diǎn),B-scan圖像可將模擬地質(zhì)的地下模型剖切面用二維表示出來(lái).

1.2 電磁波在界面的反射系數(shù)與相位變化

電磁波被探地雷達(dá)的天線激發(fā)信號(hào)射入介質(zhì)界面上會(huì)發(fā)生反射和透射現(xiàn)象,其原理類似于光入射到水面和玻璃面.在GPR中利用反射器的傳統(tǒng)技術(shù)可以將速度走時(shí)分析與CMP數(shù)據(jù)采集方法相結(jié)合,但確定雷達(dá)波的速度則需用反射波的傳播時(shí)間來(lái)確定.其方法存在缺陷,電磁波進(jìn)入介面后的回波信號(hào)深度有局限性,不能確定最低反射面以下介質(zhì)的速度.電磁波從天線中釋放出其電磁信號(hào)可能會(huì)產(chǎn)生極化,極化類型由天線的設(shè)計(jì)確定.天線的設(shè)計(jì)可以是線性的,圓形的或橢圓形的.絕大多數(shù)的探地雷達(dá)發(fā)出的波是線性的極化,線性的極化又可分為垂直極化和平行極化.極化類型決定電磁波如何反射,電磁波在界面的反射系數(shù)可定義為反射波的振幅(或能量)Er與入射波的振幅Ei之比.需要注意的是,垂直極化和平行極化的反射系數(shù)有所不同,實(shí)際探地雷達(dá)數(shù)據(jù)采集多以垂直極化為主,垂直極化的反射系數(shù)可表示為R⊥[8]:

R⊥=ErEi=ε1cosθ1-ε2-ε1sin2θ1ε1cosθ1+ε2-ε1sin2θ1

(6)

式中,ε1,ε2分別對(duì)應(yīng)第一和第二層介質(zhì)的相對(duì)介電常數(shù);θ1是入射角度.

電磁波入射到不同介質(zhì)時(shí),一部分會(huì)發(fā)生反射現(xiàn)象,另一部分電磁波則會(huì)發(fā)生透射現(xiàn)象透過(guò)分界面.電磁波在界面的透射系數(shù)可定義為Et是透射波的振幅(或能量)與Ei是入射波的振幅之比.

T⊥=EtEi=2ε1cosθ1ε1cosθ1+ε2-ε1sin2θ1

(7)

式中,T⊥是垂直極化的透射系數(shù).

理想完美電導(dǎo)體(Perfect Electric Conductor,PEC)不發(fā)生透射現(xiàn)象,是因?yàn)镻EC導(dǎo)體的相對(duì)介電常數(shù)無(wú)窮大,當(dāng)?shù)诙咏橘|(zhì)為PEC時(shí),上式中分母無(wú)窮大,分子是除PEC材料以外的其他非導(dǎo)體任意材料,透射系數(shù)T⊥趨近于零.

反射系數(shù)和透射系數(shù)之間滿足公式:

1+R⊥=T⊥

(8)

電磁波的頻率由發(fā)射的波源決定,探地雷達(dá)中多選擇Richer型波源,頻率不變.電磁波從天線釋放出信號(hào)波,信號(hào)波在同種介質(zhì)中傳播波長(zhǎng)不發(fā)生變化.當(dāng)電磁波從介質(zhì)進(jìn)入其他介質(zhì)中時(shí)波長(zhǎng)會(huì)減小,振幅也會(huì)因阻尼效應(yīng)減小;當(dāng)電磁波靠近PEC時(shí),電磁波不會(huì)進(jìn)入PEC中而是會(huì)將電磁波反射.共享相同頻率的多個(gè)采樣道電磁波相互作用產(chǎn)生相位變化,相位變化是由兩個(gè)或多個(gè)采樣道電磁波波形的振幅波峰和波谷位置關(guān)系確定,通過(guò)距離、時(shí)間、度數(shù)測(cè)量相位.介面中某一點(diǎn)的波阻抗公式為:

η=EH=(με)0.5

(9)

式中,η是電場(chǎng)與磁場(chǎng)的復(fù)振幅之比定義為波阻抗,Ω.

2 模型建模與結(jié)果分析

因本文主要討論的是探地雷達(dá)的高頻脈沖電磁波入射到地下的電磁波相位及極化方向,而設(shè)置兩層介質(zhì)界面后第二層界面的入射波是來(lái)自于第一層界面的透射波,其信號(hào)波動(dòng)不太明顯,故擬構(gòu)建兩組濕沙土樣為背景介質(zhì)的模型(如圖1,圖2所示),并分別設(shè)置2組不同尺度、位置和深度的金屬和空洞異常體.

圖1 小尺度圓柱異常體正演模型,圓柱直徑均為75 mmFig.1 The forward model of small-scale cylindrical anomalies at the cylinder diameter of 75 mm

圖2 大尺度圓柱異常體正演模型,圓柱直徑均為300 mmFig.2 The forward model of large-scale cylindrical anomalies at the cylinder diameter of 300 mm

2組模型均由表層0.03 m的空氣和0.6 m深的濕沙構(gòu)成;異常體均為圓柱形異常體,左側(cè)2個(gè)異常體采用PEC金屬電導(dǎo)體,右側(cè)2個(gè)異常體均為空洞.模型一中異常體直徑均為75 mm,深度依次為400 mm,350 mm,400 mm,350 mm;模型二中異常體直徑均為300 mm,深度依次為350 mm,250 mm,350 mm,250 mm.在gprMax正演模擬中,選擇激發(fā)源均為主頻為900 MHz的Ricker子波,激發(fā)源與接受源間距均為50 mm.

不同深度與直徑的兩組管線模型的gprMax正演結(jié)果如圖3,圖4所示.

圖3 小尺度圓柱異常體正演剖面,圓柱直徑均為75 mmFig.3 The forward section of small-scale cylindrical anomaly at the diameter of 75 mm

圖4 大尺寸圓柱異常體正演剖面,圓柱直徑均為300 mmFig.4 The forward section of large-scale cylindrical anomaly at the diameter of 300 mm

電磁波傳播到存在介電常數(shù)差異的界面時(shí)會(huì)在相應(yīng)位置會(huì)呈現(xiàn)波形,紅色與藍(lán)色分別代表正向和負(fù)向振幅,圖中顏色越深代表其電場(chǎng)強(qiáng)度越強(qiáng).正演模型的關(guān)鍵過(guò)程如下(由于兩組原理類似故僅對(duì)圖4過(guò)程進(jìn)行描述),其中左側(cè)PEC異常體上方的電磁波傳播大致可描述為2個(gè)過(guò)程.

過(guò)程1 以圖4為例.探地雷達(dá)從天線發(fā)射的高頻電磁脈沖(波形約可用Ricker子波表示),其波形為一組由正向-負(fù)向-正向-負(fù)向振幅組成的子波(圖中由紅-深藍(lán)-深紅-淺藍(lán)表示,見(jiàn)2 ns處直達(dá)波),通過(guò)較薄的空氣層傳入濕沙介質(zhì)中.電磁波在濕沙介質(zhì)中傳播至PEC金屬異常體表面,在剖面約時(shí)間深度6 ns處可見(jiàn)明顯的雙曲線(參見(jiàn)紅色箭頭),其顏色依次是淺藍(lán)-深紅-深藍(lán)-淺紅表示.通過(guò)對(duì)比直達(dá)波可知極性發(fā)生了反轉(zhuǎn).

過(guò)程2 電磁波經(jīng)PEC異常體時(shí)反射到濕沙界面和空氣界面之間,一部分電磁波透射到空氣介層,另一部分又反射到PEC異常體表面處(多次波).因?yàn)橐徊糠蛛姶挪ㄍ干涞搅丝諝饨橘|(zhì)中,圖4中黑色箭頭處電磁波產(chǎn)生衰減,但其振幅依次為微紅-淺藍(lán)-淺紅-微藍(lán),隨后再次傳播的多次波因振幅衰減在圖中已無(wú)法識(shí)別.圖2中左二的藍(lán)色圓柱完美電導(dǎo)體由于位置距離空氣界面層較遠(yuǎn),所以圖中可識(shí)別的電磁波只有一次.上述兩過(guò)程中的反射系數(shù)R⊥,透射系數(shù)T⊥以及反射與透射電場(chǎng)強(qiáng)度之比Er∶Et參見(jiàn)表1.

對(duì)于右側(cè)空洞異常體,電磁波傳播可大致描述為3個(gè)過(guò)程：

過(guò)程1 電磁波由空氣界面透過(guò)到濕沙介面之后運(yùn)動(dòng)到左測(cè)空心圓柱異常體頂部,有一部分透射到了空洞空氣介質(zhì)管內(nèi),另一部分則反射到了濕沙層頂部,回波被接收天線接收,見(jiàn)紅色箭頭處雙曲線同相軸,按接收時(shí)間順序其電場(chǎng)強(qiáng)度由淺紅-藍(lán)-紅-淺藍(lán)表示.

過(guò)程2 之前一部分電磁波透射到空洞內(nèi)部后,電磁波運(yùn)動(dòng)到空洞底部后一部分電磁波經(jīng)反射,回波被接收天線接收,見(jiàn)藍(lán)色箭頭位置,其顏色變化依次是微紅-淺藍(lán)-淺紅-微藍(lán).

過(guò)程3 當(dāng)電磁波傳播到濕沙與表層空氣界面時(shí),一部分經(jīng)反射又傳播到空洞頂部,部分電磁波經(jīng)反射后又回到地表被接收天線接收,在圖4黑色箭頭處產(chǎn)生電磁波,電場(chǎng)強(qiáng)度顏色變化依次是微紅-微藍(lán)-微紅-微藍(lán).左二與左四情況相同均因埋深較深且沒(méi)有進(jìn)行衰減補(bǔ)償處理,所以后續(xù)多次波不明顯.

通過(guò)圖片對(duì)比發(fā)現(xiàn)電磁波在濕沙界與PEC界面顏色不發(fā)生改變,正負(fù)極振幅不發(fā)生改變.電磁波由濕沙界面到圓柱空心管或由空心管到濕沙界面時(shí),顏色發(fā)生改變,正負(fù)極振幅發(fā)生相位變化.其原因是完美電導(dǎo)體沒(méi)有介電常數(shù),不發(fā)生透射,只有反射,極化方向不發(fā)生改變.空心管則是由高介電常數(shù)向低介電常數(shù)透射電磁波,極化方向發(fā)生變化,磁場(chǎng)強(qiáng)度會(huì)由衰減.反射系數(shù)與透射系數(shù)的模型正演關(guān)系見(jiàn)表1,表2.

表1 金屬異常體反射系數(shù)與透射系數(shù)的模型正演Table 1 Forward model of reflection coefficient and transmission coefficient of metal abnormal bodies

表2 空心異常體反射系數(shù)與透射系數(shù)的模型正演Table 2 Forward model of reflection coefficient and transmission coefficient of hollow abnormal bodies

3 結(jié)語(yǔ)

本文結(jié)合gprMax進(jìn)行探地雷達(dá)正演模擬,探討電磁波在地下介質(zhì)中的運(yùn)動(dòng)路徑和相位變化規(guī)律,通過(guò)gprMax對(duì)不同材質(zhì)、不同深度、不同大小的地下異常體進(jìn)行正演模擬并對(duì)比,就此算出關(guān)鍵反射界面的反射系數(shù)和投射系數(shù).通過(guò)已有實(shí)驗(yàn)數(shù)據(jù)的介電常數(shù)進(jìn)行仿真模擬地下金屬異常和空洞,對(duì)比驗(yàn)證了其結(jié)果與經(jīng)典電磁學(xué)中反射系數(shù)、投射系數(shù)以及相位變化的計(jì)算結(jié)果一致,證實(shí)了結(jié)合gprMax正演模擬的地下目標(biāo)體解釋的可行性.