1 引言

砂泥岩薄交互层是典型的低阻油气储层，其电导率存在非常明显的各向异性^[1-2]，如果能够同时测量出纵横向电阻率，对于该类储层的正确评价具有非常重要的意义.多分量感应是新一代测井技术，能够同时提供地层纵横向电阻率信息^[3-4]，与各种常规感应技术(双感应，阵列感应)相比，多分量感应不仅含有方向与井轴平行的发射和接收线圈(共轴线圈系)，而且还增加了方向与井轴垂直的发射和接收线圈(共面线圈系)，整个仪器由三个相互正交的发射线圈(T_X､T_Y和T_Z)和接收线圈(R_X､R_Y和R_Z)组成，利用这种三轴发射和接收系统，能够同时测量9个不同分量的感应磁场，提供地层纵横向电阻率､井眼相对倾角等有用信息.

由于共面线圈系的电磁场不再具有轴对称性，往往需要采用三维数值模拟技术研究多分量感应仪器的响应^[5-8]，其运算效率往往较低.因此，效率更高的解析法^[9-11]与混合法^[12-18]也得到了较广泛研究与应用.在文献^[15-16]中，我们已经用模式匹配法研究建立了一套多分量感应测井的正演与反演算法，在不含井眼和侵入带的水平层状各向异性模型上取得了满意效果.然而，当地层模型存在井眼或侵入带时，由于共面线圈系产生的电场存在着径向分量，电流密度的连续性会在柱状电阻率分界面上产生积累面电荷.当用模式匹配算法对多分量感应测井响应进行数值模拟时，由于刚度矩阵是非对称的复矩阵，矩阵元素的微小变化往往对特征值和特征向量均产生较大影响，如果忽略这些积累面电荷对电磁场的影响将得不到正确的数值模拟结果.

本文将对文献^[13-16]中的模式匹配算法进行改进，以便建立水平层状非均质各向异性地层(含有井眼和侵入带)中多分量感应测井响应快速数值模拟算法.首先利用电导率导数在柱状界面上的奇异性，引入两个附加奇异微分算子，用于描述井壁以及侵入带与原状地层分界面上的积累面电荷对电磁场的影响，然后再应用模式匹配算法求解含有附加奇异算子的偏微分方程，得到水平层状非均质各向异性地层中电磁场半解析解以及多分量感应响应有效算法.在此基础上，通过数值模拟结果研究考察井眼泥浆电阻率､侵入深度以及层厚等变化对阵列多分量感应测井响应的影响特征.

2 理论

图 1是水平层状各向异性非均质地层，地层总层数是N，其中，第n层原状地层(侵入带)横向和纵向电导率分别为σ_HT，n和σ_VT，n(σ_HX，n和σ_VX，n)(n=1，2，…，N)，侵入半径和水平界面位置分别是r_xo，n(n=1，2，…，N)和d_n(n=2，3，…，N)，井眼半径和井眼泥浆电导率分别为a和σ_md.图中还建立了相应的直角坐标系xyz，其中，z轴与井轴重合，三轴单位发射线圈用表示且与三个坐标轴方向完全一致，三个不同源距的三轴接收线圈用(R_xk，R_yk和R_zk，k=1，2，3)表示，构成一个简单的阵列多分量感应测井仪.多分量感应测井响应通过求解如下的Maxwell方程得到^[13-15]:

(1)

方程(1)中，已假定三轴发射线圈${{\hat{e}}^{p}}$(p=x，y，z)具有单位磁矩，且均匀分布在半径为ρ_s 的小圆环上，其目的是保证井轴周围的电磁场仍然可解，δ 是Dirac-δ 函数.E_p和H_p表示发射线圈${{\hat{e}}^{p}}$产生的电磁场，珔σ 是电导率张量，μ 是磁导率，Ω 为发射信号的角频率，且忽略位移电流影响.

由于电导率是分片常数函数，第n层中的纵横向电导率函数可表示为^[18]

(2)

