孔 倩, 李 鵬, 李 晶

( 1.華北電力大學 數(shù)理學院,涿州 071003 2.中國石油 東方地球物理公司,保定 072751 )

二維位場延拓的無網(wǎng)格方法研究

孔 倩1, 李 鵬1, 李 晶2

( 1.華北電力大學 數(shù)理學院,涿州 071003 2.中國石油 東方地球物理公司,保定 072751 )

位場向上延拓可歸納為求解函數(shù)所滿足的Laplace方程，將無網(wǎng)格Galerkin(EFG)方法推廣到二維位場延拓的數(shù)值計算中，詳細論述了EFG方法的基本原理和具體實施過程，建立了無網(wǎng)格Galerkin法求解的離散方程，并進行了數(shù)值求解。同時與有限差分(FD)法的數(shù)值結(jié)果進行了比較，兩種方法的求解結(jié)果基本吻合。介紹了EFG方法的應用實例，理論模型和實例的數(shù)值結(jié)果表明，EFG方法在處理二維位場延拓問題時是有效的且具有實現(xiàn)簡單的特點。

位場延拓; 無網(wǎng)格Galerkin方法; 有限差分法; 正演計算

0 引言

在地球物理中，位場延拓是許多位場數(shù)據(jù)處理和正演計算中非常重要的步驟。二維位場的向上延拓通常會求解橢圓型的偏微分方程，已有的數(shù)值延拓方法包括有限差分法、有限元法、邊界元法、積分方程法[1-9]。鄧正棟等[1]采用有限差分法對三維穩(wěn)定地電場進行了正演模擬，計算得出了地電異常體引起的異常電位場；尤淼等[2]針對二維位場延拓問題，探討了有限差分方程的建模和計算過程，有限差分法求解過程簡明而直觀，但求解速度受網(wǎng)格剖分大小的影響，且數(shù)值穩(wěn)定性和收斂性得不到保證[4]；郭榮文[5]從位場理論和有限元方法角度探討了磁場向上延拓與模擬的有限元算法問題。有限元法可以處理復雜的幾何形狀，劃分網(wǎng)格時更具有靈活性和適應性，比較容易地計算彎曲的物性界面和起伏地形，但輸入輸出的數(shù)據(jù)較多，對計算機的要求很高，且不適于處理無界區(qū)域[4]。電法勘探的邊界元方法以徐世浙[7-9]的研究成果在國際上得到很高的評價，采用邊界單元法對起伏地形二維位場向上延拓進行了數(shù)值模擬[7]，并模擬了二維地形對大地電磁場的影響[8]。邊界元法只在研究區(qū)域的邊界上剖分網(wǎng)格，計算量減少，但是對于邊界元算法得到的大型矩陣方程往往是滿秩的病態(tài)方程組[4]。因此針對位場延拓遇到的偏微分方程應研究有效的求解方法，有著非常重要的意義。

無網(wǎng)格方法是近年來迅速興起的一種數(shù)值方法, 該方法不需要生成網(wǎng)格，而是按照一些任意分布的坐標點構(gòu)造插值函數(shù)離散控制方程，可方便地模擬各種復雜形狀的流場。研究表明[10-13]，無網(wǎng)格Galerkin (Element Free Galerkin，EFG) 方法具有收斂速度快，精度較高,計算穩(wěn)定等優(yōu)點，是無網(wǎng)格方法中較成熟的一種方法。因此，我們將EFG方法推廣到二維位場延拓的數(shù)值計算中，討論EFG方法模擬位場延拓問題的可行性。

1 基礎(chǔ)介紹

建立一個三維的笛卡兒坐標系，將坐標系的xoy平面定為觀測平面， 而其z軸向上為正。假定T(x,y,0)代表觀測平面上的實際觀測位場。于是，位場的向上延拓就是從已知的位場T(x,y,0)，推導出任意高度z(z>0)的位場T(x,y,z)，這就是根據(jù)位場的拉普拉斯方程式來確定著名的Dirichlet問題的解：

▽2T(x,y,z)=0

(1)

式中的變量z是觀測平面上的延拓距離。

(2)

對式(2)擬采用無網(wǎng)格方法進行數(shù)值模擬，邊界采用罰函數(shù)法處理。

