复杂地形结构下瞬变电磁三维正演模拟−以大西洋洋中脊TAG热液区为例

聂佐夫，陶春辉,2,3*，沈金松，朱忠民

(1.上海交通大学 海洋学院，上海 200240；2.自然资源部第二海洋研究所，浙江 杭州 310012；3.自然资源部海底科学重点实验室，浙江 杭州 310012；4.中国石油大学 (北京) 地球物理学院，北京 102249)

1 引言

自1972年，多国政府和大学科学家合作开展美国国家海洋和大气管理局（National Oceanic and Atmospheric Administration,NOAA） TAG（Trans-Atlantic Geotraverse）项目以来，位于26°N大西洋洋中脊的TAG热液区一直是研究热液活动的焦点[1−2]。热液系统循环过程中形成的多金属硫化物，是一种极具经济价值的矿产资源，近年已成为海底资源研究的一大热点[3]。目前，已发现的海底硫化物矿体的规模大小差异较大，由几十米到几百米不等，因此需要一种行之有效的探测方法来确定海底多金属硫化物堆积体的位置、深度、规模，进而有效地评价其资源量及其经济价值，为后续开采工作打下基础。根据前人在有热液活动的TAG丘体[4]的钻探结果以及在已停止热液活动的区域[5]的钻探取芯结果分析可知，硫化物矿体围岩为导电性差的新鲜或蚀变玄武岩，且其导电性与岩石孔隙中的海水盐度相关。而硫化物矿体中存在较高含量的金属元素，它们与围岩存在显著的导电性差异，这为海底电磁方法探测热液多金属硫化提供了物质基础。

经过半个世纪的发展，海洋电磁在理论和装备上都有了长足的进步。20世纪80、90年代，在海洋瞬变电磁研究中比较有代表性的是Edwards研究团队，1986年Edwards和Chare[6]计算了上覆导电海水的均匀半空间电偶极子瞬变电磁二次场响应，这一项研究成为了海洋瞬变电磁法勘探的理论基础。次年，Cheesman等[7]、Everect等[8]在此基础上总结出水平电偶极子和磁偶极子两种设备十分适合用来开展瞬变电磁近底勘探的结论，随后进一步开展了海试工作，论证了其在海底电导率勘探上的可行性。最近几年，海洋电磁法探测海底热液多金属硫化物堆积体已成为研究的焦点。2012年，Swidinsky等[9]结合实际海底矿床勘探建立模型，推导同点装置一维层状介质正演公式，结果分析表明，在一定半径线圈及电流条件下，同点探测系统是确定导电矿床厚度及其上覆盖层厚度的有效工具。2015年，德国GEOMAR设计研制了一套瞬变电磁重叠回线装置MARTEMIS，并在地中海海域已知硫化物堆积体进行了初步测试，在已知的矿体附近，探测到了瞬变电磁异常的存在[10]。2010年，经自然资源部第二海洋研究所组织、湖南五维地质科技有限公司和北京先驱高技术开发公司共同研发出了海洋瞬变电磁重叠回线装置，并于南大西洋13.2°S热液区开展试验，探测到明显异常[11]。这一装置在2012年的中国大洋26航次于大西洋洋中脊TAG热液区进行了探测，得到热液区内TAG丘体的瞬变电磁响应，并利用视电阻率反演结果，将海底多金属硫化物矿体近似看作两层空间中的T型异常体[12−13]。最近几年在正演方面，国内对于有限元方法的研究也取得了一定进展。2016年，李瑞雪等[14]采用全空间矢量有限元法，对硫化物矿体进行矩形单元剖分，从频率域出发，计算出深海热液硫化物矿体的瞬变电磁响应，且与物理模拟结果具有良好的一致性。2017年，Zhang等[15]直接在时间域求解，实现了基于自适应有限元法的航空瞬变电磁三维正演模拟，并分析了地形因素对航空瞬变电磁响应的影响。2019年，殷长春等[16]实现了时间域海洋电磁三维正演模拟，改善了网格质量，提高了计算精度，并分析了在地形隆起与下凹时，海底有无块状导体时的近底电磁响应。2020年，彭荣华等[17]首次考虑海洋可控源电磁探测中的各向异性因素，利用三维正演模拟方法研究海底硫化物堆积体响应，发现各向异性因素会引起海洋可控源电磁结果失真从而不利于准确的数据解释。

