自21世纪初Baker Hughes公司成功研制出第一台三维感应测井仪器3DEX以来，有关三维感应测井响应的数值模拟、资料处理以及井下仪器的改进等已成为非常重要的研究课题 (Rabinovich et al., 2007).目前，三维有限差法 (沈金松，2003；Davydycheva et al., 2009；杨守文等，2009) 和三维有限体积法 (Hou et al., 2006；张烨等, 2012；王浩森等, 2016) 等数值方法、2.5维数值模式匹配算法 (Wang et al., 2012; 汪宏年等，2012) 以及解析方法 (Liu and Chew, 1992; Wang et al., 2008a, 2009; 魏宝君等，2009; 姚东华等，2010) 均已成为解决不同情况下三维感应测井数值模拟问题的重要方法.在仪器研制方面，目前在石油工业界已经出现了4种不同类型的三维感应测井仪器，例如Baker Hughes公司3DEX (Zhdanov et al., 2001)、Schlumberger公司的RT Scanner (Davydycheva et al., 2009)、Halliberton公司的3D-Induction Xaminer (Hou et al., 2013) 以及中国石油集团测井有限公司的TDIT (朱天竹等，2017).其中，中国石油集团测井有限公司 (CPL) 研制的三维感应测井仪器由3个相互正交的发射器、4个源距相对较短的轴向接收器和4个源距相对较长的三分量接收器组成，通过3个不同工作频率以及不同源距的单轴与三轴接收器可以同时提供96条视电导率曲线.从理论上说，利用这些测井曲线能够提取地层横向与纵向电导率、层厚以及地层倾角等信息，可以有效解决各向异性储层的评价问题.

然而，由于共面线圈系受井眼环境影响很大以及视电导率值与地层电导率间存在着较强的非线性关系，导致三维感应资料处理仍然面临诸多困难 (Wang et al, 2008b; Wang, 2011; 洪德成等，2013; Yang et al., 2014；周建美等, 2014).数值模拟结果以及井场实际测井资料均证明，多分量阵列感应测井资料不仅与地层厚度、地层横向与纵向电导率以及地层倾角与方位角等相关，而且整个测井响应同样受到井眼半径、泥浆电导率、仪器偏心距与仪器自旋角等因素的影响.在测井过程中仪器往往会发生偏心与自旋，然而，人们过去对偏心问题比较重视却往往忽略了仪器自旋对测井响应的影响.最近，我们在对多分量感应的实测资料进行处理时发现，同一仪器在同一口井中进行连续两次测量所得到的多分量视电导率曲线有时会存在较大差异.按照常规测井质量评价方法，这样的曲线应该是不合格的 (重复性差).在排除了仪器的稳定性和测量误差的原因后，基本可以断定曲线重复性差是由于仪器在下井作业过程中的偏心和自旋造成的.因此，必须分析研究仪器在偏心和自旋情况下的响应特征，解决测井曲线重复性差的问题.

本文拟利用2.5维数值模式匹配算法建立一套仪器偏心与旋转情况下TDIT测井响应的数值模拟方法，并通过数值模拟结果，对仪器偏心无自旋角、仪器偏心且存在自旋角两种情况下TIDI的响应进行对比分析，总结自旋角对TDIT响应产生影响的规律，为TDIT实测资料自旋角校正预处理和反演提供理论基础.

图 1a是包含井眼的水平层状非均质TI地层模型与仪器偏心示意图，其中，井眼半径和泥浆电导率分别为a和σ_m，地层n－1和地层n间的界面深度为d_n，地层n中横向与纵向电导率分别为σ_{h, n}和σ_{v, n}，仪器偏心距为ρ_Ecc.为简单起见，这里忽略了芯棒和侵入带的存在.图 1b是三维感应测井仪器 (TDIT) 线圈系结构示意图，整个仪器包含一个三轴发射线圈、4个短源距的轴向接收线圈 (A06、A09、A12、A15)、4个长源距的三轴接收线圈系 (A21、A27、A39、A72).图 1c是三轴发射 (三轴接收类似) 线圈系的结构示意图，其中，x和y方向的共面发射与接收线圈由对称分布在仪器芯棒两边的两个大小相同的线圈组合在一起形成的 (图 1c)，其目的是便于共面线圈系的安装并保证仪器居中情况下x和y方向的两个共面电磁响应相同.在进行多分量阵列感应电磁数值模拟时，当井眼与周围地层电导率反差较大时，如果仅仅用一个位于仪器芯棒中心的x和y方向的单一磁偶极子近似模拟图 1c所示的组合型共面线圈系的电磁响应，往往会产生较大的误差.为了取得高精度的多分量电磁响应结果，必须考虑到线圈系的这种实际分布对测井响应的影响.z方向线圈则是将导线直接缠绕在仪器芯棒上形成轴向线圈.为降低直藕信号影响，在发射线圈和每对接收线圈之间加装一组屏蔽线圈系，且屏蔽线圈缠绕方向与主接收线圈相反，L_R和L_B分别代表接收线圈系和屏蔽线圈系到发射线圈系的距离 (图 1d).仪器采用3个不同的工作频率依次进行测量，其低、中、高频分别为26.254、52.512、105.024 kHz.

