煤矸石粒度分布检测中基于凸包分析的图像分割研究

孙颖,田志辉,张雷,2*,李丽锋,赵保国,马维光,2,董磊,2,尹王保,2*,肖连团,2,贾锁堂,2

(1.山西大学 激光光谱研究所,量子光学与光量子器件国家重点实验室,山西 太原 030006;2.山西大学 极端光学协同创新中心,山西 太原 030006;3.山西河坡发电有限责任公司,山西 阳泉 045000)

0 引言

煤矸石是在掘进、开采和洗煤过程中排出的固体废料,是目前我国排放量最大的工业固体废弃物之一,但同时也是可利用的宝贵资源,能用于铺路、烧制砖瓦、供暖、发电等。在应用于循环流化床(CFB)锅炉发电时,对煤矸石的粒度分布要求严格,如果无法满足锅炉的设计要求,将会造成锅炉负荷的不稳定[1]。因此,实时检测煤矸石的粒度分布并及时反馈调整破筛机参数,对于CFB锅炉的稳定运行至关重要。

传统的煤矸石粒度分布检测主要采用筛分法,即先将样品经系列标准筛分离成若干个粒级,再分别称重后获得粒度分布。该方法虽然准确,但耗时太长(1～2 h),严重滞后于工业电力生产。相比之下,基于机器视觉图像处理的粒度分布检测方法具有操作简便、实时在线的优势[2],目前国际上已出现了商用化产品,如加拿大的WipFrag、法国的VisioRock[3-4]等,但售价昂贵、维护费用高,难以在我国电厂普遍推广。因此,开发低成本的、具有自主知识产权的机器视觉粒度分布检测技术,实现对电厂皮带煤矸石粒度分布的实时在线检测,对于提升我国煤矸石的固废利用率以及提高低热值煤的发电效率意义重大。

目前,基于机器视觉的粒度分布检测的关键技术瓶颈是如何实现颗粒图像的精确分割。针对该问题,国际上已开展了大量的相关研究,如英国诺丁汉大学的Koh等人[5-6]利用多闪光成像的方法捕获在不同照明方向时的样品图像,通过阴影得到颗粒周围的边缘信息,从而减少图像分割的误差,但该方法需要按序对样品进行多角度投影,难以应用于皮带上的移动样品;瑞典吕勒奥理工大学的Thurley等人[7-8]利用激光三角测距的机器视觉方法采集颗粒的三维数据,实现颗粒图像的精确分割,但需额外增加激光器、高速相机等设备,成本较高。除了以上技术方案外,人们也在不断改进分割算法,如加拿大阿尔伯塔大学Saha等人[9]使用Snake-PCA算法精确分割图像颗粒,评估油砂的粒度分布;加拿大曼尼托巴大学的Mebatsion等人[10]提出了基于傅里叶级数近似的分割算法,利用边界曲率成功分割了接触颗粒;加拿大劳伦森大学的Ko等人[11-12]采用了基于神经网络的分割算法,通过提取均匀性信息来建立非线性模型,得到了更加可靠的矿石粒度分布数据。

我们提出了采用基于Hessian矩阵[13]的多尺度滤波器提取颗粒边缘,并结合距离变换分水岭和凸包分析来快速、精确地对煤矸石图像进行二次颗粒分割。

1 实验装置

基于机器视觉的粒度分布检测实验装置结构如图1所示,采用了配有8 mm焦距镜头的工业相机(大恒,分辨率1 964×2 612,曝光时长5 000 μs)和LED灯条(发光功率100 W,灯珠210颗),置于装置的顶部。待测煤矸石样品均匀平铺于实验装置底部,样品表面距相机镜头0.5 m,所采集的图像由USB 3.0数据线传送至计算机进行分析。自行开发的分析软件前面板如图2所示,主体基于LabVIEW平台编制,同时调用了Matlab的图像处理工具箱。软件界面包括采集的原始图像、颗粒边缘提取图像、分割图像、颗粒识别图像和粒度分布图。煤矸石样品采集于格盟国际山西河坡发电有限责任公司的350 MW超临界热电联产1#机组的输煤皮带。

Fig.1 Schematic of the experimental set-up for coal gangue particle size distribution analysis图1 煤矸石粒度分布检测实验装置结构

Fig.2 Front panel of the coal gangue particle size distribution analysis written by LabVIEW图2 基于LabVIEW的煤矸石粒度分析界面

2 图像处理

2.1 图像预处理

为了提高原始图像的对比度,改善图像质量,在特征提取之前,对煤矸石图像进行预处理。先利用Matlab的图像处理工具箱将彩色图像转化为灰度图像,再使用自适应直方图均衡化方法增强图像局部对比度。该算法通过计算图像中多个局部区域的直方图,重新分配图像亮度,改善图像局部对比度。预处理直接影响后续的处理效果,是图像处理中的关键一步。

2.2 图像分割

实验采用了组合式的图像分割方法,即先用基于Hessian矩阵的多尺度滤波器提取预处理后的颗粒边缘,再利用距离变换分水岭和凸包分析分别进行一次和二次分割。具体的图像处理方法如下。

