一种砂体不连续界线方向自适应识别方法与流程

本发明涉及地震属性分析领域，特别是一种砂体不连续界线方向自适应识别方法。

背景技术：

在油气勘探开发中由于储层非均质性的影响，储层内部渗流单元的划分将变得十分重要。对于河流相砂体储层而言，其储量单元内部为不同尺度、级次的砂体复合体，这些相邻复合体之间的界面对于流体渗流具有不同程度的遮挡、阻流作用，即“渗流屏障”。因此，对储层非均质性的研究可以认识到储层中油、气、水的分布规律，通过地震属性分析来刻画储层中影响流体流动的不连续界线，并综合地质、测井及岩石物理等资料作为储层预测的约束条件，对油气储层的几何特征、地质特征及储层的岩性、形态、物性及含油气性等进行预测研究，将成为一种提高油气采收率、改善剩余油挖潜及注采井网部署的有效技术手段。目前，为了提高储层预测的精确度，特别是预测岩性界线、断层和砂体不连续界线等，各种数学方法被逐渐引入地震属性分析技术中，以期望识别更小尺度的界线信息来得到更精确的储层预测结果。

然而，在实际工作中，地震数据或多或少的包含有各种线性的、非线性的噪音，单纯的线性方法难以有效的在保证一定分辨率的前提下去除这些干扰，进而对边缘拾取产生一定的影响，导致储层预测难以获得最准确的预测结果。常规的压制噪声方法包括：均值滤波、中值滤波、截断均值滤波、多步中值滤波和高斯滤波等，均存在信噪比与分辨率之间的矛盾，进而对后续的边缘检测产生重要影响，甚至可能因为噪声的干扰产生伪边缘。

目前主要的边缘检测方法包括roberts算子、sobel算子、prewitt算子、log算子和canny算子等，这些方法都是线性的微分方法，通过高阶导数来放大数据之间差异，进而达到预测边缘的目的。然而，这些方法无可避免的一个缺点是在放大数据之间差异的同时也放大了干扰信息，因此对噪声敏感。尤其在地震属性分析领域中，因为地震属性包含了较多的干扰信息，从而使得应用效果不好，无法达到目前所要求的储层预测精度。

近年来，数学形态学作为一种非线性的新方法，在压制噪声的过程中，提高信噪比的同时能够最大限度的保证分辨率，具有较好的边缘保持能力；在边缘检测的过程中，能够得到更加连续、具有更多方向信息、并且具有一定增强效果的边缘预测结果。1964年，j.serra从理论上首次介绍了形态学的数学表达式，以此奠定了数学形态学的理论基础。将数学形态学引入地震属性分析技术中作为一种工具，可以从地震属性中提取描绘砂体展布形状的有用分量，如岩性边缘、断层和砂体不连续界线等，以此达到精确预测储层的目的。数学形态学的基本思想是通过一个具有特定形状和尺度的结构元素来扫描地震属性并与地震属性产生作用。数学形态学建立在集合论的基础之上，其基本操作包括膨胀运算、腐蚀运算、开运算和闭运算。以这些基本操作为基础，通过一定的运算组合方式，就构成了一种复杂的预测砂体不连续界线的方向自适应识别方法。

技术实现要素：

为解决现有技术中存在的问题，本发明提供了一种砂体不连续界线方向自适应识别方法，解决了以往识别方法的伪边界拾取率高，预测的真实结果可信度低和预测精度低的问题。

本发明采用的技术方案是，一种砂体不连续界线方向自适应识别方法，方法包括如下步骤：

步骤s1：输入地震属性并以数学形态学方法拾取界线信息；

步骤s2：根据马氏距离准则，分别计算结构元素对应敏感方向界线的权系数；

步骤s3：根据权系数对结构元素所对应的界线进行加权融合。

优选地，结构元素包括3×3尺度的十字形、矩形、0°方向、45°方向、90°方向和135°方向的结构元素，se1为3×3尺度的十字形，se2为矩形，se31为0°方向结构元素，se32为45°方向结构元素，se33为90°方向结构元素，se34为135°方向结构元素。

优选地，步骤s1包括如下子步骤：

步骤s11：令结构元素“se31”、“se32”、“se33”、“se34”分别对应0°、45°、90°、135°四个方向的结构元素；

步骤s12：令输入的地震属性表示为f，则fθse表示为当结构元素se的原点位于(x,y)处时，用结构元素se在(x,y)处对待处理地震属性f的腐蚀运算，表示为：

表示为当结构元素se的原点位于(x,y)处时，用结构元素se在(x,y)处对待处理地震属性f的膨胀运算，表示为：

步骤s13：数学形态学方法拾取界线信息的表达式为：

其中，k＝1，2，3，4。

优选地，步骤s2包括以下子步骤：

