李俊杰,嚴(yán)家斌

(1.浙江省水利水電勘測(cè)設(shè)計(jì)院,浙江杭州310002;2.中南大學(xué)地球科學(xué)與信息物理學(xué)院,有色資源與地質(zhì)災(zāi)害探查湖南省重點(diǎn)實(shí)驗(yàn)室,湖南長(zhǎng)沙410083)

有限元-點(diǎn)插值耦合法大地電磁二維正演模擬

李俊杰1,嚴(yán)家斌2

(1.浙江省水利水電勘測(cè)設(shè)計(jì)院,浙江杭州310002;2.中南大學(xué)地球科學(xué)與信息物理學(xué)院,有色資源與地質(zhì)災(zāi)害探查湖南省重點(diǎn)實(shí)驗(yàn)室,湖南長(zhǎng)沙410083)

點(diǎn)插值法(PIM)作為一種典型的全域弱式無(wú)網(wǎng)格法,該方法在地質(zhì)建模時(shí)將物性加載到只與坐標(biāo)有關(guān)的高斯積分點(diǎn)上,因此處理復(fù)雜模型時(shí)較常規(guī)網(wǎng)格方法便利,但缺點(diǎn)是計(jì)算效率低。將有限元法(FEM)與PIM耦合,形成FE-PIM,用于大地電磁二維正演模擬。利用Galerkin法代入插值法構(gòu)造的形函數(shù)并結(jié)合高斯積分公式推導(dǎo)了大地電磁二維無(wú)網(wǎng)格化總體矩陣表達(dá)式,簡(jiǎn)述了背景網(wǎng)格積分與邊界條件的加載技術(shù),理論模型的數(shù)值計(jì)算驗(yàn)證了FE-PIM算法的正確性、高效性及其在處理復(fù)雜模型上的便利性。

點(diǎn)插值法;全域弱式無(wú)網(wǎng)格法;大地電磁;有限元-點(diǎn)插值法

大地電磁二維正演多采用有限差分法(finite deference method,FDM)[1-2]、有限元法(finite element method,FEM)[3-6]等網(wǎng)格方法。對(duì)于地質(zhì)體相對(duì)規(guī)則、地形簡(jiǎn)單的模型,FDM利用較小的網(wǎng)格就能獲得較高的計(jì)算精度,效率高;而對(duì)于復(fù)雜的地質(zhì)體,就需要精細(xì)網(wǎng)格才能完成高精度的模擬。FEM適用于復(fù)雜邊界與復(fù)雜電性結(jié)構(gòu)模型的計(jì)算,但高維問(wèn)題中非結(jié)構(gòu)化三角單元生成程序?qū)崿F(xiàn)困難。無(wú)單元Galerkin法(element-free Galerkin method,EFGM)[7]是一種基于節(jié)點(diǎn)的全域弱式無(wú)網(wǎng)格方法,其形函數(shù)構(gòu)造不依賴(lài)預(yù)定義的網(wǎng)格單元,選擇離散的節(jié)點(diǎn)便可求Helmholtz方程的解。相對(duì)于低階(low-order)有限元分析,EFGM只需少量的節(jié)點(diǎn)就可獲得光滑的解,且場(chǎng)分量的導(dǎo)數(shù)同樣光滑連續(xù),對(duì)于大地電磁輔助分量的求解非常有利,可以獲得與主分量(水平電場(chǎng)或磁場(chǎng))相同的計(jì)算精度,對(duì)于節(jié)點(diǎn)的參數(shù)化處理也較單元分析容易。EFGM能處理網(wǎng)格方法難以處理的問(wèn)題[8-10],有廣闊的應(yīng)用前景。李俊杰等[11]綜述了無(wú)網(wǎng)格法進(jìn)展并推導(dǎo)了大地電磁二維問(wèn)題的無(wú)網(wǎng)格化系統(tǒng)矩陣表達(dá)式;賈曉峰等[12-13]將EFGM用于疊前地震模擬,并討論了幾種吸收邊界方法在地震模擬中與EFGM的結(jié)合,研究結(jié)果表明，EFGM精度較FDM精度高，且具有較好的穩(wěn)定性;王月英[14]分析了基函數(shù)階次對(duì)EFGM用于地震波正演計(jì)算精度的影響;馮德山等[15]、戴前偉等[16]將EFGM用于雷達(dá)波場(chǎng)的二維正演模擬,詳細(xì)推導(dǎo)了探地雷達(dá)無(wú)單元法波動(dòng)方程,并提出了適用于EFGM雷達(dá)正演模擬的透射吸收邊界、Sarma吸收邊界的改進(jìn)方法。嚴(yán)家斌等[17]將EFGM成功應(yīng)用于大地電磁二維正演模擬,并分析了支持域無(wú)量綱尺寸與高斯點(diǎn)數(shù)量對(duì)無(wú)網(wǎng)格法正演計(jì)算精度與效率的影響。無(wú)網(wǎng)格點(diǎn)插值法(point interpolation method,PIM)[18]是另一種簡(jiǎn)單高效的全域弱式無(wú)網(wǎng)格法,與EFGM的區(qū)別主要在于形函數(shù)采用過(guò)點(diǎn)插值的方法進(jìn)行構(gòu)造,邊界條件加載便利。李俊杰等[19]將PIM用于復(fù)雜模型的MT二維正演模擬,凸顯了PIM較常規(guī)網(wǎng)格方法在地質(zhì)建模上的優(yōu)越性。

我們將PIM與高效的FEM相耦合,形成有限元-點(diǎn)插值法(finite element-point interpolation method,FE-PIM),從大地電磁場(chǎng)的二維正演理論出發(fā),采用Galerkin法推導(dǎo)了FE-PIM系統(tǒng)矩陣表達(dá)式,數(shù)值計(jì)算結(jié)果驗(yàn)證了FE-PIM的有效性及其在計(jì)算復(fù)雜模型時(shí)高效、便利的特性。

1 大地電磁二維變分問(wèn)題

研究地下二維電性結(jié)構(gòu),取y軸垂直向上,z軸為走向,x軸向右且與y軸、z軸垂直,求解域Ω為矩形,4個(gè)頂點(diǎn)依次以A,B,C,D順時(shí)針編號(hào),Γ1為地質(zhì)體邊界。大地電磁二維正演滿(mǎn)足(1)式所示的變分問(wèn)題[20]:

(1a)

uAB=1

(1b)

δF(u)=0

(1c)

TE模式為:

u=Ez

(2)

TM模式為:

u=Hz

(3)

式中:ω為角頻率;μ為磁導(dǎo)率;σ為電導(dǎo)率;ε為介電常數(shù);Ez,Hz分別為電場(chǎng)及磁場(chǎng)平行于異常體走向的水平分量。

2 大地電磁二維變分問(wèn)題的無(wú)網(wǎng)格解法

PIM以積分點(diǎn)支持域內(nèi)的場(chǎng)節(jié)點(diǎn)利用多項(xiàng)式基插值構(gòu)造形函數(shù),背景網(wǎng)格用于計(jì)算總體矩陣表達(dá)式中包含的求解域及求解域邊界的積分項(xiàng)。圖1 為PIM支持域、高斯點(diǎn)、背景網(wǎng)格與場(chǎng)節(jié)點(diǎn)示意圖。由于網(wǎng)格積分多采用高斯積分法,故積分點(diǎn)也稱(chēng)高斯點(diǎn)。

圖1 全域弱式無(wú)網(wǎng)格法支持域、高斯點(diǎn)、背景網(wǎng)格與場(chǎng)節(jié)點(diǎn)圖示

2.1 支持域

如圖1所示,支持域形狀常為矩形或圓形,對(duì)于任一高斯點(diǎn)XQ,支持域尺寸d可表示為:

(4)

