邱稚鵬，李展輝，李墩柱，黃清華＊

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

2 美國加州理工學(xué)院地球與行星科學(xué)系，美國加州帕薩蒂納 91125

1 引 言

在野外瞬變電磁（Transient electromagnetics，TEM）的測(cè)量過程中，測(cè)區(qū)的地形不可避免地存在起伏，有時(shí)甚至起伏很劇烈.這就導(dǎo)致實(shí)際的測(cè)量數(shù)據(jù)除包含地下結(jié)構(gòu)的有效信息之外，還會(huì)摻雜著地形的影響.考察地形對(duì)瞬變電磁場(chǎng)影響的正演模擬對(duì)實(shí)際測(cè)量、反演方法的研究以及測(cè)量數(shù)據(jù)的解釋有著重要的指導(dǎo)意義.

現(xiàn)有的關(guān)于瞬變電磁正演的研究包括有限元法（Finite－element method，F(xiàn)EM）［1－5］、積 分 方 程 法（Integral equations method，IEM）［6－10］、時(shí)域有限差分（Finite difference time domain，F(xiàn)DTD）方法［11－22］等，其中，帶地形瞬變電磁場(chǎng)正演模擬的研究相對(duì)較少［1－2，6－9，12－14］.FEM 進(jìn)行帶地形 TEM 三維正演 時(shí) 可以很好地刻畫地形，但其求解的過程需要解大型矩陣，耗費(fèi)大量時(shí)間.IEM的主要計(jì)算量與異常體或地形的規(guī)模、形狀有關(guān)，當(dāng)?shù)匦屋^簡(jiǎn)單時(shí)，IEM的計(jì)算量很小，可以高效地進(jìn)行計(jì)算.但若要計(jì)算較為復(fù)雜的地形時(shí)，IEM會(huì)明顯受到限制.

Hohmann與Wang［15］提出了一種基于修正的Du Fort－Frankel格式［18］的有限差分法，該方法可以顯式地得到任意時(shí)間點(diǎn)的電磁場(chǎng)分量；但他們沒有考慮如何計(jì)算帶地形情況下的電磁場(chǎng).Endo與Noguchi［12］提出了垂向網(wǎng)格變換算法，通過將帶地形的物理域的場(chǎng)映射到平坦的計(jì)算域，實(shí)現(xiàn)了對(duì)地形的近似；該算法在地表采用向上延拓作為邊界條件，但在地形起伏較大時(shí)并不能獲得準(zhǔn)確的邊界條件；該文還計(jì)算了簡(jiǎn)單地形起伏的TEM場(chǎng)，但計(jì)算區(qū)域比較小，模型的設(shè)計(jì)僅限于單一模型.Mitsuhata［1］利用有限元法分析了二維長(zhǎng)偏移距瞬變電磁法（Long offset TEM，LOTEM）在電偶源激發(fā)的情況下地形的影響，認(rèn)為對(duì)于遠(yuǎn)離地形的測(cè)點(diǎn)，地形的影響不顯著.唐新功等［6－9］利用積分方程法進(jìn)行的LOTEM數(shù)值模擬表明，在電偶源情況下，地形的影響是一種“本地”的影響，并建議實(shí)際測(cè)量時(shí)選擇較平坦或開闊的地形［8］；但該方法對(duì)地形刻畫采用的是階梯近似，所以其地形刻畫比較粗糙.李墩柱［14］將非正交網(wǎng)格算法［19－20］引入到瞬變電磁場(chǎng)中，在地表引入了空氣層，避免了地形起伏給向上延拓帶來的問題；該算法在由逆變分量計(jì)算協(xié)變分量時(shí)，采用某點(diǎn)周圍四個(gè)點(diǎn)的等權(quán)投影的平均來近似該點(diǎn)的對(duì)應(yīng)分量，當(dāng)網(wǎng)格的劃分存在一定梯度時(shí)，該近似可能帶來一定的誤差.本文針對(duì)李墩柱算法中存在的上述問題，提出了一種采用非等權(quán)投影的改進(jìn)算法，在一定程度上提高了該算法的精度.

