基于共形網格技術的ATEM三維數值模擬

嵇艷鞠　任桂瑩　關珊珊　黎東升

摘? ?要：共形網格技術可減小傳統有限差分方法帶來的階梯狀網格近似誤差. 將共形網格技術運用到無源麥克斯韋方程組的推導中，采用線性加權平均方法，求得棱邊電場的等效電導率，更新了傳統時域有限差分方法的電場迭代方程;研究了具有解析式目標體的共形參數確定方法，利用投影思想得到共形網格位置坐標，進而求得共形參數. 通過與有限元方法計算得出的電磁響應進行比較，驗證了本文方法的正確性;共形網格技術在保證計算效率的條件下，與階梯狀網格近似方法相比提高了曲面三維異常體的計算精度;需要共形處理的電場越多，共形網格技術的效果越明顯，精度提高越多. 該方法為開展地下復雜異常體的航空電磁三維正反演奠定了基礎.

關鍵詞：航空電磁法;共形網格技術;三維曲面異常體;時域有限差分方法;線性加權平均

中圖分類號：P631? ? ? ? ? ? ? ? ? ? ? ? ? ? 文獻標志碼：A

Abstract：Conform mesh technique is usually adopted to eliminate the step approximation error caused by the traditional Finite Difference Time Domain（FDTD） method. In this paper，the conformal mesh technique was introduced into the derivation of passive Maxwell's equations. To update the electric field iterative equation of the traditional Finite Difference Time Domain method，a linear weighted average method was used to obtain the equivalent conductivity of the edge electric field. Using the projection theory，the position coordinates of the conformal mesh were calculated，and then the conformal parameters of the target body were obtained. The electromagnetic responses of this method are consistent with that of the finite element method. Compared with the step approximation method，the conformal mesh technique shows a higher computational accuracy when calculating a three-dimensional anomaly with a curved surface. The more the electric field is required to be processed conformally，the more obvious the effect of the conformal mesh technique is，the accuracy is also improved greatly as a result. This method sets the fundamental to calculate complex underground anomalies of airborne electromagnetic.

Key words：airborne electromagnetic surveys;conformal mesh technique;three-dimensional curved anomaly surface;Finite Difference Time Domain（FDTD） method;linear weighted average

航空瞬變電磁方法（Airborne Transient Electromagnetic Method，ATEM）采用飛行器搭載電磁系統對大地的二次場進行快速測量，通過分析二次場獲得地下電阻率分布，從而了解地質體的構造. 相比于地面瞬變電磁法，航空瞬變電磁具有探測速度快，范圍廣，可在短時間內獲得較多信息等特點，更適用于森林、沼澤、割裂地形、山區等地形復雜區域的勘探. 航空瞬變電磁法在國內外發展較為迅速，在淺覆蓋區尋找金屬礦、隱蔽洞體等典型目標體探測方面起到了重要的作用[1-4]. 常用的瞬變電磁三維正演模擬方法包括有限差分法、積分方程法、有限元方法等，其中有限差分法具有數學表達簡單、直觀高效等特點，因此發展得較為成熟. Wang等人[5]采用有限差分法進行瞬變電磁三維數值模擬. Commer等[6]2004年利用并行有限差分方法實現了電性源瞬變電磁三維數值模擬，2015年對時間步長進行了改進以減少計算時間[7]. Li等人[8]基于有限差分法研究了含水結構的瞬變電磁三維數值模擬. 隨著瞬變電磁探測技術的發展，地下復雜三維異常體結構的精細化勘探得到關注，采用傳統的時域有限差分（Finite Difference Time Domain，FDTD）方法進行地下不規則異常體的數值模擬時，常常會由于邊界處網格的階梯狀網格近似導致較大誤差，故如何處理復雜異常體的邊界問題以及提高三維不規則體的計算精度，成為地球物理勘探中的熱點問題. Cao等人[9]提出了自適應有限差分方法，采用局部細化的方式提高了三維不規則異常體的模擬精度，但增加了計算的復雜性. Chang等人[10]采用非均勻網格方式對煤礦充水區域的全波形瞬變電磁響應進行了數值模擬，但網格剖分較難. 徐巖等人[11]在瞬變電磁三維正演模擬時，將起伏地形或不規則體表面網格進行了精細剖分，雖然提高了計算精度，但會增加計算量.

