探地雷達(dá)有限元數(shù)值模擬及衰減特性探討

陳承申 ，馮德山

(1.鐵道第三勘察設(shè)計(jì)院集團(tuán)有限公司，天津 300251;2.中南大學(xué) 地球科學(xué)與信息物理學(xué)院，長沙 410083)

0 引言

探地雷達(dá)(Ground penetrating radar,GPR)正演模擬一直是該領(lǐng)域理論研究的熱點(diǎn)[1-2]，通過對典型地質(zhì)模型正演計(jì)算結(jié)果的分析，加深了對GPR反射剖面的認(rèn)識與理解[3-4]，對GPR資料解釋有指導(dǎo)作用。在GPR正演模擬中，時(shí)域有限差分法[5-7](FDTD)因簡單靈活而被廣泛采用，但它不適合復(fù)雜的物性分界面；有限單元法[8](Finite Element Method，F(xiàn)EM)不需要計(jì)算內(nèi)部邊界條件，適用于物性復(fù)雜問題、求解過程規(guī)范化，具有廣泛的適應(yīng)性；沈飚[9]進(jìn)行了GPR簡單模型正演；底青云、王妙月[10-11]推導(dǎo)了含衰減項(xiàng)的探地雷達(dá)波有限元方程，實(shí)現(xiàn)了復(fù)雜形體的GPR波的FEM正演模擬和偏移處理；由于探地雷達(dá)波在地層中傳播能量會迅速衰減[12]，在進(jìn)行采集參數(shù)設(shè)置及數(shù)據(jù)處理時(shí)必須考慮到衰減媒質(zhì)對探地雷達(dá)波的影響，Tiao等[13]利用離散時(shí)域Galerkin進(jìn)行色散介質(zhì)中的GPR正演；劉四新、曾昭發(fā)[14]利用時(shí)域有限差分法，對探地雷達(dá)波在頻散介質(zhì)中傳播進(jìn)行了數(shù)值模擬。

作者通過對雷達(dá)波在衰減媒質(zhì)中傳播特性的研究，定義了衰減比的概念及其計(jì)算公式，使用帶衰減項(xiàng)的探地雷達(dá)波方程，對非衰減媒質(zhì)和衰減媒質(zhì)的中心電磁脈沖源模型模擬計(jì)算，這對雷達(dá)波在介質(zhì)中的傳播方式有更深刻的認(rèn)識[15]。

1 GPR波動方程推導(dǎo)

Maxwell方程組描述了電磁場的運(yùn)動學(xué)和動力學(xué)規(guī)律[16]，基于電磁波理論，高頻電磁波在媒質(zhì)中的傳播規(guī)律也應(yīng)服從Maxwell方程組，其方程組的微分形式可以表示為[17-18]：

▽×E=-?B/?t

(1)

▽×H=J+?D/?t

(2)

▽·D=ρ

(3)

▽·B=0

(4)

式中B為磁感應(yīng)強(qiáng)度(T)；E為電場強(qiáng)度(V/m)；D為電位移(C/m2)；H為磁場強(qiáng)度(A/m)；ρ為電荷密度(C/m3)；J為電流密度(A/m2)。

本構(gòu)方程為式(5)：

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

(5)

式中μ為磁導(dǎo)率(H/m);ε為介電常數(shù)(F/m);σ為電導(dǎo)率(S/m)。

把方程式(5)代入方程式(1)，并求旋度得到式(6)

(6)

考慮到式(7)：

▽×▽×A=-▽2A+▽(▽·A)

2.3.1 理論考核。以出科考核方式（選擇題與病案分析題）對學(xué)生進(jìn)行綜合評價(jià)，主要考察以下幾個(gè)方面：對急性心肌梗死這一嚴(yán)重胸痛疾病的理論知識的掌握，如病理學(xué)分析、診斷及鑒別診斷、治療；實(shí)踐患者的病史采集、體格檢查、疾病診斷及鑒別診斷、治療方案。結(jié)合帶教老師對學(xué)生平時(shí)表現(xiàn)的評定，在成績中予以一定的比例（20%）。即總成績=理論成績（80%）+平時(shí)表現(xiàn)成績（20%）。規(guī)定“優(yōu)秀”≥85分，“良好”≥75分，“及格”≥60分，“不及格”＜60分。

(7)

由于ε、μ、σ為坐標(biāo)的漸變函數(shù)，它們的空間導(dǎo)數(shù)是可以忽略的，則在式(7)中右邊的第二項(xiàng)▽(▽·A)=0

整理得式(8)。

(8)

假設(shè)雷達(dá)波激勵(lì)源如電場源SE或磁場源SH，代入到式(8)中得到式(9)。

(9)

同理對式(2)中第二式兩邊求旋度，可得式(10)。

(10)

式(9)與式(10)表明，磁場H和電場E及其分量滿足相同的微分方程，可以得到時(shí)空域GPR波的二維有限元方程：

(11)

(12)

2 衰減特性分析

目前商用探地雷達(dá)大多使用調(diào)幅脈沖激勵(lì)源，其子波形式為：

