基于数字图像处理的颗粒流厚度动态提取方法研究

吴 越，李 坤，程谦恭,2,3，王玉峰，龙艳梅，姜润昱，宋 章，刘 毅

（1.西南交通大学地质工程系，四川 成都 610031；2.西南交通大学高速铁路运营安全空间信息技术国家地方联合工程实验室，四川 成都 611756；3.高铁线路工程教育部重点实验室，四川 成都610031；4.中铁二院工程集团有限责任公司，四川 成都 610031）

流动型滑坡[1−2]是高山峡谷区普遍存在的一类地质灾害，常见的流动型滑坡包括高速远程滑坡、泥石流、火山碎屑流等。长期以来，流动型滑坡的运动机理都是地质灾害领域研究的热点问题[3]。由于这一类灾害的发生往往具有突发性、瞬时性，对于其运动过程难以进行直接观测，通常采用物理模型试验的方法对其启动-运动-停积的全过程进行再现[4−6]。随着监测方法不断完善，得以从试验中提取更多、更全面的运动学参数。其中，颗粒流厚度及其演化是碎屑流物理模型试验中重点关注的要素之一。颗粒流厚度的变化可以直观反映颗粒流形态演化过程[7]，同时，颗粒流厚度也可用于计算颗粒流的其他运动学参数[8−10]。物理模型试验中，颗粒流厚度定义为底板到底板以上观察到的大量颗粒的最高点之间的垂直距离，不包括表面不连续、松散的颗粒[11]。

根据流体运动特征的观测方法，可将颗粒流厚度监测分为拉格朗日法和欧拉法两种。目前颗粒流试验中厚度的动态提取主要以欧拉法为基础。Ancey[12]采用超声波传感器获取一定面积内的颗粒流厚度的平均值；Forterre等[7]采用光吸收法，通过试验中获取的光照强度与标定的光照强度曲线进行对比间接获取颗粒流厚度，但由于标定过程与运动过程中颗粒密实度不同，最终获取的颗粒流厚度值比真实厚度值偏小；Russell等[13]和Saingier等[14]通过激光扫描方法获取了连续高精度的颗粒流厚度，但该方法不能排除离散颗粒；Zhou等[8]、Rognon等[10]和Takagi等[15]利用激光测距传感器测量颗粒流厚度，但该方法无法排除离散颗粒；Sanjadp等[16]和Augenstein等[17]通过点探针的方法获取连续的颗粒流厚度，其厚度不受离散颗粒影响，但该方法在测量过程中探头会对颗粒流表面形态产生破坏，使测量厚度偏小；Bryant等[11]通过对高速相机影像进行手动测量获取颗粒流厚度，结果准确，但工作量大，采样频率较低。

随着图像采集设备与计算机水平的不断提升，学科跨领域的实践，数字图像处理技术被逐渐应用于工程地质领域并展现出广阔的应用前景[18−20]。苗得雨等[21]提出了基于Matlab环境下土体SEM图像处理方法；吴凯等[22]使用Image-Pro Plus图像处理软件对压实黄土孔隙微观结构进行了定量评价；张嘎等[23]在触面试验中通过数字图像技术，获得了土颗粒的运动过程；彭双麒等[24]运用PCAS图像处理软件获取白格滑坡-颗粒流堆积体形态系数；White等[25]和Stanier等[26−27]基于MATLAB编制了GeoPIV程序，以此获取颗粒流的速度场与位移场。PIV分析中颗粒流速度场的范围变化反映了颗粒流厚度的变化，但获取的厚度只是实际厚度的近似值。因此，可通过数字图像处理方法进一步处理高速相机图像序列，以获取精确的颗粒流厚度。通过数字图像处理的图像，相对于原图像可以更直接地向读者解释图像中的信息[28]，更便于机器的自动理解。

本文通过开展颗粒流斜槽试验，采用高速相机对颗粒流的运动过程进行监测并获取了颗粒流通过监测断面时的图像序列，在此基础上，通过自适应中值滤波、二值化、图像腐蚀和种子填充等[29]图像处理方法对采集的图像进行处理分析，并通过编制计算机程序实现了颗粒流厚度的动态提取，形成一套完整的基于数字图像处理获取颗粒流厚度的方法。该方法减少监测传感器，节约人力，使图像处理技术在颗粒流物理模型试验中得到进一步应用。通过改变监测位置，该方法可同时应用于颗粒流运动过程中其他形态参数的动态提取，因此具有一定的应用前景和实践意义。对于实际滑坡，该方法对于滑坡演化过程研究具有一定价值。

1 斜槽试验方案

1.1 试验装置及监测

