馮 凱， 秦 策， 李 論， 郭家松， 蔡盼盼

(1.中國煤炭科工集團(tuán)西安研究院有限公司, 西安 710077；2.河南理工大學(xué) 物理與電子信息學(xué)院，焦作 454000；3.成都理工大學(xué) 地球勘探與信息技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室, 成都 610059)

0 引言

大地電磁測深(MT)是一種利用天然交變電磁場來研究地下電性結(jié)構(gòu)的地球物理勘探方法，由蘇聯(lián)學(xué)者Tikhono[1]和法國學(xué)者Cagniard[2]于20世紀(jì)50年代首次提出。其具有勘探范圍廣、深度大、不受高阻屏蔽，對低阻分辨率高等優(yōu)點(diǎn)?，F(xiàn)已廣泛應(yīng)用到礦產(chǎn)、工程、石油等勘探領(lǐng)域，并取得了不錯的效果。然而在實(shí)際MT工作時，由于野外地質(zhì)結(jié)構(gòu)的復(fù)雜，采集數(shù)據(jù)量龐大，給大地電磁反演解釋工作帶來了很大的諸多困難。正演是反演得基礎(chǔ)，因此，對復(fù)雜的地質(zhì)模型進(jìn)行精度高，速度快正演計算是MT數(shù)值模擬研究的重點(diǎn)。

有限單元法是50年代首先在彈性力學(xué)中發(fā)展起來的方法，主要適用于物性參數(shù)復(fù)雜分布的區(qū)域，但計算量非常大，隨著計算機(jī)的發(fā)展，有限單元法在解決各個工程領(lǐng)域的許多數(shù)學(xué)物理問題中，得到了廣泛地應(yīng)用，成為一種高效通用的計算方法。自Congon[3]于1971年首次將有限元引入大地電磁中，陳小斌，劉云等[4-5]相繼對二三維進(jìn)行模擬。之后許多學(xué)者[6-8]采用自適應(yīng)非結(jié)構(gòu)三角網(wǎng)格有限元方法計算電磁場。非結(jié)構(gòu)三角網(wǎng)格能夠真實(shí)地模擬各種復(fù)雜地質(zhì)結(jié)構(gòu)，該方法利用Delaunay三角劃分原理[9]對目標(biāo)區(qū)域進(jìn)行細(xì)密的網(wǎng)格剖分，而背景區(qū)域則選擇粗糙網(wǎng)格剖分。自適應(yīng)有限元對計算結(jié)果進(jìn)行誤差評估[10]，誤差較大的單元網(wǎng)格自動細(xì)化直至有限元計算結(jié)果滿足給定的誤差范圍。自適應(yīng)非結(jié)構(gòu)有限元方法因其強(qiáng)大的適應(yīng)性，越來越受到地球物理研究者的青睞，其方法技術(shù)也日臻完善[11]。在此基礎(chǔ)上，筆者采用非結(jié)構(gòu)化三角網(wǎng)格對研究區(qū)域進(jìn)行剖分，利用三次插值技術(shù)實(shí)現(xiàn)大地電磁的自適應(yīng)二維正演模擬。通過對國際標(biāo)準(zhǔn)模型COMMEMI-2D1的模擬，分別比較二次插值與本文用到的三次插值自適應(yīng)的網(wǎng)格數(shù)量和數(shù)值模擬結(jié)果，證明了本文算法的精確性與高效性，這對于后續(xù)的反演解釋工作具有借鑒意義。

1 方法原理

1.1 大地電磁邊值問題

大地電磁的二維地電模型中Ex和Hx滿足的偏微分方程為：

▽·(τ▽u)+λu=0

(1)

對于TE極化模式u=Ex，τ=1/iμω，λ=δ-iωε。對于TM極化模式u=Hx，=1/(δ-iωε)，λ=iωε。TE極化模式的上邊界條件離地面足夠遠(yuǎn)，取該處為u=1，而TM極化模式取在地面為u=1。其他邊界條件兩者一致，分別是下邊界滿足，左右邊界滿足。求解此微分方程便可得到大地電磁的地下響應(yīng)。

1.2 伽遼金法的加權(quán)余量表達(dá)式

有限單元法解題的基本思想是將研究區(qū)域剖分為許多小單元，在單元內(nèi)進(jìn)行插值[2]。筆者采用伽遼金加權(quán)余量法將微分方程轉(zhuǎn)化為所對應(yīng)的弱變分形式，取權(quán)重因子ω(x,y)，方程(1)兩邊同時乘以權(quán)重因子并積分可得式(2)。

(2)

公式(2)的左邊第一項(xiàng)，利用矢量用算公式▽·(φA)=▽φ·A+φ▽A和積分的變換公式，并帶入邊界條件可得：

