多辐射源地空瞬变电磁响应三维数值模拟研究

李貅， 胡伟明， 薛国强

1 长安大学地质工程与测绘学院， 西安 710054 2 中国科学院矿产资源研究重点实验室， 中国科学院地质与地球物理研究所， 北京 100029 3 中国科学院大学地球与行星科学学院， 北京 100049 4 中国科学院地球科学研究院， 北京 100029 5 西北有色地质矿业集团有限公司， 西安 710054

0 引言

电磁法通过观测天然或人工源激发的电磁响应获取地下电性信息(何继善和薛国强，2018;底青云等，2019).传统地面电磁法观测方式在沙漠、戈壁、山地和湿地等地区无法高效地开展工作，航空电磁法可发挥替代作用(图1).航空电磁法根据发射源布设方式可分为全航空电磁法和地空电磁法，其中全航空电磁法采集效率高，但由于受到载重限制，发射功率有限.地空电磁法将发射系统放置于地面，接收系统布设在飞行平台上(图2)，充分地结合了传统地面电磁法和全航空电磁法的优点，可实现探测深度和工作效率之间的平衡.

图1 航空电磁探测系统示意图Fig.1 Schematic diagram of airborne electromagnetic detection system

图2 多辐射源地空电磁探测系统示意图Fig.2 Schematic diagram of multi-radiation source semi-airborne electromagnetic detection system

地空瞬变电磁探测系统通常采用单一接地线源向地下发射电磁场信号，并根据地面长偏移瞬变电磁法的模式进行数据处理和解释(Mogi et al.，1998，2009；嵇艳鞠等，2013；李肃义等，2013).嵇艳鞠等(2013)采用吉林大学研发的地空瞬变电磁系统在江苏省如东县开展了海水侵蚀探测，并在内蒙古锡林郭勒盟巴彦宝利格盆地开展了地下水资源探测.Wang等(2013)引入基于多分辨分析的小波分析方法，克服飞行中的基线漂移问题.Ji等(2016)利用小波阈值法在最大限度上剔除数据中的白噪声，实现非平稳噪声的识别与剔除.李貅等(2015)利用TEM三维拟地震解释技术，为地空电磁成像解释打下基础.Wu等(2019)提出一种基于小波神经网络的噪声处理方法.

国内学者对多辐射场源下地空电磁法的数据处理与反演开展了探索性研究.张莹莹等(2015)提出一种基于多源的地空瞬变电磁方法，并给出多分量全域视电阻率定义.张莹莹等(2016)推导了多源地空瞬变电磁的视纵向电导和视深度计算公式.Zhou等(2016a,b)解决了多源发射系统的野外装置布设问题.Zhou等(2018)提出一种张量探测方法，以提升基于一对正交电性源的地空电磁系统对3D目标体的探测精度.本文采用三维矢量有限元法对两个不同地质体进行多辐射源地空瞬变电磁响应开展模拟研究，分析了多辐射源在不同辐射方向、不同飞行高度的电磁响应分布特征，从而推进多辐射源地空瞬变电磁理论方法和技术体系的建立.

1 瞬变电磁场三维矢量有限元法原理

1.1 边值问题

时间域麦克斯韦方程组为(李贺，2016)：

(1)

(2)

其中，E为电场强度，H为磁场强度，σ为介质电导率，μ0为介质的磁导率，Js为源电流密度.对(1)式两端取旋度：

(3)

将(2)式代入(3)式中得：

(4)

整理可得：

(5)

根据电磁场理论，在无源区的两种导电介质界面，电磁场满足如下四个边界条件：

(6)

式中，n为两种介质界面处的单位法向分量，B为磁感应强度，D为电位移矢量.

对于无穷远边界，电场和磁场在无穷远边界上的切向分量为零，即：

(7)

1.2 矢量变分方程

根据加权余量法，电场控制方程的加权余量方程为

(8)

其中，f为矢量基函数.根据矢量分析恒等式、高斯公式以及边界条件，可将(8)式写为

=0，

(9)

(9)式即为矢量有限元法的矢量变分方程.

1.3 矢量有限元法

