0 引言

海洋可控源电磁法(Marine Controlled Source Electromagnetic，MCSEM)近年来迅速发展成为一种与地震勘探互补的海洋油气勘探技术.经过国内外学者的共同努力，相应的MCSEM正演技术得到了较大发展，包括一维正演(Constable and Srnka, 2007；刘长胜等，2010；刘颖和李予国, 2015)、二维正演(Unsworth et al., 1993；Li and Key, 2007；沈金松和孙文博，2009；Key and Ovall, 2011；刘颖和李予国，2011)和三维正演(汤井田等, 2007, 2014；陈桂波等，2009；Schwarzbach et al., 2011；杨波等，2012；Puzyrev et al., 2013;殷长春等，2014；蔡红柱等，2015；Farquharson and Farquharson, 2015).尽管目前为止，已有这些相对成功的研究工作，然而，传统的海洋可控源电磁正演大多数都没有考虑岩石成分、结构对于电磁场的影响，或者只是从频率相关性出发(Manning and Athavale, 1986;Sternberg, 1991；Ulrich and Slater, 2004；何展翔和王绪本，2007；Veeken et al., 2009；Clennell et al., 2010；郭宁宁等，2011)，研究岩石导电特性.例如，通过引入各种复电阻率模型可以研究相应的复电阻率模型参数对于海洋可控源电磁场的影响，即激发极化效应.但是这种方式无法直接解释岩矿石中的矿物成分(如储层砂岩中的黄铁矿)是如何影响岩矿石在低频条件下的电导率特性.在海洋可控源电磁法勘探中，含有黄铁矿化的储水、油气砂层是常见的探测对象，因此，研究储水、储油气砂层中的黄铁矿或者其他矿物成分对海洋可控源电磁场的影响规律对于提高相应的数据解释的精度具有显著的理论意义和实践价值.而一个能够将沉积岩的导电特性与其组成成分合理连接起来的岩石物理模型是提高可控源电磁法数据解释精度的关键因素之一.目前常用的等效介质岩石物理模型方程主要有Archie's方程(Archie, 1942)、Hanai-Bruggeman(HB)方程(Bruggeman, 1935; Hanai, 1960a, b ; Bussian, 1983)和Asami方程(Asami, 2002).但是这些方程仅仅针对两相介质有效，对于更为接近实际情况的多相介质而言，这些方程都无法适用.基于此，Han等人基于Asami方程结合微增方式，对Asami方程进行了进一步扩展而提出了适用于多相介质的微增(Incremental)岩石物理模型(Han et al., 2015), 该模型有效地将岩、矿石的一些基本结构及岩石学特征在宏观尺度进行表示，其参数与岩、矿石的物理属性相关联，为研究多相复合岩、矿石的导电特性提供了一种定量的分析方法.

因此，本文在前人工作的基础上，引入多相微增模型实现了海洋可控源电磁三维非结构化矢量有限元数值模拟研究，并全面研究各种参数(黄铁矿含量, 孔隙度、黄铁矿电导率、粒长宽比和含水饱和度)对含有黄铁矿化储层的导电性以及海洋可控源电磁场各分量的影响规律，以期能够进一步充实海洋可控源电磁法勘探理论基础.

1 Incremental模型的引入 1.1 Asami方程

Asami方程是Asami 2007年根据Maxwelle-Wagner理论以及经过进行大量标本测试后提出的两相等效介质模型方程，公式为

(1)

其中，σ为等效介质的电导率，σ_a是背景介质的电导率，σ_g为夹杂质的电导率，ϕ为夹杂质的体积百分数，ϕ=1-ϕ为等效介质的孔隙度大小.L_x、L_y、L_z表示颗粒x、y、z方向的极化因子，表达式为

对于椭球体夹杂质而言：当长短轴的比例α < 1时，有：

(2)

1.2 多相Incremental模型

Asami等效电导率方程在颗粒纵横比(a=1、Lz=1/3)情况下，可以简化成著名的Hanai-Bruggeman(HB)等效电导率方程(Bruggeman, 1935；Hanai, 1960a, b ；Bussian, 1983).若进一步忽略夹杂质的电导率σ_g，那么Asami等效电导率方程可进一步简化成Archie′s方程(Archie，1942).然而，这些等效介质模型都是针对两相介质而已，不能适用于多相介质(如，富含黄铁矿的含油气的储水沙层).所以Asami等人提出了一种适合更加符合真实情况的多相介质的微增模型(Asami，2002；Han et al., 2015)正好填补了此部分空白.

