一种岩土体材料数字图像的分割阈值确定方法与流程

本发明涉及图像分割技术领域，尤其涉及一种岩土体材料数字图像的分割阈值确定方法。

背景技术：

随着微、细观成像技术的发展，研究人员可以利用该技术对岩土体的孔隙结构形态特征进行直接观测。对于sem扫描电镜技术，由于其较低的测试成本和良好的成像效果，被逐渐应用于各个领域。

使用sem扫描电镜技术对各种岩土体材料的表面微观孔隙结构进行直接观测，通过定量化的分析，对岩土体材料渗透率进行预测。在对sem图像定量化分析的过程中，最关键的一步就是对图像进行二值化。二值化的主要目的是将sem图像中的岩土体材料的孔隙结构与表面土骨架结构区分开来，提取出其中的孔隙结构，而孔隙结构是决定岩土体材料渗透率大小的主要因素，因此如何准确的提取孔隙结构就是数字图像法中二值化过程最为关键的技术。而二值化处理过程中，图像的分割阈值确定决定着孔隙结构提取二值化的准确性，因此，确定图像的分割阈值确定显得极为重要。

近几十年来，国内外众多学者已经提出了大量的图像二值化分割阈值确定计算方法。但由于数字图像的复杂性，以及图像分割问题依赖于具体应用领域的具体情况，因此，直到现在为止没有一种通用的二值化分割阈值确定算法。

中国专利公布号：cn102841220a；公布日：2012年12月26日，公开了一种基于孔隙率的黏土扫描电镜照片图像分割方法，包括步骤有：采用土力学常规实验方法测得黏土试样的干密度，算出干密度的三分之二次幂；实验室制备黏土试样的微结构样品；采用zd-a3冻干仪将制作的样品冻干并抽真空，真空度达到到10^-6pa时进行扫描电镜，观察样品的微结构并获得扫描照片；采用图像分析软件提取两个阈值对应的二值化数据；对比两个阈值条件下抽取到颗粒体和孔隙率的参数特征，算出孔隙体面积比；当孔隙体面积比等于干密度的三分之二次幂时，停止逼近，此时的阈值即为与所测干密度相对应的图像分割阈值确定。有益效果是采用该方法实现了三维参数向二维化的科学转变，排除了人为主观因数的干扰，使得土体微结构的研究更加准确。但是该发明的不足之处在于：该发明中leicaqwin图像分析软件对sem照片提取出的阈值，虽然进行了人为确认，但是在进行人为确认的过程中，涉及多次运算，过程繁琐的同时还易出错，从而不利于长期使用。

技术实现要素：

发明目的：针对现有岩土体材料数字图像二值化中分割阈值的确定过程复杂易出错的问题，本发明提出一种岩土体材料数字图像的分割阈值确定方法。

技术方案：为实现本发明的目的，本发明所采用的技术方案是：

一种岩土体材料数字图像的分割阈值确定方法，所述方法包括如下步骤：

s1：获取岩土体材料的sem图像的灰度直方图曲线；

s2：根据所述灰度直方图曲线，确定分割阈值t的取值范围；

s3：获取所述灰度直方图曲线的二阶导数；

s4：根据所述灰度直方图曲线的二阶导数和分割阈值t的取值范围，确定所述分割阈值t的大小。

更进一步地，所述步骤s1获取图像的灰度直方图曲线之前还包括：读取所述岩土体材料的sem图像，获取所述sem图像中每个像素的灰度级i、所述每个像素的灰度级i所对应的总像素数量ni。

更进一步地，所述步骤s1获取图像的灰度直方图曲线具体如下：

s1.1：确定所述岩土体材料灰度图像中每个像素的灰度级i；

s1.2：根据如下公式，获取所述岩土体材料灰度图像的灰度直方图曲线上各个点：

其中，i表示灰度级，n表示图像像素总数，ni表示图像中所有灰度级为i的像素的总像素个数，l表示灰度级的种类数；

s1.3：根据所述p(i)，拟合获取岩土体材料灰度图像的灰度直方图曲线。

更进一步地，所述步骤s2确定分割阈值t的取值范围具体如下：

步骤s2.1：根据所述灰度直方图曲线，确定所述灰度直方图曲线中的峰值数目；

步骤s2.2：根据所述峰值数目，确定所述岩土体材料的结构；

步骤s2.3：获取所述峰值对应像素的灰度级imax；

步骤s2.4：根据所述岩土体材料的结构和峰值对应像素的灰度级imax，确定所述分割阈值t的取值范围。

