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

孔 倩, 李 鵬, 李 晶

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

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

孔 倩1, 李 鵬1, 李 晶2

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

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

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

0 引言

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

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

1 基礎(chǔ)介紹

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

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

(1)

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

(2)

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

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

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

(3)

(4)

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

使式(4)中的J在駐點(diǎn)x達(dá)到最小，則可得到待定系數(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 模型求解

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

(7)

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

KT=q

(8)

其中：

(9)

(10)

4 方法應(yīng)用

4.1 理論模型1

有一水平地面測(cè)線，包含15個(gè)測(cè)點(diǎn)的場(chǎng)值信息。其中，地面邊界上2號(hào)點(diǎn)～14號(hào)點(diǎn)磁場(chǎng)100nT，其余兩點(diǎn)磁場(chǎng)為0nT(圖1)。

圖1 模型1中測(cè)點(diǎn)的場(chǎng)值信息
Fig.1 Field value in model 1

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

(11)

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

(12)

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

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

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

4.2 理論模型2

有一水平地面測(cè)線，同樣15個(gè)測(cè)點(diǎn)信息。地面邊界上2號(hào)點(diǎn)～9號(hào)點(diǎn)磁場(chǎng)為100nT，10號(hào)點(diǎn)～14號(hào)點(diǎn)磁場(chǎng)為10nT,其余兩點(diǎn)磁場(chǎng)為0nT,如圖4所示。

圖4 模型2中測(cè)點(diǎn)的場(chǎng)值信息
Fig.4 Field value in model 2

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

4.3 實(shí)測(cè)數(shù)據(jù)計(jì)算

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

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

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

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

圖8 磁場(chǎng)空間等值線Fig.8 Magnetic field space contour

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

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

圖9 實(shí)測(cè)總強(qiáng)度磁異常進(jìn)行向上延拓結(jié)果Fig.9 The upward continuation result of the observed data

5 結(jié)論

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

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

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

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

[1] 鄧正棟,關(guān)洪軍,聶永平,等.穩(wěn)定地電場(chǎng)三維有限差分正演模擬[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] 尤淼，魯霞，張健.有限差分法在二維位場(chǎng)延拓的應(yīng)用[J].工程地球物理學(xué)報(bào)，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] 楊曦，潘和平．井間電磁場(chǎng)時(shí)域有限差分?jǐn)?shù)值模擬[J]．地球物理學(xué)進(jìn)展，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] 胡博，岳建華，鄧帥奇．邊界元算法在電法勘探正演中的應(yīng)用綜述[J]．地球物理學(xué)進(jìn)展，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] 郭榮文.磁場(chǎng)的有限元模擬與延拓研究[D].長(zhǎng)沙：中南大學(xué)，2007. GUO R W. Finite element simulation of magnetic field and extension research[D].Changsha: Central South University, 2007. (In Chinese)

[6] 宋滔,羅勇,王宇航.有限單元法在二維位場(chǎng)延拓中的應(yīng)用[J].工程地球物理學(xué)報(bào), 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] 徐世浙,樓云菊. 起伏地形二維位場(chǎng)上延與換算的邊界單元法[J]. 物探化探計(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] 徐世浙,王慶乙，王軍.用邊界單元法模擬二維地形對(duì)大地電磁場(chǎng)的影響[J].地球物理學(xué)報(bào)，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] 徐世浙．位場(chǎng)延拓的積分-迭代法[J].地球物理學(xué)報(bào)，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ú)網(wǎng)格Galerkin法與有限元耦合新算法[J].應(yīng)用數(shù)學(xué)與力學(xué)，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áng)潔. 線性定常對(duì)流占優(yōu)對(duì)流擴(kuò)散問(wèn)題的無(wú)網(wǎng)格解法 [J]. 力學(xué)季刊, 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)東，趙國(guó)群，欒貽國(guó). 無(wú)網(wǎng)格方法中本質(zhì)邊界條件的處理[J]. 力學(xué)季刊，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月物探化探計(jì)算技術(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

國(guó)家自然科學(xué)基金項(xiàng)目(11274111);河北省自然科學(xué)基金(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