不失一般性，以四相微增模型为例，图 1为四相微增模型介质形成过程中的第一次微增过程，首先假定其中三种矿物的初始浓度分别为C₂、C₃、C₄，将它们分成N块，即每小块的浓度分别为C₂/N、C₃/N、C₄/N.接着，以水作为初始背景介质，将每种成分介质(C₂/N，C₃/N，C₄/N)根据一定的顺序依次添加，形成一次等效介质，这就完成了一次微增过程，然后将该次形成的等效介质作为下一次微增过程的初始介质，继续依次添加各种介质成分，直到循环添加完所有介质成分，即完成N次.而在每添加一种介质成分时，就运用一次Asami方程来计算介质的等效电导率，最后共进行3N次Asami方程计算.其中关键参数是各种物质成分的初始浓度值(C₂，C₃，C₄)，因为只有知道了每种微增的物质成分的初始浓度时，才能计算每次微增过程中的‘背景’介质的含量，从而利用公式(1)计算出相应的等效电导率值.例如，在计算出C₂情况下，在每次微增phase2物质时，方程(1)右端项ϕ=1-C₂/N(为‘背景’介质的含量)，接着将有关参数代入方程(1)即可计算出在添加phase2物质后的等效电导率.类似的，微增phase3物质和phase4物质亦如此，不同的是，方程(1)的右端项分别变为1-C₃/N和1-C₄/N.为了计算各种介质的初始浓度值(C₂，C₃，C₄)，首先给定相应的体积百分数(V₂, V₃, V₄)以及微增总数N, 通过附录所给出的方程(A1—A11)，来进行逆推计算求得各种介质的初始浓度值，具体可参照附录部分.以此类推，最终形成多相微增模型.本文为了简便起见，只研究海洋电磁法勘探中水，黄铁矿、石英、以及油气这四相介质组合形成的海洋沉积岩的成分、结构对电导率特性以及海洋可控源电磁场的影响特征.

图 2为在不同孔隙度下，分别改变黄铁矿的电导率、含水饱和度、黄铁矿的体积百分数、黄铁矿的颗粒纵横比的等效介质的电导率曲线.图 2a的岩石参数即σ_w=3.3、σ_qz=0.001、σ_gas=0.0001、纵横比a_qz=0.2、a_gas=0.2、S_w=0.5、V_py=0.02、σ_py=1.5:1.5:15；图 2b的具体参数为σ_py=15、S_w=0.05:0.5:1，其余参数同图 2a；图 2c的具体参数为V_py=0.002:0.02:0.1、σ_py=15，其余参数与图 2a一致；图 2d的具体参数a_py=0.02:0.2:1、σ_py=15，其余同图 2a.

从图 2a可知，在同一孔隙度下，改变黄铁矿的电导率，等效电导率变化不大; 随着孔隙度的增大，等效电导率增大；从图 2b可知，在同一含水饱和度下，孔隙度越大，等效电导率越大; 在同一孔隙度下，含水饱和度越大，等效电导率越大；从图 2c可知，在同一孔隙度下, 黄铁矿的体积百分数越大，等效电导率越大，近似成线性关系; 在同一黄铁矿的体积百分数下, 孔隙度越大，等效电导率值越大.从图 2d可知，同一孔隙度下，黄铁矿颗粒纵横比越大，等效电导率值变化不大; 在同一纵横比下，孔隙度越大，等效电导率越大.

2 正演理论 2.1 控制方程

在地球物理电磁勘探领域，通常采用低频的电磁信号，忽略位移电流，采用正弦谐变时间因子e^iωt，Maxwell方程可写为

(3)

(4)

其中，ω为角频率，μ₀为真空中的磁导率.J_s为源中电流分布，σ为电导率.

为了避免源点的奇异性，采用二次场算法，即将总场分解为一次场和二次场之和，电导率分解为背景电导率和异常电导率之和，公式为

(5)

(6)

根据(3)—(6)式，我们可以进一步推导出异常场所满足的微分方程为

(7)

通过方程(7)可发现，源的部分有关背景电场.可以通过计算全空间或者半空间或者层状介质的解析解为背景场.方程(7)可通过有限差分法、有限单元法、积分方程法、边界单元等方法求解，在求解出二次场电场后，通过以下方程，便可求出二次场磁场，进而求出总场，公式为

(8)

2.2 边界条件

当边界距离目标区域足够远时，可将二次场近似为0，因此，为了简便，本文采用狄利克雷边界条件：

(9)

式中，n为边界面上的法向单位矢量.

2.3 矢量有限元分析

文中采用非结构四面体单元网格进行空间离散，选用Whitney型矢量基函数(Jin，2002)，公式为

(10)

其中，L_i1^e为单元e的第i条边的第1个节点的权函数，L_i2^e单元e的第i条边的第2个节点的权函数，l_i^e为单元e第i条边的长度.

