海洋可控源电磁法 (Marine Controlled Source Electromagnetic, 简称MCSEM) 在海洋领域广泛地应用于洋中脊、大陆边缘地质构造研究以及海底油气资源和天然气水合物储层探测, 成效显著.

相比于传统的假设地下介质为各向同性的三维海洋可控源电磁反演技术，自然界中，地下介质的电导率常常表现为各向异性，以海底油气资源勘探为例，据估计，世界上大约30%的油气资源存在于岩性裂隙地层和泥砂岩薄互层中, 而这两种地层的宏观电阻率常常表现为各向异性.各向异性是在深入分析和解释实际电磁资料时必须加以考虑的因素，否则就会导致不容忽视的误差甚至错误的判断.目前，海洋可控源电磁油气勘探中, 薄储层电性各向异性特征常表现为垂直对称轴的横向各向同性 (Transverse Isotropy medium with Vertical symmetry axis, VTI).因此，开展VTI介质各向异性研究对认知地球介质结构和勘探开发复杂油气藏均具有重要的意义 (Constable，2010).目前，MCSEM各向异性正演算法中，广泛应用的是Li和Dai (2011)提出的MCSEM任意各向异性介质二维有限元正演模拟.随着三维各向异性MCSEM技术的逐步成熟，为解决低感应数问题，Liu和Yin (2013)研究了任意各向异性介质条件下的散度校正技术, 刘云鹤和殷长春 (2013)进行了三维航空电磁各向异性正演模拟.殷长春等 (2014)通过对MCSEM电场多分量响应及分布特征和各向同性情况的对比分析，总结了三维电各向异性MCSEM响应的影响规律和识别方法.Cai等 (2014)利用矢量有限元法进行VTI介质的三维MCSEM正演模拟.随着并行机节点计算能力的提高和内存可扩展性的进步，采用直接法求解离散方程逐渐得到重视，Li和Dai (2011)利用自适应有限元方法和直接法模拟了复杂地形MCSEM的二维响应.Da等 (2012)给出了基于多波前法的三维矢量有限元数值模拟结果.周建美等 (2014)通过利用直接法高效模拟了各向异性地层中三维MCSEM响应.Schwarzbach等 (2011)利用高阶有限元方法和直接法模拟了复杂地形MCSEM的三维响应.韩波等 (2015)采用并行直接法进行正演，模拟了海洋条件下多场源CSEM的响应.

基于非线性化反演算法的各向异性MCSEM反演中，Ramananjaona和MacGregor (2010)实现了VTI各向异性地层中MCSEM的2.5维反演.Sasaki (2013)采用高斯-牛顿算法实现了VTI各向异性地层中MCSEM的三维反演.Newman等 (2010)利用非线性共轭梯度反演算法实现了VTI各向异性地层中MCSEM的三维成像.Wiik等 (2011, 2013) 实现了层状各向异性地层中MCSEM的三维对比反演，并且还将层状各向异性地层中的三维对比源反演理论推广到MCSEM和海洋大地电磁的联合反演.Ray和Key (2012)实现了可控源电磁各向异性介质条件下的一维贝叶斯反演.

本文首先介绍了基于VTI介质的三维有限差分正演算法，针对MCSEM多源发射的问题，离散线性系统的求解采用并行直接法求解库 (MUMPS).其次，为提高三维反演速度，反演采用基于计算一阶导数和函数信息近似二阶导数的有限内存拟牛顿算法 (L-BFGS)(Haber, 2005).对于各向异性反演过程中，电阻率参数可能会出现负值、极大值或极小值等不符合一般物理常识的情况，引入基于电阻率约束的各向异性算法，提高解的可靠性.最后，利用VTI介质合成数据，分别进行电阻率各向异性覆盖层和电阻率各向异性高阻层的三维反演，并讨论:(1) 并行直接法的计算效率；(2) 各向异性覆盖层的约束反演效果；(3) 数据覆盖对电阻率各向异性高阻层反演结果的影响.

设电磁场时谐因子为e^－iωt, 非均匀各向异性介质中电流源的MCSEM三维模拟对应的二次场偏微分方程为 (Weiss and Constable, 2006)：

(1)

其中，

其中，E_p为一次场值，E_S为二次散射场值，μ₀为真空中的磁导率，σ₀为背景电导率，Δσ为异常电导率，σ_h为水平电导率，σ_v为垂向电导率，Δσ=σ－σ₀为模型电导率与参考模型电导率之差，σ₀通常假设为水平层状模型 (刘长胜等，2010).

