何建設(shè)，李俊杰＊，嚴(yán)家斌

（1．浙江省水利水電勘測設(shè)計(jì)院，杭州 310002；2．中南大學(xué)地球科學(xué)與信息物理學(xué)院，長沙 410083）

大地電磁二維正演中的無網(wǎng)格局部徑向基點(diǎn)插值法

何建設(shè)1，李俊杰1＊，嚴(yán)家斌2

（1．浙江省水利水電勘測設(shè)計(jì)院，杭州 310002；2．中南大學(xué)地球科學(xué)與信息物理學(xué)院，長沙 410083）

無單元Galerkin法作為較成熟的一種無網(wǎng)格方法，已成功應(yīng)用于有限元法觸及的領(lǐng)域，還解決了如大變形、裂紋擴(kuò)展及高速?zèng)_擊等網(wǎng)格方法較難處理的問題，但其最大的缺陷在于系統(tǒng)方程的離散需借助背景網(wǎng)格，因此該方法并非真正意義上的無網(wǎng)格方法。無網(wǎng)格局部徑向基點(diǎn)插值法采用子域法構(gòu)造系統(tǒng)方程，加權(quán)殘量只要求在局部積分域消除，大大降低了對(duì)背景網(wǎng)格的依賴，向真正的無網(wǎng)格邁進(jìn)了一大步．這里將此方法用于大地電磁二維正演，介紹了該方法的基本原理；從大地電磁二維邊值問題出發(fā)，利用子域法推導(dǎo)了與之對(duì)應(yīng)的無網(wǎng)格局部弱式系統(tǒng)方程，并用高斯積分將其離散化；論述了局部徑向基點(diǎn)插值法較無單元Galerkin法及有限元法的優(yōu)缺點(diǎn)；最后通過二維模型的計(jì)算驗(yàn)證了算法的有效性。

局部徑向基點(diǎn)插值法；大地電磁；無單元Galerkin法；有限元法

0 前言

正演是研究大地電磁的基礎(chǔ)，高維問題的大地電磁場響應(yīng)不存在解析解，為此須借助數(shù)值方法。有限單元法、有限差分法、積分方程法作為地球物理電磁法領(lǐng)域常用的網(wǎng)格方法，其最大缺陷在于求解復(fù)雜模型時(shí)網(wǎng)格生成困難。無網(wǎng)格法利用加權(quán)殘量法構(gòu)造求解問題的弱式（積分）形式或配點(diǎn)（微分）形式，并通過支持域內(nèi)的節(jié)點(diǎn)構(gòu)造形函數(shù)及建立離散系統(tǒng)方程，是一種基于節(jié)點(diǎn)的數(shù)值方法。因其無網(wǎng)格、精度高、模型加載便利等優(yōu)點(diǎn)，成為了計(jì)算力學(xué)領(lǐng)域的研究熱點(diǎn)，在地球物理學(xué)領(lǐng)域也有少量應(yīng)用，報(bào)導(dǎo)的文獻(xiàn)主要集中于無網(wǎng)格全域弱式Galerkin法［1－5］。

局部徑向基點(diǎn)插值法（local radial point inter-polation method，LRPIM）［6－7］作為無網(wǎng)格弱式方法的一種較全域弱式無網(wǎng)格法（如無單元Galerkin法（Element-free Galerkin method，EFGM）［8－10］、點(diǎn)插值法（point interpolation method，PIM）［10－11］）其最大的改進(jìn)在于數(shù)值積分只在局部積分域上進(jìn)行，不涉及背景網(wǎng)格的生成，因而向真正的無網(wǎng)格法邁進(jìn)了一大步。這里介紹了LRPIM的基本原理及大地電磁二維問題的LRPIM求解過程，通過兩個(gè)二維模型的計(jì)算，證明了算法的有效性及其在處理復(fù)雜模型問題上的優(yōu)越性。

1 局部徑向基點(diǎn)插值法相關(guān)參數(shù)

LRPIM利用位于積分點(diǎn)支持域內(nèi)的場節(jié)點(diǎn)構(gòu)造形函數(shù)，數(shù)值積分在局部積分域上進(jìn)行，因此先介紹支持域及積分域的概念。圖1為LRPIM支持域、積分域、積分點(diǎn)與場節(jié)點(diǎn)示意圖，由于局部積分域積分常選用高斯積分法，故積分點(diǎn)又稱高斯積分點(diǎn)或高斯點(diǎn)。

圖1 LRPIM支持域、積分域、積分點(diǎn)與場節(jié)點(diǎn)示意圖Fig．1 Support domain，quadrature domain，gauss points and nodes in LRPIM

