1 引言

地层内部由于受非均匀应力、热对流、温差以及物质迁移等多种因素影响，往往会引起地层电性、磁性以及热力学性质的各向异性.在石油测井中，电性各向异性的探测与研究对储层解释和评价有着十分重要的意义.当前各种常规电测井仪器只能测量到某一方向的电阻率分量，这给各向异性储层的正确解释和评价带来了困难^[1].多分量感应测井（MCIL）仪器利用三个相互正交的发射线圈（T_x、T_y和T_z）和三个相互正交的接收线圈（R_x、R_y和R_z）同时测量九个磁场分量（H_xx、H_yy、H_zz、H_xy、H_yx、H_xz、H_zx、H_yz和H_zy），其测量结果正好构成一个完整的并矢Green函数，这将更加有利于解释和评价各向异性地层中的测井资料^[2~8].然而在整个地球物理领域，目前对层状各向异性问题的研究主要采用如下两种模型：具有垂直对称轴的横向同性（TIV）模型^[9~14]和具有水平对称轴的横向同性（又称方位各向异性Azimuthal Anisotropy）模型^[15~20].这两种模型都只有一个旋转对称轴，人们利用传输线理论（TLM）、Hertz势理论以及模式匹配算法（NMM）等已经对横向同性模型中的多分量感应测井响应进行了深入研究^[6~8].但是对于更为一般的完全各向异性地层中的测井响应，现有的报道却非常少.为充分利用多分量感应测井仪器提供的丰富电磁场资料解决更为复杂的地质问题，本文以具有双轴对称性的正交各向异性（Orthogonal Anisotropy）^{[21, 22]}介质为基本模型，通过传播矩阵算法研究并建立了层状正交各向异性地层中多分量感应测井响应的数值模拟算法.最后通过数值模拟结果分别考察了各向异性系数、仪器长度、发射源频率以及井眼倾角变化时，多分量感应测井响应的特征和变化规律.

传播矩阵法作为一种处理层状介质中波传播问题的有效手段，已被广泛应用在光学、等离子物理学以及地球物理学等领域^[23~27].目前该算法的推导过程较为繁琐，不利于理解和掌握.本文通过以下方法对其推导过程进行了改进：首先将均匀介质中上行波和下行波模式的求解表示成一阶常微分方程的边值问题，使电磁波求解更加容易；然后在求解层状地层中电磁波模式的解析表达式时，利用叠加原理和边界条件给出了更为简洁的广义反射系数和电磁波振幅的递推公式，避免了由多个矩阵直接相乘导致的大指数项问题，整个算法更容易在计算机上实现.

2 理论与方法 2.1 Maxwell方程的模式分解

假定地层模型由N个水平层状正交各向异性地层组成（图 1），各个地层的界面位置用向量d_n（n=2，3，…，N）表示，电导率用对角张量σ=diag（σ_xn，σ_yn，σ_zn）（n=1，2，…，N）表示，其中σ_xn，σ_yn和σ_zn分别表示第n个地层三个主轴方向上的电导率分量.用M=表示位于r_s处三个相互正交的磁偶极子，E和H分别对应电、磁场强度，磁导率μ为常数，电磁场的时间变化因子为exp（-iωt），ω为发射源角频率.在忽略位移电流的情况下，将频率空间域中的Maxwell方程组转换到频率波数域：

(1)

由于地层电导率函数σ与坐标变量x和y无关，对方程（1）中的电磁场关于变量x和y进行Fourier变换

由直角坐标系下旋度算子“ Δ × ”的张量形式和傅氏变换微分定理，频率波数域中的旋度算子可表示为

将上式代入方程（1）中，引入水平慢度变量和并进行一系列处理后，可得频率波数域中的电磁场水平分量满足如下方程：

(2)

其中是由频率波数域中的电磁场水平分量组成的4阶列矢量；S=iωμ 为源矢量；I为4×4阶单位矩阵，上角标“T”表示转置；为4×4阶系数矩阵，可以表示成4个2×2阶的分块矩阵：

