0 引言

电测井是评价地层的重要手段, 在复杂油气识别中起着重要作用.当地层存在电性各向异性时, 相比较常规地层, 其电测井响应更为复杂(姚东华等, 2010; 洪德成和杨善德, 2011; Gao et al., 2013; Zhu et al., 2018).目前, 在各向异性地层电测井响应计算与应用研究中, 广泛采用横向各向同性(transverse isotropy, TI)模型, 指地层只存在水平与垂直电导率不同的情况, 也称单轴各向异性(Løseth and Ursin, 2007; Zhong et al., 2008; 邓少贵等, 2006).当地层发育裂缝时(Georgi et al., 2008; Tang et al., 2007), 如图 1a, 或沉积过程中横向上岩性和物性等变化时(Gao et al., 2002; Li et al., 2017; Zhou et al., 2017), 如图 1b, 地层电导率沿x、y、z三个方向均不同, 即显示双轴各向异性(biaxial anisotropy)(Hou et al., 2015; Nekut, 1994; He et al., 2016).对于页岩油气、致密油气等一些复杂油气藏, 可能存在由于沉积或成岩作用引起的层理各向异性, 又存在裂缝发育引起的诱导各向异性.因此, 基于各向同性、甚至单轴各向异性介质模型的电测井理论和应用研究, 难以有效满足这类复杂油气藏的评价需要, 亟待开展基于双轴各向异性介质模型的测井响应和评价研究.

双感应或双侧向一般只提供深、浅两种探测深度, 且只反映了地层的宏观导电特性(邓少贵等, 2018; 张庚骥, 1984; Anderson et al., 1995; Li and Wang, 2016); 阵列感应测井可以提供三种纵向分辨率、五种探测深度的电阻率信息, 但缺少分量测量, 依然难于有效识别薄互层和评价各向异性(Wang et al., 2008; Zhang et al., 2007).为适应各向异性复杂油气藏勘探开发需求, 国内外石油公司陆续推出多分量感应测井仪(王磊等, 2015; Davydycheva et al., 2009), 如Schlumberger公司的Rt-Scanner和Halliburton公司的3D-Induction Xaminer等.多分量感应测井仪可同时测量9个磁场分量, 各分量对倾角、方位角、各向异性等不同地层参数有不同的敏感性, 利用反演处理可获取地层电阻率、各向异性等信息(Wang H M et al., 2006; 汪宏年等, 2008; Hong et al., 2014; Wang L et al., 2018a, b), 但其处理方法主要基于单轴各向异性介质模型.因此, 为了适应各向异性复杂油气藏评价要求, 需要基于双轴各向异性介质模型, 分析多分量感应测井的响应特征, 提升其应用效果和能力.

双轴各向异性介质的电测井响应数值模拟主要采用解析法和有限差分法.Davydycheva和Wang(2011)利用有限差分法研究了双轴各向异性介质多分量感应测井响应特征, 但有限差分算法速度慢(Chew, 1999; Anderson, 1982; 谭茂金等, 2007), 不适用于测井资料快速反演; Nekut(1994)提出了双轴各向异性介质中磁偶极子源时域解析解, 其主要考虑仪器在零源距下的时域信息; Gao等(2002)给出了双轴各向异性介质多分量感应测井频域解析解, 但未考虑仪器旋转角等问题, 而实际中测井仪器可能与地层成任意夹角放置; Yuan等(2010)使用解析法进一步考虑了任意旋转角的问题, 但对于电磁场的高维积分问题, 特别是在倾角较高时, 其积分速度和精度仍有待提高.

本文基于均匀无限厚双轴各向异性介质模型, 推导并简化多分量感应测井空间域解析式, 深入挖掘谱域电磁场的周期特性与振荡特性, 给出一种高效、精确的三重积分方法, 并对比分析不同纵横向各向异性条件多分量感应测井响应特征.

1 双轴各向异性介质多分量感应测井解析解 1.1 双轴各向异性介质多分量感应测井基本理论

空间域中, 时谐源满足e^－iωt时, 麦克斯韦方程组可表示为：

(1)

(2)

其中, μ₀为真空磁导率; 为地层介电常数张量, 频率较低时可忽略其影响为地层电导率张量.如图 2a所示, 双轴各向异性介质电导率张量表示为

(3)

一般而言, 若不考虑线圈的尺寸及结构时, 多分量感应测井仪发射线圈可等效为磁偶极子源.对不同发射方向的磁偶极子源, 其在接收线圈处激发的电场和磁场可分别由并矢矩阵和表示：

(4)

(5)

