基于共形网格技术的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.