電性參數(shù)分塊連續(xù)變化二維MT有限元數(shù)值模擬

劉 云，王緒本

油氣藏地質(zhì)及開發(fā)工程國家重點實驗室，成都理工大學(xué)，成都 610059

電性參數(shù)分塊連續(xù)變化二維MT有限元數(shù)值模擬

劉 云，王緒本＊

油氣藏地質(zhì)及開發(fā)工程國家重點實驗室，成都理工大學(xué)，成都 610059

為了易于模擬野外復(fù)雜地形和地下任意形狀地電體模型，將有限元單元網(wǎng)格設(shè)計為三角單元；并考慮到野外實際勘探中，地球介質(zhì)的電性參數(shù)均是連續(xù)變化的情況，單元內(nèi)的場值和電性參數(shù)被設(shè)計為雙線性變化；推導(dǎo)出二維起伏地形條件下大地電磁法有限元數(shù)值模擬算法；根據(jù)單元節(jié)點主場值和線性插值形函數(shù)間的關(guān)系，計算出單元節(jié)點的輔助場值；在二維起伏地形情況下，定義TE、TM模式視電阻率和阻抗相位.4個模型的計算的結(jié)果與解析法的均方根誤差小于1%，地形模擬與前人的計算結(jié)果相符，模擬傾斜界面異常體，能有效的反映出其異常形態(tài).

大地電磁法，起伏地形，三角單元網(wǎng)格，電性參數(shù)分塊連續(xù)變化，有限元

1 引 言

大地電磁法（MT）已廣泛應(yīng)用于地球物理探測各個領(lǐng)域.我國的山區(qū)面積比較大，在起伏地形條件下開展MT工作時，地形的影響產(chǎn)生電磁場畸變，這給資料處理和解釋工作帶來很大困難，因此有必要研究帶地形的MT數(shù)值模擬算法.

目前，在二維起伏地形和地電模型條件下，國內(nèi)外學(xué)者采用過多種有限元剖分方式的二維MT有限元數(shù)值模擬.在有限元電性參數(shù)（如電導(dǎo)率、傳播系數(shù) 等）分 塊 均 勻 的 前 提 下，Wanna maker[1]、Mauriello等[2]采用矩形三角網(wǎng)格剖分、雙線性插值有限元進行數(shù)值模擬，陳小斌[3－4]、謝飛[5]用自適應(yīng)地形三角網(wǎng)格剖分、雙線性插值有限元數(shù)值模擬，周熙襄[6]、王緒本[7]、劉云等[8]采用自適應(yīng)地形四邊形網(wǎng)格剖分、雙二次插值有限元數(shù)值模擬，對斜坡型地形產(chǎn)生的地形影響做模擬和分析；認為電性參數(shù)分塊均勻的模擬方式，在單元塊內(nèi)的電性參數(shù)分布是均勻的，單元塊與單元塊之間的電性參數(shù)分布則不一定是連續(xù)的.但是在野外實際中，巖石、礦物的電性參數(shù)往往是變化的，因此采用電性參數(shù)分塊連續(xù)變化的數(shù)值模擬方法，更符合野外實際.徐世浙、李予國等[9－10]采用矩形網(wǎng)格剖分，進行電性參數(shù)分塊連續(xù)變化、雙線性插值、雙二次插值等有限元數(shù)值模擬；對于矩形網(wǎng)格剖分，如果在有限元矩形單元塊內(nèi)或單元節(jié)點處填入空氣介質(zhì)，從而形成臺階型地形模擬地形影響[11－12]，是有其局限性的[4，13].

阮百堯[14]在起伏地形二維直流電阻率法正演數(shù)值模擬中，采用矩形內(nèi)三角網(wǎng)格剖分、電導(dǎo)率分塊連續(xù)變化、雙線性插值有限元數(shù)值模擬.本文在此基礎(chǔ)上，導(dǎo)出了起伏地形三角單元網(wǎng)格剖分電性參數(shù)分塊連續(xù)變化的二維MT有限元數(shù)值模擬方法.

2 二維MT變分問題