图 1

Fig. 1

图 1 包含井眼的水平层状地层模型及TDIT结构 (a) 包含井眼的水平层状地层模型；(b) TDIT结构；(c) 三轴向线圈系结构；(d) 三轴向收发模块. Fig. 1 Schematic diagram of TDIT in horizontally layered TI formation with a borehole (a) Horizontally layered TI formation with borehole; (b) Schematic diagram of TDIT tool; (c) Structure of triaxial coil array; (d) Transmitter and receiver system of triaxial coil array.

本节将在利用2.5维模式匹配算法给出水平层状非均质各向异性地层中磁流源并矢Green函数的半解析解的基础上，根据线圈系的真实结构，应用叠加原理与旋转矩阵进一步给出仪器偏心以及存在自旋角情况下多分量阵列感应测井响应的计算方法.

在忽略仪器芯棒影响的情况下，图 1a中电导率的空间分布是轴对称的且具有如下的解析表达式为 (Wang et al., 2012):

(1)

其中，是地层n中的电导率函数，是地层n中的横向与纵向电导率组成的对角张量, 是阶跃函数，上标“±”表示为左右极限，是单位张量.为了计算仪器偏心情况下三维感应测井的真实响应，需要确定偏心磁流源并矢Green函数.为此，在柱坐标系下引入3个单位正交磁偶极子源激发的电磁场为

(2)

其中，和H^p(r, r_T)=(H_ρ^p, H_θ^p, H_z^p)^T分别是单位磁偶极子源产生的电场和磁场强度，将3个单位正交发射源产生的磁场组合在一起形成偏心磁流源磁场并矢Green函数是磁流源位置，μ为地层磁导率且假定为常数，ω=2πf为角频率，f为仪器工作频率且忽略位移电流影响.

由于电导率空间分布的轴对称性，方程 (2) 是典型的2.5维电磁正演问题，利用Fourier级数展开式，可以将该方程的求解问题转换为一系列轴对称问题 (Liu and Chew, 1992；Wang et al., 2012；汪宏年等，2012；朱天竹等，2017)，并利用数值模式匹配算法得到轴对称问题的半解析解，并将这些轴对称半解析解组合在一起得到方程 (2) 中磁场强度的半解析为

(3)

其中：

(4)

它是用数值模式匹配算法得到的地层n中磁场的第l阶谐变分量；

是径向基函数组成的2M或2M－1阶行向量，2M是模式匹配算法中径向节点的个数；是2M×2M(当l=±1时) 或2(M－1)×2(M－1)(当l≠±1时) 阶矩阵，而是M×2M或 (M－1)×2(M－1) 阶矩阵，A_l^p(n)(ρ_T, z_T) 是水平电场的l阶谐变分量的振幅；

和是电磁场在纵向上的解析解，是界面d_n+1上的广义反射系数，是地层n中的本征值方根组成的对角矩阵.

将方程 (3) 中柱坐标系下3个单位正交发射源磁场转化到直角坐标系中，可以得到直角坐标系下磁流源磁场并矢Green函数为

(5)

该式右端磁场强度的第一个下标代表接收线圈的法向方向，第二个下标代表发射线圈的法向方向.

