井与井之间砂体界面的确定主要是根据单井上划分的层界面, 运用一定的数学算法插值^[1], 或者以单井确定的地层界面约束地震资料进行层位标定^[2], 然后根据地震资料同相轴追踪形成^[3], 进而形成砂体顶面或底面构造图^[4]。单纯利用单井资料采用插值算法形成的构造图, 由于插值算法不同, 生成的构造面也会有些差别^[5-6]。重磁电方法能够较好反映构造精度, 但其探测深度相对较浅, 无法满足深层油气藏构造界面刻画^[7]。地震资料能够很好地追踪构造面^[8], 也是当前构造面解释最为有效的技术手段; 然而地震资料由于纵向分辨率的原因, 导致井间砂体界面或多或少存在一些误差, 最终造成解释的构造面也会存在误差。随着勘探程度的不断提高, 勘探技术和手段的进步, 水平井的数量增长较快, 它对于提高单井产量和原油采收率^[9]非常重要, 是油田高效开发最重要的技术之一^[10]。油田进行水平井设计和导向时, 精细构造图是提高钻井成功率非常重要的参考资料^[11]。当水平井钻完后, 由于水平井测井资料隐含丰富的地质信息, 其携带的信息是能够用于修正构造面的。利用水平井资料补充和完善传统井间砂体界面确定方法上存在的不足或减小地震构造解释存在的误差, 对于精细油藏地质建模有着非常重要的意义。当前水平井测井系列日趋丰富和完善, 对于地层边界探测能力也逐渐增强^[12], 但国内大多水平井主要为电缆测井的常规9条曲线(绝大多数不包含密度曲线)。在常规曲线中, 电阻率曲线较其他曲线探测深, 笔者以双感应测井为例, 从Maxwell方程出发^[13], 建立双感应三维数值模拟算法, 利用交互式正反演方法^[14-15]实现水平井感应测井确定砂体界面, 弥补地震刻画砂体界面精度上的不足, 从而提高精细三维地质建模精度。

1 水平井感应测井数值正演算法

水平井地层, 由于地层与井眼不再是轴对称关系^[16], 传统二维计算方法不再实用。感应测井中, 发射的电流都是交流电, 发射线圈产生的时变电磁场为时谐电磁场, 时间因子可表示为e^iωt, 所以时谐电磁场的Maxwell方程组^[17]为

$ \nabla \times \mathit{\boldsymbol{E}} = - {\rm{i}}\omega \mu \mathit{\boldsymbol{H}}, $

(1)

$ \nabla \times \mathit{\boldsymbol{H}} = \sigma \mathit{\boldsymbol{E}} + {\rm{i}}\omega \varepsilon \mathit{\boldsymbol{E}} + \mathit{\boldsymbol{J}}, $

(2)

$ \nabla \cdot \mathit{\boldsymbol{D}} = \rho , $

(3)

$ \nabla \cdot \mathit{\boldsymbol{B}} = 0, $

(4)

$ \nabla \cdot \mathit{\boldsymbol{J}} = - {\rm{i}}\omega \rho . $

(5)

式中, E为电场强度, V/m; H为磁场强度, A/m; D为电位移矢量, C/m²; B为磁感应强度, T; J为发射电流密度; A/m²; ρ为电荷密度, C/m³; μ为磁导率, H/m; ε为介电常数, F/m; ω为角频率, rad/s; σ为电导率, S/m。

引入矢量磁位A, B=∇×A, 满足库伦规范∇·A=0, 由此可得到波动方程。

根据Maxwell方程可以知道, 空间无源区域, 电导率为σ(x, y, z)的地层中, 其电场强度E满足波动方程:

$ {\nabla ^2}\mathit{\boldsymbol{E}} + {k^2}\mathit{\boldsymbol{E}} = - {\rm{i}}\omega \mu {\mathit{\boldsymbol{J}}_s}. $

(6)

