陳桂波，潘軼群，張燁

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

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

1 各向異性分層介質(zhì)電磁場(chǎng)DGF

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

圖1 各向異性分層介質(zhì)地電模型

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

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

1.1 均勻各向異性介質(zhì)DGF閉合解析計(jì)算

按Micalski（1997）所給的積分表達(dá)式，我們經(jīng)過(guò)數(shù)學(xué)推導(dǎo)，得到了均勻各向異性介質(zhì)中四種電磁場(chǎng)DGF的閉合解析解：

（1）電流源電場(chǎng)DGF解析解

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

（2）磁流源電場(chǎng)DGF解析解

其中

以上給出了電場(chǎng)DGF閉合解析解，相應(yīng)的磁場(chǎng)DGF可由對(duì)偶原理得到［10］，限于篇幅，這里不給出其具體表達(dá)式。

1.2 DGF的Sommerfeld積分形式

根據(jù)Micalski（1997）的結(jié)果，層狀各向異性介質(zhì)中的DGF由下列幾種形式的Sommerfeld積分組成：其中為傳輸線Green函數(shù)，文獻(xiàn)［5］中給出了8種傳輸線Green函數(shù)的詳細(xì)表達(dá)式，這里不再贅述。特別地，當(dāng)ρ→0時(shí)，由Bessel函數(shù)的漸近性，式（22）-（24）的Sommerfeld積分有的變?yōu)?，有的變?yōu)槠胀ǖ陌霟o(wú)窮積分

2 各向異性分層介質(zhì)DGF的混合算法

2.1 單一方法計(jì)算Sommerfeld積分的不足

（1）解析方法

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

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

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

（3）濾波方法

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

（4）鏡像方法

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

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

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

圖2 使用連分式展開(kāi)與歐拉變換方法計(jì)算GMH三個(gè)主分量對(duì)比

歸結(jié)上面的分析討論，我們得到一種計(jì)算各向異性分層介質(zhì)DGF的單點(diǎn)混合算法（SPHM）：并矢Green函數(shù)

2.2 DGF的混合算法

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

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

3 DGF混合算法應(yīng)用

圖3 應(yīng)用MPHM與SPHM法計(jì)算三層介質(zhì)中的DGF

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

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

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

本文中所有數(shù)值計(jì)算均是在同一臺(tái)主頻3.0GHz計(jì)算機(jī)上完成的。

4 結(jié)論

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

［1］Michiels B，F(xiàn)ostier 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，F(xiàn)oster 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］陳桂波.各向異性分層介質(zhì)中電磁場(chǎng)并矢Green函數(shù)的一種高效算法［J］.物理學(xué)報(bào)，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.