三维四面体单元的边与节点的关系如图 3所示.假定第i条边上的电场为E_i^e，那么四面体单元内的电场则可表示为

(11)

运用伽里金法，方程(7)的弱解形式可表示为

(12)

根据第一矢量格林函数，我们可以得到：

(13)

考虑到区域内边界是自动满足和外边界是人为规定的情况下，方程(13)中的面积分可以忽略，即：

(14)

将方程(14)代入方程(7)中得：

(15)

因为：

(16)

(17)

方程(15)可以重新写为

(18)

单元e的加权残余矩阵方程可写为

(19)

其中：

(20)

(21)

我们可以重新定义：

(22)

(23)

通过总体合成所有单元的有限元方程，便可得到总体有限元方程，公式为

(24)

其中：

(25)

对于线性方程组(24)的求解，我们利用预处理的IDR(s)迭代算法进行求解，系数矩阵采用CSR格式进行压缩存储.

2.4 线性方程组的求解

线性方程组的求解的速率是影响正演效率的关键因素之一，通常而言，主要有两种方式求解方程组.一种是直接解法，另一种是迭代求解方法.与迭代法相比，直接法求解大型方程组具有精度高、稳定性好的优点，但其内存需求大、求解效率低，在一般的微机上实现较为困难.因此，对于大型、稀疏的线性方程组的求解，大多是采用Krylov子空间迭代法.值得注意的是，采用迭代法对系数矩阵条件数较差的方程进行求解时，预条件因子能够明显提高收敛速率.本文也借鉴此思路，引入Induced Dimension Reduction (IDR(s))迭代算法结合不完全LU分解的预条件因子对线性方程组进行求解，并将其与常用的BICGSTAB、QMR等迭代算法的求解效率进行对比.具体计算流程如表 1所示.

3 算例 3.1 算法验证

为了验证文中算法的正确性，设计一个四层海洋沉积模型，其中一层为含有黄铁矿、油气、水、砂岩这四相介质的油气层.几何图形如图 4所示，模型参数如表 2所示.采用沿着x方向的水平电偶极子作为激发源，坐标为(0，0，980).激发频率为0.1 Hz，发射电流为1 A.接收站位于500:200:10000.图 4b为文中计算的数值解与采用汉克尔积分算法求解的准解析解对比曲线, 图 6为数值解与准解析解的相对误差曲线.从图 5可知，数值解与准解析解比较吻合，从图 6可以看到，相对误差最大为1.2%，对比结果验证了本文算法的正确性和有效性.图 7为IDR(s)迭代算法与常用的BICGSTAB、QMR迭代收敛曲线对比，从图 7上可知，IDR(s)算法比BICGSTAB和QMR算法稳定且收敛更快.

3.2 算例2

为了分析不同岩石结构参数对海洋可控源电磁场影响规律，设计一个含有黄铁矿、水、油气、石英四相介质的油藏模型，分别改变油气藏的岩石结构参数进行正演计算.油藏的规模大小为5 km×5 km×1 km，中心位于(2000，0，2500).表 3列出了海洋油藏模型的Incremental模型参数, 几何示意图如图 8a所示.采用COMSOL软件直接生成的非结构化网格剖分如图 8b所示，在激发源和测线附近同样进行了局部加密处理.激发源为x方向的水平电偶极子，位置在(-2000，0，980)处，发射电流为1 A，因篇幅有限，设定发射频率为0.1 Hz.接收点电极位于z=1000 m处.分析y=0，z=1000处的电磁场二次场各分量响应特征.

(1) 孔隙度对可控源电磁场的影响规律

为了分析不同的岩石孔隙度对海洋可控源电磁场的影响规律，设置不同的油气藏的岩石孔隙度进行正演计算，具体参数如表 3所示.

图 9为取不同孔隙度下的电磁场二次场虚、实分量曲线，从图 9上可知，随着孔隙度越大，含水饱和度一定的情况下，电磁场二次场各虚、实分量异常变化幅度明显减小.

(2) 黄铁矿电导率对可控源电磁场的影响规律

为了分析不同的颗粒电导率对海洋可控源电磁场的影响规律，设置不同的油气藏的颗粒电导率进行正演计算，黄铁矿电导率σ_py分别取0.15、1.5、15，孔隙度Φ=0.4.其余参数同表 3所示.

图 10为取不同黄铁矿下的电磁场二次场虚、实分量曲线，从图 10上可知，随着黄铁矿电导率的增大，电磁场二次场虚、实分量异常幅度先是明显减小，然后异常形状发生倒转.因此，可知，储水、储油沙层中的低阻黄铁矿化对于海洋可控源电磁场影响比较大.

(3) 黄铁矿体积百分数对于可控源电磁场的影响规律

