李俊杰，嚴家斌

(1.浙江省水利水電勘測設(shè)計院，浙江杭州310002;2.中南大學地球科學與信息物理學院，湖南長沙410083)

正演模擬是研究大地電磁的基礎(chǔ)手段。當?shù)叵码娦越Y(jié)構(gòu)非一維時,大地電磁場響應沒有解析表達式，必須借助正演數(shù)值模擬方法。作為大地電磁正演計算的常用網(wǎng)格方法，有限差分法(FDM)[1]、積分方程法(IEM)[2]及有限元法(FEM)[3-5]均有著各自的優(yōu)缺點：有限差分法實現(xiàn)過程直接，但無法處理復雜地球物理模型;積分方程法只須對異常體進行剖分和求積，不涉及微分方法中的吸收邊界等復雜問題，在三維電磁數(shù)值模擬研究中具有快速、方便等特點，但同樣無法應對地下電阻率很復雜時的計算;有限元法適用于復雜物性分布和復雜邊界形狀的計算，其最大缺陷在于求解復雜模型時網(wǎng)格生成困難。

無網(wǎng)格法作為網(wǎng)格方法的重要補充和發(fā)展,是近十多年來興起的一類數(shù)值計算新方法。其中，無單元Galerkin法(element-free Galerkin method，EFGM)[6]作為一種基于節(jié)點的無網(wǎng)格方法，已被成功用于地震波場[7-10]、雷達波場[11-12]及大地電磁場[13-15]響應的計算。相關(guān)文獻的研究成果表明，EFGM具有較常規(guī)網(wǎng)格方法精度高、計算復雜模型便利的特性。然而，EFGM形函數(shù)不滿足克羅內(nèi)克δ函數(shù)性質(zhì)(Ni(X)=δij)，邊界條件加載較困難。無網(wǎng)格點插值法(MPIM)[16]是一種簡單高效的無網(wǎng)格方法，較EFGM最大的改進在于形函數(shù)采用插值方法構(gòu)造，邊界條件可精確加載。此方法在彈性力學領(lǐng)域取得了很好的應用效果[16]，但在地球物理學領(lǐng)域的應用研究尚未見報導。

本文將無網(wǎng)格點插值法(MPIM)應用于大地電磁二維正演，介紹該方法的近似原理，推導大地電磁二維變分問題的無網(wǎng)格化離散過程;通過多個二維理論模型的MPIM，EFGM和FEM正演計算結(jié)果的對比，分析MPIM大地電磁二維正演數(shù)值模擬方法的特性和應用效果。

1 大地電磁二維變分問題

當?shù)叵码娦越Y(jié)構(gòu)為二維時，取走向為z軸，x軸與z軸垂直，y軸垂直向上；求解域Ω為矩形區(qū)域，4個頂點依次以A,B,C,D順時針編號,Γ1為地質(zhì)體的邊界(圖1)。

當平面電磁波以任何角度入射地面時，地下介質(zhì)中的電磁波總是以平面波形式幾乎垂直地向下傳播，滿足[17]:

圖1 大地電磁二維正演求解區(qū)域a 電場極化模式(TE模式); b 磁場極化模式(TM模式)

(2)

對于TM模式:

(3)

式中:ω為角頻率;μ為磁導率;σ為電導率;ε為介電常數(shù)。

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

MPIM利用位于積分點支持域內(nèi)的場節(jié)點構(gòu)造形函數(shù)，需用一組背景網(wǎng)格將求解區(qū)域離散以便建立離散系統(tǒng)方程，如圖2所示。由于背景網(wǎng)格的積分常選用高斯積分法，故積分點又稱高斯積分點或高斯點。

2.1 支持域

常用的支持域形狀有圓形與矩形兩種(圖2)，

圖2 無網(wǎng)格點插值法的背景網(wǎng)格、支持域、積分點與場節(jié)點

對于任一高斯點XQ，其支持域尺寸d由下式確定:

(4)