图3 单元中节点与棱边的关系(引自：李贺，2016)Fig.3 Relationship between nodes and edges in an element (from Li H, 2016)

采用六面体单元进行网格剖分，在单元中棱边与节点的关系如图3所示.在六面体单元中，将电场切向分量赋予各个单元的棱边上.采用Whitney型插值基函数(李贺，2016)，使得插值基函数的散度为0，而旋度不等于0.因此，在计算区域存在介质不连续性时，电场强度的切向分量连续，法向分量不连续，可以有效地避免“伪解”的现象.

在(9)式中，含有电场对时间的偏导数项，采用向后差分的方式进行离散，即：

(10)

将(10)式代入(9)式，通过单元分析，可将(9)式写成矩阵形式：

(11)

其中：

(12)

·NidV，

(13)

(14)

即可求得电磁响应的垂直分量.

1.4 源的加载

图4 多辐射源加载示意图(参考：孙怀凤，2013)Fig.4 Schematic diagram of multi-radiation source loading(reference to Sun H F, 2013)

将电流密度直接施加到与电场的水平分量重合的单元棱边上，如图4所示.一般瞬变电磁场的正演均采用矩形脉冲进行激发，实际工作中，发射机不可能做到矩形脉冲发射，为了符合实际情况，在正演计算中，本文采用梯形波作为激发源.

2 多辐射源地空瞬变电磁法响应模拟

为了对比多辐射场源和单辐射场源的电磁响应，设计两个不同的三维地质模型，并利用自主开发的瞬变电磁三维模拟软件对其单辐射源和多辐射源下的电磁响应特征进行了模拟分析.地空电磁系统所采集的为垂直磁场分量或其感生电动势.

2.1 模型一模拟结果

模型一：半空间电阻率为100 Ωm；异常体尺寸为200 m×200 m×30 m，电阻率1 Ωm，顶板埋深100 m；测线位于y=0 m处，飞行高度分别为10 m、30 m、50 m；分别计算在单源和双源激发情况下的电磁响应.

(1)单一地质体、单一辐射源电磁响应特征

单辐射场源与单一三维模型立体图如图5a所示，其中源1长300 m，位于(-150，-150，0)至(150，-150，0)，供电电流15 A.沿测线(y=0 m)在不同飞行高度10 m、30 m、50 m的感生电动势多测道图分别如图5b、c、d所示，由图可知垂直分量多测道图呈单峰分布，10 m飞行高度的感生电动势幅值较大，而50 m飞行高度的感生电动势幅值较小.

(2)单一地质体、多辐射源磁场响应特征

当两个辐射源的电流方向相反时的地电模型如图6a所示，源1与以上模型一致，源2位于(150，150，0)至(-150，150，0)，供电电流均为15 A.图6b、c、d分别为沿测线(y=0)在不同飞行高度10 m、30 m、50 m的感生电动势多测道图，由图可知感生电动势多测道图呈单峰分布，10 m飞行高度的感生电动势幅值较大，而50 m飞行高度的感生电动势幅值较小.同时，通过与图5计算结果比较，可以得知：当两个辐射源以相反方向的电流进行发射时，所激发的能量在空间具有一定的叠加作用，地下目标体的响应特征虽然彼此相似，但是两个电流方向相反且相互平行的辐射源电磁响应幅值明显大于单一发射源的情况，从而大大加强了探测信号强度.

2.2 模型二模拟结果

模型二：半空间电阻率为100 Ωm，异常体为两个直立板目标体，尺寸均为300 m×200 m×50 m，顶板埋深均为100 m，电阻率均为1 Ωm，两板间水平距离500 m，关于y轴对称放置(如图7a所示)；电源长均为500 m，电源距O点500 m，供电电流15 A.取飞行高度20 m，分别计算在单源、双源和四源激发情况下的电磁响应.

(1)两个地质体单辐射源垂直磁场分量响应特征

单辐射源激发模型如图7a所示，激发源1长500 m，由(-250，-500，0)至(250，-500，0)，供电电流15 A，飞行高度20 m.图7b所示的垂直磁场平面分布图，说明由单源激发时，由垂直分量响应不能区分这两个目标体.

(2)两个地质体、两个辐射源磁场垂直分量响应特征

