基于三维点云数据的集料颗粒棱角性量化评价

郝雪丽 孙朝云 耿方圆 李伟 裴莉莉 张欣

(长安大学 信息工程学院，陕西 西安 710064)

棱角性作为集料颗粒的形态特征之一，显著影响热拌沥青混合料的高温稳定性，是决定沥青混合料路用性能的重要因素。近年来，国内外学者对集料颗粒棱角性的研究取得了一定的成果。刘振清等[1]提出了基于二维图像的半径法和梯度法棱角性表征方法，并分析了棱角性指标与路面抗剪性能和颗粒指数之间的相关关系，提出了一种关于集料棱角指标的技术标准。王火明等[2]对粗集料的棱角性评价方法进行实验研究，提出了基于松装密度、振装密度和流动时间3种评价棱角方法的技术标准，并研究了针片状颗粒和浑圆颗粒与棱角性的关系。Zhang等[3]通过对粗集料的形状、尺寸、棱角性和表面纹理特性的研究，提出了基于二维图像的集料棱角性和表面纹理的综合评价指标AI。Tafesse等[4]对颗粒的棱角进行图像分析，提出了一种新的图像分析方法，并将其应用于不同的岩石样本中，该方法与已有的3种方法得到的样本棱角性排序均不同，需要继续优化设计。Yan等[5]针对颗粒形状凹凸特性，提出了一种三维固体颗粒生成算法，通过控制参数生成集料的凸/凹多面体，并引入球形度和棱角性指数对集料的形状进行评价。Yang等[6]提出了一种使用表面体素梯度向量的累计变化来量化集料棱角性，使用形态学操作前后体素的相对损失量化集料表面纹理的方法。汪海年等[7]采用数字图像处理技术，对粗集料的图像级配特征进行研究，提出了将二维数字图像级配转换为三维机械筛分级配的修正方法。张德育等[8]采用离散元方法进行沥青混合料三维离散元虚拟单轴蠕变试验，并将三维虚拟试验结果与二维虚拟试验及室内试验结果进行了比较。Su等[9]提出了一种新的量化颗粒棱角的方法，将傅里叶级数分析应用于粒子的二维图像，通过数值积分从重建的图形中计算梯度棱角指数，再对棱角指数进行评价。由此可见，目前国内外对集料棱角性的研究仍存在一些局限性，如对集料颗粒棱角性的测定缺少合理的量化评价方法，同时集料棱角性的研究大多基于二维图像进行分析处理。三维数据的优势在逐渐被认知到，因此有必要进一步开展集料三维形态量化的研究。

本文提出了两种基于三维点云数据的集料颗粒棱角性量化评价方法：基于Sobel-Feldman卷积和基于集料表面法线的集料颗粒棱角性量化方法，并将计算结果与AIMS的棱角性量化结果进行了比较，进而确定出更为合理的一种评价方法。

1 集料颗粒和三维点云数据的选取

1.1 集料颗粒选取

本文选取的是人工集料，即是岩石等经过人工机械破碎后，经筛分得到的表面粗糙且具有棱角的颗粒碎块。实验初步采用粒径为9.5～13.2、13.2～16.0、16.0～19.0 mm的石灰石、玄武岩和花岗岩集料颗粒，其中每种岩性、每档集料各50颗作为实验样本。此外，为了更好地研究集料颗粒的三维形态特征，本实验对筛分得到的集料颗粒进一步优化筛选，共挑选了5种几何形状(六面体、五面体、四面体、片状、细长状)的颗粒，如图1所示。

图1 集料分类

1.2 三维点云数据选取及采集

本文采用Gocator 3100系列双目智能传感器，搭建基于Gocator的三维图像采集系统，包括相机、光源、硬件设备等。使用Gocator采集设备对这三种岩石各档位进行三维数据采集。

如图2所示，通过以太网束连接到PC客户端和Gocator网络进行配置，完成集料颗粒三维图像采集的前期工作。

图2 Gocator 3100三维数据扫描整体图

三维数据获取具体分为以下3个方面：

(1)结构光调制器(SLM)使用蓝色LED灯，以快速的序列向待测物体投射一系列高分辨率/高对比度的光，如图3(a)所示。

