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

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

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 有色資源與地質(zhì)災(zāi)害探查湖南省重點(diǎn)實(shí)驗(yàn)室, 長(zhǎng)沙 410083 4 自然資源部覆蓋區(qū)深部資源勘查工程技術(shù)創(chuàng)新中心， 合肥 230001 5 成都大學(xué)建筑與土木工程學(xué)院, 成都 610106

0 引言

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

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

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

空解現(xiàn)象產(chǎn)生的本質(zhì)原因?yàn)樵诳諝夂偷厍蛑胁捎昧送瑯拥目刂品匠?，由于空氣和地下具有完全不同的電性特?電導(dǎo)率對(duì)比度為無(wú)窮大)，采用同樣的控制方程必然導(dǎo)致病態(tài)特征.為解決此問(wèn)題，本文提出了一種新的空氣與地球分離混合公式.解決思路為：首先，將求解區(qū)域分解為空氣、地下兩個(gè)獨(dú)立的區(qū)域，分別采用不同的控制方程；然后，利用地空分界面上電磁場(chǎng)連續(xù)性條件，進(jìn)行界面兩側(cè)的耦合與信息交換，保證電磁場(chǎng)解的唯一性.基于該解決思路，空氣層中磁場(chǎng)H分解為電矢量勢(shì)T和磁標(biāo)量勢(shì)φ，利用Coulomb規(guī)范條件，使得矢量勢(shì)T和標(biāo)量勢(shì)φ滿足具有Laplace結(jié)構(gòu)的勢(shì)場(chǎng)方程.地下區(qū)域中，采用雙旋度結(jié)構(gòu)的電場(chǎng)方程.在地下區(qū)域采用電場(chǎng)雙旋度方程的原因在于：地下介質(zhì)的電導(dǎo)率大于零，這使得雙旋度系統(tǒng)在地下區(qū)域避免了零空間問(wèn)題且有良好的穩(wěn)定性.

最后，我們計(jì)算了國(guó)際標(biāo)準(zhǔn)測(cè)試模型COMMEMI 3D-1、DTM-1與梯形山模型的大地電磁響應(yīng)，將計(jì)算結(jié)果與其他算法進(jìn)行了對(duì)比，來(lái)驗(yàn)證本文算法的正確性，以及處理起伏地表問(wèn)題的能力.我們還進(jìn)行了系統(tǒng)條件數(shù)與特征值的計(jì)算和對(duì)比分析，來(lái)檢驗(yàn)本文算法對(duì)電場(chǎng)雙旋度系統(tǒng)病態(tài)性的改進(jìn).

1 節(jié)點(diǎn)矢量混合有限元系統(tǒng)

利用大地電磁的場(chǎng)源位于求解區(qū)域外的一般性假設(shè)，在求解區(qū)域內(nèi)(見(jiàn)圖1)，大地電磁問(wèn)題的電磁場(chǎng)響應(yīng)滿足頻率域Maxwell方程組(選取時(shí)諧因子為e-iω t)：

(1)

其中，E為電場(chǎng)矢量，H為磁場(chǎng)矢量，i為純虛數(shù)單位，ω為角頻率，μ為磁導(dǎo)率，

(2)

其中，J為電流密度，σ為電導(dǎo)率，ε為介電常數(shù)，

(3)

(4)

圖1 MT地電模型示意圖Ω0表示空氣區(qū)域，Ω1表示地下區(qū)域，Γ表示地空邊界，Γ0表示空氣外部邊界，Γ1表示地球外部邊界，n0為空氣區(qū)域外邊界的外法向向量，n1為地球區(qū)域外邊界的外法向向量，n01為地空邊界上空氣一側(cè)的外法向向量，方向由空氣指向地球，n10為地空邊界上地球一側(cè)的外法向向量，方向由地球指向空氣，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 空氣層電磁場(chǎng)積分弱形式

(5)

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

(6)

其中H0表示空氣中的磁場(chǎng).引入標(biāo)量函數(shù)φ的負(fù)梯度來(lái)表示公式(6)中的矢量部分，可得到：

(7)

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

即：

(8)

將公式(8)中的矢量勢(shì)T的雙旋度結(jié)構(gòu)進(jìn)行化簡(jiǎn)，強(qiáng)加Coulomb規(guī)范條件來(lái)唯一化矢量勢(shì)T：

(9)

(10)

(11)

綜上，空氣區(qū)域的控制方程為

(12)

