基于多元回归模型的钻头前方地层速度预测方法*

郭书生 朱江梅 王世越 赵前华 王晓飞

(1. 中海石油(中国)有限公司海南分公司 海南海口 570311； 2. 中海油田服务股份有限公司油田技术事业部 河北廊坊 065201)

地震波速度是在地震资料的采集、处理、解释和评价研究中的关键参数，地震波速度求取的准确与否直接影响着地震勘探资料在应用中的可信度。随着油气勘探程度的不断提高，勘探区域地质条件变得越来越复杂、速度横纵向变化大等问题，对高精度速度预测提出了挑战。为此，各种各样的地震速度预测技术应运而生。这些地震速度预测方法按照数据基础分为三类：①地震资料为数据基础的叠加速度分析[1-3]、偏移速度分析[4-6]和层析速度反演[7-8]等速度预测方法；②以测井资料为数据基础的克里金估计[9-10]、随机模拟[11]和随机反演[12-13]等速度预测方法；③以井震资料结合为数据基础的井震匹配[14]、井震约束[15]等速度预测方法。

本文借鉴“ look ahead VSP ”思路[16]，也就是基于VSP资料进行速度预测。由于VSP技术是将井、震融于一体的一种可靠地球物理勘探手段，可为预测高精度的垂向地震速度提供更有力支撑。而基于VSP资料的钻前预测研究主要集中在钻前压力预测与深度预测方面[17-18]，而这些研究归根还是有赖于高精度钻前速度的预测，国内外许多学者对此进行了有关研究。如Khaitan等提出的基于前瞻VSP调查理论的前瞻VSP反演方法[19]；Esmersoy等提出的通过反演地震波走时估计速度的钻前速度实时预测方法[20]；张科 等提出的根据经验公式、局部构造形态做控制结合实际钻井计算钻前速度的综合速度分析方法[21]和一些基于射线追踪理论和VSP走廊叠加反演波阻抗估算速度的常规方法[22-23]。但是这些方法都是侧重于利用单一因素来预测钻头前方地层速度，没有关注VSP地震资料的多属性高频先验信息。

基于以上分析，本文将从以下两个思路展开研究：一方面，利用高分辨率VSP走廊叠加，作为高精度预测的基础；另一方面，针对现有预测技术主要关注单一因素来预测钻头前方地层速度的缺陷，提出了使用多元回归模型量化刻画客观因素对地层速度的影响，从而为高精度预测提供更好的算法基础。具体工作步骤可归纳为：首先，对基于VSP走廊叠加资料的多属性信息构建多元回归模型；其次，反演预测中频地层速度，最后再将中频地层速度等因素特征作为输入，再通过粒子滤波法[24-25]来求解该多元回归模型的状态参数，从而有效地预测钻头前方地层速度。以南海莺歌海盆地L工区L1井为目标研究区块，应用提出的方法对钻头前200 m的地层速度进行预测，以验证本文方法的有效性和稳健性。

1 算法的提出

1.1 问题描述与模型

基于以上数据条件，我们将钻头前方地层速度预测本质问题转化为有监督机器学习的回归问题，已钻地层速度数据作为回归模型的标签数据。由于涉及多个特征量的输入，所以该回归问题又属于多元回归问题的范畴，如图1所示，三维坐标空间表示由输入属性张成的空间，曲面则表示拟合得到的回归平面。

图1 多元回归示意图Fig.1 Schematic diagram of the multiple regression model

对于多元回归模型，核心是用系统状态方程对问题进行建模。根据系统状态方程，这里有两方面内容：一是将待估计的钻头前方地层速度建模为一阶马尔科夫链，也就是待估计速度序列之间为时间相关，可以用逼近函数f(·)表示；二是输入的多个特征序列和待预测速度之间为空间相关，可以用逼近函数g(·)表示。时间相关和空间相关可以用以下系统方程和状态测量方程进行刻画，具体表达式如下：

状态转移方程为

xk=f(xk-1)+rk

(1)

状态测量方程为

yk=g(xk)+vk

(2)

式(1)、(2)中：xk代表系统是tk时刻的地层速度，值得注意对于时刻tk而言此处为向量；yk代表叠加走廊和各种属性，包括反演结果;rk、vk分别代表过程噪声和测量噪声。对应于本研究，考虑建模过程是存在误差的，rk是时变建模过程中存在的误差，主要是线性或者非线性的方式逼近导致产生的误差。vk是空间关系建模所导致的误差，因为非线性函数g(·)需要泰勒展开线性化求解，而线性化则会省略多次项引入误差，已钻井段地层速度和VSP走廊叠加数据测量过程中存在的误差。

1.2 基于反演的中频地层速度预测

钻头前方未知的地层速度因其本身的复杂性，所以，我们需要分阶段来逼近预测。我们先通过反演获得钻头前方地层速度的中频信息，本质上中频信息带宽与走廊叠加资料的带宽接近。具体来说VSP走廊叠加反演需要在正演基础上才能够进行，当然正演模型在进行输入数据井震匹配时也是需要的。对于VSP走廊叠加正演模型，本文采用了常用的褶积模型进行构建。该模型定义如下：单道VSP走廊叠加S(t)可以建模为地震子波W(t)与反射系数R(t)卷积的结果：

S(t)=W(t)*R(t)

(3)

由于ΔZi远小于阻抗Zi，因此反射系数可近似用式(4)表示：

(4)

式(4)中：Zi=ρV，ΔZi为不同深度之间阻抗变化。在实际砂泥岩的预测过程中，因为密度ρ的变化范围较小，所以可以将密度ρ定义为常值，这样就可以利用Zi求取速度地层速度V。因此，式(3)可以表示为

