有耗色散地質(zhì)介質(zhì)中電磁波傳播特性的FDTD計(jì)算分析

郭 江 曹俊興 何曉燕

(成都理工大學(xué)信息工程學(xué)院,成都610059)

電磁法是重要的地球物理探測(cè)方法之一。大部分地質(zhì)介質(zhì),特別是土壤,對(duì)電磁波而言是有耗色散介質(zhì)。但限于問題的復(fù)雜性,在一般的電磁法數(shù)據(jù)處理解釋中通常并不考慮介質(zhì)色散的影響。為提高電磁探測(cè)數(shù)據(jù)的解釋精度,需要將色散的影響納入考慮。為在電磁探測(cè)數(shù)據(jù)的處理中考慮色散的影響,首先需要發(fā)展色散影響的定量計(jì)算方法,定量分析色散對(duì)電磁波場(chǎng)的影響。Teixeira等(1998,2000)[1,2]和 Wang等(2000)[3]先后研究發(fā)展了基于時(shí)域有限差分法(FDTD)和完全匹配層(PML)吸收邊界的非均勻色散良導(dǎo)介質(zhì)中高頻電磁波場(chǎng)數(shù)值計(jì)算方法。自Yee于1966年提出求解麥克斯韋方程的FDTD格式[4]以來,經(jīng)過半個(gè)多世紀(jì)的發(fā)展,目前已成為使用最廣的電磁波場(chǎng)數(shù)值計(jì)算方法之一[5～7]。但用FDTD方法處理色散介質(zhì)中的電磁場(chǎng)分布時(shí),由于介質(zhì)的介電常數(shù)是頻率的函數(shù),其本構(gòu)關(guān)系在時(shí)域成為卷積關(guān)系,給直接計(jì)算帶來困難,需采用遞推式來計(jì)算時(shí)域卷積[1,2],這無疑將增加計(jì)算存儲(chǔ)量和過程。本文用推廣的德拜介質(zhì)參數(shù)方程描述色散介質(zhì)的本構(gòu)關(guān)系,采用將此本構(gòu)關(guān)系和場(chǎng)量關(guān)系直接變換為時(shí)域微分方程的方法[6],導(dǎo)出了有耗媒色散介質(zhì)中電磁場(chǎng)的時(shí)域有限差分計(jì)算格式。新的計(jì)算格式適用于分析一般色散介質(zhì)的電磁特性。避免了時(shí)域卷積關(guān)系的計(jì)算。采用完全匹配層(PML)吸收邊界[8～11]處理邊界條件,模擬分析了介質(zhì)色散對(duì)電磁波傳播的影響,認(rèn)為一些地質(zhì)介質(zhì)的色散會(huì)對(duì)電磁波的傳播產(chǎn)生不可忽視的影響。

1 色散有耗介質(zhì)中場(chǎng)分量的FDTD差分式的推導(dǎo)

1.1 基本微分方程

各向同性色散介質(zhì)中的Maxwell旋度方程組可以寫為(僅考慮介質(zhì)的介電常數(shù)為色散的情況,設(shè)介質(zhì)的磁導(dǎo)率為常數(shù))：

介質(zhì)的本構(gòu)關(guān)系為：

根據(jù)德拜關(guān)系式我們將有耗色散介質(zhì)的介電常數(shù)表示為[6]

將式(4)代入式(3)可得在頻域中有

考慮到對(duì)任意時(shí)變場(chǎng)有

可得到(5)式在時(shí)域?yàn)?/p>

方程(1)、(2)和(6)構(gòu)成了完整的電磁場(chǎng)時(shí)域方程組。當(dāng)σ=0時(shí),該方程組表示無電導(dǎo)損耗的色散介質(zhì)的場(chǎng)方程;當(dāng)時(shí),由式(2)、(3)、(6)有

其中 ：ε=ε0ε∞。

因此,方程(1)、(7)、(8)是不計(jì)介質(zhì)色散時(shí),一般電導(dǎo)損耗介質(zhì)的場(chǎng)方程。