支持域與積分域情況如圖1所示，常用的支持域及積分域形狀有圓形與矩形兩種，對(duì)于任一高斯點(diǎn)Χ，其支持域尺寸d及積分域尺寸dq由式（1）確定。

式中：α與αq為支持域及積分域的無量綱尺寸，它們用于控制實(shí)際支持域與積分域的大小，它們都是對(duì)LRPIM計(jì)算精度有顯著影響的重要參數(shù)［7］；dc為位于高斯點(diǎn)Χ附近的平均結(jié)點(diǎn)間距，可由式（2）確定。

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

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

2 局部徑向基點(diǎn)插值算法

取走向?yàn)閆軸，X軸與Z軸垂直，Y軸垂直向上，求解域?yàn)榫匦螀^(qū)域，四個(gè)頂點(diǎn)依次以ABCD順時(shí)針編號(hào)，Γ1為地質(zhì)體與周圍介質(zhì)的邊界。

當(dāng)平面電磁波以任何角度入射地面時(shí)，地下介質(zhì)中的電磁波總以平面波形式幾乎垂直地向下傳播，當(dāng)?shù)叵码娦越Y(jié)構(gòu)為二維時(shí)，滿足如下的邊值問題［12］：

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

2．1 局部徑向基點(diǎn)插值弱式方程

采用LRPIM構(gòu)造的系統(tǒng)方程，加權(quán)殘量要求在局部積分域Ωq上消除即

式（7）中：W是以節(jié)點(diǎn)為中心的權(quán)函數(shù)，與全域Galerkin弱式法不同，LRPIM權(quán)函數(shù)W與試函數(shù)u一般取自不同的函數(shù)空間，權(quán)函數(shù)取四次樣條函數(shù)。式（7）可展開為式（8）的形式。

式（9）在邊界Γ上代入式（4）所述的邊界條件可得式（10）。

式（10）即為與大地電磁二維邊值問題對(duì)應(yīng)的LRPIM系統(tǒng)方程。

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

求解式（10）需先將其離散，將場量u表示為節(jié)點(diǎn)處場量值與形函數(shù)之積的形式有

式（12）中的Φ為用徑向基點(diǎn)插值法［6－7，13－15］構(gòu)造的形函數(shù)，其近似原理如式（13）。

式中：Ri（x，y）為MQ（multi－quadrics）函數(shù)；p（X）為用Pascal三角形確定的基函數(shù)；a與b是系數(shù)向量；n為RBF的個(gè)數(shù)，m為單項(xiàng)式的個(gè)數(shù)；dc為與節(jié)點(diǎn)間距有關(guān)的特征長度，當(dāng)節(jié)點(diǎn)均勻分布時(shí)可取dc＝；x與y為高斯點(diǎn)位置的坐標(biāo)，xi與yi為支持域內(nèi)節(jié)點(diǎn)的坐標(biāo)；αc與q為對(duì)LRPIM計(jì)算精度有較大影響的形狀參數(shù)，一般通過數(shù)值試驗(yàn)獲得，文中取αc＝2．0，q＝0．5。

式（11）采用支持域內(nèi)n個(gè)場節(jié)點(diǎn)編號(hào)，對(duì)于求解域內(nèi)所有場節(jié)點(diǎn)還應(yīng)有一個(gè)用于將局部節(jié)點(diǎn)矩陣組裝成總體剛度矩陣的總體編號(hào)體系，即從1到N編號(hào)，因此式（11）變?yōu)槭剑?4）。

值得說明的是，LRPIM總體矩陣K的加載與無網(wǎng)格全域Galerkin弱式法及傳統(tǒng)有限單元法不同，在FEM與全域Galerkin弱式法中，單元矩陣或節(jié)點(diǎn)矩陣是被系統(tǒng)地累加進(jìn)總體矩陣的，而在LRPIM中，節(jié)點(diǎn)矩陣是按行擺放而形成總體矩陣的。

2．3 局部積分域數(shù)值積分

式（14）中K表達(dá)式中包含的局部積分，可利用高斯積分法求解有

式中：ng與ngΓ分別為局部積分域及邊界域上的高斯點(diǎn)數(shù)目；ωi為第i個(gè)高斯點(diǎn)XQi的權(quán)系數(shù)；與為對(duì)應(yīng)的雅可比（Jacobian）矩陣。

LRPIM構(gòu)造的總體矩陣K具有帶狀、稀疏的性質(zhì)，但卻是非對(duì)稱的。非對(duì)稱性主要是因?yàn)長RPIM表達(dá)式采用不同的函數(shù)構(gòu)造試函數(shù)和權(quán)函數(shù)［16］?？傮w矩陣的不對(duì)稱性增加了計(jì)算耗時(shí)，這也是LRPIM的主要缺陷。