本试验采用常规的颗粒流物理模型试验装置（图1）。该装置主要分为物料区、斜槽段及堆积槽段三个区段，分别用以对应滑坡原型的源区、运动区及堆积区。其中，斜槽段和堆积槽段通过同轴装置连接，可实现斜槽倾角调节（本试验中倾角为45°）；物料供给区与斜槽段通过挡板分隔，试验开始前，关闭挡板并将物料置于挡板之后，通过快速开启挡板实现碎屑颗粒的快速释放。模型各部分的尺寸如图1（a）所示，其中，斜槽宽度为20 cm，堆积槽宽度为22.2 cm。模型槽两侧壁均采用可透视的有机玻璃板，以从侧向监测颗粒流运动及堆积过程。

为实现颗粒流厚度的动态测量，通过侧向局部监测的方法进行颗粒流运动过程图像采集。图像序列的获取采用AcutEye高速摄影系统（16M-148CXP）（帧速率1 000 fps，图像像素1 280×800 px）。高速相机布设于挡板以下132.5 cm处，在该位置处，颗粒流接近完全展开的状态，且不受颗粒堆积过程的影响，监测区范围为14 cm×10 cm，如图1（a）所示。

图1 试验装置Fig.1 Diagram of the experimental apparatus

1.2 试验工况设计

本试验采用棱角状石英砂颗粒作为试验材料，根据土工试验标准筛筛分孔径，共设定5种不同的粒径条件（GS1—GS5），如图2所示。各粒径范围及不同粒径条件下测定所得物理力学参数如表1所示；另一方面，为对比不同体积条件下颗粒流厚度的变化特征，在每种粒径条件下设定3种体积：V1表示2 L，V2表示5 L，V3表示8 L，共开展15组（GS1V1、GS1V2、······、GS5V3）工况的颗粒流试验。

图2 各粒径影像Fig.2 Images of particle sizes

表1 试验工况设置及颗粒物理参数Table 1 Test condition setting and physical properties of the granular material

2 数字图像处理方法

图3为同一拍摄范围内不同时刻、不同工况颗粒流的高速相机影像，其中，图3（a）展示了GS3V2条件下颗粒流不同位置通过监测区的图像。从图中可以看出，同一颗粒流不同位置其厚度、颗粒运动状态均不相同，自颗粒流前缘向尾部，颗粒流厚度先增大后减小，伴随厚度的变化，颗粒流态也呈现前缘和尾部较为离散、主体部分较为密集的状态。图3（b）—（d）为9种不同工况条件下颗粒流主体部分通过监测区域时的高速相机图像。从图中可以看出，对于同一体积、不同粒径的颗粒流，其厚度近似相等，随着粒径的增大，颗粒流表面颗粒逐渐呈现离散的状态；另一方面，当颗粒粒径相同时，随着颗粒流体积的增加，颗粒流厚度显著增大，且表面离散颗粒逐渐减少。由此可见，颗粒流厚度的变化是影响颗粒流态及颗粒间相互作用的关键因素之一。

图3 颗粒流高速相机影像Fig.3 Granular flow behavior recorded by the high-speed camera

由于拍摄范围较小，且拍摄方向垂直于斜槽侧壁，因此，通过适当调节拍摄角度，可以使侧向拍摄的颗粒流厚度反映颗粒流在y方向的最大厚度值。同时，通过在图像中某一位置选取一个垂直于颗粒流运动方向的纵剖面（采样剖面），并测量颗粒流不同时刻流经该剖面的高度值，即可得到颗粒流厚度的整体变化趋势。由于单次试验中所获取的图像序列数量庞大，且图像中颗粒流与背景灰度差异显著，因此可通过编写计算机程序来实现不同时刻颗粒流厚度的提取。为更加精确地获取颗粒流主体的厚度变化情况，需采用一系列图像处理方法对原始图像进行处理。常用的数字图像方法包括图像变换、图像编码压缩、图像增强和复原、图像分割、图像描述和图像识别等。本文采用的具体数字图像方法包括自适应中值滤波、图像二值化、图像腐蚀和种子填充。首先通过中值滤波消除颗粒流内部的椒盐噪声，再依次通过图像二值化、图像腐蚀及种子填充的方法生成由颗粒流边界包围的闭合区域，该区域内，颗粒流上边界定义为与下部颗粒仍保持接触的最上部颗粒的最高点，因此不包括表面不连续、松散的颗粒，下边界则定义为斜槽面与颗粒的接触面位置处，进而通过在采样剖面测量上下边界之间的像素个数来获取颗粒流流经该采样剖面的厚度变化。数字图像处理的流程如图4所示。根据高速相机拍摄的帧速率，本试验颗粒流厚度的采样频率为1 000 Hz。颗粒流厚度的初始值为颗粒流前缘通过采样剖面时的厚度值。

图4 数字图像处理流程与相对应处理结果影像Fig.4 Flow chart of digital image processing and the corresponding processing result image

2.1 自适应中值滤波

