NSST域融合FREAK及全方向相似度的泡沫崩塌率检测

廖一鹏 张进 陈诗媛 王卫星

(福州大学 物理与信息工程学院，福建 福州 350108)

浮选是工业上从矿石中提炼出目标矿物的选矿方法，生产过程中不断的向浮选槽底部输送充足的空气，矿物颗粒和各种化学药剂在搅拌机的旋转搅拌下，可浮性高的矿物颗粒随着气泡上浮到矿浆表面并形成泡沫层，从而将目标矿物与脉石分离。浮选生产工艺受多种物理化学因数的影响，研究人员发现泡沫表面的动态特征与浮选生产工艺指标密切相关[1]，气泡的崩塌程度可有效反映矿物含量[2]。目前，浮选厂主要通过肉眼观察气泡破碎程度进行相应生产变量调节，这种方式没有客观标准，机器视觉技术已在浮选自动化生产中得到成功应用[3-4]，可将其应用于气泡崩塌率的在线检测。

浮选气泡的稳定性影响了回收率、精矿品位等生产指标性能[5-6]，精确地提取气泡崩塌率对于浮选生产指标预测模型的建立至关重要，但是气泡时刻发生着位移、形变、破裂等变化，而且采集的泡沫图像无背景参照，导致崩塌率在线检测困难。目前，国内外主要采用气泡跟踪和特征点匹配两种检测方法：Morar等[7]通过分水岭分割提取相邻图像帧的各个气泡，将后一帧图像的气泡中心点映射到前一帧图像的分割图，计算分割气泡包含的中心点数进行崩塌判断，该方法受气泡运动形变的影响，小气泡的中心点易产生偏移，且准确度受限于分水岭分割结果；Jahedsaravani等[8]先对相邻图像帧进行位移校正再求差，将差值图像中的高亮区域判为相应的崩塌气泡，检测精度受气泡流动和形变的影响，位移较大或形变气泡的亮点和高亮边缘易被误检。文献[9-10]中将SIFT(Scale-invariant Feature Transform)算法应用到气泡速度和崩塌率的检测中，对相邻两帧泡沫图像进行SIFT特征点检测及匹配，根据匹配结果估算崩塌率。该方法的检测精度受限于匹配阈值的选取，存在一定误差，且SIFT算法实时性不强。

鉴于上述分析，可先通过相邻两帧图像的特征点匹配大致估算崩塌位置和范围，然后在崩塌范围内对各个气泡进一步跟踪，判断是否产生崩塌。但是SIFT匹配算法效率低，而现有气泡跟踪方法受分割精度、运动变化的影响大，有待进一步完善解决。2012年，Alahi等[11]在CVPR大会上提出了快速视网膜关键点FREAK(Fast Retina Keypoint)算法，采用类似人体视网膜的二进制描述子。该算法定位精度高、计算快的优点使其在实时性要求高的场合得到广泛应用[12-13]。FREAK可解决SIFT 在实时性方面的不足，但是采用加速分割算法检测特征点，鲁棒性有待进一步提高。近几年，多尺度几何分析的发展为浮选泡沫图像处理开辟了新的研究方向[14-15]，非下采样Shearlet变换(NSST)[16]具有多尺度多方向、平移不变性的变换特性，而且具备运算效率高、分解方向不受限制等特质。可将NSST引入到崩塌气泡的检测中，在高频子带估算崩塌位置和范围，然后在低频子带的崩塌范围内进一步定位崩塌气泡，以提高算法整体的鲁棒性。

综上所述，文中提出一种在NSST域融合改进FREAK匹配及全方向相似度的崩塌率检测方法。相邻两帧泡沫图像进行NSST分解，通过Otsu算法分割低频子带的气泡亮点；将多尺度高频子带分为多个内层和外层，在各个内层通过方向模极大值检测提取边缘区域的兴趣点；接着在本层和上下层通过非极大值抑制提取特征点，改进FREAK采样模型进行特征点描述及匹配；根据前一帧分割亮点周围的匹配点分布密度提取潜在崩塌气泡，对各潜在崩塌气泡及后一帧相应范围的分割亮点采用形状特征码进行表征；然后通过全方向相似度计算进一步确定崩塌气泡，最后根据崩塌气泡的提取结果计算崩塌率。