其中, E_ij、H_ij的下标i和j分别表示发射线圈和接收线圈的法向.

当仪器与地层坐标系z向的相对夹角α不为0°时, 需利用旋转矩阵R将接收线圈测量响应转换至仪器坐标系进行计算：

(6)

(7)

式中, T表示向量或矩阵的转置, R可表示为

(8)

其中β、γ分别为方位角和旋转角, 如图 2b所示.

1.2 双轴各向异性介质多分量感应测井响应解析解

对空间域中电场E(x, y, z)、磁场H(x, y, z)和偶极子源M(x, y, z)等分别沿x、y和z进行三重傅里叶变换, 可得谱域所对应的象函数、：

(9)

将式(2)两端求旋度, 并将式(1)代入, 可得：

(10)

式中.根据式(9), 并由麦克斯韦方程组可得：

(11)

(12)

式(11)中, Ω为系数矩阵, 可表示为

(13)

根据各方向源的不同特征, 可将式(11)展开以获得谱域中电场的解析式, 再由式(12)即可求得相应的磁场表达式.例如, 当源位于r′=(x′, y′, z′)且沿x轴方向发射时(以电场为例, 磁场不做赘述)：

(14)

将其整理后可得：

(15)

对式(15)两端进行三重傅里叶逆变换, 即可得空间域电场的解析式：

(16)

引入, 可将其转换到柱坐标系下：

(17)

2 谱域内电磁场积分简化方法

由式(17)可知, 场的解析式包含三重积分.其核函数在各方向均存在振荡性, 而振荡性函数积分一直是困扰科学与工程计算的难题, 故需对其进行简化以降低积分难度.

2.1 围线积分法

核函数关于无穷积分的振荡性随ζ的增大而增强, 这使得其高效、精确积分非常困难.经研究发现, 可利用围线积分法对其进行降维处理, 即将三重积分转化为两重积分.如式(17)中, 注意到：

(18)

可知区间内存在±ζ_o和±ζ_e四个奇异点.由场的收敛性可知, 当z－z′>0时, 积分只在虚部大于0的极点处(表示为ζ_o^>和ζ_e^>)有值, 如式(19);当z－z′ < 0时, 积分只在虚部小于0的极点处(－ζ_o^>和－ζ_e^>)有值, 如式(20).

(19)

(20)

2.2 核函数周期特性

关于定积分的核函数为三角函数, 利用其周期特性将积分区间压缩以提高计算效率.以沿x方向发射的磁偶极子源为例(其余项详见附录A)：

(21)

(22)

(23)

其中：

(24)

(25)

(26)

(27)

需要注意的是, 几种特殊情况下可进一步简化积分计算.当仪器位于地层坐标系xoz平面内时, 有：

(28)

同理仪器位于地层坐标系yoz平面内时：

(29)

仪器轴与地层坐标系z轴平行时：

(30)

经简化后, 积分区间压缩至原来的四分之一, 相应的计算量将降低至原来的四分之一.

3 核函数的快速计算与误差分析 3.1 核函数的振荡性

与TI介质中场的表达式类似(Howard, 2000; Key, 2012; Johansen and Sørensen, 1979), 核函数均存在振荡性, 且其振荡性与场分量及倾角等有关.设R_x=2 Ωm, R_y=4 Ωm, R_z=8 Ωm, β=γ=0°, 以H_zz实部为例, 不同倾角下双重积分核函数在k-ψ平面上的展布如图 3.其中, 图 3a、图 3d、图 3g、图 3j横坐标表示H_zz核函数实部关于ψ的积分区间, 单位为弧度; 纵坐标表示H_zz核函数实部关于k的积分区间, 单位为1;色标表示H_zz核函数实部值, 单位为1;其余各图可依此类推.核函数沿k方向积分可得磁场与ψ的关系, 如图 3b、图 3e、图 3h和图 3k所示; 核函数沿ψ方向积分可得磁场与k的关系, 如图 3c、图 3f、图 3i和图 3l所示.可知, 当倾角增大时, 双重积分核函数沿ψ和k方向的振荡性均增强, 沿k方向的振荡区间增大.因此, 采用固定的积分区间与高斯节点时, 将难以协调计算精度与速度.

3.2 自适应积分区间截断

由于不同场分量振荡性不同, 对不同核函数在不同倾角下的半无穷积分, 可截断积分区间, 在确保核函数收敛性的同时, 压缩其积分区间与节点, 以保证计算精度并提高计算效率.积分区间的自适应截断是通过预先设定一个截止因子, 当积分计算到几个点的幅值连续小于截止因子时停止积分, 当前区间即为截断区间.不同分量在不同倾角下应有不同的截止因子.例如, 在较大区间内(0~200)找到核函数最大幅值点, 取该幅值万分之一为截止因子, 再以此点为左端点通过二分法进行区间选取.图 4是H_zz分量的相位与幅度在倾角为60°时截止因子取最大幅值不同倍数的相对误差示意图.

