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

(1.浙江省水利水電勘測(cè)設(shè)計(jì)院，浙江杭州310002;2.中南大學(xué)地球科學(xué)與信息物理學(xué)院，湖南長(zhǎng)沙410083)

正演模擬是研究大地電磁的基礎(chǔ)手段。當(dāng)?shù)叵码娦越Y(jié)構(gòu)非一維時(shí),大地電磁場(chǎng)響應(yīng)沒(méi)有解析表達(dá)式，必須借助正演數(shù)值模擬方法。作為大地電磁正演計(jì)算的常用網(wǎng)格方法，有限差分法(FDM)[1]、積分方程法(IEM)[2]及有限元法(FEM)[3-5]均有著各自的優(yōu)缺點(diǎn)：有限差分法實(shí)現(xiàn)過(guò)程直接，但無(wú)法處理復(fù)雜地球物理模型;積分方程法只須對(duì)異常體進(jìn)行剖分和求積，不涉及微分方法中的吸收邊界等復(fù)雜問(wèn)題，在三維電磁數(shù)值模擬研究中具有快速、方便等特點(diǎn)，但同樣無(wú)法應(yīng)對(duì)地下電阻率很復(fù)雜時(shí)的計(jì)算;有限元法適用于復(fù)雜物性分布和復(fù)雜邊界形狀的計(jì)算，其最大缺陷在于求解復(fù)雜模型時(shí)網(wǎng)格生成困難。

無(wú)網(wǎng)格法作為網(wǎng)格方法的重要補(bǔ)充和發(fā)展,是近十多年來(lái)興起的一類數(shù)值計(jì)算新方法。其中，無(wú)單元Galerkin法(element-free Galerkin method，EFGM)[6]作為一種基于節(jié)點(diǎn)的無(wú)網(wǎng)格方法，已被成功用于地震波場(chǎng)[7-10]、雷達(dá)波場(chǎng)[11-12]及大地電磁場(chǎng)[13-15]響應(yīng)的計(jì)算。相關(guān)文獻(xiàn)的研究成果表明，EFGM具有較常規(guī)網(wǎng)格方法精度高、計(jì)算復(fù)雜模型便利的特性。然而，EFGM形函數(shù)不滿足克羅內(nèi)克δ函數(shù)性質(zhì)(Ni(X)=δij)，邊界條件加載較困難。無(wú)網(wǎng)格點(diǎn)插值法(MPIM)[16]是一種簡(jiǎn)單高效的無(wú)網(wǎng)格方法，較EFGM最大的改進(jìn)在于形函數(shù)采用插值方法構(gòu)造，邊界條件可精確加載。此方法在彈性力學(xué)領(lǐng)域取得了很好的應(yīng)用效果[16]，但在地球物理學(xué)領(lǐng)域的應(yīng)用研究尚未見(jiàn)報(bào)導(dǎo)。

本文將無(wú)網(wǎng)格點(diǎn)插值法(MPIM)應(yīng)用于大地電磁二維正演，介紹該方法的近似原理，推導(dǎo)大地電磁二維變分問(wèn)題的無(wú)網(wǎng)格化離散過(guò)程;通過(guò)多個(gè)二維理論模型的MPIM，EFGM和FEM正演計(jì)算結(jié)果的對(duì)比，分析MPIM大地電磁二維正演數(shù)值模擬方法的特性和應(yīng)用效果。

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

當(dāng)?shù)叵码娦越Y(jié)構(gòu)為二維時(shí)，取走向?yàn)閦軸，x軸與z軸垂直，y軸垂直向上；求解域Ω為矩形區(qū)域，4個(gè)頂點(diǎn)依次以A,B,C,D順時(shí)針編號(hào),Γ1為地質(zhì)體的邊界(圖1)。

當(dāng)平面電磁波以任何角度入射地面時(shí)，地下介質(zhì)中的電磁波總是以平面波形式幾乎垂直地向下傳播，滿足[17]:

圖1 大地電磁二維正演求解區(qū)域a 電場(chǎng)極化模式(TE模式); b 磁場(chǎng)極化模式(TM模式)

