CN108957530B - 一种基于地震相干体切片的裂缝自动检测方法 - Google Patents
一种基于地震相干体切片的裂缝自动检测方法 Download PDFInfo
- Publication number
- CN108957530B CN108957530B CN201810499447.6A CN201810499447A CN108957530B CN 108957530 B CN108957530 B CN 108957530B CN 201810499447 A CN201810499447 A CN 201810499447A CN 108957530 B CN108957530 B CN 108957530B
- Authority
- CN
- China
- Prior art keywords
- value
- crack
- pixel
- neighborhood
- algorithm
- 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.)
- Active
Links
- 230000001133 acceleration Effects 0.000 title claims abstract description 81
- 238000012360 testing method Methods 0.000 title claims abstract description 22
- 238000001914 filtration Methods 0.000 claims abstract description 42
- 230000002146 bilateral effect Effects 0.000 claims abstract description 21
- 238000000034 method Methods 0.000 claims abstract description 20
- 238000012545 processing Methods 0.000 claims abstract description 14
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 claims abstract description 11
- 230000000877 morphologic effect Effects 0.000 claims abstract description 6
- 230000008439 repair process Effects 0.000 claims abstract description 5
- 238000000605 extraction Methods 0.000 claims description 8
- 230000003321 amplification Effects 0.000 claims description 6
- 238000011835 investigation Methods 0.000 claims description 6
- 238000002372 labelling Methods 0.000 claims description 6
- 238000010606 normalization Methods 0.000 claims description 6
- 238000003199 nucleic acid amplification method Methods 0.000 claims description 6
- 230000002708 enhancing effect Effects 0.000 claims description 5
- 208000035126 Facies Diseases 0.000 claims description 4
- 239000006002 Pepper Substances 0.000 claims description 4
- 230000007812 deficiency Effects 0.000 claims description 4
- 238000006243 chemical reaction Methods 0.000 claims description 3
- 230000001427 coherent effect Effects 0.000 claims description 3
- 238000013461 design Methods 0.000 claims description 3
- 238000010845 search algorithm Methods 0.000 claims description 3
- 239000000463 material Substances 0.000 abstract description 3
- 238000004364 calculation method Methods 0.000 abstract description 2
- 238000010586 diagram Methods 0.000 description 11
- 238000001514 detection method Methods 0.000 description 4
- 239000011435 rock Substances 0.000 description 4
- 230000010261 cell growth Effects 0.000 description 2
- 230000007797 corrosion Effects 0.000 description 2
- 238000005260 corrosion Methods 0.000 description 2
- 238000011161 development Methods 0.000 description 2
- 238000006073 displacement reaction Methods 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 238000003325 tomography Methods 0.000 description 2
- 239000002699 waste material Substances 0.000 description 2
- 238000004458 analytical method Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000015572 biosynthetic process Effects 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/288—Event detection in seismic signals, e.g. microseismics
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/30—Analysis
- G01V1/307—Analysis for determining seismic attributes, e.g. amplitude, instantaneous phase or frequency, reflection strength or polarity
Landscapes
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Environmental & Geological Engineering (AREA)
- Geology (AREA)
- Acoustics & Sound (AREA)
- General Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Geophysics (AREA)
- Business, Economics & Management (AREA)
- Emergency Management (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
本发明公开了一种基于地震相干体切片的裂缝自动检测方法，涉及石油勘探图像处理技术领域，本发明包括如下步骤：S1、中值滤波去噪；S2、归一化处理；S3、双边滤波与拉普拉斯滤波处理；S4、连通分量标记去除孤立噪声；S5、形态学操作处理；S6、漫水填充算法修复2D地震相干体切片；S7、FPA算法提取骨架；S8、MyPadding函数增强裂缝连续性；S9、连通分量标记；S10、最小二乘直线拟合方法处理；S11、分类显示，按照不同的方位角的范围对裂缝进行分类，得到最终的2D地震相干体切片数据的裂缝自动检测结果，本发明能够直接利用地震相干体切片数据提取出裂缝信息，实现裂缝的方位角计算，极大提高地震切片裂缝解释效率，节约大量人力物力成本。
Description
技术领域
本发明涉及石油勘探图像处理技术领域，更具体的是涉及一种基于地震相干体切片的裂缝自动检测方法。
背景技术
目前，随着国内外油气勘探开发技术的提高，在世界范围内易开采的大型油气田基本已被开发的背景下，人们对复杂岩性、裂缝型油气田的开发越来越重视。地震相干体切片技术就是为适应这一趋势而发展起来的地震资料精细解释技术，它是指沿着三维地震数据体的某一个方向，以平面或曲面方式抽取的，具有地球物理意义的二维空间属性载体。
地震相干体切片数据的裂缝最基本的特征是由其结构形态决定的。由于地质运动，地下岩石受力产生断裂，如果构造应力较大，岩石沿断裂面发生明显位移，则称这种断裂为断层；如果构造应力较小，岩石受力后断开并沿断裂面无显著位移，则称这种断裂构造为裂缝，在地震相干体切片数据中往往需要同时提取断层与裂缝，所以本发明中的所有裂缝自动检测既指受构造应力岩石有明显位移的断层构造，又指无明显位移的裂缝构造。
现有的地震相干体切片解释的一大难点在于缺少数学化的描述，对地震相干体切片数据解释的好与坏往往依赖于地震解释人员的经验与直觉，而优秀的经验丰富的地震解释人员是非常稀少的；同时，地震相干体切片解释的另一大难点在于其是一个非常耗时、工作量巨大的任务，目前主流的做法是地震解释人员根据相干体切片数据人工一条条进行划线，并且在实际工程中由于人力的不足，地震解释人员往往只挑选典型的、重要的地震相干体切片进行解释，这就造成了大量的信息浪费，解释结果不够精确的问题。
发明内容
本发明的目的在于：为了解决现有的地震相干体切片解释缺少数学化的描述，依赖于地震解释人员的经验和直觉，并且仅对典型的的地震相干体切片进行解释，造成信息浪费、解释结果不精确的问题，本发明提供一种基于地震相干体切片的裂缝自动检测方法。
本发明为了实现上述目的具体采用以下技术方案：
一种基于地震相干体切片的裂缝自动检测方法，其特征在于，包括如下步骤：
S1、中值滤波去噪
输入2D地震相干体切片数据，对输入的2D地震相干体切片数据进行中值滤波，去除由于相干算法不足、数据采集误差等造成的椒盐噪声；
S2、归一化处理
利用离差标准化方法对中值滤波后的数据进行归一化处理，使中值滤波后的2D地震相干体切片数据变化范围在[0,-1]之间；
S3、双边滤波与拉普拉斯滤波处理
利用双边滤波方法处理归一化后的2D地震相干体切片数据，增强裂缝形态特征，保持裂缝连续性；然后再利用拉普拉斯滤波锐化双边滤波后的2D地震相干体切片，增强裂缝信息；并对拉普拉斯滤波后的2D地震相干体切片进行二值化处理；
二值化的阈值由两部分组成，其一为利用Otsu大津算法自动计算的阈值Threshold_Otsu，其二是用户输入的调整阈值Threshold_Adjust，所述二值化的阈值为上述两阈值的叠加，即
Threshold＝Threshold_Otsu+Threshold_Adjust；
S4、连通分量标记去除孤立噪声
采用四连通分量标记算法对二值化后的2D地震相干体切片的所有前景像素点进行连通分量标记，统计每个标记值的前景像素点个数，去除前景像素点个数小的连通分量，以抑制孤立噪声干扰；
S5、形态学操作处理
运用形态学闭操作处理去除了孤立噪声的二值化后的2D地震相干体切片数据，以增强裂缝的连续性和完整性；
S6、漫水填充算法修复2D地震相干体切片
由于利用二值化只能提取到裂缝的边缘轮廓信息，裂缝内部会形成大量孔洞，所以利用漫水填充算法填充形态学闭操作处理后的2D地震相干体切片中的孔洞，对2D地震相干体切片进行修复；
S7、FPA算法提取骨架
利用FPA算法对修复后的2D地震相干体切片的连通区域进行骨架提取，得到单像素宽度的裂缝；
S8、MyPadding函数增强裂缝连续性
利用MyPadding函数对单像素宽度的裂缝缺失的数据值和形态较大的波动进行修复，以增强裂缝的连续性；
S9、连通分量标记
利用八连通分量标记算法对增强了连续性的裂缝进行连通分量标记；
S10、最小二乘直线拟合方法处理
利用最小二乘直线拟合方法拟合每个八连通分量，得到拟合直线的倾角，结合实际地震资料方位角得到每条裂缝的方位角，所述实际地震资料方位角的坐标系不同于常用的直角坐标系，而是以大地正北为0°，顺时针旋转依次增加的坐标系。
S11、分类显示
按照不同的方位角的范围对裂缝进行分类，得到最终的2D地震相干体切片数据的裂缝自动检测结果。
进一步的，所述S1中用5×5的二维模板对输入的2D地震相干体切片数据进行中值滤波。
进一步的，所述S2中利用离差标准化方法对中值滤波后的数据进行归一化处理，具体包括如下步骤：
S2.1、输入n个数据x1,x2,...,xn，找出这n个数据中的最小值xmin和最大值xmax；
S2.2、根据公式
其中1≤i≤n，得到归一化后的数据y1,y2,...,yn，此时y1,y2,...,yn∈[0,1]。
进一步的，所述S3中用3×3的双边滤波模板对归一化后的2D地震相干体切片数据进行处理，且双边滤波所用空域高斯函数的标准差的默认值为0.5，值域高斯函数的标准差的默认值为0.05。
进一步的，所述S3中用3×3的拉普拉斯滤波模板对双边滤波后的2D地震相干体切片进行处理，且拉普拉斯滤波所用中心系数默认值为-8。
进一步的，所述S4具体包括如下步骤：
S4.1、默认2D地震相干体切片数据的所有前景像素点为未标记状态，从左往右、从上到下遍历每个前景像素点，检查每个前景像素点左邻域和上邻域是否被标记，如果没有，则赋予该前景像素点一个新的标记值；如果有，则将左邻域和上邻域中最小的标记值赋予该前景像素点，其中标记值是从1开始，每次以1为增幅，顺序递增生成的；
S4.2、从左往右、从上到下遍历S4.1中每个带标记值的前景像素点，用与每个前景像素点分别四连通的前景像素点中的最小的标记值替换该前景像素点的标记值；
S4.3、统计具有相同标记值的前景像素点的个数，根据预先设定的阈值Threshold_Count，删除前景像素点个数小于该阈值的标记值。
进一步的，所述阈值Threshold_Count的默认值为30。
进一步的，所述S5具体包括如下步骤：
S5.1、采用3×3的矩形结构元膨胀抑制了孤立噪声的二值图；
S5.2、采用3×3的矩形结构元腐蚀S7.1膨胀的结果。
进一步的，所述S6具体包括如下步骤：
S6.1、在2D地震相干体切片的二值图的下方和右方分别扩展一行和一列，得到扩展图，并将扩展的行和列的像素值都赋值为0；
S6.2、以扩展图左上角(1,1)像素点为起始点，运用宽度优先搜索算法搜索扩展图，将所有搜索到的像素点的像素值均标记为-1；
S6.3、遍历扩展图，如果像素点的像素值为-1，则将其设置为1，否则将其设置为0；
S6.4、截取扩展之前的二值图所在区域，得到利用漫水填充算法填充了孔洞的二值图。
进一步的，所述S7具体包括如下步骤：
S7.1、考察填充了孔洞的二值图的每一个前景像素点，判断是否存在同时满足条件(a)、(b)、(c)和(d)的前景像素点，若存在，则将其标记下来，执行S7.2，否则，执行S7.5；
S7.2、删除S7.1标记的前景像素点的像素值，执行S7.3；
S7.3、考察填充了孔洞的二值图的每一个前景像素点，判断是否存在同时满足条件(a)、(b)、(e)和(f)的前景像素点，若存在，则将其标记下来，执行S7.4，否则，执行S7.5；
S7.4、删除S7.3标记的前景像素点的像素值，完成一轮迭代，执行S7.1；
S7.5、迭代结束，得到骨架提取的结果，即单像素宽度的裂缝；
所述条件(a)、(b)、(c)、(d)、(e)和(f)具体为：
在3×3的邻域中，设中心像素点为P，则其正上方像素点为P1，并且按顺时针方向旋转，依次排布P2,P3,...,P8，则有：
(a)2≤B(P)≤6；
(b)A(P)＝1；
(c)P1×P3×P5＝0；
(d)P3×P5×P7＝0；
(e)P1×P3×P7＝0；
(f)P1×P5×P5＝0；
其中，A(P)表示P1,P2,...,P8,P1序列中“01”对出现的个数，B(P)表示P的八邻域窗口中1的个数。
进一步的，所述S8利用MyPadding函数对单像素宽度的裂缝缺失的数据值和形态较大的波动进行修复，具体为：
S8.1、设计8个3×3模板，其中：
第一个模板除了中心像素点的上邻域和下邻域像素值为1，其余像素点的像素值均为0；
第二个模板除了中心像素点的左邻域和右邻域像素值为1，其余像素点的像素值均为0；
第三个模板除了中心像素点的左上邻域和右下邻域像素值为1，其余像素点的像素值均为0；
第四个模板除了中心像素点的左下邻域和右上邻域像素值为1，其余像素点的像素值均为0；
第五个模板除了中心像素点的左邻域和上邻域像素值为1，其余像素点的像素值均为0；
第六个模板除了中心像素点的上邻域和右邻域像素值为1，其余像素点的像素值均为0；
第七个模板除了中心像素点的下邻域和右邻域像素值为1，其余像素点的像素值均为0；
第八个模板除了中心像素点的左邻域和下邻域像素值为1，其余像素点的像素值均为0；
S8.2、遍历单像素裂缝的每个像素点，以每个像素点为中心像素点，取每个像素点所在3×3的邻域分别与所述8个模板做逻辑与操作，若逻辑与值为真，则将该中心像素点的像素值设为1。
进一步的，所述S9具体为：
S9.1、默认经过S8修复的裂缝的前景像素点状态为未标记，从左到右、从上到下遍历增强了连续性的裂缝的所有前景像素点，检查该前景像素点的左邻域、上邻域、左上邻域和右上邻域是否被标记，如果没有，则赋予该前景像素点一个新的标记值；如果有，则将左邻域、上邻域、左上邻域和右上邻域中最小的标记值赋予该前景像素点，其中标记值是从1开始，每次以1为增幅，顺序递增生成的；
S9.2、从左到右、从上到下遍历S9.1中每个带标记值的前景像素点，用与每个前景像素点八连通的像素点中的最小标记值替换该前景像素点的标记值。
本发明的有益效果如下：
1、本发明能够直接利用2D地震相干体切片数据提取出切片中存在的裂缝，精确直观显示工区裂缝发育情况，克服了传统地震解释方法不能很好解释裂缝信息的缺点。
2、本发明利用直线拟合方法能够计算裂缝的方位角，并按照不同的方位角的范围对裂缝进行分类显示，能够直观显示工区裂缝分布的情况。
3、本发明通过双边滤波、孔洞填充、MyPadding函数模板匹配与骨架提取等数字图像处理算法，能够较好抑制2D地震相干体切片的噪声，同时保持裂缝的连续性与完整性，提取出完整、合理、连续性好的裂缝。
4、本发明利用计算机自动检测2D地震相干体切片中的裂缝，实现了裂缝检测的自动化，可以节约大量人力物力成本，实现较高的地震切片解释精度。
5、本发明能够移植到C++平台下运行，并且也能够移植到其他平台，适用范围广。
附图说明
图1是本发明的方法流程图。
图2是2D地震相干体切片示意图。
图3是中值滤波后的切片示意图。
图4是双边滤波后的切片示意图。
图5是拉普拉斯滤波后的切片示意图。
图6是二值化后的切片示意图。
图7是连通分量标记去除孤立噪声后的切片示意图。
图8是形态学闭操作后的切片示意图。
图9是漫水填充算法填充孔洞后的切片示意图。
图10是FPA算法骨架提取后的切片示意图。
图11是MyPadding函数增强裂缝连续性后的切片示意图。
图12是2D地震相干体切片数据的裂缝自动检测结果示意图。
具体实施方式
为了本技术领域的人员更好的理解本发明，下面结合附图和以下实施例对本发明作进一步详细描述。
实施例1
如图1到图12所示，本实施例提供一种基于地震相干体切片的裂缝自动检测方法，包括如下步骤：
S1、中值滤波去噪
输入2D地震相干体切片数据，对输入的2D地震相干体切片数据进行中值滤波，去除由于相干算法不足、数据采集误差等造成的椒盐噪声；
所述S1中用5×5的二维模板对输入的2D地震相干体切片数据进行中值滤波；
S2、归一化处理
利用离差标准化方法对中值滤波后的数据进行归一化处理，使中值滤波后的2D地震相干体切片数据变化范围在[0,-1]之间；
所述S2中利用离差标准化方法对中值滤波后的数据进行归一化处理，具体包括如下步骤：
S2.1、输入n个数据x1,x2,...,xn，找出这n个数据中的最小值xmin和最大值xmax；
S2.2、根据公式
其中1≤i≤n，得到归一化后的数据y1,y2,...,yn，此时y1,y2,...,yn∈[0,1]；
S3、双边滤波与拉普拉斯滤波处理
中值滤波虽然去除了切片中的椒盐噪声，但此时的裂缝连续性较差，裂缝形态不完整，利用双边滤波方法处理归一化后的2D地震相干体切片数据，既可以平滑图像，增强裂缝形态特征，又能不破坏图像边缘信息，保持裂缝连续性；
用3×3的双边滤波模板对归一化后的2D地震相干体切片数据进行处理，且双边滤波所用空域高斯函数的标准差的默认值为0.5，值域高斯函数的标准差的默认值为0.05；
然后再利用拉普拉斯滤波锐化双边滤波后的2D地震相干体切片，增强裂缝信息；
用3×3的拉普拉斯滤波模板对双边滤波后的2D地震相干体切片进行处理，本实施例中拉普拉斯滤波所用中心系数默认值为-8，并且另外提供-4和-32的中心系数供用户选择；
并对拉普拉斯滤波后的2D地震相干体切片进行二值化处理；
二值化的阈值由两部分组成，其一为利用Otsu大津算法自动计算的阈值Threshold_Otsu，其二是用户输入的调整阈值Threshold_Adjust，所述二值化的阈值为上述两阈值的叠加，即
Threshold＝Threshold_Otsu+Threshold_Adjust；
S4、连通分量标记去除孤立噪声
采用四连通分量标记算法对二值化后的2D地震相干体切片的所有前景像素点进行连通分量标记，统计每个标记值的前景像素点个数，去除前景像素点个数小的连通分量，以抑制孤立噪声干扰；
具体包括如下步骤：
S4.1、默认2D地震相干体切片数据的所有前景像素点为未标记状态，从左往右、从上到下遍历每个前景像素点，检查每个前景像素点左邻域和上邻域是否被标记，如果没有，则赋予该前景像素点一个新的标记值；如果有，则将左邻域和上邻域中最小的标记值赋予该前景像素点，其中标记值是从1开始，每次以1为增幅，顺序递增生成的；
S4.2、从左往右、从上到下遍历S4.1中每个带标记值的前景像素点，用与每个前景像素点分别四连通的前景像素点中的最小的标记值替换该前景像素点的标记值；
S4.3、统计具有相同标记值的前景像素点的个数，根据预先设定的阈值Threshold_Count，删除前景像素点个数小于该阈值的标记值，本实施例中所述阈值Threshold_Count的默认值为30；
S5、形态学操作处理
运用形态学闭操作处理去除了孤立噪声的二值化后的2D地震相干体切片数据，以增强裂缝的连续性和完整性；
具体包括如下步骤：
S5.1、采用3×3的矩形结构元膨胀抑制了孤立噪声的二值图；
S5.2、采用3×3的矩形结构元腐蚀S5.1膨胀的结果；
S6、漫水填充算法修复2D地震相干体切片
由于利用二值化只能提取到裂缝的边缘轮廓信息，裂缝内部会形成大量孔洞，所以利用漫水填充算法填充形态学闭操作处理后的2D地震相干体切片中的孔洞，对2D地震相干体切片进行修复，具体包括如下步骤：
S6.1、在2D地震相干体切片的二值图的下方和右方分别扩展一行和一列，得到扩展图，并将扩展的行和列的像素值都赋值为0；
S6.2、以扩展图左上角(1,1)像素点为起始点，运用宽度优先搜索算法搜索扩展图，将所有搜索到的像素点的像素值均标记为-1；
S6.3、遍历扩展图，如果像素点的像素值为-1，则将其设置为1，否则将其设置为0；
S6.4、截取扩展之前的二值图所在区域，得到利用漫水填充算法填充了孔洞的二值图；
S7、FPA算法提取骨架
利用FPA算法对修复后的2D地震相干体切片的连通区域进行骨架提取，得到单像素宽度的裂缝，具体包括如下步骤：
S7.1、考察填充了孔洞的二值图的每一个前景像素点，判断是否存在同时满足条件(a)、(b)、(c)和(d)的前景像素点，若存在，则将其标记下来，执行S7.2，否则，执行S7.5；
S7.2、删除S7.1标记的前景像素点的像素值，执行S7.3；
S7.3、考察填充了孔洞的二值图的每一个前景像素点，判断是否存在同时满足条件(a)、(b)、(e)和(f)的前景像素点，若存在，则将其标记下来，执行S7.4，否则，执行S7.5；
S7.4、删除S7.3标记的前景像素点的像素值，完成一轮迭代，执行S7.1；
S7.5、迭代结束，得到骨架提取的结果，即单像素宽度的裂缝；
所述条件(a)、(b)、(c)、(d)、(e)和(f)具体为：
在3×3的邻域中，设中心像素点为P，则其正上方像素点为P1，并且按顺时针方向旋转，依次排布P2,P3,...,P8，则有：
(a)2≤B(P)≤6；
(b)A(P)＝1；
(c)P1×P3×P5＝0；
(d)P3×P5×P7＝0；
(e)P1×P3×P7＝0；
(f)P1×P5×P5＝0；
其中，A(P)表示P1,P2,...,P8,P1序列中“01”对出现的个数，B(P)表示P的八邻域窗口中1的个数；
S8、MyPadding函数增强裂缝连续性
利用MyPadding函数对单像素宽度的裂缝缺失的数据值和形态较大的波动进行修复，以增强裂缝的连续性，具体包括如下步骤：
S8.1、设计8个3×3模板，其中：
第一个模板除了中心像素点的上邻域和下邻域像素值为1，其余像素点的像素值均为0；
第二个模板除了中心像素点的左邻域和右邻域像素值为1，其余像素点的像素值均为0；
第三个模板除了中心像素点的左上邻域和右下邻域像素值为1，其余像素点的像素值均为0；
第四个模板除了中心像素点的左下邻域和右上邻域像素值为1，其余像素点的像素值均为0；
第五个模板除了中心像素点的左邻域和上邻域像素值为1，其余像素点的像素值均为0；
第六个模板除了中心像素点的上邻域和右邻域像素值为1，其余像素点的像素值均为0；
第七个模板除了中心像素点的下邻域和右邻域像素值为1，其余像素点的像素值均为0；
第八个模板除了中心像素点的左邻域和下邻域像素值为1，其余像素点的像素值均为0；
S8.2、遍历单像素裂缝的每个像素点，以每个像素点为中心像素点，取每个像素点所在3×3的邻域分别与所述8个模板做逻辑与操作，若逻辑与值为真，则将该中心像素点的像素值设为1；
S9、连通分量标记
利用八连通分量标记算法对增强了连续性的裂缝进行连通分量标记，具体包括如下步骤：
S9.1、默认经过S8修复的裂缝的前景像素点状态为未标记，从左到右、从上到下遍历增强了连续性的裂缝的所有前景像素点，检查该前景像素点的左邻域、上邻域、左上邻域和右上邻域是否被标记，如果没有，则赋予该前景像素点一个新的标记值；如果有，则将左邻域、上邻域、左上邻域和右上邻域中最小的标记值赋予该前景像素点，其中标记值是从1开始，每次以1为增幅，顺序递增生成的；
S9.2、从左到右、从上到下遍历S9.1中每个带标记值的前景像素点，用与每个前景像素点八连通的像素点中的最小标记值替换该前景像素点的标记值；
S10、最小二乘直线拟合方法处理
利用最小二乘直线拟合方法拟合每个八连通分量，得到拟合直线的倾角，结合实际地震资料的方位角得到每条裂缝的方位角；
S11、分类显示
按照不同的方位角的范围对裂缝进行分类，得到最终的2D地震相干体切片数据的裂缝自动检测结果，如图12所示，本实施例中黑色裂缝代表方位角范围为(0°，-90°]，白色裂缝代表方位角范围为(90°，-180°]。
本实施例主要用在油气勘探开发领域，为三维地质建模、储层建模以及地震相带分析等提供基础数据，不仅能够直接利用2D地震相干体切片数据提取出裂缝信息，而且能够实现裂缝的方位角计算，能够很好地解决噪声干扰严重、形态特征不完整的切片数据提取裂缝中的裂缝检测不准确、裂缝连续性差的问题；相比于传统由地震解释人员手工解释，本实施例能够极大提高地震切片裂缝解释效率，节约大量人力物力成本。
以上所述，仅为本发明的较佳实施例，并不用以限制本发明，本发明的专利保护范围以权利要求书为准，凡是运用本发明的说明书及附图内容所作的等同结构变化，同理均应包含在本发明的保护范围内。
Claims (10)
1.一种基于地震相干体切片的裂缝自动检测方法，其特征在于，包括如下步骤：
S1、中值滤波去噪
输入2D地震相干体切片数据，对输入的2D地震相干体切片数据进行中值滤波，去除由于相干算法不足和数据采集误差造成的椒盐噪声；
S2、归一化处理
利用离差标准化方法对中值滤波后的数据进行归一化处理，使中值滤波后的2D地震相干体切片数据变化范围在[0,-1]之间；
S3、双边滤波与拉普拉斯滤波处理
利用双边滤波方法处理归一化后的2D地震相干体切片数据，增强裂缝形态特征，保持裂缝连续性；然后再利用拉普拉斯滤波锐化双边滤波后的2D地震相干体切片数据，增强裂缝信息；并对拉普拉斯滤波后的2D地震相干体切片数据进行二值化处理；
S4、连通分量标记去除孤立噪声
采用四连通分量标记算法对二值化后的2D地震相干体切片数据的所有前景像素点进行连通分量标记，统计每个标记值的前景像素点个数，去除前景像素点个数小的连通分量，以抑制孤立噪声干扰；
S5、形态学操作处理
运用形态学闭操作处理去除了孤立噪声的二值化后的2D地震相干体切片数据，以增强裂缝的连续性和完整性；
S6、漫水填充算法修复2D地震相干体切片
利用漫水填充算法填充形态学闭操作处理后的2D地震相干体切片数据，对2D地震相干体切片数据进行修复；
S7、FPA算法提取骨架
利用FPA算法对修复后的2D地震相干体切片数据的连通区域进行骨架提取，得到单像素宽度的裂缝；
S8、MyPadding函数增强裂缝连续性
利用MyPadding函数对单像素宽度的裂缝缺失的数据值和形态较大的波动进行修复，以增强裂缝的连续性；
S9、连通分量标记
利用八连通分量标记算法对增强了连续性的裂缝进行连通分量标记；
S10、最小二乘直线拟合方法处理
利用最小二乘直线拟合方法拟合每个八连通分量，得到拟合直线的倾角，结合实际地震资料方位角得到每条裂缝的方位角；
S11、分类显示
按照不同的方位角的范围对裂缝进行分类，得到最终的2D地震相干体切片数据的裂缝自动检测结果。
2.根据权利要求1所述的一种基于地震相干体切片的裂缝自动检测方法，其特征在于，所述S2中利用离差标准化方法对中值滤波后的数据进行归一化处理，具体包括如下步骤：
S2.1、输入n个数据x1,x2,...,xn，找出这n个数据中的最小值xmin和最大值xmax；
S2.2、根据公式
其中1≤i≤n，得到归一化后的数据y1,y2,...,yn，此时y1,y2,...,yn∈[0,1]。
3.根据权利要求1所述的一种基于地震相干体切片的裂缝自动检测方法，其特征在于：所述S3中用3×3的双边滤波模板对归一化后的2D地震相干体切片数据进行处理，且双边滤波所用空域高斯函数的标准差的默认值为0.5，值域高斯函数的标准差的默认值为0.05。
4.根据权利要求1所述的一种基于地震相干体切片的裂缝自动检测方法，其特征在于：所述S3中用3×3的拉普拉斯滤波模板对双边滤波后的2D地震相干体切片进行处理，且拉普拉斯滤波所用中心系数默认值为-8。
5.根据权利要求1所述的一种基于地震相干体切片的裂缝自动检测方法，其特征在于，所述S4具体包括如下步骤：
S4.1、默认2D地震相干体切片数据的所有前景像素点为未标记状态，从左往右、从上到下遍历每个前景像素点，检查每个前景像素点左邻域和上邻域是否被标记，如果没有，则赋予该前景像素点一个新的标记值；如果有，则将左邻域和上邻域中最小的标记值赋予该前景像素点，其中标记值是从1开始，每次以1为增幅，顺序递增生成的；
S4.2、从左往右、从上到下遍历S4.1中每个带标记值的前景像素点，用与每个前景像素点分别四连通的前景像素点中的最小的标记值替换该前景像素点的标记值；
S4.3、统计具有相同标记值的前景像素点的个数，根据预先设定的阈值Threshold_Count，删除前景像素点个数小于该阈值的标记值。
6.根据权利要求5所述的一种基于地震相干体切片的裂缝自动检测方法，其特征在于：所述阈值Threshold_Count的默认值为30。
7.根据权利要求1所述的一种基于地震相干体切片的裂缝自动检测方法，其特征在于，所述S6具体包括如下步骤：
S6.1、在2D地震相干体切片的二值图的下方和右方分别扩展一行和一列，得到扩展图，并将扩展的行和列的像素值都赋值为0；
S6.2、以扩展图左上角(1,1)像素点为起始点，运用宽度优先搜索算法搜索扩展图，将所有搜索到的像素点的像素值均标记为-1；
S6.3、遍历扩展图，如果像素点的像素值为-1，则将其设置为1，否则将其设置为0；
S6.4、截取扩展之前的二值图所在区域，得到利用漫水填充算法填充了孔洞的二值图。
8.根据权利要求1所述的一种基于地震相干体切片的裂缝自动检测方法，其特征在于，所述S7具体包括如下步骤：
S7.1、考察填充了孔洞的二值图的每一个前景像素点，判断是否存在同时满足条件(a)、(b)、(c)和(d)的前景像素点，若存在，则将其标记下来，执行S7.2，否则，执行S7.5；
S7.2、删除S7.1标记的前景像素点的像素值，执行S7.3；
S7.3、考察填充了孔洞的二值图的每一个前景像素点，判断是否存在同时满足条件(a)、(b)、(e)和(f)的前景像素点，若存在，则将其标记下来，执行S7.4，否则，执行S7.5；
S7.4、删除S7.3标记的前景像素点的像素值，完成一轮迭代，执行S7.1；
S7.5、迭代结束，得到骨架提取的结果，即单像素宽度的裂缝；
所述条件(a)、(b)、(c)、(d)、(e)和(f)具体为：
在3×3的邻域中，设中心像素点为P，则其正上方像素点为P1，并且按顺时针方向旋转，依次排布P2,P3,...,P8，则有：
(a)2≤B(P)≤6；
(b)A(P)＝1；
(c)P1×P3×P5＝0；
(d)P3×P5×P7＝0；
(e)P1×P3×P7＝0；
(f)P1×P5×P5＝0；
其中，A(P)表示P1,P2,...,P8,P1序列中“01”对出现的个数，B(P)表示P的八邻域窗口中1的个数。
9.根据权利要求1所述的一种基于地震相干体切片的裂缝自动检测方法，其特征在于，所述S8利用MyPadding函数对单像素宽度的裂缝缺失的数据值和形态较大的波动进行修复，具体为：
S8.1、设计8个3×3模板，其中：
第一个模板除了中心像素点的上邻域和下邻域像素值为1，其余像素点的像素值均为0；
第二个模板除了中心像素点的左邻域和右邻域像素值为1，其余像素点的像素值均为0；
第三个模板除了中心像素点的左上邻域和右下邻域像素值为1，其余像素点的像素值均为0；
第四个模板除了中心像素点的左下邻域和右上邻域像素值为1，其余像素点的像素值均为0；
第五个模板除了中心像素点的左邻域和上邻域像素值为1，其余像素点的像素值均为0；
第六个模板除了中心像素点的上邻域和右邻域像素值为1，其余像素点的像素值均为0；
第七个模板除了中心像素点的下邻域和右邻域像素值为1，其余像素点的像素值均为0；
第八个模板除了中心像素点的左邻域和下邻域像素值为1，其余像素点的像素值均为0；
S8.2、遍历单像素裂缝的每个像素点，以每个像素点为中心像素点，取每个像素点所在3×3的邻域分别与所述8个模板做逻辑与操作，若逻辑与值为真，则将该中心像素点的像素值设为1；否则，不操作。
10.根据权利要求1所述的一种基于地震相干体切片的裂缝自动检测方法，其特征在于，所述S9具体为：
S9.1、默认经过S8增强了连续性的裂缝的前景像素点状态为未标记，从左到右、从上到下遍历增强了连续性的裂缝的所有前景像素点，检查每个前景像素点的左邻域、上邻域、左上邻域和右上邻域是否被标记，如果没有，则赋予该前景像素点一个新的标记值；如果有，则将左邻域、上邻域、左上邻域和右上邻域中最小的标记值赋予该前景像素点，其中标记值是从1开始，每次以1为增幅，顺序递增生成的；
S9.2、从左到右、从上到下遍历S9.1中每个带标记值的前景像素点，用与每个前景像素点八连通的像素点中的最小标记值替换该前景像素点的标记值。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810499447.6A CN108957530B (zh) | 2018-05-23 | 2018-05-23 | 一种基于地震相干体切片的裂缝自动检测方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810499447.6A CN108957530B (zh) | 2018-05-23 | 2018-05-23 | 一种基于地震相干体切片的裂缝自动检测方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN108957530A CN108957530A (zh) | 2018-12-07 |
| CN108957530B true CN108957530B (zh) | 2019-08-23 |
Family
ID=64499363
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201810499447.6A Active CN108957530B (zh) | 2018-05-23 | 2018-05-23 | 一种基于地震相干体切片的裂缝自动检测方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN108957530B (zh) |
Families Citing this family (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109584240B (zh) * | 2018-12-20 | 2022-05-03 | 成都理工大学 | 滑坡后缘裂缝位移图像识别方法 |
| CN111580156B (zh) * | 2019-02-18 | 2022-12-02 | 中国石油天然气股份有限公司 | 一种地震零值切片自动拾取方法及系统 |
| CN111751871B (zh) * | 2019-03-29 | 2023-02-07 | 中国石油天然气集团有限公司 | 三维地震沿层相干切片处理方法、装置及电子设备 |
| CN114252913A (zh) * | 2020-09-25 | 2022-03-29 | 中国石油天然气股份有限公司 | 平面断层信息的识别方法及装置 |
Citations (13)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN1875378A (zh) * | 2003-11-12 | 2006-12-06 | 英国电讯有限公司 | 图像中的对象检测 |
| CN101738639A (zh) * | 2008-11-24 | 2010-06-16 | 中国石油天然气集团公司 | 提高岩石裂缝参数计算精度的方法 |
| CN102053277A (zh) * | 2009-10-30 | 2011-05-11 | 中国石油化工股份有限公司 | 一种利用地震资料进行检测储层裂缝发育方向的方法 |
| CN102305945A (zh) * | 2011-06-20 | 2012-01-04 | 电子科技大学 | 一种线性噪声消除方法 |
| CN103064114A (zh) * | 2011-10-18 | 2013-04-24 | 中国石油化工股份有限公司 | 一种裂缝储层的表征方法及装置 |
| CN103076623A (zh) * | 2011-10-25 | 2013-05-01 | 中国石油化工股份有限公司 | 一种基于叠前相干的裂缝检测方法 |
| CN103278864A (zh) * | 2013-05-10 | 2013-09-04 | 中国石油天然气股份有限公司 | 洞缝型储层的地质特征参数及分布的测定方法及装置 |
| CN103592690A (zh) * | 2013-10-24 | 2014-02-19 | 长江大学 | 基于电成像测井孔隙度谱信息自动识别储层裂缝的方法 |
| CN103927526A (zh) * | 2014-04-30 | 2014-07-16 | 长安大学 | 一种基于高斯差分多尺度边缘融合的车辆检测方法 |
| CN104570083A (zh) * | 2013-10-29 | 2015-04-29 | 中国石油天然气集团公司 | 基于多维地震属性的地质体自动识别方法 |
| CN103135131B (zh) * | 2011-11-28 | 2015-05-20 | 中国石油化工股份有限公司 | 一种针对裂缝性储层预测的解释装置 |
| CN105425299A (zh) * | 2015-11-17 | 2016-03-23 | 中国石油天然气股份有限公司 | 确定地层裂缝分布的方法和装置 |
| CN107918153A (zh) * | 2018-01-10 | 2018-04-17 | 成都理工大学 | 一种地震信号相干性高精度检测方法 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US9389326B2 (en) * | 2011-03-23 | 2016-07-12 | Global Ambient Seismic, Inc. | Methods, systems and devices for near-well fracture monitoring using tomographic fracture imaging techniques |
-
2018
- 2018-05-23 CN CN201810499447.6A patent/CN108957530B/zh active Active
Patent Citations (13)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN1875378A (zh) * | 2003-11-12 | 2006-12-06 | 英国电讯有限公司 | 图像中的对象检测 |
| CN101738639A (zh) * | 2008-11-24 | 2010-06-16 | 中国石油天然气集团公司 | 提高岩石裂缝参数计算精度的方法 |
| CN102053277A (zh) * | 2009-10-30 | 2011-05-11 | 中国石油化工股份有限公司 | 一种利用地震资料进行检测储层裂缝发育方向的方法 |
| CN102305945A (zh) * | 2011-06-20 | 2012-01-04 | 电子科技大学 | 一种线性噪声消除方法 |
| CN103064114A (zh) * | 2011-10-18 | 2013-04-24 | 中国石油化工股份有限公司 | 一种裂缝储层的表征方法及装置 |
| CN103076623A (zh) * | 2011-10-25 | 2013-05-01 | 中国石油化工股份有限公司 | 一种基于叠前相干的裂缝检测方法 |
| CN103135131B (zh) * | 2011-11-28 | 2015-05-20 | 中国石油化工股份有限公司 | 一种针对裂缝性储层预测的解释装置 |
| CN103278864A (zh) * | 2013-05-10 | 2013-09-04 | 中国石油天然气股份有限公司 | 洞缝型储层的地质特征参数及分布的测定方法及装置 |
| CN103592690A (zh) * | 2013-10-24 | 2014-02-19 | 长江大学 | 基于电成像测井孔隙度谱信息自动识别储层裂缝的方法 |
| CN104570083A (zh) * | 2013-10-29 | 2015-04-29 | 中国石油天然气集团公司 | 基于多维地震属性的地质体自动识别方法 |
| CN103927526A (zh) * | 2014-04-30 | 2014-07-16 | 长安大学 | 一种基于高斯差分多尺度边缘融合的车辆检测方法 |
| CN105425299A (zh) * | 2015-11-17 | 2016-03-23 | 中国石油天然气股份有限公司 | 确定地层裂缝分布的方法和装置 |
| CN107918153A (zh) * | 2018-01-10 | 2018-04-17 | 成都理工大学 | 一种地震信号相干性高精度检测方法 |
Non-Patent Citations (2)
| Title |
|---|
| "利用相干体技术研究碳酸盐岩裂缝特征";万学鹏 等;《断块油气田》;20070531;第14卷(第3期);第43-45页 * |
| "综合利用测井—地震方法识别火成岩裂缝";黄悍东 等;《石油地球物理勘探》;20151031;第50卷(第05期);第942-950页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN108957530A (zh) | 2018-12-07 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN108957530B (zh) | 一种基于地震相干体切片的裂缝自动检测方法 | |
| Mallast et al. | Derivation of groundwater flow-paths based on semi-automatic extraction of lineaments from remote sensing data | |
| Thannoun | Automatic extraction and geospatial analysis of lineaments and their tectonic significance in some areas of Northern Iraq using remote sensing techniques and GIS | |
| CN104459768B (zh) | 一种基于可视化的三维空间目标地质体追踪方法 | |
| US20160188956A1 (en) | Identifying matching properties between a group of bodies representing a geological structure and a table of properties | |
| GB2502681A (en) | Color model transform and edge detection for seismic data | |
| CN106971156A (zh) | 一种基于面向对象分类的稀土开采区遥感信息提取方法 | |
| CN116152461B (zh) | 地质建模方法、装置、计算机设备及计算机可读存储介质 | |
| CN109738947A (zh) | 一种圈定砂岩型铀矿找矿远景区的物化探组合方法 | |
| CN103969683B (zh) | 一种三维地震解释中基于约束的批量拾取层位面的方法 | |
| CN109872393B (zh) | 一种基于地上、地下地质信息的三维地质数据处理方法 | |
| Prabhakaran et al. | An automated fracture trace detection technique using the complex shearlet transform | |
| Mohammadpour et al. | Automatic lineament extraction method in mineral exploration using CANNY algorithm and hough transform | |
| CN103035006A (zh) | 一种LiDAR辅助下基于LEGION的高分辨率航空影像分割方法 | |
| Gong et al. | Roof-cut guided localization for building change detection from imagery and footprint map | |
| Everett et al. | Remote sensing and GIS enable future exploration success | |
| Wylie Jr et al. | Well-log tomography and 3-D imaging of core and log-curve amplitudes in a Niagaran reef, Belle River Mills field, St. Clair County, Michigan, United States | |
| Mallast et al. | Semi-automatic extraction of lineaments from remote sensing data and the derivation of groundwater flow-paths. | |
| Li et al. | Fracture extraction from FMI based on multiscale mathematical morphology | |
| CN108805154A (zh) | 一种基于空间聚类的地质断层识别方法 | |
| WO2022245869A1 (en) | Utilization of geologic orientation data | |
| Kurtzman et al. | Improving fractured carbonate-reservoir characterization with remote sensing of beds, fractures, and vugs | |
| Refayee et al. | Fault and fracture detection in unconventional reservoirs: A Utica shale study | |
| CN107621655A (zh) | 基于DoS滤波的三维数据断层增强方法 | |
| Haryono et al. | A comparison of lineament and fracture trace extraction from LANDSAT ETM+ panchromatic band and panchromatic aerial photograph in Gunungsewu karst area, Java-Indonesia |
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 | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |