空間-波數(shù)域三維大地電磁場(chǎng)積分方程法數(shù)值模擬

戴世坤, 陳輕蕊*, 凌嘉宣, 李昆

1 中南大學(xué)地球科學(xué)與信息物理學(xué)院, 長(zhǎng)沙 410083 2 中南大學(xué)有色金屬成礦預(yù)測(cè)與地質(zhì)環(huán)境監(jiān)測(cè)教育部重點(diǎn)實(shí)驗(yàn)室, 長(zhǎng)沙 410083 3 西南石油大學(xué)地球科學(xué)與技術(shù)學(xué)院, 成都 610500

0 引言

電磁法勘探廣泛應(yīng)用于探測(cè)地殼物質(zhì)結(jié)構(gòu)、普查石油天然氣、煤田、地?zé)岷蛯ふ业叵滤徒饘俚V產(chǎn)等資源中.頻率域電磁法可分為天然場(chǎng)源和人工源法.電磁法正演方法主要有有限差分法、有限單元法、有限體積法和積分方程法四種.有限差分法原理簡(jiǎn)單，是以差分和差商代替求導(dǎo)的一種數(shù)值方法，目前應(yīng)用得比較多的是Yee(1966)創(chuàng)立的交錯(cuò)網(wǎng)格，能夠很好的處理場(chǎng)在介質(zhì)內(nèi)部的不連續(xù)性，在求解大地電磁各向異性介質(zhì)中的場(chǎng)(殷長(zhǎng)春等，2014；薛帥等，2017)、帶地形(Mackie et al.，1993，1994；陳伯舫等，1998；董浩等，2014)和算法并行(Tan et al.，2006；李焱等，2012；秦策等，2017)等方面都有突破， Varilsüha和Candansayar(2018)對(duì)比了基于不同規(guī)范下大地電磁控制方程的有限差分正演精度和速度，總結(jié)了direct EM, the Coulomb-gauged, the ungauged, the Lorenz, 和the axial-gauged formulations這四種方法在高頻和低頻迭代計(jì)算時(shí)的一些規(guī)律.羅威等(2019)基于交錯(cuò)網(wǎng)格有限差分法研究了球坐標(biāo)系下的大地電磁三維正演；有限單元法是依據(jù)變分原理或加權(quán)余量法推導(dǎo)出的與定解問(wèn)題相等價(jià)的積分弱解形式，是求解邊值問(wèn)題中數(shù)學(xué)理論最為完備的數(shù)值方法，在地球物理中應(yīng)用的較多的可分為節(jié)點(diǎn)有限元和矢量有限元法.節(jié)點(diǎn)有限元(Zyserman and Santos，2000；Mitsuhata and Uchida，2004；Farquharson and Miensopust，2011；馮建新等，2012；Grayver and Bürg，2014)不滿足電場(chǎng)法向分量不連續(xù)性，需要進(jìn)行散度校正.矢量有限元可以確保不同物性邊界處切線方向場(chǎng)的連續(xù)性，且自動(dòng)滿足無(wú)源區(qū)散度為零的條件，不需要散度校正，可以克服傳統(tǒng)節(jié)點(diǎn)有限元出現(xiàn)偽解的不足，由于這個(gè)特性，矢量有限元(Liu et al.，2008；顧觀文等，2014；Ren et al.，2014；Kordy et al.，2016；Prihantoro et al.，2016)在大地電磁數(shù)值模擬中得到了很好的應(yīng)用，已逐漸代替節(jié)點(diǎn)有限元，成為地球物理電磁場(chǎng)數(shù)值模擬的標(biāo)準(zhǔn)方式(湯井田等，2015).有限體積法又稱控制體積法，先將整個(gè)計(jì)算區(qū)域離散成一系列不重疊的控制網(wǎng)格，然后將微分方程在每一個(gè)控制體積進(jìn)行體積分，單元合成得到線性方程組.有限體積法(Haber and Ascher, 2000; Haber and Ruthotto, 2014；Weiss，2013；Jahandari and Farquharson，2014；陳輝等，2018)是有限差分和有限單元法的一種中間產(chǎn)物，近十幾年才發(fā)展比較迅速，相對(duì)于有限單元法，計(jì)算精度略低但效率更高.上述三種算法需要對(duì)整個(gè)計(jì)算區(qū)域進(jìn)行剖分，對(duì)于大規(guī)模使得形成的矩陣方程大、要求解的未知量個(gè)數(shù)多，對(duì)計(jì)算機(jī)的性能要求高.

