梁建剛， 甕晶波

(1.中國地質(zhì)調(diào)查局天津地質(zhì)調(diào)查中心，天津 300170;2.中南大學，湖南 長沙 410083)

2.5 維人工源時間域電磁場數(shù)值模擬(正演)問題近些年發(fā)展很快(Leppin，1992;Everett et al.，1993;Mitsuhata，2000;Meng et al.，1999)，這是因為較之于純二維的數(shù)值模擬，2.5 維方法考慮到了場源的三維性質(zhì)，更貼近于野外實際地質(zhì)情況，另一方面，較之真正的三維方法，2.5 維方法大大減少了計算量，得以在普通性能的微機上順利實現(xiàn)(王華軍等，2003;熊彬等，2006;熊彬，2005;金建銘，1998)。

時間域的電磁場問題，實際上是一個空間加時間的四維邊值問題。

瞬變電磁2.5 維有限單元法正演二次場算法中以均勻大地或?qū)訝畲蟮氐乃沧冸姶彭憫?yīng)為一次場，其間異常體的響應(yīng)為二次場(純異常)，二次場的場源為異常體處的一次場，從而避免了計算階躍脈沖的響應(yīng)。其公式推導(dǎo)及實現(xiàn)思路是:首先對電磁場作拉氏變換，將時間域問題變成拉氏域(或稱頻率s 域)中的問題，實現(xiàn)降維(消除時間變量t)，考慮到二維地電條件下ε、μ 和σ 與坐標軸y 無關(guān)，沿y 軸對電磁場作傅氏變換，轉(zhuǎn)換成(s，m)域中的邊值問題。在(s，m)域中實現(xiàn)電磁場的求解，然后再經(jīng)拉氏逆變換、傅氏逆變換后重新回到時間空間域。

編寫程序的過程中，一次場的求取和傅氏逆變換兩次用到余弦變換，拉氏逆變換則用普遍接受的G-S 變換得以實現(xiàn)。線性方程則在證明剛度矩陣對陣性的基礎(chǔ)上，經(jīng)LDLT實現(xiàn)。

1 (s，m)域中求解

1.1 變分問題

用e，h 表示(s，m)域中的電場和磁場，再經(jīng)過邊界問題的推導(dǎo)，得到瞬變電磁2.5 維正演等價的泛函分析(王華軍等，2003):

式中，up= ep，y，vp= - ihp，yu = es，y，v = hs，y，=

e，y，hp，y，es，y，hs，y分別為拉式傅氏(s，m)域沿走向方向的一次場，二次場分量，k 為介質(zhì)的傳播系數(shù)，m 為沿走向的傅氏域波數(shù)，s 為拉式空間變量，ε為介質(zhì)的介電常數(shù)，μ 為磁導(dǎo)率，σ、δσ 為介質(zhì)的電導(dǎo)率及異常電導(dǎo)率(異常體與均勻介質(zhì)電導(dǎo)率之差)，ηu、ηv為區(qū)域邊界系數(shù)，i 為復(fù)數(shù)單位。

1.2 有限單元法

用有限單元法求解變分問題，將(1)式離散化后變?yōu)?

1.3 矩形剖分條件下單元分析:

單元積分1:

單元積分2:

單元積分3:

單元積分4:

單元積分5:

為了使總的剛度矩陣對稱，需要對反對稱矩陣做適當線性變換:

單元積分6:

單元積分7:

單元積分8:

單元積分9:

單元積分10:

單元積分11:

1.4 總體合成

對上述求得的子剛度矩陣先由4 ×4 的矩陣擴展為8 ×8 的矩陣，然后合成總體矩陣為:

再者，各單元擴展后的列向量W，Wp是相同的，所以各單元積分相加時只需將ke和kpe相加，即:

1.5 求變分

對(15)式求變分，并令其為零，δJ(u，v)=δWTKW + WTKδW + δWTKpWp= 0 考慮到K 的對稱性，有δWTKW = WTKδW，并且δWT≠0，故而有

在(16)式中，K、Kp由單元分析求得，Wp是異常體處的一次電磁場，由解析式求得。求解方程組可得到W，即二次電場、磁場分量。

這樣就把變分問題轉(zhuǎn)化為求解線性方程組，這是有限單元法的精髓。

1.6 線性方程組的求解

