基于图像处理的道砟颗粒二维廓形追踪提取算法

李 朋,刘鉴兴,肖杰灵,王 平

(1.西南交通大学高速铁路线路工程教育部重点实验室,成都 610031； 2.西南交通大学土木工程学院,成都 610031)

随着我国铁路提速与高速化发展，有砟轨道道床病害日益突出。通过离散元法来模拟有砟道床细观、宏观力学行为与劣化机制，是发展和完善高速铁路有砟轨道结构设计与运营维护技术的一项重要研究手段[1-3]。其中，精确提取道砟颗粒廓形是保障仿真结果可靠性的重要基础。除服务于离散元仿真分析外，准确的道砟颗粒廓形还可以应用于评估散体道床道砟堆积体的物理特性，如堆积密实度、道砟棱角系数、道床级配等[4]。

在提取颗粒的几何形态时，通常采用数码相机、激光扫描或CT设备等计算机视觉测量技术。其中，由于数码相机的技术复杂度及成本低于激光扫描等方式，因此其被广泛用于获取颗粒的二维特征[5-6]。然而，在运用数码相机拍摄的数字图像来提取道砟颗粒二维廓形时，使用的传统边缘检测方法，如Sobel算子、Roberts算子等虽然具有实现简单，运算速度快等优点，但其抗干扰能力差，边缘宽度往往比实际大，且边缘像素孤立或分小段连续[7-9]。不仅降低了道砟颗粒廓形的提取精度，而且孤立和分断的无序点云也不利于后期道砟形态表征的评估。一些学者提出了不同的边界检测改良算法[10-12]，但这些算法大都针对特定的情况和物体，完全套用这些算法，亦不能获得理想的道砟颗粒二维廓形。

为提高道砟颗粒二维廓形的连续性和准确性，同时得到有序的廓形数据点列以弥补无序数据点云在分析道砟形态特征时的不足，本文针对性地提出一种基于图像处理的道砟二维廓形追踪提取算法。算法特点在于抗外界环境干扰能力高，能够一次性处理多个道砟颗粒，提取的道砟颗粒二维廓形不仅连续性好，而且检测出来的边缘是单像素宽。通过“追踪”的思想，使提取的廓形数据形成有序点列。通过将算法提取的廓形数据与游标卡尺测量的实际道砟廓形数据对比分析，发现两种粒径的差值均在1 mm之内。考虑测量误差，可以肯定本文算法的有效性与可靠性。

1 搭建消影拍照台

随意地将道砟放置于多变的自然环境中，由于光照条件不易控制，会有很多复杂的噪声影响，降低图像的拍摄质量。为减小不良因素对图像的影响，同时避免繁琐的图像预处理步骤，需要搭建一个特殊的消影拍照台，该拍照台由数码相机、白色背景板、透明玻璃板和LED灯组等构件组成，如图1所示。

图1 消影拍照台

2 标定与图像预处理

在拍摄道砟颗粒的二维图像之前，要进行标定以获得相机的内外参数，进而获得道砟颗粒的尺寸等信息。本文采用了经典的棋盘格标定法进行标定。首先制作了棋盘格标定板，并拍摄了20张其不同角度的图像。不同角度拍摄所得道砟颗粒的图像有较大区别。针对二维廓形的研究，需考察道砟颗粒具有代表性角度下拍摄的图像。使道砟颗粒从距地面0.5 m自由落下，待其在平坦地面上稳定后，将其移至拍照台进行拍摄。此后运用MATLAB标定工具箱获得了相机的内外参数，用以计算得出道砟颗粒的真实尺寸。

拍摄完标定板后，将6个道砟颗粒均布于白色背景板上(后期实践证明，选择与道砟颗粒颜色差异大的背景色即可，如绿色、蓝色等)，每批次摆放的道砟颗粒大致在同一位置。

