CN106842299B - 一种基于地震属性的裂缝定量化预测的方法 - Google Patents

Abstract

本发明提供了一种基于地震属性的裂缝定量化预测的方法。该方法包括：获取目的层的三维地震数据和测井数据，基于这两项数据获得目的层的构造解释；对三维地震数据进行去噪处理；对去噪处理后的三维地震数据体进行相干属性计算，得到相干属性体；对去噪处理后的三维地震数据体进行曲率属性计算，得到曲率属性体；对相干属性体进行值域校正；基于曲率属性体以及校正后的相干属性体，计算得到裂缝预测属性体；根据裂缝预测属性体对地下裂缝的空间展布信息进行预测。本发明提供的技术方案能够快速准确地定量化描述地下裂缝空间展布信息。

Description

一种基于地震属性的裂缝定量化预测的方法

技术领域

本发明涉及一种基于地震属性的裂缝定量化预测的方法，属于地球物理勘探技术领域。

背景技术

自近代油气勘探开始以来，已逾百年历史。油气勘探目标已由早期的寻找构造圈闭，到近年来寻找岩性圈闭。随着勘探的不断深入，人们发现在地下许多地层中有裂缝发育，裂缝发育的地层往往储层也较为发育，相比裂缝不发育的储层还伴随着油气高产。研究认为，裂缝不仅是渗流通道，也是重要的储集空间，是影响储层发育的重要因素之一，所以准确预测地下裂缝空间展布是油气勘探的重要研究内容。

由于地下构造复杂，裂缝的空间展布形态千变万化，对其准确预测难度较大。常用的地震裂缝预测技术包括相干、曲率、蚂蚁追踪、体积膨胀系数和叠前方位各向异性分析等技术。

相干技术是目前应用最广泛的裂缝预测技术，通过地震道组合分析空间不连续性来预测裂缝，可以较好刻画裂缝空间展布趋势，该技术到目前为止已经发展了三代。第一代技术主要基于归一化互相关算法，由Bahorich等人提出；第二代技术基于多道相似性测量算法，通过构建协方差矩阵计算地震道之间的相似性，进一步增加了抗噪性和稳定性，由Marfurt等人提出；第三代技术基于本征结构算法，通过矩阵的特征结构计算相似性，在第二代算法基础上提高了横向分辨率。相干算法发展至今可以较好预测各种尺度裂缝形态的空间展布，但只能进行定性预测，未能达到定量化描述。

曲率技术作为地震几何属性的一种，近几年来在构造识别和地震解释方面得到迅速发展和应用。曲率在数学上是用于度量曲线上某一点弯曲程度的数值，表明曲线偏离直线的程度。曲率越大表示曲线的弯曲程度越大。对于直线而言，不管为水平形态或者倾斜形态，其曲率始终为零。曲率属性用于描述地质体的几何变化，与地震反射体的弯曲程度相对应，对岩层的弯曲、褶皱和裂缝、断层等反应敏感，是寻找地层裂缝空间展布的有效手段。由于裂缝密度信息与地层的弯曲程度相关，所以可以利用曲率属性定量化描述裂缝密度信息。但是由于岩层弯曲和褶皱这些非裂缝因素也可以引起曲率响应，所以利用曲率属性预测裂缝得到的结果往往多于地下实际存在的裂缝。

其余方法也各有优缺点。蚂蚁追踪属性基于蚂蚁算法刻画地下裂缝空间展布，可以定性描述中小尺度裂缝的空间展布形态，但空间展布趋势不清晰。体积膨胀系数通过对比地层原始沉积形态体积和现今形态体积预测地下裂缝发育情况。叠前方位各向异性分析技术通过分析不同方位角地震属性的差异预测地下裂缝发育情况。体积膨胀系数和叠前方位各向异性分析技术都可以定量化预测裂缝发育情况，但裂缝空间展布预测精度较低。

可见，上述单项技术都无法准确地定量化预测裂缝的空间展布信息。