由于电磁波在海水中快速衰减，若计算域的边界离异常区域足够远，我们可以认为二次场已充分衰减，则可以在网格边界上应用Dirichlet边界条件，将每个边界面上的网格单元的切向分量的值赋为0(Um and Alumbaugh, 2007)，即

(2)

三维MCSEM正演中，直接解法的优点是矩阵分解、求近似逆都与方程右端项无关，分解的结果可重复用于其他右端项，方程右端项仅与一次场以及场源相关.因此，这种解法的优点是对于有多个激发源的正演问题只需将矩阵做一次分解.

为验证三维数值解的正确性，利用三维程序计算近似二维模型E_x分量的幅值及相位响应，并与Scripps海洋研究所公开的基于自适应有限元2.5维海洋可控源电磁程序MARE2DEM进行比较.测试模型如下：空气电阻率为2×10⁸ Ωm；海水电阻率0.3 Ωm，深度1 km；覆盖层为VTI各向异性介质，X方向从-2.5 km到2.5 km，覆盖层的顶面深度为0.2 km，层厚度0.4 km，横向电阻率ρ_h=0.5 Ωm，垂向电阻率ρ_v=2 Ωm, Y方向从-100 km到100 km；高阻薄层为各向同性，X方向从-1.25 km到1.25 km，厚度为0.2 km，电阻率为50 Ωm，Y方向从-100 km到100 km；背景电阻率为1 Ωm.如图 1所示二维各向异性正演模型，X方向的电偶极子源放置于 (0 km, 0 km, 0.95 km) 处，观测测线 (Y=0 km) 为0.5 km间距 (从-4 km到4 km).网格剖分为92×56×34，在观测区域和源所在区域，网格是均匀的，间距为0.1 km.覆盖层的顶面深度为0.2 km.收发距之间采用先大后小的策略，以因子3等比增长添加边界网格.

图 1

Fig. 1

图 1 二维各向异性正演模型 (a) 横向电阻率模型; (b) 垂向电阻率模型. Fig. 1 2D anisotropic earth for CSEM modeling (a) Horizontal resistivity model; (b) Vertical resistivity model.

利用此网格剖分计算频率为0.25 Hz的幅值及相位响应，结果如图 2所示.通过与数值解的对比，幅值和相位的平均相对误差小于1%.机器配置为:CPU为Intel Core I7-4770(4核心8线程，主频3.4 GB).内存为16 GB.操作系统为Ubuntu14.04.计算频率为0.25 Hz时一次正演所消耗的时间为669 s.

图 2

Fig. 2

图 2 三维可控源电磁正演结果验证 (a) 幅值相对误差; (b) 相位相对误差. Fig. 2 Modeling results for 3D CSEM (a) Amplitude relative error; (b) Phase relative error.

为获得较真实的反演模型，除了提高反演技术手段外，尽可能施加先验信息作用于反演过程是缩小“解空间”、减少解非唯一性和提高解可靠性的有效办法.借鉴 (Commer and Newman, 2009) 在MCSEM各向同性介质反演中采用的对数变换方式，将基于电阻率不等式约束的反演问题变换为无约束最优化问题.按照Tikhonov正则化理论, MCSEM反演问题可归结为求解以下目标函数的极小值问题:

(3)

ψ(m) 为总的目标函数，ψ_d(m)=‖V(d-F(m))‖₂²为数据目标函数项，ψ_mh(m_h)=‖W(m_h)-(m_h)₀‖₂², ψ_mv(m_v)=‖W(m_v)-(m_v)₀‖₂²分别为模型水平和垂直方向约束目标函数项，ψ_m(m) 选用拉普拉斯算子.λ_h和λ_v为正则化因子，使得数据目标函数项和模型约束目标函数项达到一种平衡.d为观测数据向量，(m_h)₀、(m_v)₀为初始模型矢量，F为模型的正演响应向量，V为数据方差矩阵，W为模型光滑矩阵.l、u分别为电导率的下、上界.目标函数确立后，反演问题就是采用L-BFGS算法极小化目标函数的问题.在L-BFGS (赵宁等，2016) 算法中，目标函数及其梯度信息是反演的关键信息.对式 (3) 的目标函数求梯度得：

(4)

根据式 (4) 中的定义，并令A为正演响应向量的雅可比矩阵，则式 (4) 可进一步表示为

