李俊杰，嚴家斌，皇祥宇

(1.浙江省水利水電勘測設計院，浙江杭州 310002； 2.中南大學 地球科學與信息物理學院有色資源與地質災害探查湖南省重點實驗室，湖南長沙 410083)

?

無單元Galerkin法大地電磁三維正演模擬

李俊杰1，嚴家斌2，皇祥宇2

(1.浙江省水利水電勘測設計院，浙江杭州 310002； 2.中南大學 地球科學與信息物理學院有色資源與地質災害探查湖南省重點實驗室，湖南長沙 410083)

無單元Galerkin法(EFGM)作為一種相對成熟的無網格方法，避免了網格剖分，其精度高，適用于復雜電導率分布和復雜邊界形狀的計算。本文將EFGM用于大地電磁三維正演，詳述了三維EFGM形函數(shù)的構造過程，從大地電磁三維變分問題出發(fā)，利用Galerkin法結合高斯積分公式推導了相應的系統(tǒng)矩陣離散表達式，簡述了邊界條件的加載技術，研究了支持域尺寸對EFGM三維正演計算精度的影響，最后通過數(shù)值計算驗證了EFGM三維算法的正確性。

無網格法 無單元Galerkin法 大地電磁 支持域尺寸

Li Jun-jie, Yan Jia-bin, Huang Xiang-yu. Three-dimensional forward modeling of magnetotellurics using the element-free Galerkin method[J]. Geology and Exploration, 2015, 51(5):0946-0952.

作為網格法數(shù)值計算的重要補充和發(fā)展，無網格法是近十多年來興起的一類數(shù)值計算新方法(Garbonetal.,2002；徐凱軍等,2006；沈金松，2003；陳輝等，2011；劉長生等, 2010；Tongetal.,2009；Yangetal.,2013；Sunetal.,2014)，其形函數(shù)的構造避免了網格的生成過程, 物性加載在只與坐標位置有關的高斯點上，因此節(jié)點規(guī)則分布下也能計算復雜模型,現(xiàn)已成為力學領域的研究熱點。無單元Galerkin法(Element-Free Galerkin Method，EFGM)是一種較成熟的擬合型無網格方法，其形函數(shù)的構造過程不存在奇異性問題,且能通過權函數(shù)的適當選取實現(xiàn)不同階次的函數(shù)近似。p型自適應過程實現(xiàn)便利，不但廣泛應用于有限元法所觸及的領域，且成功解決了某些網格方法難以處理的問題(Belytschkoetal., 1994；Lietal., 2012；Liuetal., 2009),在地震波場(賈曉峰等,2006)、雷達波場(馮德山等,2013)、大地電磁場(李俊杰等, 2014a,2014b, 2014c, 2015；嚴家斌等，2014；Wittkeetal.,2014)及直流電法(李俊杰等,2015)二維響應的計算上也取得了成功。EFGM具有精度高、處理復雜模型便利等優(yōu)勢,然而有關EFGM在大地電磁三維正演中的應用還未見報導。

本文嘗試將EFGM應用于MT三維正演，介紹了三維移動最小二乘(MLS)近似原理，推導了MT三維變分問題的EFGM總體矩陣表達式，簡述了含背景網格的三維高斯積分技術，最后通過數(shù)值計算驗證了EFGM算法的正確性。

1 大地電磁三維變分問題

當?shù)叵码娦越Y構為三維時，取向東為x軸，x軸與y軸垂直，z軸垂直向下,求解域V為六面體區(qū)域，上下表面的4個頂點分別以A、B、C、D與E、F、G、H順時針編號(圖1)。

大地電磁三維正演滿足式(1)的變分問題(Tongetal., 2009)：

(1)

圖1 大地電磁三維模型Fig. 1 Three-dimensional MT model

2 大地電磁三維變分問題的無網格解法

2.1 支持域

EFGM利用位于高斯點支持域內的場節(jié)點構造形函數(shù)，支持域為高斯點附近人為劃定的一個區(qū)域，其概念相當于FEM的單元,常用的支持域形狀有球體與長方體兩種。對于任一高斯點XQ，其支持域尺寸d由式(2)確定：

d=αdc

