橫向同性分層介質(zhì)三維電磁模擬的高階DTA算法

陳桂波，張燁，潘軼群

（1.長(zhǎng)春理工大學(xué) 理學(xué)院，長(zhǎng)春 130022；2.長(zhǎng)春市農(nóng)業(yè)機(jī)械研究院，長(zhǎng)春 130062）

矩量法是進(jìn)行電磁場(chǎng)三維數(shù)值模擬的一種常用而有效的方法，在農(nóng)業(yè)遙感、探地雷達(dá)、資源勘探等領(lǐng)域具有十分廣泛的應(yīng)用［1，2］。該方法首先建立一個(gè)關(guān)于散射體內(nèi)電場(chǎng)或電流分布的第二類弗雷德霍姆積分方程，然后將積分方程離散形成一個(gè)矩陣方程，求解該矩陣方程可以得到散射體內(nèi)的電場(chǎng)分布，進(jìn)而得到空間任意位置的電磁場(chǎng)分布［3］。

由于矩量法一般只針對(duì)有限大小的散射體進(jìn)行離散，因此與需要全空間離散的有限差分與有限元法相比，矩量法具有離散數(shù)目少、求解精度高等優(yōu)點(diǎn)。但該方法在使用過程中也存在一個(gè)不容忽視的問題，由于積分方程的離散化系數(shù)矩陣為復(fù)數(shù)致密矩陣，存儲(chǔ)和計(jì)算該矩陣需要大量的計(jì)算機(jī)內(nèi)存和計(jì)算時(shí)間，致使矩量法在實(shí)際應(yīng)用中往往受到較大的限制。為了克服這個(gè)問題，一些學(xué)者提出了積分方程的近似解法，可達(dá)到滿意的計(jì)算精度并且具有占用內(nèi)存少、計(jì)算效率高等固有優(yōu)點(diǎn)，在近年來(lái)逐漸發(fā)展成為計(jì)算電磁學(xué)領(lǐng)域中的一個(gè)研究熱點(diǎn)。其中由量子力學(xué)理論發(fā)展而來(lái)的Born近似是最常見的一種近似方法［4］，該方法雖然形式簡(jiǎn)單、易于編程，但使用條件極為苛刻，往往要求其有頻率較低以及散射體尺寸較小等弱散射條件。Habashy等人在Born近似的基礎(chǔ)上提出了擴(kuò)展Born近似算法［5］，該方法通過引入散射張量將散射體內(nèi)電場(chǎng)與入射電場(chǎng)聯(lián)系起來(lái)，其計(jì)算精度和適用范圍都要優(yōu)于Born近似，但該方法存在一個(gè)明顯不足就是當(dāng)散射體距發(fā)射源較近以及散射體內(nèi)電場(chǎng)分布極不均勻時(shí)計(jì)算精度較差；Zhdanov等人分別提出了QL和QA近似［6，7］，并進(jìn)一步給出了相應(yīng)的高階級(jí)數(shù)展開形式；Song等人基于局部非線性近似原理提出了一種依賴于入射源的對(duì)角張量近似（DTA）算法［8］，數(shù)值實(shí)例證明該方法相對(duì)于前幾種方法具有更高的計(jì)算精度和更廣泛的適用性，并且已經(jīng)成功用于層狀介質(zhì)的三維反演成像［9］。但是，上述幾種方法的應(yīng)用目前還都只局限于各向同性介質(zhì)，而關(guān)于各向異性介質(zhì)條件下相應(yīng)的近似算法研究卻鮮有報(bào)道。各向異性介質(zhì)的電磁模擬與特性分析是在實(shí)際中經(jīng)常會(huì)遇到的問題，因此本文將相對(duì)較新、性能較好的DTA算法推廣應(yīng)用到各向異性分層介質(zhì)（為簡(jiǎn)單起見，本文只考慮橫向同性分層介質(zhì)），為進(jìn)一步提高DTA求解積分方程的計(jì)算精度，本文提出一種滿足壓縮映射條件的高階DTA級(jí)數(shù)算法。

文中首先給出求解積分方程的DTA算法基本原理，然后基于壓縮算子理論得到一種總是收斂的高階DTA級(jí)數(shù)解，最后通過數(shù)值實(shí)例對(duì)比檢驗(yàn)了本文算法的有效性。

1 DTA近似基本原理

圖1 三維模型示意圖

DTA算法用于求解積分方程的基本原理是引入一個(gè)與入射源有關(guān)的對(duì)角張量將散射體內(nèi)的散射電場(chǎng)與入射電場(chǎng)聯(lián)系起來(lái)，即：

其中 Γ=diag(γx,γy,γz)為對(duì)角張量，其各個(gè)元素表達(dá)式可利用文獻(xiàn)［8］中的局部非線性近似方法來(lái)得到：