將介質(zhì)的電導(dǎo)損耗在介質(zhì)的介電常數(shù)中考慮,會(huì)給色散介質(zhì)中 FDTD公式的建立帶來方便。它使我們可以將介質(zhì)參數(shù)ε*(ω)和頻域場(chǎng)分量直接變換為時(shí)域的場(chǎng)方程,以便進(jìn)一步得到時(shí)域場(chǎng)分量的有限差分式,避免了由于時(shí)域卷積形式的出現(xiàn)所需的用迭代法求卷積值而增加的計(jì)算存儲(chǔ)量和過程。

1.2 FDTD差分方程式

用以上導(dǎo)出的完整的時(shí)域場(chǎng)方程(1)、(2)和(6),在直角坐標(biāo)系中取Yee元胞[4]為空間電磁場(chǎng)FDTD離散的基本單元,設(shè) Δx=Δ y=Δ z=δ,運(yùn)用中心有限差分方程式,離散(1)、(2)、(6)式,可得到9個(gè)場(chǎng)分量(Ex,Ey,Ez,Dx,Dy,Dz,Hx,)的差分方程式,其中x分量為：

其中：

1.3 單軸各向異性完全匹配吸收層中場(chǎng)量的差分格式

Gedney各向異性PM L介質(zhì)中場(chǎng)方程仍為麥克斯韋方程的形式[10,11]

其中

