基于数字图像处理的浮选泡沫速度特征提取及分析

唐朝晖，刘金平，桂卫华，阳春华

(中南大学 信息科学与工程学院，湖南 长沙，410083)

摘  要：以浮选泡沫图像序列为对象，研究浮选泡沫图像序列速度特征提取方法，分析泡沫速度特征与浮选性能间的关系。提出一种抗尺度快速变化和具有旋转不变性的模板匹配算法，利用宏块跟踪技术对浮选泡沫图像序列的泡沫速度特征进行估计，得出像素级的泡沫运动速度参量。然后，采用二维拉格朗日曲面插值方法提取亚像素位移，得出精确的亚像素级位移参量。结果表明：在浮选过程中，减少浮选泡沫运动速度的紊乱程度能减小已粘附在泡沫上的矿物粒子的脱附率，进而提高浮选精矿品位，降低尾矿中有用矿物含量。

关键词：浮选；泡沫图像；模板匹配；旋转不变性；亚像素配准

中图分类号：TP391.41          文献标志码：A         文章编号：1672-7207(2009)06-1616-07

Froth bubbles speed characteristic extraction and analysis based on digital image processing

TANG Zhao-hui, LIU Jin-ping, GUI Wei-hua, YANG Chun-hua

(School of Information Science and Engineering, Central South University, Changsha 410083, China)

Abstract: Based on the study of the bubble image sequences, the characteristics of the forth velocities were extracted, and the relationship between the froth speed characteristic and the froth flotation performance was discussed. Hence, a new fast template matching algorithm was presented. This algorithm is invariant to arbitrary rotation revolving and scaling, and using macro-block tracking method to calculate the froth motion velocity, the integral pixels velocity characteristic can be obtained; then, the sub-pixels velocity parameters was extracted using the two-dimension Lagrange surface interpolation and the full precise sub-pixel scale froth bubbles motion velocity was gained. The results show that the decrease of froth velocity turbulence degree in the flotation production is able to decrease the detach ratio of the mineral particles on the bubbles and then to improve the concentrate grade and decrease the ratio of minerals in tailings.

Key words: flotation; froth images; template matching; rational invariant; sub-pixels registration

浮选是一个极其复杂的物理化学过程，难以建立准确的实时控制模型^[1]。在实际的矿物浮选生产中，有经验的操作工人主要通过观察浮选槽表面泡沫视觉信息来完成浮选过程操作。浮选泡沫表面视觉特征参数一般表现在2个方面：一是基于单帧图像的泡沫参数，称为静态参数，比如泡沫的大小、形状、颜色特征、纹理特征^[2]等；另一方面是基于图像序列的泡沫图像特征参数，称为动态特征，一般包括泡沫运动速度、泡沫稳定度等。分析泡沫图像的静态和动态特征，发现它们与浮选工艺参数密切相关^[3]，如泡沫颜色特征反映所携带的矿物的种类和品位；泡沫的大小和形状反映泡沫矿化程度；而浮选泡沫的速度特征与浮选产量、浮选的精矿品位等密切相关。计算机技术的快速发展促进了数字图像处理技术在浮选控制中的应用，并取得较好的浮选效果^[4-6]。因而，研究浮选泡沫图像处理、开发融合了数字图像处理的浮选控制技术，提高浮选性能，达到提高精矿品位指标、降低尾矿品位指标及药剂消耗，从而降低生产消耗的目的，具有较高的市场价值和企业应用价值。国内外许多学者研究了重金属浮选泡沫图像静态特征的提取方法^[7-9]，但对于其动态特征的提取方法并没有详细说明和分析。要建立准确的矿物浮选视觉控制模型，了解浮选泡沫的动态特征和静态特征很重要。动态特征中的稳定性提取是在泡沫速度特征提取的基础上进行的^[5]，所以，动态特征的提取主要是对泡沫图像序列进行速度估计。为此，本文作者总结目前常用的运动估计算法，提出一种快速的可以用于在线检测的浮选泡沫图像运动特征提取方法，它采用具有旋转不变性和抗尺度变化的模板匹配算法，通过实时跟踪泡沫图像序列中的泡沫块来计算泡沫的运动速度。该方法具有亚像素水平的估算精度，并具有抗浮选泡沫任意角度旋转和尺度变化的特点。最后，分析浮选泡沫速度特征和浮选主要性能参数即尾矿矿物含量和总精矿品位间的关系，以便为矿物浮选的高效、稳定生产提供指导。

