基于自適應有限元正演的大地電磁法三維反演算法研究

秦策， 劉幸飛， 王緒本， 孫衛斌， 趙寧*

1 河南理工大學物理與電子信息學院， 河南焦作 454000 2 成都理工大學地球物理學院， 成都 610059 3 中國石油集團東方地球物理公司綜合物化探處， 河北涿州 072750

0 引言

大地電磁法是實際工作中應用非常廣泛的一種電磁勘探方法，在深部電性結構、礦產、油氣和地熱資源勘探等領域發揮了巨大的作用(趙國澤等, 2007).相對于一、二維反演，三維在消除假異常等方面具有很大優勢(Siripunvaraporn, 2012).近年來，隨著計算機硬件能力和方法理論的進步，同時三維數據采集也已逐步普及，三維反演應用越來越廣泛(Dong et al., 2020; 楊文采等, 2020).反演需要使用正演來計算電磁場響應和靈敏度矩陣，三維正演是三維反演的基礎，其計算精度越高，正演響應和靈敏度矩陣的精度也越高，正演計算速度越快，反演的效率也越高.在眾多三維正演方法中，交錯網格有限差分法有著理論簡單、易于實現、計算量小等優點，是目前在三維反演中使用最多的正演方法(Siripunvaraporn et al., 2005; 胡祥云等, 2012; Kelbert et al., 2014; 董浩等, 2014; 秦策等, 2017; 阮帥等, 2020).但是，結構化六面體網格只能對地形進行近似，影響了對地形的處理效果.另外，反演中采用同套正演和反演網格，且網格不能自適應變化，這帶來了反演網格設置的問題.若網格過密，則增加了反演的非唯一性；若網格過稀，則無法保證正演響應和靈敏度矩陣的計算精度，影響了反演的可靠性.在實際工作中，通常會使用不同的反演網格做多次反演，大大增加了數據處理的工作量和難度.

最近十年來，有限元法在電磁法的三維正演中得到了迅速的發展.非結構化網格(四面體和形變六面體)能夠很好地模擬復雜地形和異常體.另外，自適應有限元法能夠對網格進行自適應加密，相比全局網格加密可以更高效地提高正演響應的精度(Ren et al., 2013; Grayver and Bürg, 2014; 殷長春等, 2017; 趙寧等, 2019).很多學者基于有限元法實現了大地電磁法的三維反演(Grayver, 2015; Usui, 2015; Kordy et al., 2016b; Jahandari and Farquharson, 2017)，也有學者將有限元法應用在了其他電磁法的三維反演中(Liu et al., 2019; Chen et al., 2020)，這些研究取得了非常好的應用效果.更進一步地，若能夠將自適應有限元正演方法應用在三維反演中，可以預見能夠得到更好的效果.我們認為主要的困難是如何處理網格自適應加密的問題.一方面，大地電磁法的觀測頻率范圍非常寬，因此希望能夠自適應地對正演網格進行加密，不同頻率使用不同的正演網格，以提高計算精度；另一方面，很多反演方法要求反演網格不能改變，反演過程中只能有一套網格，這與正演網格的自適應加密產生了矛盾.為解決此問題，我們設計了正演網格和反演網格相分離的策略，即保持反演網格不變，正演網格從反演網格出發進行自適應加密，以確保正演響應和偏導數計算的精確性，同時避免反演網格的過度參數化.

另外，有限元正演中使用的網格類型也是很多學者關心的問題.理論上，任何非結構化網格均可以模擬復雜異常體和地形.在實際應用中，由于四面體網格更容易生成，在電磁法領域應用最為廣泛(Ren et al., 2013; Usui, 2015; Jahandari and Farquharson, 2017).也有一些學者利用非結構化的六面體網格得到了較好的結果(Grayver, 2015; Kordy et al., 2016a).至于哪一種網格更優，目前還沒有定論.Cifuentes和Kalbag(1992)在結構問題的模擬結果中表明六面體網格的精度和穩定性相對四面體網格較高，Tadepalli等(2011)在生物力學中的模擬也得到了類似的結果，而在Ramos和Sim?es(2006)的股骨模擬中，四面體網格表現出了更大的優勢.我們認為該問題與具體的領域相關，有待進一步研究.在本文的反演算法中，我們使用的網格分離策略需要保證加密前后的網格具有層級關系，因此選擇使用非結構化六面體網格并利用八叉樹方式對其進行加密.

