基于模式识别技术的铁路道砟模型库建立及离散元模拟

贺 勇，张宗堂，樊宝杰，柳光磊，张巨峰

(1.长沙市规划勘测设计研究院，湖南 长沙 410007；2.湖南科技大学 岩土工程稳定控制与健康监测湖南省重点实验室，湖南 湘潭 411201；3.湖南科技大学 资源环境与安全工程学院，湖南 湘潭 411201；4.湖南省地质矿产勘查开发局四一四地质队，湖南 益阳 413000)

1 概述

有砟铁路铺设了大量的道砟颗粒[1-3]。传统的颗粒形状参数测量法通过手动测量或与标准图表视觉比较来确定参数[4]。这些传统方法存在很多局限性，例如易受噪声，振动的影响[5]，手动测量和视觉比较法测量非常主观且缓慢[6]。因此，需要自动、快速、准确的铁路道砟颗粒测量方法来代替传统的测量法。

当颗粒彼此不接触的情况下，获取图像并进行图像分析相对容易。RASCHKE、BANTA、ALTUHAFI等[7-9]都使用了这种分离方法。但在实际情况下，铁路道砟颗粒都是彼此接触的，因此难以对其进行颗粒分离。

在接触颗粒的图像中，计算铁路道砟颗粒形状更加困难，首先需要在图像中识别出颗粒，并准确定义其边界。此外，考虑三维铁路道砟颗粒，道砟颗粒间不仅彼此接触，而且会有颗粒重叠，在相机视图中存在遮挡。因此，只有小部分道砟颗粒可以完整地显示二维投影。只有图像中完整地显示出连续轮廓的颗粒才可用于形状测量。因此，在ZHENG等[10]提出的一种半自动方法中，仍然需要科研人员在道砟颗粒的图像中挑选出完整轮廓颗粒。

本文基于模式识别，提出了一种自动识别完整轮廓颗粒的算法，省去了人力工作，对道砟颗粒形状评估结果进行统计分析。建立了颗粒形状库，用于在数值模型中随机重建具有所需形状值的铁路道砟颗粒。基于所提出的道砟颗粒形状库和离散元方法，对真实形状的道砟颗粒进行了双轴压缩试验，研究了道砟颗粒细长度对其力学特性的影响。

2 道砟颗粒的自动获取

为节省时间和成本，并考虑到铁路道砟样本数量多、精度要求低的特点，本文采用了低精度的摄影设备(例如手机和数码相机)进行拍摄，获得了多个铁路道砟颗粒图像，如图1所示。本文提取的目标是道砟颗粒的不规则形状。

图1 道砟图像示例

由于算法的简单性、高精度性和高计算效率，在本文中使用经过改进的Viola-Jones算法[11]来自动识别道砟颗粒。

整个过程主要分为4个主要步骤：训练集准备，训练集图像的缩放和旋转，Adaboost和Cascade训练分类器，基于滑动窗口方法识别颗粒。

2.1 训练集准备

首先，考虑到道砟颗粒的最常见形状，其长宽比基本为1：1。拍摄1 206张道砟照片，对其中的一些道砟图像进行裁剪、缩放并作为模型的训练数据集。每个裁剪后的图像都是原图像中的一个小方块。其中，包含道砟颗粒完整轮廓的图像被标记为正图像，而其他图像被标记为负图像。所有的裁剪图像的大小均为32×32像素。将这些图像作为输入数据，计算机可以进行模型训练和学习来获得模式识别中的Haar-Like特征[12]。并将特征结果存储在分类器中。试样图像如图2所示。

(a)正图像

2.2 训练集图像的缩放和旋转

大多数铁路道砟颗粒的长宽比为1:1至2:1。对于不处于该比例的道砟，很难通过拍摄的方式来大量收集图像数据，但这些高长宽比的道砟颗粒是真实存在的，所以分类器也需要检测到这些道砟。因此，需要一种新的获取大量图像的方法。为了快速获得大量可用于训练的细长道砟颗粒图像，本文将道砟长宽比为1:1的现有图像进行数字拉伸[13]。通过调整之前的1:1训练图像(32×32像素)的大小，生成1.5:1训练图像(48×32像素)，2:1训练图像(64×32像素)，3:1训练图像(96×32)像素)和4:1训练图像(128×32像素)，如图3所示。

(a)1.5:1图像

