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

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

0 引言

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

最近十年來(lái)，有限元法在電磁法的三維正演中得到了迅速的發(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; 殷長(zhǎng)春等, 2017; 趙寧等, 2019).很多學(xué)者基于有限元法實(shí)現(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ù)見(jiàn)能夠得到更好的效果.我們認(rèn)為主要的困難是如何處理網(wǎng)格自適應(yīng)加密的問(wèn)題.一方面，大地電磁法的觀測(cè)頻率范圍非常寬，因此希望能夠自適應(yīng)地對(duì)正演網(wǎng)格進(jìn)行加密，不同頻率使用不同的正演網(wǎng)格，以提高計(jì)算精度；另一方面，很多反演方法要求反演網(wǎng)格不能改變，反演過(guò)程中只能有一套網(wǎng)格，這與正演網(wǎng)格的自適應(yīng)加密產(chǎn)生了矛盾.為解決此問(wèn)題，我們?cè)O(shè)計(jì)了正演網(wǎng)格和反演網(wǎng)格相分離的策略，即保持反演網(wǎng)格不變，正演網(wǎng)格從反演網(wǎng)格出發(fā)進(jìn)行自適應(yīng)加密，以確保正演響應(yīng)和偏導(dǎo)數(shù)計(jì)算的精確性，同時(shí)避免反演網(wǎng)格的過(guò)度參數(shù)化.

另外，有限元正演中使用的網(wǎng)格類型也是很多學(xué)者關(guān)心的問(wèn)題.理論上，任何非結(jié)構(gòu)化網(wǎng)格均可以模擬復(fù)雜異常體和地形.在實(shí)際應(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)，目前還沒(méi)有定論.Cifuentes和Kalbag(1992)在結(jié)構(gòu)問(wèn)題的模擬結(jié)果中表明六面體網(wǎng)格的精度和穩(wěn)定性相對(duì)四面體網(wǎng)格較高，Tadepalli等(2011)在生物力學(xué)中的模擬也得到了類似的結(jié)果，而在Ramos和Sim?es(2006)的股骨模擬中，四面體網(wǎng)格表現(xiàn)出了更大的優(yōu)勢(shì).我們認(rèn)為該問(wèn)題與具體的領(lǐng)域相關(guān)，有待進(jìn)一步研究.在本文的反演算法中，我們使用的網(wǎng)格分離策略需要保證加密前后的網(wǎng)格具有層級(jí)關(guān)系，因此選擇使用非結(jié)構(gòu)化六面體網(wǎng)格并利用八叉樹(shù)方式對(duì)其進(jìn)行加密.

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

1 三維正演算法

1.1 控制方程

取時(shí)諧因子為eiω t，大地電磁法中電場(chǎng)和磁場(chǎng)所滿足的偏微分方程為

(1)

(2)

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

(3)

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

n×E=n×E0,

(4)

其中n是邊界網(wǎng)格面上的外法向向量，E0是邊界上一維介質(zhì)的電場(chǎng)響應(yīng)，可以通過(guò)解析方法計(jì)算(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)化的六面體單元.同時(shí)，為滿足電場(chǎng)的連續(xù)性條件，將形函數(shù)定義在單元的棱邊上(Nédélec, 1986)，如圖1所示.在(3)式兩邊同時(shí)點(diǎn)乘矢量形函數(shù)V，并對(duì)全區(qū)域積分:

(5)

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

(6)

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

(7)

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

(8)

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

Ke=s,

(9)

任意點(diǎn)P處的電場(chǎng)值可由公式

(10)

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

(11)

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

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

2 三維反演算法

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

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

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

(12)

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

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

mk+1=mk+αkpk,

(13)

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

2.2 L-BFGS方法

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

sk=mk+1-mk,

(14)

yk=gk+1-gk,

(15)

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

pk=-Bkgk.

(16)

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

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

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

(17)

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

(18)

和曲率條件

(19)

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

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

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

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

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

3 數(shù)值算例

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