本文的第1節介紹了三維自適應有限元正演方法，包括線性方程組的求解算法和面向目標的后驗誤差估計方法.第2節給出了本文中使用的L-BFGS反演算法以及網格分離策略.第3節中通過兩個正演算例驗證了正演算法的正確性和對地形的模擬精度，并重點對一個三維地形模型的正演數據進行了不同方法的反演，驗證了本文網格分離策略的優勢和處理地形問題的有效性.

1 三維正演算法

1.1 控制方程

取時諧因子為eiω t，大地電磁法中電場和磁場所滿足的偏微分方程為

(1)

(2)

其中σ是介質的電導率，ω是角頻率，μ是磁導率，其值為4π×10-7.對(1)式兩邊求旋度，并將(2)式帶入，可得介質中電場所滿足的二階矢量亥姆霍茲方程:

(3)

在邊界上施加Dirichlet邊界條件，即

n×E=n×E0,

(4)

其中n是邊界網格面上的外法向向量，E0是邊界上一維介質的電場響應，可以通過解析方法計算(Cagniard, 1953).

圖1 非結構化六面體單元上的矢量形函數定義Fig.1 The definition of vector shape function on unstructured hexahedral element

1.2 自適應矢量有限單元法

使用數值方法求解上述偏微分方程，需要將求解區域進行離散化.為能夠模擬起伏地形和復雜異常體，本文使用非結構化的六面體單元.同時，為滿足電場的連續性條件，將形函數定義在單元的棱邊上(Nédélec, 1986)，如圖1所示.在(3)式兩邊同時點乘矢量形函數V，并對全區域積分:

(5)

其中V∈H(curl;Ω).H(curl;Ω)是Hilbert空間上的旋度平方可積函數空間，其定義為

(6)

根據矢量恒等式和分部積分原理，式(5)可轉換為

(7)

式(7)即為與式(3)所等價的泛函形式.將式(7)在區域中的每個單元進行積分并求和，可得

(8)

各單元上的積分可以由高斯數值積分方法進行計算(Venkateshan and Swaminathan, 2014).最終可得如下線性方程組

Ke=s,

(9)

任意點P處的電場值可由公式

(10)

計算.根據法拉第定律，點P處的磁場可表示為

(11)

以上兩式中，ei是點P所在單元中第i個棱邊上形函數的插值系數，Vi(P)是點P處第i個棱邊上的形函數值.進一步可由電磁場值計算得到觀測點處的阻抗張量響應.

在有限元法中，除了單元上形函數的階數，數值解的精度基本取決于網格的大小.在三維情況下，全局網格加密會急劇地增加問題規模，是非常不經濟的.本文使用基于面向目標后驗誤差方法的自適應網格加密來更有效地改善數值解的精度.附錄B中給出了后驗誤差估計理論和自適應網格加密方法.

2 三維反演算法

2.1 目標函數

在反演中，我們的目標是獲取一個模型，其正演響應能夠在一定程度上擬合觀測數據，同時又符合實際的地質規律.根據正則化反演理論(Tong et al., 2018)，反演目標函數可表示為

φ(m)=(d-F)TV-1(d-F)+λmTLTLm,

(12)

上式中第一項為數據擬合項，衡量著數據擬合程度，第二項為模型約束項，λ為正則化因子，控制著模型約束項在目標函數中的比重，d為觀測數據向量，m為待反演模型向量，F為模型的正演響應向量，V為數據方差矩陣，L為拉普拉斯算子的離散形式.

目前常用的反演方法大多派生自牛頓法，通過迭代法求目標函數的極小值，迭代形式為

mk+1=mk+αkpk,

(13)