三维感应测井仪器在井下测量过程中，不可避免的会发生偏心和自旋，不仅改变了仪器在井中的位置还影响了三轴线圈系的指向.图 2是仪器偏心情况下出现的两种可能的情况:仪器偏心但自旋角为零 (图 2a)、仪器偏心且自旋角不为零 (图 2b)，其中，ρ_Ecc和ρ_tool分别是仪器偏心距和仪器芯棒半径，仪器自旋角用φ表示.在图 2中建立了一个平面坐标系Oxy，其中，x轴位于井轴和仪器轴的垂直连线且坐标原点位于井轴，y轴也垂直于井轴，而z轴与井轴平行且垂直向下.当仪器偏心且自旋时，3个正交发射线圈的方向可以分别表示为T_x=(cosφ sinφ 0)^T、T_y=(－sinφ cosφ 0)^T和T_z=(0 0 1)^T.由于三分量接收线圈与发射线圈是平行的，3个正交接收线圈方向为R_x=T_x、R_y=T_y和R_z=T_z.此外，由于x和y方向的共面发射与接收线圈对称贴附在半径为ρ_tool的仪器芯棒两边，因此，x方向的两个发射线圈位置为r_{Tx, 1}=(ρ_Ecc－ρ_toolcosφ, －ρ_toolsinφ, z_T) 和R_{Tx, 2}=(ρ_Ecc+ρ_toolcosφ, ρ_toolsinφ, z_T)，y方向的两个发射线圈位置为r_{Ty, 1}=(ρ_Ecc－ρ_toolsinφ, －ρ_toolcosφ, z_T) 和r_{Ty, 2}=(ρ_Ecc－ρ_toolsinφ, ρ_toolcosφ, z_T).同样地，x方向的两个接收线圈位置为r_{Rx, 1}=(ρ_Ecc－ρ_toolcosφ, －ρ_toolsinφ, z_R) 和r_{Rx, 2}=(ρ_Ecc+ρ_toolcosφ, ρ_toolsinφ, z_R)，y方向的两个接收线圈位置为r_{Ry, 1}=(ρ_Ecc+ρ_toolsinφ, －ρ_toolcosφ, z_R) 和r_{Ry, 2}=(ρ_Ecc－ρ_toolsinφ, ρ_toolcosφ, z_R)，其中，z_T和z_R分别是发射与接收线圈所在的深度位置.在数值模拟过程中，为了尽可能减少Green函数的计算次数以便提高数值模拟效率，z方向发射线圈分别用与x、y方向的共面发射线圈位置相同的4个垂直方向磁流源近似，即z方向发射线圈的位置为r_{Tz, α}=r_{Tx, α}, r_{Tz, α+2}=r_{Ty, α}, α=1, 2；同样地，z方向接收线圈假定与x、y方向的4个共面接收线圈位置相同，这时z方向接收线圈的位置为r_{Rz, α}=r_{Rx, α}, r_{Rz, α+2}=r_{Ry, α}, α=1, 2.

图 2

Fig. 2

图 2 井下仪器偏心与自旋情况下发射线圈位置与方向 (a) 自旋角φ=0；(b) 自旋角φ≠0. Fig. 2 Schematic diagram of position and directions for transmitting coils at the condition of eccentricity and spin (a) Spin angle φ=0; (b) Spin angle φ≠0.

根据线圈系的这种真实分布情况下，利用叠加原理以及方程 (5) 给出的磁流源并矢Green函数的半解析表达式，可以给出发射线圈T_p(p=x, y, z)在接收线圈R_q(q=x, y, z) 上的磁场强度计算公式为

(6)

从公式 (6) 可以看出，只需要利用4个不同发射位置r_{Tx, 1}、r_{Tx, 2}、r_{Ty, 1}和r_{Ty, 2}上并矢Green函数在4个不同接收位置r_{Rx, 1}、r_{Rx, 2}、r_{Ry, 1}和r_{Ry, 2}上的值就能够计算出多分量发射与接收线圈上磁场强度的9个分量，有效提高了计算效率.TDIT收发线圈系统由一个三分量发射线圈、一个三分量屏蔽线圈和主接收线圈组成 (如图 1d)，因此，根据复合线圈系结构计算出复合线圈系上的磁场张量，并结合不同分量的仪器常数得到视电导率张量为 (朱天竹等，2017)：

(7)

其中，W_pq和H_pq(p、q=x、y、z) 分别为TDIT的仪器系数和合成磁场张量.

本节首先通过对比三层模型上2.5维模式匹配算法与三维有限体积算法的数值结果，检验2.5维模式匹配算法的有效性与计算效率；然后，通过一个包含井眼的7层水平层状各向异性地层模型，考察仪器偏心且自旋角φ=0°和φ=30°两种情况下三维感应测井的完整96条测井曲线，通过对比分析，了解仪器在偏心且不同自旋角情况下，不同源距的各个测井分量的响应特征与变化规律.

为了考察直井偏心情况下模式匹配算法的计算精度，选择一个包含井眼的三层水平层状TI地层模型，其井眼半径为0.1 m，井眼泥浆电阻率为0.5 Ωm，仪器沿x方向偏心且偏心距为0.05 m，并假定仪器自旋角为零，井眼外为三层水平分层TI地层，其横向电阻率分别为2.0、20.0、5.0 Ωm，纵向电阻率分别为6.0、80.0、12.0 Ωm，中间层厚度为2 m.TDIT采用3个不同频率、4个短源距轴向接收系统和4个长源距三轴向接收系统，可以同时采集96条不同阵列、不同频率和不同分量的测井曲线.因为计算不同阵列、不同频率的视电导率张量的核心算法类似，为了清楚起见，我们仅挑选出低频情况下阵列A21中的9个分量进行比较.图 3是仪器在偏心情况下，分别利用有限体积 (FV) 和2.5D数值模式匹配算法 (2.5D NMM) 计算的视电导率9分量的对比结果.从图中可以看出，两种方法计算的视电导率各分量均吻合得非常好，证明了利用NMM计算多分量阵列感应测井响应的有效性.利用FV同时计算上述地层模型中96条测井曲线大约需要13 h，而利用NMM算法只需2 h，说明NMM算法具有较高的计算效率.此外，从图 3中可以看出，当仪器沿x轴发生偏移且自旋角为零时，主分量σ_{a, xx}、σ_{a, yy}、σ_{a, zz}和交叉分量σ_{a, xz}、σ_{a, zx}对地层模型具有明显的响应，而剩下的4个交叉分量σ_{a, xy}、σ_{a, yx}和σ_{a, yz}、σ_{a, zy}恒为零.