自适应中值滤波可用于去除椒盐（脉冲）噪声，平滑非脉冲噪声，并减少图像边界失真[29]。自适应中值滤波算法是先选取一个a×b的矩形窗口Sab进行噪声监测，算法中包含两个进程A和B，其中，Smax为矩形窗口Sab允许的最大尺寸，算法流程如图5所示。

图5 自适应中值滤波流程图Fig.5 Flow chart of adaptive median filtering

其中进程A包含两个算法，分别为：

进程B包含两个算法，分别为：

式中：Zxy—坐标为（x,y）处的灰度值；

Zmin—矩形窗口中最小灰度值；

Zmed—矩形窗口中灰度值的中值；

Zmax—矩形窗口中最大灰度值。

自适应中值滤波是在图像二值化以前对原图像进行的预处理，以此来增强图像与背景的对比度，去除由颗粒孔隙与杂色颗粒造成的噪声干扰，从而得到可直接用于分析颗粒运动特征量的图像。从图4（a）可以看出，由于颗粒之间的孔隙灰度值与背景接近，染色颗粒灰度值较低，在图像中颗粒流内部会形成与背景灰度接近的暗斑，这些暗斑在图像二值化时会在颗粒流中产生空洞，从而对厚度的提取造成影响，这种暗斑即可视为图像中的椒盐噪声。通过自适应中值滤波可以有效地去除椒盐噪声造成的干扰，并能较好地保留颗粒流的边界特性，提高剖面高度提取算法的精确度，通过自适应中值滤波处理后的图像如图4（b）所示。

2.2 图像二值化

图像二值化可以通过适当的阈值反映出感兴趣的图像轮廓，使图像轮廓特征更加明显，且二值化后图像灰度值只有0或255，可以起到简化图形的作用。图像二值化是先选定一个灰度值作为阈值Zt，通过对比Zxy与Zt，若ZxyZt，则Zxy值替换为255。

通过图像二值化使图像中颗粒流与背景对比更明显，同时二值化之前图像灰度值分布在0~255之间。通过适当的阈值将图像的像素点的灰度值设置为0或255，阈值的大小需要根据光照强度进行调整，若光照强度一定，可以设定一个相同的阈值。图像二值化可以使图像数据量大大减少，提高计算机图像处理效率，结果如图4（c）所示。

2.3 图像腐蚀

图像腐蚀是数学形态学中的一种最基本的形态学算子[29]，可以快速将图形边界灰度值处理为背景灰度值。图像腐蚀技术可以通过式（5）表示。

式中：z—所有点的集合；

DΘE—E对D的腐蚀；

Dc—D的补集；

∅—空集。

通过图像腐蚀技术可以快速提取图像中颗粒流的轮廓，具体操作步骤为：选定适当大小的结构元素，对二值化图像进行图像腐蚀处理，得到失去边界的图像，通过原二值化图像P1与经过图像腐蚀图像P2做差得到图像P3。图像腐蚀掉的边界图像P3即为颗粒流轮廓图像，如图4（d）所示。

2.4 种子填充

种子填充又称边界填充算法，是数字图像处理方法中常用的填充算法。基本思想是从封闭区间内部一点开始，由内逐步向外在各个方向进行搜索，对搜索区域进行判断，不是边界则给定填充颜色，若是边界则停止填充。

基于上一步图像腐蚀获取的图像中颗粒流闭合轮廓，选定颗粒流取样剖面位于底板上方轮廓内一点作为种子点，给定填充颜色灰度值，通过种子填充算法对各像素点逐个填充，直到抵达颗粒流轮廓边界时停止填充颜色，从而实现对颗粒流主体部分填充的效果，使颗粒流主体部分与背景分离，种子填充后所得图像如图4（e）所示。

3 结果与讨论

3.1 试验结果分析

根据以上方法，基于MATLAB平台编写了颗粒流厚度动态提取程序，进行不同试验条件下颗粒流厚度的自动提取，得到的厚度（h）随时间（t）的演化曲线如图6所示。从图中可以看出，从颗粒流前缘开始h迅速增大，然后达到一个厚度变化较小的相对稳定状态，之后h逐渐减小直至到达一个临界值，随后h开始第二次上升之后逐渐减小到零，表明颗粒流已完全通过监测区域。另一方面，从图6可以看出，颗粒流厚度减小到临界值以后，相邻时间间隔的厚度值波动幅度显著增大，说明颗粒流表面开始出现显著的颗粒离散化，对应于图4中的颗粒流尾部区段。而在临界值以前，颗粒流厚度变化较为连续，颗粒流主要呈密集状态，因此该区段可视为颗粒流的主体部分。

