许汪歆, 袁天辰, 杨 俭

(上海工程技术大学 城市轨道交通学院, 上海 201620)

我国铁路运营里程位于世界前列.随着运营时间的增加,轨道故障时有发生.轨道结构在列车动载荷和环境温度等因素反复作用下,空吊等病害日益突显.目前轨道系统多采用安全检测车巡检并辅以定期人工检修,但这些检修方法可能存在漏检,将严重威胁我国高铁线路的运营安全.因此,如何有效地对轨道结构进行故障诊断和预警是保证铁路运行稳定的重要研究方向.

近年来,基于智能算法的结构健康监测系统已逐渐应用到机械设备的状态监测与故障诊断中[1].研究人员结合振动信号冗余与数据稀疏化理论做了一些探索.温江涛等[2]采用信号稀疏表示理论,实现信号稀疏表示.刘翠响等[3-4]针对稀疏字典较难重构的问题,使用奇异值分解获得过完备字典,实现信号的稀疏化表示.但现有的故障诊断方法大多基于信号处理的方法[5],侧重结合监督学习方法对机械设备进行故障诊断[6-7],没有将信号稀疏表示与故障诊断结合.

基于无监督学习的机械设备故障诊断是对数据进行特征提取[8].原始信号虽然能够最完整地表达机械设备(轨道结构)的全部信息,但在一定程度上存在冗余,不能较好地实现特征提取.而轨道结构振动过程中产生的数据集也存在样本类标签未知、各数据集分布不平衡、非线性、多工况、瞬态等诸多问题.针对上述问题,本文提出一种新的轨道结构智能诊断方法.通过将轨道结构振动信号稀疏化表示,降低振动信号“冗余性”;通过将稀疏化表示的振动信号数据集样本转化为多个特征指标,将特征指标两两组合,实现从不同维度识别轨道结构故障,为轨道结构故障智能诊断提供参考.

1 车辆—轨道耦合动力学计算模型

1.1 轨道结构振动方程

结合车辆—轨道耦合垂向振动结构断面模型,以速度v沿轨道运行,则有下列条件成立.

1) 根据钢轨微段垂向受力平衡,可得钢轨振动方程为

(1)

式中:EI为钢轨抗弯刚度;Zr(x,t)为钢轨振动位移变量;Frsi(t)为轨枕支点反力;xi为各轨枕支点坐标,xi=ils(i=1～N);p(t)为单侧车轮轮轨垂向作用力.

2) 根据模型中第i号轨枕受力状态,可得其振动方程为

(2)

式中:Msi为单位长度轨枕质量;Cpi、Cbi分别为轨下垫层、道床阻尼;Kpi、Kbi分别为轨下垫层、道床刚度;Yh(xi)为振型.

3) 根据第i号离散道床块受力状态,可得其平衡振动方程为

KbiZsi(t)-CwiZb(i+1)(t)-

KwiZb(i-1)(t)=0

(3)

图1 车辆—轨道耦合垂向振动结构断面模型Fig.1 Structural section model of vehicle-track coupling vertical vibration

式中:Mbi为单位长度道床参振质量;Cfi、Cwi分别为轨下路基及道床剪切阻尼;Kfi、Kwi分别为轨下路基及道床剪切刚度.

轨下基础结构服役状态对轨道结构动力性能影响的仿真工况见表1.表中,故障设置轨枕空吊、道床翻浆、道床板结以及正常4种工况,列车运行速度分别为200、150和100 km/h,故障点设置在第100号轨枕处,进行仿真计算.

表1 轨枕出现不同故障时数值模拟计算工况Table 1 Numerical simulation calculation conditions of different faults for track sleeper

1.2 数值仿真结果及分析

仿真计算中,改变轨道结构的质量、刚度和阻尼矩阵以及列车运行速度、轨道不平顺等参数,可以模拟研究不同轨道工况下轨道结构振动响应的变化规律.在分析车辆通过轨道故障结构路段的耦合振动时,选取轨枕加速度作为动力学响应指标分析.为模拟轨道结构振动响应工况,设定采样频率为10 kHz,轨道振动响应时间历程为0.8～1.7 s,采样点数为N=9 000,满足振动响应分析要求,如图2所示.图3为不同行车速度下轨道结构故障处轨枕加速度峰值直方图.

图2 100号轨枕振动加速度时程曲线Fig.2 Curve of vibration acceleration with time history for No.100 track sleeper

图3 97～103号轨枕对应不同故障加速度峰值Fig.3 Acceleration peeks of No.97 to 103 track sleeper for different faults