拍摄的道砟颗粒图像无需进行图像预处理，而且由于拍摄时使用的背景与道砟颗粒对比度很大，也不需要进行灰度变换。如有特殊情况，本文推荐采用的主要图像预处理方法有：图像平滑、图像去噪、图像锐化等[13]。进行预处理的目的是对不符合要求的图像进行某些适当的变换，从而使获得的道砟颗粒图像清晰并且满足后期处理的要求。

3 道砟颗粒二维廓形追踪提取算法

算法主要包括颜色空间转换与二值化处理、获取道砟边界、初始化、追踪循环等过程。其流程如图2所示。

图2 算法流程

3.1 颜色空间的转换与二值化处理

为了方便分辨图片中的像素是属于背景色还是属于道砟颗粒，注意到明度最能反映道砟颗粒与背景的差异，且RGB颜色空间不方便描述颜色之间的距离，因而考虑将RGB颜色空间转化到Lab颜色空间[14-16]。两者之间的转化用XYZ颜色空间传递，即先将RGB颜色空间转化到XYZ颜色空间，再将XYZ颜色空间转化到Lab颜色空间。本文采用一种通用的空间转化算法[17]，具体参数取值通过试验获得。

在将RGB颜色空间转化到Lab颜色空间后，设置Lab颜色空间的明度(即L、a、b中的L矩阵)阈值，将提取的明度矩阵转化为二值矩阵，用来表征道砟颗粒与背景色。基于上述操作，可快速通过梯度算子求出二值矩阵的梯度，从而得到二维图像中道砟颗粒的边界像素点集。

3.2 追踪提取二维道砟廓形

在一幅二维图像中为每一个道砟选取一个位于道砟颗粒轮廓线以内的点作为起点(该点应大致位于道砟颗粒中心)。因前文要求每批次道砟颗粒摆放位置大致相同，所以后续不同批次的道砟颗粒都可以用第一次选点作为起点，这样可提高数据处理的速度。

3.2.1 初始化

(1)寻找起始点

用50×50的像素矩阵(本文结合二维相片中道砟颗粒的边界宽度，选取50×50的像素矩阵。可根据实际情况调整，如设置40×40或60×60等不同形式的像素矩阵)从起点向正上方移动，找到道砟颗粒的边界，即判断像素矩阵框内是否存在道砟颗粒边界像素点。具体方法如下：

在Lab颜色空间中，边界像素点的值为非0数(一般接近1)，其余空白像素点部分的值均为0。所以在50×50的像素矩阵范围内求和

(1)

其中，Ni为第i个像素点的值。当C非零时，就可以说明在此范围内存在至少一个像素点，从而将矩阵中心移动至此，该像素点称为0号中心像素点；当C为零时，则说明在此范围内不存在像素点，则继续移动像素框。

(2)初次拟合

此步骤为追踪的关键，目的在于运用最小二乘法原理确定像素矩阵沿道砟廓形的移动方向，如图3所示。

图3 初次拟合示意

已知通过颜色空间转换和梯度算子求得的道砟颗粒边界像素点集为

(2)

将像素矩阵中的道砟颗粒边界像素点集令为[A]，即

(3)

设直线拟合系数为(k，m)，直线表达式为y=kx+m，则有

(4)

(5)

(3)初次步进

在求得直线后，以10个像素为步长(为保障提取精度及效率，10个像素为宜)将中心像素点按照直线任意方向移动(“任意”只适用于初次步进)，此时选取的方向线与水平正方向形成夹角，如图4所示为正方向情况。

图4 步进示意

移动以后判断像素矩阵框内是否存在道砟颗粒边界像素点，若存在，则将50×50的像素矩阵中心移至距离最近的道砟颗粒边界像素点上。

3.2.2 追踪

此步骤是本文算法的核心，与初始化步骤中的内容大致相同，但需做一些改动，具体内容如下。

(1)拟合

运用最小二乘法原理，求出移动后新像素矩阵内廓形像素点的拟合直线，根据此直线确定下一次步进方向。

(2)步进