通过对比图6中不同体积条件下的厚度变化曲线，可以看出：粒径相同时，颗粒流体积越大，整体厚度越大，颗粒流持续时间越长，主要表现为主体区段运动时间的延长，而离散段持续时间近似相等。颗粒体积相同时，不同粒径颗粒流主体部分持续时间较为接近，且厚度演化趋势近于一致，而对于粒径较大的颗粒，离散段持续时间显著减小。此外，随着粒径的增大和体积的减小，颗粒流主体部分相邻时间间隔的厚度波动幅度逐渐增大，表明颗粒流表面更加离散，该变化趋势与图3（b）—（d）中所反映的不同工况条件下的颗粒流流态变化趋势一致。

图6 各工况颗粒流厚度曲线Fig.6 Variations in granular flow thickness with time for the tests

3.2 误差分析

为分析本文中所提出的数字图像处理方法的合理性与准确性，通过将上述结果对比人工测读的颗粒流厚度值对该方法所产生的误差进行分析。针对每一组试验所采集的图像序列，每隔10张选取一张图像进行厚度的人工读取。人工读取的位置选取与图像处理中的同一列像素，以图像中斜槽表面与颗粒接触面为底边界，自下而上等间距（10 px）添加辅助点（图7），通过辅助点目测读取像素个数，并将其换算为实际厚度。通过对比图像处理获取的颗粒流厚度与目测读取获取的颗粒流厚度，可得到各组试验误差比例图（图8）（采用目测读取厚度减图像处理厚度）。从图8中可以看出，图像处理获取的颗粒流厚度能够较为准确地与手动测量的颗粒流厚度相匹配，但也存在个别误差较大的数据。还可以看出，同一体积工况下，误差比例相近；同一粒径工况下，体积增大，误差比例有一定增加，但相对趋势相同。由此可知粒径、体积等因素对颗粒流厚度误差产生的影响较小。

图7 颗粒流厚度目测读取Fig.7 Granular flow thickness which is manually measured

图8 各试验工况下误差比例Fig.8 Error ratio for the tests

为查明产生较大误差原因，将每组试验图像识别的颗粒流厚度与目测读取的颗粒流厚度进行对比，如图9所示（对比差异等于目测读取厚度减去图像处理厚度）。可以看出，图像处理获取的颗粒流厚度在颗粒流主体部分与实际厚度吻合度较高，而在颗粒流尾部二者的吻合程度相对较低，在某些工况下呈现出较大的差异。说明该方法对于密集流态的颗粒流具有良好的适用性。为分析颗粒流尾部差异较大的原因，对颗粒流尾部的原始图像与经过一系列处理后的图像进行了对比分析，结果如图10所示。可以看出，在对颗粒流进行侧向拍摄时，离散段的颗粒流颗粒在图像中只能呈现二维形态。在图像处理过程中，由于图像处理处于二维平面，而实际颗粒流斜槽试验是在三维空间中进行，三维空间降维到二维平面，本处于不同空间位置的部分离散颗粒降维之后在图像上呈现出重叠的现象，造成颗粒连续的视觉假象，图像处理由此将不连续的颗粒处理为连续的颗粒，从而导致种子填充区域远大于实际颗粒区域（图10），使获取的颗粒流厚度大于实际高度。另外，结合图8中可以发现，存在较大计算误差的均为目测读取小于图像处理，由此可推断较大误差的产生均为图像处理时识别的颗粒流边界过大，从而导致图像处理获取的颗粒流厚度偏大。

图9 数字图像处理厚度曲线与目测读取厚度曲线对比Fig.9 Comparison between the digital image processing thickness curve and the manual measurement thickness curve

图10 图像处理误差图像Fig.10 Image processing error image

4 结论与展望

（1）本文提出了一种基于高速相机获取颗粒流剖面图像，并通过自适应中值滤波、图像二值化、图像腐蚀及种子填充等一系列图像处理方式获取动态颗粒流厚度的方法。

（2）通过颗粒流的厚度剖面，可以清晰地看出厚度剖面整体趋势为先增高后降低，再升高之后降低，能够反映出颗粒流的流态化形态，且后一次增高时对比图片发现飞溅颗粒明显，由此可以推断后一次增高是由于颗粒流离散段进入了拍摄区域。

（3）利用图像处理获取的颗粒流厚度，能够准确地反映颗粒流密实段的厚度，能够较准确地反映颗粒流离散段的厚度，实现了颗粒流厚度获取的自动化，提高了效率，降低了劳动强度，精度可靠，为快速获取大量颗粒流厚度提供了新的方法。

基于数字图像处理的方式可以获取动态颗粒流厚度，由此类推，可进一步将该方法推广用于颗粒流试验形态学参数获取，如堆积体的堆积轮廓、堆积面积及堆积体扩散速率的计算等。同时，在实际滑坡中，该方法可以用于滑坡演化过程影像的处理，通过处理结果获取滑坡前缘扩散速率、覆盖范围等运动特征参数。