更进一步地，步骤s2.2确定所述岩土体材料的结构具体如下：

如果所述灰度直方图曲线只有一个峰值，所述岩土体材料的结构为含孔隙结构；

如果所述灰度直方图曲线有两个峰值，所述岩土体材料的结构为含裂缝结构。

更进一步地，步骤s3获取所述灰度直方图曲线的二阶导数之前包括：获取所述灰度直方图曲线的一阶导数，所述一阶导数为:

其中，i表示灰度级，ni表示图像中每个像素的灰度级i所对应的总像素数量。

更进一步地，确定所述孔隙结构的分割阈值t具体如下：

sa4.1：所述灰度直方图曲线只有一个峰值，所述分割阈值t的取值范围为：

t＜imax

其中，imax为峰值对应像素的灰度级；

sa4.2：获取所述灰度直方图曲线的二阶导数，所述二阶导数为：

其中，i表示灰度级，ni表示图像中每个像素的灰度级i所对应的总像素数量；

sa4.3：根据所述二阶导数，确定在分割阈值t的取值范围内二阶导数的最大值aimax；

sa4.4：确定所述二阶导数最大值aimax所对应像素的灰度级it,所述分割阈值t为：

t＝it

其中，it为二阶导数最大值aimax所对应像素的灰度级。

更进一步地，确定所述裂缝结构的分割阈值t具体如下：

sb4.1：所述灰度直方图曲线有两个峰值，所述分割阈值t的取值范围为：

imax1＜t＜imax2

其中，imax1为第一个峰值对应像素的灰度级，imax2为第二个峰值对应像素的灰度级；

sb4.2：获取所述灰度直方图曲线的二阶导数，所述二阶导数为：

其中，ni表示图像中所有灰度级为i的像素的总像素个数，i表示灰度级；

sb4.3：根据所述二阶导数，确定在分割阈值t的取值范围内二阶导数的最大值aimax；

sb4.4：确定所述二阶导数最大值aimax所对应像素的灰度级it,所述分割阈值t为：

t＝it

其中，it为二阶导数最大值aimax所对应像素的灰度级。

有益效果：与现有技术相比，本发明的技术方案具有以下有益技术效果：

(1)本发明先是通过待测岩土体材料的sem图像的灰度直方图曲线，确定分割阈值的取值范围，再通过灰度直方图曲线的二阶导数确定分割阈值的具体大小，通过不断的缩小范围直至精确取值，能够更进一步地保证分割阈值取值的准确性；

(2)本发明根据待测岩土体材料自身的sem图像进行分析，能够确保分割阈值更加满足待测岩土体材料数字图像二值化的需求，从而准确的从数字图像中将岩土体材料的孔隙或裂缝结构与表面土骨架结构区分开来；

(3)本发明为后续基于数字图像的岩土体材料的深入研究提供了准确的分割阈值，为准确的提取岩土体材料的孔隙或裂缝结构，提供了有效的技术支持。

附图说明

图1是本发明的流程示意图；

图2是膨润土sem图像；

图3是膨润土颗粒基本单元的示意图；

图4是sem扫描岩土体材料原理剖面示意图；

图5是灰度直方图曲线对应岩土体结构示意图；

图6是不同分割阈值压实膨润土二值图；

图7是sem图像随分割阈值增大的二值化提取过程；

图8是压裂后煤样sem图像；

图9是压裂后煤样灰度直方图曲线；

图10是膨润土的灰度直方图曲线；

图11是膨润土的灰度直方图曲线二阶导数曲线；

图12是膨润土的最终二值化图像；

图13是压裂后的煤样的灰度直方图曲线二阶导数曲线；

图14是压裂后的煤样的最终二值化图像。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述。其中，所描述的实施例是本发明一部分实施例，而不是全部的实施例。因此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。

参考图1，本实施例提供了一种岩土体材料数字图像的分割阈值确定方法，本方法计算得到的分割阈值能够有效地提取出sem图像中的孔隙结构或裂缝结构，同时也符合岩土体材料孔隙及裂缝的实际分布情况。从而为基于数字图像法的岩土体材料细观机理研究提供了良好的技术支持。

参考图4，图4为扫描电子显微镜通过用聚焦电子束扫描样品的表面来产生样品表面的图像，其中扫描电子显微镜英文为：scanningelectronmicroscope，可以缩写为sem，简称扫描电镜，是一种电子显微镜。

使用sem电镜扫描得到的岩土体材料表面微观结构图像为灰度图像，其中灰度图像的灰度级从0到255，总共有256种取值，同时sem图像中每个像素点颜色的深浅都代表着一个灰度级。

