CN104422959A - 一种检测储层边界的曲率属性方法 - Google Patents

Abstract

本发明提供了一种检测储层边界的曲率属性方法，属于油田勘探领域。本方法包括：(1)对叠后地震资料进行常规处理，得到叠后地震数据体的主测线与联络测线排列图；(2)在所述叠后地震数据体的主测线与联络测线排列图中，以一个地震数据点为中心圆点，获取该中心圆点的沿层振幅；(3)计算该中心圆点的导数值；(4)计算该中心圆点的振幅曲率；(5)依次对叠后地震数据体中的每一个地震数据点重复步骤(2)至步骤(4)，得到所有地震数据点的振幅曲率，所有地震数据点的振幅曲率即构成用于储层边界检测的曲率属性体。利用本发明可以获取更加精确的边界表述特征，因此可以获得更好的复杂边界储层的预测效果。

Description

一种检测储层边界的曲率属性方法

技术领域

本发明属于油田勘探领域，具体涉及一种检测储层边界的曲率属性方法，针对实际油气田勘探中存在的储层边界问题，利用高精度叠后地震资料主测线(联络测线)的振幅导数结果计算二次曲率属性。

背景技术

曲率属性作为地震属性中的一种，已经为研究碳酸盐岩地区缝洞型储层密集发区带及岩溶储层发育的最有效方法。传统的曲率计算方法多用于描述复杂储层的连续性以及延伸性。但在很多有效区带内圈定的曲率属性带有成片的模糊特征，不能精确描述缝洞型储层的边界特征。也就是说，一般的曲率属性都是描述断裂或者裂缝储层的局部突变点，展现的是储层在整个三维空间的内部连续性，而很少被用于来描述储层的边界特征。

发明内容

本发明的目的在于解决上述现有技术中存在的难题，提供一种检测储层边界的曲率属性方法，在高精度叠后地震资料纵横测线的振幅导数的基础上，开展振幅导数曲率的计算方法，可以发掘叠后地震资料中的曲率属性的边界特征，挖掘出比传统的曲率算法更为准确的储层边界特征，因此可以获得更好储层预测效果。

本发明是通过以下技术方案实现的：

一种检测储层边界的曲率属性方法，所述方法包括以下步骤：

(1)对叠后地震资料进行常规处理，得到叠后地震数据体的主测线与联络测线排列图；

(2)在所述叠后地震数据体的主测线与联络测线排列图中，以一个地震数据点为中心圆点，获取该中心圆点的沿层振幅；

(3)计算该中心圆点的导数值；

(4)计算该中心圆点的振幅曲率；

(5)依次对叠后地震数据体中的每一个地震数据点重复步骤(2)至步骤(4)，得到所有地震数据点的振幅曲率，所有地震数据点的振幅曲率即构成用于储层边界检测的曲率属性体。

所述步骤(2)是这样实现的：

在叠后地震数据体的主测线与联络测线排列图中，假定该中心圆点的主测线的序号为I，联络测线的序号为J，则该中心圆点的振幅值为A0；

中心圆点的右边相邻主测线上的点，即主测线的序号为I，联络测线的序号为J+1的点的振幅值为A1，主测线的序号为I，联络测线的序号为J+2的振幅值为A2；

中心点下方沿着联络测线方向的点，即主测线的序号为I+1，联络测线的序号为J的点的振幅值为A3，主测线的序号为I+2，联络测线的序号为J的点的振幅值为A4；(这5个振幅值都是从叠后地震数据中直接读出来的)

从叠后地震数据体中读取线间距和道间距；所述线间距是两条相邻主测线之间的间距；道间距是两条相邻联络测线之间的间距；假设线间距和道间距均为M米。

所述步骤(3)是这样实现的：

用a表示中心圆点的振幅沿主测线方向变化的二阶导数，b表示中心圆点的振幅沿联络测线方向变化的二阶导数，c表示中心圆点的振幅沿主测线、联络测线变化的偏导数，d表示中心圆点的振幅沿主测线方向变化的一阶导数，e表示中心圆点的振幅沿联络测线变化的一阶导数，f和g为求取二阶导数过程中使用的中间结果，计算公式如下：

$d=\frac{A1-A0}{M}$

$e=\frac{A3-A0}{M}$

$f=\frac{A2-A1}{M}$

$g=\frac{A4-A3}{M}$

$a=\frac{f-d}{M}$

