定量而精确的颗粒速度测量是获取颗粒流动特性、高效组织工质流动的前提，是研究与优化系统结构的设计基础，在能源、动力、化工与环境等领域具有广泛而重要的工程应用价值^［1−2］。

基于数字图像处理的颗粒检测方法是图像处理技术的重要分支，相比于其他测量方法图像法具有直观、可靠、简便等优点^［3−4］，尤其适用于稀疏颗粒相多参数在线测量^［5−6］。单帧长曝光轨迹图像法^［7−8］通过适当延长工业相机的曝光时间以获得离散颗粒的单帧“拖影”图像，即运动轨迹，通过轨迹宽度和长短同时获取颗粒粒径和运动速度信息，从而可以降低对测量系统帧率的要求。已经有许多学者对图像法运动场测量的方法及其应用做了大量研究，例如：吴学成等^［9］对煤粉颗粒进行了基于轨迹图像的速度和粒径测量；王若琳等^［10］采用单帧运动模糊图像处理测量了颗粒粒度信息；张超等^［11］利用轨迹图像实现了喷嘴液滴粒径和速度测量系统构建；冯明亮等^［12］从轨迹识别的角度，对单帧长曝光轨迹法测速上限的影响因素进行了理论分析和实验研究。但是，轨迹法的测量精度还有进一步提升的空间，尤其需要从图像处理方面做进一步的研究以获得更准确的结果。

本文基于颗粒运动中的成像过程，即轨迹成像法，首先，通过理论分析，提出基于轨迹图像灰度分布的颗粒测速方法，推导颗粒速度测量的理论公式；其次，基于仿真图像的处理分析，证明上述方法的可行性和可靠性；最后，通过实验研究，对基于该方法的颗粒速度测量误差进行分析，并与前期工作中采用的直接二值化结合Regionprops函数或Radon变换测速方法的误差进行对比研究，为利用背光照明式成像系统测量运动颗粒提供参考和依据。

1 基于轨迹图像灰度分布的颗粒测速原理 1.1 单帧长曝光成像原理

图1为典型的背光式单帧长曝光图像法测量系统示意图。光源直接朝向成像系统照明，物体（如颗粒）对光线形成遮挡从而在相机中投影成像，即形成背景亮物体暗的“剪影”图像。当在曝光时间t内颗粒平行于成像平面有运动时，像场的光能量（或能量的遮挡效应）于曝光时间t内在图像传感器上的累积，产生电荷信号，经过数模转换，最终由计算机量化为数字图像。因此，颗粒图像实际上是颗粒在曝光时间内经过的所有位置的像的叠加，如图2所示，运动模糊的图像正蕴含了物体的速度信息。

本文对处于运动状态下的颗粒进行如下三点假设：①运动颗粒为圆形，直径为D；②由于曝光时间很短，故可以假设颗粒在曝光时间内为近似的匀速直线运动，速度为v；③颗粒在垂直于成像平面方向上速度为0，即为二维运动。

基于上述假设，颗粒粒径等于轨迹图像的宽度；颗粒在曝光时间t内通过的位移L等于S和D的差，则颗粒在曝光时间内的平均速度为

$\qquad v = \frac{L}{t} = \frac{{S - D}}{t}$

(1)

常规速度测量的方法即通过图像二值化后识别轨迹的长轴长度和短轴长度，已知相机的曝光时间，基于式（1）获得速度。获取S和D的图像处理算法主要有两种：①Matlab软件自带的图像处理工具箱中的Regionprops函数；②图像在一个特定角度下的径向线方向投影得到轨迹图像的Radon变换函数。Regionprops函数是Matlab软件中一个重要的图像分析函数，它可以快速度量图像区域的像素数、重心、等效圆直径等参数^［13］。学者们利用函数中的等效短轴和等效长轴这两个参数分别来近似计算颗粒轨迹图像的宽度和长度，从而进一步提取出颗粒的粒径和速度信息。Radon变换函数由奥地利数学家约翰·雷登于1917年提出，其数学意义是平面内函数f（x，y）沿特定直线的线积分，在二维图像中，一幅图像的Radon变换结果就是这幅图像在一个特定角度下的径向线方向的投影^［14］。但由成像原理导致轨迹边缘灰度与背景灰度接近，且存在系统噪声，难以合理选取阈值大小，二值化操作常常带来较大误差。因此，本文拟从图像灰度分布的角度，探索颗粒速度测量的方法，避免或减小二值化过程带来的误差。