为了分析不同的黄铁矿的体积百分数对海洋可控源电磁场的影响规律，设置不同的油气藏中的黄铁矿体积百分数进行正演计算，黄铁矿体积百分数V_py分别取0.02、0.2、0.4.其余参数同表 3所示.

图 11为取不同黄铁矿的体积百分数下的电磁场二次场虚、实分量曲线，从图 11上可以看到，随着黄铁矿的体积百分数越大，电磁场二次场虚、实分量的变化趋势与在一定的体积百分数下，不断增大黄铁矿的电导率的情况下是相似的，即异常幅度先减弱，然后接着异常形状发生倒转.

(4) 含水饱和度可控源电磁场的影响规律

为了分析不同的含水饱和度对海洋可控源电磁场的影响规律，设置不同的油气藏中的含水饱和度进行正演计算，含水饱和度S_w分别取0.05、0.4、0.8，孔隙度Φ=0.5.其余参数同表 3所示.

图 12为取不同含水饱和度的电磁场二次场虚、实分量曲线，从图 12上可知，在孔隙度一定的情况下，随着含水饱和度的增大，电磁场二次场虚、实分量异常幅度均减小，并且变化幅度比较大.

(5) 颗粒纵横比对可控源电磁场的影响规律

为了分析不同的颗粒纵横比对海洋可控源电磁场的影响规律，设置不同的油气藏中的黄铁矿的颗粒纵横比进行正演计算，黄铁矿的颗粒纵横比a_py分别取0.03、0.3、0.9，Φ=0.5，S_w=0.5，其余参数同表 3所示.

图 13为取不同黄铁矿颗粒纵横比时的电磁场二次场虚、实分量曲线，从图 13上可知，随着颗粒纵横比的增大，电磁场二次场虚、实分量的异常幅度变大，但是异常变化幅度与孔隙度，含水饱和度、黄铁矿电导率和体积百分数而言，颗粒纵横比对电磁场的影响相比小一些.

4 结论

本文通过引入Incremental模型计算地下多相介质的等效电导率，结合非结构化矢量有限元法实现了频率域海洋可控源电磁法三维数值模拟研究，并分析了地下介质的孔隙度、含水饱和度等参数对于可控源电磁场二次场的影响特征，可得结论：

(1) 引入多相微增模型表征多相岩、矿石的导电特性，是可行的，有效的，可为在宏观尺度上进一步了解岩、矿石的成分、结构对其导电特性的影响提供一种新的途径.

(2) 岩石孔隙度在一定的情况下，含水饱和度、黄铁矿电导率、体积百分数、颗粒纵横比越大，电磁场二次场变化异常幅度变化比较大.在含水饱和度一定的情况下，岩石孔隙度越大，电磁场二次场异常幅度减弱.因此在进行数据解释时，若加以考虑这方面的因素，有助于进一步提高数据解释的精度.

(3) 为了更深入了解孔隙度、含水饱和度等参数对MCSEM探测的影响，还需要在反演过程中考虑这些参数.因此，基于多相Incremental模型的海洋可控源电磁法三维反演将是本文后续重要工作之一.

致谢

感谢两名匿名审稿人对本文提出宝贵的修改意见，才使得文章得以进一步完善.

附录  多相介质在微增过程中的体积百分数计算公式

如图 1所示，在第一次微增过程中(step1)，加入浓度为C₂/N的第一种物质成分时，各相介质的体积百分数计算公式为

(A1)

其中，V_W为水的体积百分数，V₂、V₃、V₄分别为初始浓度是C₂、C₃、C₄的三种物质成分.

继续加入浓度为C₃/N的第二种物质成分时(step2)，各相介质的体积百分数计算公式为

(A2)

继续加入浓度为C₄/N的第三种物质成分时(step3)，各相介质的体积百分数计算公式为

(A3)

在完成第一次微增后的等效介质作为下一次微增过程的初始背景介质，在第i次微增过程中(2 < i < N)，各相介质的体积百分数计算公式为

(A4)

Step2：

(A5)

Step3：

(A6)

值得注意的时，在每一次微增过程中，各种组成成分体积百分数都满足关系为

(A7)

其中i=1:N, j=1, 3.

根据公式(A1)—(A7), 在给定微增总数N的以及水等四种物质成分的体积百分数(V_W, V₂, V₃, V₄)前提下，最终可获得方程组为

(A8)

(A9)

(A10)

(A11)

通过求解方程(A8)—(A11), 便可得到三种物质的初始浓度值(C₂, C₃, C₄).然后在每次微增过程中，通过初始浓度可以计算每种物质成分在每次微增时的体积百分数，进而求出每次微增过程中的背景介质的体积百分数，再回代入公式(1)中去，便可得到每次微增过程的等效电导率值.