$b=\frac{g-e}{M}$

$c=\frac{e-d}{M}.$

所述步骤(4)是这样实现的：

采用以下公式计算该中心圆点的振幅曲率：

高斯曲率Kg=4ab-c*c

最大正曲率Kpos=(a+b)+sqrt((a-b)*(a-b)+c*c)

最小负曲率Kneg=(a+b)-sqrt((a-b)*(a-b)+c*c)

倾角曲率Kdip=2.0*(a*d*d+b*e*e+c*d*e)/((d*d+e*e)*dt+1.e-25)。

与现有技术相比，本发明的有益效果是：本发明对高精度叠后地震数据的振幅进行了二次改造，形成针对不同测线方向的导数值，然后用该导数值替换原有的地震振幅数据，来计算该振幅点的曲率属性，该振幅属性可以发掘叠后地震资料中的边界特征、计算的振幅曲率相比传统的曲率算法可以获取更加精确的边界表述特征，因此可以获得更好的复杂边界储层的预测效果。

附图说明

图1是不同曲率定义示意图。

图2是宏面元中心点以及相邻点地震道示意图。

图3是本发明方法的步骤框图。

图4-1是实施例中的最大曲率。

图4-2是实施例中的高斯曲率。

图4-3是实施例中的最大正曲率。

图4-4是实施例中的最小负曲率。

图4-5是实施例中的倾角曲率。

具体实施方式

下面结合附图对本发明作进一步详细描述：

本发明利用叠后高精度地震资料的振幅，开展不同测线方向上的导数计算，得到振幅导数的结果；在该振幅导数的结果上开展曲率计算，可以更好的描述不同类型储层的边界特征。

(1)导数公式及曲率的定义

①一阶导数

设y=f(x)在点x₀处的某个邻域有定义，当自变量x在点x₀处有增量Δx时，函数有相应的增量Δy，

当Δx→0时，的极限存在，称这个极限为函数y=f(x)在点x₀处的导数，记作f′(x₀)或即 $f^{\′}\left(x_0\right)=y^{\′}{|}_{x=x_0}=\underset{\mathrm{\Δx}\→0}{\lim }\frac{\mathrm{\Δy}}{\mathrm{\Δx}}=\underset{\mathrm{\Δx}\→0}{\lim }\frac{f(x_0+\mathrm{\Δx})-f\left(x_0\right)}{\mathrm{\Δx}}---\left(1\right)$

②二阶导数

函数y=f(x)的导数y′=f′(x)仍是x的可导函数时，则称y′=f′(x)的导数为做函数y=f(x)的二阶导数，记作

y″，f"(x)或 $\frac{d^{2}y}{d{x}^2},---\left(2\right)$

即y″=(y′)′，f"(x)=[f′(x)]′ $\frac{d^{2}y}{d{x}^2}=\frac{d}{\mathrm{dx}}\left(\frac{\mathrm{dy}}{\mathrm{dx}}\right).$ 或 $\frac{d^{2}f\left(x\right)}{d{x}^2}=\frac{d}{\mathrm{dx}}\left(\frac{\mathrm{df}\left(x\right)}{\mathrm{dx}}\right)---\left(3\right)$

相应地，把y=f(x)的导数f′(x)叫做函数y=f(x)的一阶导数.类似地，二阶导数的导数，叫做三阶导数，三阶导数的导数叫做四阶导数，…，一般地，(n-1)阶导数的导数叫做n阶导数，分别记作

y"′，y⁽⁴⁾，…，y⁽ⁿ⁾或 $\frac{d^{3}y}{d{x}^3},\frac{d^{4}y}{d{x}^4},\·\·\·,\frac{d^{n}y}{d{x}^n}.---\left(4\right)$

③偏导数

设函数z=f(x，y)在点(x₀，y₀)的某一邻域内有定义，当y固定在y₀而x在x₀处有增量Δx时，相应地函数有增量

f(x₀+Δx,y₀)-f(x₀，y₀).       (5)

如果极限

$\underset{\mathrm{\Δx}\→0}{\lim }\frac{f(x_0+\mathrm{\Δx},y_0)-f(x_0,y_0)}{\mathrm{\Δx}}---\left(6\right)$