而电磁场的垂直分量与水平分量满足如下关系

(3)

方程（2）是关于矢量V（z）的一阶线性微分方程组，其解可分解成上行波和下行波模式的组合.首先利用本征值理论将矩阵A表示为

(4)

其中是由矩阵A的四个本征值组成的对角矩阵，Λ=diag（λ₁，λ₂），L是由矩阵A的本征向量组成的4×4阶矩阵，可分解成4个2×2阶的分块矩阵^{[28, 29]}：

(5)

其逆矩阵

(6)

本征值和本征向量的具体求解见附录A.

将（4）式代入方程（2）中，令电磁波模式W（z）=，发射源=，经整理后可得上行波u和下行波d分别满足下面的方程：

(7)

其中

(8)

2.2 电磁波在均匀正交各向异性介质中的模式

在均匀正交各向异性介质中，由于不存在反射界面，方程（7）中由Σ_s^-和Σ_s⁺产生的上行波u和下行波d分别满足条件和，这里的z_s^-和z_s⁺分别表示z_s的左右极限.进一步利用积分公式和-z_s）dz=1，可将非齐次方程（7）转化为如下的齐次方程定解问题：

(9)

求解方程（9），可得均匀正交各向异性介质中上行波和下行波的解析表达式：

(10)

2.3 电磁波在层状正交各向异性地层中的模式

对于层状正交各向异性地层，用表示第n个地层的系数矩阵，其对应的本征向量和本征值矩阵分别由和表示.设发射源Σ_s位于中间的第k层，由叠加原理可知，该层中的电磁波是由发射源Σ_s产生的入射波与上、下层界面上的反射波的叠加.即在区域z_k < z < z_s中，电磁波模式可表示为

(11)

而在区域z_s < z < z_k+1中，电磁波模式的表达式为

(12)

其中A和B均为2阶列向量，分别表示下界面和上界面上反射波的振幅.

利用波场传播理论^[30]可知，A等于入射到界面z=z_k+1上的整个下行波的反射

(13)

而B等于入射到界面z=z_k上的整个上行波的反射

(14)

其中和分别表示界面z_k和z_k+1的广义反射系数，它们均为2×2阶矩阵.

求解方程（13）和（14）得

(15)

(16)

其中

(17)

将（15）~（17）式代入（11）和（12）式，经整理可得发射源所在层的电磁波模式表达式

(18)

其中为2×2阶的单位矩阵，

(19)

均为2阶列向量，表示该层中电磁波的振幅.

当源Σ_s位于顶层时，由于没有上界面，（11）和（12）式中的B=0.介质1中的电磁波模式可简化为

(20)

其中.

当源Σ_s位于底层时，只需令（11）和（12）式中的A=0，即可得到类似于上面的电磁波模式简化表达式.

对于发射源Σ_s以外的其他地层，电磁波模式也具有类似（18）式的解析形式

(21)

其中A_n⁺和A_n^-分别表示位于发射源所在层（k）下部和上部的各个地层中的电磁波振幅.当n=1时，假设；当n=N时，假设.

2.4 电磁波振幅和广义反射系数的递推公式

下面进一步确定（18）~（21）式中的电磁波振幅以及广义反射系数.首先考虑n≥k时，由（21）式可将第n+1层中的电磁波模式表示为

显然，在界面z=z_n+1上，第n+1层中的下行波是该层中的上行波的反射与第n层中的下行波的透射的叠加：

同样，第n层中的上行波是该层中的下行波的反射与第n+1层中的上行波的透射的叠加：

求解上面两式可得电磁波振幅A_n+1⁺和广义反射系数的递推公式：

(22)

(23)

其中和以及和分别表示界面z=z_n+1的反射系数和透射系数（具体的计算公式见附录B）.

采用类似的方法可得n≤k时，电磁波振幅A_n-1^-和广义反射系数的递推公式：

(24)

(25)

其中.

2.5 多分量感应测井响应计算

