一种岩石块度尺寸分布的测量方法与流程

本发明属于岩石块度检测技术领域，尤其涉及一种岩石块度尺寸分布的测量方法，适用于运动的或者爆堆的岩石块度图像的测量和分析。

背景技术：

在采石和矿业工程中，对岩石块度的尺寸分布的测量是非常重要的。石料就是自然岩块和爆破及机械破碎的岩块的混合体，主要是用于建筑，公路，铁路和大坝等。为了判断石料的质量，对石料颗粒的尺寸和形状参数进行估算是必要的。石料的尺寸分布不但是用来评估产品质量的一个数据，而且还是调整破碎机或爆破生产的重要信息，例如：在采石生产中，调节破碎机间隙及在矿业工程中调节打孔的孔径等。破碎机通常被设定用来生产某个严格指定的相对较窄尺寸范围内的石料，比如从16mm到30mm。通常破碎机操作的一个主要指标就是平均尺寸。在自动破碎控制系统中，从实时系统发回的包括平均石料尺寸的反馈信号，就显示了流水线上破碎过程的实际进展。在实际应用中，从破碎机出来的破碎颗粒在一条传送带上传输，在其上方放置一个ccd摄像头向下拍摄，然后用图像处理、分割和分析对获取的图像中的颗粒进行测量。

在矿业生产中，块度尺寸分布是衡量整个采矿生产成本的关键，几乎可以影响整个采矿成本的30％。如果块度尺寸过低，岩石爆破的成本就会提高；相反，如果尺寸过大，不仅要进行二次爆破会增加费用，而且货车装载、卡车及火车运输的费用都会增高。因此，岩石爆破的平均尺寸是使矿业生产达到最优化的一个重要因素，是控制整个采矿生产过程的重要信息。

人工测量劳动强度大，采样量有限，而且很难进行实时测量和分析，从而影响生产效率。图像处理可以克服此弊病，但是岩石块度图像是最复杂的颗粒图像，原因是：岩石块度表面粗糙，有明显的三维体信息，不同的颜色及相互粘连和重叠等会给图像的分割带来很大的困难。

技术实现要素：

针对上述现有技术的缺点，本发明的目的在于提供一种岩石块度尺寸分布的测量方法，有利于快速准确的检测出岩石块度图像中目标的尺寸分布。

为达到上述目的，本发明采用如下技术方案予以实现。

一种岩石块度尺寸分布的测量方法，所述方法包括如下步骤：

步骤1，获取岩石块度的彩色rgb图像，将所述彩色rgb图像转换为灰度图像；

步骤2，对所述灰度图像进行分数阶积分平滑滤波，得到平滑后的图像；

步骤3，对所述平滑后的图像进行边界扫描，得到二值边界图像；

步骤4，对所述二值边界图像进行骨架提取，得到目标骨架图像；

步骤5，对所述目标骨架图像进行三次循环的短线去除及断线连接，得到结果图像；

步骤6，对所述结果图像进行毛刺去除及断点连接，得到标号图像；所述标号图像中包含多个目标块度；

步骤7，获取所述目标标号图像中每个目标块度的最小外接矩形，并计算每个目标块度的面积和长度；所述每个目标块度的面积为其最小外接矩形的面积，所述每个目标块度的长度为其最小外接矩形的长度；

步骤8，对所有目标块度的面积和长度进行统计，得到岩石块度尺寸分布。

本发明技术方案的特点和进一步的改进为：

(1)步骤1中，

获取岩石块度的彩色rgb图像f(x，y)rgb，将所述彩色rgb图像f(x，y)rgb转换为灰度图像f(x，y)＝max(f(x，y)r，f(x，y)g，f(x，y)b)；

其中，(x，y)表示彩色rgb图像中像素点的位置，其与灰度图像中像素点的位置一一对应，f(x，y)r表示像素点(x，y)处的r值，f(x，y)g表示像素点(x，y)处的g值，f(x，y)b表示像素点(x，y)处的b值。

(2)步骤2中，

对所述灰度图像进行分数阶积分平滑滤波，得到平滑后的图像f(x，y)＝(f(x，y)*h)/8，其中，f(x，y)表示灰度图像，h表示对所述灰度图像进行分数阶积分平滑滤波的滤波器；

对所述灰度图像进行分数阶积分平滑滤波的滤波器h为5×5的滤波器，具体为：

(3)步骤3中，

对所述平滑后的图像f(x，y)进行canny边界扫描，得到二值边界图像g(x，y)，边界为白色。

(4)步骤5中，对所述目标骨架图像g1(x，y)进行三次循环的短线去除及断线连接，得到结果图像g2(x，y)，具体包括：