式中，k为传播常数, 满足k²=iωμσ(x, y, z), 式(6)右边由发射线圈产生, 由于发射线圈的半径相对于求解空间位置非常小, 可以直接将源带入式(6)求解出三维空间电场强度。采用背景场(E^b)和散射场(E^s)结合的办法^[18], 总场E为

$ \mathit{\boldsymbol{E}} = {\mathit{\boldsymbol{E}}^{\rm{b}}} + {\mathit{\boldsymbol{E}}^{\rm{s}}}. $

(7)

式(6)可分解为以下两个方程:

$ {\nabla ^2}{\mathit{\boldsymbol{E}}^{\rm{b}}} - \left( {{\omega ^2}\mu \varepsilon - {\rm{i}}\omega \mu {\sigma _{\rm{b}}}} \right){\mathit{\boldsymbol{E}}^{\rm{b}}} = - i\omega \mu {\mathit{\boldsymbol{J}}_s}. $

(8)

$ {\nabla ^2}{\mathit{\boldsymbol{E}}^{\rm{s}}} + \left( {{\rm{i}}\omega \mu {\sigma _{\rm{r}}}} \right){\mathit{\boldsymbol{E}}^{\rm{s}}} = - {\rm{i}}\omega \mu {\sigma _{\rm{r}}}{\mathit{\boldsymbol{E}}^{\rm{b}}}. $

(9)

式中, σ_b为背景电导率; σ_r为扰动电导率。式(8)可以用解析方法求取, 式(9)采用有限元素法^[19]计算, 得到背景场和散射场大小后, 带入式(7)便可求得整个求解空间电场强度。

2 水平井双感应测井响应分析

依据上述算法思想, 用程序语言建立双感应三维快速数值模拟正演算法, 用于考察双感应在水平井地层中的测井响应。为明确双感应测井能够反映多远的地层界面以及双感应测井对地层界面的敏感性, 本文中主要分析水平井中双感应测井在地层界面附件的响应特征以及不同地层厚时的响应特征。

2.1 不同地层厚度响应特征

采用地层模型为上中下三层介质模型^[20], 如图 1所示, 上下棕色地层为围岩, 中间黄色地层为目的层(砂岩), 仪器(灰色部分)居于目的层中部, 与地层平行, 不考虑井眼尺寸及钻井液影响, 设砂岩电阻率R_t为20 Ω·m, 厚度H从1 m变化到20 m, 上下围岩无限厚, 围岩电阻率R_s分别为1、2、4、6、10、30、50、100、200 Ω·m, 模拟不同地层厚度H情况下水平井中双感应测井的响应特征。图 2为感应测井随地层厚度变化时的模拟结果。从图 2可以看出, 当地层厚度H逐渐增加时, 不论围岩是高阻(R_s＞R_t)还是低阻(R_s＜R_t), 双感应视电阻率值都逐渐接近目的层(砂岩)真实值R_t; 当地层较薄时, 双感应视电阻率偏离目的层真实值, 地层越薄, 围岩电阻率与目的层电阻率相差越大, 双感应视电阻率偏离目的层真实值越远, 表明地层越薄, 受围岩影响越大。水平井中, 通常采用探测深度来描述仪器探测边界的能力, 从图 2可以看出, 均匀地层中, 深感应在地层厚度大于10 m(图 2(a)), 中感应在地层厚度大于8 m(图 2(b))时, 双感应视电阻率值基本不随围岩电阻率和地层厚度变化而变化, 说明此时双感应基本不受上下围岩影响, 即双感应测井能够反映8 m内地层界面的信息; 但按照双感应几何因子80%贡献率, 深感应在地层厚度大于4 m, 中感应在地层厚度大于2 m时, 测量结果基本不受上下围岩影响, 即实际测量过程中双感应对2 m内探测到的界面信息较为可靠。

2.2 层界面响应特征

