瞬變電磁快速三維正演

李瑞雪，王鶴, 2，席振銖, 2，蔣歡，劉愿愿

?

瞬變電磁快速三維正演

李瑞雪1，王鶴1, 2，席振銖1, 2，蔣歡1，劉愿愿1

(1. 中南大學(xué)地球科學(xué)與信息物理學(xué)院，湖南長沙，410083；2. 中南大學(xué)海洋礦產(chǎn)探測(cè)技術(shù)與裝備研究所，湖南長沙，410083)

為了提高瞬變電磁法三維正演計(jì)算速度，一方面，采用矢量有限單元法與積分方程法相混合，使24維單元矩陣的階數(shù)降低到12維，減少大型稀疏矩陣的計(jì)算時(shí)間，同時(shí)增強(qiáng)邊界條件的連續(xù)性，避免出現(xiàn)偽解；另一方面，引進(jìn)Flow-through Hankel transform快速計(jì)算技術(shù)，解決不同收發(fā)距的貝塞爾函數(shù)重復(fù)計(jì)算的問題，進(jìn)一步節(jié)省計(jì)算時(shí)間。研究結(jié)果表明：引入矢量有限單元法的混合法和Flow-through Hankel transform快速計(jì)算技術(shù)不但能保證3維計(jì)算精度，而且能提高計(jì)算速度。

瞬變電磁法；3維正演；Flow-through Hankel transform；矢量有限單元法

瞬變電磁法以其靈敏度高、探測(cè)深度大、分辨率高以及適應(yīng)性強(qiáng)的優(yōu)點(diǎn)，已經(jīng)迅速發(fā)展為地球物理勘探方法中常用方法，但其正、反演技術(shù)進(jìn)展緩慢，特別是3維計(jì)算速度不能滿足實(shí)際工程應(yīng)用的要求[1?2]。瞬變電磁3維正演方法主要有有限差分法、有限單元法和積分方程法[3?4]。其中，有限單元法和有限差分法計(jì)算速度慢，占用內(nèi)存大[5?7]；積分方程法只適合模擬簡(jiǎn)單模型[8?9]。BAKR等[10]使用迭代混合解法求解擴(kuò)散電磁場(chǎng)，該方法先設(shè)定1個(gè)邊界場(chǎng)的初始值，求出內(nèi)部場(chǎng)，然后根據(jù)傅里葉變換求出新的邊界場(chǎng)，一直重復(fù)這個(gè)過程，直到計(jì)算結(jié)果誤差在規(guī)定范圍內(nèi)為止。這種算法適用于異常體與圍巖的電導(dǎo)率差很小的模型，但當(dāng)電導(dǎo)率差達(dá)200倍時(shí)，這種算法精度不高。CHUNG等[11]提出了1種基于求解二次磁場(chǎng)的正演方法，克服了基于總場(chǎng)求解只適合于電導(dǎo)率差小的模型的缺點(diǎn)。GUPTA等[12?13]將節(jié)點(diǎn)有限元法和積分方程法結(jié)合進(jìn)行瞬變電磁3維正演，在保證計(jì)算精度的基礎(chǔ)上成倍提高了計(jì)算速度，但會(huì)導(dǎo)致偽解，并引起邊界條件不連續(xù)。在瞬變電磁的數(shù)值模擬中，麥克斯韋方程推導(dǎo)的電磁場(chǎng)的表達(dá)式是含0階或者1階貝塞爾函數(shù)的積分形式，因?yàn)樨惾麪柡瘮?shù)不收斂，求解這類積分必須用到漢克爾變換。國內(nèi)外學(xué)者針對(duì)這方面進(jìn)行了大量研究。GUPTASARMA等[14]求解了0階漢克爾變換的61點(diǎn)和120點(diǎn)濾波系數(shù)、1階漢克爾變換的47點(diǎn)和140點(diǎn)濾波系數(shù)，提高了計(jì)算速度和精度。但瞬變電磁法二次場(chǎng)的計(jì)算將異常體作為場(chǎng)源，導(dǎo)致收發(fā)距在每個(gè)測(cè)點(diǎn)都不同，因此，需要重復(fù)進(jìn)行漢克爾變換求解貝塞爾函數(shù)，耗費(fèi)大量的計(jì)算時(shí)間。RAICHE[15]使收發(fā)距與漢克爾濾波系數(shù)同步對(duì)數(shù)等間隔變化，避免了重復(fù)計(jì)算貝塞爾函數(shù)。為此，本文作者一方面將矢量有限單元法和積分方程法結(jié)合實(shí)現(xiàn)瞬變電磁3維正演中大型矩陣的降維快速計(jì)算，另一方面運(yùn)用Flow-through漢克爾變換快速計(jì)算技術(shù)，解決不同收發(fā)距貝塞爾函數(shù)重復(fù)計(jì)算的問題，進(jìn)一步節(jié)約計(jì)算時(shí)間。為了表達(dá)方便，將該方法稱為改進(jìn)混合 算法。