3.3 高斯节点数选取

两重积分核函数振荡性不一致, 故对其选择不同的节点数以提高计算效率.定积分可得到精确的结果, 故先计算定积分再计算半无穷积分以提高精度.如图 5是H_xx、H_xz与H_zz分量的相位与幅度在不同倾角下随ψ节点数变化的误差示意图.由图可知, 当节点增多时, 区间内节点数达到饱和, 且倾角越小时所需节点数越少; 同时, 不同场分量饱和节点数不同.

各分量半无穷积分区间不同, 故分别选取不同的节点数.如图 6是H_xx、H_xz与H_zz分量相位与幅度在不同倾角下, 误差随积分节点数变化示意图.由图可知, 不同场分量在不同倾角下饱和节点数不同.当节点数趋于饱和时, 继续减小误差需增大积分区间并增加节点数.若精度要求不高, 区间内节点未饱和即能达到预期效果.如表 1, 表示倾角为0°、45°、85°时, 在保证不同场分量相位与幅度相对误差均小于0.1%情况下, 所需高斯节点数与相应的积分区间.

显而易见, 计算耗时随倾角不同而不同.经简化表达式、优取积分区间与节点数, 在相对误差小于0.1%时, 即使倾角为85°其计算耗时仍小于0.02 s.如图 7是将模型简化后所得H_zz分量与精确解(Moran and Gianzero, 1979)之间的对比图, 其中图 7a、图 7b是各向同性(R=2 Ωm)对比图, 图 7c、图 7d为单轴各向异性(R_x=2 Ωm, R_y=8 Ωm)对比图, 结果相对误差小于10^-5.

4 双轴各向异性介质多分量感应测井响应特征分析

以同轴、共面分量为实例, 针对电阻率较高与较低的地层, 深入探讨其多分量感应测井的响应特征.假定分别表征地层横向、纵向各向异性.假设β=γ=0°, ω=20 kHz, L=40×2.54 cm, 考虑到油藏电阻率的差异性, 令R_y分别为2 Ωm、50 Ωm, 在不同倾角下, 不同λ_yx、λ_yz的同轴、共面分量视电阻率取对数后的响应特征分别如图 8Ⅰ、图 8Ⅱ所示, 其中横坐标为λ_yz, 纵坐标为λ_yx.

由图 8Ⅰ、Ⅱ可知, 在两种背景电阻率情况下, 其视电阻率响应特征相似：倾角为0°时, R_axx不受横向各向异性影响, R_azz不受纵向各向异性影响, R_ayy同时受横向、纵向各向异性影响, 且纵向各向异性影响更大.随倾角增大, 横向各向异性对R_ayy的影响逐渐增大, 纵向各向异性对其影响逐渐减小, 当背景电阻率较大时, 上述现象更为明显.随倾角增大, 横向各向异性对R_azz的影响逐渐减小, 对R_axx影响逐渐增大, 纵向各向异性对R_azz的影响逐渐增大, 对R_axx影响逐渐减小, 这种现象与背景电阻率无关.倾角接近90°时, R_axx只受横向各向异性影响, R_azz只受纵向各向异性影响, 而R_ayy受横向各向异性影响更大.

综上可知, 倾角较小时, 可用R_axx、R_ayy反映地层纵向各向异性信息, 用R_azz表征地层横向各向异性; 倾角较大时, R_axx、R_ayy反映地层横向各向异性信息, R_azz表征地层纵向各向异性.

5 结论

根据电磁场积分核函数周期性, 针对定积分压缩其积分区间, 可将双轴各向异性介质多分量感应测井解析解计算速度提升至少4倍; 根据电磁场积分核函数振荡特性, 针对半无穷积分自适应截取积分区间, 并根据积分区间优选节点数, 可保证计算速度的同时将误差降低至0.1%.

基于本文提供的双轴各向异性介质多分量感应解析算法, 模拟了倾斜双轴各向异性地层多分量感应测井响应, 得到了同轴、共面分量响应与纵横向各向异性系数、倾角的关系, 有助于理解复杂油气藏电磁响应特征及其各向异性评价.

附录A

沿y方向发射的磁偶极子源：

(A1)

(A2)

(A3)

沿z方向发射的磁偶极子源：

(A4)

(A5)

(A6)