建立两层地层模型, 假设地层无限厚, 井眼与地层界面平行, 如图 3(a)所示, 仪器(图 3(a)灰色部分)从目的层(R_t=20 Ω·m)逐渐靠近并进入围岩(围岩电阻率R_s分别为1、4、8、30、40 Ω·m)。从图 3(b)可以看到, 当围岩电阻率比较低时(相当于目的层高阻), 从高阻地层逐渐进入低阻地层时, 视电阻率并不是逐渐降低, 而是在地层界面(ds)附近先升高后降低, 在界面处形成一个尖角(图 3(b), R_s=1 Ω·m黑线); 而当围岩电阻率较高时(相当于目的层低阻), 从目的层进入围岩, 电阻率逐渐升高, 界面处基本上没有变化(图 3(b), R_s=40 Ω·m红线)。

当地层电阻率相对低时(目的层高阻), 仪器在低阻层内(即围岩内)距离地层界面大于1 m时, 双感应视电阻率基本接近地层真实值, 表明此时双感应测量值基本不受上覆围岩(目的层)影响(图 3(b), 下部), 而当目的层电阻率相对较高(围岩低阻), 仪器在目的层内距离地层界面大于4 m时, 双感应视电阻率基本接近地层真实值, 表明此时双感应不受下覆围岩影响(图 3(b), 上部)。以上分析表明, 地层界面上下地层电阻率相差越大, 双感应对地层界面反应越灵敏; 仪器所处地层电阻率高于围岩, 仪器能够反应的界面距离越远, 反之, 能够反应的界面距离越近。

3 反演算法

双感应实际测量值可以表示为一个N维的列向量d, 同样, 地层参数(如地层厚度、地层电阻率、地层个数等)也可以表示为M维的列向量m, 即

$ \mathit{\boldsymbol{d}} = {\left[ {{d_1},{d_2}, \cdots ,{d_N}} \right]^{\rm{T}}}, $

(10)

$ \mathit{\boldsymbol{m = }}{\left[ {{m_1},{m_2}, \cdots ,{m_M}} \right]^{\rm{T}}}. $

(11)

此时地层参数m和测井响应y的关系表示如下:

$ \mathit{\boldsymbol{y}} = \mathit{\boldsymbol{f}}\left( {\mathit{\boldsymbol{m}},\mathit{\boldsymbol{d}}} \right). $

(12)

上述关系式相当于N个函数, 它们之中每一个并不是线性关系, 为此在求解时是通过寻找一组最优解使得方程的平方和误差最小:

$ \mathit{\boldsymbol{ \boldsymbol{\varPhi} }} = \sum\limits_{i = 1}^N {{{\left| {{\mathit{\boldsymbol{f}}_i} - {\mathit{\boldsymbol{y}}_i}} \right|}^2}} < {e_p}. $

(13)

其中e_p为给定的误差范围, 式(13)等价为求$\mathit{\boldsymbol{ \boldsymbol{\varPhi} }} = \sum\limits_{i = 1}^N {f_i^2} $的极小值, 定义函数:

$ \mathit{\boldsymbol{ \boldsymbol{\varPhi} }} = \frac{1}{2}\mathit{\boldsymbol{f}}{\left( x \right)^{\rm{T}}}\mathit{\boldsymbol{f}}\left( x \right). $

(14)

若f在求解区域上可微, 即Φ对x求导:

$ \mathit{\boldsymbol{g}}\left( x \right) = \mathit{\boldsymbol{Df}}{\left( x \right)^{\rm{T}}}\mathit{\boldsymbol{f}}\left( x \right) = 0. $

(15)

其中

$ \mathit{\boldsymbol{Df}}{\left( x \right)^{\rm{T}}} = \left[ {\begin{array}{*{20}{c}} {\frac{{\partial {y_1}}}{{\partial {x_1}}}}&{\frac{{\partial {y_2}}}{{\partial {x_1}}}}& \cdots &{\frac{{\partial {y_m}}}{{\partial {x_1}}}}\\ \vdots&\vdots &{}& \vdots \\ {\frac{{\partial {y_1}}}{{\partial {x_n}}}}&{\frac{{\partial {y_2}}}{{\partial {x_n}}}}& \cdots &{\frac{{\partial {y_m}}}{{\partial {x_n}}}} \end{array}} \right]. $

将f(x)在x_i点线性展开:

