劉　堯，李延清，陶鵬飛，胡尊平

（新疆維吾爾自治區(qū)地質(zhì)調(diào)查院，新疆烏魯木齊830000）

瞬變電磁三維數(shù)值模擬進展

劉堯*，李延清，陶鵬飛，胡尊平

（新疆維吾爾自治區(qū)地質(zhì)調(diào)查院，新疆烏魯木齊830000）

介紹了瞬變電磁法三維數(shù)值模擬的研究現(xiàn)狀，對目前流行的幾種正演方法進行了詳細(xì)介紹，并對未來的發(fā)展方向進行了展望。20世紀(jì)80年代以來，隨著計算機技術(shù)的快速發(fā)展，瞬變電磁正演研究已經(jīng)由早期的一維模型的解析方法，轉(zhuǎn)向了高維數(shù)值計算方法。目前許多學(xué)者在從事三維的數(shù)值計算，尤其是針對任意三維復(fù)雜模型的求解更是成為目前的研究熱點。20世紀(jì)80年代，進行時間域三維瞬變電磁法數(shù)值計算主要采用積分方程法，現(xiàn)已很少使用；之后則主要采用：有限單元法、時域有限差分法；另外還有近年來興起的SLDM算法。目前瞬變電磁三維正演，雖然已獲得長足發(fā)展，但還是處于起步階段，大量問題亟待解決?？偟膩碚f應(yīng)該在以下方面取得突破：數(shù)值算法的改進、計算效率的大幅提高、邊界條件的復(fù)雜化、以及如何消除地形對瞬變電磁場的影響。

瞬變電磁法；三維數(shù)值模擬；有限單元法；有限差分法；積分方程法

瞬變電磁法（Transient Electromagnetic Method，簡稱TEM）屬于時間域主動源（人工源）電磁法，利用布設(shè)在地表的不接地回線向地下發(fā)射一次脈沖磁場（磁源），或利用接地線源直接向地下發(fā)送一次脈沖電流場（電源），使地下導(dǎo)電介質(zhì)因電磁感應(yīng)現(xiàn)象而產(chǎn)生感應(yīng)渦旋電流，從而形成隨時間變化的二次電場和磁場；在一次場的脈沖間歇，測量二次場隨時間的衰減［1-2］。TEM觀測的是純二次場，可在近區(qū)測量，裝置靈活多樣，并且對低阻反應(yīng)靈敏、體積效應(yīng)小、分辨率高、操作簡單、工作效率高。所以近些年來，瞬變電磁法正日益凸顯出其在地質(zhì)勘探中的重要地位，且已被廣泛應(yīng)用于礦產(chǎn)資源勘探、水文地質(zhì)調(diào)查和工程勘查等領(lǐng)域［3-4］。

瞬變電磁三維數(shù)值模擬研究，可以幫助我們清楚地認(rèn)識瞬變電磁場在地下介質(zhì)內(nèi)的傳播規(guī)律，為時間域電磁法的深入研究提供了必要的理論基礎(chǔ)，同時對指導(dǎo)野外瞬變電磁勘探數(shù)據(jù)處理解釋具有重要的理論和現(xiàn)實意義。因此，對三維瞬變電磁響應(yīng)進行數(shù)值模擬研究很有必要。

1　瞬變電磁正演研究現(xiàn)狀

在瞬變電磁法應(yīng)用日益廣泛的同時，其反演解釋技術(shù)卻相對滯后。由于瞬變電磁場的正、反演計算需要耗費大量的計算機內(nèi)存和時間，個人PC機幾乎無法完成這項任務(wù)，嚴(yán)重限制了瞬變電磁法理論研究的進一步深入，使得瞬變電磁法長期停留在一維、2.5維正、反演研究階段。并且在相當(dāng)長的時間里，都著重于如何獲得電磁場問題的解析解，但事實上能夠用解析解來表征的電磁場問題非常有限。之后，雖然發(fā)展了一些數(shù)值計算方法，但是受當(dāng)時計算條件的限制，這些方法均無法充分發(fā)揮其作用［5］。進入21世紀(jì)后，計算機技術(shù)得到快速發(fā)展，使得我們可以進一步借助計算機的強大計算能力，來研究三維模型的瞬變電磁響應(yīng)，提高瞬變電磁法的應(yīng)用水平和解釋精度。

20世紀(jì)末期，隨著計算機技術(shù)的快速發(fā)展，瞬變電磁正演研究由早期的一維模型的解析方法，逐步轉(zhuǎn)向了高維數(shù)值計算方法。目前國內(nèi)外許多學(xué)者在從事二維與三維的數(shù)值計算，尤其是針對任意三維復(fù)雜模型的求解更是成為目前的研究熱點［6-7］。