在实际生产过程中，研究人员往往面对叠后地震数据。如何利用叠后地震数据准确地定量化预测地下裂缝的空间展布信息成为地震储层描述综合研究亟待解决的问题。

发明内容

为解决上述技术问题，本发明的目的在于提供一种基于地震属性的裂缝定量化预测的方法，该方法能够利用叠后地震数据快速准确地定量化描述地下裂缝空间展布信息，为储层评价及油气勘探提供技术支撑。

为达到上述目的，本发明提供了一种基于地震属性的裂缝定量化预测的方法，其包括：

获取目的层的三维地震数据和测井数据，基于这两项数据获得目的层的构造解释；

对三维地震数据进行去噪处理，得到去噪处理后的三维地震数据体；

对去噪处理后的三维地震数据体进行相干属性计算，得到相干属性体；

对去噪处理后的三维地震数据体进行曲率属性计算，得到曲率属性体；

对相干属性体进行值域校正，得到校正后的相干属性体；

基于曲率属性体以及校正后的相干属性体，计算得到裂缝预测属性体；

根据裂缝预测属性体和目的层的构造解释对地下裂缝进行预测；

对相干属性体进行值域校正包括：

设定门槛值，相干属性体中采样点的值小于所述门槛值时，则确定该采样点进行值域校正后的值为1；相干属性体中采样点的值大于或等于所述门槛值时，则确定该采样点进行值域校正后的值为0；

所述门槛值的取值为0≤所述门槛值≤1；

所述裂缝预测属性体的计算公式如式13所示

Z_i＝E_i′×Q_i 式13

在式13中，i＝1，2，…，M；Z_i为裂缝预测属性体中的采样点的值；E′_i为校正后的相干属性体中对应位置的采样点的值；Q_i为曲率属性体中对应位置的采样点的值。

在上述方法中，优选地，所述三维地震数据为保幅处理后的偏移叠加地震数据；所述测井数据为地震工区内钻穿目的层的测井数据；井-震标定后进行地震数据目的层构造解释，其是在前述的三维地震数据和测井数据基础上得到的，可以用于后期对目的层段开展裂缝解释提供地层层段的顶底界面。

与只能定性预测裂缝的相干属性体、蚂蚁追踪属性体相比，本发明提供的技术方案通过将可以较好预测裂缝空间展布形态的地震相干属性体和能够定量化预测裂缝密度信息的曲率属性体进行融合，并进一步消除了非裂缝信息的干扰，从而可以利用叠后地震数据快速准确地定量化描述地下裂缝空间展布信息，最终获得的预测结果相对于可以定量化预测裂缝的体积膨胀属性体和叠前方位各向异性分析技术而言，预测精度更高；与此同时，与可以定量化预测裂缝的曲率属性体而言，本发明提供的技术方案消除了非裂缝信息的干扰，预测结果更加准确，为储层评价及油气勘探提供了技术支撑。

在上述方法中，优选地，对三维地震数据进行去噪处理时，所采用的去噪方法包括中值滤波法或均值滤波法。

在上述方法中，优选地，对三维地震数据进行去噪处理时包括：

在三维地震数据体中定义一个时窗，该时窗包含2j+1个样本点；此处样本点的个数需要小于地震数据体纵向时间域样本点的个数，且样本点的个数为奇数，故此处设定为2j+1个样本点；

采用中值滤波法或均值滤波法计算所述时窗中的第j+1个样本点的值；其中，

所述中值滤波法的计算公式如式1所示

在式1中，i＝1，2，…，N；为采用中值滤波法获得的时窗中第j+1个样本点的数值；

所述均值滤波法的计算公式如式2所示

在式2中，i＝1，2，…，N；为采用均值滤波法获得的时窗中第j+1个样本点的数值。

在上述方法中，优选地，对去噪处理后的三维地震数据体进行相干属性计算，得到相干属性体包括以下过程：

过程1，将去噪处理后的三维地震数据体作为一个分析窗口，该分析窗口中有N道地震数据，K个采样点，则所述去噪处理后的三维地震数据体表示为式3所示的矩阵D：

