基于非結構化網格的瞬變電磁2.5維有限元正演模擬

趙曉博，朱自強，李建慧，彭凌星

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

基于非結構化網格的瞬變電磁2.5維有限元正演模擬

趙曉博，朱自強，李建慧，彭凌星

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

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

非結構化網格;Delaunay三角化;瞬變電磁;有限單元法

0 前言

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

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

目前，國內、外對于人工源電磁場的2.5維有限元數值模擬已做了大量研究，并取得了豐碩的成果，但基本上是基于結構化網格實現的［3～10］。結構化網格在人工源電磁場問題的數值模擬中，有下面兩個缺點:

(1)結構化網格的幾何適應能力通常較差，因此對于較復雜計算區域，很難控制網格節點的分布情況，并會產生較大的離散誤差，從而影響有限元計算精度。

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

而相對于結構化網格，非結構化網格對復雜區域邊界和約束情形有很強的適應能力，可以更好地適應不規則計算區域和快速實現網格生成自動化，在網格生成中也只需要對異常區和測點處做加密處理。因此，利用非結構化的網格剖分，可以提高有限元計算的精度，同時也方便對不規則異常體的模擬。

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

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

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

在這里:ub=eb，y

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

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

解此線性方程組，便可求得頻率波數域中的各節點u、v值。

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

2 測點感應電動勢計算

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

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

根據式(3)，可求得異常場形成的感應電動勢εa(t)。

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

根據式(4)，可計算出我們平常測量中所關心的總場感應電動勢ε(t)。

其中 εb(t)為背景場形成的感應電動勢。

3 非結構化網格法

網格可分為結構化和非結構化兩類。結構化網格，每個內部節點都被相同數目的單元所包含;而非結構化網格中，包含每個內部節點的數目是不相同的。

Delaunay三角化方法是目前應用較為廣泛的非結構化網格生成方法。Delaunay三角化方法是將平面上的一組已經給定的點聯結成三角形，并具有以下特點:

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

(2)可以覆蓋整個計算區域。

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

Delaunay三角化方法具有以下優點:

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

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

而一般的三角網格并不具備這些優點［13、14］。

Delaunay算法分為下述兩個步驟。

3.1 邊界網格的初始化

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

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

3.2 內部結點的自動生成

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

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

4 正演算例

4.1 H型模型

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

圖2(見下頁)為采用非結構化網格剖分的有限單元解與解析解的對比圖。

表1為非結構化網格計算效果與熊彬［9］采用的一個矩形網格中細分為四個三角形網格的計算效果作的對比。由表1可以看出，在取得大約同樣計算精度的情況下，非結構化網格所需單元數和節點數，遠遠少于結構化網格所需，因而大大提高了計算效率。

表1 H型模型非結構化網格與結構化網格計算效果對比Tab. 1 The computational efficiency comparison of unstructuredgrid and structured grid for H － model

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

4.2 低阻圓柱體模型

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

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

圖4 為測點位于圓柱體正上方時，此模型對應的非結構化網格剖分結果圖。在有限元計算中，與結構化網格只需一次剖分不同的是，采用非結構化網格當測點移動時，網格需要重新剖分。為了提高計算精度，作者在測點和異常體附近均作了網格加密。

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

圖4 測點位于圓柱體中心地面投影點處時模型非結構化網格剖分圖Fig. 4 Unstructured mesh of the model when the measuringpoint at cylindrical center's projection point on theground

5 結論

利用非結構化網格，是模擬復雜二維地質體瞬變電磁響應的一種新工具。利用非結構化網格，可允許數值模型中包含復雜二維地質體，這是傳統的結構化網格所不能做到的。因此，作者在本文的研究對提高現有瞬變電磁2.5維有限單元算法應用范圍有很好的意義。

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

通過非結構化三角形網格代替傳統的結構化網格，對地電斷面做剖分，對異常區和測點附近做加密處理，在擬合地電斷面的情況下，盡可能地減少了不要的網格的生成。通過文中程序計算表明，在計算精度相當的情況下，采用非結構化網格比結構化網格所需單元數和節點數大大減少，因此可以大幅減少計算資源，提高了計算效率。

［1］ 呂國印.瞬變電磁法的現狀與發展趨勢［J］.物探化探計算技術，2007，29(增刊):111.

［2］ 薛國強，李貅，底青云.瞬變電磁法正反演問題研究進展［J］.地球物理學進展，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］.地球物理學報，2003，46(6):855.

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

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

［11］任政勇，湯井田.基于局部加密非結構化網格的三維電阻率法有限元數值模擬［J］.地球物理學報，2009，52(10):2627.

［12］湯井田，王飛燕.基于非結構化網格的2.5D直流電阻率模擬［J］.物探化探計算技術，2008，30(5):413.

［13］陳建軍.非結構化網格生成及其并行化的若干問題研究［D］.杭州:浙江大學博士論文，2006

［14］肖根如，甘衛軍，陳為濤.地應變計算Delaunay三角網在MATLAB與GMT環境下的相互轉換［J］.大地測量與地球動力學，2010，30(3):122.

［15］李金銘.地電場與電法勘探［M］.北京:地質出版社，2005

P631.3+25

A

1001—1749(2011)05—0517—05

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

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