石 明, 馮德山, 李開鵬，王 珣

(1.中南大學(xué) 地球科學(xué)與信息物理學(xué)院，湖南 長沙，410083；2. 有色資源與地質(zhì)災(zāi)害探查湖南省重點(diǎn)實(shí)驗(yàn)室, 湖南 長沙，410083；3. 貴州省有色金屬和核工業(yè)地質(zhì)勘查局物化探總隊(duì)，貴州 都勻，558004)

?

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

石 明1,3, 馮德山1,2*, 李開鵬3，王 珣1,2

(1.中南大學(xué) 地球科學(xué)與信息物理學(xué)院，湖南 長沙，410083；2. 有色資源與地質(zhì)災(zāi)害探查湖南省重點(diǎn)實(shí)驗(yàn)室, 湖南 長沙，410083；3. 貴州省有色金屬和核工業(yè)地質(zhì)勘查局物化探總隊(duì)，貴州 都勻，558004)

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

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

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

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

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

1 三維大地電磁邊值問題

以e-iwt表示諧變場(chǎng)的時(shí)間因子，大地電磁滿足的Maxwell方程組可表示為：

(1)

式中：ω為角頻率；μ為介質(zhì)的磁導(dǎo)率；σ為電導(dǎo)率；ε為介電常數(shù).聯(lián)立兩式，消去H，并假定巖石和空氣中的μ為常數(shù)得：

(2)

式(2)為電場(chǎng)矢量E所應(yīng)滿足的微分方程,式中k為波數(shù),k2=iωμσ+ω2εμ.將電場(chǎng)E寫成分量形式為：

E=Exex+Eyey+Ezez．

(3)

在三維情況下，Ex，Ey和Ez是相互聯(lián)系的，無法將其中一個(gè)分量單獨(dú)分離出來.根據(jù)廣義變分原理，將微分方程和邊界條件都考慮在內(nèi)，利用加權(quán)余量法[9]推導(dǎo)得到的三維大地電磁邊值問題相關(guān)的加權(quán)余量方程為：

(4)

2 矢量FEM網(wǎng)格離散及單元分析

矢量FEM與常規(guī)節(jié)點(diǎn)FEM求解過程非常類似，僅僅在區(qū)域剖分、插值基函數(shù)選取及單元分析方面略有不同，主要步驟包括：(ⅰ) 區(qū)域剖分；(ⅱ) 選用矢量插值基函數(shù)；(ⅲ) 單元分析；(ⅳ) 總體合成；(ⅴ) 多元函數(shù)求偏導(dǎo)數(shù)；(ⅵ) 解線性代數(shù)方程組.

圖1 矢量有限元三維區(qū)域剖分圖

圖2 矢量六面體單元中的節(jié)點(diǎn)與棱邊編號(hào)

通過分配常切向場(chǎng)分量給單元的每一條邊，單元內(nèi)x,y,z3個(gè)方向上場(chǎng)分量可分別表示為：

(5)

其中：

(6)

矢量FEM定義的基函數(shù)具有零散度和非零旋度的特性.而且，構(gòu)成離散單元各個(gè)小平面上的切向場(chǎng)僅由組成小平面的棱邊上的切向場(chǎng)決定，不僅保證了穿越棱邊的切向場(chǎng)連續(xù)，而且保證了穿越離散單元表面時(shí)切向場(chǎng)的連續(xù)性.用六面體單元對(duì)整個(gè)區(qū)域進(jìn)行剖分.單元的坐標(biāo)變換關(guān)系為：

x=xc+ξ·a/2,

y=yc+η·b/2,

z=zc+ζ·c/2．

(7)

式中：xc，yc，zc是子單元中心的坐標(biāo)，微分關(guān)系為：

(8)

將式(4)中的區(qū)域積分分解為各單元積分之和：

(9)

由插值基函數(shù)形式，可得到

(10)

δEx,δEy,δEz分別為電場(chǎng)Ex,Ey,Ez的變分.將式(10)帶入，可得到式(9)中的第一項(xiàng)積分為：

∫e×δE·×EdΩ=∫eδET·(×NT)×

(11)

(12)

將形函數(shù)代入并進(jìn)行積分求解，可以得到式(9)第一積分項(xiàng)的三維單元?jiǎng)偠染仃囅禂?shù)：

其中N1，N2，N3均為經(jīng)過三重積分后得到的4×4的單元矩陣，結(jié)果如下：

將形函數(shù)代入式(9)的第二項(xiàng)單元積分中：

(13)

式中：

K11=K22=K33=∫ek2NNTdΩ，

