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

(1.浙江省水利水電勘測設(shè)計院,浙江杭州310002;2.中南大學(xué)地球科學(xué)與信息物理學(xué)院,有色資源與地質(zhì)災(zāi)害探查湖南省重點實驗室,湖南長沙410083)

有限元-點插值耦合法大地電磁二維正演模擬

李俊杰1,嚴(yán)家斌2

(1.浙江省水利水電勘測設(shè)計院,浙江杭州310002;2.中南大學(xué)地球科學(xué)與信息物理學(xué)院,有色資源與地質(zhì)災(zāi)害探查湖南省重點實驗室,湖南長沙410083)

點插值法(PIM)作為一種典型的全域弱式無網(wǎng)格法,該方法在地質(zhì)建模時將物性加載到只與坐標(biāo)有關(guān)的高斯積分點上,因此處理復(fù)雜模型時較常規(guī)網(wǎng)格方法便利,但缺點是計算效率低。將有限元法(FEM)與PIM耦合,形成FE-PIM,用于大地電磁二維正演模擬。利用Galerkin法代入插值法構(gòu)造的形函數(shù)并結(jié)合高斯積分公式推導(dǎo)了大地電磁二維無網(wǎng)格化總體矩陣表達(dá)式,簡述了背景網(wǎng)格積分與邊界條件的加載技術(shù),理論模型的數(shù)值計算驗證了FE-PIM算法的正確性、高效性及其在處理復(fù)雜模型上的便利性。

點插值法;全域弱式無網(wǎng)格法;大地電磁;有限元-點插值法

大地電磁二維正演多采用有限差分法(finite deference method,FDM)[1-2]、有限元法(finite element method,FEM)[3-6]等網(wǎng)格方法。對于地質(zhì)體相對規(guī)則、地形簡單的模型,FDM利用較小的網(wǎng)格就能獲得較高的計算精度,效率高;而對于復(fù)雜的地質(zhì)體,就需要精細(xì)網(wǎng)格才能完成高精度的模擬。FEM適用于復(fù)雜邊界與復(fù)雜電性結(jié)構(gòu)模型的計算,但高維問題中非結(jié)構(gòu)化三角單元生成程序?qū)崿F(xiàn)困難。無單元Galerkin法(element-free Galerkin method,EFGM)[7]是一種基于節(jié)點的全域弱式無網(wǎng)格方法,其形函數(shù)構(gòu)造不依賴預(yù)定義的網(wǎng)格單元,選擇離散的節(jié)點便可求Helmholtz方程的解。相對于低階(low-order)有限元分析,EFGM只需少量的節(jié)點就可獲得光滑的解,且場分量的導(dǎo)數(shù)同樣光滑連續(xù),對于大地電磁輔助分量的求解非常有利,可以獲得與主分量(水平電場或磁場)相同的計算精度,對于節(jié)點的參數(shù)化處理也較單元分析容易。EFGM能處理網(wǎng)格方法難以處理的問題[8-10],有廣闊的應(yīng)用前景。李俊杰等[11]綜述了無網(wǎng)格法進(jìn)展并推導(dǎo)了大地電磁二維問題的無網(wǎng)格化系統(tǒng)矩陣表達(dá)式;賈曉峰等[12-13]將EFGM用于疊前地震模擬,并討論了幾種吸收邊界方法在地震模擬中與EFGM的結(jié)合,研究結(jié)果表明，EFGM精度較FDM精度高，且具有較好的穩(wěn)定性;王月英[14]分析了基函數(shù)階次對EFGM用于地震波正演計算精度的影響;馮德山等[15]、戴前偉等[16]將EFGM用于雷達(dá)波場的二維正演模擬,詳細(xì)推導(dǎo)了探地雷達(dá)無單元法波動方程,并提出了適用于EFGM雷達(dá)正演模擬的透射吸收邊界、Sarma吸收邊界的改進(jìn)方法。嚴(yán)家斌等[17]將EFGM成功應(yīng)用于大地電磁二維正演模擬,并分析了支持域無量綱尺寸與高斯點數(shù)量對無網(wǎng)格法正演計算精度與效率的影響。無網(wǎng)格點插值法(point interpolation method,PIM)[18]是另一種簡單高效的全域弱式無網(wǎng)格法,與EFGM的區(qū)別主要在于形函數(shù)采用過點插值的方法進(jìn)行構(gòu)造,邊界條件加載便利。李俊杰等[19]將PIM用于復(fù)雜模型的MT二維正演模擬,凸顯了PIM較常規(guī)網(wǎng)格方法在地質(zhì)建模上的優(yōu)越性。

