颗粒物质分选过程中的速度和角度分布的测量与研究

陶子寅， 李 然， 陈 泉， 修文正， 华云松， 杨 晖

（上海理工大学 光电信息与计算机工程学院，上海 200093）

颗粒物质在自然界、日常生活及生产技术中普遍存在，是地球上存在最多、最为人们所熟悉的物质类型之一[1]。与流体运动相比，颗粒的运动非常复杂，兼具流体和固体的某些特征，但与它们又都不相同，颗粒流与流体重要的区别之一是粒料中存在与流体不同的混合与分离现象，特别是不同性质颗粒的分离和分层[2]。由于不同颗粒物质的混合与分聚对于社会生产中的诸如粮食干燥和储存、药品加工等生产工艺有着重要的影响，因此对于滚筒内颗粒的混合与分选过程的研究具有很重要的理论和实际意义：合理地应用颗粒的混合与分离作用能够节省能源与加工费用，同时也能在一定程度上促进颗粒力学的发展[3]。

针对不同颗粒所组成的颗粒体系的混合与分聚现象，近年来国内外学者作了许多研究。Kwapinska等[4]通过离散元模拟方法研究了旋转速度、转筒直径、装载量等参数对颗粒径向混合的影响。赵永志等[5]利用三维三方程线性弹性–阻尼离散单元模型讨论了滚筒转速、颗粒装载率等参数对薄滚筒内二元S型颗粒体系分层的影响，当转速较高时，滚筒内形成大颗粒在外、小颗粒在内、具有圆形界面的月亮模式，当转速较低时形成具有波浪形界面的花瓣模式。Arntz等[6]通过玻尔兹曼方程提出了一种通过混合熵来量化颗粒体系混合特性的方式。欧阳鸿武等[7]系统地研究颗粒混合物体系在半填充转鼓中的花瓣斑图形成过程，探讨规则花瓣斑图的形成条件。

本文在上述研究内容的基础上，针对二维转筒，创新性地将光流法和图像法应用于颗粒体系的表面流动层速度与倾斜角检测上，同时结合混合熵这一量化颗粒体系分选程度的指标，重点探究在不同转速下颗粒体系的混合熵、表面流动层平均速度和倾斜角三者之间的关系。

1 实验方法

1.1 实验装置

图1为本实验中所采用的实验装置，由二维滚筒、可调光源、步进电机和高速面阵CCD相机组成。滚筒采用有机玻璃材料加工制成，内径R=140 mm，筒壁厚度为5 mm，两个壁面之间的间隙d=10 mm。滚筒与步进电机相连，通过输出电压的方波频率来驱动电机，进而控制滚筒转速。实验过程中滚筒的转速范围为0～0.12 rad/s。将转筒与高速面阵CCD相机（3 F04 M型号）放置在同一水平面上，对转筒侧壁面进行拍照，获得颗粒体系的正视图。相机分辨率为 2320×1720，单位像元尺寸长度为7.4 μm，采样频率设置为330帧/s。转筒内填充的颗粒材料为两种颜色的电镀玻璃球形颗粒，绿色颗粒的均值粒径d1=2.6 mm，红色颗粒的均值粒径d2=1.25 mm。将两种颗粒以体积比1∶1进行填充，转筒内颗粒材料的填充度（颗粒体系体积与二维滚筒总体积的比值）为38%。

图 1 二维滚筒实验装置图Fig.1 Experimental device diagram for two-dimensional drum

图 2 表面流动层的区域划分示意图Fig. 2 Schematic diagram of the regional division on the surface flow layer

图2为二维滚筒内的颗粒体系流动层的位置划分示意图。图中划分了10个区域，每个区域长度为170像素，宽度为32像素，相邻两个区域中心位置的间隔为200像素，并分别命名为Ⅰ～Ⅹ；以流动层的中点位置为原点o建立直角坐标系，平行于表面流动层的方向为坐标x轴，垂直于流动层的方向为坐标y轴。

1.2 流动层颗粒的速度测量算法

目前光流法大致有两种计算方法：稠密光流和稀疏光流[7]。其区别在于稠密光流是将每个像素与速度关联；稀疏光流则是先指定一组具有明显特征的像素点进行匹配。因此稀疏光流相比于稠密光流不仅结果会更加稳定，相应的运算量也会减少许多。本文所采用的测速算法即是一种基于金字塔分层的Lucas-Kanade光流法，属于稀疏光流这一类，相应的特征像素点的选取采用Shi-Tomasi角点检测算法。

