海洋可控源电磁法勘探包括海底高阻油气资源勘探^[1-5]、深海热液硫化物矿勘探^[6]、海底工程检测^[7]等, 能够有效进行油气层识别和海上油气储层定量评价^[1-4], 2000年以来得到了大力发展^[8-9]。针对不同水深和油气藏规模, 海洋电磁法发展了多种不同的电磁勘探装置^[10-16], 如频率域深拖拽发射-接收装置、频率域浅拖拽发射-接收装置、时间域深拖拽发射-接收装置、时间域浅拖拽发射-接收装置、时间域海底垂直发射-垂直接收装置等。国内外学者对不同海洋电磁法勘探装置的高阻目标体探测能力进行了很多对比研究, 如SHANTSEV等^[17]、GUO等^[18]比较了频率域深拖拽和浅拖拽装置在不同水深情况下对海底高阻目标体的探测能力, 建议在深水域采用深拖拽装置, 浅水域采用浅拖拽装置。GOLDMAN等^[19]比较了频率域深拖拽水平发射-固定阵列接收装置、时间域海底垂直发射-垂直接收装置、时间域回线源装置对海底高阻目标体的探测能力, 建议在深水域采用频率域方法, 浅水域以及横向小目标体探测采用时间域方法。MACGREGOR等^[10]对各类装置探测海底高阻目标体的优缺点进行了较为系统的总结。

深海热液硫化物矿^[20]广泛分布于深海的大洋扩张脊和火山构造带等区域, 具有显著的经济开采价值。日益增长的资源需求推动着海底矿产勘探的发展, 海洋可控源电磁法在深海热液硫化物矿勘探中的应用得到越来越多的关注^{[6, 21-23]}。深海热液硫化物矿是一种典型的低阻目标体, 相应的深海低阻目标体电磁法探测主要采用时间域海底中心回线或重叠回线装置^{[6, 21-22]}。而目前广泛应用的海洋电磁法各类装置^[11-16]主要应用于海底高阻油气探测, 关于深海低阻目标体探测能力的讨论和分析较少。本文研究分析了目前广泛应用的海洋电磁法装置对深海低阻目标体(以热液硫化物矿为例)的探测能力。考虑深海区域的目标体探测采用浅拖拽装置接收信号时需要经过两次厚的海水层吸收, 会导致电磁信号幅值过小^[17], 同时考虑到海底操作的效率, 本文主要针对工作效率较高的频率域深拖拽发射-接收装置^[12]和时间域深拖拽发射-接收装置^[14]进行讨论。通过计算典型的深海三维低阻目标体模型的电场响应, 分析这两类海洋电磁法勘探装置对深水域低阻目标体的探测能力。

1 正演理论

三维频率域正演电磁响应采用耦合势有限体积法计算^[24]。频率域电磁法满足的控制方程为：

$ \nabla \times \nabla \times E + {\rm{i}}\omega {\mu _0}\sigma \mathit{\boldsymbol{E}} = - {\rm{i}}\omega {\mu _0}\mathit{\boldsymbol{J}} $

(1)

式中:E为频率域电场强度矢量, J为外加电性源, ω为圆频率, μ₀为真空磁导率, σ为地层电导率。引入电场的矢势A和标势φ, 将电场分解为无散场和无旋场之和：

$ \mathit{\boldsymbol{E}} = \mathit{\boldsymbol{A}} + \nabla \varphi $

(2)

控制方程转换为关于矢势与标势的混合Helmholtz方程：

$ {\nabla ^2}{A_x} + {\rm{i}}\omega {\mu _0}\left[ {\sigma \left( {{A_x} + \frac{{\partial \varphi }}{{\partial x}}} \right)} \right] = - {\rm{i}}\omega {\mu _0}{J_x} $

(3a)

$ {\nabla ^2}{A_y} + {\rm{i}}\omega {\mu _0}\left[ {\sigma \left( {{A_y} + \frac{{\partial \varphi }}{{\partial y}}} \right)} \right] = - {\rm{i}}\omega {\mu _0}{J_y} $

(3b)

$ {\nabla ^2}{A_z} + {\rm{i}}\omega {\mu _0}\left[ {\sigma \left( {{A_z} + \frac{{\partial \varphi }}{{\partial z}}} \right)} \right] = - {\rm{i}}\omega {\mu _0}{J_z} $