(2)

對(duì)于TM模式:

(3)

式中:ω為角頻率;μ為磁導(dǎo)率;σ為電導(dǎo)率;ε為介電常數(shù)。

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

MPIM利用位于積分點(diǎn)支持域內(nèi)的場(chǎng)節(jié)點(diǎn)構(gòu)造形函數(shù)，需用一組背景網(wǎng)格將求解區(qū)域離散以便建立離散系統(tǒng)方程，如圖2所示。由于背景網(wǎng)格的積分常選用高斯積分法，故積分點(diǎn)又稱高斯積分點(diǎn)或高斯點(diǎn)。

2.1 支持域

常用的支持域形狀有圓形與矩形兩種(圖2)，

圖2 無(wú)網(wǎng)格點(diǎn)插值法的背景網(wǎng)格、支持域、積分點(diǎn)與場(chǎng)節(jié)點(diǎn)

對(duì)于任一高斯點(diǎn)XQ，其支持域尺寸d由下式確定:

(4)

式中:α為支持域的無(wú)量綱尺寸，它用于控制實(shí)際支持域的大小，是對(duì)無(wú)網(wǎng)格法計(jì)算精度影響很大的一個(gè)參數(shù)[16]。dc為位于高斯點(diǎn)XQ附近的平均結(jié)點(diǎn)間距，可由下式確定：

(5)

式中:A為預(yù)估的支持域面積;n為包含在A中的節(jié)點(diǎn)數(shù)。對(duì)于節(jié)點(diǎn)均勻分布的情況，dc為節(jié)點(diǎn)間距。本文采用矩形支持域，故有兩個(gè)方向的支持域尺寸，即

(6)

式中:dcx與dcy分別為橫、縱向節(jié)點(diǎn)間距;αx與αy為對(duì)應(yīng)的支持域無(wú)量綱尺寸。為便于程序設(shè)計(jì)，研究中常取αx=αy=α，本文取α=1.0。

2.2 點(diǎn)插值法近似

求解域Ω中任意一點(diǎn)X處的場(chǎng)變量u(X)的近似表達(dá)式為

(7)

式中:pT(X)為二維空間坐標(biāo)XT=[xy]的基函數(shù);a是系數(shù)向量;m為單項(xiàng)式的個(gè)數(shù)。p(X)可用Pascal三角形確定，對(duì)于一維(1D)和二維(2D)空間，其線性基函數(shù)分別為

pT(X)=1xm=2p=1

(1D)

pT(X)=[1xy]m=3p=1

(2D)

二次基函數(shù)分別為

pT(X)=[1xx2]m=3p=2

(1D)

pT(X)=[1xyx2xyy2]

m=6p=2

(2D)

完備的p階多項(xiàng)式一般可由下式表示

基函數(shù)階次的增加不能顯著提高無(wú)網(wǎng)格法的近似精度，反而降低了計(jì)算效率[8]，一般取線性基即可。將(7)式寫(xiě)成矩陣的形式:

(8)

式中:U=[u1u2…un]T;a=[a1a2…an]T;P的展開(kāi)式為

(9)

傳統(tǒng)PIM中局部支持域內(nèi)的節(jié)點(diǎn)數(shù)總是等于基函數(shù)個(gè)數(shù)(m=n),因此(9)式是一方陣的形式，系數(shù)向量a可通過(guò)矩陣的逆運(yùn)算得到,即

(10)

將(10)式代回(7)式得

(11)

式中，ΦT(X)即為計(jì)算點(diǎn)X在支持域內(nèi)的PIM形函數(shù),其表達(dá)式為

ΦT(X)=pT(X)P-1=

[φ1(X)φ2(X) …φn(X)]

(12)

因?yàn)椴捎枚囗?xiàng)式基函數(shù)，PIM形函數(shù)的求導(dǎo)運(yùn)算較容易，有

(13)