其中，下标“P= H，V"分别表示横向与纵向方向，是Heaviside函数，上标“±"分别对应变量的上下极限.对方程(2)两端求导，则横向与纵向电导率导数为

(3)

(3) 式表明在井壁以及侵入带和原状地层交界面上，电导率导数是奇异的.

2.1 水平发射线圈的电磁场

利用直角坐标系(xyz)中，单位正交向量${{\hat{e}}^{x}}$和${{\hat{e}}^{y}}$与柱坐标系(ρ，θ，z)中单位正交向量${{{\hat{e}}}_{\rho }}$和${{{\hat{e}}}_{\theta }}$的关系${{\hat{e}}^{x}}$=cosθ${{{\hat{e}}}_{\rho }}$ -sinθ${{{\hat{e}}}_{\theta }}$ 和${{\hat{e}}^{y}}$=sinθ${{{\hat{e}}}_{\rho }}$ +cosθ${{{\hat{e}}}_{\theta }}$，并与欧拉公式cosθ=0.5(e^iθ+e^-iθ)和sinθ=-0.5i(e^iθ -e^-iθ)结合，可以将直角坐标系中的两个单位向量${{\hat{e}}^{x}}$和${{\hat{e}}^{y}}$转换到柱坐标系中，并可展开成如下的Fourier级数:

(4)

其中，M_l^x=0.5(1 il)^T 和M_l^y=-0.5il(1 il)^T分别是单位向量${{\hat{e}}^{x}}$和${{\hat{e}}^{y}}$在柱坐标系中的Fourier展开式系数(l阶谐变分量).

将(1)式中水平磁流源${{\hat{e}}^{x}}$和${{\hat{e}}^{y}}$产生的电磁场也转换到柱坐标系中，并将电磁场Ep= (E_ρ^p，E_θ^p ，E_z^p)^T和Hp= (H_ρ^p，H_θ^p ，H_z^p)^T(p=x，y)也展开成如下形式的Fourier级数:

(5)

将(4)式和(5)式代入(1)式，并对方程进行分解，得到水平电场的l阶谐变分量E_{s, l}^p= (σ_HE_ρ，l^p，E_θ，l^p)^T 满足的方程(见附录):

(6)

其中${{{\tilde{L}}}_{E,L}}$是2×2阶的微分算子，定义为

(7)

(7) 式中的第二项为附加奇异算子，起因于(3)式中电导率导数的奇异性，主要反映柱状界面上积累面电荷对电场的影响.

(8)

为保证方程(6)在井轴附近仍然可解并保证其解存在和唯一，需要如下的边界条件^{[9， 15， 17]}:

(9)

其中，ρ_MN 是足够小的数(如10^-6m)，ρ_MX 是径向上的最大求解范围.

磁场强度水平方向和垂直方向的l阶谐变分量H_s，l^p= (H_ρ，l^p，H_θ，l^p)^T 和H_z，l^p与E_s，l^p= (σ_HE_ρ，l^p，E_θ，l^p)^T间的关系方程^[13-14]:

(10)

(11)

是2×2阶的微分算子且右端的第二算子也是奇异微分算子.

2.2 垂向发射线圈的电磁场

垂向单位发射线圈${{{\hat{e}}}^{z}}$是轴对称的，其产生的电场E^z中只有一个非零分量E_θ，0^z､磁场H^z中仅有两个非零分量H_ρ，0^z，H_z，0^z，用同样方法可以从(1)式中得到电场E_θ，0^z 满足的方程:

(12)

对应的边界条件^{[9， 15， 17]}:

(13)

磁场的两个非零分量H_ρ，0^z,H_z，0^z与电场E_θ，0^z的关系为^[13-14]

(14)

由于(12)~(14)式的求解过程与(6)~(11)式完全类似且前者更为简单，因此下面仅给出(6)~(11)式的求解过程.

2.3 完全柱状介质中的电磁场

