0 引言

电阻率测井中, 砂泥岩薄交互层、部分裂缝层等储层通常被视为一维层状横向同性(TI)地层, 其电导率由横向分量(与层界面平行)和垂直分量(与层界面垂直)表征.由于TI模型突出了砂泥岩薄交互层电导率的主要特征, 且模型较为简单、其电磁响应也容易进行数值模拟, 因而得到了极其广泛的研究与应用(沈金松, 2004;孙向阳等, 2008; Zhong et al., 2008;陈桂波等, 2009).然而对于许多更为复杂环境中形成的沉积地层(例如流水沉积或者风蚀沉积), 这时在沉积地层中会出现交错层, 其地层横向电导率平面(即层理面)与地层界面不再平行, 该类地层在宏观上将会展现很强的各向异性特征, TI模型将不再适用.研究表明, 在感应测井的资料解释方面, 利用各向异性信息可提高储层评价和测井资料的解释精度(汪宏年等, 2008;刘云鹤等, 2018).为了满足复杂地层的评价需要, 在常规TI模型的基础上引入层理面相对于地层坐标系的方位角和倾角(即层理方位角和层理倾角), 可建立适用性更广的一维交错地层模型(Anderson et al., 2001; Wang et al., 2017).传统的感应测井仪器只能提供共轴方向的视电导率, 而多分量感应测井(MCIL)仪器能同时提供地层横纵向电导率、地层边界以及井眼倾角等信息(Barber et al., 2004; Wang, 2011), 将对各向异性地层评价发挥积极作用.此外, MCIL响应的交叉分量可以在相对较短的源距获得更远的探测距离(Hartmann et al., 2014), 更有利于复杂地层中测井资料的解释和评价.由三个正交发射线圈和三个正交接收线圈组成的MCIL仪器能够同时提供磁场张量的九个分量, 其资料处理与解释方法获得了较为广泛的研究(Wang et al., 2012; Yang et al., 2014;朱天竹等, 2017;白彦等, 2018).

目前, 电磁场数值模拟方法主要包括有限差分法(FD)、有限体积法(FV)、有限元法(FE)和数值模式匹配法(NMM)等数值方法(Davydycheva et al., 2003;张烨等, 2012;孙向阳等, 2008;林蔺等, 2017; Wang et al., 2020), 以及主要针对特定模型的解析法(Wang et al., 2008).从理论上说，三维数值模拟技术能够求解任意复杂地层中的电磁场, 应用范围更广, 但其耗时往往很长且对计算机的内存有更高的要求.而解析法往往速度快、精度高, 是对其他数值方法计算效果进行验证的有力工具.此外, 层状介质中背景场的计算也往往采用一维模型中的解析解(沈金松, 2004; Streich, 2009).因此, 一维层状模型的解析解在理论研究和实际应用上都具有重要意义.Michalski和Mosig(1997)利用传输线法(TLM)推导了水平层状TI介质中电流源和磁流源的并矢格林函数, 陈桂波等(2009)在此基础上实现了格林函数的快速计算; Yuan等(2010)考察了无限大双轴各向异性地层中的测井响应; Løseth和Ursin(2007)及罗鸣和李予国(2015)实现了层状任意各向异性模型的电磁场模拟并用于海洋可控源正演; 姚东华等(2010)和He等(2016)模拟了层状双轴各向异性地层中的多分量感应测井响应, 但两者都没有考虑地层坐标系与电导率主轴坐标系不一致的情形.

本文中, 我们应用传播矩阵法(PMM)在姚东华等(2010)的基础上给出了适用性更广的一维层状交错地层中电磁场的解析公式.编写了Fortran计算机程序, 并利用数值模拟结果对多个不同模型上MCIL仪器的测井响应特征、电导率的层理方位角和倾角以及井眼倾角变化对测井响应的影响进行了系统考察.

1 方法与理论 1.1 交错地层模型及Maxwell方程