總體上講，以往針對(duì)地形影響的研究大都局限于對(duì)山峰、山谷、斜坡等地形類型的簡(jiǎn)單分析，對(duì)決定地形影響效果的各種因素的討論并不充分.為此，本文利用我們改進(jìn)的非正交網(wǎng)格算法開展了瞬變電磁場(chǎng)的數(shù)值模擬研究，探討了地形的深度、寬度和相對(duì)源的距離等參數(shù)對(duì)瞬變電磁場(chǎng)的影響，獲得了一些初步認(rèn)識(shí).

2 非正交網(wǎng)格FDTD

2.1 非正交網(wǎng)格

空間中的任意矢量都可以按照非正交的一般坐標(biāo)系（圖1）進(jìn)行展開.若選取該坐標(biāo)系的三分量作為基矢量，則可以得到任意矢量的展開式（以電場(chǎng)為例）：

圖1 一般坐標(biāo)系逆變基矢量與協(xié)變基矢量［20］Fig.1 Covariant vectors and contravariant vector in a general coordinate system［20］

其中，

協(xié)變基與逆變基有如下關(guān)系：

gij與gij稱為度量張量的分量.

同樣的，任意矢量也可對(duì)逆變基矢量進(jìn)行展開.將協(xié)變基定義為非正交網(wǎng)格的棱，則該坐標(biāo)系中線積分與面積分可分別表示如下：

2.2 非正交網(wǎng)格FDTD

對(duì)于麥克斯韋方程組，散度方程可以利用邊界條件自然滿足，實(shí)際計(jì)算中主要是對(duì)兩個(gè)旋度方程進(jìn)行計(jì)算.網(wǎng)格的劃分采用Yee網(wǎng)格［21］，電場(chǎng)定義在網(wǎng)格的棱心和整數(shù)時(shí)間點(diǎn)上，磁場(chǎng)各分量定義在

在非正交坐標(biāo)系中，相應(yīng)于上述協(xié)變基矢量Ai，存在一套相應(yīng)的逆基，稱為逆變基矢量Ai，其定義如下：網(wǎng)格的面心和半整數(shù)時(shí)間點(diǎn)，考察積分形式的安培定律和法拉第定律，在各向同性的情況下可以表示為：

將式中的線積分和面積分利用式（5），（6）表示.對(duì)時(shí)間的一階導(dǎo)數(shù)采用中心差分近似，可以得到

其中，Ji，j，k，JHi，j，k分別表示電場(chǎng)和磁場(chǎng)相應(yīng)位置協(xié)變基矢量圍成的晶格體積.最后，為完成電場(chǎng)與磁場(chǎng)的迭代，需要將等式左邊的電場(chǎng)與磁場(chǎng)的逆變分量轉(zhuǎn)換為協(xié)變分量.考慮到Y(jié)ee網(wǎng)格的任意點(diǎn)只定義了一個(gè)分量，相應(yīng)的另外兩個(gè)分量需要由該點(diǎn)周圍點(diǎn)的分量平均近似得到［14］：

其中g(shù)′ii＝Ai·AHi.Ai，AHi分別表示電場(chǎng)與磁場(chǎng)相應(yīng)位置的基矢量.非正交網(wǎng)格的穩(wěn)定性條件為［20］：

其中，sup代表上限（upper－bound）.

同時(shí)，因?yàn)樗沧冸姶艌?chǎng)近似為擴(kuò)散場(chǎng)，為了保證方程中的波動(dòng)不會(huì)影響到擴(kuò)散項(xiàng)，還需要保證最大時(shí)間步長(zhǎng)［15］：

其中α的選取一般小于0.2.

2.3 改進(jìn)的非正交網(wǎng)格FDTD

在上一小節(jié)介紹的非正交網(wǎng)格FDTD方法中，由逆變分量計(jì)算協(xié)變分量時(shí)，不同方向分量間的投影采用了等權(quán)投影作為近似.這在網(wǎng)格的劃分存在梯度時(shí)，可能帶來一定的誤差.為此，本文進(jìn)行以下改進(jìn)，將上述等權(quán)投影改為相應(yīng)位置實(shí)際的不等權(quán)投影，即通過每一點(diǎn)場(chǎng)的分量乘以相應(yīng)位置上的度量張量來計(jì)算相應(yīng)的投影，從而減小在帶地形模型中由于網(wǎng)格劃分存在梯度所帶來的誤差.這樣，式（11）、（12）可修正如下：

2.4 初始條件與邊界條件