(3c)

$ \nabla \cdot \sigma \left( {\mathit{\boldsymbol{A}} + \nabla \varphi } \right) = - \nabla \cdot \mathit{\boldsymbol{J}} $

(3d)

式中:A_x, A_y和A_z为矢势A在直角坐标系中的3个分量; J_x, J_y和J_z为电性源J在直角坐标系中的3个分量。选择足够大的计算区域Ω, 则在区域外边界∂Ω上的电磁场将非常小, 因此可以采用简单的截断边界条件：

$ \begin{array}{l} \mathit{\boldsymbol{n}} \times \mathit{\boldsymbol{A}}\left| {_{\partial \mathit{\Omega }}} \right. = 0\\ \mathit{\boldsymbol{n}} \cdot \mathit{\boldsymbol{A}}\left| {_{\partial \mathit{\Omega }}} \right. = 0\\ \;\;\;\;\varphi \left| {_{\partial \mathit{\Omega }}} \right. = 0 \end{array} $

(4)

其中, n为边界单位矢量。为了求解控制方程(3)在边界条件(4)下的数值解, 采用Yee氏交错网格对求解区域进行空间离散, 利用高斯定理对旋度和散度算子进行离散处理, 发射源采用直接离散方法, 最终得到关于矢势A和标势φ的离散控制方程：

$ \mathit{\boldsymbol{FX}} = \mathit{\boldsymbol{Q}} $

(5)

式中:系数矩阵F为大型非对称稀疏复矩阵, X为离散网格上的未知数, Q为离散的源项。采用直接法求解器PARDISO ^[25]求解该离散方程, 得到三维频率域正演电磁响应。

三维时间域正演电磁响应采用基于隐式时间步迭代的拟态有限体积法计算。忽略位移电流, 回线源瞬变电磁法对应的时间域Maxwell方程为：

$ \frac{{\partial \mathit{\boldsymbol{b}}}}{{\partial t}} = - \nabla \times \mathit{\boldsymbol{e}} $

(6a)

$ \nabla \times \mu _0^{ - 1}\mathit{\boldsymbol{b}} - \sigma \mathit{\boldsymbol{e}} = \mathit{\boldsymbol{s}} $

(6b)

式中:e是时间域电场强度矢量, b是时间域磁感应强度矢量, t是时间, s是时间域外加源项。采用简单的自然边界条件

$ \begin{array}{l} \mathit{\boldsymbol{b}} \times \mathit{\boldsymbol{n}} = 0\\ \mathit{\boldsymbol{n}} \times \mathit{\boldsymbol{e}} = 0 \end{array} $

(7)

则回线源瞬变电磁法在t=0时刻空间中只存在稳定的磁场分布, 即初始条件为:

$ \mathit{\boldsymbol{b}}\left( {t = 0} \right) = {\mathit{\boldsymbol{b}}^0} $

(8)

式中, b⁰为t=0时刻空间中的磁场分布。将(8)式采用弱形式表示, 并采用Yee氏交错网格进行空间离散, 利用积分形式的斯托克斯定理处理电场旋度的离散, 得到控制方程空间离散的矩阵表示。采用关断源, 并采用无条件稳定的欧拉后向差分格式进行时间步离散, 得到最终的离散控制方程为：

$ \left( {\mathit{\boldsymbol{C}}_{{\rm{URL}}}^{\rm{T}}{\mathit{\boldsymbol{M}}_\mu }{\mathit{\boldsymbol{C}}_{{\rm{URL}}}} + \Delta {t^{ - 1}}{\mathit{\boldsymbol{M}}_\sigma }} \right){\mathit{\boldsymbol{e}}^n} = \Delta {t^{ - 1}}\mathit{\boldsymbol{C}}_{{\rm{URL}}}^{\rm{T}}{\mathit{\boldsymbol{M}}_\mu }{\mathit{\boldsymbol{b}}^{n - 1}} $

(9a)

$ {\mathit{\boldsymbol{b}}^n} = {\mathit{\boldsymbol{b}}^{n - 1}} - \Delta t{\mathit{\boldsymbol{C}}_{{\rm{URL}}}}{\mathit{\boldsymbol{e}}^n} $

(9b)

