殷長春，盧永超，劉云鶴，張 博，齊彥福，蔡 晶

吉林大學(xué)地球探測科學(xué)與技術(shù)學(xué)院，長春 130026

0 引言

隨著我國經(jīng)濟(jì)的快速增長，礦產(chǎn)資源的消耗逐漸增多，地質(zhì)條件有利地區(qū)的礦產(chǎn)資源已近枯竭，目標(biāo)轉(zhuǎn)向廣大無人區(qū)。由于這些地區(qū)施工困難，常規(guī)的地面電磁法很難滿足勘探需求，航空電磁法勘探因此成為重要的勘查技術(shù)手段[1]。三維積分方程正演是航空電磁數(shù)值模擬的重要技術(shù)，它通過對(duì)異常區(qū)域進(jìn)行剖分，建立異常區(qū)域內(nèi)電場滿足的第二類Fredholm積分方程，通過直接求解或利用迭代技術(shù)求解離散后的矩陣方程得到空間電磁場。在小規(guī)模異常體模擬中，系數(shù)矩陣占用存儲(chǔ)小，求解速度快。然而，對(duì)于大型異常體，由于形成的線性方程組系數(shù)矩陣的密實(shí)性，計(jì)算耗時(shí)多、存儲(chǔ)量大，限制了積分方程算法的應(yīng)用。

為了避免大型線性方程組的求解，提高計(jì)算效率，學(xué)者提出許多近似方法。典型的方法包括擴(kuò)展Born近似[2]、準(zhǔn)線性近似(quasi-linear or QL)[3]、準(zhǔn)線性級(jí)數(shù)(QLSE)[4]、準(zhǔn)解析(QA)近似和準(zhǔn)解析級(jí)數(shù)(QASE)[5]等。然而，上述近似方法只有在滿足其特定的使用條件下才能獲得理想的計(jì)算結(jié)果。通過對(duì)QL技術(shù)的改進(jìn)，Zhdanov等[6]實(shí)現(xiàn)了局部準(zhǔn)線性近似(LQL)的頻率域電磁正反演，并在對(duì)加拿大拉布拉多沃伊西灣鎳-銅-鈷礦航空電磁數(shù)據(jù)反演中取得了良好的效果。然而，該算法需要通過求解最小值問題來獲得電反射系數(shù)，計(jì)算效率低。劉永亮等[7]提出了快速準(zhǔn)線性近似方法，并將其應(yīng)用于三維頻譜激電反演研究，取得了很好的效果。針對(duì)多源問題的海洋可控源電磁正演模擬，Ueda等[8]結(jié)合積分方程(IE)和QL技術(shù)，提出了基于多重網(wǎng)格的準(zhǔn)線性(MGQL)近似算法。MGQL算法的關(guān)鍵是假設(shè)背景場與異常場之間存在線性關(guān)系?；贛GQL技術(shù)，人們首先在粗網(wǎng)格上求解積分方程得到總場和異常場，進(jìn)而求解粗網(wǎng)格上的電反射系數(shù)，并利用該系數(shù)在細(xì)網(wǎng)格上做線性插值得到離散后的電反射系數(shù)。在細(xì)網(wǎng)格上，利用電反射系數(shù)與背景場求解異常場，最后在細(xì)網(wǎng)格上獲得接收點(diǎn)處的場值。與全積分方程法相比，對(duì)于大型異常體和網(wǎng)格剖分較多的情況，MGQL通過粗化網(wǎng)格可以大量節(jié)省計(jì)算時(shí)間。本文將該技術(shù)引入到多源航空電磁三維正演模擬中，由于反射系數(shù)計(jì)算在粗網(wǎng)格上進(jìn)行，而細(xì)網(wǎng)格上的場值通過插值獲得，因此可極大提高多源航空電磁的計(jì)算效率。

本文在MGQL算法的基礎(chǔ)上，利用計(jì)算電磁學(xué)中的處理技術(shù)[9]優(yōu)化粗網(wǎng)格上積分方程的計(jì)算過程，將發(fā)射源首次發(fā)射時(shí)的格林函數(shù)存儲(chǔ)為Toeplitz矩陣形式，避免了不同源格林函數(shù)的重復(fù)計(jì)算，進(jìn)而利用快速傅里葉變換(FFT)技術(shù)實(shí)現(xiàn)矩陣與向量的乘積[10]，提升MGQL的計(jì)算效率。MGQL算法網(wǎng)格粗化與計(jì)算精度有關(guān)，因此在精度允許的范圍內(nèi)，調(diào)節(jié)粗網(wǎng)格大小可有效降低存儲(chǔ)量、加快計(jì)算速度。最后，通過將本文算法結(jié)果和已發(fā)表的典型模型計(jì)算結(jié)果對(duì)比驗(yàn)證精度，并通過與開源代碼MarcoAir對(duì)比驗(yàn)證本文MGQL算法的有效性。