由于岩土体材料的孔隙结构或裂缝结构距离电子束发射口较远，因此最终的sem图像中代表孔隙或裂隙结构的像素点呈现出的灰度级一般在0～90；而岩土体材料的表面土骨架结构上的凸起颗粒距离较近，因此最终图像中代表岩土体材料表面土骨架结构的像素点呈现出的灰度级一般在150～255；同时其中岩土体材料的表面土骨架结构占大多数，具体大于50％，且处于同一平面中，因此其灰度级一般在90～150，对应像素数量最多。

在本实施例中，具体地讲，本分割阈值确定方法具体包括如下步骤：

步骤s1：通过matalb代码读取待测岩土体材料的sem图像，获取sem图像中每个像素的灰度级i以及每个像素的灰度级i所对应的总像素数量ni。

步骤s2：获取待测岩土体材料的sem图像的灰度直方图曲线，具体过程如下：

步骤s2.1：根据每个像素的灰度级i、每个像素的灰度级i所对应的总像素数量ni以及如下公式，获取待测岩土体材料灰度图像的灰度直方图曲线上的各个点：

其中，i表示灰度级，n表示图像像素总数，ni表示图像中所有灰度级为i的像素的总像素个数，l表示灰度级的种类数；

步骤s2.2：根据步骤s2.1中的各个p(i)，拟合获取岩土体材料灰度图像的灰度直方图曲线。

步骤s3：根据灰度直方图曲线，确定分割阈值t的取值范围，具体过程如下：

步骤s3.1：根据灰度直方图曲线，确定灰度直方图曲线中的峰值数目；

步骤s3.2：根据峰值数目，确定待测岩土体材料的结构，具体为：

如果灰度直方图曲线中只有一个峰值，待测岩土体材料的结构为含孔隙结构；

如果灰度直方图曲线中有两个峰值，待测岩土体材料的结构为含裂缝结构；

步骤s3.3：根据峰值数目，获取相应峰值对应像素的灰度级imax；

步骤s3.4：根据待测岩土体材料的结构、相应峰值对应像素的灰度级imax，确定分割阈值t的取值范围。

值得注意的是，由于待测岩土体材料结构的不确定，分割阈值t的取值范围的确定也是不固定的，由结构对应峰值的数目决定。

步骤s4：获取灰度直方图曲线的二阶导数，具体过程如下：

步骤s4.1：获取灰度直方图曲线的一阶导数，一阶导数为:

其中，i表示灰度级，ni表示图像中每个像素的灰度级i所对应的总像素数量；

步骤s4.2：获取灰度直方图曲线的二阶导数，二阶导数为：

其中，i表示灰度级，ni表示图像中每个像素的灰度级i所对应的总像素数量。

步骤s5：根据灰度直方图曲线的二阶导数和分割阈值t的取值范围，确定分割阈值t的大小，具体过程如下：

步骤s5.1：根据步骤s4.2中的各个ai，拟合获取二阶导数曲线；

步骤s5.2：确定在二阶导数曲线中，分割阈值t取值范围内的最大值aimax，同时明确最大值aimax所对应像素的灰度级it,分割阈值t为：

t＝it

其中，it为二阶导数最大值所对应像素的灰度级。

实施例1含孔隙结构的岩土体材料

本实施例提供了一种岩土体材料数字图像的分割阈值确定方法，在本实施例中，具体地讲，待测岩土体材料选择膨润土。其中膨润土数字图像的分割阈值确定方法，具体如下：

步骤sa1：参考图2，图2为膨润土的sem图像，通过matalb代码读取膨润土sem图像，获取膨润土sem图像中每个像素的灰度级i以及每个像素的灰度级i所对应的总像素数量ni。

步骤sa2：获取膨润土sem图像的灰度直方图曲线，具体过程如下：

步骤sa2.1：根据每个像素的灰度级i、每个像素的灰度级i所对应的总像素数量ni以及如下公式，获取膨润土灰度图像的灰度直方图曲线上的各个点：

其中，i表示灰度级，n表示图像像素总数，ni表示图像中所有灰度级为i的像素的总像素个数，l表示灰度级的种类数。

步骤sa2.2：根据步骤sa2.1中的各个p(i)，拟合获取膨润土灰度图像的灰度直方图曲线。参考图10，图10为膨润土的灰度直方图曲线，从膨润土的灰度直方图曲线中，能够直观的看到膨润土sem图像中每个灰度级i的像素数量的占比情况。

