劉長(zhǎng)生 ，湯井田，任政勇 ，馮德山

(1. 中南大學(xué) 信息物理工程學(xué)院，湖南 長(zhǎng)沙，410083；2. 長(zhǎng)沙航空職業(yè)技術(shù)學(xué)院 計(jì)算機(jī)系，湖南 長(zhǎng)沙，410014；3. 瑞士聯(lián)邦理工學(xué)院(EIH) 地球物理系，瑞士 蘇黎士，CH8092)

自 Coggon[1]提出地球物理電磁中的有限單元算法(Finite element method, FEM)以來(lái)，F(xiàn)EM開(kāi)始在電磁勘探領(lǐng)域得到廣泛應(yīng)用[2-5]。Badea等[6]利用節(jié)點(diǎn)型的線性有限單元模擬了可控源音頻大地電磁法；Mitsuhata等[7]基于電場(chǎng)標(biāo)量勢(shì)和磁場(chǎng)矢量勢(shì)的 T-Ω公式，利用線性的向量有限單元計(jì)算了三維大地電磁模型；Nelson等[8-10]利用非結(jié)構(gòu)化網(wǎng)格來(lái)解決網(wǎng)格剖分的幾何離散化誤差問(wèn)題，從而計(jì)算了三維電磁矢量有限元模擬；阮百堯等[11]利用節(jié)點(diǎn)型有限單元實(shí)現(xiàn)了三維地電斷面的正演等；王燁等[12]采用基于棱邊的矢量有限元方法計(jì)算了三維大地電磁模型。目前，地球物理電磁FEM計(jì)算存在1個(gè)問(wèn)題：數(shù)值計(jì)算結(jié)果的精度完全取決于初始模型的離散化網(wǎng)格。對(duì)于簡(jiǎn)單的電磁模型，基于經(jīng)驗(yàn)可以得到較優(yōu)化的離散網(wǎng)格；而對(duì)于復(fù)雜 3D地電模型，電磁場(chǎng)的形態(tài)和變化趨勢(shì)復(fù)雜，僅憑人為經(jīng)驗(yàn)難以得到優(yōu)化網(wǎng)格。為了提高復(fù)雜區(qū)域的電磁場(chǎng)精度，網(wǎng)格可以進(jìn)行加密(h型自適應(yīng)[13-14])，另外也可以提供局部擬合電磁場(chǎng)的形狀多項(xiàng)式函數(shù)的階次p(p型自適應(yīng)[15-16])。相對(duì)于p型自適應(yīng)策略，h型自適應(yīng)策略更容易嵌入現(xiàn)有的代碼之中，從而縮短動(dòng)運(yùn)算時(shí)間。Key等[17]利用h型自適應(yīng)網(wǎng)格的思想對(duì)二維大地電磁模型的模擬進(jìn)行了研究；Li等[18]將這一概念引入到二維可控源電磁法之中。對(duì)于三維模擬，本文作者利用非結(jié)構(gòu)化網(wǎng)格，提出基于后驗(yàn)誤差的h型自適應(yīng)加密的3D大地電磁模擬算法。

1 三維MT矢量有限元數(shù)學(xué)模型

三維 MT 模型的一般地電結(jié)構(gòu)如圖 1所示，其中：0σ，1σ和2σ分別為空氣層、地下地層和地下埋藏的異常體的電導(dǎo)率。

圖1 三維MT模型的一般地電結(jié)構(gòu)Fig.1 MT model of three-dimensional structure of general earth electricity

在空氣層和地下地層區(qū)域，電場(chǎng)滿足下面微分控制方程[19]：

式中：E為電場(chǎng)強(qiáng)度；i為虛數(shù)；ω為角頻率；μ為自由空間磁導(dǎo)率；σ為MT模型的電導(dǎo)率。

采用理論可靠和執(zhí)行簡(jiǎn)單的Direlet邊界條件[20]，并采用水平地層地電結(jié)構(gòu)在邊界上的電磁場(chǎng)近似作為邊界條件。

式中：n為外邊界的外單位法向量；E0為外邊界上給定的已知電場(chǎng)強(qiáng)度。

方程(1)和(2)定義了MT問(wèn)題的邊值問(wèn)題。借助于矢 量 恒 等 式Gauss分部積分公式，其對(duì)應(yīng)的變分公式為：

式中：H(curl)為 Hilbert空間向量， H (curl)=；Ω 為求解區(qū)域；L2(Ω)為平方可視的Hilbert空間；b為雙線性表達(dá)式；f為單線性表達(dá)式。

此處采用矢量有限元法求解方程(3)所表示的電場(chǎng)分布。為此，需要把整個(gè)求解區(qū)域離散成一系列四面體單元。它與節(jié)點(diǎn)型有限元不同之處在于：四面體單元中電場(chǎng)強(qiáng)度E由每條邊的切向方向定義。考慮任意的四面體單元 e，并把該單元中電場(chǎng)的近似表達(dá)式表示為：