2　三維瞬變電磁正演研究方法

20世紀(jì)80年代之后，時間域三維瞬變電磁法數(shù)值計算主要采用積分方程法、有限單元法、時域有限差分法和SLDM算法［8］。

（1）積分方程法。SanFilipo&Hohmann（1985）首次利用時域積分方程法計算了由矩形回線激發(fā)的均勻半空間中三維模型的瞬變電磁響應(yīng)。之后，Newman&Hohmann（1986）利用積分方程法計算了層狀大地的頻率域電磁響應(yīng)，然后通過余弦變換得到時間域響應(yīng)。Zhdanov（2006）給出了求解三維模型瞬變電磁響應(yīng)的積分方程法。Swidinsky&Edwards（2009）利用積分方程法模擬了由長接地導(dǎo)線激發(fā)的瞬變電磁響應(yīng)。Cox（2010）首次采用積分方程法對航空電磁數(shù)據(jù)進行三維數(shù)值模擬［9］。

（2）有限單元法。Kuo&Cho（1980）利用時域有限單元法對低阻礦體進行了數(shù)值模擬，取得了良好的效果。Pridmore&Hohmann（1981）采取有限單元法對三維電磁探測進行了求解。Sugeng（1998）使用矢量有限單元法對回線源在電性差異較大的地電模型中激發(fā)的瞬變電磁場進行了三維數(shù)值模擬。Stalnaker（2004）利用節(jié)點有限單元法對帶地形的磁性源瞬變電磁場響應(yīng)進行了三維數(shù)值模擬。Epov&Shurina（2007）利用矢量有限單元法對層狀介質(zhì)的瞬變電磁響應(yīng)進行了三維模擬。Borner（2008）提出一種基于傅氏變換和Krylov子空間的新型高效的有限元算法［11］。Um（2010）提出一種結(jié)合有限差分方法和矢量有限元方法的新算法，用來模擬海洋瞬變電磁響應(yīng)。李建慧（2011）年利用矢量有限元方法對中心回線瞬變電磁響應(yīng)進行了成功的數(shù)值模擬［8］。此外，譜有限元法、高階自適應(yīng)有限元法等也被應(yīng)用到了瞬變電磁三維數(shù)值模擬中［7］。

（3）時域有限差分法。Oristaglio&Hohmann（1984）提出基于Dufort-Frankel法的有限差分法，用來研究二維地質(zhì)體的瞬變響應(yīng)［11］。Adhidjaja&Hohmann（1989）試圖對Oristaglio&Hohmann的算法進行改進，以用于求解三維瞬變電磁問題，但是沒有成功。后來，Wang &Hohmann（1993）采用改進的Dufort-Frankel法給出了求解三維瞬變電磁場的時域有限差分算法，成功地將Oristaglio&Hohmann提出的二維算法移植到了三維問題的求解中來［12］。Endo&Noguchi（1999）以 Wang &Hohmann的算法為基礎(chǔ)，利用坐標(biāo)變換的方法解決了帶地形的回線源瞬變電磁三維正演問題。Commer&Newman（2004）針對電性源長偏移瞬變電磁問題，對Wang&Hohmann的算法進行改進，用數(shù)值模擬的方法計算關(guān)斷時間之前的電場值，以此代替解析解作為初值，并引入并行算法，使得有限差分法可以模擬復(fù)雜三維模型的瞬變電磁響應(yīng)。史紅蓓（2009）對上述方法進行了改進，利用軟源施加的方式成功引進初始場，克服了Wang&Hohmann方法以解析解來充當(dāng)初始場的局限性［13］。孫懷鳳（2013）考慮了關(guān)斷時間，得到了全時域的瞬變電磁響應(yīng)［7］。

（4）SLDM算法。瞬變場進行數(shù)值模擬時，在對計算區(qū)域離散網(wǎng)格剖分之后，需要利用有限差分法來求解一個關(guān)于Maxwell方程組的大型常微分方程。無論是顯式離散，還是隱式離散，均需要進行時間步迭代，使得效率偏低。因此Druskin（1999）基于Krylov子空間投影建立了針對各向異性三維模型的譜Lanczos迭代分解法（SLDM）。后來諸多學(xué)者對SLDM算法做出了改進，其中B?rner（2008）基于頻率域有限元離散和Krylov子空間投影技術(shù)建立了模型簡化方法；Zaslavsky （2011）結(jié)合最優(yōu)子空間投影技術(shù)給出了一般化的SLDM方法，并考慮了IP效應(yīng)［14］。