求解線性方程組KU＝0還需加載邊界條件，常用邊界條件的加載方法有罰函數(shù)法［4－5］和Lagrange乘子法［8］。由于LRPIM形函數(shù)滿足Kronecker delta函數(shù)特性（Ni（X）＝δij），邊界條件也可直接加載［7］。

3 正演計(jì)算

為了研究算法的應(yīng)用效果，計(jì)算了二維模型（圖2）：模型一為正方形，背景電阻率為1 000 Ω·m，異常體電阻率為100Ω·m，邊長為400m，異常體頂部到地面的距離為800m；模型二為圓形，半徑為200m；模型三與四為橢圓，長半軸為300 m，短半軸為200m，前者為直立橢圓，后者為水平橢圓。計(jì)算時(shí)采用1 681（41×41）個(gè)場節(jié)點(diǎn)，節(jié)點(diǎn)等間距分布于問題域。

圖3為頻率f＝100Hz時(shí)模型一的數(shù)值計(jì)算結(jié)果，由圖3可知，LRPIM計(jì)算結(jié)果與FEM及另兩種全域Galerkin弱式無網(wǎng)格方法（EFG、PIM）一致，視電阻率與視相位曲線均很好地反映出了方形異常體的存在，證明了LRPIM求解大地電磁二維問題的有效性。

無網(wǎng)格法是一類基于節(jié)點(diǎn)的數(shù)值算法，其模型參數(shù)的加載基于高斯積分點(diǎn)而非單元，高斯點(diǎn)的位置可由坐標(biāo)確定，該特性使其在復(fù)雜模型的構(gòu)建上較常規(guī)網(wǎng)格方法方便。圖4為模型二的無網(wǎng)格法計(jì)算結(jié)果，由圖4可知，三種方法的異常曲線只在里程0km附近有細(xì)微的差別，均較好地反映出了異常體的存在，LRPIM異常幅值較其余兩種無網(wǎng)格全域弱式法大。

圖5為模型三與四的視電阻率計(jì)算結(jié)果，LRPIM與EFG、PIM計(jì)算結(jié)果基本一致，但曲線更窄，里程0km附近水平橢圓異常幅值較直立橢圓大（圖5）。

圖2 二維理論模型Fig．2 Two－dimensional theoretical models

圖3 頻率f＝100Hz時(shí)模型一數(shù)值計(jì)算結(jié)果Fig．3 Numerical solutions of model 1when frequency is 100Hz

4 結(jié)論

1）將LRPIM應(yīng)用于大地電磁二維正演，介紹了無網(wǎng)格法中相關(guān)參數(shù)（支持域、積分域、高斯積分點(diǎn)）的基本概念。

2）從大地電磁二維邊值問題出發(fā)，采用子域加權(quán)殘量法結(jié)合高斯積分公式，推導(dǎo)出了與之對(duì)應(yīng)的LRPIM系統(tǒng)方程及總體矩陣表達(dá)式。

3）截面方形二度體模型的LRPIM計(jì)算結(jié)果與FEM一致，驗(yàn)證了LRPIM在大地電磁二維正演中的有效性。復(fù)雜地質(zhì)模型電磁響應(yīng)的計(jì)算一般需采用非結(jié)構(gòu)化網(wǎng)格，此種網(wǎng)格的生成算法較復(fù)雜。作者采用節(jié)點(diǎn)規(guī)則分布的LRPIM計(jì)算了圓與橢圓二維模型的電磁響應(yīng)，計(jì)算結(jié)果較好地反應(yīng)出了異常體的存在，體現(xiàn)了LRPIM在計(jì)算復(fù)雜模型上的優(yōu)越性。

圖4 頻率f＝100Hz時(shí)圖模型二數(shù)值計(jì)算結(jié)果Fig．4 Numerical solutions of model 2when frequency is 100Hz

圖5 頻率f＝100Hz時(shí)圖模型三與模型四的視電阻率數(shù)值計(jì)算結(jié)果Fig．5 Numerical solutions of apparent resistivity for model 3and model 4when frequency is 100Hz

［1］賈曉峰，胡天躍，王潤秋．復(fù)雜介質(zhì)中地震波模擬的無單元法［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．（In Chinese）

［2］王月英．地震波正演模擬中無單元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．（In Chiness）

［3］馮德山，王洪華，戴前偉．基于無單元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．（In Chinese）

