CN103604946B - 一种面向矿物泡沫浮选的嵌入式泡沫流速测量系统及方法 - Google Patents

Abstract

本发明公开了一种面向矿物泡沫浮选的泡沫流速测量系统及方法。本发明构建了嵌入式浮选泡沫采集平台，以工业摄像机获取矿物浮选过程中泡沫图像为基础，针对采集到的图像在DSP处理器内完成图像的预处理，利用当前帧和参考帧中图像块的相关性分析方法提取了泡沫速度特征。本发明所阐述的嵌入式解决方案通过分析泡沫图像得到泡沫速度特征，作为一种生产状态重要的参考指标，可实时反映生产状况，具有体积小、成本低、功耗低、可靠性高，特别适用于选矿厂等环境恶劣、难以测量泡沫流速的工业现场。

Description

一种面向矿物泡沫浮选的嵌入式泡沫流速测量系统及方法

技术领域

本发明属于嵌入式图像处理领域，具体涉及一种面向矿物泡沫浮选的泡沫流速测量系统及测量方法。

背景技术

矿物浮选是矿物加工中应用最广泛的一种选矿方法，通过极其复杂的物理化学反应过程进行矿化起泡。应用矿物浮选是为了提高原矿品位，达到矿物还原冶炼的要求。多年来矿物浮选主要靠技术工人通过观察浮选槽中泡沫的特征，手动调节浮选操作，而调节过程由于人员的轮换以及调节的主观性，不能对泡沫特征进行准确的判断，从而导致操作的主观性和随意性较大，使得浮选过程难以处于最优生产状态。

浮选过程中产生形态和大小不等的各种气泡来携带矿粒，浮选泡沫是包含有气、液、固的三相复杂体，流动性强，容易破裂。浮选过程中泡沫流速是一个重要的指标。泡沫流速能够直接反映矿物浮选过程中矿物泡沫的状态及液位的信息，而矿物泡沫状态与矿物的浮选产品品位直接相关，泡沫的液位则直接决定了工况正常与否以及所需矿物的有效回收。由于泡沫位于浮选过程中矿浆表面，含有大量细小的矿物颗粒，传统的流速传感器只能测量矿浆的流速，很难直接测量泡沫的流速，采用机器视觉方法是可行的方案。浮选过程往往环境较恶劣，各种噪音、电磁波、震动和粉尘等污染使得传统的计算机很难应用在浮选流速测量上，而嵌入式DSP具有实时性好、抗干扰能力强和能耗低的优点，非常适合用来设计专用的泡沫流速测量系统。

发明内容

本发明的目的在于解决矿物浮选过程中泡沫流速难以用传统流速传感器测量的问题，提出了一种用于矿物浮选过程中的嵌入式泡沫流速测量系统及方法。

本发明的技术方案是，一种面向矿物泡沫浮选的嵌入式泡沫流速测量系统，包括嵌入式视频处理模块、光源、工业摄像机、交换机、计算机和浮选槽，所述的光源设置于浮选槽上方，所述的工业摄像机与光源平行安装且与光源高度相同，所述的嵌入式视频处理模块、交换机和计算机依次通信串联。

一种使用上述装置的面向矿物泡沫浮选的嵌入式泡沫流速测量方法，包括以下步骤：

步骤一：光源照射泡沫表面，工业摄像机拍摄PAL制式视频，并将视频传输至视频解码单元；

步骤二：视频解码单元将模拟图像信息转换成数字图像数据，并存储到SDRAM中；

步骤三：对存储到SDRAM中的图像数据进行处理，得出泡沫流速；

步骤四：将得到的泡沫流速通过网口传输到上位机数据库中。

所述的一种方法，所述的步骤三的步骤包括：

步骤(1)：将SDRAM存储的当前帧和前一帧图像进行灰度化处理，并将灰度化后的两帧图像存储在SDRAM中；

步骤(2)：对转换后的灰度图像采用直方图均衡化增强图像对比度，统计图像中每个灰度级的像素出现的概率并通过对其均匀处理的变换函数对图像进行灰度变换，最终得到直方图均衡化后的图像；

步骤(3)：选取当前帧的前一帧图像为参考帧，将参考帧和当前帧图像以像素为基本单位划分为64*64大小子区域，分别对参考帧和当前帧的子区域在水平和垂直方向灰度投影并对相同子区域进行相关运算，得到当前帧所有子区域相对于参考帧水平和垂直方向的运动矢量；在所有子区域运动矢量的基础上，根据六参数仿射模型和最小二乘法计算当前帧相对于参考帧的全局运动矢量；由计算得到的水平和垂直两个方向的运动矢量能够计算出最终的泡沫移动速度，连续计算一秒钟内所有25帧图像中泡沫的速度，求取平均速度即为流速。

