電導率連續(xù)變化大地電磁無網(wǎng)格法數(shù)值模擬*

黃廷哲，嵇艷鞠，2?，黃婉玉，關珊珊，姜曜

(1.吉林大學 儀器科學與電氣工程學院，吉林 長春 130026；2.地球信息探測儀器教育部重點實驗室，吉林 長春 130026)

電導率連續(xù)變化大地電磁無網(wǎng)格法數(shù)值模擬*

黃廷哲1，嵇艷鞠1，2?，黃婉玉1，關珊珊1，姜曜1

(1.吉林大學 儀器科學與電氣工程學院，吉林 長春 130026；2.地球信息探測儀器教育部重點實驗室，吉林 長春 130026)

針對地下空間目標體多為電性參數(shù)連續(xù)變化且分布形態(tài)復雜等特點，結(jié)合無網(wǎng)格法物性參數(shù)加載方便、計算精度較高、自適應分析便利等優(yōu)勢，基于徑向基函數(shù)構(gòu)造了連續(xù)導電模型下的形函數(shù)，推導了等價線性方程組，進行了電導率連續(xù)變化的二維大地電磁數(shù)值模擬，給出了形狀參數(shù)最優(yōu)值.通過電導率連續(xù)變化的水平層狀模型驗證了算法的正確性，計算結(jié)果均方根相對誤差不超過0.36%，精度優(yōu)于有限元法.討論了電導率連續(xù)變化的水平模型與均勻分塊模型的電磁響應差異，對連續(xù)變化的地塹模型和不同傾斜角度下油藏注水模型進行了數(shù)值模擬.研究表明：連續(xù)變化模型和均勻分塊模型差異明顯，在實際反演解釋中采用連續(xù)變化模型有利于提高反演的精度；TM模式的觀測方式對于異常體的傾斜分布具有更好的分辨能力；無網(wǎng)格法避免了復雜的模型輸入和網(wǎng)格的生成，更適合計算電性參數(shù)連續(xù)變化和復雜分布的異常體響應，將成為復雜電性和分布形態(tài)下電磁探測高精度數(shù)值模擬新方法.

連續(xù)變化；大地電磁；無網(wǎng)格法；徑向基點插值

大地電磁測深法(Magnetotelluric Sounding，MT)是通過在地面觀測相互正交的電場和磁場分量來研究地球內(nèi)部電性結(jié)構(gòu)的一種電磁探測方法，在實際探測中得到了廣泛的應用[1-2].在電性參數(shù)分塊均勻的前提下， Mauriello[3]采用矩形三角網(wǎng)格剖分、雙線性插值有限元進行了數(shù)值模擬;Kerry等[4]采用非結(jié)構(gòu)化網(wǎng)格實現(xiàn)了任意復雜地電模型的大地電磁自適應有限元模擬.電性參數(shù)分塊均勻的模擬方式認為在單元塊內(nèi)電性參數(shù)是均勻分布的，但是在實際的深部地球探測中，由于地溫梯度變化以及壓實作用的存在，巖礦石電性參數(shù)往往是連續(xù)變化的，例如在煤田和油氣探測中，含煤巖層和飽和油藏的巖石物理學參數(shù)就表現(xiàn)出了連續(xù)變化的性質(zhì)[5]，分塊均勻模型只是連續(xù)變化模型的特殊表現(xiàn)形式，采用電性參數(shù)連續(xù)變化的剖分方式更符合電磁勘探的實際情況.另外，某些初步反演方法得到的電性參數(shù)也是連續(xù)變化的，連續(xù)導電模型的正演可為此類模型的反演打下基礎，因此有必要對電導率連續(xù)變化的模型進行數(shù)值模擬.劉云等[6]對連續(xù)變化的模型進行了正演計算，但都是采用將電性參數(shù)賦予單元節(jié)點上的有限元方法，在子單元內(nèi)電性參數(shù)采用插值構(gòu)造，并未真正反映真實的電導率連續(xù)變化的情況.

