陳承申 ，馮德山

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

0 引言

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

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

1 GPR波動方程推導

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

▽×E=-?B/?t

(1)

▽×H=J+?D/?t

(2)

▽·D=ρ

(3)

▽·B=0

(4)

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

本構方程為式(5)：

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

(5)

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

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

(6)

考慮到式(7)：

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

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

(7)

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

整理得式(8)。

(8)

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

(9)

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

(10)

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

(11)

(12)

2 衰減特性分析

目前商用探地雷達大多使用調幅脈沖激勵源，其子波形式為：

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

(13)

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

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

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

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

2.1 非衰減媒質

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

MЁ+KE=SE

(14)

(15)

圖3分別給出非衰減媒質帶衰減項和不帶衰減項模擬計算結果。

2.2 衰減媒質

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

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

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

為了驗證在GPR正演模擬計算時，當σ?εω時，衰減項(σ/ε)(?E/?t)和(σ/ε)(?H/?t)的作用幾乎可以忽略，當地下介質導電性較強甚至很強的時候，電導率σ值大，σ?εω不成立，衰減項(σ/ε)(?E/?t)和(σ/ε)(?H/?t)的作用不可以忽略，作者自定義了探地雷達電磁波衰減比的概念及其計算公式。

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

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

(16)

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

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

4 結論

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

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

參考文獻：

[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ā)，徐波. 三維頻散介質中地質雷達信號的FDTD數值模擬[J]. 吉林大學學報:地球科學版, 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]. 地球物理學報, 1995, 36(4):499-506．

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

[10] 底青云，王妙月. 雷達波有限元仿真模擬[J]. 地球物理學報，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ā).頻散介質中地質雷達波傳播的數值模擬[J].地球物理學報,2007,50(1): 320-326.

[15] 陳承申. 探地雷達二維有限元正演模擬[D].長沙:中南大學，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.