复杂地形瞬变电磁三维正演模拟与地形效应分析

屈少波，朱 姣，姜志海，李毛飞，何治隆

(1.中国矿业大学 矿业工程学院，江苏 徐州 221116；2.山西宁武榆树坡煤业有限公司，山西 忻州 036700；3.中国矿业大学 资源与地球科学学院，江苏 徐州 221116)

煤矿突水灾害精准超前预报预警一直是煤矿防治水工作追求的目标。常用的矿井瞬变电磁法通常在巷道内施工，采用多匝小回线作为发射线圈，由于自感、互感现象严重，并且巷道内金属干扰过多，导致矿井瞬变电磁信号受干扰严重[1-2]。近年来学者们将目光转移到传统的地面瞬变电磁法。瞬变电磁法即为典型的时间域电磁法，可采用不接地回线和接地导线作为发射源，即磁性源与电性源系统[3-4]。我国地域广袤，山地、丘陵、高原等复杂地形众多，而地形是影响电磁场分布的主要因素之一[5-7]，基于平坦地形的瞬变电磁数据处理与解释方法与实际施工勘探具有较大的差异，大大降低了其结果的可信度，所以将地形的影响纳入瞬变电磁法的研究中具有现实意义[8-11]。

地形对短偏移距瞬变电磁法[12](Short Offset Transient Electromagnetic Method，SOTEM))和长偏移距瞬变电磁法(Long Offset Transient Electromagnetic Method，LOTEM)的探测效果均有较大的影响，通常接收点处地形对观测信号的影响大于发射源处，并且地形对电场E 的影响更为严重[13]。除地面装置外，航空瞬变电磁系统中的地形效应也广泛存在，其对接收信号的影响程度与电性参数、飞行高度与轨迹等因素密切相关[14-15]。针对地形效应，一方面可以对装置改进，如采用矩形小回线或地面-巷道装置等[16]，另一方面可对数据（如感应电动势或视电阻率参数等）进行地形校正等处理[17]。总体而言，地形效应的特征与去除仍然是瞬变电磁法勘探的重要课题，也是矿井中精确定位隐蔽突水致灾源的决定性因素之一。

三维电磁数值模拟是研究地形影响的有效途径，也是三维反演的基础[18-19]。随着计算机硬件设备与软件性能的提升，瞬变电磁法的研究已经从一维走向三维[20-22]。目前主流的三维数值模拟方法主要包括积分方程法、有限差分法以及有限单元法等[23-27]。其中，四面体网格的非结构有限元法，可对复杂地形的细节变化进行精细刻画，拥有巨大优势[28-30]。瞬变电磁数值模拟大致分为两种思路：一种是在频域计算出电磁场分布，接着由频时转换求得时域电磁场[31]；一种是直接在时域求解Maxwell 方程组[32]，可有效减小频时转换带来的误差，从而提高数值模拟精度。本文采用第二种策略。本文基于矢量有限元法，以典型的山峰、山谷模型为例分别讨论电性源和磁性源装置系统发射源和接收测点处受地形影响的规律与原因。进而利用煤层采掘工程平面图获得实际地形数据，并以此建立复杂地形三维正演模型，最后通过数据分析得出起伏地形的瞬变电磁响应特征。

1 矢量有限元正演模拟理论

从麦克斯韦方程出发，利用伽辽金加权余量法，w以作为权函数，则单元电场扩散方程的残差Re可表示为：

将计算区域内所有单元合并成全局系统，令加权残差R=0，表示为矩阵形式：

对式(2)中的时间导数项采用二阶后退欧拉格式离散[13]。这样，只需求解式(2)中的线性方程就可得到计算区域所有棱边的电场值，进而得到接收点位置的电磁场参数。

2 理论模型三维正演模拟

为探究起伏地形对地面瞬变电磁响应的影响规律，本文采用有限元三维正演技术直接在时域求解Maxwell方程组得到电磁场，获得起伏地形条件下地面瞬变电磁三维正演响应。综合考虑实际地质环境，野外地形多数可简化为连绵起伏的山峰山谷构造，因此，本文在三维模型的设计过程中，选择山峰、山谷这两类较为典型的起伏地形模型进行地面瞬变电磁法的数值模拟。图1为山峰和山谷模型的参数和正演网格，对发射源及接收测点附近进行局部网格加密。

图1 理论模型与网格示意Fig.1 Theoretical models and their grids

2.1 电性源瞬变电磁系统正演模拟