图 1给出了交错层组成的水平层状地层模型、倾斜井眼轨迹以及由三个相互正交的发射线圈(T_x, T_y和T_z)与三个相互正交的接收线圈(R_x, R_y和R_z)组成的MCIL仪器结构示意图.在图 1a中还给出了地层坐标系xyz: x轴与y轴平行于地层界面, z轴垂直于地层界面向下.地层模型的层数为N, 层界面深度记为z_n, n=2, 3, …, N, 第n层中的电导率张量用横向电导率σ′_{h, n}和垂直电导率σ′_{v, n}以及层理方位角α_n和倾角β_n表示.图 2给出了模型第n层中电导率主轴坐标系x′_ny′_nz′_n与地层坐标系的相对关系:x′_ny′_n平面表示层理面, 此平面内任意方向的电导率都为σ′_{h, n}, z′_n轴垂直于此平面, 沿z′_n轴的电导率为σ′_{v, n}.将x′_ny′_n平面和xy平面的交线确定为y′_n轴, 则可以把y轴与y′_n轴的夹角定义为层理方位角α_n(0≤α_n＜2π), z轴与z′_n轴的夹角定义为层理倾角β_n(－π/2≤β_n≤π/2).进一步利用xyz与x′_ny′_nz′_n之间的复合旋转矩阵：

(1)

可以得到第n层中电导率张量在地层坐标系的表达式：

(2)

图 3中展示了几个具体模型的层理面示意图和对应的层理方位角和倾角：(a)α=0°, β=0°; (b)α=0°, β=45°; (c)α=45°, β=45°; (d)α=90°, β=45°.其中(a)是一个传统TI模型, (b)、(c)和(d)是交错模型, 可以看作是模型(a)分别绕α=0°、α=45°和α=90°的交线轴旋转45°而来.

在地层坐标系xyz下, 位于r_s=(x_s, y_s, z_s)处的磁偶极子发射源M=(M_x, M_y, M_z)激发的电磁场E和H满足频率域中的Maxwell方程(时间变化因子为表示角频率):

(3)

其中, μ₀是真空中的磁导率, 表示复电导率张量, ε₀是真空介电常数, 表示单位矩阵.

1.2 上下模式波的分解及其基础解

为了求解方程(3), 首先对电磁场E(r)和H(r)进行二维Fourier变换:

(4)

并引入波数域下的旋度算子:

(5)

将式(4)和(5)代入方程(3)中, 同时消去和分量, 经整理得到如下常微分方程(Løseth and Ursin, 2007;姚东华等, 2010):

(6)

其中是由电磁场水平分量组成的四阶列向量. 是4×4的系数矩阵, 可以分块化为四个子矩阵, 具体表达式为

和

其中p_x=k_x/ω和p_y=k_y/ω表示水平慢度.源向量S是一个四阶列向量并可分解为, 其中.一旦得到了电磁场的水平分量, 则电磁场垂直分量和可通过下式获得:

(7)

根据矩阵理论, 通过对矩阵的相似对角化可以实现方程组(6)的解耦, 得到新的对角化微分方程组(Løseth and Ursin, 2007;姚东华等, 2010):

(8)

其中对角矩阵由的四个特征值组成, 可以表示为分块对角矩阵, 这里的和分别由正虚部和负虚部特征值组成.而向量w(z)= 表示模式波, 表示新的源向量, 这里是的4×4阶特征矩阵.根据矩阵的性质, 其特征矩阵以及逆矩阵可表示为如下分块矩阵的形式(Løseth and Ursin, 2007):

(9)

如果令, 其中u和d分别代表上行波和下行波, 则方程(8)可以进一步分解成关于上行波和下行波的两个方程:

(10)

和

(11)

对于无限大均匀介质, 由于不存在界面的反射, 在源的上方区域只存在上行波, 源的下方区域只存在下行波.分别在对应区域求解(10)和(11), 可以得到无限大均匀介质中上行波和下行波的基础解(姚东华等, 2010):

(12)

1.3 在1D层状地层中的解

对于水平层状地层模型(图 1a), 第n层的系数矩阵为, 相应的特征值矩阵和特征向量矩阵分别为和.本节首先给出模式波在含源层的形式解, 再给出其余的不含源的各层的形式解并进一步推导相邻层模式波振幅的递推公式, 据此可求得任意层模式波的解.

1.3.1 含源层的解

当源位于第n(1 < n < N)层时, 即z_n < z_s < z_n+1, 该层中的模式波可以看作是由源直接激发的入射波和上下边界面反射波的叠加.利用式(12)的基础解, 在源上方和源下方区域的模式波可分别表示成:

(13)

和

(14)

其中U_n和D_n是二阶向量, 分别表示从下边界z_n+1反射而来的上行波的振幅和从上边界z_n反射而来的下行波的振幅.根据波传播理论, 可以利用广义反射系数矩阵将U_n和D_n分别表示为(Chew, 1995;姚东华等, 2010):

(15)

和

(16)

其中和分别是边界面z_n+1和z_n处的广义反射系数矩阵(见附录A).利用(15)和(16)式求解U_n和D_n得到

(17)

和

(18)

其中和是2×2传播矩阵.将式(17)和(18)代回(13)和(14)式经化简可得模式波的表达式:

(19)

和

(20)

其中A_n^－和A_n⁺分别表示第n层的上行波和下行波的振幅, 其具体表达式为

(21)

当发射源位于顶层或者底层时, 通过类似推导也可得到模式波的解, 见附录B.

1.3.2 其他不含源地层中的解及递推公式

类似于Chew(1995)给出的结果, 不含源层的模式波也可表示为式(19)和(20)的形式, 而电磁场水平分量在层界面z_n处的连续性条件将给出局部反射系数矩阵和透射系数矩阵 (见附录A).在此基础上, 本节给出相邻层之间模式波振幅递推公式.

对于z_s以下的各层, 在边界z=z_n+1处, 第n+1层的下行波等于该层的上行波在界面处的反射波与第n层下行波的透射波的叠加, 利用(19)和(20)给出的模式波的形式解可以得到如下方程:

(22)

类似地, 第n层的上行波等于该层下行波的反射加上n+1层上行波的透射, 由此得到另一个方程:

(23)

联立(22)和(23), 求解得到在源的下方各层中振幅的递推公式:

(24)

以及广义反射系数矩阵的递推公式:

(25)

并令.

类似地, 可以推导出在源z_s上方各层中振幅的递推公式:

(26)

以及广义反射系数的递推公式：

(27)

并令.

1.4 空间域中的电磁场及感应测井视电导率

利用上面求得的模式波w(z)的解并结合公式, 得到电磁场水平分量表达式:

(28)

进一步应用式(7), 可以得到垂直分量和的解.根据上面的推导结果, 可以分别给出三个坐标轴方向的单位磁偶极子产生的磁场, 对其进行组合得到波数域下的磁场并矢格林函数:

(29)

其中表示p方向磁偶极子产生的磁场强度在q方向的分量.再利用Fourier逆变换可以得到空间域中的磁场并矢格林函数

(30)

为了计算(30)右端的积分, 首先用足够大的正方形有限区域[－k_max, k_max]×[－k_max, k_max]代替原积分的二维无穷积分区域, 然后借助二维Gauss-Legendre积分公式计算积分值.一维n阶Gauss-Legendre积分公式为

(31)

其中x_i和w_i分别表示积分节点和权重, 而对于定义在正方形区间[－k_max, k_max]×[－k_max, k_max]上的函数f(x, y), 通过逐次对x和y积分并利用(31)可以得到二维Gauss-Legendre积分公式：

(32)

其中x_i和w_i分别表示x方向的积分节点和权重, y_j和v_j则对应于y方向的积分节点和权重.

在倾角为θ的倾斜井眼中, 可以通过坐标转化将地层坐标系中的磁场并矢Green函数转换到仪器坐标系中:

(33)

其中, 为地层坐标系与仪器坐标系间的旋转矩阵:

(34)

最后, 根据仪器坐标系中磁场强度的各分量H_qp^Tool, 计算出视电导率的各个分量(Zhdanov et al., 2001)：

(35)

其中仪器常数各分量分别为K_xx=K_yy=K_xy=K_yx=ωμ₀/8πL, K_zz=ωμ₀/4πL, K_xz=K_zx=K_yz=K_zy=ωμ₀/16πL, 这里L表示仪器源距.

2 数值结果

本节首先利用TLM(陈桂波等, 2009)和三维有限体积法(3DFV)(张烨等, 2012)与本文的PMM算法的数值结果进行对比, 以验证PMM算法的有效性.之后利用PMM算法模拟多层模型和无限大均匀模型中MCIL的测井响应, 用于系统考察存在交错层理的情况下层理方位角和倾角以及井眼倾角等参数变化对MCIL响应的影响.为简单起见, 固定MCIL仪器的源距为L=1.016 m, 工作频率为f=20 kHz.