存在，则称此极限为函数z=f(x，y)在点(x₀，y₀)处对x的偏导数，记作 $\frac{\∂z}{\∂x}{|}_{\begin{array}{c}x=x_0\\ y=y_0\end{array}},\frac{\∂f}{\∂x}{|}_{\begin{array}{c}x=x_0\\ y=y_0\end{array}}{,z}_x{|}_{\begin{array}{c}x=x_0\\ y=y_0\end{array}},$ 或f_x（x₀,y₀)    .(7)

④曲率的含义

对于三维空间的反射面形态来讲，不同切平面与曲面相交就可能在计算点形成无数条交线(如图2所示)，所构成的曲率也就有无限多个。数学上已证明法平面与曲面相交之交线(曲线)上计算出的曲率子集是最为有效的，并且其中两两正交的一组曲率属性能够最佳反映曲面的起伏形态：如极大/极小曲率(Kmax/Kmin)、最大正/最小负曲率、走向/趋向曲率(Ks/Kd)，并利用这些正交曲率分量的组合来反映曲面的变形形态、变形强度及变形大小，如平均曲率、高斯曲率、弯曲度，以及反映曲面三维起伏形态——穹隆、脊梁、沟谷、地堑的形态指数等参数。利用这组不同的曲率属性就可以实现对三维起伏面变形的详细描述。图1为曲率属性的定义图，用来说明不同的曲率，比如图中Ks，Kd以及Kmax，Kmin之间的关系。

⑤高精度曲率体的计算方法

反射为倾平界面下，反射面可表示为：

T(x，y)=Dx+Ey+F

反射为曲界面下，反射面可表示为：

T(x，y)=Ax²+Bxy+Cy²+Dx+-Ey+F      (8)

如将分析点一直选为坐标原点x=y=t=0处，就有如下的几何属性计算结果。

曲率属性

沿x，y方向的曲率正交，均为法曲率：

Kx=2A/(1+D²)^3/2    (9)

Ky=2C/(1+E²)^3/2    (10)

由于扫描是在各个方向进行的，故可以获得一簇法曲率。依据各种层面曲率属性的定义，可以获得反映不同特性的曲率属性体(附图2)。

⑥曲率属性系列的计算及地质意义

由上方法原理可见，由于曲率属性为二维属性，两两正交的一组法曲率属性能够最佳反映曲面的起伏形态及变形强度，故可从中分离出一组不同的曲率属性系列：

a.极大/极小曲率(Kmax/Kmin)：不同方向曲率计算中出现的两个正交的最大及最小的曲率分量，也被称为主曲率，它是计算其它曲率的基础。

该组曲率属性对界定断层及断层的几何形态、描述断层的大小非常有效。以这个属性表示的断层表现为正曲率值和负曲率值彼此相邻。另外，曲率值也确定了断层的错断方向，正的曲率值代表上升盘，负的曲率值代表下降盘。两者强度愈大，断裂规模愈大。正是这种相反的曲率值使得判定断层的断错方向成为一种有用的工具。

b.平均曲率Km：极大曲率与极小曲率的平均；反映各个方向的平均变形强度；

Km=(Kmax+Kmin)/2=(Kx+Ky)/2   (11)

c.高斯曲率Kg：高斯曲率(Gauss，1827)也称为全曲率。定义为主曲率的积：

Kg=Kmax*Kmin     (12)

该定理表明层面的等距弯曲不会改变层面上任一点的高斯曲率。简言之如果层面以某种形式褶皱了，但没有断开、拉伸或压缩，，那么高斯曲率仍保持为常量。平面在各个方向上的曲率都是零；而一个柱体或锥体，在极小曲率方向的曲率为零，在极大曲率方向上的曲率为非零，故曲率仍为零。因此很多形态不能单独用高斯曲率加以区分，还需要平均曲率信息加以辅助。

d.最大正/最小负曲率(Kpos/Kneg)：不同方向曲率计算中出现的正值之最大及负值之最小的两个正交的曲率分量，分离突起与凹陷部位，分别描述发生弯曲的强度与宽窄大小；Kpos定义为地形突出部分(相对地下埋深方向来讲，与地貌学中的方向正好相反，即在地貌学中Kpos与Kneg互换)，为突起(背斜)顶部曲率图，与应力、应变量成正比；Kneg为凹陷(向斜)曲率图，也与应力、应变量成正比；用这种方式引入曲率不仅放大了层面中的断层及裂缝，也放大了一些小的线性构造，能得到地质体的边缘显示，但不包含形态的信息。解释中需根据实际地质情况判别裂缝主要发育于背斜顶部与向斜坳部、或两者都发育来妁情选取；

