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

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

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

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

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

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

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

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

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

圖1 理論模型

1 邊值問(wèn)題及其變分問(wèn)題

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

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

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

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

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

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

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

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

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

3 數(shù)值模擬計(jì)算

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

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

圖3 理論模型

圖4 TE、TM極化模式下0點(diǎn)電阻率對(duì)比圖

圖5 TE極化極化模式下0點(diǎn)相位對(duì)比圖

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

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

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

4 結(jié) 論

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

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

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

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

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

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

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

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

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

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

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

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

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

［11］老大中.變分法基礎(chǔ)［M］.北京：國(guó)防工業(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ú)網(wǎng)格數(shù)值模擬的并行算法及并行實(shí)現(xiàn)研究［D］.合肥：中國(guó)科學(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

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

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