式中:C_URL为旋度算子的离散形式, M_μ和M_σ分别为磁导率μ₀和电导率σ离散形成的矩阵, Δt为迭代时间步长, eⁿ和bⁿ分别为第n次时间步迭代得到的电场和磁场。先通过解析求得初始场b⁰, 然后通过求解时间步迭代的线性方程组(9), 即可得到不同时刻的电磁场响应。本文采用直接法求解器PARDISO^[25]求解线性方程组(9), 同时, 为了保证计算精度和效率, 选取分段等间隔的时间步长。

2 探测能力分析

频率域和时间域深拖拽发射-接收装置探测深海三维低阻目标体模型设置参考文献[6]中的一维热液硫化物矿模型。当海水深度大于1000m时认为是深海域^[8], 可以忽略空气层, 因此可以设置如图 1所示的正演模型。该模型上半空间为海水层, 电导率为3S/m; 海底地层电导率为1S/m, 海底下方20m存在一个电导率为20S/m的三维矩形低阻目标体, 其x方向长度为500m, y方向长度为500m, z方向厚度为30m。采用海底深拖拽水平发射-水平接收装置, 发射源与接收机位于同一高度, 操作方便, 距离海底上方5m, 发射源水平方向距离海底下方三维目标体边界100m(图 1)。

2.1 频率域响应分析

在探测埋深较大的高阻油气藏时, 一般采用较低的发射频率(0.1~10.0Hz)。由于本文研究的低阻目标体埋深较浅, 因此采用相对高的发射频率(100Hz)^[11]。深海域海洋电磁法在100Hz时的噪声水平为2×10^－15V/Am^2[11]。发射源位置水平方向距离海底下方三维目标体边界100m, 如图 1所示, 随船拖拽的多个接收机与发射机位于同一高度。图 2为不同偏移距时频率域深拖拽水平发射-接收装置的响应振幅和归一化振幅。图 2a中黑色圆圈为采用开源软件Dipole 1D^[26]计算的不含异常体的一维模型的解析解, 红色曲线为采用三维耦合势有限体积算法计算的不含异常体的一维模型的数值解, 蓝色曲线为采用三维耦合势有限体积算法计算的三维模型(如图 1所示)的数值解, 黑色虚线为噪声水平; 图 2b为存在3D低阻体的电场响应相对不含低阻体的半空间模型的归一化振幅。图 2a中黑色圆圈与红色曲线重合, 说明本文算法的计算结果是有效的。图 2a中红色曲线和蓝色曲线对比可知, 该装置对于海底低阻目标体存在显著异常, 最小归一化振幅小于0.02, 最大异常响应位于偏移距500m处(图 2b), 其主要的异常响应均大于噪声水平, 说明该装置能够有效地探测海底低阻异常。

为了更清晰地理解频率域深拖拽水平发射-接收装置对海底低阻体的探测原理, 绘制了含有低阻目标体和不含低阻目标体两种情况三维模型的电流和电场分布图。图 3为不含低阻异常体的半空间模型的电场幅值(即$\sqrt {E_x^2 + E_y^2 + E_z^2} $)及电流分布图(等值线表示电场分布, 箭头表示电流方向, 下同), 图 4为不含低阻异常体的一维模型的x方向电场幅值|E_x|分布图; 图 5为含有低阻异常体的三维模型的电场幅值及电流分布图, 图 6为三维模型的x方向电场幅值|E_x|分布图, 白色矩形为目标体所在区域。由图 3和图 5可以清晰地看出, 电性源的空间电流分布既存在垂直分量, 也存在水平分量, 因此该装置既对高阻目标体灵敏(典型的海底高阻油气藏勘探), 也对海底低阻目标体灵敏。图 5说明, 由于三维目标体电阻率低, 导致其所在区域的能量衰减更快, 三维低阻体周围区域的电场值减小, 电场的x分量也相应减小(图 6)。对比图 4和图 6可知, 低阻体的存在导致电流分布发生畸变, 从而在偏移距500m区域出现一个极小的归一化电场E_x, 这与图 2中电场响应随偏移距的变化相对应。

2.2 时间域响应分析