发射源和接收测点分别位于地形上(以山峰为例)的电性源瞬变电磁模型如图2 所示。计算域范围为3 km×3 km×0.6 km，背景电阻率为100 Ω·m，激励源为600 m 的长导线，发射电流为1 A。测线分布范围：x=[-250 m,250 m]，y=[-250 m,250 m]，线距与点距均为25 m，共计21×21=441 个测点，采样时间为10-6～10-2s，共有对数等间隔采样点31 个。山峰、山谷地形底半径为150 m，谷底与山顶高度均为50 m。

图2 电性源系统发射源或接收测点在山峰上Fig.2 Electrical-source TEM models with the transmitting source or receiving point located at a mountain peak

2.1.1 接收测点在起伏地形上

当发射源位于水平地表，断电后整个观测区域内不同时刻(10-6、10-5、10-4、10-3、10-2s)的感应电动势(dBz/dt)剖面如图3 所示。水平地表条件下的响应如图3a-图3e 所示，由电磁感应原理可知，供电导线中的电流突然断开后，会在附近产生与导线源同向的感应电流，同时在四周形成回流电流。随着时间的推移，感应电流开始向下扩散，并且范围逐渐扩大。由此，导线两侧形成2 个响应中心，其中心位置随着时间逐渐向远离导线方向移动。

图3 不同时刻水平地表、接收测点在山峰/山谷模型的感应电动势(dBz/dt)剖面Fig.3 Profiles showing the induced electromotive force (dBz/dt) at different times in the horizontal surface model and models with the receiving point at a mountain peak/valley

当发射源位于水平地表，接收点在山峰(图3f-图3j)和山谷(图3k-图3o)地形上时，观测响应发生了改变。早期时间道，山峰表现为“吸引作用”，出现一个响应高值中心(图3f)，同时在峰顶另一侧出现反号现象。随后，反号现象迅速消失，但山峰导致的正响应中心依然存在(图3g)。当响应中心经过地形时，山峰表现为“排斥作用”，形成一个响应低值中心(图3h-图3j)。山谷的响应规律则与山峰地形完全相反，早期响应特征表现为“排斥作用”(图3k)，同样存在变号现象，中晚期响应表现为“吸引作用”(图3m-图3o)。

图4 为单条测线上各测点的衰减曲线。水平地表模型中，x=0 m 测线(图4a)上的测点沿中心呈对称分布，两侧测点响应基本重合。而y=0 m 测线(图4d)上的测点随着偏移距的增大早期响应数值逐渐减小。山峰地形影响主要集中在早期时间道，并且山峰顶部测点dBz/dt曲线上凸，但在山峰底部测点曲线则呈下凹趋势。另外，考虑偏移距的影响将其与水平地表模型对比可发现，在山峰范围内右侧测点产生反号现象(图4e)。山峰区域测点由左到右响应数值先增大再减小，直至出现正负交界面，随后负响应数值也先增大后减小。山谷模型响应与山峰模型相反，二者正演结果与图3 相吻合。

图4 水平地表、接收测点在山峰/山谷模型电磁响应衰减曲线Fig.4 Electromagnetic response attenuation curves derived from the horizontal surface model and models with the receiving point at a mountain peak/valley

由此可见，地形效应的产生包括多种因素：山峰相当于水平地表上方的一个低阻体，山谷则可视为埋藏水平地表下方的一个高阻体。突然断电后，山峰会在两侧产生与半空间的主电流相反的散射电流场，早期过强的散射电流场导致感应电压的变号现象。在山峰中表现为靠近源一侧产生同向的感应电压，远离源导线的一侧产生反向的感应电压[13]。同样的道理，山谷的响应也存在变号现象，源导线两侧感应电压的方向与山峰恰好相反。其次，由电流密度分布可知，早期低阻体内激励的涡旋场较大，使得接收到的感应电动势随之增大，反之高阻体中接收的感应电动势则减小。山峰位置对应响应信号增强，山谷对应响应信号减弱，表现为响应参数曲线的变化与地形起伏形态一致。另外，频率较高的电磁信号穿透深度较浅，感应涡流产生的二次场在较早的时刻就被接收到，因此，早期瞬变电磁信号中高频成分占主导，主要反映较浅位置的地下介质，地形的特征在早期能够被捕捉到。反之，晚期瞬变电磁信号中低频成分占主导，主要反映较深位置的背景介质，受地形影响较小。最后，起伏地形与深部结构之间存在复杂的耦合效应，导致接收信号发生改变。

2.1.2 发射源在起伏地形上