采用多項(xiàng)式作為基函數(shù)時(shí)，PIM插值形式簡(jiǎn)單，并且對(duì)于規(guī)則分布的場(chǎng)節(jié)點(diǎn)的插值具有較高的精度。此外PIM形函數(shù)滿足克羅內(nèi)克δ函數(shù)性質(zhì)，因此邊界條件加載較容易。

然而(9)式的求解是以P非奇異為前提的，當(dāng)多項(xiàng)式基選用不合理或場(chǎng)節(jié)點(diǎn)分布不恰當(dāng)時(shí)，P可能呈現(xiàn)高度病態(tài)甚至呈現(xiàn)奇異性，雖然存在隨機(jī)移點(diǎn)法[16]、局部坐標(biāo)變換法[18]等處理P奇異性的方法，但這些方法均或多或少存在不足之處，限制了多項(xiàng)式PIM的應(yīng)用。

2.3 無(wú)網(wǎng)格離散系統(tǒng)方程的構(gòu)造

定義了無(wú)網(wǎng)格形函數(shù)，便可以構(gòu)造(1)式的離散系統(tǒng)方程，將變分符號(hào)代入(1)式有

(14)

u=[φ1φ2…φn][u1u2…un]T

(15)

式中:Φ為形函數(shù)矩陣;n為支持域內(nèi)的節(jié)點(diǎn)數(shù);u為支持域內(nèi)n個(gè)節(jié)點(diǎn)的場(chǎng)向量。對(duì)(15)式進(jìn)行變分運(yùn)算有

(16)

將(16)式和(15)式代入(14)式，得

(17)

(17)式采用支持域內(nèi)n個(gè)場(chǎng)節(jié)點(diǎn)編號(hào)。對(duì)于求解域內(nèi)所有場(chǎng)節(jié)點(diǎn),還應(yīng)有一套總體編號(hào)體系，用于將局部節(jié)點(diǎn)矩陣組裝成總體剛度矩陣。當(dāng)節(jié)點(diǎn)I和節(jié)點(diǎn)J位于不同支持域時(shí)，其相應(yīng)的被積函數(shù)為0，此時(shí)可得到按求解域節(jié)點(diǎn)編號(hào)的總體編號(hào)體系方程

(18)

(18)式可簡(jiǎn)寫(xiě)為

(19)

由于δUT是任意的，故(19)式成立的條件為

(20)

(20)式即為無(wú)網(wǎng)格點(diǎn)插值法構(gòu)造的系統(tǒng)方程。

2.4 求解域背景網(wǎng)格積分

(20)式中的K包含對(duì)求解域Ω與求解域邊界Γ的積分，為計(jì)算這些積分，須將求解域離散成一組背景網(wǎng)格(圖2)，總體積分可表示成這些單元積分之和的形式，每個(gè)單元的積分利用高斯積分法求解,即

(21)

背景網(wǎng)格積分是MPIM中最重要的數(shù)值計(jì)算問(wèn)題之一，總積分點(diǎn)數(shù)一般取支持域內(nèi)場(chǎng)節(jié)點(diǎn)數(shù)量的2～8倍，但高斯點(diǎn)數(shù)目的增加并不能顯著提高無(wú)網(wǎng)格法的計(jì)算精度，反而極大地降低了計(jì)算效率[15]，因此本文每個(gè)背景單元僅采用4(2×2)個(gè)高斯點(diǎn)，每個(gè)邊界單元采用2個(gè)高斯點(diǎn)。求解線性方程組KU=0還需加載邊界條件，可采用與FEM相同的罰函數(shù)法[16]加載。先將剛度矩陣中相應(yīng)的對(duì)角元素KII變?yōu)棣罧II(α為懲罰系數(shù)，其值可取104～1010),然后用邊界上的uI值取代方程組右端的零向量即可。

3 正演計(jì)算及效果分析

3.1 計(jì)算精度

圖3 模型一

