蘇 洲，胡文寶，朱 毅 （油氣資源與勘探技術(shù)教育部重點實驗室（長江大學(xué)），湖北荊州434023）

二維大地電磁正演中無網(wǎng)格算法研究

蘇 洲，胡文寶，朱 毅 （油氣資源與勘探技術(shù)教育部重點實驗室（長江大學(xué)），湖北荊州434023）

無網(wǎng)格法是近幾年來發(fā)展的一種新的基于變分原理的數(shù)值計算方法，由于在計算形函數(shù)中不需要劃分網(wǎng)格，在力學(xué)、電磁學(xué)等領(lǐng)域得到了廣泛的研究?；跓o網(wǎng)格法在大地電磁勘探正演中的應(yīng)用進行了研究，首先對無網(wǎng)格法的基本原理進行了闡述，并利用廣義變分原理推導(dǎo)出了相應(yīng)的離散方程，編制了相應(yīng)的程序。最后通過理論模型的計算結(jié)果檢驗了該算法的正確性。

無網(wǎng)格法；大地電磁法；移動最小二乘擬合（MLS）；有限元法

大地電磁測深法（MT）是一種以巖石的電性差異為基礎(chǔ)和前提，利用天然交變電磁場研究地球電性結(jié)構(gòu)的勘探方法［1～3］。在解決MT正演問題時，有限差分法［4］和有限元法［5，6］是主要的數(shù)值計算方法；盡管這些方法取得了巨大成功，但是這些方法都是基于網(wǎng)格的數(shù)值計算方法，也有其自身的局限性［7］，如網(wǎng)格部分計算成本高、計算精度依賴于單元剖分的形狀和大小等。產(chǎn)生上述問題的一個主要原因是這些方法都是基于網(wǎng)格剖分來建立離散系統(tǒng)方程的，為克服這些困難，無網(wǎng)格方法（Meshless）作為有限元法（FEM）的一種重要補充具有重要的研究意義。

無網(wǎng)格法（Meshless）是近10年來興起的一種數(shù)值計算方法［7］，它是在有限元法的基礎(chǔ)上發(fā)展起來的，但是不同于有限元；無網(wǎng)格法的近似函數(shù)是建立在一系列離散節(jié)點上的，然后在這些節(jié)點上利用形函數(shù)求出其近似值，代入相應(yīng)的變分問題得到系統(tǒng)方程。

無網(wǎng)格法已在力學(xué)［8］、油藏滲流問題［9］、電磁學(xué)［10］等領(lǐng)域得到了廣泛的研究，然而在求解地球物理正演問題中的理論和應(yīng)用研究鮮少出現(xiàn)。為此，筆者利用無網(wǎng)格法對二維MT正演做了相關(guān)研究。首先從麥克斯韋方程組出發(fā)，推導(dǎo)出二維介質(zhì)中邊值所對應(yīng)的廣義變分問題；然后用節(jié)點離散求解區(qū)域，得到建立在離散節(jié)點上的無網(wǎng)格形函數(shù)，求解變分問題得到邊值問題對應(yīng)的離散方程，解方程后代入視電阻率計算公式計算地面上待求點的視電阻率。

圖1 理論模型

1 邊值問題及其變分問題

大地電磁場可以看作是從高空垂直入射到地面的平面電磁波［1］，基于這種假設(shè)可以推導(dǎo)二維情況下大地電磁測深法所對應(yīng)的邊值問題［2］，如圖1所示：

式（1）所對應(yīng)的變分問題為求如下泛函I（u）的極值：

2 二維大地電磁正演中無網(wǎng)格算法研究

微分方程及其邊界條件與變分問題是等價的，因而求解微分方程可以等價為求解變分問題［11］。其基本思想為：求解變分問題的最終目的是得到在求解區(qū)域中N個離散節(jié)點的值，假設(shè)待求函數(shù)在N個節(jié)點有準(zhǔn)確值，任意點的函數(shù)值可以通過形函數(shù)來表示，將其代入變分問題中，求泛函極小即可得到關(guān)于N個節(jié)點值的方程。因而在求解變分問題時，首先要構(gòu)造形函數(shù)。