(14)

求解式(14)，可以得到式(9)第二項(xiàng)的單元?jiǎng)偠染仃囅禂?shù)為：

(15)

當(dāng)單元的邊界面5678落在區(qū)域底面邊界EFGH上時(shí)，邊界積分

(16)

式中：

K11=∫χNNTdΓ，

(17)

這里Ni，Nj為四邊形單元的矢量基函數(shù)，不是六面體單元的基函數(shù).求解式(17)，可以得到式(9)第3積分項(xiàng)的三維矢量FEM單元?jiǎng)偠染仃囅禂?shù)：

(18)

在一個(gè)單元內(nèi)，將三項(xiàng)積分得出的矩陣K1e，K2e，K3e相加，再將單元的系數(shù)矩陣擴(kuò)展成由全體棱邊組成的系數(shù)矩陣，再將各單元相加，對(duì)(9)式進(jìn)行離散、擴(kuò)展，將頂界面ABCD上的邊界條件代入線性方程組后，解方程組，就可以得到各個(gè)棱邊的Ex,Ey,Ez的值.根據(jù)式(1)中第一式，進(jìn)行有限元求解即可得到相應(yīng)的Hx，Hy，Hz的值.

3 三維矢量FEM大地電磁正演模擬

為了驗(yàn)證三維矢量FEM程序的正確性和精度，選取國際通用COMMEMI 3D-1 標(biāo)準(zhǔn)測(cè)試模型，模型示意圖如圖3所示.

圖3中異常體大小為1 000 m×2 000 m×2 000 m，異常體上界面據(jù)地面250 m，電阻率為0.5·m，背景電阻率為100·m，MT測(cè)線分別布置在x，y坐標(biāo)軸上，范圍為0 km～2.5 km.其中x，y，z方向網(wǎng)格數(shù)為 41×41×39(8 層空氣層)，計(jì)算區(qū)域?yàn)?0 000 m×25 000 m×35 000 m， 空氣層厚度為 5 000 m，空氣層電阻率設(shè)為 1015·m .

圖3 COMMEMI 3D-1標(biāo)準(zhǔn)測(cè)試模型

圖4為f=10 Hz時(shí)XY模式與YX模式下矢量FEM正演模擬所得到的視電阻率曲線.VFEM3D表示三維矢量FEM計(jì)算結(jié)果，黑色豎線為COMMEMI提供的誤差值.由圖中可見，f=10 Hz時(shí)矢量FEM的計(jì)算結(jié)果與COMMEMI提供的結(jié)果盡管存在細(xì)小誤差，但相差不大，且在誤差允許范圍之內(nèi)，表明在該頻率下矢量FEM的計(jì)算結(jié)果可信.在圖4(a)XY模式下矢量FEM計(jì)算結(jié)果較COMMEMI偏小，特別是第一個(gè)點(diǎn)已經(jīng)超出誤差棒的范圍.這一現(xiàn)象普遍存在，推斷為矢量FEM后處理精度不足引起的.而圖4(b)YX模式下VFEM3D結(jié)果與COMMEMI結(jié)果吻合相當(dāng)好.

x/mm

x/mm

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

x/mm

x/mm

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

圖6 COMMEMI 3D-1A模型單個(gè)低阻異常體張量阻抗等值線圖

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

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

圖7 COMMEMI 3D-1A模型單個(gè)高阻異常體張量阻抗等值線圖

4 結(jié) 論

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

2) 設(shè)置三維COMMEMI 3D-1模型進(jìn)行矢量FEM的計(jì)算，模擬結(jié)果與COMMEMI提供的數(shù)據(jù)擬合效果很好，驗(yàn)證了矢量有限元程序的正確性.通過對(duì)比高低阻異常體的張量阻抗，分析了不同異常下張量阻抗的特點(diǎn)，進(jìn)一步認(rèn)識(shí)了MT的響應(yīng)特性.

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

LIU jian-xin, TONG Xiao-zhong, GUO Rong-wen,etal. Magnetotelluric sounding exploration: data processing inversion and interpretation [M].Beijing: Science Press, 2012:1-12.(In Chinese)

[2] 底青云，王若. 可控源音頻大地電磁數(shù)據(jù)正反演及方法應(yīng)用[M]. 北京：科學(xué)出版社，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，等. 大地電磁法三維交錯(cuò)采樣有限差分?jǐn)?shù)值模擬[J]. 地球物理學(xué)報(bào)，2003，46(5)：705-711.