設(shè)計(jì)一個(gè)3層介質(zhì)模型(模型一，圖3)，驗(yàn)證MPIM計(jì)算大地電磁場(chǎng)響應(yīng)的高精度特性。該模型第1層電阻率ρ1=500Ω·m,層厚h1=1km;第2層電阻率ρ2=2000Ω·m,層厚h2=3km;第3層電阻率ρ3=100Ω·m,層厚h3=4km。場(chǎng)節(jié)點(diǎn)等間距分布于問(wèn)題域，對(duì)于TE模式，使用3321(41×81)個(gè)場(chǎng)節(jié)點(diǎn)，3200(40×80)個(gè)背景網(wǎng)格，橫、縱向節(jié)點(diǎn)間距均為200m;對(duì)于TM模式，使用1681(41×41)個(gè)場(chǎng)節(jié)點(diǎn)，1600(40×40)個(gè)背景網(wǎng)格，節(jié)點(diǎn)間距與TE模式相同。分別采用EFGM，MPIM和FEM進(jìn)行正演計(jì)算與對(duì)比，有限元法節(jié)點(diǎn)的分布與無(wú)網(wǎng)格法一致。

表1和表2分別為EFGM，MPIM和FEM計(jì)算的視電阻率及視相位數(shù)值相對(duì)誤差。相對(duì)誤差值的大小對(duì)應(yīng)著數(shù)值方法計(jì)算精度的高低，相對(duì)誤差值越大，精度越低，反之則精度越高。由表1可見(jiàn),3種數(shù)值方法計(jì)算精度相當(dāng)，且精度均隨頻率的增高而降低，精度變化范圍約為3×10-10～4×10-3;當(dāng)頻率f<10-3Hz時(shí)，MPIM精度約高出EFGM一倍，大于此頻點(diǎn)時(shí)MPIM精度與EFGM相當(dāng)。視相位數(shù)值計(jì)算精度變化規(guī)律與視電阻率類似，精度約為2×10-10～4×10-4(表2)。

3.2 有效性

為了進(jìn)一步說(shuō)明MPIM求解MT問(wèn)題的有效性，設(shè)計(jì)圖4所示的二維模型(模型二):圍巖電阻率ρ1=1000Ω·m，方形二度異常體電阻率ρ2=100Ω·m，方形截面邊長(zhǎng)L=400m，異常體頂部到地面的距離h=800m，場(chǎng)節(jié)點(diǎn)的分布與計(jì)算3層介質(zhì)模型時(shí)相同(圖4a)。分別采用EFGM，MPIM和FEM進(jìn)行正演計(jì)算與對(duì)比，圖5為頻率f=100Hz時(shí)模型二的數(shù)值計(jì)算結(jié)果曲線，由圖5可見(jiàn)，MPIM計(jì)算結(jié)果與EFGM和FEM基本一致，視電阻率曲線極小值與視相位曲線極大值的橫坐標(biāo)(X=0)正好對(duì)應(yīng)于方形截面中心在地面的投影，證明了二維情況下MPIM的正確性。

圖6為模型二的FEM，EFGM與MPIM視電阻率計(jì)算結(jié)果斷面圖，可見(jiàn)3種數(shù)值方法計(jì)算結(jié)果仍然一致且對(duì)稱性良好。TE模式下，在100Hz頻點(diǎn)附近，MPIM和EFGM與FEM斷面圖存在細(xì)微差異;TM模式下，形函數(shù)采用插值法構(gòu)造的數(shù)值方法MPIM與FEM計(jì)算結(jié)果一致，較形函數(shù)采用擬合法構(gòu)造的EFGM在10-4～10-3Hz低頻段存在輕微差異。

表1 3層模型視電阻率數(shù)值解的相對(duì)誤差

表2 3層模型視相位數(shù)值解的相對(duì)誤差

圖4 二維理論模型a 模型二; b 模型三; c 模型四; d 模型五

圖5 頻率f=100Hz時(shí)圖模型二的數(shù)值計(jì)算結(jié)果a 視電阻率數(shù)值計(jì)算結(jié)果; b 視相位數(shù)值計(jì)算結(jié)果

3.3 計(jì)算效率