角点为图像中轮廓边界的相交点[8]，也被定义为在任意方向的一个微小变动都会引起灰度很大变化的点。现有的角点检测算法主要有基于边缘轮廓和基于灰度图像两类。其中由于基于灰度的角点检测算法避免了在提取轮廓时存在的误差，因而应用更为广泛[8]。本文所采用的是Shi-Tomasi角点检测算法，这是一种针对Harris角点检测[9]的改进算法，获取Harris角点中的强角点。

Harris角点检测算法是一种直接基于灰度图像的角点检测算法，其主要原理是对于强度信息为I(x,y)的图像，在像素点 (x,y)处将一个局部的小窗口w(x,y)沿任意方向移动一段微小位移(u,v)，w(x,y)为滑动窗口函数，通过引入灰度变化函数E(x,y)来表征图像的自相关性。Harris对于灰度变化函数E(x,y)的定义如下：

对位移后的I(x+u,y+v)进行泰勒展开，可得

O表示泰勒展开式中的高阶项。忽略高阶无穷小项，可将E(x,y)化简为二次型：

对于实对称矩阵M，通过对其对角化处理后可以得到它的两个特征值 λ1和 λ2，由于这两个特征值表征了图像在x轴和y轴的表面曲率，因此可以通过这两个特征值的大小来判别图像中的角点、边缘以及平坦区域。Harris在此基础上给出了角点响应函数（CRF）的定义：

式中， D et(M)表示矩阵M的行列式，Det(M)=λ1λ2； T r(M)表示矩阵M的迹， T r(M)=λ1+λ2。当目标点所求得的CRF值大于阈值时，Harris认为其为角点。Shi-Tomasi在Harris研究的基础上进行了改进，只选取矩阵M的两个特征值中较小的一个与阈值进行比较，如果大于阈值，即可以得到强角点。

图3为通过Shi-Tomasi算法检测到的区域角点结果。检测到的角点将作为后续Lucas-Kanade稀疏光流法的特征点基础。

图 3 Shi-Tomasi角点检测算法处理混合颗粒体系示意图Fig.3 Schematic diagram of Shi-Tomasi corner detection algorithm for processing the mixed particle system

光流场的计算最初由Horn和Schunck提出[10]，其约束方程的建立基于两个假设：相邻图像之间的灰度基本保持不变以及相邻图像之间的时间间隔很小，从而时间的变化不会引起位置的剧烈变化。假定某一时刻t图像上某一像素点 (x,y)所对应的灰度值为I(x,y,t)，经过一个很小的时间间隔Δt后，运动后的像素点所对应的灰度值变为I(x+Δx,y+Δy,t+Δt)，根据灰度基本保持不变的假设可得

将运动后像素点处的灰度值通过泰勒展开，并忽略高阶项后可得

式 中 ，Ix，Iy和It分 别 表 示 灰 度 信 号 在x，y和t上 的偏导。

仅有一个光流约束方程是无法确定u和v这两个光流分量的，必须引入其他的附加约束条件。于是Lucas和Kanade在此基础上引入了新的假设：假设在一个小的空间邻域 Ω 上运动矢量保持恒定。通过空间一致性的假设，可以考虑借助领域内的n个 像素点来建立n个光流约束方程用来估计光流分量。下面以一个 3 ×3的领域窗口为例，可以建立9个方程：

对于上式的超定方程A[u,v]T=-B，可以采用最小二乘法来求解 [u,v]T的最优解：

由于Lucas-Kanade光流法在计算过程中假设了相邻图像之间的位置不会发生剧烈的变化，因而在物体运动速度较快时，该假设不成立，从而导致后续的计算产生较大偏差。为了解决上述弊端，本文采用由Jean-Yves Bouguet提出的一种基于金字塔分层的改进Lucas-Kanade光流法[11]。该方法的主要原理是将图像的长宽每次缩放为原来的一半，通过对图像的缩放分层来起到减小像素位移的目的，同时金字塔这种结构也可以起到逐层分解、自上而下修正运动量的作用。

计算时先计算出图像金字塔的顶层（即分辨率最低的图片）的光流，将得到的结果代入相邻两层图像的递归公式中计算出下一层图像的猜测光流向量gL-1：

式中：dL表示像素点在第L层迭代运算中的剩余光流向量，通过标准的Lucas-Kanade光流算法迭代得到；gL表示第L层的猜测光流向量，并且假设顶层的g=[0,0]T。如此循环就能通过顶层的光流向量计算得到底层的猜测光流向量g0和剩余光流向量d0，从而得到最终的光流结果，即

式中，Lm表示金字塔的最大层数。

金字塔式的Lucas-Kanade光流法可以带来很高的像素位移增益，从而可以在不使用大窗口的情况下处理运动幅度较大的图片。

1.3 二元颗粒体系的倾斜角测量方法