步骤sa3：根据灰度直方图曲线，确定分割阈值t的取值范围，具体过程如下：

步骤sa3.1：根据灰度直方图曲线，确定灰度直方图曲线中的峰值数目，在本实施例中，参考图10可知，灰度直方图曲线中的峰值只有一个，因此膨润土的结构为孔隙结构。

根据压汞法及岩心覆压孔渗自动测试系统测试得到岩土体的孔隙率小于30％，其与表面土骨架相比来说，孔隙结构占小部分。在本实施例中，具体地讲，根据压汞法及岩心覆压孔渗自动测试系统测试得到干燥压实膨润土的孔隙率为20％～30％，对于膨润土sem图像来说，其中的孔隙结构应该在25％左右。参考图3可知，膨润土的孔隙结构又分有几个不同的层次，通过sem电子镜扫描观察到的膨润土孔隙结构属于表观孔隙结构，故通过数字图像法得到的孔隙率应该小于20％。

步骤sa3.2：根据步骤sa3.1可知，膨润土的灰度直方图曲线中只有一个峰值，因此在本实施例中，只需要确定一个峰值对应像素的灰度级imax。

步骤sa3.3：根据膨润土的结构、相应峰值对应像素的灰度级imax，确定分割阈值t的取值范围，因为膨润土的结构为孔隙结构，且其灰度直方图曲线只有一个峰值，故膨润土数字图像的分割阈值t的取值范围为：

t＜imax

其中，imax为峰值对应像素的灰度级。

参考图5，由于所有含孔隙结构的岩土体材料的灰度直方图都是单峰曲线，所以膨润土的结构为孔隙结构。其中由于表面土骨架占大多数，因此峰值处对应像素一定属于表面土骨架。值得注意的是，通过压汞法及岩心覆压孔渗自动测试系统测试得到干燥压实膨润土的孔隙率小于30％，因此，膨润土的数字图像中代表孔隙结构的像素也肯定小于30％，故膨润土中的表面土骨架结构超过70％，因此灰度直方图中峰值部分对应像素肯定是代表表面土骨架。显然地，孔隙结构对应的灰度级一般在0～90，而剩下的在90～225的灰度级对应的则为表面土骨架结构，故代表孔隙结构像素的灰度级小于代表表面土骨架像素的灰度级。

参考图6，图6为随着分割阈值t从0逐渐增大到255时，压实膨润土二值图的变化过程。其中，图6中黑色部分就是根据不同的分割阈值提取出的孔隙结构，显然有一些是非常不合理的。但是通过这个变化过程，可以观察到随着分割阈值t的逐渐增大，孔隙率也在逐渐增大。

参考图7，图7为二值化图随着分割阈值t的变化过程，其中二值化图像孔隙率的增大实际上就是图中黑色像素数量随着分割阈值t的增大而增大。当分割阈值t从孔隙结构的灰度级到达表面土骨架灰度级时，由于表面土骨架结构占大多数，故图像中每个像素的灰度级i所对应的总像素数量ni会存在一个突变，因此，找到这个从孔隙结构到表面土骨架结构突变的这个点对应的灰度级i，即确定分割阈值t的大小。

步骤sa4：获取灰度直方图曲线的二阶导数，具体过程如下：

步骤sa4.1：获取灰度直方图曲线的一阶导数，一阶导数为:

其中，i表示灰度级，ni表示图像中每个像素的灰度级i所对应的总像素数量。

步骤sa4.2：获取灰度直方图曲线的二阶导数，二阶导数为：

其中，i表示灰度级，ni表示图像中每个像素的灰度级i所对应的总像素数量。

步骤sa5：根据灰度直方图曲线的二阶导数ai和分割阈值t的取值范围，确定分割阈值t的大小，具体过程如下：

步骤sa5.1：根据步骤sa4.2中的各个ai，拟合获取二阶导数曲线，参考图11，图11为膨润土的灰度直方图曲线的二阶导数曲线。

步骤sa5.2：确定在分割阈值t取值范围内，二阶导数曲线中的最大值aimax，同时确定最大值aimax所对应像素的灰度级it,分割阈值t为：

t＝it

其中，it为二阶导数最大值所对应像素的灰度级。

参考图12，图12为膨润土最终的二值化图像。其中分割阈值t＝it，故在膨润土的灰度直方图曲线中，灰度级i小于分割阈值t处所对应的范围为膨润土的孔隙结构，灰度级i大于分割阈值t处所对应的范围为膨润土的表面土骨架结构。