本发明的技术效果在于，通过嵌入式设备测得的泡沫流速，有效的解决了矿物浮选过程中泡沫流速难以测定的问题，测量精度能够达到工业现场要求。本发明解决了传统流速传感器难以测量矿物浮选过程中泡沫流速的问题，具有体积小、成本低的特点。本发明采用平台为DM642的嵌入式系统，能够很好的在现场具有较强震动、酸雾腐蚀、电磁干扰的环境下工作，保证了测量的准确性和可靠性。本发明针对嵌入式系统的特点设计了专用的算法，即分块投影匹配算法。投影算法较传统的块匹配算法减少了算法的时间消耗，而分块求取运动矢量并根据摄像机模型求解全局矢量的方法提高了计算的精确度，能够满足现场泡沫流速测量中对实时性和测量精度的要求。

附图说明

图1是矿物浮选流速测量方法硬件结构示意图；

图2是嵌入式视频处理模块的硬件结构图；

图3是泡沫流速测量算法流程图；

其中1为嵌入式视频处理模块、2为光源、3为工业摄像机、4为同轴电缆、5为路由器、6为上位计算机、7为泡沫层、8为浮选槽、9为泡沫槽。

具体实施方式

矿物浮选泡沫流速测量系统硬件结构如图1所示，主要由嵌入式视频处理模块1、光源2、工业摄像机3、同轴电缆4、用于通信数据交换传输的路由器5和上位计算机6构成。浮选槽是矿浆与药物混合反应发生的场所，细化泡沫，使得矿物粘合泡沫之上，浮到矿浆面并形成矿化泡沫，其表面是泡沫产生最多的区域，具有良好的隔离性。泡沫槽是盛放由刮板在浮选槽表面刮出的矿化泡沫的器皿，能够进一步提高矿浆中含矿量，为下一步精选提供品位更高的矿浆。泡沫槽是盛放浮选槽表面刮出的矿化泡沫的容器。工业CCD摄像机3垂直安装在浮选槽上方，与嵌入式视频处理模块通过视频输出线相连接。嵌入式视频处理模块选用TI公司的DM642作为核心处理芯片，视频解码部分选用A/D视频输入芯片TVP5158，配置总线选用芯片PCF8584作为I2C总线控制器，SDRAM存储器采用型号为HY57V283220T的内存芯片。板上还提供McASP与PCI总线用于芯片之间数据的传输。摄像机与浮选槽8表面泡沫的垂直高度为70-110cm。系统光源2选用200W的高频荧光灯，与摄像机3的水平距离为5-10cm，保证摄像机采集图像所需的充足稳定的光照。摄像机3实时采集泡沫层7的图像并传输到嵌入式视频处理模块上的TVP5158视频解码单元，并通过解码芯片将模拟视频信号转换为YUV的数字图像信号，存储在嵌入式视频处理模块的SDRAM中。DM642的核心处理芯片根据流速测量算法对SDRAM中存储的多帧数字图像进行处理，最终得到当前浮选槽8中矿浆的流速。嵌入式视频处理模块利用以太网将计算得到的实时矿浆流速数据通过同轴电缆4传输到选矿厂的内部路由器5上，最终传输到终端处理计算机6并存储到数据库中。

泡沫流速测量方法具体实现如下：

对嵌入式视频处理模块上型号为HY57V283220T的SDRAM中存储的多帧泡沫图像，提取当前帧和前一帧两帧图像，对两帧图像进行灰度化处理。由于工业摄像机采集的视频信号由嵌入式视频处理模块上的视频解码单元TVP5158转换后的格式为YUV，而亮度信号Y可以直接表征每一帧图像的像素的灰度值，从而直接存储单帧图像像素的Y值完成对图像的灰度化处理。DM642内部硬件结构如图2所示，视频解码单元TVP5158解码后得到的Y值通过视频传输接口传输到DM642内部的以太网控制器中，以太网控制器通过EMIF接口将单帧Y值数据存储在SDRAM内存中，DM642中的C64x内核将要计算的数据提取到高速二级缓存及一级缓存中并将计算完成后的结果存储在SDRAM内。泡沫流速的计算需在DM642内核中以流水线方式完成若干次数据处理，最终计算结果存储在SDRAM中。