根据图3的仿真结果,不同工况下轨枕加速度峰值明显不同.当轮轨冲击作用力沿着钢轨传递到轨枕,轨枕空吊的加速度产生较大变化.轨枕空吊引起的轨枕加速度增加表明轨道结构故障区段的轨枕间距和轨下支承刚度已经发生改变.

由于轮轨作用力随机性和周期性并存,并且车轨耦合动力学具有非线性特点,上述试验结果具有一定的离散性,仅从动力学角度无法较为全面地诊断出轨道结构故障.因此,本文提出基于密度聚类算法的振动信号处理和特征指标提取的故障诊断方法,对轨道结构进行故障分类.

2 轨道振动信号稀疏表示

稀疏表示是由Coifman和Wickerhauser于1992年提出的概念[10].稀疏表示是压缩感知与信号重构的基础,因此有必要研究信号的稀疏表示.轨道振动信号稀疏表示可以有效地解决低功耗无线传输以及压缩感知信号重构在机械设备故障诊断领域中的应用.然而,轨道振动时变信号并不具备稀疏性,因此直接将原始振动数据通过低功耗无线传感节点传输是不可行的.研究信号稀疏表示问题的本质是研究对稀疏字典ψ的设计.利用字典冗余特性可以更好地捕捉信号本质特征.而绝大多数时变信号可以在某个特定字典下稀疏表示,因此在轨道结构振动信号传输之前,可以将其稀疏化.对稀疏字典的研究,主要包括完备字典及基于字典学习理论的过完备字典[11].基于挑选的字典,利用稀疏表示对轨道结构振动信号进行建模是可行的.寻找稀疏表示下最优基的构造涉及较多理论,本文选择标准正交基矩阵作为轨道振动信号稀疏表示的字典,即正交基的分解是唯一的.信号稀疏表示流程如图4所示.

图4 稀疏表示理论Fig.4 Sparse representation theory

2.1 基于完备字典的信号稀疏表示

稀疏表示字典按照矩阵ψ中元素个数不同分为完备字典和过完备字典.字典中元素按照关键信息进行编码,完备字典中的标准正交基矩阵具有明确的解析表达,本文选取离散余弦变换(Discrete Cosine Transform,DCT)以及离散傅里叶变换(Discrete Fourier Transform,DFT)[12].

对于N点一维序列{x(n):0,1,…,N-1},基于离散余弦变换的信号稀疏表达式为

(4)

其中,k=0,1,2,…,N-1.

将上式单位化后,得到矩阵X=Cx,其中,时域信号x=[x(0),x(1),…,x(N-1)]T;频域向量X=[X(0),X(1),…,X(N-1)]T;标准正交基矩阵C第i行第j列元素可表示为

(5)

其中, 1≤i,j≤N.

因此可以使用CT作为离散余弦变换的字典矩阵.基于离散傅里叶变换的信号稀疏表达式为

(6)

2.2 完备字典对轨道结构振动信号稀疏表示性能分析

根据表1轨道结构工况,结合上述正交基矩阵以及轨道结构振动响应,对99号轨枕时域振动信号稀疏分解性能进行分析,行车速度为150 km/h.通过仿真分析,得到轨道结构振动信号在两种正交稀疏化字典下对振动数据的信号稀疏效果,如图5所示.200和100 km/h行车速度下,轨道结构振动信号稀疏效果与图5一致.表2给出3种行车速度及4种轨道结构状态下稀疏系数与最佳稀疏系数. 由表2可知,4种轨道结构状态都存在稀疏系数,可知轨道振动信号存在冗余性,而事实上,信号稀疏表示在理论上是欠定线性系统,最佳稀疏系数的求解具有不确定性.因此,有必要在今后的研究工作中进一步研究信号稀疏的最优化问题.

图5 150 km/h速度下两种正交稀疏字典的稀疏效果Fig.5 Sparse effect of two orthogonal sparse dictionaries under speed of 150 km/h

表2 轨道结构在不同工况下的稀疏系数与最佳稀疏系数Table 2 Sparse coefficient and the optimal sparse coefficient of track structure under different working condition

2.3 轨道振动信号特征指标提取

由于原始振动信号具有冗余性,直接采用原信号进行特征提取不能较好地实现密度聚类的轨道结构故障分类.根据轨道结构振动信号x,将N=9 000个振动信号进行稀疏化表示,取不同轨道结构、不同位置的振动信号,转换成振动信号特征指标,用以描述轨道结构的服役状态,进而对轨道结构进行故障诊断.根据文献[13-14]对所选取特征指标定义,选取能够体现指标绝对水平的平均值、均方根值、方根幅值3个有量纲指标以及体现指标概率密度分布的波形指标、峰值指标、脉冲指标、裕度指标、偏斜度指标、峭度指标等6个无量纲特征参数,每一个特征指标代表一个维度,特征指标两两组合,可以实现从不同维度观测轨道结构的故障表现形式.以150 km/h行车速度为例,提取的特征指标样本能较好地表示轨道结构故障,见表3.

