李展輝，黃清華

北京大學(xué)地球與空間科學(xué)學(xué)院地球物理學(xué)系，北京 100871

1 引言

時域有限差分（Finite difference time domain，F(xiàn)DTD）是瞬變電磁法（Transient electromagnetic，TEM）正演常用方法之一.Goldman和Stoyer（1983）首次使用隱式FDTD模擬三維軸對稱模型，隨后Oristaglio和Hohmann（1984）使用修改的Du Fort-Frankel格式模擬了TEM二維模型，Wang和Hohmann（1993）在Du Fort-Frankel格式的基礎(chǔ)上發(fā)展了一種TEM三維FDTD正演方法，為后續(xù)部分學(xué)者借用（孫懷鳳等，2013；許洋鋮等，2012；Commer and Newman，2004；Endo and Noguchi，2002）.FDTD對TEM的三維正演是在有限模型空間內(nèi)模擬能擴散到更大區(qū)域的瞬變場，因而需要引入截斷邊界條件.正演中通常采用狄利克雷邊界條件（Dirichlet boundary condition，DBC）（Wang and Hohmann，1993），并設(shè)置邊界足夠遠.DBC在擴散場未抵達或剛抵達邊界時效果較好，但隨著延遲時間的推移，反射場將對原擴散場產(chǎn)生較大的干擾.要消除這種干擾，則需要更大的模型，這將增加內(nèi)存和計算量.因而一些學(xué)者開始研究其他的吸收邊界來代替DBC，比如Mur吸收邊界（秦臻和胡文寶，2004；姬金祖，2008；岳建華，2007）以及在其基礎(chǔ)上修改而來的新吸收邊界條件（楊海燕和岳建華，2009）.但是他們并沒有對吸收邊界的效果作詳細的討論，也沒有給出與理論解的對比.李墩柱（2010）將電磁波正演中應(yīng)用最普遍的完全匹配層（Perfectly Matched Layer，PML）吸收邊界（Berenger，1994，1996；李展輝等，2009；Chew and Weedon，1994）應(yīng)用到TEM半空間模型正演中，然而文中模型的水平范圍為［-3000m，3000m］，半空間電導(dǎo)率為0.01S·m-1，但是最長延遲時間僅為8ms，渦流電場形成的“煙圈”還遠未到達邊界，晚期吸收效果并不明確，在參數(shù)設(shè)置上也僅僅說明需要多次嘗試.本文深入研究了PML在TEM正演中的吸收效果，并最終采用了其中的一種——復(fù)頻率參數(shù)完全匹配層（Complex Frequency Shifted Perfectly Matched Layer，CFS-PML）（Kuzuoglu and Mittra，1996；Roden and Gedney，2000；姜永金等，2004；Berenger，2012）.首先通過推導(dǎo)擴散場在常規(guī)PML內(nèi)部的平面波解的形式，指出了常規(guī)PML在TEM正演中失效的原因，然后給出了擴散場在CFS-PML內(nèi)部的平面波解，研究了CFS-PML各項參數(shù)對吸收效果的影響，并給出了CFS-PML的參數(shù)選擇方案.最后將CFSPML應(yīng)用到全空間和半空間模型中以驗證CFS-PML的有效性.同時為了提高半空間模型的正演效率，本文采用了李展輝和黃清華（2011）提出的新方法，而全空間模型正演則采用普遍使用的Du Fort-Frankel格式.

2 理論推導(dǎo)

2.1 含PML的Maxwell頻域方程組

頻率域帶PML的旋度方程可以表示為（Chew and Weedon，1994）

其中

常規(guī)PML中（Chew and Weedon，1994），

其中κu為網(wǎng)格延拓因子，由于本文采用非等距網(wǎng)格，κu可直接包含于網(wǎng)格中.因此在下文中所有的κu值都設(shè)置為1.令ε′＝εΔ-jσ／ω，可得

2.2 常規(guī)PML和CFS-PML的吸收性質(zhì)

設(shè)式（5）的平面波解為

定義復(fù)波速

可得到頻散方程

對于常規(guī)PML，式（6）可寫為

其中

而TEM的低頻近似滿足

那么復(fù)波速可以近似為

將式（13）代入到式（10）中，并且考慮在x方向吸收，設(shè)置 （σx，σy，σz）＝ （σx，0，0），可得

設(shè)