經(jīng)驗證，所求得的線性方程組的矩陣是對稱的稀疏矩陣。求解這樣的方程組時，若完全套用一般格式，則不僅要存貯大量的零元素，而且還要讓它們參加運算，在存貯量、計算量方面定會造成巨大浪費。本文按照徐世浙(1994)中的做法，等帶寬存貯稀疏剛度矩陣，并用書中所附LDLT程序進行求解(在橫向剖分40個單元，縱向剖分30個單元的情況下半帶寬為66)。

2 拉氏逆變換、傅氏逆變換及其關(guān)鍵點

經(jīng)上述有限元計算求解出域的電磁場值，必須進行拉氏逆變換、傅氏逆變換才能轉(zhuǎn)換成時間空間域的電磁場值。

2.1 拉氏逆變換

考慮到Gaver-Stehfest 變換是純實數(shù)運算，而且只需對較少的拉氏變量作計算，是一種計算較快的算法，本文直接用它來作拉氏逆變換，實現(xiàn)頻域到時域的轉(zhuǎn)換(羅延鐘等，2000;熊彬，2005)。

設(shè)拉氏空間的像函數(shù)為F(s)，實空間的像原函數(shù)為f(t)則:

其中N 是取決于計算機位數(shù)的正整數(shù)。本文中取N= 12。

2.2 傅氏逆變換及一次場計算

由于拉氏逆變換算法對精度要求高，必須先做拉氏逆變換再做傅氏逆變換(圖1)。傅氏逆變換主要完成從波速域向三維(x，y，z)空間的轉(zhuǎn)換，計算公式為(肖明順，2008;昌彥君等，1995):

在瞬變電磁2.5 維正演中，上式有三個顯著的意義:

(1)當y ≠0 時，即計算非主剖面電磁場，也就是旁測剖面的電磁場值，這是2.5 維算法區(qū)別于純2 維算法的顯著特征之一。

(2)當z ＞0 時，即計算地下介質(zhì)中的電磁場值。這為井中瞬變電磁的正演提供了途經(jīng)。

(3)上式中的參數(shù)對應(yīng)著電磁法中的收發(fā)矩參數(shù)，對于瞬變電磁法的中心回線裝置，收發(fā)矩x= 0。

因此，完整的2.5 維瞬變電磁正演中的傅氏逆變換應(yīng)該囊括以上三種情況。本文中關(guān)心的主剖面、地表上、零收發(fā)距的瞬變電磁場的參數(shù)為:x =0，y = 0，z = 0。

考慮到位于主剖面的回線源所產(chǎn)生的垂向磁分量在空間域中相對于y 具有偶對稱性，對回線中心的主剖面上y = 0 的hz，s作傅氏逆變換時，正弦項為零，余弦項為1，積分內(nèi)項變?yōu)楹撕瘮?shù)本身。

這樣(19)式就變?yōu)?

加之在進行有限元計算前，一次場的求解過程中，對給定的m 和發(fā)射回線的半邊長b，sin(mb)為常量，而傅氏逆變換又是在相同波數(shù)

序列中進行的，所以在求解一次場ey，hy時不計算sin(mb)，而把它放在傅氏逆變換中。這樣原本偶對稱的hz，s貌似變成了奇對稱。經(jīng)這么處理，瞬變電磁2.5 維有限元正演中的傅氏逆變換就轉(zhuǎn)化為正弦變換。

(20)式進一步變?yōu)?

編程時直接使用(21)式，用正弦變換實現(xiàn)傅氏逆變換。

而二次場計算所涉及的一次場計算應(yīng)用公式文中直接用王華軍等(2003)的計算公式:

其中

(22)、(23)、(24)三式均含有正弦或余弦高振蕩函數(shù)，用正余弦變換進行求解(王華軍，2004)。

至此，程序?qū)崿F(xiàn)中應(yīng)用到兩次正弦變換(一次場求取與傅氏逆變換)。

2.3 波數(shù)選擇

本文中傅氏域波數(shù)的選擇，是從兩個方面考慮的:一方面，有限元計算時選擇的波數(shù)跟進行傅氏逆變換時的波數(shù)是一樣的，本文中的傅氏逆變換最后轉(zhuǎn)化為正弦變換;另一方面，根據(jù)Hz 隨波數(shù)變化情況，即不同波數(shù)范圍對積分的貢獻程度。綜合上述兩個方面，程序?qū)崿F(xiàn)時波數(shù)范圍的選擇方案為，21個波數(shù)指數(shù)等間隔取為10-8～1 范圍內(nèi)正弦變換的節(jié)點enΔ/k，在進行傅氏逆變換之前用三次樣條插值法指數(shù)等間隔插值成161個波數(shù)及場值。

