陳桂波，張燁，潘軼群

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

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

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

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

1 DTA近似基本原理

圖1 三維模型示意圖

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

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

這里Ebξ為入射電場Eb的ξ分量，而EBξ為矢量

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

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

2 高階DTA級數(shù)

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

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

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

其中

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

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

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

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

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

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

圖4 散射體為低阻時的散射場

圖5 目標體為高阻時的散射場

4 結(jié)論

本文提出了一種計算橫向同性分層介質(zhì)電磁散射問題的高階DTA級數(shù)算法。通過與現(xiàn)有方法對比發(fā)現(xiàn)，本文算法對于散射體為低阻和高阻時均具有較高的計算精度，而傳統(tǒng)Born近似以及DTA在目標體為低阻時計算精度較低。引入壓縮映射的高階DTA級數(shù)收斂很快，只需要較少的階數(shù)就可以接近于嚴格方法的計算結(jié)果，而計算效率卻比嚴格方法有了較大提高。本文方法可用于橫向同性分層介質(zhì)中各種電磁響應特征與規(guī)律的系統(tǒng)考察與分析，也可為進行復雜介質(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］陳桂波，汪宏年，姚敬金，等.用積分方程法模擬各向異性地層三維電性散射體的電磁響應［J］.地球物理學報，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ì)中電磁場并矢格林函數(shù)的一種高效算法［J］.物理學報，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.