為PML的吸收系數(shù)，APML越小，代表吸收越強烈.從式（15）中可以看出，隨著角頻率ω的減小，APML也會隨著減小，而當ω→0時，亦有APML→0.這意味著PML對低頻場的吸收隨著頻率的降低而更為強烈，過強的吸收在數(shù)值計算中會引起大的異常（Berenger，2000，1996）.

對于CFS-PML，經(jīng)過類似的推導(dǎo)之后可以得到

從式（16）可以看出，αx的出現(xiàn)使得CFS-PML避免了在ω接近0時出現(xiàn)吸收非常強烈的現(xiàn)象.

設(shè)置CFS-PML的一個參考頻率：

把式（17）代入式（16），得

當f?fαx時，式（18）近似為式（14），CFS-PML退化成常規(guī)PML；當f?fαx時，式（18）近似為

這僅僅是x方向的網(wǎng)格延拓.

瞬變場的優(yōu)勢頻率會隨著時間逐漸變低.當優(yōu)勢頻率足夠低并進入PML內(nèi)部時，PML會對這些優(yōu)勢頻率產(chǎn)生過強的吸收進而使整個瞬變場產(chǎn)生異常.而CFS-PML能避免這種情況，但是σx／αx這個比例非常重要，從式（19）可以看出，對于f?fαx頻率部分只是一個網(wǎng)格延拓因子，如果過大同樣會引起數(shù)值異常.為此我們可以設(shè)置另一個參考頻率：

2.3 含CFS-PML的Maxwell時域方程組

時域方程組可以用卷積的形式表達出來（Roden and Gedney，2000），以Ex為例：

其中ˉsy（t）和ˉsz（t）分別是1／sy（ω）和1／sz（ω）的拉普拉斯逆變換：

其中δ（t）是狄拉克函數(shù)，ζu（t）可以表述如下：

其中u（t）為單位階躍函數(shù).那么最終式可以表示為

其中

其他分量可以類似地表達出來.

2.4 CFS-PML參數(shù)選取研究

其中L為PML的厚度.假設(shè)場沿x方向傳播，X＝x.定義CFS-PML吸收系數(shù)為

式中

式（28）中存在兩個可控參數(shù)：fαx和L.吸收系數(shù)A越小代表額外吸收越強烈.很明顯，吸收系數(shù)A是CFS-PML厚度L的遞減函數(shù)，L越大，吸收越強烈.下面將對參數(shù)fαx選取作具體的分析.

假設(shè)介質(zhì)是均勻的，且電導(dǎo)率為σ＝1×10-3S·m-1，場源到PML內(nèi)邊界的距離為L1＝1400m，PML厚度為L2＝600m，關(guān)注的最大時間為tmax＝1ms，那么抵達邊界并穿過PML返回場內(nèi)的最低頻率為

圖1 不同參考頻率下CFS-PML對擴散場的額外吸收系數(shù)比較Fig.1 The additional absorbing coefficient of CFS-PML under different reference frequencies

3 數(shù)值模擬結(jié)果及討論

本文除巷道帶異常體模型外，所有正演均采用相同的網(wǎng)格：網(wǎng)格數(shù)量為128×128×128，中央網(wǎng)格大小為Δxmin＝Δymin＝Δzmin＝10m，從中央到四周逐漸增加網(wǎng)格步長，最大處為Δxmax＝10Δxmin，Δymax＝10Δymin，Δzmax＝7Δzmin.圖2描述了網(wǎng)格的主要特征.CFS-PML的層數(shù)為8層，已經(jīng)包含在上述網(wǎng)格中.

圖2 瞬變電磁三維正演網(wǎng)格示意圖Fig.2 A sketch map of the mesh for TEM 3Dforward modeling

3.1 均勻全空間

為了方便與理論解對比，我們采用了均勻全空間模型.同時為了探討CFS-PML對不同介質(zhì)的適應(yīng)性，我們選取了兩種不同的介質(zhì)：σ＝1×10-2S·m-1和σ＝1×10-3S·m-1.以水平正方形線框作為發(fā)射線框，線框邊長L＝70m，垂直方向位于z＝0處，水平方向位于x-y平面的正中央.采用脈沖形式的發(fā)射電流，最終比對垂直磁場Hz.

