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

?

瞬變電磁快速三維正演

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

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

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

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

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

1 理論基礎

1.1 降維快速計算原理

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

其中：()，()，，和()分別為電場強度、磁場強度、磁導率、電導率和電流密度。

對式(1)求散度，并將式(2)代入，同時將電場分解為一次電場和二次電場，可得二次電場的擴散方程：

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

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

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

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

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

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

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

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

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

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

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

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

1.2 Flow-through快速漢克爾變換

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

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

式中：為采樣起始點；為采樣間隔；=1， 2，…，。

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

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

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

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

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

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

程序偽代碼如下：

!n是1和2的中間值

=1

DO=0，2

n=

計算并存儲

計算(1)

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

EXIT

END DO

DO=0?1,1, ?1

n1=

計算并存儲

計算(1)

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

EXIT

END DO

=

DO=n1?1,1, ?1

n=

計算并存儲

計算(k)

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

EXIT

END DO

！求取積分

DO=2，?1

DO=n，n

計算(k)

END DO

END DO

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

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

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

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

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

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

1—改進混合算法解；2—解析算法解。

表1 不同時刻下改進混合算法誤差

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

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

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

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

圖2 直立薄板模型

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

3 3維快速計算算例

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

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

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

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

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

4 結論

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

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

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

參考文獻：

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

[2] 李貅. 瞬變電磁測深的理論與應用[M]. 西安: 陜西科學技術出版社, 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] 薛國強, 李貅, 底青云. 瞬變電磁法正反演問題研究進展[J]. 地球物理學進展, 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] 李建慧, 朱自強, 曾思紅, 等. 瞬變電磁法正演計算進展[J]. 地球物理學進展, 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]. 長沙: 中南大學地球科學與信息物理學院, 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] 孫懷鳳, 李貅, 李術才, 等. 考慮關斷時間的回線源激發(fā)TEM三維時域有限差分正演[J]. 地球物理學報, 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] 張輝, 李桐林, 董瑞霞. 體積分方程法模擬電偶源三維電磁響應[J]. 地球物理學進展, 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] 金建銘. 電磁場有限元法[M]. 西安: 西安電子科技大學出版社, 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

大洋“十二五”重大項目(DY125-11-R-03)；國家自然科學基金資助項目(41304090)；深圳市戰(zhàn)略新興產(chǎn)業(yè)發(fā)展專項資金資助項目(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)

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