$ \mathit{\boldsymbol{f}}\left( x \right) = \mathit{\boldsymbol{Df}}\left( {{x^i}} \right)\left( {x - {x^i}} \right) + \mathit{\boldsymbol{f}}\left( {{x^i}} \right). $

(16)

带入式(15)得到的解记为xⁱ⁺¹, 则有

$ {x^{i + 1}} = {x^i} - {\left[ {\mathit{\boldsymbol{Df}}{{\left( {{x^i}} \right)}^{\rm{T}}}\mathit{\boldsymbol{Df}}\left( {{x^i}} \right)} \right]^{ - 1}}\mathit{\boldsymbol{Df}}{\left( {{x^i}} \right)^{\rm{T}}}\mathit{\boldsymbol{f}}\left( {{x^i}} \right). $

(17)

令G(x)=Df(x)^TDf(x), 将式(15)带入式(17), 则式(17)可变为

$ {x^{i + 1}} = {x^i} - \mathit{\boldsymbol{G}}{\left( {{x^i}} \right)^{ - 1}}\mathit{\boldsymbol{g}}\left( {{x^i}} \right). $

(18)

如果G(x)正定, 方程(18)改写为

$ {\mathit{\boldsymbol{w}}_i} = - \mathit{\boldsymbol{G}}{\left( {{x^i}} \right)^{ - 1}}\mathit{\boldsymbol{g}}\left( {{x^i}} \right) = - \mathit{\boldsymbol{G}}{\left( {{x^i}} \right)^{ - 1}}\nabla \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}\left( {{x^i}} \right), $

(19)

则

$ \mathit{\boldsymbol{w}}_i^{\rm{T}}\nabla \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}\left( {{x^i}} \right) = - \mathit{\boldsymbol{w}}_i^{\rm{T}}\nabla \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}\left( {{x^i}} \right){\mathit{\boldsymbol{G}}_i} < 0. $

(20)

上式表明w_i与-∇Φ(xⁱ)方向一致, 为防止G(xⁱ)奇异或病态, 引入一阻尼因子γ,

$ \mathit{\boldsymbol{\hat G}}{\left( {{x^i}} \right)^{ - 1}} = \mathit{\boldsymbol{G}}\left( {{x^k}} \right) + \gamma . $

(21)

当阻尼因子γ＞0时, G(x^k)是对称正定矩阵^[21], 因此有

$ {\mathit{\boldsymbol{w}}_i}\left( \gamma \right) = - \mathit{\boldsymbol{G}}{\left( {{x^i}} \right)^{ - 1}}\nabla \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}\left( {{x^i}} \right). $

(22)

式(22)即为函数Φ在(xⁱ)处的下降方向。于是可构造迭代算法:

$ {x^{i + 1}} = {x^i} - \mathit{\boldsymbol{\hat G}}{\left( {{x^i}} \right)^{ - 1}}\nabla \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}\left( {{x^i}} \right). $

(23)

式(23)称为阻尼最小二乘法, 当γ＞0时, 总Φ(xⁱ⁺¹)＜Φ(xⁱ)可保证, 故{xⁱ}总是收敛的, 但当γ较大时, {xⁱ}收敛速度下降, 而当γ太小, 则收敛域较小, 因此选择合适的γ较为困难, 通常原则上γ的取值应为较小的正数如0.01, 然后通过程序自动调节γ, 确定搜索方向^[19], 求响应函数f的解, 具体算法如图 4所示。

为验证反演算法可靠性, 以图 3(a)地层模型为基础, 开展两组地层模型反演处理; 地层模型中, 目的层电阻率R_t为1 Ω·m, 围岩电阻率R_s为20 Ω·m; 对仪器位于目的层中以及仪器位于围岩中两种情况进行反演, 具体模型参数如表 1所示。

