基于空氣地球分離策略的大地電磁三維正演

李袁傲, 任政勇,2,3,4*, 湯井田,2,3,4, 吳啟紅

1 中南大學地球科學與信息物理學院, 長沙 410083 2 中南大學有色金屬成礦預測與地質環境監測教育部重點實驗室, 長沙 410083 3 有色資源與地質災害探查湖南省重點實驗室, 長沙 410083 4 自然資源部覆蓋區深部資源勘查工程技術創新中心， 合肥 230001 5 成都大學建筑與土木工程學院, 成都 610106

0 引言

大地電磁法(magnetotellurics，MT)是一種利用天然交變電磁場研究地球內部電性結構的地球物理方法(Cagniard,1953)，被廣泛運用于礦產資源勘察、巖石圈電性結構研究以及石油天然氣勘探等領域(Chave and Jones, 2012).當前，大地電磁數據采集趨向于三維方式，采集區域地質情況復雜且常具有起伏地形，因此需要研發高效的大地電磁數據三維反演方法.反演的迭代過程中需要多次調用正演算法來計算模型更新量，正演響應的準確性保證了更新模型的正確性.因此，研發快速高精度的大地電磁三維正演算法是大地電磁數據解釋的研究熱點.

大地電磁正演算法主要基于積分方程法，有限差分法和有限單元法.積分方程法只離散化異常體，適合于較小異常體模型的計算.由于受限于簡單模型才有格林函數，積分方程法一般只能用于簡單模型電磁場響應的計算與分析(Hohmann，1975; Lajoie and West，1976; Zhdanov et al.，2006; 王若等，2009).有限差分法采用差分來近似微分，具有方法簡單、編程容易等優點，在積分方程法之后迅速受到了普及，被廣泛用于三維地電磁感應問題的求解 (Mackie et al.，1994; Newman and Alumbaugh，2002; Siripunvaraporn et al.，2002; 沈金松，2003; 譚捍東等，2003; Streich，2009; 李展輝等，2009; Dong and Egbert，2019) .但是，由于差分近似一般適用于直線、不適合于復雜折線和曲線，有限元差分方法在處理任意起伏地形時存在一定的困難.有限元法采用在小區域內利用簡單多項式來近似復雜未知場的思想，在小區域幾何形狀的選取上具有任意性，可以精確逼近三維復雜地質情況以及起伏地形，從而計算出高精度的電磁場響應，因此有限單元法逐漸成為了當前地電磁感應正演問題的主流算法.

Coggon(1971)率先基于有限單元法計算了地球電磁場響應.隨后，地電磁感應問題的有限元計算方法在網格離散化技術、線性方程組求解技術、并行加速技術、后處理方法與技術等方面取得了長足進步(Raiche，1994; Avdeev，2005; B?rner，2010; Newman，2014).與此相對的是，在微分控制方程最優化設計方向的研究成果較少.目前，地電磁感應問題的有限元計算一般采用電場方程、磁場方程和勢場方程.利用電場方程，Nam等(2007)和Farquharson和Miensopust (2011)進行了大地電磁問題的有限元計算；Ren等(2013)將面向目標自適應有限元法應用到大地電磁三維正演計算中.利用磁場方程，Siripunvaraporn等(2002)和Franke等(2007)進行了大地電磁三維場分布的計算.勢場方程方面，Haber(2000)基于A-φ勢實現了大地電磁三維有限元模擬；Mitsuhata和Uchida(2004)實現了基于T-Ω勢的大地電磁三維有限元模擬.

空解現象產生的本質原因為在空氣和地球中采用了同樣的控制方程，由于空氣和地下具有完全不同的電性特征(電導率對比度為無窮大)，采用同樣的控制方程必然導致病態特征.為解決此問題，本文提出了一種新的空氣與地球分離混合公式.解決思路為：首先，將求解區域分解為空氣、地下兩個獨立的區域，分別采用不同的控制方程；然后，利用地空分界面上電磁場連續性條件，進行界面兩側的耦合與信息交換，保證電磁場解的唯一性.基于該解決思路，空氣層中磁場H分解為電矢量勢T和磁標量勢φ，利用Coulomb規范條件，使得矢量勢T和標量勢φ滿足具有Laplace結構的勢場方程.地下區域中，采用雙旋度結構的電場方程.在地下區域采用電場雙旋度方程的原因在于：地下介質的電導率大于零，這使得雙旋度系統在地下區域避免了零空間問題且有良好的穩定性.

