一种基于地震相干体切片的裂缝自动检测方法与流程

本发明涉及石油勘探图像处理技术领域，更具体的是涉及一种基于地震相干体切片的裂缝自动检测方法。

背景技术

目前，随着国内外油气勘探开发技术的提高，在世界范围内易开采的大型油气田基本已被开发的背景下，人们对复杂岩性、裂缝型油气田的开发越来越重视。地震相干体切片技术就是为适应这一趋势而发展起来的地震资料精细解释技术，它是指沿着三维地震数据体的某一个方向，以平面或曲面方式抽取的，具有地球物理意义的二维空间属性载体。

地震相干体切片数据的裂缝最基本的特征是由其结构形态决定的。由于地质运动，地下岩石受力产生断裂，如果构造应力较大，岩石沿断裂面发生明显位移，则称这种断裂为断层；如果构造应力较小，岩石受力后断开并沿断裂面无显著位移，则称这种断裂构造为裂缝，在地震相干体切片数据中往往需要同时提取断层与裂缝，所以本发明中的所有裂缝自动检测既指受构造应力岩石有明显位移的断层构造，又指无明显位移的裂缝构造。

现有的地震相干体切片解释的一大难点在于缺少数学化的描述，对地震相干体切片数据解释的好与坏往往依赖于地震解释人员的经验与直觉，而优秀的经验丰富的地震解释人员是非常稀少的；同时，地震相干体切片解释的另一大难点在于其是一个非常耗时、工作量巨大的任务，目前主流的做法是地震解释人员根据相干体切片数据人工一条条进行划线，并且在实际工程中由于人力的不足，地震解释人员往往只挑选典型的、重要的地震相干体切片进行解释，这就造成了大量的信息浪费，解释结果不够精确的问题。

技术实现要素：

本发明的目的在于：为了解决现有的地震相干体切片解释缺少数学化的描述，依赖于地震解释人员的经验和直觉，并且仅对典型的的地震相干体切片进行解释，造成信息浪费、解释结果不精确的问题，本发明提供一种基于地震相干体切片的裂缝自动检测方法。

本发明为了实现上述目的具体采用以下技术方案：

一种基于地震相干体切片的裂缝自动检测方法，其特征在于，包括如下步骤：

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地震相干体切片数据提取出裂缝信息，而且能够实现裂缝的方位角计算，能够很好地解决噪声干扰严重、形态特征不完整的切片数据提取裂缝中的裂缝检测不准确、裂缝连续性差的问题；相比于传统由地震解释人员手工解释，本实施例能够极大提高地震切片裂缝解释效率，节约大量人力物力成本。

以上所述，仅为本发明的较佳实施例，并不用以限制本发明，本发明的专利保护范围以权利要求书为准，凡是运用本发明的说明书及附图内容所作的等同结构变化，同理均应包含在本发明的保护范围内。