模型一, 仪器位于低阻目的层内, 当地层界面距离仪器2 m时, 反演地层界面误差为13.2 cm, 误差相对小于2 m距离时较大, 由于目的层低阻, 能够反应的界面距离较近, 当距离大于2 m时, 由于反应界面不灵敏, 导致反演结果误差相差较大; 模型二, 仪器位于高阻围岩内, 当地层界面距离仪器3.5 m时, 反演地层界面误差为12 cm, 比模型一中2 m距离误差还小, 这是由于地层高阻时能够反应的界面距离较远, 这与上文分析一致。图 5表明两组模型反演界面距离与实际模型设定的距离基本都在45°线分布, 整体上误差较小, 反演算法可以用于实际资料处理。

4 地质应用实例

根据以上算法, 实际资料处理时, 具体步骤如下:

(1) 水平井大斜度段TVD校正;

(2) TVD校正后的曲线与导眼井或邻井对比, 确定水平井钻遇目的层;

(3) 根据曲线特征, 将导眼井或邻井目的层段上下10 m(根据仪器探测深度以及轨迹走向和地层的构造情况调节)地层细分层, 并标上砂体界面位置;

(4) 反演确定每一细分层的地层电阻率;

(5) 根据地震剖面确定砂体界面的展布情况、走向等;

(6) 调节地层模型和水平井井眼轨迹的相互关系、进出层位置;

(7) 反演出砂体界面沿井眼轨迹的走向。

以某油田J区块一口水平井为例。图 6为该油田J区块根据井标定后的地震解释形成的砂体Ta顶面构造图, 该构造图主要根据地震构造解释得出, 由于解释时并没有考虑到水平井信息, 使得解释成果后, 该区有9口水平井在构造面上部, 即没有出现在砂体Ta中, 而是在上覆泥岩层中穿行, 与实际情况不符合。鉴于该区水平井覆盖率较大, 可以充分利用水平井信息来标定地震解释的井间砂体顶面构造的展布情况。以其中一口水平井XX25H井为例, 说明解释过程。

XX25H井为从南向西北钻探的一口水平井(图 6), 按照上述流程, 首先进行TVD校正, 确定该井水平段钻探主要目的层位。图 7为水平井大斜度段TVD校正后与导眼井对比图(其中，SP为自然电位曲线，GR为伽马曲线，RILD为深感应电阻率，RILM为中感应电阻率，RFOC为微球聚焦电阻率，AC为声波时差曲线)。从图 7可以很明确地看出, 该水平井主要钻探层段为对应导眼井的4 603.5~4 638.05 m段, 即图 7中蓝色层(层号4), 图 7左图第三道为导眼井双感应电阻率和微球聚焦电阻率, 电性曲线上可以很明显看出4号砂体为低阻油藏, 4号层电阻率比上覆3号泥岩电阻率小, 油层电阻率约为1~1.5 Ω·m, 导眼井4号砂体顶界面Ta在4 603.5 m。

为便于反演砂体顶界面, 将导眼井砂体顶界面4 603.5 m上下15 m井段地层根据曲线变化, 采用活化函数法^[22]进行分层(该分层目的是为了反演而细分), 分层确定每一小层边界后, 结合原始测量曲线、井眼尺寸及钻井液电阻率, 并利用本文中反演算法对感应电阻率进行反演得到每一细分层的电阻率, 如图 8所示, 图 8第5道台阶状曲线为反演后的每一细分层的地层电阻率(FRT), 每一台阶表示细分层的界面, 第6道为地层模型, 颜色深浅表示地层电阻率高低。据此, 完成导眼井地层模型建立, 并将该模型用于水平井作为初始模型。