共形網格技術是近年來新興的一種數值模擬技術，它基于有限差分方法，對目標體邊界進行特殊處理，以減小傳統FDTD方法所產生的階梯狀網格近似誤差，提高了有限差分方法的計算精度，在不增加網格數量的基礎上，提高了計算效率. Ilinca等人[12]將共形網格技術應用于材料物理學中，解決了材料復雜界面的熱傳遞問題. Wang等人[13]提出了一種改進的時域有限差分共形方案，可快速預測雷達截面以及復雜三維結構的感應電流分布. Guo等人[14]利用改進的共形FDTD網格劃分算法模擬復雜三維模型. 武超[15]將共形網格技術應用于無線電物理方面，研究了復雜目標的電磁特性計算. 范宜仁等人[16]將共形網格技術應用于隨鉆電磁波測井響應數值模擬中，為隨鉆快速反演提供支持. Nicolas等人[17]提出了一種基于六面體網格的共形方案，并在材料邊緣處進行特定分裂，實現了細化和共形. Wei等人[18]擴展了共形時域有限差分方法，將該方法應用于地震中的彈性波計算，并取得了很好的結果. Cabello等人[19]將共形網格技術應用于雷達散射截面的計算中，取得了較好的效果. 但在三維地質體探測中，如何采用共形網格技術更好地擬合曲面異常體的邊界問題還有待研究.

本文在傳統FDTD方法的基礎上，結合共形網格技術開展了航空瞬變電磁的數值模擬研究. 基于介質參數等效的思想，更新了傳統FDTD方法的電場表達式;對于具有解析式的典型目標體，討論了基于解析式建模的FDTD共形網格生成過程，采用共形網格技術實現了地下曲面異常體的三維數值模擬，提高了傳統FDTD方法對曲面異常體進行數值模擬時的精度.

1? ?共形網格技術

利用常規FDTD方法進行數值模擬時，Yee元胞僅可以賦一個電導率值，在模擬曲面異常體邊界時勢必會帶來由于階梯狀網格近似而產生的誤差. 圖1為球形異常體二維截面階梯狀網格近似圖，介質1表示目標體內部，介質2表示目標體外部，其中陰影區域表示位于目標體表面附近的網格，為了減小誤差，在附近網格處，采用介質參數等效技術，更新傳統FDTD迭代方程，這種技術叫做共形網格技術，而目標體附近的網格被稱作共形網格.

在計算時，從曲面異常體的邊界出發，確定共形網格，基于線性加權平均思想，重新計算邊界共形網格處電場的等效電導率值，將等效電導率代入相應位置的電場表達式，進行更精確的數值模擬[20]. 其中，對于具有解析式的異常體模型，可直接根據其解析式以及整個計算區域網格的坐標，得到具體共形網格的坐標位置以及需要共形處理的電場所在棱邊位置，為異常體表面的共形計算提供參數.

2? ?基于共形網格技術的電場迭代方程推導

圖2為航空瞬變電磁法的工作原理示意圖. 當發射線圈中通入發射電流時，在整個空間產生一次場，發射電流的瞬間變化在大地中形成渦流并產生二次場，接收線圈接收到二次場后以感應電壓的形式記錄下來，文中的數值模擬采用發射線圈與接收線圈共中心的方式.

圖8（a）為階梯狀網格近似和共形網格技術處理后的電磁響應與有限元方法的電磁響應對比圖. 由于有限元方法采用四面體模擬地下異常體，故該方法可精確剖分三棱柱模型，因此將有限元方法對三棱柱模型的數值模擬作為基準解. 由圖8（a）可知，3種方法得出的曲線一致性較好. 從圖8（b）中可以看出，在1 ms之前，兩種處理方法的相對誤差相差不大，而在1 ms后，可看出利用共形網格技術處理后的相對誤差小于階梯狀網格近似方法處理的相對誤差，減小了傳統FDTD方法在目標剖分過程中引入的階梯狀網格近似誤差，驗證了共形網格技術的有效性. 由于在對三棱柱模型進行階梯狀網格近似時，采用的最小網格為10 m，因此兩者之間的誤差相差不大，當增大最小剖分網格尺寸后，經共形網格技術處理后的數值模擬優勢會更加明顯. 圖8（c）為空中垂直磁場的等值線分布.

4.2? ?球形異常體模型的數值模擬

采用與4.1節中同樣的系統模型參數，建立

4 060 m × 4 060 m × 2 950 m的包含吸收邊界層的計算區域，大地被剖分為111×111×77個網格. 初始場及時間步長的選取方式與4.1節相同. 在發射源的正下方設置球形異常體模型，半徑為100 m，埋深為100 m，剖分異常體的網格步長為10 m，計算區域示意圖如圖9（a）所示. 以球形異常體的最大橫截面為例，作如圖9（b）所示的球形異常體截面階梯狀網格近似示意圖，并將與球形異常體邊界有交點的棱邊處的電場進行共形處理. 其中，異常體電導率為1 S/m，圍巖電導率為0.01 S/m. 球形異常體模型基準解采用5 m×5 m×5 m的小網格進行剖分. 計算區域不變，將球形異常體附近網格剖分成5 m×5 m×5 m的小網格，此時大地被剖分為163×163×160個網格，將此時的電磁響應作為基準解進行分析. 圖10（a）為階梯狀網格近似和共形網格技術處理后的電磁響應與基準解的電磁響應對比圖. 圖10（b）為兩種處理方法的相對誤差曲線，從圖中可以看出在0.4 ms之前，兩種處理方法的誤差相差不大，而在0.4 ms后，共形處理的相對誤差明顯小于階梯狀網格近似處理的相對誤差，共形網格技術在處理復雜異常體邊界問題時具有明顯優勢，相對誤差減小了1.65%. 圖10（c）為空中垂直磁場的等值線分布圖.

