基于非正交網(wǎng)格的帶地形三維瞬變電磁場模擬

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

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

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

1 引 言

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

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

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

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

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

2.1 非正交網(wǎng)格

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

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

其中，

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

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

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

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

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

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

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

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

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

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

同時，因為瞬變電磁場近似為擴散場，為了保證方程中的波動不會影響到擴散項，還需要保證最大時間步長［15］：

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

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

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

2.4 初始條件與邊界條件

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

其中，

對式（17）求t的一階數(shù)值導數(shù)即可以得到電場的初始響應，再求旋度即可得到磁場的初始響應.

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

3 算法驗證

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

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

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

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

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

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

圖3 帶虛擬的200m山谷高斯地形界面的理論解以及非正交網(wǎng)格FDTD算法的數(shù)值模擬結(jié)果的對比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)格算法與理論解誤差對比；（c）高斯界面處2ms兩種非正交網(wǎng)格算法與理論解誤差對比；（d）高斯界面處3ms兩種非正交網(wǎng)格與理論解誤差對比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

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

4.1 地形深度的影響

為了考察地形深度變化對瞬變電磁場的影響，我們令地形寬度特征參數(shù)l1和l2均為400m，圖5a給出了深度分別為h1＝100m，h2＝200m，h3＝300m和h4＝400m的山谷模型的剖面示意圖，源位于點（0，0，0），地形模型在地表的中心點相對源的距離｜x1｜為1600m，且在x軸負方向，地下介質(zhì)電導率為10－2S／m，空氣中電導率近似取為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中的山谷模型計算得到的地形深度h對瞬變電磁場的歸一化影響隨時間的演化圖（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

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

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

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

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

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

4.2 地形寬度的影響

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

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

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

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

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

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

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

圖9 y軸方向地形寬度dy對瞬變電磁場在x軸負半軸的歸一化影響隨時間的演化（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 考察地形寬度影響，不同時刻的地形整體的歸一化影響Fig.10 Temporal variation of the total topographical effect normalized by a half－space model while the width of topography changes

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

圖12 考慮地形中心與源距離變化，不同時刻的地形整體的歸一化影響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)格可以有效地模擬帶地形的瞬變電磁場，與以往帶地形瞬變電磁正演中采納的梯度網(wǎng)格和垂向變換網(wǎng)格相比，非正交網(wǎng)格更利于刻畫精細的地形，從而提高算法精度.本文提出的采用非等權(quán)投影的改進非正交網(wǎng)格FDTD算法，進一步提高了帶地形瞬變電磁場三維正演算法的精度.

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

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

（

）

［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ì)大學，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］ 唐新功，胡文寶，嚴良俊.地塹地形對長偏移距瞬變電磁測深的影響研究.工程地球物理學報，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］ 唐新功，胡文寶，嚴良俊.瞬變電磁法對存在山谷地形時的多個異常體的探測能力研究.地震地質(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］ 唐新功，胡文寶，嚴良俊.斜坡和垂直斷層地形的瞬變電磁響應.石油天然氣學報，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ì)大學，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)格模擬帶地形的瞬變電磁場［碩士論文］.北京：北京大學，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ù)值計算.地球物理學報，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)斷時間的回線源激發(fā)TEM三維時域有限差分正演.地球物理學報，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ì)大學，2008.Chen D D.Study of three－dimensional forward of TEM［Master′s thesis］（in Chinese）.Beijing：China University of Geosciences，2008.