式中:α為支持域無(wú)量綱尺寸,其大小影響無(wú)網(wǎng)格法的計(jì)算精度[17];dc為支持域內(nèi)高斯點(diǎn)XQ的平均結(jié)點(diǎn)間距,有:

(5)

式中:A為支持域面積;n為包含在A中的節(jié)點(diǎn)總數(shù)。本文采用矩形支持域,x,y方向的支持域尺寸為：

(6)

式中:dcx與dcy分別為x,y方向的節(jié)點(diǎn)間距;αx與αy為對(duì)應(yīng)的無(wú)量綱尺寸。為便于程序設(shè)計(jì),一般令αx=αy=α,本文取α=1。

2.2 FE-PIM離散系統(tǒng)方程的構(gòu)造

PIM計(jì)算效率低但處理復(fù)雜模型便利[19],FEM求解高效,由于PIM與FEM系統(tǒng)矩陣采用相同的方法構(gòu)成,故PIM與FEM具有很好的耦合條件。

如圖2所示,FE-PIM的基本原理是將求解域Ω劃分為實(shí)線(xiàn)表示的FEM區(qū)域Ω1及虛線(xiàn)表示的PIM背景單元區(qū)域Ω2,異常體置于Ω2,無(wú)窮遠(yuǎn)邊界單元設(shè)為Ω1。將變分符號(hào)δ代入(1)式并將FEM區(qū)域單元離散化,有:

(7)

式中:e表示FEM的子單元;CD為求解域底邊界;Γe為邊界單元。

圖2 FE-PIM計(jì)算模型二有限單元與背景網(wǎng)格分布圖示

將FEM區(qū)域的積分項(xiàng)即δF1(u)離散,有:

(8)

(9)

式中:U為全部節(jié)點(diǎn)的場(chǎng)值向量;KFEM為FEM區(qū)域計(jì)算的總體系統(tǒng)矩陣。求得δF1(u)的離散表達(dá)式后,將δF2(u)離散,將u表示為PIM形函數(shù)向量與場(chǎng)節(jié)點(diǎn)向量之積的形式,有

(10)

式中:Φ為用PIM構(gòu)造的形函數(shù)矩陣[18-19];n為支持域內(nèi)的節(jié)點(diǎn)數(shù);u為支持域內(nèi)n個(gè)節(jié)點(diǎn)的場(chǎng)向量值。對(duì)(10)式進(jìn)行變分運(yùn)算,有:

(11)

將(10)式與(11)式代入(1)式的第1項(xiàng)并引入總體編號(hào)體系,總體矩陣最終變?yōu)?

(12)

式中:M1與Mn表示PIM區(qū)域第一個(gè)節(jié)點(diǎn)與最后一個(gè)節(jié)點(diǎn)的編號(hào);KPIM為背景單元域PIM計(jì)算的總體系統(tǒng)矩陣。KPIM表達(dá)式中包含對(duì)求解域Ω2的積分,積分利用高斯積分法求解,有:

(13)

(14)

由于δUT的任意性,(14)式成立的條件變?yōu)?

(15)

(15)式即為用FE-PIM構(gòu)造的系統(tǒng)矩陣,由于FEM與PIM構(gòu)造的總體矩陣KFEM及KPIM均具有帶狀、稀疏、對(duì)稱(chēng)的性質(zhì),故K也具備該特性。

線(xiàn)性方程組KU=0的求解需加載(1)式所示的邊界條件,PIM與FEM邊界條件均可直接加載,計(jì)算時(shí)先找出邊界節(jié)點(diǎn)在總體剛度矩陣中對(duì)應(yīng)的KII,KII對(duì)應(yīng)的行與列設(shè)置為0,其值設(shè)置為1,最后將方程右端項(xiàng)對(duì)應(yīng)的值PI設(shè)置為邊界上的預(yù)設(shè)場(chǎng)值即可。PIM邊界條件可精確加載,該特性是PIM較EFGM最大的優(yōu)勢(shì)。

3 正演計(jì)算