如圖1所示，假設(shè)三維地質(zhì)模型的電、磁等參數(shù)沿構(gòu)造方向（x方向）上無變化，這時可將三維模型問題轉(zhuǎn)化為二維模型問題進行處理.此時，與二維MT邊值問題相應(yīng)的變分問題[13]是

圖1 區(qū)域Ω有限元網(wǎng)格剖分Fig.1 Division of region Ω with FEM grids

式（1）中，u為Ex或Hx，Ω為二維研究區(qū)域.

對于TE模式（Ex型）：上邊界AB處離地面足夠遠.則

對于TM模式（Hx型）：上邊界AB在地面上.則

其中，σ和μ分別為介質(zhì)電導(dǎo)率和磁導(dǎo)率，ω為圓頻率，ε為介電常數(shù).下邊界CD處離目標區(qū)域足夠遠，這時σ為CD 處均勻介質(zhì)電導(dǎo)率；左、右邊界AC、BD離目標區(qū)域足夠遠.

3 有限單元法

用有限單元法求解式（1）的二維 MT變分問題，具體步驟如下：

3.1 網(wǎng)格剖分

對于TE模式，在高頻時，空氣介質(zhì)中位移電流的影響不可忽略[3]，因此研究區(qū)域為空氣和地下；對于TM模式，場值Hx在空氣中近似為常數(shù)（在y和z方向的偏導(dǎo)數(shù)近似為零），即不考慮空氣介質(zhì)中位移電流的影響，此時AB邊界設(shè)在地面上，研究區(qū)域為地下.

如圖1所示，研究區(qū)域Ω用三角單元進行網(wǎng)格剖分.如圖1a所示，當模擬水平地形時，考慮到地下介質(zhì)體形狀的任意性，首先對研究區(qū)域進行矩形網(wǎng)格剖分，再在每一個矩形網(wǎng)格內(nèi)剖分出4個三角單元網(wǎng)格.這樣，一方面，避免了三角網(wǎng)格過于尖銳的情況；另一方面，可利用三角形的斜邊模擬地形線和任意介質(zhì)體傾斜的界面；如圖1b所示，當模擬起伏地形時，在圖1a網(wǎng)格剖分的基礎(chǔ)上，地空邊界采用沿實際地形線進行網(wǎng)格剖分.

3.2 線性插值

如圖2所示，假設(shè)三角單元e內(nèi)電磁場u和電性參數(shù)τ，λ，k呈線性變化，則在每一個單元中：

圖2 三角單元eFig.2 Triangle element e

這里，ui為三角單元節(jié)點處的電磁場，τi、λi、ki為三角單元節(jié)點處的電性參數(shù)，（i）、1、2和3為三角單元的三角節(jié)點號.是關(guān)于y和z的 線 性 插 值 形 函 數(shù) （i＝1，2，3），其 中Δ ＝是三角單元的面積，其中

其中，（y1，z1）、（y2，z2）和（y3，z3）是三角單元節(jié)點號的坐標.

3.3 單元分析

在整個研究區(qū)域Ω中，將泛函（1）式離散化，表示為所有三角單元e的線性組合，即

3.3.1 單元分析1

u對y求偏導(dǎo)數(shù)，有

3.3.2 單元分析2

因此，對于（2）式中第二個積分項，有

其中，c11＝ （6 22）λ，c12＝ （2 21）λ，c13＝（2 12）λ，c21＝ （2 21）λ，c22＝ （2 62）λ，c23＝ （1 22）λ，c31＝ （2 12）λ，c32＝ （1 22）λ，c33＝ （2 26）λ，λ＝ （λ1λ2λ3）T.

3.3.3 單元分析3

3.4 系數(shù)矩陣總體合成及求變分

將研究區(qū)域Ω內(nèi)每個單元系數(shù)矩陣ke按照總體節(jié)點號進行擴展，得到，相加得到總體系數(shù)矩陣K，即