2.1 算法验证

图 4同时给出了三层地层模型中由TLM和PMM计算得到的常规TI模型上的MCIL响应以及用3DFV和PMM计算得到的含交错层理的模型上MCIL响应间的对比, 简单起见只考虑直井的情形.该模型的顶层和底层是电导率相同且不含交错层理的TI层, 其横向和垂直电导率分别是0.5 S·m^－1和0.2 S·m^－1, 中间层是厚度为3m的交错层, 横向和垂直电导率分别是0.1 S·m^－1和0.05 S·m^－1, 而层理方位角与和倾角分别选取两组不同的值α, β=0°和α=0°, β=60°, 分别对应于常规TI模型和交错模型.图 4中包含MCIL响应的三个主分量(σ_{a, xx}, σ_{a, yy}和σ_{a, zz})和两个交叉分量(σ_{a, xz}和σ_{a, zx}), 其中由3DFV和PMM计算得到的交错地层模型(α=0°, β=60°)的响应分别用离散的三角符号和虚线表示, 而TLM和PMM计算得到的TI地层模型(α, β=0°)响应分别用离散的星型和实线表示, 结果显示在TI模型中PMM与TLM的计算结果完全重合, 在交错模型下PMM与3DFV结果的主分量也符合较好, 并且交叉分量除了界面附近外也基本重合, 证明了PMM算法的有效性.此外, 数值结果还显示TI模型中交叉分量σ_{a, xz}和σ_{a, zx}均等于零, 而在由PMM与3DFV计算得到的交错模型中间层上的交叉分量响应非常明显, 产生该现象的原因是直井情况下TI模型的电导率相对于仪器轴的空间分布是完全对称的, 而在α=0°, β=60°的交错地层上, 电导率分布不再是对称的.同时还可以看到, 共面主分量σ_{a, xx}和σ_{a, yy}在界面附近出现了负响应, 这是因为地层电导率的不连续性会导致界面处的面电荷积累, 进而影响场的分布.

2.2 交错层的层理方位角α=0°时倾角β变化对响应的影响

本节将考察在交错层的层理方位角α=0°的情况下, 由于层理倾角β与井眼倾角θ变化对多分量测井响应的影响.为此, 设计出由各向同性地层与交错地层交互组合成的九层模型(见图 5a), 中间七个地层的厚度均为3 m, 上下围岩与中间三个各向同性地层电导率均相等并假定为0.1 S·m^－1, 而四个交错地层具有相同的横向电导率σ′_h=5 S·m^－1和垂直电导率σ′_v=0.2 S·m^－1, 从上到下各交错层的倾角依次为β=0°, 15°, 30°, 60°.同时还假定有三个不同倾角θ=30°, 0°, －30°的井眼, 图 5(b—f)同时给出了三个井眼中的MCIL响应的三个主分量(σ_{a, xx}, σ_{a, yy}和σ_{a, zz})和两个交叉分量(σ_{a, xz}和σ_{a, zx})的正演结果, 其他四个分量σ_{a, xy}, σ_{a, yx}, σ_{a, yz}, σ_{a, zy}都等于零, 没有展示其结果.

首先分析MCIL在直井中的响应, 从图可以看出在四个不同层理倾角的交错地层中, 响应结果差异十分明显, 在β=0°的常规TI层中, 交叉分量σ_{a, xz}和σ_{a, zx}均等于零, 且共面主分量σ_{a, xx}和σ_{a, yy}完全相同, 而随着从上到下的各交错层的层理倾角β不断增大, 相应的共面分量σ_{a, xx}和σ_{a, yy}会单调增大, 而共轴分量σ_{a, zz}单调减小.在θ=30°的斜井情况下, 层理倾角β=0°的常规TI层中的交叉分量σ_{a, xz}和σ_{a, zx}不再为零, 并且在四个交错层中, 三个主分量都大致在β=30°的交错层中取得极值, 而交叉分量的绝对值在该层取得极小值, 产生该现象的主要原因是该层的层理面正好与井轴垂直.而在θ=－30°的斜井情况下, 在β=60°的交错层三个主分量取得极值, 而交叉分量的绝对值取得极小值, 而这种情况正好是层理面与井轴平行的情形.