TAG热液区拥有丰富的地质和矿体产状的背景资料，是目前研究程度最高的活动热液区之一，在TAG丘体开展瞬变电磁三维正演研究，分析空间电磁场的分布特征和接收装置的响应特征，有利于我们了解仪器状态、海底地形以及矿体分布对瞬变电磁二次场响应的影响。结合该区域开展的海洋瞬变电磁探测资料，可以检验瞬变电磁方法在探测海底多金属硫化物堆积体方面的可行性与有效性。

同时，由于目前海洋瞬变电磁实际数据处理仍然以一维反演拼接成拟二维剖面为主[18]，但利用一维层状介质模型的瞬变电磁响应来解释复杂地形情况下的响应会存在一定假象[19]，因此有必要采用高维的反演方法，并充分考虑地形起伏、探测装置测线轨迹以及装置的姿态等因素，从而得到更为合理的电导率剖面。因此，本研究可为未来海底多金属硫化物堆积体电导率剖面的精确反演工作打下基础。

2 TAG热液区地质背景

TAG热液区位于大西洋洋中脊26°08′N附近，热液区面积约为 5 km×5 km，水深范围为 3 430～3 670 m（图1）。由1处存在热液活动的TAG丘体和多处已停止热液活动的多金属硫化物丘体组成[2]。其中，TAG丘体存在高温热液喷口（温度最高达363℃），已形成直径为200 m，高为50 m的环形热液产物堆积体，活动黑烟囱体高约12 m[20]。大洋钻探计划对TAG丘体进行了一系列钻孔采样[4]，在丘体周围的5个区域（TAG-1至TAG-5）共计钻孔17个，最深钻孔达125 m。通过对钻孔中取回的岩芯进行一系列地球化学分析，在海山堆积体内发现了高品位的Cu、Zn、Ag和Au多金属硫化物[21]，并得到 TAG 丘体的内部结构[4,21−22]：丘体表层为硫化物角砾，下伏富含金属元素的黄铁矿、二氧化硅与硬石膏角砾，最深处为金属含量稍低的脉状矿化硅化玄武岩。对TAG热液区内多处硫化物丘体采集的硫化物进行物性测量工作，样品包括硫化物和作为参照的玄武岩，测得的部分结果如表1所示[23−24]。可见，即使硫化物孔隙度低于玄武岩时，也具有更高的电导率，说明硫化物并不符合基于电解质导电的阿尔奇定律[25]。硫化物因其导电形式更类似于半导体从而具有高电导率[26]。

表1 TAG热液区样品及部分陆地玄武岩样品的电导率测量结果[23−24]Table 1 Conductivity measurement results of sulfide samples from the TAG hydrothermal field and various onshore basalt samples[23−24]

图1 大西洋洋中脊TAG热液区高分辨率地形图（a）以及TAG丘体内部地质结构（b）Fig.1 High-resolution bathymetry map of the TAG hydrothermal field at the Mid-Atlantic Ridge (a) and the geological structure of the TAG mound (b)

3 瞬变电磁响应的数值模拟方法

对于瞬变电磁场计算，我们采用了库伦规范准则下磁矢量势−电标量势与自适应有限元相结合的方法[27]，在此仅给出算法相关的基本关系。在各向同性海水和海底地层中，瞬变电磁场满足如下Maxwell方程：

式中，E为电场强度（单位：V/m）；H为磁场强度（单位：A/m）；J为传导电流密度（单位：A/m2）；Je为外加电流密度（单位：A/m2）；µ和 σ 分别为介电常数和电导率。

由公式（1）和公式（2）可得到电场强度E的双旋度方程，公式为

利用亥姆霍兹分解理论，矢量场E可变为

式中，A为矢量势；ϕ为标量势，将公式（4）代入公式（3），得到下式：

同时，将公式（4）代入连续性方程∇×σE=−∇×Je可得公式：

为求解方程（5）和方程（6），考虑各向同性介质的Dirichlet边界条件，对于求解域Ω和闭合边界Γ，可得：

式中，n为单位法向量。本文采用Galerkin加权余值法导出方程（5）、方程（6）的弱解形式，并用 Delaunay非结构化四面体单元剖分，四面体单元内各边与节点的矢量位与标量位可表示为[16]

图2 四面体单元自由度（修改自参考文献 [17]）Fig.2 Degrees of freedom associated with a tetrahedral element (modified from reference [17])

建立有限元方程如下求解。得到 A ,ϕ：

式中，系数矩阵元素和源项的计算关系可参考文献 [17,27−28]，而对于时间步的离散使用后向欧拉法[29−30]。相应地，磁场和电场的计算关系为

4 TAG丘体电性参数的三维建模