在本文數(shù)值模擬中，采用磁偶極子源作為激勵(lì)源，便于與相應(yīng)的理論解進(jìn)行對(duì)比檢驗(yàn).為避免磁偶極子源的奇異性，數(shù)值模擬中常使用某一時(shí)刻的初始場(chǎng)來代替源的影響.模擬脈沖信號(hào)下的電磁場(chǎng)，在柱坐標(biāo)下，磁矩為M的磁偶極子產(chǎn)生的電場(chǎng)的階躍響應(yīng)為［22］

其中，

對(duì)式（17）求t的一階數(shù)值導(dǎo)數(shù)即可以得到電場(chǎng)的初始響應(yīng)，再求旋度即可得到磁場(chǎng)的初始響應(yīng).

在上邊界，由于地形起伏的影響，采用向上延拓邊界條件有一定的困難.在本研究中直接將空氣層作為模型的一部分，用電導(dǎo)率為10－6S／m的高阻層來近似.數(shù)值模擬表明這種近似并不會(huì)對(duì)結(jié)果造成明顯的影響.由于空氣層的引入會(huì)導(dǎo)致空氣層中的波動(dòng)傳播很遠(yuǎn)，因此，對(duì)于外邊界，我們采用卷積完全匹配層（CPML）作為吸收邊界.

3 算法驗(yàn)證

為了驗(yàn)證算法的有效性，我們?cè)诰鶆虬肟臻g模型（圖2，紅實(shí)線為地下均勻介質(zhì)與空氣的界面）的地下介質(zhì)內(nèi)部引入虛擬的高斯地形界面（圖2中黑實(shí)線所示），在圖2中品紅實(shí)線以下的均勻半空間采納非正交網(wǎng)格.整個(gè)計(jì)算區(qū)域（包括品紅線以上網(wǎng)格與底部網(wǎng)格）均采用非正交網(wǎng)格算法進(jìn)行計(jì)算.非正交網(wǎng)格算法在正交網(wǎng)格中自動(dòng)退化為常規(guī)的FDTD算法.在網(wǎng)格的底部，網(wǎng)格的劃分是由非正交網(wǎng)格逐漸向最底層吸收邊界的正交網(wǎng)格過渡得到.

考慮到實(shí)際測(cè)量時(shí)是記錄z方向的電動(dòng)勢(shì)（即磁場(chǎng)Hz分量的響應(yīng)），本文只針對(duì)Hz分量的結(jié)果進(jìn)行討論.圖3給出了上述模型在虛擬地形界面處不同時(shí)刻的理論解（藍(lán)實(shí)線）以及李墩柱的非正交FDTD算法［14］（紅實(shí)線）和本文提出的改進(jìn)的非正交FDTD算法（淡藍(lán)實(shí)線）的數(shù)值模擬結(jié)果.從圖3可以看出，無論是李墩柱算法還是本文的改進(jìn)算法的數(shù)值模擬結(jié)果都與理論解非常接近，兩種數(shù)值算法的結(jié)果也很一致，表明非正交FDTD算法可以有效地模擬帶地形的瞬變電磁場(chǎng).

為了進(jìn)一步考察本文提出的改進(jìn)算法的效果，我們對(duì)圖3的黑色矩形區(qū)域進(jìn)行放大，如圖4a所示.結(jié)果顯示，與李墩柱算法的結(jié)果（紅實(shí)線）相比，本文改進(jìn)算法的結(jié)果（淡藍(lán)實(shí)線）更接近于理論解（藍(lán)實(shí)線），與理論解的誤差分析結(jié)果（圖4（b—d））也表明改進(jìn)算法比李墩柱算法在精度方面有一定的提高.

4 數(shù)值計(jì)算結(jié)果與討論

為簡(jiǎn)單起見，本文統(tǒng)一用高斯型的地形來近似各種地形模型，定義垂直方向的z坐標(biāo)軸正方向向下，0點(diǎn)在地表，地形上任意點(diǎn)（x，y）的z坐標(biāo)為：

圖2 帶虛擬地形的均勻半空間模型的網(wǎng)格劃分剖面圖Fig.2 Gridding profile of a half－space model with a virtual topographical interface