無網(wǎng)格方法是伴隨著工程計算領域的深入發(fā)展應運而生的一種新興數(shù)值計算方法，其近似函數(shù)僅建立在一系列離散節(jié)點上，克服了傳統(tǒng)方法對于網(wǎng)格的依賴，解決了動態(tài)斷裂力學、金屬體積成形等常規(guī)方法難以模擬的問題[7-8].在計算電磁學領域,Cingoski等[9]首次將無單元Galerkin法(Element-Free Galerkin method,EFGM)應用于電磁學計算中，成功解決了電壓互感器問題.Ho等[10]總結(jié)了無網(wǎng)格法在計算電磁學中取得的初步成果.馮德山等[11]研究了無單元Galerkin法的二維探地雷達正演模擬，結(jié)果表明在相同的節(jié)點條件下EFGM比矩形剖分的有限單元法精度要高.李俊杰等[12]研究了大地電磁中無網(wǎng)格與有限元的耦合算法.目前無網(wǎng)格法主要用于電導率分塊均勻模型的電磁數(shù)值模擬[13]，而且形函數(shù)大多數(shù)采用移動最小二乘法構(gòu)造，導致本質(zhì)邊界條件的加載較為困難，對于實際勘探中廣泛存在的電性參數(shù)連續(xù)變化的情況還未見到有相關成果發(fā)表.

本文提出用無網(wǎng)格徑向基點插值法(Radial Point Interpolation Method,RPIM)模擬電導率連續(xù)變化的大地電磁模型，推導了對應等價線性方程組，將物性參數(shù)的加載在一個個離散的點上，較好地模擬了復雜異常體的形狀和電性參數(shù).基于Matlab軟件平臺編制了相應程序并計算了典型的電導率連續(xù)變化模型，結(jié)果表明無網(wǎng)格法能夠有效地計算復雜模型下的電磁場響應，視電阻率和相位的斷面圖均較好地反映了電性參數(shù)連續(xù)變化的趨勢，同時無網(wǎng)格法具有高精度計算和自適應分析便利的特點，有望成為復雜媒質(zhì)和復雜構(gòu)造下的電磁法數(shù)值模擬一種新的強有力方法.

1 無網(wǎng)格法正演理論

1.1 連續(xù)導電媒質(zhì)形函數(shù)的構(gòu)造

多項式基函數(shù)插值法被廣泛運用于各類數(shù)值方法中，為了避免傳統(tǒng)的多項式基函數(shù)插值所引起的奇異性，針對地下空間的連續(xù)導電媒質(zhì)問題，采用徑向基函數(shù)(RBF)構(gòu)造無網(wǎng)格法形函數(shù).由于無網(wǎng)格法不依賴于預定義的單元構(gòu)造形函數(shù)，因此引入了支持域的概念，支持域一般為矩形或者圓形，以計算點為中心，如圖1所示.

圖1 無網(wǎng)格法節(jié)點離散示意圖Fig.1 Nodes discrete sketch of meshfree method

在求解域Ω內(nèi)，任意一節(jié)點XT=(x,y)的場變量u(X)在局部支持域Ωs的徑向基點插值表達式為[14]：

Us=R0a+Pmb

(1)

節(jié)點場值向量Us為:

Us={u1u2…un}T

(2)

RBF的力矩矩陣為:

(3)

多項式力矩矩陣為:

(4)

式中:n為RBF的個數(shù);m為多項式基函數(shù)的個數(shù)；ai和bj為待定常數(shù).Ri(rk)中的rk的表達式為:

(5)

本文選用的RBF為復合2次(MQ)函數(shù)，多項式基函數(shù)的個數(shù)為3，MQ函數(shù)其表達式為：

(6)

其中αc和q為形狀參數(shù)；ri是計算點與節(jié)點之間的距離.

然而式(1)中有(n+m)個變量，添加m個約束條件：

(7)

聯(lián)立式(1)和式(7)可得到如下矩陣方程：

(8)

求解式(8)可得到：

(9)

可將式(9)重寫為：

u(x) =RT(x)a+PT(x)b=

(10)

式中的RPIM形函數(shù)可表示為：

{φ1(x)φ2(x)…φn(x)φn+1(x)…φn+m(x)}

(11)

1.2 MT邊值問題無網(wǎng)格法離散

(12)

圖2 二維大地電磁模擬示意圖Fig.2 Sketch map of 2D MT simulation

大地電磁兩種極化模式分別如下：

(13)

對于TE模式，計算區(qū)域包含空氣層；對于TM模式，研究區(qū)域為地下空間，上邊界AB在地面上.求取泛函F(u)對場量u的偏導數(shù)，將u(x)=ΦT(x)U帶入，令泛函F(u)對場量u的偏導數(shù)為0，從而得到如下矩陣表達式：

