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

(1.信息工程學(xué)院 長(zhǎng)安大學(xué)， 西安 710064；2.淺層探測(cè)與測(cè)井研究實(shí)驗(yàn)室 休斯敦大學(xué)，德克薩斯州 美國(guó) 77004)

0 引言

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

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

1 TLM原理

1.1 Maxwell方程組

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

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

(1)

(2)

(3)

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

定義

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

(4)

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

1.2 傳輸線網(wǎng)絡(luò)方程與Maxwell方程間的關(guān)系

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

(5)

(6)

(7)

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

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

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

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

V=EzΔz

(8)

Ix=-ΔyHy

(9)

Iy=ΔxHx

(10)

(11)

(12)

(13)

(14)

Is=I(t)

(15)

1.3 傳輸線中的節(jié)點(diǎn)表達(dá)式

為了計(jì)算傳輸線中每個(gè)節(jié)點(diǎn)的電壓與電流，我們?cè)谶@里定義節(jié)點(diǎn)的阻抗為C，由圖3可知，此節(jié)點(diǎn)的阻抗可以定義為：

C=Cx+Cy+Cs

(16)

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

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

(17)

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

(18)

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

(19)

(20)

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

(21)

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

(22)

(23)

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

(24)

同理可得到傳輸線上的導(dǎo)納特性，即

(25)

(26)

(27)

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

Y∞=G

(28)

1.4 傳輸線網(wǎng)絡(luò)中的節(jié)點(diǎn)關(guān)系

在均勻介質(zhì)中，波的傳播過(guò)程[9]如圖4所示，等效為電壓波在傳輸線網(wǎng)格上傳播。圖4(a)中為傳播過(guò)程的初始狀態(tài)，設(shè)想輸入一個(gè)幅度為1V的沖激脈沖，它們?cè)谀┒诵纬傻碾妷悍瓷湎禂?shù)為(1/2，結(jié)果就形成了圖4(b)的結(jié)果，依次類推，圖4(c)為迭代兩次的結(jié)果，圖4表明了電壓波在傳輸線矩陣上的傳播過(guò)程。

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

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

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

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

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

(29)

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

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

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

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

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τ時(shí)的傳輸電壓Vi可由t=(p-1)τ時(shí)刻各點(diǎn)的反射電壓獲得。由公式(30)可以得到t=pτ時(shí)的反射電壓Vr。依此類推，我們可以迭代計(jì)算出不同時(shí)刻各個(gè)節(jié)點(diǎn)的電壓狀態(tài)。

2 算法實(shí)現(xiàn)

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

(1)設(shè)置模型。

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

(3)計(jì)算時(shí)間迭代，根據(jù)各個(gè)節(jié)點(diǎn)的入射電壓及反射電壓計(jì)算出節(jié)點(diǎn)總電壓。

(4)迭代結(jié)束，輸出計(jì)算結(jié)果。

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

3 數(shù)值算例

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

圖7 測(cè)試模型Fig．7 Test model

圖8 雷達(dá)入射脈沖Fig．8 Radar pulse

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

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

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

圖10 第一層介質(zhì)電導(dǎo)率不變，改變介電常數(shù)分別為4、6、8、16，得到的正演模擬結(jié)果圖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ì)的電阻率，通過(guò)仿真發(fā)現(xiàn)，在介電常數(shù)不變的情況下，電導(dǎo)率越大，物質(zhì)對(duì)電磁波的吸收越大，電磁波衰減地越快。

圖11 第一層介質(zhì)介電常數(shù)不變，改變電導(dǎo)率分別為0.005、0.05、0.5、1，得到的正演模擬結(jié)果圖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 結(jié)論

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

參考文獻(xiàn)：

[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].地球物理學(xué)報(bào), 2012,55(11): 3774-3785.

[3] 鄧世坤,王惠濂.探地雷達(dá)圖像的正演合成與偏移處理[J].地球物理學(xué)報(bào),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.