1 泡沫图像帧及NSST多尺度变换

为在线检测气泡崩塌率，在浮选槽上方安装高清工业CCD摄像机和光源，气泡在刮板作用下朝溢流口方向移动，采集的气泡图像帧如图1所示。前两帧图像较为稳定，第3帧图像的气泡产生了崩塌。前两帧图像中的气泡有一定的位移，每个气泡因反射外部光源而在气泡顶部区域出现高亮点，亮点形状与气泡相似，其大小与气泡成正比。当气泡产生崩塌时，相应的亮点消失，产生新的气泡和亮点，而周围气泡受影响也产生崩塌或形变，前后帧图像在崩塌区域的信息产生急剧变化。

非下采样Shearlet变换[16]是对Shearlet变换的优化改进，通过非下采样金字塔滤波器组(NSLP)去掉了下采样环节，使图像具有平移不变性、抑制伪吉布斯效应，同时通过剪切滤波器(SF)实现对高频子带的多方向分解。Easley等[17]在合成小波理论的基础上通过对基本小波函数进行扩张、剪切以及平移变换，并结合仿射系统理论构建剪切波基函数，其对二维图像处理的表达式为：

(1)

式中：x表示输入二维图像的像素；ψ表示合成仿射系统中的元素，ψ∈L2(R)2；MHY(ψ)是一个合成小波集的剪切波；k、l和v分别为多尺度变换的尺度、方向和平移的系数；H为控制多尺度分解的各向异性膨胀矩阵；Y为控制多方向分解的剪切矩阵。泡沫图像分解流程如图1所示，采用NSLP进行多尺度分解，图像k次分解后，得到1个低频图像和k尺度高频子带，然后采用SF对各高频子带进行l级方向分解，得到大小与原图像相同的2l个方向子带。分解后的低频子带图像保留了气泡的轮廓信息，边缘不明显，去除了大量噪声，气泡亮点对比度高，有利于亮点分割；分解后的高频子带包含了气泡的纹理、边缘细节和噪声点，适合于匹配特征点的检测。

图1 泡沫图像帧及多尺度变换

2 多尺度域气泡崩塌率检测

2.1 高频子带改进FREAK检测及过程描述

P(i，j)=

(2)

(∀α∈[-1，1]，∀β∈[-1，1])

FREAK算法是基于人类视网膜成像原理而提出的特征描述方法，借鉴人眼的中央视觉细胞与周缘视觉细胞的关系，构建了快速视网膜关键点采样模型，该模型由7层同心圆环构成，每层都带有6个采样点，离中心附近的采样点密集，有利于提取更多的细节信息，而外围的采样点稀疏，便于捕捉更多的轮廓信息。为提高在线检测的实时性，对FREAK采样模型进行简化。如图2所示，将原有7层同心圆环的采样模型简化为4层，减少采样点数，提高效率，也能减少误匹配点数。

图2 简化的FREAK采样模型

FREAK描述符是通过采样点对的强度比较级联而成的一个二进制数据串，假设某采样点的描述符为F：

(3)

(4)

HF(F1，F2)=F1⊕F2

(5)

对于任何特征点p，通过式(5)计算p与各待匹配点的汉明距离HF，找到汉明距离最近的待匹配点p1和次近的待匹配点p2，并将p1、p2与p的汉明距离分别记为HF1和HF2，如果满足HF1/HF2<0.6，则接受p1为匹配点。

2.2 低频图像亮点分割及全方向相似度计算