3.1 正演算例

3.1.1 DTM1模型

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

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

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

表1 DTM1模型自適應(yīng)加密過(guò)程中自由度數(shù)目及計(jì)算時(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é)者對(duì)此模型進(jìn)行了模擬(Miensopust et al., 2013).(0 km,0 km)位于三個(gè)異常體交界處的上方，其響應(yīng)最為奇異.圖6中展示了本文的計(jì)算結(jié)果與Mackie等(1993)的有限差分結(jié)果、Nam等(2007)等的有限元結(jié)果的對(duì)比.從圖中可以看出，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)勢(shì)之一是可以精確地模擬起伏地形.一般來(lái)說(shuō)，四面體單元的適應(yīng)性最強(qiáng)，可以模擬任意復(fù)雜的地形.在本文中，我們使用非結(jié)構(gòu)化六面體單元，通過(guò)對(duì)單元進(jìn)行形變也可模擬起伏地形.通過(guò)對(duì)一個(gè)二維山脊地形(Wannamaker et al., 1986)進(jìn)行模擬來(lái)驗(yàn)證程序?qū)Φ匦翁幚淼恼_性.該模型的背景電阻率為100 Ωm，模型的中間位置包含一平臺(tái)狀的地形，如圖7所示.計(jì)算了頻率為2 Hz時(shí)的響應(yīng).圖8中展示了初始網(wǎng)格和經(jīng)過(guò)5次自適應(yīng)加密得到的最終網(wǎng)格.初始網(wǎng)格非常稀疏，最終網(wǎng)格中，網(wǎng)格得到了充分加密.與DTM1模型類似，測(cè)點(diǎn)附近網(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.

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

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

3.2 反演算例

3.2.1 簡(jiǎn)單三維模型

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

圖10 簡(jiǎn)單三維模型示意圖Fig.10 The diagram of simple 3D model

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

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

3.2.2 三維地形模型

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

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

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

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

圖13 三維地形模型交錯(cuò)網(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 三維地形模型交錯(cuò)網(wǎng)格反演結(jié)果Fig.14 Inversion results of 3D topography model using staggered grids

圖15 三維地形模型交錯(cuò)網(wǎng)格反演過(guò)程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)格較密時(shí)反演非唯一性過(guò)強(qiáng)，反演結(jié)果并不理想.我們進(jìn)一步使用本文的方法對(duì)地形數(shù)據(jù)進(jìn)行反演，目標(biāo)區(qū)域的網(wǎng)格單元尺寸為2500 m，并對(duì)地表附近的網(wǎng)格進(jìn)行了加密和形變以適應(yīng)地形.反演初始模型為100 Ωm的均勻半空間，同時(shí)我們也進(jìn)行了不包含地形的反演.反演結(jié)果如圖16所示.

圖16 三維地形模型反演結(jié)果(a) 真實(shí)模型； (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.

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

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

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

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

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

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

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

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

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

(A1)

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

(A2)

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

(A3)

Py=c,

(A4)

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

Bz=d.

(A5)

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

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

記有限元解為E，單元上的后驗(yàn)誤差可表示為

(B1)

(B2)

(B3)

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

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

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

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

(B4)

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

附圖 B1八叉樹(shù)局部網(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個(gè)網(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ù)不同，則會(huì)產(chǎn)生懸掛點(diǎn).如附圖B2a所示，細(xì)網(wǎng)格中紅色棱邊與相鄰粗網(wǎng)格中較長(zhǎng)的棱邊部分重合，藍(lán)色棱邊與相鄰粗網(wǎng)格的面相交，破壞了有限元解的切向連續(xù)性，須對(duì)其施加約束.附圖B2b是附圖B2a中網(wǎng)格交界面的平面圖，與自由度x1、x2關(guān)聯(lián)的棱邊的長(zhǎng)度是與x0關(guān)聯(lián)的棱邊的長(zhǎng)度的一半，它們之間應(yīng)滿足的條件為

(B5)

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

(B6)