為了避免偽數(shù)值解，節(jié)點(diǎn)型有限元法必須加入額外的罰項(xiàng)，進(jìn)而加大了理論難度。而矢量有限元法由于其矢量形狀函數(shù)的散度為0，電流密度的無(wú)散度要求自然得到滿足；并考慮到矢量形狀函數(shù)僅僅沿四面體邊的切向連續(xù)，因此，電場(chǎng)的切向連續(xù)性條件得到滿足，法向不連續(xù)性條件被保留，故基于矢量形狀函數(shù)的四面體單元能有效地避免出現(xiàn)偽解。

MT 模型求解的大型復(fù)線性方程為：

其中：U為定義在各個(gè)四面體單元的邊上電場(chǎng)切向未知數(shù)向量。矩陣A和B由各個(gè)單元的貢獻(xiàn)所得[21]：

式中：Bi僅僅在位于外Dirichlet邊界面上的邊不為0。

2 基于殘差的后驗(yàn)誤差估計(jì)

采用Nédélec四面體線性單元，其求解過(guò)程可以表述為：定義Tk為計(jì)算區(qū)域Ω的一系列四面體網(wǎng)格剖分單元，F(xiàn)k為網(wǎng)格中的一系列三角面；定義Tk上的有限元空間Uk為[21]：

在Tk上，有限元計(jì)算可表示為：求解Ek∈Uk，使其滿足