在式3中，d_KN为第N道的第K个样点值；

过程2，按式4所示的公式计算上述分析窗口中的三维数据体的协方差矩阵C

在式4中，为矩阵D中的第k行向量其表示数据体的第k个采样点的集合；d_k为的转置，d_k＝[d_k1,d_k2,…,d_kN]^T；

过程3，假设λ_k是协方差矩阵C的第k个特征值，λ₁是其中最大的特征值，则相干属性的计算公式如式5所示，按照式5所示的公式计算得到相干属性体；

在式5中，λ_k是协方差矩阵C的第k个特征值，λ_m是其中最大的特征值。

在上述方法中，优选地，所述对去噪处理后的三维地震数据体进行曲率属性计算，得到曲率属性体包括：

将去噪处理后的三维地震数据体先转化为倾角数据体；

根据倾角数据体计算得到数据体内任意反射点的曲率。

在上述方法中，优选地，将去噪处理后的三维地震数据体转化为倾角数据体包括以下过程：

过程1，计算三维地震数据体的瞬时频率；

过程2，计算三维地震数据体的瞬时波数；

过程3，根据得到的三维地震数据体的瞬时频率和瞬时波数，计算三维地震数据体的倾角。

在上述方法中，优选地，所述计算三维地震数据体的瞬时频率包括：已知三维地震数据体在空间任意点是时间标量u(t,x,y)，根据式6所述的计算公式，计算得到瞬时频率ω

式6中，u^H(t,x,y)为所述时间标量u关于时间t的希尔伯特变换。

在上述方法中，优选地，所述计算三维地震数据体的瞬时波数包括：按式7所示的公式，分别计算三维地震数据体沿x方向和y方向的瞬时波数

式7中，k_x为x方向的瞬时波数，k_y为y方向的瞬时波数。

在上述方法中，优选地，根据得到的三维数据体的瞬时频率和瞬时波数，利用式8所示的公式计算三维地震数据体的倾角dip(px，qy)：

式8中，p_x、q_y分别为x方向和y方向的倾角分量。

在上述方法中，三维空间中某点的曲率值通过与其相邻道和样点的倾角值拟合出的空间区面方程进行计算；优选地，根据倾角数据体计算得到数据体内任意反射点的曲率包括以下过程：

过程1，根据式9所示的最小二乘逼近原理，对方程两边求微分后，代入倾角分量p_x、q_y，得到式10所示的方程系数：

z(x,y)＝ax²+by²+cxy+dx+ey+f 式9

过程2，根据式10所示的方程系数，得到式11-1所示的最大正曲率K_pos的计算公式

K_pos＝(a+b)+[(a-b)²+c²]^1/2 式11-1

或者，得到式11-2所示的最大负曲率K_neg的计算公式

K_neg＝(a+b)-[(a-b)²+c²]^1/2 式11-2。

在上述方法中，优选地，按照式12所示的公式对相干属性体E进行值域校正包括：相干属性体中采样点的值最大为1，表示完全无相关性；最小值为0，表示有极强的相关性。设定门槛值，所述门槛值为目的层裂缝发育的临界值，相干属性体中采样点的值大于或等于门槛值时，判定此处无裂缝发育；相干属性体中采样点的值小于该值时，判定此处有裂缝发育。在对相干属性体进行值域校正时，相干属性体中采样点的值小于所述门槛值时，则确定该采样点进行值域校正后的值为1，表示该样点处有裂缝发育；相干属性体中采样点的值大于或等于所述门槛值时，则确定该采样点进行值域校正后的值为0；

在式12中，E_i为相干属性体E的任意采样点；E′_i为值域校正后的相干属性体E′对应位置的采样点；E_max为门槛值。值域校正后的相干属性体E′中采样点E′_i值为1的位置表示该样点处有裂缝发育。

在上述方法中，优选地，所述门槛值E_max的取值为0≤E_max≤1。

在本发明提供的技术方案中，所述曲率属性体K可以为最大正曲率K_pos、最大负曲率为K_neg或者其它的曲率属性体。