当滚筒内的二元颗粒体系混合达到稳定状态后，利用高速CCD相机抓拍二元滚筒内的颗粒体系正视图，并将抓拍到的图像进行灰度化等预处理，随后再利用边缘轮廓提取以及曲线的拟合得到最终的表面曲线轮廓，从而进一步测量出此时堆积颗粒的倾斜角 θ[12-14]。图4为在 ω =0.031 rad/s时通过图像法所测得的颗粒体系倾斜角示意图。

图 4 滚筒颗粒体系倾斜角示意图Fig.4 Schematic diagram of the angle of slope in the particle system

1.4 颗粒混合程度的熵

针对颗粒物质随时间变化的混合或是分离程度如何量化这一问题，早期有学者们提出过一些不同的方法[15-24]。而在2010年时，Arntz等[6]提出了一种基于混合熵的方法，该方法有着统计学的理论基础，相较于先前的其他方法也很容易应用于不同类型的分离过程。于是借用此方法（即总体熵）来对不同转速下的颗粒分选程度进行量化。首先在二维滚筒的正视图上划分1 2 ×12个网格单元，用黑色直线标注出来，如图5所示。在定义网格之后，通过玻尔兹曼公式计算每个网格下的混合熵s(k)，计算公式如下：

式中，xa(k)和xb(k)分别指代在网格k中颗粒a所占的数量比以及颗粒b所占的数量比。在求解出每个网格下的混合熵后，通过加权求和的方式计算得出在某一时刻t下的总体熵S(t)，计算公式如下

式中：N表示整个二维滚筒内的总颗粒数；n(k)表示单个网格k中的颗粒个数。

根据Arntz等[6]的研究结果：在同一转速下的两种颗粒所组成的颗粒体系，它的总体熵最终会在一个值附近波动并趋于稳定。因此在本实验中，通过改变转速，选择在不同转速下颗粒体系达到稳定状态时的总体熵作为该转速下的颗粒分选程度的量化指标。

图 5 二维滚筒的正视图Fig. 5 Front view of the two-dimensional drum

2 结果与讨论

2.1 颗粒体系的总体熵

图6是不同转速下处于稳定状态的二维滚筒颗粒体系的总体熵的变化趋势，以及分别选取转速为0.023，0.065和0.102 rad/s下所拍摄的颗粒体系正视图（分别对应于图 6中的（a），（b）和（c）），图6（d）中横轴表示二维滚筒的转速，纵轴表示颗粒体系的总体熵。在0～0.051 rad/s的转速范围内，颗粒的运动状态表现为离散雪崩，颗粒体系的总体熵随着转速的增加逐渐减小。由图6（a）可以看到，颗粒体系的分选状态表现为花瓣模式（petals pattern），并且随着转速的增加逐渐转向月亮模式（moon pattern）。0.051～0.079 rad/s的转速范围内，颗粒处于由间歇性的雪崩向连续雪崩转变的过渡状态，颗粒体系的总体熵在定值0.614的附近波动。而在0.079～0.120 rad/s的转速范围内，颗粒的运动状态表现为连续雪崩，颗粒体系的总体熵随着转速的增加逐渐增大。值得注意的是，对于图6的（b）和（c），滚筒内的颗粒体系在这两个转速下所表现出的分选状态从拍摄的图像上很难看出差异，但从总体熵这个量化指标上看却存在着较为明显的差异，因此可以通过总体熵这个指标来进一步地描述颗粒分选过程中的细微变化。

2.2 表面流动层的速度特征

本实验在测量出混合总体熵的变化趋势后，进一步地对不同转速下的颗粒体系的表面流动层的速度进行研究，并最终希望将颗粒速度与之前所研究的总体熵相对应。

图 6 二维滚筒内颗粒体系总体熵的变化曲线Fig. 6 Variation curve of the total entropy of the particle system in the two-dimensional drum

实验分别选取了 0.058，0.065 和 0.072 rad/s 3 个转速，且在保证颗粒体系的分选过程达到稳定状态的情况下对本文所划分的10个区域内的速度分布进行测量研究，速度的测量方法遵循前文所介绍的Lucas-Kanade光流法，同时对于3个转速下的区域速度的空间分布进行横向对比。图7为在0.058，0.065和0.072 rad/s这3个转速下，区域Ⅱ内的速度分布图。从图7可以看出，当分选过程达到稳定状态时，颗粒体系的表面流动层上的速度呈现出正态分布，且标准差都在0.035左右。因而后续对速度的空间分布进行横向对比时，选择用每个区域内的平均速度作为该区域位置下的流动层速度。