1.2 轨迹图像灰度分布

根据运动轨迹成像原理，对圆点轨迹图像的灰度分布进行分析。假设有一不透光圆点半径为R，以速度v由左向右水平运动，成像系统放大倍率为M，图像传感器的像元大小为P，圆点在曝光时间t内形成的理想轨迹如图3（a）所示。假设运动轨迹长度大于直径（实验时可通过调节曝光时间以满足该条件），可将其分成Ⅰ和Ⅱ两个区域分别进行研究，其中区域Ⅰ为轨迹的两端[见图3（b）]，区域Ⅱ为轨迹的中间部分[见图3（c）]。

由于运动轨迹具有对称性，以轨迹几何中心为原点，以颗粒运动方向为x轴建立像素坐标系，取轨迹的左上部分区域进行分析。以区域Ⅰ为例，在背光成像方式下，区域Ⅰ中坐标点（x，y）在成像过程中被圆点遮光的时间t_z为

$\qquad {t_{\textit{z}}} = \dfrac{L}{v} = \dfrac{{\sqrt {{{(MR)}^2} - {{(Py)}^2}} + \dfrac{{vtM}}{2} - Px}}{v}$

(2)

式中，L为该像素被圆点遮光时间下，圆点运动的轨迹长度。由于像素点相对灰度（该像素点灰度与背景灰度的差值） $\Delta G_{\rm{I}}$ 与该点被遮光的时间t_z成线性关系，即

$ \Delta {G_{\rm I}} = k{t_{\textit{z}}} = k\left( {\frac{{\sqrt {{{(MR)}^2} - {{(Py)}^2}} - Px}}{v} + \frac{{Mt}}{2}} \right)$

(3)

式中，k为系统常数，与光源强度和系统响应有关。

同理推导得到区域Ⅱ中的图像灰度值分布为

$\qquad \Delta {G_{\rm II}} = k\frac{{2\sqrt {{{(MR)}^2} - {{(Py)}^2}} }}{v}$

(4)

可见，区域Ⅰ的灰度在同一y值下与x成线性关系，且线性系数为v的函数；而区域Ⅱ的灰度与x无关，只与y有关。本文基于区域Ⅰ的灰度分布进行速度测量。当y = 0时，即轨迹区域Ⅰ长轴线上的像素点灰度分布为

$\qquad \Delta {G_{{\rm I}0}} = k\left( {\frac{{MR - Px}}{v} + \frac{{Mt}}{2}} \right)$

(5)

拟合曲线 $\Delta {G_{{\rm I}0}} {\text{～}} x$ ，得到斜率 $ - \dfrac{{Pk}}{v}$ 。已知像元大小P，通过系统标定得到系数k，即可获得速度信息v。可见，该图像处理过程与圆点大小R和曝光时间t均无关。关于系数k的标定，则可以在已知某一圆点大小R和曝光时间t的条件下，基于式（5）进行 $\Delta {G_{{\rm I}0}} {\text{～}} \left( {\dfrac{{MR - Px}}{v} + \dfrac{{Mt}}{2}} \right)$ 的标定获得，后文中将对实验标定结果作详细阐述。

1.3 图像处理流程

针对基于灰度分布拟合的圆点速度测量过程，对圆点运动轨迹图像的信息提取算法进行研究，形成的图像处理流程如图4所示。

具体步骤为：

（1）图像预处理：读取圆点轨迹图像I_P[如图5（a）所示]和背景图像I_b，对I_b−I_P的图像经过维纳去噪后得到图像I，如图5（b）所示；

（2）采用Otsu算法计算整幅图像阈值，并进行图像二值化，二值化图像边界位置如图5（b）中黑色闭合曲线所示，若一幅图像中有多条轨迹则进行轨迹图的分割；

