1 引言

深海热液硫化物矿广泛分布于大洋扩张脊、弧后扩张中心以及热点火山构造带内，其经济价值在近几十年内被深刻认识.目前已知的典型硫化物矿热液区范围巨大，分布着储量可观的硫化物丘体和硫化物矿堆积体，矿藏丰富(邓希光, 2007; Богданов, 2007).当前，日益减少的陆地金属矿产资源以及迅速增长的金属资源需求推动着深海热液硫化物矿探测的发展，如何在热液区内探测深海热液硫化物矿床成为重要问题.

近年来，应用于深海矿产勘察的海洋瞬变电磁法受到广泛关注(Cheesman et al., 1990; Yu et al., 1997; Danielsen et al., 2003; Börner et al., 2015; Swidinsky et al., 2015).深海热液硫化物矿的瞬变电磁响应特征是应用该方法探矿的基础，已有学者对其做了研究.Cheesman等(1987)计算了海洋双半空间模型多种装置的瞬变电磁响应，刘长胜和林君(2006)通过海洋层状模型的频率域响应表达式推导了中心回线装置和重叠回线装置的瞬变电磁响应表达式，Swidinsky等(2012)在上述研究基础上计算了深海洋底导电层的中心回线装置及重叠回线装置瞬变电磁响应；周胜等(2012)运用全空间理论计算方法模拟了深海热液硫化物矿床层状模型重叠回线装置瞬变电磁响应；胡俊华等(2013)采用余弦变换多项式近似法对海洋层状模型的中心回线装置瞬变电磁响应进行了一维正演计算.Jang和Kim (2015)进行了中心回线瞬变电磁法探测深海热液矿床的一维正反演研究.以上模拟计算，把海底硫化物矿体简化为层状模型.然而，Evans和Everett (1994), Богданов, (2007), Swidinsky等(2012)分别对Galapagos洋脊热液矿化堆，大洋低速扩张脊热液成矿系统，俾斯麦海块状硫化物矿堆的研究表明，深海地球物理环境以洋底面为界分为海水与玄武岩基底两个半空间，硫化物在洋底面堆积成矿，洋底面以下为热液喷出通道及蚀变带，硫化物矿体整体呈锥形或透镜状，是全空间中的三维体，席振铢等(2016)通过分析研究大西洋洋中脊以及西南印度洋洋中脊热液硫化物矿实测数据将其简化为三维目标体.因而，采用三维正演计算是模拟深海热液硫化物矿瞬变电磁响应的有效手段.

为提高计算速度，增强计算实用性，本文将深海热液硫化物矿简化为低阻长方体进行三维数值模拟，应用全空间矢量有限元法计算其瞬变电磁响应.首先，通过双半空间模型的计算验证全空间矢量有限元法的正确性；其次，建立深海热液硫化物矿的三维电性模型，应用全空间矢量有限元法计算其瞬变电磁响应，总结响应特征；最后，通过物理实验验证深海热液硫化物矿瞬变电磁响应数值模拟结果的正确性.

2 计算方法

如图 1所示，深海硫化物矿体、海水和围岩的电磁边界复杂，Γ_in是海水与洋底的界面，Γ_a是热液硫化物矿的外边界，Γ_out是计算区域外边界.在这些边界上，电场的法向分量不连续，切向连续.传统的节点有限元计算方法需要所有场分量连续，因而不能获得精确解.不同于节点有限元法，矢量有限元法将场分量作为整体进行插值求得更精确的解(Jin, 2002; Sugeng, 1998; Cai et al., 2014; 黄威等, 2016).

2.1 边值问题

海洋瞬变电磁法的实施采用拖曳式工作方式，拖曳高度距海底仅50 m (周胜等，2012).由于海水的高电导率特性，需考虑海水层内电磁场的扩散.同时，由于深海热液硫化物矿多产出于水深达1.5~3.7 km的大洋扩张脊和弧后扩张中心等(Herzig and Hannington, 1995)，而根据趋肤深度公式δ≈503(Nabighian, 1987)，即便是0.1 Hz的低频信号，其趋肤深度约为918 m，可近似认为扩散至海面也已衰减殆尽.因此，本文海洋瞬变电磁的正演计算将海水层视为无限大半空间，在全空间内求解.