1 MGQL正演理論

根據(jù)積分方程理論，在水平層狀各向同性介質(zhì)中，對(duì)麥克斯韋方程組進(jìn)行變換，可得電場E與磁場H的如下求解方程[11-12]：

[Δσ(r)E(r)]dV+Eb(r')，

(1)

[Δσ(r)E(r)]dV+Hb(r')。

(2)

E(r)=Eb(r)+Ea(r)。

(3)

根據(jù)Zhdanov等[3]提出的準(zhǔn)線性近似理論，假設(shè)背景場和異常場之間存在線性關(guān)系，即

(4)

(5)

(6)

式中，EB(rj)表示異常場的Born近似[13]，其表達(dá)式為

(7)

利用最小化求解方法，由式(6)可得電反射系數(shù)的求解表達(dá)式為

(8)

Ea(rc)=E(rc)-Eb(rc)。

(9)

(10)

因此，式(4)可以寫成獨(dú)立的標(biāo)量式，即

；l=x,y,z。

(11)

必須指出的是，依據(jù)公式(4)中的假設(shè)，反射系數(shù)計(jì)算雖然簡單，但是如果背景場的某一個(gè)分量為零，將導(dǎo)致反射系數(shù)無法獲取。為此，本文根據(jù)Gao等[15]提出的針對(duì)各向異性介質(zhì)三維電磁感應(yīng)測井的模擬算法，假設(shè)異常場和層狀半空間背景場絕對(duì)值存在如下線性關(guān)系：

Ea(r)≈λ(r)·|Eb(r)|。

(12)

式中，向量λ(r)=(λx,λy,λz)。式(12)表明，在粗網(wǎng)格上，由于|Eb(rc)|≠0，則可以得到電反射系數(shù)的計(jì)算公式為

；l=x,y,z。

(13)

利用式(13)計(jì)算λ(rc)之后，通過線性插值計(jì)算細(xì)網(wǎng)格上的λ(rf)，則由式(12)可得單元rf中心處的異常場為

Ea(rf)≈λ(rf)·|Eb(rf)|。

(14)

而總場由式(3)給出，最后利用式(2)可計(jì)算觀測點(diǎn)處的電磁場[16-19]。

2 模型試驗(yàn)和結(jié)果分析

2.1 精度驗(yàn)證

針對(duì)Newman等[20]給出的三維航空電磁模型，我們將本文算法的結(jié)果與其給出的全積分方程結(jié)果進(jìn)行對(duì)比，以驗(yàn)證MGQL算法的精確性。模型如圖1所示，粗網(wǎng)格剖分?jǐn)?shù)為2×10×10(x×y×z，下同)個(gè)，細(xì)網(wǎng)格剖分?jǐn)?shù)為4×10×20個(gè)。 測線沿著x方向，發(fā)射頻率為900 Hz，系統(tǒng)飛行高度20 m。三維異常體的幾何參數(shù)和埋深如圖1所示。

圖2a、b給出利用本文算法和Newman文中 IE計(jì)算的電磁響應(yīng)對(duì)比結(jié)果，可見MGQL的磁場Hx、Hz實(shí)虛分量曲線與Newman的結(jié)果吻合很好。圖2c、d給出的相對(duì)誤差均保持在5%以內(nèi)，說明本文算法具有很高的精度。

據(jù)文獻(xiàn)[20]。圖1 三維正演模型Fig.1 Models for 3D modeling

2.2 單個(gè)異常體模型

航空電磁勘探在地形條件比較復(fù)雜的區(qū)域應(yīng)用較多，目標(biāo)體上方的覆蓋層會(huì)對(duì)異常響應(yīng)產(chǎn)生影響[21-22]。因此，本節(jié)采用的模型與圖1類似，只是在地表加一層電阻率為10 Ω·m、厚度為10 m覆蓋層，異常體為30 m×180 m×70 m, 電阻率為1 Ω·m，埋深40 m，粗網(wǎng)格剖分網(wǎng)格數(shù)為2×10×10個(gè)，細(xì)網(wǎng)格剖分網(wǎng)格數(shù)為4×10×20個(gè)。測線沿x方向，飛行高度30 m，頻率900 Hz。

