大地電磁三維矢量有限元正演模擬

石明+馮德山+李開鵬+王珣

摘 要：從Maxwell方程出發(fā)，開展了三維大地電磁場所滿足的邊值問題研究，利用加權余量法導出了三維大地電磁有限元方程.介紹了三維矢量有限元六面體網格剖分方式、插值基函數選取，推導了三維大地電磁矢量有限元正演的單元剛度系數矩陣及離散格式.編制了三維矢量有限元大地電磁正演的Matlab程序.三維COMMEMI 3D-1模型的視電阻率曲線與國際通用的標準測試數據能很好地擬合，驗證了作者編寫的矢量有限元正演程序的正確性.通過對高、低阻異常體的阻抗張量形態(tài)分析，說明張量阻抗等值線圖能用以大致判斷異常體特性，豐富了大地電磁響應特征的表達方式.

關鍵詞：矢量有限元；大地電磁；正演模擬；張量阻抗

中圖分類號：P631 文獻標識碼：A

文章編號：1674-2974（2016）10-0119-07

Abstract：Starting from the Maxwell equations， this article studied the boundary conditions of 3D MT. By using the weighted residual method， we derived the three-dimensional MT finite element equation. The three-dimensional vector finite element hexahedral meshing mode was introduced and the basis functions were selected. Then we derived the three-dimensional magnetotelluric vector finite element stiffness coefficient matrix and discrete format. A three-dimensional vector finite element magnetotelluric forward Matlab program was done. The apparent resistivity curve of the dimensional COMMEMI 3D-1 model matches the international standard test data， which proves the correctness of 3D magnetotelluric forward program. With the analysis of high and low resistivity anomalies， it shows that tensor impedance map can roughly determine the anomaly characteristics， which enriches the magnetotelluric response characteristics of expression.

Key words：vector finite element method；magnetotelluric； forward modeling； impedance tensor

大地電磁（MT）是以電離層激發(fā)的天然交變電磁場為場源，在地表觀測相互正交的電場、磁場分量來獲取地電構造信息的一種重要地球物理勘探方法[1].MT不需要龐大的發(fā)射源設備，只需采用比較輕便的接收設備，野外工作方便、成本低，被廣泛應用于地殼和上地幔電性結構的研究，在石油天然氣勘探、礦產資源勘探、工程與環(huán)境普查等領域，發(fā)揮著舉足輕重的作用[2-9].可以預見，三維MT勘探技術是地球物理中深層領域的研究熱點及今后MT的發(fā)展趨勢，而三維MT正演是理解MT勘探物理現象并認識地質體電磁響應規(guī)律的有效手段，顯然尤其重要.

盡管矢量FEM擁有諸多優(yōu)點，但在地球物理的電磁法正演領域中，其應用并不多見，尚需要進一步完善.目前的研究主要包括：Yoshimura 等 [10]開展了矢量FEM的MT響應數值模擬，并將矢量FEM的計算結果與交錯網格FDM的計算結果進行了對比；Mitsuhata等 [11]利用矢量FEM和節(jié)點FEM耦合的方法對三維MT數值模擬；Nam [12]采用不規(guī)則六面體矢量FEM直接計算電場，研究了起伏地形下MT的電阻率和相位的變化規(guī)律；劉長生等[13]將完全非結構化四面體單元引入到矢量有限元中，實現了三維大地電磁h-型自適應矢量有限元正演；王燁[14]開展了高頻率大地電磁法矢量有限元正演，并采用改進的威爾金森方法求解大型病態(tài)方程組，提高了迭代速度；顧觀文等[15]開展了矢量有限元法MT三維地形數值模擬，研究了地形起伏下三維阻抗張量的變化規(guī)律；楊軍等[16]采用非結構四面體單元的三維矢量FEM實現了海洋可控源電磁數值模擬；蘇曉波等[17]采用規(guī)則六面體單元的三維矢量FEM實現了大地電磁數值模擬，并對網格剖分的重要性進行了研究.

在前人基礎上，作者推導了三維大地電磁矢量FEM正演的離散形式，應用矢量FEM算法計算了三維COMMEMI 3D-1國際模型[18]的MT視電阻率及模型張量阻抗，研究了高低阻異常體的電磁響應特性，有效地指導了MT的資料解釋.

圖5（a），（b）分別為f=0.1 Hz時XY模式與YX模式下矢量FEM正演視電阻率曲線.分析圖5（a），（b）可知，兩幅圖中的矢量FEM曲線與COMMEMI所提供的數據都能夠很好地吻合，說明無論是在低頻還是高頻部分，應用矢量FEM開展三維大地電磁正演，都具有較高的精度，同時也驗證了矢量FEM算法及程序的正確性.