图5 为源在地形上，接收点位于水平地表时整个观测区域内不同时刻电磁扩散剖面。每列时间相同，即10-6、10-5、10-4、10-3、10-2s。可以看出发射源位于地形上对dBz/dt响应等值线形态并没有明显影响，但是响应幅值增大。其中山峰地形使得早期响应增强(图5f-图5j)，而山谷地形则削弱了响应(图5k-图5o)。同样的，测点的衰减曲线具有相同规律(图6)，同模型测点电磁响应衰减曲线形态基本一致，拐点位置也基本相同。区别在于源在山峰和源在山谷模型的响应相较水平地表型(图5a-图5e)更大，曲线整体上移。这是由于发射源在地形之上，相对于水平地表二者的相对高程和电磁波传播路径发生改变。发射源与接收测点的布置均是以水平距离为准，实际的发射源导线长度大于水平距离，响应曲线整体得到抬升。另外，山峰和山谷地形在早期响应增强的幅值也有所差异，是由于山峰地形与导线源的距离虽然增大，但空气作为中间介质对电磁波吸收更弱，早期信号相对增强。而山谷地形与导线源的距离同样增大，但中间介质为导电大地，早期响应幅值降低。

图5 不同时刻水平地表、发射源在山峰/山谷模型的电磁扩散剖面Fig.5 Profiles showing electromagnetic diffusion at different times in the horizontal surface model and models with transmitting source at mountain peak/valley

图6 水平地表、发射源在山峰/山谷模型电磁响应衰减曲线Fig.6 Electromagnetic response attenuation curves derived from the horizontal surface model and models with the transmitting source at a mountain peak/valley

2.2 磁性源瞬变电磁系统正演模拟

磁性源系统模型与电性源相似，发射源是长宽均为420 m 的正方形回线，坐标满足x=[-210 m,210 m]，y=[-210 m,210 m]，发射电流为1 A。采用面积性观测方式进行测量，测线分布为：x=[-300 m,300 m]，y=[-300 m,300 m]，点距25 m，共25×25=625 个测点。设计山峰模与山谷模型如图7 所示，将源或接收测点置于地形之上。

图7 磁性源系统发射源导线或接收测点在山峰上Fig.7 Magmatic-source TEM models with the transmitting source or receiving point at a mountain peak

2.2.1 接收测点在起伏地形上

图8 为源在水平地表时断电后整个观测区域内不同时刻(10-6、10-5、10-4、10-3、10-2s)电磁扩散剖面。在水平地表情况下(图8a-图8e，图9a)早期涡流紧挨发射回线，响应与发射回线形状相同，随着电流环向外、向下扩散，逐渐变为圆形扩散。由于回线中心位置磁力线垂直，垂直磁场分量最大。当接收测点在山峰(图8f-图8o)的响应扩散剖面等值线形态相似。区别在于接收点在山峰上时，回线中心位置测点处的早期响应较水平地形表现为响应“增强”的特征，山峰越高信号强度越大；中、晚期信号较水平地形表现为响应“减弱”的特征，山峰越高信号强度越小(图8f-图8j，图9b)。接收点在山谷上时响应变化则与山峰相反(图8k-图8o，图9c)。磁性源系统的发射线圈可以视为4 条首尾相接的电性源导线，磁性源系统地形影响产生的原因与电性源相同，也是由于收发装置相对高程的改变和地形所引起的介质电性差异造成。

图8 不同时刻水平地表、发射源在山峰/山谷模型电磁扩散剖面Fig.8 Profiles showing electromagnetic diffusion at different times in the horizontal surface model and models with the transmitting sources at a mountain peak/valley

图9 水平地表、接收测点在山峰/山谷模型电磁响应衰减曲线(y=0 m 测线)Fig.9 Electromagnetic response attenuation curves (y=0 m) derived from the horizontal surface model and models with the receiving point at a mountain peak/valley

2.2.2 2 条发射源导线均在起伏地形上

图10 为源导线在地形上时断电后整个观测区域内不同时刻电磁扩散剖面。图10a-图10e 为水平地表模型；图10f-图10j 中发射源左侧和右侧导线均在山峰；图10k-图10o 中发射源左侧导线在山峰、右侧在山谷；图10p-图10t 中为发射源左侧导线在山谷、右侧在山峰；图10u-图10y 中发射源左侧和右侧导线均在山谷。当2 条源导线均在起伏地形上时，相较水平地表模型dBz/dt响应数值整体变大，并且电磁扩散剖面等值线形态发生改变。起伏地形模型中，山峰与山谷下方早期时间道响应数值均变小。随着时间的推移，中晚期山峰对应等值线低值中心，山谷对应等值线高值中心。由于模型对称，两侧均为山峰或山谷时，响应等值线呈对称分布。晚期等值线图山峰位置出现低值中心，山谷位置则对应高值中心。选取图10 中1 条测线，即y=0 m，且x=[-300 m,0 m]，得到测线上每个测点的电磁衰减曲线(图11)。可以看出，相较水平地表模型，2 条源导线在地形上时各模型dBz/dt响应衰减曲线形态相似，但曲线均整体上移。由于晚期响应参数的数值较小，测点衰减曲线无法显示地形的影响。