(5)

通过互易定理计算式 (5) 中雅可比矩阵的转置与一个向量的乘积，即A(m)^TV^－1(d－F(m))，从而构建目标函数的梯度信息.

反演过程中数据拟合误差 (RMS) 的定义：N_d为反演数据个数, F为模型的正演响应向量, W_d为数据方差矩阵, d为反演数据.

(6)

模型设置如下：空气电阻率2×10⁸ Ωm，海水电阻率0.3 Ωm，海水深度1 km；覆盖层为VTI各向异性的棱柱体，距离海底0.2 km，模拟区域为5×2×0.4 km³，横向电阻率ρ_h=0.5 Ωm，垂向电阻率ρ_v=2 Ωm；高阻薄层为各向同性的棱柱体，模拟区域为2.5×2×0.2 km³，电阻率为50 Ωm；背景电阻率为1 Ωm，如图 3所示.测区共85个测点, 范围为[-4 km, -2 km, 1 km]×[4 km, 2 km, 1 km]，测线间距1 km，测点间距0.5 km.网格剖分同2.2节的网格剖分形式，电偶极子源沿X方向放置于测区，发射频率为1 Hz.反演采用海底为均匀半空间1 Ωm的初始模型，观测同线方向 (inline), λ_h和λ_v均为1.合成数据加入5%的随机噪声.反演内存消耗为11.27 GB.

图 3

Fig. 3

图 3 各向异性覆盖层模型 (a) 横向电阻率模型; (b) 垂向电阻率模型. Fig. 3 The anisotropic model of covering layer (a) Horizontal resistivity model; (b) Vertical resistivity model.

图 4为VTI介质模型的合成数据分别进行各向同性反演 (表 1，a项) 和各向异性反演 (表 1，b项) 的RMS对比图.从图中可知，各向异性反演具有更小的拟合差.图 5为发射频率为1Hz时，不同发射源个数的反演耗时.三维正演中发射源每改变一次激发的位置 (或方位)，所形成的线性系统方程只改变右端项，而不改变系数矩阵.从图中可知，85个发射源反演耗时比17个发射源反演耗时增加了约3倍，基于并行直接法的MCSEM非常适合于海洋电磁所特有的多场源问题，场源的增加只会带来一部分的额外计算量.

图 4

Fig. 4

图 4 各向异性与各向同性收敛曲线对比 Fig. 4 The RMS convergence rate for isotropic and anisotropic inversion

表 1(Table 1)

表 1 电阻率各向异性算例 Table 1 Inversion cases used in anisotropic model 序号 合成数据 反演类型 a VTI各向异性 各向同性反演 b VTI各向异性 各向异性 (无约束) c VTI各向异性 各向异性 (电阻率约束)

表 1 电阻率各向异性算例 Table 1 Inversion cases used in anisotropic model

图 5

Fig. 5

图 5 不同发射源个数的反演耗时 Fig. 5 The consuming time of different source numbe

模型设置如4.1数值算例.图 6a(表 1，a项)，对三维VTI介质各向异性模型的响应数据进行各向同性反演，从XZ、YZ反演结果剖面分析，其中，电阻率各向异性低阻覆盖层和高阻层的位置、电阻率值与实际情况有很大的偏差，会造成电阻率各向异性覆盖层下高阻层的误判.图 6b、6c(表 1b、c项)，结果表明，无约束各向异性反演对于各向异性覆盖层和下面的高阻薄层都有很好的反映，但对各向异性覆盖层和高阻薄层边界圈定不足.基于不等式约束的各向异性反演不仅对异常体电阻率值有很好的反映，并可以圈定各向异性覆盖层及其高阻薄层的位置和边界.对于电阻率各向异性覆盖层模型，通过表 1的对比，基于不等式约束的各向异性反演能够获得很好的垂向电阻率信息，可以更好的识别高阻目标层.

图 6

Fig. 6

图 6 电阻率各向异性覆盖层反演结果 (a) 各向异性的合成数据，各向同性反演; (bl) 无约束反演结果 (横向电阻率); (b2) 无约束反演结果 (垂向电阻率)；（cl) 电阻率约束反演结果 (横向电阻率); (c2) 电阻率约束反演结果 (垂向电阻率). Fig. 6 Inversion results of anisotropic model for covering layer (a) Inline+broadside, isotropic inversion; (bl) Inline, unrestraint inversion (Horizontal resistivity); (b2) Inline, unrestraint inversion (Vertical resistivity); (cl) Inline+broadside, restraint inversion (Horizontal resistivity); (c2) Inline+broadside, restraint inversion (Vertical resistivity).