其中pk為搜索方向，控制著模型的修正方向，αk為步長，控制著模型在搜索方向上的修正大小.在眾多反演方法中，非線性共軛梯度法(NLCG)(Rodi and Mackie, 2001)和有限內存的BFGS方法(L-BFGS)(Byrd et al., 1994)只需計算目標函數值及其梯度，計算量較小，比較適合三維反演.其中，L-BFGS相對NLCG具有更快的收斂速度，同時確定步長所需的搜索次數也較少(秦策, 2018).經綜合考慮，本文在反演中使用L-BFGS方法.

2.2 L-BFGS方法

在L-BFGS方法中，只需存儲最近l次迭代中模型的修正量sk和梯度的修正量yk，其中

sk=mk+1-mk,

(14)

yk=gk+1-gk,

(15)

因此僅需要保存2l個向量，占用的內存空間較少.每一次反演迭代中，使用{s1,s2,…,sk}和{y1,y2,…,yk}近似計算Hessian矩陣，記近似Hessian矩陣的逆為Bk，搜索方向pk可表示為

pk=-Bkgk.

(16)

Bk的計算方法可參考Nocedal和Wright(2006).

在計算步長αk時，目標函數應獲得足夠的下降，同時計算量也不能太大.理想情況下，步長αk應是一元函數

f(αk)=φ(mk+αkpk)

(17)

的全局極小值點.但是在實際中，計算其局部極小也需要多次計算目標函數.更可行的策略是通過不精確的線搜索來獲得滿足條件的步長αk，既能使目標函數獲得充分的下降，又花費盡量小的計算代價.最常用的條件是Wolfe條件(Nocedal and Wright, 2006)，即充分下降條件

(18)

和曲率條件

(19)

其中c1、c2為常數，且0

2.3 網格分離策略

在自適應有限元方法中，正演網格會自適應地根據后驗誤差估計值進行加密，不同頻率得到的最終網格也不相同.同時，L-BFGS算法也要求反演過程中網格不發生變化.為解決此問題，我們使用將正演網格與反演網格相分離的策略.

圖2 網格分離策略示意圖Fig.2 The schematic diagram of mesh separation strategy

算法1反演過程中梯度計算流程1:確定反演網格剖分T2:令gk=03:forf=1,…,Nfdo4: 令Tf0=T5: fori=1,…,Nmaxdo6: 對Tfi-1進行自適應加密,得到Tfi7: 令Tf=Tfi8: end for9: 在Tf上計算梯度gfk10: 令gk=gk+Dgfk11:end for

3 數值算例

基于本文的正演和反演算法，我們使用C++程序設計語言開發了正反演程序.程序設計過程中使用了開源的有限元程序庫deal.II(Bangerth et al., 2007; Arndt et al., 2021).本文的所有算例均在河南理工大學高性能計算中心的集群上完成，計算節點配備了Intel Xeon E5 2680 V4 CPU，包含14個處理器核心，主頻為2.4 GHz，內存128 GB.

3.1 正演算例

3.1.1 DTM1模型

為驗證本文正演算法的正確性，我們使用標準模型DTM1(Dublin Test Model 1)(Miensopust et al., 2013)對程序進行測試.該模型的背景電阻率為100 Ωm，其中包含了三個電阻率差異非常大的塊狀異常體，異常體的位置、尺寸和電阻率如圖3所示.

圖3 DTM1模型示意圖，圖片修改自Miensopust等(2013)Fig.3 Sketch of DTM1 model, figure revised from Miensopust et al.(2013)

我們計算了10-4Hz至10 Hz范圍內的21個頻率的響應.正演過程中，初始網格中單元數為6498，自由度數為21640，每次加密約10%的網格，經過10次自適應加密，表1中給出了自適應加密過程中自由度的變化和計算時間.圖4是全局網格加密和自適應加密過程中歸一化后驗誤差的變化趨勢，可以看到隨著網格的加密，誤差穩步下降.自適應網格加密時誤差下降的速度更快，意味著可以用較小的計算量獲得更高的計算精度.圖5展示了頻率分別是10 Hz和0.01 Hz時的自適應網格加密結果，可以看到網格得到了充分的加密，頻率為10 Hz時淺部網格加密的較多，而頻率為0.01 Hz時深部的網格更加稠密.這與我們對電磁波傳播規律的認識是一致的，高頻時電磁波衰減較快，因此淺部的網格加密較多；低頻時電磁波衰減慢，深部的網格也需要加密.另外，由于電場穿過介質分界面是不連續的，所以網格在電阻率變化劇烈的地方加密的較多.從網格的自適應加密結果來看，本文使用的后驗誤差估計方法是有效的，能夠較好地反映數值解的誤差分布.同時，由于我們使用了面向目標的后驗誤差估計方法，觀測點附近的網格都得到了加密，可以大幅度提高觀測點處響應的精度.

