各向異性分層介質電磁場并矢Green函數的混合算法

陳桂波，潘軼群，張燁

（1.長春理工大學 理學院，長春 130022；2.長春市農業機械研究院，長春 130062）

積分方程法是求解層狀介質三維電磁正演問題的一種常見方法，該方法的難點在于如何計算層狀介質的并矢Green函數（Dyadic Green’s Function，DGF），許多學者對此類問題作了深入研究，但多數只局限于各向同性介質模型［1-3］，而針對各向異性介質的研究則相對較少［4，5］。對于相對比較簡單的各向異性介質，如橫向同性各向異性介質，其層狀結構中的DGF一般也是由Sommerfeld積分組成，該積分的被積函數一般呈現出較強的振蕩特性與慢衰減特性，從而導致數值計算非常緩慢，這往往成為實際應用中的瓶頸問題。因此，國內外許多學者對Sommerfeld積分的快速求解問題作了研究［6-9］，并開發出相應的快速算法，歸結起來主要有：解析方法，改進的直接積分方法，濾波方法和復鏡像方法。盡管存在以上多種快速算法，但各種方法都有其獨特的適用條件，而且計算精度和計算效率對于介質參數、場源距離等因素的影響有很大程度的依賴關系，也就是說，沒有哪一種方法能夠滿足所有條件下Sommerfeld積分的快速計算。因此，開發一種能夠綜合運用各種方法的優勢并且根據參數條件自動選擇最優計算方案以達到最佳計算效果的混合算法，具有十分重要的理論與實際意義。

1 各向異性分層介質電磁場DGF

考慮圖1所示的N層各向異性介質（為簡單起見，本文只考慮單軸各向異性介質），各個地層的層界面位置為dn(n=2,3,…,N)、復介電常數和 磁 導 率 為 對 角 張 量diag(μnh,μnh,μnv) (n=1,2,…,N)。其中分別表示第n層橫向和縱向復介電常數，而μnh,μnv分別表示第n層橫向和縱向磁導率。復介電常數ε*n由公式 ε*np=εnp+(iωε0ρnp)-1(p=h,v)確定，其中 εnp為相對介電常數，ρnh,ρnv分別表示第n層橫向和縱向電阻率。

圖1 各向異性分層介質地電模型

對于空間中任意的激發源分布(J,M)，假定時諧因子為 eiωt，i=-1，其產生的電場和磁場滿足Maxwell方程組：

如果引入層狀介質中的DGF，E,H解可表示為［5］：

1.1 均勻各向異性介質DGF閉合解析計算

按Micalski（1997）所給的積分表達式，我們經過數學推導，得到了均勻各向異性介質中四種電磁場DGF的閉合解析解：

（1）電流源電場DGF解析解

其中α,β=x,y且α≠β，

（2）磁流源電場DGF解析解

其中

以上給出了電場DGF閉合解析解，相應的磁場DGF可由對偶原理得到［10］，限于篇幅，這里不給出其具體表達式。

1.2 DGF的Sommerfeld積分形式

根據Micalski（1997）的結果，層狀各向異性介質中的DGF由下列幾種形式的Sommerfeld積分組成：其中為傳輸線Green函數，文獻［5］中給出了8種傳輸線Green函數的詳細表達式，這里不再贅述。特別地，當ρ→0時，由Bessel函數的漸近性，式（22）-（24）的Sommerfeld積分有的變為0，有的變為普通的半無窮積分

2 各向異性分層介質DGF的混合算法

2.1 單一方法計算Sommerfeld積分的不足

（1）解析方法

解析方法一般利用數學推導給出Sommerfeld積分的解析表達式，由于其計算速度較快并且易于對物理現象的直觀解釋和理解，解析方法在均勻半空間以及一些特殊模型問題的計算與分析中占有重要地位。Sommerfeld、Foster和Wait等人在這方面做過深入的研究［6，12，13］，尤其是Sommerfeld經過復雜數學推導得到的Sommerfeld恒等式，至今仍在許多領域被廣泛使用。解析方法的不足是應用范圍十分狹窄，僅限于均勻介質、半空間以及一些較特殊的模型。