相比之下，積分方程法有如下優(yōu)點(diǎn)：(1)只需對(duì)異常體進(jìn)行剖分，計(jì)算量小，占用內(nèi)存少；(2)積分方程的解析解具有半解析解的精度，常常被用于測(cè)試新開(kāi)發(fā)算法的精度；(3)基于積分方程的反演算法精度和效率高(任政勇等，2017)，所以開(kāi)發(fā)高精度高效率的積分方程正演算法仍具有一定的研究?jī)r(jià)值.Raiche(1974)、Hohmann(1971，1975)及Wannamaker和Hohmann(1984)為了積分方程技術(shù)在三維電磁模擬中的應(yīng)用做了很多基礎(chǔ)性的工作；Berdichevsky和Zhdanov(1984)介紹了大地電磁場(chǎng)譜域滿足的方程和表達(dá)式；Wannamaker(1991)對(duì)大地電磁積分方程法的精度和效率進(jìn)行了初探；陳久平等(1990)，殷長(zhǎng)春和樸化榮(1994)分別對(duì)半空間、兩層及多層大地介質(zhì)中的3D異常體的電磁場(chǎng)積分方程法正演進(jìn)行了研究.鮑光淑等(1999)進(jìn)行了均勻半空間頻率域三維電磁響應(yīng)的數(shù)值模擬,深入研究了均勻?qū)щ姲肟臻g頻域三維電磁散射問(wèn)題；Zhdanov和Fang(1997)，Hursán和Zhdanov(2002)使用積分方程法模擬了3D地電結(jié)構(gòu)的電磁響應(yīng)，改進(jìn)了格林函數(shù)算子；Zhdanov等(2006)提出了一種新的積分方程(IE)方法，用于非均勻背景電導(dǎo)率(IBC)復(fù)雜結(jié)構(gòu)的三維電磁(EM)建模.美國(guó)猶他大學(xué)(2006)推出了一個(gè)基于積分方程法的三維正演軟件INTEM3D，該軟件可以對(duì)水平層狀介質(zhì)中的三維地電結(jié)構(gòu)用積分方程法進(jìn)行頻率域電磁模擬；后來(lái)學(xué)者研究了(Farquharson et al.，2006；Zhdanov et al.，2007)用于大電導(dǎo)率對(duì)比模型三維電磁建模的新方法；徐凱軍等(2006)提出利用結(jié)合連分式展開(kāi)的高斯求積代替常規(guī)的快速漢克爾變換方法，進(jìn)行了半空間大地電磁正演模擬；張輝等(2006)用直接求解體積分方程的方法模擬了電偶源激發(fā)時(shí)均勻?qū)щ姲肟臻g頻率域三維電磁響應(yīng)，張量格林函數(shù)采用差分近似的方法求解，結(jié)合連分式展開(kāi)的高斯求積求解含有貝塞爾積分，取得了較好的計(jì)算效果；阮百堯等(2007)提出了一種用邊界元法計(jì)算大地電磁場(chǎng)三維地形影響的數(shù)值模擬方法；張羅磊等(2010)基于MNS技術(shù)進(jìn)行了大地電磁三維正演模擬，提出在空間-波數(shù)域計(jì)算張量格林函數(shù)，大大提高了計(jì)算效率；陳桂波等(2009)采用積分方程法進(jìn)行了各向異性地層電磁場(chǎng)三維數(shù)值模擬研究，分析了各向異性對(duì)三維電磁響應(yīng)的影響特征和規(guī)律；Zaslavsky等(2011)采用有限差分和積分方程混合的方法(稱為有限差分積分方程法)，降低了每次迭代的計(jì)算成本；任政勇等(2017)采用四面體單元、解析的并矢Green函數(shù)奇異積分表達(dá)式，借助于PARDISO高性能并行直接求解器，進(jìn)行了三維大地電磁積分方程正演.