Ku=b

(14)

式中：u={u1u2… un}T；b=[0 0 … 0]；

K=K1-K2+K3=

∫ΩλΦiΦjdΩ+∫CDτkΦiΦjdΩ

為了方便數(shù)值積分，將求解域Ω劃分為若干積分單元，在求解域內(nèi)和邊界積分均采用高斯積分實現(xiàn).上邊界AB是本質(zhì)邊界，其場值始終是1個單位，為了加載本質(zhì)邊界，將系數(shù)矩陣K的對角線元素修正為αKii，對應的右端項b修正為αKii.其中α是一較大的數(shù)，這樣使系數(shù)矩陣只在兩處發(fā)生了變化，本質(zhì)邊界條件處理較為簡單.

Krylov子空間方法被認為是一種有效的求解大型稀疏線性方程組方法[16].基于Krylov子空間的迭代方法收斂速度快，求解精度高，而且穩(wěn)定性好.Krylov子空間方法常常與預處理技術結(jié)合求解線性方程組，本文采用有容許誤差的不完全LU分解作為預條件子和穩(wěn)定雙共軛梯度法(BICGSTAB)進行求解，容許誤差為1×10-10，最大迭代次數(shù)為1 000次，對方程組(14)求解即可得到各個節(jié)點處的場值.

(15)

式中：h是垂直方向上等間距兩個節(jié)點的距離.

得到輔助場值之后，可以通過式(16)計算視電阻率和相位：

(16)

2 模型計算及分析

2.1 RPIM最優(yōu)形狀參數(shù)選取

徑向基點插值無網(wǎng)格法采用RBF構(gòu)造形函數(shù)，RBF中的形狀參數(shù)對計算精度有直接影響，其最優(yōu)值需要通過數(shù)值試驗獲得[17].建立ρ=100 Ω·m的均勻半空間模型，表1分別為不同支持域無量綱尺寸α，不同高斯點個數(shù)n，不同q值，以及不同αc下的視電阻率計算結(jié)果.選取計算頻率為1 Hz，計算機硬件平臺為CPU Intel core i7-4770 3.40 GHz，RAM 8 GB，軟件編程平臺為MATLAB 2015a.

由表1可以看出，支持域無量綱尺寸對于計算結(jié)果有較大影響，α的最優(yōu)值介于1.1～1.6之間.高斯點個數(shù)n越大，計算結(jié)果越精確，但耗時增加較多，4點高斯即可滿足精度要求.MQ函數(shù)中的形狀參數(shù)q和αc具有較好的普適性，αc的取值范圍約在[-45,45]之間，q的取值范圍約為[-2.3,0.9]，在此范圍內(nèi)計算結(jié)果均保持較高的計算精度，超出此范圍形函數(shù)的構(gòu)造將出現(xiàn)奇異性.綜合以上分析，本文選取的形狀參數(shù)如下：支持域無量綱尺寸α為1.1，高斯點數(shù)量為4個，q=0.5，αc=1.

表1 不同形狀參數(shù)下的視電阻率計算結(jié)果Tab.1 Apparent resistivity calculation results under different shape parameters

2.2 連續(xù)變化算法驗證

建立如圖3所示的水平層狀模型,求解域橫向規(guī)模為20 km，縱向規(guī)模為30 km.第1層及第3層電阻率是均勻的，第1層電阻率為100 Ω·m，厚度為2 km，底層電阻率為200 Ω·m.中間層含有一個電導率連續(xù)變化的水平層，其電阻率值隨著深度線性增大，深度為4 km.節(jié)點均勻分布于求解域內(nèi)，一共201×301個節(jié)點，節(jié)點間距100 m.背景網(wǎng)格的頂點和節(jié)點重合.

圖3 電導率連續(xù)變化層狀模型Fig.3 Continuous conductivity layered model

為了驗證算法的正確性，將圖3所示的連續(xù)導電媒質(zhì)層分為19層按照層狀模型求出解析解，在相同的節(jié)點分布下取地面中心點的測點為研究對象，將無網(wǎng)格法數(shù)值解與解析解、有限元法數(shù)值解進行對比,視電阻率和阻抗相位結(jié)果如表2和表3所示.