式中:α為支持域的無量綱尺寸，它用于控制實際支持域的大小，是對無網(wǎng)格法計算精度影響很大的一個參數(shù)[16]。dc為位于高斯點XQ附近的平均結(jié)點間距，可由下式確定：

(5)

式中:A為預估的支持域面積;n為包含在A中的節(jié)點數(shù)。對于節(jié)點均勻分布的情況，dc為節(jié)點間距。本文采用矩形支持域，故有兩個方向的支持域尺寸，即

(6)

式中:dcx與dcy分別為橫、縱向節(jié)點間距;αx與αy為對應的支持域無量綱尺寸。為便于程序設(shè)計，研究中常取αx=αy=α，本文取α=1.0。

2.2 點插值法近似

求解域Ω中任意一點X處的場變量u(X)的近似表達式為

(7)

式中:pT(X)為二維空間坐標XT=[xy]的基函數(shù);a是系數(shù)向量;m為單項式的個數(shù)。p(X)可用Pascal三角形確定，對于一維(1D)和二維(2D)空間，其線性基函數(shù)分別為

pT(X)=1xm=2p=1

(1D)

pT(X)=[1xy]m=3p=1

(2D)

二次基函數(shù)分別為

pT(X)=[1xx2]m=3p=2

(1D)

pT(X)=[1xyx2xyy2]

m=6p=2

(2D)

完備的p階多項式一般可由下式表示

基函數(shù)階次的增加不能顯著提高無網(wǎng)格法的近似精度，反而降低了計算效率[8]，一般取線性基即可。將(7)式寫成矩陣的形式:

(8)

式中:U=[u1u2…un]T;a=[a1a2…an]T;P的展開式為

(9)

傳統(tǒng)PIM中局部支持域內(nèi)的節(jié)點數(shù)總是等于基函數(shù)個數(shù)(m=n),因此(9)式是一方陣的形式，系數(shù)向量a可通過矩陣的逆運算得到,即

(10)

將(10)式代回(7)式得

(11)

式中，ΦT(X)即為計算點X在支持域內(nèi)的PIM形函數(shù),其表達式為

ΦT(X)=pT(X)P-1=

[φ1(X)φ2(X) …φn(X)]

(12)

因為采用多項式基函數(shù)，PIM形函數(shù)的求導運算較容易，有

(13)

采用多項式作為基函數(shù)時，PIM插值形式簡單，并且對于規(guī)則分布的場節(jié)點的插值具有較高的精度。此外PIM形函數(shù)滿足克羅內(nèi)克δ函數(shù)性質(zhì)，因此邊界條件加載較容易。

然而(9)式的求解是以P非奇異為前提的，當多項式基選用不合理或場節(jié)點分布不恰當時，P可能呈現(xiàn)高度病態(tài)甚至呈現(xiàn)奇異性，雖然存在隨機移點法[16]、局部坐標變換法[18]等處理P奇異性的方法，但這些方法均或多或少存在不足之處，限制了多項式PIM的應用。

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

定義了無網(wǎng)格形函數(shù)，便可以構(gòu)造(1)式的離散系統(tǒng)方程，將變分符號代入(1)式有

(14)

u=[φ1φ2…φn][u1u2…un]T

(15)

式中:Φ為形函數(shù)矩陣;n為支持域內(nèi)的節(jié)點數(shù);u為支持域內(nèi)n個節(jié)點的場向量。對(15)式進行變分運算有

(16)

將(16)式和(15)式代入(14)式，得

(17)

(17)式采用支持域內(nèi)n個場節(jié)點編號。對于求解域內(nèi)所有場節(jié)點,還應有一套總體編號體系，用于將局部節(jié)點矩陣組裝成總體剛度矩陣。當節(jié)點I和節(jié)點J位于不同支持域時，其相應的被積函數(shù)為0，此時可得到按求解域節(jié)點編號的總體編號體系方程

(18)

(18)式可簡寫為

(19)

由于δUT是任意的，故(19)式成立的條件為

(20)

(20)式即為無網(wǎng)格點插值法構(gòu)造的系統(tǒng)方程。