为了研究(6)~(11)式的快速求解方法，首先假定模型中不含水平界面(即完全柱状介质)，并用有限长区间[ρ_MN，ρ_MX]逼近半无限区间，并将其剖分成M(为偶数)个小区间､剖分节点位置为ρ_α，α=1，2，…，M+1.选用二次Newton插值函数作为径向节点ρ_α上的基函数φα(ρ)，(α = 1，…，M)，基函数定义如下^[19]:

在内边界节点ρ=ρ₁ =ρ_MN(α=1)上:φ₁(ρ)=；

在其他奇数节点ρ =ρ_2k-1(α =2k-1)上:

在偶数节点ρ =ρ2k(α =2k)上:

利用Galerkin法，将方程(6)中的水平电场的l阶谐变分量展开成如下形式^[14-15]:

(15)

这里的是已知的2×2M阶基函数矩阵，S₁(ρ)= (φ₁(ρ)，…，φ_M-1(ρ)，φ_M(ρ))^T 是由径向节点上的基函数构成的M维列向量；E_l，α 是未知的2M维列向量，Q_l，α(z，z_s)和c_l，α ^p分别是未知函数与未知常数.将(15)式代入(6)式，两边同乘以矩阵S(ρ)ρ 后，关于ρ 在区间[ρ_MN，ρ_MX]上进行积分:

(16)

其中，${{{\tilde{B}}}_{l}}$和${\tilde{D}}$分别是2M×2M阶刚度矩阵和质量矩阵，$\tilde{F}_{l}^{p}$是2M阶列向量，分别定义如下:

(17)

通过分离变量，从(16)式可以得到如下三个独立方程:

(18)

(19)

(20)

求解(18)式，得到2M个本征值κ_l，α²和相应的本征向量E_l，α，(α =1，…，2M).方程(19)的解:

(21)

进一步求解方程(20)确定向量C_l^p= (c_l，1^p，c_l，2^p，…，c_l，2M^p)^T是本征向量矩阵.最后得到完全柱状介质中水平电场的l阶谐变分量的半解析解.

(22)

其中，=diag(κ_l，1，κ_l，2，…，κ_l，2M)是2M个本征值κ_l，α²的平方根组成的对角矩阵且Re(κ_l，α)>0.

2.4 水平层状非均质地层中的电磁场

利用文献^[14-15]的结果，第n层非均质地层中的p-方向的磁偶极子产生的水平电场的l阶谐变分量可表示为如下形式:

(23)

其中，A_l^p(n)±(ρ_s, z_s)是电场振幅，

(24)

(24) 式中符号表示各个地层界面上的广义反射系数.

将(23)式代入(10)式中，经整理，得第n层中的水平和垂直磁场的l阶谐变分量(H_s，l^p和H_z，l^p):

(25)

(26)

其中，，而矩阵和通过如下方程确定:

(27)

(28)

最后，得到由x-方向的单位水平磁偶极子${{\hat{e}}^{x}}$产生的磁场强度:

(29)

y-方向水平磁偶极子${{\hat{e}}^{y}}$产生的磁场强度，由(17)式中的第三个公式得:

(30)

设r_R= (ρ_R，θ_R，z_R)为接收线圈位置，将柱坐标系下的解(29)式转换到直角坐标系中，得磁场强度主分量H_xx和H_yy:H_xx=H_ρ^x(r_R；r_s)cosθ_R-H_θ^x(r_R；r_s)sinθ_R和H_yy= H_ρ^y(r_R；r_s)sinθ_R+H_θ^y(r_R；r_s)cosθ_R.

采用类似方法求解(12)~(14)式，可以确定垂直磁偶极子${{{\hat{e}}}^{z}}$产生的磁场强度H_z(r_R；r_s)，其主分量记为H_zz= Hzz(r_R；r_s).

2.5 多分量感应视电导率

利用磁场的三个主分量H_xx､H_yy和H_zz，并按照公式^[4]:和，计算出多分量感应仪器的测井响应，其中，L是发射线圈和接收线圈间距离.

3 数值计算结果