设置第一阈值，第二阈值以及第三阈值，所述第一阈值小于所述第二阈值，所述第二阈值小于所述第三阈值；

采用第一阈值过滤掉所述目标骨架图像中的噪声线条和噪声点，并进行边界间隙的连接，得到第一中间图像；

采用第二阈值过滤掉所述第一中间图像中的噪声线条和噪声点，并进行边界间隙的连接，得到第二中间图像；

采用第三阈值过滤掉所述第二中间图像中的噪声线条和噪声点，并进行边界间隙的连接，得到第三中间图像，作为结果图像g2(x，y)。

(5)步骤8中，对所有目标块度的面积和长度进行统计，得到岩石块度尺寸分布，具体包括：

计算所有目标块度的面积之和，并获取所有目标块度的长度中的最大值imax和最小值imin，将岩石块度尺寸分布划分为10个尺寸等级，相邻两个尺寸等级之间的间隔为：(imax-imax)/10；

根据岩石块度尺寸分布的10个尺寸等级，得到每个目标块度所属的尺寸等级；分别将属于同一个尺寸等级的目标块度的面积进行累加，得到每个尺寸等级包含的目标块度的总面积；

根据每个尺寸等级包含的目标块度的总面积与所有目标块度的面积之和的比值，得到每个尺寸等级的累积百分比；

以10个尺寸等级为x轴，每个尺寸等级的累积百分比为y轴，绘制岩石块度尺寸分布的二维分布图，得到岩石块度尺寸分布。

本发明提供的一种通过提取不完整矿岩块度轮廓来估算块度尺寸分布的方法，通过将彩色图像转换为灰度图像；进行分数阶积分平滑滤波；对平滑后的图像进行canny边界扫描；进行骨架提取；进行三次循环的短线去除及断线连接；进行毛刺去除及最后的断点连接；最小外接矩形获取及块度长度的测量；尺寸分布统计和计算，有利于快速准确地检测出复杂图像中颗粒的尺寸分布，该方法适用于实时在线的运动颗粒视频图像处理。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明实施例提供一种岩石块度尺寸分布的测量方法的流程示意图；

图2为本发明实施例提供最小外界矩形的示意图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本发明实施例提供一种岩石块度尺寸分布的测量方法，如图1所示，所述方法包括如下步骤：

步骤1，输入rgb彩色岩石块度图像f(x，y)rgb：目前市面上的相机及视频相机都是彩色的，所以以彩色图像为例，图像可以来自运输带上的运动矿岩块度图像也可以是来自矿山爆堆或是产状斗里的图像。

步骤2，将彩色图像转换为灰度图像：由于岩石料堆或在运输带上都呈暗色居多，故选择rgb三频道里最亮的点作为灰度图像相应的点，这样能增强图像的亮度，便于后续的处理，故此如下采用公式获得灰度图像f(x，y)＝max(f(x，y)r，f(x，y)g，f(x，y)b)，f(x，y)r，f(x，y)g，f(x，y)b分别为同一场景的红绿蓝三频道的图像。

步骤3：进行分数阶积分平滑滤波：去除所述颗粒图像中的噪声，该过程可表示为：

由差分表达式，便可以得到分数阶微分的差分系数为：

a0＝1，a1＝v；a2＝＝(v-v²)/2；a3＝(v³-3v²+2v)/6

为了充分利用每一个邻域像素的作用，模板中的每一个像素权值都不应该为零，而权值的大小根据邻域点与中心点的距离及分数阶微分结果级数展开式中的项数系数来确定。

考虑到距离的计算较为繁杂，为了简化这种繁杂的计算，微小的距离差异可以忽略不计，如：距离(两相邻像素的斜线距离)。因此有：

a1-2＝(3v-v²)/4；a2-3＝(v³-6v²+5v)/12；a3-4＝-(v⁴-10v³+23v²-14v)/48

根据上述分析，5x5模板的滤波器可以写成(a＝2)：

f(x，y)＝(f(x，y)*h)/8，其中h表示平滑滤波器。

步骤4：对平滑后的图像进行canny边界扫描：得到二值边界图像g(x，y)，边界为白色，具体过程如下：

对平滑后的图像f(x，y)进行canny边界扫描后的梯度图像为g(x，y)，高低阈值采取最大熵值得方法求得。

使用类之间的最大信息交叉熵来获得双阈值，最后通过一个粗略道路追踪程序跟踪目标边缘。使用贝叶斯和交叉熵理论确定梯度幅值图像的阈值。图像被划分成两大类：目标(o)及背景(b)，一幅图像应该有两个正态分布，其中的参数可以从灰度图像的直方图得到：