圖3 帶虛擬的200m山谷高斯地形界面的理論解以及非正交網(wǎng)格FDTD算法的數(shù)值模擬結(jié)果的對(duì)比Fig.3 Comparison of the analytical solution and the numerical results of the non－orthogonal FDTD method of a half－space model with a virtual topographical interface of a 200mdepth Gaussian－type subsurface

圖4 （a）圖3中矩形區(qū)域的放大圖；（b）高斯界面處1ms兩種非正交網(wǎng)格算法與理論解誤差對(duì)比；（c）高斯界面處2ms兩種非正交網(wǎng)格算法與理論解誤差對(duì)比；（d）高斯界面處3ms兩種非正交網(wǎng)格與理論解誤差對(duì)比Fig.4 （a）The enlarged figure of the black rectangular region in Fig.3；Error comparison between Li′s FDTD method and modified FDTD method on the Gaussian－type subsurface（b）1ms，（c）2ms，（d）3ms

其中，點(diǎn)（x1，y1）為地形在地表的中心位置，h為地形深度（見圖5a，即地形在點(diǎn)（x1，y1）的z值，正值表示山谷，負(fù)值表示山峰），寬度特征參數(shù)l1（見圖5a）和l2（見圖5b）分別為描述沿水平方向x和y的地形寬度特征的參數(shù).

4.1 地形深度的影響

為了考察地形深度變化對(duì)瞬變電磁場(chǎng)的影響，我們令地形寬度特征參數(shù)l1和l2均為400m，圖5a給出了深度分別為h1＝100m，h2＝200m，h3＝300m和h4＝400m的山谷模型的剖面示意圖，源位于點(diǎn)（0，0，0），地形模型在地表的中心點(diǎn)相對(duì)源的距離｜x1｜為1600m，且在x軸負(fù)方向，地下介質(zhì)電導(dǎo)率為10－2S／m，空氣中電導(dǎo)率近似取為10－6S／m.

圖5 （a）在半空間地表帶高斯型山谷地形模型的剖面圖；（b）半空間地表高斯型山谷地形俯視圖Fig.5 （a）Vertical section of a half－space model with a Gaussian－type topography on the ground surface；（b）Cross－section of a Gaussian－type topographical model

圖6 基于圖5a中的山谷模型計(jì)算得到的地形深度h對(duì)瞬變電磁場(chǎng)的歸一化影響隨時(shí)間的演化圖（a）h＝100m；（b）h＝200m；（c）h＝300m；（d）h＝400m.Fig.6 Evolution of normalized influence due to the varied height，hof a topographical model

由于帶地形情形下的場(chǎng)可以看作半空間的場(chǎng)與地形影響的疊加，為了定量描述地形對(duì)瞬變電磁場(chǎng)的影響，我們將帶地形模型在地表的數(shù)值模擬結(jié)果減去半空間模型的相應(yīng)結(jié)果，并計(jì)算其與半空間模型背景場(chǎng)的比值.作為一個(gè)例子，我們針對(duì)圖5a所示的山谷模型進(jìn)行了數(shù)值模擬，圖6給出了地形深度對(duì)瞬變電磁場(chǎng)的歸一化影響隨時(shí)間演化的結(jié)果.圖中地形的寬度由品紅線標(biāo)出，五角星指示源的位置.

圖6表明，隨著地形深度增加，地形對(duì)場(chǎng)的影響的強(qiáng)度和范圍也逐漸變大，但影響的范圍集中在地形附近.圖中折線的產(chǎn)生是由于在“煙圈”附近的值很小，所以場(chǎng)的變化與背景場(chǎng)的比值很大，等值線變化比較密集.因此，折線的軌跡可以表示“煙圈”傳播的軌跡.

為了進(jìn)一步定量考察地形對(duì)瞬變電磁場(chǎng)的整體影響，我們基于圖5a所示帶地形模型計(jì)算得到的磁場(chǎng)Hz分量，并選擇地形所在的x軸負(fù)半軸作為考察對(duì)象，按照以下公式計(jì)算地形整體影響相對(duì)于半空間模型理論解的歸一化相對(duì)誤差rn＝其中，Hzm和Hzt分別代表帶地形模型數(shù)值解和半空間模型理論解在x軸負(fù)半軸的值.圖7給出了計(jì)算得到的地形深度的整體歸一化影響隨時(shí)間的變化曲線.