MPIM中模型參數(shù)的加載基于高斯積分點(diǎn)而非單元，高斯點(diǎn)的位置可由坐標(biāo)確定，該特性使其在復(fù)雜地質(zhì)模型構(gòu)建上較常規(guī)網(wǎng)格方法方便。為表明MPIM的這一優(yōu)勢(shì)，采用與前文相同的節(jié)點(diǎn)分布計(jì)算了如圖4b至圖4d所示的二維模型:模型三為圓形截面低阻體，電阻率ρ2=100Ω·m，圍巖電阻率ρ1=1000Ω·m，低阻體截面半徑R=200m，異常體頂部到地面的距離h=800m;模型四與模型五為出露于地表的條帶狀低阻體，異常體截面呈平行四邊形，下底長(zhǎng)a=400m，上下底間距d=600m，異常體右邊界與地面的夾角分別呈45°與30°，其電阻率及圍巖電阻率與模型二相同。

圖7為模型三的MPIM和EFGM視電阻率計(jì)算結(jié)果,由圖7可見(jiàn)，TE模式異常橫向拉長(zhǎng)，TM模式異?？v向拉長(zhǎng)，兩種模式下數(shù)值計(jì)算結(jié)果均很好地反映了低阻異常體的存在。

圖6 模型二視電阻率數(shù)值計(jì)算結(jié)果a TE模式FEM視電阻率; b TM模式FEM視電阻率; c TE模式EFGM視電阻率; d TM模式EFGM視電阻率; e TE模式MPIM視電阻率; f TM模式MPIM視電阻率

圖7 模型三視電阻率數(shù)值計(jì)算結(jié)果a TE模式EFGM視電阻率; b TM模式EFGM視電阻率; c TE模式MPIM視電阻率; d TM模式MPIM視電阻率

圖8為模型四與模型五的TE模式MPIM和EFGM計(jì)算結(jié)果，橫坐標(biāo)為里程，縱坐標(biāo)為以2為底的對(duì)數(shù)工作頻率。由圖8可見(jiàn)，①模型四與模型五的視電阻率斷面呈現(xiàn)左右不對(duì)稱，低阻區(qū)域呈“上窄下寬”狀分布，模型四視電阻率變化范圍均為100～1060Ω·m，異常特征顯著;②低阻區(qū)域呈傾斜條帶狀分布，傾向與地質(zhì)體的傾向相同;③模型五視電阻率變化范圍均為100～1160Ω·m，其斷面較模型四低阻條帶呈現(xiàn)的傾角更大，且EFGM異常區(qū)斷面較MPIM更為寬泛，如圖8c和圖8d 所示)。

以上模型試算分析可見(jiàn)EFGM和MPIM的高精度特性及其在計(jì)算復(fù)雜地質(zhì)模型上的優(yōu)越性。雖然無(wú)網(wǎng)格法與有限元法均達(dá)到精度要求，但在計(jì)算耗時(shí)上卻存在較大差異,表3列出了3種數(shù)值算法計(jì)算模型二的耗時(shí)，工作頻率為10-4～103Hz，共17個(gè)頻點(diǎn)。由表3可見(jiàn)：①3種數(shù)值算法的計(jì)算速度由快到慢依次為FEM，MPIM，EFGM，EFGM計(jì)算耗時(shí)約為FEM的10倍;②MPIM速度較EFGM快，TM模式下效率提高8%，TE模式下提高5%。

無(wú)網(wǎng)格法(MPIM和EFGM)計(jì)算的耗時(shí)主要花費(fèi)在形函數(shù)的構(gòu)建及數(shù)值積分上[15]，計(jì)算效率的提升是無(wú)網(wǎng)格法廣泛應(yīng)用需解決的首要問(wèn)題，無(wú)網(wǎng)格方法與高效網(wǎng)格方法的耦合或許可成為一個(gè)可行的方向。

圖8 TE模式模型四與模型五視電阻率無(wú)網(wǎng)格法計(jì)算結(jié)果a 模型四EFGM視電阻率計(jì)算結(jié)果; b模型四MPIM視電阻率計(jì)算結(jié)果; c 模型五EFGM視電阻率計(jì)算結(jié)果; d 模型五MPIM視電阻率計(jì)算結(jié)果