其中，p(g/i)为i条件下的概率密度，t是阈值，g是灰度值，μ和σ分别为正太分布的均值和方差。

两个类别的方差估计如下：

其中，目标类的先验概率是背景类的先验概率是因此它们的类间平均值分别是知其中，t是阈值，g是灰度值，l为灰度上界。

后验概率由贝叶斯概率公式获得：

通过不同区域中的像素的最大后验概率获得最优阈值。

基于单个像素后验概率的类间交叉熵为：

据此，可获得类间差。为了简化计算，以灰度值g替换像素灰度尺度s，这样，用灰度直方图替换了概率分布。它可以被重新写为：

式中，l为灰度上界，t是一个灰度阈值。

要取得基于最大交叉熵类之间的最佳的阈值t^*，可以通过一个搜索操作进行：

是g的梯度向量，

是平滑图像g的梯度图像。平滑参数σgauss即所谓的滤波器尺寸参数，最后根据所求的t值得到二值边界图像g(x，y)。

为了避免canny的双边界和不连续性，这里附加了动态阈值处理。那么在离散的栅格中，通过边界的细线化，边缘像素为：在二值图像中，0就表示边缘像素，非0则表示非边缘像素，记边界图像为εg(x，y)，或者更精确地表示为εg(x，y；σgauss)。边缘密度总是根据这个边缘图像εg的ne值计算出来的。

步骤5，对二值边界图像g(x，y)进行骨架提取：模拟烧草模型的过程，由目标的边界开始逐渐向内演化，逐步搜索到目标中轴的位置。其主要思想是均匀地逐层剥掉目标的边界点，剩下最里层保留(否则会影响连通性)的部分就得到了目标的骨架，从而得到骨架图像g1(x，y)。

步骤6，对骨架图像g1(x，y)进行三次循环的短线去除及断线连接：用第一阈值3滤掉尽可能多的噪声线条及噪声点，再缝合边界间隙，然后目标细线化。而后用第二阈值5滤掉尽可能多的噪声线条及噪声点，缝合边界间隙。最后用第三阈值7滤掉尽可能多的噪声线条及噪声点，再缝合边界间隙，得到结果图像g2(x，y)。

步骤7，对结果图像g2(x，y)进行毛刺去除及最后的断点连接：为了去除毛刺，进行交叉点检测，然后去除交叉点将一个目标分成若干段，再进行短线删除，在讲交叉点放回原处连接了两段长线，在每个目标的两端搜索附近是否有其他线段，若有，进行断点连接，完成所有目标操作后，做图像标号，得到标号图像g3(x，y)。

步骤8，最小外接矩形获取及块度长度的测量：在标号图像g3(x，y)中，对每一个目标进行最小外接矩形获取及块度长度的测量，具体过程如下：

对于上述处理后的标号图像g3(x，y)，定义每个目标的矩公式为：

由此得到质心

再根据下式得到的θ值，可以确定图像的主轴与次轴方向，主轴即为通过物体质心的直线，次轴方向即为通过物体的质心同时垂直于主轴线的方向。

如图2所示，在图形上做两条与主轴方向平行的直线l1和直线l2，将这两条直线分别沿垂直于主轴的方向即次轴方向不断平移，直到两条平移直线与图形轮廓线最外面的点相交，即下一次平移时直线与图形将无交点。同理，做两条平行于次轴方向的直线l3和直线l4，将这两条直线分别沿垂直于次轴的方向即主轴方向不断向外平移，直到两条平移直线与图形轮廓线最外面的点相交为止，最后这四条相交的直线构成一个矩形来确定ferret的边界。

按上面的法则，对标号图像g3(x，y)中的k个目标(块度)，计算每一个目标的面积sj长度ij，j＝1，2，...，k。

步骤9，尺寸分布统计和计算：在长度ij(j＝1，2，...，k)里找出最大值imax和最小值imin，将岩块分成10个等级，等级间隔为：(imax-imax)/10，对每一个等级的岩块，累加所有岩块的面积sj，最后以面积为权重，计算块度尺寸的累计分布代替实际的矿岩块度尺寸分布。在xy直角坐标系中，x轴是尺寸，y轴表示块度的累积百分比，当x轴上在最大尺寸时，y轴上的累积百分比是100％。

本领域普通技术人员可以理解：实现上述方法实施例的全部或部分步骤可以通过程序指令相关的硬件来完成，前述的程序可以存储于计算机可读取存储介质中，该程序在执行时，执行包括上述方法实施例的步骤；而前述的存储介质包括：rom、ram、磁碟或者光盘等各种可以存储程序代码的介质。

以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以所述权利要求的保护范围为准。