為了求解公式(12)，采用Galerkin節(jié)點(diǎn)有限元法(Jin，2002)，先假設(shè)余量rT為

(13)

余量rφ為

(14)

采用線性插值基函數(shù)，四面體單元中任意點(diǎn)處的矢量勢(shì)T以及標(biāo)量勢(shì)φ可由節(jié)點(diǎn)基函數(shù)表示為

(15)

其中Txi，Tyi，Tzi分別為矢量勢(shì)T的三個(gè)分量，Li為節(jié)點(diǎn)基函數(shù)，在空氣區(qū)域中滿足以下條件：

?Ω0Li·rTdv=0,

?Ω0Li·rφdv=0.

(16)

利用Green恒等式：

(17)

(18)

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

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

(19)

寫成分量形式為

(20)

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

1.2 地球電磁場(chǎng)積分弱形式

地下區(qū)域中，電場(chǎng)滿足雙旋度方程：

(21)

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

(22)

采用棱邊基函數(shù)，四面體單元中的電場(chǎng)可表示為

(23)

其中E1i為對(duì)應(yīng)棱邊上的電場(chǎng)切向分量，Ni為棱邊基函數(shù)，計(jì)算區(qū)域內(nèi)滿足以下條件：

?Ω0Ni·REdv=0.

(24)

(25)

其中a,b分別為矢量函數(shù)，則

(26)

(27)

(28)

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

1.3 空氣-地球耦合系統(tǒng)

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

(29a)

(29b)

(29c)

根據(jù)矢量公式：

(30)

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

(31)

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

(32)

(33)

在地空分界面上施加如下電磁場(chǎng)連續(xù)性條件:

(34)

(35a)

(35b)

(35c)

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

(36)

其中：

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

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

2 數(shù)值試驗(yàn)

本小節(jié)所有計(jì)算均在CPU主頻為2.30 GHz，RAM超過(guò)200 GB的小型AMAX工作站上完成.

2.1 COMMIME 3D-1模型

3D-1模型是國(guó)際標(biāo)準(zhǔn)測(cè)試COMMEMI項(xiàng)目(Zhdanov et al., 1997)的一個(gè)模型.它是均勻半空間中嵌套一個(gè)導(dǎo)電棱柱，其頂埋深為250 m，背景電阻率為100 Ωm，棱柱電阻率為0.5 Ωm，尺寸為1000 m×2000 m×2000 m，空氣電阻率為1016Ωm,計(jì)算區(qū)域尺寸為70 km×70 km×140 km.

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

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

圖2 3D-1模型視電阻率與相位，測(cè)試頻率為0.1 Hz, Curl-Curl system代表了自適應(yīng)有限元結(jié)果(Ren et al., 2013),Air-Earth decomposition 代表了本文計(jì)算的結(jié)果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模型視電阻率與相位相對(duì)誤差(高精度自適應(yīng)有限元計(jì)算結(jié)果為參考值(Ren et al., 2013))，測(cè)試頻率為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表明兩個(gè)系統(tǒng)的最大特征值分布具有相似性.圖5表明地空分離公式系統(tǒng)有更多的非零或模長(zhǎng)更大的最小特征值.當(dāng)一個(gè)系統(tǒng)的特征值全部位于復(fù)數(shù)平面的左半部分時(shí)，系統(tǒng)是收斂的.與此相反，當(dāng)特征值全部位于復(fù)數(shù)平面的右半部分時(shí)，系統(tǒng)則是發(fā)散的.從圖4與圖5可以直觀的看出，空氣地球分離公式的特征值更多地集中在復(fù)數(shù)平面的左半部分且距離原點(diǎn)更遠(yuǎn)，因而空氣地球分離公式相對(duì)于電場(chǎng)雙旋度公式有更好的穩(wěn)定性.根據(jù)條件數(shù)的定義(最大特征值與最小特征值之比)，地空分離系統(tǒng)擁有更小的條件數(shù).為了定量研究地空分離系統(tǒng)的條件數(shù)，計(jì)算了不同頻率下的條件數(shù)(從0.001 Hz，0.01 Hz，…，10 Hz共六個(gè)對(duì)數(shù)等間隔分布的頻點(diǎn)).測(cè)試結(jié)果如圖6所示，結(jié)果表明混合系統(tǒng)的條件數(shù)相對(duì)于雙旋度系統(tǒng)均降低，下降可達(dá)3個(gè)數(shù)量級(jí).