3.1 算法驗(yàn)證

PIM與FEM在一維介質(zhì)大地電磁(MT)正演中具有較高的精度[19],FE-PIM為兩者的耦合,先通過(guò)圖3所示的COMMEMI 2D-0模型[21]來(lái)驗(yàn)證文中算法的正確性。此模型異常體頂部與地表重合,縱向規(guī)模50km,橫向規(guī)模20km,電阻率為1Ω·m,其左、右兩側(cè)電阻率分別為10Ω·m和2Ω·m,深度大于50km區(qū)域?yàn)槌瑢?dǎo)層,電阻率為0,求解域TE模式取120km×110 km,TM模式取120km×70km,場(chǎng)節(jié)點(diǎn)橫縱向等間距分布,節(jié)點(diǎn)間距為1km。

圖3 COMMEMI 2D-0模型

表1列出了PIM與FEM計(jì)算COMMEMI 2D-0模型的視電阻率結(jié)果。從表1可以看出,兩種方法計(jì)算結(jié)果與COMMEMI解析解基本吻合,驗(yàn)證了數(shù)值算法的有效性。此外,表1中FEM計(jì)算結(jié)果也與文獻(xiàn)[22]一致。

表1 COMMEMI 2D-0模型數(shù)值方法視電阻率計(jì)算結(jié)果

3.2 有效性及優(yōu)勢(shì)

為進(jìn)一步驗(yàn)證FE-PIM有效性及其優(yōu)勢(shì),設(shè)計(jì)了如圖4所示的二維地質(zhì)模型:模型二異常體為截面方形二度體,異常體邊長(zhǎng)400m,埋深800m;模型三為水平橢圓,長(zhǎng)半軸300m,短半軸200m,埋深800m;模型四為等腰直角三角形地壘,斜邊長(zhǎng)800m;模型五為出露于地表的條帶狀低阻體,異常體截面呈平行四邊形,下底長(zhǎng)400m,上、下底間距800m,異常體右邊界與地面的夾角呈45°。

模型背景電阻率1000Ω·m,異常體電阻率100Ω·m。場(chǎng)節(jié)點(diǎn)等間距分布于求解域,TE模式3321(41×81)個(gè)場(chǎng)節(jié)點(diǎn),TM模式1681(41×41)個(gè)場(chǎng)節(jié)點(diǎn),橫、縱向節(jié)點(diǎn)間距均為200m,FEM與FE-PIM采用相同的節(jié)點(diǎn)分布。

圖5顯示了頻率為100Hz時(shí)PIM與FEM模型二視電阻率計(jì)算結(jié)果。從圖5可以看出,求解域兩側(cè)邊界的視電阻率值約為1000Ω·m,低阻異常在里程0附近最顯著,TE模式視電阻率值約780Ω·m,TM模式約920Ω·m。兩種模式下FE-PIM,PIM與FEM的計(jì)算結(jié)果基本一致,驗(yàn)證了FE-PIM二維算法的正確性。

圖4 二維地質(zhì)模型

圖5 頻率為100Hz時(shí)模型二視電阻率計(jì)算結(jié)果

圖6為不同模型的FE-PIM視電阻率計(jì)算結(jié)果。從圖6a和圖6b可以看出,TE模式低阻異常中心與模型中心一致,能明顯地反映出橢圓低阻體的存在,但異常規(guī)模比實(shí)際低阻體大,視電阻率為740～1080Ω·m。TM模式低阻異常呈直立狀并無(wú)限向下延伸,異常的水平寬度反映了實(shí)際低阻體水平方向的尺寸,視電阻率為800～1060Ω·m。從圖6c可以看出,模型四山脊位置在視電阻率斷面圖上顯示為高阻極大值區(qū)域,兩側(cè)的山谷區(qū)則呈現(xiàn)低阻極小值,且頻率越高此特性越明顯,視電阻率變化范圍為900～1950Ω·m,此現(xiàn)象表明地形對(duì)大地電磁測(cè)深法影響規(guī)律與直流電阻率法存在差異。斷面圖高阻區(qū)的橫向規(guī)模約800m,與山脊橫向規(guī)模相當(dāng),驗(yàn)證了FE-PIM的有效性。從圖6d 可以看出,模型五的視電阻率斷面圖左右非對(duì)稱(chēng),低阻區(qū)域呈“上窄下寬”條帶狀傾斜分布,傾向與地質(zhì)體相同,視電阻率為100～1100Ω·m,較好地反映出了異常體的物性及空間賦存狀態(tài)。