2.2.1 基于多尺度Hessian矩阵的图像边缘提取

为获取图像中的颗粒信息,首先要提取图像中的颗粒边缘。Hessian 矩阵是一个多元函数的二阶偏导数构成的实对称矩阵,描述了函数的局部曲率,常被用于医学血管图像[14]、路面裂缝[15]等线状结构检测。将一幅二维图像各像素点的强度值定义为函数I(x,y),则相应的Hessian矩阵可表示为:

(1)

式中,Ixx(x,y)、Iyy(x,y)分别为I(x,y)在x方向、y方向上的二阶偏微分,Ixy(x,y)、Iyx(x,y)分别为I(x,y)在x、y方向和y、x方向上的混合偏微分,反映该点的局部二阶结构。Hessian矩阵边缘检测滤波器是计算某一像素点的Hessian矩阵的两个特征值λ1和λ2,它们分别代表该像素点的最大和最小曲率。假设|λ1|≤|λ2|,则Hessian矩阵的两个特征值有如下规律:如果该点属于灰度分布均匀的图像背景点时,则λ1、λ2的绝对值均很小;如果该点是孤立点,则在该点的任意方向曲率都很大,即λ1、λ2的绝对值均偏大;而当该点属于图像边缘点时,沿边缘方向曲率非常小,因此有|λ1|≈0,而垂直于边缘方向曲率最大,则λ2的绝对值偏大。因此为提取图像中的颗粒边缘,构造如下函数:

λmax=max(|λ1|,|λ2|) .

(2)

根据该函数计算图像中每个点的Hessian矩阵的特征值,并选择绝对值大的特征值λmax作为所求边缘图像在该点处的强度值,最后所有点的λmax值共同组成颗粒的边缘图像。

由于采集的图像中颗粒间隔的宽度不一,构造多尺度的边缘检测滤波器更适用于灵活地提取颗粒边缘。因此,我们将上述Hessian矩阵边缘检测滤波器与高斯函数相结合[16],构造多尺度边缘滤波器。根据二阶导数的卷积性质可知,与高斯函数结合后每点的Hessian矩阵中各元素均可写为:

∂2I(x,y;σ)=I(x,y)⊗σ2∂2G(x,y;σ) .

(3)

式中,G为高斯函数,可表示为:

(4)

其中,σ为多尺度滤波器的尺度因子,通过对每个(x,y)点迭代尺度因子σ,得到最大响应值f(x,y)作为所求边缘图像在该点处的强度值,响应输出函数f(x,y)可表达为:

f(x,y)=max(λmax(σ)),σ∈[σmin,σmax] .

(5)

迭代尺度因子的变化范围为[σmin,σmax],为获得更好的增强效果,记实际颗粒间隙范围为[d0,d1],尺度因子σ合适的尺度范围[14]选取为[d0/4,d1/4],则σmin=d0/4,σmax=d1/4。迭代步长为(σmin-σmax)/n,n为正整数,n值越大,则迭代次数越多,边缘的细节提取也越多,计算量也越大。对于给定尺度σ的响应函数,σ越小,对宽度细小的边缘的提取效果越好,σ越大,对宽度较大的边缘的提取效果越好。

2.2.2 基于距离变换分水岭的一次分割

用双阈值将所获边缘图像转换为二值图,其中高阈值通过最常用的最大类间方差法获得,低阈值设置为高阈值的一半;再进行取反、填充、开、闭运算等形态学图像处理,消除颗粒表面纹理产生的伪边缘;最后使用距离变换,即计算每个像素到最近非零值像素的距离,并对距离图像进行分水岭变换[17]。分水岭算法将图像看作测地学上的拓扑表面,通过寻找图像中代表颗粒和边缘的汇水盆地和脊线,从而提取堆积颗粒的边界。因此实验使用分水岭算法完成对颗粒图像的一次分割。

2.2.3 基于凸包分析的二次分割

一次分割图像仍会存在少量欠分割颗粒,这些欠分割的特征参数凸度率[18]均小于已分割颗粒,因此实验设计了一种凸包分析算法对图像进行二次分割。每个颗粒的凸度率可以表示为:

(6)

式中,S为颗粒原始面积,Sc为凸包面积。凸包[19-20]是包含了颗粒的边缘和边缘内所有点的最小多边形,常用于事物形状的表征,因此凸包面积Sc大于或等于颗粒原始面积S,凸度率范围在0和1之间。在进行实验时发现,需进行二次分割的颗粒边缘凹陷较明显,颗粒原始面积远小于它的凸包面积,凸度率的值较低。实验通过检查每个区域的特征(即凸度率),根据凸包区域和原始区域之间的差异(即如果区域低于某一阈值,则该区域需要二次分割),来判别需要二次分割的区域。实际操作中确定凸度率阈值ε0,对阈值以下的颗粒进行如下二次分割:求取每个原始颗粒的凸包图像,将原始颗粒与凸包图像相减,找到凸包中不属于原始颗粒的区域,消除过小的连通区域后,由于这些较小的区域通常导致无效的分割,得到剩余的多个连通片段的二值图;求得剩余连通区域的凸包点,计算连通区域间的所有凸包点间距;在原始图像上用暗线连接最小间距的两个凸包点,完成对颗粒图像的二次分割。