表1 DTM1模型自適應加密過程中自由度數目及計算時間Table 1 Number of DoFs and computational time for DTM1 model using adaptive mesh refinement

圖4 DTM1模型10 Hz和0.01 Hz自適應網格加密歸一化誤差收斂速度Fig.4 Normalized estimated errors using adaptive mesh refinement for the DTM1 model for frequency 10 Hz and 0.01 Hz

圖5 DTM1模型第10次自適應加密結果(a)10 Hz；(b)0.01 Hz.Fig.5 Adaptive mesh refinement result of DTM1 model

已有很多學者對此模型進行了模擬(Miensopust et al., 2013).(0 km,0 km)位于三個異常體交界處的上方，其響應最為奇異.圖6中展示了本文的計算結果與Mackie等(1993)的有限差分結果、Nam等(2007)等的有限元結果的對比.從圖中可以看出，Zxy和Zyx的視電阻率和相位響應吻合良好，這證明了本文所采用方法的正確性.

圖6 DTM1模型(0 km, 0 km)處的響應Fig.6 The response of DTM1 model at (0 km, 0 km)

3.1.2 地形模型

非結構化網格的優勢之一是可以精確地模擬起伏地形.一般來說，四面體單元的適應性最強，可以模擬任意復雜的地形.在本文中，我們使用非結構化六面體單元，通過對單元進行形變也可模擬起伏地形.通過對一個二維山脊地形(Wannamaker et al., 1986)進行模擬來驗證程序對地形處理的正確性.該模型的背景電阻率為100 Ωm，模型的中間位置包含一平臺狀的地形，如圖7所示.計算了頻率為2 Hz時的響應.圖8中展示了初始網格和經過5次自適應加密得到的最終網格.初始網格非常稀疏，最終網格中，網格得到了充分加密.與DTM1模型類似，測點附近網格的加密次數更多，解的精度得到了保證.

圖7 二維山脊地形模型示意圖Fig.7 The diagram of 2D ridge topography model

圖8 二維山脊地形模型網格自適應加密結果(a) 初始網格； (b) 最終網格.Fig.8 Adaptive mesh refinement result of 2D ridge topography model(a) Initial mesh； (b) Final mesh.

使用開源的二維自適應有限元正演程序MARE2DEM(Key, 2016)對此模型進行了模擬，并和我們的計算結果進行對比，如圖9所示.可以看到，兩者吻合良好，這驗證了我們使用的非結構化六面體網格在處理地形問題時的正確性.另外，從響應中也可發現，地形的影響非常大.因此，在地形起伏情況下，其影響是不能忽略的，必須加以處理，否則會對反演結果產生嚴重的干擾.我們將在反演算例部分對地形的處理方法進行討論.

圖9 二維山脊地形模型響應(頻率為2 Hz)Fig.9 Response of 2D ridge topography model (frequency is 2 Hz)

3.2 反演算例

3.2.1 簡單三維模型

我們首先通過對一個簡單模型進行反演來驗證算法的效率.此模型的背景電阻率為100 Ωm，其中包含4個塊狀異常體(2個高阻異常體和2個低阻異常體)，電阻率分別為10 Ωm和1000 Ωm，異常體的尺寸為10 km×10 km×5 km，相鄰異常體的間距為5 km，異常體的埋深為2.5 km，如圖10所示.

圖10 簡單三維模型示意圖Fig.10 The diagram of simple 3D model