選擇的測線為z＝0和y＝0平面的交線.我們定義在CFS-PML與正常區(qū)域交接處Hz易號的時間點為早期、晚期的分界點.首先對σ＝1×10-2S·m-1的三個早期時間點進行對比，根據(jù)牛之鏈（2007）中式（2-33）計算得到早晚期分界延遲時間約為4.5ms，那么首先選取t＝0.35ms，1.23ms和2.64ms三個時間點.圖3為在上述三個時間點測線上分別使用CFS-PML、DBC得到的Hz值與理論值的對比.可以看出，在早期三者幾乎一致.這是因為擴散場尚未或者剛剛抵達邊界，邊界條件的區(qū)別未能體現(xiàn)出來.因此在下文中，本文將著重擴散晚期的對比.

圖4為t＝6.80ms，12.50ms，和20.00ms時分別使用CFS-PML，DBC得到的Hz值和理論解的對比.可以看出t＝6.8ms時使用DBC的結(jié)果在邊界處已經(jīng)出現(xiàn)了較大的異常，隨后這種異常開始擴散到場的中央，結(jié)果越來越偏離理論解.而使用CFS-PML吸收邊界時，不管在中央還是在邊界，Hz值依然和理論解符合得很好.CFS-PML的吸收作用在晚期非常顯著.

圖5為σ＝1×10-3S·m-1情況下不同時刻使用CFS-PML、DBC得到的Hz值和理論解之間的對比，同時加入了CFS-PML中設(shè)置fαx＝∞所得的結(jié)果.σ＝1×10-3S·m-1模型中場的擴散速度要明顯快于σ＝1×10-2S·m-1的情況，但是兩種情況下的場在一定的條件下具有相似性（王華軍，2008）.根據(jù)這種相似性，t＝2ms時場的形態(tài)類似于σ＝1×10-2S·m-1時t＝20ms情況下場的形態(tài)，因此本文將σ＝1×10-3S·m-1情況下最晚的一個時間點設(shè)定為t＝2ms.從圖5可以看出，在晚期CFS-PML的優(yōu)勢依然非常明顯，而且經(jīng)過適當設(shè)置后的CFS-PML也明顯優(yōu)于設(shè)置fαx＝∞的CFS-PML.

圖3 全空間σ＝1×10-2S·m-1模型中，t＝0.35ms，1.23ms和2.64ms時使用CFS-PML吸收邊界、DBC的Hz分量和理論解的對比Fig.3 The comparison among the Hzcomponents obtained from CFS-PML，DBC，and the analytical solutions at the times of 0.35ms，1.23ms，and 2.64ms of a whole space model withσ＝1×10-2S·m-1

圖4 全空間σ＝1×10-2S·m-1模型中，t＝6.80ms，12.50ms和20.00ms時使用CFS-PML吸收邊界條件、DBC的Hz分量和理論解的對比Fig.4 The comparison among the Hzcomponents obtained from CFS-PML，DBC，and the analytical solutions at the times of 6.80ms，12.50ms，and 20.00msof a whole space model withσ＝1×1 0-2 S·m-1

圖5 全空間σ＝1×10-3S·m-1模型中，t＝0.680ms和2.00ms時使用CFS-PML吸收邊界條件、DBC的Hz分量和理論解的對比，對比中還加入了在CFS-PML中設(shè)置fαx＝ ∞ 所得的Hz分量Fig.5 The comparison among the Hzcomponents obtained from CFS-PML，DBC，and the analytical solutions at the times of 0.68ms and 2.00ms of a whole space model withσ＝1×10-3S·m-1.The Hz component obtained from CFS-PML with fαx＝ ∞is also included for reference

為顯示CFS-PML內(nèi)部場的情況，我們將σ＝1×10-3S·m-1模型中2ms時CFS-PML內(nèi)外的場都描繪出來并與理論解進行對比，如圖6所示.可以看出，CFS-PML內(nèi)部的場與理論解截然不同.CFS-PML內(nèi)部出現(xiàn)了Hz的正負變化，意味著渦流的存在，而理論解并沒有Hz符號的變化.CFS-PML雖然不能將擴散場完全吸收，但是其在吸收的同時以某種方式將渦流限制在了CFS-PML內(nèi)部，使得正常區(qū)域的擴散場不受影響.

3.2 均勻半空間