为分析总结不同层理倾角和井眼倾角对MCIL响应影响的变化特征, 这里引入θ_rel=|θ－β|用于描述井轴与层理面之间的相对倾角, θ_rel的大小直接影响了交错层上MCIL响应.例如, 在层理面倾角β=15°的交错层中, 直井和θ=30°的斜井两种情况对应的θ_rel均为15°(见图 6), 图 5中的结果清楚显示了这两种情况相对应的主分量响应值几乎完全重合.类似地, θ=0°的直井在β=0°的TI层的响应值与θ=30°的斜井在β=30°交错层中的主分量响应也几乎完全相同.此外, 交叉分量响应的绝对值取得极小值时对应着井轴和层理面垂直或者平行的情形, 即θ_rel=0°或者θ_rel=90°.需要指出的是, 这些规律只适用于在交错层中间部分的响应值, 在层边界附近上, 由于受到围岩影响, 共面分量(σ_{a, xx}和σ_{a, yy})以及交叉分量(σ_{a, xz}和σ_{a, zx})都容易出现“犄角”, 而各种情况下σ_{a, zz}分量的曲线都较为光滑.

为了进一步考察层理倾角对磁场空间分布的影响, 考虑中间层为交错层的三层模型.该模型的上下围岩是电导率σ=0.1 S·m^－1的各向同性地层, 而中间的交错层厚度为5 m, 横向和垂直电导率分别为5 S·m^－1和0.2 S·m^－1.发射源位于中间层的中间处, 图 7给出了交错层的层理倾角分别为β=0°, 45°(层理面在图中用灰线表示), 发射源(以T_z为例)倾角分别为θ=0°, 45°共计四种组合的磁场强度虚部的空间分布图.对比图 7a和图 7d不难发现, 尽管发射源的方向不同, 但由于发射源和层理面的相对倾角θ_rel都等于0°, 源附近的磁场空间分布几乎完全相同.在图 7b和图 7c中也可以看到源附近相似的磁场分布, 因为相对倾角都为θ_rel=45°.图中同时给出了沿发射源方向的直线上数个观察点处的磁场方向, 可以看到图 7a中各磁场方向与发射源严格地同向, 这说明直井情况下TI模型的交叉分量为零, 而图 7d中在层中间附近的各磁场方向与发射源方向近似同向, 说明θ_rel=0°时交错模型的各交叉分量响应的绝对值较小, 而在边界附近由于围岩的响应较大, 磁场方向偏离发射源方向的程度较大.

2.3 层理方位角α和倾角β同时变化对响应的影响

下面将综合考查层理方位角α和倾角β变化对MCIL响应的影响, 首先采用与图 7相似的三层模型, 上下围岩是电导率为0.1 S·m^－1的各向同性地层, 中间层为厚度为5 m交错层, 其横向和垂直电导率分别为5 S·m^－1和0.2 S·m^－1, 层理倾角固定为β=45°不变, 而层理方位角α分别取0°, 45°和90°三个不同的值(见图 8j).

图 8(a—i)是直井情况下由三个不同层理方位角模型对应的MCIL的九个分量的响应, 由图可看出:当层理面与仪器的T_y轴平行时(α=0°), 与y相关的交叉分量σ_{a, xy}, σ_{a, yx}, σ_{a, yz}, σ_{a, zy}都等于零, 层理面与仪器的T_x轴平行时(α=90°), 与x相关的交叉分量σ_{a, xy}, σ_{a, yx}, σ_{a, xz}, σ_{a, zx}都等于零, 而当层理面与仪器T_x或T_y轴都不平行时(α=45°)视电导率的所有分量都不等于零.利用(2)式我们也可以求得这三个模型的中间层在地层坐标系中的电导率张量, 分别为

容易看出视电导率各分量是否为零与中间层的电导率各分量有简单的对应关系.此外, α=45°和α=90°模型的层理面分别可看作α=0°模型的层理面绕z轴旋转45°和90°而来, 因此σ_{a, zz}响应保持一致.