1 理論基礎(chǔ)

1.1 降維快速計(jì)算原理

在各向同性的3維均勻介質(zhì)中，忽略位移電流和磁流密度，電磁場(chǎng)滿足麥克斯韋方程[13]：

其中：()，()，，和()分別為電場(chǎng)強(qiáng)度、磁場(chǎng)強(qiáng)度、磁導(dǎo)率、電導(dǎo)率和電流密度。

對(duì)式(1)求散度，并將式(2)代入，同時(shí)將電場(chǎng)分解為一次電場(chǎng)和二次電場(chǎng)，可得二次電場(chǎng)的擴(kuò)散方程：

其中：()，()和分別為二次電場(chǎng)、一次電場(chǎng)和電導(dǎo)率異常。

對(duì)異常體和邊界圍巖區(qū)域進(jìn)行網(wǎng)格剖分，剖分網(wǎng)格為線性等參六面體，在每個(gè)單元內(nèi)二次電場(chǎng)的近似解為[16]

則單元內(nèi)的余量為

引入伽遼金法，令單元內(nèi)余量加權(quán)為0，得

將式(6)代入式(7)，考慮單元內(nèi)全部12個(gè)形函數(shù)，則在每個(gè)單元內(nèi)式(7)可寫為

其中：為12×12階單元矩陣；為12×1階列向量；為與源有關(guān)的矩陣。

對(duì)所有單元方程式(8)進(jìn)行組合得方程組：

其中：為稀疏矩陣；1為與源有關(guān)的向量。將二次電場(chǎng)()分為異常體內(nèi)部二次電場(chǎng)和邊界二次電場(chǎng)，則式(9)變?yōu)?/p>

其中：為對(duì)權(quán)矢量和場(chǎng)矢量采用內(nèi)部棱邊形函數(shù)得到的剛度矩陣；為對(duì)權(quán)矢量采用內(nèi)部棱邊形函數(shù)、對(duì)場(chǎng)矢量采用邊界棱邊形函數(shù)得到的剛度矩陣。和有類似的含義。

將式(10)寫為總場(chǎng)的形式，并引入積分方程法[13]，可得

其中：G1和G2分別為一次電場(chǎng)和二次電場(chǎng)的格林函數(shù)。式(11)方程組階數(shù)通常是幾百階，可采取直接解法求解，求得一次電場(chǎng)強(qiáng)度。最后，根據(jù)

可求出二次磁場(chǎng)強(qiáng)度。其中：()為層狀介質(zhì)中的一次磁場(chǎng)強(qiáng)度；G2為二次磁場(chǎng)的格林函數(shù)。

1.2 Flow-through快速漢克爾變換

考慮含貝塞爾函數(shù)的積分：

其中：J0()為0階第一類貝塞爾函數(shù)；為源點(diǎn)和接收點(diǎn)的距離，且

式中：為采樣起始點(diǎn)；為采樣間隔；=1， 2，…，。

由式(14)可見：()只與收發(fā)距有關(guān)。但瞬變電磁法中二次場(chǎng)的計(jì)算將異常體看作發(fā)射源，因此，每個(gè)測(cè)點(diǎn)的收發(fā)距都不同，常規(guī)的快速漢克爾變換在每個(gè)收發(fā)距下都得進(jìn)行1次，浪費(fèi)了大量重復(fù)計(jì)算時(shí)間。Flow-through 漢克爾變換[15]提供了1種新的計(jì)算思路。令漢克爾變換的系數(shù)和收發(fā)距都呈對(duì)數(shù)等間隔分布，式(13)將變?yōu)?/p>