最后，我們計算了國際標準測試模型COMMEMI 3D-1、DTM-1與梯形山模型的大地電磁響應，將計算結果與其他算法進行了對比，來驗證本文算法的正確性，以及處理起伏地表問題的能力.我們還進行了系統條件數與特征值的計算和對比分析，來檢驗本文算法對電場雙旋度系統病態性的改進.

1 節點矢量混合有限元系統

利用大地電磁的場源位于求解區域外的一般性假設，在求解區域內(見圖1)，大地電磁問題的電磁場響應滿足頻率域Maxwell方程組(選取時諧因子為e-iω t)：

(1)

其中，E為電場矢量，H為磁場矢量，i為純虛數單位，ω為角頻率，μ為磁導率，

(2)

其中，J為電流密度，σ為電導率，ε為介電常數，

(3)

(4)

圖1 MT地電模型示意圖Ω0表示空氣區域，Ω1表示地下區域，Γ表示地空邊界，Γ0表示空氣外部邊界，Γ1表示地球外部邊界，n0為空氣區域外邊界的外法向向量，n1為地球區域外邊界的外法向向量，n01為地空邊界上空氣一側的外法向向量，方向由空氣指向地球，n10為地空邊界上地球一側的外法向向量，方向由地球指向空氣，n為約定的地空邊界上法向向量，方向由地球指向空氣.Fig.1 Schematic diagram of MT geoelectric modelΩ0 represents the air area, Ω1 represents the underground area, Γ represents the air-earth interface，Γ0 represents the outer boundary of air，Γ1 represents the outer boundary of earth，n0 is the outer normal vector of the outer boundary of the air region, n1 is the outer normal vector of the outer boundary of the earth region, n01 is the outer normal vector of the air side on the air-earth interface, and the direction is from the air to the earth, n10 is the outer normal vector on the side of the earth on the air-earth interface, and the direction is from the earth to the air, n is the normal vector on the agreed air-earth interface, with the direction from the earth to the air.

1.1 空氣層電磁場積分弱形式

(5)

將公式(5)代入公式(2)中得到：

(6)

其中H0表示空氣中的磁場.引入標量函數φ的負梯度來表示公式(6)中的矢量部分，可得到：

(7)

將公式(7)代入公式(2)并取旋度得到：

即：

(8)

將公式(8)中的矢量勢T的雙旋度結構進行化簡，強加Coulomb規范條件來唯一化矢量勢T：

(9)

(10)

(11)

綜上，空氣區域的控制方程為

(12)

為了求解公式(12)，采用Galerkin節點有限元法(Jin，2002)，先假設余量rT為

(13)

余量rφ為

(14)

采用線性插值基函數，四面體單元中任意點處的矢量勢T以及標量勢φ可由節點基函數表示為

(15)

其中Txi，Tyi，Tzi分別為矢量勢T的三個分量，Li為節點基函數，在空氣區域中滿足以下條件：

?Ω0Li·rTdv=0,

?Ω0Li·rφdv=0.

(16)

利用Green恒等式：

(17)

(18)

綜合公式(15)—(18)，給出公式(12)的積分弱形式：

高溫環境下GH3536基體軟化，導致其粘著趨勢上升，存在外部載荷時易發生塑性變形，在磨球滑動過程中以剝層的方式不斷脫落的磨屑被膠合在一起，形成層狀物(見圖5a)并逐漸被撕裂。另外，基體磨痕的局部區域可見沿滑動方向的犁溝，這可能是由于部分磨屑在磨損過程中嵌入到基體中，在滑動中推擠基體使之塑性流動并犁出溝槽，從而引起磨粒磨損。由圖5b可見，NiAlW涂層的磨痕比基體更為光滑和平整，沒有明顯的單個顆粒沉積物脫落和涂層整塊粘著或撕裂現象。涂層的磨損機制以微切削為主，其表面的單個顆粒沉積物凸起部分首先被切削，在外界的反復撞擊與擠壓下發生破碎，最后脫落并造成涂層的磨損[15-16]。

(19)

寫成分量形式為

(20)

公式(19)與公式(20)即為空氣中的控制方程積分弱形式.

1.2 地球電磁場積分弱形式

地下區域中，電場滿足雙旋度方程：

(21)

采用Galerkin矢量有限元法(Jin，2002)，求解公式(21).先假設余量RE為

(22)

采用棱邊基函數，四面體單元中的電場可表示為

(23)

其中E1i為對應棱邊上的電場切向分量，Ni為棱邊基函數，計算區域內滿足以下條件：

?Ω0Ni·REdv=0.

(24)

(25)

其中a,b分別為矢量函數，則

(26)

(27)

(28)

公式(28)即為地下介質區域控制方程的積分弱形式.