图8为在3个不同的转速下，颗粒速度沿着表面流动层的分布情况，图中横轴表示各区域内的平均速度矢量距离流动层中心点的相对位置，纵轴表示表面流动层的平均速度。由图8可以看出，颗粒速度的空间分布在不同转速下呈现出近似相同的变化趋势，靠近坐标原点的区域由于颗粒经过加速的过程而导致速度较大，但总体的速度变化比较平稳；位于表面流动层两侧靠近边缘的区域由于颗粒在该处会发生碰撞等现象，颗粒之间的交换更多而导致速度较小，速度的变化较大；于此同时，3根曲线都在区域Ⅴ处产生了一个向下的突变，这与同种颗粒的速度空间分布[25]相比存在不同。通过对图8的分析可以认为在区域Ⅴ处发生速度上的突变之前，流动层的速度变化是比较理想的，并且靠近坐标原点附近的区域相比于两侧的区域，其颗粒速度更能反映在当前转速条件下的流动层颗粒速度的真实情况，因而在接下来的分析中选择突变之前的区域Ⅵ作为主要的测量区域。

图 7 区域Ⅱ内的速度分布直方图Fig. 7 Histogram of the velocity distribution in region II

图 8 表面流动层的速度空间分布随转速的变化曲线Fig. 8 Spatial distribution of the velocity on the surface flow layer at different rotating speeds

2.3 颗粒速度、倾斜角与总体熵的关系

图9是在不同转速条件下，颗粒运动的平均速度与堆积颗粒倾斜角的变化曲线，图中横轴表示二维滚筒的转速，左侧纵轴表示区域Ⅵ下的流动层平均速度，右侧纵轴表示颗粒体系的倾斜角。从图中可以看出，随着转速的增加，颗粒运动速度与倾斜角的变化趋势一致，呈现出一个先增大后平稳再增大的变化过程。雪崩过程中的颗粒运动速度与倾斜角呈正比例关系，这是因为在转筒内颗粒流中，当堆积颗粒的倾斜角增大时，颗粒床的重心位置也抬高，使得雪崩前颗粒的重力势能增加，导致崩塌过程中颗粒的运动速度变大。在0～0.051 rad/s范围内，转筒内颗粒的运动状态为离散雪崩，颗粒的运动速度与倾斜角随着转速的增加呈线性增大。而在0.079～0.120 rad/s的范围内，转筒内颗粒的运动状态为连续雪崩，此时颗粒的运动速度与倾斜角也随着转速的增加呈线性增大。值得注意的是，在0.051～0.079 rad/s的范围内，平均速度与倾斜角为定值，分别为0.085 m/s和32.85°。这可能是由于转筒内混合颗粒处于过渡运动状态（过渡运动状态是由间歇性的雪崩向连续雪崩转变的中间状态）导致的。

图 9 颗粒运动平均速度与倾斜角随转速的变化曲线Fig. 9 Average speed and angle of slope of the particle motion at different rotating speeds

对比图9与图6发现，转速大于0.051 rad/s时，颗粒的运动速度、倾斜角与颗粒体系的总体熵的变化趋势一致，在0.051～0.079 rad/s的范围内运动速度、倾斜角与总体熵三者均为定值，而在0.079～0.120 rad/s的范围内三者均随着转速的增加而增大。这说明了颗粒的运动速度和倾斜角影响着颗粒体系的总体熵。另一方面，也可以根据颗粒的运动速度或堆积颗粒的倾斜角去判断总体熵的变化，即颗粒体系的分选程度。

3 结 论

本文在滚筒转速为0～0.120 rad/s，滚筒内二维颗粒体系填充度为38%的条件下，采用Lucas-Kanade稀疏光流法与图像法观察均值粒径分别为2.6 mm和1.25 mm的规则球形颗粒所组成的颗粒系统的分选过程。通过测量颗粒体系的总体熵、颗粒体系表面流动层的平均速度与倾斜角，得到如下结论：

a. 针对二维滚筒内混合颗粒体系的分选程度，可以通过图像法与总体熵相结合的方法进行判别，当图像法观察不出明显的差异（如本实验过程中转速为 0.065 rad/s以及 0.102 rad/s时），可以通过测量颗粒体系的总体熵来区别分选过程中的细微变化。

b. 颗粒体系表面流动层的平均速度与倾斜角的正比例关系，说明了在崩塌过程中平均速度受倾斜角的影响。

c. 颗粒体系的总体熵、表面流动层平均速度和倾斜角与转速之间的关系，在转速大于0.051 rad/s时三者的变化趋势相一致，说明在混合颗粒的运动状态转为连续雪崩的条件下，颗粒体系的总体熵，即分选程度与平均速度和倾斜角有关，可以借助平均速度和倾斜角来判断总体熵的变化。