其中：為常數(shù)；為對(duì)數(shù)等間隔分布的收發(fā)距；w為漢克爾濾波系數(shù)；x為漢克爾濾波系數(shù)橫坐標(biāo)值。

從式(15)可以看出漢克爾變換的計(jì)算時(shí)間主要受2個(gè)方面的影響：1) 漢克爾系數(shù)的節(jié)點(diǎn)數(shù)；2) 函數(shù)值的計(jì)算次數(shù)。其中，第2個(gè)方面的影響遠(yuǎn)大于第1個(gè)方面的影響，因此，可以通過對(duì)函數(shù)值的重復(fù)利用來減少計(jì)算時(shí)間。故令收發(fā)距對(duì)數(shù)等間隔分布，在每個(gè)數(shù)量級(jí)取15個(gè)點(diǎn)進(jìn)行計(jì)算。令最小的收發(fā)距1比實(shí)際測(cè)量中的最小收發(fā)距小，最大的收發(fā)距k比實(shí)際測(cè)量中的最大收發(fā)距大。因?yàn)?i>x也是對(duì)數(shù)等間隔分布的，定義

則式(15)變?yōu)?/p>

對(duì)1個(gè)固定的收發(fā)距k來說，定義它所需的最小漢克爾系數(shù)為n，它所需的最大漢克爾系數(shù)為n。顯而易見，計(jì)算最小收發(fā)距1時(shí)需計(jì)算最大漢克爾系數(shù)n，計(jì)算最大收發(fā)距k時(shí)需計(jì)算最小漢克爾系數(shù)n，故在實(shí)際程序中需用n到n的漢克爾系數(shù)。因此，

為了滿足所有收發(fā)距范圍，令收發(fā)距從10?1m到104m呈對(duì)數(shù)等間隔分布。在計(jì)算過程中，首先對(duì)于最小收發(fā)距1，求取對(duì)應(yīng)濾波系數(shù)n1到n的函數(shù)值并存儲(chǔ)，同時(shí)計(jì)算此時(shí)的積分值。其次，對(duì)于最大收發(fā)距k，求濾波系數(shù)n到n1對(duì)應(yīng)的函數(shù)值并存儲(chǔ)，同時(shí)求取積分值。在計(jì)算其他收發(fā)距時(shí)，直接調(diào)用已經(jīng)儲(chǔ)存的函數(shù)值，節(jié)約計(jì)算時(shí)間。

程序偽代碼如下：

!n是1和2的中間值

=1

DO=0，2

n=

計(jì)算并存儲(chǔ)

計(jì)算(1)

IF (|Re()|<|Re()|*10?6or |Im()|<|Im()|*10?6)

EXIT

END DO

DO=0?1,1, ?1

n1=

計(jì)算并存儲(chǔ)

計(jì)算(1)

IF (|Re()|<|Re()|*10?6or |Im()|<|Im()|*10?6)

EXIT

END DO

=

DO=n1?1,1, ?1

n=

計(jì)算并存儲(chǔ)

計(jì)算(k)

IF (|Re()|<|Re()|*10?6or |Im()|<|Im()|*10?6)

EXIT

END DO

！求取積分

DO=2，?1

DO=n，n

計(jì)算(k)

END DO

END DO

對(duì)于實(shí)際的任意收發(fā)距，判斷出它在k和k+1之間，通過3次樣條插值算出對(duì)應(yīng)的函數(shù)值。這種計(jì)算方法避免了不同收發(fā)距下的函數(shù)值的重復(fù)計(jì)算，提高了計(jì)算效率，節(jié)約了計(jì)算時(shí)間和內(nèi)存。

2 數(shù)值計(jì)算精度與速度比對(duì)

2.1 數(shù)值計(jì)算精度比對(duì)

檢驗(yàn)數(shù)值計(jì)算精度的最好途徑是與解析解比對(duì)。瞬變電磁法均勻半空間具有解析解，不失一般性，選擇電阻率為100 Ω?m均勻大地進(jìn)行比對(duì)。