两个目标体水平距离500 m，模型参数不变，如图8a所示，激发源1与图7模型一致，激发源2长500 m，位于(250，500，0)至(-250，500，0)，供电电流15 A.当利用双源同时激发时，取飞行高度30 m，在图8b所示的垂直磁场平面图中，出现了两个异常体的外边界轮廓，但是两个异常体之间的边界仍然显得比较模糊.

(3) 两个地质体、四个辐射源磁场垂直分量响应特征

两个目标体模型参数不变，如图9a所示，利用四个激发源同时进行激发(源长均为500 m，源1、源2与以上模型一致，源3位于(500，-250，0)至(500，250，0)，源4位于(-500，250，0)至(-500，250，0)，供电电流均为15A，飞行高度为30 m)，可以明显看到两个目标体分异明显.当利用四个辐射源激发时，由图9b所示的垂直磁场平面分布图可知，两个目标体的异常对称，异常幅值较高，其内边界的清晰度提高.通过与图8所示的计算结果相比，可以认为四个发射源的垂直磁场响应比两个发射源的响应效果更好，所反映的地质目标更真实.

3 结论

地空电磁法兼具探测深度大和探测效率高的优点，对于深部资源调查具有重要意义.在目前的实际应用中，地空电磁法一般单源激发，导致所接收到的信号较弱，难以获得高精度解释结果.本文基于自主研发的三维矢量有限元法正演模拟软件对多辐射源地空瞬变电磁响应进行三维模拟，当采用单源激发时，由于仅在一个方向激发，场的幅值和分辨率受到限制，获取的地质体的信息不全面，或者目标体并不能得到有效的识别.当采用多辐射场源作为地空电磁法的发射源时，能够获得不同角度的电磁场的辐射信息，从而获得比单源激发更高的分辨率.特别是四个源同时进行激发时，对于两个目标体有很好的识别能力.利用源的排列及电流方向等因素对信号影响的差异，合理布设电性源，可以达到加大勘探深度，提高多个目标体的分辨能力的目的.

图5 不同飞行高度感生电动势多测道图(a) 模型三维立体图； (b) 飞行高度10 m感生电动势多测道图； (c) 飞行高度30 m感生电动势多测道图； (d) 飞行高度50 m感生电动势多测道图.Fig.5 Multi-channel diagram of the dBz/dt with different flight heights(a) 3D view of model ； (b) Multi-channel diagram of dBz/dt with flight altitude 10 m； (c) Multi-channel diagram of dBz/dt with flight altitude 30 m; (d) Multi-channel diagram of dBz/dt with flight altitude 50 m.

图6 不同飞行高度感生电动势多测道图(a) 模型三维立体图； (b) 飞行高度10 m感生电动势多测道图； (c) 飞行高度30 m感生电动势多测道图； (d) 飞行高度50 m感生电动势多测道图.Fig.6 Multi-channel diagram of the dBz/dt with different flight heights(a) Model overhead view; (b) Multi-channel diagram of dBz/dt with flying height 10m; (c) Multi-channel diagram of dBz/dt with flying height 30m; (d) Multi-channel diagram of dBz/dt with flying height 50 m.

图7 双地质体、单一辐射源模型和垂直磁场分量平面图(a) 模型三维立体图； (b) 飞行高度20 m垂直磁场分量平面图.Fig.7 Plane diagram of the Bz of two geological bodies with a single radiation source(a) 3D view of model； (b) Plane diagram of Bz with flight altitude 20 m.

图8 双地质体、双辐射源模型和垂直磁场分量平面图(a) 模型三维立体图； (b) 飞行高度30 m垂直磁场分量平面图.Fig.8 Plane diagram of the Bz of two geological bodies with two radiation sources(a) 3D view of model； (b) Plane diagram of Bz with flight altitude 20 m.

图9 双地质体、四个辐射源模型和垂直磁场分量平面图(a) 模型三维立体图； (b) 飞行高度20 m垂直磁场分量平面图.Fig.9 Plane diagram of the Bz of two geological bodies with four radiation sources(a) 3D view of model； (b) Plane diagram of Bz with flight altitude 30 m.