(2)技术实现步骤

技术实现流程如图3所示，步骤如下：

①叠后地震资料的常规处理，处理后得到的还是叠后地震资料。不过是更能反映地层构造特征的高精度叠后地震数据资料。

叠后地震资料的常规处理资料的预处理主要针对前期的常规地震保幅处理工作，主要包括道编辑、带通滤波、真振幅恢复、静校正、速度分析、剩余静校正、地表振幅一致性补偿、反褶积及动校正等。

②叠后地震数据的沿层振幅

如图2所示的叠后地震数据主测线与联络测线排列的示意图(一般把沿着采集地震数据主方向的测线叫做主测线，垂直主测线的叫做联络测线。步骤1中得到的叠后地震数据当沿着主测线和联络测线方向可以被称作不同测线方向的数据，实际上还是一个数据，而且，对于本方法来说，由于主测线和联络测线的地震数据具有对称性，所以主测线和联络测线的不同划分不会影响最终的计算结果)。图中圆点为某一个时间层位上的地震数据点，假定图中中心圆点的线号(即主测线的序号)为I，道号(即联络测线的序号)为J(线号和道号可以分别从该点的叠后地震数据的道头中获取)该点的振幅值为A0；中心点右边相邻主测线的点(各相邻主测线之间的间距是相等的，各相邻联络测线之间的间距也是相等的)，其振幅依次为A1(线号为I，道号为J+1)和A2(线号为I，道号为J+2)；中心点下方沿着联络测线方向的点，其振幅分别为A3(线号为I+1，道号为J)和A4(线号为I+2，道号为J)。叠后地震数据的排量都是循序渐进的，依次加1。