本文采用了李展輝和黃清華（2011）介紹的方法在半空間模型中引入空氣層.由于空氣層中波速遠大于地下場的擴散速度，由式（10）可知，針對地下介質(zhì)設(shè)計的CFS-PML在空氣層中對波的吸收基本為零，因此在空氣層中需要在CFS-PML介質(zhì)內(nèi)設(shè)置一定的電導(dǎo)率進行吸收.半空間模型中，地下介質(zhì)的電導(dǎo)率設(shè)置為σ＝1×10-2S·m-1，空氣層中電導(dǎo)率設(shè)置為1×10-6S·m-1，并在靠近邊界處逐漸增加.地空邊界設(shè)置在z＝0處，激勵源的位置與形式均與全空間模型一致.

圖7顯示了使用CFS-PML吸收邊界、DBC所得到的Hz分量和理論解之間的對比.從圖中可以看出，t＝6.80ms時使用DBC的Hz分量在兩端已經(jīng)偏離了理論解.隨著時間的推移，誤差越來越大，在20ms時，誤差已經(jīng)在一個量級以上.而使用CFS-PML的Hz與理論解的吻合度大幅優(yōu)于使用DBC的情況，但是也有一定的誤差，而且這種誤差也隨時間的推移緩慢地增加.這可能是CFS-PML不能對空氣層中的場產(chǎn)生有效的吸收造成的，越到晚期越明顯.圖8詳細展示了線圈中心Hz相對誤差隨時間變化.從圖中可以看出，線圈中心相對誤差隨著時間逐漸增加，在20ms時誤差達到了13%.實際應(yīng)用中需要根據(jù)精度要求調(diào)整模型的大小以控制相對誤差.

圖6 全空間σ＝1×10-3S·m-1模型CFS-PML內(nèi)外Hz分量展示，并與理論解比較以體現(xiàn)CFS-PML內(nèi)部場特征Fig.6 The Hzcomponent in and out the CFS-PML of a whole space model withσ＝1×10-3S·m-1 at 2ms

3.3 巷道模型

全空間瞬變電磁法典型的應(yīng)用之一為煤礦巷道內(nèi)地下水的探測（李宇等，2012；于景 等，2011；郭純等，2006）.模型如圖9所示.地下水所在的位置電導(dǎo)率較高，因此模型中設(shè)定圍巖電導(dǎo)率為0.05S·m-1，異常體電導(dǎo)率為0.5S·m-1.巷道內(nèi)通常很狹小，因此發(fā)射線圈設(shè)定為3m×3m，并采用重疊回線的方式記錄感應(yīng)電壓.由于該模型不存在理論解，為了獲取類似理論解的參考解，我們將原網(wǎng)格x，y方向各增加60個網(wǎng)格，在z方向增加40個網(wǎng)格.總網(wǎng)格數(shù)量達到188×188×168，以保證在我們關(guān)心的時間窗內(nèi)，邊界反射場沒有進入觀測區(qū)域.考慮到模型尺寸問題，最大的延遲時間設(shè)置為0.6ms，因為根據(jù)擴散深度公式（牛之璉，2007）

圖7 半空間模型中，t＝6.80ms，12.50ms和20.00ms時使用CFS-PML吸收邊界條件、DBC的Hz分量和理論解的對比Fig.7 The comparisons among the Hzcomponents obtained from CFS-PML，DBC，and the analytical solutions at the times of t＝6.80ms，12.50ms，and20.00ms of a half space model

圖8 回線中心Hz的相對誤差隨時間分布圖Fig.8 The relative errors of Hzvarying over time at the center of the rectangular loop

計算可得此時的擴散深度為124m，已經(jīng)超出了本模型范圍.式（31）中ρ為電阻率，t為延遲時間.結(jié)果如圖10所示.從圖中可以看出CFS-PML的結(jié)果與參考解一直符合得很好，最終的相對誤差不超過6%.

圖9 巷道內(nèi)含異常體模型，模型關(guān)于y＝0對稱Fig.9 Underground roadway model with an anomalous body，which is symmetric about the y＝0plane

圖10 巷道模型的重疊回線測量垂直感應(yīng)電動勢結(jié)果圖，左下角的小圖為相對于參考解的相對誤差Fig.10 Vertical emf of the underground roadway model.The relative errors to the reference solution are shown in the embedded figure

4 結(jié)論