1  宏块跟踪

目前常用的基于数字图像处理的运动值估计算法主要有基于空间域的块运动分析方法和基于频域的相位相关法。基于块运动分析的方法有光流法^[10]和宏块跟踪法^[11]。光流法虽然能得出比宏块跟踪法更精细的结果，但由于其对光照条件比较敏感，且其计算复杂性高，算法开销大，不适合实时系统的特征提取。基于频率域的相位相关法^[12]对光照噪声有抗干扰性，但其必须进行傅里叶变换，所以，其在原理、计算上较复杂。宏块跟踪法具有算法简单的特点。为了满足现场实时处理要求，针对运动估算算法的特点，本文采用基于宏块跟踪的方法对泡沫图像序列运动特征进行估计。宏块跟踪法是指从图像序列中相邻帧中前一帧的某一位置选出一子块(称为宏块)作为模板，再在当前帧中搜索其最佳的匹配位置，利用其在前后2帧中的位置之差来进行运动参量估计。把前一帧分为M×N块，分别搜索其在当前帧的最佳匹配位置就能估计出全帧图像各子块的运动速度。宏块跟踪中最重要的2个要素是模板匹配原则和搜索策略。

1.1  模板匹配准则

现在常用的匹配准则有最小绝对距离差(SAD)，最小均方误差统计(MSE)、互相关函数(NCCF)等^[13]。由于归一化互相函数具有抗白噪声能力强优点，所以，泡沫图像的块跟踪可以选择该准则，其计算方法如下：

对式(1)可以进行相应的简化以减少计算代价，但只化简式(1)还不能从根本上解决宏块跟踪的计算代价，因为宏块跟踪需要用模板在当前帧设定的搜索区域内进行搜索，如何减少搜索次数是提高计算速度的关键。

1.2  搜索策略

要在当前帧中找到从前一帧中的选定宏块的最佳匹配位置，必须进行相应搜索，其最简单的方法当然是全搜索(FS)，即在图像搜索区域的每个位置进行匹配，但这样会耗费大量时间在不相关的位置进行相关计算，不适合实时检测。故本文所用的搜索策略是分级三步搜索法^[14]，其原则是采用图像金字塔分层匹配准则，从粗到精，沿着最小误差(最大互相关系数)的方向进行匹配搜索，即先在分辨率小的图像上搜索最佳匹配位置，再在该位置周围的小区域进行高分辨率搜索。

图1所示为分级三步搜索示意图，它表示一个简单的两级搜索匹配。首先，利用图像金字塔原则，将图像模板和待搜索图像都在低分辨率(比如为实际图像分辨率的一半)下搜索。搜索采用三步法进行，先搜索中心位置(标号为0)的周围8个位置(图中标号1，其距搜索中心0的距离为搜索区域半径的一半，这里称为步长)，找到最大互相关位置(图中用方框标出)；然后，以该位置为搜索中心再搜索其周围的8个位置(图中标号为2，每个位置距其搜索中心的步长为上次步长的一半)，以此类推直到搜索步长为0为止。然后，用高分辨率的图像和模板在以上一级搜索得到的最大互相关位置为中心的位置进行小区域的类似操作。假设搜索区域为S_M×S_M，采用全搜索要进行S_M×S_M次图像归一化互相关计算，如果采用上面的两级三步搜索法将只需要8log₂S_M次搜索，极大地减少了相关计算次数，提高搜索速度。如果采用更多级的搜索，那么，其搜索次数还可以减少。

(a) 低分辨率；(b) 高分辨率

图1  分级三步搜索法示意图