對(duì)于電阻率為100 Ω?m的均勻半空間，劃分出長×寬×高為100 m×100 m×100 m的區(qū)域進(jìn)行均勻網(wǎng)格剖分，分別是方向(?50 m，50 m)，方向(?50 m，50 m)，方向(?50 m，?150 m)。這3個(gè)方向各有5個(gè)節(jié)點(diǎn)，在地面布1條測(cè)線，從(?200 m，0 m)到(200 m，0 m)，共11個(gè)測(cè)點(diǎn)。發(fā)射線圈采用100 m×100 m的方形回線，電流為1 A。計(jì)算機(jī)為ThinkPadE430，i5-3210M處理器，4G內(nèi)存，計(jì)算時(shí)間為390 s(11個(gè)測(cè)點(diǎn)，每個(gè)測(cè)點(diǎn)125個(gè)網(wǎng)格)。計(jì)算結(jié)果見圖1，具體誤差見表1。

從圖1可見：程序模擬的均勻半空間曲線和解析解得出的曲線形態(tài)基本一致，也基本重合。由表1可見：除早期的幾個(gè)點(diǎn)外，電流歸一化感應(yīng)電動(dòng)勢(shì)相對(duì)誤差基本在10%以內(nèi)，滿足精度要求，可以用來計(jì)算更復(fù)雜的地電模型。

1—改進(jìn)混合算法解；2—解析算法解。

表1 不同時(shí)刻下改進(jìn)混合算法誤差

2.2 數(shù)值計(jì)算速度比對(duì)

為了驗(yàn)證改進(jìn)混合算法計(jì)算速度，通過模擬直立薄板模型響應(yīng)與2.5維有限單元法、3維積分方程法進(jìn)行比對(duì)。

直立薄板模型如圖2 所示。在電阻率為100 Ω?m的均勻半空間下有電阻率為1 Ω?m的異常體，異常體走向?yàn)榉较?，范圍??200 m，200 m)，沿方向范圍為(?15 m，15 m)，方向范圍為(50 m，200 m)，可以認(rèn)為是直立薄板。在地面布測(cè)線1條，從(?300 m，0 m)到(300 m，0 m)，共25個(gè)測(cè)點(diǎn)，發(fā)射線圈采用100 m×100 m的方形回線，電流為1 A。假設(shè)接收線圈為1個(gè)點(diǎn)，響應(yīng)結(jié)果和計(jì)算時(shí)間見圖3。

由圖3可見，該算法與2.5維有限單元法和積分方程法得到的曲線形態(tài)基本一致，在異常體區(qū)域曲線完全重合?；旌纤惴ǖ挠?jì)算時(shí)間是2.5維有限單元法的10%，與積分方程法相比，計(jì)算速度提高了28%，體現(xiàn)了該算法的快速性。

圖2 直立薄板模型

1—改進(jìn)混合算法，20 min；2—2.5維有限元法，3.0 h；3—3維積分方程法，28 min。

3 3維快速計(jì)算算例

電阻率為100 Ω?m的均勻半空間中地表下方有1個(gè)長方體，其電阻率為1 Ω?m，走向?yàn)榉较颍唧w范圍為(?200 m，200 m)；方向范圍為(?100 m，100 m)，方向范圍為(50 m，250 m)。發(fā)射線圈采用100 m×100 m的方形回線，電流為10 A，接收線圈為1個(gè)點(diǎn)。布置7條測(cè)線，測(cè)線沿方向，7條測(cè)線的坐標(biāo)在(?300 m，300 m)內(nèi)均勻分布；每條測(cè)線13個(gè)測(cè)點(diǎn)，在(?300 m，300 m)內(nèi)均勻分布，具體模型如圖4所示。為了保證計(jì)算精度，對(duì)上述異常體進(jìn)行均勻網(wǎng)格剖分，共剖分為8×5×5個(gè)單元，網(wǎng)格剖分如圖5所示。

圖4 3維長方體模型示意圖

圖5 3維長方體網(wǎng)格剖分示意圖