3 程序流程圖

實現(xiàn)2 維模型多測點剖面的瞬變電磁2.5 維有限單元法正演需要有四重循環(huán)，從里層到外層分別為:①測點;②s(拉氏逆變換);③m(波數(shù)或傅氏逆變換);④t(采樣時間)。具體的程序流程如圖1。

4 正演結(jié)果

4.1 模型1(低阻)

為了驗證2.5 維程序的正確性及精度，用2.5維程序計算了有解析解的一維模型。

圖1 程序流程圖Fig.1 The flow chart of procedure

模型1 的參數(shù)為:ρ1= 100 Ω·m，h1= 100 m，ρ2= 20 Ω·m，h2= 50 m，ρ3= 100 Ω·m。中心回線裝置的發(fā)射回線邊長。

從圖2 可以看出，純異常(二次場)有限元解與解析解曲線形態(tài)完全相同，但幅值有一定的誤差。其中在異常的峰值點相對誤差為5.49%。

4.2 模型2

模型的參數(shù)為:ρ1= 100 Ω·m，h1= 100 m，ρ2= 20 Ω·m，h2= 50 m，ρ3= 100 Ω·m。中心回線裝置的發(fā)射回線邊長。

從圖3 可以看出，純異常(二次場)有限元解與解析解曲線形態(tài)完全相同，但幅值有一定的誤差，且在時間軸方向，有限元解與解析解之間有大約一個采樣間隔的時間位移，兩條曲線的其中在異常的峰值點相對誤差為0.7%。

圖2 模型1 二次場比較圖Fig.2 The second field compare chart of model 1

在 Inter(R)core(TM)i5-2540 M CPU，2.60 GHz，4.00G 內(nèi)存的筆記本電腦上，剖分30 ×40個單元格，異常體5 ×40個單元格，計算16個波數(shù)，41個時間需要602 s。

圖3 模型2 二次場比較圖Fig.3 The second field compare chart of model 2

綜合模型1 和模型2，可以看到有限元解與解析解曲線形態(tài)完全相同，在時間軸方向均有一定的時間位移(有限元解滯后于解析解)，經(jīng)分析認為是在兩種介質(zhì)的界面，解析解可以實現(xiàn)瞬間過渡，而有限元則對應(yīng)于兩個不同物性特征的剖分單元，可以在單元剖分時在界面處剖分加密，但永遠無法實現(xiàn)瞬間過渡?；谕瑯釉颍诩兌萎惓５耐砥诜聪虍惓?，有限元解滯后于解析解，且幅值相對較小。

4.3 模型3

模型的參數(shù)為:均勻半空間ρ1= 100 Ω·m，異常體h1= 100 m，厚度為10 m，寬100 m(x 為450 ～550 m，中心位置x =500 m)，ρ2=10 Ω·m(模擬采空區(qū)，如圖4)。中心回線裝置的發(fā)射回線邊長2b =50 m，測量范圍x 為350 ～650 m。

圖4 模型3 地電斷面示意圖Fig.4 The geological-electricity sketch map of model 3

圖5 為異常中心左側(cè)11個測點的二次場的場值，從圖上可以看到，隨著測量點由遠到近逐漸接近異常體中心，二次場的異常范圍逐漸變小，二次場峰值逐漸變大，峰值對應(yīng)的時間逐漸變小(其中異常中心位置二次場極值出現(xiàn)在32 μs)。這一現(xiàn)象可以用波場傳播的角度解釋，隨著由遠到近電磁波遇到異常體的距離也逐漸減小，響應(yīng)時間逐漸提前，異常效應(yīng)逐漸變大。

圖5 各測點二次場圖Fig.5 The second field chart of different measuring points

圖6 為21個測點的總場的衰減曲線?？梢詤⒄請D5 解釋圖6 中的各種現(xiàn)象，從上圖可以看到，高阻圍巖條件下的低阻異常體模型，總場衰減曲線大致經(jīng)過如下三個階段:

