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

吉林大學(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ù)，它通過對異常區(qū)域進(jìn)行剖分，建立異常區(qū)域內(nèi)電場滿足的第二類Fredholm積分方程，通過直接求解或利用迭代技術(shù)求解離散后的矩陣方程得到空間電磁場。在小規(guī)模異常體模擬中，系數(shù)矩陣占用存儲小，求解速度快。然而，對于大型異常體，由于形成的線性方程組系數(shù)矩陣的密實性，計算耗時多、存儲量大，限制了積分方程算法的應(yīng)用。

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

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

1 MGQL正演理論

根據(jù)積分方程理論，在水平層狀各向同性介質(zhì)中，對麥克斯韋方程組進(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)可以寫成獨立的標(biāo)量式，即

；l=x,y,z。

(11)

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

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

(12)

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

；l=x,y,z。

(13)

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

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

(14)

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

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

2.1 精度驗證

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

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

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

2.2 單個異常體模型

航空電磁勘探在地形條件比較復(fù)雜的區(qū)域應(yīng)用較多，目標(biāo)體上方的覆蓋層會對異常響應(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個，細(xì)網(wǎng)格剖分網(wǎng)格數(shù)為4×10×20個。測線沿x方向，飛行高度30 m，頻率900 Hz。

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

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

2.3 多個異常體模型

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

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

3 計算效率分析

本文提出的MGQL算法的最大優(yōu)點是在保證計算精度的前提下達(dá)到快速計算的目標(biāo)。為了分析MGQL算法的計算效率，本文將其與開源軟件MarcoAir (Version 2.3.1)及全積分方程算法結(jié)果作對比。本文設(shè)計與圖1類似的模型，異常體大小為30 m×200 m×100 m，電阻率1 Ω·m，埋深30 m，飛行高度30 m, 頻率900 Hz，測線沿x方向。通常，全矩陣存儲積分方程法需要占用大量的內(nèi)存空間，用于存儲格林系數(shù)矩陣。受內(nèi)存限制，全矩陣存儲積分方程法計算網(wǎng)格數(shù)有限。使用MGQL技術(shù)能夠減少網(wǎng)格數(shù)、節(jié)約存儲，可實現(xiàn)多網(wǎng)格數(shù)剖分的正演模擬。本文的計算在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ù)的對稱性和塊迭代算法提升積分方程的計算速度[22]。當(dāng)異常體剖分單元較少時，直接解法優(yōu)勢較大；當(dāng)剖分單元較多時，迭代法計算速度較快。我們對比時選用MarcoAir迭代法求解線性方程組。本文MGQL算法中將格林函數(shù)存儲為Toeplitz矩陣形式，同時利用FFT實現(xiàn)矩陣與向量乘積，并采用穩(wěn)定雙共軛梯度法(BICGSTAB)求解線性方程組[23-24]，以進(jìn)一步加快計算速度。假設(shè)航空電磁系統(tǒng)發(fā)射源和接收機(jī)共移動60次。表1給出在相同計算精度條件下，本文算法計算時間與MarcoAir的對比結(jié)果。由表1可見：當(dāng)網(wǎng)格數(shù)較少時，MarcoAir求解速度快于MGQL算法；隨著模型剖分網(wǎng)格數(shù)增加，MGQL算法計算速度遠(yuǎn)遠(yuǎn)超過MarcoAir。例如，當(dāng)網(wǎng)格數(shù)增加到10×20×20時，MGQL的計算速度是MarcoAir近54倍。另外，對于計算設(shè)備的有限存儲空間和時間，MarcoAir可以計算的網(wǎng)格數(shù)不超過20×40×40個，而對于MGQL算法來說，可以計算的網(wǎng)格數(shù)為50×100×100。由于MGQL在粗網(wǎng)格上求解異常場，計算量極大減少，使得該算法在計算效率上擁有很大的優(yōu)勢。

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

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

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

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

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為計算時間；*|*|*|*表示計算時間大于24 h。

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

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

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

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

4 結(jié)論

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

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

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

[1] 殷長春，張博，劉云鶴，等.航空電磁勘查技術(shù)發(fā)展現(xiàn)狀及展望[J].地球物理學(xué)報，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é)報，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] 張金會，孫建國. 三維直流電場積分方程中奇異性的近似處理[J]. 吉林大學(xué)學(xué)報(地球科學(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é)報(地球科學(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] 張博，殷長春，劉云鶴，等. 起伏地表頻率域/時域航空電磁系統(tǒng)三維正演模擬研究[J]. 地球物理學(xué)報，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.