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

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

0 引言

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

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

1 TLM原理

1.1 Maxwell方程組

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

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

(1)

(2)

(3)

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

定義

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

(4)

其中δ為單位脈沖函數(shù)。

1.2 傳輸線網(wǎng)絡方程與Maxwell方程間的關系

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

(5)

(6)

(7)

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

圖2 網(wǎng)格為M×N的傳輸線網(wǎng)絡Fig．2 M×N-transmission line network

圖3 網(wǎng)格單元等效電路圖Fig．3 Equivalent circuit of the TLM cell

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

V=EzΔz

(8)

Ix=-ΔyHy

(9)

Iy=ΔxHx

(10)

(11)

(12)

(13)

(14)

Is=I(t)

(15)

1.3 傳輸線中的節(jié)點表達式

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

C=Cx+Cy+Cs

(16)

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

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

(17)

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

(18)

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

(19)

(20)

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

(21)

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

(22)

(23)

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

(24)

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

(25)

(26)

(27)

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

Y∞=G

(28)

1.4 傳輸線網(wǎng)絡中的節(jié)點關系

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

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

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

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

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

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

(29)

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

節(jié)點反射電壓定義如下：

[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上的反射系數(shù)，S21=Ty代表分支2與分支1之間的傳輸系數(shù)。以此類推，矩陣中對角線上的元素分別代表不同分支上的反射系數(shù)，其余元素則為不同分支間的反射系數(shù)。

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

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)

反射系數(shù)U和傳輸系數(shù)W可表示如下：

(42)

(43)

(44)

(45)

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

(46)

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

2 算法實現(xiàn)

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

(1)設置模型。

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

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

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

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

3 數(shù)值算例

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

圖7 測試模型Fig．7 Test model

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

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

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

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

圖10 第一層介質(zhì)電導率不變，改變介電常數(shù)分別為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

同樣條件下，改變第一層介質(zhì)的電阻率，通過仿真發(fā)現(xiàn)，在介電常數(shù)不變的情況下，電導率越大，物質(zhì)對電磁波的吸收越大，電磁波衰減地越快。

圖11 第一層介質(zhì)介電常數(shù)不變，改變電導率分別為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方法模擬探地雷達電磁波在均勻多層介質(zhì)模型中的傳播，通過對比本文算法與時域有限差分法模擬結果可以看出，兩種算法計算結果吻合非常好，這表明了TLM算法用于雷達波模擬的正確性和有效性。同時將TLM算法用于對已有模型的定性分析，結果表明：介電常數(shù)和電導率對雷達電磁波有很大影響，隨著介電常數(shù)和電導率的增大，雷達波的衰減越快，介質(zhì)層對電磁波吸增加。

參考文獻：

[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.