計(jì)算所用電腦為ThinkPadE430，i5-3210M處理器，4G內(nèi)存，計(jì)算時(shí)間為3 430 s，平均1條測(cè)線(13個(gè)測(cè)點(diǎn)，每個(gè)測(cè)點(diǎn)324個(gè)單元)運(yùn)算時(shí)間為490 s，計(jì)算結(jié)果如圖6所示。因?yàn)楫惓ｓw埋深較淺，因此，在=0.5 ms時(shí)，長方體的響應(yīng)已經(jīng)十分明顯，隨著時(shí)間的延長，二次場(chǎng)逐步衰減，探測(cè)深度逐漸加深，異常體響應(yīng)越來越弱，直到=25 ms時(shí)異常響應(yīng)完全消失。由圖6可知：程序正演結(jié)果所得異常與物理模型中異?；緦?duì)應(yīng)，異常范圍與物理模型異常范圍也基本符合，進(jìn)一步驗(yàn)證了程序的正確性和快速性，也說明了3維正演的必要性。

圖中數(shù)值單位：nT/s

4 結(jié)論

1) 將矢量有限單元法和積分方程法相混合，降低了單元矩陣的階數(shù)，減少了大型稀疏矩陣的計(jì)算時(shí)間，增強(qiáng)了邊界條件的連續(xù)性，提高了3維計(jì)算精度。

2) 引進(jìn)Flow-through Hankel transform快速計(jì)算技術(shù)，解決了不同收發(fā)距的貝塞爾函數(shù)重復(fù)計(jì)算浪費(fèi)時(shí)間的問題，提高了計(jì)算速度。

3) 改進(jìn)混合算法不但能夠保證3維計(jì)算精度，而且提高了計(jì)算速度，促進(jìn)了瞬變電磁3維正演技術(shù)的發(fā)展。

參考文獻(xiàn)：

[1] 牛之鏈. 時(shí)間域電磁法原理[M]. 長沙: 中南大學(xué)出版社, 2007: 1?6. NIU Zhilian. Theory of time domain electromagnetic method[M]. Changsha: Central South University Press, 2007: 1?6.

[2] 李貅. 瞬變電磁測(cè)深的理論與應(yīng)用[M]. 西安: 陜西科學(xué)技術(shù)出版社, 2002: 1?4. LI Xiu. The theory and application of transient electromagnetic sounding[M]. Xi’an: Shanxi Science and Technology Press, 2002: 1?4.

[3] 薛國強(qiáng), 李貅, 底青云. 瞬變電磁法正反演問題研究進(jìn)展[J]. 地球物理學(xué)進(jìn)展, 2008, 23(4): 1165?1172. XUE Guoqiang, LI Xiu, DI Qingyun. Research progress in TEM forward modeling and inversion calculation[J]. Progress in Geophysics, 2008, 23(4): 1165?1172.

[4] 李建慧, 朱自強(qiáng), 曾思紅, 等. 瞬變電磁法正演計(jì)算進(jìn)展[J]. 地球物理學(xué)進(jìn)展, 2012, 27(4): 1393?1400. LI Jianhui, ZHU Ziqiang, ZENG Sihong, et al. Progress of forward computation in transient electromagnetic method[J]. Progress in Geophysics, 2012, 27(4): 1393?1400.

[5] EPOV M I, SHURINA E P, NECHAEV O V. 3D forward modeling of vector field for induction logging problems[J]. Russian Geology and Geophysics, 2007, 48(9): 770?774.

[6] 李建慧. 基于矢量有限單元法的大回線源瞬變電磁法三維數(shù)值模擬[D]. 長沙: 中南大學(xué)地球科學(xué)與信息物理學(xué)院, 2011: 3?7. LI Jianhui. 3D numerical simulation for transient electromagnetic field excited by large source loop based on vector finite element method[D]. Changsha: Central South University School of Geosciences and Info-Physics, 2011: 3?7.

[7] 孫懷鳳, 李貅, 李術(shù)才, 等. 考慮關(guān)斷時(shí)間的回線源激發(fā)TEM三維時(shí)域有限差分正演[J]. 地球物理學(xué)報(bào), 2013, 56(3): 1049?1064.SUN Huaifeng, LI Xiu, LI Shucai, et alThree-dimensional FDTD modeling of TEM excited by a loop source considering ramp time[J]. Chinese Journal of Geophysics, 2013, 56(3): 1049?1064.