这5个振幅组成了中心圆点计算振幅曲率所需要的原始地震数据点。同时从叠后地震数据中可以读取到的是相邻主测线和联络测线的间距。(主测线和联络测线的间距简称线间距和道间距(即两条相邻主测线之间的间距为线间距，两条相邻联络测线之间的间距为道间距)，一般情况下，线间距和道间距是相等的，这里假设线间距和道间距均为M米。

③叠后地震数据点的导数计算

为了简洁明了，用a表示振幅沿主测线方向变化的二阶导数，b表示振幅沿联络测线方向变化的二阶导数，c表示沿主、联络测线变化的偏导数，d表示振幅沿主测线方向变化的一阶导数，e表示沿联络测线变化的一阶导数，f和g为求取二阶导数过程中使用的中间结果。

a，b，c，d，e的计算公式方法为：

$d=\frac{A1-A0}{M}$

$e=\frac{A3-A0}{M}$

$f=\frac{A2-A1}{M}$

$g=\frac{A4-A3}{M}$

$a=\frac{f-d}{M}$

$b=\frac{g-e}{M}$

$c=\frac{e-d}{M}$

上面公式中的M为线间距和道间距。

④叠后地震数据的振幅曲率

在a，b，c，d，e几个参数确定后，改造后的用于检测储层边缘的几个曲率属性可以采用以下公式计算：

最大曲率Kmax=Max(Kx，Ky)，其中Kx和Ky分别代表x方向和y方向的曲率值，Kmax取这两个值中最大的那个值；

高斯曲率Kg=4ab-c*c

最大正曲率Kpos=(a+b)+sqrt((a-b)*(a-b)+c*c)

最小负曲率Kneg=(a+b)-sqrt((a-b)*(a-b)+c*c)

倾角曲率Kdip=2.0*(a*d*d+b*e*e+c*d*e)/((d*d+e*e)*dt+1.e-25)

⑤用于储层边缘检测的曲率属性体

通过步骤2到步骤4，，可以得到叠后地震数据的中心点在该层位上的振幅曲率值，依次以各个点作为中心点，然后重复上面的步骤，就可以得到整个层位面上的曲率值，然后再计算各个时间层位上的值，就可以得到整个叠后地震数据的三维振幅曲率属性。必须指出的是：由于面元以及采集系统的限制，为了保证中心点周边面元地震道的个数足够，通常不计算地震工区边上的几条线。同时由于时窗的限制，在整个地震道上进行计算时，也不计算地震道整个采样时间上的前面半个时窗以及后半个时窗内的数据。

图4-1至图4-5为塔河油田某区振幅曲率属性图，五种曲率中，振幅最大曲率突出的是最大变化量，着重刻画裂缝以及缝洞发育的边界；振幅高斯曲率突出的是沿圆弧面的变化特征，对于预测溶洞的边界，该属性效果较好；最大正曲率和最小负曲率对应一个振幅变化的正梯度和负梯度；倾角曲率是沿最大倾角方向的振幅变化。他们分别着重刻画的是缝洞发育边界的突变。从图中的储层预测结果可以看到，各个不同属性对于复杂缝洞型储层的各个边界显示了较好的检测效果。

本发明针对实际油田勘探中面临的复杂储层边界检测问题，采用的方法是针对实际叠后地震资料开展具有振幅导数信息的三维边缘检测曲率体计算研究，对叠后中心点地震道展开展一阶、二阶以及偏导数计算、然后利用导数计算的结果开展三维曲率属性计算，并最终将该步骤推广到整个工区，得到整个地震数据的三维振幅曲率属性计算结果，进而对工区发育的储层边界进行预测。

上述技术方案只是本发明的一种实施方式，对于本领域内的技术人员而言，在本发明公开了应用方法和原理的基础上，很容易做出各种类型的改进或变形，而不仅限于本发明上述具体实施方式所描述的方法，因此前面描述的方式只是优选的，而并不具有限制性的意义。

Claims (4)

1.一种检测储层边界的曲率属性方法，其特征在于：所述方法包括以下步骤：

(1)对叠后地震资料进行常规处理，得到叠后地震数据体的主测线与联络测线排列图；

(2)在所述叠后地震数据体的主测线与联络测线排列图中，以一个地震数据点为中心圆点，获取该中心圆点的沿层振幅；

(3)计算该中心圆点的导数值；

(4)计算该中心圆点的振幅曲率；

(5)依次对叠后地震数据体中的每一个地震数据点重复步骤(2)至步骤(4)，得到所有地震数据点的振幅曲率，所有地震数据点的振幅曲率即构成用于储层边界检测的曲率属性体。

2.根据权利要求1所述的检测储层边界的曲率属性方法，其特征在于：所述步骤(2)是这样实现的：

在叠后地震数据体的主测线与联络测线排列图中，假定该中心圆点的主测线的序号为I，联络测线的序号为J，则该中心圆点的振幅值为A0；

中心圆点的右边相邻主测线上的点，即主测线的序号为I，联络测线的序号为J+1的点的振幅值为A1，主测线的序号为I，联络测线的序号为J+2的振幅值为A2；

中心点下方沿着联络测线方向的点，即主测线的序号为I+1，联络测线的序号为J的点的振幅值为A3，主测线的序号为I+2，联络测线的序号为J的点的振幅值为A4；

从叠后地震数据体中读取线间距和道间距；所述线间距是两条相邻主测线之间的间距；道间距是两条相邻联络测线之间的间距；假设线间距和道间距均为M米。

3.根据权利要求2所述的检测储层边界的曲率属性方法，其特征在于：所述步骤(3)是这样实现的：

用a表示中心圆点的振幅沿主测线方向变化的二阶导数，b表示中心圆点的振幅沿联络测线方向变化的二阶导数，c表示中心圆点的振幅沿主测线、联络测线变化的偏导数，d表示中心圆点的振幅沿主测线方向变化的一阶导数，e表示中心圆点的振幅沿联络测线变化的一阶导数，f和g为求取二阶导数过程中使用的中间结果，计算公式如下：

$d=\frac{A1-A0}{M}$

$e=\frac{A3-A0}{M}$

$f=\frac{A2-A1}{M}$

$g=\frac{A4-A3}{M}$

$a=\frac{f-d}{M}$

$b=\frac{g-e}{M}$

$c=\frac{e-d}{M}.$

4.根据权利要求3所述的检测储层边界的曲率属性方法，其特征在于：所述步骤(4)是这样实现的：

采用以下公式计算该中心圆点的振幅曲率：

高斯曲率Kg=4ab-c*c

最大正曲率Kpos=(a+b)+sqrt((a-b)*(a-b)+c*c)

最小负曲率Kneg=(a+b)-sqrt((a-b)*(a-b)+c*c)

倾角曲率Kdip=2.0*(a*d*d+b*e*e+c*d*e)/((d*d+e*e)*dt+1.e-25)。