如圖3所示，采用矩形內(nèi)三角網(wǎng)格剖分方式，比較矩形網(wǎng)格剖分來說，增加了一個中間節(jié)點5，大大增加了計算量.但是節(jié)點5只與該矩形內(nèi)4個角點有關(guān)系，而與其它不相鄰的矩形網(wǎng)格沒有任何直接聯(lián)系，所以在單元分析中，運用高斯消元法事先可消去這個節(jié)點.節(jié)點5只是一個虛設(shè)的節(jié)點，并不包含在總體節(jié)點數(shù)中.到解方程結(jié)束后，亦可根據(jù)矩形網(wǎng)格內(nèi)4個角點的已知場值和節(jié)點系數(shù)之間的關(guān)系，直接計算出節(jié)點5的場值.這樣，一方面不增加節(jié)點總數(shù)，節(jié)省了計算量，另一方面每個三角單元仍具有

圖3 矩形內(nèi)節(jié)點間的關(guān)系Fig.3 Relationship of nodes in the rectangle

各自的不同物性.具體做法為：節(jié)點5的消除，將ki，j（i，j ＝ 1，2，3，4，5）替 換 為節(jié) 點 5 的 恢 復(fù)，u5＝

4 輔助場、視電阻率和阻抗相位的計算

4.1 輔助場的計算

輔助場是通過計算主場值沿某方向的方向?qū)?shù)得到的[17－19].當各節(jié)點主場值求出后，由3.2節(jié)可知，單元內(nèi)場值u可表示為則u對y求偏導(dǎo)數(shù)，得

其中a1、a2、a3，是與三角單元的三頂點坐標有關(guān)的常數(shù)，所以在單元內(nèi)，u的偏導(dǎo)數(shù)是常數(shù).

同理可得

在大地電磁法的數(shù)值模擬中，常以地面三角單元邊上任意點或三角單元的節(jié)點為實際測量點，此時要考慮到相鄰單元間場值偏導(dǎo)數(shù)的不連續(xù)性問題[13].對于彼此相鄰的單元、相同單元邊或節(jié)點上場值的偏導(dǎo)數(shù)，可以在所屬不同單元中予以分別計算，再取平均值.本文以地面三角單元的節(jié)點為實際測量點.不難發(fā)現(xiàn)，對于TE模式，需要計算同一測點的地表和地表上一空氣層所有單元的場值偏導(dǎo)數(shù)，再取平均值；而對于TM模式，則只需求出同一測點的地表所有單元的場值偏導(dǎo)數(shù)，再取平均值即可.

4.2 視電阻率定義和阻抗相位的計算

對于TE模式：u＝Ex，輔助場在二維起伏地形條件下，Ex、Hy在野外實際中均可測得，因此，視電阻率的定義仍采用普遍形式的阻抗定義方式，即

5 模型計算

為驗證本文正演算法的有效性，首先與水平層狀一維連續(xù)介質(zhì)理論模型的解析法、以及前人已有地形模型做計算分析，之后用該算法對傾斜界面異常體模型做計算分析.

在以下所有模型的計算中，有限元網(wǎng)格的剖分方法如下：如圖1所示，在橫向網(wǎng)格中，目標區(qū)域采用等間隔網(wǎng)格，間隔距等于點距.左、右稀疏網(wǎng)格數(shù)各為18個，網(wǎng)格大小等于點距乘以系數(shù)，其中18個系數(shù)依次分別為1，1，1，1，1，2，2，2，2，4，4，4，8，8，16，32，64，128；在縱向網(wǎng)格中，空氣稀疏網(wǎng)格坐標（TM模式無空氣網(wǎng)格，單位為m）14個，依次分別為100000，50000，10000，5000，1000，500，200，100，55，35，25，15，10，5.地下網(wǎng)格的剖分，從1～10層的網(wǎng)格大小（單位為 m）依次分別為5，5，5，5，10，10，10，15，15，20，第1層網(wǎng)格大小的默認值為5m.當然，也可以任意設(shè)置第1層網(wǎng)格大小，這樣，從1～10層的網(wǎng)格大小則相應(yīng)地等比例放大或縮小.從11層及以下的網(wǎng)格剖分為等比網(wǎng)格，該層的網(wǎng)格大小等于上一層網(wǎng)格大小乘以等比系數(shù)，等比系數(shù)一般為1.1或1.2.這樣，橫向網(wǎng)格數(shù)則取決于測點數(shù)的多少，縱向網(wǎng)格數(shù)則取決于探測深度的選擇；當模擬起伏地形時，上、下地表附近還要加上地形網(wǎng)格，地形網(wǎng)格數(shù)等于所有測點的高程數(shù)（對于重復(fù)的高程，只計1個），地形網(wǎng)格大小等于相鄰高程之差.