[8] SWIDINSKY A, EDWARDS R N. The transient electromagnetic response of a resistive sheet straight forward but not trivial[J]. Geophysics Journal International, 2009, 179(3): 1488?1498.

[9] 張輝, 李桐林, 董瑞霞. 體積分方程法模擬電偶源三維電磁響應(yīng)[J]. 地球物理學(xué)進(jìn)展, 2006, 21(2): 386?390. ZHANG Hui, LI Tonglin, DONG Ruixia. Modeling 3-D electromagnetic responses of the electric dipole using volume integral equation method[J]. Progress in Geophysics, 2006, 21(2): 386?390.

[10] BAKR S A, MANNSETH T. An approximate hybrid method for electromagnetic scattering from an underground target[J]. Geoscience and Remote Sensing, IEEE Transactions on Geoscience and Remote Sensing, 2013, 51(1): 99?107.

[11] CHUNG Y, SON J S, LEE T J, et al. Three-dimensional modelling of controlled-source electromagnetic surveys using an edge finite-element method with a direct solver[J]. Geophysical Prospecting, 2014, 62(6): 1468?1483.

[12] GUPTA P K, BENNETT L A, RAICHE A P. Hybrid calculations of the three-dimensional electromagnetic response of buried conductors[J]. Geophysics, 1987, 52(3): 301?306.

[13] GUPTA P K, RAICHE A P, SUGENG F. Three-dimensional time-domain electromagnetic modelling using a compact finite-element frequency-stepping method[J]. Geophysical Journal International, 1989, 96(3): 457?468.

[14] GUPTASARMA D, SINGH B. New digital linear filters for Hankel J0and J1transforms[J]. Geophysical Prospecting, 1997, 45(5): 745?762.

[15] RAICHE A. A flow-through Hankel transform technique for rapid, accurate Green’s function computation[J]. Radio Science, 1999, 34(2): 549?555.

[16] 金建銘. 電磁場(chǎng)有限元法[M]. 西安: 西安電子科技大學(xué)出版社, 1998: 176?182. JIN Jianming. The finite element method in electromagnetics[M]. Xi’an: Xi Dian University Press, 1998: 176?182.

(編輯 陳燦華)

Fast 3D forward modeling of transient electromagnetic

LI Ruixue1, WANG He1, 2, XI Zhenzhu1, 2, JIANG Huan1, LIU Yuanyuan1

(1. School of Geosciences and Info-Physics, Central South University, Changsha 410083, China;2. Marine Mineral Exploration Technology and Equipment Institute, Central South University, Changsha 410083, China)

In order to enhance the speed of 3D transient electromagnetic forward modeling, two new techniques were employed. One is combining vector finite-element method with integral equation method to decrease the order of elemental matrix from 24 to 12, which can save much time of the computing large sparse matrix, keep the continuity of boundary conditions, and avoid spurious solutions; the other is adopting flow-through Hankel transform technique to resolve the exhausted time problem through calculating. Bessel function repeatedly while resetting the offsets every time. The results show that the two new techniques not only keep the accuracy of 3D TEM forward technique but also speed up the calculation.

transient electromagnetic method; 3D forward modeling; Flow-through Hankel transform; vector finite-element method

10.11817/j.issn.1672-7207.2016.10.026

P319.1+2

A

1672?7207(2016)10?3477?06

2015?10?12；

2015?12?24

大洋“十二五”重大項(xiàng)目(DY125-11-R-03)；國家自然科學(xué)基金資助項(xiàng)目(41304090)；深圳市戰(zhàn)略新興產(chǎn)業(yè)發(fā)展專項(xiàng)資金資助項(xiàng)目(CXZZ20120618165608947)(Project(DY125-11-R-03) supported by the Ocean Major during “Twelfth Five Year Plan”; Project(41304090) supported by the National Natural Science Foundation of China; Project(CXZZ20120618165608947) supported by Strategic Emerging Industry Development Special Foundation of Shenzhen)

席振銖，教授，博士生導(dǎo)師，從事瞬變電磁方法與技術(shù)研究；E-mail：xizhenzhu@163.com