(5)

式(5)中：St表示时刻t的走廊叠加资料;wt表示时刻t的子波。在实践中，反演问题的目的是根据S和G的知识给出X的估计。然而，由于这个问题的不适定性，方程大部分是使用最小平方正则化，其可表示如下：

(6)

f(x)=‖x‖1

这里‖·‖1为L1范数，该目标函数的求解有许多求解方法，比如共轭梯度[26]等，鉴于篇幅不在这里赘述。

1.3 完整预测算法

为了进一步提升地层速度预测效果，需要将反演的中频地层速度和属性数据一起输入多元回归模型，再利用粒子滤波算法来求解多元回归模型。作为典型的回归模型，它可以具体分成2个阶段，第1个阶段是训练阶段，也就是已知地层速度情况下，训练1个多元回归模型；第2个阶段是预测阶段，输入新的反演速度和属性数据，预测新的地层速度。

在粒子滤波算法中，对于一个动态非线性离散时间系统，其状态空间模型通常写为：

状态转移方程为

xk,1=Axk-1,1+rk,1

(7)

状态测量方程为

yk=g1(xk,1)+vk,1

(8)

式(7)、(8)中：A表示将公式(1)中f(·)线性化，具体通过历史数据拟合得到，所谓历史数据是指已钻探段采集得到的速度；xk,1=[dk,hk]T代表在第k个时刻地层速度预测所对应的模型参数；rk,1、vk,1分别代表参数状态转移的过程噪声和地层速度的测量噪声。

在粒子滤波算法中，系统状态的概率分布估计是通过一组具有权重的随机样本(称为粒子)给出的。对于一个粒子数是ND的粒子滤波算法，具体计算步骤如下:

图2 粒子滤波器示意图[24]Fig.2 Schematic diagram of particle filter[24]

在利用粒子滤波估计出速度预测的多元回归模型的状态参数后，将根据式(7)、(8)可直接对钻头前方地层速度进行预测。

当然，在不同的地质环境下，影响钻头前方地层速度的客观因素并不仅仅只有反演结果和叠加走廊差分衍生属性，本文提出的方法作为一个开放结构，在未来研究工作中，更多样和密切的客观因素可以以系统状态方程的特征输入到我们的方法中，进一步提升钻头前方地层速度的预测精度。

1.4 具体算法步骤

本文提出的算法具体步骤如下：

输入：已钻井段地层速度X，VSP走廊叠加Y

2) 求解公式(6)的中频地层速度x；

3) 将以上地震属性和中频地层速度输入回归模型，利用粒子滤波器进行求解；

②根据粒子的重要度权值重采样；

输出：yk

2 应用实例分析

本文以莺歌海盆地L工区L1井为目标研究区块，该井为高温超高压井，目的层温度为179～188 ℃。该井在钻到3 834.2 m后，根据邻井的信息，认为前方将会钻遇超高压地层的台阶面，希望将当前井段的中完在超高压地层的台阶面之上，越接近超高压地层的台阶面越好，因此决定测量中途VSP以标定层位深度及预测下部井段的地层信息。VSP资料采集的最大深度为3 826 m。图3显示了基于多元回归模型的钻头前方地层速度预测的结果。图中第2列为录井岩性剖面(后续实钻)及已知井段通过VSP测量时深关系转换的地层速度，第3列是垂直地震走廊叠加道；第4列是常规地震的井旁道；第5列是后续实钻后测量的伽马测井曲线。第6列是后续实钻后测量的声波测井曲线(洋红色)及利用本文方法预测的获得速度(黑色)和常规波阻抗反演获得的速度(绿色)。通过图3可以看出，首先垂直地震走廊叠加道(第3列)的分辨率明显比常规地震的井旁道(第4列)高，有助于钻头前方的砂体的识别。其次，在第6列反演预测地层速度和已钻井段实测地层速度在已钻井段(3 826 m以上)非常一致。未钻井段(3 826 m以下)经过后续实钻表明，此次预测的速度趋势与实测数据相近。

图3 垂直地震走廊叠加预测速度与实钻结果对比Fig.3 Vertical seismic data corridor stack section predicted velocity vs. actual drilling

根据此结果对目的层进行预测，见表1钻前预测部分。本井实时钻前深度预测结果与后续实钻结果有很好的一致性，对目的层T31，T31-A1深度的预测和实钻深度误差0 m，对目的层T31-B1的预测和实钻深度差为-3 m，对目的层T31-C1的预测和实钻深度差为12 m。对于目的层T31-C1，从录井资料和伽玛曲线分析，该层为砂泥岩互层，波阻抗变化比较频繁，而VSP资料的纵向分辨率不及电缆测井，所以对急速变化的波阻抗特征反应不及时，加上预测距离增大等因素的影响，导致预测误差较大。

表1 目的层预测深度与实钻结果对比Table 1 Predicted results of the depth of the target layer

3 结论

本文提出了基于多元回归模型的钻头前方地层速度预测方法，该方法通过建立非线性多元回归模型，充分利用了VSP资料中的多属性的信息，也即叠加走廊差分衍生属性(比如振幅导数、瞬时振幅导数等)，并使用非线性粒子滤波法对模型进行求解。通过在莺歌海盆地L工区L1井目标研究区块的实际应用表明，本文方法反演预测的速度比常规波阻抗反演获得的速度频率获得了提高，同时速度变化的趋势与实测数据更为相近，验证了本文方法的有效性。