浮选气泡图像中，每个气泡顶部出现高亮区域，亮点形状与气泡相似，其大小与气泡成正比，可对气泡亮点进行分割，将亮点分割代替气泡分割。浮选气泡原图受噪声、高亮边缘带的干扰，分割结果存在大量噪声孤立点和线形边缘带，个别气泡亮点因边缘带而发生粘合，不利于后续处理。NSST分解后的低频子带图像去除了大量噪声和边缘细节，气泡亮点对比度高，分割效果好。如果气泡产生崩塌，前后两帧图像在崩塌区域的亮点个数和形状产生变化，可通过前后帧气泡的形状相似度进行崩塌判断。

形状是高级别的视觉信息，具有较强的视觉表征性，能从语义上表达检测目标的内容。带符号三角形面积形状描述法[19]具有平移、旋转、抗扭曲和缩放不变性，同时包含局部与全局特征，适合于对动态运动的气泡亮点的形状描述。气泡亮点区域较小，形状都类似椭圆，为进一步提高亮点的识别精度，在原来形状描述法的基础上，文中采用一组平行于主方向的等距线将亮点区域分割成若干份，如图3所示。

图3 形状分割示意图

计算每个分割份的形状复杂度，得到一组复杂度特征，沿主方向旋转分割线于不同的角度，得到不同方向的复杂度特征，使得复杂度特征带有方向性，具有更高的抗形变能力和识别精度。

亮点形状相似性检测时，待检测的亮点可能发生过旋转，导致各分割线的方向不一致，为了确保相同的分割方向，将惯性主轴方向作为形状的主方向[20]，该方向不随形状的转换而改变。从主方向开始，每次间隔特定角度θ对形状进行分割，N-1条等距平行线将形状分割为N等份，假设分割形状包含H个轮廓点，任意3个连续轮廓点(xh-t，yh-t)、(xh，yh)、(xh+t，yh+t)组成的带符号三角形面积ST为

(6)

式中，h∈[1，H]、t∈[1，(H-1)/2]，当t从 1 增加到[(H-1)/2]时，可得到不同尺度的ST，通过所有轮廓点在各尺度下的最大ST和最小ST计算各个分割份的形状复杂度CS：

(7)

根据式(7)计算第m个方向第n个分割形状的复杂度CS(m，n)，然后计算各方向下的复杂度特征DC(m)及全方向复杂度特征DC，计算如下：

(8)

(9)

根据各个方向上的复杂度特征差异程度DV(m)计算全方向复杂度特征差异程度DV：

(10)

最后，两个待检测形状的全方向相似度ρ为

ρ=1-DV

(11)

2.3 算法实现流程及步骤

综上分析，文中浮选泡沫崩塌率检测算法流程如图4所示，具体实现步骤如下。

步骤1摄像机实时采集浮选槽表面的泡沫图像，并将图像传输给后台处理主机。

步骤2抽取时间间隔为Δt的连续两帧图像It和It+1，对It和It+1进行NSST多尺度分解，分别得到1个低通子带图像和k个尺度高频子带，各尺度高频子带再分解为l个方向子带。

步骤3把It和It+1的k个尺度高频子带分为多个内层和外层，在各个内层的边缘区域提取兴趣点，然后对各个兴趣点通过式(2)在本层和上下层通过非极大值抑制提取特征点。

步骤4对It和It+1的各特征点采用改进的FREAK采样模型进行描述，根据式(3)-(5)计算匹配点间的汉明距离HF。对于It的任何特征点p，找到汉明距离最近的待匹配点p1和次近的待匹配点p2，并将p1、p2与p的汉明距离分别记为HF1和HF2，如果满足HF1/HF2<0.6，则接受p1为匹配点。将匹配结果中一定长度和斜率范围的点作为RANSAC算法[21]的内点，进一步剔除误匹配点，并计算匹配点间的平均位移d。

步骤5对It和It+1的低频子带图像进行Otsu分割，去除分割结果中的线状干扰区域。统计It中各分割亮点的大小Sq和总个数Q，根据式(12)计算各个气泡的半径Rq，其中SI为图像的大小。