不同的方法構(gòu)造的形函數(shù)不同［12］：有限元法中，其形函數(shù)是通過單元的一組固定節(jié)點利用插值技術(shù)構(gòu)造的；而在無網(wǎng)格法中，任意點的場變量是由該點局部支持域中的一組場節(jié)點近似表示的。該次研究采用移動最小二乘近似法來構(gòu)造無網(wǎng)格形函數(shù)。設(shè)求解域Ω內(nèi)，待求場函數(shù)u在N個節(jié)點處的場值ui＝u（xi）（i＝1，2，…，N）是已知的，基于這些已知節(jié)點構(gòu)造待求函數(shù)u在求解域Ω上的近似uh（x），設(shè)表達式為：

式中，x為計算點為x的鄰域中的結(jié)點；a（x）為待求系數(shù)為基函數(shù)；I表示基函數(shù)的個數(shù)。在求解域Ω內(nèi)，N個節(jié)點處定義權(quán)函數(shù)。筆者分別對三次和四次樣條權(quán)函數(shù)進行研究。

設(shè)計算點x的定義域Ωx包括N個節(jié)點，近似函數(shù)在節(jié)點＝xi處的加權(quán)離散平方和為：J＝

本質(zhì)邊界條件在式（2）中并未得到滿足，因而采用廣義變分原理［12］將場函數(shù)滿足的本質(zhì)邊界條件引入泛函，使有約束條件的變分問題變成無約束條件的變分問題。構(gòu)造修正泛函常用的方法有拉格朗日乘子法和罰函數(shù)法［12］。該次研究采用拉格朗日乘子法得到修正泛函I＊（u，珔λ）：

在有限元中，求解域Ω被離散成一系列單元，積分問題可以轉(zhuǎn)化為對各單元積分的和。在無網(wǎng)格法中，求解域Ω是用節(jié)點離散表示的，不存在網(wǎng)格。因此在無網(wǎng)格法中需要特殊的積分方案。用規(guī)則網(wǎng)格覆蓋求解域Ω，將對求解域Ω的積分轉(zhuǎn)化為對各規(guī)則單元的積分之和，在每個格子中使用高斯積分［12］，如圖2所示。通過上述方法得到式（6），求解得各節(jié)點處的值，其值為各節(jié)點的HY，EY值，代入視電阻率計算公式［2］即可。

3 數(shù)值模擬計算

對于上述求解過程，利用Matlab編寫了二維MT正演的無網(wǎng)格算法程序，并與有限元方法做了對比研究，根據(jù)圖3所示模型進行了計算。結(jié)果見圖4～7。

圖2 無網(wǎng)格法中背景網(wǎng)格單元積分點及其節(jié)點圓形影響域

圖3 理論模型

圖4 TE、TM極化模式下0點電阻率對比圖

圖5 TE極化極化模式下0點相位對比圖

由圖4、5可見，無網(wǎng)格法和有限元法有相似的計算結(jié)果，證明了二維情況下算法的正確性；由圖6、7可見，二維等值線可以很好地反映異常體特征。此外，由于無網(wǎng)格法形函數(shù)構(gòu)造不需要網(wǎng)格信息，因而對于求解區(qū)域節(jié)點的設(shè)置沒有限制，相對于有限元法和有限差分法簡單了許多；在計算過程中發(fā)現(xiàn)權(quán)函數(shù)和節(jié)點影響域?qū)τ嬎憔扔兄苯拥挠绊懀辉谟嬎阈史矫?，由于無網(wǎng)格法需要在每個積分點重新計算其鄰域中所包含的節(jié)點數(shù)，需要更多的節(jié)點信息，因而需要更長的計算時間。現(xiàn)有的資料中已經(jīng)出現(xiàn)過無網(wǎng)格并行算法［13］，為了將無網(wǎng)格法用于實際，研究地球物理無網(wǎng)格算法的并行計算將是今后研究的方向之一。

圖6 TM極化模式下視電阻率和相位等值線圖

圖7 TE極化模式下視電阻率和相位等值線圖

4 結(jié) 論

