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

嵇艷鞠，任桂瑩，關(guān)珊珊?，黎東升

（1.吉林大學(xué) 地球信息探測儀器教育部重點實驗室，吉林長春130026；2.吉林大學(xué) 儀器科學(xué)與電氣工程學(xué)院，吉林 長春130026；3.吉林大學(xué)國家地球物理探測儀器工程技術(shù)研究中心，吉林長春130026）

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

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

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

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

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

圖1 球形異常體二維截面階梯狀網(wǎng)格近似圖Fig.1 Spherical two-dimensional section step approximation

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

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

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

圖2 航空瞬變電磁工作原理圖Fig.2 Schematic diagram of airborne transient electromagnetic working principle

在用有限差分方法對三維目標(biāo)體進行數(shù)值模擬時，基于Maxwell旋度方程進行差分離散，式（1）和（2）為時域Maxwell旋度方程.常規(guī)FDTD方法在數(shù)值模擬時，需要將整個計算區(qū)域離散成規(guī)則的六面體，即近似為圖3所示的Yee元胞.電場和磁場在時間上相差半個步長交替取樣，Maxwell旋度方程離散后推導(dǎo)得到顯式差分方程.擬合復(fù)雜模型的曲面邊界時，在處理介質(zhì)共形中，只涉及到電導(dǎo)率的變化，因此只需更新Maxwell方程組的電場迭代方程[21].

式中：e（r，t）為電場強度；b（r，t）為磁感應(yīng)強度；h（r，t）為磁場強度；j（r，t）為傳導(dǎo)電流密度；ε（r）為介電常數(shù).

圖3 FDTD中的Yee元胞Fig.3 The Yee cell of FDTD

基于共形網(wǎng)格技術(shù)，將傳統(tǒng)電場迭代方程中的電導(dǎo)率進行等效替代.圖4為共形網(wǎng)格中電導(dǎo)率的等效.圖4（a）為三維FDTD共形網(wǎng)格示意圖，網(wǎng)格被目標(biāo)體分為介質(zhì)1和介質(zhì)2兩部分，由于介質(zhì)不同，電導(dǎo)率存在差異，因此，需求出共形網(wǎng)格中各電場位置的等效電導(dǎo)率.圖4（b）為三維共形網(wǎng)格在xoy面上的投影，介質(zhì)1的電導(dǎo)率為σ1，介質(zhì)2的電導(dǎo)率為σ2，介質(zhì)1與介質(zhì)2的分界面將一條棱邊分為兩部分，因電場分布在棱邊上，故相應(yīng)棱邊上的等效電導(dǎo)率值由分界面所截長度的線性加權(quán)平均得到.式（3）為A邊上x方向電場的等效電導(dǎo)率.

以x方向為例，將式（3）代入傳統(tǒng)電場迭代方程中，得式（4）所示更新后的電場方程.y、z方向同理可得.

圖4 共形網(wǎng)格中電導(dǎo)率的等效Fig.4 Equivalent of conductivity in a conformal grid

3 共形參數(shù)確定方法

通過分析三維共形網(wǎng)格可知，在利用共形網(wǎng)格技術(shù)進行數(shù)值模擬時，其關(guān)鍵是得到需要共形的網(wǎng)格位置以及網(wǎng)格共形參數(shù)lz1和lz2.本文以球形異常體表面的共形網(wǎng)格為例，分析說明具有解析式目標(biāo)體的共形參數(shù)確定方法.

若交點z滿足式（5），則可判斷z所在的網(wǎng)格編號為k：

設(shè)zmin所在的網(wǎng)格編號為kmin，zmax所在的網(wǎng)格編號為kmax，lz1k為z方向編號為k的網(wǎng)格線被目標(biāo)體所截的內(nèi)部長度，lz2k為z方向編號為k的網(wǎng)格線被目標(biāo)體所截的外部長度，由此可得z方向的共形網(wǎng)格坐標(biāo)及共形參數(shù)，如式（6）和式（7）所示.

圖5 球形異常體共形網(wǎng)格參數(shù)示意圖Fig.5 Schematic diagram of spherical conformal mesh parameters

同理，可分別得到x方向、y方向的共形網(wǎng)格坐標(biāo)及共形參數(shù).

4 模型計算及分析

4.1 精度驗證

