基于非結構化網格的三維大地電磁自適應矢量有限元模擬

劉長生 ，湯井田，任政勇 ，馮德山

(1. 中南大學 信息物理工程學院，湖南 長沙，410083；2. 長沙航空職業技術學院 計算機系，湖南 長沙，410014；3. 瑞士聯邦理工學院(EIH) 地球物理系，瑞士 蘇黎士，CH8092)

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

1 三維MT矢量有限元數學模型

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

圖1 三維MT模型的一般地電結構Fig.1 MT model of three-dimensional structure of general earth electricity

在空氣層和地下地層區域，電場滿足下面微分控制方程[19]：

式中：E為電場強度；i為虛數；ω為角頻率；μ為自由空間磁導率；σ為MT模型的電導率。

采用理論可靠和執行簡單的Direlet邊界條件[20]，并采用水平地層地電結構在邊界上的電磁場近似作為邊界條件。

式中：n為外邊界的外單位法向量；E0為外邊界上給定的已知電場強度。

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

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

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

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

MT 模型求解的大型復線性方程為：

其中：U為定義在各個四面體單元的邊上電場切向未知數向量。矩陣A和B由各個單元的貢獻所得[21]：

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

2 基于殘差的后驗誤差估計

采用Nédélec四面體線性單元，其求解過程可以表述為：定義Tk為計算區域Ω的一系列四面體網格剖分單元，Fk為網格中的一系列三角面；定義Tk上的有限元空間Uk為[21]：

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

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

定義數值解的誤差為：

經推導，適合于MT模型的后驗誤差估計為：

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

3 自適應策略

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

在四面體網格Tk上，利用有限元方法求出Ek，根據文獻[22]中的誤差收斂理論，得出當前網格上的誤差估計，并得到單元上的局部誤差估計指示值：

相應地，全局誤差指示值和最大誤差估計指示值分別為：

當計算單元誤差超過給定的誤差限后，采用相應的加密策略，具體算法流程如圖2所示。

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

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

4 COMMEMI 3D-1 MT模型數值模擬算例

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

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

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

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

注意到測線附近的單元并沒有隨著自適應加密過程的運行而對網格顯著加密。這一現象導致測線附近數值精度提高(見圖3)：視電阻率的平均誤差從初始網格上的3.49%減低為第12次網格上的0.87%，但是，隨著網格的再加密，視電阻率平均誤差始終保持在0.50%～1.00%；相位曲線的平均誤差從初始網格上的1.23%減低為第20次時的0.21%。但是，隨著網格的再加密，相位曲線平均誤差始終保持在0.20%～0.50%。

5 結論

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

(2) 推導出基于殘差的三維大地電磁矢量有限元后驗誤差估計公式，為三維大地電磁自適應矢量有限元數值模擬的實現奠定了基礎。

(3) 在完全非結構化四面體單元剖分及優化基礎上，結合三維大地電磁矢量有限元后驗誤差估計公式，提出了基于非結構化網格的三維大地電磁h型自適應矢量有限元計算策略，保證了對復雜大地電磁模型數值計算的精度和可靠性。

(4) 電磁模型的復雜性總體上不影響方法的收斂性，且可以達到預期的計算精度。可見，h型自適應矢量有限元可以保證復雜模型的計算精度和速度，具有廣闊的應用前景。

[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] 李大潛. 有限元素法在電法測井中的應用[M]. 北京: 石油工業出版社, 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] 羅延鐘, 張桂青. 電子計算機在電法勘探中的應用[M]. 武漢:武漢地質學院出版社, 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]. 北京: 科學出版社,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] 阮百堯, 熊彬, 徐世浙. 三維地電斷面電阻率測深有限元數值模擬[J]. 地球科學, 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] 王燁. 基于矢量有限元的高頻大地電磁法三維數值模擬[D].長沙: 中南大學信息物理工程學院, 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.