采用时间域深拖拽发射-接收装置探测深海三维低阻热液硫化物矿目标体类似于长偏移距瞬变电磁法^[27-29], 存在最优化偏移距选取问题。本文通过对比不同偏移距情况下一维低阻模型的瞬变电磁响应, 得到近似最优化的偏移距。图 7为一维低阻层状模型, 分别计算偏移距为50, 100, 300, 500, 800m的模型地层响应, 如图 8所示。其中实线为不含低阻层的半空间地层响应, 虚线为含有低阻层的地层响应。由图 8可见, 当偏移距为300m时, 低阻层的异常响应最大。当偏移距逐渐大于300m或者小于300m时, 低阻层的异常响应都表现为逐渐减小, 可知图 7所示的一维低阻层状模型最优化偏移距在300m左右。因此, 本文在随后计算图 1所示三维模型的瞬变电磁响应时, 偏移距选取300m。

采用三维有限体积正演算法计算图 1所示三维模型的正演响应。图 9为时间域深拖拽水平发射-接收装置的电场响应振幅和归一化振幅, 其中图 9a中黑色圆圈为采用开源软件Dipole1D结合正余弦数字滤波算法^[30]计算得到的不含低阻异常体的一维模型的解析解, 红色曲线为采用三维有限体积算法计算的不含低阻异常体的一维模型的数值解, 蓝色曲线为采用三维有限体积算法计算的三维低阻模型的数值解, 黑色虚线为噪声水平^[31]; 图 9b为存在三维低阻体的电场响应相对不含低阻异常体的一维模型的归一化振幅。图 9a中黑色圆圈与红色曲线重合, 说明本文算法的计算结果是有效的。图 9a中红色曲线和蓝色曲线对比可知, 时间域深拖拽水平发射-接收装置对于海底低阻目标体存在明显异常, 最小归一化振幅为0.55, 最大异常响应位于5×10^-3s时刻。其主要异常响应均大于噪声水平, 说明该装置能够有效探测海底低阻异常。

为了清晰地理解时间域深拖拽水平发射-接收装置对海底低阻体的探测原理, 绘制了含有低阻目标体和不含低阻目标体的三维模型电流和电场分布图。图 10为不含低阻异常体的模型在10^-5s时刻x方向电场E_x幅值分布图, 图 11为含有低阻异常体的模型在10^-5s时刻x方向电场E_x幅值分布图; 图 12为不含低阻异常体的模型在10^-2s时刻x方向电场E_x幅值分布图, 图 13为含有低阻异常体的模型在10^-2s时刻x方向电场E_x幅值分布图。由图 10和图 11可知, 在10^-5s时刻, 三维低阻异常体的存在虽然能够导致其所在区域场值减小, 从而导致空间中的电场分布产生畸变, 但由于能量主要集中在发射源附近, 因此三维低阻异常体的存在对于300m接收点处的电场分布影响甚微, 表现为图 9中10^-5s时刻存在三维低阻异常体和不存在三维低阻异常体两种情况的电场幅值|E_x|相同。随着时间增加, 能量逐渐向外扩散, 三维低阻异常体对接收点信号的影响逐渐增大, 表现为图 9中归一化振幅随时间增加逐渐减小, 到5×10^-3s时刻三维低阻异常体对接收点信号的影响达到最大, 此后随着时间增加, 三维低阻异常体的影响逐渐减小, 对应接收点处的归一化振幅逐渐增大。对比图 12和图 13可知, 在10^-2s时刻, 三维低阻异常体的存在导致电场幅值|E_x|的分布产生较大畸变, 三维低阻异常体对能量的强吸收作用导致接收点处的电场幅值低于不含三维低阻异常体时的电场幅值, 对应图 9中10^-2s时刻电场的归一化振幅小于1.0。

3 结束语

本文研究了频率域和时间域深拖拽发射-接收装置对于深海三维低阻目标体的探测能力。通过计算典型的深水域三维低阻目标体模型的电场响应、归一化电场响应和空间中电场分布特点, 分析了两种不同海洋电磁方法的应用效果, 取得以下成果：

对于本文所计算的三维模型, 频率域响应的最大异常位于偏移距500m处, 时间域响应的最大异常位于偏移距300m处。频率域和时间域深拖拽发射-接收装置均能够在小偏移距观测到深海三维低阻目标体的电场异常, 频率域响应的最小归一化振幅小于0.02, 时间域响应的最小归一化振幅为0.55, 频率域电磁法的异常幅值要明显大于时间域电磁法。

本文研究内容能够为海洋电磁法勘探, 特别是小目标体和深海热液硫化物矿勘探提供一定的参考。