f(t)=t2e-α tsinω0t

(13)

式中ω0為GPR天線的主頻；α值的大小影響著電磁波在傳播過程中能量衰減的速率，這里α=0.93ω0。圖1為主頻100MHz的GPR子波圖。

圖1 100 MHz電磁脈沖子波Fig．1 100MHz electromagnetic pulse wavelet

為了說明雷達(dá)波在不同介質(zhì)中的傳播特性，建立模型如圖2所示，9.0 m×9.0 m大小的正方形均勻模型，單元大小為0.1 m×0.1 m，并把中心頻率為100 MHz的電磁脈沖激勵(lì)源置于模擬區(qū)域的正中心，采樣間隔為0.5 ns。

圖2 均勻模型示意圖Fig．2 Sketch map of homogeneous model

2.1 非衰減媒質(zhì)

設(shè)置模型介質(zhì)為干花崗巖，其中相對介電常數(shù)為ε=5，電導(dǎo)率為σ=10-7S/m，對于探地雷達(dá)電磁波，其頻率ω很高，有10-7=σ?εω=0.027 8，衰減項(xiàng)(σ/ε)(?E/?t)或者(σ/ε)(?H/?t)的作用幾乎可以忽略，此時(shí)式(9)與式(10)變?yōu)榧儾▌臃匠?，GPR有限元波動方程變?yōu)槭?14)與式(15)。

MЁ+KE=SE

(14)

(15)

圖3分別給出非衰減媒質(zhì)帶衰減項(xiàng)和不帶衰減項(xiàng)模擬計(jì)算結(jié)果。

2.2 衰減媒質(zhì)

模型設(shè)置為浸水花崗巖介質(zhì)，其中電導(dǎo)率為σ=0.01 S/m，相對介電常數(shù)為ε=7，當(dāng)媒質(zhì)為良導(dǎo)體或者近似于良導(dǎo)體時(shí)，σ較大，0.01=σ?εω=0.038 9不成立，衰減項(xiàng)(σ/ε)(?E/?t)或者(σ/ε)(?H/?t)的作用不能忽略，此時(shí)探地雷達(dá)波在這種介質(zhì)中傳播時(shí)，就會被媒質(zhì)吸收和發(fā)生頻散(圖4)。

圖3 GPR波在干花崗巖中傳播時(shí)的全波場快照Fig．3 The full field snapshot of the GPR wave transmitted in the dry granite(a)帶衰減項(xiàng)5 ns時(shí)刻全波場快照；(b)不帶衰減項(xiàng)5 ns時(shí)刻全波場快照；(c)帶衰減項(xiàng)15 ns時(shí)刻全波場快照；(d)不帶衰減項(xiàng)15 ns時(shí)刻全波場快照；(e)帶衰減項(xiàng)25 ns時(shí)刻全波場快照；(f) 不帶衰減項(xiàng)25 ns時(shí)刻全波場快照；(g)帶衰減項(xiàng)35 ns時(shí)刻全波場快照；(h) 不帶衰減項(xiàng)35 ns時(shí)刻全波場快照

3 干花崗巖模型與浸水花崗巖模型正演結(jié)果對比

為了驗(yàn)證在GPR正演模擬計(jì)算時(shí)，當(dāng)σ?εω時(shí)，衰減項(xiàng)(σ/ε)(?E/?t)和(σ/ε)(?H/?t)的作用幾乎可以忽略，當(dāng)?shù)叵陆橘|(zhì)導(dǎo)電性較強(qiáng)甚至很強(qiáng)的時(shí)候，電導(dǎo)率σ值大，σ?εω不成立，衰減項(xiàng)(σ/ε)(?E/?t)和(σ/ε)(?H/?t)的作用不可以忽略，作者自定義了探地雷達(dá)電磁波衰減比的概念及其計(jì)算公式。

概念 1 GPR電磁波在介質(zhì)中傳播時(shí)，可以計(jì)算得到t時(shí)刻帶衰減項(xiàng)的磁場 (或電場)峰值，設(shè)該峰值為a,不帶衰減項(xiàng)的磁場(或電場)峰值，設(shè)該峰值為b，b與a這兩個(gè)峰值的差值與不帶衰減項(xiàng)的磁場(或電場)峰值b的比值，這個(gè)比值稱為GPR電磁波衰減比[15]。

GPR電磁波衰減比可由式(16)計(jì)算得到[15]

(16)

利用GPR二維有限元正演模擬程序,對非衰減介質(zhì)(σ?εω成立，干花崗巖相關(guān)參數(shù)模型)和衰減介質(zhì)(σ?εω不成立，浸水的花崗巖相關(guān)參數(shù)模型)分別進(jìn)行數(shù)值模擬，并得出模擬結(jié)果(如表1所示)，截取電磁脈沖發(fā)射后的某些時(shí)刻的電場峰值，將這些值代入計(jì)算公式(16)，可以得到不同時(shí)刻不同介質(zhì)的衰減比。