(12)

步骤6气泡崩塌可能带动周围气泡产生崩塌或形变，使得崩塌气泡以及外围气泡区域的匹配点数急剧下降。因此，对It中的各分割亮点，统计以分割亮点的质心为圆心、以2Rq为半径的圆形区域内的匹配点数，如果匹配点数为零，则该气泡为潜在崩塌气泡。

步骤7对It中的各个潜在崩塌气泡，在以相应亮点的质心为圆心和以2(Rq+d)为半径的圆形区域内搜索It+1中的所有分割亮点，根据式(8)-(11)计算潜在崩塌气泡亮点与It+1中圆形区域内所有分割亮点的全方向相似度ρ，如果ρ都小于阈值Tρ，则该气泡为崩塌气泡。

步骤8最后根据提取的崩塌气泡数量B及相应亮点大小Sb，计算崩塌率：

(13)

图4 气泡崩塌率检测流程图

3 实验结果与分析

为验证文中方法的有效性，以福建金东矿业股份有限公司的铅矿浮选槽作为实验对象，实验的硬件平台为Intel(R)Core(TM)i5-4570 CPU@3.20 GH 、4.00 GB(RAM)，软件运行环境为Windows 7 Matlab 2014a。通过大量实验对提出方法进行验证，给出了各实验步骤的结果及分析，并与现有方法进行了结果比较分析。

在NSST域的多尺度高频子带提取特征点，采用改进的FREAK采样模型进行特征点描述及匹配，借助RANSAC算法进一步剔除误匹配点，实验结果及比较如图5所示。抽取连续两帧分辨率为256×256的气泡图像It和It+1进行NSST分解，得到1个低频子带图像和5个高频子带图像，把It和It+1的2、4尺度高频子带定为内层，其他高频子带定为外层。首先在两个内层通过模极大值检测提取兴趣点，然后在本层和上下层通过非极大值抑制提取的特征点如图5(b)、(c)、(f)和(g)所示，融合后的特征点如图5(d)和(h)所示，特征点主要集中在边缘区域。对各特征点采用改进的FREAK采样模型进行描述及匹配，RANSAC算法剔除误匹配点后的结果如图5(l)所示。SIFT、BRISK、FREAK算法的匹配结果如图5(i)、(j)和(k)所示，匹配点对较多，但都存在误匹配点。为客观评价算法性能，分别采用4种算法对相同的300对图像进行实验，各项指标平均值统计如表1所示：SIFT算法的匹配效率最高，而且正确匹配率达到98.73%，但是运算效率很低；BRISK算法运算效率有所提高，但是匹配效率降低，而且出现大量误匹配点对，正确匹配率为89.09%；FREAK算法的运算效率大大提高，但是匹配效率和正确匹配率有待进一步提高；文中方法检测的特征点较少，但是匹配效率高，去除了所有误匹配点，而且运算效率满足在线检测的要求，整体性能最好。

为评估本文形状全方向相似度检测方法的性能，使用MPEG-7 CE-1 Part B图库和气泡亮点分割图像进行测试，旋转角θ=60°，分割份N=5，部分实验结果如图6所示，第1列为待检测的原始形状，前两个为MPEG-7 CE-1 Part B图库的形状，第3个为分割的气泡亮点形状。第2列是对原始形状缩放后的相似形状，计算得到的相似度ρ都在95%以上，具有尺度不变性；第3列是对原始形状旋转一定角度后的相似形状，计算得到的相似度ρ都在93%以上，具有旋转不变性；第4列是产生细小形变的相似形状，计算得到的相似度ρ保持在85%～90%之间，具有抗细小形变的能力；第5列和第6列是有较大差异的形状，随着差异程度的加大，相似度急剧下降，对差异形状具有较强的识别能力。浮选气泡亮点受光照或运动变化的影响，前后帧图像中同一气泡的亮点形状可能产生缩放、旋转、细小形变等，当气泡崩塌后，相应的气泡亮点消失，产生多个新的小气泡及亮点，形状上有一定的差异。因此，文中提出的形状相似度检测方法适用于崩塌气泡的判断，通过大量实验验证，形状相似度的判断阈值Tρ取0.860 0可取得最佳效果。