從表中可以看出，無網(wǎng)格法的數(shù)值解與解析解、有限元數(shù)值解結(jié)果高度一致，其中無網(wǎng)格法視電阻率和阻抗相位與解析解的均方根相對誤差不超過0.36%，而有限元法均方根相對誤差不超過0.74%，體現(xiàn)了無網(wǎng)格法高精度計算的特點.由于此處無網(wǎng)格法采用了均勻節(jié)點分布，在高頻時趨膚深度較淺，因此導致在高頻時計算精度略低于有限元法，計算時間約為有限元法的10倍左右，但無網(wǎng)格法的物性參數(shù)加載在離散的點上，避免了復雜的模型輸入文件和網(wǎng)格的生成，計算結(jié)果的均方根相對誤差更小，表明無網(wǎng)格法可以有效地進行較高精度下的連續(xù)導電媒質(zhì)模型下大地電磁響應計算.

2.3 連續(xù)變化模型和均勻分塊模型對比

為了說明連續(xù)變化模型分析的必要性，將圖3所示的連續(xù)變化層的電導率取平均，并和原來的模型進行對比分析.第1層和第3層電阻率與圖3相同，中間層電阻率取為圖3電阻率的平均值，厚度仍為4 km，其余參數(shù)均不變.連續(xù)導電媒質(zhì)模型與一維層狀模型的地面中心點處的視電阻率和相位對比結(jié)果如圖4所示.從圖中可以看出，兩者形態(tài)基本一致，但在一些頻點上差異較大，其中視電阻率的最大差異可達9.8%，相位的最大差異可達7.6%，這表明連續(xù)變化模型和均勻分塊模型的電磁響應差異較為明顯，如果實際地電模型的電導率為連續(xù)變化，而依然采用傳統(tǒng)的分塊均勻模型進行解釋，將會造成較大的偏差.

表2 視電阻率數(shù)值解和解析解對比Tab.2 Comparison of apparent resistivity between numerical solution and analytical solution

表3 相位數(shù)值解和解析解對比Tab.3 Comparison of phase between numerical solution and analytical solution

(a)視電阻率

(b)相位圖4 連續(xù)變化模型與分塊均勻模型結(jié)果對比Fig.4 Comparison between continuous medium model and homogeneous model

2.4 連續(xù)變化地塹模型計算

分別建立如圖5所示的地塹連續(xù)變化模型，第1層電阻率10 Ω·m，厚度為2 km，中間1層為過渡層，其電阻率隨深度線性增大，第2層電阻率的變化為:

(17)

式中：ρ1為第1層電阻率；ρ3為第3層電阻率；H為當前的層厚度；h1,h2分別為第1層和第2層的層厚.在地塹存在的區(qū)域，h2為4 km，其他區(qū)域h2為2 km.第3層為均勻媒質(zhì)層，電阻率為100 Ω·m.

圖5 地塹連續(xù)變化模型Fig.5 Continuous conductivity graben model

(a)視電阻率

(b)相位圖6 模型5(a)地面中心點測深曲線Fig.6 Ground center point sounding curve of model 5(a)

圖5(b)模型的地下連續(xù)導電媒質(zhì)交界面為起伏變化，用傳統(tǒng)的方法已較難模擬.自適應有限元方法雖然可以對該模型進行模擬，但該方法需要經(jīng)過多次迭代及誤差估計，提高局部精度的同時卻降低了整體的收斂速度，而無網(wǎng)格法脫離了網(wǎng)格的限制，只需離散的點便可構(gòu)造形函數(shù)，可以非常方便地對該模型進行模擬.圖7為圖5(b)模型的無網(wǎng)格法視電阻率計算結(jié)果.

(a)TE模式

(b)TM模式圖7 模型5(b)視電阻率模擬結(jié)果Fig.7 Resistivity simulation results of model 5(b)

從模擬結(jié)果中可以看出，起伏地形在視電阻率擬斷面圖較為明顯，較為清晰地反映出了所設計模型電阻率連續(xù)變化的趨勢.由于切向電場的連續(xù)性，TM模式更加清晰地反映了媒質(zhì)交界處的起伏地形情況.模擬結(jié)果符合TE模式垂向分辨率較高，而TM模式橫向分辨率較高的特點.

2.5 實際油藏注水模型計算

在石油地球物理勘探中，由于受到含油飽和度、含油高度、孔隙結(jié)構(gòu)以及油水密度差等因素影響，油藏的電導率往往表現(xiàn)出連續(xù)變化的趨勢.為了模擬實際中的電導率連續(xù)變化模型以及進一步突出無網(wǎng)格法計算復雜模型便利的優(yōu)勢，建立如圖8所示的飽和油藏區(qū)注水模型，注水區(qū)呈橢圓形分布，電導率分布函數(shù)為：

