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

秦策， 劉幸飛， 王緒本， 孫衛(wèi)斌， 趙寧*

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

0 引言

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

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

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

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

1 三維正演算法

1.1 控制方程

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

(1)

(2)

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

(3)

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

n×E=n×E0,

(4)

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

圖1 非結(jié)構(gòu)化六面體單元上的矢量形函數(shù)定義Fig.1 The definition of vector shape function on unstructured hexahedral element

1.2 自適應(yīng)矢量有限單元法

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

(5)

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

(6)

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

(7)

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

(8)

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

Ke=s,

(9)

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

(10)

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

(11)

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

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

2 三維反演算法

2.1 目標(biāo)函數(shù)

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

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

(12)

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

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

mk+1=mk+αkpk,

(13)

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

2.2 L-BFGS方法

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

sk=mk+1-mk,

(14)

yk=gk+1-gk,

(15)

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

pk=-Bkgk.

(16)

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

在計算步長αk時，目標(biāo)函數(shù)應(yīng)獲得足夠的下降，同時計算量也不能太大.理想情況下，步長αk應(yīng)是一元函數(shù)

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

(17)

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

(18)

和曲率條件

(19)

其中c1、c2為常數(shù)，且0

2.3 網(wǎng)格分離策略

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

圖2 網(wǎng)格分離策略示意圖Fig.2 The schematic diagram of mesh separation strategy

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

3 數(shù)值算例

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

3.1 正演算例

3.1.1 DTM1模型

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

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

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

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

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

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

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

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

3.1.2 地形模型

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

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

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

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

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

3.2 反演算例

3.2.1 簡單三維模型

我們首先通過對一個簡單模型進(jìn)行反演來驗證算法的效率.此模型的背景電阻率為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個頻率(對數(shù)等間隔地分布在10-3～10 Hz范圍內(nèi))的阻抗張量響應(yīng)，并在數(shù)據(jù)中加入2%的高斯噪聲，對數(shù)據(jù)進(jìn)行了反演.初始數(shù)據(jù)擬合差為16.7，經(jīng)過18次迭代下降至0.97，反演結(jié)果如圖11所示.從圖中可以看到，四個塊狀異常體的電阻率和形態(tài)都被反演出來，且與真實模型吻合良好.驗證了本文反演算法的收斂速度和反演程序的正確性.

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

3.2.2 三維地形模型

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

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

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

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

圖13 三維地形模型交錯網(wǎng)格剖分示意圖(a) 水平網(wǎng)格尺寸500 m； (b) 水平網(wǎng)格尺寸1000 m； (c) 水平網(wǎng)格尺寸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 三維地形模型交錯網(wǎng)格反演結(jié)果Fig.14 Inversion results of 3D topography model using staggered grids

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

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

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

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

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

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

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

圖18 包含地形的反演過程最后一次迭代中的正演網(wǎng)格(a) 反演網(wǎng)格； (b) 0.1 Hz時的正演網(wǎng)格； (c) 10 Hz時的正演網(wǎng)格.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 結(jié)論

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

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

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

附錄A 復(fù)系數(shù)線性方程組求解算法

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

(A1)

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

(A2)

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

(A3)

Py=c,

(A4)

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

Bz=d.

(A5)

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

附錄B 面向目標(biāo)后驗誤差估計方法

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

(B1)

(B2)

(B3)

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

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

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

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

(B4)

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

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

附圖B2 非協(xié)調(diào)網(wǎng)格示意圖(a) 由相鄰網(wǎng)格加密次數(shù)不同產(chǎn)生的非協(xié)調(diào)網(wǎng)格； (b) (a) 中左側(cè)4個網(wǎng)格與右側(cè)網(wǎng)格相交界面.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.

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

(B5)

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

(B6)