(a)图像It

(b) It高频尺度2特征检测

(c) It高频尺度4特征检测

(d) It特征点检测结果

(e)图像It+1

(f) It+1 高频尺度2特征检测

(g) It+1 高频尺度4特征检测

(h) It+1 特征点检测结果

(i) SIFT算法结果

(j) BRISK算法结果

(k)FREAK算法结果

(l)改进FREAK算法结果

表1 特征点检测与匹配结果比较

Table 1 Comparison of feature points detection and matching results

方法特征点匹配点匹配效率/%误匹配点正确匹配比/%时间/sSIFT(723,758)31643.71498.732.2354BRISK(1032,1056)38537.314289.090.1207FREAK(985,997)36136.644886.700.0325本文 (696,754)29141.810100.000.2138

为验证本文崩塌率检测方法的精度，选取产生崩塌的前后两帧256×256图像It和It+1进行实验，并与现有文献方法进行结果比较分析，实验结果如图7所示。提取的前后两帧图像It和It+1如图7(a)所示。It+1中有部分气泡产生了崩塌，在NSST域的高频子带采用改进的FREAK算法提取的特征点及匹配结果如图7(b)所示。崩塌区域的

图6 形状相似度检测结果及比较

图像信息发生巨大变化，使得崩塌气泡及周围不存在匹配点，根据匹配结果计算两帧图像的平均位移d=4.12。直接采用Otsu算法分割的气泡亮点如图7(c)所示，存在大量噪声孤立点和线形边缘带，个别分割亮点区域发生粘合，不利于后续处理。分解后的低频子带图像如图7(d)所示，去除了大量噪声和边缘细节，气泡亮点对比度高。分割后的气泡亮点如图7(e)所示，有效去除了噪声孤立点和边缘带，亮点区域完全分离，根据分割结果估算各个气泡的半径Rq。将高频子带检测的匹配点对映射到低频子带的亮点分割图像上，如图7(f)所示。对It中的各分割亮点，搜索以亮点质心为圆心、以2Rq为半径的圆形区域内的匹配点数，如果无匹配点，则该气泡为潜在崩塌气泡，检测结果如图7(g)的黄色区域。对It中的各个潜在崩塌气泡，以2(Rq+d)为半径的圆形区域内搜索It+1中的分割亮点，并计算两者间的全方向相似度ρ，如果ρ都小于阈值Tρ，则确定该气泡为崩塌气泡。崩塌气泡提取结果如图7(g)的绿色区域，与实际崩塌气泡较为吻合。采用人工选取崩塌气泡计算得到的崩塌率Pb=26.53%，文中方法根据崩塌气泡数量B及相应亮点大小Sb，通过式(11)计算得到崩塌率Pb=29.39%，误差较小。

(a)图像It和It+1

(b)It和It+1的特征点及匹配结果

(c)It和It+1的Otsu分割结果

(d)It和It+1的低频子带图像

(e)It和It+1的低频子带分割结果

(f)It和It+1的气泡亮点与匹配点

潜在崩塌气泡 崩塌气泡检测结果

气泡中心点检测 崩塌气泡检测结果

差值图像 崩塌气泡检测结果

It的匹配结果 It+1的匹配结果