σ(x,z)=σ0+Aexp(-(0.000 1x2+

0.000 01(z-1 500)2))

(18)

其中背景電導率為0.01 S·m，A=0.99.為了研究不同傾斜角度下的大地電磁響應，分別對注水區(qū)橢圓異常旋轉(zhuǎn)0°，30°，60°和90°，不同旋轉(zhuǎn)角度下的電阻率分布情況如圖8所示.

(a)旋轉(zhuǎn)0°

(b) 旋轉(zhuǎn)30°

(c)旋轉(zhuǎn)60°

(d) 旋轉(zhuǎn)90°圖8 飽和油藏注水區(qū)模型Fig.8 Saturated reservoir water flooding model

該模型油藏注水區(qū)呈橢圓形，并且電導率呈非線性分布，采用傳統(tǒng)的有限元方法已較難模擬.但無網(wǎng)格法不需要利用預定義的網(wǎng)格信息，可以方便地對該模型進行異常體形狀的模擬及電性參數(shù)的加載.采用和圖5相同的無網(wǎng)格參數(shù)，縱向采用非均勻節(jié)點分布，為了突出顯示連續(xù)導電媒質(zhì)模型的響應，分別計算相同旋轉(zhuǎn)角度下的分塊均勻模型的視電阻率，并與連續(xù)導電媒質(zhì)模型下的響應相除，圖9為在1 Hz頻率下的視電阻率剖面曲線.

(a)TE模式

(b)TM模式圖9 不同傾斜角度下的視電阻率剖面曲線Fig.9 Apparent resistivity curves under different inclination angle

從圖9中可以看出，TE模式對異常體的傾向不敏感，在不同傾斜角度下的視電阻率剖面曲線幾乎重合，而TM模式對異常體的傾向較為敏感，在旋轉(zhuǎn)角度為0°和90°時呈現(xiàn)了較好的對稱性，對不同的傾斜角度具備了一定的分辨能力.這是由于在TM模式下，磁場分量平行于走向，存在著與構(gòu)造走向垂直的電場分量，非均勻媒質(zhì)接觸面在電場的作用下充電而產(chǎn)生了多余的附加電荷，導致TM極化模式受異常體的影響較為明顯.在TE極化模式中，電場分量僅僅沿走向方向，不存在與走向垂直的電場，在不均勻媒質(zhì)交界面不產(chǎn)生電荷積累，因此TE模式對不同傾斜角度下的不均勻體幾乎表現(xiàn)不出分辨能力.

3 結(jié) 論

本文針對野外探測中廣泛存在的電性參數(shù)連續(xù)變化的實際情況以及實際反演的需要，提出將無網(wǎng)格法應用到電磁法數(shù)值模擬中，基于徑向基函數(shù)構(gòu)造了連續(xù)導電模型下的形函數(shù)，推導了等價線性方程組，進行了電導率連續(xù)變化的二維大地電磁數(shù)值模擬，通過數(shù)值試驗給出了RPIM形狀參數(shù)最優(yōu)值.主要結(jié)論如下：

1)電導率連續(xù)變化水平模型的無網(wǎng)格法計算結(jié)果均方根相對誤差不超過0.36%，算法精度優(yōu)于有限元法.連續(xù)變化模型和均勻分塊模型差異明顯，其中視電阻率的最大差異可達9.8%，在實際反演解釋中采用連續(xù)導電媒質(zhì)模型有利于提高反演的精度.TM極化模式對于異常體的傾斜狀況更加敏感.

2)基于徑向基函數(shù)的無網(wǎng)格法形函數(shù)構(gòu)造脫離了網(wǎng)格的限制，本質(zhì)邊界條件加載非常容易，計算精度相對較高，物性參數(shù)的加載在一個個離散的計算點上，可通過形狀參數(shù)來調(diào)整支持域的大小，以實現(xiàn)在模型的特征區(qū)精細化模擬，特別適合連續(xù)變化等復雜形態(tài)地質(zhì)構(gòu)造的計算，自適應分析便利，具有較強的實用性，將成為復雜電性和分布形態(tài)下電磁探測高精度數(shù)值模擬新方法.