Fig.1  Schematic diagram of three-step hierarchical searching

2  抗尺度变化和旋转不变性的快速浮选泡沫图像序列运动估计

直接用图像的灰度进行互相关匹配运算，只能计算出无角度、无尺度变化的图像序列的运动特征参量。然而，在机械搅拌式浮选机浮选过程中，每个浮选槽中都有搅拌机在不停地搅拌，以增加矿物粒子和泡沫的接触机会，这样，浮选槽中的泡沫不可避免的发生形变即角度旋转和尺度变化。对发生了几何变化的泡沫图像序列进行模板匹配，会影响宏块跟踪准确   性。因此，有必要找出一种更好的匹配算法进行宏块跟踪，这种方法要能够迅速而准确地估计出发生角度旋转和尺度形变的泡沫图像的运动参量。为此，本文作者研究一种快速的抗尺度变化和旋转不变性的模块匹配方法。

2.1  具有旋转和尺度不变性的快速模块匹配法

2.1.1  2个基于圆投影的统计量

为了在模板匹配过程中处理目标旋转和尺度变化的情况，引入圆投影作为预处理过程，它将二维的图像转变成2个一维向量。设模板为T(x, y)，大小为    M×N。首先，将其从直角坐标T(x, y)转换到极坐标LPT(ρ, θ)(其中：ρ表示半径(以图像的中心为原点，ρ≤)，θ为角度轴，θ∈[1, 2, 3, …, n])。

然后，定义2个新的一维统计量S_T(ρ)和G_T(θ)。它们的计算公式如下：

也就是说，定义圆投影向量S_T(r)为沿着半径方向到模板中心距离为r的圆周上像素灰度之和；定义径向投影向量G_T(r)为角度等于θ时半径方向的像素灰度之和。由于S_T(r)记录了一组同心圆圆周上的灰度特征值，所以，它具有旋转不变性；而G_T(θ)记录了径向的灰度特征值，所以，它对含尺度缩放变换的图像对具有较强的相关性。圆投影概念如图2所示。

图2  圆投影概念

Fig.2  Circular projection conception

通过将传统的模板进行极坐标变换(对要进行搜索的当前帧也分块并进行同样的转换，计算出对应的圆投影向量S_F(ρ)和径向投影向量G_F(θ))，然后，用S_T(ρ)和G_T(θ)这2个分别具有旋转不变性和尺度不变性的投影向量进行相关性计算，就能获得准确的匹配位置，从而进行正确的宏块跟踪。具有旋转和尺度变化的图像对的统计向量S_T(ρ)和G_T(θ)对比曲线如图3所示。通过提取图像中泡沫子块的S_T(ρ)和G_T(θ)统计向量的一维变化趋势也可以看出，它们具旋转和尺度不变性。利用这2个统计量进行的相关计算将获得抗任意角度旋转和抗尺度变化的精确匹配结果。

(a) 模板图像；(b) 模板图像经放大1.2倍并旋转15?的参考图像；(c) 模板图像和参考图像的S_T(ρ)对比曲线；(d) 模板图像和参考图像的G_T(θ)对比曲线

图3 具有旋转和尺度变化的图像对的统计向量S_T(ρ)和G_T(θ)对比曲线

Fig.3 Curves comparison of statistic S_T(ρ) and G_T(θ) of image-pair with revolution and scaling

2.1.2  相关性计算

在进行模板匹配时，仍然采用相关函数进行模块匹配。假设：

，

。

它们分别代表模板和待匹配图子图的圆投影向量。S_T(ρ)和S_f(ρ)之间的相关性计算表示为：

统计向量G_T(θ)和G_F(θ)的相关性C_TFθ的计算也采用类似方法。最终的宏块相似度可以度量为：

2.2  二维拉格朗日插值方法的亚像素位移提取