文献[7]中对It和It+1进行分水岭分割，将It+1各气泡中心点映射到It分割图，通过It中气泡包含的中心点个数进行崩塌判断，结果如图7(h)所示。受崩塌和运动形变的影响，It+1中小气泡中心点易被挤入It中的大气泡而产生误判，且受光照影响，大气泡区域的过分割、小气泡区域的欠分割也会降低检测精度，检测的Pb=32.82%，存在一定误差。文献[8]中先对It和It+1进行位移校正再求差，将差值图像中的高亮区域判为崩塌气泡，结果如图7(i)所示。由于受整体运动不均匀性的影响，It和It+1的大气泡亮点区域、高亮边缘区域易发生偏移，使得差值图像中的高亮区域增多，检测的Pb=35.52%，误判的崩塌气泡较多。文献[9-10]中通过It和It+1的SIFT匹配结果估算崩塌率，结果如图7(j)所示。当崩塌判断阈值取0.12，计算的Pb=40.83%。该方法受限于判断阈值的选取，误差较大，实用性不强。综上分析，文中方法受运动形变、运动不均匀的影响小，能准确提取各个崩塌气泡，检测精度高，实用性强。

将文中方法进行工业现场试验，以福建金东矿业股份有限公司梅仙浮选厂的铅矿浮选槽作为试验对象，采集2019年3月11日至2019年3月29日期间的正常浮选、过浮选和欠浮选3种工况下的泡沫视频图像进行实验。各工况下选取1 800对图像进行测试，以专家人工标注方法计算的崩塌率为标准，采用崩塌率Pb、误差绝对值Eb和运行时间Tb对实验结果进行评估，并与现有文献方法进行比较分析，3项指标的平均值统计结果如表2所示。实验数据表明：正常浮选下泡沫流动速度中等且方向稳定，文献[7]和文献[8]方法受不均匀流动和运动形变影响，文献[9-10]方法受限于判断阈值的选取，存在一定的误差，而文中方法受影响小，检测精度较高；过浮选下泡沫流动缓慢而造成堆积，相邻两帧图像的信息变化较小，文献[7]和文献[8]方法受运动影响小，而文献[9-10]方法的匹配点数上升，检测精度都有一定提高；欠浮选下泡沫快速运动且变化大，运动方向不稳定，导致文献[7]、文献[8]和文献[9-10]方法的检测精度急剧下降，而文中方法仍保持较高的精度。在运算效率方面：文献[7]和文献[8]方法的运算效率较高，且受工况变化的影响小，各工况下的运算时间波动较小；文献[9-10]方法和文中方法采用基于特征点匹配，在各工况下的匹配点数量急剧变化，导致各工况下的运行时间波动较大。文献[9-10]方法的实时性较差，文中方法提高了运算效率，平均运行时间为0.403 1 s，满足在线检测的需求，且各工况下保持较高的检测精度。综上分析，本文方法检测的崩塌率能定量表征各工况下浮选气泡的稳定性，且在不同浮选工况下均表现出良好的鲁棒性，满足浮选生产在线检测和动态变化需求。

表2 气泡崩塌率检测结果客观评价

4 结论

针对浮选表面气泡受不均匀流动、运动形变影响而导致崩塌率在线检测困难的问题，提出了一种在NSST域结合改进快速视网膜关键点匹配及全方向相似度计算的气泡崩塌率检测方法。在多尺度高频子带借鉴BRISK算法思想提取边缘区域的特征点，改进FREAK采样模型并用于特征点描述及匹配，改善了FREAK算法的整体鲁棒性；基于形状全方向复杂度特征计算的全方向相似度检测具有抗缩放、旋转、细小形变的特点，而对差异形状具有较强的识别能力；结合前一帧亮点周围匹配点分布密度与前后两帧亮点的全方向相似度计算的方法能有效提取崩塌气泡。实验结果表明：改进的FREAK算法较SIFT、BRISK、FREAK算法的匹配效果好，崩塌率检测方法受泡沫不均匀流动、运动形变的影响小，能有效提取出各个崩塌气泡，检测精度较现有方法有较大提高，能定量表征各生产工况下气泡的稳定性，且在不同工况下均表现出良好的鲁棒性，满足浮选生产在线检测和动态变化需求。