本文利用高精度多波束测深数据构建TAG丘体的三维地形图[31]，为得到海底以下的电导率分布信息，参考大洋钻探及地球化学探测结果，建立如图3所示的三维地电模型。模型核心部分为400 m×400 m×300 m的三维结构，其中TAG丘体的尺寸约为180 m×180 m×150 m。模型核心部分外加厚 1 km 的边界层以确保电磁波衰减，并在层边界添加Dirichlet边界条件。根据钻孔数据确定的TAG丘体内部的电导率分布如图3c和图3d所示，丘体内部高品位的黄铁矿电导率为60 S/m，丘体表层的硫化物角砾电导率为30 S/m，丘体深部为电导率稍低的脉状矿化硅化玄武岩，电导率约为20 S/m。正演模拟中海水电导率设为3 S/m（海水温度为2～3℃），围岩玄武岩电导率设为0.1 S/m。

在瞬变电磁正演模拟中，发射源采用圆形重叠回线装置，半径为10 m。发射线圈中施加20 A的阶跃电流，为观测丘体上方的瞬变电磁响应，接收关断后感应电动势 ∂ B/∂t 随时间的衰减，接收间距为10 m，测线如图3a中红线所示。

图3 TAG 丘体所在区域三维电导率分布模型Fig.3 3D conductivity model for the TAG mound

模型建立后，采用非结构化网格（自由四面体）对三维模型进行离散化。为了保证正演模拟的精度，应对发射源和地形变化大、电导率产生差异的界面附近使用局部更为精细的网格。为提高运算效率，根据源的位置对网格进行局部剖分，利用局部网格替代网格数量庞大的全局网格，以减少系统矩阵的自由度。采用PARDISO直接求解器来求解方程组，该求解器是目前最快的线性稀疏矩阵求解方法之一[32]。

5 电磁场响应的特征分析

5.1 电磁场时空分布特征

为了解电磁场在该地质模型中的传播规律，确定有效信号的观测范围，笔者正演计算了瞬变响应在电流关断后的时空分布。在TAG丘体模型西南角离底20 m处建立重叠回线装置，线圈半径为10 m，电流为20 A，通过正演得到TEM脉冲响应∂ Bz/∂t在海底地形上的扩散过程如图4所示。当发射电流关断后，变化的磁场会在硫化物丘体内部产生感应涡流，它会随观测时间不断向下和向外扩散且幅值衰减形成“烟圈效应”[33–34]，由于硫化物矿体与围岩的电导率差异巨大，大多数涡流集中于TAG丘体区域。

图4 发射线圈关闭后4个时间步长上TAG丘体地形处∂ Bz/∂t扩散图Fig.4 Forward model of the diffusion of ∂Bz/∂t across the bathymetry of the TAG mound area at four timesteps after transmitter turn-off

为验证瞬变电磁法对高导硫化物的探测能力，本文同时也对纯玄武岩基底的海底进行了模拟，即把原模型中具有导电性的硫化物角砾、黄铁矿和脉状矿化硅化玄武岩都替换为高阻玄武岩，将其正演响应视为背景场，并将含硫化物的异常场与背景场做归一化，对∂ Bz/∂t归一化后取对数后的结果如图5所示，该结果可以指示瞬变电磁响应对高导体的可探测性。在 1×10−5s、3.162 3×10−4s和 1×10−3s时刻，线圈上方的海水因电导率一致而没有产生异常，但下方TAG丘体处的高导体导致了明显的 ∂ Bz/∂t差异。同时由于线圈位于丘体西南侧（剖面左侧），涡流系统经过高导体后，导致整个剖面右侧的高阻围岩响应也出现了明显的偏差，而两次正演模拟在剖面右侧的差异很小。

图5 发射线圈关断后，由背景场归一化后的异常场，其中第二列图为西南−东北方向的竖直剖面Fig.5 Forward model of the normalized ∂Bz/∂t for a homogeneous seafloor model vs.the model which includes the conductive target after transmitter turned off,graphs at the second column are the SW-NE sections of the TAG mound area

5.2 实测瞬变电磁响应的拟合与影响因素分析

5.2.1 TAG 热液区瞬变电磁数据采集

海上瞬变电磁测量采用重叠回线装置[12]，线圈为1.95 m×0.75 m 的矩形多匝回线，线圈匝数 40 匝，等效面积为 58.8 m2，发送电流为 20 A 的双极性方波，占空比为50%，方波频率为0.625 Hz。海上作业时，科考船以 1～2 kn（0.5～1 m/s）的速度拖曳瞬变电磁拖体匀速前进，通过调节光缆的长度使拖体保持在距海底上方 30～50 m。