2 EFG方法的基本原理[11-12]

EFG方法采用Galerkin變分原理對控制方程進行空間離散，對未知函數(shù)使用移動最小二乘法進行近似。由移動最小二乘法，待求函數(shù)u(x)在點x的鄰域內(nèi)可近似為式(3)。

(3)

(4)

式中：w(x-xI)為節(jié)點的權(quán)函數(shù)；n為節(jié)點數(shù)目。

使式(4)中的J在駐點x達到最小，則可得到待定系數(shù)ai(x)滿足如下方程組

(5)

所以a(x)=A-1(x)B(x)u

式中：A=pTW(x)p，B=pTW(x)，u=(u1,u2,……,un)T， 將a(x)代入式(2)，則

(6)

3 模型求解

對式(2)采用上述EFG方法，取φi、φj為形函數(shù)ΦI(x)，邊界條件采用罰函數(shù)方法[10]處理，α為罰因子，得到相應的變分形式為式(7)。

(7)

將式(7)，表述成下列矩陣形式：

KT=q

(8)

其中：

(9)

(10)

4 方法應用

4.1 理論模型1

有一水平地面測線，包含15個測點的場值信息。其中，地面邊界上2號點～14號點磁場100nT，其余兩點磁場為0nT(圖1)。

圖1 模型1中測點的場值信息
Fig.1 Field value in model 1

將模型1中x方向水平地面長度定為300 m，高度定為600 m。將上述問題歸結(jié)為如下Laplace方程的求解：

(11)

用分離變量法求出的上述邊界條件下的解析解：

(12)

同時，分別采用EFG方法和有限差分(FiniteDifference,FD)方法求解上述方程。EFG方 法求解空間等間距均勻布置15×29個節(jié)點，四個本質(zhì)邊界條件節(jié)點上的罰因子取為1015。圖2給出了模型1中采用EFG方法得出的整個求解區(qū)域的磁場曲面，驗證了無網(wǎng)格方法在求解位場延拓問題中的可行性。圖3給出了模型1中分別采用EFG方法和FD方法得到的磁場空間誤差圖。由圖2、圖3可見，相對于FD方法，EFG方法的誤差較小，并且相對于FD方法，EFG方法基于點的近似不需要網(wǎng)格的初始劃分，減少了計算難度。

圖2 EFG方法的數(shù)值解磁場曲面Fig.2 Magnetic field surface by EFG method

圖3 模型1中EFG方法和FD方法的磁場空間誤差Fig.3 Magnetic field space error by EFG and FD method in model 1(a) EFG方法磁場空間誤差圖；(b)FD方法磁場空間誤差圖

4.2 理論模型2

有一水平地面測線，同樣15個測點信息。地面邊界上2號點～9號點磁場為100nT，10號點～14號點磁場為10nT,其余兩點磁場為0nT,如圖4所示。

圖4 模型2中測點的場值信息
Fig.4 Field value in model 2

將模型2中x方向水平地面長度定為300 m，高度定為600 m。 圖4上部空間等間距均勻布置15×29個節(jié)點，四個本質(zhì)邊界條件節(jié)點上的罰因子取為1015。圖5給出了模型2中分別采用EFG方法和FD方法得到的磁場空間等值線。橫軸、縱軸分別代表x、z軸的節(jié)點編號。由圖4可見， EFG方法和有限差分方法的求解的等值線的形態(tài)十分相似，并且相對于FD方法，EFG方法的等值線更為光滑。

4.3 實測數(shù)據(jù)計算

圖6給出了一個某地地面上的實測總強度磁異常曲線[14]，該數(shù)據(jù)包含45個測點的場值信息。由于此地區(qū)淺部覆蓋有一層磁性較強并且很不均勻的玄武巖，因此磁場出現(xiàn)劇烈的跳動，在該曲線上有兩個局部異常。為了消除地區(qū)玄武巖的干擾，對該數(shù)據(jù)進行延拓，采用EFG方法進行數(shù)值求解。將求解區(qū)域Ω擴展到無窮遠邊界，下邊界的場值采用測量的邊界場值按余弦衰減的方式衰減至零[6]。

圖5 EFG方法和FD方法磁場空間等值線Fig.5 Magnetic field space contour by EFG and FD method in model 2(a) EFG方法磁場空間等值線；(b) FD方法磁場空間等值線