計算了21個頻率(對數等間隔地分布在10-3～10 Hz范圍內)的阻抗張量響應，并在數據中加入2%的高斯噪聲，對數據進行了反演.初始數據擬合差為16.7，經過18次迭代下降至0.97，反演結果如圖11所示.從圖中可以看到，四個塊狀異常體的電阻率和形態都被反演出來，且與真實模型吻合良好.驗證了本文反演算法的收斂速度和反演程序的正確性.

圖11 簡單三維模型反演結果Fig.11 Inversion result of simple 3D model

3.2.2 三維地形模型

在前文的正演地形算例中，我們看到，大地電磁響應受地形影響非常嚴重，因此在處理實測數據時，必須對地形進行處理.一些學者提出了地形校正的方法(Nam et al., 2008)，即根據地形的響應特征，將地形的干擾從觀測數據中分離出去，再對不包含地形影響的數據進行反演.另外一種思路是不對數據進行任何處理，將地形包含在初始模型中進行反演.已有研究表明，即使使用臺階狀的網格近似地形，也可以提高對地形的模擬精度，一定程度上減弱地形對反演結果的影響 (董浩等, 2014; 余輝等, 2019; 顧觀文和李桐林, 2020).

我們建立了如圖12所示的地形模型，通過對該地形模型的正演合成數據進行反演來討論地形數據的處理方法.模型的背景電阻率為100 Ωm，包含2個塊狀異常體，電阻率分別為10 Ωm和1000 Ωm，尺寸為10 km×10 km×5 km.兩個塊狀異常體上方各有一個平臺狀的正地形，高度為2 km.觀測區域范圍為22.5 km×37.5 km，覆蓋了整個地形區域，測點間距2.5 km，如圖12b中十字符號所示.觀測頻率共21個，對數等間隔地分布在10-3Hz至10 Hz范圍內.

圖12 三維地形模型示意圖Fig.12 The diagram of 3D topography model

使用本文的自適應正演方法對此模型進行計算，并在阻抗張量響應中加入了2%的高斯噪聲，得到地形模型的合成數據.我們首先使用近年來在學術界應用比較廣泛的三維反演程序ModEM(Kelbert et al., 2014)對此數據進行反演.ModEM使用的是交錯網格，對地形的近似程度取決于地形附近網格單元的大小，因此我們使用了三種不同尺度的網格.在地形區域，縱向的網格單元尺寸設置為100 m，橫向網格單元尺寸分別選擇500 m、1000 m和2500 m.如圖13所示，三種尺度的網格都能在一定程度上對地形進行模擬，單元尺寸越小，對地形的近似程度越高，但同時也會導致區域內網格數目急劇增長.圖14展示了使用不同尺寸網格的反演結果，圖15是反演過程中RMS的收斂過程.可以看到，單元尺寸為2500 m時的反演結果很好地恢復了低阻和高阻異常體，但是背景中有虛假異常.經過50次迭代RMS只下降到約2.7，這是由于對地形的近似比較粗糙，不能很好地擬合數據.單元尺寸為1000 m和500 m時對地形的近似比較好，經過約30次迭代數據即得到了很好的擬合，低阻異常體的位置和電阻率都得到了較好的反映.但是，在結果中幾乎看不到高阻異常體.我們推測主要有兩方面的原因，一方面由于大地電磁法本身對高阻體不靈敏，另一方面可能是因為反演參數過多增大了反演的非唯一性.

圖13 三維地形模型交錯網格剖分示意圖(a) 水平網格尺寸500 m； (b) 水平網格尺寸1000 m； (c) 水平網格尺寸2500 m.Fig.13 Staggered grids of 3D topography model(a) Cell size 500 m； (b) Cell size 1000 m； (c) Cell size 2500 m.

圖14 三維地形模型交錯網格反演結果Fig.14 Inversion results of 3D topography model using staggered grids

圖15 三維地形模型交錯網格反演過程RMS收斂曲線Fig.15 Convergence curve of RMS in the inversion process of 3D topography model using staggered grids