表3 150 km/h速度下轨枕不同故障下指标间部分样本分布(N=9 000)Table 3 Partial sample distribution among indexes under different faults of sleeper under speed of 150 km/h

3 基于密度聚类算法的轨道结构故障诊断

具有噪声的基于空间密度的聚类(DBSCAN)算法是通过数据样本分布的紧密程度确定类别划分,从数据样本密度的角度来考察样本之间的可连续性,并基于可连接样本不断扩展聚类类别以获得最终聚类效果.通过ε,min Pts来刻画样本分布的紧密程度.密度聚类算法对聚类后的类别形状没有限制,因此可以描述任意形状类别.而类别间距离刻画了不同类别间属性,通过对轨道结构振动数据特征形成、特征提取、特征选择和故障诊断,可识别出不同轨道结构故障,表征轨道结构服役状态.

3.1 k-距离图、ε、min Pts确定

DBSCAN要选择一种距离度量,在二维平面讨论数据分区问题.给定特征参数集P={p(i);i=0,1,…,n},对于任意点p(i),计算每个点与其他所有点之间的欧几里得距离,按照降序排序,得到k-距离图,如图6所示.通过输入指标数据,按DBSCAN启发式方法寻找k-距离图排序的凹陷点(阈值点),确定每个分区的邻域半径ε值.min Pts为以对象p为核心对象时其ε邻域中包含对象最小数目.

(7)

式中:dij为数据集中所有距离;dc为一个截断距离,由dij进行降序排列得到.对于每一个数据节点i,都能找到所有比该节点密度大的节点j,选取其中最小的dij,记为δi,则

图6 不同工况下k-距离图Fig.6 k-distance diagram under different working conditions

(8)

根据式(5)和式(6),本文选取min Pts=4.以列车运行速度150 km/h为例,特征参数点密集处属于同类密度半径,反映在k-距离图中为平缓曲线,如d、e、f,h段急剧变化,因此确定ε=0.4[15].

3.2 基于DBSCAN的轨道结构故障诊断

DBSCAN通过形成具有“同质”数据之间的类别,在轨道结构数据输入时,算法会判断每一个数据点是否满足核心对象条件,进而判断是否有必要形成新的类别[16].根据表3中150 km/h工况下特征指标选取,对轨道结构故障的4种工况聚类,进行故障诊断,仿真结果如图7所示.

图中虚线椭圆表示轨道结构聚类所属类别.不同类别间距离用d1,d2表示.特征指标代表维度,不同特征指标之间两两组合,可以实现从不同维度对聚类结果观测,进而判别出轨道结构故障.从图中可以看出,图7(a)完全可以识别出4个类别,对于正常工况与道床板结而言,在图7(b)局部放大图中可以看出这两者之间是可以区分的.图7(c)与7(d)比较可见,道床翻浆与轨枕空吊在 不同特征指标中平面分布位置不同.正常与道床板结工况在图7(e)能区分开.以上充分说明,根据不同特征指标之间的组合,从不同维度观测,可以判断出轨道结构故障类别划分.

图7 150 km/h速度时不同特征指标聚类结果Fig.7 Clustering results of different characteristics indexes at a speed of 150 km/h

3.3 聚类有效性评价指标

基于轨道结构振动数据集D={x1,x2,…,xn},通过聚类给出类别划分C={C1,C2,…,Cm},则聚类有效性评价指标可以表示为

(9)

(10)

dcen(Ci,Cj)=dist(μi,μj)

(11)

基于密度聚类的轨道结构聚类有效性评价指标函数,道床板结与正常工况这两个类别的dmin(Ci,Cj)与dcen(Ci,Cj)几乎相等,这也进一步表明,道床板结与正常工况这两类几乎属于同一个类别.对于道床翻浆与轨枕空吊,无论是在何种特征指标下,dmin(Ci,Cj)与dcen(Ci,Cj)数值不相等,这表明这两个类别类间分离性较大.

4 结 论

本文研究密度聚类算法下轨道结构故障诊断问题.通过将原始振动信号进行稀疏表示、特征提取,在密度聚类算法下可识别出轨道结构故障,为轨道结构故障诊断无损化、智能化提供理论参考.