廖一鹏,陈诗媛,杨洁洁,王志刚,王卫星
(1. 福州大学 物理与信息工程学院,福建 福州350108;2. 福建金东矿业股份有限公司,福建 三明365101)
浮选生产是工业上从矿石中提炼出金属矿物的一种有效选矿方法,通过加入特定的化学药剂将金属从矿石中分离,药剂添加的好坏直接影响最终生产指标的高低。目前,浮选厂主要通过人眼观察泡沫表面状态的变化进行相应的药剂量调节,这种方式主观性强且实时性差。 近年来,研究人员通过机器视觉技术提取泡沫表面特征来进行加药状态预测,取得了一定成果:Zhang等[1]通过分水岭算法提取泡沫尺寸分布特征,并建立非线性预测模型实现加药状态的有效识别;Ai 等[2]提出一种基于泡沫形状加权大小分布特征的加药状态预测方法,进一步提高了预测准确度;Tang 等[3]提出一种基于泡沫大小累积分布特征的加药状态预测方法,减少了噪声干扰,提高了稳定性;Li 等[4]提出一种基于泡沫图像深度学习的加药故障状态检测方法,在大数据集条件下具有较高的识别精度。以上方法主要通过提取泡沫尺寸、形态特征进行加药状态预测,Cao 等[5]通过试验发现泡沫表面的流动特征与药剂添加量之间有较强的相关性,在不同加药状态下呈现不同的流动特性,增加泡沫流动特征作为加药状态识别的驱动,将有助于提高识别精度。
目前,国内外学者主要研究泡沫静态特征的提取,而对流动特征提取方法的研究较少,浮选过程中泡沫时刻发生着崩塌、合并、新产生等变化,而且浮选槽表面雾气大、粉尘多,采集的图像受光照和噪声影响大,流动特征检测困难。陈良琴等[6]结合气泡亮点跟踪和相位相关法估计气泡的流动速度,该方法采用Otsu 分割气泡亮点,分割精度受气泡形变和光照的影响,而且相位相关法计算获取的是两帧图像的相对位移,不能精确描述泡沫的局部运动特征;WANG 等[7]将SIFT(Scale Invariant Feature Transform)算法应用到泡沫流动特征的检测中,根据匹配结果计算流动特征,提取的流动特征在各加药状态下有一定的区分度,但SIFT 计算复杂,实时性不强;Nakhaei等[8-9]采用像素点跟踪技术计算连续帧间的平均流动速度,该方法运行效率高,但是像素值在流动过程容易受噪声和光照影响而产生变化,导致像素点跟踪出错。
Roblee 等[10]提 出 的ORB(Oriented Fast and Rotated Brief)算法具有旋转不变性和运算效率高的优点,适用于实时性要求高的应用场合,解决了SIFT 在实时性方面的不足,但该算法提取的特征点容易出现簇集现象,而且抗噪声和尺度变换的能力差,有待改善。近几年,多尺度几何分析的发展为图像处理领域开辟了新的研究方向 ,非 下 采 样Shearlet 变 换[11](NSST,Nonsubsampled Shearlet Transform)具有多尺度多方向特性,分解后的图像具有平移不变性,而且具备运算效率高、分解方向不受限制等特质,可将NSST 引入到ORB 算法的改进,提高算法整体的鲁棒性。
鉴于上述分析,本文提出一种基于NSST 域改进ORB 算法的泡沫流动特征提取方法,并用于加药状态识别。在NSST 域先采用尺度相关性去除噪声,在各个内层通过方向模极大值检测提取边缘区域的兴趣点,然后在本层和上下层通过非极大值抑制进一步提取特征点,解决特征点的簇集现象和尺度不变性。采用多尺度BRIEF(Binary Robust Independent Elementary Features)描述子对特征点进行描述,结合泡沫流动趋势动态调整搜索的匹配区域,以提高匹配效率和精度,根据匹配结果计算泡沫流动特征。 最后,构建行列自编码极限学习机对泡沫形态、尺寸分布特征和流动特征进行融合,然后通过自适应随机森林算法进行加药状态识别。
NSST 是对剪切波变换的优化改进[11],去掉了下采样环节使图像具有平移不变性和抑制伪吉布斯效应,泡沫图像通过k 级分解,得到1个低频图像和k 个不同尺度的高频子带图像,各高频子带图像可再分解为多个方向子带,且所有分解图像大小相同。泡沫图像的3 级NSST 多尺度分解如图1 所示,分解后的低频子带保留泡沫的轮廓信息,而泡沫纹理、边缘细节和噪声点被分解到各高频方向子带,可在多尺度高频子带进行噪声去除和特征点检测。
图1 泡沫图像NSST 多尺度分解Fig. 1 Foam image NSST multiscale decomposition
泡沫图像在光源的作用下,每个气泡表面呈现高亮区,流动过程中受光照影响大,而气泡边缘像素点受光照影响小,可将NSST 域多尺度边缘带作为特征点检测区域,以减少光照影响,同时解决特征点的簇集现象。在进行边缘系数检测前,先进行噪声去除处理,随着分解尺度的加深,噪声系数在尺度间是弱相关的,会迅速衰减,边缘系数强相关而保持稳定[12],可根据尺度间系数相关强弱将噪声点从边缘中分离。定义像素点( i,j ) 在k 尺 度l 方 向 子 带 上 的 系 数 为clk( i,j ),尺 度 相 关 系 数 为Corrkl(i,j),K个 尺 度 间 的 系 数为 第k尺 度 上 第l方 向 子 带的 系 数 能 量 ,是 归 一 化 处 理 ,则 像 素 点(i,j) 的 尺度相关系数为:
如果该点的原系数大于尺度相关系数,则该点为噪声点,将该系数值置零,否则系数值保持不变:
为保证尺度不变性,在NSST 域对ORB 特征点检测方法进行改进,借鉴BRISK(Binary Robust Invariant Scalable Keypoints)[13]特 征 点 检 测方法,特征点检测在尺度空间上进行,把NSST多尺度高频子带分为多个内层和外层,在各个内层的边缘区域提取兴趣点,然后在本层和上下层通过非极大值判断进行特征点检测:计算内层上(i,j) 点在各个方向分量的模值,假设l方向的模值mod [clk(i,j) ] 最大,以该点为中心,计算半径为3 的圆上16 个系数点在l方向的模值,如果16个模值至少有9 个的值小于mod [clk(i,j) ],则该点为潜在的兴趣点。然后该兴趣点在本层和上下层通过非极大值抑制进行特征点判断,同层中该点的l方向的模值mod [clk(i,j) ]必须要大于其余八邻域中的点,上下层对应的2×9 个点的l方向的模值也必须小于mod [clk(i,j) ],特征点P(i,j)的检测公式如式(3)所示:
其中:m∈[-1,1 ],n∈[-1,1 ]。
为使特征点描述具有旋转不变性,必须先规定特征点的主方向,ORB 将特征点与邻域质心的夹角定义为主方向[14],特征点邻域图像矩定义为:
其中:I(i,j)为像素点(i,j)的灰度值,mpq为p+q阶矩。特征点邻域图像质心为:
特征点与质心的夹角为特征点的主方向:
获取特征点主方向后,本文采用多尺度系数进行特征描述,在特征点周围选取M对随机点,将M对随机点旋转σ 角度后,对每对点的低频子带梯度以及高频子带系数进行比较,将比较结果组合起来作为描述子,定义比较准则:
其中:I(a) 和I(b) 为随机点对a和b的低频子带梯度,P(a) 和P(b) 为全方向的K个尺度系数乘积之和。通过M对随机点的比较,构成M维的二进制特征描述子:
每个特征点对应的描述符B均为M位的二进制比特串,匹配时采用异或操作计算对应描述符B1和B2之间的汉明距离HD,用于衡量两个特征点间的近似度,值越大代表相似度越小,反之越大:
浮选槽表面的泡沫受刮板的作用其流动方向和速度在短时间内相对稳定,相对于前一帧图像中的特征点,后一帧图像的匹配点应落于沿泡沫流动方向的一定范围内,可根据这特性缩小匹配的搜索区域,假设vx和vy分别为前两帧图像所有匹配点的水平速度和竖直速度的均值,其合成速度为v。定义搜索区域为以当前特征点p为圆心,以r为半径的圆,搜索区域的半径公式为:
其中:nr为区域扩大系数,初始值为1,在搜索区域中找到与该点汉明距离最近的关键点p1和近次的关键点p2,并将p1,p2与p 的汉明距离分别记为HD1和HD2,如果满足HD1/ HD2<0.5,则接受p1为匹配点,反之依据式(11)循环扩大搜索范围直至找到匹配点或搜索至全图。最后计算匹配点对间线段的长度和斜率,将一定长度和斜率范围的点作为RANSAC[15]内点,剔除误匹配点,提高匹配精度。
假设时间间隔为△t 的连续两帧图像It和It+1共有N 对匹配点,它们在It和It+1的位置分别水平流动速度vnx和垂直流动速度vny,通过式(13)计算平均水平流动速度vx和平均垂直流动速度vy:
通过式(14)计算当前的泡沫平均流动速度和方向,假设前一时刻检测的水平流动速度为vx′,垂直流动速度为vy′,通过式(15)计算当前的水平流动加速度ax和垂直流动加速度ay,通过式(16)计算流动速度和方向的无序度,-v 和方向-θ为前一段时间的平均流动速度和方向。
假设摄像机镜头到目标物体的物距为S,镜头焦距为f,摄像机的光学放大倍率为μ,焦距和物距的关系为:
式中:vT为泡沫表面流动速度实际值,由此可推出泡沫流动速度与摄像机高度及镜头焦距的关系式为:
根据精矿品位和尾矿品位的等级,将浮选加药状态分为正常、过量、欠量、故障等4 种状态[4]。泡沫表面的形态、尺寸分布特征是加药状态评判的重要依据,研究发现泡沫表面流动特征与药剂添加状态也有较强的相关性[5],本文综合考虑两类特征作为加药状态评判的依据。黄凌霄等[16]在多尺度域提取泡沫的亮点个数、平均面积、标准差、椭圆率和分形维数作为等效形态、尺寸特征,提取的特征与泡沫类别的相关性较高,本文采用这5 个特征作为泡沫表面的形态、尺寸分布特征,结合提取的6 个流动特征作为加药状态评判的依据。
随机森林(Random Forest,RF)[17]是基于多决策树的集成学习算法,对异常值不敏感,抗噪声能力强,具有很好的分类性能。为了提升决策树的多样性,降低树之间的相似性,避免产生过拟合现象,先对形态、尺寸分布特征和流动特征进行融合和维度扩展。自编码极限学习机[18]令输入与输出相等,通过前馈神经网络来完成高层次的特征提取,具有较高的特征抽取和学习能力:
式中:β 是隐含层和输出层间的连接权值矩阵,X表示与输入相等的输出值,H 表示隐层单元的输出矩阵。随机生成输入与隐层间的权重和偏置可得隐层输出矩阵H,通过式(20)求得输出权重矩阵β:
经过训练后,隐层节点H 的输入是输出权重的转置βT,通过式(21)对原始特征进行高维特征抽取:
为实现特征融合和高维扩展,如图2 所示,先通过两个自编码极限学习机分别对形态、尺寸分布特征和流动特征进行学习和抽取,得到N1维特征向量H1和N2维特征向量H2,然后通过式(22)将向量H1和向量H2的转置进行叉乘,将两类特征融合在一起,得到N3(N3= N1× N2)维特征向量E,有助于扩大不同决策树间样本采样和特征选取的差异程度,提升决策准确率。
随机森林采用有放回抽样的方式随机从训练样本中抽取Nt个样本子集建立Nt棵决策树,在决策树的创建过程中,从所有特征中随机选取比例为Pm的特征变量创建决策树,Nt和Pm在模型中都是超参数,不同组合对最终的分类效果产生重要影响。本文以测试样本集的识别精度作为适应度,采用量子菌群算法[19]对两个超参数进行优化。浮选加药状态识别的具体步骤如下:
Step 1:采集不同加药状态下浮选槽表面的泡沫图像,并将视频图像分类存储。
Step 2:对不同加药状态下的泡沫视频图像,在NSST 域提取泡沫的6 个等效形态、尺寸分布特征,以及5 个流动特征。
Step 3:对不同加药状态下视频图像提取的等效形态、尺寸分布特征和流动特征,分别通过自编码极限学习机进行学习和抽取,然后将两组特征行列矩阵相乘进行特征融合,创建样本数据集。
Step 4:通过样本数据集训练随机森林多分类器,并采用量子菌群算法[19]对Nt和Pm两个超参数进行自适应优化。
Step 5:实时采集浮选槽表面的泡沫图像,提取泡沫的形态、尺寸分布特征和流动特征,然后通过行列自编码极限学习机对特征进行抽取和融合,最后采用训练的自适应随机森林进行决策分类。
为验证本文方法的性能,以福建金东矿业股份有限公司的铅矿浮选槽作为实验对象,在精选槽泡沫表面溢流口上方80 cm 处安装工业CCD摄像机及光源,实验硬件平台为Intel(R)Core(TM)i7-9700 CPU@3 GHz ,8 GB(RAM),软件运行环境为Windows 7 Matlab 2014a。
在NSST 域进行ORB 算法改进实验,并与SIFT 算法、ORB 算法的实验结果比较,结果如图3 所示。抽取间隔△t=0. 5 s 的连续两帧256×256 图像It和It+1,NSST 分解得到1 个低频图像和5 个尺度高频子带,把第2,4 尺度高频子带定为内层,而第1,3,5 尺度高频子带定为外层。两个内层提取的特征点如图3(b),(c),(f)和(g)所示,融合后得到特征点检测结果如图3(d)和(h)所示,大部分特征点都均匀分布在泡沫边缘区域,有效解决了簇集现象。对各特征点采用多尺度BRIEF 描述子进行特征描述,并结合前两帧图像检测的流动速度(33. 99 pixel/s)和方向(21. 54°)动态调整搜索的匹配区域,匹配结果如图3(i)所示,因循环扩大搜索的匹配范围而出现了误匹配点,通过RANSAC 算法进一步剔除了所有误匹配点。传统SIFT 算法的检测结果如图3(k)所示,匹配点对较多且分布均匀,但出现了误匹配点。传统ORB 算法的检测结果如图3(l)所示,匹配点对较少且簇集,也出现了少量误匹配点。
为验证本文改进ORB 算法的抗噪声和尺度变换性能,分别对100 对叠加均值为0、不同方差的高斯白噪声图像和100 对不同尺度比例的图像进行匹配实验,并与近年所提出的改进SIFT[20]和ORB[14]算法进行结果比较分析,平均匹配精度(正确匹配点对数与总匹配点对数的比率)和运行时间如表1 所示,在噪声方差为10%、尺度比为1∶2 情况下:改进的SIFT 较原SIFT 的运算时间较少了一半,抗噪声和尺度变换性能较好,平均匹配精度较高,但是运算效率有待进一步提高;改进的ORB 算法的抗噪声和尺度变换性大大提升,而且保持较高的运算效率;本文改进ORB 算法的抗噪声和尺度变换能力大大提升,抗噪声性能优于改进的SIFT 和ORB 算法,抗尺度变换性能与改进SIFT 相当,但是运算效率高于改进的SIFT 算法。当噪声方差增加到30%、尺度比增大到1∶8 情况下:3 种算法的匹配性能都大幅度下降,但是本文改进ORB 算法保持高于85% 的匹配精度,具有较优的抗噪声和尺度变换能力。
图3 改进ORB 的实验结果及比较Fig. 3 Improved ORB algorithm experimental results and comparison
表1 匹配精度及运行时间比较Tab. 1 Comparison of matching accuracy and running time
为验证本文流动特征检测的精度,从采集的图像库中,人工截取实际位移为(20,20)的两帧256×256 图像It和It+1进行仿真实验,计算速度大 小v、方 向θ、相 对 误 差Ev=|v-v0|/v0×100% 和Eθ=|θ-θ0|/θ0× 100%(v0和θ0为 实际值),以及运行时间,并与现有文献方法进行比较分析,实验结果如图4 所示,结果数据统计如表2 所示:文献[6]方法的实验结果如图4(f),对大津阈值分割的气泡图像进行相位相关计算,运算效率高,但获取的是两帧图像的整体相对位移,检测精度低;文献[7]方法的匹配结果和速度矢量图如图4(g)所示,速度矢量图上的带箭头线段代表各个特征点的速度矢量,检测范围广且分布均匀,但是SIFT 算法的运算效率低。文献[8-9]方法的像素跟踪结果和速度矢量图如4(h)所示,跟踪错误率低、检测效率高,但是跟踪点较少、分布不均匀。本文方法的实验结果如图4(e)所示,匹配点多且分布较为均匀,Ev=0. 012 0%,Eθ=0. 045 3%,检测精度高,而且运算效率较文献[7]方法有较大提高。
表2 特征检测结果及比较Tab. 2 Feature detection results and comparison
为验证本文方法的抗噪声和光照影响的性能,对图像叠加均值为0、方差为20% 的高斯白噪声,如图5(a)和(b)所示,对It+1通过伽马校正进行亮度非线性调整,如图5(c)所示。采用不同方法对两种图像进行实验,结果如图5(d)~图5(g)所示,数据统计如表2 所示:文献[6]方法的功率谱峰值图如5(d)所示,获取的是两帧图像的整体相对位移,运算效率高,但气泡分割结果受噪声和光照变化影响大,导致检测精度不高;文献[7]方法的速度矢量图如5(e)所示,检测范围广且分布均匀,抗噪性能好,但是运算效率低,光照影响下的检测精度不高;文献[8-9]方法的像素跟踪矢量图如5(f)所示,检测效率高,但是跟踪点较少、分布不均匀,在噪声和光照变化的干扰下跟踪容易出错,导致检测精度低。本文方法的结果如图5(g)所示,匹配点多且分布均匀,在噪声和光照变化影响下匹配点减少但分布均匀,相对误差Ev和Eθ都较小,仍然保持较高的检测精度,抗噪声和光照影响能力强。
图4 特征检测效果及比较Fig. 4 Feature detection effect and comparison
图5 抗噪声和光照影响性能Fig. 5 Noise and light resistance affect performance
为了使相机镜头不受其气泡飞溅矿浆的影响,浮选表面离相机镜头的距离为80 cm,匹配镜头为CMP 百万级工业镜头M0814-MP,调节镜头的焦距位8 mm。 采集铅矿精选槽在正常加药、欠量、过量、故障等4 种状态下的泡沫视频图像,然后进行泡沫形态、尺寸分布特征和流动特征提取,以平行刮板转轴的方向为X轴方向、流溢口的方向为Y轴方向,实验结果如图6 所示:正常加药下泡沫流动速度保持在20~50 pixel/s 之间且变化缓慢,流动方向为75°~100°,大部分朝着流溢口的方向流动,加速度较小且变化不大,流动速度和方向有一定无序度;加药过量下泡沫流动速度处在5~25 pixel/s,流动缓慢且曲线呈现周期性的上下波动,流动方向稳定在90°附近,加速度较小且波动不大,但速度呈现极大的无序性;欠药量下泡沫流动速度在40~60 pixel/s 之间且上下波动较大,流动方向稳定在70°左右,加速度较大且变化幅度大,速度和方向的无序度小;故障状态下泡沫流动速度在50~70 pixel/s之间,流动快且上下波动较大,流动方向在50°~70°之间来回变化,加速度较大且变化幅度大,流动方向呈现较大的无序性。实验数据表明:提取的流动特征与实际加药状态下的泡沫流动趋势吻合,实现了流动特征的定量描述,能准确表征各加药状态下泡沫的流动特性,且在各状态下具有一定的区分度,可作为加药状态预测的有效特征。
图6 不同加药状态下的流动特征提取结果Fig. 6 Extraction results of flow characteristics under different dosing states
采集4 种加药状态下各2 500 组泡沫图像,其中2 000×4 组作为训练集,其他500×4 组作为测试集,提取各组图像的5 个形态、尺寸分布特征和6 个流动特征,然后通过行列自编码极限学习机对这些特征学习和抽取,设置自编码极限学习机的隐含节点与输入节点数相同(N1=6,N2=5),两组特征行列矩阵相乘后得到30 维的特征向量。为验证本文行列自编码极限学习机的性能,分别采用原始的11 个特征和融合后的30 个特征对随机森林进行训练和测试,参数Nt取值范围:100~900,Pm取值范围:0. 1~0. 9,在取值范围内调整两个参数进行实验,两种方法的测试准确率曲面图如图7(a),(b)所示:本文通过行列自编码极限学习机的特征学习和抽取,将两类特征进行融合及维度扩展,有助于扩大不同决策树间样本采样和特征选取的差异程度,避免产生过拟合现象,整体准确率大大提升,分类精度受Nt和Pm的影响而产生的波动小,稳定性好。 为验证本文加药状态识别效果,训练过程中采用量子菌群算法自适应获取随机森林算法的最优参数Nt=503,Pm=0. 6,各加药状态的识别结果如表3 所示,本文方法能有效识别4 种不同的加药状态,平均识别准确率达到97. 85%。
图7 测试准确率曲面图Fig. 7 Surface diagram of test accuracy
表3 加药状态识别结果Tab. 3 Identification results of dosing status
为验证本文方法抗噪声和光照变化影响性能,对测试集图像叠加均值0、方差5% ~40% 的高斯白噪声,采用伽马校正对图像非线性亮度调整模拟光照变化影响,测试结果如表4 所示:当高斯白噪声方差<25% 时,保持较高的识别精度,抗噪声性能好;当伽马校正系数在0. 8~1. 3 之间,保持较高的识别精度,抗光照变化影响能力强;当噪声方差>40% 时,或伽马校正系 数>1. 5、<0. 7 时,泡沫形态、尺寸和流动特征检测误差大而导致识别精度急剧下降。
表4 噪声与光照变化测试结果Tab. 4 Test results of noise and light change
为进一步验证本文方法的优势,采用相同的数据集对现有文献方法进行试验,统计各方法的平均识别精度、标准差以及平均运行时间,实验结果如表5 所示:文献[1]采用分水岭分割方法提取泡沫大小分布特征,建立Hammerstein-Wiene非线性模型对泡沫加药状态进行预测,运算效率高,识别精度达到90% 以上,但是提取的特征单一,且识别精度受分割精度影响大,识别精度有待进一步提高,稳定性差;文献[2]采用形状加权大小分布来表征泡沫特征,通过最小二乘支持向量机对泡沫加药状态进行识别,进一步提高了识别精度,标准差也较小;文献[3]基于时间序列图像提取泡沫大小累积分布特征,采用形状加权大小分布来表征泡沫特征,建立XGBoost 非线性模型实现泡沫加药状态识别,具有较高的识别精度,稳定性好,但是运算效率低;文献[4]构建CNN 网络提取泡沫图像特征,通过训练的SVM进行浮选加药状态判别,大规模数据集训练下具有较高的识别率,但是在等同的小样本数据集下识别精度较低;本文方法的流动特征提取精度高,构建行列自编码极限学习机对泡沫形态、尺寸分布特征和流动特征进行融合,然后通过自适应随机森林对加药状态识别,扩大了决策树间样本采样和特征选取的差异程度,泛化性能好,平均识别率达97. 85%,标准差为1. 01%,识别精度和稳定性较文献[1,2,4]方法有较大提升,与文献[3]方法的识别精度相当,但提高了运算效率。
表5 不同方法的识别效果Table 5 Recognition effect of different methods
本文提出了一种基于NSST 域改进ORB 算法的泡沫流动特征提取及加药状态识别方法。对各高频方向子带采用尺度相关性去除噪声系数,有效提高了算法的抗噪性能。在多尺度域结合方向模极大值检测和非极大值抑制在边缘区域提取特征点,采用多尺度BRIEF 描述子对特征点描述,解决了特征点的簇集现象和尺度不变性,使特征点描述具有旋转不变性,改进的ORB算法的抗噪声和尺度变换性能增强。构建行列自编码极限学习机对泡沫形态、尺寸分布特征和流动特征进行融合,然后通过自适应随机森林对加药状态识别,扩大了决策树间样本采样和特征选取的差异程度,以减少过拟合现象,提升决策准确率。本文方法实现了泡沫流动特征的定量描述,受噪声和光照影响小,检测精度和运算效率较现有方法有较大提高,能准确地表征不同加药状态下泡沫表面的流动特性,并应用于加药状态识别,平均识别精度达97. 85%,较现有文献方法有较大提升,为后续的加药量优化控制奠定基础。