0 引言

传统的单井感应测井技术有着良好的纵向分辨率，能够对井眼及其周围地层的地质特性进行精细描述，但固有的横向探测能力不足的缺点导致其无法对储层构造、裂缝发育方向等重要信息进行探测.井间电磁测井作为电法测井技术的前沿，从根本上解决了单井技术的横向探测能力不足问题，实现了“井孔”与“井间”的信息平衡(栗建军, 2004).井间电磁测井技术的基本原理是将发射线圈置于一口井中，向地层中激发频率为5 Hz~1 kHz的交变电磁场，由电磁感应作用使地层中形成感应涡流产生二次感应电磁场，利用置于另一口井中的接收线圈获取二次场信号来分析井间电导率二维乃至三维构造信息.近二十年来，井间电磁测井技术应用及设备研制方面的研究取得重大进展，如EMI公司的XBH-2000 (曾文冲等, 2001; 臧德福等, 2013)和斯伦贝谢公司的Deep-Look(Zaki et al., 2009; 臧德福等, 2013)都已成功地应用于勘探生产实践中.

砂泥岩薄交互层、裂缝性储层(沈金松, 2009)作为典型的电导率各向异性储层亟待被开发利用，但采用传统轴向线圈的井间测井仪器难以对各向异性储层进行探测.例如，斯伦贝谢公司在利用Deep-Look井间电磁成像系统对巴西油田勘探时发现，当井眼发生倾斜时，地层电导率各向异性会导致井间电磁响应曲线发生“畸变”现象，使测井数据解释产生严重偏差(Beer et al., 2010).目前，对于地层各向异性的测量主要采用多分量感应测井技术，其原理是利用三轴正交发射与接收线圈，同时测量单井中每个深度点上磁场的九个分量，进而提供井眼附近地层的横纵电导率信息(Kennedy et al., 2001).若能将三轴线圈系应用在井间电磁测井中，深入研究电导率各向异性储层中多分量井间电磁响应特征，对正确解释井间电磁测井资料极为重要.

目前，数值仿真是分析井间电磁响应特征的有效手段.近些年，国内外学者在井间电磁响应的数值正演算法及模拟计算方面开展较深入研究，并取得一定成果.Zhdanov建立了不含井眼条件下三维井间电磁响应的区域准线性近似方法，并以此为基础建立快速三维反演算法(Zhdanov and Yoshioka, 2003).Donadille利用二维有限元方法模拟了裂缝性地层的井间电磁响应，同时指出裂缝方向和位置可作为井间电磁数据资料解释的辅助信息(Donadille and Al-Ofi, 2012).杨曦利用FDTD算法模拟基于电偶极子的二维轴对称模型井间电磁响应(杨曦和潘和平, 2008).常欣莉等人利用体积分方程方法模拟了倾斜及层状各向异性地层的井间电磁响应，重点分析了井眼倾角和背景层各向异性系数对测井响应的影响(常欣莉等, 2015).目前，井间电磁数值模拟研究大多是二维或不含井眼的简单三维模型，且对各向异性介质中的电磁响应模拟研究相对较少.目前井间电磁响应的三维频域正演研究面临以下两方面问题：首先是低感应数问题，井间电磁测井属于低频测井，若采用直接离散电磁场的数值算法，会导致数值迭代算法收敛缓慢(Weiss and Newman, 2003).其次是计算规模庞大问题，在井间电磁数值计算中既要保证空间求解域的尺度又要兼顾井眼周围的计算精度势必会大幅增加离散网格数和计算时间.解决低感应数问题的有效手段是采用低感应数预处理技术，其本质是将电场分解成一个无源场(矢势)和无旋场(标势)叠加形式，建立关于矢势与标势的耦合方程组，其离散方程组系数矩阵具有更小的条件数，因此利用迭代算法求解离散方程组具有强收敛性.该预处理技术适用于有限差分、有限体积以及有限元等数值方法，并已应用于多分量感应测井(张烨等, 2012; Hou et al., 2006; Sun and Zhao, 2008)、随钻测井(王浩森等, 2016; Novo et al., 2010)和可控源大地电磁(周建美等, 2014; 陈辉等, 2016; Haber et al., 2000)的数值计算中.对于计算规模问题通常采用非均匀网格技术来降低对网格数量的需求.此外，采用并行计算技术可有效提升数值方法的计算速度.

