施 莹， 庄 哲， 林建辉

(西南交通大学牵引动力国家重点实验室 成都, 610031)

引 言

轮对轴承是高速列车走行部的重要部件，它既承担了由列车自重及负载合成的垂向力，还承担了轮轨间横向力。在高速列车长期运行过程中，这些外力的联合作用使轮对轴承有可能发生如剥落、点蚀和烧损等故障，这些故障一旦出现，极易在短时间内恶化，最终危及高速列车的运行安全。因此，研究高速列车轮对轴承的故障诊断方法对保障列车安全运行具有重要意义。

轮对轴承故障检测的方法主要包括油脂(液)监测、声发射监测、轴温监测及振动监测[1]。然而，油脂(液)的监测方法不适合轴承的实时故障监测；声发射技术因轴承故障引起的瞬态弹性波衰减现象使其无法得到广泛使用；轴温监测是目前较常用的检测手段：当轴箱温度超过事先设定的阈值时，轴温监测装置发出警报，但该方法只对严重故障敏感，无法诊断早期故障[2-3]。轴箱振动信号不但易于监测，还包含丰富的轴承状态信息。因此，振动监测将成为轮对轴承故障，特别是早期故障的有效监测手段。

轴承故障引起的脉冲常与车轮踏面缺陷或轮对不圆时引起的脉冲混杂在一起，从这些非线性和非平稳的具备调制特性、强噪声干扰的信号中提取有用故障特征一直是一个挑战[4]。

为了有效地应对挑战，众多学者提出了许多先进的信号处理技术，包括经验模态分解[5](empirical mode decomposition,简称EMD)，小波变换(wavelet transform,简称WT)，并成功地应用于旋转机械的故障诊断中[6-7]。EMD非常适合于分析非线性非平稳信号，并能将信号分解为一组本征模函数(intrinsic mode function,简称IMF)和残余分量，已成功地应用于故障检测、诊断等领域的研究。EMD的缺点(如缺乏理论基础，噪声的敏感性，模态混叠等问题[8])限制了其进一步广泛应用。为了解决这些问题，聚合经验模态分解[9](ensemble empirical mode decomposition,简称EEMD)等算法被陆续提出，但以上问题并没有从根本上得到解决。WT是另外一种常用的时频分析方法。相较于离散小波变换(discrete wavelet transformation,简称DWT)的快速算法，连续小波变换(continuous wavelet transform,简称CWT)巨大的计算成本和冗余系数给实际应用带来了巨大困难，小波包变换(wavelet packet transform,简称WPT)的提出使信号的高频部分分解性能得到了较大提升，然而DWT和WPT的分解效果受小波基函数影响较大。

基于交替方向乘子法[10](alternating direction method of multipliers,简称ADMM)的卷积稀疏表示于2015年提出，该方法为轴承故障诊断的特征提取提供了良好前景。然而，CSR提取脉冲的性能受到其惩罚因子的影响较大：当惩罚因子过大时，提取的脉冲小于真正的脉冲数；当惩罚因子过小时，提取的脉冲大于真正的脉冲数，这便会给故障的精确诊断带来极大困难。

流形学习是一种非线性降维方法，目前广泛应用的流形学习算法一般分为两种:基于全局的流形学习分析方法以及基于局部的流形学习分析方法。基于全局的分析方法主要有等距映射[11]和最大方差展开算法[12](maximum variance unfolding,简称MVU)等。基于局部的分析方法主要有局部切空间分析[13](local tangent space alignment,简称LTSA)、局部线性嵌入[14](locally linear embedding,简称LLE)和拉普拉斯特征映射算法[15](Laplacian eigenmaps,简称LE)等。

笔者提出了一种基于卷积稀疏表示、希尔伯特变换和流形学习的新的故障检测方法。首先，通过在不同惩罚因子下的卷积稀疏表示提取不同稀疏特征的脉冲；然后，针对提取的一系列脉冲进行希尔伯特变换，构造脉冲包络空间；最后，利用等距映射流形学习算法对脉冲包络空间求解低维本征包络，以实现故障诊断。