但是上述方法最終都?xì)w結(jié)于大型線性方程組的求解，三維電磁場(chǎng)空間域數(shù)值模擬計(jì)算量大，存儲(chǔ)要求高，限制了現(xiàn)有方法的大規(guī)模應(yīng)用.針對(duì)這一問(wèn)題，本文提出一種空間-波數(shù)域三維電磁場(chǎng)數(shù)值模擬方法，該方法利用電磁場(chǎng)積分方程是一個(gè)三維卷積的特性，沿水平方向進(jìn)行二維傅里葉變換，將三維空間域卷積問(wèn)題轉(zhuǎn)換為多個(gè)不同波數(shù)的空間垂向一維積分問(wèn)題，一維積分相對(duì)獨(dú)立，并行度好，由此大大減少了計(jì)算量和存儲(chǔ)需求；保留垂向?yàn)榭臻g域，將一維積分垂向可離散為多個(gè)單元積分之和，每個(gè)單元采用二次形函數(shù)表征散射電流變化，得出單元積分的解析表達(dá)式，計(jì)算精度高、效率高，垂向單元網(wǎng)格靈活，可以準(zhǔn)確模擬任意復(fù)雜模型，兼顧計(jì)算精度與計(jì)算效率；利用壓縮算子，將電磁場(chǎng)采用迭代法求解，占用內(nèi)存小，計(jì)算速度快.該方法充分利用快速傅里葉變換的高效性和一維形函數(shù)積分的高精度特性，實(shí)現(xiàn)了三維電磁場(chǎng)高效高精度數(shù)值模擬.與常規(guī)積分方程法不同，本文方法適用于地下任意復(fù)雜介質(zhì)，計(jì)算效率高；且相較于其他算法，隨著計(jì)算規(guī)模增大，算法優(yōu)勢(shì)越明顯.

1 三維電磁積分方程

1.1 基本原理

E(r)=Eb(r)+?vG(r,r′)Δσ(r′)E(r′)dv,

(1)

式中r(x,y,z)表示觀測(cè)點(diǎn)位置，r′(xs,ys,zs)表示異常體位置，v為異常體體積，E(r)表示r處的總電場(chǎng)，Eb(r)為r處的背景電場(chǎng)，式中G(r,r′)為r′處的源在r處的電場(chǎng)張量格林函數(shù)，式中Δσ是異常導(dǎo)電率，Δσ(r′)=σ(r′)-σb(r′)，σ(r′)是異常體的導(dǎo)電率，σb(r′)是背景導(dǎo)電率.

設(shè)散射電場(chǎng)的表達(dá)式為

(2)

式中Es(r)為散射電場(chǎng)，J(r′)為散射電流.

電場(chǎng)與磁場(chǎng)之間的關(guān)系為

(3)

將式(2)進(jìn)行二維傅里葉變換，可得

(4)

式(4)采用迭代法求解.定義一個(gè)線性算子：

G(·)=G(Δσ·E).

(5)

式(1)寫(xiě)為

E=Eb+G(Δσ·E).

(6)

對(duì)于式(6)，理論上可以采用迭代法逐次逼近進(jìn)行求解，即可以寫(xiě)為

E(n+1)=Eb+G(Δσ·E(n)).

(7)

E(n+1)和E(n)分別表示第(n+1)和n次計(jì)算得的總場(chǎng)，n≥0，由泛函分析中的Banach定理可知，要使得式(7)收斂的條件是

‖G(Δσ·(E(n+1)-E(n)))‖<κ‖Δσ·(E(n+1)-E(n))‖,

(8)

式中κ<1，‖·‖表示2-范數(shù)，使得上式成立的條件是

‖G‖<1.

(9)

基于Singer(1995)、Pankratov等(1997)，Zhdanov和Fang(1997)和Avdeev等(2002)提出的一種波恩級(jí)數(shù)法， Hursán和Zhdanov(2002)從能量不等式出發(fā)，對(duì)算子G進(jìn)行了修正，構(gòu)造了滿足式(9)的壓縮算子，該算子能使得迭代穩(wěn)定收斂(Gao and Torres-Verdin,2006).其迭代格式如下:

E(n)=Eb+G(Δσ·E(n-1)),

(10)

E(n)=αE(n)+βE(n-1),

(11)

式(11)中左端項(xiàng)E(n)是通過(guò)壓縮算子更新的第n迭代的總場(chǎng)，它將作為第(n+1)次迭代計(jì)算的初值，式中α,β是與背景電導(dǎo)率σb、異常體電導(dǎo)率與背景電導(dǎo)率的差Δσ有關(guān)的張量:

(12)

(13)

利用水平方向二維傅里葉變換，將散射電流滿足的三重卷積轉(zhuǎn)化為多個(gè)空間域垂向的一維積分，不同波數(shù)之間的積分相互獨(dú)立，并行性好；垂向一維積分離散為多個(gè)單元積分之和，每個(gè)單元采用二次形函數(shù)表征散射電流，求得單元積分的解析表達(dá)式，垂向網(wǎng)格剖分靈活，能兼顧計(jì)算精度和計(jì)算效率；用迭代求解電磁場(chǎng)，占用內(nèi)存小，計(jì)算速度快.

1.2 波數(shù)域格林函數(shù)

張量格林函數(shù)表示點(diǎn)源(單元偶極子源)的響應(yīng)，是一個(gè)3×3的矩陣，在直角坐標(biāo)系中表示為