(2)两台摄像机从不同的视角捕捉反射光的信息，如图3(b)所示。

(3)传感器根据立体相关原理从光模式中产生三维点云数据，如图3(c)所示。

(a)结构光投射物体示意图

(c)三维点云图

左摄像机

右摄像机

Gocator输出的图像一般为16 bit的RGB图，需要使用Halcon13软件对图像数据进行三通道拆分处理，集料三维图像示例如图4所示。

图4 集料三维图像

2 集料颗粒棱角性评价规范

2.1 集料颗粒的棱角性特征

集料颗粒的形态特征主要包括粒形、棱角以及表面纹理，如图5所示。集料颗粒的棱角反映集料颗粒宏观局部的变化状态，表示集料颗粒表面尖锐程度或者凹下、凸起的程度。集料颗粒棱角性和沥青混合料的强度有着紧密的关联，可以使集料间产生较大的摩擦力，以防止集料的流动，从而使得路面结构更加稳定、抗滑。

图5 集料颗粒形态图

2.2 AIMSⅡ集料颗粒棱角性评价方法

现有规范中的集料颗粒棱角性评价方法存在需人工操作、效率较低、误差较大等问题，而AIMSⅡ[10- 12]中提供的方法是现行比较公认的集料颗粒棱角性量化方法，因此本文将其作为比较的标准。

2.2.1 AIMSⅡ集料颗粒棱角性评价原理

AIMSⅡ集料图像采集系统由图像采集硬件设备和计算机组成。如图6所示，图像采集硬件设备是一个封闭系统，包括高分辨率的相机、光源和放置集料的托盘。

AIMSⅡ量化棱角是通过计算集料颗粒表面轮廓凸起的尖锐程度的变化率来表示集料的棱角性。梯度棱角(GA)是量化棱角的主要指标，其取值范围为0～104，越接近圆形的集料颗粒的GA值越小，圆形的集料颗粒的梯度棱角值为0。计算集料

图6 AIMS样品摆放和试验系统内部结构图

二维图像轮廓相邻点之间的梯度向量倾角的变化量，即可得到集料的梯度棱角值，GA的计算公式为

(1)

式中：θi为集料颗粒在轮廓图像上第i点位置的梯度向量值(如图7所示)，n为集料颗粒轮廓图像边缘点的个数。

图7 集料颗粒梯度向量

2.2.2 AIMSⅡ集料颗粒棱角性结果分析

使用AIMSⅡ测量软件检测粒径为9.5、13.2、16.0 mm石灰石集料的四面体、五面体、六面体集料颗粒，其梯度棱角GA如图8所示。从图中可以看出：粒径为16.0 mm的集料相较于粒径为 9.5 mm的集料，GA值整体上更平缓，即随着粒径的增大，GA值同比减小，棱角性变好；相同粒径下六面体集料的GA值几乎都是小于四面体和五面体集料的，说明四面体集料的棱角性最丰富，六面体集料更光滑且更稳定。

图8 石灰石集料颗粒的AIMS统计结果

3 基于Sobel-Feldman卷积的集料颗粒棱角性量化方法

3.1 Sobel-Feldman卷积算法的基本原理

为了避免噪声的影响，在构造梯度算子时还要考虑平滑处理，这样能在滤除噪声的同时保留集料的棱角。Sobel-Feldman算子包括两个操作：垂直于导数方向的平滑和导数方向的差分，即该算子包含横向和纵向两组卷积因子，分别将图像与之进行卷积，便可得出横向和纵向的差分近似值。平滑算子可以从帕斯卡三角形[13]中得到，运算表达式如下：h(-1)=1，h(0)=2，h(1)=1，h′(-1)=-1,h′(0)=0,h′(1)=1。

根据以上算子，在二维Sobel-Feldman运算中，x和y方向的梯度矩阵可以表示为

(2)

则x和y方向的卷积核分别为

如果A为输入图像，Gx和Gy分别为横向和纵向的输出图像，则有

(3)

(4)

即输出图像G为

(5)

