基于傳輸線法模擬探地雷達在多層介質模型中的傳播

李 婷，劉 策，2，郭 晨，張 勇

(1.信息工程學院 長安大學， 西安 710064；2.淺層探測與測井研究實驗室 休斯敦大學，德克薩斯州 美國 77004)

0 引言

探地雷達是一種近幾十年發展起來的一種地下目標的有效探測手段，通過分析檢測到的雷達數據，可以獲得測量物質的有效信息。正演數值模擬可以有效地模擬電磁波在地下結構中傳播，是解釋探地雷達數據的重要方法。探地雷達數值模擬方法主要包括時域有限差分法(Finite difference time domain method)[1]、有限單元法(Finite Element Method)[2]、射線追蹤法[3]、積分方程法[4]等，其中FDTD是應用最廣泛的方法，適合物理參數分布簡單或幾何特征規則的模型。事實上，與FDTD方法相對應的另一種時域方法即傳輸線矩陣法[5-7]，同樣具備FDTD方法的優點。

傳輸線法是由Johns和Beurle[7]提出的一種時間域分析電磁場方法。TLM方法基于Huygens波的傳播原理，通過研究時間和空間上的離散波在不同導波結構中的傳播情況，來獲得導波結構的種種傳輸特性。該方法最初應用于波導的不連續性及散射問題的分析，隨著計算機技術的發展，近年來TLM被廣泛應用在微波、數字電路仿真、電磁散射、雷達探測模擬、石油測井等具有波傳播特征的領域中。作者應用TLM方法進行均勻多層介質的正演模擬，通過典型激勵得到了合成的反射波，并與已有的FDTD算法進行比較，驗證了TLM方法的正確性和適用性。

1 TLM原理

1.1 Maxwell方程組

TLM方法是依據傳輸線矩陣中的電壓電流方程與Maxwell方程組的類比性，用傳輸線矩陣網格中的電壓和電流脈沖來模擬某區域中的電磁場傳播。如圖1所示，電磁脈沖發射出的電磁波垂直穿透介質層表面，脈沖探頭處的電池強度用I(t)表示，其中I(t)隨著時間的變化而發生變化。對TM模而言，此時區域的電磁場關系，場分量為Ez、Hx、Hy, Maxwell方程可表示如下：

圖1 位于(x0,y0)處的電磁脈沖源Fig．1 An electromagnetic-pulse source locate on the (x0,y0)

(1)

(2)

(3)

其中μ、σ和ε分別為磁導率、電導率和介電常數；Jsz為發射源的電流密度。

定義

Jsz=I(t)δ(x-x0)δ(y-y0)

(4)

其中δ為單位脈沖函數。

1.2 傳輸線網絡方程與Maxwell方程間的關系

在TLM方法中，均勻的二維空間可以用傳輸線來模擬。我們可以將圖1中的介質層劃分成網格為M×N的傳輸線網絡(圖2)。它實際上是一個半波系統，每隔Δl的距離都聯接著一對長為Δl/2的開路線(圖3)。則傳輸線方程可表示如下：

(5)

(6)

(7)

其中Lx和Ly分別為x、y軸上電感；C為電容；G為電阻器；Is是單位為安培電流源。

圖2 網格為M×N的傳輸線網絡Fig．2 M×N-transmission line network

圖3 網格單元等效電路圖Fig．3 Equivalent circuit of the TLM cell

通過對比式(1)-式(3)與式(5)-式(7)發現，傳輸線網絡方程與二維電磁場波動方程的關系如下：

V=EzΔz

(8)

Ix=-ΔyHy

(9)

Iy=ΔxHx

(10)

(11)

(12)

(13)

(14)

Is=I(t)

(15)

1.3 傳輸線中的節點表達式

為了計算傳輸線中每個節點的電壓與電流，我們在這里定義節點的阻抗為C，由圖3可知，此節點的阻抗可以定義為：

C=Cx+Cy+Cs

(16)

其中Cx、Cy分別為x、y軸上的電容值；Cs為并聯電容的電容值。

在傳輸線矩陣中，電壓或電流沿傳輸線矩陣傳播，其傳播速度v與傳輸線上的電容與電感有關，則x軸上的傳播速度可以定義為：

(17)

傳輸線x軸方向上節點間的傳播時間可表示為：

(18)

同理，y軸上節點間的傳播時間可表示為：

(19)

(20)

其中ε0、μ0分別為真空中的介電常數和磁導率；h為傳輸線長度，由文獻[8]中可知h需滿足以下條件：