本节首先通过柱状介质模型中多分量感应响应的解析解^[9]与本文的NMM 解的对比，检验NMM的有效性，然后，通过数值结果考察阵列多分量感应测井响应.假定阵列多分量感应仪由一个三分量的发射线圈和三个源距分别为0.8､1.2 m 和1.8 m三分量接收器组成(见图 1)，发射线圈的工作频率为50kHz和200kHz.假定井眼半径a=0.1m, 由于在垂直井眼中σ_A^XX=σ_A^YY，所以本节中的数值结果仅包含两个主分量σ_A^XX和σ_A^ZZ.在数值模拟过程中，方程(6)~(14)中的径向求解区间的范围为ρ_MN =10^-6m 和ρ_MX =60m, 总节点数M+1=119，其中，节点ρ_α(α = 1，2，…，40)按等间距增加且间距为0.005m, 而节点ρ_α(α=40，42，…，M+1)按对数等间距增加.

图 2 是一个完全柱状各向异性地层中，NMM和解析法^[9]计算得到的不同源距条件下多分量感应的两个主分量σ_A^XX 和σ_A^ZZ 以及归一化后的三个电场强度(E_YX，E_XY和E_YZ)的实部对比结果.该模型的纵(横)向电阻率R_VT(R_HT)是90(20)Ωm, 泥浆电阻率R_MD 分别为0.1和5Ωm, 且ρs=0.005m､接收线圈位置为(xR，yR)= (0.01m, 0).图 2a和图 2b分别是σ_A^XX 和σ_A^ZZ 结果的对比，图 2c和图 2d是归一化的电场强度分量E_YX､E_XY和E_YZ的实部对比情况，其中实线和虚线是NMM 数值结果而各个离散符号对应于解析法结果，结果清楚显示NMM 和解析法两种方法所得结果吻合得非常好.

对于含有井眼和侵入带的水平层状非均质各向异性地层中NMM 算法的检验，在文献^{[8， 17]}中已经与三维有限体积差分算法的数值结果进行了详细对比，故这里从略.

为考查层厚和泥浆电阻率变化对阵列多分量感应测井影响，并了解仪器的纵向分辨率大小，假定模型厚度分别为1､1.5､2､3m 和4m 的高阻层与厚度均为3m 的低阻层组成的交互层，上下围岩和低阻层的横向､纵向电阻率为1/(σ_HT，σ_VT)= (2，6)Ωm, 而高阻层横向､纵向电阻率为1/(σ_HT，σ_VT)= (25，85)Ωm.图 3是泥浆电阻率分别为Rmd=0.1､5.0Ωm时阵列多分量感应测井仪器的正演结果.数值结果清晰显示，源距､工作频率､层厚以及泥浆电阻率变化等，均对测井响应产生明显的影响.如果细心观测数值结果可以发现几个有趣的现象，首先分析图 3a中低阻泥浆情况下共面线圈系的主响应σ_A^XX ，不难看出，源距增加会导致整个视电导率增大而频率增加导致视电导率减小，且在低阻层的层界面附件出现大量的负响应；σ_A^XX 对不同厚度的高阻层均有明显反映，且在厚度2m 以下的薄层上，由于受围岩影响较大，其视电导率明显大于厚层上的视电导率，此外，高阻地层上几乎没有负响应.图 3c中高阻泥浆对应σ_A^XX 响应，其变化规律与图 3a中低阻泥浆情况有很大不同，高阻泥浆对应的视电导率σ_A^XX 明显大于低阻泥浆时的视电导率，并且源距增加会导致整个视电导率降低.特别需要指出的是，不论是在低阻还是在高阻泥浆情况下，其视电导率曲线σ_A^XX 均难以直接指示出地层电导率的大小，且共面线圈系的纵向分辨率最多为1 m.图 3b和3d是低阻和高阻泥浆情况下共轴线圈系响应σ_A^ZZ ，共轴线圈系的响应特征要简单得多，从曲线的变化形态能够较易识别各个地层上电导率的相对大小，且源距(或频率)增加，会导致电导率地层上的视电导率减小､电阻率地层的视电导率增加，电导率泥浆和电阻率泥浆下的测井响应有一定的差别，但没有共面线圈系情况明显.图 3中含有314个深度点，在2.5GHz双核笔记本电脑上完成两个不同电阻率泥浆情况下的所有计算，使用的CPU 时间为4340.8+4472.4s.