5.1 一維連續(xù)介質(zhì)模型

如圖4所示，三層一維連續(xù)介質(zhì)地電斷面模型，第一層電阻率為100Ωm、層厚度為1000m；第二層為連續(xù)介質(zhì)層，介質(zhì)電阻率隨深度線性減小，其變化關(guān)系如圖中所示，層厚度為1000m；第三層電阻率為1Ωm；測點區(qū)域為0～4km，點距為100m，共41個測點，頻率范圍為103～10－3Hz，以10為底的對數(shù)間隔采樣，共51個頻點.

圖4 三層地電斷面模型Fig.4 Schematic diagram for geoelectric section of model 3 layers

在一維層狀介質(zhì)的解析法中，將中間的連續(xù)介質(zhì)層剖分為10層進行模擬.在二維有限元數(shù)值模擬中，取測線中心處第21號點測點為研究對象.則解析法與二維數(shù)值解法的結(jié)果對比如圖5所示.

圖5 有限元數(shù)值解與解析法結(jié)果比較Fig.5 Comparison between analytical and numerical solutions by FEM

在TE、TM視電阻率曲線圖和阻抗相位曲線圖上，解析法和數(shù)值解曲線形態(tài)、值的大小基本一致，這表明本文方法能夠有效地模擬連續(xù)介質(zhì)模型.正演結(jié)果的數(shù)據(jù)統(tǒng)計表明，視電阻率、阻抗相位的解析法結(jié)果和數(shù)值模擬結(jié)果均方根誤差均小于1%.本文的數(shù)值模擬結(jié)果與解析法結(jié)果之間的誤差分析主要為：一方面，與有限元數(shù)值模擬的網(wǎng)格剖分有關(guān)系，如高頻網(wǎng)格和低頻網(wǎng)格的剖分方式是不一樣的；另一方面，解析法是基于電性參數(shù)分層均勻的計算方式，而本文的數(shù)值模擬方法是基于電性參數(shù)分塊連續(xù)變化的模擬.

5.2 地形模型

如圖6下部所示為圓柱凹陷地形模型，Wanna maker等[1]和徐世浙等[20]在文獻中對這個模型做過模擬.測點區(qū)域為0～0.4km，點距為10m，共41個測點，圓半徑為50m，頻率為0.01Hz，計算TM模式下視電阻率.

本文數(shù)值模擬結(jié)果如圖6上部曲線圖所示，對照文獻[1，20]中模擬結(jié)果，在0.01Hz頻點處，TM模式下，三種模擬方法與地形測點對應(yīng)的視電阻率起伏形態(tài)和值的大小基本一致.

圖6 圓柱凹陷地形模型與TM模式計算結(jié)果Fig.6 Topography of columnar sunken model and Result of TM

圖7所示為山脊地形模型，在文獻[1，20]中對這個模型做過模擬.測點區(qū)域為0～4km，點距為50m，共81個測點，起伏落差為100m，山底寬度為2.4km，頻率為10Hz.

圖7 山脊地形模型Fig.7 Chine terrain model

本文數(shù)值模擬結(jié)果如圖8所示.對照文獻[1，20]中的模擬結(jié)果，在10Hz頻點處，TE、TM模式下，三種模擬方法與地形測點對應(yīng)的視電阻率、阻抗相位起伏形態(tài)和值的大小基本一致.

