基于COMSOL的起伏地形三维大地电磁数值模拟

李忠乾，李长伟,2，宋卫文，周 翔，伍师莹，熊 彬

(1.桂林理工大学地球科学学院，广西桂林 541006；2.广西隐伏金属矿产勘查重点实验室，广西桂林 541006)

0 引言

大地电磁是一种天然源低频电磁勘探方法，广泛应用于资源勘探、工程勘察，地质调查、环境和水资源调查，以及深部地质结构探测等领域(付光明等,2019；仇根根等,2020)。数值模拟对于大地电磁的方法研究、数据解释等具有重要意义。随着计算技术、电磁勘探技术的发展，电磁法的数值模拟已经有了相当多的研究成果，简单地质模型的正演问题已经发展完善，但对于复杂地质模型而言，需要进一步探索求解的技术方法。地形是引起电磁响应畸变的一个重要因素(Nam et al.，2008；薛国强等，2016)，模拟起伏地形下的电磁响应，对地形进行校正，对于提高大地电磁的数据解释质量有重要意义。近年来，国内外有许多学者对此进行了相关研究。Han et al.(2009)对比了有限差分法与有限单元法处理三维MT问题的优缺点；Ting and Hohmann(2012)利用积分方程成功实现小规模三维地质体的正演模拟；杨波等(2012)应用有限体积法模拟了一个带地形的海底三维频率与可控源正演，得到海底地形变化对于CSEM电场分量的变化要大于磁场分量的变化；李芸芸(2014)利用有限单元法进行了起伏地形井地的数值模拟，总结了不同埋深对于低高阻异常响应影响规律；李俊杰等(2015)利用Galerkin无网格法研究了复杂电导率分布和复杂边界形状的地电模型正演；赵越等(2017)研究了三维起伏地形条件下航空瞬变电磁响应特征；Lin et al.(2017)研究了三维地形对可控源音频大地电磁数据的影响；李静和等(2018)采用了积分方程法研究了起伏地形频率域可控源电磁二维快速正反演模拟；冯德山等(2018)利用不规则四边形的插值方式实现了二维起伏地形的正演模拟；钟苏美等(2018)采用有限差分法研究起伏地形下可控源音频大地电磁三维数值模拟，提出了新的起伏地形三维正演方程的边界条件；马炳镇(2018)采用三维时间域有限差分法模拟了起伏地形地面瞬变电磁法三维正演；殷长春等(2019)采用自适应非结构有限元算法求解起伏地形时间域海洋电磁三维正演；Ansari et al.(2020)利用多种软件将真实地质体构建出来，进行四面体网格剖分，实现了三维大地电磁正演；顾观文和李桐林(2020)研究了基于矢量有限元的带地形大地电磁三维正演；严家斌等(2020)研究了基于电场矢量有限元三维高频大地电磁倾子响应；Li et al.(2020)基于改进的有限元算法实现了大地电磁的三维正演；万文武等(2020)基于谱元法实现了低频平面电磁波电磁场数值模拟。在上述这些数值模拟方法中，对于复杂模型，积分方程法在求解格林函数时存在困难；有限差分法实现简单，但无法用非结构化网格进行剖分；有限体积法难以实现高阶精度；有限单元法可采用非结构化网格剖分，能灵活地进行局部网格加密，求解过程规范，对于复杂地质地形的模拟有天然优势。