s

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

復(fù)雜地質(zhì)模型電磁響應(yīng)的計(jì)算一般需采用非結(jié)構(gòu)化網(wǎng)格，此種網(wǎng)格的生成算法較復(fù)雜。本文將MPIM應(yīng)用于大地電磁二維正演數(shù)值模擬，介紹了MPIM的基本原理及大地電磁二維變分問(wèn)題的無(wú)網(wǎng)格化求解過(guò)程;采用節(jié)點(diǎn)規(guī)則分布的無(wú)網(wǎng)格法(MPIM和EFGM)計(jì)算出了柱形二度體及傾斜條帶狀二度體模型的電磁響應(yīng)。理論模型的MPIM，EFGM和FEM正演計(jì)算結(jié)果的對(duì)比分析表明:MPIM適用于大地電磁正演計(jì)算，其計(jì)算精度很高，數(shù)值算例表明精度最高可達(dá)10-9數(shù)量級(jí);MPIM與EFGM精度相當(dāng)，但MPIM的計(jì)算效率更高，本文17個(gè)頻點(diǎn)下的數(shù)值計(jì)算表明TM模式下效率提高8%，TE模式下提高5%。研究結(jié)果驗(yàn)證了MPIM大地電磁二維正演數(shù)值模擬方法的有效性，體現(xiàn)了MPIM在處理復(fù)雜地質(zhì)模型正演問(wèn)題上的優(yōu)越性。

參 考 文 獻(xiàn)

[1] 譚捍東，余欽范，Booker J，等.大地電磁法三維交錯(cuò)采樣有限差分?jǐn)?shù)值模擬[J].地球物理學(xué)報(bào),2003,46(5):705-711

Tan H D,Yu Q F,Booker J,et al.Magnetotelluric three-dimensional modeling using the staggered-grid finite difference method[J].Chinese Journal of Geophysics (in Chinese),2003,46(5):705-711

[2] Wannamaker P E.Advances in three-dimensional magnetotelluric modeling using integral equations[J].Geophysics,1991,56(11):1716-1728

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

[7] 賈曉峰,胡天躍.滑動(dòng)最小二乘法求解地震波波動(dòng)方程[J].地球物理學(xué)進(jìn)展，2005，20(4):920-924

Jia X F,Hu T Y.Solving seismic wave equation by moving least squares (MLS) method[J].Progress in Geophysics,2005,20(4):920-924

[8] 賈曉峰,胡天躍,王潤(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

[9] 賈曉峰,胡天躍,王潤(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 model-

ing and imaging by using element-free method[J].Progress in Geophysics,2006，21(1):11-17

[10] 王月英.地震波正演模擬中無(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

[11] 馮德山,王洪華,戴前偉.基于無(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

[12] 戴前偉,王洪華.基于隨機(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

[13] 蘇洲,胡文寶,朱毅.二維大地電磁正演中無(wú)網(wǎng)格算法研究[J].石油天然氣學(xué)報(bào),2012,34(5):87-90

Su Z,Hu W B,Zhu Y.Meshfree method used in two-dimensional magnetotelluric forwarding[J].Journal of Oil and Gas Technology,2012,34(5):87-90

[14] 李俊杰,嚴(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(in Chinese),2014,29(1):452-461

[15] 李俊杰.基于弱式無(wú)網(wǎng)格法的大地電磁二維正演[D].長(zhǎng)沙:中南大學(xué)地球科學(xué)與信息物理學(xué)院,2014

Li J J.Magnetotelluric two-dimensional forward by weak-form meshfree methods[D].Changsha:Central South University of Geosciences and Info-Physics (in Chinese),2014

[16] 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

[17] 徐世浙.地球物理中的有限單元法[M].北京:科學(xué)出版社，1994:229-235

Xu S Z.Finite element method for geophysics[M].Beijing:Science Press,1994:229-235

[18] Liu G R.Meshfree methods:moving beyong the finite element method[M].USA:CRC Press,2002:123-126