为了在任何方向上都能识别细长道砟颗粒，原始图像以15°的增量逆时针旋转，如图4所示。对于每个旋转的图像，滑动窗依次扫描整个图像。由于某些道砟颗粒倾斜放置，滑动窗法可能无法检测到这些道砟。但是如果旋转图像，这些道砟将转换为竖直方向上，易于识别。将图像旋转15°、30°、45°、60°、75°、90°、105°、120°、135°、150°、165°后，可以获得13张旋转的图像。图4给出了0°、30°、60°的简单示例。

(a)0°

合并这些图像的所有结果。在每个图像中，假如成功识别出道砟，则添加一个道砟颗粒的边界框。然后，在合并时所有边界框会被叠加，但由于旋转，同一道砟会生成一些不同的边界框。此时，将图像旋转到原点0°位置，这相当于只旋转了滑动窗。接着选择面积最小的边界框作为此道砟的最终边界框，如图5所示。

(a)所有旋转合并结果

2.3 Adaboost和Cascade训练分类器

为了比较数据集中的正负图像，将所有标记的图像用来训练模型，得到正负图像之间的特征差异。首先，使用Adaboost和Cascade方法[14]。然后通过滑动窗口方法[15]，即可成功识别道砟颗粒。

Adaboost是一种分类算法，可以将一些弱分类器整合为一个强分类器。而且该强分类器可用于Cascade方法。对于m个正像(yi=1)和n个负像(yi=0)，具有m+n个实例(xi，yi)的训练数据，其中xi表示图像，y表示图像的正负。对于某些给定的值T，该算法迭代t=1，2，…，T次。对于数据实例的权重wi，将每个实例初始化为wi=1。在每次迭代中，将使用相同的分类方法(如决策树)，但数据集的权重会在每次迭代后会发生变化，增加所有错误分类实例的权重，因此，在下一次迭代中会更注意这些错误分类的点。在第t次迭代中，使用当前权重wi构建函数ft并在训练数据集上计算误差errt，如果errt=0(即该方法在训练数据集上非常准确)或errt≥0.5，算法终止；否则，对于每个正确分类的实例，本文保持其权重不变，而对于每个错误分类的实例，将权重乘一个因子(1-errt)/errt。最后，返回分类器的线性组合，其中第t个分类器的输出权重为log(1-errt/errt)。将分类器用于测试一个实例，每一个分类器ft将输出2类的概率分布。将这些概率按相关乘数的比例进行组合。详细算法如下所示。

图6 Adaboost算法

Cascade方法的训练过程分为多个阶段。在每个阶段，通过Adaboost算法训练一个强分类器。此外，需要预先给定3个参数：阶段个数s，检测率d和假阳性率f。其中d表示正确检测的最小百分比，f表示错误识别为正的负图像的最大百分比。每个阶段中，需在Adaboost算法中搜索最小的T，以使由Adaboost训练的强分类器C(x)的检测率大于d且假阳性率小于f。

在如图7所示的假定示例中，将100个正图像和100个负图像示例作为输入数据。每个阶段的原始参数为s=3，d=0.99，f=0.25。然后，在每个阶段中，通过Adaboost训练一个具有适当T的强分类器，并分别假设T为3、10、35。对于第一阶段，模型有200个输入图像，且分类器满足检测率大于d且假阳性率小于f的要求。由于参数d=0.99，因此至少100×0.99=99个图像被正确分类为正，而f=0.25，因此至多100×0.99=99个图像被正确分类为正，而f=0.25，因此至多100×0.25=25被错误分类。因此，分类器将输出124个正图像，其中99个是真实颗粒，将124个图像作为阶段二的输入。对于输出为负的76个图像，则直接拒绝。然后对第二阶段进行类似的计算，由于只有99张正图像，因此检测率为98/99=0.99，假阳性率为6/25=0.24，可以满足要求。然后，在阶段三中，可以得到98个图像的输出，其中包含97个颗粒和1个非颗粒，检测率和假阳性率分别为97/98=0.99和1/6=0.17。因此，在s=3个阶段之后，可得到总检测率D=d3=0.97和假阳性率F=f3=0.015。

图7 Cascade算法的一个假定示例

2.4 基于滑动窗口方法识别颗粒

通过之前得到的分类器，可以获得包含真实颗粒的图像。接下来需要识别完整轮廓道砟颗粒。给定一个铁路道砟颗粒新图像，使用滑动窗口技术，在每个扫描位置，将滑动窗口内的区域输入到检测器中。检测器将在窗口区域内找到并提取特征。基于这些特征，分类器最终确定该图像是正图像还是负图像。如果为正图像，则在扫描位置添加一个边界框，从而将其标识为包含完整轮廓颗粒的区域。具体说明如下。