本文基于耦合势方法在低频电磁勘探数值仿真中的优点，首先建立含井眼的横向同性地层中井间电磁测井响应的耦合势有限差分离散方程组，并利用不完全LU分解预处理技术和稳定双共轭梯度法实现对离散方程组的快速求解.然后，通过与有限元软件COMSOL计算结果对比，验证算法有效性.最后，利用数值仿真结果重点研究三维横向同性地层中多分量井间电磁测井响应特征，并系统考察地层各向异性、发射源位置以及宏观异常体等因素对井间电磁测井响应不同分量的影响.本文数值仿真结果可为各向异性地层中井间电磁测井资料的解释以及测井仪器的改进提供重要数据支撑和理论基础.

1 研究方法

本文首先建立井间电磁测井模型所满足的Helmholtz方程及边界条件，然后利用耦合势的Yee氏交错网格对模型求解区域进行离散，并结合预处理技术、迭代算法以及多核并行技术实现对超大型稀疏线性方程组的快速运算.

1.1 频率域耦合势Helmholtz方程及其边界条件

非均质各向异性地层中频域Maxwell方程表示为

(1a)

(1b)

ω是发射线圈圆频率，对于各向同性的非铁磁性介质，磁导率用真空磁导率μ₀表示，ε代表介电常数，I是单位张量，Mδ(r-r_s)表示位于r_s处的磁偶极子源.在横向同性地层中σ为3×3对角张量，若考虑到位移电流效应，引入复电导率：

其中σ_h^*，σ_v^*分别是地层的水平电导率和垂直电导率.

将(1b)代入(1a)，得到Helmholtz方程：

(2)

在井间电磁测井发射源频率范围下(＜1 kHz)，|ωμ₀σ^*|≪∇× ∇×E，则(2)式近似为

直接离散∇×∇×E，会得到相对病态的离散化系数矩阵，这给迭代算法的收敛稳定性带来严重困难，形成低频电磁勘探数值计算中经常遇到的低感应数问题.

解决低感应数问题的有效方法是利用势函数描述电磁场，引入矢势A和标势ϕ，将电场标势为矢势与标势梯度叠加的形式:

(3)

由于矢势描述电磁场的不唯一性，需要对矢势A施加规范条件.本文采用库伦规范条件：∇·A=0.于是(3)式的物理意义是将电场表示为无源感应电场A和有源电场∇ϕ叠加的形式.由(1a)可得，磁场强度：H=i∇×A/ωμ₀.

将式(3)代入式(2)，得到基于矢势和标势的Helmholtz方程：

并利用∇×∇×A=-∇²A，得

(4)

为了得到一个完整并且对角占优的线性方程组，我们对(4)两边求散度得到

(5)

其中∇·[∇×Mδ(r-r_s)]=0.联立(4)和(5)建立耦合势线性方程组：

(6)

由于内存空间和计算资源的限制，实际数值计算只能在有限大小的六面体区域Ω内进行，在求解域外边界∂Ω施加截断边界条件.本文采用完美电导体边界条件(PEC condition)，即

由此，得到矢势、标势的截断边界条件：

(7)

1.2 三维耦合势有限体积离散

首先基于耦合势的Yee氏交错网格对求解域进行区域分解，规定离散分量在Yee氏交错网格的位置遵循的原则是：矢势分量位于网格面的中心点上，离散标势定义在网格的体中心处.然后以离散分量为中心建立四种类型控制体积单元，每种控制体积单元对应的区域和体积为

在离散的过程中认为控制体积V_{i, j, k}内部电导率是均匀的，然而在体积单元之间的边界上电导率是不连续的，甚至可能反差很大.于是V_{i, j, k}内部等效电导率表示为

其中i指h, v, 代表水平或垂直方向，n为控制体积内部介质的种类数.

下面以(4)式中的x分量方程离散为例介绍有限体积离散方法，首先对其在控制体积V_{i+1/2, j, k}上做积分平均：

(8)

由高散度定理，将(8)式等号左侧的第一个体积分可转化为如下的面积分，利用基于体积元上六个面的中值积分公式和一阶差分公式，将其离散为如下形式：

对第二个体积分进行直接离散，得

其中〈σ_h〉_{i+1/2, j, k}是两体积元V_{i, j, k}和V_{i+1, j, k}边界上的水平电导率平均值，为符合电导率串联法则，该电导率用相邻体积元等效电导率的调和平均值确定：

(4) 式y，z分量方程离散方法与x分量一致，值得注意的是，对于z分量离散方程，电导率为地层的垂直电导率.方程(8)右端的源项离散结果为

对于(5)式在控制体积V_ijk上积分平均：

并根据高斯散度定理得到

其中J=σ^*(A+ ∇ϕ)为传导电流密度矢量，于是(5)式可以离散成

其中：

