潘林冬， 李予國,2*， 趙 慧

(1.中國海洋大學(xué) 海洋地球科學(xué)學(xué)院，青島 266100;2.海底科學(xué)與探測技術(shù)教育部重點實驗室，青島 266100；3.國防科技大學(xué),長沙 410073)

0 引言

大地電磁測深法(MT)是地球物理學(xué)中的一種重要探測方法。它根據(jù)電磁感應(yīng)原理研究天然場源在大地中激勵的電磁場分布，并由觀測到的電磁場研究地球電性結(jié)構(gòu)。大地電磁測深法在地殼和上地幔地質(zhì)構(gòu)造研究等領(lǐng)域應(yīng)用廣泛[1-3],在石油天然氣勘探、地?zé)岷偷叵滤Y源調(diào)查等領(lǐng)域也發(fā)揮著重要作用[4-5]。1990年以來，大地電磁測深法已成功地應(yīng)用到大洋中脊和大陸架地質(zhì)構(gòu)造研究中[6-8]。

巖礦石在形成過程中受濕度、含水率等多種因素的影響，這些因素將導(dǎo)致巖礦石的電導(dǎo)率連續(xù)變化[9]。電導(dǎo)率連續(xù)變化介質(zhì)大地電磁場正演方法研究始于20世紀80年代。Kao 等[10-12]提出了層狀介質(zhì)電阻率隨深度線性變化或指數(shù)變化時大地電磁場計算方法;徐世浙等[13]利用有限元法實現(xiàn)了一維分層電阻率連續(xù)變化介質(zhì)大地電磁場正演。近年來，電導(dǎo)率分塊連續(xù)變化介質(zhì)二維和三維電磁場模擬研究也取得了一些進展。李予國等[14]用有限元法模擬了電導(dǎo)率分塊連續(xù)變化二維大地電磁場響應(yīng)；戴前偉等[15]討論了電導(dǎo)率分塊線性變化二維高頻大地電磁場有限元正演模擬；阮百堯等[16]提出了電導(dǎo)率連續(xù)變化三維電阻率測深有限元模擬方法；李勇等[17]用有限元法模擬電導(dǎo)率分塊連續(xù)變化三維可控源電磁場響應(yīng)。

眾所周知，結(jié)構(gòu)網(wǎng)格難以精確模擬海底地形起伏和傾斜界面等復(fù)雜地質(zhì)構(gòu)造和地下復(fù)雜形狀異常體。近年來，基于非結(jié)構(gòu)網(wǎng)格剖分技術(shù)的自適應(yīng)有限元方法在大地電磁場和海洋可控源電磁場數(shù)值模擬中取得了很好的應(yīng)用效果[18-22]。趙慧等[23]研究了非結(jié)構(gòu)三角單元網(wǎng)格自適應(yīng)有限元海洋大地電磁場二維正演算法，它能夠真實地模擬海底地形起伏、傾斜界面等復(fù)雜地質(zhì)構(gòu)造和地下復(fù)雜形狀異常體，并能夠克服海水與海底圍巖巨大電阻率差異引起的數(shù)值模擬困難，提供可靠穩(wěn)定的數(shù)值解。在自適應(yīng)網(wǎng)格細化過程中，采用后驗誤差估計指導(dǎo)網(wǎng)格細化，只對具有較大局部誤差的網(wǎng)格進行細化，從而減少了計算內(nèi)存需求和節(jié)約了計算時間。我們將以上算法推廣到二維電導(dǎo)率連續(xù)變化介質(zhì)中，通過在有限元網(wǎng)格節(jié)點上進行電導(dǎo)率賦值，并且利用單元塊屬性區(qū)分開不同區(qū)域塊體同一邊界上的電導(dǎo)率賦值，使得它能夠模擬電導(dǎo)率連續(xù)變化介質(zhì)二維大地電磁場響應(yīng)。

1 基本理論

1.1 大地電磁場的邊值問題及有限單元法

考慮二維大地電磁模型，假設(shè)在x軸方向電導(dǎo)率保持不變，z軸垂直地表指向地下。平行走向電磁場分量滿足如下偏微分方程[24]

▽·(τ▽u)+λu=0

(1)

其中：ω為角頻率；σ為電導(dǎo)率(S/m)；介質(zhì)的磁導(dǎo)率μ0=4π×10-7H/m。