這里Ebξ為入射電場(chǎng)Eb的ξ分量，而EBξ為矢量

的ξ分量。利用式（3），可將散射體內(nèi)的電場(chǎng)寫為：

再由方程（1）和（2）就可以得到空間任意位置的電磁場(chǎng)。

2 高階DTA級(jí)數(shù)

方程（1）也可以寫成算子方程的形式：

根據(jù)積分方程理論，方程（7）的解可由迭代法得到：

由上式容易看出，若E(0)=Eb，則E(1)即為傳統(tǒng)Born近似解，k＞1時(shí)方程（1）的解可以展開成高階Born級(jí)數(shù)的形式；若，則 E(1)為DTA解，同理在這種情況下k＞1時(shí)方程（1）的解也可以展開成高階級(jí)數(shù)的形式，稱之為高階DTA級(jí)數(shù)。根據(jù)算子理論可知，應(yīng)用迭代法求解方程（7）只有在A為壓縮算子，也就是的L2范數(shù)小于1的條件下才是收斂的。遺憾的是，算子A只有在散射體尺寸較小以及電導(dǎo)率差異較小等弱散射條件下才是壓縮算子。為了克服這個(gè)問題，利用文獻(xiàn)［11］中介紹的線性變換方法，在方程（1）兩端同時(shí)左乘一個(gè)對(duì)角張量，可以得到關(guān)于參量E的算子方程：

其中

由文獻(xiàn)［3］可知，算子 C是 L2(Ω)上的壓縮算子，也就是說(shuō)算子方程（10）滿足壓縮映射條件，進(jìn)而應(yīng)用迭代法將方程（10）的解展開成高階級(jí)數(shù)對(duì)于任意電導(dǎo)率分布等參數(shù)條件下總是收斂的。為了驗(yàn)證這一點(diǎn)，考慮具有圖2所示結(jié)構(gòu)的三層橫向同性介質(zhì)模型，中間層厚度為20m，其余兩層均為半空間。在圖示直角坐標(biāo)系中，大小為30m×120m×90m的散射體位于第三層介質(zhì)中，其中心坐標(biāo)為（0，0，75）。假定第一層介質(zhì)為空氣，第二層介質(zhì)的橫向和縱向電導(dǎo)率分別為0.1S/m和0.02S/m，第三層介質(zhì)的橫向和縱向電導(dǎo)率分別為0.01S/m和0.005S/m。在（－10m，0，0）處有一個(gè)沿x方向單位大小的電偶極子，發(fā)射頻率為1000Hz，接收點(diǎn)位于（10m，0，0）。分別采用基于方程（7）的傳統(tǒng)迭代法和基于方程（10）的壓縮迭代法計(jì)算了散射電場(chǎng)x分量的相對(duì)誤差（這里表示應(yīng)用嚴(yán)格矩量法的計(jì)算結(jié)果［3］），圖3（a）～（c）給出了散射體的電導(dǎo)率為0.1S/m、0.05S/m、0.0005S/m時(shí)兩種方法計(jì)算的相對(duì)誤差隨迭代次數(shù)的變化曲線?？梢钥闯?，傳統(tǒng)迭代法的計(jì)算結(jié)果均不收斂，尤其是在散射體電導(dǎo)率大于分層介質(zhì)電導(dǎo)率時(shí)，散射體電導(dǎo)率越高，發(fā)散則越嚴(yán)重，而壓縮迭代法在三種情況下的計(jì)算結(jié)果都是收斂的，可見本文引入的壓縮算子對(duì)于各種情況下積分方程迭代收斂性能的提高作用是非常明顯的。

圖2 三層橫向同性介質(zhì)模型

3 數(shù)值計(jì)算結(jié)果

圖3 傳統(tǒng)迭代法與壓縮迭代法的相對(duì)誤差

考慮一個(gè)與圖2具有相同幾何參數(shù)的三層介質(zhì)模型，并且依舊假定第一層介質(zhì)為空氣，而第二層和第三層介質(zhì)的橫向電導(dǎo)率分別為0.1S/m和0.01S/m，縱向電導(dǎo)率分別為0.02S/m和0.005S/m。單位大小的垂直磁偶極子可沿x軸任意移動(dòng)，接收點(diǎn)也位于x軸上，發(fā)射源與收收點(diǎn)之間的距離為150m，并且選擇發(fā)射源與接收點(diǎn)的中心作為記錄點(diǎn)，發(fā)射頻率為1000Hz。分別采用傳統(tǒng)Born近似、DTA近似、本文的高階DTA級(jí)數(shù)方法以及嚴(yán)格矩量法［3］計(jì)算了散射磁場(chǎng)的垂直分量分布，圖4和圖5分別給出了當(dāng)散射體電導(dǎo)率為0.1S/m和0.001S/m時(shí)的計(jì)算結(jié)果。從圖4中曲線可以看出，本文的高階DTA級(jí)數(shù)法收斂速度較快，3階展開就幾乎可以達(dá)到精確結(jié)果，并且DTA近似（0階DTA級(jí)數(shù)）的計(jì)算精度要高于傳統(tǒng)Born近似。從計(jì)算效率來(lái)看，在同一臺(tái)計(jì)算機(jī)上，嚴(yán)格矩量法計(jì)算圖中31個(gè)節(jié)點(diǎn)需要3分鐘左右，而3階DTA近似則小于1分鐘，可見本文算法在保證計(jì)算精度的同時(shí)可使積分方程的求解效率得到有效提高。