模型设置如下：空气电阻率为2×10⁸ Ωm；海水电阻率0.3 Ωm，深度1 km；高阻薄层为VTI各向异性的棱柱体，距离海底0.6 km，2.5×2×0.2 km³，横向电阻率ρ_h=30 Ωm，垂向电阻率ρ_v=150 Ωm；背景电阻率为1 Ωm，如图 7所示.多测线共85个测点，范围为[-4 km, -2 km, 1 km]×[4 km, 2 km, 1 km]，测线间距1 km，测点间距0.5 km.测线从 (-4 km, 0 km, 1 km) 到 (4 km, 0 km, 1 km)，共17个测点.电偶极子源沿X方向放置于测区，发射频率为1 Hz.网格剖分同2.2节的网格剖分形式，反演采用的合成数据为理论数据加入5%的高斯噪音，反演采用海底为均匀半空间1 Ωm的初始模型，λ_h和λ_v均为1.发射频率为1 Hz.

图 7

Fig. 7

图 7 高阻层模型 (a) 横向电阻率模型; (b) 垂向电阻率模型. Fig. 7 The high resistance reservoir model (a) Horizontal resistivity model; (b) Vertical resistivity model.

图 8a(表 2，a项)，三维VTI介质各向异性模型的响应数据进行各向同性介质反演，从XZ、YZ反演结果剖面分析，高阻层的位置和电阻率值与实际情况有一定的偏差.图 8b、8c(表 2b、c项)，结果表明，Inline数据对各向异性高阻层有一定的反映，但对异常体的位置反映不准确.Inline和broadside数据不仅对异常体电阻率值有很好的反映，并圈定了高阻薄层的位置.对于4.2节电阻率各向异性高阻层模型，通过表二的对比，Inline和broadside数据的VTI各向异性反演效果更佳.

图 8

Fig. 8

图 8 电阻率各向异性高阻层反演结果 (a) 各向异性的合成数据，各向同性反演; (b1) 横向电阻率反演结果 (Inline); (b2) 垂向电阻率反演结果 (Inline); (c1) 横向电阻率反演结果 (Inline和broadside); (c2) 垂向电阻率反演结果 (Inline和broadside). Fig. 8 Inversion results of anisotropic model for high resistance layer (a) Inline+broadside, isotropic inversion; (b1) Inline, unrestraint inversion (Horizontal resistivity); (b2) Inline, unrestraint inversion (Vertical resistivity); (c1) Inline+broadside, restraint inversion (Horizontal resistivity); (c2) Inline+broadside, restraint inversion (Vertical resistivity).

表 2(Table 2)

表 2 电阻率各向异性高阻层算例 Table 2 Inversion cases used in anisotropic high resistance reservoir model 序号 合成数据 反演类型 a Inline和broadside 各向同性反演 b Inline VTI各向异性约束反演 c Inline和broadside VTI各向异性约束反演

表 2 电阻率各向异性高阻层算例 Table 2 Inversion cases used in anisotropic high resistance reservoir model

本文实现了基于纵向各向异性 (VTI) 介质的频率域海洋可控源电磁三维反演.针对海洋可控源电磁多源发射的问题，正演调用大规模并行矩阵直接求解 (MUMPS) 求解离散线性系统，提高了多源问题的计算效率.三维反演采用基于电阻率不等式约束的有限内存BFGS (L-BFGS) 算法，缩小了“解空间”，提高了解的可靠性.利用VTI介质各向异性合成数据，分别进行了电阻率各向异性覆盖层和电阻率各向异性高阻层的三维反演，针对前者，高阻层的反演结果受各向异性覆盖层的影响较大，采用Inline方向的观测系统，各向异性反演能够获得很好的垂向电阻率信息，可以更好的识别高阻目标层.针对后者，Inline和broadside数据不仅对异常体电阻率值有很好的反映，并圈定了高阻薄层的位置.

本文的研究由中石油东方地球物理公司 (BGP) 资助，感谢BGP国家“千人计划”专家余刚，综合物化探处何展翔总工、孙卫斌副总工、刘雪军副总工的帮助与支持.感谢中国海洋大学李予国教授对反演模型设置的交流和犹他大学Alexander Gribenko教授对三维正演计算速度提出的建议.