1 基本原理

1.1 卷积稀疏表示脉冲提取

传统的信号稀疏表示理论基于正交线性变换，通过将信号展开为一系列同形态基函数线性加权和的形式。然而，在工程实际中大部分信号都不是单一形态的成分，因而很难用单一的形态基函数稀疏表示。因此，就需要利用两种或两种以上的不同形态基函数对信号加以表示，其基本思想是用超完备的冗余基函数原子库取代正交的基函数库，原子库可以由任意基函数构成，从而最优地匹配原始信号的结构。从原子库中寻求信号的最佳线性原子组合来表示信号，被称为信号的稀疏表示。

卷积稀疏表示理论将信号表示为一系列原子与稀疏系数的卷积和，其定义[10]为

(1)

其中：sk∈Rn为第k项待分析的振动信号；n为信号长度；dm∈Rp为一系列原子；p为原子长度；xk,m∈Rn-p+1为与给定信号sk和原子dm相关的稀疏系数；λ∈R+为惩罚因子；K为给定信号sk的个数；M为原子dm的个数。

当给定一组信号和原子，式(1)中的稀疏系数可由平移不变稀疏编码[10](shift invariant sparse coding,简称SISC)得到，该方法将信号表示为多个原子与系数的卷积形式，自适应地学习原子与系数达到对原信号的稀疏表达，其表达式为

(2)

式(2)中，不同稀疏系数xk,m针于不同信号sk进行解耦，稀疏系数xm表示为

(3)

当信号sk及稀疏系数xk,m已知，式(1)中的原子dm可由平移不变性字典学习[10](shift invariant dictionary learning,简称SIDL)得到，该算法为变量交替更新算法，即轮流更新x和d，表达式为

(4)

所以，dm=1∀m。

(5)

1.2 希尔伯特脉冲包络提取

关于轴承故障的特征信息包含在提取的脉冲的包络中，提取的脉冲的希尔伯特变换定义[16]为

(6)

(7)

(8)

(9)

1.3 脉冲包络空间流形降维

流形学习将式(9)所示的高维脉冲包络空间IESm映射至一个低维本征脉冲包络中。文中采用等距映射流形算法进行降维计算。

基于等距映射流形方法是一种非迭代、非线性的数据降维算法，其基本原理是：将任意两个数据点之间的测地线距离(也就是两点间的最短距离)当成流形的几何表述，通过多维缩放算法(multiple dimensional scaling,简称MDS)算法保持这个点与点之间的最短距离，形成局部测地距离和低维非线性特征空间的映射关系，进而获取最佳的低维流形结果，还可以满足数据测地距离在低维流形空间维持不变的特点。

等距映射流形算法的优化函数[11]如下

(10)

τ(DG)-τ(Dγ)L2

(11)

那么Isomap的最优化解Y=(y1,y2,…,ym)可以通过τ(DG)的最大特征值表示的向量得到。

1.4 估测原子类型数量

在式(1)中，原子的数目是不确定的，而是基于轴承故障振动特性来确定[17]。脉冲可以用共振频率和阻尼参数两个参数来完全描述，共振频率可以定义不同的原子，因此利用原子的主频率来估计原子数，考虑到轴承的故障特性，在轮对轴承故障检测中，一般预先设定为4，然后通过CSR的迭代得到M的准确值。每个原子的主频率通过傅里叶变换提取fm(m=1,2,…,M)。如果任何两个主频率的差大于频率分辨率fsp-1(fs表示采样频率)，则表明有原子；否则M=M-1，循环执行CSR并且计算主频率直到任何两个主频率的差大于频率分辨率。输出M的终值作为原子数目估测值。

1.5 基于卷积稀疏表示及流形学习故障诊断总体框架