圖2 MGQL算法與Newman IE(1995)正演結(jié)果對(duì)比Fig.2 Comparison of MGQL results from this paper with those of Newman (1995)

圖3給出含覆蓋層與不含覆蓋層的單個(gè)異常體航空電磁響應(yīng)對(duì)比。由3圖可見：覆蓋層不影響異常響應(yīng)曲線形態(tài)，但是異常響應(yīng)幅值發(fā)生明顯變化，含覆蓋層異常體電磁響應(yīng)明顯大于不含覆蓋層的電磁響應(yīng)，說明良導(dǎo)覆蓋層增強(qiáng)了異常體的電磁響應(yīng)(背景值得到提升)。

2.3 多個(gè)異常體模型

在三維航空電磁正演計(jì)算中，多個(gè)異常體之間存在電磁耦合。利用Zhdanov等[14]提出的不均勻背景電導(dǎo)率法(IBC)，當(dāng)存在多個(gè)異常體時(shí)，對(duì)每個(gè)異常體使用積分方程進(jìn)行計(jì)算，然后通過耦合迭代求解電磁響應(yīng)，但這種方法計(jì)算速度較慢。本文采用MGQL算法，由于積分是在粗糙網(wǎng)格上進(jìn)行的，耦合迭代速度極大地提升。這種方法對(duì)網(wǎng)格數(shù)巨大的模型求解優(yōu)勢(shì)更為明顯。為驗(yàn)證MGQL算法對(duì)多個(gè)異常體模擬的有效性，我們?cè)O(shè)計(jì)如圖4所示的模型，頻率為900 Hz，粗網(wǎng)格剖分?jǐn)?shù)為4×10×10個(gè)，細(xì)網(wǎng)格剖分?jǐn)?shù)為8×10×20個(gè)，測線沿著x方向，飛行高度30 m。

為了方便對(duì)比，本文分別計(jì)算存在單個(gè)異常體和多個(gè)異常體時(shí)電磁響應(yīng)。圖5a、b給出存在單個(gè)異常體時(shí)Hx和Hz響應(yīng)，而圖5c、d給出存在兩個(gè)異常體時(shí)Hx和Hz響應(yīng)曲線。由圖5可以看出，由于兩個(gè)異常體之間存在電磁耦合，單個(gè)異常體電磁響應(yīng)的疊加明顯大于兩個(gè)異常體組合的電磁響應(yīng)；說明異常體之間的二次耦合不可忽視。

3 計(jì)算效率分析

本文提出的MGQL算法的最大優(yōu)點(diǎn)是在保證計(jì)算精度的前提下達(dá)到快速計(jì)算的目標(biāo)。為了分析MGQL算法的計(jì)算效率，本文將其與開源軟件MarcoAir (Version 2.3.1)及全積分方程算法結(jié)果作對(duì)比。本文設(shè)計(jì)與圖1類似的模型，異常體大小為30 m×200 m×100 m，電阻率1 Ω·m，埋深30 m，飛行高度30 m, 頻率900 Hz，測線沿x方向。通常，全矩陣存儲(chǔ)積分方程法需要占用大量的內(nèi)存空間，用于存儲(chǔ)格林系數(shù)矩陣。受內(nèi)存限制，全矩陣存儲(chǔ)積分方程法計(jì)算網(wǎng)格數(shù)有限。使用MGQL技術(shù)能夠減少網(wǎng)格數(shù)、節(jié)約存儲(chǔ)，可實(shí)現(xiàn)多網(wǎng)格數(shù)剖分的正演模擬。本文的計(jì)算在Intel(R) Xeon(R) CPU E5-2667 v3@3.20 Hz處理器上進(jìn)行，內(nèi)存128 GB，操作系統(tǒng)為64位Win7系統(tǒng)。