圖6 不同模型的視電阻率FE-PIM計(jì)算結(jié)果

FE-PIM延續(xù)PIM地質(zhì)建模便利特點(diǎn)的同時(shí)也兼顧了計(jì)算效率,表2列出了17個(gè)頻點(diǎn)PIM,FE-PIM與FEM的計(jì)算耗時(shí)。從表2可以看出,PIM耗時(shí)約為FEM耗時(shí)的8～9倍,FE-PIM耗時(shí)相對(duì)PIM耗時(shí)呈著降低,無(wú)網(wǎng)格區(qū)域(10×4)的FE-PIM與PIM相比，TE模式下效率約提高87%,TM模式下效率約提高83%,隨著FEM區(qū)域的增大,FE-PIM耗時(shí)逐漸接近FEM耗時(shí)。

表2 模型二不同數(shù)值方法計(jì)算耗時(shí)

4 結(jié)束語(yǔ)

1) 將FE-PIM應(yīng)用于大地電磁二維正演模擬,通過(guò)把求解區(qū)域劃分為有限元區(qū)域及背景單元區(qū)域,利用Galerkin法和高斯積分公式導(dǎo)出了大地電磁二維FE-PIM耦合的系統(tǒng)方程離散表達(dá)式,介紹了邊界條件直接加載的技巧,COMMEMI 2D-0 模型的數(shù)值計(jì)算驗(yàn)證了文中PIM與FEM算法的正確性。

2) 采用節(jié)點(diǎn)規(guī)則分布的FE-PIM計(jì)算了截面為方形、橢圓、平行四邊形及三角地壘二維模型的大地電磁場(chǎng)響應(yīng),視電阻率斷面圖較好地反映了地質(zhì)體產(chǎn)生的電磁異常響應(yīng)及空間賦存狀態(tài),驗(yàn)證了耦合算法的正確性,凸顯了FE-PIM處理復(fù)雜地質(zhì)模型的優(yōu)越性。

3) 比較分析了PIM,FE-PIM及FEM的計(jì)算效率,PIM計(jì)算耗時(shí)明顯高于FEM計(jì)算耗時(shí);FE-PIM計(jì)算耗時(shí)受異常體區(qū)域的影響較大,當(dāng)異常體區(qū)域較小時(shí),FE-PIM計(jì)算耗時(shí)遠(yuǎn)小于PIM計(jì)算耗時(shí)而略高于FEM計(jì)算耗時(shí),但FE-PIM綜合了FEM計(jì)算高效及PIM處理復(fù)雜模型便利的優(yōu)點(diǎn),是一種優(yōu)越的數(shù)值方法。

[1] 胡祥云,霍光譜,高銳,等.大地電磁各向異性二維模擬及實(shí)例分析[J].地球物理學(xué)報(bào),2013,56(12):4268-4277 Hu X Y,Huo G P,Gao R,et al.The magnetotelluric anisotropic two-dimensional simulation and case analysis[J].Chinese Journal of Geophysics,2013,56(12):4268-4277

[2] 董浩,魏文博,葉高峰,等.基于有限差分正演的帶地形三維大地電磁反演方法[J].地球物理學(xué)報(bào),2014,57(3):939-952 Dong H,Wei W B,Ye G F,et al.Study of three-dimensional magnetotelluric inversion including surface topography based on finite-difference method[J].Chinese Journal of Geophysics,2014,57(3):939-952