2.4 求解域背景網(wǎng)格積分

(20)式中的K包含對求解域Ω與求解域邊界Γ的積分，為計算這些積分，須將求解域離散成一組背景網(wǎng)格(圖2)，總體積分可表示成這些單元積分之和的形式，每個單元的積分利用高斯積分法求解,即

(21)

背景網(wǎng)格積分是MPIM中最重要的數(shù)值計算問題之一，總積分點數(shù)一般取支持域內(nèi)場節(jié)點數(shù)量的2～8倍，但高斯點數(shù)目的增加并不能顯著提高無網(wǎng)格法的計算精度，反而極大地降低了計算效率[15]，因此本文每個背景單元僅采用4(2×2)個高斯點，每個邊界單元采用2個高斯點。求解線性方程組KU=0還需加載邊界條件，可采用與FEM相同的罰函數(shù)法[16]加載。先將剛度矩陣中相應的對角元素KII變?yōu)棣罧II(α為懲罰系數(shù)，其值可取104～1010),然后用邊界上的uI值取代方程組右端的零向量即可。

3 正演計算及效果分析

3.1 計算精度

圖3 模型一

設(shè)計一個3層介質(zhì)模型(模型一，圖3)，驗證MPIM計算大地電磁場響應的高精度特性。該模型第1層電阻率ρ1=500Ω·m,層厚h1=1km;第2層電阻率ρ2=2000Ω·m,層厚h2=3km;第3層電阻率ρ3=100Ω·m,層厚h3=4km。場節(jié)點等間距分布于問題域，對于TE模式，使用3321(41×81)個場節(jié)點，3200(40×80)個背景網(wǎng)格，橫、縱向節(jié)點間距均為200m;對于TM模式，使用1681(41×41)個場節(jié)點，1600(40×40)個背景網(wǎng)格，節(jié)點間距與TE模式相同。分別采用EFGM，MPIM和FEM進行正演計算與對比，有限元法節(jié)點的分布與無網(wǎng)格法一致。

表1和表2分別為EFGM，MPIM和FEM計算的視電阻率及視相位數(shù)值相對誤差。相對誤差值的大小對應著數(shù)值方法計算精度的高低，相對誤差值越大，精度越低，反之則精度越高。由表1可見,3種數(shù)值方法計算精度相當，且精度均隨頻率的增高而降低，精度變化范圍約為3×10-10～4×10-3;當頻率f<10-3Hz時，MPIM精度約高出EFGM一倍，大于此頻點時MPIM精度與EFGM相當。視相位數(shù)值計算精度變化規(guī)律與視電阻率類似，精度約為2×10-10～4×10-4(表2)。

3.2 有效性

為了進一步說明MPIM求解MT問題的有效性，設(shè)計圖4所示的二維模型(模型二):圍巖電阻率ρ1=1000Ω·m，方形二度異常體電阻率ρ2=100Ω·m，方形截面邊長L=400m，異常體頂部到地面的距離h=800m，場節(jié)點的分布與計算3層介質(zhì)模型時相同(圖4a)。分別采用EFGM，MPIM和FEM進行正演計算與對比，圖5為頻率f=100Hz時模型二的數(shù)值計算結(jié)果曲線，由圖5可見，MPIM計算結(jié)果與EFGM和FEM基本一致，視電阻率曲線極小值與視相位曲線極大值的橫坐標(X=0)正好對應于方形截面中心在地面的投影，證明了二維情況下MPIM的正確性。

圖6為模型二的FEM，EFGM與MPIM視電阻率計算結(jié)果斷面圖，可見3種數(shù)值方法計算結(jié)果仍然一致且對稱性良好。TE模式下，在100Hz頻點附近，MPIM和EFGM與FEM斷面圖存在細微差異;TM模式下，形函數(shù)采用插值法構(gòu)造的數(shù)值方法MPIM與FEM計算結(jié)果一致，較形函數(shù)采用擬合法構(gòu)造的EFGM在10-4～10-3Hz低頻段存在輕微差異。