图10 不同时刻水平地表、发射源在山峰/山谷地形模型电磁扩散剖面Fig.10 Profiles showing electromagnetic diffusion at different times in the horizontal surface model and models with the transmitting source at a mountain peak/valley

图11 水平地表、发射源在山峰/山谷模型电磁响应衰减曲线(y=0 m 测线)Fig.11 Electromagnetic response attenuation curves (y=0 m) derived from the horizontal surface model and models with the transmitting sources at a mountain peak/valley

综上所述，无论是磁性源还是电性源系统，地形对瞬变电磁法勘探的影响均较大，尤其是早期信号。在野外实测数据处理中需要充分考虑地形因素，复杂地形模型的整体响应结果可以看作是各简单地形的组合，本文主要讨论的山峰与山谷地形可为复杂模型的研究提供参考。

3 基于实际地形的电磁正演模拟

山西宁武县榆树坡煤矿实际地貌如图12 所示，地表为黄土覆盖，经长期冲刷切割，沟谷纵横、梁峁绵延，地形较为复杂。俯瞰工区整体地势(图13a)为中部高东西两侧低、南高北低，东部发育一南北向沟谷。梁垣近东西、东南向连绵展布，沟谷多呈西北向延伸。图13b为采掘工程平面图，可将图13b 中等高线数据进行抽样插值并转化为二维等值线(图13c)，以获得三维起伏地表参数(图13d)，并据此建立三维地形模型进行正演计算。

图12 测区野外地貌Fig.12 Field landforms in the study area

图13 实际地形三维正演建模步骤Fig.13 Steps for 3D forward modeling of the complex field terrains in the study area

研究区地势最高点位于东北部，最大相对高差达130 m。在三维地形数据的基础上设置观测和接收装置进行起伏地形正演建模，如图14、图15 所示。模型计算域为1 700 m×1 700 m×1 500 m，扩边域为80 km×80 km×80 km。

图14 野外实际地形建模Fig.14 Modeling of field terrains in the study area

图15 三维正演模型及其网格Fig.15 3D forward model and its grids

研究区域地形崎岖复杂，为提高探测分辨率，克服地形影响，同时采用电性源装置和大回线源装置，共设计3 组模型(表1)：(1) 电性源：源导线长600 m，共4 条平行的接收测线，点距10 m，每条测线共计51 个测点；(2) 磁性源A 与磁性源B：发射线框400 m×600 m，共4 条接收测线，点距10 m，每条测线共有26 个测点。4 条测线的切片如图16 所示，源导线和接收点附近网格均进行局部加密。由图16 可知，4 条测线地形起伏较为相似，均呈现出左高右低的特征，并且地形面网格较为精细，随着深度增大网格逐渐粗糙。

表1 地形模型装置参数Table 1 Parameters for devices in the terrain models

图16 起伏地形模型中4 条测线网格切面Fig.16 Grid sections of four survey lines in the rugged terrain model

本文地形起伏模型考虑的是源导线完全贴合地形的情况，需要在剖分完毕的网格文件中搜索源导线的坐标与高程，并形成多条首尾相连的导线作为正演模拟的输入信息文件。此模型中，电性源导线由124 条线段组成，磁性源A 和磁性源B 的线圈共有452 条线段。采用阶跃波作为发射波形，采样时间为10-6～10-2s，共有对数等间隔采样点31 个，电流为1 A。空气电阻率为108Ω·m，背景电阻率为100 Ω·m。

4 结果与分析

下文将对这2 种装置系统的响应结果进行详细讨论。野外实际地形的响应分析中将时间参数分为两段，即t1：1×10-6～3.981×10-5s 和t2：5.412×10-5～1×10-2s。多测点、多时间道系统响应图中，t1和t2时间的响应分开绘图示意，而单点、多时间测点衰减曲线中，以蓝色色块标识t2时间段内响应。

4.1 电性源瞬变电磁系统正演响应