我們將PIM與高效的FEM相耦合,形成有限元-點插值法(finite element-point interpolation method,FE-PIM),從大地電磁場的二維正演理論出發(fā),采用Galerkin法推導(dǎo)了FE-PIM系統(tǒng)矩陣表達(dá)式,數(shù)值計算結(jié)果驗證了FE-PIM的有效性及其在計算復(fù)雜模型時高效、便利的特性。

1 大地電磁二維變分問題

研究地下二維電性結(jié)構(gòu),取y軸垂直向上,z軸為走向,x軸向右且與y軸、z軸垂直,求解域Ω為矩形,4個頂點依次以A,B,C,D順時針編號,Γ1為地質(zhì)體邊界。大地電磁二維正演滿足(1)式所示的變分問題[20]:

(1a)

uAB=1

(1b)

δF(u)=0

(1c)

TE模式為:

u=Ez

(2)

TM模式為:

u=Hz

(3)

式中:ω為角頻率;μ為磁導(dǎo)率;σ為電導(dǎo)率;ε為介電常數(shù);Ez,Hz分別為電場及磁場平行于異常體走向的水平分量。

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

PIM以積分點支持域內(nèi)的場節(jié)點利用多項式基插值構(gòu)造形函數(shù),背景網(wǎng)格用于計算總體矩陣表達(dá)式中包含的求解域及求解域邊界的積分項。圖1 為PIM支持域、高斯點、背景網(wǎng)格與場節(jié)點示意圖。由于網(wǎng)格積分多采用高斯積分法,故積分點也稱高斯點。

圖1 全域弱式無網(wǎng)格法支持域、高斯點、背景網(wǎng)格與場節(jié)點圖示

2.1 支持域

如圖1所示,支持域形狀常為矩形或圓形,對于任一高斯點XQ,支持域尺寸d可表示為:

(4)

式中:α為支持域無量綱尺寸,其大小影響無網(wǎng)格法的計算精度[17];dc為支持域內(nèi)高斯點XQ的平均結(jié)點間距,有:

(5)

式中:A為支持域面積;n為包含在A中的節(jié)點總數(shù)。本文采用矩形支持域,x,y方向的支持域尺寸為：

(6)

式中:dcx與dcy分別為x,y方向的節(jié)點間距;αx與αy為對應(yīng)的無量綱尺寸。為便于程序設(shè)計,一般令αx=αy=α,本文取α=1。

2.2 FE-PIM離散系統(tǒng)方程的構(gòu)造

PIM計算效率低但處理復(fù)雜模型便利[19],FEM求解高效,由于PIM與FEM系統(tǒng)矩陣采用相同的方法構(gòu)成,故PIM與FEM具有很好的耦合條件。

如圖2所示,FE-PIM的基本原理是將求解域Ω劃分為實線表示的FEM區(qū)域Ω1及虛線表示的PIM背景單元區(qū)域Ω2,異常體置于Ω2,無窮遠(yuǎn)邊界單元設(shè)為Ω1。將變分符號δ代入(1)式并將FEM區(qū)域單元離散化,有:

(7)

式中:e表示FEM的子單元;CD為求解域底邊界;Γe為邊界單元。

圖2 FE-PIM計算模型二有限單元與背景網(wǎng)格分布圖示

將FEM區(qū)域的積分項即δF1(u)離散,有:

(8)

(9)

式中:U為全部節(jié)點的場值向量;KFEM為FEM區(qū)域計算的總體系統(tǒng)矩陣。求得δF1(u)的離散表達(dá)式后,將δF2(u)離散,將u表示為PIM形函數(shù)向量與場節(jié)點向量之積的形式,有

(10)

式中:Φ為用PIM構(gòu)造的形函數(shù)矩陣[18-19];n為支持域內(nèi)的節(jié)點數(shù);u為支持域內(nèi)n個節(jié)點的場向量值。對(10)式進(jìn)行變分運(yùn)算,有:

(11)

將(10)式與(11)式代入(1)式的第1項并引入總體編號體系,總體矩陣最終變?yōu)?

(12)

式中:M1與Mn表示PIM區(qū)域第一個節(jié)點與最后一個節(jié)點的編號;KPIM為背景單元域PIM計算的總體系統(tǒng)矩陣。KPIM表達(dá)式中包含對求解域Ω2的積分,積分利用高斯積分法求解,有:

(13)

(14)

由于δUT的任意性,(14)式成立的條件變?yōu)?

(15)

(15)式即為用FE-PIM構(gòu)造的系統(tǒng)矩陣,由于FEM與PIM構(gòu)造的總體矩陣KFEM及KPIM均具有帶狀、稀疏、對稱的性質(zhì),故K也具備該特性。

線性方程組KU=0的求解需加載(1)式所示的邊界條件,PIM與FEM邊界條件均可直接加載,計算時先找出邊界節(jié)點在總體剛度矩陣中對應(yīng)的KII,KII對應(yīng)的行與列設(shè)置為0,其值設(shè)置為1,最后將方程右端項對應(yīng)的值PI設(shè)置為邊界上的預(yù)設(shè)場值即可。PIM邊界條件可精確加載,該特性是PIM較EFGM最大的優(yōu)勢。

3 正演計算

3.1 算法驗證

PIM與FEM在一維介質(zhì)大地電磁(MT)正演中具有較高的精度[19],FE-PIM為兩者的耦合,先通過圖3所示的COMMEMI 2D-0模型[21]來驗證文中算法的正確性。此模型異常體頂部與地表重合,縱向規(guī)模50km,橫向規(guī)模20km,電阻率為1Ω·m,其左、右兩側(cè)電阻率分別為10Ω·m和2Ω·m,深度大于50km區(qū)域為超導(dǎo)層,電阻率為0,求解域TE模式取120km×110 km,TM模式取120km×70km,場節(jié)點橫縱向等間距分布,節(jié)點間距為1km。

圖3 COMMEMI 2D-0模型

表1列出了PIM與FEM計算COMMEMI 2D-0模型的視電阻率結(jié)果。從表1可以看出,兩種方法計算結(jié)果與COMMEMI解析解基本吻合,驗證了數(shù)值算法的有效性。此外,表1中FEM計算結(jié)果也與文獻(xiàn)[22]一致。

表1 COMMEMI 2D-0模型數(shù)值方法視電阻率計算結(jié)果

3.2 有效性及優(yōu)勢

為進(jìn)一步驗證FE-PIM有效性及其優(yōu)勢,設(shè)計了如圖4所示的二維地質(zhì)模型:模型二異常體為截面方形二度體,異常體邊長400m,埋深800m;模型三為水平橢圓,長半軸300m,短半軸200m,埋深800m;模型四為等腰直角三角形地壘,斜邊長800m;模型五為出露于地表的條帶狀低阻體,異常體截面呈平行四邊形,下底長400m,上、下底間距800m,異常體右邊界與地面的夾角呈45°。

模型背景電阻率1000Ω·m,異常體電阻率100Ω·m。場節(jié)點等間距分布于求解域,TE模式3321(41×81)個場節(jié)點,TM模式1681(41×41)個場節(jié)點,橫、縱向節(jié)點間距均為200m,FEM與FE-PIM采用相同的節(jié)點分布。

圖5顯示了頻率為100Hz時PIM與FEM模型二視電阻率計算結(jié)果。從圖5可以看出,求解域兩側(cè)邊界的視電阻率值約為1000Ω·m,低阻異常在里程0附近最顯著,TE模式視電阻率值約780Ω·m,TM模式約920Ω·m。兩種模式下FE-PIM,PIM與FEM的計算結(jié)果基本一致,驗證了FE-PIM二維算法的正確性。

圖4 二維地質(zhì)模型

圖5 頻率為100Hz時模型二視電阻率計算結(jié)果

圖6為不同模型的FE-PIM視電阻率計算結(jié)果。從圖6a和圖6b可以看出,TE模式低阻異常中心與模型中心一致,能明顯地反映出橢圓低阻體的存在,但異常規(guī)模比實際低阻體大,視電阻率為740～1080Ω·m。TM模式低阻異常呈直立狀并無限向下延伸,異常的水平寬度反映了實際低阻體水平方向的尺寸,視電阻率為800～1060Ω·m。從圖6c可以看出,模型四山脊位置在視電阻率斷面圖上顯示為高阻極大值區(qū)域,兩側(cè)的山谷區(qū)則呈現(xiàn)低阻極小值,且頻率越高此特性越明顯,視電阻率變化范圍為900～1950Ω·m,此現(xiàn)象表明地形對大地電磁測深法影響規(guī)律與直流電阻率法存在差異。斷面圖高阻區(qū)的橫向規(guī)模約800m,與山脊橫向規(guī)模相當(dāng),驗證了FE-PIM的有效性。從圖6d 可以看出,模型五的視電阻率斷面圖左右非對稱,低阻區(qū)域呈“上窄下寬”條帶狀傾斜分布,傾向與地質(zhì)體相同,視電阻率為100～1100Ω·m,較好地反映出了異常體的物性及空間賦存狀態(tài)。