為了減少截?cái)噙吔缫鸬姆瓷湔`差,ση和χη取如下漸變形式：

應(yīng)注意的是,在吸收邊界的不同交疊區(qū)(14)式具有不同的形式。在圖1所示的區(qū)域1,完全匹配介質(zhì)界面垂直x軸,此處sy=sz=1;在區(qū)域2,sx==1;在區(qū)域3=1;在區(qū)域4,=1;在區(qū)域 5,=1;在區(qū)域6,取(14)形式(即若PML介質(zhì)與計(jì)算區(qū)域的分界面垂直 η軸則僅sη≠1,如 1和2區(qū)域;若PML介質(zhì)與計(jì)算區(qū)域有兩個(gè)分界面分別垂直 x軸或y軸,則≠1,≠1,如區(qū)域5;其他區(qū)域以此類推)。在區(qū)域1,由(12)、(14)式可得

圖1 各向異性完全匹配層區(qū)域示意圖Fig.1 sketch of the region of the anisotropic PM L

借助輔助變量

由(15)式的第二和第三式,通過頻域到時(shí)域的轉(zhuǎn)換后,容易得到Dy和Dz的時(shí)域推進(jìn)計(jì)算公式,其y分量為：

由(15)式第一式,通過頻域到時(shí)域的轉(zhuǎn)換,采用常用的FDTD離散方法可以得到的時(shí)域推進(jìn)計(jì)算公式,與(10)式相同。從D′x到Ex的計(jì)算則采用以下方法,首先將式(17)轉(zhuǎn)換到時(shí)域,再利用中心差分格式離散可得到的二階精度的差分表達(dá)式

使用(16)、(11)、(18)式,可完成從H(經(jīng)D)到E的時(shí)域推進(jìn)計(jì)算。從E到H的推進(jìn)計(jì)算可由(13)式用類似方法得到。

對(duì)交角區(qū)域6,(12)式為

為便于將(19)式過渡到時(shí)域FDTD格式,也為了避免方程中出現(xiàn)三階以上的時(shí)間偏導(dǎo)數(shù),引入輔助變量

2 數(shù)值計(jì)算及分析

2.1 介質(zhì)色散對(duì)電磁波場(chǎng)傳播特性影響分析

考慮二維色散介質(zhì)中Blackman-Harris脈沖源激勵(lì)的電磁場(chǎng)的分布,設(shè)入射波上限頻率為1.5×108Hz,脈沖寬度為2.5×10-8s,空間離散采用均勻網(wǎng)格剖分,網(wǎng)格尺寸取Δx=Δy=δ=0.166 m,時(shí)間步長取 Δ t=3.9×10-10s。計(jì)算空間為136×136個(gè)離散單元,場(chǎng)源設(shè)置在(68,68)網(wǎng)格處,外邊界采用Gedney各向異性PML吸收邊界截?cái)唷M L的電導(dǎo)率表達(dá)式中

為研究介質(zhì)色散對(duì)電磁波傳播特性的影響,我們?nèi)∮?jì)算域中距場(chǎng)源20δ處的點(diǎn)P處的電場(chǎng)(68,48)為分析對(duì)象。圖3展示了點(diǎn)P處隨時(shí)間的變化。分析圖3所示的曲線,可知色散介質(zhì)中由于色散的影響電磁波的傳播位相落后于非色散介質(zhì)(即在時(shí)間上有滯后),脈沖形狀發(fā)生畸變,寬度有較明顯的展寬,幅度減小。所以介質(zhì)的色散會(huì)使電磁脈沖展寬、位相滯后、幅度減小,這在目標(biāo)識(shí)別中是必須考慮的問題。

為考察色散介質(zhì)電導(dǎo)率變化對(duì)電磁波傳播的影響,我們?cè)?36×136的計(jì)算域,模擬計(jì)算了不同電導(dǎo)率土壤(參數(shù)=6×10-10s;電導(dǎo)率分別為 σ1=0,0.001 6,0.008,0.016 S/m)中電磁波的傳播。場(chǎng)源仍用Blackman-Harris脈沖,置于計(jì)算域中心(68,68)。場(chǎng)域內(nèi)點(diǎn)P(68,63)處Ez隨時(shí)間的變化曲線示于圖4。圖4表明色散介質(zhì)電導(dǎo)率的變化只會(huì)改變脈沖的幅度,對(duì)其傳播過程中的位相、波形不會(huì)產(chǎn)生明顯的影響。

2.2 分層介質(zhì)中電磁脈沖的傳播

圖2 色散土壤(A)和非色散土壤(B)中不同時(shí)刻的波場(chǎng)(Ez)快照Fig.2 distributions at different time steps in the media with dispersion(A)and without dispersion(B)

地質(zhì)雷達(dá)等電磁探測(cè)所涉及的地質(zhì)介質(zhì)可以用分層介質(zhì)模型來近似。為研究介質(zhì)色散對(duì)電磁探測(cè)的影響,我們計(jì)算了三層色散與非色散介質(zhì)中的反射場(chǎng)。三層介質(zhì)之頂層(層1)為空氣;層2(模擬土壤)的電性參數(shù)為層3的電性參數(shù)為0.000 4 S/m。計(jì)算空間剖分為310×235個(gè)網(wǎng)格單元,層1與層2的分界面位于 z=157δ,層 2與層3的分界面位于 z=187δ,場(chǎng)源置于(155,155),觀測(cè)點(diǎn)設(shè)在點(diǎn)P(160,156)處。計(jì)算結(jié)果如圖5中的虛線所示。圖5中的實(shí)線是層2作為非色散有耗介質(zhì)時(shí)的反射場(chǎng)=0.001 6 S/m)。圖5曲線上的第一個(gè)波峰是自由空間與層2界面的反射波,第二個(gè)波峰為層2與層3界面的反射波。圖5表明,將土壤作為色散介質(zhì)與非色散介質(zhì),反射電磁波的差別是比較大的。當(dāng)將土壤(層2)作為色散介質(zhì)時(shí),空氣和土壤界面的反射波幅要較作為非色散介質(zhì)時(shí)的反射波幅為大,但相位相同;電磁波脈沖穿過色散介質(zhì)時(shí),其速度會(huì)降低,旅行時(shí)會(huì)增加。因此,電磁探測(cè)數(shù)據(jù)的處理應(yīng)將介質(zhì)的色散納入考慮。

圖3 介質(zhì)色散對(duì)Ez隨時(shí)間變化的影響Fig.3 Curves showing the effect of the dispersion on Ez

圖4 色散介質(zhì)電導(dǎo)率變化對(duì)Ez隨時(shí)間變化的影響Fig.4 Curves showing the effect of the conductivity change on Ez

圖5 觀測(cè)點(diǎn)反射場(chǎng)分量Ey隨時(shí)間的變化Fig.5 Curves showing the Eychange with time at the receiver

3 結(jié)論

a.用推廣的德拜介質(zhì)參數(shù)方程描述一般有耗色散介質(zhì)的本構(gòu)關(guān)系,采用將介質(zhì)本構(gòu)關(guān)系、場(chǎng)量關(guān)系直接變換為時(shí)域微分方程的方法,導(dǎo)出了有耗介質(zhì)中電磁場(chǎng)的時(shí)域有限差分計(jì)算公式,避免了時(shí)域卷積關(guān)系的出現(xiàn),簡化了色散介質(zhì)中電磁場(chǎng)問題的處理。

b.將Gedney單軸各向異性完全匹配邊界條件應(yīng)用于有耗色散介質(zhì)中電磁場(chǎng)的模擬計(jì)算,導(dǎo)出了吸收層中的FDTD計(jì)算公式。數(shù)值計(jì)算結(jié)果表明,該邊界條件的吸收效果良好。

c.色散介質(zhì)會(huì)使電磁脈沖在傳播過程中波形展寬、位相滯后、幅度減小。因此,電磁探測(cè)數(shù)據(jù)的處理解釋應(yīng)將介質(zhì)色散的影響納入考慮。

[1]T EIXEIRA F L,CHEW W C,STRAKA M,et al.Finite-difference time-domain simulation of ground penetrating radar on dispersive,inhomogeneous,and conductive soils[J].IEEE Trans on Geoscience and Sensing,1998,36(6)：1928-1936.

[2]TEIXEIRA F L,CHEW W C.Finite-difference computation of transient electromagnetic waves for cylin-drical geometries in complex media[J].IEEE Trans on Geoscience and Sensing,2000,38(4)：1530-1543.

[3]WANG T,ORISTAGLIO M L.3-D simulation of GRP surveys over pipes in dispersive soils[J].Geophysics,2000,65(5)：1560-1568.

[4]YEE K S.Numerical solution initial boundary problems involving Maxwell's equations in isotropic media[J].IEEE Trans on AP,1966,14(4)：302-307.

[5]TAFLOVE A.Computational Electrodynamics：The Finite-Difference Time-Domain Method[M].Boston：Artech House,1995.

[6]高本慶.時(shí)域有限差分法FDTD Method[M].北京：國防工業(yè)出版社,1995.

[7]葛德彪,閆玉波.電磁波時(shí)域有限差分方法[M].西安：西安電子科技大學(xué)出版社,2003.

[8]BERENGER J P.A perfectly matched layer for the absorption of electromagnetic waves[J].J Comput Phys,1994,114(2)：185-200.

[9]SACKS Z S,KINGSLAND D M,LEE D M,et al.A perfectly matched anisotropic absorber for use as an absorbing boundary condition[J].IEEE T rans Antennas Propagat,1995,AP-43(12)：1460-1463.

[10]GEDNEY S D.An anisotropic perfectly matched layer absorbing media for the truncation of FDTD lattices[J].IEEE Trans Antennas Propagat,1996,AP-44(12)：1630-1639.

[11]孔繁敏,李康,劉新,等.波動(dòng)方程FDTD算法的PM L吸收邊界條件的實(shí)現(xiàn)與驗(yàn)證[J].微波學(xué)報(bào),2004,20(1)：1-5.

[12]王禹,袁乃昌.色散媒質(zhì)中ADI-FDTD的PML[J].電子與信息學(xué)報(bào),2005,27(10)：1677-1680.

[13]馬孝義,馬建倉.土壤水分介電測(cè)量的頻率上限分析[J].水土保持研究,2002,9(2)：82-86.