瞬變電磁2.5維有限元正演的程序?qū)崿F(xiàn)過(guò)程

梁建剛， 甕晶波

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

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

時(shí)間域的電磁場(chǎng)問(wèn)題，實(shí)際上是一個(gè)空間加時(shí)間的四維邊值問(wèn)題。

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

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

1 (s，m)域中求解

1.1 變分問(wèn)題

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

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

e，y，hp，y，es，y，hs，y分別為拉式傅氏(s，m)域沿走向方向的一次場(chǎng)，二次場(chǎng)分量，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 有限單元法

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

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

單元積分1:

單元積分2:

單元積分3:

單元積分4:

單元積分5:

為了使總的剛度矩陣對(duì)稱，需要對(duì)反對(duì)稱矩陣做適當(dāng)線性變換:

單元積分6:

單元積分7:

單元積分8:

單元積分9:

單元積分10:

單元積分11:

1.4 總體合成

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

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

1.5 求變分

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

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

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

1.6 線性方程組的求解

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

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

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

2.1 拉氏逆變換

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

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

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

2.2 傅氏逆變換及一次場(chǎng)計(jì)算

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

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

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

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

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

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

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

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

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

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

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

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

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

其中

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

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

2.3 波數(shù)選擇

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

3 程序流程圖

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

4 正演結(jié)果

4.1 模型1(低阻)

為了驗(yàn)證2.5 維程序的正確性及精度，用2.5維程序計(jì)算了有解析解的一維模型。

圖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ā)射回線邊長(zhǎng)。

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

4.2 模型2

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

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

圖2 模型1 二次場(chǎng)比較圖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個(gè)單元格，異常體5 ×40個(gè)單元格，計(jì)算16個(gè)波數(shù)，41個(gè)時(shí)間需要602 s。

圖3 模型2 二次場(chǎng)比較圖Fig.3 The second field compare chart of model 2

綜合模型1 和模型2，可以看到有限元解與解析解曲線形態(tài)完全相同，在時(shí)間軸方向均有一定的時(shí)間位移(有限元解滯后于解析解)，經(jīng)分析認(rèn)為是在兩種介質(zhì)的界面，解析解可以實(shí)現(xiàn)瞬間過(guò)渡，而有限元?jiǎng)t對(duì)應(yīng)于兩個(gè)不同物性特征的剖分單元，可以在單元剖分時(shí)在界面處剖分加密，但永遠(yuǎn)無(wú)法實(shí)現(xiàn)瞬間過(guò)渡?；谕瑯釉?，在純二次異常的晚期反向異常，有限元解滯后于解析解，且幅值相對(duì)較小。

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ā)射回線邊長(zhǎng)2b =50 m，測(cè)量范圍x 為350 ～650 m。

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

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

圖5 各測(cè)點(diǎn)二次場(chǎng)圖Fig.5 The second field chart of different measuring points

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

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

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

圖6 測(cè)量剖面總場(chǎng)衰減曲線圖Fig.6 The total field decline chart of measuring profile

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

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

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

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

圖7 8.192 ms 二次場(chǎng)剖面圖Fig.7 The second field profile at 8.192 ms

5 結(jié)論

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

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

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

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

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

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

牛之璉.2007.時(shí)間域電磁法原理［M］. 長(zhǎng)沙:中南大學(xué)出版社:37-48.

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

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

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

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

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

徐世浙.1994 地球物理中的有限單元法［M］.北京:科學(xué)出版社: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.