2012年的大洋26航次第2航段对海上瞬变电磁装置进行了测试，结合TAG热液区矿体位置对测线进行设计，利用深海瞬变电磁双拖体装置采集海底多金属硫化物矿体的响应数据，共完成了6条测线[13]。其中EM02测线由西南向东北从TAG丘体正上方经过，并经过热液停止活动的MIR区，测线总长为2.67 km。该测线数据在这两处已知区域均发现明显的瞬变电磁异常，具有良好的可靠性，因此本文以EM02测线为例来说明实际数据采集中存在的影响因素。

测线位置及测量得到的瞬变响应如图6所示，图6b和图6c分别对应TAG丘体和MIR区上方的瞬变电磁剖面，该数据已做中值滤波处理，并利用源归一化来剔除线圈大小与发射电流对响应的影响。在此基础上，为消除仪器本身对瞬变电磁响应的影响从而突出异常的相对变化，利用实测与模拟的海水背景响应对实测数据进行归一化处理。可以看出，该测线所测得的电压幅值存在一定的震荡，在已知的硫化物矿体上方的瞬变电磁响应有明显的隆起，且TAG丘体处的响应呈“双峰”状。为探究其原因，有必要对数据采集过程中可能影响实测数据的各个因素进行研究。

图6 TAG丘体与MIR区周边多波束地形图及大洋26航次瞬变电磁EM02测线（a）、TAG丘体瞬变电磁归一化响应平剖图（b）和MIR区瞬变电磁归一化响应平剖图（c）Fig.6 The multibeam bathymetry of the area contains the TAG mound and the MIR zone with the EM02 survey line from COMRA 26th Cruise (a),normalized section of TEM response at the TAG mound (b),and normalized section of TEM response at the MIR zone (c)

5.2.2 瞬变电磁响应拟合与分析

本节模拟了重叠回线装置的近底拖曳测量方式，离底高为20 m时得到的瞬变电磁响应以及高导异常点对应位置如图7所示。高导异常响应与无明显异常的背景响应在5 ms前没有明显区别，因为早期响应大多贡献自仪器附近的海水。在5 ms后，二者开始产生差异，深部高导体的存在导致了更强的二次场，从而使异常响应衰减变慢。将存在异常的响应所在的位置投点至平面图可以发现，高导异常出现的位置均在TAG丘体上方，没有异常的背景响应则采集自远离丘体的测点上。晚期20 ms处异常最为明显的响应来自于TAG丘体的山顶处，此处矿体厚度最大。

图7 瞬变电磁重叠回线系统接收线圈响应（a）和高导异常所在测点位置（b）Fig.7 Transient electromagnetics response for coincident loop system (a) and the locations that consist high conductivity anomalies (b)

我们将正演模拟结果与TAG热液区实测瞬变电磁数据进行对比。经过预处理后的部分实测数据如图 8 所示，采集有效时窗范围为 3～30 ms，30 ms后的数据因信噪比过低需被剔除。黑色曲线为已知TAG丘体的瞬变电磁响应，浅红色与灰色曲线分别为离底高度20 m时异常最大与最小的正演模拟结果。可知TAG丘体的三维正演结果与该区域的实测数据能够很好地拟合，且实测数据表明，已有效探测到了高电导率异常。

图8 TAG 热液区部分实测瞬变电磁数据Fig.8 Part of the collected TEM data at the TAG hydrothermal field

为研究仪器离底高度变化对海底多金属硫化物探测效果的影响，把离底高度分别设置为40 m、60 m和80 m来对比其二次场响应特征。为方便比较，绘制电压平剖图如图9所示。综合TAG丘体地形可以看出，由于地形起伏及仪器拖曳方式改变了线圈和高导矿体之间的耦合关系，使早期响应的异常形态发生畸变。而对于晚期响应，当矿体厚度逐渐变大时，响应的幅值也变大。对比 20 m，40 m、60 m 与 80 m的正演模拟结果，可知仪器离TAG丘体越远，异常幅值越小，且早期畸变响应也变小。理论上而言，即使是离底高度为80 m时，晚期的瞬变电磁响应也能探测到丘体的高导异常，但实际上，实测数据是存在噪声干扰的。以中国大洋26航次在TAG热液区采集到的瞬变电磁数据为例，实测数据的有效时窗范围为 1～30 ms，30 ms后的数据由于幅值很小，受到的噪声干扰也就随之变大，往往被剔除。图中的方框虚线即代表30 ms时的背景场平剖曲线，即有效信号的范围在该曲线以上，此时可以发现，离底高度80 m的平剖曲线已几乎与背景场重合，无法有效反映出TAG丘体的高导异常。