1.3 空氣-地球耦合系統

結合公式(19)和公式(28)以及圖1中所定義的邊界外法向向量，空氣地球耦合系統的積分弱形式為

(29a)

(29b)

(29c)

根據矢量公式：

(30)

公式(29a)可以進一步改寫為

(31)

將公式(29)中與地空界面有關的項移到等式左邊，且在地空界面上有n=n10=-n01，可以得到：

(32)

(33)

在地空分界面上施加如下電磁場連續性條件:

(34)

(35a)

(35b)

(35c)

公式(35a)—(35c)中的右端項均可由1D層狀介質的解析解給出，將公式(35)寫成矩陣形式：

(36)

其中：

公式(35)和公式(36)即為空氣地球分離策略的有限元離散方程.

其中i,j=x,y(x≠y)，Zij為阻抗，ρij為視電阻率，φij為相位，E1i和H1j分別為電場E1、磁場H1在對應方向上的分量.

2 數值試驗

本小節所有計算均在CPU主頻為2.30 GHz，RAM超過200 GB的小型AMAX工作站上完成.

2.1 COMMIME 3D-1模型

3D-1模型是國際標準測試COMMEMI項目(Zhdanov et al., 1997)的一個模型.它是均勻半空間中嵌套一個導電棱柱，其頂埋深為250 m，背景電阻率為100 Ωm，棱柱電阻率為0.5 Ωm，尺寸為1000 m×2000 m×2000 m，空氣電阻率為1016Ωm,計算區域尺寸為70 km×70 km×140 km.

計算網格包含148576個四面體單元.傳統雙旋度方程共173335個自由度，本文提出的地空分離公式共168978個自由度，頻率為0.1 Hz.我們將測試結果(Air-Earth decomposition)與高精度自適應有限元計算結果(Curl-Curl system)(Ren et al., 2013)進行對比，圖2展示了視電阻率和相位的對比，圖3展示了這兩種方法的相對誤差.空氣地球分離公式計算用時 19.0 s，電場雙旋度公式計算用時 21.8 s.

圖2和圖3顯示混合公式系統與雙旋度公式系統能夠獲得一致的視電阻率與相位曲線.視電阻率值從測線端點到y=0處逐漸降低，過y=0之后逐漸增大，曲線具有對稱性.視電阻率最大相對誤差0.34%，相位最大相對誤差0.18%，驗證了算法的正確性和精度.為了進一步研究混合系統的數值特征，我們計算并對比了特征值分布(如圖4和圖5).

圖2 3D-1模型視電阻率與相位，測試頻率為0.1 Hz, Curl-Curl system代表了自適應有限元結果(Ren et al., 2013),Air-Earth decomposition 代表了本文計算的結果Fig.2 Comparison of calculated apparent resistivity and phase on the 3D-1 model among our solution (denoted by Air-Earth decomposition) and the reference solution (Curl-Curl system, Ren et al., 2013), at a frequency of 0.1 Hz

圖3 3D-1模型視電阻率與相位相對誤差(高精度自適應有限元計算結果為參考值(Ren et al., 2013))，測試頻率為0.1 HzFig.3 Relative error of apparent resistivity and phase by comparing to the reference solution (Curl-Curl system,Ren et al., 2013) for the 3D-1 model at a frequency of 0.1 Hz

圖4表明兩個系統的最大特征值分布具有相似性.圖5表明地空分離公式系統有更多的非零或模長更大的最小特征值.當一個系統的特征值全部位于復數平面的左半部分時，系統是收斂的.與此相反，當特征值全部位于復數平面的右半部分時，系統則是發散的.從圖4與圖5可以直觀的看出，空氣地球分離公式的特征值更多地集中在復數平面的左半部分且距離原點更遠，因而空氣地球分離公式相對于電場雙旋度公式有更好的穩定性.根據條件數的定義(最大特征值與最小特征值之比)，地空分離系統擁有更小的條件數.為了定量研究地空分離系統的條件數，計算了不同頻率下的條件數(從0.001 Hz，0.01 Hz，…，10 Hz共六個對數等間隔分布的頻點).測試結果如圖6所示，結果表明混合系統的條件數相對于雙旋度系統均降低，下降可達3個數量級.

圖4 3D-1模型上地空分離公式和雙旋度電場公式對應的前1000個最大特征值分布的對比分析，測試頻率為0.1 HzFig.4 The distribution of 1000 maximum eigenvalues for two methods (one using the Air-Earth decomposition,and one using the Curl-Curl system) on the 3D-1 model, at a frequency of 0.1 Hz