对灰度化处理后的两帧图像进行直方图均衡化处理。直方图均衡化是使原直方图变换为具有均匀密度分布的直方图，然后按该直方图调整原图像的一种图像处理技术，能够增强图像全局对比。变换函数T(r)是原图像的累积分布函数并满足条件：(1)T(r)在0≤r≤1范围内是单值单增的函数；(2)对0≤r≤1有0≤T(r)≤1。设N为待处理帧图像中像素的总数，n_k为第k级灰度的像素数，r_k表示第k个灰度级，则该帧灰度级为r_k的像素出现概率P(r_k)的计算方法见公式(1)，对其直方图均衡化后的像素的灰度值s_k见公式(2)。直方图均衡化后的图像数据存储在SDRAM中的特定内存区域，不同帧图像存储的内存段不同，存储的图像数据仅限于最新两帧。

$P\left(r_k\right)=\frac{n_k}{N}(k=\mathrm{0,1,2},...,L-1)---\left(1\right)$

$s_k=T\left(r_k\right)=\underset{j=0}{\overset{k}{\mathrm{\Σ}}}{P}_r\left(r_j\right)=\underset{j=0}{\overset{k}{\mathrm{\Σ}}}\frac{n_j}{N}---\left(2\right)$

为了提高DM642处理的实时性，数据的存储采用定点方式，采用移位方式替代除法运算，优化图像处理速度，根据处理图像大小动态分配SDRAM中存储空间。

基于分块投影匹配的泡沫流速测量算法具体流程如图3所示，主要分为图像分块、分块图像投影、分块运动矢量计算、全局运动矢量计算及求取25帧泡沫流速平均值几个步骤组成。将直方图均衡化后的两帧图像进行分区，区域大小取决于处理帧图像大小，对于工业摄像机摄取的PAL制式视频单帧图像大小为720*576,选取子区域大小为64*64将其划分为若干个图像块。完成子区域划分后，对每一个子区域做水平和垂直两个方向的灰度投影，具体投影变换计算如公式(3)所示。

公式(3)中，Y_k(i,j)表示第k个子区域图像上(i,j)处的灰度值，Col_k(j),Row_k(i)分别表示第k个子区域图像第j列、第i行所有像素的灰度值之和。灰度累加计算在嵌入式视频处理模块的核心处理芯片上完成，灰度值以定点数表示。计算得到的临时数据暂时存储在嵌入式视频处理模块的片上RAM中，提高了计算的速度。子区域投影计算后的数据分别保存在SDRAM中的特定内存段中，当前帧与前一帧相同子区域之间的运动矢量根据子区域水平和垂直方向的投影曲线进行相关计算，根据两条相关曲线的谷点确定当前帧相对于前一帧的运动矢量。水平方向互相关运算如公式(4)所示，垂直方向互相关运算如公式(5)所示。

$\begin{array}{c}{\mathrm{Col}}_k\left(j\right)=\underset{i}{\mathrm{\Σ}}{Y}_k(i,j)\\ {\mathrm{Row}}_k\left(i\right)=\underset{j}{\mathrm{\Σ}}{Y}_k(i,j)\end{array}---\left(3\right)$

$\begin{array}{cc}B\left(i\right)=\underset{j=1}{\overset{d}{\mathrm{\Σ}}}\[G_{t+k}(j+i+1)-G_t(m+j){\]}^2& 1\≤i\≤2m+1---\left(4\right)\end{array}$

$\begin{array}{cc}C\left(j\right)=\underset{j=1}{\overset{e}{\mathrm{\Σ}}}\[G_{t+k}(j+i-1)-G_t(m+i){\]}^2& 1\≤j\≤2m+1---\left(5\right)\end{array}$

其中，公式(4)中，G_t(j)、G_t+k(j)分别为第t帧和t+k帧对应子区域第j列的投影，即第j列所有像素灰度值的和，B(j)为第t帧和t+k帧对应子区域第j列的投影值的方差，j的取值范围为当前投影子区域的列数,即第1列到第d列。公式(5)中，G_t(i)、G_t+k(i)分别为第t帧和第t+k帧对应子区域的第i行的投影，即第i行所有像素灰度值的和，C(j)为第t帧和第t+k帧对应子区域的第i行的投影值的方差,i的取值范围为当前投影子区域的行数，即第1列到第e列。公式(4)、(5)中的m为最大偏移量的值，设当I_min、J_min分别为B(i)、C(j)最小时i、j的取值，则第t帧子区域相对于第t+k帧在水平、垂直方向的位移分别为公式(6)、(7)所示。