(3)

圖1 三次插值自然坐標(biāo)Fig.1 Cubic interpolation natural coordinates

將式(3)帶入式(2)可得：

(4)

1.3 三次插值有限元分析

設(shè)u是單元中的三次函數(shù)：

u(x,y)=a1x3+a2x3+a3x3+a4x3+

a5x3+a6x3++a7x3

+a8x3+a9x3+a10

(5)

式(5)有10個系數(shù)。除三角的頂點(diǎn)外，在每邊上各取二到三分點(diǎn)，再取三角形的重心c,共10個點(diǎn)(如圖1)。根據(jù)面積坐標(biāo)的定義，易推出邊界上的三分點(diǎn)的面積坐標(biāo)。以jm邊上的三分點(diǎn)p、q為例：

P點(diǎn)：Li=0,Lj=2/3,Lm=1/3;q點(diǎn)：Li=0,Lj=1/3,Lm=2/3。

其他邊上的三分點(diǎn)的面積坐標(biāo)與此類似，三角形重心的面積坐標(biāo)是C點(diǎn)：Li=1/3,Lj=1/3,Lm=1/3。

單元上的三次插值的函數(shù)可表示為:

u=Niui+Njuj+Nmum+Npup+

Nquq+Nrur+Nsus+Nlul+

Nuuu+Ncuc

(6)

其中，Ni、…、Nc是形函數(shù)，它們與面積坐標(biāo)的關(guān)系是：

頂點(diǎn)：Ni=1/2(3Li-1)(3Li-2)Li

邊界點(diǎn)：Np=9/2LjLm(3Lj-1);Nq=9/2LjLm(3Lm-1)

中心點(diǎn)：Nc=27LiLjLm

顯然，這些形函數(shù)是x、y的三次函數(shù)，所以式(6)中的u也是形x、y的三次函數(shù)。

1.4 單元分析和總體矩陣的合成

i=1,2,…，n

(7)

其中基函數(shù)滿足：

(8)

由式(8)可得式子(4)第一項(xiàng)的積分值，同理可算出第二項(xiàng)和第三項(xiàng)的積分。最后單元的積分可表示為：

(9)

通常只需要存儲非零元素，本文采用定帶寬存儲方式，利用LU分解法求解系數(shù)矩陣組成的線性方程組，最終得到全區(qū)域上的節(jié)點(diǎn)場值u，至此有限元的求解全部結(jié)束。

1.5 h-型的自適應(yīng)算法

自適應(yīng)有限元可以使網(wǎng)格自動剖分，在滿足所設(shè)要求的情況下完成計算。而網(wǎng)格劃分有①h-型；②p-型；③h-p型三種形式。h-型劃分是保證各單元內(nèi)的形函數(shù)階數(shù)不變，并逐步加密網(wǎng)格，使結(jié)果逼近精確解；p-型劃分是保證各單元的幾何尺寸不變的同時，逐步增加各單元的階數(shù)，從而逼近精確解；h-p型劃分是在h-型基礎(chǔ)上再進(jìn)行p-型劃分，效果較好但編程困難，很難實(shí)現(xiàn)。筆者采取h-型的自適應(yīng)有限元，先對給定的模型劃分相對較粗糙的網(wǎng)格；然后通過DWR法計算局部誤差，并根據(jù)誤差來指向性的加密指定區(qū)域的網(wǎng)格；根據(jù)加密的網(wǎng)格生成新的網(wǎng)格；再次計算局部網(wǎng)格的誤差直到誤差小于所給定的誤差閾值。最終所生成的網(wǎng)格便是我們所需要的網(wǎng)格剖分，在此網(wǎng)格下進(jìn)行的模擬即能滿足精度又能提高計算速度。

2 算例

2.1 COMMEMI-2D1自適應(yīng)有限元模擬

COMMEMI-2D1模型(圖2)是一個眾多學(xué)者參與的國際合作項(xiàng)目，不同的研究學(xué)者采用不同的模擬方法對此模型進(jìn)行計算，并與發(fā)布的參考結(jié)果作對比來驗(yàn)證算法的有效性。為了方便與參考解作對比，我們采用10 Hz頻率進(jìn)行數(shù)值模擬，與Zhdanov et al[12]發(fā)布的結(jié)果進(jìn)行對比來證明本文算法的有效性。

首先對COMMEMI-2D1模型進(jìn)行粗網(wǎng)格剖分(圖3)，其中研究區(qū)域?yàn)閳D3中的粉紅色部分，剖分區(qū)域長為40 km, 寬為8 km，電阻率為100 Ω·m；深藍(lán)色區(qū)域?yàn)榭諝鈹U(kuò)展層；紅色區(qū)域?yàn)樽筮吔绾拖逻吔绲臄U(kuò)展層；淺藍(lán)色區(qū)域?yàn)?.5 Ω·m的低阻矩形塊，矩形塊位于地表以下250 m處，長為1 km，寬為2 km。引用國際流行的開源工具Triangle對模型進(jìn)行非結(jié)構(gòu)化網(wǎng)格的剖分，初始網(wǎng)格剖分為

