大地電磁二維正演中的無網格局部徑向基點插值法

何建設，李俊杰＊，嚴家斌

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

大地電磁二維正演中的無網格局部徑向基點插值法

何建設1，李俊杰1＊，嚴家斌2

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

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

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

0 前言

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

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

1 局部徑向基點插值法相關參數

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

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

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

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

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

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

2 局部徑向基點插值算法

取走向為Z軸，X軸與Z軸垂直，Y軸垂直向上，求解域為矩形區(qū)域，四個頂點依次以ABCD順時針編號，Γ1為地質體與周圍介質的邊界。

當平面電磁波以任何角度入射地面時，地下介質中的電磁波總以平面波形式幾乎垂直地向下傳播，當地下電性結構為二維時，滿足如下的邊值問題［12］：

式中：ω為角頻率；μ為磁導率；σ為電導率；ε為介電常數。

2．1 局部徑向基點插值弱式方程

采用LRPIM構造的系統方程，加權殘量要求在局部積分域Ωq上消除即

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

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

式（10）即為與大地電磁二維邊值問題對應的LRPIM系統方程。

2．2 無網格離散系統方程的構造

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

式（12）中的Φ為用徑向基點插值法［6－7，13－15］構造的形函數，其近似原理如式（13）。

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

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

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

2．3 局部積分域數值積分

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

式中：ng與ngΓ分別為局部積分域及邊界域上的高斯點數目；ωi為第i個高斯點XQi的權系數；與為對應的雅可比（Jacobian）矩陣。

LRPIM構造的總體矩陣K具有帶狀、稀疏的性質，但卻是非對稱的。非對稱性主要是因為LRPIM表達式采用不同的函數構造試函數和權函數［16］。總體矩陣的不對稱性增加了計算耗時，這也是LRPIM的主要缺陷。

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

3 正演計算

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

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

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

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

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

圖3 頻率f＝100Hz時模型一數值計算結果Fig．3 Numerical solutions of model 1when frequency is 100Hz

4 結論

1）將LRPIM應用于大地電磁二維正演，介紹了無網格法中相關參數（支持域、積分域、高斯積分點）的基本概念。

2）從大地電磁二維邊值問題出發(fā)，采用子域加權殘量法結合高斯積分公式，推導出了與之對應的LRPIM系統方程及總體矩陣表達式。

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

圖4 頻率f＝100Hz時圖模型二數值計算結果Fig．4 Numerical solutions of model 2when frequency is 100Hz

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

［1］賈曉峰，胡天躍，王潤秋．復雜介質中地震波模擬的無單元法［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］．地球物理學進展，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法探地雷達正演模擬［J］．地球物理學報，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］嚴家斌，李俊杰．無網格法在大地電磁正演計算中的應用［J］．中南大學學報：自然科學版，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］李俊杰，嚴家斌．無網格點插值法大地電磁二維正演數值模擬［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］李俊杰，嚴家斌．無網格法進展及其在地球物理學中的應用［J］．地球物理學進展，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］李俊杰．基于弱式無網格法的大地電磁二維正演［D］．長沙：中南大學，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］戴前偉，張富強，楊坤平，等．電導率分塊線性變化的二維高頻大地電磁法有限元正演模擬［J］．物探化探計算技術，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］李瑩，夏茂輝，董凱．無網格局部徑向點插值法求解Helmholtz方程［J］．鄭州大學學報：理學版，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］嚴默非，楊慶，常文潔．局部徑向基點插值法模擬裂紋尖端附近應力場［J］．武漢理工大學學報，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

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

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

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