MarcoAir屬于早期研發(fā)的積分方程算法，它利用格林函數(shù)的對(duì)稱性和塊迭代算法提升積分方程的計(jì)算速度[22]。當(dāng)異常體剖分單元較少時(shí)，直接解法優(yōu)勢(shì)較大；當(dāng)剖分單元較多時(shí)，迭代法計(jì)算速度較快。我們對(duì)比時(shí)選用MarcoAir迭代法求解線性方程組。本文MGQL算法中將格林函數(shù)存儲(chǔ)為Toeplitz矩陣形式，同時(shí)利用FFT實(shí)現(xiàn)矩陣與向量乘積，并采用穩(wěn)定雙共軛梯度法(BICGSTAB)求解線性方程組[23-24]，以進(jìn)一步加快計(jì)算速度。假設(shè)航空電磁系統(tǒng)發(fā)射源和接收機(jī)共移動(dòng)60次。表1給出在相同計(jì)算精度條件下，本文算法計(jì)算時(shí)間與MarcoAir的對(duì)比結(jié)果。由表1可見：當(dāng)網(wǎng)格數(shù)較少時(shí)，MarcoAir求解速度快于MGQL算法；隨著模型剖分網(wǎng)格數(shù)增加，MGQL算法計(jì)算速度遠(yuǎn)遠(yuǎn)超過MarcoAir。例如，當(dāng)網(wǎng)格數(shù)增加到10×20×20時(shí)，MGQL的計(jì)算速度是MarcoAir近54倍。另外，對(duì)于計(jì)算設(shè)備的有限存儲(chǔ)空間和時(shí)間，MarcoAir可以計(jì)算的網(wǎng)格數(shù)不超過20×40×40個(gè)，而對(duì)于MGQL算法來說，可以計(jì)算的網(wǎng)格數(shù)為50×100×100。由于MGQL在粗網(wǎng)格上求解異常場，計(jì)算量極大減少，使得該算法在計(jì)算效率上擁有很大的優(yōu)勢(shì)。

圖3 不含覆蓋層和含覆蓋層單個(gè)異常體航空電磁響應(yīng)Fig.3 AEM responses for an anomalous body in a homogeneous earth with and without overburden

圖4 多個(gè)異常體三維正演模型Fig.4 Multiple anomalous bodies for 3D modeling

a，b.單個(gè)異常體存在時(shí)；c，d.多個(gè)異常體存在時(shí)。圖5 單個(gè)異常體和多個(gè)異常體的航空電磁響應(yīng)對(duì)比Fig.5 Comparison of AEM responses for single anomalous body and multiple anomalous bodies

表1MGQL和MarcoAir不同網(wǎng)格剖分的計(jì)算時(shí)間對(duì)比

Table1ComputationtimeforMGQLandMarcoAirfordifferentgrids

Nx×Ny×NztMGQL/mintMarcoAir/min5×5×50.020.015×10×100.110.5610×20×201.3775.6120×40×4022.36****30×60×60112.38****50×100×1001318.17****100×100×100********

注：N為網(wǎng)格數(shù)；t為計(jì)算時(shí)間；*|*|*|*表示計(jì)算時(shí)間大于24 h。

對(duì)于表1設(shè)計(jì)的模型，我們進(jìn)一步在MGQL與全積分方程法(全矩陣存儲(chǔ)和直接矩陣向量乘積算法)之間進(jìn)行計(jì)算效率對(duì)比，結(jié)果如圖6所示。對(duì)比結(jié)果表明，MGQL算法明顯具有計(jì)算速度快、占用內(nèi)存少的優(yōu)勢(shì)。隨著剖分網(wǎng)格增加，全積分方程計(jì)算時(shí)間和占用內(nèi)存呈指數(shù)增加；但對(duì)MGQL算法，隨網(wǎng)格數(shù)增加，其計(jì)算時(shí)間和占用內(nèi)存增長緩慢，凸顯出MGQL算法的高效性。

圖6 MGQL算法與全積分方程法計(jì)算時(shí)間和占用內(nèi)存對(duì)比Fig.6 Comparison of MGQL with classic IE from calculation time and memory requirement