电磁场可通过解麦克斯韦方程组求得，取时间因子e^iωt，全空间的麦克斯韦方程为

(1)

(2)

其中, E是电场，H为磁场，ω为角频率，μ为磁导率，ε为介电常数，σ为电导率，J为源电流密度.对公式(1)取旋度，利用公式(2)消去磁场H，得电场的矢量波动方程：

(3)

上式中的电场是包含一次场和二次场的总场，总场的直接求解需考虑场源的空间奇异性，为模拟场源附近电场的快速变化，对模型进行剖分时须在场源附近做网格细化，导致待求量的增加.为解决上述问题，本文求解二次电场.将上式中的电场E分解：

(4)

式中，E^p表示一次电场，E^s表示二次电场，则电场矢量波动方程可分解为

(5)

(6)

式中，Δσ表示电导率异常，是异常体电导率与背景电导率的差值.对于一次电场E^p，计算其在两个半空间背景模型中的解析解.将本文要模拟的重叠回线装置发射线圈分成多段作电偶极源考虑，计算每段电偶极源TE极化模式的电场(Nabighian, 1987)，并沿线圈积分求得一次电场.

深海环境下，海水电导率为3 S·m^-1，玄武岩电导率取0.01 S·m^-1，热液金属硫化物矿电导率取10 S·m^-1，本文采用0.1~10⁵ Hz的频域响应进行时频转换，以上参数取值范围内，ω²με << ωμσ，因此忽略位移电流做准静态近似，待求解的二次电场矢量波动方程简化为扩散方程：

(7)

如图 1所示的计算区域，无穷远边界Γ_out上二次电场衰减为0，满足齐次狄利克雷边界条件：

(8)

内部介质分界面Γ_in及异常体界面Γ_a上，电场切向分量相等：

(9)

表示分界面单位法向向量，E₊^s、E_{_}^s表示界面不同侧的二次电场.

2.2 矢量有限元分析

本文采用的矩形单元如图 2所示，图中数字表示棱边局部编号.

根据矩形单元矢量基函数表达式(Jin, 2002)，单元内的二次电场可由矢量基函数表示为

(10)

E_i^s是棱边上的待求二次电场值，由于矢量基函数N_i满足

(11)

因此对于公式(10)表示的二次电场有

(12)

上式证明矢量基函数表示的电场自然满足零散度条件，从而保证求解结果二次可微，避免非物理解的出现.由公式(7)表示的二次电场扩散方程的残数为：

(13)

将单元内二次电场表达式(10)代入公式(13)并应用伽辽金方法(Jin，2002)可得单元内残数R^e，其矩阵形式为

(14)

其中[A^e]，[B^e]，[C^e]为12×12阶单元矩阵：

(15)

(16)

(17)

对所有单元组合单元方程(14)式，令与每一棱边相关的残数加权积分为0，得方程组：

(18)

其中，K为N×N阶整体矩阵，e、b为N阶向量，N是剖分模型总棱边数.

对于位于外边界上的棱边，其全局编码为n，令

(19)

(20)

实现狄利克雷条件的强加.由于矢量基函数N_i只在第i边上有切向分量，因此由式(10)表示的二次电场既保证了在棱边上的切向连续性，也保证了在单元界面上的切向连续性.

采用共轭梯度法求解方程组(18)得二次电场，二次磁场通过对式(10)应用法拉第定律求得：

(21)

以上求解出频率域内二次磁场H^s，对其进行Fourier反变换求得二次磁场时间域响应(考夫曼等, 1987).

3 模型计算

首先通过与Swidinsky双半空间模型计算结果进行对比验证程序正确性，然后将深海热液硫化物矿简化为双半空间中的低阻长方体，研究其瞬变电磁响应特征.本文计算均在CPU为2.5 GHz，RAM为4GB的个人计算机上完成.

3.1 双半空间模型

图 3为Swidinsky双半空间模型，海水电导率取3 S·m^-1，洋底半空间电导率分别取1、3 S·m^-1、10、30 S·m^-1和100 S·m^-1.