表1 3層模型視電阻率數(shù)值解的相對誤差

表2 3層模型視相位數(shù)值解的相對誤差

圖4 二維理論模型a 模型二; b 模型三; c 模型四; d 模型五

圖5 頻率f=100Hz時圖模型二的數(shù)值計算結(jié)果a 視電阻率數(shù)值計算結(jié)果; b 視相位數(shù)值計算結(jié)果

3.3 計算效率

MPIM中模型參數(shù)的加載基于高斯積分點而非單元，高斯點的位置可由坐標確定，該特性使其在復雜地質(zhì)模型構(gòu)建上較常規(guī)網(wǎng)格方法方便。為表明MPIM的這一優(yōu)勢，采用與前文相同的節(jié)點分布計算了如圖4b至圖4d所示的二維模型:模型三為圓形截面低阻體，電阻率ρ2=100Ω·m，圍巖電阻率ρ1=1000Ω·m，低阻體截面半徑R=200m，異常體頂部到地面的距離h=800m;模型四與模型五為出露于地表的條帶狀低阻體，異常體截面呈平行四邊形，下底長a=400m，上下底間距d=600m，異常體右邊界與地面的夾角分別呈45°與30°，其電阻率及圍巖電阻率與模型二相同。

圖7為模型三的MPIM和EFGM視電阻率計算結(jié)果,由圖7可見，TE模式異常橫向拉長，TM模式異?？v向拉長，兩種模式下數(shù)值計算結(jié)果均很好地反映了低阻異常體的存在。

圖6 模型二視電阻率數(shù)值計算結(jié)果a TE模式FEM視電阻率; b TM模式FEM視電阻率; c TE模式EFGM視電阻率; d TM模式EFGM視電阻率; e TE模式MPIM視電阻率; f TM模式MPIM視電阻率

圖7 模型三視電阻率數(shù)值計算結(jié)果a TE模式EFGM視電阻率; b TM模式EFGM視電阻率; c TE模式MPIM視電阻率; d TM模式MPIM視電阻率

圖8為模型四與模型五的TE模式MPIM和EFGM計算結(jié)果，橫坐標為里程，縱坐標為以2為底的對數(shù)工作頻率。由圖8可見，①模型四與模型五的視電阻率斷面呈現(xiàn)左右不對稱，低阻區(qū)域呈“上窄下寬”狀分布，模型四視電阻率變化范圍均為100～1060Ω·m，異常特征顯著;②低阻區(qū)域呈傾斜條帶狀分布，傾向與地質(zhì)體的傾向相同;③模型五視電阻率變化范圍均為100～1160Ω·m，其斷面較模型四低阻條帶呈現(xiàn)的傾角更大，且EFGM異常區(qū)斷面較MPIM更為寬泛，如圖8c和圖8d 所示)。

以上模型試算分析可見EFGM和MPIM的高精度特性及其在計算復雜地質(zhì)模型上的優(yōu)越性。雖然無網(wǎng)格法與有限元法均達到精度要求，但在計算耗時上卻存在較大差異,表3列出了3種數(shù)值算法計算模型二的耗時，工作頻率為10-4～103Hz，共17個頻點。由表3可見：①3種數(shù)值算法的計算速度由快到慢依次為FEM，MPIM，EFGM，EFGM計算耗時約為FEM的10倍;②MPIM速度較EFGM快，TM模式下效率提高8%，TE模式下提高5%。

無網(wǎng)格法(MPIM和EFGM)計算的耗時主要花費在形函數(shù)的構(gòu)建及數(shù)值積分上[15]，計算效率的提升是無網(wǎng)格法廣泛應用需解決的首要問題，無網(wǎng)格方法與高效網(wǎng)格方法的耦合或許可成為一個可行的方向。

圖8 TE模式模型四與模型五視電阻率無網(wǎng)格法計算結(jié)果a 模型四EFGM視電阻率計算結(jié)果; b模型四MPIM視電阻率計算結(jié)果; c 模型五EFGM視電阻率計算結(jié)果; d 模型五MPIM視電阻率計算結(jié)果