1.3 离散方程组建立及求解

式(4)和式(5)的离散结果进行排列和整理，得到关于矢势A和标势ϕ的离散化线性代数方程组：

通过分析(4)式离散方程发现，对于每个离散矢势分量(以A_z(i, j, k+1/2)为例)，一共有8个离散分量与其相关(其中包括6个离散矢势分量和2个标势分量)，并且所有相关分量均以A_z(i, j, k+1/2)为中心对称分布，如图 1a所示.同样地，(5)式中每个标势分量都有12个相关离散分量，围绕其对称分布，如图所示.由此可以推断出：系数矩阵中前三行子矩阵中最多包含9个非零元，第四行子矩阵最多包含13个非零元，如图 1b所示.因此，系数矩阵中的非零元个数远小于矩阵尺度，是一个超大型稀疏非对称矩阵，Mumps和Pardiso等直接算法由于内存需求量太大无法求解离散方程，本文采用不完全LU分解的稳定双共轭梯度法BICGSTAB对其进行求解(Saad, 1996).在编程的过程中利用CSR压缩存储格式对稀疏矩阵进行压缩存储.此外，采用Openmp多核并行编程技术(Hermanns, 2002)提升算法的计算速度.

2 数值仿真结果及讨论 2.1 井间电磁模型建立及网格划分

多分量井间电磁测量原理如图 2所示，为了能够研究地层电导率的各向异性特征，本文假设发射线圈T和接收线圈R均为三轴正交指向，并且线圈中心与井轴重合，发射线圈频率选择范围在0.1~1 kHz, 并根据电磁场强度在有损介质中的衰减特性，将计算区域选择为最低发射频率的三倍趋肤深度，其中σ_min为模型中电导率的最小值.以保证在截断边界上场值接近于0.图 2中圆点代表数据采集时发射线圈和接收线圈的位置点.

为了提升数值计算效率，本文采用非均匀直角网格对求解域进行划分，基于井间电磁灵敏度分析结果(Spies and Habashy, 1995)，在井眼及发射-接收线圈附近的高灵敏度区域采用小尺寸均匀网格，网格长度在水平方向与竖直方向分别为0.1 m和0.5 m.而在低灵敏度区域网格尺寸呈指数增加，如图 3所示.

2.2 算法验证

为验证耦合势有限差分算法在计算三维井间电磁响应中的有效性，与商用有限元软件COMSOL的数值仿真结果进行对比.验证模型选取与图 2模型相类似的含井眼但不含异常体的三层地层模型，其中间层为层厚40 m的水平横向同性层，层边界位置为120 m和160 m, 水平和垂直电导率分别为σ_h=0.1 S·m^-1，σ_v=3 S·m^-1.第一层和第三层为电导率同性地层，电导率为σ_h=σ_v=1 S·m^-1，井眼电导率为σ_bh=0.001 S·m^-1，发射-接收井距为100 m.假设发射线圈工作频率为100 Hz，位于中间层121 m靠近层边界处，其线圈指向为水平x方向.在仿真计算中求解域范围是400 m×400 m×250 m，网格个数为140×140×180个.

图 4中散点线和实线分别代表由有限体积和有限元算法计算得到磁场虚部.可以看出三个接收方向的磁场曲线几乎完全重合，经计算两者平均偏差低于0.5%，由此验证本文算法在三维井间电磁模型中的有效性.此外对比三轴磁场分量大小可以看出，H_xy分量要比H_xx和H_xz小4~5个量级.

2.3 多分量井间电磁响应特征

对于z方向发射线圈，在接收井中磁场虚部的三个分量如图 5a所示.在井轴所在的x-z平面上，没有y方向磁场分量，因此响应分量H_zy=0 A·m^-1.磁场H_zz分量曲线比较平缓，在深度为100 m附近曲线达到最大值，而不是发射源位置120 m处出现最大值.产生这一现象的原因是：中间层的水平电导率小于两侧地层的电导率, 在两侧地层离发射源深度较近的区域内产生的感应电流较强，因此磁场峰值发生向高电导率层偏移的现象.磁场H_zx分量对层边界非常敏感，在层边界处出现两个拐点，这是由感应电流在界面上不连续造成的.根据电磁场的边值关系式，感应电场在分界面上是连续的，但由于电导率的不连续，根据欧姆定律的微分形式J_i=σE_i可知，感应电流在分界面上不连续，进而造成了二次场的不连续.对于x方向发射源，三轴接收线圈的磁场如图 5b所示.磁场的H_xy分量近似为零，H_xx、H_xz分量不为零.其中H_xz的曲线形态类似正弦型.由于发射源的位置处于地层一和地层二的交界面附近，不是地层模型的中心位置，因此正负响应的强度和深度范围均不对称.H_xx分量主要表现为正响应，其曲线在层边界处也出现拐点，说明H_xx分量对层边界也比较敏感.最后，y方向发射线圈对应的磁场三分量如图 5c可见，H_yx和H_yz近似为零，H_yy的曲线比较平滑且没有负响应，但对层边界的敏感程度不如分量H_xx，只能在曲线中分辨出160 m深度的层边界.