上述反演結果表明，使用較粗的網格很難擬合數據，而網格較密時反演非唯一性過強，反演結果并不理想.我們進一步使用本文的方法對地形數據進行反演，目標區域的網格單元尺寸為2500 m，并對地表附近的網格進行了加密和形變以適應地形.反演初始模型為100 Ωm的均勻半空間，同時我們也進行了不包含地形的反演.反演結果如圖16所示.

圖16 三維地形模型反演結果(a) 真實模型； (b) 包含地形的反演結果； (c) 不包含地形的反演結果.Fig.16 Inversion result of 3D topography model(a) The true model； (b) The inversion result with topography； (c) The inversion result without topography.

可以看到，在包含地形的反演中，兩個異常體的尺寸、位置和電阻率都得到了較好的反映，模型背景也較為干凈.相對地，在不包含地形時，反演結果較差，高阻體基本沒有反演出來，反演結果中也出現了許多虛假異常.另外，低阻體的形狀不準確，且位置發生了下移.我們認為主要的原因是正地形產生的低電阻率異常掩蓋了高阻體的響應，導致高阻體未反演出來，同時也增強了低阻體的響應，導致其位置下移.

圖17展示了兩種情況下數據擬合差隨迭代次數的變化趨勢，在包含地形的情況下，經過23次迭代數據擬合差由10.3下降至0.98，而在不包含地形的情況下，經過50次迭代數據擬合差由20.5下降至2.3，且在后20次迭代中幾乎沒有下降，可以預見繼續迭代下去也不會進一步下降.這說明在不包含地形的情況下，很難尋找到一個模型能夠擬合地形數據，反演過程陷入局部極小.反演結果中的虛假異常體可能是迭代過程中為強行擬合地形數據所產生的“噪聲”.

圖17 三維地形模型反演過程RMS收斂曲線Fig.17 Convergence curve of RMS in the inversion process of 3D topography model

在反演初始模型中，我們使用了較為稀疏的網格.在計算梯度過程中，正演網格由反演網格出發自適應加密5次，每次加密約10%的網格.圖18展示了反演網格和頻率為0.1 Hz和10 Hz時包含地形的反演過程中最后一次迭代的正演網格.反演網格中單元數目為14400，經過加密，單元數目分別為442954和500375.從圖中可以看到，兩個頻率的正演網格都得到了充分的加密，且10 Hz的正演網格淺部加密的較多，而0.1 Hz的正演網格深部加密的較多，與前文中DTM1模型的結果是一致的.這也說明了網格分離策略的優勢，不同頻率的正演使用不同的網格，從而保證所有頻率的計算精度.

圖18 包含地形的反演過程最后一次迭代中的正演網格(a) 反演網格； (b) 0.1 Hz時的正演網格； (c) 10 Hz時的正演網格.Fig.18 The forward mesh in the last iteration of the inversion process with topography(a) Inversion mesh； (b) Forward mesh of 0.1 Hz； (c) Forward mesh of 10 Hz.

4 結論

本文基于自適應有限元算法，實現了大地電磁法的三維正演程序，并通過網格分離策略將其應用到了的大地電磁法的三維反演中.在正演中，我們使用了非結構化六面體網格和面向目標的后驗誤差估計方法，計算精度較高且能夠模擬起伏地形，兩個正演算例驗證了處理復雜模型和地形的有效性.反演中，通過使用兩套網格，將反演網格和正演網格分離.加密的正演網格保證了正演響應和靈敏度的計算精度，保證了反演的可靠性，同時較為稀疏的反演網格也不會增多反演參數個數.最后通過對三維地形模型的反演討論了地形數據的處理方法.本文的方法具有一定的通用性，也可用于其他電磁法的三維反演中.

本文還有很多需要改進的地方.首先，在反演過程中，我們沒有進行反演網格的加密.主要有兩方面的原因，一方面，我們使用的L-BFGS方法要求反演過程中反演參數個數不能變化，否則會破壞反演過程中存儲的輔助信息的一致性；另一方面，對于實測數據，我們并不知道需要在哪些位置加密反演網格以提高分辨率.Grayver(2015)使用模型的空間導數來加密反演網格，在合成數據的反演中顯示了良好的效果，可以提高反演結果中特定位置的分辨率，但目前尚不清楚該方法是否適用于復雜的實測數據.另外，本文只對理論模型進行了試算，還未對實測數據進行測試.后續將針對這兩方面進一步開展研究工作.

