空間-波數域三維大地電磁場積分方程法數值模擬

2022-06-02 01:14:50戴世坤陳輕蕊凌嘉宣李昆

地球物理學報 2022年6期

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

1 中南大學地球科學與信息物理學院, 長沙 410083 2 中南大學有色金屬成礦預測與地質環境監測教育部重點實驗室, 長沙 410083 3 西南石油大學地球科學與技術學院, 成都 610500

0 引言

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

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

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

1 三維電磁積分方程

1.1 基本原理

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

(1)

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

設散射電場的表達式為

(2)

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

電場與磁場之間的關系為

(3)

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

(4)

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

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

(5)

式(1)寫為

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

(6)

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

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

(7)

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

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

(8)

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

‖G‖<1.

(9)

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

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

(10)

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

(11)

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

(12)

(13)

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

1.2 波數域格林函數

張量格林函數表示點源(單元偶極子源)的響應，是一個3×3的矩陣，在直角坐標系中表示為

(14)

式中G的第一列表示x方向單位電偶極子源產生電場，第二列表示y方向單位電偶極子源產生的電場，第三列表示z方向單位電偶極子源產生的電場.求解張量格林函數可以轉化為求解x、y、z三個方向電偶極子源產生的電場問題.

本文采用的電導率背景模型為均勻層狀介質模型，推導均勻層狀介質波數域張量格林函數的基本思路是：從Maxwell方程組出發，引入矢量位和標量位，將方程轉換成矢量位和標量位的方程，再帶入洛倫茲規范條件，得到矢量位所滿足的Helmholtz方程.對亥姆霍茲方程進行x、y方向的傅里葉變換，將三維偏微分方程降為一維常微分方程，帶入邊界條件，得到位的波數域表達式，利用波數域位與場之間的關系，求解波數域場的表達式，分別求解x、y、z三個方向電偶極子源產生的電磁場得波數域表達式，即可合成求得波數域的張量格林函數.推導過程和表達式詳見附錄A.

根據傅里葉變換的微分性質，將推導空間域張量格林函數的思路應用到波數域張量格林函數的推導，波數域張量格林函數表達式中無復雜的漢克爾積分，將空間域奇異點的計算轉化為零波數的計算，方法簡單，計算量大大減少，效率高；能大大提升積分方程的正演計算速度.

1.3 一維積分計算

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

(15)

2 算法

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

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

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

3 算例

3.1 算法驗證

將平面波作為一次場，進行大地電磁場數值模擬，用美國猶他大學開發的基于積分方程法的三維正演軟件INTEM3D的計算結果為參照，驗證方法的正確性.測試的計算機為Intel(R) Core(TM) i7-6700HQ CPU主頻為2.60 GHz,內存為16 GB、64位win10系統，算法在Microsoft Visual Studio 2015開發平臺上運行.

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

設迭代終止的條件為：相鄰兩次迭代所有節點電場模的總和的相對誤差小于10-4.表達式為

(16)

式中|En|表示第n次迭代的總場的模，|En+1|表示第n+1次迭代的總場的模.

圖3和圖4分別是頻率為1 Hz本文算法(space wavenumber domain integral electromagnetic method, SWIEM)的數值解與軟件INTEM3D計算的地面視電阻率和相位等值線圖，從圖中可看出，數值解和傳統積分方程解的等值線吻合程度高，ρxx、ρxy、ρyx和ρyy分量的相對均方根誤差(Wu, 2016)分別為1.12%、0.052%、0.062%、1.13%，φxx、φxy、φyx和φyy分量的相對均方根誤差分別為3.12%、0.0024%、0.0011%、5.13%，ρxx、ρyy、φxx和φyy數量級小，舍入誤差稍大，ρ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不同頻率視電阻率和相位對應的相對誤差曲線.圖中，不同頻率視電阻率和相位相對誤差均小于0.5%. 通過對比不同頻率的電阻率和相位相對誤差，進一步表明本文空間波數域電磁三維積分方程數值模擬方法對天然大地電磁場的適應性、穩定性和可靠性.

圖5 地面y=0 km不同頻率視電阻率和相位的相對誤差曲線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節低頻驗證模型，分別改變異常體大小(異常體頂面埋深不變，計算頻率為10 Hz)、異常體頂面深度(異常體大小不變，計算頻率為1 Hz)和計算頻率，記錄電場三分量達到迭代終止條件所需的迭代次數，結果如表1、表2和表3所示.

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

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

表3 計算頻率與迭代次數Table 3 Frequencies and iteration times

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

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

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

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

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

3.3 計算效率

利用3.1節的低頻模型，設頻率為10 Hz，進行數值模擬正演，記錄不同網格剖分個數算法迭代一次的耗時和占用內存，探究算法的計算效率.

表6是不同網格采用標準FFT迭代一次的耗時與內存占用，圖6為其相應的變化曲線圖，結合圖、表發現，隨著網格個數的增多，迭代一次所用時間和所占內存近似呈線性增長，當網格數為200×200×100(4080501個節點)時，本文算法串行迭代計算一次的時間約為4 s，正演計算總耗時約50 s，內存占用約1.3 G，計算速度快，占用內存少，在普通的筆記本上也能高效計算.

表6 不同規模串行迭代一次耗時與內存占用Table 6 Time consumption and memory consumption of one iteration for the different grids

圖6 不同網格迭代一次所需時間和內存Fig.6 Time and memory required to iterate through different grids

3.4 復雜模型適應性

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

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

表7 本文算法計算效率Table 7 Computational efficiency of the SWIEM

表8 AGMG-GCR算法不同網格計算效率Table 8 Computational efficiency of the AGMG-GCR

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

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

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

圖9 不同數值模擬方法的計算結果對比Fig.9 Comparison of the results of different numerical simulation methods

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

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

4 結論

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

本文算法在不同波數一維積分計算和不同深度節點電磁場正、反傅里葉變換均具有高度并行性，采用并行計算，計算效率將進一步提升.此外，本文僅研究大地電磁場三維數值模擬，若將背景場改為均勻層狀介質可控源電磁場，即可進行可控源電磁場三維數值模擬.

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

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

附錄A 波數域張量格林函數公式推導

將空間-波數域電場張量格林函數的求解轉化為x,y,z三個方向的電偶源產生的場，頻率域Maxwell為