s

4 結(jié)束語

復雜地質(zhì)模型電磁響應的計算一般需采用非結(jié)構(gòu)化網(wǎng)格，此種網(wǎng)格的生成算法較復雜。本文將MPIM應用于大地電磁二維正演數(shù)值模擬，介紹了MPIM的基本原理及大地電磁二維變分問題的無網(wǎng)格化求解過程;采用節(jié)點規(guī)則分布的無網(wǎng)格法(MPIM和EFGM)計算出了柱形二度體及傾斜條帶狀二度體模型的電磁響應。理論模型的MPIM，EFGM和FEM正演計算結(jié)果的對比分析表明:MPIM適用于大地電磁正演計算，其計算精度很高，數(shù)值算例表明精度最高可達10-9數(shù)量級;MPIM與EFGM精度相當，但MPIM的計算效率更高，本文17個頻點下的數(shù)值計算表明TM模式下效率提高8%，TE模式下提高5%。研究結(jié)果驗證了MPIM大地電磁二維正演數(shù)值模擬方法的有效性，體現(xiàn)了MPIM在處理復雜地質(zhì)模型正演問題上的優(yōu)越性。

參 考 文 獻

[1] 譚捍東，余欽范，Booker J，等.大地電磁法三維交錯采樣有限差分數(shù)值模擬[J].地球物理學報,2003,46(5):705-711

Tan H D,Yu Q F,Booker J,et al.Magnetotelluric three-dimensional modeling using the staggered-grid finite difference method[J].Chinese Journal of Geophysics (in Chinese),2003,46(5):705-711

[2] Wannamaker P E.Advances in three-dimensional magnetotelluric modeling using integral equations[J].Geophysics,1991,56(11):1716-1728

[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)格的三維大地電磁自適應矢量有限元模擬[J].中南大學學報(自然科學版),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] Belytschko T,Lu Y Y,Gu L.Element-free Galerkin methods[J].International Journal for Numerical Methods in Engineering,1994,37(2),229-256

[7] 賈曉峰,胡天躍.滑動最小二乘法求解地震波波動方程[J].地球物理學進展，2005，20(4):920-924

Jia X F,Hu T Y.Solving seismic wave equation by moving least squares (MLS) method[J].Progress in Geophysics,2005,20(4):920-924

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

[9] 賈曉峰,胡天躍,王潤秋.無單元法用于地震波波動方程模擬與成像[J].地球物理學進展，2006，21(1):11-17

Jia X F,Hu T Y,Wang R Q.Wave equation model-

ing and imaging by using element-free method[J].Progress in Geophysics,2006，21(1):11-17

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

[11] 馮德山,王洪華,戴前偉.基于無單元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

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

[13] 蘇洲,胡文寶,朱毅.二維大地電磁正演中無網(wǎng)格算法研究[J].石油天然氣學報,2012,34(5):87-90

Su Z,Hu W B,Zhu Y.Meshfree method used in two-dimensional magnetotelluric forwarding[J].Journal of Oil and Gas Technology,2012,34(5):87-90

[14] 李俊杰,嚴家斌.無網(wǎng)格法進展及其在地球物理學中的應用[J].地球物理學進展,2014,29(1):452-461

Li J J,Yan J B.Developments of meshless method and applications in geophysics[J].Progress in Geophysics(in Chinese),2014,29(1):452-461

[15] 李俊杰.基于弱式無網(wǎng)格法的大地電磁二維正演[D].長沙:中南大學地球科學與信息物理學院,2014

Li J J.Magnetotelluric two-dimensional forward by weak-form meshfree methods[D].Changsha:Central South University of Geosciences and Info-Physics (in Chinese),2014

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

[17] 徐世浙.地球物理中的有限單元法[M].北京:科學出版社，1994:229-235

Xu S Z.Finite element method for geophysics[M].Beijing:Science Press,1994:229-235

[18] Liu G R.Meshfree methods:moving beyong the finite element method[M].USA:CRC Press,2002:123-126