由各个地层中的电磁波模式W（z）解析表达式和方程V（z）=LW（z），可得频率波数域中的电磁场水平分量

(26)

再通过（3）式即可确定频率波数域中的电磁场垂直分量和.进一步取（8）式中的M_x、M_y和M_z分别等于1，将这三个相互正交的磁偶极子产生的磁场组合起来，即构成频率波数域中的磁场并矢Green函数

通过如下Fourier逆变换

(27)

即可得到频率空间域中的磁场并矢Green函数.由于（27）式是二维无限平面上的广义积分，为了保证计算精度和运算效率，采用二维Patterson自适应求积算法结合有限连分式展开技术^[31~33]专门开发了相应的计算软件.

为了考察不同倾斜井眼情况下的多分量感应测井响应，需要将上面地层坐标系XYZ下的计算结果转化到仪器坐标系X'Y'Z'^[2].设井眼（即Z'轴）与Z轴的夹角为α，OX'Z'平面与OXY平面垂直且与OXZ平面的夹角为θ.利用坐标旋转可得仪器坐标系X'Y'Z'与地层坐标系XYZ间的旋转矩阵

这样，仪器坐标系下的磁场并矢Green函数可表示为

(28)

多分量感应测井主要是通过测量二次场产生的感生电动势^[6]来探测地层电导率的空间分布和不同方向上的电导率大小，所以其响应特征主要与磁场并矢Green函数的虚部（）有关.

3 数值结果

本节采用前面提出的PMM算法具体考察正交各向异性地层中多分量感应仪器的测井响应，用和分别表示水平各向异性系数和垂直各向异性系数.当m=1.0时，介质为垂直横向同性；当m=n时，介质为方位各向异性.

3.1 算法验证

为检验PMM算法的有效性和计算精度，在TIV模型上对比给出PMM算法和TLM算法的数值模拟结果.令σ_x=σ_y=σ_h，σ_z=σ_v，地层模型参数见表 1.图 2给出了垂直井眼情况下发射源频率f=20kHz，仪器长度L=1.5 m时，由两种不同算法计算得到的磁场强度主分量曲线（H_xx、H_yy和H_zz）.观察该图可以看到两者吻合的非常好，说明用PMM算法计算测井响应是可行的，并且其计算结果的精度也是较高的.

3.2 各向异性系数变化对测井响应的影响

为考察各向异性对测井响应的影响，建立如下模型，地层模型参数见表 2.图 3给出了垂直井眼情况下，发射源频率f=20kHz，仪器长度L=1.0m，垂直各向异性系数n²=1.0，水平各向异性系数m²=0.25，0.5，1.0，2.0，4.0时，测井响应的计算结果.从图中可以看到：磁场强度y轴响应H_yy不受水平各向异性系数m²变化的影响，而x轴响应H_xx随m²的增大而增大，z轴响应H_zz随m²的增大而减小，并且测井响应在低阻层受水平各向异性系数变化的影响比在高阻层明显；此外，磁场强度的x轴和y轴响应在层界面附近较为复杂，出现“犄角”现象，而z轴响应相对简单，无“犄角”现象出现.“犄角”现象的产生，主要是由于层界面附近积累了大量面电荷导致的，随着水平各向异性系数的减小，x轴响应的“犄角”逐渐消失，说明层界面上的面电荷随地层y轴电导率的增大而减少.同样，图 4给出了水平各向异性系数m²=1.0，垂直各向异性系数n²=0.25，0.5，1.0，2.0，4.0时，测井响应的计算结果.观察图 4发现：磁场强度z轴响应H_zz不受垂直各向异性系数n²变化的影响，图形较为简单，而x轴和y轴响应完全相同，且随垂直各向异性系数的增大而减小，图形较为复杂，出现“犄角”现象.

3.3 仪器长度和发射源频率变化对测井响应的影响