用上面提出的具有旋转不变性的模板匹配算法进行的宏块跟踪来估算泡沫图像序列的运动速度只能估算出整数像素的位移，然而，由于在实际的现场图像采集中，摄像头的轻微抖动(现场由于搅拌机等的运行，振动较大)或者图像分辨率的局限性可能使泡沫的位移不是整数像素，而在应用中要计算速度的紊乱程度，因此，泡沫速度的计算精度至关重要。为了获得更为精确、有用的泡沫图像运动特征，有必要对泡沫图像亚像素位移进行估计。目前，能进行亚像素估计的算法主要有插值法、差分迭代法、双傅里叶变换和基于梯度的方法等^[15]。所有这些方法中，要满足精确性和实时性要求是关键。受文献[16]的启发，引入基于拉格朗日插值法的曲面拟合进行亚像素配准，具体操作如下：对已获得的整像素匹配点及其周围8点的相关系数矩阵，拉格朗日插值法插值为以下形式的四次曲面：

该非线性方程组可以通过Newton迭代法^[17]求出其极值点位置(x_max, y_max)，迭代的初值估计为(0, 0)。(x_max, y_max)为具有亚像素精度的宏块跟踪位置，利用这个位置和上一帧取得模板的位置之差，就可以得到具有亚像素精度的位移向量。另外，若将相关系数取自然对数后再进行插值，则结果更精确^[18]。因此，在具体插值之前，都先对相关系数取自然对数，再进行迭代计算。

3  浮选泡沫图像速度特征提取及结果分析

3.1  速度特征提取准确性测试

为了验证具有亚像素位移图像的速度估计的精确性，采用计算机自动生成已知像素位移的图像进行实验。实验步骤如下：对大小为256×256的泡沫图像上采样并平移整数像素(x, y)后再采样，这样，就能产生具有亚像素位移的实验图像，同时，图像也进行不同角度的旋转。表1所示为其中4次亚像素配准实验结果。从表1可见：采用本文算法能取得高精度的泡沫位移。

表1  亚像素配准实验结果

Table 1  Experiment results of sub-pixel registration

3.2  浮选泡沫图像序列速度特征分析

以某选矿厂轻金属浮选泡沫图像序列为研究对象，在2008-04-02对浮选泡沫图像序列进行速度特征实时提取。在9:00，使用工业摄像机在粗选槽采集一段视频并进行速度提取，其所提取的图像速度特征如图4所示。在16:00采集同一槽所得的浮选泡沫图像序列，其提取的泡沫图像速度特征如图5所示。

图4  某选矿厂4月2日9:00粗选槽浮选泡沫图像速度特征曲线

Fig.4  Flotation froth image sequences velocity curve of rougher cell at 9:00 on April 2

图5  某选矿厂4月2日16:00粗选槽浮选泡沫图像速度特征曲线

Fig.5  Flotation froth image sequences speed characteristic curve of rougher cell at 16:00 on April 2

由于所选用的拍摄帧率为15帧/s，所以，每帧代表的时间是1/15 s)，根据相机的焦距、物距计算出每个像素代表宽为0.041 5 mm，高为0.041 5 mm。从图4和图5可以看出：浮选泡沫的速度以某个值为中心的周围浮动，并且具有一定的周期性。由于矿浆的流动、矿粒和泡沫的接触主要由叶轮搅拌引起，它的周期特性与搅拌机搅拌的频率等有直接联系。浮选泡沫运动的速度可以表示为：

为了能比较不同速度均值下浮选泡沫运动的紊乱程度，定义

在实际分析中，常采用一段时间的平均紊乱程度来衡量浮选泡沫运动速度的紊乱程度。选泡沫运动的紊乱程度是搅拌叶轮雷诺数的具体体现^[19]，矿浆叶轮雷诺数定义为：

叶轮雷诺数的变化直接影响泡沫区矿粒和气泡的稳定漂浮，较小的叶轮雷诺数会使矿物粒子与泡沫不能充分接触，但太大的叶轮雷诺数会减弱矿物粒子在泡沫上的黏附程度，反而减少了泡沫中的矿物含量。从式(10)可以看出，由于叶轮直径和叶轮转速基本固定，因此，当叶轮雷诺数过大时，可以通过降低矿浆密度或者增加矿浆黏度(实际上就是控制加药量)来调整叶轮雷诺数，从而减小浮选泡沫速度的紊乱程度。