3)相比于分塊均勻的正演模擬，許多有效的反演方法(如 Bostick，Occam 等 )都是假設地下導電媒質(zhì)是連續(xù)的，而無網(wǎng)格法為連續(xù)模型正演計算提供了新的思路，在電磁探測領域?qū)⒕邆鋸V闊的應用前景.

致謝：吉林大學李桐林教授，陳漢波博士及中國海洋大學盧杰博士、潘林冬碩士對本文的研究和寫作提供了指導和幫助，審稿專家對本文提出了寶貴意見和建議，在此表示衷心的感謝.

[1] 顧觀文,吳文鸝,李桐林.大地電磁場三維地形影響的矢量有限元數(shù)值模擬[J].吉林大學學報:地球科學版,2014,44(5):1678-1686.

GU Guanwen,WU Wenli,LI Tonglin.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)

[2] 湯井田，李晉，肖曉,等.基于數(shù)學形態(tài)濾波的大地電磁強干擾分離方法[J].中南大學學報:自然科學版，2012,43(6):2215-2221.

TANG Jintian,LI Jin,XIAO Xiao,etal.Magnetotelluric sounding data strong interference separation method based on mathematical morphology filtering[J].Journal of Central South University:Science and Technology,2012,43(6):2215-2221.(In Chinese)

[3] MAURIELLO P,PATELLA D.Principles of probability tomography for natural-source electromagnetic induction fields[J].Geophysics，1999,64(5):1403-1417.

[4] KERRY K,CHESTER W.Adaptive finite-element modeling using unstructuredgrids:The 2D maggnetotelluric example[J].Geophysics，2006,71(6):291-299.

[5] ГРЕЧУХИН В В著，蔡柏林等譯.應用地球物理方法研究含煤建造[M].北京:地質(zhì)出版社,1987:51-60.

ГРЕЧУХИН В В.Translated by CAI Bailin，etal.Application of geophysical methods studying coal bearing building[M].Beijing:Geology Press，1987:51-60.(In Chinese)

[6] 劉云，王緒本.電性參數(shù)分塊連續(xù)變化二維MT有限元數(shù)值模擬[J].地球物理學報，2012,55(6):2079-2086.

LIU Yun,WANG Xuben.The FEM for modeling 2-D MT with continuous variation of electric parameters within each block [J].Chinese Journal of Geophysics，2012，55(6):2079-2086.(In Chinese)

[7] 龍述堯,張國虎.基于MLPG法的動態(tài)斷裂力學問題[J].湖南大學學報:自然科學版,2012,39(11):41-45.

LONG Shuyao,ZHANG Guohu.An analysis of the dynamic fracture problem by the meshless local Petrov-Galerkin method[J].Journal of Hunan University:Natural Sciences,2012,39(11):41-45.(In Chinese)

[8] 鄭剛，伍素珍，李光耀，等.金屬體積成形過程的無網(wǎng)格RPIM方法分析[J].湖南大學學報:自然科學版,2010,37(10):41-46.

ZHENG Gang,WU Suzhen ,LI Guangyao,etal.Numerical simulation of bulk forming processes by radial point interpolation method(RPIM) [J].Journal of Hunan University:Natural Sciences，2010，37(10):41-46.(In Chinese)

[9] CINGOSKI V，MIYAMOTO N，YAMASHIT A.Hybrid element-free Galerkin-finite element method for electromagnetic field computations[J].IEEE Transactions Magnetics,2000,36(4):1543-1547.

[10]HO S L，YANG S，MACHADO J M,etal.Application of a meshless method in electromagnetics[J].IEEE Transactions on Magnetics，2001,37(5):3198-3202.

[11]馮德山,王洪華,戴前偉.基于無單元Galerkin法探地雷達正演模擬[J].地球物理學報,2013，56(1):298-308.

FENG Deshan,WANG Honghua,DAI Qianwei.Forward simulation of ground penetrating radar based on the element-free Galerkin method [J].Chinese Journal of Geophysics,2013,56(1):298-308.(In Chinese)

[12]李俊杰,嚴家斌.大地電磁二維正演中的有限元-徑向基點插值法[J].中國有色金屬學報,2015,25(5):1314-1324.

LI Junjie,YAN Jiabing.Magnetotelluric two-dimensional forward by finite element-radial point interpolation method[J].The Chinese Journal of Nonferrous Metals,2015，25(5):1314-1324.(In Chinese)

