趙曉博，朱自強(qiáng)，李建慧，彭凌星

(中南大學(xué) 地球科學(xué)與信息物理學(xué)院，湖南長(zhǎng)沙 410083)

基于非結(jié)構(gòu)化網(wǎng)格的瞬變電磁2.5維有限元正演模擬

趙曉博，朱自強(qiáng)，李建慧，彭凌星

(中南大學(xué) 地球科學(xué)與信息物理學(xué)院，湖南長(zhǎng)沙 410083)

利用Delaunay三角化這種網(wǎng)格非結(jié)構(gòu)化方法，通過編程實(shí)現(xiàn)了二維模型的非結(jié)構(gòu)化三角形網(wǎng)格剖分，并編寫了中心回線法瞬變電磁2.5維有限元正演程序。與前人計(jì)算結(jié)果對(duì)比，在取得相同計(jì)算精度的情況下，與結(jié)構(gòu)化網(wǎng)格相比，非結(jié)構(gòu)化網(wǎng)格所需網(wǎng)格和節(jié)點(diǎn)數(shù)量大大減少，計(jì)算效率更高。通過將非結(jié)構(gòu)化網(wǎng)格法引入到瞬變電磁2.5維正演模擬中，實(shí)現(xiàn)了對(duì)復(fù)雜二維地電模型的有限元數(shù)值模擬，提高了現(xiàn)有有限元算法的應(yīng)用范圍。

非結(jié)構(gòu)化網(wǎng)格;Delaunay三角化;瞬變電磁;有限單元法

0 前言

瞬變電磁法是一種時(shí)間域電磁勘探方法，已廣泛應(yīng)用于能源、礦產(chǎn)、水文、工程、環(huán)境等領(lǐng)域。目前，瞬變電磁法資料解釋水平較低，基本上停留在一維階段，而二維或三維瞬變電磁反演在解釋技術(shù)方面仍處于探索研究階段，離真正實(shí)用階段仍有相當(dāng)大的距離［1、2］。由于正演是反演、資料解釋的基礎(chǔ)，所以開展瞬變電磁法的正演研究很有必要。

2.5 維正演模擬是在實(shí)際中應(yīng)用較多的一種正演模型，與三維模型相比，2.5維正演模擬只需對(duì)截面而不是整個(gè)體積作離散處理，節(jié)省了計(jì)算量。同時(shí)，與一維、二維模型相比，2.5維正演模擬又能較好的近似地質(zhì)情況。

目前，國(guó)內(nèi)、外對(duì)于人工源電磁場(chǎng)的2.5維有限元數(shù)值模擬已做了大量研究，并取得了豐碩的成果，但基本上是基于結(jié)構(gòu)化網(wǎng)格實(shí)現(xiàn)的［3～10］。結(jié)構(gòu)化網(wǎng)格在人工源電磁場(chǎng)問題的數(shù)值模擬中，有下面兩個(gè)缺點(diǎn):

(1)結(jié)構(gòu)化網(wǎng)格的幾何適應(yīng)能力通常較差，因此對(duì)于較復(fù)雜計(jì)算區(qū)域，很難控制網(wǎng)格節(jié)點(diǎn)的分布情況，并會(huì)產(chǎn)生較大的離散誤差，從而影響有限元計(jì)算精度。

(2)結(jié)構(gòu)化網(wǎng)格很難處理網(wǎng)格的疏密過渡。而在電磁場(chǎng)的有限元正演模擬中，要做一些微分運(yùn)算，需要對(duì)局部做網(wǎng)格加密處理，這時(shí)使用結(jié)構(gòu)化網(wǎng)格就必須將網(wǎng)格劃分得足夠小，從而使節(jié)點(diǎn)數(shù)大大增加，增加了不必要的計(jì)算量，影響了有限元計(jì)算速度［11、12］。