1）無網(wǎng)格數(shù)據(jù)結(jié)構(gòu)簡單，只需要節(jié)點信息，脫離了網(wǎng)格的限制，計算結(jié)果具有高階連續(xù)性。

2）在構(gòu)造無網(wǎng)格形函數(shù)時，權(quán)函數(shù)的選擇和節(jié)點影響域?qū)τ嬎憬Y(jié)果有著較大的影響，對于具體模型計算前，要根據(jù)數(shù)值試驗確定上述參數(shù)，筆者采用4次樣條權(quán)函數(shù)，影響域權(quán)重取2.3。

3）由于無網(wǎng)格法在每個計算點處要重新計算相應(yīng)的形函數(shù)值，因而其計算復(fù)雜度較高。

［1］樸化榮.電磁測深法原理［M］.北京：地質(zhì)出版社，1990.1～50.

［2］徐世浙.地球物理中的有限單元法［M］.北京：科學(xué)出版社，1994.50～164.

［3］胡建德.電法勘探中的數(shù)學(xué)模型［J］.數(shù)學(xué)的實踐與認識，2004，34（2）：26～30.

［4］譚捍東，余欽范，John Booker，等.大地電磁法三維交錯采樣有限差分?jǐn)?shù)值模擬［J］.地球物理學(xué)報，2003，46（5）：705～710.

［5］陳樂壽.有限元法在大地電磁測深正演計算中的應(yīng)用與改進［J］.石油物探，1981，20（3）：84～103.

［6］陳小斌，張翔，胡文寶.有限元直接迭代算法在MT二維正演計算中的應(yīng)用［J］.石油地球物理勘探，2000，35（4）：487～496.

［7］陶文銓，吳學(xué)紅，戴艷俊.無網(wǎng)格數(shù)值求解方法［J］.中國機電工程學(xué)報，2010，30（5）：1～10.

［8］張雄，劉巖.無網(wǎng)格法［M］.北京：清華大學(xué)出版社，2004.14～103.

［9］李玉坤，姚軍，黃朝琴，等.油藏滲流問題的無網(wǎng)格法分析［J］.中國石油大學(xué)學(xué)報（自然科學(xué)版），2003，31（2）：95～104.

［10］茅昕光，林鶴云.電磁場數(shù)值分析的無單元Galerkin方法［J］.東南大學(xué)學(xué)報，2003，33（1）：31～35.

［11］老大中.變分法基礎(chǔ)［M］.北京：國防工業(yè)出版社，2007.11～104.

［12］Liu G R，Gu Y T.An introduction to meshfree methods and their programming［M］.Berlin：Springer，2005.54～114.

［13］曾清紅.無網(wǎng)格數(shù)值模擬的并行算法及并行實現(xiàn)研究［D］.合肥：中國科學(xué)技術(shù)大學(xué)，2006.

［編輯］ 龍 舟

87 Meshfree Method Used in Two－dimensional Magnetotelluric Forwarding

SU Zhou，HU Wen－bao，ZHU Yi

（AuthorsAddress：Key Laboratory of Exploration Technologies for Oil and Gas Resources（Yangtze University），Ministry of Education，Jingzhou434023，Hubei，China）

As a new numerical computational method based on the principle of variation，meshfree methods has maintained a rapid development in recently years.It has widely applied in studying the problem of the mechanics and electromagnetics because of avoiding the onerous mesh generation.The application of meshfree method in magnetotelluric forwarding was studied.First，the fundamental principle of meshfree method was described and the particular meshless calculating formulas were deduced through the generalized variation principle.The program is verified by comparison with the analytic response of a symmetrical ladders model and with finite difference results of laterally inhomogeneous model.

meshfree method；magnetotelluric；moving least－square method；FEM

book=112,ebook=112

P631.325

A

1000－9752（2012）05－0087－04

2012－02－10

國家“973”規(guī)劃項目（2007CB209607）；國家自然科學(xué)基金項目（40727001）。

蘇洲（1986－），男，2009年大學(xué)畢業(yè)，碩士生，現(xiàn)主要從事電磁測深正演方法的研究工作。