TAN Han-dong, YU Qin-fan, JOHN B,etal. Magnetotelluric three-dimensional modeling using the staggered-grid finite difference method[J]. Chinese Journal of Geophysics, 2003, 46(5):705-711.(In Chinese)

[4] 陳輝，鄧居智，譚捍東,等. 大地電磁三維交錯(cuò)網(wǎng)格有限差分?jǐn)?shù)值模擬中的散度校正方法研究[J]. 地球物理學(xué)報(bào)，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] 李焱，胡祥云，楊文采,等. 大地電磁三維交錯(cuò)網(wǎng)格有限差分?jǐn)?shù)值模擬的并行計(jì)算研究[J]. 地球物理學(xué)報(bào)，2012，55(12)：4036-4043.

LI Yan, HU Xiang-yun, YANG Wen-cai,etal. 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]. 物探化探計(jì)算技術(shù)，2009，31(5)：421-425.

XU Ling-hua, TONG Xiao-zhong, LIU Jian-xin,etal. 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è)，吾守艾力,等. 復(fù)雜三維介質(zhì)的大地電磁正演模擬[J]. 地球物理學(xué)進(jìn)展，2012，27(5)：1981-1988.

LIANG Sheng-xian, ZHANG Sheng-ye, WU Shouaili,etal. Magnetotelluric forward modeling in complex three- dimensional media [J]. Progress in Geophys, 2012, 27(5): 1981-1988. (In Chinese)

[9] 金建銘.電磁場(chǎng)有限元方法[M]. 西安：西安電子科技大學(xué)出版社, 1998:164-177.

JIN Jian-ming. Finite element method of electromagnetic field [M]. Xi’an: Xi’an 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 theT-Ωfinite-element method [J].Geophysics, 2004, 69(1):108-119.

[12]NAM M J, KIM H J, SONG Y,etal. 3D magnetotelluric modelling including surface topography [J]. Geophysical Prospecting, 2007, 55(2): 277-287.

[13]劉長生，湯井田，任政勇,等. 基于非結(jié)構(gòu)化網(wǎng)格的三維大地電磁自適應(yīng)矢量有限元模擬[J]. 中南大學(xué)學(xué)報(bào):自然科學(xué)版，2010，41(5)：1855-1860.

LIU Chang-sheng, TANG Jing-tian, REN Zheng-yong,etal. 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]王燁. 基于矢量有限元的高頻大地電磁法三維數(shù)值模擬[D]. 長沙：中南大學(xué)地球科學(xué)與信息物理學(xué)院,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]顧觀文，吳文鸝，李桐林. 大地電磁場(chǎng)三維地形影響的矢量有限元數(shù)值模擬[J]. 吉林大學(xué)學(xué)報(bào):地球科學(xué)版,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]楊軍，劉穎，吳小平. 海洋可控源電磁三維非結(jié)構(gòu)矢量有限元數(shù)值模擬[J]. 地球物理學(xué)報(bào)，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].地球物理學(xué)進(jìn)展，2015，30(4)：1772-1778.

SU Xiao-bo, LI Tong-lin, ZHU Cheng,etal. 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,etal. 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.

Three-dimensional Magnetotelluric Forward Modeling Using Vector Finite Element Method

SHI Ming1,3，F(xiàn)ENG De-shan1,2?，LI Kai-peng3，WANG Xun1,2

(1．School of Geosciences and Info-Physics，Central South Univ，Changsha, Hunan 410083, China； 2. Key Laboratory of Non-ferrous Resources and Geological Detection of Hunan Province, Changsha, Hunan 410083, China； 3. Geophysical and Geochemical Prospecting Team, Non-ferrous Metals and Nuclear Industry Geological Exploration Bureau of Guizhou Province, Duyun, Guizhou 558004, China)

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.

vector finite element method；magnetotelluric; forward modeling; impedance tensor

1674-2974(2016)10-0119-07

2016-03-28

國家自然科學(xué)基金資助項(xiàng)目(41574116), National Natural Science Foundation of China(41574116)；中南大學(xué)創(chuàng)新驅(qū)動(dòng)項(xiàng)目(2015CX008)；中南大學(xué)教師研究基金資助項(xiàng)目(2014JSJJ001); 中南大學(xué)升華育英人才計(jì)劃(2012)；湖湘青年創(chuàng)新創(chuàng)業(yè)平臺(tái)培養(yǎng)對(duì)象共同資助項(xiàng)目(2013)

石 明(1969-)，男，湖南益陽人，中南大學(xué)博士研究生，高級(jí)工程師

?通訊聯(lián)系人，E-mail: fengdeshan@126.com

P631

A