[3] Key K,Weiss C.Adaptive finite-element modeling using unstructured grids:the 2D magnetotelluric example[J].Geophysics,2006,71(6):291-299

[4] 劉長(zhǎng)生,湯井田,任政勇,等.基于非結(jié)構(gòu)化網(wǎng)格的三維大地電磁自適應(yīng)矢量有限元模擬[J].中南大學(xué)學(xué)報(bào)(自然科學(xué)版),2010,41(5):1855-1859 Liu C S,Tang J T,Ren Z Y,et al.Three-dimension magnetotellurics modeling by adaptive edge finite-element using unstructured meshes[J].Journal of Central South University(Science and Technology),2010,41(5):1855-1859

[5] Ren Z Y,Tang J T.3D direct current resistivity modeling with unstructured mesh by adaptive finite-element method[J].Geophysics,2010,75(1):7-17

[6] 柳建新,孫麗影,童孝忠,等.基于自適應(yīng)有限元的起伏地形MT二維正演模擬[J].地球物理學(xué)進(jìn)展,2012,27(5):2016-2023 Liu J X,Sun L Y,Tong X Z,et al.Undulating terrain 2D MT forward modelling with adaptive finite element[J].Progress in Geophysics,2012,27(5):2016-2023

[7] Belytschko T,Lu Y Y,Gu L.Element-free Galerkin methods[J].International Journal for Numerical Methods in Engineering,1994,37(2):229-256

[8] Belytschko T,Lu Y Y,Gu L.Fracture and crack growth by element-free Galerkin methods[J].Modelling and Simulation in Materials Science and Engineering,1994,2(3):519-534

[9] Li D M,Bai F N,Cheng Y M,et al.A novel complex variable element-free Galerkin method for two-dimensional large deformation problems[J].Computer Methods in Applied Mechanics and Engineering,2012,233-236(1):1-10

[10] Liu L C,Dong X H,Li C X.Adaptive finite element-element-free Galerkin coupling method for bulk metal forming processes[J].Journal of Zhejiang University-Science A,2009,10(3):353-360

[11] 李俊杰,嚴(yán)家斌.無(wú)網(wǎng)格法進(jìn)展及其在地球物理學(xué)中的應(yīng)用[J].地球物理學(xué)進(jìn)展,2014,29(1):452-461 Li J J,Yan J B.Developments of meshless method and applications in geophysics[J].Progress in Geophysics,2014,29(1):452-461

[12] 賈曉峰,胡天躍,王潤(rùn)秋.復(fù)雜介質(zhì)中地震波模擬的無(wú)單元法[J].石油地球物理勘探,2006,41(1):43-48 Jia X F,Hu T Y,Wang R Q.A mesh-free algorithm of seismic wave simulation in complex medium[J].Oil Geophysical Prospecting,2006,41(1):43-48

[13] 賈曉峰,胡天躍,王潤(rùn)秋.無(wú)單元法用于地震波波動(dòng)方程模擬與成像[J].地球物理學(xué)進(jìn)展,2006,21(1):11-17 Jia X F,Hu T Y,Wang R Q.Wave equation modeling and imaging by using element-free method[J].Progress in Geophysics,2006,21(1):11-17

[14] 王月英.地震波正演模擬中無(wú)單元Galerkin法初探[J].地球物理學(xué)進(jìn)展,2007,22(5):1539-1544 Wang Y Y.Study of element-free Galerkin method in the seismic forward modeling[J].Progress in Geophysics,2007,22(5):1539-1544

[15] 馮德山,王洪華,戴前偉.基于無(wú)單元Galerkin法探地雷達(dá)正演模擬[J].地球物理學(xué)報(bào),2013,56(1):298-308 Feng D S,Wang H H,Dai Q W.Forward simulation of ground penetrating radar based on the element-free Galerkin method[J].Chinese Journal of Geophysics,2013,56(1):298-308