本文將CFS-PML應(yīng)用到TEM三維FDTD正演中，通過適當?shù)卦O(shè)置CFS-PML參數(shù)保證了數(shù)值計算的穩(wěn)定性和良好的吸收效果.全空間和半空間數(shù)值結(jié)果表明CFS-PML在瞬變場早期與傳統(tǒng)的狄利克雷邊界條件的結(jié)果一致，晚期則明顯優(yōu)于狄利克雷邊界條件的結(jié)果，并且能有效地節(jié)省計算時間和存儲空間.但是針對地下擴散場設(shè)計的CFSPML并不能有效地吸收空氣中的電磁場，導(dǎo)致半空間模型結(jié)果到后期也產(chǎn)生一定的誤差，并且隨著時間緩慢的增長.實際應(yīng)用中需要根據(jù)TEM正演的最大延遲時間以及對誤差的容忍度適當調(diào)整模型的大小，以確保誤差在要求之內(nèi).本文最后將CFSPML應(yīng)用到與實際情況貼近的巷道模型中，采用重疊回線測量感應(yīng)電動勢，并與參考解進行對比和誤差分析，再一次證明了CFS-PML的有效性.

Bérenger J P.1994.A perfectly matched layer for the absorption of electromagnetic waves.Journal of Computational Physics，114（2）：185-200.

Bérenger J P.1996.Three-dimensional perfectly matched layer for the absorption of electromagnetic waves.Journal of Computational Physics，127（2）：363-379.

Bérenger J P.1996.Perfectly matched layer for the FDTD solution of wave-structure interaction problems.IEEE Transactions on Antennas and Propagation，44（1）：110-117.

Bérenger J P.2000.Numerical reflection of evanescent waves by PMLs：origin and interpretation in the FDTD case.Expected consequences to other finite methods.Int.J.Numer.Model，13（2-3）：103-114.

Bérenger J P.2012.An optimized CFS-PML for wave-structure interaction problems.IEEE Transactions on Electromagnetic Compatibility，54（2）：351-358.

Chew W C，Weedon W H.1994.A 3Dperfectly matched medium from modified Maxwell′s equations with stretched coordinates.Microwave and Optical Technology Letters，7（13）：599-604.

Commer M，Newman G.2004.A parallel finite-difference approach for 3Dtransient electromagnetic modeling with galvanic sources.Geophysics，69（5）：1192-1202.

Endo M，Noguchi K.2002.Three-dimensional modeling considering the topography for the case of the time-domain electromagnetic method.Amsterdam：Elsevier Science BV，35：85-107.

Goldman M M，Stoyer C H.1983.Finite-difference calculations of the transient field of an axially-symmetric earth for vertical magnetic dipole excitation.Geophysics，48（7）：953-963.

Guo C，Liu B Z，Bai D H.2006.Prediction of water disasters ahead of tunneling in coal mine using continuous detection by UWTEM.Seismology and Geology（in Chinese），28（3）：456-462.

Ji J Z，Gao X，Liu Z H.2008.Simplified format application of 2D two-rank Mur absorbing boundary condition.Journal of Shenyang Institute of Aeronautical Engineering（in Chinese），25（2）：10-13.

Jiang Y J，Mao J J，Cai S L.2004.Implementation and analysis of the perfectly matched layer media with CFS for the PSTD method.Journal of Microwaves（in Chinese），20（4）：36-39.

Kuzuoglu M，Mittra R.1996.Frequency dependence of the constitutive parameters of causal perfectly matched anisotropic absorbers.IEEE Microwave and Guided Wave Letters，6（12）：447-449.

Li D Z.2010.Numerical Simulation of Transient Electromagnetic Fields（in Chinese）［Master′s thesis］.Beijing：Peking University.

Li Y，Tong B，Guo C F，etal.2012.Study on application and denoising of mine full space transient electromagnetic method.Coal Engineering（in Chinese），（11）：29-31.

Li Z H，Huang Q H，Wang Y B.2009.A 3-D staggered grid PSTD method for borehole radar simulations in dispersive media.Chinese J.Geophys.（in Chinese），52（7）：1915-1922.

Li Z H，Huang Q H.2011.A 3Dcurvilinear staggered-grid finite difference time-domain method with an impulse source for transient electromagnetic method.Chinese Geophysics（in Chinese），240.

Niu Z L.2007.Principle of Time Domain Electromagnetic Method（in Chinese）.Changsha：Central South University Press.