为进一步考察井眼倾角对测井响应的影响, 图 9是θ=30°的斜井中正演模拟结果, 其地层模型参数与图 8完全相同.由于井眼倾斜, MCIL仪器与层理面的相对位置关系相对图 8发生了变化(图 9j), 响应结果的各个分量发生了明显变化.但α=0°模型的层理面仍然和仪器的T_y轴平行, 其结果也与图 8中α=0°的结果相似:与y相关的各交叉分量σ_{a, xy}, σ_{a, yx}, σ_{a, yz}, σ_{a, zy}都等于零.而α=45°和α=90°的层理面与仪器的T_x和T_y轴都不再平行, 因此视电导率的所有分量都不为零.由于斜井的仪器坐标系与地层坐标不再重合, 因此类似图 8中视电导率响应与中间层在地层坐标系下电导率的简单对应关系也不再成立.

前面的几个正演模拟结果均是在固定层理方位角α或者倾角β的情况下得到的.为了考察α和β同时变化对MCIL响应的影响, 考虑一个横向和垂直电导率分别为5 S·m^－1和0.2 S·m^－1的无限大交错地层, 只考虑直井的情形.图 10在极坐标上展示了视电导率的九个分量随α和β连续变化的模拟结果, 极坐标的极角和极径分别代表α和β.由图可看出主分量σ_{a, xx}和σ_{a, yy}在方位角α方向相差90°, 此外完全一致, 在|β| < 45°时对α的变化不敏感.而主分量σ_{a, zz}仅与β相关, β为0°和90°时分别取得极大值与极小值.对交叉分量而言, σ_{a, xz}, σ_{a, zx}, σ_{a, yz}, σ_{a, zy}表现出相似的特征:在整个αβ平面内取得两个极大值和两极小值, 取极值时|β|≈60°, 而σ_{a, xy}, σ_{a, yx}与其他交叉分量显著不同, 整个αβ平面内取得四个极大值和四个极小值, 取极值时|β|=90°.

4 结论

本文应用PMM算法模拟了一维层状交错地层模型中任意朝向的磁偶极子激发的电磁场, 并通过与TLM和3DFV模拟结果对比验证了该算法的有效性.之后利用多个具体模型考察了井眼倾角和交错层的层理方位角α以及倾角β对MCIL响应的影响, 并给出常规TI地层模型和交错地层模型中磁场的空间分布图.模拟结果显示当层理方位角α=0°时, 井轴与层理面的相对倾角θ_rel很大程度上决定了响应结果.对于更一般的情形, α和β以及井眼倾角θ都对MCIL响应有显著的影响.此外, 相较TI模型, 交错地层中的视电导率响应结果更为复杂, 往往会出现更多的非零交叉分量, 这增大了资料处理和解释的难度, 基于交错模型的反演算法将是日后有待研究的课题.最后需要指出的是对于水平井或者井眼倾角接近90°的大斜度井, 由于被积函数的高度振荡性本文所使用Gauss-Legendre积分方法将不再适用, 需要适应性更强的积分技巧.

附录A 透射系数矩阵和局部反射系数矩阵

利用电磁场水平分量在界面z=z_n(n=2, 3, …, N)处的连续性条件可以确定透射系数矩阵 和局部反射系数矩阵假设一个振幅为A_n－1⁺的下行入射波从第n－1层透过界面z_n射入第n层.不考虑其他界面的影响, 则界面z_n以上及以下区域的模式波可以表示为

(A1)

和

(A2)

考虑到电磁场水平分量在界面z=z_n处的连续性条件b_n－1(z_n^－)=b_n(z_n⁺), 这里z_n^－和z_n⁺分别表示z_n的左右极限, 再利用公式和式(A1)、(A2), 可以得到:

(A3)

再将式(9)的结果代入式(A3), 求解和得到如下表达式:

(A4)

同理, 假设一个上行入射波从第n层透过界面z_n射入第n－1层, 经过类似的推导最终可得到:

(A5)

附录B 顶层或底层含发射源时的模式波

当发射源在顶层时(n=1), 不存在来自于上界面的反射, 该层的模式波可以表示为

(B1)

和

(B2)

其中表示此层的下行波的振幅.

类似地, 当发射源在底层时(n=N), 不存在来自于下界面的反射, 该层的模式波可以表示为

(B3)

和

(B4)

其中表示该层上行波的振幅.