将上述导眼井建立的地层模型用于水平井, 首先根据测井曲线特征, 初步调整地层模型, 通过曲线对比, 确定水平井从4 723 m进入砂体顶界面Ta, 然后, 根据地震剖面确定地层的大体倾向, 最后采用本文所述反演流程算法, 反演出砂体Ta顶界面沿井眼轨迹的走向, 图 9为最终成果图。图 9从上到下, 第1道蓝色为自然伽马(GR)曲线, 橙色为自然电位(SP)曲线, 第2道为深度道, 第3道为原始测量深感应电阻率曲线(蓝色)与正演模拟的深感应电阻率曲线(RILD_M, 红色虚线), 第4道为原始测量中感应电阻率曲线(蓝色)与正演模拟的中感应电阻率曲线(RILM_M, 红色虚线), 第5道为声波时差(AC, 蓝色)曲线和中子(CNL, 黑色)曲线, 第6道为井眼轨迹和地层相互关系, 红色为实钻井眼轨迹, 沿井眼轨迹红色虚线为双感应有效探测范围, 蓝色实线即为反演出的砂体顶界面Ta, 通过响应特征分析, 水平井中感应测井探测范围为8 m, 但井筒2 m范围内(即井沿井眼轨迹上下红色虚线框内)的影响最大, 故认为2 m内的反演结果较为可靠, 大于2 m的层界面主要根据2 m内细分层的层界面趋势调整, 第3道和第4道正演的双感应曲线是根据井眼轨迹与地层关系, 结合实际测量仪器参数正演得出, 正演结果与实测结果基本一致, 说明反演构造界面可靠。

图 9第6道黑色实线为在三维地质模型中, 沿井眼轨迹切一条剖面, 得到的砂体顶面构造线, 该构造线为地震刻画的砂体顶界线, 可以看出与蓝色实线(反演后的Ta界线)存在较大差别, 井眼从4 723 m接触到砂体界面开始到水平段结束, 地震解释构造线TVD深度整体要低于反演后的构造线, 即地震解释该井完全在3号泥岩层穿行, 然而该井初期日产油86.5 t, 显然水平段在砂体里穿行, 地震解释结果与实际情况不符。反演后Ta界线与地震解释Ta界线TVD差如图 10所示, 二者相差1.5~3.5 m(正数表示反演后的砂体界面TVD深度小于地震解释构造线TVD深度), 这些差异主要是由测井与地震的分辨尺度不同造成的, 由于测井的精度要高于地震精度, 因此利用该方法反演的砂体界面可靠性较高。

对研究区其他水平井做相同处理, 利用水平井反演砂体界面约束后的砂体顶界面构造图如图 11所示, 对比图 11和图 6发现, 整体上, 二者的构造趋势一致, 这是由于图 10反演的基础也是地震剖面约束的整体走向和倾向, 但由于图 11反演过程中井间砂体构造面充分利用水平井的测井信息, 因此在局部地方与图 6存在差别, 如图 11中A处、B处都是分开的, 而图 6相应地方为连续, 且B处构造形态存在较大差别, 图 11中C处为3个分散的小高点, 而在图 6中对应地方为两个大的高点。反演后, 该区原9口水平井全部在层内, 由此可以得出, 充分利用水平井测井信息反演后的构造图可以提高油田精细地质建模精度, 为后期开发方案部署提供坚实依据。

5 结论

(1) 数值模拟结果表明, 水平井中双感应测井能够反映2 m内地层界面信息。

(2) 地层界面上下地层电阻率数值相差越大, 双感应测井对地层界面信息反应越灵敏; 当仪器所处地层电阻率比围岩高, 则双感应能够反应的地层界面越远。

(3) 利用水平井感应测井资料进行构造解释, 可以辅助弥补地震由于精度不高造成的误差, 为油藏精细地质建模提供更为可靠的参数。

(4) 该方法基于砂体横向变化不大的储层, 将砂体界面作为唯一输出反演参数, 因此对于油藏内砂体横向变化大、非均质性强的地区, 并不一定实用, 但反演思想和流程是可以借鉴的。

(5) 水平井测井资料隐含丰富的地质信息, 对于水平井覆盖率较高的地区, 本文中所述方法可以达到精细砂体构造面的目的, 然而感应测井探测深度相对较浅, 探测范围外的砂体界面可信度降低。当前针对随钻电磁波探边仪器研发力度加大, 远探测方位随钻仪器的研制及其解释方法的应用将是水平井未来发展的重要方向, 这将加强水平井测井资料地质解释应用程度。

致谢 感谢中国石化石油勘探开发研究院赵磊博士的帮助。