圖7 不同時(shí)刻的地形深度整體的歸一化影響Fig.7 Temporal variation of the total topographical effect normalized by a half－space model while the height of topography changes

沿x軸負(fù)半軸各點(diǎn)所受地形影響隨時(shí)間演化結(jié)果（圖6）表明，在早期，地形對(duì)垂直磁場(chǎng)的影響主要集中在地形附近.圖7對(duì)地形整體效應(yīng)的計(jì)算結(jié)果則顯示，地形對(duì)場(chǎng)整體的影響在早期相對(duì)較小，這反映了場(chǎng)的主要能量此時(shí)還沒有傳過來；隨著“煙圈”到達(dá)地形處，地形對(duì)場(chǎng)產(chǎn)生的整體影響相對(duì)較大，但其隨時(shí)間增加的幅度逐漸趨于平緩，在一定程度上反映了地形的影響也隨著場(chǎng)的衰減而減弱.圖8給出了地形最低點(diǎn)處場(chǎng)值隨時(shí)間的變化曲線.可以看出隨著地形深度的增加，該點(diǎn)處場(chǎng)偏離半空間結(jié)果的幅度也相應(yīng)地增加，這與地形整體的歸一化影響相一致.

4.2 地形寬度的影響

由于本研究采用的是高斯型地形模型，其寬度的改變可以通過寬度特征參數(shù)l1和l2來控制.為方便起見，我們定義高斯型地形的寬度d為地形的兩側(cè)的深度減小為最大深度1／e處的距離，即x和y軸方向的寬度分別滿足dx＝2l1和dy＝2l2（參見圖5b的模型俯視圖）.為了考察地形寬度對(duì)瞬變電磁場(chǎng)的影響，作為具體的算例，保持源的位置、地形中心位置以及x軸方向的寬度等參數(shù)與圖5a中400m山谷模型一致（亦即l1不變），令l2的取值從400m變化到700m（步長(zhǎng)為100m），即相應(yīng)y軸方向的寬度由800m變化到1400m（步長(zhǎng)為200m）.

類似圖6的做法，我們將帶地形模型在地表的數(shù)值模擬結(jié)果減去半空間模型的理論解，并相對(duì)半空間的理論解進(jìn)行歸一化.圖9給出了y軸方向地形寬度變化模型對(duì)瞬變電磁場(chǎng)在x軸負(fù)半軸（x軸方向地形寬度由品紅線標(biāo)出）的歸一化影響隨時(shí)間的演化結(jié)果.

圖9表明，盡管地形的y軸方向?qū)挾扔?00m變化到1400m，地形對(duì)x軸負(fù)半軸的磁場(chǎng)垂直分量的影響范圍和幅度并沒有顯著的變化.類似圖7的做法，圖10給出了地形y軸方向?qū)挾茸兓瘜?duì)x軸負(fù)半軸磁場(chǎng)垂直分量的整體影響情況，結(jié)果顯示，與地形深度相比，上述影響依舊不太顯著.

4.3 地形與源距離變化的影響

考慮到源與地形中心的距離可能也是影響瞬變電磁場(chǎng)的相關(guān)因素之一，我們比較了相同地形（深度、寬度）條件下，源與地形中心距離改變時(shí)，地形模型對(duì)瞬變電磁場(chǎng)的影響.作為具體的算例，我們選擇深度為400m的山谷模型，水平方向的兩個(gè)寬度特征參數(shù)均為400m，令地形中心與源的距離分別為1600m，1720m，1840m和1960m.

從圖11可以看出，“煙圈”抵達(dá)地形處的時(shí)間隨著源與地形相對(duì)位置的增加而增加，地形的影響范圍在“煙圈”到達(dá)地形處以前，主要集中在地形附近，幾乎不隨源與地形中心相對(duì)距離的變化而明顯變化，但受影響的強(qiáng)度隨之有所改變.隨著“煙圈”抵達(dá)地形處，地形對(duì)瞬變電磁場(chǎng)的影響范圍和強(qiáng)度都發(fā)生了變化.

圖8 不同地形深度的地形最低點(diǎn)處Hz分量隨時(shí)間變化曲線Fig.8 Temporal variation of component Hzat the lowest point of topography while the height of topography changes