通過對以上兩例均勻大地起伏地形模型的模擬，得到了與前人研究相符的結(jié)論.分塊均勻法認為地面或地下單元內(nèi)節(jié)點之間的電性參數(shù)是均勻分布的，而本文的分塊連續(xù)法則認為地面或地下單元內(nèi)節(jié)點之間的電性參數(shù)是連續(xù)變化的（空氣介質(zhì)除外）.顯然，分塊均勻法是分塊連續(xù)法的特殊形表現(xiàn)形式，分塊連續(xù)法更具有普遍性.用三角網(wǎng)格模擬起伏地形、分塊連續(xù)變化的剖分方式且更加符合野外實際地質(zhì)情況.

5.3 傾斜界面異常體模型

如圖9所示為“W”字形傾斜界面異常體模型，背景電阻率為100Ωm，異常體電阻率為20Ωm，其上頂面距地面100m，下底面距地面500m，異常體上頂面外寬4800m，下底面外寬2400m.測點區(qū)域為0～8km，點距為100m，共81個測點，頻率為103～10－1Hz，以10為底的對數(shù)間隔采樣，共41個頻點.

圖9 傾斜界面異常體模型Fig.9 Model of slopeing interface abnormity body

有限元數(shù)值模擬網(wǎng)格剖分為：最大探測深度選擇為10km，等比網(wǎng)格系數(shù)選擇為1.1；TE模式，網(wǎng)格數(shù)為116×66（橫向網(wǎng)格數(shù)為116，左、右稀疏網(wǎng)格數(shù)各為18；縱向網(wǎng)格數(shù)為66，其中包括空氣稀疏網(wǎng)格數(shù)為14）；TM模式，網(wǎng)格數(shù)為116×52（橫向網(wǎng)格數(shù)同TE模式；縱向網(wǎng)格數(shù)為52，沒有空氣網(wǎng)格）.大型稀疏矩陣的存儲方式均按文獻[13，15]中變帶寬方式存儲，用不帶平方根的Cholesky分解法求解線性方程組[13，16].

本文數(shù)值模擬方法的結(jié)果剖面圖如圖10所示，計算每一個頻點的時間為0.78s（TE和TM模式同時計算）.在TE、TM兩種極化模式下，視電阻率、阻抗相位剖面上都出現(xiàn)明顯低阻異常，“W”形異常形態(tài)清晰可見，表明采用三角網(wǎng)格模擬任意傾斜界面異常體的方法是有效的.

圖10 （a）TE視電阻率剖面（單位：Ωm）；（b）TE阻抗相位剖面（單位（°））；（c）TM視電阻率剖面（單位：Ωm）；（d）TM阻抗相位剖面（單位（°））Fig.10 TE profile of apparent resistivity（a）；TE profile of impedance phase（b）；TM profile of apparent resistivity（c）；TM profile of impedance phase（d）

6 結(jié) 論

本文提出三角單元網(wǎng)格剖分、電性參數(shù)分塊連續(xù)變化二維MT有限元數(shù)值模擬方法，比較矩形單元剖分、電性參數(shù)分塊均勻的模擬方法，這樣既不增加節(jié)點，節(jié)省了計算時間，又能很好模擬野外實際地形和任意地電體模型.通過對一維連續(xù)介質(zhì)模型模擬計算，與解析法計算結(jié)果對比，驗證了本文算法的可靠性；對前人地形模型模擬計算，結(jié)果與前人結(jié)果相符；對任意傾斜界面異常體模型進行模擬，結(jié)果表明本方法是有效的.同時指出，本方法仍然是考慮各向同性介質(zhì)情況，對各向異性介質(zhì)分塊連續(xù)變化的數(shù)值模擬方法將是本文的后續(xù)研究工作.

致 謝 在本文程序的編制、有限元網(wǎng)格剖分方法以及輔助場的計算中，陳小斌研究員和作者進行多次坦誠交流和討論，并提出許多寶貴意見，謹在此表示衷心感謝.

（References）

[1] Wanna maker P E，Stodt J A，Rijofi L.Two－dimensional topographic responses in magnetotellurics modeled using finite elements.Geophys.，1986，51（11）：2131－2144.

[2] Mauriello P，Patella D.Principles of probability tomography for natural－source electromagnetic induction fields.Geophys.，1999，64（5）：1403－1417.

