引　言

颗粒物质是由大量固体颗粒堆积而成的离散系统，广泛存在于自然界与工业生产中^[1-5]。颗粒物质的崩塌是其中一种常见的运动模式，常见于雪崩、泥石流等自然灾害中^[6-7]。颗粒崩塌可以划分为表面崩塌颗粒流与内部静态堆积颗粒。内部颗粒总体上可以认为是相对静止的，但其中会有少数颗粒发生局部、微小的非仿射位移，导致局部堆积结构的变化，这种局部的结构上的变化也被称为颗粒重排^[8-9]。内部颗粒的局部堆积结构变化很大程度上影响着整个系统的动力学行为，研究其在崩塌前后的分布、变化特征可以更好地对崩塌行为进行预测。

目前研究人员在已知颗粒具体位置的模拟系统中已经提出了许多描述堆积结构的方法，例如传统方法的Voronoi图用来对堆积颗粒的体积进行划分，新提出的Cell网格方法将堆积颗粒按照接触颗粒的中心点沿着逆时针方向连接起来围成的区域进行划分，保留了接触颗粒的连通性信息，在颗粒统计力学中发挥着重要作用^[10-12]，这些方法的有效性已被广泛证明。但在实验系统中，如何获得每个颗粒的位置依然是一个难题。图像法在二维圆盘的堆积结构测量中已被广泛使用，华云松等^[13]通过连通域特征提取二维圆盘颗粒的空间坐标，并通过Cell量化了其在二维筒仓中的几何堆积结构。Seguin等^[14]使用颗粒追踪以及光弹性技术考察了二维压缩系统中接触力网络的变化。但对三维紧密堆积的颗粒系统来说，获得其堆积结构等颗粒微观特性的实验手段是有限的。Brodu等^[15]使用激光片光照射浸没在折射率匹配的浴液中的浸透荧光染料的水凝胶颗粒，采集颗粒系统横截面图像，通过一组横截面图像叠加重建出颗粒的表面。但该方法对颗粒材质以及颗粒系统限制较大且时间分辨率较小，难以捕捉到堆积结构的变化。Born等^[16]使用散射方法，用太赫兹电磁辐射在颗粒介质中的传输特性量化颗粒堆积结构，但该方法不够直观且还处于理论阶段，距离实际使用还有一段距离。

在三维颗粒系统成像困难的情况下，可以考虑用表面颗粒的二维信息表征三维系统的特征。如在颗粒流测速领域，研究人员使用粒子图像测速（particle image velocimetry，PIV）方法测量表面颗粒的二维流场^[17]，从而一定程度上反映了三维颗粒流的运动特征。三维滚筒中的颗粒系统沿滚筒轴向大致是均匀分布的，滚筒侧面颗粒的堆积结构一定程度上也可以反映整个系统的堆积结构。但之前在PIV中对颗粒表面图像的处理不涉及具体颗粒的识别，如何准确识别表面单层颗粒，滤除堆积孔隙中露出的内部颗粒的干扰，还需要进一步优化。

本文使用Voronoi图研究了三维滚筒中的颗粒体系在崩塌周期中侧面颗粒的堆积结构。首先使用图像法对滚筒侧面表层颗粒进行识别，使用具有颗粒反光形状特征的特征模板去除第二层颗粒的影响；再使用Voronoi图来表征滚筒侧面表层颗粒的结构特征，以此来反映三维滚筒内部颗粒的结构；然后研究崩塌周期中重排区域堆积结构的变化，并通过图像差分的方法估算出崩塌规模，分析其与堆积结构之间的相关性。

1 实　验

旋转滚筒内的颗粒崩塌相对容易控制与测量，并且可以持续产生周期性的崩塌流，其已成为研究颗粒崩塌流动力学的常用模型系统^[18-19]。实验中使用低速旋转的三维滚筒为后续研究提供周期性的崩塌流，如图1（a）所示。滚筒中颗粒崩塌系统可以划分为表面崩塌流（主动层）以及内部静态堆积的区域（被动层），如图1（b）所示。颗粒速度梯度为0的区域为主被动层交界处，根据之前对该实验系统的测量，距离自由表面垂直深度大于10 mm的区域为被动层^[20]，此处自由表面的角度取上、下休止角的均值。