圖6為應用矢量FEM正演計算COMMEMI3D-1模型得到的張量阻抗.由圖可見，10 Hz與0.1 Hz兩個頻率下的張量阻抗形態(tài)基本一致，10 Hz的數值較0.1 Hz要大.對比圖中4個不同的張量阻抗，可以發(fā)現，圖6（a），（d），（e）和（h）中兩個頻率下的Zxx與Zyy分為四瓣，且阻抗值較小，四瓣的中心反映了異常體的邊界，而圖6（b），（c），（f）和（g）中的Zxy與Zyx阻抗值較大，反映了入射場的特性.根據張量阻抗理論可知，當構造為二維構造時，Zxx和Zyy為零，即當異常體走向方向越長，Zxy與Zyx越小，Zxy與Zyx差異也越大.由此，張量阻抗分解后，無需做反演即可以判斷出異常體的簡單特性.

為了進一步認識大地電磁的響應特性，對比高、低阻異常體張量阻抗的不同，在圖3中COMMEMI3D-1測試模型的基礎上，僅將低阻異常體改為1 000 Ω·m高阻異常體.其他參數均與國際模型相同.應用三維矢量FEM開展三維高阻異常體模型的張量阻抗研究.

圖7為應用三維矢量FEM正演的10 Hz大地電磁張量阻抗圖.分析圖7（a）與圖7（d）可知，高阻異常體張量阻抗中的Zxx與Zyy同樣分為四瓣，且阻抗值較小，其四瓣的中心反映了異常體的邊界.由于異常體x方向與y方向的比值為1∶2，圖7（b）中的Zxy與圖7（c）中的Zyx差異較大.對比高低阻異常10 Hz時的阻抗相位Zxx，雖然兩者都為四瓣，但是阻抗值正負值的分布正好相反，低阻異常體四瓣的中心向外輻射，幅值變小趨于0；而高阻異常體四瓣的中心向外輻射，幅值變小趨于0之后會發(fā)生反轉之后再次趨于0.Zyy具有相同的規(guī)律.對比Zxy和Zyx，高阻異常體中心僅出現一個閉合異常形態(tài)，而低阻異常體則形態(tài)更為復雜.

4 結 論

1） 介紹了三維矢量有限元區(qū)域剖分方式，對矢量FEM插值基函數以及單元插值方式進行了闡述，應用Galerkin算法，推導了三維矢量FEM大地電磁方程離散格式，編制了矢量FEM三維MT的Matlab模擬程序.

2） 設置三維COMMEMI 3D-1模型進行矢量FEM的計算，模擬結果與COMMEMI提供的數據擬合效果很好，驗證了矢量有限元程序的正確性.通過對比高低阻異常體的張量阻抗，分析了不同異常下張量阻抗的特點，進一步認識了MT的響應特性.

參考文獻

[1] 柳建新，童孝忠，郭榮文，等. 大地電磁測深勘探：資料處理反演與解釋[M]. 北京：科學出版社，2012：1-12.

LIU jian-xin， TONG Xiao-zhong， GUO Rong-wen， et al. Magnetotelluric sounding exploration： data processing inversion and interpretation [M].Beijing： Science Press， 2012：1-12.（In Chinese）

[2] 底青云，王若. 可控源音頻大地電磁數據正反演及方法應用[M]. 北京：科學出版社，2008：1-8.

DI Qing-yun， WANG Ruo. Controlled source audio magnetotelluric data inversion and methods application [M].Beijing： Science Press， 2008： 1-8.（In Chinese）

[3] 譚捍東，余欽范，JOHN B，等. 大地電磁法三維交錯采樣有限差分數值模擬[J]. 地球物理學報，2003，46（5）：705-711.

TAN Han-dong， YU Qin-fan， JOHN B， et al. Magnetotelluric three-dimensional modeling using the staggered-grid finite difference method[J]. Chinese Journal of Geophysics， 2003， 46（5）：705-711.（In Chinese）

[4] 陳輝，鄧居智，譚捍東，等. 大地電磁三維交錯網格有限差分數值模擬中的散度校正方法研究[J]. 地球物理學報，2011，54（6）：1649-1659.

CHEN Hui， DENG Ju-zhi， TAN Han-dong， et，al. Study on divergence correction method in three-dimensionalmagnetotelluric modeling with staggered-grid finite difference method [J]. Chinese Journal of Geophysics， 2011， 54（6）： 1649-1659. （In Chinese）