［4］嚴(yán)家斌，李俊杰．無網(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．（In Chinese）

［5］李俊杰，嚴(yán)家斌．無網(wǎng)格點(diǎn)插值法大地電磁二維正演數(shù)值模擬［J］．石油物探，2014，53（5）：617－626．

LI J J，YAN J B．Magnetotelluric two－dimensional forward numerical modeling by meshfree point interpolation method［J］．Geophysical Prospecting for Petroleum，2014，53（5）：617－626．（In Chinese）

［6］LIU G R，GU Y T．A local radial point interpolation method（LRPIM）for free vibration analyses of 2－D solids［J］．Journal of Sound and Vibration，2001，246 （1）：29－46．

［7］LIU G R，YAN L，WANG J G，Gu Y T．Point interpolation method based on local residual formulation u-sing radial basis functions［J］．Structural Engineering and Mechanics，2002，l4（6）：7l3－732．

［8］BELYTSCHKO T，LU YY，GU L．Element－free Galerkin methods［J］．International Journal for Numerical Methods in Engineering，1994，37（2），229－256．

［9］李俊杰，嚴(yán)家斌．無網(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．（In Chinese）

［10］李俊杰．基于弱式無網(wǎng)格法的大地電磁二維正演［D］．長沙：中南大學(xué)，2014．

LI J J．Magnetotelluric two－dimensional forward by weak－form meshfree methods［D］．Changsha：Central South University，2014．（In Chinese）

［11］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．

［12］戴前偉，張富強(qiáng)，楊坤平，等．電導(dǎo)率分塊線性變化的二維高頻大地電磁法有限元正演模擬［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 conduc-tivity within each block［J］．Computing Techniques for Geophysical and Geochemical Exploration，2012，34 （5）：552－558．（In Chinese）

［13］WANG J G，LIU G R．A point interpolation meshless method based on radial basis functions［J］．International Journal for Numerical Methods in Engineering，2002，54（11）：1623－1648．

［14］李瑩，夏茂輝，董凱．無網(wǎng)格局部徑向點(diǎn)插值法求解Helmholtz方程［J］．鄭州大學(xué)學(xué)報(bào)：理學(xué)版，2012，44 （4）：26－30．

LI Y，XIA M H，DONG K．A meshless local radial point interpolation method for the Helmholtz equation ［J］．Journal of Zhengzhou University：Natural Science Edition，2012，44（4）：26－30．（In Chinese）

［15］嚴(yán)默非，楊慶，常文潔．局部徑向基點(diǎn)插值法模擬裂紋尖端附近應(yīng)力場［J］．武漢理工大學(xué)學(xué)報(bào)，2010，32 （5）：141－145．

YAN M F，YANG Q，CHANG W J．Local radialbasis point interpolation method for crack tip stress field［J］．Journal Of Wuhan University of Technology，2010，32（5）：141－145．（In Chinese）

［16］ATLURI S N，KIM H G，CHO J Y．A critical assessment of the truly meshless local Petrov－Galerkin （MLPG）and Local Boundary Integral Equation （LBIE）methods［J］．Computational Mechanics，1999，24（5）：348－372．

Two－dimensional forward of magnetotelluric using meshless local radial point interpolation method

HE Jian-she1，LI Jun-jie1＊，YAN Jia-bin2
（1．Zhejiang Design Institute of Water Conservancy and Hydroelectric Power，Hangzhou 310002，China；2．School of Geosciences and Info－Physics，Central South University，Changsha 410083，China）

Element－free Galerkin method（EFGM）as a relatively mature meshless method not only has been successfully applied in the field of where the finite element method（FEM）used，but also addressed many problems which grid method is difficult to deal with such as large－deformation，fracture and crack growth and high speed impact．Its biggest drawback is that discrete process of system equation involves the background grid generation therefore cannot be called the true sense of meshless method．Subdomain method is used to construct the system equation in meshless local radial point interpolation method （LRPIM），the weighted residual only be demanded eliminating in local integral domain greatly reduces the dependence on the background grid therefore LRPIM forwards a major step to the true meshless method．This paper devotes two－dimensional magnetotelluric forward by LRPIM and its basic principle is introduced．Weak form equation is derived by subdomain method corresponding to the magnetotelluric two－dimensional boundary value problem then the equation is discrete by gauss integral method．The advantages and disadvantages of LRPIM are analyzed comparing to EFGM and FEM．At last，the validity of the algorithm is verified by two－dimensional model calculations．

local radial point interpolation method；magnetotelluric；element-free Galerkin method；finite element method

P 631．3

A

10．3969／j．issn．1001-1749．2015．03．01

1001-1749（2015）03-0267-06

2014-08-01 改回日期：2014-10-13

國家自然科學(xué)基金（40874055）；湖南省自然科學(xué)基金（07JJ5065）

何建設(shè)（1970－），男，高級(jí)工程師，主要從事水電勘測與電磁法正演研究，E-mail：751141742＠qq．com。

＊通信作者：李俊杰（1989－），碩士，主要從事地球物理場無網(wǎng)格化正演研究，E-mail：838885421＠qq．com。