此时需保证步进方向与初次步进时相同(顺时针或逆时针)，但确定步进方向时的线为直线而非射线，会有两个相反方向。观察到像素框移动步长小时，道砟颗粒的二维轮廓梯度改变量不会太大，所以前一次产生的“α1”和本次拟合直线与水平正方向产生的夹角“α2”应相差不会太大，即如图5所示，其中

α3=α2+180°

(6)

图5 α2和α3示意

基于上述条件，可设定夹角参数θ来判断50×50的像素矩阵沿顺时针(或逆时针)运动的方向。其中

θ=|α1-α2|<180°

(7)

(3)校准

步进后，若在50×50的像素矩阵范围内存在道砟颗粒边界像素点，则将像素矩阵中心移至距离最近的道砟颗粒边界像素点上，如图6所示。移动到距离最近点的方法如下。

设50×50的像素矩阵中心所在位置坐标为(a，b)，像素矩阵范围内存在的任意一个道砟颗粒边界像素点的坐标为(xj，yj)，两点之间的距离可以表示为

(8)

利用计算机可以快速计算出minLj，则校准后的点可表示为

(a′，b′)=(a，b)+minLj·δ

(9)

其中，δ表示单位方向向量，则定义移动的方向。

图6 校准示意

3.2.3 输出有序二维廓形数据点列

进行下一次移动，每次移动后判断新的边界位置是否与起点边界位置重合或接近(小于5像素)，即是否回到0号中心像素点，如果回到0号中心像素点则停止移动并输出所记录的点列，将每一次的中心像素点坐标连线，绘制拟合的道砟颗粒二维廓形图；否则继续追踪。在第一颗道砟廓形寻找完成后，由于前期将各批道砟均摆放在相同位置，便可定义各颗道砟50×50的像素矩阵的初始点，由此循环上述“初始化—追踪”步骤，即可一次性获取6颗道砟的二维廓形。

3.3 二维廓形点列数据优化

通过上述廓形提取算法获取道砟颗粒二维廓形表示为

C={(xi，yi)|i=1，2，…，n}

(10)

n的取值根据道砟大小及相机精度不同而变化，本文中n在100至300内波动。

如图7所示，由于直接获取的廓形数据由大量密集点列组成，而在一定范围内，大多数点对于廓形变化的表征作用并不明显，造成了点列的冗杂。点列数据越多，对后期仿真效率的影响也就越大，因此针对一些不需要精细道砟颗粒二维廓形的情况，可对数据进行压缩优化。优化的目的是在不影响道砟颗粒整体形状的情况下，对边界节点个数进行精简，也即忽略部分纹理、变化率较小的表面波动，而将道砟颗粒的形状表达出来。下面说明优化方法。

图7 数据优化示意

设置一个角度变化阈值，每积累超过该阈值则引入一个新的点，在小于该阈值时则减少点数以省略不关注的细节。以此达到用最优的方式记录廓形边界数据的目的。

步骤1：设定一个角度控制阈值，该阈值用于表征边界廓形的变化剧烈程度。阈值越小，对廓形的描述越精确，但廓形数据的数量就相对较多；角度越大，对廓形的描述则越粗糙，廓形数据的数量就相对较少。

步骤2：定义角度记录变量β，用于存储角度累计值。

步骤3：定义角度增量βi，用于表示相邻两点之间的角度增量值；初始β=0。

步骤4：选取任意一点为迭代初始点，计算输入的廓形点列中相邻两个点之间的角度βi

(11)

判断β=β+βi是否超过角度控制阈值ρ

(12)

超过则记录当前节点坐标，并重置β=0；否则累计角度记录变量β，再进行下一个点的计算。

步骤5：每次计算后判断是否已经遍历了所有节点，即所计算的点与迭代初始点坐标信息是否相同。相同则停止计算，输出记录点，形成新的多边形廓形，否则继续循环此计算。

步骤6：设求出的经过压缩的点为

P={(xj，yj)|k1，k2，k3,…，kN}

(13)

则在遍历处理前的所有点后，所求得的解的相邻两个点之间应该满足下列条件

(14)

4 实例分析及可靠性验证

4.1 实例分析