下面說明MGQL算法如何通過粗細(xì)網(wǎng)格結(jié)合，實(shí)現(xiàn)航空電磁響應(yīng)高效三維數(shù)值模擬。采用的模型與表1相同，航空電磁系統(tǒng)發(fā)射源和接收機(jī)共移動(dòng)60次，異常體離散為2×10×10個(gè)粗網(wǎng)格和4×10×20個(gè)細(xì)網(wǎng)格。我們分別采用如下3種方法：1)在粗網(wǎng)格上應(yīng)用積分方程法(BICGSTAB-FFT)；2)在精細(xì)網(wǎng)格上應(yīng)用積分方程法；3)在粗細(xì)多重網(wǎng)格上應(yīng)用MGQL算法。圖7是3種方法的計(jì)算的Hz響應(yīng)對(duì)比。由7圖可以看出：方法1)的結(jié)果與方法3)的計(jì)算結(jié)果差距較大，而方法2)和3)的響應(yīng)曲線吻合度高。3種方法的計(jì)算時(shí)間分別為13.12、69.61和16.06 s。很顯然，本文算法在保證計(jì)算精度的同時(shí)，大大提升了航空電磁數(shù)據(jù)正演模擬的計(jì)算效率。

圖7 3種方法計(jì)算的Hz響應(yīng)對(duì)比Fig.7 Comparison of Hz responses calculated by three methods

4 結(jié)論

本文提出的MGQL算法有效地解決了積分方程計(jì)算時(shí)占用內(nèi)存大、計(jì)算速度慢的問題。通過移動(dòng)場源航空電磁三維模型試驗(yàn)，并與已發(fā)表的算法結(jié)果對(duì)比，表明響應(yīng)誤差均小于5%。另外本文算法模擬多個(gè)異常體也取得很好的結(jié)果，說明MGQL算法在滿足精度要求的條件下能夠求解大型復(fù)雜異常體模型，是解決多源航空電磁問題的有效算法。

在相同模型保證計(jì)算精度的前提下，本文算法相比于MarcoAir算法和全積分方程法，網(wǎng)格數(shù)越多，相對(duì)計(jì)算時(shí)間越少、占用存儲(chǔ)越小，說明利用矩陣Toeplitz性質(zhì)和BICGSTAB-FFT有效地提升了計(jì)算效率。探究精算精度與網(wǎng)格的關(guān)系發(fā)現(xiàn)，采用組合網(wǎng)格可以使用最小的內(nèi)存和計(jì)算時(shí)間獲得高的計(jì)算精度。因此，未來我們將嘗試?yán)肕GQL算法進(jìn)行多源航空電磁數(shù)據(jù)的快速成像和反演。

致謝: 感謝開源軟件MarcoAir的作者，感謝吉林大學(xué) “千人計(jì)劃”電磁研究團(tuán)隊(duì)全體成員。

[1] 殷長春，張博，劉云鶴，等.航空電磁勘查技術(shù)發(fā)展現(xiàn)狀及展望[J].地球物理學(xué)報(bào)，2015，58(8): 2637-2653.

Yin Changchun，Zhang Bo，Liu Yunhe，et al. Review on Airborne EM Technology and Developments [J]. Chinese Journal Geophysics，2015，58(8): 2637-2653.

[2]Habashy T M，Groom R W，Spies B R. Beyond the Born and Rytov Approximations: A Nonlinear Approach to Electromagnetic Scattering [J]. J Geophys Res，1993，98(B2): 1759-1775.

[3] Zhdanov M S，F(xiàn)ang S. Quasi-Linear Approximation in 3-D Electromagnetic Modeling [J]. Geophysics，1996，61(3): 646-665.

[4] Zhdanov M S，F(xiàn)ang S. Quasi-Linear Series in Three-Dimensional Electromagnetic Modeling [J]. Radio Science，1997，32(6): 2167-2188.

[5]Zhdanov M S，Dmitriev V I，F(xiàn)ang S，et al. Quasi-Analytical Approximations and Series in 3D Electromagnetic Modeling [J]. Geophysics，2000，65(6): 1746-1757.

[6] Zhdanov M S，Tartaras E. Three-Dimensional Inver-sion of Multitransmitter Electromagnetic Data Based on the Localized Quasi-Linear Approximation [J]. Geophys J Int，2002，148 (3): 506-519.

[7] 劉永亮，李桐林，胡英才，等. 快速擬線性近似方法及三維頻譜激電反演研究[J]. 地球物理學(xué)報(bào)，2015，58(12): 4709-4717.

Liu Yongliang，Li Tonglin，Hu Yingcai，et al. Fast Quasi-Linear Approximation and the Three-Dimensional Spectrum of Induced Polarization Inversion Study [J]. Chinese Journal Geophysics，2015，58(12): 4709-4717.

[8]Ueda T，Zhdanov M S. Fast Numerical Modeling of Multitransmitter Electromagnetic Data Using Multigrid Quasi-Linear Approximation [J]. IEEE Transactions on Geoscience and Remote Sensing，2006，44(6): 1428-1434.