最后，为考查阵列多分量感应仪器的探测深度，假定上下围岩和导电率地层横向､纵向电阻率均为1/(σ_HT，σ_VT)= (2，6)Ωm, 导电率地层厚度均为3m, 5个电阻率地层厚度均3.5m, 其原状地层和侵入带横向､纵向电阻率均为1/(σ_HT，σ_VT)= (25，85)Ωm和1/(σ_XT，σ_VX)= (12，35)Ωm 且从上到下其各地层的侵入半径分别是0.2､0.4､0.6m､0.9m和1.2m.图 4是泥浆电阻率分别为R_md=0.1和5.0Ωm时，阵列多分量感应测井响应的数值模拟结果.从图 4b中低阻泥浆情况下共面线圈系的主响应σ_A^XX 曲线同样可以看出，源距增加会导致整个视电导率增大而频率增加导致视电导率减小，且在低阻层的层界面附件出现大量的负响应；高阻地层上负响应明显降少，且不同侵入深度地层上的视电导率曲线的差异也很小，从曲线形态差异和其视电导率的大小难以判断侵入深度变化对于测井响应的影响，不过，如果对视电导率采用对数表示，在侵入半径1.2m 的地层上，仍然可以看出侵入带影响.图 4d是高阻泥浆情况下共面线圈系的主响应σ_A^XX 曲线，同样地，高阻泥浆对应的视电导率σ_A^XX 明显大于低阻泥浆时的视电导率，并且源距增加会导致整个视电导率降低.图 4c和4e给出的低阻和高阻泥浆情况下的共轴线圈系响应σ_A^ZZ ，对侵入深度的变化有非常明显响应.图 4中含有394个深度点，在2.5GHz双核笔记本电脑上完成两个不同电阻率泥浆情况下的所有计算，使用的CPU 时间为6478.0+6597.6s.

4 结论

本文用模式匹配算法建立了水平层状非均质各向异性地层中阵列多分量感应测井响应快速算法.利用电导率导数在柱状界面的奇异性，得到两个附加奇异微分算子，能够有效描述积累面电荷对电磁场的影响，使模式匹配算法能够快速准确地模拟非均质各向异性地层中阵列多分量感应的测井响应.

数值计算结果显示，源距､工作频率､地层厚度､侵入半径以及泥浆电阻率等变化，均对共面线圈系的主响应σ_A^XX 产生明显影响.由于共面线圈系的空间响应函数与共轴线圈系的空间响应函数存在非常大的差别，共面线圈系主响应对泥浆电阻率的变化非常敏感，低阻泥浆情况下的视电导率σ_A^XX 变化特征与高阻泥浆对应的视电导率σ_A^XX 存在着很大差异；低阻泥浆视电导率明显小于高阻泥浆的视电导率；由于积累面电荷影响，在低阻层界面附件往往出现大量的负响应，而高阻层中的负响应明显较少；在不同侵入深度的高阻层上，由于共面线圈系的视电导率非常小，从曲线形态差异和视电导率大小难以判断侵入深度变化对于测井响应的影响.与共面线圈系响应相比，共轴线圈系响应σ_A^ZZ 要简单得多，从测井曲线的变化形态较易识别各个地层上电导率的相对大小，泥浆电阻率变化对共轴线圈系响应有一定影响，但没有共面线圈系明显.

附录

将l阶谐变电场强度的旋度公式Δ×E_l=代入Maxwell方程中得

(A1)

从(A1)中提取出柱坐标系中(ρ，θ)两个水平方向对应的偏微分方程:

(A2)

其中(A2)右端的下标“s"表示旋度的水平分量.进一步利用0得

(A3)

将(A3)代人(A2)，并经整理得到方程(6)和(7).