圖9 y軸方向地形寬度dy對(duì)瞬變電磁場(chǎng)在x軸負(fù)半軸的歸一化影響隨時(shí)間的演化（a）dy＝800m；（b）dy＝1000m；（c）dy＝1200m；（d）dy＝1400m.Fig.9 Evolution of normalized influence due to a topographical model with a varied width（dy）along the y－axis

圖10 考察地形寬度影響，不同時(shí)刻的地形整體的歸一化影響Fig.10 Temporal variation of the total topographical effect normalized by a half－space model while the width of topography changes

圖12給出了地形對(duì)瞬變電磁場(chǎng)的整體影響，結(jié)果反映源與地形相對(duì)位置的變化導(dǎo)致地形對(duì)瞬變電磁場(chǎng)的影響隨時(shí)間而變化，不再是單調(diào)的影響.在早期，地形對(duì)場(chǎng)整體的影響相對(duì)有限，但隨著場(chǎng)逐漸到達(dá)地形處，地形對(duì)場(chǎng)的整體影響逐漸增大，且隨地形與源相對(duì)距離的增加而減小，這在一定程度上反映了地形對(duì)場(chǎng)整體的影響隨著距離的增加而推遲.

圖12 考慮地形中心與源距離變化，不同時(shí)刻的地形整體的歸一化影響Fig.12 Temporal variation of the total topographical effect of a topographical model with a varied distance from the source

5 結(jié) 論

本文的數(shù)值模擬研究表明，非正交網(wǎng)格可以有效地模擬帶地形的瞬變電磁場(chǎng)，與以往帶地形瞬變電磁正演中采納的梯度網(wǎng)格和垂向變換網(wǎng)格相比，非正交網(wǎng)格更利于刻畫精細(xì)的地形，從而提高算法精度.本文提出的采用非等權(quán)投影的改進(jìn)非正交網(wǎng)格FDTD算法，進(jìn)一步提高了帶地形瞬變電磁場(chǎng)三維正演算法的精度.

利用我們提出的改進(jìn)非正交網(wǎng)格FDTD算法，對(duì)不同地形模型進(jìn)行了數(shù)值模擬，結(jié)果表明，地形的深度、寬度以及地形與源的相對(duì)位置等都可能影響瞬變電磁場(chǎng)的測(cè)量結(jié)果，但影響程度不一.瞬變電磁場(chǎng)的地形響應(yīng)與地形的深度存在正相關(guān)，而寬度的影響可能與測(cè)線方向和寬度變化方向有關(guān)，地形與源的相對(duì)位置變化對(duì)瞬變電磁場(chǎng)的影響隨時(shí)間并非單調(diào)變化，但隨著“煙圈”的抵達(dá)時(shí)間推后，地形對(duì)瞬變電磁場(chǎng)歸一化影響效果也相應(yīng)地延遲.

致 謝 感謝兩位審稿專家提出建設(shè)性修改意見.

（

）

［1］ Mitsuhata Y.2－D electromagnetic modeling by finite－element method with a dipole source and topography.Geophysics，2000，65（2）：465－475.

［2］ 肖明順.帶地形的瞬變電磁2.5維有限元數(shù)值模擬研究［碩士論文］.北京：中國地質(zhì)大學(xué)，2008.Xiao M S.Study on 2.5Dnumerical modeling for transient electromagnetic method by finite element method with topography［Master′s thesis］（in Chinese）.Beijing：China University of Geosciences，2008.

［3］ Sugeng F.Modeling the 3DTDEM response using the 3D full－domain finite－element method based on the hexahedral edge－element technique.ExplorationGeophysics，1998，29（4）：615－619.

［4］ Li J H，Zhu Z Q，Liu S C，et al.3Dnumerical simulation for the transient electromagnetic field excited by the central loop based on the vector finite－element method.Journalof GeophysicsandEngineering，2011，8（4）：560－567.

［5］ Xiong B.2.5Dforward for the transient electromagnetic response of a block linear resistivity distribution.Journalof GeophysicsandEngineering，2011，8（1）：115－121.

［6］ 唐新功，胡文寶，嚴(yán)良俊.地塹地形對(duì)長(zhǎng)偏移距瞬變電磁測(cè)深的影響研究.工程地球物理學(xué)報(bào)，2004，1（4）：313－317.Tang X G，Hu W B，Yan L J.Graben topographyic effects to the long offset transient electromagnetic responses.ChineseJournalofEngineeringGeophysics（in Chinese），2004，1（4）：313－317.

