CN109948291A - 一种砂体不连续界线方向自适应识别方法 - Google Patents
一种砂体不连续界线方向自适应识别方法 Download PDFInfo
- Publication number
- CN109948291A CN109948291A CN201910256924.0A CN201910256924A CN109948291A CN 109948291 A CN109948291 A CN 109948291A CN 201910256924 A CN201910256924 A CN 201910256924A CN 109948291 A CN109948291 A CN 109948291A
- Authority
- CN
- China
- Prior art keywords
- structural element
- boundary line
- recognition methods
- sand body
- follows
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
Landscapes
- Geophysics And Detection Of Objects (AREA)
Abstract
本发明公开了一种砂体不连续界线方向自适应识别方法，方法包括如下步骤：步骤S1：输入地震属性并以数学形态学方法拾取界线信息；步骤S2：根据马氏距离准则，分别计算结构元素对应敏感方向界线的权系数；步骤S3：根据权系数对结构元素所对应的界线进行加权融合；解决了以往识别方法的伪边界拾取率高，预测的真实结果可信度低和预测精度低的问题。
Description
技术领域
本发明涉及地震属性分析领域，特别是一种砂体不连续界线方向自适应识别方法。
背景技术
在油气勘探开发中由于储层非均质性的影响，储层内部渗流单元的划分将变得十分重要。对于河流相砂体储层而言，其储量单元内部为不同尺度、级次的砂体复合体，这些相邻复合体之间的界面对于流体渗流具有不同程度的遮挡、阻流作用，即“渗流屏障”。因此，对储层非均质性的研究可以认识到储层中油、气、水的分布规律，通过地震属性分析来刻画储层中影响流体流动的不连续界线，并综合地质、测井及岩石物理等资料作为储层预测的约束条件，对油气储层的几何特征、地质特征及储层的岩性、形态、物性及含油气性等进行预测研究，将成为一种提高油气采收率、改善剩余油挖潜及注采井网部署的有效技术手段。目前，为了提高储层预测的精确度，特别是预测岩性界线、断层和砂体不连续界线等，各种数学方法被逐渐引入地震属性分析技术中，以期望识别更小尺度的界线信息来得到更精确的储层预测结果。
然而，在实际工作中，地震数据或多或少的包含有各种线性的、非线性的噪音，单纯的线性方法难以有效的在保证一定分辨率的前提下去除这些干扰，进而对边缘拾取产生一定的影响，导致储层预测难以获得最准确的预测结果。常规的压制噪声方法包括：均值滤波、中值滤波、截断均值滤波、多步中值滤波和高斯滤波等，均存在信噪比与分辨率之间的矛盾，进而对后续的边缘检测产生重要影响，甚至可能因为噪声的干扰产生伪边缘。
目前主要的边缘检测方法包括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)是来自笛卡尔乘积Z2的一个整数对，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表示边缘拾取方法，表示膨胀运算，Θ表示腐蚀运算。
Claims (5)
1.一种砂体不连续界线方向自适应识别方法，其特征在于，方法包括如下步骤：
步骤S1：输入地震属性并以数学形态学方法拾取界线信息；
步骤S2：根据马氏距离准则，分别计算结构元素对应敏感方向界线的权系数；
步骤S3：根据权系数对结构元素所对应的界线进行加权融合。
2.根据权利要求1所述的砂体不连续界线方向自适应识别方法，其特征在于，所述结构元素包括3×3尺度的十字形、矩形、0°方向、45°方向、90°方向和135°方向的结构元素，所述SE1为3×3尺度的十字形，所述SE2为矩形，所述SE31为0°方向结构元素，所述SE32为45°方向结构元素，所述SE33为90°方向结构元素，所述SE34为135°方向结构元素。
3.根据权利要求2所述的砂体不连续界线方向自适应识别方法，其特征在于，所述步骤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。
4.根据权利要求1所述的砂体不连续界线方向自适应识别方法，其特征在于，所述步骤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表示该点全方位马氏距离和。
5.根据权利要求1所述的砂体不连续界线方向自适应识别方法，其特征在于，所述步骤S3的加权融合的表达式为：
式中，ωk表示结构元素的权值，Ek表示边缘拾取方法，表示膨胀运算，Θ表示腐蚀运算。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910256924.0A CN109948291A (zh) | 2019-04-01 | 2019-04-01 | 一种砂体不连续界线方向自适应识别方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910256924.0A CN109948291A (zh) | 2019-04-01 | 2019-04-01 | 一种砂体不连续界线方向自适应识别方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN109948291A true CN109948291A (zh) | 2019-06-28 |
Family
ID=67013258
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201910256924.0A Pending CN109948291A (zh) | 2019-04-01 | 2019-04-01 | 一种砂体不连续界线方向自适应识别方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN109948291A (zh) |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112255673A (zh) * | 2020-09-27 | 2021-01-22 | 中国石油天然气股份有限公司 | 一种基于地震反演的砂体顶界面自动追踪方法 |
| CN113362312A (zh) * | 2021-06-11 | 2021-09-07 | 中海石油(中国)有限公司 | 一种储层不连续界限方位及宽度的检测方法 |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104715472A (zh) * | 2014-04-17 | 2015-06-17 | 中国石油化工股份有限公司 | 基于数学形态学的岩性体自动追踪方法 |
| CN109254320A (zh) * | 2018-10-18 | 2019-01-22 | 中国海洋石油集团有限公司 | 基于图像处理方法的地震属性优化和砂体叠置区预测方法 |
-
2019
- 2019-04-01 CN CN201910256924.0A patent/CN109948291A/zh active Pending
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104715472A (zh) * | 2014-04-17 | 2015-06-17 | 中国石油化工股份有限公司 | 基于数学形态学的岩性体自动追踪方法 |
| CN109254320A (zh) * | 2018-10-18 | 2019-01-22 | 中国海洋石油集团有限公司 | 基于图像处理方法的地震属性优化和砂体叠置区预测方法 |
Non-Patent Citations (1)
| Title |
|---|
| 贺萌: "基于自适应形态学的边缘检测及应用", 《中国优秀硕士学位论文全文数据库·信息科技辑》 * |
Cited By (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112255673A (zh) * | 2020-09-27 | 2021-01-22 | 中国石油天然气股份有限公司 | 一种基于地震反演的砂体顶界面自动追踪方法 |
| CN112255673B (zh) * | 2020-09-27 | 2024-03-26 | 中国石油天然气股份有限公司 | 一种基于地震反演的砂体顶界面自动追踪方法 |
| CN113362312A (zh) * | 2021-06-11 | 2021-09-07 | 中海石油(中国)有限公司 | 一种储层不连续界限方位及宽度的检测方法 |
| CN113362312B (zh) * | 2021-06-11 | 2024-03-08 | 中海石油(中国)有限公司 | 一种储层不连续界限方位及宽度的检测方法 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Igbinosa | Comparison of edge detection technique in image processing techniques | |
| Shafiabadi et al. | Identification of reservoir fractures on FMI image logs using Canny and Sobel edge detection algorithms | |
| CN109948291A (zh) | 一种砂体不连续界线方向自适应识别方法 | |
| CN109254320A (zh) | 基于图像处理方法的地震属性优化和砂体叠置区预测方法 | |
| CN103901469B (zh) | 地震数据的恢复方法 | |
| Li et al. | Pixel-level detection and measurement of concrete crack using faster region-based convolutional neural network and morphological feature extraction | |
| CN111025409A (zh) | 一种水淹层评价方法、装置及存储介质 | |
| Zhen et al. | Design and evaluation of an FFT-based space-time image velocimetry (STIV) for time-averaged velocity measurement | |
| CN105353409B (zh) | 一种用于抑制全波形反演震源编码串扰噪音的方法和*** | |
| CN106772599B (zh) | 一种计算地层横波速度的方法及装置 | |
| CN111859709A (zh) | 一种含水层结构变异转移概率的地质统计模拟方法及装置 | |
| CN110716239B (zh) | 一种基于电成像的测井砂砾岩体岩性精细评价方法 | |
| CN107621655B (zh) | 基于DoS滤波的三维数据断层增强方法 | |
| CN104062681B (zh) | 一种基于分数阶导数的地震层位追踪预处理方法 | |
| CN109212609A (zh) | 基于波动方程延拓的近地表噪音压制方法 | |
| CN113900141B (zh) | 油气分布预测方法及装置 | |
| Rao et al. | A review on edge detection technique in image processing techniques | |
| Lee et al. | Total variation-based image noise reduction with generalized fidelity function | |
| CN106772596A (zh) | 一种确定叠前时间偏移速度场的方法及装置 | |
| CN110097133A (zh) | 一种用于盾构机大数据的预处理方法 | |
| CN113325474B (zh) | 生物礁判别方法 | |
| CN108318936A (zh) | 一种地层划分处理方法和装置 | |
| CN110555221A (zh) | 一种计算区域地层升降幅度的方法及装置 | |
| CN116630811B (zh) | 一种河流提取方法、装置、终端设备和可读存储介质 | |
| CN112324422B (zh) | 一种电成像测井缝洞识别方法、***及孔隙结构表征方法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| RJ01 | Rejection of invention patent application after publication | ||
| RJ01 | Rejection of invention patent application after publication |
Application publication date: 20190628 |