本发明的有益效果：

1)本发明提供的技术方案创造性地将地震相干属性体和地震曲率属性体进行了融合，在裂缝空间展布位置赋予可以定量化描述裂缝密度信息的曲率数值，从而可以定量化描述地下裂缝空间展布信息；

2)本发明提供的技术方案设定了利用相干属性体预测裂缝空间展布的门槛值，从而确定裂缝空间展布的具体位置，在与曲率属性体融合过程中将由地层弯曲和褶皱引起的干扰响应有效去除。

附图说明

图1为本发明实施例提供的一种基于地震属性的裂缝定量化预测方法的流程图；

图2为Z段地层裂缝的原始地震剖面；

图3为基于本发明实施例提供的方法获得的Z段地层裂缝预测平面；

图4为基于本发明实施例提供的方法获得的Z段地层裂缝预测剖面；

图5为基于相干算法获得的Z段地层裂缝预测平面；

图6为基于相干算法获得的Z段地层裂缝预测剖面；

图7为基于曲率算法获得的Z段地层裂缝预测平面；

图8为基于曲率算法获得的Z段地层裂缝预测剖面。

具体实施方式

为了对本发明的技术特征、目的和有益效果有更加清楚的理解，现对本发明的技术方案进行以下详细说明，但不能理解为对本发明的可实施范围的限定。

以下实施例以四川盆地Z段地层裂缝为研究对象。

实施例

本实施例提供了一种基于地震属性的裂缝定量化预测的方法(流程如图1所示)，其包括以下步骤：

步骤101，获取保幅处理后的偏移叠加地震数据(基于该偏移叠加地震数据可以获得Z段地层裂缝的原始地震剖面，如图2所示)和地震工区内钻穿目的层的测井数据；井-震标定后进行地震数据目的层构造解释，其主要用于后期对目的层段开展裂缝解释提供地层层段的顶底界面。

步骤102，对三维地震数据进行去噪处理，得到去噪处理后的三维地震数据体：

去噪处理主要去除随机噪音，以提高信噪比，本实施中采用中值滤波法进行去噪处理，处理过程如下所述：

在三维地震数据中定义一个时窗，该时窗包含2j+1个样本点，其中，2j+1小于地震数据体纵向时间域样本点的个数；则该时窗中第j+1个样本点的值为：

中值滤波法：

式1中：i＝1，2，…，N；为采用中值滤波法获得的时窗中第j+1个样本点的数值；

步骤103，对去噪处理后的三维地震数据体进行相干属性计算，得到相干属性体E：

过程1，将去噪处理后的三维地震数据体作为一个分析窗口，该分析窗口中有N＝9道地震数据，K＝9个采样点，则去噪处理后的三维地震数据体可表示为式3所示的矩阵D：

式3中，d_KN为第N道的第K个样点值；

过程2，矩阵D中的第k行向量表示数据体的第k个采样点的集合，则整个去噪处理后的三维地震数据体的协方差矩阵C可表示为：

过程3，设λ_k是协方差矩阵C的第k个特征值，其中λ₁是最大的特征值，则相干属性的计算公式如式5所示，按照式5所示的计算公式可以计算得到相干属性体E：

在该步骤中，基于相干算法对Z段地层裂缝进行预测可以获得图5和图6所示的结果；其中，图5和图6分别为基于相干算法获得的Z段地层裂缝预测平面和剖面。

步骤S104，对去噪处理后的三维地震数据体进行曲率属性计算，得到曲率属性体K：

1)将去噪处理后的三维地震数据体转化为倾角数据体

过程1，计算三维地震数据体的瞬时频率

已知三维地震数据体在空间任意点的时间标量u(t,x,y)，则瞬时频率ω可以表示为

式6中：u^H(t,x,y)为所述时间标量u关于时间t的希尔伯特变换；

过程2，计算三维地震数据体的瞬时波数

三维地震数据体沿x方向和y方向的瞬时波数k_x、k_y分别为：