（2）基于直接積分的方法

直接積分法利用計算機程序，通過不斷的外延積分進行Sommerfeld積分的高精度計算，其受參數影響較小，計算精度較高，但由于Sommerfeld積分被積函數的性質，該積分在積分路徑上收斂緩慢。雖然一些學者對直接積分法作了改進，引入非線性變換手段來加速積分收斂［14，15］，但往往仍難以滿足需要調用大量Sommerfeld積分的三維電磁模擬的計算。

（3）濾波方法

數字濾波技術誕生于上實際70年代初，由Ghosh首先使用［16］，之后Johnsen，Koefoed，Anderson等人對該方法作了改進和發展［17-19］。該方法通過輸入函數（核函數）與濾波因子在給定長度下褶積求和從而快速得到積分近似值。由于避免了Bessel函數的計算和數值積分，使Sommerfeld積分的計算速度得到大大提高。但問題在于該方法對參數范圍較為敏感，尤其是在橫向距離ρ較小或較大的情況，由于核函數與Bessel函數相比變化劇烈，此時數字濾波方法的精度較差。

（4）鏡像方法

鏡像方法源于幾何光學理論，在地電領域最早由Wait在計算半空間磁偶極子的輻射時引入［20］。而后，Lindell曾利用連續復鏡像法計算Green函數［21］。近年來，離散復鏡像法［22，23］在電磁學領域的應用越來越廣泛，該方法利用現有的Bessel積分恒等式結合復指數序列逼近得到Sommerfeld積分的閉合解。其優點是計算速度很快，但適用條件較高，往往要求工作頻率限定在某一范圍內且源點和場點位于同一介質層中。

由前面的分析和討論可知，Sommerfeld積分的各種計算方法都有其特殊的適用條件，而沒有哪一種方法是一勞永逸的。因此，本文擬從參數條件出發，開發一種能夠綜合運用各種方法的優勢并根據參數范圍選擇最優計算方案的高效混合算法。這里需要說明的是，本文所研究的算法主要是針對其在積分方程中的應用，而在利用積分方程法求解三維電磁問題時，一般計算量較大的是源點和場點位于同一地層或相鄰地層時的DGF。由電磁理論，這時可將DGF分為對應均勻無限大介質中的直射項和反映地層邊界條件的反射項。為了提高計算精度與效率，直射項采用前面推導出的DGF的解析式計算；而對于反射項，其所包含的Sommerfeld積分往往沒有解析解，一般需要用數值方法來計算。

橫向距離ρ趨近于零時DGF的計算是在三維電磁正演模擬中經常遇到的一個問題，但在以往相關文獻中卻很少被提及。本文在前面式（21）和式（25）已經給出了當ρ→0時直射項的解析解和反射項中Sommerfeld積分的表達式及計算方案，在精度允許范圍內，一般ρ＜10-5即可以使用ρ→0時的結果來近似。當10-5≤ρ≤10-3或ρ≥103時，由于數字濾波方法對ρ很小或很大時特別敏感，若采用數字濾波方法可能會得到較大誤差甚至錯誤結果，而基于直接積分的方法則受參數影響較小并且可以得到較高的精度，這里我們采用一種新近開發的使用高階窗函數改良Sommerfeld積分被積函數的性質并結合非線性變換手段加速積分收斂的高效直接積分算法［20］。經過研究多種加快級數收斂速度的變換方法我們發現，使用連分式展開具有更好的收斂速度，圖2為相同條件下分別使用連分式展開和歐拉變換［14］計算了均勻介質中磁流源磁場DGF的三個主分量，在精度相當的情況下兩種方法所用計算時間分別為2.81、4.15秒，因此我們以高階窗函數結合連分式展開法（HWCFM）作為改進的直接積分方法。當10-3＜ρ＜103時，由于我們需要計算的四種場型DGF的Sommerfeld積分的被積函數均不存在性質特別差的情況，因此采用數字濾波方法計算Sommerfeld積分一般可以得到較高的計算精度與效率。由于離散復鏡像法的適用條件過于苛刻，因此我們在混合算法中暫不引入離散復鏡像法。