式中：)(2ΩL∈f；)(2ΩL?∈g。

定義數(shù)值解的誤差為：

經(jīng)推導(dǎo)，適合于MT模型的后驗(yàn)誤差估計(jì)為：

式中：C為一依賴于問(wèn)題的常數(shù)；ηT和ηF分別為四面體單元和三角面的局部誤差。

3 自適應(yīng)策略

采用h型自適應(yīng)有限元對(duì)MT模型進(jìn)行數(shù)值模擬，網(wǎng)格剖分采用約束的 Delaunay四面體(DT)非結(jié)構(gòu)化方法及相應(yīng)的加密算法。對(duì)于任意模型，只要給出初始網(wǎng)格和全局誤差控制值，利用h型自適應(yīng)迭代策略，計(jì)算機(jī)即可全自動(dòng)求解，無(wú)需人工干預(yù)。

在四面體網(wǎng)格Tk上，利用有限元方法求出Ek，根據(jù)文獻(xiàn)[22]中的誤差收斂理論，得出當(dāng)前網(wǎng)格上的誤差估計(jì)，并得到單元上的局部誤差估計(jì)指示值：

相應(yīng)地，全局誤差指示值和最大誤差估計(jì)指示值分別為：

當(dāng)計(jì)算單元誤差超過(guò)給定的誤差限后，采用相應(yīng)的加密策略，具體算法流程如圖2所示。

圖2 具體算法流程圖Fig.2 Flow chart of calculation method

算法中影響自適應(yīng)加密過(guò)程的終止因子為1個(gè)給定小正數(shù)ε和加密密度因子β。通常來(lái)說(shuō)，ε小于1，而β的最佳值為0.5～0.7。對(duì)于不同的問(wèn)題，正數(shù)ε也相差較大，ε[1∈×10-10, 0.1]。另外，在實(shí)際計(jì)算過(guò)程中，還引入最大迭代次數(shù)n和最大內(nèi)存消耗M來(lái)控制自適應(yīng)終止，一般取n=10，M=500 MB。

4 COMMEMI 3D-1 MT模型數(shù)值模擬算例

為了加強(qiáng)對(duì)比分析，選取國(guó)際通用的標(biāo)準(zhǔn)測(cè)試模型COMMEMI 3D-1 model進(jìn)行模擬[23]。其中：六面體異常體的電阻率為 0.5 Ω·m，圍巖的背景電阻率為100 Ω·m，其對(duì)比度達(dá)1∶200。求解的空間維數(shù)如下：x方向長(zhǎng)度為100 km，y方向長(zhǎng)度為100 km，z方向上空氣層厚度為50 km，大地層厚度為50 km。MT測(cè)線分布在x和y坐標(biāo)軸上，長(zhǎng)度范圍均為0～3 km。測(cè)點(diǎn)個(gè)數(shù)為121，工作頻率為f=0.1 Hz。

自適應(yīng)控制參數(shù)為：ε=0.01，β=0.7，n=25，M=500 MB。實(shí)際迭代步數(shù)為20步，內(nèi)存最大為128.5 MB。初始網(wǎng)格為 10 553個(gè)節(jié)點(diǎn)，65 896個(gè)單元和153 060條邊；自適應(yīng)加密的第12步網(wǎng)格規(guī)模為13 238個(gè)節(jié)點(diǎn)，79 658個(gè)單元和186 854條邊；最后一步的網(wǎng)格為19 538個(gè)節(jié)點(diǎn)，113 281個(gè)單元和267 798條邊。為了計(jì)算方便，測(cè)試平臺(tái)選用DELL D620筆記本電腦進(jìn)行數(shù)值模擬，計(jì)算頻點(diǎn)分別是32，16，8，4，2，1，0.5，0.25和0.1 Hz，總耗時(shí)為317.461 4 s。其中：初次網(wǎng)格是采用tetgen軟件進(jìn)行剖分，自動(dòng)加密采取二分加密，時(shí)間均少于 1 s。自動(dòng)頻點(diǎn)計(jì)算最大耗時(shí)56.811 0 s，最小耗時(shí)23.319 5 s，平均耗時(shí)35.273 5 s。

圖3 頻率為0.1 Hz、含不均勻立方體MT模型的地表測(cè)點(diǎn)上的相位和視電阻率曲線Fig.3 Phase and resistivity curves containing uneven surface cube measuring point when frequency is 0.1 Hz

通過(guò)模擬得到 3D-1模型的自適應(yīng)加密網(wǎng)格的后驗(yàn)誤差分布。以0.1 Hz頻點(diǎn)模擬情況進(jìn)行分析：在初始網(wǎng)格上的單元誤差較大，平均誤差為74.48%，最大誤差為 165.00%。因此，較大的誤差分布驅(qū)使設(shè)計(jì)的自適應(yīng)加密策略運(yùn)作，并加密誤差大的單元。這一結(jié)果在第12和20次網(wǎng)格上可以明顯看出。在第12次網(wǎng)格上，單元平均誤差減至3.33%，最大誤差減至6.12%；在第 20次網(wǎng)格上，單元平均誤差減至 0.93%,最大誤差減至3.98%，其平均誤差小于給點(diǎn)的1%，因此，加密過(guò)程終止。數(shù)值模擬結(jié)果與Sasaki針對(duì)本模型所得計(jì)算結(jié)果較吻合(見(jiàn)圖3中的Sasaki(SFD))[7]。

注意到測(cè)線附近的單元并沒(méi)有隨著自適應(yīng)加密過(guò)程的運(yùn)行而對(duì)網(wǎng)格顯著加密。這一現(xiàn)象導(dǎo)致測(cè)線附近數(shù)值精度提高(見(jiàn)圖3)：視電阻率的平均誤差從初始網(wǎng)格上的3.49%減低為第12次網(wǎng)格上的0.87%，但是，隨著網(wǎng)格的再加密，視電阻率平均誤差始終保持在0.50%～1.00%；相位曲線的平均誤差從初始網(wǎng)格上的1.23%減低為第20次時(shí)的0.21%。但是，隨著網(wǎng)格的再加密，相位曲線平均誤差始終保持在0.20%～0.50%。

5 結(jié)論

(1) 根據(jù)電場(chǎng)微分控制方程和邊界條件，采用基于矢量形狀函數(shù)的四面體單元，其矢量形狀函數(shù)的散度為0，電流密度的無(wú)散度要求就被自然滿足，因此，能有效地避免偽解的出現(xiàn)，從而建立了三維MT矢量有限元數(shù)學(xué)模型。

(2) 推導(dǎo)出基于殘差的三維大地電磁矢量有限元后驗(yàn)誤差估計(jì)公式，為三維大地電磁自適應(yīng)矢量有限元數(shù)值模擬的實(shí)現(xiàn)奠定了基礎(chǔ)。

(3) 在完全非結(jié)構(gòu)化四面體單元剖分及優(yōu)化基礎(chǔ)上，結(jié)合三維大地電磁矢量有限元后驗(yàn)誤差估計(jì)公式，提出了基于非結(jié)構(gòu)化網(wǎng)格的三維大地電磁h型自適應(yīng)矢量有限元計(jì)算策略，保證了對(duì)復(fù)雜大地電磁模型數(shù)值計(jì)算的精度和可靠性。

(4) 電磁模型的復(fù)雜性總體上不影響方法的收斂性，且可以達(dá)到預(yù)期的計(jì)算精度?？梢?jiàn)，h型自適應(yīng)矢量有限元可以保證復(fù)雜模型的計(jì)算精度和速度，具有廣闊的應(yīng)用前景。

[1] Coggon J H. Electromagnetic and electrical modeling by the finite element method[J]. Geophysics, 1971, 36(1): 132-145.

[2] Pridmore D F, Hohmanns G W, Ward S H, et al. An investigation of finite-element modeling of electric and electromagnetic data in three dimensions[J]. Geophysics, 1981, 46(7): 1009-1024.