假定有限元模擬區(qū)域離二維不均勻體足夠遠，以致不均勻體產(chǎn)生的異常電磁場在其邊界上衰減為零。在模擬區(qū)域左、右邊界上可以采用Dirichlet邊界條件，其電磁場表達式在文獻[23]中已經(jīng)給出。

二維大地電磁場邊值問題與下列泛函極值問題等價[23-24]。

(2)

其中：Ω為模擬區(qū)域；Γ為模擬區(qū)域Ω的邊界。

筆者采用三角網(wǎng)格有限單元法求解上述變分問題，將模擬區(qū)域Ω剖分成ne個三角形單元。由此，方程(2)的積分就轉(zhuǎn)換為各個單元積分之和：

(3)

其中：Γe表示底邊界CD上的線段；e表示三角單元。

在以往的數(shù)值模擬中，假定單元內(nèi)的電性參數(shù)(τ和λ)為常數(shù)，即在單元內(nèi)電性是均勻的，電性參數(shù)保持不變，但是連續(xù)變化的電性參數(shù)可能更符合實際情況。在我們的算法中，假定單元內(nèi)電性參數(shù)(τ和λ)與電磁場u均是線性變化的，且可以近似為：

(4)

其中：ui、τi和λi分別是三角單元中第i(i=1,2,3)個節(jié)點上u、τ和λ的函數(shù)值；Ni是三角單元的形函數(shù)，其表達式如式(5)所示。

(5)

a1=z2-z3,b1=y3-y2,c1=y2z3-y3z2

a2=z3-z1,b2=y1-y3,c2=y3z1-y1z3

(6)

a3=z1-z2,b3=y2-y1,c3=y1z2-y2z1

方程(3)中第一個單元面積分可以寫成式(7)。

(7)

其中：ue=(u1,u2,u3)T；K1e為單元矩陣

(8)

如果令式(8)中,τ1=τ2=τ3=τ,即得電性參數(shù)為常數(shù)情形的單元矩陣。

方程(3)中第二個單元面積分可以表示為式(9)。

(9)

(10)

如果令式(10)中,λ1=λ2=λ3=λ,即得電性參數(shù)為常數(shù)情形的單元矩陣K2e。

圖1 線性插值三角單元Fig.1 Atriangular element

(11)

其中：K3e是一個3×3的對稱矩陣，其非零元素為：

(12)

其中：

δ11=12τ1k1+3τ1k2+3τ2k1+2τ2k2

δ21=3τ1k1+2τ1k2+2τ2k1+3τ2k2

(13)

δ22=2τ1k1+3τ1k2+3τ2k1+12τ2k2

令式(13)中τ1=τ2=τ3=τ，k1=k2=k3=k，即得到CD邊上電性參數(shù)為常數(shù)時單元矩陣K3e的元素。

計算得到每個單元的系數(shù)矩陣后，將其擴展和合并成總體系數(shù)矩陣為式(14) 。

(14)

于是得到泛函的離散形式為式(15)。

(15)

對方程(15)求變分得

KU=0

(16)

方程(16)的系數(shù)矩陣是對稱的，但是奇異的。換句話說，需要代入邊界條件，線性方程組才是唯一可解的。代入邊界條件后，方程(16)變?yōu)槭?17)。

(17)

求解方程(17)即可得到電磁場平行走向分量。

由上述有限單元法得到的平行走向電場或磁場值的精度與網(wǎng)格緊密相關(guān)，合理可靠的離散網(wǎng)格設(shè)計是獲得高精度數(shù)值模擬結(jié)果的關(guān)鍵。對于簡單的地電模型，基于經(jīng)驗可以得到優(yōu)化離散網(wǎng)格，而對于復(fù)雜模型，僅憑研究者的經(jīng)驗難以取得優(yōu)化網(wǎng)格。新近發(fā)展起來的自適應(yīng)有限元法能夠自動細化網(wǎng)格,并提供可靠的數(shù)值解,該方法所涉及的主要技術(shù)是后驗誤差估計和自適應(yīng)網(wǎng)格細化。

1.2 非結(jié)構(gòu)網(wǎng)格和自適應(yīng)有限元模擬技術(shù)