在地球物理电磁法正演方面，目前还很少见到界面友好、功能足够强大、通用化强的专业软件。COMSOL是一个基于有限元方法的多物理场仿真软件，有一些学者尝试利用COMSOL进行地球物理电磁法数值模拟。如王小龙等(2011)以点电源和线电源为例研究了基于的COMSOL的直流电法正演；Butler and Sinha(2012)利用COMSOL实现了电阻率法、激发激化法、电磁法等地球物理正演；冯硕等(2013)实现了基于COMSOL的井地电阻率法正演；陈小燕(2014)研究了COMSOL软件的大地电磁正演模拟；唐杰和熊彬(2018)利用COMSOL，基于非结构化网格剖分，实现了二维大地电磁正演；李付龙等(2019)基于COMSOL实现了大地电磁二维正反演；郝耀军等(2020)将COMSOL应用于掘进巷道直流电法超前探测研究。上述这些成果大都是二维模拟，而三维模拟以及考虑起伏地形的模拟成果比较少。本文研究了基于COMSOL的起伏地形大地电磁三维正演，将起伏地形参数化为曲面，建立几何结构，赋予结构体电性，进行非结构化网格剖分，求解线性方程组得到电磁响应，利用模拟结果进行地形校正，获得了较好的效果。

1 大地电磁三维有限元正演

设角频率为ω，时间因子为e-iωt，大地电磁正演的基本方程为：

∇×E=-iωμH

(1)

∇×H=(σ-iωε)E

(2)

其中H为磁场强度(A/m),E为电场强度(V/m)，ω是频率(rad/s)，μ是介质的磁导率(H/m)，σ是导电率(S/m)，ε是介电常数(F/m)。

将(1)式代入(2)得到电场E满足的微分方程：

(3)

其中k2=iωμσ+ω2εμ≈iωμσ。

设初始大地电磁场是平面波场，初始电场E的极化方向沿X轴，选取足够大的六面体区域，三维不均匀体产生的异常电磁场在边界区域Γ上为0，如图1所示(徐世浙，1994)，电磁场的边界条件是：

图1 三维大地电磁模型示意图Fig.1 Schematic diagram of three dimensional magnetotelluric model

(1)ABCD面上，

Ex=1,Ey=0,Ez=0

(4)

(2)在四个垂直边界上，电磁场的传播方向垂直向下，与边界面的法向垂直，即

E×H⊥Γ

② 网点：网点清晰，角度准确，不出重影。印刷品50%网点的增大值范围一般为10%～20%，最大不超过25%。线划粗细变形率不超过15%。

(5)

(3)在EFGH面上，电磁场按指数规律向下传播：

Ex=ce-kz,Ey=0,Ez=0

(6)

在求解出电场分布后，利用法拉第电磁感应定律求得磁场:

(7)

视电阻率和相位由下面的公式定义：

(8)

(9)

2 COMSOL Multiphysics建模

COMSOL是以有限元法为基础，通过求解偏微分方程来实现真实物理现象的仿真，其强大的内插函数以及参数化曲面能够精确地构建出复杂几何结构，并且可以灵活实现局部网格精化，广泛应用于各个领域的科学研究以及工程计算。COMSOL中物理场有多个模块可供选择，我们选取AC/DC板块来研究大地电磁正演。在COMSOL(本文使用版本为COMSOL Multiphysics 5.6试用版)中构建出复杂地电模型，然后进行网格剖分及方程组求解，具体实现流程如图2所示。

图2 COMSOL仿真流程图Fig.2 Flow chart of COMSOL simulation

为建立起伏地形，将高程数据导入到COMSOL建立模型，如图3所示。设置计算区域为长方体，将起伏地表进行参数化(张亮等，2014；Butler and Zhang，2016；何紫兰等，2020；刘江涛等，2021)，利用参数化将长方体截取为起伏地形模型，删除多余实体(图4)，再进行非结构化网格剖分，如图5所示。

图3 起伏地形示意图Fig.3 Sketch showing undulating terrain

图4 参数化曲面Fig.4 Parametized surface

图5 非结构化网格剖分Fig.5 Unstructured grid

3 验证COMSOL仿真正确性以及频率的选择