模型中没有地形起伏和复杂结构地质体，采用易于实现的矩形单元对模型进行剖分.剖分网格在海水-洋底界面处加密，剖分后各方向节点数21×21×14，总节点数6174，总棱边数17493，总单元数5200.洋底半空间电导率不同取值情况下计算时间分别为1459.89、970.81、665.69、279.88、211.25 s.图 4为双半空间模型的归一化磁场衰减曲线以及误差图，其中符号为本文计算结果，线条为Swidinsky计算结果，二者吻合良好，最大误差2.97%，平均误差0.31%，说明本文算法正确.

3.2 三维模型

将深海热液硫化物矿简化为位于洋底面上的低阻方形板状体，并将海水及洋底视为均匀半空间，建立如图 5所示的三维电性模型，海水电导率取3 S·m^-1，玄武岩电导率取0.01 S·m^-1，热液硫化物矿电导率取10 S·m^-1.图 5所示的深海热液硫化物矿三维电性模型异常体结构规则，采用矩形单元对模型进行剖分，剖分网格在异常体边界及洋底面处加密，图 6为网格剖分示意图，红色区域为异常体，东西走向及南北走向延伸200 m，厚度50 m，计算区域2000 m×2000 m×2000 m，剖分后各方向节点数47×47×46，总节点数101614，总棱边数298309，总单元数95220.发射线圈及接收线圈采用10 m×10 m的方形回线，激发电流10 A，上升沿0.2 ms，持续时间400 ms，关断时间为0.2 ms.测线沿东西方向，线距50 m，共9条测线，每条测线9个测点，点距50 m，主测线中心测点位于异常体中心，共计81个测点，总计算耗时263586.15 s，平均耗时3254.15 s.

图 7是主测线瞬变电磁响应多道剖面图，为直观地分析深海热液硫化物矿的瞬变电磁响应特征，取81个测点不同时间道的响应值绘制3D瞬变电磁响应图，如图 8所示.

由图 7及图 8可见，异常体以外区域，响应值低且平缓，变化不大；在异常体边界处，响应值由异常体外部至内部陡然增大，至一定幅值后趋向平缓，形成形似异常体的高响应值平台，响应异常形态极好地反映了异常体形态，二者近似一致.对比早期、中期、晚期的观测信号，早期响应幅值更高，且更清晰的反映异常体形态，信号幅值随时间衰减，且形态偏向平滑，此时响应异常形态由方形向圆形渐变.此变化主要由于随时间推移，电磁信号主要来自深部地层，出露于洋底的硫化物矿响应影响逐渐减小，响应异常形态趋向模糊.

4 物理实验

为验证数值模拟结果，根据物理模拟相似性准则设计并进行物理实验(Nabighian, 1987).实验仪器采用湖南五维地质科技公司的MTEM-08海洋瞬变电磁仪，发射线圈及接收线圈采用直径为10 cm的圆形线圈，激发电流波形为双极性矩形方波，发射频率6.25 Hz.实验采用120 cm×120 cm×85 cm水槽模拟全空间，水槽下方模拟高阻洋底半空间，水槽底部中心放置20 cm×20 cm×5 cm铜板模拟低阻深海热液硫化物矿.水槽内注入盐水模拟海水，观测线距5 cm，点距5 cm，拖曳高度5 cm，即发射及接收线圈沿铝板上顶面移动.图 9是取所有观测点不同时间道数据绘制的电流归一化磁场响应图，由图可见，物理实验采集的数据受到一定干扰，存在误差，但整体规律与数值模拟趋势一致，说明了数值模拟结果的正确性.

5 结论

(1)  采用全空间矢量有限元法，对硫化物矿体采用矩形单元模型剖分，应用Galerkin法推导有限元方程，对频率域响应进行Fourier反变换的计算方案，能够有效地计算出深海热液硫化物矿矿体的瞬变电磁响应.

(2)  通过简化的双半空间模型解析解和三维模型物理模拟分别与矢量有限元数值模拟进行比对，证明了矢量有限元算法计算全空间电磁环境是可行的.

(3)  对于海水、矿体以及围岩复杂电磁边界，应用全空间矢量有限元法模拟计算深海热液硫化物矿瞬变电磁响应异常与物理模拟响应趋势一致，而且计算方法简单精确，异常幅值明显，边界清晰.