［7］ 唐新功，胡文寶，嚴(yán)良俊.瞬變電磁法對(duì)存在山谷地形時(shí)的多個(gè)異常體的探測(cè)能力研究.地震地質(zhì)，2005，27（2）：316－323.Tang X G，Hu W B，Yan L J.The detectability of transient electromagnetic method to multiple 3Dbodies with valley topography.SeismologyandGeology（in Chinese），2005，27（2）：316－323.

［8］ Tang X G，Hu W B，Yan L J.Topographic effects on long offset transient electromagnetic response.AppliedGeophysics，2011，8（4）：277－284.

［9］ 唐新功，胡文寶，嚴(yán)良俊.斜坡和垂直斷層地形的瞬變電磁響應(yīng).石油天然氣學(xué)報(bào)，2011，33（10）：57－60.Tang X G，Hu W B，Yan L J.The transient electromagnetic responses of slope and vertical fault topography.Journalof OilandGasTechnology（in Chinese），2011，33（10）：57－60.

［10］ Zhdanov M S，Lee S K，Yoshioka K.Integral equation method for 3Dmodeling of electromagnetic fields in complex structures with inhomogeneous background conductivity.Geophysics，71（6）：G333－G345.

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

［12］ Endo M，Noguchi K.Chapter 6three－dimensional modeling considering the topography for the case of time－domain electromagnetic method.MethodsinGeochemistryand Geophysics，2002，35：85－107.

［13］ 肖懷宇.帶地形的瞬變電磁法三維數(shù)值模擬研究［碩士論文］.北京：中國地質(zhì)大學(xué)，2006.Xiao H Y.Three－dimensional numerical modeling considering the topography of TEM［Master′s thesis］（in Chinese）.Beijing：China University of Geosciences，2006.

［14］ 李墩柱.用正交網(wǎng)格模擬帶地形的瞬變電磁場(chǎng)［碩士論文］.北京：北京大學(xué)，2010.Li D Z.TEM simulation with topography using nonorthogonal grid FDTD［Master′s thesis］（in Chinese）.Beijing：Peking University，2010.

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

［16］ 許洋鋮，林君，李肅義等.全波形時(shí)間域航空電磁響應(yīng)三維有限差分?jǐn)?shù)值計(jì)算.地球物理學(xué)報(bào)，2012，55（6）：2105－2114.Xu Y C，Lin J，Li S Y，et al.Calculation of full－waveform airborne electromagnetic response with three－dimension finitedifference solution in time－domain.ChineseJ.Geophys.（in Chinese），2012，55（6）：2105－2114.

［17］ 孫懷鳳，李貅，李術(shù)才等.考慮關(guān)斷時(shí)間的回線源激發(fā)TEM三維時(shí)域有限差分正演.地球物理學(xué)報(bào)，2013，56（3）：1049－1064.Sun H F，Li X，Li S C，et al.Three－dimensional FDTD modeling of TEM excited by a loop source considering ramp time.ChineseJ.Geophys.（in Chinese），2013，56（3）：1049－1064.

［18］ Du Fort E C，F(xiàn)rankel S P.Stability conditions in the numerical treatment of parabolic differential equations.MathematicalTablesandOtherAidstoComputation，1953，7（43）：135－152.

［19］ Holland R.Finite－difference solution of Maxwell′s equations in generalized nonorthogonal Coordinates.IEEETransactionson NuclearScience，1983，30（6）：4589－4591.

［20］ Lee J F，Palandech R，Mittra R.Modeling three－dimensional discontinuities in waveguides using nonorthogonal FDTD algorithm.IEEETransactionsonMicrowaveTheoryand Techniques，1992，40（2）：346－352.

［21］ Yee K S.Numerical solution of initial boundary value problems involving Maxwell′s equations in isotropic media.IEEETransactionsonAntennasandPropagation，1966，14（3）：302－307.

［22］ 陳丹丹.瞬變電磁法三維正演研究［碩士論文］.北京：中國地質(zhì)大學(xué)，2008.Chen D D.Study of three－dimensional forward of TEM［Master′s thesis］（in Chinese）.Beijing：China University of Geosciences，2008.