2.3 粒度分布计算

为了获取颗粒的粒度分布,需要对分割后二值图中的颗粒进行识别与分析。利用Matlab图像处理工具箱识别二值图像中代表颗粒的连通区域,提取颗粒的平均Feret直径[21],从而计算相应的粒度分布结果。

3 实验结果与分析

利用本煤矸石粒度分布检测实验装置采集的煤矸石原始图像如图3a。为了获得颗粒的粒度分布信息,首先对图像进行预处理,即对原始图像进行自适应直方图均衡化处理,如图3b,图像的对比度相比原始图像有所增强,更有利于后续的边缘提取;然后利用多尺度Hessian矩阵边缘检测滤波器提取颗粒边缘,通过实验测试和计算,选取尺度因子σ范围是[1.5, 9.5],步长是1,所获边缘图像如图3c;最后对图像继续进行阈值处理、形态学开、闭运算以及一次和二次分割,所获可识别颗粒的分割情况如图3d和3e。

a) original image; b) pre-processed image; c) multi-scale filtered edge image;d) segmentation image; e) identified particle imageFig.3 Image processinga)原始图像;b)预处理图像;c)多尺度滤波后边缘图像;d)分割图像;e)识别颗粒图3 图像处理

颗粒的凸包二次分割过程如图4所示。其中,图4a是原始颗粒的二值图,将其与相应凸包图像(图4b)相减,得到包含多个连通区域的二值图像,消除过小区域后得到图4c。求得剩余连通区域的凸包点,计算区域间的所有凸包点间距,将最小间距的两凸包点在原始图像上用暗线相连,实现对粘连颗粒的分割,结果如图4d。在凸包分析中,图像中的欠分割颗粒的凸度率均小于0.75,因此设置阈值ε0=0.75。当某颗粒凸度率小于该阈值时,就对其进行二次分割。图4e、4f为图片中同一区域的一次、二次分割的局部放大对比图,其中经二次分割的颗粒圈在了红色方框中。

a) original image; b) convex hull image; c) subtracted image; d) segmentation results;e) partially enlarged image of the primary split; f) partially enlarged image of the secondary splitFig.4 Convex hull segmentationa)原始图像;b)凸包图像;c)差值图像;d)分割结果；e)一次分割的局部放大;f)二次分割的局部放大图4 凸包分割

为验证本粒度分布图像处理方法的精确性和重复性,实验采用了六组(a、b、c、e、f、g组)粒度分布不同的煤矸石样品,将各组样品颗粒平铺于实验检测装置底部,每组样品采集10张图像(每次采集前都将颗粒进行重新随机分散),并用分水岭和凸包分析分割算法进行一次和二次分割。实验所获在<1 mm、1～3 mm、3～5 mm、5～8 mm、>8 mm粒级上的颗粒粒度分布与标准筛分法所测结果的对比如图5。值得说明的是,该处的筛分法是指采用标准筛对样品进行逐级筛分后,再用图像识别法对各级分散颗粒进行计数。由图5可知,基于机器视觉的图像处理方法基本能客观反映样品的粒度分布趋势,与常规的一次分割相比,基于凸包分析的二次分割方法将各粒级粒度分布的平均绝对测量误差由3.78%降至2.23%,同时通过对比标准偏差(图中竖直短线)可知粒度分布的平均不确定度也由2.27%降至1.91%。可见,采用基于凸包分析的二次分割方法能够有效提高粒度分布的检测精度和重复性。此外,经测试本算法的平均分析时长为34 s,能够进一步应用于工业现场实现在线检测。

Fig.5 Comparison of particle size distribution results of three samples by traditional sieving, primary split and secondary split图5 分别利用筛分法、一次分割、二次分割所测六三组样品的颗粒粒度分布对比

4 结论

为了实现对低热值煤电厂中输煤皮带上煤矸石粒径的在线检测,以及时反馈指导调整破筛机工艺参数,提高矸石利用率和发电效率,课题组前期在实验室搭建了基于机器视觉的煤矸石粒度检测实验装置,开发了用于预测粒度分布的图像分析方法,重点解决图像颗粒分割的关键问题。该方法中使用多尺度Hessian矩阵边缘检测滤波器提取图像中的颗粒边缘,并结合距离变换分水岭算法分割图像中堆积的煤矸石颗粒,同时创新性地引入了基于凸包分析的二次分割算法,进一步校正了欠分割颗粒区域。基于三组煤矸石样品的实验测试结果表明,与常规的单次分割相比,基于凸包分析的二次分割能将颗粒粒度分布的平均绝对测量误差和平均不确定度分别由3.78%和2.27%降至2.23%和1.91%,有效地提高了粒度分布检测的精确性和重复性。本研究为下一步开发应用于工业现场的皮带式煤矸石粒度分布在线检测设备奠定了前期基础。