图9 不同离底高度下的电压平剖图（a）和对应测线轨迹（b）Fig.9 Section of TEM response under different survey altitudes (a) and the corresponding survey lines (b)

在实际近底勘探应用中，瞬变电磁重叠回线仪器会因拖曳过程及洋流扰动而晃动。为研究仪器姿态对瞬变电磁响应的影响，在线圈模型建立时引入倾角，以X轴方向为例，当线圈位于TAG丘体正上方20 m及80 m时，分别计算线圈倾角为45°和60°时的接收线圈响应，结果如图10所示。可以看出由于激发的一次场随着线圈发生倾斜，导致二次场产生变化，同时接收线圈的倾斜也使得实测数据产生干扰。对于离底高度为20 m的高导矿体响应，这一干扰主要发生在信号的早期（小于 3 ms）及晚期（大于 70 ms）；而对于离底高度为80 m的响应，姿态变化产生的干扰减弱，且集中在30 ms之后。只有当仪器离底高度远大于仪器探测深度时，所测得的响应可等效为将仪器置于无限空间海水中的响应，此时线圈的姿态变化则不会对所测响应产生影响。因此，线圈的离底高度和海底地质体电导率分布都会使其在姿态不稳时的响应产生干扰。

图10 线圈在不同离底高度下，沿X轴倾角变化时的接收线圈响应曲线Fig.10 Transient electromagnetic method response under different loop angle along the X-axis and different altitude

在上文中，我们发现正演模拟的早期响应存在明显的畸变，这一现象在实测数据中也有体现。为更好地了解这一现象产生的原因，有必要进一步探究地形起伏以及测线布置对瞬变电磁响应的影响。设置3组不同的正演模型：起伏地形+水平测线、水平地形+水平测线以及水平地形+起伏测线。其中，水平地形指的是将原TAG丘体区域地形转换为水平地形，地下高导体形态与之前模型一致；起伏地形为原TAG丘体区域；水平测线的离底高度为20 m，起伏测线为原模型测线。分别进行正演模拟，得到不同的正演模拟结果（图11）。

图11 不同地形及测线条件下的瞬变电磁响应平剖图Fig.11 Section of TEM response under different topography and survey line conditions

起伏测线条件下，当探测过程不受热液区崎岖地形影响时，可以发现其与起伏地形下的平剖图相比，早期响应更为平缓；此时的畸变信号来自仪器离底高度的变化，因为其改变了重叠回线装置以下的电导率分布。对起伏地形进行水平测线观测，可以发现早期响应不再发生畸变，其幅值变化可以更好地反映海底良导体的厚度，但固定深度的观测方式会使丘体两侧水深更深的矿体响应发生衰减，需要结合晚期响应才可圈定异常范围。

6 结论

本文采用有限元法实现了海洋瞬变电磁响应的三维正演模拟。通过多波束、钻探和地质资料建立了复杂地形与结构的TAG丘体三维电导率模型，并计算得到了电磁场的时空分布，正演响应结果与TAG热液区的实测数据能够很好地拟合。研究了仪器不同离底高度与姿态条件下的瞬变电磁响应规律，以及地形因素和测线设置对瞬变电磁响应的影响。总结得到以下认识：（1）电流关断后的磁场在富含多金属硫化物的TAG丘体内产生涡流，其接收线圈响应特征与高阻体响应特征具有明显差异，因此瞬变电磁重叠回线装置可以有效探测海底多金属硫化物堆积体；（2）瞬变电磁的有效探测深度与所测数据的时窗有关，以TAG热液区实测瞬变电磁数据的信噪比为例，应保证瞬变电磁重叠回线装置在拖曳过程中离底高度不高于60 m，才能有效探测到高导体异常；（3）复杂的热液区地形以及仪器拖曳过程中高度的变化会严重影响瞬变电磁早期响应，因此在后续数据处理过程中应结合精确的海底地形数据以及超短基线定位数据，以防止对实测数据的错误解释；或者在保证离底高度不大于探测深度的前提下，采用水平测线的拖曳方式以降低早期响应的扰动，简化后续处理反演流程；（4）拖曳过程中线圈姿态变化会干扰瞬变电磁的响应，所以有必要记录探测装置的姿态，来更好地分析所探测异常的来源。本文验证了瞬变电磁法在海底硫化物勘探中良好的应用前景，同时考虑了复杂地形、测线轨迹、仪器姿态因素对数据的扰动。在实际反演过程中，若已知并引入这些参数，将能够为复杂地形与结构条件下的热液区硫化物堆积体的精确三维反演提供支撑。