(21)

由公式(11)、式(12)、式(19)、式(20)可以推導出Cx、Cy的表達式：

(22)

(23)

最終Cs可以由公式(13)、式(16)、式(22)、式(23)得到，即

(24)

同理可得到傳輸線上的導納特性，即

(25)

(26)

(27)

圖3中的并聯電容可以由半無限均勻無損耗傳輸線表示，即可得到

Y∞=G

(28)

1.4 傳輸線網絡中的節點關系

在均勻介質中，波的傳播過程[9]如圖4所示，等效為電壓波在傳輸線網格上傳播。圖4(a)中為傳播過程的初始狀態，設想輸入一個幅度為1V的沖激脈沖，它們在末端形成的電壓反射系數為(1/2，結果就形成了圖4(b)的結果，依次類推，圖4(c)為迭代兩次的結果，圖4表明了電壓波在傳輸線矩陣上的傳播過程。

我們將電壓波在傳輸線圖絡上的傳播具體化，以便形象地解釋節點電壓間的關系。圖5中各節點電路形式與圖3等效相同，分支(1)、分支(3)分別為X軸上的線電阻，分支(2)、分支(4)分別為y軸上的線電阻，分支(5)為開路節點，分支(6)等效為半無限均勻無損耗傳輸線。

圖4 電壓波在傳輸線網絡上的擴散過程Fig．4 The voltage wave diffusion on the transmission line network (a)初始狀態；(b)第一次迭代；(c)第二次迭代

設V(m,n,p)為t=p(時的網格(m,n)處的節點電壓，Vr(m,n,k,p) 、Vi(m,n,k,p)為第k個分支上的反射電壓和傳輸電壓，節點分支如圖5所示。

圖5 單元節點(m,n)與其相鄰節點Fig．5 Node(m,n ) and its adjacent node

在二維介質中，任意單元節點(m,n)都與四條TLM線相鄰，每條TLM線上存在正向波和反向波，而每條線上的正向波和反向波又分別來自四條TLM線上波的貢獻。則每個節點的電壓都可由式(29)表示。

(29)

其中Vr(m,n,k,p)和Vi(m,n,k,p)分別定義為5X1的矩陣，其中k=1,2,…。

節點反射電壓定義如下：

[Vr(m,n,k,p)]=[S][Vi(m,n,k,p)]

(30)

式中S為5×5的矩陣，其表示為：

(31)

而矩陣中的 Γ和T可定義成如下形式：

(32)

(33)

(34)

Tx=1+Γx

(35)

Ty=1+Γy

(36)

Tc=1+Γc

(37)

矩陣元素S11=Γy表示分支1上的反射系數，S21=Ty代表分支2與分支1之間的傳輸系數。以此類推，矩陣中對角線上的元素分別代表不同分支上的反射系數，其余元素則為不同分支間的反射系數。

每個節點的傳輸電壓由兩部分組成，一部分是來自該節點的反射電壓，另一部分為相鄰4個節點上的傳輸電壓。 每一個分支的傳輸電壓可表示如下：

Vi(m,n,1,p)=U1Vr(m,n,1,p-1)+

W1Vr(m,n-1,3,p-1)

(38)

Vi(m,n,2,p)=U2Vr(m,n,2,p-1)+

W2Vr(m-1,n,4,p-1)

(39)

Vi(m,n,3,p)=U3Vr(m,n,3,p-1)+

W3Vr(m,n+1,1,p-1)

(40)

Vi(m,n,4,p)=U4Vr(m,n,4,p-1)+

W4Vr(m+1,n,2,p-1)

(41)

反射系數U和傳輸系數W可表示如下：

(42)

(43)

(44)

(45)

Wk=1-Uk,k=1,…,4

(46)

由公式(38) - 公式(41) 中可知，t=pτ時的傳輸電壓Vi可由t=(p-1)τ時刻各點的反射電壓獲得。由公式(30)可以得到t=pτ時的反射電壓Vr。依此類推，我們可以迭代計算出不同時刻各個節點的電壓狀態。

2 算法實現

由第二節TLM的原理及其計算迭代公式可知，只要知道任意t=(p-1)τ時刻，各節點上的電壓大小與方向，就可以得到t=pτ時刻網絡中各個節點上的場值。關于TLM的具體實現過程如圖6所示，我們將TLM的基本算法步驟總結如下：

(1)設置模型。

(2)初始化節點反射電壓，其中Vr(m,n,k,0)=0 表示在任意網格t=0時的反射電壓為0，t=0時脈沖源處的電壓為Vr(m0,n0,k,0)=Vs

(3)計算時間迭代，根據各個節點的入射電壓及反射電壓計算出節點總電壓。

(4)迭代結束，輸出計算結果。

圖6 TLM算法流程圖Fig．6 Flow chart of TLM

3 數值算例

首先通過電磁波在多層介質中的傳播問題來驗證TLM算法的正確性。如圖7所示， 假設雷達發射端與接收端相距0.5 m，測試模型共有三層介質層，各層介電常數與電阻率分別為4、0.005；6、0.01；8、0.02。雷達發射脈沖采用中心頻率為1 GHz的正弦脈沖，如圖8所示。

圖7 測試模型Fig．7 Test model

圖8 雷達入射脈沖Fig．8 Radar pulse

通過對比圖9中TLM方法與FDTD方法的正演模擬結果，可以發現，TLM算法與FDTD方法可以達到同樣的正演效果。

圖9 TLM與FDTD正演模擬結果Fig．9 Forward modeling result of TLM and FDTD

改變圖7模型中第一層介質的介電常數值，其余參數保持不變。分別將介電常數改為4、6、8、16。通過TLM方法模擬電磁波在介質層的傳播，可以發現，在電導率不變的情況下，介質層介電常數的增加，將會導致電磁波間的衰減越來越大，而反射脈沖的強度越來越低。

圖10 第一層介質電導率不變，改變介電常數分別為4、6、8、16，得到的正演模擬結果圖Fig．10 The forwarding modeling result: the conductivity of 1st-layer is the same as the test model, changed dielectric constant to 4,6,8 and 16,respectively

同樣條件下，改變第一層介質的電阻率，通過仿真發現，在介電常數不變的情況下，電導率越大，物質對電磁波的吸收越大，電磁波衰減地越快。

圖11 第一層介質介電常數不變，改變電導率分別為0.005、0.05、0.5、1，得到的正演模擬結果圖Fig．11 The forwarding modeling result: conductivity of the 1st layer is the same as the test model, changed conductivity to 0.005、0.05、0.5 and 1,respectively

4 結論

作者應用二維TML方法模擬探地雷達電磁波在均勻多層介質模型中的傳播，通過對比本文算法與時域有限差分法模擬結果可以看出，兩種算法計算結果吻合非常好，這表明了TLM算法用于雷達波模擬的正確性和有效性。同時將TLM算法用于對已有模型的定性分析，結果表明：介電常數和電導率對雷達電磁波有很大影響，隨著介電常數和電導率的增大，雷達波的衰減越快，介質層對電磁波吸增加。

參考文獻：

[1] TING LI, CE LIU, CHEN GUO ,et al. Research on absorbing boundary condition in FDTD method [C]// Computer Science and Network Technology 2012 2nd International Conference.Changchun,China, 2012: 1574-1577.

[2] 馮德山,陳承申,王洪華. 基于混合邊界條件的有限單元法GPR正演模擬[J].地球物理學報, 2012,55(11): 3774-3785.

[3] 鄧世坤,王惠濂.探地雷達圖像的正演合成與偏移處理[J].地球物理學報,1993,3(4):528-536.

[4] SIMSEK E, LIU JIANGUO, LIU QINGHUO. A Spectral Integral Method (SIM) for Layered Media [J]. IEEE Transactions on antennas and propagation,2006,54(6): 3878- 3884.

[5] LIU CE, SHEN, LIANG C. Numerical Simulation of Subsurface Radar for Detecting Buried Pipes[J].IEEE Transactions on Geoscience and Remote Sensing, 1991,29(5):795- 798.

[6] HOEFER W J R.The transmission line matrix method theory and applications[J].IEEE Transactions on Microwave Theory and Techniques,1985,33(10):882-893.

[7] JOHNS P B,BEURLE R L. Numerical solution of 2-dimensional scattering problems using a transmission-line matrix[J]. Proceedings of the Institution of Electrical Engineers,1971,118(9): 1203-1208.

[8] LIU C, SHEN L.C. Response of Electromagnetic-Pulse Logging Sonde in Axially-Symmetrical Formation[J]. IEEE Transactions on Geoscience and Remote Sensing,1991,29:214-221.

[9] CHRISTOPOULOS C.The Transmission-Line Modeling Method TLM[M]. IEEE/OUP Series on Electromagnetic Wave Theory, 1995.