(14)

式中G的第一列表示x方向單位電偶極子源產(chǎn)生電場(chǎng)，第二列表示y方向單位電偶極子源產(chǎn)生的電場(chǎng)，第三列表示z方向單位電偶極子源產(chǎn)生的電場(chǎng).求解張量格林函數(shù)可以轉(zhuǎn)化為求解x、y、z三個(gè)方向電偶極子源產(chǎn)生的電場(chǎng)問(wèn)題.

本文采用的電導(dǎo)率背景模型為均勻?qū)訝罱橘|(zhì)模型，推導(dǎo)均勻?qū)訝罱橘|(zhì)波數(shù)域張量格林函數(shù)的基本思路是：從Maxwell方程組出發(fā)，引入矢量位和標(biāo)量位，將方程轉(zhuǎn)換成矢量位和標(biāo)量位的方程，再帶入洛倫茲規(guī)范條件，得到矢量位所滿足的Helmholtz方程.對(duì)亥姆霍茲方程進(jìn)行x、y方向的傅里葉變換，將三維偏微分方程降為一維常微分方程，帶入邊界條件，得到位的波數(shù)域表達(dá)式，利用波數(shù)域位與場(chǎng)之間的關(guān)系，求解波數(shù)域場(chǎng)的表達(dá)式，分別求解x、y、z三個(gè)方向電偶極子源產(chǎn)生的電磁場(chǎng)得波數(shù)域表達(dá)式，即可合成求得波數(shù)域的張量格林函數(shù).推導(dǎo)過(guò)程和表達(dá)式詳見(jiàn)附錄A.

根據(jù)傅里葉變換的微分性質(zhì)，將推導(dǎo)空間域張量格林函數(shù)的思路應(yīng)用到波數(shù)域張量格林函數(shù)的推導(dǎo)，波數(shù)域張量格林函數(shù)表達(dá)式中無(wú)復(fù)雜的漢克爾積分，將空間域奇異點(diǎn)的計(jì)算轉(zhuǎn)化為零波數(shù)的計(jì)算，方法簡(jiǎn)單，計(jì)算量大大減少，效率高；能大大提升積分方程的正演計(jì)算速度.

1.3 一維積分計(jì)算

式(4)中，垂向積分為空間域，將zs方向的一維積分離散為多個(gè)單元之和，波數(shù)域張量格林函數(shù)的指數(shù)項(xiàng)表達(dá)式中存在變量zs，將波數(shù)域格林函數(shù)中與z，zs無(wú)關(guān)的系數(shù)提取到積分外，剩余積分的表達(dá)式可歸納為通式

(15)

2 算法

利用迭代法求電場(chǎng)，計(jì)算電場(chǎng)的流程如圖1所示，利用如圖所示流程求得電場(chǎng)數(shù)值后，利用式(3)可求得磁場(chǎng).

圖1 空間波數(shù)域算法流程圖Fig.1 Flow chart of the space-wavenumber domain algorithm to calculate electric field

圖2 模型平面示意圖Fig.2 Model schematic plane

3 算例

3.1 算法驗(yàn)證

將平面波作為一次場(chǎng)，進(jìn)行大地電磁場(chǎng)數(shù)值模擬，用美國(guó)猶他大學(xué)開(kāi)發(fā)的基于積分方程法的三維正演軟件INTEM3D的計(jì)算結(jié)果為參照，驗(yàn)證方法的正確性.測(cè)試的計(jì)算機(jī)為Intel(R) Core(TM) i7-6700HQ CPU主頻為2.60 GHz,內(nèi)存為16 GB、64位win10系統(tǒng)，算法在Microsoft Visual Studio 2015開(kāi)發(fā)平臺(tái)上運(yùn)行.

模型XOY平面投影如圖2所示，背景為均勻半空間介質(zhì)，上半空間為空氣，導(dǎo)電率σ1=10-12S·m-1，下半空間導(dǎo)電率σb=0.01 S·m-1，將頻率分別設(shè)為0.01 Hz、1 Hz、100 Hz和10000 Hz，進(jìn)行大地電磁場(chǎng)三維數(shù)值模擬，計(jì)算范圍x方向-1000～500 m，y方向-1000～500 m，剖分網(wǎng)格節(jié)點(diǎn)個(gè)數(shù)101×101×101，三個(gè)方向均勻剖分，Δx、Δy均為10 m，異常體范圍x方向-100～100 m，y方向-200～200 m，異常體導(dǎo)電率σ=0.1 S·m-1.基于趨膚深度的考慮，頻率為0.01 Hz、1 Hz和100 Hz的z方向計(jì)算范圍設(shè)為0～1000 m,異常體范圍z方向200～400 m，Δz為10 m；當(dāng)頻率為10000 Hz時(shí)，z方向計(jì)算范圍設(shè)為0～400 m,異常體范圍z方向40～120 m，Δz為4 m.傅里葉變換采用4個(gè)高斯點(diǎn)的Gauss-FFT法(Wu and Tian, 2014)實(shí)現(xiàn).