圖5 3D-1模型上地空分離公式和雙旋度電場公式對應的后1000個最小特征值分布的對比分析，測試頻率為0.1 HzFig.5 The distribution of 1000 minimum eigenvalues for two methods (one using the Air-Earth decomposition,and one using the Curl-Curl system) on the 3D-1 model at a frequency of 0.1 Hz

圖6 3D-1模型地空分離公式和雙旋度電場公式的條件數對比Fig.6 Comparison of condition numbers for two methods (one using the Air-Earth decomposition, and one using the Curl-Curl system) on the 3D-1 model at different frequencies

2.2 梯形山模型

為了測試對地形的適應性，計算了梯形山模型(Nam et al., 2007).如圖7所示，梯形山頂部為450 m×450 m的正方形，底部為2000 m×2000 m的正方形且高為450 m.地形電阻率為100 Ωm，空氣的電阻率為1016Ωm，測線方向為x=0方向，計算區域尺寸為25 km×25 km×50 km.計算區域包含186168個四面體單元，雙旋度系統共227241個自由度，地空分離系統共219957個自由度，測試頻率為2 Hz.地空分離公式和雙旋度電場公式的模型視電阻率與相位曲線對比如圖8所示，相應的相對誤差見圖9.空氣地球分離公式計算用時 54.26 s，電場雙旋度公式計算用時 59.66 s.

圖7 梯形山模型Fig.7 The trapezoidal hill model

圖8 梯形山模型視電阻率與相位，測試頻率為2 Hz，Curl-Curl system代表了自適應有限元結果(Ren et al., 2013),Air-Earth decomposition代表了本文計算的結果Fig.8 Comparison of calculated apparent resistivity and phase on the trapezoidal hill model among our solution (denoted by Air-Earth decomposition) and the reference solution (Curl-Curl system, Ren et al., 2013) at a frequency of 2 Hz

測試結果顯示：視電阻率最大相對誤差0.25%，相位最大相對誤差0.50%.在視電阻率曲線中，ρxy在梯形山基底部位的電阻率稍高于圍巖的電阻率，而在起伏地形以及山頂部分則低于圍巖電阻率，ρyx則呈現兩端高中間低的電阻率變化趨勢，這是由于電流的地形效應引起的.xy模式中的電流沿y方向(垂直測線方向)傳播并在山體的側面上向山頂彎曲，使得沿測線方向山腳處的電流發散而山脊和山頂處電流聚集，導致了山腳的電阻率高于圍巖而山脊和山頂的電阻率低于圍巖.在yx模式中，電流沿x軸方向(沿測線方向)傳播，電流的傳播路徑發生彎曲且都在x=0方向上匯聚，因而在yx模式中，ρyx低于圍巖電阻率.

為了研究混合系統與雙旋度系統的數值特征，我們計算了特征值分布，測試結果如圖10和圖11.與3D-1模型類似，圖10表明兩個系統具有相似的最大特征值分布，圖11表明地空分離系統擁有模長更大的最小特征值，且地空分離系統的特征值更加集中在復數平面的左半部分且距離原點更遠.

圖10 梯形山模型上地空分離公式和雙旋度電場公式對應的前1000個最大特征值分布的對比分析，測試頻率為1 HzFig.10 The distribution of 1000 maximum eigenvalues for two methods (one using the Air-Earth decomposition,and one using the Curl-Curl system) on the trapezoidal hill model at a frequency of 1 Hz

圖11 梯形山模型上地空分離公式和雙旋度電場公式對應的后1000個最小特征值分布的對比分析，測試頻率為1 HzFig.11 The distribution of 1000 minimum eigenvalues for two methods (one using the Air-Earth decomposition,and one using the Curl-Curl system) on the trapezoidal hill model at a frequency of 1 Hz

我們在這個模型中選取0.01 Hz，0.1 Hz，…，100 Hz共五個頻點測試了系統的條件數.如圖12所示，模型的條件數分布與COMMEMI 3D-1模型的條件數分布特征一致，條件數有1～3個數量級的改善，在10 Hz頻率下兩系統條件數差異最大.

圖12 梯形山模型地空分離公式和雙旋度電場公式的條件數對比Fig.12 Comparison of condition numbers of the Air-Earth decomposition and the Curl-Curl system on the trapezoidal hill model at different frequencies

2.3 DTM-1模型

該模型由嵌套在均勻半空間中的三個異常體組成，其模型示意圖如圖13，其中ρ1=10 Ωm,ρ2=1 Ωm,ρ3=10000 Ωm，介質電阻率最大對比度可達10000∶1，模型外邊界距離地面250 km.計算區域尺寸為250 km×250 km×500 km，測線沿y=0方向分布.