由于颗粒的最小像素约为50×50，因此，首先使用50×50的滑动窗口扫描整个图像。在每个窗口中，将裁剪图输入分类器。由于先前的分类器都是在32×32图像的基础上训练的，因此，大于32×32的裁剪图像应首先缩小为32×32，并利用Cascade算法中的C1(x)分类器进行测试。如果输出为正，则进入Cascade方法的阶段2。否则，将拒绝负图像，且滑动窗口将以1个像素移动到下一个位置。如果裁剪图像成功通过所有s个分类器，则它将有很大概率是包含完整轮廓道砟颗粒的正图像。然后，将在此位置的滑动窗口标记为该颗粒的边界框。接下来，继续移动滑动窗口直到扫描完整个图像。滑动窗口在水平方向和垂直方向上的移动增量均为1个像素，并且移动路径如图8(a)所示。

当50×50滑动窗口成功扫描整个图像时，在横纵两个方向上将滑动窗口的大小加5个像素，如图8(b)所示。并重复之前的滑动过程，直到滑动窗口扩展到最大尺寸300×300像素为止。由于大多数裁剪后的图像都会被分类为负图像并被拒绝，因此只有极少数图像可以成功通过所有s个分类器，故该算法计算效率较高。

(a)最小滑动窗口

当滑动窗口足够大时，滑动窗口的几次移动都可能包含同一个颗粒，从而产生多个边界框。在这种情况下，取平均值作为一个新的单个窗口，如图9所示。

(a)颗粒的多个边界框

3 形状分析和颗粒库的建立

在分类器识别出完整轮廓颗粒后，本文通过一系列计算几何算法来计算道砟颗粒形状参数。主要通过3个几何形状评价指标来描述颗粒的形状：细长度、棱角圆度和粗糙度，如图10所示。

图10 道砟颗粒形状示例

3.1 几何形状的评价指标

图10展示了道砟颗粒的几何形状评价指标。对于理想的圆形或等边多边形道砟颗粒，其长宽比等于1。对于细长度(EI)，其式如下：

EI=S/L

(1)

式中：S为最小主轴方向上的颗粒直径；L为最大主轴方向上的颗粒直径。

此外，棱角圆度(Rd)描述了颗粒表面棱角处接近圆的程度。

(2)

式中：Rinsc是最大内切圆的半径；Rcirc是位于质心的外接圆的半径。

粗糙度(Rg)反映了叠加在颗粒表面拐角处的形状不规则性。粗糙度计算如下。首先，需要分解颗粒轮廓。根据极性系统中的角度θ和极性直径d，将二维颗粒扩展为与θ和d相对应的连续点。然后，使用“局部加权回归平滑”(LOESS)方法将这些离散点转换为平滑曲线/曲面，得到以下式(3)：

(3)

式中：n是二维颗粒轮廓上的点数。yi是点i的y轴坐标。yi-LOESS是平滑后i的y轴坐标。

3.2 统计形状分析

首先，分析在颗粒库中获得的二维颗粒的形状指标信息。每个形状指标的统计分布如图11所示。

由图11可见，棱角圆度和粗糙度的分布相对狭窄，近似于偏态分布；而细长度的分布较广，近似于正态分布。

(a)棱角圆度的分布频率

3.3 建立颗粒库

为了便于检索和制样，将道砟颗粒的形状指标与坐标一起存储，存储格式为[数字，轮廓，细长度，棱角圆度，粗糙度]。

从颗粒库进行制样的步骤如下：①输入所需形状参数的范围；②然后，根据这些值的范围，可以自动搜索出符合要求的二维轮廓；③选择可视化的二维轮廓并将其输出到“dxf”文件中。

颗粒库建立完成后，选择满足特定要求的道砟颗粒二维轮廓将非常方便，以方便进行铁路道砟的离散元建模，从而为后续研究道砟颗粒的几何信息与力学性能之间的关系奠定基础。

4 双轴剪切行为仿真

4.1 双轴剪切试验的数值模拟

首先，构建5个具有不同细长度的试样(EI=0.3、0.45、0.6、0.75、0.9)。试样的初始大小为9 m×18 m。每个试样包含大约5 000个颗粒。本文利用线弹性模型表示颗粒之间的接触法则。通过移动试样的4个刚性边界(见图12)并确保边界上的力恒定在100 kPa(数值伺服控制的应力边界条件)来压缩试样。为了使试样在剪切过程中处于准静态(为了使试样的加载过程为准静态加载)，剪切应变率应该非常小。因此，引入惯性参量如下：

(a)初始试样

(4)

4.2 双轴剪切试验结果

4.2.1应力比随轴向应变的变化