設(shè)迭代終止的條件為：相鄰兩次迭代所有節(jié)點(diǎn)電場(chǎng)模的總和的相對(duì)誤差小于10-4.表達(dá)式為

(16)

式中|En|表示第n次迭代的總場(chǎng)的模，|En+1|表示第n+1次迭代的總場(chǎng)的模.

圖3和圖4分別是頻率為1 Hz本文算法(space wavenumber domain integral electromagnetic method, SWIEM)的數(shù)值解與軟件INTEM3D計(jì)算的地面視電阻率和相位等值線圖，從圖中可看出，數(shù)值解和傳統(tǒng)積分方程解的等值線吻合程度高，ρxx、ρxy、ρyx和ρyy分量的相對(duì)均方根誤差(Wu, 2016)分別為1.12%、0.052%、0.062%、1.13%，φxx、φxy、φyx和φyy分量的相對(duì)均方根誤差分別為3.12%、0.0024%、0.0011%、5.13%，ρxx、ρyy、φxx和φyy數(shù)量級(jí)小，舍入誤差稍大，ρxy、φxy、φyx和ρyx誤差小于1‰，表明算法正確.

圖3 地面視電阻率等值線圖Fig.3 The profile contour map of the apparent resistivity on the ground

圖4 地面相位等值線圖Fig.4 The profile contour map of the phase on the ground

圖5為y=0 km不同頻率視電阻率和相位對(duì)應(yīng)的相對(duì)誤差曲線.圖中，不同頻率視電阻率和相位相對(duì)誤差均小于0.5%. 通過(guò)對(duì)比不同頻率的電阻率和相位相對(duì)誤差，進(jìn)一步表明本文空間波數(shù)域電磁三維積分方程數(shù)值模擬方法對(duì)天然大地電磁場(chǎng)的適應(yīng)性、穩(wěn)定性和可靠性.

圖5 地面y=0 km不同頻率視電阻率和相位的相對(duì)誤差曲線Fig.5 The relative errors of apparent resistivity and phase for different frequencies in line y=0 km on the ground

3.2 迭代算法收斂性

利用3.1節(jié)低頻驗(yàn)證模型，分別改變異常體大小(異常體頂面埋深不變，計(jì)算頻率為10 Hz)、異常體頂面深度(異常體大小不變，計(jì)算頻率為1 Hz)和計(jì)算頻率，記錄電場(chǎng)三分量達(dá)到迭代終止條件所需的迭代次數(shù)，結(jié)果如表1、表2和表3所示.

表1 異常體大小與迭代次數(shù)Table 1 Sizes of anomalies and iteration times

表2 異常體頂面深度與迭代次數(shù)Table 2 Depths of anomalies and iteration times

表3 計(jì)算頻率與迭代次數(shù)Table 3 Frequencies and iteration times

綜合表1—3可知，分別改變異常體大小、頂面埋深和頻率大小，電場(chǎng)三分量達(dá)到計(jì)算精度的迭代次數(shù)相差不大，說(shuō)明異常體大小、頂面埋深和頻率對(duì)算法的收斂速度影響較小，幾乎不影響算法的迭代收斂性.

而背景導(dǎo)電率與異常體導(dǎo)電率的差異對(duì)算法收斂的速度影響大.研究導(dǎo)電率差異對(duì)算法收斂速度的影響，保持異常體大小和埋深不變，設(shè)計(jì)與背景導(dǎo)電率不同差異異常體，記錄其迭代次數(shù).利用3.1節(jié)模型，設(shè)場(chǎng)源為x方向極化，頻率為1 Hz，分別設(shè)高阻和低阻異常體的導(dǎo)電率為背景導(dǎo)電率5倍，10倍，50倍，100倍和1000倍，分別記錄電場(chǎng)三分量達(dá)到迭代終止條件所需的迭代次數(shù).