圖13 DTM-1模型結構Fig.13 Structure of DTM-1 model

計算區域包含289178個四面體單元，雙旋度系統共343978個自由度，地空分離系統共343615個自由度.選取測試頻率為0.01 Hz進行計算，空氣地球分離公式計算用時 70.47 s，電場雙旋度公式計算用時 71.86 s.模型的電磁響應曲線如圖14所示，相對誤差如圖15所示，結果顯示兩種系統的響應曲線吻合度十分高，視電阻率最大相對誤差0.20%，相位最大相對誤差0.43%.與梯形山模型視電阻率曲線分析類似，該模型的視電阻率的變化趨勢依舊可以根據電流的傳播方向來分析.由于ρ3=10000 Ωm高于圍巖電阻率，使得通過異常體的電流會向周圍發散，則在x=[-25000 m,0]的范圍內電流向測點聚集，反應在視電阻率曲線上則是這一段范圍內的視電阻率呈現下降趨勢.相反，ρ1=10 Ωm,ρ2=1 Ωm低于圍巖電阻率，電流向異常體聚集，所以在x=[0,25000 m]范圍內測點處的電流發散，視電阻率升高.

圖14 DTM-1模型視電阻率與相位，測試頻率為0.01 Hz，Curl-Curl system代表了自適應有限元結果(Ren et al., 2013),Air-Earth decomposition代表了本文計算的結果Fig.14 Comparison of calculated apparent resistivity and phase on the DTM-1 model among our solution (denoted by Air-Earth decomposition) and the reference solution (Curl-Curl system, Ren et al., 2013) at a frequency of 0.01 Hz

圖15 DTM-1模型電磁響應相對誤差(高精度自適應有限元計算結果為參考值(Ren et al., 2013))，測試頻率為0.01 HzFig.15 Relative error of our apparent resistivity and phase by comparing to the reference solution (Curl-Curl system,Ren et al., 2013) for the DTM-1 model at a frequency of 0.01 Hz

同樣，我們也在該工作頻率下測試了兩個系統的特征值分布，如圖16和圖17所示,DMT-1模型的特征值分布特征與3D-1模型以及梯形山模型的測試結果相似，空氣地球分離公式系統矩陣的特征值更多地分布在復數平面左半部分距離原點更遠處.我們在這個基礎上選擇了從0.001 Hz，0.01 Hz，…，100 Hz共六個對數等間隔分布的頻點進行了條件數的測試，結果如圖18所示.根據圖18可以得知，混合系統對雙旋度系統的病態性做出了改進，在1 Hz的工作頻率下，條件數的改善最明顯，可達3個數量級.

圖16 DTM-1模型上地空分離公式和雙旋度電場公式對應的前1000個最大特征值分布的對比分析，測試頻率為0.01 HzFig.16 The distribution of 1000 maximum eigenvalues for two methods (one using the Air-Earth decomposition,and one using the Curl-Curl system) on the DTM-1 model at a frequency of 0.01 Hz

圖17 DTM-1模型上地空分離公式和雙旋度電場公式對應的后1000個最小特征值分布的對比分析，測試頻率為0.01 HzFig.17 The distribution of 1000 minimum eigenvalues for two methods (one using the Air-Earth decomposition,and one using the Curl-Curl system) on the DTM-1 model at a frequency of 0.01 Hz

圖18 DTM-1模型地空分離公式和雙旋度電場公式的條件數對比Fig.18 Comparison of condition numbers of the Air-Earth decomposition and the Curl-Curl system on the DTM-1 model at different frequencies

3 結論

我們推導了一種新的基于空氣地球分離策略的大地電磁混合計算公式，在空氣中采用Laplace算子取代雙旋度算子，將雙旋度系統分解為Helmholtz-Curl耦合系統，有效改善了系統矩陣在低頻條件下的病態性特征.

模型測試驗證了新系統的正確性與精度.在相同網格的條件下，本文算法的自由度數低于電場雙旋度方程的自由度數.在計算速度上，空氣地球分離公式稍高于電場雙旋度公式.特征值分布與條件數分布表明新系統具有更好的特征值分布與更低的條件數，在數值上更加穩定，對雙旋度系統的病態性做出了有效的改進，因而新系統應更適用于低頻條件下的大地電磁正演工作.

本研究基于空氣地球分離思路推導了一種具有良好條件數的混合公式，為尋求大地電磁最優化控制方程開啟了新的求解思路.基于新思路，我們期待更多、性能更佳的大地電磁混合公式.