通過討論分析，從表1中數(shù)據(jù)可以看出,干花崗巖模型GPR電磁波在45 ns時(shí),衰減比僅為0.005%，而在浸水花崗巖模型中，GPR電磁波在15ns時(shí)，其GPR電磁波衰減比就已經(jīng)為58.158%，到45 ns時(shí)更是高達(dá)94.877%，這說明了衰減項(xiàng)在非衰減媒質(zhì)GPR正演模擬時(shí)，是可以不用考慮其吸收作用，相反在衰減媒質(zhì)GPR正演模擬時(shí)，卻是必須研究的一項(xiàng)重要內(nèi)容。

4 結(jié)論

通過對兩種性質(zhì)相反的介質(zhì)模型的探地雷達(dá)正演模擬計(jì)算，說明了探地雷達(dá)電磁波在衰減介質(zhì)中傳播時(shí)能量會被介質(zhì)迅速吸收，衰減項(xiàng)的作用明顯；在非衰減介質(zhì)中傳播時(shí)能量被吸收的很少，衰減項(xiàng)是可以忽略不計(jì)的。作者通過給出的GPR電磁波衰減比概念，為理解雷達(dá)波在不同介質(zhì)中傳播的特性，提供了一種新的思路及方法。因此在探地雷達(dá)正演模擬計(jì)算時(shí)，必須考慮衰減項(xiàng)的衰減吸收作用，這樣的計(jì)算結(jié)果才更貼近實(shí)際情況。

表1 干花崗巖與浸水花崗巖模型正演結(jié)果對比表

參考文獻(xiàn)：

[1] FENG D S, DAI Q W. GPR numerical simulation of full wave field based on UPML boundary condition of ADI-FDTD[J]. NDT&E International, 2011,44(6): 495-504.

[2] SHAARI A, AHMAD R S, CHEW T H. Effects of antenna-target polarization and target-medium dielectric contrast on GPR signal from non-metal pipes using FDTD simulation[J]. NDT & E International, 2010, 43(5): 403-408.

[3] JAMES IRVING, ROSEMARY KNIGHT. Numerical modeling of ground-penetrating radar in 2-D using MATLAB[J] .Computers & Geosciences, 2006,32 (9): 1247-1258．

[4] GIANNOPOULOS A. Modelling ground penetrating radar by GprMax[J]. Construction and Building Materials , 2005, 19(10): 755-762．

[5] 劉四新，曾昭發(fā)，徐波. 三維頻散介質(zhì)中地質(zhì)雷達(dá)信號的FDTD數(shù)值模擬[J]. 吉林大學(xué)學(xué)報(bào):地球科學(xué)版, 2006, 36(1)：123-127．

[6] BERGMANN T，ROBERTSSON J O.A，HOLLIGER K. Numerical properties of staggered finite-difference solutions of Maxwell’s equations for ground-penetrating radar modeling[J]. Geophysical Research Letters, 1996, 23(1):45-48．

[7] CARCIONE，JOSE M. Radiation patterns for 2-D GPR forward modeling[J]. Geophysics, 1998, 63(2):424-430．

[8] 王妙月，郭亞曦，底青云. 二維線性流變體波的有限元模擬[J]. 地球物理學(xué)報(bào), 1995, 36(4):499-506．

[9] 沈飚. 探地雷達(dá)波波動方程研究及其正演模擬[J]. 物探化探計(jì)算技術(shù)，1994，16(1)：29-33．

[10] 底青云，王妙月. 雷達(dá)波有限元仿真模擬[J]. 地球物理學(xué)報(bào)，1999,42(6):818-825．

[11] DI QING-YUN，WANG MIAO-YUE. Migration of ground-penetrating radar data with a finite-element method that considers attenuation and dispersion[J]. Geophysics, 2004, 69(2):472-477．

[12] RICHARD D, MILLER.Vertical resolution of a seismic survey in stratigraphic sequences less than loom deep in Sortheastern Kansas[J]. Geophysics,1995,60(2):423-430．

[13] LU TIAO, CAI WEI, ZHANG PING-WEN. Discontinuous Galerkin Time-Domain Method for GPR Simulation in Dispersive Media[J]. IEEE Transactions on Geoscience and Remote Sensing, 2005, 43(1):72-80．

[14] 劉四新,曾昭發(fā).頻散介質(zhì)中地質(zhì)雷達(dá)波傳播的數(shù)值模擬[J].地球物理學(xué)報(bào),2007,50(1): 320-326.

[15] 陳承申. 探地雷達(dá)二維有限元正演模擬[D].長沙:中南大學(xué)，2011.

[16] YEE K S. Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media[J].IEEE Transactions on Antennas and Propagation,1996,14 (3):302-307.

[17] TAFLOVE A. Re-inventing electromagnetics:supercomputing solution of Maxwell’s equations via direct time integration on space grids[J]. Computing Systems in Engineering, 1992, 3(1):153-168.

[18] ALVAREZ G B, LOULA A F D, DUTRA DO CARMO E G,et al . A discontinuous finite element formulation for the Helmholtz equation. Comput[J]. Methods Appl. Mech. Engrg., 2006,195:4018-4035.