基于卷积稀疏表示及流形学习的轴承故障诊断流程如图1所示。

图1 基于卷积稀疏表示及流形学习的故障诊断流程Fig.1 Flowchart of the fault detection based on CSR and the manifold

具体步骤包括：

1) 准备分析信号；

2) 估测使用主频分析的原子类型数目；

3) 根据式(1)～(5)提取不同稀疏特性的脉冲；

4) 根据式(9)构建脉冲包络空间；

5) 针对脉冲包络空间展开流形降维，求解低维本征脉冲包络；

6) 计算本征10包络的包络谱判断轴承故障。

2 仿真验证

为了证明所提方法的有效性，设置如下仿真信号

(12)

其中：L为模拟脉冲数目，为39；Al为l的故障脉冲振幅，为5×10-9m；Tp为故障特征频率fp的倒数；fp设为49.1 Hz；β为结构阻尼系数，为1 200 Ns/m；ω为固有频率，为2 000 Hz。

图2表示仿真信号原始时域波形及信噪比为-7 dB的含白噪声信号波形。

图2 模拟信号Fig.2 Simulation signals

图3(a)表示基于CSR-HT-Isomap方法提取的低维本征脉冲包络，其中，惩罚因子λ=6∶1∶15，Isomap的近邻点个数[11]取10，可以发现，该方法成功地提取了全部39条脉冲包络。图3(b)表示该低维本征脉冲包络的包络谱，可以发现该方法将故障特征频率fp及其20次谐波均成功提取，证明了笔者所提方法的有效性。

图3 CSR-HT-Isomap故障特征提取Fig.3 Fault feature extraction by the CSR-HT-Isomap method

现比较所提方法与传统的基于EEMD和WPT包络空间的Isomap特征提取效果，其中，EEMD中添加的白噪声的幅值选择为0.2σ，分解次数为400。WPT中母小波为Db8小波。以IMF和多尺度包络空间为基础的本征脉冲包络如图4所示。IMF包络空间发现10次谐波的故障特征频率，多尺度包络空间发现5次谐波，而且通过基于IMF和多尺度包络空间的包络幅值和包络谱幅值均小于图3。这些表明，笔者所提方法要优于IMF包络和多尺度包络空间。

图4 本征包络流形提取Fig.4 The intrinsic envelope extracted by the manifold with Isomap and their envelope spectra

3 台架实验验证

为了获得高速列车轮对轴承振动信号，采用高速列车轮对跑合实验台进行台架实验，如图5所示。跑合实验台由电机、驱动轮、加载装置、轮对及轴箱组成，实验中，对轮对轴承做了人工伤处理，设置了外圈故障与滚柱故障，外圈故障的深度、宽度皆为1 mm，长度为5 mm，滚柱故障的深度、宽度皆为1 mm。如图6所示，加速度传感器安置于轴箱上，轴承参数见表1。

图5 轮对跑合实验台及加速度传感器Fig.5 The test bench and measurement sensor

图6 高速列车轮对轴承人工伤Fig.6 Man-made fault photos of the wheelset of high speed train

滚子数量滚子直径/mm节圆直径/mm接触角/rad转速/(r·s-1)1926.91800.157 120.57

3.1 外圈故障检测

当轮对轴承外圈存在图6(a)所示故障时，轴箱振动信号如图7所示。

图7 外圈故障轴箱振动信号Fig.7 Vibration signals of bearing outer race fault

外圈的故障特征频率fBPFO表示为

(13)

其中:Nb为滚子数量;Bd为滚子直径;Pd为节圆直径;φ为接触角;fw为转速；经计算fBPFO=166.6 Hz。

利用CSR-HT-Isomap方法提取的本征脉冲包络如图8所示，其中图8(b)成功地发现了外圈的故障特征频率及其7阶谐波。

图8 CSR-HT-Isomap故障特征提取Fig.8 Fault feature extraction by the CSR-HT-Isomap method