[3] 李大潛. 有限元素法在電法測(cè)井中的應(yīng)用[M]. 北京: 石油工業(yè)出版社, 1980: 15-67.LI Da-qian. The application of finite element method in electric well logging[M]. Beijing: Petroleum Industry Press, 1980:15-67.

[4] 羅延鐘, 張桂青. 電子計(jì)算機(jī)在電法勘探中的應(yīng)用[M]. 武漢:武漢地質(zhì)學(xué)院出版社, 1987: 51-99.LUO Yan-zhong, ZHANG Gui-qing. Application of electronic computer in electrical prospecting[M]. Wuhan: Press of Wuhan College of Geology, 1987: 51-99.

[5] 徐世浙. 地球物理中的有限單元法[M]. 北京: 科學(xué)出版社,1994: 23-28.XU Shi-zhe. Finite element method in geophysics[M]. Beijing:Science Press, 1994: 23-28.

[6] Badea E A, Everett M E, Newman G A, et al. Finite-element analysis of controlled-source electromagnetic induction using Coulomb-gauged potentials[J]. Geophysics, 2001, 66(3):786-799.

[7] Mitsuhata Y, Uchida T. 3D magnetotelluric modeling using the T-Ω finite-element method[J]. Geophysics, 2004, 69(1):108-119.

[8] Nelson E M. Advances in 3D Electromagnetic finite element modeling[C]//Proceedings of the IEEE 1997 Particle Accelerator Conference. Vancouver, Canada, 1998: 1837-1840.

[9] Sugeng F. Modeling the 3D TDEM response using the 3D full-domain finite-element method based on the hexahedral edge-element technique[J]. Exploration Geophysics, 1998, 29(4):615-619.

[10] Yoshimura R, Oshiman N. Edge-based finite element approach to the simulation of geoelectromagnetic induction in a 3-D sphere[J]. Geophysical Research Letters, 2002, 29(3):1039-1045.

[11] 阮百堯, 熊彬, 徐世浙. 三維地電斷面電阻率測(cè)深有限元數(shù)值模擬[J]. 地球科學(xué), 2001, 26(1): 73-77.RUAN Bai-yao, XIONG Bin, XU Shi-zhe. Finite element method of modeling resistivity sounding on 3D geoelectic section[J]. Earth Science, 2001, 26(1): 73-77.

[12] 王燁. 基于矢量有限元的高頻大地電磁法三維數(shù)值模擬[D].長(zhǎng)沙: 中南大學(xué)信息物理工程學(xué)院, 2008: 35-62.WANG Ye. A study of 3D high frequency magnetotellurics modeling by edge-based finite element method[D]. Changsha:Central South University. School of Info-physics and Geometrics Engineering, 2008: 35-62.

[13] Zienkiewicz O C, Zhu J Z. Adaptive and mesh generation[J]. Int J Num Meth Eng, 1991, 32(4): 783-810.

[14] Tani K, Yamada T. H-version adaptive finite element method using edge element for 3D non-linear magnetostatic problems[J].IEEE Transactions on Magnetics, 1997, 33(2): 1756-1759.

[15] Szabo B A. Mesh design for the p-version of the finite element method[J]. Computer Methods in Applied Mechanics and Engineering, 1986, 55(1): 181-197.

[16] WEN Dai-gang, JIANG Ke-xun. P-version adaptive computation of FEM[J]. IEEE Transactions on Magnetics, 1994, 30(5):3515-3518.

[17] Key K, Weiss C. Adaptive finite-element modeling using unstructured grids: The 2D magnetotelluric example[J].Geophysics, 2006, 71(6): 291-299.

[18] LI Y G, Key K, Constable S. An adaptive finite element modeling of 2-D marine controlled-source electromagnetic fields[C]//18th IAGAWG 1.2 Workshop on Electromagnetic Induction in the Earth. El Vendrell, Spain, 2006: S3-S14.

[19] Harrington R F. Time harmonic electromagnetic fields[R]. New York: McGraw-Hill Book Co, 1961: 37-85.

[20] Nam M J, Kim H J, Song Y, et al. 3D magnetotelluric modeling including surface topography[J]. Geophysical Prospecting, 2007,55(2): 277-287.

[21] LIU Chang-sheng, REN Zheng-yong, TANG Jing-tian, et al.Three-dimensional magnetotellurics modeling using edge-based finite element unstructured meshes[J]. Applied Geophysics, 2008,5(3): 170-180.

[22] CHEN Zhi-ming, WANG Long, ZHENG Wei-ying. An adaptive multilevel method for time-harmonic Maxwell equations with singualarities[J]. Siam J Sci Comput, 2007, 29(1): 118-138.

[23] Zhdanov M S, Varentsov I M, Weaver J T, et al. Methods for modeling electromagnetic fields Results from COMMEMI—the international project on the comparison of modeling methods for electromagnetic induction[J]. Journal of Applied Geophysics,1997, 37(3/4): 265-271.