[16] 戴前偉,王洪華.基于隨機(jī)介質(zhì)模型的GPR無(wú)單元法正演模擬[J].中國(guó)有色金屬學(xué)報(bào),2013,23(9):1-9 Dai Q W,Wang H H.Element free method forward modeling of GPR based on random medium model[J].The Chinese Journal of Nonferrous Metals,2013,23(9):1-9

[17] 嚴(yán)家斌,李俊杰.無(wú)網(wǎng)格法在大地電磁正演計(jì)算中的應(yīng)用[J].中南大學(xué)學(xué)報(bào)(自然科學(xué)版),2014,45(10):3513-3520 Yan J B,Li J J.Magnetotelluric forward calculation by meshless method[J].Journal of Central South University(Science and Technology),2014,45(10):3513-3520

[18] Liu G R,Gu Y T.A point interpolation method for two-dimensional solids[J].International Journal for Numerical Methods in Engineering,2001,50(4):937-951

[19] 李俊杰,嚴(yán)家斌.無(wú)網(wǎng)格點(diǎn)插值法大地電磁二維正演數(shù)值模擬[J].石油物探,2014,53(5):617-626 Li J J,Yan J B.Magnetotelluric two-dimensional for-

ward numerical modeling by meshfree point interpolation method[J].Geophysical Prospecting for Petroleum,2014,53(5):617-626

[20] 戴前偉,張富強(qiáng),楊坤平,等.電導(dǎo)率分塊線(xiàn)性變化的二維高頻大地電磁法有限元正演模擬[J].物探化探計(jì)算技術(shù),2012,34(5):552-558 Dai Q W,Zhang F Q,Yang K P,et al.The finite element method for modeling 2-D high-frequency electromagnetic with continuous variation of conductivity within each block[J].Computing Techniques for Geophysical and Geochemical Exploration,2012,34(5):552-558

[21] 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):133-271

[22] 許建榮.起伏地形條件下大地電磁測(cè)深二維正反演研究及應(yīng)用[D].長(zhǎng)沙:中南大學(xué)地球科學(xué)與信息物理學(xué)院,2010 Xu J R.Research and applications of 2-D MT forward modeling and inversion with topography[D].Changsha:Central South University of Geosciences and Info-Physics,2010

(編輯：顧石慶)

Magnetotelluric two-dimensional forward modelling by finite element-point interpolation coupling method

Li Junjie1,Yan Jiabin2

(1.ZhejiangDesignInstituteofWaterConservancyandHydroelectricPower,Hangzhou310002,China; 2.KeyLaboratoryofNonferrousResourcesandGeologicalHazardDetection,SchoolofGeosciencesandInfo-Physics,CentralSouthUniversity,Changsha410083,China)

Point interpolation method (PIM) is a typical global weak-form meshfree numerical calculation method.The physical properties of PIM are loaded on Gauss integral points which are only related with coordinates during geological modeling.Therefore,PIM is more convenient while dealing complex model than conventional grid method,but the former computation efficiency is low.We couple FEM and PIM into finite element-point interpolation method (FE-PIM) for magnetotelluric two-dimensional forward.Firstly,magnetotelluric two-dimensional discrete system matrix is deduced through substituting the shape function constructed by interpolation method and combining Galerkin method with Gauss integral formula.Then,the principle of background grid integral and the loading technique of boundary conditions are briefly characterized.Finally,the effectiveness and high efficiency of the FE-PIM algorithm and the convenience for complex models are proved by several numerical models.

point interpolation method,global weak-form meshfree method,magnetotelluric,finite element-point interpolation method

2014-12-02;改回日期：2015-02-15。

李俊杰(1989—),碩士,助理工程師,現(xiàn)主要從事大地電磁無(wú)網(wǎng)格化正演模擬研究工作。

國(guó)家自然科學(xué)基金項(xiàng)目(40874055)和湖南省自然科學(xué)基金項(xiàng)目(14JJ2012)共同資助。

P631

A

1000-1441(2015)04-0477-08

10.3969/j.issn.1000-1441.2015.04.015