通常，为了提高效率，输出图像也可以使用绝对值近似计算，即

|>G|=|>GX|+|>GY|

(6)

计算梯度的方向

(7)

得到集料的输出图像后，根据集料中心点，遍历集料边缘上所有像素点，每10°选取一个边缘像素点(如图9所示)，计算该点以及相邻点的梯度方向差值，遍历求和得到集料的棱角性指数(AI)：

(8)

式中，θα为某点梯度的方向。

图9 梯度像素取值图

3.2 优化的Sobel-Feldman卷积算法

一般Sobel-Feldman操作是使用3×3大小的卷积核进行卷积，但这样的尺寸对边缘变化敏感，因为AI值是基于表面像素梯度的变化，而梯度是根据周围的像素确定的。本文选择增大卷积核，这样能够考虑更多的周围像素，从而使得梯度的变化变小，计算结果更准确。Sobel-Feldman卷积算法流程图如图10所示。

图10 Sobel-Feldman卷积算法流程图

根据帕斯卡三角形可得到5×5卷积核如下：

同理,可得到该卷积核下的棱角性指数。

在3.1节集料颗粒棱角性的量化只针对边缘点的一个像素值，为了使棱角量化结果更加精确，这里采用相邻的两个像素来计算，选取边缘点附近的两个像素的平均梯度方向作为该点的梯度方向，如图11(b)所示，然后继续对集料进行旋转，得到集料的棱角性指数。

3.3 基于Sobel-Feldman卷积的集料颗粒棱角性评价

针对一个像素下不同卷积内核计算集料棱角性

图11 集料颗粒的梯度方向

指数AI，不同粒径下不同形状的石灰岩集料颗粒的AI统计结果如表1所示。由表中可以看出，3×3卷积核下计算出来的AI结果比5×5卷积核下的计算结果均偏大，说明较大的卷积核对周围情况不敏感，即卷积核越大，计算结果越趋于稳定。

表1 不同卷积核下集料棱角性指数统计结果

进一步使用相邻两像素的平均梯度值代替某一点的梯度值对Sobel操作进行优化，AI结果如表2所示。由表中可以看出，使用单个像素和使用相邻两个像素梯度计算的AI间没有显著差异，当像素点在集料颗粒的突变锐点上时，使用单个像素梯度计算的AI值会大于使用两个相邻像素梯度计算的AI值。

表2 不同计算方式下集料棱角性指数统计结果

如图12所示，考虑到AI值的稳定性，一般情况下会优先考虑使用相邻两像素梯度的均值来计算集料的AI。

图12 不同位置的梯度方向示意图

综合上述分析结果可知：

(1)扩大卷积内核后的Sobel-Feldman卷积相比3×3内核，AI值更低，这是由于较小的内核对周围的变化更敏感。

(2)相邻两像素法比单个像素法更稳定。这是由于曲面的细节、尖锐点上的像素会与相邻像素在梯度方向上产生巨大的差异。

(3)与常规算法相比，改进后的Sobel-Feldman卷积对图像分辨率的敏感性降低，对不同方法得到的图像鲁棒性更强。

4 基于表面法线的集料颗粒棱角性量化方法

表面法线是三维模型的重要指标，为了量化集料的棱角，可通过计算集料表面的法线，对法线方向进行聚类来得到集料表面个数[14]，从而分析出集料表面个数与集料棱角间的关系。物体表面法线示意图如图13所示。

图13 表面法线示意图

4.1 主元分析法理论