而相對(duì)于結(jié)構(gòu)化網(wǎng)格，非結(jié)構(gòu)化網(wǎng)格對(duì)復(fù)雜區(qū)域邊界和約束情形有很強(qiáng)的適應(yīng)能力，可以更好地適應(yīng)不規(guī)則計(jì)算區(qū)域和快速實(shí)現(xiàn)網(wǎng)格生成自動(dòng)化，在網(wǎng)格生成中也只需要對(duì)異常區(qū)和測(cè)點(diǎn)處做加密處理。因此，利用非結(jié)構(gòu)化的網(wǎng)格剖分，可以提高有限元計(jì)算的精度，同時(shí)也方便對(duì)不規(guī)則異常體的模擬。

本文作者通過編程實(shí)現(xiàn)了二維模型的非結(jié)構(gòu)化網(wǎng)格剖分，針對(duì)實(shí)際生產(chǎn)中常用的中心回線模型，編寫了瞬變電磁2.5維問題的有限元正演程序。通過程序計(jì)算表明，在取得相同計(jì)算精度的情況下，與結(jié)構(gòu)化網(wǎng)格相比，非結(jié)構(gòu)化網(wǎng)格所需網(wǎng)格和節(jié)點(diǎn)數(shù)量大大減少，計(jì)算效率更高。另外，采用非結(jié)構(gòu)化網(wǎng)格，容易實(shí)現(xiàn)復(fù)雜二維地電模型的數(shù)值模擬，也可以提高現(xiàn)有有限元算法的應(yīng)用范圍。

1 2.5維瞬變電磁法的有限元法

對(duì)于一個(gè)如圖1所示的中心回線瞬變電磁2.5維問題，它是一個(gè)關(guān)于 (x，y，z，t)的四維電磁場(chǎng)問題。對(duì)時(shí)間變量t作拉普拉斯變換，對(duì)走向y方向做傅里葉變化，可將此連續(xù)的四維電磁場(chǎng)問題降維為離散的頻率波數(shù)域中的關(guān)于x、z方向的二維問題。通過一系列的公式推導(dǎo)，可將頻率波數(shù)域中的中心回線法瞬變電磁場(chǎng)2.5維邊值問題表示如下［8］:

在這里:ub=eb，y

其中 Ω為求解區(qū)域;G為求解區(qū)域外邊界;eb，y、hb，y與 ea，y、ha，y分別為頻率波數(shù)域中電、磁場(chǎng)的走向分量的背景場(chǎng)、異常場(chǎng);m為離散的波數(shù)集;s為離散的頻率集;k為介質(zhì)的傳播系數(shù);ε、μ、σ分別為介質(zhì)的介電常數(shù)、磁導(dǎo)率和電導(dǎo)率，δσ=σσb為相對(duì)異常電導(dǎo)率。

對(duì)式(1)通過區(qū)域剖分(文中采用三角形網(wǎng)格)、線性插值、單元分析、總體合成、求變分等步驟，便可得到關(guān)于變分問題式(1)的線性方程組，詳細(xì)過程可參見文獻(xiàn)［9］，這里不再贅述。

解此線性方程組，便可求得頻率波數(shù)域中的各節(jié)點(diǎn)u、v值。

圖1 中心回線瞬變電磁2.5維問題示意圖Fig. 1 The schematic diagram of 2. 5 － D TEM witha central loop

2 測(cè)點(diǎn)感應(yīng)電動(dòng)勢(shì)計(jì)算

通過有限元計(jì)算求得的u、v值，為頻率波數(shù)域中異常體走向分量的電場(chǎng)、磁場(chǎng)異常場(chǎng)的場(chǎng)值，可根據(jù)式(2)求得頻率波數(shù)域中z分量的磁場(chǎng)異常場(chǎng)場(chǎng)值［8］:

再對(duì)ha，z反拉氏變換和反傅氏變化之后，便可求得給定點(diǎn)空間異常場(chǎng)分量的瞬變響應(yīng)Ha，z(t)。

根據(jù)式(3)，可求得異常場(chǎng)形成的感應(yīng)電動(dòng)勢(shì)εa(t)。

式中 S為接收線圈的等效面積;μ0為真空中的磁導(dǎo)率。

根據(jù)式(4)，可計(jì)算出我們平常測(cè)量中所關(guān)心的總場(chǎng)感應(yīng)電動(dòng)勢(shì)ε(t)。