圖6 實測數(shù)據(jù)曲線Fig.6 Curve of the measured data

圖7 求解區(qū)域Fig.7 Computation domain

圖8 磁場空間等值線Fig.8 Magnetic field space contour

圖8給出了采用EFG方法得到的磁場空間等值線。由圖8可見,通過位場向上延拓,基本消除了淺部玄武巖這一磁性體的干擾,突出了深部磁性體所產(chǎn)生的磁場。

為了分析在不同高度的磁異常的衰減關(guān)系,本文主要選取了三個不同的高度進行分析,分別是向上延拓1個測點距，5個測點距和20個測點距，對不同高度的磁異常強度進行曲線對比,如圖9所示。由圖9可見，隨著高度的增加,異常體磁場衰減很快, 異常體產(chǎn)生磁場的異常也就越不明顯。當向上延拓20個測點距時，已分辨不出上述的兩個局部異常了,這與實際情況相符。EFG方法的數(shù)值模擬結(jié)果對研究區(qū)域場有重要作用,這也是向上延拓的必要性所在。

圖9 實測總強度磁異常進行向上延拓結(jié)果Fig.9 The upward continuation result of the observed data

5 結(jié)論

采用EFG方法，對不同模型和實例的二維位場延拓問題進行了數(shù)值求解，同時與有限差分(FD)法的數(shù)值結(jié)果進行了比較。數(shù)值結(jié)果表明：

1)EFG方法在數(shù)值求解過程中，只需輸入節(jié)點信息和邊界條件,減少了有限差分法中復雜的網(wǎng)格生成過程，相對于FD方法，采用EFG方法均能很好地處理二維位場延拓問題，驗證了無網(wǎng)格方法在求解位場延拓問題中的可行性。

2)采用EFG方法進行實例求解，數(shù)值模擬結(jié)果與實際形態(tài)趨勢基本一致，實現(xiàn)了位場的向上延拓,表明了方法求解位場延拓問題的有效性。

3)無網(wǎng)格方法需要通過背景網(wǎng)格進行高斯積分，并且由于邊界條件的特殊處理等問題,使得無網(wǎng)格法目前的計算速度尚低于FD法，因此EFG法在位場延拓問題中的應用還需要更深入的研究。

[1] 鄧正棟,關(guān)洪軍,聶永平,等.穩(wěn)定地電場三維有限差分正演模擬[J].石油物探，2001,40(1):107-114. DENG ZH D,GUAN H J, NIE Y P，et al. 3D finite difference forward modeling of stable geoeleclrical field[J].Geophysical Prospecting for Petroleum,2001,40(1):107-114. (In Chinese)

[2] 尤淼，魯霞，張健.有限差分法在二維位場延拓的應用[J].工程地球物理學報，2010,7(3):269-273. YOU M, LU X, ZHANG J. Application of FDM to two-dimensional potential field extension[J]. Chinese Journal of Engineering Geophysics, 2010,7(3):269-273. (In Chinese)

[3] 楊曦，潘和平．井間電磁場時域有限差分數(shù)值模擬[J]．地球物理學進展，2008，23(2)：573-682． YANG X, PAN H P. The simulation of cross-hole electromagnetic fields using FDTD method[J]. Progress in geophysucs, 2008, 23(2): 573-682．(In Chinese)

[4] 胡博，岳建華，鄧帥奇．邊界元算法在電法勘探正演中的應用綜述[J]．地球物理學進展，2010，25(3):1024-1030. HU B, YUE J H, DENG SH Q. Review on the forward modelling by the boundary element method in electric exploration[J]. Progress in Geophys, 2010，25(3):1024-1030. (In Chinese)

[5] 郭榮文.磁場的有限元模擬與延拓研究[D].長沙：中南大學，2007. GUO R W. Finite element simulation of magnetic field and extension research[D].Changsha: Central South University, 2007. (In Chinese)

[6] 宋滔,羅勇,王宇航.有限單元法在二維位場延拓中的應用[J].工程地球物理學報, 2012,9(2):134-138. SONG T, LUO Y, WANG Y H. Application of FEM to two-dimensional potential field extension[J]. Chinese Journal of Engineering Geophysics, 2012,9(2):134-138. (In Chinese)