[9] 周后型.大型電磁問題快速算法的研究[D].南京：東南大學(xué)，2002.

Zhou Houxing. Inverstigations on Fast Algorithms for Electrically Large Eletromagnetic Problems [D]. Nanjing: Southeast University, 2002.

[10] Zhang Z Q，Liu Q H. Three-Dimensional Weak-Form Conjugate- and Biconjugate-Gradient FFT Methods for Volume Integral Equations [J]. Microwave and Optical Technology Letters，2001，29(5): 350-356.

[11] Hohmann G W. Three-Dimesional Induced Polariza-tion and Electromagnetic Modeling [J]. Geophysics，1975，40(2): 309-324.

[12]Wannamaker P E，Hohmann G W，Sanfilipo W A. Electromagnetic Modeling of Three-Dimensional Using Integral Equations [J]. Geophysics，1984，49(1): 60-74.

[13] Born M. Optic [M]. New York: Springer，1933.

[14] Zhdanov M S，Lee S K，Yoshioka K. Integral Equation Method for 3D Modeling of Electromagnetic Fields in Complex Structures with Inhomogeneous Background Conductivity [J]. Geophysics，2006，71(6): G333-G345.

[15] Gao G Z，Torres-Verdin C，F(xiàn)ang S. Fast 3D Modeling of Borehole Induction Measurements in Dipping and Anisotropic Formations Using a Novel Approximation Technique [J]. Petrophysics，2004，45(4): 335-349.

[16] 張金會(huì)，孫建國. 三維直流電場積分方程中奇異性的近似處理[J]. 吉林大學(xué)學(xué)報(bào)(地球科學(xué)版)，2009，39(5): 923-928.

Zhang Jinhui，Sun Jianguo. Treatment of Singularity in Integration Equations for 3D DC Electrical Field [J]. Journal of Jilin University (Earth Science Edition)，2009，39(5): 923-928.

[17] Xiong Z. Electromagnetic Fields of Electric Dipoles Embedded in a Stratified Anisotropic Earth [J]. Geophysics，1989，54(12): 1643-1646.

[18] Xiong Z. Electromagnetic Modeling of 3-D Structures by the Method of System Iteration Using Integral Equations [J]. Geophysics，1992，57(12): 1556-1561.

[19] 張輝，李桐林，董瑞霞. 基于電偶源的體積分方程法三維電磁反演[J]. 吉林大學(xué)學(xué)報(bào)(地球科學(xué)版)，2006，36(2): 284-288.

Zhang Hui，Li Tonglin，Dong Ruixia. 3D Electromagnetic Inversion by Volume Integral Equation Method Based on Current Dipole Source [J]. Journal of Jilin University (Earth Science Edition)，2006，36(2): 284-288.

[20] Newman G A，Alumbaugh D L. Frequency-Domain Modeling of Airborne Electromagnetic Responses Using Staggered Finite Differences [J]. Geophysical Prospecting，1995，43: 1021-1042.

[21] 黃威，殷長春，賁放，等. 頻率域航空電磁三維矢量有限元正演模擬[J]. 地球科學(xué)，2016，41(2): 331-342.

Huang Wei，Yin Changchun，Ben Fang，et al. 3D Forward Modeling for Frequency AEM by Vector Finite Element [J]. Earth Science，2016，41(2): 331-342.

[22] 張博，殷長春，劉云鶴，等. 起伏地表頻率域/時(shí)域航空電磁系統(tǒng)三維正演模擬研究[J]. 地球物理學(xué)報(bào)，2016，59(4): 1506-1520.

Zhang Bo，Yin Changchun，Liu Yunhe，et al. 3D Modeling on Topographic Effect for Frequency-/Time-Domain Airborne EM Systems [J].Chinese Journal Geophysics，2016，59(4): 1506-1520.

[23] 王德智.基于積分方程技術(shù)的三維電磁法正演模擬研究[D].長春：吉林大學(xué)，2015.

Wang Dezhi. Three-Dimensional EM Forward Modeling Based on Integral Equation Method [D]. Changchun: Jilin University，2015.

[24] Chew W C，Tong M S，Hu B. Integral Equation Me-thods for Electromagnetic and Elastic Waves [M]. San Rafeal: Morgan & Claypool Publishers，2009.