采用与3.2节相同的地层模型，将电阻率地层的水平各向异性系数（m²）和垂直各向异性系数（n²）分别取为0.5和2.0.图 5给出了垂直井眼情况下发射源频率f=20kHz，仪器长度L=1.0，2.0，3.0m时，测井响应的计算结果.可以看到：随着仪器长度的增大，响应曲线的纵向分辨率逐渐降低.图 6又给出了垂直井眼情况下，仪器长度L=1.0 m，发射源频率f=20，30，50kHz时，测井响应的计算结果.观察发现：发射源频率的增大将导致测井响应的增大，这是由于趋肤效应增强造成的.

3.4 井眼倾角变化对测井响应的影响

仍采用3.3小节的地层模型，图 7给出了发射源频率f=20kHz，仪器长度L=1.0 m，井眼倾角α=0°，30°，60°时，测井响应的计算结果.从图中可以明显看到：当井眼倾角不为零时，垂直磁偶极子M_z在水平x方向以及水平磁偶极子M_x在垂直z方向上均产生负响应，即磁场强度分量H_xz和H_zx不再为零；响应曲线除H_zz外，H_xx、H_yy、H_xz和H_zx均在层边界附近出现“犄角”现象，且“犄角”现象在高阻层比在低阻层更明显；随着井眼倾角的增大，磁场强度分量H_xx和H_yy不断增大，H_zz、H_xz和H_zx不断减小.以上结果表明：在井眼倾斜时，测井响应是地层三个主轴方向电导率的综合反映.这就增大了实际资料处理和解释的难度，要想从测井曲线中提取出三个主轴方向上的电导率，需要借助特殊的反演手段.

4 结论

本文通过传播矩阵理论推导出层状正交各向异性地层中多分量感应测井响应的正演模拟公式.将均匀介质中的电磁波求解问题表示成一阶常微分方程的边值问题，使电磁波模式的求解更加容易；由递推公式给出多层介质的广义反射系数和电磁波振幅，避免了计算多个矩阵的直接相乘；采用二维Patterson自适应求积算法结合有限连分式展开技术计算傅氏逆变换，提高了运算效率和计算精度.通过考察不同各向异性系数、不同井眼倾角以及仪器长度和工作频率变化时的多分量感应测井响应特征发现：水平各向异性系数的变化不改变y轴的测井响应，垂直各向异性系数的变化不改变z轴的测井响应；仪器长度的增大导致测井响应纵向分辨率的降低；发射源频率的增大又会造成趋肤效应的增强；井眼倾角的变化不但影响测井响应的三个主轴分量，还将产生x轴和z轴方向的负响应.最后需要指出的是，本文只考虑了分层均匀介质的正演模拟，对于含有井眼和侵入带的测井响应正演，以及提取地层真实电导率的反演技术，将是日后有待研究的新课题.

附录A本征值及本征向量的求解

矩阵A的本征值由如下方程确定

(A1)

这涉及到求解关于本征值λ的一元四次方程^[34].对于正交各向异性介质，由于矩阵A是斜对角的分块矩阵，其本征值具有如下的简单解析表达式

(A2)

(A3)

(A4)

其中

选取归一化的本征向量矩阵^{[28, 29]}

(A5)

将（A2~A5）式代入方程（4），经一系列运算得

其中

附录B单层界面的反射系数和透射系数

对于只有一个界面的2层模型，设层界面位置z=z₁，与介质1和介质2对应的系数矩阵的本征值和归一化本征向量矩阵分别为和.首先假设在介质1中有一向下入射的单位电磁波，遇到界面z=z₁后发生反射和透射.此时介质1和介质2中的电磁波模式可表示为

(B1)

(B2)

其中R_1，2和T_1，2分别表示界面z=z₁的反射系数和透射系数，它们都是2×2阶矩阵.

利用电磁场水平分量V（z）=LW（z）在界面z=z₁上连续，结合方程（B1）和（B2）可得反射系数R_1，2和透射系数T_1，2满足下面的方程

(B3)

求解上式得

(B4)

同理，假设在介质2中有一个向上入射的单位电磁波I'，遇到界面z=z₁后发生反射和透射，可得另一组反射系数R_2，1和透射系数T_2，1的表达式

(B5)