由于pH值发生波动(实际是加药量引起的)引起浮选槽内的矿浆浓度发生变化，而电动搅拌机的搅拌频率基本没变化。从表2所示的速度特征统计结果可以看出：在16:00的泡沫速度紊乱程度大于9:00时的泡沫运动紊乱程度，从而也增大了矿物粒子的脱附浮选泡沫的概率，结果是尾矿中的金属含量变大，精矿品位降低，从而使16:00时的浮选效率降低。

表2  所提取的速度特征统计比较

Table 2  Statistic comparison of extracted velocity vector

对该选矿厂轻金属矿物浮选进行长时间观察，分析浮选泡沫速度的紊乱程度和浮选工艺参数之间的关系，发现浮选泡沫速度紊乱程度和浮选主要的性能指标即尾矿(金属氧化物/氧化硅)指标和精矿(金属氧化物/氧化硅)指标的关系如图6所示。从图6可以看出：

(a) 尾矿；(b) 精矿

图6  浮选泡沫速度紊乱程度和浮选性能间的关系

Fig.6  Relationship between froth velocity turbulence degree and flotation performance

浮选泡沫的紊乱程度增加使矿物回收率降低，矿物的品位下降，表现为尾矿金属氧化物含量升高和总精矿金属氧化物含量下降。可见，在矿物矿浮选过程中，要提高浮选效率，在满足矿物粒子和浮选泡沫充分接触时，要尽可能降低浮选泡沫速度的紊乱程度，也就是尽量减少已黏附在泡沫上的矿物粒子的脱附概率，从而有效地提高浮选效率。

4  结  论

a. 在浮选机的搅拌下，浮选泡沫不可避免地发生角度旋转和尺度变化。抗尺度变化和具有旋转不变性的浮选泡沫图像速度特征提取方法，解决了这种具有角度旋转和尺度变化的浮选泡沫图像的速度特征精确提取的问题。

b. 在轻金属浮选生产中，尽量减少泡沫运动速度的紊乱程度，也就是说，在保证矿物粒子和泡沫充分接触的前提下尽量减少已黏附在泡沫上的矿物粒子的脱附概率。

c. 采用该速度特征提取方法能实时采集矿物浮选泡沫图像，提取速度特征，观察速度紊乱程度，为浮选生产提供指导，保证矿物浮选的优化稳定运行，从而提高浮选生产效率。该方法可推广到不同的矿物浮选生产中，提高金属矿物浮选的自动化水平。

参考文献：

[1] 张 勇, 王 峰, 潘学军, 等. 粗糙集理论在阳离子反浮选控制中的应用[J]. 中南工业大学学报: 自然科学版, 2003, 34(4): 368-372.

ZHANG Yong, WANG Feng, PAN Xue-jun, et al. Rough set theory and its application in cation anti-flotation control[J]. Journal of Central South University: Natural Science, 2003, 34(4): 368-372.

[2] 刘文礼, 路迈西, 王 勇, 等. 媒泥浮选泡沫图像纹理特征的提取及泡沫状态的识别[J]. 化工学报, 2003, 54(6): 832-836.

LIU Wen-li, LU Mai-xi, WANG Yong, et al. Extraction of textural feature and recognition of coal flotation froth[J]. Journal of Chemical Industry and Engineering, 2003, 54(6): 832-836.

[3] 何桂春, 黄开启. 浮选指标与浮选泡沫数字图像关系研究[J]. 金属矿山, 2008(8): 96-101.
HE Gui-chun, HUANG Kai-qi. Study of the relation between flotation indexes and froth digital images[J]. Metal Mine, 2008(8): 96-101.

[4] 阳春华, 杨尽英, 牟学民, 等. 基于聚类预分割和高低精度距离重构的彩色浮选泡沫图像分割[J]. 电子与信息学报, 2008, 30(6): 1286-1290.
YANG Chun-hua, YANG Jin-ying, MOU Xue-min, et al. A segmentation method based on clustering per-segmentation and high-low scale distance reconstruction for color froth image[J]. Journal of Electronic & Information Technology, 2008, 30(6): 1286-1290.