觀察圖5曲線可以發(fā)現(xiàn)，幾種方法計(jì)算結(jié)果之間的差別沒有圖4曲線那樣明顯，這是由于當(dāng)散射體電導(dǎo)率為0.001S/m時(shí)電導(dǎo)率差異相對(duì)于分層介質(zhì)的電導(dǎo)率來(lái)說(shuō)較小，產(chǎn)生的散射較弱，導(dǎo)致散射體內(nèi)入射場(chǎng)在總場(chǎng)中所占有的比例較大，因此各種方法的近似程度均比目標(biāo)體電導(dǎo)率為0.1S/m時(shí)要高。在實(shí)際計(jì)算中，可以根據(jù)散射體的性質(zhì)來(lái)選擇較少階數(shù)的DTA展開，以便在保證計(jì)算精度的同時(shí)進(jìn)一步提高計(jì)算速度。

圖4 散射體為低阻時(shí)的散射場(chǎng)

圖5 目標(biāo)體為高阻時(shí)的散射場(chǎng)

4 結(jié)論

本文提出了一種計(jì)算橫向同性分層介質(zhì)電磁散射問題的高階DTA級(jí)數(shù)算法。通過與現(xiàn)有方法對(duì)比發(fā)現(xiàn)，本文算法對(duì)于散射體為低阻和高阻時(shí)均具有較高的計(jì)算精度，而傳統(tǒng)Born近似以及DTA在目標(biāo)體為低阻時(shí)計(jì)算精度較低。引入壓縮映射的高階DTA級(jí)數(shù)收斂很快，只需要較少的階數(shù)就可以接近于嚴(yán)格方法的計(jì)算結(jié)果，而計(jì)算效率卻比嚴(yán)格方法有了較大提高。本文方法可用于橫向同性分層介質(zhì)中各種電磁響應(yīng)特征與規(guī)律的系統(tǒng)考察與分析，也可為進(jìn)行復(fù)雜介質(zhì)中電磁資料的快速三維反演成像提供一套有效正演算法。

［1］Chew W C，Jin J M，Michielsen E，et al.Fast and efficient algorithm in computational electromagnetics［C］.Norwood：Artech House，2001.

［2］Zhdanov M S.Geophysical electromagnetic theory and methods［C］.Amsterdam：Elsevier，2009.

［3］陳桂波，汪宏年，姚敬金，等.用積分方程法模擬各向異性地層三維電性散射體的電磁響應(yīng)［J］.地球物理學(xué)報(bào)，2009，52（8）：2174-2181.

［4］Chew W C.Waves and fields in inhomogeneous media［C］.New York：Van Nostrand Reinhold，1990.

［5］Habashy TM，Groom R W，Spies B R.Beyond the Born and Rytov approximations：anonlinearapproach to electromagnetic scattering［J］.Journalof Geophysical Research，1993，98（B2）：1759-1775.

［6］Zhdanov M S，F(xiàn)ang S.Quasi-linear approximation in 3D electromagnetic modeling［J］.Geophysics，1996，61（3）：646-665.

［7］Zhdanov M S，Dmitriev V I，F(xiàn)ang S，et al.Quasi-analytical approximations and series in electromagnetic modeling［J］.Geophysics，2000，65（6）：1746-1757.

［8］SongLP，LiuQH.Anewapproximationto three-dimensional electromagnetic scattering ［J］.IEEE Geoscience and Remote Sensing Letters，2005，2（2）：238-242.

［9］Song L P，Liu Q H.Fast three-dimensional electromagnetic nonlinear inversion in layered media with a novelscattering approximation［J］.Inverse Problems，2004，20（6）：171-194.

［10］陳桂波，汪宏年，姚敬金，等.水平層狀各向異性介質(zhì)中電磁場(chǎng)并矢格林函數(shù)的一種高效算法［J］.物理學(xué)報(bào)，2009，58（3）：1608-1618

［11］Pankratov O V，Kuvshinov A V，Avdeev D B.High-performance three-dimensional electromagnetic modeling using modified Neumann series：Anisotropic earth［J］.J Geomag Geoelectr，1997（49）：1541-1547.