圖6 不同模型的視電阻率FE-PIM計算結(jié)果

FE-PIM延續(xù)PIM地質(zhì)建模便利特點的同時也兼顧了計算效率,表2列出了17個頻點PIM,FE-PIM與FEM的計算耗時。從表2可以看出,PIM耗時約為FEM耗時的8～9倍,FE-PIM耗時相對PIM耗時呈著降低,無網(wǎng)格區(qū)域(10×4)的FE-PIM與PIM相比，TE模式下效率約提高87%,TM模式下效率約提高83%,隨著FEM區(qū)域的增大,FE-PIM耗時逐漸接近FEM耗時。

表2 模型二不同數(shù)值方法計算耗時

4 結(jié)束語

1) 將FE-PIM應(yīng)用于大地電磁二維正演模擬,通過把求解區(qū)域劃分為有限元區(qū)域及背景單元區(qū)域,利用Galerkin法和高斯積分公式導(dǎo)出了大地電磁二維FE-PIM耦合的系統(tǒng)方程離散表達(dá)式,介紹了邊界條件直接加載的技巧,COMMEMI 2D-0 模型的數(shù)值計算驗證了文中PIM與FEM算法的正確性。

2) 采用節(jié)點規(guī)則分布的FE-PIM計算了截面為方形、橢圓、平行四邊形及三角地壘二維模型的大地電磁場響應(yīng),視電阻率斷面圖較好地反映了地質(zhì)體產(chǎn)生的電磁異常響應(yīng)及空間賦存狀態(tài),驗證了耦合算法的正確性,凸顯了FE-PIM處理復(fù)雜地質(zhì)模型的優(yōu)越性。

3) 比較分析了PIM,FE-PIM及FEM的計算效率,PIM計算耗時明顯高于FEM計算耗時;FE-PIM計算耗時受異常體區(qū)域的影響較大,當(dāng)異常體區(qū)域較小時,FE-PIM計算耗時遠(yuǎn)小于PIM計算耗時而略高于FEM計算耗時,但FE-PIM綜合了FEM計算高效及PIM處理復(fù)雜模型便利的優(yōu)點,是一種優(yōu)越的數(shù)值方法。

[1] 胡祥云,霍光譜,高銳,等.大地電磁各向異性二維模擬及實例分析[J].地球物理學(xué)報,2013,56(12):4268-4277 Hu X Y,Huo G P,Gao R,et al.The magnetotelluric anisotropic two-dimensional simulation and case analysis[J].Chinese Journal of Geophysics,2013,56(12):4268-4277

[2] 董浩,魏文博,葉高峰,等.基于有限差分正演的帶地形三維大地電磁反演方法[J].地球物理學(xué)報,2014,57(3):939-952 Dong H,Wei W B,Ye G F,et al.Study of three-dimensional magnetotelluric inversion including surface topography based on finite-difference method[J].Chinese Journal of Geophysics,2014,57(3):939-952

[3] Key K,Weiss C.Adaptive finite-element modeling using unstructured grids:the 2D magnetotelluric example[J].Geophysics,2006,71(6):291-299