过程3，根据三维地震数据体的瞬时频率ω和瞬时波数k_x、k_y，计算三维地震数据体的倾角dip(p_x,q_y)

式8中：p_x、q_y分别为x方向和y方向的倾角分量。

2)根据倾角数据体计算得到数据体内任意反射点的曲率：

三维空间中某点的曲率值通过与其相邻道和样点的倾角值拟合出的空间曲面方程进行计算，本实施例设定x、y方向以该点为中心两边各取4道，根据式9所示的最小二乘逼近原理方程

z(x, y)＝ax²+by²+cxy+dx+ey+f 式9

对式9所示的方程两边求微分后代入倾角分量p_x、q_y，从而得到式10所示的方程系数

则对描述本地区目的层段裂缝构造最有效的最大正曲率K_pos可表示为：

K_pos＝(a+b)+[(a-b)²+c²]^1/2 式11-1。

在该步骤中，基于曲率算法对Z段地层裂缝进行预测可以获得图7和图8所示的结果；其中，图7和图8分别为基于曲率算法获得的Z段地层裂缝预测平面和剖面。

步骤S105，对相干属性体E的值进行值域校正，以确定裂缝的空间展布形态：

相干属性体E中样点的值最大为1，表示完全无相关性，最小为0，表示有极强的相关性；给定门槛值E_max＝0.99，当相干属性体E中样点的值小于E_max时，将该样点进行值域校正后的值设定为1；当相干属性体E中样点的值大于或等于E_max时，将该样点进行值域校正后的值设定为0，具体表示过程如式12所示

在式12中，E_i为相干属性体E的任意采样点，E′_i为值域校正后的相干属性体E′对应的采样点，E_max为相干属性体E采样点确定裂缝是否发育的门槛值；

值域校正后的相干属性体E′中样点值为1的位置表示该样点处有裂缝发育。

步骤S106，将值域校正后的相干属性体E′与曲率属性体K进行属性融合得到裂缝预测属性体Z(如式13所示)，该裂缝预测属性体Z具有预测灵敏，能够实现定量化预测的优点

Z_i＝E_i′×Q_i 式13

式13中，i＝1，2，…，M；Z_i为裂缝预测属性体Z的任意采样点的值，E′_i为值域校正后相干属性体E′对应位置的采样点的值，Q_i为最大正曲率属性体对应位置的采样点的值。

步骤S107，利用最终得到的裂缝预测属性体Z和目的层的构造解释对地下裂缝的空间展布信息进行预测，获得预测结果(如图3和图4所示)，从图中可以看出，四川盆地Z段地层岩性为灰岩，脆性大易于产生裂缝。

对图3-图8进行对比分析，可以发现，基于相干算法得到的Z段地层裂缝预测结果只能定性预测裂缝空间展布位置，基于曲率算法得到的Z段地层裂缝预测结果通过数值相对大小可以定量化预测裂缝空间展布信息，但是图中还包含一些非裂缝的干扰信息，而基于本发明实施例的方法得到的Z段地层裂缝预测结果不仅定量化预测裂缝空间展布信息，而且非裂缝干扰信息基本被消除(去除了由于地层弯曲和褶皱引起的干扰响应)。

综上所述，本发明提供的一种基于地震属性的裂缝定量化预测方法通过将可以较好预测裂缝空间展布形态的地震相干属性体和能够定量化预测裂缝密度信息的曲率属性体进行融合，从而可以利用叠后地震数据快速准确地定量化描述地下裂缝空间展布信息，为储层评价及油气勘探提供技术支撑。

Claims (9)

1.一种基于地震属性的裂缝定量化预测的方法，其包括：

获取目的层的三维地震数据和测井数据，基于这两项数据获得目的层的构造解释；

对三维地震数据进行去噪处理，得到去噪处理后的三维地震数据体；

对去噪处理后的三维地震数据体进行相干属性计算，得到相干属性体；

对去噪处理后的三维地震数据体进行曲率属性计算，得到曲率属性体；

对相干属性体进行值域校正，得到校正后的相干属性体；