Oristaglio M L，Hohmann G W.1984.Diffusion of electromagnetic fields into a two-dimensional earth-a finite-difference approach.Geophysics，49（7）：870-894.

Qin Z，Hu W B.2004.Modified Mur absorbing boundary conditions for lossy medium FDTD calculation.Journal of Jianghan Petroleum Institute（in Chinese），26（2）：76-77.

Roden J A，Gedney S D.2000.Convolution PML（CPML）：An efficient FDTD implementation of the CFS-PML for arbitrary media.Microwave and Optical Technology Letters，27（5）：334-339.

Sun H F，Li X，Li S C，etal.2013.Three-dimensional FDTD modeling of TEM excited by a loop source considering ramp time.Chinese J.Geophys.（in Chinese），56（3）：1049-1064.

Wang H J.2008.Time domain transient electromagnetism all time apparent resistivity translation algorithm.Chinese J.Geophys.（in Chinese），51（6）：1936-1942.

Wang T，Hohmann G W.1993.A finite-difference，time-domain solution for 3-dimensional electromagnetic modeling.Geophysics，58（6）：797-809.

Xu Y C，Lin J，Li S Y，etal.2012.Calculation of full-waveform airborne electromagnetic response with three-dimension finitedifference solution in time-domain.Chinese J.Geophys.（in Chinese），55（6）：2105-2114.

Yang H Y，Yue J H.2009.Application of absorbing boundary condition in whole-space computation of transient electromagnetic response.Journal of China University of Mining ＆Technology（in Chinese），38（2）：263-268.

Yu J C，Liu Z Q，Liao J J，etal.2011.Application of full space transient electromagnetic method to mine water prevention and control.Coal Science and Technology（in Chinese），39（9）：110-113.

Yue J H，Yang H Y，Hu B.2007.3Dfinite difference time domain numerical simulation for TEM in-mine.Progress in Geophysics（in Chinese），22（6）：1904-1909.

附中文參考文獻

郭純，劉白宙，白登海.2006.地下全空間瞬變電磁技術(shù)在煤礦巷道掘進頭的連續(xù)跟蹤超前探測.地震地質(zhì)，28（3）：456-462.

姬金祖，高旭，劉戰(zhàn)合.2008.二維Mur二階吸收邊界條件簡化格式的應(yīng)用.沈陽航空工業(yè)學(xué)院學(xué)報，25（2）：10-13.

姜永金，毛鈞杰，柴舜連.2004.CFS-PML邊界條件在PSTD算法中的實現(xiàn)與性能分析.微波學(xué)報，20（4）：36-39.

李墩柱.2010.地形起伏下的瞬變電磁場數(shù)值模擬［碩士論文］.北京：北京大學(xué).

李宇，童兵，郭昌放等.2012.礦井全空間瞬變電磁法應(yīng)用及其消噪研究.煤炭工程，（11）：29-31.

李展輝，黃清華，王彥賓.2009.三維錯格時域偽譜法在頻散介質(zhì)井中雷達模擬中的應(yīng)用.地球物理學(xué)報，52（7）：1915-1922.

李展輝，黃清華.2011.瞬變電磁正演之帶激勵源的三維曲線交錯網(wǎng)格時域有限差分方法.中國地球物理，240.

牛之璉.2007.時間域電磁法原理.長沙：中南大學(xué)出版社.

秦臻，胡文寶.2004.用于有耗介質(zhì)FDTD計算的改進的Mur吸收邊界條件.江漢石油學(xué)院學(xué)報，26（2）：76-77.

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

王華軍.2008.時間域瞬變電磁法全區(qū)視電阻率的平移算法.地球物理學(xué)報，51（6）：1936-1942.

許洋鋮，林君，李肅義等.2012.全波形時間域航空電磁響應(yīng)三維有限差分數(shù)值計算.地球物理學(xué)報，56（6）：2105-2114.

楊海燕，岳建華.2009.吸收邊界條件在全空間瞬變電磁計算中的應(yīng)用.中國礦業(yè)大學(xué)學(xué)報，38（2）：263-268.

于景邨，劉振慶，廖俊杰等.2011.全空間瞬變電磁法在煤礦防治水中的應(yīng)用.煤炭科學(xué)技術(shù)，39（9）：110-113.

岳建華，楊海燕，胡搏.2007.礦井瞬變電磁法三維時域有限差分數(shù)值模擬.地球物理學(xué)進展，22（6）：1904-1909.