表4和表5分別為不同導(dǎo)電率差異的低阻異常體和高阻異常體的電場(chǎng)三分量達(dá)到計(jì)算精度的迭代次數(shù).表中，高阻和低阻異常的迭代次數(shù)都隨著異常導(dǎo)電率與背景導(dǎo)電率差異的增大而增多；高阻達(dá)到計(jì)算精度的迭代次數(shù)比低阻少，收斂更快；因一次場(chǎng)為x方向極化，所以電場(chǎng)Ex分量的收斂慢，迭代次數(shù)比其他兩個(gè)分量多，Ey分量收斂最快；若將一次場(chǎng)改為y方向極化，則Ey分量收斂慢，Ez次之，Ex分量收斂最快.

表4 低阻異常與迭代次數(shù)Table 4 Low-resistivity anomalies and iteration times

表5 高阻異常與迭代次數(shù)Table 5 High-resistivity anomalies and iteration times

3.3 計(jì)算效率

利用3.1節(jié)的低頻模型，設(shè)頻率為10 Hz，進(jìn)行數(shù)值模擬正演，記錄不同網(wǎng)格剖分個(gè)數(shù)算法迭代一次的耗時(shí)和占用內(nèi)存，探究算法的計(jì)算效率.

表6是不同網(wǎng)格采用標(biāo)準(zhǔn)FFT迭代一次的耗時(shí)與內(nèi)存占用，圖6為其相應(yīng)的變化曲線圖，結(jié)合圖、表發(fā)現(xiàn)，隨著網(wǎng)格個(gè)數(shù)的增多，迭代一次所用時(shí)間和所占內(nèi)存近似呈線性增長(zhǎng)，當(dāng)網(wǎng)格數(shù)為200×200×100(4080501個(gè)節(jié)點(diǎn))時(shí)，本文算法串行迭代計(jì)算一次的時(shí)間約為4 s，正演計(jì)算總耗時(shí)約50 s，內(nèi)存占用約1.3 G，計(jì)算速度快，占用內(nèi)存少，在普通的筆記本上也能高效計(jì)算.

表6 不同規(guī)模串行迭代一次耗時(shí)與內(nèi)存占用Table 6 Time consumption and memory consumption of one iteration for the different grids

圖6 不同網(wǎng)格迭代一次所需時(shí)間和內(nèi)存Fig.6 Time and memory required to iterate through different grids

3.4 復(fù)雜模型適應(yīng)性

利用陳輝等(2018)設(shè)計(jì)的復(fù)雜模型，如圖7所示，計(jì)算范圍22.6 km×30.3 km×20.6 km，背景采用均勻半空間介質(zhì)，上半空間為空氣，下半空間背景電阻率為100 Ωm，計(jì)算頻率為0.01 Hz，三個(gè)異常體電阻率分別為ρ1=10 Ωm，ρ2=1 Ωm，ρ3=1000 Ωm.圖7為該模型在YOZ剖面和XOY平面示意圖.水平x、y方向均勻剖分，z方向非均勻剖分.以x極化方向的電磁場(chǎng)為例，記錄不同網(wǎng)格算法迭代一次耗時(shí)與計(jì)算總耗時(shí)，并對(duì)比文獻(xiàn)(陳輝等，2018)中AGMG-GCR算法的計(jì)算效率.表7為本文算法不同網(wǎng)格的計(jì)算效率.表8為陳輝等(2018)中AGMG-GCR算法的效率，因本文與文獻(xiàn)(陳輝等，2018)中的算法平臺(tái)和程序運(yùn)行環(huán)境一致，能進(jìn)行對(duì)比.

圖7 復(fù)雜模型示意圖Fig.7 Schematic diagram of the complex model

表7 本文算法計(jì)算效率Table 7 Computational efficiency of the SWIEM

表8 AGMG-GCR算法不同網(wǎng)格計(jì)算效率Table 8 Computational efficiency of the AGMG-GCR

圖8為不同算法計(jì)算不同未知量個(gè)數(shù)時(shí)x極化方向電磁場(chǎng)時(shí)的耗時(shí)，圖(a)表示迭代一次的耗時(shí)，圖(b)表示計(jì)算總耗時(shí).由圖8可以看出，兩種方法的計(jì)算耗時(shí)隨著未知量個(gè)數(shù)的增加均呈現(xiàn)出線性增長(zhǎng)的規(guī)律.在相同的算法平臺(tái)和程序運(yùn)行環(huán)境下，對(duì)比AGMG-GCR方法，本文算法效率高，且隨著網(wǎng)格個(gè)數(shù)的增多，本文算法耗時(shí)的優(yōu)勢(shì)越來(lái)約明顯，待求解未知量約為4000000個(gè)時(shí)，AGMG-GCR算法迭代一次耗時(shí)約3 s，本文串行算法約為0.5 s，速度快5倍；待求解未知量約為6000000個(gè)時(shí)，AGMG-GCR方法計(jì)算x極化方向電場(chǎng)總耗時(shí)約1000 s，本文串行算法耗時(shí)約100 s，速度快9倍.相比陳輝等(2018)中計(jì)算效率最高的AGMG-GCR方法，本文算法耗時(shí)更短，為大規(guī)模問(wèn)題的大地電磁三維數(shù)值模擬提供重要的技術(shù)保障.