致謝本文的研究過程中使用了河南理工大學高性能計算中心的計算設備，特此表示感謝.也感謝審稿專家百忙之中審閱我們的論文并提出寶貴建議.

附錄A 復系數線性方程組求解算法

令K=Kr+iKi,e=er+iei,s=sr+isi，式(9)可改寫為

(A1)

為了保證其對稱性，將上式中的第二行兩端同時乘以-1，可得

(A2)

式(A2)中矩陣階數為式(9)中矩陣階數的2倍，與式(9)完全等價.為方便起見，我們將式(A2)中的系數矩陣記為A.構造分塊對角矩陣：

(A3)

Py=c,

(A4)

其中c是任意向量.由于P具有分塊對角特性，上式可以轉換為求解兩次如下方程：

Bz=d.

(A5)

我們使用共軛梯度法來求解式(A5)，并使用輔助空間算法作為預條件(Hiptmair and Xu, 2007).

附錄B 面向目標后驗誤差估計方法

記有限元解為E，單元上的后驗誤差可表示為

(B1)

(B2)

(B3)

其中，he和hf分別是單元e的外接球和面f的外接圓的直徑，nf表示面上的外法向向量，[·]表示相鄰單元交界面處的跳躍量.

給定有限元解E，計算式(B2)需要在每個單元上對偏微分方程的殘差進行積分.計算式(B3)中的跳躍量需要對[·]中的式子分別在相鄰的兩個單元交界面上進行積分，再計算其差值.上述積分可使用高斯數值積分方法來計算，即在每個單元上的8個積分點處計算待積分函數值，再乘以權重并求和(Venkateshan and Swaminathan, 2014).

在大地電磁法中，通常主要關心觀測點所在位置解的精度.我們使用面向目標的自適應加密策略來針對性地提高測點處解的精度(Ren et al., 2013; 殷長春等, 2017).即通過在觀測點處放置一個點源來構造對偶問題，使用對偶問題解的后驗誤差估計值對原后驗誤差估計值進行加權.由于點源的奇異性很強，所以它的后驗誤差估計值能夠有效地區分對觀測點處精度影響較大的單元，使用加權后的誤差估計值指導網格加密可以快速提高觀測點處解的精度.記對偶問題的解為ED，則面向目標的后驗誤差估計值可表示為

ηgo=ηe(E)ηe(ED).

(B4)

由式(B4)得到的后驗誤差估計值可以指導正演計算過程中的局部網格加密.對于六面體網格，通常使用八叉樹的方式對其進行加密，即把一個六面體單元的每條邊都一分為二，連接各邊中點、面中心點和單元中心點，可得到八個單元，如附圖B1所示.

附圖 B1八叉樹局部網格加密示意圖Appendix Fig.B1 The schematic diagram of octree local mesh refinement

附圖B2 非協調網格示意圖(a) 由相鄰網格加密次數不同產生的非協調網格； (b) (a) 中左側4個網格與右側網格相交界面.Appendix Fig.B2 The schematic diagram of non-conforming mesh(a) Non-conforming mesh generated by different refinements of adjacent cells； (b) Intersection of 4 left cells and right cell.

需要注意的是，若相鄰單元的加密次數不同，則會產生懸掛點.如附圖B2a所示，細網格中紅色棱邊與相鄰粗網格中較長的棱邊部分重合，藍色棱邊與相鄰粗網格的面相交，破壞了有限元解的切向連續性，須對其施加約束.附圖B2b是附圖B2a中網格交界面的平面圖，與自由度x1、x2關聯的棱邊的長度是與x0關聯的棱邊的長度的一半，它們之間應滿足的條件為

(B5)

與自由度y0關聯的藍色棱邊與相鄰單元的面相交，y0應等于與其同方向的兩個自由度(y1、y2)的平均值，即

(B6)