圖11為球形異常體電磁響應切片圖，異常體由于埋深較淺，對源的影響在早期較為明顯，隨著時間的推移，電磁響應逐漸向下、向外擴散，且響應逐漸減小.

4.3? ?橢球形異常體模型的數值模擬

系統模型參數的選取方式均與4.1節相同. 計算區域及網格剖分數與4.2節相同. 在發射源的正下方放置一個半長軸a=150 m，半短軸b=50 m，埋深為100 m的橢球形異常體，采用基于解析式建模的共形網格技術求解橢球形異常體的時間域電磁響應. 圖12（a）為計算區域示意圖，圖12（b）為橢球截面階梯狀網格近似示意圖. 其中，異常體電導率為1 S/m，圍巖電導率為0.01 S/m. 將橢球形異常體附近網格剖分成5 m×5 m×5 m的小網格，網格剖分數同4.2節中的描述，以此時的電磁響應作為基準解進行分析. 圖13（a）為階梯狀網格近似和共形網格技術處理后的電磁響應與基準解的電磁響應對比圖. 圖13（b）為兩種處理方法的相對誤差曲線圖，從圖中可以看出在0.5 ms之前，兩種處理方法的誤差相差不大，而在0.5 ms后共形處理的相對誤差明顯小于階梯狀網格近似處理的相對誤差，其中相對誤差減小了4.64%. 與圖10（b）中曲線進行比較，在處理橢球形異常體時，共形網格處理的相對誤差相對于階梯狀網格近似處理的相對誤差減小較為明顯. 通過分析圖9（b）、圖12（b）可知，需要共形處理的棱邊電場越多，共形處理的效果越明顯，精度就越高.

圖14為橢球形異常體電磁響應切片圖. 通過比較球形與橢球形異常體在t = 0.5 ms，y = 1 980 m的電磁響應切片圖，在電磁場擴散的過程中，能夠反映出三維異常體的近似形狀，并且橢球形異常體的擴散速度要高于球形異常體的電磁場擴散速度. 圖15為t = 7 ms時，z分別為150、200、250 m時的橢球形異常體感應電動勢等值線分布圖.

本文分別計算了三棱柱、球形異常體、橢球形異常體作為三維異常體時的時間域電磁響應，在相同運行環境下，計算效率如表1所示. 其中三棱柱模型利用共形網格技術以及階梯狀網格近似方法的運行時間分別為2 460 s和2 214 s;球形異常體模型利用共形網格技術以及階梯狀近似方法的運行時間分別為2 743 s和1 447 s;橢球形異常體模型利用共形網格技術以及階梯狀網格近似方法的運行時間分別為2 909 s和1 564 s. 由于三棱柱模型的共形參數確定只涉及到二維平面，故相對容易獲得，因此與階梯狀網格近似方法的運行時間相比較相差不多，耗時為階梯狀網格近似處理的1.11倍. 球形異常體模型和橢球形異常體模型為具有解析式的目標體，需根據解析式確定共形參數，因此在電場和磁場的每一次迭代計算中都需要判斷x、y、z三個方向共形網格的位置，即異常體表面與網格相交位置，以及交點所分相應棱邊的內外邊長. 并將共形參數代入迭代式中進行計算，因此與三棱柱模型相比要更耗時，與階梯狀網格近似相比，其時間比值分別為1.90倍和1.86倍. 由此可知，共形網格技術可在保證計算效率的同時提高計算精度.

5? ?結? ?論

本文基于傳統FDTD方法，結合了共形網格技術，采用介質參數等效思想，求得了線性加權平均電導率，更新了電場迭代方程;研究了有解析式的三維目標體共形網格技術;采用三棱柱模型驗證了共形網格技術的有效性，提高了地下曲面三維異常體的數值模擬精度. 主要結論如下：

1）在不改變網格大小及數量的基礎上，共形網格技術僅對電導率進行線性加權平均處理，并將計算求得的等效電導率代入傳統FDTD電場表達式，更新了電場迭代方程. 共形網格技術在保證計算效率的條件下，與階梯狀網格近似方法相比提高了曲面三維異常體的計算精度.

2）可通過三維坐標投影法確定出共形網格位

置，進而判斷出所截網格棱邊位置，以確定出曲面共形參數.