[7] 徐世浙,樓云菊. 起伏地形二維位場上延與換算的邊界單元法[J]. 物探化探計算技術(shù),1987,9(3): 195-199. XU SH Z, LOU Y J. Application of the boundary element method to the continuation upward and transformation of 2-D potential field on relief topography[J]. Computing techniques for Geophysical and Geochemical Exploration, 1987,9(3): 195-199. (In Chinese)

[8] 徐世浙,王慶乙，王軍.用邊界單元法模擬二維地形對大地電磁場的影響[J].地球物理學報，1992,35(3):380-388. XU SH Z, WANG Q Y, WANG J. Modeling 2-D terrain effect on MT by the boundary element method [J]. Chinese Journal of Geophysics, 1992,35(3):380-388.(In Chinese)

[9] 徐世浙．位場延拓的積分-迭代法[J].地球物理學報，2006，49(4)：1176-1182. XU SH Z. The integral-iteration method for continuation of potential fields[J]. Chinese Journal of Geophysics, 2006, 49(4):1176-1182.(In Chinese)

[10]BELYTSCHKO T, LU Y Y, GU L. Element-free Galerkin methods [J]. International journal for numerical methods in fluids, 1994, 37(2): 229-256.

[11]趙光明，宋順成. 無網(wǎng)格Galerkin法與有限元耦合新算法[J].應用數(shù)學與力學，2005，26(8):899-904. ZHAO G M，SONG SH C. New Algorithm of Coupling Element-Free Galerkin and Finite Element Method[J]. Applied Mathematics and Mechanics, 2005,26(8):899-904. (In Chinese)

[12]張小華, 歐陽潔. 線性定常對流占優(yōu)對流擴散問題的無網(wǎng)格解法 [J]. 力學季刊, 2006, 27(2): 220-226. ZHANG X H, OUYANG J. The element free Galerkin method for steady convection dominated convection- diffusion problem [J]. Chinese Quarterly of Mechanics, 2006, 27(2): 220-226. (In Chinese)

[13]王衛(wèi)東，趙國群，欒貽國. 無網(wǎng)格方法中本質(zhì)邊界條件的處理[J]. 力學季刊，2002，23(4)：521-527. WANG W D, ZHAO G Q, LUAN Y G. Treatment of essential boundary conditions for element-free Galerkin method[J]. Chinese Quarterly of Mechanics,2002,23(4):521-527. (In Chinese)

[14]譚承澤，郭紹雍．磁法勘探教程[M ]．北京:地質(zhì)出版社,1984． TAN CH Z，GUO SH Y. Magnetic prospecting guide[M].Beijing: Geological Publishing House,1984. (In Chinese)

第39卷 第2期2017年3月物探化探計算技術(shù)COMPUTING TECHNIQUES FOR GEOPHYSICAL AND GEOCHEMICAL EXPLORATIONVol.39 No.2Mar. 2017

Research on element free Galerkin method for 2-D potential field extension

KONG Qian1, LI Peng1, LI Jing2

(1.School of Applied Mathematics and Physics, North China Electric Power University, Zhuozhou 071003, China; 2.BGP,CNPC, Baoding 072751, China)

The upward continuation of potential field can be induced to determine the function for satisfying the Laplace equation. The Element Free Galerkin (EFG) method is extended and employed for numerical computations of 2-D potential field extension. The principle of the EFG method and its numerical implementation are described. The discrete equation of EFG method is built for 2-D potential field extension. The results of our method are compared with that of finite difference (FD) method. Simulation results show that the proposed method and FD method have the same performances. In addition, a practical case is given in this paper. The theoretic and real numerical results show that EFG method is efficient and easy-implementation in solving the problem of potential field extension.

potential field extension; the element free Galerkin method; the finite difference method; forward calculation

2015-12-15 改回日期：2016-02-15

國家自然科學基金項目(11274111);河北省自然科學基金(F2015502014)

孔倩(1981-), 女，博士，主要研究偏微分方程數(shù)值解法，E-mail:qiankongkong @126.com。

1001-1749(2017)02-0155-06

P 631.3

A

10.3969/j.issn.1001-1749.2017.02.01