为验证COMSOL仿真的正确性和精度，我们采用国际标准测试模型COMMEMI 3D-1进行测试。该模型包括电阻率为100 Ω·m的围岩，及其中赋存的一个电阻率为0.5 Ω·m的低阻体，低阻体模型大小为1 km×2 km×2 km，埋深为0.25 km，模型总尺寸为4 km×4 km×4 km。在XY模式中设置上界面的磁场强度Hy为1000 A/m；在YX模式下设置上界面的磁场强度Hx为1000 A/m。取X轴方向为测线方向，取0.1 Hz下XY模式的COMSOL模拟结果与IE方法计算出的视电阻率进行对比，如表1和图6所示。两种方法的结果比较表明，COMSOL仿真结果与IE的计算结果相吻合，二者的最大相对误差在4.29%，说明利用COMSOL软件进行大地电磁模拟是可行的，可以进行下一步复杂地形的大地电磁正演。为了后续选取一个合适频率来进行研究，在求解国际标准测试模型COMMEMI 3D-1时，频率分别取1000 Hz、100 Hz、10 Hz、1 Hz、0.1 Hz、0.01 Hz、0.001 Hz，得到的正演结果如图7所示。由图7a、b可以看出，无论是XY模式还是YX模式，随着频率的降低，在x为[-5 km,5 km]处低阻体的异常响应变得越来越明显，但是当频率降低到0.001 Hz时，在曲线的两侧视电阻率值明显要小于背景场电阻率，所以我们在进行频率值的选取时不能一味追求低频率。表2为各个频率下正演所需时间、内存以及自由度，电脑配置为Intel(R) Core(TM) i7-6005U CPU 2.5G Hz RAM-8GB。为保证结果准确，在频率处于[10 Hz,1000 Hz]情况下网格加密，其自由度非常大，所以计算时间久，占用内存也很高。在频率[0.001 Hz,1 Hz]区间，网格剖分密度一致，所占物理内存与虚拟内存较为接近。在内存够用的情况下，应优先选择耗时最短的频率，由表2可见0.01 Hz时耗时最短。综上，后续研究中均取0.1 Hz、0.01 Hz两个频率进行研究。

图6 IE的结果和COMSOL的模拟值对比Fig.6 Comparison of IE results and COMSOL simulation

表1 IE的结果和COMSOL的模拟值对比Table 1 Comparison of IE results and COMSOL simulation values

图7 多个频率视电阻率曲线Fig.7 Apparent resistivity curves of multiple frequenciesa-XY模式；b-YX模式a-XY mode；b-YX mode

表2 各个频率正演所需时间、内存及自由度Table 2 Time，memory and degrees of freedom required by forward modeling of each frequency

4 模型算例

模型一：为研究起伏地形对正演模拟结果的影响，首先构建一个不含异常体的起伏地形模型，地层的电阻率值设置为 100 Ω·m，如图8a所示。设置计算区域为70 km×70 km×60 km，对计算区域进行非结构化网格剖分，如图8b所示。各段测线在X轴方向两端端点分别为：(-35km，0，-13.5 km)，(-25 km，0，-11 km)；(-25 km，0，-11 km)，(-18 km，0，-9 km)；(-18 km，0，-9 km)，(-10 km，0，-13 km)；(-10 km，0，-13 km)，(-4 km，0，-17 km)；(-4 km，0，-17 km)，(-2 km，0，-17.5 km)；(-2 km，0，-17.5 km)，(2 km，0，-15.5 km)；(2 km，0，-15.5 km)，(12 km，0，-9 km)；(12 km，0，-9 km)，(17 km，0，-10 km)；(17 km，0，-10 km)，(27 km，0，-16 km)；(27 km，0，-16 km)，(35 km，0，-19 km)。

图8 无异常体起伏地形模型Fig.8 Undulating terrain model without anomaly bodya-无异常体起伏地形模型；b-沿测线非结构化网格剖分切面a-an undulating terrain without abnormal body；b-section of unstructured grid along survey line