野外实际地形的发射源与接收点均在起伏地表上方，电性源4 条测线的多测点、多时间道系统电磁响应在图17 中展示。可以看出，t1时间段内(图17a-图17d)复杂地形使得dBz/dt曲线产生严重畸变，观测响应变得十分复杂。响应曲线与地形起伏形态基本一致，即接收点处高程大于水平地表时则响应增强，低于水平地表时则被削弱。随时间推移逐渐，响应曲线变化与地形起伏形态相反并且开始出现明显的跳变，这是由于电性源系统受地形影响，不同接收测点响应参数发生变号所导致。t2时间段内(图17e-图17h)早期仍然受到地形影响，尤其L3 和L4 测线5.412×10-5～1×10-4s 之间曲线光滑与地形起伏相反，但随着时间推移各测点响应值趋于一致。图18 中蓝色区域为t2时间段内的衰减曲线。4 测线上测点整体形态相似，t1时间测线出现变号现象，t2时间晚期响应则基本重合。

图17 实际地形模型电性源多测点、多时间道系统电磁响应Fig.17 Electromagnetic response of an electrical-source system with multiple survey points and time channels in the real field terrain model

图18 实际地形模型电性源单点、多时间道系统电磁响应衰减曲线Fig.18 Electromagnetic response attenuation curves of an electrical-source system with single survey point and multiple time channels in the real field terrain model

4.2 磁性源瞬变电磁系统正演响应

磁性源瞬变电磁法的响应曲线与电性源不同，多测点、多时间道系统电磁响应(图19、图20)中磁性源装置响应信号均为正值，不存在变号现象。磁性源A 的接收测线对应图16 中每条地形剖面的前半段(点1-点26)，磁性源B 则对应后半段(点26-点51)，可以看出t1时间内(图19a-图19d、图20a-图20d)早期响应曲线与实测地形吻合良好，晚期则相反。t2时间内(图19e-图19h、图20e-图20h)多测点的dBz/dt早期响应曲线光滑，形态与地形起伏相反，晚期则逐渐趋于水平。此外，单点、多时间道系统的衰减曲线中(图21、图22)各测点t1时间响应幅值差别较大，与地形起伏的形态规律相似，t2时间曲线则重合。

图19 实际地形模型磁性源A 多测点、多时间道系统电磁响应Fig.19 Electromagnetic response of a magnetic-source (A) system with multiple survey points and time channels in the real field terrain model

图20 实际地形模型磁性源B 多测点、多时间道系统电磁响应Fig.20 Electromagnetic response attenuation curves of a magnetic-source (A) system with single survey point and multiple time channels in the real field terrain model

图21 实际地形模型磁性源A 单点、多时间道系统电磁响应衰减曲线Fig.21 Electromagnetic response of a magnetic-source (B) system with multiple survey points and time channels in the real field terrain model

图22 实际地形模型磁性源B 单点、多时间道系统电磁响应衰减曲线Fig.22 Electromagnetic response attenuation curves of a magnetic-source (B) system with single survey point and multiple time channels in the real field terrain model

由此可得，电性源和磁性源系统在t1时间的多测点、多时间道系统曲线能够明显识别地形的起伏变化，并保留更多的细节信息。而t2时间内的早期时间道也能得出地形的大致变化。

5 结论

a.利用基于二阶后退欧拉时间离散的矢量有限元法进行理论模型数值模拟，系统研究了电性源和磁性源瞬变电磁系统在复杂地形(山峰、山谷)条件下的响应特征。电性源系统会出现响应变号现象，响应强度与地形高程呈正相关；磁性源系统无变号现象，但响应规律与电性源类似。

b.提取实际地形数据，采用四面体网格建模技术获得与实际地形高度吻合的三维正演模型。模拟结果表明：无论电性源或磁性源，早期响应曲线能够清晰反映地形细节，曲线形态与起伏地形一致；而晚期响应则逐渐趋于平缓，主要反映背景介质的影响。

c.地形效应的主要成因是由于电磁波传播路径改变导致地下介质电性分布差异。利用四面体网格正演模型能精细刻画复杂地形，模拟结果可为复杂地形区域的瞬变电磁数据采集、处理及解释提供理论指导。

d.研究结果对提高复杂地形环境下瞬变电磁法探测精度具有重要意义，是实现煤矿突水致灾源精准超前预报预警的理论基础。

致谢：本论文受国家自然科学基金、江苏省自然科学基金、江苏省“卓博计划”资助。感谢审稿专家和编辑对论文撰写的帮助及宝贵意见。

符号注释：

E(t) 为t时刻的电场；j(t) 为电流密度；J为源项矩阵；M为质量矩阵；R为一个单元内6 条棱边的权重函数；Re为第e个单元的加权残差；S为刚度矩阵；t为观测时间；V为单元体积；w为权函数；μ为磁导率；σ为电导率；Ωe为第e个单元域。