為了驗證共形網(wǎng)格技術(shù)的有效性，建立5 980 m×5 980 m×3 910 m的包含吸收邊界條件的計算區(qū)域，大地被剖分為127×127×63個非均勻網(wǎng)格.采用發(fā)射與接收共中心的測量裝置，發(fā)射線圈位于整個計算區(qū)域的中心，距地面高度為30 m，發(fā)射波形為負(fù)階躍波，電流大小為300 A，其中三維異常體為三棱柱模型，位于發(fā)射線圈的正下方，截面為等腰三角形，腰長為，底邊長為120 m，三棱柱高度為300 m，三維異常體埋深為250 m，異常體電導(dǎo)率為1 S/m，圍巖電導(dǎo)率為0.01 S/m.圖6為三棱柱模型設(shè)置示意圖.

圖6 三棱柱模型設(shè)置示意圖Fig.6 Schematic diagram of the setting of the triangular prism model

在數(shù)值計算時，采用加載初始場的方式代替實際探測中的發(fā)射源，采用數(shù)字濾波算法解決計算式中含有的一階貝塞爾函數(shù)積分問題，之后利用線性濾波方法進行頻時轉(zhuǎn)換，得到時間域的初始場值[23-24].

為了驗證采用共形網(wǎng)格技術(shù)進行數(shù)值模擬的有效性，將階梯狀網(wǎng)格近似結(jié)果及共形網(wǎng)格處理后結(jié)果與有限元方法進行了對比.三棱柱模型截面電導(dǎo)率階梯狀網(wǎng)格近似賦值示意圖如圖7（a）所示，其中虛線內(nèi)部近似為異常體電導(dǎo)率值.在利用共形網(wǎng)格技術(shù)處理時，需要將共形網(wǎng)格內(nèi)部電場處的電導(dǎo)率作等效，截面上需共形處理的電場位置示意圖如圖7（b）所示，以三棱柱截面左側(cè)為例，箭頭所在位置表示x、y方向需共形的電場位置坐標(biāo).

圖7 三棱柱模型截面共形處理示意圖Fig.7 Schematic diagram of a conformal mesh of a triangular prism model cross section

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

圖8 三棱柱模型計算結(jié)果圖Fig.8 Calculation results of the triangular prism model

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

圖9 球形異常體模型設(shè)置示意圖Fig.9 Schematic diagram of the sphere model setting

圖10 球形異常體模型計算結(jié)果圖Fig.10 Calculation results of sphere model

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

圖11 球形異常體電磁響應(yīng)切片圖Fig.11 Electromagnetic response slice diagram of spherical anomalous body

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

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

圖12 橢球形異常體模型設(shè)置示意圖Fig.12 Schematic diagram of the ellipsoid model setting

圖13 橢球形異常體模型計算結(jié)果圖Fig.13 Calculation results of ellipsoid model

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

本文分別計算了三棱柱、球形異常體、橢球形異常體作為三維異常體時的時間域電磁響應(yīng)，在相同運行環(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倍.球形異常體模型和橢球形異常體模型為具有解析式的目標(biāo)體，需根據(jù)解析式確定共形參數(shù)，因此在電場和磁場的每一次迭代計算中都需要判斷x、y、z三個方向共形網(wǎng)格的位置，即異常體表面與網(wǎng)格相交位置，以及交點所分相應(yīng)棱邊的內(nèi)外邊長.并將共形參數(shù)代入迭代式中進行計算，因此與三棱柱模型相比要更耗時，與階梯狀網(wǎng)格近似相比，其時間比值分別為1.90倍和1.86倍.由此可知，共形網(wǎng)格技術(shù)可在保證計算效率的同時提高計算精度.

圖14 橢球形異常體電磁響應(yīng)切片圖Fig.14 Electromagnetic response slice diagram of ellipsoidal anomalous body

圖15 t=7 ms時橢球形異常體感應(yīng)電動勢等值線分布圖Fig.15 Contour distribution of ellipsoidal anomalous electromagnetic response at t=7 ms

表1 共形網(wǎng)格技術(shù)及階梯狀網(wǎng)格近似計算效率對比表Tab.1 The calculation efficiency comparison of conformal technology and stepped grid approximation

5 結(jié)論

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

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

2）可通過三維坐標(biāo)投影法確定出共形網(wǎng)格位置，進而判斷出所截網(wǎng)格棱邊位置，以確定出曲面共形參數(shù).

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

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