拍摄照片：本文选用了遮光纸箱、铝合金架、数码相机、A3大小的白纸、透明玻璃板以及LED光源等组成消影拍照台。利用蓝牙遥控器控制数码相机拍摄道砟颗粒的照片。

颜色空间转换：将照片信息导入Matlab软件，进行从RGB颜色空间到Lab颜色空间的转化并求取廓形像素点集。

追踪提取二维廓形：根据算法编写程序对转换后的图片进行道砟颗粒二维廓形追踪提取工作。图8中蓝色的闭合曲线为提取的6个道砟颗粒二维廓形。可以看出提取的廓形与道砟颗粒拟合度良好。

二维廓形点列数据优化：为方便叙述，现用提取的二维廓形图中的第一行第一列道砟颗粒轮廓单独说明。编写程序，对提取的道砟颗粒二维廓形进行数据优化处理。图中红线为优化廓形结果，红色圆圈为最终保留的轮廓数据。

图8 实例分析

4.2 算法可靠性验证

除通过观察提取的廓形与道砟颗粒外形拟合度来了解算法可靠性外，本文在像素坐标系下获得道砟颗粒廓形数据后，通过标定得到的相关信息，反算出廓形在真实坐标系下的尺寸，再利用游标卡尺量取真实颗粒的最大粒径与反算结果进行比较，旁证本文算法的有效性和可靠性。

通过上述二维廓形追踪提取算法得到图8示例颗粒的轮廓数据点共186组，同时求得其中两点之间最大距离为1 007.542 1像素距离。由标定反算得到像素距离与真实距离之间的关系为

(15)

由此可以得到像素坐标系下最大粒径反算结果为

L=距离×S=1 007.542 1×0.094 0=94.709 0 mm游标卡尺测得的道砟真实粒径为95.26 mm。通过对比，可知获得的道砟廓形结果与真实粒径相差0.551 0 mm。

随后用游标卡尺测量了多个道砟的最大粒径，两种粒径的差值均在1 mm之内，考虑到测量误差的存在，可以认为提取的道砟颗粒廓形数据有效可靠。部分测量结果和误差分析结果如表1所示。

此外，棱角指数[18]是定量描述颗粒几何特征的最为关键的参数，本文采用Rao C等[19]提出的方法计算颗粒的棱角指数。基于图8所提取出的实际廓形，将其轮廓简化为n边形。把区间[0°，180°]分为18个长度为10°的等值区间，计算n边形相邻顶点对应的内角变化值在这18个区间上的概率分布，按下式计算棱角度A

表1 部分颗粒测量结果及误差分析 mm

(16)

式中，e为区间E=[e，e+10]的起点值；P(e)为多边形相邻内角变化值α=|αi+1αi|位于区间E内概率。

用式(17)求三个投影图棱角度的加权平均值，即为棱角指数AI。

(17)

式中，k为投影图编号；S(k)为投影图k的面积。

表2列出了Masad E等[20]测定的典型碎石的棱角指数值和本文中所采取算法计算的棱角指数，可看出同种碎石的棱角指数相近，进一步旁证本文算法的有效性和可靠性。

表2 典型碎石的棱角指数范围

5 结论

针对道砟颗粒的二维廓形，提出了基于图像识别的道砟颗粒二维廓形追踪提取算法。通过实例，分析了算法的精度及耗时，得出以下结论。

(1)RGB颜色空间不方便描述颜色之间的距离，Lab颜色空间更加便于分辨图片中的像素是属于背景色还是属于道砟颗粒，且明度最能反映道砟颗粒与背景的差异。

(2)改变像素矩阵和移动步长的大小等参数，可以对道砟颗粒二维廓形数据提取的精度和时间进行控制。精度越高，则耗时越多，反正亦然。

(3)本算法提取的二维廓形数据为有序的点列，而非无序的点云数据。

(4)提取出来的廓形与实际廓形的最大粒径在考虑测量误差及精度的条件下二者数值基本相同。本算法适用于提取道砟颗粒的二维廓形。