3　目前研究中亟需解決的問題

目前瞬變電磁三維正演，雖然已獲得長足發(fā)展，但總的來說還是處于起步階段，大量問題亟待解決。

3.1滿足復(fù)雜3D模型的亞網(wǎng)格剖分技術(shù)

對于復(fù)雜3D模型，均勻網(wǎng)格剖分技術(shù)已不能很好地反映模型邊界處的場值變化。為了精細(xì)地反映場的變化，需要加密網(wǎng)格，利用非均勻網(wǎng)格剖分方案來進行模型離散。然而目前普遍采用的時域有限差分法，計算過程中必須滿足CFL穩(wěn)定性條件，使得迭代時間步不能太大，大大降低了迭代效率。近年來快速發(fā)展的交錯逆網(wǎng)格時域有限差分法（ADI-FDTD），其時間步長不需要滿足CFL條件，可根據(jù)精度要求適當(dāng)增大時間步間隔以提高迭代效率。

3.2計算效率的提高

瞬變電磁數(shù)值模擬計算量大，耗時長，要想其快速發(fā)展，必須進一步提高計算效率，可以考慮使用CPU+ GPU異構(gòu)并行算法，不僅可以提高程序的運算效率，還能解決內(nèi)存需求大的問題。

3.3吸收邊界條件

另外，要施加合適的邊界條件，以使得復(fù)雜三維模型在網(wǎng)格截斷處的電磁波反射最小化，比較好的解決方法就是設(shè)立Mur吸收邊界條件和完全匹配層吸收邊界。

3.4地形對瞬變電磁響應(yīng)的影響

最后需要對帶地形的瞬變電磁響應(yīng)進行研究，因為實際工作中，瞬變電磁測量都是在地形起伏的地區(qū)進行測量，我們必須深入了解地形對瞬變電磁的影響，才能在數(shù)據(jù)處理時進行地形校正。

［1］樸化榮.電磁測深法原理［M］.北京：地質(zhì)出版社，1990.

［2］牛之鏈.時間域電磁法原理［M］.長沙：中南大學(xué)出版社，2007.

［3］蔣邦遠(yuǎn).實用近區(qū)磁源瞬變電磁法勘探［M］.北京：地質(zhì)出版社，1998.

［4］李貅.瞬變電磁測深的理論與應(yīng)用［M］.西安：陜西科學(xué)技術(shù)出版社，2002.

［5］王長清.瞬變電磁場——理論和計算［M］.北京：北京大學(xué)出版社，2011.

［6］薛國強，李貅，底青云.瞬變電磁法正反演問題研究進展［J］.地球物理學(xué)進展，2008，23（4）：1165-1172.

［7］孫懷鳳，李貅，李術(shù)才，等.考慮關(guān)斷時間的回線源激發(fā)TEM三維時域有限差分正演［J］.地球物理學(xué)報，2013，56（3）：1049-1064.

［8］李建慧.瞬變電磁法正演計算進展［J］.地球物理學(xué)進展，2012，27（4）：1393-1400.

［9］Cox L H，Zhdanov M S.Advanced Computational Methods of Rapid and Rigorous 3D Inversion of Airborne Electromagnetic Data［J］.Communications in Computational Physics，2008，3（1）：160-179.

［10］B?rner R，Ernst O G，Spitzer K.Fast 3-D Simulation of Transient Electromagnetic Fields by Model Reduction in the Frequency Domain Using Krylov Subspace Projection［J］.Geophysical Journal International，2008，173（3）：766-780.

［11］ Oristaglio M L，Hohmann G W.Diffusion of Electromagnetic Fields into a Two-dimensional Earth：A Finite-difference Approach［J］.Geophysics，1984，49（7）：870-894.

［12］Wang.T，Hohmann.G.W.A finite-difference，Time-domain Solution for Three-dimensional Electromagnetic Modeling［J］Geophysics，1993（58）：797-809.

［13］史紅蓓.導(dǎo)電媒質(zhì)中三維似穩(wěn)瞬態(tài)電磁場的FDTD數(shù)值模擬［D］.江蘇大學(xué)，2009.

［14］ Zaslavsky M，Druskin V，Knizhnerman L.Solution of 3D Time-domain Electromagnetic Problems Using Optimal Subspace Projection［J］.Geophysics，2011，76（6）：F339-F351.

P59

A

1004-5716（2016）03-0179-03

2015-12-29

劉堯（1989-），男（漢族），陜西渭南人，現(xiàn)從事地球物理正反演研究及勘探工作。