[3] 陳小斌，張翔，胡文寶.有限元直接迭代算法在MT二維正演計算中的應(yīng)用.石油地球物理勘探，2000，35（4）：487－496.Chen X B，Zhang X，Hu W B.Application of finiteelementdirect iteration algorithm to MT 2－D forward computation.Oil Geophysical Prospecting （in Chinese），2000，35（4）：487－496.

[4] 陳小斌.MT二維正演計算中地形影響的研究.石油物探，2000，39（3）：112－120.Chen X B.On the research of the influence of terrain to MT 2D forward computation.Geophysical Prospecting for Petroleum （in Chinese），2000，39（3）：112－120.

[5] 謝飛.帶地形大地電磁場二維有限元數(shù)值模擬研究.北京：中國地質(zhì)大學(xué)，2006.Xie F. Two－dimension Topographic Numeric Simulation Research in Magnetotellurics Using Finite Elements（in Chinese）.Beijing：China University of Geosciences，2006.

[6] 周熙襄，鐘本善.電法勘探數(shù)值模擬技術(shù).成都：四川科學(xué)技術(shù)出版社，1986.Zhou X X，Zhong B S，et al.Numerical Modeling for Electrical Exploration （in Chinese）.Chengdu：Sichuan Science and Technology Press，1986.

[7] 王緒本，李永年，高永才.大地電磁測深二維地形影響及其校正方法研究.物探化探計算技術(shù)，1999，21（4）：328－332.Wang X B，Gao Y N，Gao Y C.Two dimentional topographic responses in magneto telluric sounding and its correction methods.Computing Techniques for Geophysical and Geochemical Exploration （in Chinese），1999，21（4）：328－332.

[8] 劉云，王緒本.大地電磁二維自適應(yīng)地形有限元正演模擬.地震地質(zhì)，2010，32（3）：382－391.Liu Y，Wang X B.FEM using adaptive topography in 2DMT forward modeling.Seismology and Geology （in Chinese），2010，32（3）：382－391.

[9] 徐世浙，于濤，李予國等.電導(dǎo)率分塊連續(xù)變化的二維MT有限元模擬（Ⅰ）.高校地質(zhì)學(xué)報，1995，1（2）：65－73.Xu S Z，Yu T，Li Y G，et al.The finite element method for modeling 2－D MT field on a geoelectrical model with continuous variation of conductivity within each block.Geological Journal of Universities （in Chinese），1995，1（2）：65－73.