现比较笔者所提的基于CSR-HT-Isomap方法与基于传统的EEMD和WPT包络空间的Isomap特征提取效果，其中，EEMD中添加的白噪声的幅值选择为信号标准差，分解次数为400。WPT中母小波为Db8小波。以IMF和多尺度包络空间为基础的本征脉冲包络流形提取如图9所示。IMF包络空间发现了故障信息中7次谐波的故障特征频率，如图9(b)所示，多尺度包络空间仅仅发现2次谐波的故障特征频率，如图9(d)所示，而且通过基于IMF和多尺度包络空间得到的包络谱幅值均小于基于CSR-HT-Isomap方法提取的本征脉冲包络空间包络谱幅值，如图8(b)、图9(b)及图9(d)所示。这些表明在提取由轴承外圈故障特征时，脉冲包络空间的表现要优于IMFs和多尺度包络空间，无论是在放大故障特征频率的谐波数方面，还是在强化本征包络谱方面，前者均优于后两者。

图9 本征包络流形提取Fig.9 The intrinsic envelope extracted by the manifold with Isomap and their envelope spectra

3.2 滚子故障检测

当轮对轴承外圈存在图6(b)所示故障时，轴箱振动信号如图10所示。

图10 滚子故障轴箱振动信号Fig.10 Vibration signals of a bearing roller fault

滚子的故障特征频率fBSF表示为

(14)

计算后，fBSF=67.9 Hz。由CSR-HT-Isomap算法提取的轴承滚子故障振动本征脉冲包络如图11所示。

图11 CSR-HT-Isomap故障特征提取Fig.11 Fault feature extraction by the CSR-HT-Isomap method

根据图12所示的滚子和滚道的运动关系。在轮对轴承滚子负载区(load zone，简称LZ)中有从脉冲B1经E1到H1有7个脉冲。在图11(a)中有7个脉冲在滚子负载区，在轮对轴承滚子非负载区域(none load zone，简称NLZ)，由于在滚子非负载区域缺乏接触力导致作用在滚子故障处的力太小而无法激发脉冲。在脉冲B1到脉冲H1这7个脉冲中，脉冲B1和脉冲H1的位置在滚子负载区的边界，当有缺陷的滚子处于脉冲E1的位置上时，在滚道和有缺陷的滚子之间出现了一个缺口，相当于脉冲F1和D1的两个正常滚子形成的对称结构造成缺口的产生。所以脉冲B1，H1和E1的包络边缘小于脉冲C1，D1，F1和G1的包络边缘。通过CSR-HT-Isomap方法获得的本征包络成功地发现了滚子故障时轴承的振动特点，同时本征包络的包络谱成功地显示了故障特征频率fBSF和它的2次谐波。仔细观察还可发现，包络谱的频率呈现周期性波动的特点，除第一部分频率峰值围绕在转频左右，其他的峰值皆围绕在滚子故障特征频率及其倍频频率附近，这为故障的识别提供了有力的证据。

图12 滚子与滚道间的运动关系Fig.12 The relative motion between the rollers and races

比较基于IMFs包络空间和多尺度包络空间的本征脉冲包络如图13所示。可以发现，其包络空间无法体现滚子故障时滚子处于不同位置时的轴承振动特点，包络谱的幅值较小，且故障特征频率淹没在噪声频率中，无法反应轴承故障状态。

图13 本征包络流形提取Fig.13 The intrinsic envelope extracted by the manifold with Isomap and their envelope spectra:

4 结束语

笔者提出了一种基于CSR-HT-Isomap相结合的故障检测方法，通过仿真数据以及台架实验数据表明：该轮对轴承故障诊断方法可以很好地提取仿真信号、轴承内圈及滚动体故障特征，通过与基于聚合经验模态分解和小波包变换的包络空间算法进行比较，证明本方法在提取本征包络、强化本征包络谱以及放大故障特征频率的谐波数方面具备较大优势，具有一定的工程应用前景。