1.1 实验装置

实验装置由高速相机、LED光源、三维滚筒与直流电机组成，如图1（a）所示。滚筒由透明有机玻璃制成，直径140 mm，长20 mm。两根金属支撑杆水平固定且相互平行，均可以旋转，两者之间通过一根传动皮带连接。旋转电机采用直流电源供电，控制其中一根支撑杆旋转，另一根支撑杆在传动带作用下以相同的角速度旋转。支撑杆上固定4个圆柱形橡胶块，滚筒放置于橡胶块上，轴线与支撑杆平行，通过橡胶块的摩擦力带动旋转。改变直流电源输出电压可以控制滚筒转速，滚筒转速设置为0.19 r/min。滚筒内填充粒径为0.95 ～1.05 mm的干燥、无黏性电镀玻璃珠，平均粒径为1 mm。滚筒装置放置在水平实验平台上，采用高速相机沿着滚筒轴向采集滚筒侧面图像，尽可能保证相机与滚筒中心处于一条水平线上。

本实验中使用的面阵高速相机为千眼狼5F04M型黑白摄像机，该相机基于CMOS图像传感器，最大内存为16 GB。在实验过程中选取了两种帧率。采集整个滚筒侧面图像时，目标主要为表面崩塌颗粒流，其运动速度较大，不需要太高的帧率。同时为了消除电机振动、转速不稳定带来的误差，相机以5000 帧/s的帧率采样，每帧图像分辨率为1984×1680像素，每10帧图像作平均，最终得到的图像序列为500 帧/s。在采集内部颗粒重排事件图像时，颗粒重排的持续时间大概在1 ms数量级，故为进一步捕捉颗粒重排过程，相机采用10000 帧/s的帧率进行采样，图像分辨率为256×1396像素。

这两种帧率下采集得到的单帧图像并没有明显差异，区别主要体现在相邻两帧图像之间的差异上：高帧率下相邻两帧图像差异小，便于观察重排过程中颗粒详细的运动轨迹；低帧率下相邻两帧图像之间差异大，直接展示出重排颗粒的起始位置与终止位置，省略了重排的具体过程，便于重排事件的统计，并且低帧率下可以采集更多的连续的崩塌周期，便于统计具有相同历史堆积结构的重排事件。在后文中，2.1节对重排细节的讨论中使用的是高帧率下采集的图像，2.2节对重排与崩塌之间相关性的讨论中使用的是低帧率下采集的图像。

同时在采集过程中用LED灯照射滚筒侧面进行补光，以提高采集图像的对比度。实验开始时，滚筒内颗粒自由堆积，自由表面处于水平状态；在滚筒匀速旋转5 min后认为系统进入稳定的崩塌周期，开始图像采集。图1（b）给出了其中一帧典型的实验图像。

1.2 颗粒堆积图像处理方法

图2给出了5000 frame/s帧率下获得的图像序列平均方法示意图。原始图像每10帧分为一组，针对一组内的10张图像的每个像素点的灰度值取均值，得到500 frame/s的图像序列。平均后的图像序列与10000 frame/s帧率下采集的图像序列一起进行图像畸变校正。

图3（a）给出了校正前的图像。由于相机或实验装置摆放位置有偏差，图像存在一定程度的畸变。选取比例尺上、下沿与滚筒圆周的交点A、B、C、D 4个点作为关键点进行透视变换。已知滚筒直径为140 mm，根据AB、CD的弦长即可计算出4点的真实坐标。透视变换的结果如图3（b）所示。根据图中比例尺可知，一个像素对应的实际长度约为0.667 mm。

滚筒侧面颗粒的堆积结构一定程度上可以反映整个系统的结构特征。测量其堆积结构首先需要对每个颗粒的位置进行识别。在滚筒侧面图像中选取被动层边长为17 mm的正方形区域为例，颗粒识别过程如图4所示。