圖8 x極化方向電磁場(chǎng)計(jì)算耗時(shí)(a) 迭代一次耗時(shí); (b) 計(jì)算總耗時(shí).Fig.8 Computational efficiency of electromagnetic field in x polarization direction(a) Computation time to iterate once; (b) Computation time to calculate total.

圖9為該模型用不同數(shù)值模擬方法計(jì)算得到的阻抗Zxy和Zyx的實(shí)部與虛部的解，測(cè)點(diǎn)位于地面y=0 m，x方向-4000～4000 m，共21個(gè)測(cè)點(diǎn)，每個(gè)測(cè)點(diǎn)間距400 m.其中FDM表示有限差分法計(jì)算的結(jié)果，是采用大地電磁三維正反演代碼ModEM(Kelbert et al., 2014; Dong and Egbert,2019)計(jì)算得到；FEM為大地電磁三維有限元直接解法的計(jì)算結(jié)果(Xiong et al.,2018)，SWIEM為本文算法的計(jì)算結(jié)果.從圖中可以看出，三種算法計(jì)算得到的阻抗實(shí)部和虛部吻合得很好.本文算法與有限元結(jié)果的最大相對(duì)誤差為1.1%，與有限差分法結(jié)果的最大相對(duì)誤差為1.5%，誤差均在可接受的范圍內(nèi)，表明本文算法能在滿足精度要求的前提下達(dá)到高效率.

圖9 不同數(shù)值模擬方法的計(jì)算結(jié)果對(duì)比Fig.9 Comparison of the results of different numerical simulation methods

圖10為該復(fù)雜模型在地面的視電阻率和相位圖，組合異常體電阻率差異大，在異常體位置異常明顯，本文算法對(duì)復(fù)雜模型的適應(yīng)性強(qiáng)，適用于大規(guī)模三維地下任意復(fù)雜模型大地電磁快速正演模擬.

圖10 地面視電阻率和相位圖Fig.10 Apparent resistivity and phase diagram of the ground

4 結(jié)論

利用空間域電磁場(chǎng)積分方程為三重卷積的特點(diǎn)，本文提出一種空間-波數(shù)域數(shù)值模擬方法，該方法借助水平方向二維傅里葉變換，將三維空間域卷積問(wèn)題轉(zhuǎn)換為多個(gè)波數(shù)之間相對(duì)獨(dú)立的空間垂向一維積分問(wèn)題，由此大大減少了計(jì)算量和存儲(chǔ)需求并且算法高度并行.一維積分垂向可離散為多個(gè)單元積分之和，每個(gè)單元采用二次形函數(shù)表征電流變化，可得出單元積分的解析表達(dá)式.保留垂向?yàn)榭臻g域，優(yōu)勢(shì)之一在于可根據(jù)實(shí)際情況靈活調(diào)整單元疏密程度，兼顧計(jì)算精度與計(jì)算效率；優(yōu)勢(shì)之二在于用形函數(shù)擬合求得積分的解析解，計(jì)算精度和效率高.該方法充分利用一維形函數(shù)積分的高效和高精度、不同波數(shù)一維積分之間相對(duì)獨(dú)立高度并行和快速傅里葉變換的高效性, 實(shí)現(xiàn)電磁場(chǎng)高效高精度的三維數(shù)值模擬.將積分方程軟件INTEM3D的數(shù)值解與本文算法的數(shù)值解對(duì)比驗(yàn)證了算法正確；異常體與背景場(chǎng)的導(dǎo)電率差異越大，算法所需迭代次數(shù)越多，高阻異常體比低阻異常體收斂速度快；隨著網(wǎng)格個(gè)數(shù)的增多，算法耗時(shí)和所占內(nèi)存近似呈線性增長(zhǎng)，與目前主流方法相比，速度快一個(gè)數(shù)量級(jí)以上，且隨著計(jì)算規(guī)模的增大，算法優(yōu)勢(shì)越明顯. 本文的正演算法為大規(guī)模高效、高精度的反演研究奠定了基礎(chǔ).