模型二：将异常体加入模型，模型包含一个电阻率值为1000 Ω·m的高阻异常体(图9a右)，一个电阻率值为0.5 Ω·m的低阻异常体(图9b左)，围岩的电阻率值为100 Ω·m，两个异常体模型大小设置为10 km×10 km×10 km，埋深为0.25 km ，高阻体上方为山谷地形，低阻体上方为山峰地形，如图9a所示。计算区域设置为70 km×70 km×60 km，取空气层的电阻率为109Ω·m，对计算区域进行非结构化网格剖分，如图9b所示。测线设置与模型一相同。

图9 起伏地形低高阻异常体模型Fig.9 Low and high resistance anomaly body model under undulating terraina-起伏地形低高阻异常体模型；b-沿测线起伏地形非结构化网格剖分切面a-low-high resistance anomaly for undulating terrain；b-section of unstructured grid along survey line with undulating terrain

取两个频率f=0.1 Hz、f=0.01 Hz，进行正演计算，视电阻率曲线如图10所示。通过图10a，10c可以看出无异常体平坦地形的视电阻率曲线与背景场吻合。由图10b可以看出在XY模式中，视电阻率曲线与地形呈镜像关系。图10(d)YX模式对地形反应不明显。通过对比图10e和10f视电阻率曲线，可以看出在XY模式下，起伏地形在X在[-15 km，-5 km]处引起视电阻率曲线振荡，在极大值处视电阻率值降低。对比图10g和10 h视电阻率曲线可以看出:在YX模式下，高阻异常体上方地形使电阻值减小，起伏地形影响相对XY模式不明显。

图10 平坦地形视电阻率曲线(a) (c)(e)(g)、起伏地形视电阻率曲线(b)(d)(f)(h)Fig.10 Apparent resistivity curves (a)，(c)，(e) and (g) for flat terrain，apparent resistivity curves (b)，(d)，(f) and (h) for undulating terraina-无异常体平坦地形视电阻率曲线XY；b-无异常体起伏地形视电阻率曲线XY；c-无异常体平坦地形视电阻率曲线YX；d-无异常体起伏地形视电阻率曲线YX；e-低高阻平坦地形视电阻率曲线XY；f-低高阻起伏地形视电阻率曲线XY；g-低高阻平坦地形视电阻 率曲线YX；h-低高阻起伏地形视电阻率曲线YXa-apparent resistivity curveof flat terrain without abnormal body in XY mode;b-apparent resistivity curve for undulating terrain without abnormal body in XY mode;c-apparent resistivity curve of flat terrain without abnormal body in YX mode;d-apparent resistivity curve of undulating terrain without abnormal body in YX mode;e-apparent resistivity curveof flat terrain with low and high resistance in XY mode;f-apparent resistivity curve for low and high resistance relief terrain in XY mode;g-apparent resistivity curve of flat terrain with low and high resistance in YX mode;h-apparent resistivity curve for low and high resistance relief terrain in YX mode

图11是有地形和无地形情况下计算得到的视电阻率和相位等值线图。由图11可以看出，在XY模式下，起伏地形引起视电阻率等值线出现条带状分布，在YX模式下地形对视电阻率的分布影响相对不明显。相位则是在YX模式下出现条带状异常，XY模式下不明显。

图11 平坦地形视电阻率和相位等值线(a)(c)(e)(g)、起伏地形视电阻率和相位等值线(b)(d)(f)(h)Fig.11 Apparent resistivity and phase contours (a)(c)(e)(g) of flat terrain apparent resistivity and phase contours (b)(d)(f)(h) of undulating terraina-平坦地形低高阻视电阻率XY模式 ；b-起伏地形低高阻视电阻率XY模式；c-平坦地形低高阻相位XY模式；d-起伏地形低高阻相位XY模式；e-平坦地形低高阻视电阻率YX模式 ；f-起伏地形低高阻视电阻率YX模式；g-平坦地形低高阻相位YX模式；h-起 伏地形低高阻相位YX模式a-low and high apparent resistivity in flat terrain XY mode；b-low and high apparent resistivity in undulating terrain in XY mode；c-low and high resistance phase in flat terrain in XY model；d-low and high resistance phase in undulating terrain in XY mode；e-low and high apparent resistivity in flat terrain in YX mode；f-low and high apparent resistivity in undulating terrain in YX mode；g-low and high resistance phase in flat terrai YX model；h-low and high resistance phase in undulating terrain in YX model