颗粒表面在LED灯光照射下的反光具有丰富的灰度信息，表层颗粒与第二层颗粒的反光具有显著差异。取10个典型表层颗粒的反光图像的平均值作为颗粒反光的特征模板，如图4（b）所示。用该模板遍历原始图像上的每一个位置，计算原始图像和特征模板的互相关系数。记原始图像为I，特征模板为B，则互相关系数可以表示为

$ F(k,l)=\sum\limits_{p,q=1}^MI(k+p,l+q)B(p,q) $

(1)

式中：M为特征模板的边长；k、l为模板在原始图像上的位置。互相关计算结果如图4（c）所示。

互相关计算结果依然保持了圆形特征，故采用霍夫变换对互相关计算结果的灰度图像进行识别。霍夫变换是将x-y平面上的点映射到参数空间的变换，参数空间是由待识别目标的形状决定的。对于圆形目标而言，其在x-y平面上的方程为

$ {r^2} = {(x - a)^2} + {(y - b)^2} $

(2)

可以看出圆有3个参数r，a和b，其中a和b分别为x和y方向上的圆心位置，r为圆的半径。圆的参数表示为

$ \left\{\begin{gathered}x=a+r\cos\alpha \\ y=b+r\sin\alpha \\ \end{gathered}\right. $

(3)

由式（3）可知，圆的参数空间为R³，即x-y平面同一个圆形上所有的点可以经过霍夫变换映射到参数空间中的一个点（a, b, r）。霍夫变换检测圆的基本算法步骤为

1）对原始图像进行Canny算子边缘检测，得到边缘检测后的二值图像；

2）创建一个大小与参数空间一致的三维累加器矩阵A[a, b, r]；

3）对于图像中的每个边缘点（x_i, y_i），画出所有通过该点并且r在指定范围内的圆，并映射到参数空间即可得到一个曲面，曲面上所有点在累加器矩阵中对应位置的值加一；

4）在累加器矩阵A中设置一个阈值搜索局部最大值，获得图像中圆的中心位置及半径。

霍夫检测结果如图4（d）所示，红色圆形即为检测到的圆形亮点。

最终根据圆心位置得到的所有颗粒中心的坐标如图4（e）所示。该方法一方面去除了第二层颗粒的反光对表层颗粒识别的干扰[图4（d）中的红色方框]，另一方面与直接阈值分割的方法相比，修正了反光点与颗粒中心位置的偏移，有效识别出表层颗粒的位置。

完成每个颗粒位置的识别后采用Voronoi细胞划分方法对其堆积结构进行表征。给定一组离散的颗粒，每个颗粒的Voronoi细胞是比任何其他颗粒更接近它的空间区域。对于颗粒s_i、s_j，作其连线的中垂线将平面划分为两个半空间，其中一个半空间中所有的点距离s_i更近。根据Voronoi细胞的定义，s_i的Voronoi细胞由所有最靠近它的点组成，所有距离s_i比s_j更近的点组成的半空间H_ij可以表示为

$ {H_{ij}} = \{ x \in M|d(x,{s_i}) \leqslant d(x,{s_j})\} $

(4)

式中，d表示空间中点x到颗粒s_i /s_j的距离。将s_i的Voronoi单元定义为所有这些半空间的交集，表示为

$ V({s_i}) = \mathop \cap \limits_j {H_{ij}} $

(5)

式中，j表示系统中所有其他颗粒。而远离中心颗粒的颗粒不影响其Voronoi细胞，Voronoi细胞的形状通常仅由附近的颗粒决定。由式（4）和式（5）提供的定义迭代地使用相交半空间切割图4中颗粒s_i周围的区域即可得到整个系统的Voronoi图，如图5所示。

图6为被动层区域Voronoi多边形面积与边数的分布，其中S表示单个Voronoi细胞的面积，Ps表示面积S的概率密度函数，n表示单个Voronoi细胞多边形的边数，N表示每种边数Voronoi细胞的个数。Voronoi多边形面积服从伽马分布，图中红色曲线为式（6）预测的伽马分布，该式是由Aste等^[21]提出的三维随机堆积颗粒系统中Voronoi体积的分布，其中V为每个颗粒Voronoi多面体的体积，V_min为整个颗粒系统内Voronoi多面体的最小体积，$ \overline V$为平均体积。该式是变量(V－V_min)的伽马分布，其特征为形状参数k和尺度参数$ \overline V $－V_min，此处用面积S代替式(6)中的体积V，实验测量结果与理论预测具有良好的一致性，验证了侧面单层颗粒的堆积结构信息可以反映三维系统的堆积结构特征。边数服从正态分布，平均边数为6，最大边数为8，最小边数为4。

$ p(V)=f(V,k)=\frac{k^k}{\Gamma(k)}\frac{(V-V_{\min})^{k-1}}{(\overline{V}-V_{\min})^k}\exp(-k\frac{V-V_{\min}}{\overline{V}-V_{\min}}) $

(6)

2 结果与讨论 2.1 颗粒重排对堆积结构的影响

实验结果显示，在崩塌周期中，滚筒中被动层区域静态堆积的颗粒并不是完全相对静止的，局部颗粒会发生一些非仿射位移，导致局部结构的变化。这种局部堆积结构的变化也被称为颗粒重排。重排颗粒的非仿射位移非常微小，从原始图像中难以直接观察到其位置变化，故以1.2节中识别出的颗粒坐标为基础，认为在0.1 ms内与其相邻颗粒相对位移超过粒径20%的颗粒参与了颗粒重排，相邻的重排颗粒组成了一个单次重排事件。滚筒被动层区域的颗粒系统具有许多不同的堆积结构，这些堆积结构具有不同的稳定性，即有些结构在外力作用下相对更容易发生重排。为了考察局部重排区域的结构特征以及重排前后非仿射位移带来的结构的变化，本节研究了不同重排区域的结构特征。首先对单次重排规模进行统计。图7给出了每次重排过程中参与重排颗粒个数的统计结果。

图7（a）给出了测量区域示意图。Li等^[20]的研究显示颗粒重排事件的分布在崩塌系统上下游之间存在明显区别，但其在被动层中的传播没有固定的方向，故测量区域的划分参考地面坐标。选取崩塌开始前自由表面最高点的水平位置作为左边界，如图7（a）中A点所示；崩塌结束时自由表面最低点的水平位置作为右边界，如图7（b）中B点所示。等间隔将其划分为七个区域，得到每个区域水平方向上宽度为17 mm，竖直方向上包含该宽度内所有颗粒。统计115个崩塌周期每个区域被动层中每次颗粒重排涉及的颗粒个数，持续时间约480 s。由于区域六、七的被动层中颗粒在表面颗粒撞击作用下堆积相对较为紧密，基本不发生颗粒重排，统计结果显示115个周期中区域六只发生了2次重排，区域七没有发生重排。可见，区域六、七的结果具有偶然性，故后续的讨论中不考虑这两个区域。区域二至区域四的被动层面积根据区域一的面积进行归一化，得到归一化后每个区域中参与重排颗粒个数分布，如图7（c）所示。

由于只有1～2个颗粒参与的颗粒重排对其邻域颗粒的堆积结构不会造成影响，可以认为其是在“颗粒牢笼”中相对自由的颗粒，不参与后续重排区域结构变化的讨论，因此统计结果从3个颗粒参与的颗粒重排开始展示。结果显示，重排涉及的范围大部分为3～6个颗粒，最大范围为10个颗粒，且出现频率大致随参与重排颗粒数量增加而减小；比较不同区域的结果可以看出，区域一内重排发生的次数与规模都是最大的，且区域一至区域五重排发生的次数与规模逐渐减小。

在参与重排颗粒个数分布结果的基础上可以对二维颗粒系统中参与重排颗粒的堆积结构进行进一步考察。对重排区域内颗粒进行二维Voronoi划分，根据颗粒的运动趋势以及Voronoi图结构可以将重排颗粒划分为三种基本团簇：线形、三角形以及四边形，认为具有相近方向位移的颗粒处于同一颗粒团簇。图8为重排前后相邻两帧图像中三种基本颗粒团簇对应的颗粒位置及Voronoi图变化。蓝色点为重排前颗粒位置，红色点为重排后颗粒位置，团簇内颗粒具有相近的运动趋势，可以认为在重排过程中团簇内颗粒作为一个整体一起运动。

由图7（b）中给出的参与重排的颗粒个数分布可知，除了只包含单个基本团簇的小范围重排（3～4个颗粒）之外，还存在许多较大范围的颗粒重排（5～10个颗粒）。对于这些较大范围的重排，可以根据不同的运动趋势将其划分为几个基本团簇。图9给出了一次由8个颗粒参与的颗粒重排过程。

如图9所示，蓝色点为重排前颗粒位置，红色点为重排后颗粒位置。将该颗粒重排区域划分为Ⅰ、Ⅱ和Ⅲ3个基本团簇，其中Ⅰ、Ⅱ为三角形，Ⅲ为四边形。区域Ⅰ作为一个整体向右下方运动，区域Ⅲ作为一个整体向左下方运动，区域Ⅱ作为区域Ⅰ、Ⅲ之间的桥梁存在，将区域Ⅰ指向右下方的压力传递到区域Ⅲ中。类似于图9中给出的划分方法，对所有颗粒重排进行基本团簇的划分，统计在115个崩塌周期被动层颗粒重排中3种基本团簇出现的次数，如表1所示。

根据表1给出的结果可知，3种基本团簇中线性团簇出现次数最多，三角形团簇出现次数最少。即线性团簇结构相对最不稳定，更容易出现在重排区域中；而三角形团簇结构相对较为稳定，较少参与颗粒重排。

对区域一中线形团簇指向的方向进行统计，由于该系统中线形团簇在重力作用下均向下运动，取线形团簇中最上端颗粒为起始点，最下端颗粒为终点，团簇的方向向量记为n，如图10中插图所示。计算方向向量n与z轴负方向的夹角θ，夹角位于z轴右侧记为正值，位于左侧记为负值。θ的分布情况如图10所示。由图10可知，大部分θ分布在0°～60°之间，即大部分线形团簇的方向指向右下方向。这说明线形团簇在右下方向上具有更高的自由度，而在接近于指向圆心的方向不易发生结构变化。

2.2 崩塌与堆积结构的相关性分析

颗粒系统内部堆积结构的变化在崩塌的发展过程中起到十分重要的作用。本小节主要分析了崩塌间隔阶段被动层堆积结构变化的剧烈程度与下一次崩塌之间的相关性。一个崩塌周期中参与崩塌颗粒的最大数量可以用来衡量该次崩塌的规模。图11给出了崩塌颗粒数量的估计过程。图11（a）、11（b）为崩塌过程中相邻两帧图像，间隔时间为2 ms。在滚筒圆周上选取5个点，根据其坐标使用最小二乘法拟合出圆心O的位置。已知滚筒以0.19 r/min的速度逆时针旋转，为消除滚筒旋转带来的颗粒位移的干扰，将图11（b）以O为旋转中心顺时针旋转3.98×10⁻⁵ rad，旋转前后的图像中心点保持不变，截取与图11(a)相同大小的部分，两者作差，差分结果的伪彩色图如图11(c)所示。

统计了100张崩塌间隔期间相对静止的颗粒系统的差分图像灰度值分布，结果显示灰度值最大为0.17。认为灰度值大于0.2的区域的颗粒堆积结构发生了变化。对差分图像进行二值化处理，如图12（a）所示；再以直径为20像素的圆形结构元素对其进行闭运算，以填补二值化图像中崩塌区域的空白，得到崩塌区域如图12（b）中白色区域所示。该区域在崩塌周期中的最大深度约为9.87 mm，与主动层最大深度基本一致，说明该阈值选择比较合理。估算参与崩塌的颗粒数量为

$ N=0.667^2\times\frac{cL\varphi}{v} $

(7)

式中：c为崩塌区域像素数；L为滚筒长度；v为单个颗粒的体积；$ \varphi $为颗粒体系的体积分数，根据之前的数值模拟结果可知约为0.58。

被动层中颗粒堆积结构变化的剧烈程度可以用发生重排的颗粒数量来衡量。重排颗粒的统计与崩塌颗粒采用相同的阈值，即认为差分图像中灰度值大于0.2的被动层颗粒发生了重排。统计22个连续的崩塌周期的最大崩塌颗粒数量以及每个崩塌阶段之前一个崩塌间隔中堆积结构发生变化的颗粒的个数，如图13（a）所示。其中，横坐标N_arr为崩塌间隔期间发生重排的颗粒个数，纵坐标N_ava为每个崩塌间隔之后紧跟的一次崩塌中崩塌颗粒的数量。可以看出两者之间具有较强的线性相关关系，计算得到相关系数R=0.8968，拟合出线性相关斜率k=20.58。不难看出，崩塌规模与该次崩塌前颗粒系统内重排颗粒个数呈正线性相关，即崩塌前静止抬升阶段重排颗粒数量越多，参与崩塌的颗粒数量越多，崩塌规模越大。崩塌颗粒个数与重排颗粒个数之间大致服从N_ava=20.58×N_arr+1428，线性拟合结果如图13(a)中直线所示。

认为差分图像中刚出现灰度值大于0.2的区域时为崩塌开始时刻，没有灰度值大于0.2的区域时为崩塌终止时刻。图13（b）给出了崩塌终止后0.1 s内发生重排的颗粒个数与该崩塌终止时刻到下个崩塌开始时刻之间的时间间隔之间的关系。纵坐标t为崩塌终止时刻到下一次崩塌开始时刻之间的时间间隔。可以看出，两者具有较强的负相关性，计算得到相关系数R=−0.8875，拟合出线性相关斜率k=−1.982。颗粒系统在静止抬升阶段发生重排的次数的增加会使得下一次崩塌的发生提前，两次崩塌之间的时间间隔t与重排颗粒个数N_arr之间的线性负相关大致服从t=−1.982×N_arr+275.9，线性拟合结果如图13(b)中直线所示。这说明颗粒重排带来的微小扰动会使还未达到上休止角的颗粒系统提前发生崩塌，重排颗粒数量越多，扰动越大，崩塌的发生越提前。

由图13的结果可知，崩塌间隔期间被动层区域颗粒堆积结构变化越剧烈，接下来的一次崩塌规模会越大，并且下一次崩塌的发生会提前。据此，可以通过崩塌前颗粒系统内部堆积结构变化的信息对崩塌行为进行一定程度的预测，为雪崩、泥石流等自然灾害的防制提供了更充足的参考依据。

3 结　论

研究了低速旋转滚筒中周期性颗粒崩塌系统中被动层的堆积结构，使用颗粒反光的特征信息有效去除了内部颗粒反光对表层颗粒识别的干扰，准确识别出滚筒侧面单层颗粒的位置；研究了滚筒侧面被动层区域颗粒的堆积结构的分布特征及其在崩塌周期中的变化；使用图像差分方法估算了崩塌规模，分析了堆积结构变化与崩塌之间的相关性。

1）被动层中局部的堆积结构变化可以划分为线形、三角形、四边形3种基本团簇，其中线形团簇最不稳定，更容易参与到局部重排中，而三角形团簇相对最稳定，更不容易参与到重排中。其中线形团簇在垂直于半径的方向上具有更高的自由度，而在接近于指向圆心的方向不易发生结构变化。

2）实验结果揭示了滚筒中颗粒被动层靠上游的位置结构相对更脆弱，更容易发生重排。这可能是由于颗粒在主动层下游区域进入被动层时具有较强的各向异性，可以支撑指向圆心方向的压缩而难以支撑垂直于半径方向的剪切作用。滚筒内颗粒在崩塌间隔期间重排次数越多意味着主力链受到破坏的程度越大，导致被动层堆积结构难以再为主动层颗粒提供足够的支撑，使得下一个崩塌的规模扩大，崩塌开始时刻提前。这一结论为人们更好地预测与防治雪崩、泥石流、山体滑坡等自然灾害提供了进一步的理论支撑。