非結(jié)構(gòu)網(wǎng)格具有更好地模擬海底地形起伏和地下復(fù)雜異常體的能力；自適應(yīng)有限元模擬技術(shù)是指在得到的初始網(wǎng)格基礎(chǔ)上，利用基于梯度恢復(fù)的后驗誤差估計對需要提高精度的網(wǎng)格進行加密，使得到的解滿足精度要求。自適應(yīng)有限元模擬涉及的主要技術(shù)問題，是后驗誤差估計和自適應(yīng)網(wǎng)格細化。

本文算法采用基于梯度恢復(fù)的后驗誤差估計方法對有限元網(wǎng)格進行評判，指導(dǎo)離散網(wǎng)格的細化，提高解數(shù)值的精度[25]。在L2范數(shù)條件下，恢復(fù)梯度與有限元解的梯度之差為給定單元上的局部誤差指示子:

ηe=‖(R-I)▽uh‖L2(e)

(18)

式中：R為一個梯度恢復(fù)算子；I為單位算子；▽uh為有限元解的梯度；恢復(fù)梯度為R▽uh。Bank等[26]提出了一種新的超收斂梯度恢復(fù)算子R=SmQh，其中Qh為L2投影算子，m為光滑算子平滑迭代次數(shù)(通常取2)。采用平滑算子Sm對Qh▽uh進行平滑，uh為分段線性的，因此▽uh是分段常數(shù)。隨著單元減小，后驗誤差估計逐漸逼近真實的梯度誤差[26]。后驗誤差估計表達式為式(19)。

ηe=‖(SmQh-I)▽uh‖L2(e)

(19)