5 地形校正

由前面的算例可以看出，起伏地形会引起视电阻率畸变，给数据解释造成干扰。考虑利用COMSOL的模拟结果进行地形校正，采用常规的 “比值法”(李金铭，2005)，即计算单一地形因素影响下的视电阻率与背景场视电阻率的比值，再用有异常体起伏模型的视电阻率除以这个比值，得到校正后的视电阻率。

取对数后写为

为验证方法的效果，构建山峰、山谷模型，求解区域设置为4 km×4 km×6 km，围岩的电阻率值为100 Ω·m，空气层的电阻率设为109Ω·m，异常体的电阻率设置为0.5 Ω·m，异常体模型大小为1 km×1 km×2 km，如图12所示。山峰测点坐标为(-2 km，0，0.1 km)，(-1.3 km，0，-1 km)，(-0.9 km，0，-1.5 km)，(-0.3 km，0，-1.8 km)，(0.1 km，0，-1.7 km)，(0.5 km，0，-1.4 km)，(1 km，0，-0.85 km),(2 km，0，0.2 km)。山谷测点坐标为(-2 km，0，-1.9 km)，(-0.7 km，0，-0.4 km)，(-0.3 km，0，-0.15 km)，(0，0，0)，(0.5 km，0，-0.2 km)，(2 km，0，-1.7 km)。

图12 山峰与山谷模型Fig.12 Peak and valley modelsa-山峰；b-山谷a-peak；b-valley

图13a，b分别为XY模式下山峰、山谷地形中异常体埋深250 m时地形校正情况。由图13可以看出，起伏地形会影响视电阻率曲线的形态，曲线两侧远端不再趋向于背景值，异常体附近视电阻率曲线出现跳跃，经过地形校正后曲线形态光滑，远端趋于背景值，异常体上方极值点更为明显。

图13 地形校正对比图Fig.13 Comparison of terrain correctiona-XY模式下埋深250 m山峰地形校正对比图；b-XY模式下埋深250 m山谷地形校正对比图a-comparison of topography correction of mountain peak at 250 m burial depth in XY mode ；b-comparison of topography correction of valley at 250 m burial depth in XY mode

6 结论

本文利用COMSOL实现了带地形MT三维正演模拟，将地形起伏数据转化为参数化曲面，再利用COMSOL进行正演模拟，将起伏地形的模拟结果与无地形起伏时的结果进行对比，总结出山谷、山峰地形对三维大地电磁响应的影响，并利用“比值法”进行地形校正，得出如下结论：

(1)本文进行了COMSOL软件模拟值与IE结果对比，证明利用COMSOL进行大地电磁三维正演是可靠的。

(2)通过数值算例分析可以看出，当X方向取为垂直地形走向时，XY模式下，起伏地形引起明显的视电阻率曲线和等值线图的畸变，对相位的影响较小；在YX模式下，起伏地形会引起相位的畸变，对视电阻率分布的影响相对不明显。

(3)通过对比地形校正前后，发现经过地形校正的数据在异常体远侧回归至背景场值，视电阻率曲线形态简单化，由异常体引起的极值更为明显。因此在地形复杂区域开展勘探工作时，采用软件模拟结果进行比值法地形改正应是可供选择的方法之一。

[附中文参考文献]

陈小燕.2014.基于COMSOL软件的大地电磁正演模拟[D].西安:长安大学：56-64.

冯德山,刘金宝,王珣.2018.含地形起伏的MT模型非规则四边形FEM正演与反演[J].中南大学学报(自然科学版),49(3):626-632.

