基于共形網(wǎng)格技術(shù)的ATEM三維數(shù)值模擬

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

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

關鍵詞：航空電磁法;共形網(wǎng)格技術(shù);三維曲面異常體;時域有限差分方法;線性加權(quán)平均

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

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

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

1? ?共形網(wǎng)格技術(shù)

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

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

2? ?基于共形網(wǎng)格技術(shù)的電場迭代方程推導

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

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

4.2? ?球形異常體模型的數(shù)值模擬

采用與4.1節(jié)中同樣的系統(tǒng)模型參數(shù)，建立

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

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

4.3? ?橢球形異常體模型的數(shù)值模擬

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

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

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

5? ?結(jié)? ?論

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

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

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

置，進而判斷出所截網(wǎng)格棱邊位置，以確定出曲面共形參數(shù).

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

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

參考文獻

[1]? ? 許洋鋮，林君，李肅義，等. 全波形時間域航空電磁響應三維有限差分數(shù)值計算[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]? ? 駱燕.航空瞬變電磁法在我國淺覆蓋區(qū)尋找多金屬礦應用前景[C]//中國地質(zhì)學會2013年學術(shù)年會論文摘要匯編——S02資源與環(huán)境地球物理勘查理論與方法技術(shù)分會場. 昆明：中國地質(zhì)學會，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]? 范宜仁，李煒，李虎，等. 基于時域有限差分亞網(wǎng)格與共形網(wǎng)格技術(shù)的隨鉆電磁波測井響應數(shù)值模擬[J]. 測井技術(shù)，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，F(xiàn)OUQUET 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.