[10] 李予國，徐世浙，劉斌等.電導(dǎo)率分塊連續(xù)變化的二維MT有限元模擬（Ⅱ）.高校地質(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 field on a geoelectrical model with continuous variation of conductivity within each block.Geological Journal of Universities （in Chinese），1996，2（4）：448－452.

[11] 樸化榮.電磁測深法原理.北京：地質(zhì)出版社，1990.Piao H R.Theory of Electro magnetics（in Chinese）.Beijing：Geology Press，1990.

[12] 曾國.大地電磁二維有限元正演數(shù)值模擬.湖南：中南大學(xué)，2008.Zeng G.FEM in 2－D Forward Modeling （in Chinese）.Hunan：Central South University，2008.

[13] 徐世浙.地球物理中的有限單元法.北京：科學(xué)出版社，1994.Xu S Z.FEM in Geophysics（in Chinese）.Beijing：Science Press，1994.

[14] 阮百堯.三角單元部分電導(dǎo)率分塊連續(xù)變化點源二維電場有限元數(shù)值模擬.廣西科學(xué)，2001，8（1）：1－3.Ruan B Y.2－D electrical modeling due to a current point by FEM with variation of conductivity within each triangular element.Guangxi Sciences（in Chinese），2001，8（1）：1－3.

[15] 阮百堯，熊 彬.大型對稱變帶寬方程組的Cholesky分解法.物探化探計算技術(shù)，2000，22（4）：361－368.Ruan B Y，Xiong B.A cholesky decomposition method for large－scale symmetrical system of equations with varying band－width.Computing Techniques for Geophysical and Geochemical Exploration （in Chinese），2000，22（4）：361－368.

[16] 劉德貴，費景高，于泳江等.FORTRAN算法匯編（第一分冊）.北京：國防工業(yè)出版社，1980.Liu D G，F(xiàn)ei J G，Yu Y J，et al.FORTRAN Algorithm Collection （Volume 1） （in Chinese）.Beijing：National Defense Industry Press，1980.

[17] 馬為.MT二維正演中輔助場計算新方法研究.北京：中國地震局，2007.Ma W.A Study of New Algorithm for the Calculation of Auxiliary Field in MT 2D Forward Modeling （in Chinese）.Beijing：China Earthquake，2007.

[18] 馬為，陳小斌，趙國澤.大地電磁測深二維正演中輔助場的新算法.地震地質(zhì)，2008，30（2）：525－533.Ma W，Chen X B，Zhao G Z.A new algorithm for the calculation of auxiliary field in MT 2D forward modeling.Seismology and Geology （in Chinese），2008，30（2）：525－533.

[19] 史明娟，徐世浙，劉斌.大地電磁二次函數(shù)插值的有限元法正演模擬.地球物理學(xué)報，1997，40（3）：421－430.Shi M J，Xu S Z，Liu B.Finite element method using quadratic element in MT forward modeling.Chinese J.Geophys.（in Chinese），1997，40（3）：421－430.

[20] 徐世浙，王慶乙，王軍.用邊界單元法模擬二維地形對大地電磁場的影響.地球物理學(xué)報，1992，35（3）：380－388.Xu S Z，Wang Q Y，Wang J.Modeling 2－D terrain effect on MT by the boundary element method.Chinese J.Geophys.（in Chinese），1992，35（3）：380－388.

The FEM for modeling 2－D MT with continuous variation of electric parameters within each block

LIU Yun，WANG Xu－Ben＊
State Key Laboratory of Oil and Gas Reservoir Geology and Exploitation，Chengdu University of Technology，Chengdu610059，China

To model arbitrarily shaped two－dimensional topography and structures in field work，triangular element grid was used in the finite element method （FEM）.In view of the fact of continuous variation of the subterranean rock－mineral electric parameters，the electromagnetic field and some electric parameters of models are designed to bilinear variation within each triangular element in our numeric modeling method，and which is developed for modeling two dimensional magnetotelluric （MT）under the field topography condition.The calculation for mulae of the auxiliary field，and the definition of apparent resistivity and impedance phase are deduced according to the relationship between the main fields of the three nodes and the linear shape function within each element.By calculating a continuous medium model and two topography models set up by other scholars to test our method，the result of our method shows a high accuracy （the mean square error is less than 1%），and the results of modeling two topography models accord with other scholar′s，too.Through modeling a sloping interface abnormity body，we find that our method can model arbitrarily complicated terrain and geoelectric bodies preferably.

MT，Topography，Triangle grids，Continuous variation of electric parameters within each block，F(xiàn)EM

10.6038/j.issn.0001－5733.2012.06.029

P631

2011－06－17，2012－05－08收修定稿

國家高技術(shù)研究發(fā)展計劃重點項目子課題（2009AA06Z108）資助.

劉云，男，1973年生，博士研究生，主要從事電磁場數(shù)值模擬及反演成像技術(shù)研究.E－mail：liuyun0077＠163.com

＊通訊作者 王緒本，男，教授、博士生導(dǎo)師，主要從事電磁場數(shù)值模擬及反演成像技術(shù)研究.E－mail：wxb5501＠yahoo.com.cn

劉云，王緒本.電性參數(shù)分塊連續(xù)變化二維 MT有限元數(shù)值模擬.地球物理學(xué)報，2012，55（6）：2079－2086，

10.6038/j.issn.0001－5733.2012.06.029.

Liu Y，Wang X B.The FEM for modeling 2－D MT with continuous variation of electric parameters within each block.Chinese J.Geophys.（in Chinese），2012，55（6）：2079－2086，doi：10.6038/j.issn.0001－5733.2012.06.029.

（本文編輯 汪海英）