2.4 地层各向异性对测井响应各分量的影响

为考察各向异性对测井响应的影响，建立三层地层模型，其中一、三层为电导率各向同性地层，地层电导率为1 S·m^-1.第二层为横向同性地层，横向电导率为0.1 S·m^-1，垂直电导率分别选取0.1 S·m^-1、0.5 S·m^-1、1.0 S·m^-1、2.0 S·m^-1、3.0 S·m^-1.从图 6(a-c)中可以看出测井响应的H_zz、H_zx和H_xz分量与地层各向异性无关.由Maxwell方程(1b)可知：对于z方向磁流源，磁场只与地层横向电导率有关，而与地层垂直电导率无关.对于x方向磁流源，磁场响应中的H_xx和H_xy分量与地层电导率的垂直分量有关，而H_xz分量与地层水平分量有关.同样地，对于y方向磁流源，磁场响应的H_yy和H_yx两个分量与地层垂直电导率有关.由此可见，H_xz、H_zx和H_zz的曲线在不同垂直电导率的条件下完全重合，如图 6(a—c)所示.而磁场响应H_xx和H_yy曲线不重合，如图 6(d、e)所示，H_yy随着地层垂直电导率的增加而升高，当垂直电导率大于2 S·m^-1后，曲线无明显上升趋势.而H_xx随着地层垂直电导率的增加而降低，且下降幅度逐渐减小.

2.5 发射源位置变化时井间电磁测井响应特征

发射源在40~200 m区间范围运动，发射点间距与测量点间距均为4 m，多分量井间电磁响应随发射-接收位置变化特征由图 7所示.H_zz分量没有负响应，当接收线圈与发射线圈位于同一深度时，H_zz分量达到极大值，且在低电阻层中H_zz分量强度要高于高电阻层.而H_zx分量和H_xz分量的正负响应以对角线近似对称分布，且在中间层边界处出现明显拐点，因此从交叉分量可以清楚分辨出层边界位置.由于H_xx分量与H_yy分量均能反映地层的垂直电阻率信息，所以在第二层中H_xx分量和H_yy分量出现较强的正响应，说明该层的垂直电阻率较低.此外，H_yy分量略强于H_xx分量.

2.6 异常体对井间电磁测井三轴主分量的影响

在三层井间模型的各向异性层中，添加高阻异常体，其长宽高均为35 m，且中心位于两井的中心位置上，异常体电导率为0.01 S·m^-1.图 8为有、无异常体两种情况下的井间电磁三主分量图.通过对比发现，H_xx分量和H_yy分量对异常体几乎没有任何反应，表明其存在井间区域分辨率不足的缺点，而H_zz分量在异常体存在深度上，相比于无异常情况，磁场强度明显升高，这说明H_zz分量在井间区域具有更好的分辨率.

3 结论

本文建立了基于横向同性地层的多分量井间电磁响应的三维耦合势有限体积差分算法.通过与有限元软件COMSOL的仿真结果进行对比验证该算法的有效性.

数值仿真结果表明：井间三轴发射与接收线圈电磁响应包含丰富的地层信息，在垂直井模型中除主分量xx，yy，zz和交叉分量xz，zx不为零外，其余交叉分量的测井响应均近似为零.在非零测井响应分量中xx和yy分量对地层各向异性比较敏感，随着地层垂直电导率的增大，xx分量响应呈下降趋势而yy分量则呈上升趋势，而主分量zz和交叉分量xz，zx对地层各向异性不敏感，因此xx和yy分量可用于有效判断地层背景电导率各向异性特征.交叉分量xz和zx对地层边界异常敏感，具有极高的纵向分辨率，因此可用于分辨层边界的位置.对于传统井间电磁测井所使用的zz分量，相比于xx，yy分量和主分量，其对于井间区域的异常体具有很强的灵敏度和更高的横向分辨率，因此可利用该特性探测井间油气藏构造或检测井间裂缝发育方向.

本文结果可以为井间电磁的三维反演提供快速正演方法，也可为井间电磁成像奠定基础.