本文算法在不同波數(shù)一維積分計(jì)算和不同深度節(jié)點(diǎn)電磁場(chǎng)正、反傅里葉變換均具有高度并行性，采用并行計(jì)算，計(jì)算效率將進(jìn)一步提升.此外，本文僅研究大地電磁場(chǎng)三維數(shù)值模擬，若將背景場(chǎng)改為均勻?qū)訝罱橘|(zhì)可控源電磁場(chǎng)，即可進(jìn)行可控源電磁場(chǎng)三維數(shù)值模擬.

值得注意的是，本文采用的壓縮算子在低阻異常體導(dǎo)電率與背景導(dǎo)電率對(duì)比度大(>100)時(shí)，迭代次數(shù)需上百次，影響算法效率；在復(fù)雜地質(zhì)條件下，本文算法需要精細(xì)剖分來(lái)擬合地下異常體，增加算法的計(jì)算量.但由于這兩個(gè)問(wèn)題也普遍存在于常規(guī)空間域數(shù)值模擬方法中，且在滿足計(jì)算精度的前提下，相比常規(guī)算法，本文算法效率高，此時(shí)新算法依然有一定的優(yōu)勢(shì).此外，本文算例的傅里葉變換采用標(biāo)準(zhǔn)FFT，網(wǎng)格均勻剖分，下一步將研究非均勻網(wǎng)格條件下的空間-波數(shù)域正演，進(jìn)一步提高算法的適用性和效率.

致謝感謝中國(guó)地質(zhì)大學(xué)(武漢)羅天涯博士后提供復(fù)雜模型有限單元法的數(shù)值解數(shù)據(jù)，感謝中南大學(xué)王旭龍博士生提供復(fù)雜模型ModEM軟件有限差分法的數(shù)值解數(shù)據(jù).感謝審稿專家和編輯提出的寶貴意見(jiàn).

附錄A 波數(shù)域張量格林函數(shù)公式推導(dǎo)

將空間-波數(shù)域電場(chǎng)張量格林函數(shù)的求解轉(zhuǎn)化為x,y,z三個(gè)方向的電偶源產(chǎn)生的場(chǎng)，頻率域Maxwell為

(A1)

(A2)

(A3)

(A4)

(A5)

電場(chǎng)與矢量位的關(guān)系式為

(A6)

(1)全空間波數(shù)域張量格林函數(shù)

設(shè)源為x方向，全空間僅存在矢量位Ax，源點(diǎn)坐標(biāo)為(xs,ys,zs)，將式(A5)進(jìn)行水平方向二維傅里葉變換，可得表達(dá)式

(A7)

(A8)

利用式(A6)在空間-波數(shù)域的表達(dá)式可求得x方向偶極源空間-波數(shù)域電場(chǎng)表達(dá)式.同理可得，y、z方向偶極源空間-波數(shù)域的表達(dá)式，此處不再贅述.

(2)層狀介質(zhì)波數(shù)域張量格林函數(shù)求解

(A9)

(A9′)

(A10)

根據(jù)(考夫曼和凱勒,1987)的推導(dǎo)，將推導(dǎo)過(guò)程轉(zhuǎn)化為在波數(shù)域求解，直接寫(xiě)出矢量位系數(shù)的表達(dá)式.

(i)x/y方向偶極源電磁場(chǎng)

水平矢量位

構(gòu)造：

(A11)

(A12)

設(shè)遞推：

(A13)

(A14)

源層系數(shù)的解為

(A15)

再利用各層矢量位的遞推關(guān)系可得各層系數(shù)表達(dá)式.

垂直矢量位

(A16)

式中X′表示X對(duì)z方向的導(dǎo)數(shù)，設(shè)

Vt=Pteut(z-zt)+Qte-ut(z-zt-1),

(A17)

直接寫(xiě)出V的解析表達(dá)式.

構(gòu)造：

(A18)

(A19)

(A20)

(A21)

源層系數(shù)的解為

(A22)

再利用各層矢量位之間的遞推關(guān)系可得各層系數(shù)表達(dá)式.

(A23)

(A24)

再利用式(A8)在空間-波數(shù)域的表達(dá)式可求得x、y方向偶極源空間-波數(shù)域電場(chǎng)表達(dá)式，此處不再贅述.

(ii)z方向偶極源電磁場(chǎng)

(A25)

(A26)

再利用各層矢量位之間的遞推關(guān)系可得各層系數(shù)表達(dá)式.

同理利用式(A8)在空間-波數(shù)域的表達(dá)式可求得z方向偶極源空間-波數(shù)域電場(chǎng)表達(dá)式，此處不再贅述.

附錄B 一維單元積分解析解

空間-波數(shù)域電磁場(chǎng)一維單元積分的表達(dá)式為

(B1)

(B2)

(B3)

式(B3)的解析解為

(B4)