图 3

Fig. 3

图 3 模式匹配法 (NMM) 和有限体积法 (FV) 的计算结果对比 Fig. 3 Comparison of methods NMM and FV

为了更系统的对TDIT 96条测井曲线进行综合分析，选取一个7层地层模型，其井眼半径为0.1 m，井眼泥浆为电阻率为100.0 Ωm的高阻泥浆.井眼外地层为水平层状TI地层，其横向电阻率分别为2.0、10.0、1.0、22.0、3.0、0.3、3.0 Ωm，纵向电阻率分别为2.0、40.0、2.0、80.0、12.0、1.5、3.0 Ωm，中间层厚分别为0.9、1.5、1.5、2.0、0.6 m.为了考察仪器偏心时，自旋角对测井响应的影响，分别计算了自旋角φ=0°和自旋角φ=30°时TDIT96条测井曲线，假定偏心距ρ_Ecc=0.05 m.

图 4是该模型在自旋角φ=0°时，TDIT的96条测井曲线.其中，图 4a1—a3是TDIT短源距收发线圈 (A06—A15) 的响应曲线，图例A06_H和A06_M分别代表高频和中频情况下阵列线圈A06的测井响应，其余图例类似.由于阵列线圈A06—A15的源距较小，因此响应曲线均具有较高的纵向分辨率.余下的图 4b1—d3中9幅子图，分别代表长源距三分量线圈系 (A21—A72) 的测井响应，图例A21_M和A21_L分别代表中频和低频情况下阵列线圈A21的测井响应，其余图例类似.9幅图分别对应公式 (7) 中视电导率张量的各个分量.从图中可以看出，除了视电导率分量σ_{a, xx}、σ_{a, yy}、σ_{a, zz}、σ_{a, xz}和σ_{a, zx}，其余4个视电导率分量σ_{a, xy}、σ_{a, yx}、σ_{a, yz}和σ_{a, zy}恒为零.

图 4

Fig. 4

图 4 TDIT 96条测井曲线 (ρEcc=0.05 m, φ=0°) Fig. 4 96 logging curves of TDIT (ρEcc=0.05 m, φ=0°)

图 5是在相同地层模型下自旋角φ=30°时TDIT的96条测井曲线.与图 4相比，测井响应曲线更加复杂，视电导率各分量全部不为零.同一个测井仪器在同样的地层模型中，仅仅由于仪器自身的自旋，其测井曲线发生了较大变化，给测井数据带来较大的不确定性.因此，为了更准确的得到地层模型参数，在对测井资料反演前，必须根据测井资料计算仪器自旋角，对测井资料进行自旋角校正.此外，由于A72阵列源距较大，其响应对层边界的响应不够灵敏，如果让A72阵列的测井响应参与反演，势必会给地层模型参数的反演带来较大的不确定度.因此，将来的反演将主要考虑A21、A27、A39阵列线圈的测井响应.

图 5

Fig. 5

图 5 TDIT 96条测井曲线 (ρEcc=0.05 m, φ=30°) Fig. 5 96 logging curves of TDIT (ρEcc=0.05 m, φ=30°)

本文利用模式匹配算法建立了水平层状非均质各向异性地层中仪器偏心且存在自旋角情况下，三维感应测井仪器响应的数值模拟方法，该算法能够提供96条不同阵列、不同频率、不同分量的视电导率曲线.数值结果表明，仪器偏心情况下，自旋角对TDIT响应的影响很大.当仪器偏心且自旋角为零时，每组阵列的张量视电导率的9个分量中只有3个主分量和两个次主分量不为零；而当仪器偏心且自旋角不为零时，每组阵列的张量视电导率的9个分量均不为零.由于仪器在井下测量时，自旋角是不可控的，在不同采样点处，其方位角可能完全随机，所以，自旋角的存在给资料的解释和处理增加了难度，如何从测量资料中提取出自旋角角并进行自旋角校正，对于简化资料的处理是非常重要的.