(2)

式中：α為支持域的無量綱尺寸，是對無網格法計算精度與效率影響很大的重要參數(shù)；dc為位于高斯點XQ附近的平均結點間距，可由式(3)確定：

(3)

式中：B為預估的支持域體積；n為包含在B中的節(jié)點數(shù)。對于節(jié)點均勻分布的情況，dc為節(jié)點間距。本文采用六面體支持域，故有三個方向的支持域尺寸，即式(4)：

(4)

式中：dcx、dcy與dcz分別為直角坐標系三個方的的節(jié)點間距；αx、αy與αz為對應的支持域無量綱尺寸。為便于程序設計，研究中常取αx=αy=αz=α，本文中取α=1.0,即支持域的大小與FEM單元相等。

2.2 移動最小二乘近似(moving least-square，MLS)

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

(5)

式中：p(X)為三維空間坐標XT=[x,y,z]的基函數(shù)，可用Pascal三角形確定(趙熙強等, 2015)，其線性基函數(shù)為pT(X)={1xyz}；m為單項式的個數(shù)；a(X)是與高斯點位置相關的系數(shù)向量。MLS近似中的節(jié)點數(shù)通常大于基函數(shù)個數(shù)(n>m),系數(shù)向量a(X)需通過如式(6)的加權殘差法求得：

(6)

式中W(X-Xi)為權函數(shù)，常用無網格權函數(shù)有三次樣條函數(shù)、四次樣條函數(shù)及指數(shù)函數(shù)(Belytschkoetal., 1994)。本文選擇四次樣條函數(shù)為權函數(shù)，表達式如式(7)：

(7)

式中：Ri=di/rω；di=|X-Xi|。di為節(jié)點Xi與高斯點X之間的距離，rω為權函數(shù)支持域尺寸。通過式(6)的運算將產生如式(8)的線性關系：

A(X)a(X)=B(X)Us

(8)

式中：Us={u1u2…un}T為支持域內所有節(jié)點的節(jié)點場函數(shù)向量，A(X)與B(X)表達式如下：

(9)

B(X)=[W1(X)p(X1)W2(X)p(X2)…Wn(X)p(Xn)]

(10)

式(9)中Wi(X)=W(X-Xi),求解式(8)可得：

a(X)=A-1(X)B(X)Us

(11)

將上式代入式(5)得：

φiui=ΦT(X)Us

(12)

ΦT(X)即為計算點X在支持域內的MLS形函數(shù)，其表達式為：

2.3 無網格離散系統(tǒng)方程的構造

將計算點處E的三分量Ex、Ey、Ez表示為形函數(shù)與場節(jié)點之積的形式，有：

(14)

式中：Φ為MLS近似構造的形函數(shù)矩陣；n為支持域內的節(jié)點數(shù)；Ex、Ey與Ez為支持域內n個節(jié)點的場向量。將式(14)代入式(1)有：

(15)

利用式(15)將式(1)最終表示為式(16)的形式：

δEKE=δE(K1+K2+K3)E=0

(16)

式(16)即為EFGM三維系統(tǒng)方程。

K的表達式中包含對求解域V與求解域邊界面Γ的積分，為計算這些積分，一般將求解域離散成一組背景網格，以K1與K2矩陣左上角的元素為例,其積分可表示成單元積分之和的形式，有

(17)

求解線性方程組KE=0還需加載邊界條件。EFGM邊界條件可用罰函數(shù)法加載，將剛度矩陣中相應的對角元素KII變?yōu)棣罧II，α為罰系數(shù)，其值可取104～1010,然后將邊界上的E值取代方程組右端的零向量即可。

3 數(shù)值計算

為了驗證EFGM算法的正確性，本文計算了若干數(shù)值模型：模型一為ρ1=1000Ω·m的均勻介質模型,用于支持域無量綱尺寸最優(yōu)參數(shù)值的試驗, 采用(11×11×66)個場節(jié)點，(10×10×65)個背影網格，節(jié)點間距為100m,空氣層厚度為500m。模型二(圖2)為二層介質模型，第一層電阻率ρ1=1000Ω·m，層厚h1=3km；第二層電阻率ρ2=10Ω·m，層厚h2=3km,空氣層厚度及節(jié)點分布情況與模型一相同；模型三(圖3)為六面體低阻模型，背景電阻率為1000Ω·m,異常體電阻率為10Ω·m,長寬高均為600m, 埋深800m，異常體中軸線的投影對應于求解域的原點，采用(41×41×46)個場節(jié)點，節(jié)點間距為200m?？諝鈱雍穸葹?000m。