圖2 COMMEMI-2D1模型Fig.2 COMMEMI-2D1 model

圖3 COMMWMI-2D1模型粗網(wǎng)格Fig.3 COMMWMI-2D1 Model coarse grid

圖4 10 HZ下的二次插值和三次插值自適應(yīng)網(wǎng)格Fig.4 10Hz Quadratic interpolation and cubic interpolation adaptive grid(a)二次插值;(b)三次插值

圖5 二次插值和三次插值模型的自適應(yīng)有限元MT響應(yīng)結(jié)果Fig.5 Adaptive finite element MT response results of quadratic interpolation and cubic interpolation results(a)TE極化模式的視電阻率響應(yīng)；(b)TM極化模式的視電阻率響應(yīng)；(b)TE極化模式的阻抗相位響應(yīng)；(d)TM極化模式的阻抗相位響應(yīng)

79個網(wǎng)格，分別采用二次插值和三次插值進(jìn)行電磁場模擬。

地表設(shè)置13個測點(diǎn)，對異常體矩形塊采用1%的誤差控制進(jìn)行自適應(yīng)計算?？梢钥闯鲈谕日`差控制的范圍內(nèi)最終計算得到的自適應(yīng)網(wǎng)格，異常體內(nèi)和電性分界面處被自動加密，三次插值的網(wǎng)格更加的稀疏(圖4)，節(jié)點(diǎn)數(shù)也遠(yuǎn)遠(yuǎn)小于二次插值的自適應(yīng)網(wǎng)格。

分別計算MT的響應(yīng)結(jié)果(圖5)與已發(fā)布的結(jié)果相比較證明本算法的有效性和可行性。TE模式和TM模式的視電阻率模擬結(jié)果與參考解非?？拷瑥亩C明本算法的正確性。由模型響應(yīng)結(jié)果可得當(dāng)?shù)叵麓嬖诘妥璁惓ｓw時，TE和TM的視電阻率和阻抗相位都有很好的響應(yīng)，但是TE模式的阻抗相位在異常體周圍顯示低值，而TM模式正好相反。

2.2 三次插值有限元的精度驗(yàn)證

為了比較三次差值的精度，在此模型的基礎(chǔ)上利用自適應(yīng)三次插值最終迭代后的網(wǎng)格(圖4b)做為比較精度的網(wǎng)格。這套網(wǎng)格的優(yōu)點(diǎn)是對于三次插值在這套網(wǎng)格上有和很好的結(jié)果，而二次插值在這套網(wǎng)格上會產(chǎn)生較大的誤差，模擬結(jié)果如圖6所示。

圖6 二次插值和三次插值模型有限元MT響應(yīng)結(jié)果Fig.6 MT response results of quadratic interpolation and cubic interpolation results(a)TE極化模式；(b)TM極化模式

由圖6可以看出，二次插值由圖中黑線表示在遠(yuǎn)離低阻矩形塊時計算的視電阻率，較三次插值由圖中綠線表示和參考接由圖中紅線表示具有較大的誤差，由于在地表以下電性分界面處網(wǎng)格比較稀疏，二次插值的單個網(wǎng)格節(jié)點(diǎn)數(shù)少于三次插值，在進(jìn)行有限元數(shù)值計算過程中所得的節(jié)點(diǎn)場值和輔助場值較真值誤差較大，所以導(dǎo)致二次插值在目標(biāo)區(qū)域以外處表現(xiàn)為波浪狀，而三次插值很好得減小了這種誤差，結(jié)果與參考解接近。

3 結(jié)論

通過對同一套粗網(wǎng)格分別利用二次插值和三次插值進(jìn)行自適應(yīng)計算，得出在同一后驗(yàn)誤差的控制下三次插值的自適應(yīng)在經(jīng)過9次的迭代過程中網(wǎng)格數(shù)量，遠(yuǎn)遠(yuǎn)少于二次插值迭代25次的自適應(yīng)網(wǎng)格。同時通過對同一套網(wǎng)格分別進(jìn)行二次插值和三次插值的有限元計算，將TE模式和TM模式的視電阻率與Zhdanov等[7]對COMMEMI-2D1模型的模擬結(jié)果進(jìn)行比較，得到三次插值的精度也是要高于二次插值的算法。驗(yàn)證了三次插值數(shù)值模擬方法計算速度快，有很高的計算精度，對今后反演解釋工作具有借鑒意義。