3）在共形網格處理中，網格棱邊與異常體邊界的交點越多，即需要共形處理的電場越多，共形網格技術的效果越明顯，精度提高越多.

本文目前開展的工作主要是針對有解析式的三維目標模型進行共形網格技術的研究，而對于無法寫出解析式的任意復雜模型，則應分以下幾點進行共形網格技術處理：1）通過3ds Max等建模軟件建立三維目標模型. 2）根據目標模型，生成.stl文件，導出三角面元數據. 3）由三角面元數據，通過編程生成Yee網格坐標以及共形參數，并導入3ds Max中進行可視化. 4）結合有限差分方法與共形網格技術進行計算，最終實現任意三維復雜模型的高精度數值模擬.

參考文獻

[1]? ? 許洋鋮，林君，李肅義，等. 全波形時間域航空電磁響應三維有限差分數值計算[J]. 地球物理學報，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 finite-difference solution in time-domain[J]. Chinese Journal of Geophysics，2012，55（6）：2105—2114.（In Chinese）

[2]? ? 駱燕.航空瞬變電磁法在我國淺覆蓋區尋找多金屬礦應用前景[C]//中國地質學會2013年學術年會論文摘要匯編——S02資源與環境地球物理勘查理論與方法技術分會場. 昆明：中國地質學會，2013：239—241.

LUO Y. Application of aviation transient electromagnetic method to search for polymetallic mines in China[C]//China Geological Society 2013 Annual Conference Abstracts Compilation—S02 Resources and Environmental Geophysical Exploration Theory and Methodology Sub-venue. Kunming：China Geological Society，2013：239—241.（In Chinese）

[3]? ? 殷長春，任秀艷，劉云鶴，等. 航空瞬變電磁法對地下典型目標體的探測能力研究[J]. 地球物理學報，2015，58（9）：3370—3379.

YIN C C，REN X Y，LIU Y H，et al. Research on the detection ability of aerial transient electromagnetic method to typical targets in underground[J]. Chinese Journal of Geophysics，2015，58（9）：3370—3379. （In Chinese）

[4]? ? 劉偉. 航空瞬變電磁法探測地下隱蔽洞體的前景淺析[J]. 兵器裝備工程學報，2017，38（6）：169—175.

[16]? 范宜仁，李煒，李虎，等. 基于時域有限差分亞網格與共形網格技術的隨鉆電磁波測井響應數值模擬[J]. 測井技術，2015，39：561—566.

FAN Y R，LI W，LI H，et al. Numerical simulation of electromagnetic wave logging response while drilling based on time domain finite difference subgrid and conformal grid technology [J]. Well Logging Technology，2015，39：561—566. （In Chinese）

[17]? NICOLAS G，FOUQUET T，GENIAUT S，et al. Improved adaptive mesh refinement for conformal hexahedral meshes[J]. Advances in Engineering Software，2016，102：14—28.

[18]? WEI S L，ZHOU J Y，ZHUANG M W，et al. A 3-D enlarged cell technique （ECT） for elastic wave modelling of a curved free surface[J]. Geophysical Journal International，2016，206：1921—1932.

[19]? CABELLO M R，ANGULO L D，ALVAREZ J，et al. A new efficient and stable 3D conformal FDTD[J]. IEEE Microwave and Wireless Components Letters，2016，26：553—555.

[20]? YU W Y，MITTRA R. A conformal finite difference time domain technique for modeling curved dielectric surfaces[J]. IEEE Microwave & Wireless Components Letters，2002，11：25—27.

[21]? 葛德彪，閆玉波. 電磁波時域有限差分方法[M]. 3版. 西安：西安電子科技大學出版社，2011：9—28.

GE D B，YAN Y B. Finite-difference time-domain method for electromagnetic waves [M]. 3rd ed. Xian：Xidian University Press，2011：9—28. （In Chinese）

[22]? 胡曉娟. 復雜目標電磁散射的FDTD及FDFD算法研究[D]. 西安：西安電子科技大學物理與光電工程學院，2007：16—17.

HU X J. Study of FDTD and FDFD algorithms for electromagnetic scattering by complex targets[D]. Xian：School of Physics and Optoelectronic Engineering，Xidian University，2007：16—17.（In Chinese）

[23]? 關珊珊. 基于GPU的三維有限差分直升機瞬變電磁響應并行計算[D]. 長春：吉林大學儀器科學與電氣工程學院，2012：19—27.

GUAN S S. Parallel calculation of 3D FDTD helicopter transient electromagnetic response based on the GPU[D]. Changchun：College of Instrument Science and Electrical Engineering，Jilin University， 2012：19—27.（In Chinese）

[24]? WANG T，HOHMANN G W. A finite-difference，time-domain solution for thtee-dimensional electromagnetic modeling[J]. Geophysics，1993，58：797—809.