圖2 模型二Fig. 2 model

圖3 模型三Fig. 3 model 3

圖4為不同支持域無量綱尺寸下模型一的EFGM計算結果。如圖4所示,當支持域無量綱尺寸α=1～1.2時，計算結果與解析解基本一致,隨著α的增大數(shù)值解開始偏離解析解，且α越大計算精度越低。數(shù)值試驗過程中還發(fā)現(xiàn)當α>1.2時求解系統(tǒng)方程組的迭代次數(shù)會隨著頻率的增高而增大，當α>1.7時相同參數(shù)設置下計算結果開始發(fā)散。α增大的同時伴隨著計算效率的降低，因此本文數(shù)值計算取α=1。

圖5為二層介質模型的EFGM數(shù)值計算結果。如圖5所示，EFGM計算結果與解析解基本一致, 驗證了算法的正確性。圖6 為頻率f=10Hz時三維低阻模型的數(shù)值計算結果，視電阻率與阻抗相位斷面圖在中部區(qū)域呈現(xiàn)極小值，極小值的兩側附近呈現(xiàn)極大值,異常形態(tài)呈直立的橢圓形,背景區(qū)域的視電阻率接近1000Ω·m，阻抗相位接近45°(圖6)，較好地反映出了地下異常體的存在，進一步凸顯了EFGM在大地電磁三維正演中的有效性。

圖4 不同支持域無量綱尺寸模型一的EFGM數(shù)值計算結果Fig. 4 Calculation results of EFGM for model 1 with different support domain dimensionless sizesa-視電阻率; b-阻抗相位a-apparent resistivity; b-impedance phase

圖5 二層介質模型的EFGM數(shù)值計算結果Fig. 5 Calculation results of EFGM for two-layer medium modela-阻抗相位;b-視電阻率a-impedance phase;b-apparent resistivity

然而,EFGM模擬三維電磁場響應計算效率較低,模型三采用的剖面方式在64位PC電腦Matlab平臺上運行耗時約2h,內存占用約450MB，計算耗時主要來源于形函數(shù)構造過程中高斯點支持域內節(jié)點的搜索及含背景網格的高斯積分,因此研究高效的積分方法是EFGM三維模擬進入實用的重要手段。

圖6 頻率為10Hz時模型三EFGM數(shù)值計算結果Fig. 6 Calculation results of model 3 by EFGM when frequency is 10Hz

4 結論

本文詳述了三維移動最小二乘(MLS)近似原理，推導了MT三維變分問題的無網格離散系統(tǒng)方程，簡述了三維背景網格高斯積分與邊界條件的罰函數(shù)法加載技術。

EFGM求解MT三維變分問題支持域無量綱尺寸最優(yōu)區(qū)間為α=1～1.2,計算精度隨α的增大而降低，α的增大同時增加了計算耗時，建議MT正演中取α=1.0。

層狀介質模型視電阻率與阻抗相位計算結果均與解析解一致，10Hz頻點地表電磁響應特征也較好地反映出了三維異常體的存在，驗證了算法的正確性。

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

Chen Hui, Deng Ju-zhi, Tan Han-dong. 2011. Study on divergence correction method in three-dimensional magnetotelluric modeling with staggered-grid finite difference method[J]. Chinese Journal of Geophysics, 54(6): 1649-1659(in Chinese with English abstract)

Feng De-shan, Wang Hong-hua, Dai Qian-wei. 2013. Forward simulation of Ground Penetrating Radar based on the element-free Galerkin method[J]. Chinese Journal of Geophysics, 56(1): 298-308(in Chinese with English abstract)

Garbon H, Zhdanov M S. 2002. Comtraction intergral equation method in three-dimensional electromagnetic modeling[J]. Radio Science, 37(6): 1089-1101