基于曲率属性体以及校正后的相干属性体，计算得到裂缝预测属性体；

根据裂缝预测属性体和目的层的构造解释对地下裂缝进行预测；

对相干属性体进行值域校正包括：

设定门槛值，相干属性体中采样点的值小于所述门槛值时，则确定该采样点进行值域校正后的值为1；相干属性体中采样点的值大于或等于所述门槛值时，则确定该采样点进行值域校正后的值为0；

所述门槛值的取值为0≤所述门槛值≤1；

所述裂缝预测属性体的计算公式如式13所示

Z_i＝E_i′×Q_i 式13

在式13中，i＝1，2，…，M；Z_i为裂缝预测属性体中的采样点的值；E_i′为校正后的相干属性体中对应位置的采样点的值；Q_i为曲率属性体中对应位置的采样点的值。

2.根据权利要求1所述的方法，其中，对三维地震数据进行去噪处理时，所采用的去噪方法包括中值滤波法或均值滤波法。

3.根据权利要求1所述的方法，其中，对去噪处理后的三维地震数据体进行相干属性计算，得到相干属性体包括以下过程：

过程1，将去噪处理后的三维地震数据体作为一个分析窗口，该分析窗口中有N道地震数据，K个采样点，则所述去噪处理后的三维地震数据体表示为式3所示的矩阵D：

在式3中，d_KN为第N道的第K个样点值；

过程2，按式4所示的公式计算上述去噪处理后的三维地震数据体的协方差矩阵C

在式4中，为矩阵D中的第k行向量

过程3，确定相干属性的计算公式，由此得到相干属性体。

4.根据权利要求3所述的方法，其中，在过程3中，所述相干属性的计算公式如式5所示

在式5中，λ_k是协方差矩阵C的第k个特征值，λ_m是其中最大的特征值。

5.根据权利要求1所述的方法，其中，对去噪处理后的三维地震数据体进行曲率属性计算，得到曲率属性体包括：

将去噪处理后的三维地震数据体转化为倾角数据体；

根据倾角数据体计算得到数据体内任意反射点的曲率。

6.根据权利要求5所述的方法，其中，将去噪处理后的三维地震数据体转化为倾角数据体包括以下过程：

过程1，计算三维地震数据体的瞬时频率；

过程2，计算三维地震数据体的瞬时波数；

过程3，根据得到的三维数据体的瞬时频率和瞬时波数，计算三维地震数据体的倾角。

7.根据权利要求6所述的方法，其中，

所述过程1包括：三维地震数据体在空间任意点是时间标量u(t,x,y)，根据式6所述的计算公式，计算得到瞬时频率ω

式6中，u^H(t,x,y)为所述时间标量u关于时间t的希尔伯特变换；

所述过程2包括：按式7所示的公式，分别计算三维地震数据体沿x方向和y方向的瞬时波数

式7中，k_x为x方向的瞬时波数，k_y为y方向的瞬时波数；

所述过程3包括：根据得到的三维数据体的瞬时频率和瞬时波数，利用式8所示的公式计算三维地震数据体的倾角：

式8中，p_x、q_y分别为x方向和y方向的倾角分量。

8.根据权利要求5或6所述的方法，其中，所述根据倾角数据体计算得到数据体内任意反射点的曲率包括：

过程1，根据式9所示的最小二乘逼近原理，对方程两边求微分后，代入倾角分量p_x、q_y，得到式10所示的方程系数：

z(x,y)＝ax²+by²+cxy+dx+ey+f 式9

过程2，根据式10所示的方程系数，得到式11-1所示的最大正曲率K_pos的计算公式

K_pos＝(a+b)+[(a-b)²+c²]^1/2 式11-1

或者，得到式11-2所示的最大负曲率K_neg的计算公式

K_neg＝(a+b)-[(a-b)²+c²]^1/2 式11-2。

9.根据权利要求1所述的方法，其中，所述三维地震数据为保幅处理后的偏移叠加地震数据，所述测井数据为地震工区内钻穿目的层的测井数据。