[5] Kaartinen J, Hatonen J, Hyotyniemi H, et al. Machine- vision-based control of zinc flotation: A case study[J]. Control Engineering Practice, 2006, 14(12): 1455-1466.

[6] Hasu V. Image analysis in mineral flotation[R]. Finland: Control Engineering Laboratory, Helsinki University of Technology, 2002.

[7] Lin B, Recke B, Knudsen J K H, et al. Bubble size estimation for flotation processes[J]. Minerals Engineering, 2008, 21(7): 539-548.

[8] Liu J J, Macgregor J F, Duchesne C, et al. Flotation froth monitoring using multiresolutional multivariate image analysis[J]. Minerals Engineering, 2005, 18(1): 65-76.

[9] Citir C, Aktas Z, Berber R. Off-line image analysis for froth flotation of coal[J]. Computers & Chemical Engineering, 2004, 28(5): 625-632.

[10] ZHANG Hai-yan, ZHOU Mai-li. Objects tracking based on optical flow in polar-log images[J]. Journal of Computational Information System, 2006, 2(2): 829-832.

[11] XU Jie-bin, PO Lai-man, Cheung C K. Adaptive motion tracking block matching algorithms for video coding[J]. IEEE Trans on Circuits and Systems for Video Technology, 1999, 9(7): 1025-1029.

[12] Balci M, Foroosh H. Subpixel estimation of shifts directly in the Fourier domain[J]. IEEE Trans on Image Processing, 2006, 15(7): 1965-1972.

[13] Wang H, Mersereau R. Fast algorithms for the estimation of motion vectors[J]. IEEE Trans on Image Processing, 1999, 8(3): 435-438.

[14] Li R, Zeng B, Lion M L. A new three-step search algorithm for block motion estimation[J]. IEEE Trans Circuits Systems for Video Technology, 1994, 4(4): 438-442.

[15] ZHANG Jun, JIN Guan-chang, MA Shao-peng, et al. Application of an improved subpixel registration algorithm on digital speckle correlation measurement[J]. Optics & Laser Technology, 2003, 35(7): 533-542.

[16] 王会峰, 冀 芳, 汪大宝. 激光基准成像测量光斑图像的亚像素检测算法[J]. 量子电子学报, 2008, 25(4): 402-406.
WANG Hui-feng, JI Fang, WANG Da-bao. A sub-pixel edge detection algorithm for laser spot in measurement system using laser beam as straight datum[J]. Chinese Journal of Quantum Electronics, 2008, 25(4): 402-406.

[17] Giusti M, Lecerf G, Salvy B. A grobner free alternative for polynomial system solving[J]. Journal of Complexity, 2001, 17(1): 154-211.

[18] 孟利波, 马少鹏, 金 观. 数字散斑相关测量中亚像素位移测量方法比较[J]. 实验力学, 2003, 18(3): 343-348.
MENG Li-bo, MA Shao-peng, JIN Guan. On the performance of the subpixel displacement estimations used in digital speckle correlation method(DSCM)[J]. Journal of Experimental Mechanics, 2003, 18(3): 343-348.

[19] 曾克文, 余永富. 浮选矿浆紊流强度对矿物浮选的影响[J]. 金属矿山, 2000(9): 17-20.
ZENG Ke-wen, YU Yong-fu. Effect of the turbulent strength of the flotation pulp on mineral flotation[J]. Metal Mine, 2000(9): 17-20.

收稿日期：2009-01-28；修回日期：2009-04-28

基金项目：国家自然科学重点基金资助项目(60634020)；国家自然科学基金资助项目(60874069)

通信作者：唐朝晖(1965-)，男，湖南长沙人，博士，教授，从事智能控制、智能信息处理、复杂工业过程故障诊断等研究；电话：0731-88879628；E-mail: zhtangcq@hotmail.com