[5] 李焱，胡祥云，楊文采，等. 大地電磁三維交錯網格有限差分數值模擬的并行計算研究[J]. 地球物理學報，2012，55（12）：4036-4043.

LI Yan， HU Xiang-yun， YANG Wen-cai， et al. A study on parallel computation for 3D magnetotelluric modeling using the staggered-grid finite difference method [J]. Chinese Journal of Geophysics， 2012， 55（12）： 4036-4043. （In Chinese）

[6] 徐凌華，童孝忠，柳建新，等. 基于有限單元法的二維/三維大地電磁正演模擬策略[J]. 物探化探計算技術，2009，31（5）：421-425.

XU Ling-hua， TONG Xiao-zhong， LIU Jian-xin， et al. Solution strategies for 2D and 3D magnetotelluric forward modeling based on the finite element method [J]. Computing Techniques for Geophysical and Geochemical Exploration， 2009， 31 （5）： 421-425. （In Chinese）

[7] MOGI T. Three-dimensional modeling of magnetotelluric data using finite- element method [J]. Journal of Applied Geophysics， 1996， 35（2）： 185-189.

[8] 梁生賢，張勝業(yè)，吾守艾力，等. 復雜三維介質的大地電磁正演模擬[J]. 地球物理學進展，2012，27（5）：1981-1988.

LIANG Sheng-xian， ZHANG Sheng-ye， WU Shouaili， et al. Magnetotelluric forward modeling in complex three- dimensional media [J]. Progress in Geophys， 2012， 27（5）： 1981-1988. （In Chinese）

[9] 金建銘.電磁場有限元方法[M]. 西安：西安電子科技大學出版社， 1998：164-177.

JIN Jian-ming. Finite element method of electromagnetic field [M]. Xian： Xian Electronic University Press， 1998：164-177.（In Chinese）

[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-1042.

[11]MITSUHATA Y， UCHIDA T. 3D magnetotelluric modeling using the T-Ω finite-element method [J].Geophysics， 2004， 69（1）：108-119.

[12]NAM M J， KIM H J， SONG Y， et al. 3D magnetotelluric modelling including surface topography [J]. Geophysical Prospecting， 2007， 55（2）： 277-287.

[13]劉長生，湯井田，任政勇，等. 基于非結構化網格的三維大地電磁自適應矢量有限元模擬[J]. 中南大學學報：自然科學版，2010，41（5）：1855-1860.

LIU Chang-sheng， TANG Jing-tian， REN Zheng-yong， et al. Three-dimension magnetotellurics modeling by adaptive edgefinite-element using unstructured meshes [J]. Journal of Central South University： Science and Technology，2010， 41（5）： 1855-1860. （In Chinese）

[14]王燁. 基于矢量有限元的高頻大地電磁法三維數值模擬[D]. 長沙：中南大學地球科學與信息物理學院，2008：35-62.

WANG Ye. A study of 3D high frequency magnetotelluric modeling by edge-based finite element method [D]. Changsha： Central South University. School of Geosciences and Info-Physics， 2008：35-62. （In Chinese）

[15]顧觀文，吳文鸝，李桐林. 大地電磁場三維地形影響的矢量有限元數值模擬[J]. 吉林大學學報：地球科學版，2014，44（5）：1678-1686.

GU Guan-wen， WU Wen-li， LI Tong-lin. Modeling for the effect of magnetotelluric 3D topography based on the vector finite-element method [J]. Journal of Jilin University： Earth Science Edition， 2014，44（5）：1678-1686. （In Chinese）

[16]楊軍，劉穎，吳小平. 海洋可控源電磁三維非結構矢量有限元數值模擬[J]. 地球物理學報，2015，58（8）：2827-2838.

YANG Jun， LIU Ying， WU Xiao-ping. 3D simulation of marine CSEM using vector finite element method on unstructured grids [J]. Chinese Journal of Geophysics， 2015， 58（8）： 2827-2838. （In Chinese）

[17]蘇曉波，李桐林，朱成，等.大地電磁三維矢量有限元正演研究[J].地球物理學進展，2015，30（4）：1772-1778.

SU Xiao-bo， LI Tong-lin， ZHU Cheng， et al. Study of three-dimensional MT forward modeling using vector finite element method [J].Progress in Geophysics， 2015，30（4）：1772-1778. （In Chinese）

[18]ZHDANOV M S， VARENTSOV I M， WEAVER J T， et al. Methods for modelling electromagnetic fields Results from COMMEMI—the international project on the comparison of modelling methods for electromagnetic induction [J]. Journal of Applied Geophysics， 1997， 37（3/4）：133-271.