其中 εb(t)為背景場(chǎng)形成的感應(yīng)電動(dòng)勢(shì)。

3 非結(jié)構(gòu)化網(wǎng)格法

網(wǎng)格可分為結(jié)構(gòu)化和非結(jié)構(gòu)化兩類。結(jié)構(gòu)化網(wǎng)格，每個(gè)內(nèi)部節(jié)點(diǎn)都被相同數(shù)目的單元所包含;而非結(jié)構(gòu)化網(wǎng)格中，包含每個(gè)內(nèi)部節(jié)點(diǎn)的數(shù)目是不相同的。

Delaunay三角化方法是目前應(yīng)用較為廣泛的非結(jié)構(gòu)化網(wǎng)格生成方法。Delaunay三角化方法是將平面上的一組已經(jīng)給定的點(diǎn)聯(lián)結(jié)成三角形，并具有以下特點(diǎn):

(1)形成的三角形互不重疊。

(2)可以覆蓋整個(gè)計(jì)算區(qū)域。

(3)每一個(gè)點(diǎn)均在不包括該點(diǎn)的三角形的外接圓外。

Delaunay三角化方法具有以下優(yōu)點(diǎn):

(1)用Delaunay三角化方法連接成的三角形，能取得最大的最小角。這意味著，對(duì)給定的一組點(diǎn)，用Delaunay三角化方法生成的三角形的邊長(zhǎng)均勻性是最好的。

(2)Delaunay三角化方法生成的網(wǎng)格，在空間上是全局優(yōu)化的結(jié)構(gòu)，能滿足高精度計(jì)算中需要網(wǎng)格盡量飽滿的要求，而且可以方便地對(duì)已經(jīng)生成的網(wǎng)格進(jìn)行加點(diǎn)和減點(diǎn)操作并保證新生成的網(wǎng)格仍然是Delaunay網(wǎng)格。

而一般的三角網(wǎng)格并不具備這些優(yōu)點(diǎn)［13、14］。

Delaunay算法分為下述兩個(gè)步驟。

3.1 邊界網(wǎng)格的初始化

生成一個(gè)包含所有給定邊界點(diǎn)集P0的大三角形T0。?P∈P0，將外接圓包含P的所有三角形找出，記下這些三角形形成的邊界，并刪除這些三角形。將點(diǎn)P與空洞邊界上每一邊連接，形成新的三角形，將新生成的三角形加入三角形集合T中，并將點(diǎn)P從邊界點(diǎn)集合P0中刪除。

重復(fù)上述操作，直到集合P0為空。將所有包含T0頂點(diǎn)的三角形從集合T中刪除。

3.2 內(nèi)部結(jié)點(diǎn)的自動(dòng)生成

在三角形集合T中，尋找外接圓半徑最大的三角形，在其外接圓心處加入點(diǎn)P。尋找外接圓包含P的所有三角形，將點(diǎn)P與刪除這些三角形后所形成的空洞邊界上每一邊連接形成新的三角形，并將新生成的三角形加入集合T中。對(duì)現(xiàn)存三角形按外接圓半徑大小排序，最大的在序列的頂上。

重復(fù)上述步驟，直到集合T中三角形的最大外接圓半徑，小于某一給定值即可終止迭代。

4 正演算例

4.1 H型模型

H型三層地電斷面模型的參數(shù)為 ρ1=100Ω·m、h1=100m、ρ2=20 Ω·m、h2=100m、ρ3=100Ω·m。采用中心回線裝置，發(fā)送回線邊長(zhǎng)50m，接收線圈等效面積50m2，供電電流1A。

圖2(見下頁(yè))為采用非結(jié)構(gòu)化網(wǎng)格剖分的有限單元解與解析解的對(duì)比圖。