δ_x=m-1-I_min(6)

δ_y=m-1-J_min(7)

得到各个子区域的运动矢量后，利用摄像机6参数仿射模型和最小二乘法计算当前帧相对于前一帧的全局运动矢量。6参数仿射模型由摄像机的参数模型推导而来。在摄像机坐标体系内，描述三维物体运动位置的变化如公式(8)所示。公式(8)中，三维物体从t_k时刻点(X,Y,Z)移动到t_k+1时刻点(X',Y',Z')，这种运动可以分解为矩阵R描述的旋转和矢量T描述的平移。r_ij(i,j=0…9)为旋转矢量矩阵R的元素，t₁,t₂,t₃分别为点(X,Y,Z)相对于点(X',Y',Z')在三方向上的平移矢量。

$\left[\begin{array}{c}x\text{'}\\ y\text{'}\\ z\text{'}\end{array}\right]=R\left[\begin{array}{c}x\\ y\\ z\end{array}\right]+T=\left[\begin{array}{ccc}{r}_{11}& r_{12}& r_{13}\\ r_{21}& r_{22}& r_{23}\\ r_{31}& r_{32}& r_{32}\end{array}\right]\left[\begin{array}{c}x\\ y\\ z\end{array}\right]+\left[\begin{array}{c}{t}_1\\ t_2\\ t_3\end{array}\right]---\left(8\right)$

图像上的点是三维物体在二维平面的投影，二维成像平面上的目标点坐标(x,y)和三维点坐标(X,Y,Z)之间的关系如公式(9)所示。

$\begin{array}{cc}(x,y)=(X,Y)& (x\text{'},y\text{'})=(X\text{'},Y\text{'})---\left(9\right)\end{array}$

联立公式(8)和公式(9)可以推出坐标点的求解函数，如公式(10)所示。

$\begin{array}{c}x\text{'}=r_{11}x+r_{12}y+(r_{13}Z+t_1)\\ y\text{'}=r_{21}x+r_{22}y+(r_{23}Z+t_2)\end{array}---\left(10\right)$

将公式(9)中的参数提取出来，可以表示为一般化形式，如公式(11)所示。

$\begin{array}{c}x\text{'}=\mathrm{ax}+\mathrm{by}+c\\ y\text{'}=\mathrm{dx}+\mathrm{ey}+f\end{array}---\left(11\right)$

当前帧和前一帧对应子区域像素点之间的运动关系如公式(12)所示。

$\begin{array}{c}{x}_i^{\′}+\mathrm{\Δ}{x}_i=x_i\\ y_i^{\′}+\mathrm{\Δ}{y}_i=y_i\end{array}---\left(12\right)$

公式(12)中i为子区域编号，Δx_i、Δy_i为当前帧第i个子区域图像相对于前一帧的水平、垂直方向运动矢量。根据公式(12)定义能量函数如公式(13)所示：

$\begin{array}{c}{E}_x=\underset{i=0}{\overset{n}{\mathrm{\Σ}}}{(a{x}_i+b{y}_i+c-x_i^{\′})}^2\\ E_y=\underset{i=0}{\overset{n}{\mathrm{\Σ}}}{(d{x}_i+e{y}_i+f-y_i^{\′})}^2\end{array}---\left(13\right)$

公式(13)中n为子区域的数量，当能量值E_x、E_y得最小值时得到的参数即为全局运动参数，由最小二乘法可得计算公式(14)：

$\frac{\∂E_x}{\∂a}=\frac{\∂E_x}{\∂b}=\frac{\∂E_x}{\∂c}=\frac{\∂E_y}{\∂d}=\frac{\∂E_y}{\∂e}=\frac{\∂E_Y}{\∂f}=0---\left(14\right)$

从而由公式(14)能够得到一组方程如公式(15)所示。

$\left\{\begin{array}{c}a\underset{i=0}{\overset{n}{\mathrm{\Σ}}}{x}_i^2+b\underset{i=0}{\overset{n}{\mathrm{\Σ}}}{y}_i{x}_i+c\underset{i=0}{\overset{n}{\mathrm{\Σ}}}-\underset{i=0}{\overset{n}{\mathrm{\Σ}}}{x}_i{x}_i^{\′}=0\\ a\underset{i=0}{\overset{n}{\mathrm{\Σ}}}{x}_i{y}_i+b\underset{i=0}{\overset{n}{\mathrm{\Σ}}}{y}_i^2+c\underset{i=0}{\overset{n}{\mathrm{\Σ}}}{y}_i-\underset{i=0}{\overset{n}{\mathrm{\Σ}}}{y}_i{x}_i^{\′}=0\\ a\underset{i=0}{\overset{n}{\mathrm{\Σ}}}+b\underset{i=0}{\overset{n}{\mathrm{\Σ}}}{y}_i+\mathrm{nc}-\underset{i=0}{\overset{n}{\mathrm{\Σ}}}{x}_i^{\′}=0\\ d\underset{i=0}{\overset{n}{\mathrm{\Σ}}}{x}_i{y}_i+e\underset{i=0}{\overset{n}{\mathrm{\Σ}}}{y}_i^2+f\underset{i=0}{\overset{n}{\mathrm{\Σ}}}{x}_i-\underset{i=0}{\overset{n}{\mathrm{\Σ}}}{y}_i{x}_i^{\′}=0\\ d\underset{i=0}{\overset{n}{\mathrm{\Σ}}}{x}_i{y}_i+e\underset{i=0}{\overset{n}{\mathrm{\Σ}}}{y}_i^2+f\underset{i=0}{\overset{n}{\mathrm{\Σ}}}{y}_i-\underset{i=0}{\overset{n}{\mathrm{\Σ}}}{y}_i{y}_i^{\′}=0\\ d\underset{i=0}{\overset{n}{\mathrm{\Σ}}}{x}_i+e\underset{i=0}{\overset{n}{\mathrm{\Σ}}}{y}_i+\mathrm{nf}-\underset{i=0}{\overset{n}{\mathrm{\Σ}}}{y}_i^{\′}=0\end{array}\right\}---\left(15\right)$

由公式(15)的方程组可以求解得到运动参数矢量{a,b,c,d,e,f}，进而求得代表泡沫流速的矢量{c,f}，求取一秒内连续25帧泡沫图像运动矢量并计算平均值，作为流速平均速度。求解的数学运算过程调用TI的数学函数库，从而大大提高求解的速度，满足实时性要求。

Claims (1)

1.一种面向矿物泡沫浮选的嵌入式泡沫流速测量方法，其特征在于，使用一种面向矿物泡沫浮选的嵌入式泡沫流速测量系统，包括嵌入式视频处理模块、用于照射浮选槽内泡沫的光源、用于拍摄浮选槽内泡沫图像的工业摄像机、上位计算机和浮选槽，所述的光源设置于浮选槽上方，所述的工业摄像机与光源平行安装且与光源高度相同，所述的嵌入式视频处理模块与计算机通信连接，所述的嵌入式视频处理模块包括处理器、视频解码单元和存储单元，所述的处理器分别与视频解码单元和存储单元通信连接，所述的视频解码单元通信连接工业摄像机，处理器通信连接至上位计算机；

包括以下步骤：

步骤一：光源照射泡沫表面，工业摄像机拍摄被照射的泡沫表面并输出PAL制式视频传输至视频解码单元；

步骤二：视频解码单元将模拟图像信息转换成数字图像数据，并存储到存储单元中；

步骤三：处理器对存储到存储单元中的图像数据进行处理，得出泡沫流速；图像数据处理步骤包括：

步骤(1)：将存储单元存储的当前帧和前一帧图像进行灰度化处理，并将灰度化后的两帧图像存储在存储单元中；

步骤(2)：对转换后的灰度图像采用直方图均衡化增强图像对比度，统计图像中每个灰度级的像素出现的概率并通过对其均匀处理的变换函数对图像进行灰度变换，最终得到直方图均衡化后的图像；

步骤(3)：选取当前帧的前一帧图像为参考帧，将参考帧和当前帧图像以像素为基本单位划分为64*64大小子区域，分别对参考帧和当前帧的子区域在水平和垂直方向灰度投影并对相同子区域进行相关运算，得到当前帧所有子区域相对于参考帧水平和垂直方向的运动矢量；在所有子区域运动矢量的基础上，根据六参数仿射模型和最小二乘法计算当前帧相对于参考帧的全局运动矢量；由计算得到的水平和垂直两个方向的运动矢量能够计算出最终的泡沫移动速度，连续计算一秒钟内所有25帧图像中泡沫的速度，求取平均速度即为流速。

步骤四：将得到的泡沫流速传输到上位计算机的数据库中。