步骤s21：根据拾取的边缘信息，计算灰度距离，其表达式为：

d0°＝d3+d4+d5+d7+d8+d9

d45°＝d2+d4+d5+d6+d8+d9

d90°＝d2+d3+d5+d6+d7+d9

d135°＝d2+d3+d4+d6+d7+d8

其中，d0°、d45°、d90°、d135°分别对应该边界点0°、45°、90°、135°四个方向的灰度距离，dk表示该点的灰度值，k＝2，3，4，5，6，7，8，9。

步骤s22：根据边界点0°、45°、90°、135°四个方向的灰度距离，计算出四个方向的结构元素的权值，

d＝d0°+d45°+d90°+d135°

其中，ω1、ω2、ω3、ω4分别表示0°、45°、90°、135°四个方向的结构元素的权值，d表示该点全方位马氏距离和。

优选地，步骤s3的加权融合的表达式为：

式中，ωk表示结构元素的权值，ek表示边缘拾取方法，表示膨胀运算，θ表示腐蚀运算。

本发明砂体不连续界线方向自适应识别方法的有益效果如下：

1.本发明提出的数学形态学砂体不连续界线检测方法，具有边缘保持能力的滤波效果，降低了伪边界的拾取率，提高了真实结果的可信度。

2.提出的数学形态学砂体不连续界线检测方法，具有四个方向的结构元素，每一个方向的结构元素对相对应方向的边缘敏感，通过马氏距离准则进行融合，融合结果保持了强边界的显示效果，并且对弱边界具有明显的增强效果，对砂体的岩性界线、断层以及砂体不连续界线的预测和识别有非常明显的优势。可与现有的地震属性分析技术结合使用，进一步提高预测精度。

附图说明

图1为本发明砂体不连续界线方向自适应识别方法的总流程框图。

图2为本发明砂体不连续界线方向自适应识别方法的结构元素示意图。

图3为本发明砂体不连续界线方向自适应识别方法的马氏距离计算模板示意图。

图4为本发明砂体不连续界线方向自适应识别方法的灰度腐蚀算法图。

图5为本发明砂体不连续界线方向自适应识别方法的灰度膨胀运算图。

具体实施方式

下面对本发明的具体实施方式进行描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。

如图1所示，一种砂体不连续界线方向自适应识别方法，方法包括如下步骤：

步骤s1：输入地震属性并以数学形态学方法拾取界线信息；

步骤s2：根据马氏距离准则，分别计算结构元素对应敏感方向界线的权系数；

步骤s3：根据权系数对结构元素所对应的界线进行加权融合。

本实施方案的结构元素包括3×3尺度的十字形、矩形、0°方向、45°方向、90°方向和135°方向的结构元素，se1为3×3尺度的十字形，se2为矩形，se31为0°方向结构元素，se32为45°方向结构元素，se33为90°方向结构元素，se34为135°方向结构元素

本实施方案的步骤s1包括如下子步骤：

步骤s11：令结构元素“se31”、“se32”、“se33”、“se34”分别对应0°、45°、90°、135°四个方向的结构元素；

步骤s12：令输入的地震属性表示为f，则fθse表示为当结构元素se的原点位于(x,y)处时，用结构元素se在(x,y)处对待处理地震属性f的腐蚀运算，表示为：

表示为当结构元素se的原点位于(x,y)处时，用结构元素se在(x,y)处对待处理地震属性f的膨胀运算，表示为：

步骤s13：数学形态学方法拾取界线信息的表达式为：

其中，k＝1，2，3，4。

本实施方案的步骤s2包括以下子步骤：

步骤s21：根据拾取的边缘信息，计算灰度距离，其表达式为：

d0°＝d3+d4+d5+d7+d8+d9

d45°＝d2+d4+d5+d6+d8+d9

d90°＝d2+d3+d5+d6+d7+d9

d135°＝d2+d3+d4+d6+d7+d8

其中，d0°、d45°、d90°、d135°分别对应该边界点0°、45°、90°、135°四个方向的灰度距离，dk表示该店的灰度值，k＝2，3，4，5，6，7，8，9。。

步骤s22：根据边界点0°、45°、90°、135°四个方向的灰度距离，计算出四个方向的结构元素的权值，

d＝d0°+d45°+d90°+d135°

其中，ω1、ω2、ω3、ω4分别表示0°、45°、90°、135°四个方向的结构元素的权值，d表示该点全方位马氏距离和。

本实施方案的步骤s3的加权融合的表达式为：

式中，ωk表示结构元素的权值，ek表示边缘拾取方法，表示膨胀运算，θ表示腐蚀运算。

本实施方案在实施时，步骤一：输入地震属性，并以数学形态学方法拾取界线信息，方法如下：分别使用以结构元素“se31”、“se32”、“se33”、“se34”为核心的界线检测算子对地震属性拾取界线；