为了研究试样的剪切强度，本文计算了试样的平均应力p和偏应力q。在二维双轴剪切试验中，有效平均应力p'和偏应力q定义如下：

p'=(σ1+σ2)/2

(5)

q=σ1-σ2

(6)

图 13(a)展示了在具有不同EI值的试样中，应力比q/p'随着轴向应变ε1的变化。初始阶段随着轴向应变的增加，所有试样的应力比迅速增加，然后达到峰值阶段，之后应力比逐渐减小并趋于稳定。当剪切达到稳定阶段时，试样的应力比与EI呈负相关关系。

4.2.2体积应变随轴向应变的变化

此外，本文研究了剪切过程中不同试样的体积变化。考虑到本次模拟的边界是刚性墙，可以根据刚性边界的位移近似计算轴向应变和体积应变。

ε1=(h0-h)/h0

εv=(v0-v)/v0

(7)

式中：h0和h为试样的初始高度和当前时刻的高度；v0和v为试样的初始体积和当前时刻的体积；负体积应变表示体积膨胀。为了研究稳定后试样的临界状态特性，将所有试样剪切到轴向应变ε1=30%。在这种形变下，基本满足临界状态的典型条件(即应力和体积随应变保持恒定)。

如图 13(b)所示，在初始阶段，体积应变随着轴向应变的增加而减小，而当轴向应变达到大约1%时，体积应变随着轴向应变的增加而增加。当轴向应变在20%到30%的范围内达到临界状态时，试样的体积变形保持恒定，这表明所有试样均表现出由剪切诱发的剪胀变形并伴随应变软化。此外，当试样达到稳定状态时，EI值越小(EI=0.3、0.45、0.6)，体积应变越大，而EI值越大(EI=0.75，0.9)，体积应变越小。因此，EI值较小的试样的剪胀性越强。

4.2.3配位数随轴向应变的变化

在模拟的每个时间步长中，计算每个颗粒的平均配位数(MCN)。图 13(c)展示了具有不同EI值的试样的平均配位数随轴向应变的变化。在剪切初期，所有试样的平均配位数随轴向应变增加。当轴向应变大于2%时，平均配位数达到峰值并开始降低。当轴向应变大于10%时，平均配位数基本稳定。每个试样的平均配位数达到稳定值后，随着EI值从0.3增加到0.9，平均配位数从8.5减少到5.2，可以发现其与EI值呈负相关关系。

4.2.4滑动接触的百分比随轴向应变的变化

滑动接触的百分比遵循Mohr-Coulomb定律。滑动系数定义为：

(8)

SP=NSC/NC×100%

(9)

式中：NSC表示颗粒系统中滑动接触的数量，NC是接触的总数。

图 13(d)展示了每个试样的SP随轴向应变的变化。所有试样的SP随轴向应变的变化基本相同。随着ε1的增加，SP迅速增加到峰值，然后逐渐下降到稳定值。

图 13(d)还表明，接触滑动率随EI值的增加而降低。当EI值从0.3增加到0.9时，SP的峰值从0.42减小到0.2左右，SP的稳定性值从18%减小到5%左右。

(a)偏应力随轴向应力的变化

4.2.5EI的摩擦角变化

本节计算了峰值状态和临界状态内摩擦角Φ。非黏性颗粒材料中的Φ定义为：

(10)

式中：σ1和σ2分别为主要和次要主应力。

如图 14所示，不同EI值的试样的剪切强度有明显的变化趋势。在临界状态下，摩擦角与EI值负相关。当EI值从0.3增加到0.9时，临界摩擦角从0.54减小到0.35左右。此外，峰值摩擦角首先随EI值的增大而增大，然后减小。

图14 摩擦角随细长度的变化

5 结论

本文使用模式识别方法获得了40 644个铁路颗粒道砟的二维图像，计算并存储了各颗粒的3个形状指标(细长度，棱角圆度和粗糙度)，建立了二维铁路道砟颗粒模型库，进行了形状指标的统计分析。然后，基于所建立的模型库，重建了二维道砟颗粒的虚拟试样。最后，研究了二维道砟颗粒的双轴剪切特性。结果表明，细长度对道砟颗粒的力学行为具有重要影响。

本文的主要创新如下：

a.提出了模式识别方法，全面获取道砟颗粒的二维投影的信息。与传统的方法相比，基于模式识别的方法具有简单、经济等优势。

b.建立了二维道砟颗粒的形态指标模型库，并通过统计方法进行了统计分析。

c.基于离散元算法研究二维道砟颗粒的力学性能。分析了细长度对剪切性能的影响，为进一步研究形状指标与道砟颗粒的力学性能之间的定性和定量关系奠定了基础。