(1)二次場異常在早期相對一次場(總場)較小，500 m 處(異常體中心位置)純二次場在32 μs時早期假異常達到最大，但這時二次場相對一次場很小，在總場曲線上幾乎沒有反應(yīng)。

(2)181 ～512 μs 看到的明顯的正異常。

圖6 測量剖面總場衰減曲線圖Fig.6 The total field decline chart of measuring profile

(3)1.024 ms 之后的貌似與一維正演縱向(時間)特征相矛盾的反?，F(xiàn)象。為此加大測量范圍并提取了8.192 ms 時刻的二次場曲線。

二次場在8.192 ms時刻均為正值，與181 μs時刻方向相同，之所以出現(xiàn)圖5 中的現(xiàn)象，是要因為圖中沒有畫出一次場的水平，而由于地電斷面橫向不均勻性，二次場在異常體中心位置較小，其最大值出現(xiàn)在距離異常中心410 m 和590 m 處。這一現(xiàn)象是1 維正演中無法看到了的，是煙圈效應(yīng)的體現(xiàn)(牛之璉，2007)。

另外曲線呈現(xiàn)出較好的對稱性，是對熊彬等(2006)的修正。

從圖7 可以看到模型3 地電條件下，8.192 ms時異常場的數(shù)量級已經(jīng)到pV，是否能夠有效探測跟儀器的靈敏度、地質(zhì)噪聲等密切相關(guān)。

圖7 8.192 ms 二次場剖面圖Fig.7 The second field profile at 8.192 ms

5 結(jié)論

(1)本文以前人推導(dǎo)的2.5 維瞬變電磁泛函公式為基礎(chǔ)，進一步梳理了矩形剖分下的有限元正演過程，對程序?qū)崿F(xiàn)過程中的一次場求取、傅氏逆變換、拉氏逆變換以及剛度矩陣對稱性及其解法進行了詳細的闡述，并運用自己編寫的程序計算了三個模型。

(2)瞬變電磁2.5 維正演較一維正演更貼近客觀事實，尤其是針對煤田巷道等模型，模型計算出的響應(yīng)可以部分解釋瞬變電磁波場的眼圈效應(yīng);相對三維正演又節(jié)約時間。

(3)瞬變電磁2.5 維正演在計算精度、運算速度上還有待進一步提高。在精度提高的情況下，可以計算地質(zhì)體的電場響應(yīng)，并與地質(zhì)噪聲、儀器噪聲、儀器靈敏度相比較，從而實現(xiàn)瞬變電磁勘察2D地質(zhì)體的可行性分析。

昌彥君，張桂青.1995.電磁場從頻率域轉(zhuǎn)換到時間域的幾種算法比較［J］.物探化探計算技術(shù)，17(3):25-29.

金建銘. 1998.電磁場有限元方法［M］. 西安:西安電子科學出版社:1-34.

羅延鐘，昌彥君. 2000.G-S 變換的快速算法［J］.地球物理學報，43(5):684.

牛之璉.2007.時間域電磁法原理［M］. 長沙:中南大學出版社:37-48.

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

王華軍. 2004. 正余弦變換的數(shù)值濾波算法［J］. 工程地球物理學報，1(4):329-335.

肖明順. 2008. 帶地形的瞬變電磁2. 5 維有限元數(shù)值模擬研究［D］.武漢:中國地質(zhì)大學:1-50.

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

熊彬. 2005.關(guān)于瞬變電磁法2.5 維正演中的幾個問題［J］. 物探化探計算技術(shù)，28(2):124-128.

徐世浙.1994 地球物理中的有限單元法［M］.北京:科學出版社:1-157.

Everett M E，Edwards R N. 1993. Transient marine electromagnetics:the 2.5D forward problem［J］. Geophys.Int，113(3):545-561.

Leppin M. 1992.Electromagnetic modeling of 3D sources over 2D inhomogeneties in the time domain［J］. Geophysics，57(2):994-1003.

Meng YL，Li W D，Zhdanov M，et al. 1999.2.5-D eletromagnetic forward modeling in the time and frequency domains using the finite element method［A］//SEG’69 annual Meeting Expanded Abstracts，Tulsa:Society of Exploration Geophysicists:1-7.

Mitsuhata Y J. 2000. 2-D electromagnetic modeling by finite element method with a dipole sourceand topography［J］. Geophysics，65(4):465-475.