方法如下：以结构元素“se1”对地震属性进行有顺序的腐蚀运算和膨胀运算，再以结构元素“se2”对得到的结果进行有顺序的膨胀运算和腐蚀运算；

在实整数集合z中，如果f(x,y)和b(x,y)这两个函数均为离散函数，则坐标(x,y)是来自笛卡尔乘积z²的一个整数对，f和b是对每个(x,y)坐标赋以灰度值(来自实整数集合r的一个实数)的函数。如果灰度级也是整数，则用z代替r。

灰度腐蚀算法，当b的原点(中心点)位于(x,y)处时，用结构元素b在(x,y)处对待处理图像f的腐蚀运算定义为待处理图像f中与结构元素b重合区域的最小值。腐蚀运算的运算符表示为“θ”，利用结构元素b腐蚀待处理图像f可表示为fθb，即:

该式指出，将结构元素的原点与待处理图像中每一个像素点的位置重合，在任意位置的腐蚀运算由包含在与b重合区域中的f的所有值中的最小值确定。

图4所示为灰度形态学腐蚀运算的一个具体说明。图中，左侧数值阵列表示一个大小为6×6个像素的待处理图像f，右侧数值阵列为腐蚀结果。结构元素b为3×3大小的窗口。使用结构元素b对待处理图像f进行腐蚀运算，表示为b的原点(中心点)与图像中每一个像素点重合时，求取窗口包含区域的待处理图像中的最小值，并将该值赋给结构元素的原点。以此方式，扫描待处理图像中的每一个像素点，便可得到最后的腐蚀结果。

当b的原点(中心点)位于(x,y)处时，用结构元素b在(x,y)处对待处理图像f的膨胀运算定义为待处理图像f中与结构元素b重合区域的最大值。膨胀运算的运算符表示为利用结构元素b膨胀灰度图灰度形态学膨胀运算像f可表示为即：

该式指出，将结构元素的原点与待处理图像中每一个像素点的位置重合，在任意位置的膨胀运算由包含在与b重合区域中的f的所有值中的最大值确定。

图5所示为灰度形态学膨胀运算的一个具体说明，图中，左侧数值阵列表示一个大小为6×6个像素的待处理图像f，右侧数值阵列为膨胀结果。结构元素b为3×3大小的窗口。使用结构元素b对待处理图像f进行膨胀运算，可以表示为b的原点(中心点)与图像中某一个像素点重合时，求取窗口包含区域的待处理图像的最大值，并将该值赋给结构元素的原点。以此方式，扫描待处理图像中的每一个像素点，便可得到最后的膨胀结果。

通过对数字图像的表示方法、集合理论、二值形态学以及灰度形态学的研究，我们介绍了数学形态学的两个基本运算：膨胀运算和腐蚀运算。基于这两种基本运算，形成了以开运算、闭运算、开-闭运算、闭-开运算等复杂的运算组合方式为特点的图像和信号处理技术。数学形态学方法可以提取和预测图像中某种具有特殊形状的表达，如边界、骨架等。

其中，结构元素“se1”为十字形结构元素，结构元素“se2”为矩形结构元素，如图2所示，其中，1表示参与计算数据的位置，0表示不参与计算数据的位置。令输入的地震属性表示为f，则fθse表示为当结构元素se的原点(中心点)位于(x,y)处时，用结构元素se在(x,y)处对待处理地震属性f的腐蚀运算，表示为：

表示为当结构元素se的原点(中心点)位于(x,y)处时，用结构元素se在(x,y)处对待处理地震属性f的膨胀运算，表示为：

最后，界线拾取方法表示为：

其中，k＝1，2，3，4。

步骤二：根据马氏距离准则，分别计算结构元素“se31”、“se32”、“se33”、“se34”对应敏感方向界线的权系数；

方法如下：在图3所示的模板中，灰度距离计算如下：

d0°＝d3+d4+d5+d7+d8+d9

d45°＝d2+d4+d5+d6+d8+d9

d90°＝d2+d3+d5+d6+d7+d9

d135°＝d2+d3+d4+d6+d7+d8

其中，d0°、d45°、d90°、d135°分别对应该边界点0°、45°、90°、135°四个方向的灰度距离，dk表示该点的灰度值，k＝2，3，4，5，6，7，8，9。

令ω1、ω2、ω3、ω4分别表示0°、45°、90°、135°四个方向的结构元素的权值，d表示该点全方位马氏距离和，则有：

d＝d0°+d45°+d90°+d135°

步骤三：按照所求得的权系数对结构元素“se31”、“se32”、“se33”、“se34”所对应的界线进行加权融合；

方法如下：

式中，ωk表示结构元素的权值，ek表示边缘拾取方法，表示膨胀运算，θ表示腐蚀运算。