实施例2含裂缝结构的岩土体材料

本实施例提供了一种岩土体材料数字图像的分割阈值确定方法，在本实施例中，具体地讲，待测岩土体材料选择压裂后的煤样。其中压裂后煤样的数字图像的分割阈值确定方法，具体如下：

步骤sb1：参考图8，图8为压裂后煤样的sem图像，通过matalb代码读取压裂后煤样的sem图像，获取压裂后煤样sem图像中每个像素的灰度级i以及每个像素的灰度级i所对应的总像素数量ni。

步骤sb2：获取压裂后煤样sem图像的灰度直方图曲线，具体过程如下：

步骤sb2.1：根据每个像素的灰度级i、每个像素的灰度级i所对应的总像素数量ni以及如下公式，获取压裂后煤样灰度图像的灰度直方图曲线上的各个点：

其中，i表示灰度级，n表示图像像素总数，ni表示图像中所有灰度级为i的像素的总像素个数，l表示灰度级的种类数。

步骤sb2.2：根据步骤sb2.1中的各个p(i)，拟合获取压裂后煤样灰度图像的灰度直方图曲线。参考图9，图9为压裂后煤样sem图像的灰度直方图曲线。

步骤sb3：根据灰度直方图曲线，确定分割阈值t的取值范围，具体过程如下：

步骤sb3.1：根据灰度直方图曲线，确定灰度直方图曲线中的峰值数目，在本实施例中，参考图9可知，灰度直方图曲线中的峰值有两个，因此压裂后煤样的结构为裂缝结构。

步骤sb3.2：根据步骤sb3.1可知，压裂后煤样的灰度直方图曲线中有两个峰值，因此在本实施例中，需要确定两个峰值分别对应像素的灰度级imax1和imax2，其中，imax1为第一个峰值对应像素的灰度级，imax2为第二个峰值对应像素的灰度级。

步骤sb3.3：根据压裂后煤样的结构、相应峰值对应像素的灰度级imax1和imax2，确定分割阈值t的取值范围，因为压裂后煤样的结构为孔隙结构，且其灰度直方图曲线有两个峰值，故压裂后煤样数字图像的分割阈值t的取值范围为：

imax1＜t＜imax2

其中，imax1为第一个峰值对应像素的灰度级，imax2为第二个峰值对应像素的灰度级。

参考图9，含裂缝结构的岩土体材料的灰度直方图曲线均为双峰曲线，故压裂后煤样的结构为含裂缝结构的岩土体材料，在灰度级i从小到大变化的过程中，图像中第一个峰值处的像素代表着裂缝处的像素，而第二个峰值处的像素代表着表面骨架处的像素，因此含裂缝结构的岩土体材料的分割阈值t介于两个峰值之间。

步骤sb4：获取灰度直方图曲线的二阶导数，具体过程如下：

步骤sb4.1：获取灰度直方图曲线的一阶导数，一阶导数为:

其中，i表示灰度级，ni表示图像中每个像素的灰度级i所对应的总像素数量。

步骤sb4.2：获取灰度直方图曲线的二阶导数，二阶导数为：

其中，i表示灰度级，ni表示图像中每个像素的灰度级i所对应的总像素数量。

步骤sb5：根据灰度直方图曲线的二阶导数ai和分割阈值t的取值范围，确定分割阈值t的大小，具体过程如下：

步骤sb5.1：根据步骤sb4.2中的各个ai，拟合获取二阶导数曲线，参考图13，图13为压裂后煤样的灰度直方图曲线的二阶导数曲线。

步骤sb5.2：确定在分割阈值t取值范围内，二阶导数曲线中的最大值aimax，同时确定最大值aimax所对应像素的灰度级it,分割阈值t为：

t＝it

其中，it为二阶导数最大值所对应像素的灰度级。

参考图14，图14为压裂后煤样最终的二值化图像。其中分割阈值t＝it，故在压裂后煤样的灰度直方图曲线中，灰度级i小于分割阈值t处所对应的范围为压裂后煤样的裂缝结构，而灰度级i大于分割阈值t处所对应的范围则为压裂后煤样的表面土骨架结构。

以上示意性的对本发明及其实施方式进行了描述，该描述没有限制性，附图中所示的也只是本发明的实施方式之一，实际的结构和方法并不局限于此。所以，如果本领域的普通技术人员受其启示，在不脱离本发明创造宗旨的情况下，不经创造性的设计出与该技术方案相似的结构方式及实施例，均属于本发明的保护范围。