圖2 使用連分式展開與歐拉變換方法計算GMH三個主分量對比

歸結上面的分析討論，我們得到一種計算各向異性分層介質DGF的單點混合算法（SPHM）：并矢Green函數

2.2 DGF的混合算法

在利用DGF研究電磁輻射與散射問題時，往往需要計算沿某一方向或不同方向上的一系列場點的函數值，此時若采用直接的單點算法，雖然其中應用了一些加快計算速度的手段，但往往還是難以滿足實際工作需要。對于一些特殊而用又常用的模型，我們充分考慮構成DGF的Sommerfeld積分被積函數的某些性質，得到一些計算速度更快且不失精度的全空間DGF的算法。

首先考慮對于固定的(z,z')，需要計算某一范圍內的橫向距離ρ對應的DGF值。當ρ≤10-3或ρ≥103時，采用前面介紹的SPHM計算。而對于實際中經常遇到的位于(10-3,103)范圍內的多個ρ對應的Sommerfeld積分的計算，我們引入Anderson（1982）提出的延遲褶積算法，選擇一個參數序列ρN+1-j=e-Δ(j-1)ρmax，j=1,2,…,N ，ρ1=ρmin，其中Δ為對數坐標下的濾波系數采樣間隔，ρmax,ρmin分別為ρ的最大值和最小值。當 j=1時，計算全部濾波系數下的離散褶積，并保存所有的輸入函數值。之后對于參數序列的每一次延遲即j=2,3,…,N，根據所選特殊參數序列的性質，輸入函數的指標j分別相對于前一次時增加一個采樣間隔，這時使用已保存的輸入函數值和濾波系數做褶積求得 j=2,3,…,N時的Sommerfeld積分，最后利用插值方法可以快速得到[ρmin,ρmax]范圍內所有的ρ對應的DGF，以上延遲褶積算法的詳細過程參見文獻［8］。本文將上述計算任意參數范圍內多點DGF的算法簡稱為MPHM。MPHM大大提高了全空間上DGF的計算效率，圖3（a）為使用MPHM與SPHM計算了三層介質電流源DGF分量，兩種方法計算圖中50個節點所用時間分別為0.61和0.96秒，可見引入MPHM使計算速度得到了較大提高。這里需要說明的是，本文SPHM和MPHM中使用的數字濾波方法都是變長度的自適應數字濾波方法［18］。

3 DGF混合算法應用

圖3 應用MPHM與SPHM法計算三層介質中的DGF

圖4 本文方法與解析方法計算半空間介質表面電偶極子的輻射電場（地層的平均電阻率為100Ω?m，f=1000Hz）

考察源偶極子位于地表以及地中時的電磁輻射問題，本文所選用模型實例均假定為非磁性的。按照Wait的方法［13］，得到位于各向異性介質表面的水平電偶極子在地表產生的電場水平分量的解析式。令沿x軸方向的單位大小的電偶極子位于地表(x',y')處，則在地表(x,y)處產生的電場水平分量為其中和λe為半空間均勻介質的橫向復介電常數和電各向異性系數，波數kh,ke的計算公式前面已經給出。圖4為分別利用本文DGF方法和解析方法計算不同各向異性系數條件下均勻半空間介質表面的電偶極子輻射場，可以看出，兩種方法計算的曲線基本一致，并且介質的各向異性對電場水平分量的實部幾乎沒有影響，但對虛部的影響卻較明顯。圖5是利用本文算法計算半空間介質中的電偶極子輻射場，與文獻［21］的結果吻合得較好。

圖5 本文方法與Xiong（1989）計算半空間介質中的偶極子輻射

本文中所有數值計算均是在同一臺主頻3.0GHz計算機上完成的。

4 結論

利用傳輸線Green函數和復變積分理論得到均勻介質中四種場型DGF的閉合解析解，并特別給出了場源橫向距離趨于零時的解。對于層狀介質中的Sommerfeld積分，對比分析各種算法的優缺點與適用條件，得到一種快速、高精度且不受參數范圍限制的DGF的混合算法，并引入延遲褶積方法進一步加快特殊幾何模型下全空間DGF的計算。研究DGF在計算層狀各向異性介質偶極子輻射場以及積分方程中的應用，并得到了一些具有重要參考價值的數值結果。本文方法可直接用于積分方程法求解三維電磁正演問題，另外也可作為研究電導率與磁導率各向異性對地下資源勘探的綜合影響的一種手段。

［1］Michiels B，Fostier J，Bogaert I，et al.Full-wave simulations of electromagnetic scattering problems with billions of unknowns［J］.IEEE Transactions on Antennas and Propagation，Publication，2015，63（2）：796-799.

［2］Wiedenmann O，Eibert T F.A domain decomposition method for boundary integral equations using a transmission condition based on the near-zone couplings［J］.IEEE Transactions on Antennas and Propagation，2014，62（8）：4105-4114.

［3］Markkanen J，Yla Oijala P.Discretization of electric current volume integral equation with piecewise linear basis functions［J］.IEEE Transactions on Antennas and Propagation，2014，62（9）：4877-4880.

［4］Xiong Zonghou.Induced-polarization and electromagnetic modeling of a three-dimensional body buried in a two-layer anisotropic earth［J］.Geophysics，1986，51（12）：2235-2246.

［5］MicalskiK A，Mosig JR.Multilayered media Green’s functions in integral equation formulations［J］.IEEE transactions on antennas and propagation，1997，45（3）：508-519.

［6］Sommerfeld A. Partialdifferentialequations in physics［A］.New York：Academic，1949.

［7］Cui T J.Fast evaluation of Sommerfeld integrals for EM wave scattering and radiation by three-dimensional buried objects［J］.IEEE Transactions on Geoscience and Remote Sensing，1999，37（3）：887-900.

［8］Anderson W L.Fast Hankel Transforms using related and lagged convolutions［J］.ACM Transactions on Mathematical Software，1982，8（4）：344-368.

［9］Aksun M I.A robust approach for the derivation of closed form green functions［J］.IEEE Trans.Microwave Theory Tech，1996，44（5）：651-658.

［10］Tai C T.Dyadic Green’s Functions in Electromagnetic Theory［A］.2nd ed.IEEE Press，New York，1993.

［11］Piessens R，Ede Doncker-Kapenga.QUADPACK：A SubroutinePackageforAutomaticIntegration［A］.New York：Springer-Verlag，1983.

［12］Campbell G A，Foster R M.Fourier Integrals for Practical Applications［A］.New York，Bell Telephone System，1931.

［13］WaitJR.Geo-electromagneticsm［M］.Academic Press，1982.

［14］Cornille P.Computation of Hankel transforms［J］.SIAM Rev，1972，14：278-285.

［15］Micalski K A.Extrapolation methods for Sommerfeld Integral Tails［J］.IEEE transactions on antennas and propagation，1998，46（10）：1405-1418.

［16］Johansen H K，Sorensen K.Fast Hankel transforms［J］.Geophysical Prospecting，1979，27（4）：876-901.

［17］KoefoedO，GhoahD P.Computationoftype curves for electromagnetic depth sounding with a horizontal transmitting coil by means of a digital linear filter［J］.Geophysical Prospecting.1972，20（2）：406-420.

［18］Anderson W L.Numerical integration of related Hankel transforms of orders 0 and 1 by adaptive digital filtering［J］.Geophysics，1979，44（7）：1287-1305.

［19］Wait J R.Image theory of a quasistatic magnetic dipole over a dissipative half-space［J］.Electron.Lett，1969，5（13）：281-282.

［20］陳桂波.各向異性分層介質中電磁場并矢Green函數的一種高效算法［J］.物理學報，2009，58（3）：1608-1618.

［21］Zonghou Xiong.Electromagnetic fields of electric dipoles embedded in a stratified anisotropic earth［J］.Geophysics，1989，54（12）：1643-1646.