圖4 3D-1模型上地空分離公式和雙旋度電場(chǎng)公式對(duì)應(yīng)的前1000個(gè)最大特征值分布的對(duì)比分析，測(cè)試頻率為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模型上地空分離公式和雙旋度電場(chǎng)公式對(duì)應(yīng)的后1000個(gè)最小特征值分布的對(duì)比分析，測(cè)試頻率為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模型地空分離公式和雙旋度電場(chǎng)公式的條件數(shù)對(duì)比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 梯形山模型

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

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

圖8 梯形山模型視電阻率與相位，測(cè)試頻率為2 Hz，Curl-Curl system代表了自適應(yīng)有限元結(jié)果(Ren et al., 2013),Air-Earth decomposition代表了本文計(jì)算的結(jié)果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

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

為了研究混合系統(tǒng)與雙旋度系統(tǒng)的數(shù)值特征，我們計(jì)算了特征值分布，測(cè)試結(jié)果如圖10和圖11.與3D-1模型類似，圖10表明兩個(gè)系統(tǒng)具有相似的最大特征值分布，圖11表明地空分離系統(tǒng)擁有模長(zhǎng)更大的最小特征值，且地空分離系統(tǒng)的特征值更加集中在復(fù)數(shù)平面的左半部分且距離原點(diǎn)更遠(yuǎn).

圖10 梯形山模型上地空分離公式和雙旋度電場(chǎng)公式對(duì)應(yīng)的前1000個(gè)最大特征值分布的對(duì)比分析，測(cè)試頻率為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 梯形山模型上地空分離公式和雙旋度電場(chǎng)公式對(duì)應(yīng)的后1000個(gè)最小特征值分布的對(duì)比分析，測(cè)試頻率為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

我們?cè)谶@個(gè)模型中選取0.01 Hz，0.1 Hz，…，100 Hz共五個(gè)頻點(diǎn)測(cè)試了系統(tǒng)的條件數(shù).如圖12所示，模型的條件數(shù)分布與COMMEMI 3D-1模型的條件數(shù)分布特征一致，條件數(shù)有1～3個(gè)數(shù)量級(jí)的改善，在10 Hz頻率下兩系統(tǒng)條件數(shù)差異最大.

圖12 梯形山模型地空分離公式和雙旋度電場(chǎng)公式的條件數(shù)對(duì)比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模型

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

圖13 DTM-1模型結(jié)構(gòu)Fig.13 Structure of DTM-1 model

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

圖14 DTM-1模型視電阻率與相位，測(cè)試頻率為0.01 Hz，Curl-Curl system代表了自適應(yīng)有限元結(jié)果(Ren et al., 2013),Air-Earth decomposition代表了本文計(jì)算的結(jié)果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模型電磁響應(yīng)相對(duì)誤差(高精度自適應(yīng)有限元計(jì)算結(jié)果為參考值(Ren et al., 2013))，測(cè)試頻率為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

同樣，我們也在該工作頻率下測(cè)試了兩個(gè)系統(tǒng)的特征值分布，如圖16和圖17所示,DMT-1模型的特征值分布特征與3D-1模型以及梯形山模型的測(cè)試結(jié)果相似，空氣地球分離公式系統(tǒng)矩陣的特征值更多地分布在復(fù)數(shù)平面左半部分距離原點(diǎn)更遠(yuǎn)處.我們?cè)谶@個(gè)基礎(chǔ)上選擇了從0.001 Hz，0.01 Hz，…，100 Hz共六個(gè)對(duì)數(shù)等間隔分布的頻點(diǎn)進(jìn)行了條件數(shù)的測(cè)試，結(jié)果如圖18所示.根據(jù)圖18可以得知，混合系統(tǒng)對(duì)雙旋度系統(tǒng)的病態(tài)性做出了改進(jìn)，在1 Hz的工作頻率下，條件數(shù)的改善最明顯，可達(dá)3個(gè)數(shù)量級(jí).

圖16 DTM-1模型上地空分離公式和雙旋度電場(chǎng)公式對(duì)應(yīng)的前1000個(gè)最大特征值分布的對(duì)比分析，測(cè)試頻率為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模型上地空分離公式和雙旋度電場(chǎng)公式對(duì)應(yīng)的后1000個(gè)最小特征值分布的對(duì)比分析，測(cè)試頻率為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模型地空分離公式和雙旋度電場(chǎng)公式的條件數(shù)對(duì)比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 結(jié)論

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

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

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