[4] 劉長生,湯井田,任政勇,等.基于非結(jié)構(gòu)化網(wǎng)格的三維大地電磁自適應(yīng)矢量有限元模擬[J].中南大學(xué)學(xué)報(自然科學(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] 柳建新,孫麗影,童孝忠,等.基于自適應(yīng)有限元的起伏地形MT二維正演模擬[J].地球物理學(xué)進(jìn)展,2012,27(5):2016-2023 Liu J X,Sun L Y,Tong X Z,et al.Undulating terrain 2D MT forward modelling with adaptive finite element[J].Progress in Geophysics,2012,27(5):2016-2023

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

[8] Belytschko T,Lu Y Y,Gu L.Fracture and crack growth by element-free Galerkin methods[J].Modelling and Simulation in Materials Science and Engineering,1994,2(3):519-534

[9] Li D M,Bai F N,Cheng Y M,et al.A novel complex variable element-free Galerkin method for two-dimensional large deformation problems[J].Computer Methods in Applied Mechanics and Engineering,2012,233-236(1):1-10

[10] Liu L C,Dong X H,Li C X.Adaptive finite element-element-free Galerkin coupling method for bulk metal forming processes[J].Journal of Zhejiang University-Science A,2009,10(3):353-360

[11] 李俊杰,嚴(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

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

[13] 賈曉峰,胡天躍,王潤秋.無單元法用于地震波波動方程模擬與成像[J].地球物理學(xué)進(jìn)展,2006,21(1):11-17 Jia X F,Hu T Y,Wang R Q.Wave equation modeling and imaging by using element-free method[J].Progress in Geophysics,2006,21(1):11-17

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

[15] 馮德山,王洪華,戴前偉.基于無單元Galerkin法探地雷達(dá)正演模擬[J].地球物理學(xué)報,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

[16] 戴前偉,王洪華.基于隨機(jī)介質(zhì)模型的GPR無單元法正演模擬[J].中國有色金屬學(xué)報,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

[17] 嚴(yán)家斌,李俊杰.無網(wǎng)格法在大地電磁正演計算中的應(yīng)用[J].中南大學(xué)學(xué)報(自然科學(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

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

[19] 李俊杰,嚴(yán)家斌.無網(wǎng)格點插值法大地電磁二維正演數(shù)值模擬[J].石油物探,2014,53(5):617-626 Li J J,Yan J B.Magnetotelluric two-dimensional for-

ward numerical modeling by meshfree point interpolation method[J].Geophysical Prospecting for Petroleum,2014,53(5):617-626

[20] 戴前偉,張富強(qiáng),楊坤平,等.電導(dǎo)率分塊線性變化的二維高頻大地電磁法有限元正演模擬[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 conductivity within each block[J].Computing Techniques for Geophysical and Geochemical Exploration,2012,34(5):552-558

[21] Zhdanov M S,Varentsov I M,Weaver J T,et al.Methods for modelling electromagnetic fields results from COMMEMI-the international project on the comparison of modelling methods for electromagnetic induction[J].Journal of Applied Geophysics,1997,37(3):133-271

[22] 許建榮.起伏地形條件下大地電磁測深二維正反演研究及應(yīng)用[D].長沙:中南大學(xué)地球科學(xué)與信息物理學(xué)院,2010 Xu J R.Research and applications of 2-D MT forward modeling and inversion with topography[D].Changsha:Central South University of Geosciences and Info-Physics,2010

(編輯：顧石慶)

Magnetotelluric two-dimensional forward modelling by finite element-point interpolation coupling method

Li Junjie1,Yan Jiabin2

(1.ZhejiangDesignInstituteofWaterConservancyandHydroelectricPower,Hangzhou310002,China; 2.KeyLaboratoryofNonferrousResourcesandGeologicalHazardDetection,SchoolofGeosciencesandInfo-Physics,CentralSouthUniversity,Changsha410083,China)

Point interpolation method (PIM) is a typical global weak-form meshfree numerical calculation method.The physical properties of PIM are loaded on Gauss integral points which are only related with coordinates during geological modeling.Therefore,PIM is more convenient while dealing complex model than conventional grid method,but the former computation efficiency is low.We couple FEM and PIM into finite element-point interpolation method (FE-PIM) for magnetotelluric two-dimensional forward.Firstly,magnetotelluric two-dimensional discrete system matrix is deduced through substituting the shape function constructed by interpolation method and combining Galerkin method with Gauss integral formula.Then,the principle of background grid integral and the loading technique of boundary conditions are briefly characterized.Finally,the effectiveness and high efficiency of the FE-PIM algorithm and the convenience for complex models are proved by several numerical models.

point interpolation method,global weak-form meshfree method,magnetotelluric,finite element-point interpolation method

2014-12-02;改回日期：2015-02-15。

李俊杰(1989—),碩士,助理工程師,現(xiàn)主要從事大地電磁無網(wǎng)格化正演模擬研究工作。

國家自然科學(xué)基金項目(40874055)和湖南省自然科學(xué)基金項目(14JJ2012)共同資助。

P631

A

1000-1441(2015)04-0477-08

10.3969/j.issn.1000-1441.2015.04.015