（3）针对单个轨迹图像，利用Radon算法确定轨迹中心点并建立坐标系；

（4）如图5（b）所示，提取颗粒轨迹的长轴长度L₁和短轴长度L₂，以长轴一端点 ${x_{\rm{A}}}$ 作为选取像素点的起始位置（ ${x_{\rm{A}}} = {L_1}/2$ ）；从该点起向坐标中心方向共选取 ${L_2}/2$ 个像素位置，作为待拟合区域，即 $\left( {{x_{\rm{B}}},{x_{\rm{A}}}} \right)$ ；并读取图像I中该区域的灰度值 $\Delta G$ （由于第一步预处理，此处已为相对灰度值）。

（5）拟合 $\Delta G {\text{～}} x$ ，得到斜率 $ - Pk/v$ ，代入图像传感器的像元大小P和标定系数k，即可计算出颗粒速度v。

2 仿真图像分析

为了对基于轨迹图像灰度分布的测速方法可行性进行验证，仿真生成不同圆点半径和不同速度下的轨迹图片，按照1.3节的图像处理流程进行处理获得速度，并与理论值进行比较，分析噪声的影响。

基于Matlab软件平台，采用Motion算子卷积清晰圆点图片形成运动模糊轨迹。以图6（a）为例，它是半径为81个像素的圆点水平运动301个像素时形成的图片，其中运动模糊前的图片背景灰度为0，颗粒静止成像的前景灰度为0.8；添加高斯噪声（噪声方差3.25）后如图6（b）所示。假设成像系统曝光时间为20 μs，放大倍率为1倍，像元大小为5.3 μm，根据图6（a）获得的标定系数k为1.026 8×10⁷。

假设圆点半径由11变化到101个像素，运动模糊轨迹长度由101变化至1 001个像素，获得不同圆点半径和不同速度下标定的k值，如图7（a）所示。可见：该系数与所选圆点大小及运动速度基本无关；在圆点较小时，由于拟合数据点数量过少（约为圆点半径对应的像素个数），标定的k值变化范围较大；在圆点较大时，若运动长度接近或小于圆点直径，由于所拟合区域将包含图3中的区域Ⅱ， $\Delta {G_{{\rm I}0}} {\text{～}} x$ 数据将偏离线性关系，上述方法将不再适用，即该方法的前提是运动轨迹足够长（这一点在实际测量中可以通过调节曝光时间来保证）。取标定系数的平均k值（1.030 6×10⁷），根据1.3节所述的处理流程，得到速度计算值，并与理论值进行比较，结果如图7（b）和（c）所示。图7（b）为不含噪声图像的测速结果误差，其值等于各k值与平均k值的相对误差；随着圆点直径的增大，测量相对误差有整体减小的趋势，这一前提是运动轨迹长度大于圆点直径；但随着轨迹长度的进一步增加，由于运动模糊导致前景图像灰度逐渐接近背景，即所拟合曲线的斜率逐渐降低，误差又逐渐增大，即测速范围存在上限，而该上限与粒径有关。图7（c）为仿真图片添加方差为3.25的高斯噪声后的测速结果。在轨迹长度较长时，拟合曲线的灰度变化可能仅有几个灰度级，噪声造成的影响非常大，极个别到100%以上。为更好地显示可测量区域的误差，此处仅显示501个运动长度以下的测速误差。可见，误差水平整体提高，但圆点直径和轨迹长度的影响趋势与图7（b）类似。

3 实验装置及参数标定 3.1 背光式图像法圆点测速实验装置

本文搭建了一套背光式图像法圆点速度测量系统，如图8所示。该系统包括工业相机（FL3−U3−13Y3M）、远心镜头（2 ×）、直流电机、光度计、转速仪和平行光源等。拍摄对象为固定在铝板扇叶上标定板中的黑色圆点。实验时标准圆点板随着铝板扇叶的转动获得转速，并被调至镜头焦平面，以保证相机镜头轴线、光源轴线及标定圆点板所拍摄区域中心位于同一水平线上，且与标定板平面垂直，避免远离光轴位置的离焦模糊现象。相机采集图片位深选为10位。值得注意的是，测量的轨迹长度一般为圆点直径的3～8倍，在此区间内运动轨迹的平均曲率仅为0.003 5，因此，圆点可认为是沿着水平方向运动，此时适用1.2节中的运动轨迹成像原理。