由于点云的数据量比较大，为了简化计算，本文使用一种比较高效的点云法线估计方法——主元分析法(PCA)。如图14所示,PCA是一种线性投影，不仅可以用来提取数据的特征分量，还可用于计算点云数据的主方向。PCA的思想是将n维特征数据映射到k(k

(1)将需要处理的数据组成n×m矩阵X。

(2)将矩阵X的每一行数据进行零均值化处理，即每个数据减去该行的均值。

(3)求出协方差矩阵。

(4)利用数学方法解出协方差矩阵的特征值及其对应的特征向量。

(5)按特征值大小顺序将特征向量从上到下按行排列成矩阵，并取前k行组成矩阵P。

(6)Y=PX即为降维到k维后的数据。

图14 PCA示意图

为了使得协方差矩阵的特征向量就是理想的K维特征，PCA必须满足两个原则：

(1)最大方差原则

一般希望投影后的投影值尽可能分散开来，这种分散，在数学上用方差来表示，代表每个元素和该段均值相差的平方和的均值，表示为

(9)

(2)最小平方误差原则

4.2 基于最小二乘局部平面拟合的法线估计

根据主成分分析法，对点云法线的求解可以通过最小二乘法来求解局部数据点的拟合平面[15]，然后以此平面的法线来表示待求数据点的法线。因此，估计集料表面法线的方案变为求解集料点云数据的局部拟合平面，并分析该邻域内协方差矩阵的特征矢量和特征值，找到特征的最小值，该最小值所对应的特征向量即为拟合平面的法向量。

对点云中的任意一点xi，选取距离点云xi最近的K个邻域点，根据最小二乘原理，可解出这些点所拟合的局部平面P，该平面的法向量记为n，d为点到拟合平面的距离，则拟合平面P可表示为

(10)

(11)

协方差矩阵C存在的特征值及特征向量为

(12)

由于本文只需要数据点法线方向的特征值和特征向量，不妨设n0=(a,b,c)，则有

(13)

进而

(14)

解得

(15)

依据最小二乘原理，当该表达式取得最小值时，其对应的特征向量即为法线的方向。集料表面法线如图15所示。

图15 集料颗粒的表面法线效果

4.3 三维集料颗粒多维数据聚类方法

由于得到集料的相邻面法线方向是不同的，且同一平面上的法线方向趋于平行，因此本文提出了一种基于法向量的三维数据聚类方法，用于判断集料颗粒所在面的个数，从而分析集料的棱角性。

聚类分析是根据数据点的相似性来对元素进行分类的，同一个面上的法线坐标信息具有相似性，故本文采用基于快速搜索与寻找密度峰值的聚类方法[16]对法线坐标数据进行聚类，它是一种新颖的聚类算法，能对任意形状、任意维度的数据进行聚类，所需参数少并且对噪声不敏感。该算法主要包括两个部分：首先通过用户输入的参数值，计算每个样本的局部密度及距离，找到样本聚类中心，并根据决策图选择合适的聚类中心；然后将剩余的其他样本点分配到与其最近的高密度邻域的集群中。具体实现过程如下。

对于每个样本数据点i，计算其局部密度ρi和它与高密度点的距离δi。这两个量都只依赖于数据点i和j之间的距离dij，假设dij满足三角不等式，则定义数据点i的局部密度ρi为

(16)

式中：当x<0时χ(x)=1，否则χ(x)=0；dc为截止距离；ρi为比dc更接近点i的点的个数。

计算点i与任何其他密度较高的点之间的最小距离为

(17)

对于局部密度最高的样本点，δi=maxjdij。值得注意的是，对于密度上的局部或全局最大值点，δi比典型的最近邻距离大得多。因此，聚类中心被认为是δi值异常大的点。

图16(a)显示了二维空间中的点，容易看出，密度最大值在点1和点10处，称这两个点为聚类中心；26、27和28属于一个集群，10、13、15、19、22属于另一个集群，其余的属于集群1。图16(b)是δi作为各点ρi的函数曲线图，又称为决策图，可以看出，点9和点10的ρ值虽然类似，但点9属于集群1而点10属于另一个集群；正如预期那样，点10具有高的δ值和相对高的ρ值，可称为聚类中心点；点26、27、28的δ值相对较高而ρ值相对较低，因为它们是孤立杂点，所以可以看作是由单个点组成的集群。

图16 二维数据点集聚类示意图

在找到集群中心之后，每个剩余的点被分配到与其最近的高密度邻域的集群中。在该算法中，为每个集群或者簇找到一个边界区域，定义为分配给该集群的点集，与属于其他集群的数据点之间的距离为dc。每个集群的边界区域中密度最高的样本的密度称为临界密度(记为ρb)，将集群中密度大于ρb的点视为集群中的一部分，反之视为噪声。

采用以上算法对集料颗粒法线进行聚类处理，粒径为9.5 mm的石灰岩四面体、五面体及六面体集料的法线聚类结果如图17所示。从图中可以看出，在选择聚类中心之后，每个就近点被分配到一个集群，经分析可知，同一集群中的点属于同一个平面。

图17 石灰岩集料颗粒的聚类决策图和聚类结果

为了进一步分析聚类个数与集料表面个数之间的关系，对石灰岩集料颗粒的聚类结果进行了统计，结果如表3所示。可以看出，四面体的聚类个数最多，五面体次之，六面体的聚类个数最少。究其原因，集料颗粒在扫描为三维点云数据时，只是针对集料的表面进行数据获取，标准四面体为三棱锥形状，扫描到的图像呈3个面；标准五面体为三棱柱形状，扫描到的图像呈2个面；标准六面体为立方体，扫描到的图像呈1个面。对集料颗粒来说，还存在很多外界因素使得扫描到的形状有差异，会对集料的聚类结果产生影响，但总体趋势一致。

为了验证分析结果，采用橡皮泥样本数据作为对比，具体方法为：将橡皮泥制作为标准试件模型，然后进行图像采集，获得其点云数据的法向量，并用聚类法对法向量进行聚类，结果如图18所示，最终得到的聚类结果依次为3、2、1，表明以上结果成立。

表3 石灰岩集料颗粒的聚类结果的统计结果

(a)四面体

(b)五面体

(c)六面体

5 集料颗粒棱角性量化结果分析

本文使用3种不同的方法(AIMS集料量化方法、基于Sobel-Feldman卷积的集料棱角性量化方法以及基于表面法线的集料棱角性量化方法)对集料的棱角性进行量化，石灰岩(粒径为9.5、13.2、16.0 mm)的AIMS与基于Sobel-Feldman卷积的集料棱角量化多项式拟合结果如图19所示。

图19 石灰岩集料的AIMS与基于Sobel-Feldman卷积的棱角性拟合结果

可以看出，对于粒径为9.5～16.0 mm石灰岩集料，基于Sobel-Feldman卷积的方法的棱角性计算结果与AIMS的棱角性量化结果具有比较显著的正相关关系，相关系数分别为0.868 2、0.905 3和0.880 2，因此可以认为，基于Sobel-Feldman卷积的集料颗粒棱角性量化评价方法测定集料颗粒棱角性是合理的。

进一步对集料使用相邻两像素的Sobel-Feldman计算方法求棱角性AI，分析不同岩性集料颗粒的棱角性趋势，结果如图20所示，其中每档粒径集料依次包含四面体、五面体、六面体集料颗粒各10颗。由图20可以看出，各档粒径石灰岩集料的棱角性指数均小于玄武岩和花岗岩，说明在本次试验中石灰岩的棱角性最好，玄武岩最差，六面体集料棱角性明显比四面体和五面体好，且随着粒径的增大，棱角性指数同比变小。由4.3节可知，在基于表面法线的聚类方法中，六面体集料颗粒的最终聚类个数最少，四面体集料颗粒的聚类个数最多。由此可知：集料越尖锐，集料棱角性指数越大，法线聚类个数越多；集料越圆滑，集料棱角性指数越小，法线聚类个数越少。

图20 集料颗粒的AI结果对比

6 结论

本文基于三维点云数据，采用不同尺度下的3种方法对不同粒径、不同岩性、不同形状的集料颗粒的棱角性进行分析。结果表明：粒径为9.5～16.0 mm石灰岩集料的Sobel-Feldman卷积法棱角性计算结果与AIMS棱角量化结果具有显著的正相关关系；随着粒径的增大，棱角性指数同比变小，且集料越尖锐，集料棱角性指数越大，法线聚类个数越多；集料越圆滑，集料棱角性指数越小，法线聚类个数越少。因此，基于Sobel-Feldman卷积的集料棱角量化方法评价集料颗粒棱角性具有更高的准确性，可为沥青混合料中关于可压实性、高温抗车辙性、低温抗裂性、抗滑性及虚拟力学模型的研究提供参考价值。