冯硕,刘得军,张颖颖,秦民君,兀凤娜.2013.基于COMSOL的井地电阻率正演研究[J].断块油气田,20(5):589-592.

付光明,黄进调,刘阳,李晓兵,陈国玉,黄艺.2019.高密度电阻率法与CSAMT法在江西会昌县坝背地热勘查中的综合探测[J].中国地质,46(4):927-936.

顾观文,李桐林.2020.基于矢量有限元的带地形大地电磁三维反演研究[J].地球物理学报，63(6):2449-2465.

郝耀军,李浩,武艺.2020.基于COMSOL掘进巷道直流电法超前探测研究[J].煤矿现代化,(6):104-106.

何紫兰,朱鹏飞,白芸,曹珂,孔维豪.2020.复杂地质体三维实体建模方法[J].地质与勘探,56(1):190-197.

李付龙,李清坤,李文帝,何好,尚延杰,李维强.2019.基于Comsol仿真软件的大地电磁法拟二维反演研究[J].勘察科学技术,(4):25-28.

李静和,何展翔,熊彬.2018.起伏地形频率域可控源电磁二维快速正反演[J].地质与勘探,54(2):325-331.

李金铭.2005.地电场与电法勘探[M].北京：地质出版社:187-189.

李俊杰,严家斌,皇祥宇.2015.无单元Galerkin法大地电磁三维正演模拟[J].地质与勘探,51(5):946-952.

李芸芸.2014.起伏地形井地电位三维有限元数值模拟[D].长沙:中南大学：27-49.

刘江涛,赵洁,吴发富.2021.结构方程模型及其在地学数据建模中的回顾与展望[J].地质力学学报,27(3):350-364.

马炳镇.2018.起伏地形下地面瞬变电磁法三维正演数值模拟研究[J].物探与化探,42(4):777-784.

仇根根,方慧,吕琴音,彭炎.2019.武夷山北段及相邻区深部电性构造与成矿分析：基于三维大地电磁探测结果[J].中国地质,46(4):775-785.

唐杰,熊彬.2018.基于COMSOL软件的非结构化网格下的二维大地电磁正演模拟[J].工程地球物理学报,15(3):347-356.

万文武,李长伟,熊彬,高磊,宋卫文.2020.基于谱元法的低频平面波电磁场数值模拟[J].桂林理工大学学报,40(4):703-712.

王小龙,冯宏,田华光,吴迪.2011.基于COMSOL Multiphysics的直流电法正演模拟[J].煤田地质与勘探,39(5):76-80.

徐世浙.1994.地球物理中的有限单元法[M].北京:科学出版社:247-248.

薛国强,闫述,陈卫营.2016.电磁测深数据地形影响的快速校正[J].地球物理学报,59(12):4408-4413.

严家斌,胡涛,林旭,赵玄,皇祥宇.2020.基于电场矢量有限元三维高频大地电磁倾子响应与感应矢量研究[J].地质与勘探，56(1):123-136.

杨波,徐义贤,何展翔,孙卫斌.2012.考虑海底地形的三维频率域可控源电磁响应有限体积法模拟[J].地球物理学报,55(4):1390-1399.

殷长春,惠哲剑,张博,刘云鹤,任秀艳.2019.起伏海底地形时间域海洋电磁三维自适应正演模拟[J].地球物理学报,62(5):1942-1953.

张亮,姚磊华,王迎东.2014.基于COMSOL Multiphysics的三维地质建模方法[J].煤田地质与勘探,42(06):14-19.

赵越,李貅,王祎鹏,郭建磊,曾友强.2017.三维起伏地形条件下航空瞬变电磁响应特征研究[J].地球物理学报，60(1):383-402.

钟苏美,林昌洪,谢裕春.2018.起伏地形下可控源音频大地电磁三维数值模拟[J].现代地质,32(2):398-405.