基于局部誤差指示子的網(wǎng)格細化，對于獲得求解區(qū)域的精確解是合適的，可得到全局精度的解。但是在很多情況下，僅僅需要求得特定位置(如測點)的精確解[18]。如果想要減小測點處數(shù)值解的誤差，需要考慮整個區(qū)域?qū)植空`差的影響，可采用對偶加權(quán)誤差估計方法[27]解決這個問題。

采用內(nèi)積形式，將方程(1)轉(zhuǎn)化成等價的變分問題，可得：

B(u,v)=0

(20)

其中：

(21)

式中：TE模式時u=Ex,k=1,q=iωu0σ；TM模式時，u=Hx,k=1/σ,q=iωu0。

定義泛函G為電磁場偏微分方程的精確解與有限元數(shù)值解的誤差u-uh的函數(shù)，則有如下對偶方程為式(22)。

B*(w,v)=G(v)

(22)

這里B*(w,v)是對偶或伴隨算子，假定B*(w,v)=B(v,w)，則有

G(u-uh) =B*(w,u-uh)

=B(u-uh,w)

=B(u-uh,w-wh)

(23)

其中：w和wh分別為對偶問題的解析解與有限元解。由式(23)可得:

(R-I)▽uhdΩ

(24)

根據(jù)式(24)，定義基于對偶加權(quán)的誤差指示子：

(25)

其中，

(26)

(27)

自適應(yīng)有限元算法的基本思路如下：①將模擬區(qū)域進行初始剖分，產(chǎn)生一個粗糙的初始網(wǎng)格，并在該初始網(wǎng)格上進行有限元正演計算得到數(shù)值解；②計算每個單元的局部誤差，選取一定比例具有較大局部誤差的單元進行網(wǎng)格細化，產(chǎn)生一個新網(wǎng)格,并在新網(wǎng)格上進行有限元正演計算；③重復(fù)以上過程，直到有限元解滿足假定的收斂精度或達到最大網(wǎng)格細化次數(shù)為止。我們設(shè)定的細化網(wǎng)格百分比為5%，收斂精度為1%。具體操作流程如圖2所示。

2 算例

2.1 一維模型

為了驗證本文算法的正確性及有效性,模擬一個三層地電模型(圖3)的大地電磁測深響應(yīng)。假定第一層和第三層的電導(dǎo)率分別為0.1S/m和0.01S/m，第一層的厚度為1 000m。第二層是厚度為500m的過渡層，其電導(dǎo)率隨深度線性增加。

首先，按文獻[24]中一維電導(dǎo)率連續(xù)變化有限元理論編寫了程序；然后，對比程序計算結(jié)果與徐世浙《有限單元法》中的結(jié)果一致；最后，我們模擬電導(dǎo)率連續(xù)變化一維模型大地電磁場響應(yīng),并與一維電導(dǎo)率連續(xù)變化有限元算法結(jié)果進行比較。

圖2 自適應(yīng)有限元正演算法流程圖Fig.2 Adaptive finite element forwarding algorithm flow chart

圖3 一維電導(dǎo)率連續(xù)變化地電模型Fig.3 A 1D model with linear variation of conductivity

圖4為筆者所述三角單元網(wǎng)格剖分自適應(yīng)有限元算法得到的大地電磁測深視電阻率曲線和相位曲線，與用文獻[13]一維有限元法及二維矩形網(wǎng)格剖分有限元方法所得計算結(jié)果的對比。由圖4可見，三種方法的計算結(jié)果非常一致。

假定視電阻率相對誤差和相位絕對誤差為：

(28)

(29)

圖4 一維層狀模型大地電磁響應(yīng)曲線Fig.4 Magretoele ctromagnetic response curve of one-dimensional layered model(a)視電阻率測深曲線圖；(b)相位曲線

圖5 等效二維地塹模型Fig.5 Equivalent two-dimensional graben model(a)二維含電導(dǎo)率連續(xù)變化層地塹模型圖；(b)二維電導(dǎo)率分塊均勻體地塹模型

2.2 二維地塹模型

圖5(a)為一個四層地塹模型。第一層是電導(dǎo)率為3.333 3S/m，厚度為1 000m的海水層。第二層為500m厚的過渡層，其電導(dǎo)率隨深度線性增加：

σ(z)=3.3333-3/500×(z-1000)

(30)

圖6 視電阻率隨周期變化曲線圖Fig.6 MT apparent resistivities for TE mode and TM mode(a)TE模式；(b)TM模式

圖7 相位隨周期變化曲線圖Fig.7 MT impedance phase for TE mode and TM mode(a)TE模式；(b)TM模式

2.3 復(fù)雜地電模型

設(shè)計一個如圖8所示的復(fù)雜二維海洋地電模型，該模型關(guān)于y=0 m對稱。海水層電導(dǎo)率3.333 3 S/m，厚度為1 000 m。在海底下方，有兩個厚度各為1 000 m的電導(dǎo)率隨深度線性變化水平地層。第一層的電導(dǎo)率隨深度線性減少，由頂部處的3.333 3 S/m減小到底部處的0.333 3 S/m。第二層的電導(dǎo)率隨深度線性增加，由頂部處的0.333 3 S/m增加到底部處的1 S/m。海底下方第三層是電導(dǎo)率為1 S/m的半空間，但其中含有一個形態(tài)較復(fù)雜的二維異常體。該異常體由上、下兩部分組成，下面部分是電阻率為10 S/m的高導(dǎo)體，而上面部分電導(dǎo)率線性連續(xù)變化。

圖9和圖10分別為圖8地電模型的初始網(wǎng)格和細化的最終網(wǎng)格，由圖9和圖10可以看出在接收站處網(wǎng)格得到了加密。

我們模擬了該模型的大地電磁測深響應(yīng)。計算頻率范圍為10-4Hz～101Hz，計算頻率點均勻分布，每個數(shù)量級取10個頻點，總共計算了51個頻率的大地電磁場響應(yīng)。接收站在-15 000 m～15 000 m范圍內(nèi)按點距為500 m布設(shè)在海底，總共61個測點。

圖8 二維復(fù)雜地質(zhì)體電導(dǎo)率連續(xù)變化地電模型Fig.8 A two-dimensional model with linear variation of conductivity

圖9 二維復(fù)雜地質(zhì)體初始網(wǎng)格Fig.9 Initial mesh of two - dimensional complex geological

圖10 二維復(fù)雜地質(zhì)體細化最終網(wǎng)格Fig.10 Two dimensional complex geological body refinement final mesh

圖11為TE模式大地電磁測深視電阻率和相位斷面圖。由圖11可以很清楚地推斷出深部高導(dǎo)體(低阻體)，也能看出海底下方兩個電導(dǎo)率連續(xù)變化的水平層。圖12為TM模式的視電阻率和相位斷面圖。對比圖11(a)和圖12(a)可知，TE模式深部視電阻率曲線是圈閉的，而TM模式視電阻率曲線在最下方并未圈閉，說明了TE模式對勾畫異常體的范圍效果更好些。

圖11 TE模式視電阻率斷面圖和相位斷面圖Fig.11 TE mode apparent resistivity pseudosection and phase pseudosection(a) 視電阻率斷面圖；(b)相位斷面圖

圖12 TM模式視電阻率斷面圖和相位斷面圖Fig.12 TM mode apparent resistivity pesudosection and phase pseudosection(a) 視電阻率斷面圖；(b)相位斷面圖

3 結(jié)論

研究了電導(dǎo)率連續(xù)變化二維介質(zhì)大地電磁場自適應(yīng)有限元正演算法。與電導(dǎo)率分塊均勻介質(zhì)數(shù)值模擬方法不同，我們在單元節(jié)點上進行電導(dǎo)率賦值，在不同區(qū)域塊體的相同邊界的節(jié)點電導(dǎo)率利用區(qū)域塊編號屬性來區(qū)分，單元內(nèi)各點處的電導(dǎo)率通過線性插值獲到，從而使得我們能夠模擬電導(dǎo)率線性連續(xù)變化介質(zhì)電磁場響應(yīng)。同時，由于采用了非結(jié)構(gòu)三角單元網(wǎng)格剖分技術(shù)，我們的算法能夠更好地模擬復(fù)雜地電模型。

致謝

感謝外審專家提出的寶貴意見，感謝劉穎博士的幫助和討論。

[1] 馬曉冰，孔祥儒. 青藏高原巖石圈熱狀態(tài)及其東西部差異[J].地球物理學(xué)進展，2001，16(3)：12-20. MA X B，KONG X R. The thermal status of Qinghai-Tibet plateau and the differen-ces between the western and the eastern plateau[J]. Progress in Geophysics，2001，16(3)：12-20.(In Chinese)

[2] 魏文博，鄧明，譚捍東，等.我國海底大地電磁探測技術(shù)研究的進展[J].地震地質(zhì)，2001，23(2)：131-137. WEI W B，DENG M，TAN H D，et al. Develo-pment of marine magnetotelluricprospe-cting technique in China[J]. Seismology a-nd Geology, 2001，23(2)：131-137. (In Chinese)

[3] 劉國棟. 我國大地電磁測深的發(fā)展[J]. 地球物理學(xué)報，1994，37(增刊1)：301-310. LIU G D. Delopment of magnetotelluric m-ethod in china[J]. Actageophysicasinic-a，1994，37(supp.1)：301-310. (In Chinese)

[4] 何展翔，賈進斗，茍量.非地震技術(shù)在油氣勘探開發(fā)中的作用[J].石油勘探與開發(fā)，2001，28(4)：70-72. HE Z X, JIA J D, GOU L. The role of non-seismic techniques in exploration and development of oil and gas[J]. Petroleum E-xploration and Development，2001，28(4)：70-72. (In Chinese)

[5] 趙建良，陳天振，張晉，等. MT法在開封凹陷地?zé)豳Y源調(diào)查中的應(yīng)用[J].物探與化探，2010，34(2)：163-166. ZHAO J L，CHEN T Z，ZHANG J,et al. The applic-ation of the MT method to the investin-gation of geothermal resources in Kaif-eng depression[J].Geophysical and Geoche-mical Exploration，2010，34(2)：163-166. (In Chinese)

[6] EVANS R L，LAW L K，LOUIS B S，et al.Theshallow porosity structure of the Eel- shelf, northern California: results of a towed electromagnetic survey[J]. Marine Geology，1999，154(1)：211-226.

[7] HOVERSTEN G H，CONSTABLE S，MORRISON H F. Marine magnetotellurics for base salt mapping: Gulf of Mexico field-test at the Gemini structure[J]. Geophysics，2000，65(5)：1476-1488.

[8] WORZEWSKI T，JEGEN K M，KOPP H，et al. Magnetotelluric image of the fluid cycle in the costa rican subduction zone[J]. Nature Geoscience，2011，4(2)：108-111.

[9] 李金銘.地電場與電法勘探[M].北京：地質(zhì)出版社，2005. LI J M. Geoelectric field and electrical exploration[M]. BeiJing:Geological Publishing Ho-use，2005. (In Chinese)

[10]KAO D, RANKIN D.Magnetotelluricrespo-nse on inhomogeneous layered earth[J]. G-eophysics，1980，45: 1973-1802.

[11]KAO D. Magnetotelluric response on vertically inhomogeneous earth[J]. J. Geophys. Res，1981，86： 3027-3038.

[12]KAO D. Magnetotelluric response on vertically inhomogeneous earth having conductivity varying exponentially with depth[J]. Geophysics，1982，47：89-99.

[13]徐世浙，劉斌.電導(dǎo)率分層連續(xù)變化的水平層的大地電磁正演[J].地球物理學(xué)報，1995，38(2)：262-268. XU S Z，LIU B. A numerical method for calculating MT field on a layered model with continuous change of conductivity in each layer[J]. Actageophysicasinica，1995，38(2)：262-268. (In Chinese)

[14]李予國，徐世浙，劉斌,等.電導(dǎo)率分塊連續(xù)變化的二維MT有限元模擬(II)[J]. 高校地質(zhì)學(xué)報，1996，2(4)： 448-452. LI Y G，XU S Z，LIU B,et al. The finite element method for modeling 2-D MT fi-eld on a geoelectrical model with con-ductivity within each block[J]. Geologic-al Journal of China Universities，1996，2(4)：448-452.

[15]戴前偉，張富強，楊坤平，等.電導(dǎo)率分塊線性變化的二維高頻大地電磁法有限元正演模擬[J].物探化探計算技術(shù)，2012，34(5)：552-559. DAI Q W，ZHANG F Q，YANG K P. et al. Th-e finite element method for modeling 2-D high-frequency electromagnetic wit-h continuous of variation of conducti-vity within each block[J]. Computing tec-hniques for geophysical and geochemical exploration，2012，34(5)：552-559. (In Chinese)

[16]阮百堯，熊彬.電導(dǎo)率連續(xù)變化的三維電阻率測深有限元模擬[J].地球物理學(xué)報，2002，45(1)：131-138. RUAN B Y，XIONG B. A finite element mo-deling of three-dimensional resistive-ty sounding with continuous conductiv-ity[J]. Chinese journal of geophysics，2002，45(1)：131-138. (In Chinese)

[17]李勇，吳小平，林品榮. 基于二次場電導(dǎo)率分塊連續(xù)變化的三維可控源電磁有限元數(shù)值模擬[J].地球物理學(xué)報，2015，58(3)：1072-1087. LI Y，WU X P，LIN P R. Three-dimensiona-l controlled source electromagnetic f-inite simulation using the secondary field for continuous variation of elec-trical conductivity within each block[J]. Chinese journal of geophysics，2015，58(3)：1072-1087. (In Chinese)

[18]KEY K，WEISS C. Adaptive finite-elemen-t using unstructured grids：The 2D magn-etotelluric example[J]. Geophysics，2006，71(6)：G291-G299.

[19]LI Y G，KEY K. 2D marine controlled so-urce electromagnetic modeling: Part 1-An adaptive finite element algorithm[J]. Geophysics，2007，72(2)：WA51-WA62.

[20]LI Y G，CONSTABLE S. 2D marine control-led source electromagnetic modeling: Part 2-The effect of bathymetry[J]. Geoph-ysics, 2007，72(2)：WA63-WA71.

[21]LI Y G，PEK J. Adaptive finite element modeling of two=dimensional magnetite-lluric fields in general anisotropic media[J]. Geophys J Int，2008，175(3)：942-954.

[22]LI Y G，DAI S K. Finite element modeli-ng of marine controlled-source electro-magnetic responses in two-dimensional dipping anisotropic conductivity struct-ures[J]. Geophys J Int，2011，185(2)： 622-636.

[23]趙慧，劉穎，李予國.自適應(yīng)有限元海洋大地電磁場二維正演模擬[J].石油地球物理勘探，2014，49(3)： 578-585. ZHAO H，LIU Y，LI Y G. Adaptive finite element forward modeling for two-dimen-sional marine magnetotelluric field[J]. OGP，2014，49(3)：578-585. (In Chinese)

[24]徐世浙.地球物理中的有限單元法[M].北京：科學(xué)出版社，1994. XU S Z.Finite element method in geoph-ysics[M]. BeiJing:Science Publishing Company，1994. (In Chinese)

[25]劉穎. 海洋可控源電磁法二維有限元正演和反演 [D].青島:中國海洋大學(xué)，2014. LIU Y. 2D finite element modeling and inversion for marine controlled-source electromagnetic fields[D]. Qindao: OUC，2014. (In Chinese)

[26]BANK R E ,XU J. Asymptotically exact a posteriori error estimators, Part I: Grids with superconvergence[J]. SIAM journal on numerical analysis, 2013,41(6):2294-2312.

[27]OVALL J S. Asymptotically exact functional error estimators based on superconvergent gradient recovery[J]. Numersische mathematic, 2006，102(3): 543-558.