3.2 参数k的标定

基于背光式图像法实验装置，首先对该实验系统进行标定，实验选取直径分别为80、100、150、200、250和400 μm共6种圆点，每种圆点运动速度分别取5.0、6.3、7.6、8.9、10.1、11.4和12.6 m·s⁻¹，拍摄并筛选出含有运动轨迹的517张实验图片进行处理，获取相应的x与 $\Delta {G_x}$ 。以 $\left( {\dfrac{{MR - Px}}{v} + \dfrac{{MT}}{2}} \right)$ 为横坐标，以 $\Delta {G_x}$ 为纵坐标绘制散点图并拟合出斜率k，结果如图9（a）所示。该图中整体趋势与仿真结果一致。当圆点直径较小时，用于拟合的数据较少，且在高速条件下相对灰度值变化较小[如图9（b）所示]，拟合误差增加；即测速范围上限随着粒径的降低而降低，如本文实验条件下，直径100 μm圆点的运动速度上限约为8.9 m·s⁻¹。粒径较大或速度较低时，拟合曲线如图9（c）所示，符合基于灰度分布拟合速度测量方法的测速范围。

上述仿真和实验结果均表明：本文提出的基于轨迹图像灰度分布的测速方法更适合使用在圆点直径较大的速度场检测中；测速范围存在上、下限，下限要保证运动长度大于颗粒直径，上限需保证拟合区域的图片灰度有较明显的变化；还需要深入研究确定上、下限的方法；此外，噪声的存在明显影响图像的质量，进而影响测速结果的准确性。图像去噪方法值得进一步深入研究。

3.3 速度测量结果与误差分析

根据图像处理算法原理结合实验图像，采用灰度分布拟合法、二值化结合Regionprops 函数法以及二值化结合Radon变换函数法等三种测速方法进行图像处理并获得速度测量误差，结果如图10所示。

由图10中可知，当圆点直径较小时，三种图像处理方法的相对误差相差无几，但随着圆点直径增加，二值化结合Regionprops函数法和Radon变换函数法相对误差大大增加，而灰度分布拟合法的测速相对误差则无明显增加，甚至有所降低。究其原因，如图11所示，黑色区域为二值化之后的轨迹图像（Radaon函数处理基于该区域），椭圆为与二值化区域具有相同标准二阶中心矩的椭圆（Regionprops函数处理基于该区域）；利用Radon函数法处理运动轨迹时，其径向线方向的投影无法解决二值化过程中的长轴损失问题，所以，其对测速结果带来的误差最大；而Regionprops函数法测速结果虽优于Radon变换函数法，究其原因是由于其等效椭圆长轴拉伸了部分条件下的二值化轨迹图像，但是不同情况下长轴损失的程度各不相同。相比较而言，本文提出的利用轨迹成像原理速度法由于避免了直接基于二值化图进行测量，测速误差可降低5%～25%。

4 结　论

本文从图像灰度分布识别角度，对单帧长曝光轨迹法测速误差进行理论和实验研究，提出一种基于轨迹图像灰度分布的测速方法。具体结论为：

（1）单帧长曝光圆点轨迹的长轴线两端区域（区域I），像素点相对灰度与该点距轨迹质心的距离在理论上呈较为严格的线性关系，且标定系数k与圆点大小R和曝光时间t均无关；

（2）在对仿真图像的分析中发现，圆点直径越大，本文提出的基于轨迹图像灰度分布的测速方法计算结果越准确；在测速范围内轨迹长度即颗粒速度对其影响不大；图像噪声使得误差整体增大，但不影响上述规律；

（3）搭建了圆点测速实验装置并通过实验比较了三种测速方法；在实验条件下，对比Otsu二值化结合Regionprops函数法以及二值化结合Radon函数法。本文提出的基于灰度分布拟合速度法测速误差可降低5%～25%，可靠性更高。