Jia Xiao-feng , Hu Tian-yue, Wang Run-qiu. 2006 . A mesh-free algorithm of seismic wave simulation in complex medium[J]. Oil Geophysical Prospecting, 41(1): 43-48(in Chinese with English abstract)

Li Dong-ming, Bai Fu-nong, Cheng Yu-min. 2012. A novel complex variable element-free Galerkin method for two-dimensional large deformation problems[J]. Computer Methods in Applied Mechanics and Engineering, 233-236(1): 1-10.

Li Jun-jie, Yan Jia-bin. 2014. Developments of meshless method and applications in geophysics[J]. Progress in Geophysics, 29(1): 452-461(in Chinese with English abstract)

Li Jun-jie, Yan Jia-bin. 2014. Magnetotelluric two-dimensional forward numerical modeling by meshfree point interpolation method[J]. Geophysical Prospecting for Petroleum, 53(5): 617-626(in Chinese with English abstract)

Li Jun-jie, Yan Jia-bin. 2014. Magnetotelluric two-dimensional forward modeling using meshless local Petrov-Galerkin method[J]. Coal Geology & Exploration, 42(6): 101-104(in Chinese with English abstract)

Li Jun-jie, Yan Jia-bin. 2015. A mesh-free method for the variational problem of the point source two-dimensional electric field[J]. Geophysical and Geochemical Exploration, 39(3): 606-609(in Chinese with English abstract)

Li Jun-jie, Yan Jia-bin. 2015. Magnetotelluric two-dimensional forward by finite element-radial point interpolation method[J]. The Chinese Journal of Nonferrous Metals, 25(5): 1314-1324(in Chinese with English abstract)

Liu Chang-sheng, Tang Jing-tian, Ren Zheng-yong. 2010. Three-dimension magnetotellurics modeling by adaptive edge finite-element using unstructured meshes[J]. Journal of Central South University (Science and Technology), 41(5): 1855-1859(in Chinese with English abstract)

Liu Lei-chao, Dong Xiang-huai, Li Cong-xin. 2009. Adaptive finite element-element-free Galerkin coupling method for bulk metal forming processes[J]. Journal of Zhejiang University-Science A, 10(3): 353-360

Shen Jin-son. 2003. Modeling of 3D electromagnetic responses in frequency domain by using staggered grid finite difference method[J]. Chinese Journal of Geophysics, 46(2): 281-288(in Chinese with English abstract)

Sun Wei-feng, Song Yan, Gong Yan-jie, Gui Li-li. 2014. Characteristics and distribution prediction of structural fissures in the Lower Cretaceous Xiagou Formation in the Qingxi oilfield[J]. Geology and Exploration, 50(6): 1181-1189(in Chinese with English abstract)

Tong Xiao-zhong, Liu Jian-xin, Xie Wei. 2009. Three-dimensional forward modeling for magnetotelluric sounding by finite element method[J]. Journal of Central South University of Technology, 16(1): 136-142

Wittke J, Tezkan B. 2014. Meshfree magnetotelluric modeling[J]. Geophysical Journal International, 198(2): 1255-1268.

Xu Kai-jun, Li Tong-lin, Zhang Hui. 2006. Three-dimensional magnetotelluric forward modeling using integral equation[J]. Northwestern Seismological Journal, 28(2): 104-107(in Chinese with English abstract)

Yan Jia-bin, Li Jun-jie. 2014. Magnetotelluric forward calculation by meshless method[J]. Journal of Central South University (Science and Technology), 45(10): 3513-3520(in Chinese with English abstract)

Yang Li-gong, Liu Ji-shun, Yin Li-jun, Liu Wei-ming, Liu Wen-heng, Guo Jun. 2013. Application of the dual-frequency IP method to rapid prospecting of mineral resources in the Lizi area of Gansu Province[J]. Geology and Exploration, 49(2): 330-336(in Chinese with English abstract)

Zhao Xi-qiang, Li Lin. 2015. Further generalization of Pascal functional matrix and application[J]. Periodical of Ocean University of China, 45(3): 136-140(in Chinese with English abstract)

[附中文參考文獻]