[13]嵇艷鞠,黃廷哲,黃婉玉,等.起伏地形下各向異性的2D大地電磁無網(wǎng)格法數(shù)值模擬[J].地球物理學報,2016,59(12):4483-4493.

JI Yanju,HUANG Tingzhe,HUANG Wanyu,etal.2D anisotropic magnetotelluric numerical simulation using meshfree method under undulating terrain[J].Chinese Journal of Geophysics,2016,59(12):4483-4493.(In Chinese)

[14]LIU G R,GU Y T著.王建明，周學軍譯.無網(wǎng)格法理論及程序設計[M].濟南：山東大學出版社，2007:45-48.

LIU G R,GU Y T.Translated by WANG Jianming,ZHOU Xuejun.An introduction to meshfree methods and their programming [M].Jinan:Shandong University Press,2007:45-48.(In Chinese)

[15]徐世浙.地球物理中的有限單元法[M].北京:科學出版社，1994:229-239.

XU Shizhe.Finite element method in geophysics[M].Beijing:Science Press,1994:229-239.(In Chinese)

[16]柳建新，蔣鵬飛，童孝忠,等.不完全LU分解預處理的BICGSTAB算法在大地電磁二維正演模擬中的應用[J].中南大學學報:自然科學版,2009,40(2):484-491.

LIU Jianxin,JIANG Pengfei,TONG Xiaozhong,etal.Application of BICGSTAB algorithm with incomplete LU decomposition preconditioning in two-dimensional magnetotelluric forward modeling[J].Journal of Central South University:Science and Technology,2009,40(2):484-491.(In Chinese)

[17]李俊杰，嚴家斌.RPIM求解點源二維變分問題的最優(yōu)形狀參數(shù)[J].物探與化探,2015,39(6):1233-1237.

LI Junjie,YAN Jiabing.Optimal shape parameters of RPIM for resolving point source two-dimensional variational problem [J].Geophysical and Geochemical Exploration,2015,39(6):1233-1237.(In Chinese)

Magnetotellurics Simulation with Continuous Variation of Conductivity Based on Meshfree Method

HUANG Tingzhe1，JI Yanju1 2?，HUANG Wanyu1，GUAN Shanshan1，JIANG Yao1

(1.College of Instrumentation and Electrical Engineering，Jilin University，Changchun 130026，China；2.Key Laboratory of Geo-Exploration Instrumentation，Ministry of Education，Jilin University，Changchun 130026，China)

In view of the problem in which most target bodies are continuously changing and the distribution is complex in the actual geological structures,2-D magnetotelluric forward modeling is simulated combined with the advantage of meshfree method where parameter loading is convenient,calculation accuracy is high and adaptive analysis is convenient.The shape function of the continuous medium model is constructed by the radial basis function,the equivalent linear equations are derived and the corresponding program is compiled.The optimal values of the shape parameters are given by numerical experiments.The correctness of meshfree method is verified by a horizontal continuous variation layered model.The root mean square error of the calculation results is not more than 0.36%,and its accuracy is superior to the finite element method.The electromagnetic response of the horizontal layer with the continuous variation of electrical conductivity and the uniform block model is discussed,and a graben terrain model and a reservoir water flooding model whose conductivity is continuously variable are calculated,respectively.The conclusions are given: the electromagnetic responses between the continuous variation model and uniform block model are significantly different,the use of the continuous variation model is helpful to improve the accuracy of inversion in the actual interpretation,and the observation method of the TM mode has better resolution for the inclination distribution of anomalous body.The meshfree method does not need the complicated mesh generation,and it is more suitable to calculate the response of continuous variation and complex distribution models.Therefore,meshfree method is expected to be a new and robust method for the numerical simulation of complex medium and complex structures.

continuous variation；magnetotelluric；meshfree method；radial point interpolation

1674-2474(2017)08-0091-09

10.16339/j.cnki.hdxbzkb.2017.08.014

2016-12-11

國家自然科學基金資助項目(41674109),National Natural Science Foundation of China(41674109);國家重大科研裝備研制資助項目(ZDYZ2012-1-05),National Major Scientific Research and Equipment Development Project(ZDYZ2012-1-05)

黃廷哲(1988-)男，河南南陽人，吉林大學博士研究生

?通訊聯(lián)系人：E-mail:jiyj@jlu.edu.cn

TM153

A