表1為非結(jié)構(gòu)化網(wǎng)格計(jì)算效果與熊彬［9］采用的一個(gè)矩形網(wǎng)格中細(xì)分為四個(gè)三角形網(wǎng)格的計(jì)算效果作的對(duì)比。由表1可以看出，在取得大約同樣計(jì)算精度的情況下，非結(jié)構(gòu)化網(wǎng)格所需單元數(shù)和節(jié)點(diǎn)數(shù)，遠(yuǎn)遠(yuǎn)少于結(jié)構(gòu)化網(wǎng)格所需，因而大大提高了計(jì)算效率。

表1 H型模型非結(jié)構(gòu)化網(wǎng)格與結(jié)構(gòu)化網(wǎng)格計(jì)算效果對(duì)比Tab. 1 The computational efficiency comparison of unstructuredgrid and structured grid for H － model

圖2 H型模型2.5維有限元法解與解析解曲線對(duì)比圖Fig. 2 The comparison of 2. 5D FEM solution and analytical solution for H － model

4.2 低阻圓柱體模型

模型如圖3所示，電阻率為100Ω·m的圍巖中有個(gè)沿一方向無限延伸的電阻率為0.5Ω·m的圓柱體，圓柱體半徑為50m，中心點(diǎn)埋深為100m。采用中心回線裝置，發(fā)送回線邊長(zhǎng)50m，接收線圈等效面積50m2，供電電流1A。

圖3 圓柱體模型圖Fig.3 Themodeldiagramofcylindricalbody

圖4 為測(cè)點(diǎn)位于圓柱體正上方時(shí)，此模型對(duì)應(yīng)的非結(jié)構(gòu)化網(wǎng)格剖分結(jié)果圖。在有限元計(jì)算中，與結(jié)構(gòu)化網(wǎng)格只需一次剖分不同的是，采用非結(jié)構(gòu)化網(wǎng)格當(dāng)測(cè)點(diǎn)移動(dòng)時(shí)，網(wǎng)格需要重新剖分。為了提高計(jì)算精度，作者在測(cè)點(diǎn)和異常體附近均作了網(wǎng)格加密。

下頁(yè)圖5(a)為低阻圓柱體模型的瞬變電磁測(cè)深剖面圖。由圖5(a)可見:①在1.0×10－6s～2.4×10－5s期間，各測(cè)道的剖面曲線都很平緩;②在3.6×10－5s～9.6×10－5s期間，各測(cè)道出現(xiàn)對(duì)稱與球頂?shù)牡凸犬惓?③在9.6×10－5s之后，各測(cè)道均呈現(xiàn)明顯的對(duì)稱與球頂?shù)膯畏瀹惓?，與物理模擬結(jié)果一致［15］。下頁(yè)圖5(b)為低阻圓柱體模型不同測(cè)點(diǎn)處的異常曲線，圖5(b)中dx表示測(cè)點(diǎn)距圓柱體中心點(diǎn)地面投影處的距離。從圖5(b)中可以看出，當(dāng)dx=0m(測(cè)點(diǎn)位于圓柱體中心正上方)時(shí)，異常曲線幅值最大。隨著dx的增大，異常幅值逐漸衰減，并逐漸趨近于均勻半空間的場(chǎng)值。

圖4 測(cè)點(diǎn)位于圓柱體中心地面投影點(diǎn)處時(shí)模型非結(jié)構(gòu)化網(wǎng)格剖分圖Fig. 4 Unstructured mesh of the model when the measuringpoint at cylindrical center's projection point on theground

5 結(jié)論

利用非結(jié)構(gòu)化網(wǎng)格，是模擬復(fù)雜二維地質(zhì)體瞬變電磁響應(yīng)的一種新工具。利用非結(jié)構(gòu)化網(wǎng)格，可允許數(shù)值模型中包含復(fù)雜二維地質(zhì)體，這是傳統(tǒng)的結(jié)構(gòu)化網(wǎng)格所不能做到的。因此，作者在本文的研究對(duì)提高現(xiàn)有瞬變電磁2.5維有限單元算法應(yīng)用范圍有很好的意義。

圖5 低阻圓柱體模型有限元法數(shù)值模擬結(jié)果Fig. 5 The FEM numerical simulation result of low resistance cylindrical body model