陳 輝, 鄧居智, 譚捍東. 2011. 大地電磁三維交錯網格有限差分數(shù)值模擬中的散度校正方法研究[J]. 地球物理學報, 54(6): 1649-1659

馮德山,王洪華,戴前偉. 2013. 基于無單元Galerkin法探地雷達正演模擬[J].地球物理學報, 56(1): 298-308

賈曉峰, 胡天躍, 王潤秋. 2006 . 復雜介質中地震波模擬的無單元法[J]. 石油地球物理勘探, 41(1): 43-48

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

李俊杰, 嚴家斌. 2014.無網格點插值法大地電磁二維正演數(shù)值模擬[J].石油物探, 53(5): 617-626

李俊杰, 嚴家斌. 2015. 大地電磁二維正演中的有限元-徑向基點插值法[J]. 中國有色金屬學報, 25(5): 1314-1324

李俊杰, 嚴家斌. 2015. 點源二維電場變分問題的無網格解法[J]. 物探與化探, 39(3): 606-609

李俊杰,嚴家斌. 2014. 無網格局部Petrov-Galerkin法大地電磁二維正演[J]. 煤田地質與勘探， 42(6)：101-104

劉長生, 湯井田, 任政勇. 2010. 基于非結構化網格的三維大地電磁自適應矢量有限元模擬[J]. 中南大學學報(自然科學版), 41(5): 1855-1859

沈金松. 2003. 用交錯網格有限差分法計算三維頻率域電磁響應[J]. 地球物理學報, 46(2): 281-288

孫維鳳, 宋 巖, 公言杰, 桂麗黎. 2014. 青西油田下溝組構造裂縫發(fā)育特征與分布預測[J]. 地質與勘探, 50(6): 1181-1189

徐凱軍, 李桐林, 張 輝, 2006. 利用積分方程法的大地電磁三維正演[J]. 西北地震學報, 28(2): 104-107

嚴家斌, 李俊杰. 2014. 無網格法在大地電磁正演計算中的應用[J]. 中南大學學報(自然科學版), 45(10): 3513-3520

楊立功, 劉繼順, 尹利君, 劉衛(wèi)明, 劉文恒, 郭 軍. 2013. 雙頻激電偶極在甘肅李子地區(qū)快速找礦中的應用[J]. 地質與勘探, 49(2): 330-336

趙熙強, 李 琳. 2015. Pascal函數(shù)矩陣的進一步推廣及應用[J]. 中國海洋大學學報, 45(3): 136-140

Three-Dimensional Forward Modeling of Magnetotellurics Using the Element-Free Galerkin Method

LI Jun-jie1, YAN Jia-bin2, HUANG Xiang-yu2

(1.ZhejiangDesignInstituteofWaterConservancyandHydroelectricPower,Hangzhou,Zhejiang310002; 2.KeyLaboratoryofNon-ferrousResourcesandGeologicalHazardDetection,SchoolofGeosciencesandInfo-Physics,CentralSouthUniversity,Changsha,Hunan410083)

The element-free Galerkin method (EFGM) is a relatively mature high-precision mesh-free method which avoids the mesh generation and has the advantages in dealing with complex conductivity distribution and geometrical boundaries. This work applies EFGM to three-dimensional forward modeling of magnetotellurism (MT). We introduce the shape function constructing process of three-dimensional EFGM and then deduce the discrete system matrix expression corresponding to the three-dimensional variational problem of MT by combining the Galerkin method and the Gauss integral formula. Then, the approach of assigning boundary conditions is briefly described. Finally, the effect of support domain size on the solution accuracy is studied and the correctness of three-dimensional EFGM is verified by numerical calculations of several models.

mesh-free method, element-free Galerkin method, magnetotelluric, support domain size

2015-01-13；

2015-06-16；[責任編輯]郝情情。

國家自然基金項目(40874055)和湖南省自然基金資助項目(14JJ2012)聯(lián)合資助。

李俊杰(1989年-)，男，2014年畢業(yè)于中南大學，獲碩士學位，助理工程師，從事大地電磁無網格化正演研究。E-mail：lijunjiecsu@163.com。

[文獻標識碼]A [文章編號]0495-5331(2015)05-0946-07