通過非結(jié)構(gòu)化三角形網(wǎng)格代替?zhèn)鹘y(tǒng)的結(jié)構(gòu)化網(wǎng)格，對(duì)地電斷面做剖分，對(duì)異常區(qū)和測(cè)點(diǎn)附近做加密處理，在擬合地電斷面的情況下，盡可能地減少了不要的網(wǎng)格的生成。通過文中程序計(jì)算表明，在計(jì)算精度相當(dāng)?shù)那闆r下，采用非結(jié)構(gòu)化網(wǎng)格比結(jié)構(gòu)化網(wǎng)格所需單元數(shù)和節(jié)點(diǎn)數(shù)大大減少，因此可以大幅減少計(jì)算資源，提高了計(jì)算效率。

［1］ 呂國(guó)印.瞬變電磁法的現(xiàn)狀與發(fā)展趨勢(shì)［J］.物探化探計(jì)算技術(shù)，2007，29(增刊):111.

［2］ 薛國(guó)強(qiáng)，李貅，底青云.瞬變電磁法正反演問題研究進(jìn)展［J］.地球物理學(xué)進(jìn)展，2008，23(4):1165.

［3］UNSWORTH M J，TRAVIS B J，CHAVE A D. Electromagneticinduction by a finite electric dipole source overa 2－D earth ［J］. Geophysics，1993，58( 2) : 198.

［4］YONFLIANG MENG，WEIDONG Li， MICHAEL S.ZHDANOV，et al. 2. 5 － D electromagnetic forwardmodeling in the time and frequency domains using the finiteelement method SEG'69 Annual Meeting ExpandedAbstracts［M］. Tulsa: Society of Exploration Aeophysicists，1999.

［5］ MITSUHATA Y. 2 － D electromagnetic modeling by finite－ element method with a dipole source and topography［J］. Geophysics，2000，65( 2) : 465.

［6］ KONG F N， JOHNSTAD S E，R?TEN T，et al. A 2.5D finite － element modeling difference method for marineCSEM modeling in stratified anisotropic media ［J］.2008， 73( 1) : F9.

［7］ XIONG B. A new formula based on the independent electricand magnetic fields for 2. 5 － D forward of TEM［Z］. The 19th international workshop on electromagneticinduction in the earth，2008: 536.

［8］ 王華軍，羅延鐘.中心回線瞬變電磁法2.5維有限單元算法［J］.地球物理學(xué)報(bào)，2003，46(6):855.

［9］ 熊彬，羅延鐘.電導(dǎo)率分塊均勻的瞬變電磁2.5維有限元數(shù)值模擬［J］.地球物理學(xué)報(bào)，2006，49(2):590.

［10］底青云，MARTYNUNSWORTH，王妙月.有限元法2.5維CSAMT數(shù)值模擬［J］.地球物理學(xué)進(jìn)展，2004，19(2):317.

［11］任政勇，湯井田.基于局部加密非結(jié)構(gòu)化網(wǎng)格的三維電阻率法有限元數(shù)值模擬［J］.地球物理學(xué)報(bào)，2009，52(10):2627.

［12］湯井田，王飛燕.基于非結(jié)構(gòu)化網(wǎng)格的2.5D直流電阻率模擬［J］.物探化探計(jì)算技術(shù)，2008，30(5):413.

［13］陳建軍.非結(jié)構(gòu)化網(wǎng)格生成及其并行化的若干問題研究［D］.杭州:浙江大學(xué)博士論文，2006

［14］肖根如，甘衛(wèi)軍，陳為濤.地應(yīng)變計(jì)算Delaunay三角網(wǎng)在MATLAB與GMT環(huán)境下的相互轉(zhuǎn)換［J］.大地測(cè)量與地球動(dòng)力學(xué)，2010，30(3):122.

［15］李金銘.地電場(chǎng)與電法勘探［M］.北京:地質(zhì)出版社，2005

P631.3+25

A

1001—1749(2011)05—0517—05

2011－02－18 改回日期:2011－06－19

趙曉博(1985－)，男，甘肅正寧人，碩士，主要從事瞬變電磁法理論研究。