基于K-奇异值分解和层次化分块正交匹配算法的滚动轴承故障诊断

张文颢 李永健 张卫华

1.西南交通大学牵引动力国家重点实验室，成都 ，610031 2.五邑大学轨道交通学院，江门，529020 3.西华大学汽车测控与安全四川省重点实验室，成都，610039

0 引言

轴箱轴承是列车转向架中的关键部件，同时也是易损部件。轴箱轴承发生故障后直接影响列车的动力学性能，严重时可导致脱轨[1],因此及时发现和诊断轴箱轴承故障具有重要意义。滚动轴承诊断方法主要有振动信号诊断方法、温度诊断方法、声音诊断方法，其中振动信号诊断方法应用最为广泛且易于实时在线监测[2]。目前主要的降噪方法有很多，如本征时间尺度分解[3]、小波阈值降噪[4]、经验模态分解[5]以及集合经验模态分解等。

近年来，基于稀疏表示理论[6]的信号处理方法在图像处理[7]、语音识别[8]、压缩感知[9]等领域得到了广泛应用。由于旋转机械故障的特征成分在整体信号中往往表现出稀疏性，因此有越来越多的学者尝试将稀疏表示理论应用于机械故障诊断领域，并取得了一定的成果。稀疏表示的基本思想是从冗余字典中选取具有最佳线性组合的原子对原信号线性表示。目前，稀疏表示理论的研究主要分为两个方面：冗余字典的构造和稀疏系数的求解。字典的构造方式主要有两种：①基于已有分析字典及其组合的构造方式；②基于字典自学习的构造方式[10]。分析字典虽然具有详细的数学描述，且目前该方法研究较为深入，但是这种基于固定数学公式的字典构造方法存在适应性不足的缺点。字典自学习的构造方法(如K-SVD、最佳方向法等)是近年来较新的研究方向，该方法能够根据原信号自适应地进行训练，从而得到匹配度更高的原子库。对于稀疏系数求解，目前经典的算法为基追[11](basis pursuit，BP)和匹配追踪[12](matching pursuit，MP)。在MP算法的基础上人们又相继提出了正交匹配追踪[13](orthogonal matching pursuit，OMP)、批处理正交匹配追踪[14](batch orthogonal matching pursuit，Batch-OMP)等算法，克服了MP算法存在过匹配的现象。文献[15]在图像处理领域提出了一种基于层次化的分块正交匹配(hierarchized block wise orthogonal matching pursuit，HBW-OMP)算法，与传统的MP算法和OMP算法相比，该方法的表示结果在小波域具有更好的稀疏性。文献[16]在HBW-OMP的基础上优化了原子的选取方法并应用于音频文件压缩处理，取得了良好的效果。本文利用层次化选择的分块正交匹配方法(hierarchized block wise optimized orthogonal matching pursuit, HBW-OOMP)方法的高稀疏性和运算速度快等优点，与K-奇异值分解(K singular value decomposition, K-SVD)字典相结合，对轴承振动信号完成重要特征提取，并进行故障诊断。

1 冗余字典构造

K-SVD是一种字典自学习的训练算法[14]，其核心思想是通过约束问题的稀疏表征和奇异值分解更新算法交替进行，最终得到自适应的冗余字典。相比其他由固定基底扩张出来的冗余字典(如DCT、小波基、Gabor字典等)，K-SVD自适应性更强，方法更为高效。K-SVD字典训练方法主要步骤如下：

(1)字典D的初始化，可直接将训练信号按列赋值到初始字典D或者采用随机数生成。

(2)采用稀疏表示系数求解算法以下式为目标函数进行系数求解：

(1)

式中，xi为原信号向量；ai为稀疏系数向量。

(3)更新字典。对字典D中任意第k列都按下式更新:

(2)

利用奇异值分解对Ek进行更新。

(4)字典更新完成后，重复步骤(2)和步骤(3)，直到满足迭代终止条件。

2 稀疏系数求取

假设采集的原信号为y∈Rn，冗余字典D∈Rn×q，根据稀疏表示理论，原信号y可表示为

(3)

式中，qi为字典D中原子;ci为稀疏表示系数;ε为信号残差。

原信号中冲击成分与冗余字典中原子的结构具有较强相关性，投影系数较大，可由少量原子线性表示；而信号中的噪声成分与原子相关性较弱，投影系数小，故需要大量原子才能表示。所以有效的原子选取方法对稀疏表示质量至关重要。

在OMP算法的基础上，文献[16]提出了一种层次化选择的分块正交匹配方法HBW-OOMP。HBW-OOMP的分块化以及层次化选择方法避免了计算复杂度及存储量的显著增加。相比直接对原信号采用OMP算法进行匹配，该方法运算速度快，表示结果稀疏度更高。HBW-OOMP的主要步骤如下。

(1)对长度为N的原信号进行分块处理得到M个块，每个块的数据长度为d且互不重叠，即

(4)

(2)分别对各个块xT[m]与冗余字典D进行内积运算，找出与各个块对应的内积最大的原子qm:

qm=|>〈D,xT[m]〉|m=1,2，…,M

(5)

(3)根据下式选出匹配度最高的块xT[n]，并对块xT[n]进行操作:

n=arg maxqm

(6)

RxT[n]=xT[n]-〈xT[n],qn〉qn

(7)

cn=|>〈xT[n],qn〉|

(8)

式中,qn为选出的原子。

xT[n]*=RxT[n]

(9)

(10)

(5)重复步骤(3)、步骤(4)，直到选取原子数等于给定值或者信号残差小于某个阈值。

(6)重构信号。

通过对全局稀疏度的限制而不是对每个块的逼近质量的控制，信号的整体表示质量会有很大提升。由于各分块没有重叠部分，所以避免了块效应。HBW-OOMP在每一次迭代中提出了一种层次化的选择块的方法，因此包含冲击成分的块会优先选出。传统方法通过设定阈值进而控制最终选取的原子个数。但是阈值的设定往往需要根据经验和多次的尝试。原子数过少导致重构信号包含较少的故障特征，原子数过多又会使得重构信号存在更多的噪声和冗余特征。本文迭代终止条件不再是信号残差小于某个阈值，而是采用信号包络谱峭度最大准则[10]，即对信号的包络谱进行峭度的求解。由于故障冲击成分存在规律性，反映在包络谱中则表现为峰值的有序性，信号的包络谱峭度越大，重构信号中的随机波动的幅值越低。与设定阈值相比，该方法能够更大程度地保留冲击成分，减少噪声和冗余信息，根据原信号自适应地确定原子数。

3 实验研究

3.1 仿真测试

为了验证本文所提方法的有效性，根据实际振动信号构造如下故障轴承仿真模型：

(11)

其中，第一部分表示由滚动轴承故障引起的冲击信号;Ta为相邻冲击的标准时间间隔；τi为随机变量，表示由于滑移导致两个冲击之间的时间滞后现象，通常为Ta的1%～2% ；Ai用来模拟信号的幅值；Dj为幅值参数；Td为时间参数；sj(t)为冲击响应函数。等式右边第二部分用来仿真由外部敲击以及电磁干扰导致的混入信号，由于混入信号的不可预知性，时间参数Td和幅值参数Dj为随机变量。等式右边第三部分是由叶片、轴等部件产生的离散谐波干扰。n(t)为高斯白噪声。

sj(t)表示如下：

si(t)=e-Bitcos(2πfit+φi)

(12)

其中,fi为冲击引起的共振频率；Bi为阻尼系数；φi为相位。

设定仿真故障频率为135 Hz，仿真原信号见图1。应用本文方法对仿真信号进行稀疏表示。首先对原信号进行K-SVD字典学习：将信号按照一定重叠率分割成300个样本，每个样本长度为32，训练后得到32×300的冗余字典。然后采用HBW-OOMP进行原子的选取以及稀疏系数的求取：将原信号分成长度为32的块，且各个块之间无重叠部分。由于迭代次数越多，重构信号越逼近原信号，导致重构信号中逐渐包含噪声信息，所以限制最大迭代次数为100。每次迭代时重构信号的包络谱峭度见图2，可以看出在原子数为42时包络谱峭度最大，故最终迭代次数为42。重构信号的时域图以及包络谱见图3、图4。本文所提方法运行时间为1.129 s。对同一段原信号采用Gabor-OMP方法进行稀疏表示，终止条件是该方法得到的残差与本文方法得到的残差一致，重构信号见图5。该方法运行时间为23.992 s，选取原子数103。可以看出本文方法在表示结果和运行时间上比传统的Gabor-OMP算法有很大改进。

图1 仿真原信号Fig.1 Simulation of the original signal

图2 不同迭代次数下重构信号的包络谱峭度(仿真信号)Fig.2 The envelop spectrum kurtosis in each step (simulation signal )

图3 本文方法重构信号(仿真信号)Fig.3 The reconstructed signal by the method proposed (simulation signal)

图4 重构信号包络谱(仿真信号)Fig.4 The envelop spectrum for the reconstructed signal (simulation signal)

图5 Gabor-OMP方法重构信号(仿真信号)Fig.5 The reconstructed signal by Gabor-OMP (simulation signal)

3.2 实验研究

实验采用的振动信号来自天马轴承厂试验台所测实验数据。表1给出了所用轴承的基本参数。设置4种轴承状态：正常、外圈故障、内圈故障和滚动体故障。采用电火花加工方法在外圈、内圈和滚动体上加工出宽0.1 mm、深0.43 mm的划痕来模拟故障。

表1 轴承参数Tab.1 The specifications of tested bearings

电机设置转速为1 169 r/min，采样频率5 120 Hz。根据故障特征频率计算公式可得到：外圈故障频率f0=141 Hz,内圈故障频率f1=193.2 Hz,滚动体故障频率f2=133 Hz。

3.2.1外圈故障特征提取

轴承外圈故障时的时域信号见图6。外圈故障状态下，故障位置固定不动，较高的转速下时域冲击信号变得非常密集。为了便于观察重构信号，取800样本点进行分析。将原信号分割成长度为20的块，并进行K-SVD字典学习得到冗余字典。通过本文方法得到的各迭代次数下重构信号包络谱峭度曲线见图7。迭代次数为23时可得到峭度最大的包络谱，故最终迭代次数为23。重构的外圈故障信号见图8，包络谱见图9。图中表明特征频率为f0=140 Hz，与理论值相符。对同一段原信号采用Gabor-OMP方法进行稀疏表示，通过使该方法得到的残差与本文方法得到的残差相一致作为终止条件，重构信号见图10。该方法运行时间7.102 s，选取原子数89。

图6 外圈故障时域信号Fig.6 Vibration signal of outer race fault

图7 不同迭代次数下重构信号的包络谱峭度(外圈故障)Fig.7 The envelop spectrum kurtosis in each step (outer race fault)

图8 外圈故障重构信号Fig.8 The reconstructed signal of outer race fault

图9 外圈故障重构信号包络谱Fig.9 The envelop spectrum of the reconstructed signal of outer race fault

图10 Gabor-OMP方法重构信号(外圈故障)Fig.10 The reconstructed signal by Gabor-OMP(outer race fault)

3.2.2轴承内圈故障特征提取

轴承内圈故障时的时域信号见图11。将原信号分割成长度为32的块，并进行K-SVD字典学习得到冗余字典。通过本文方法得到的各迭代次数下重构信号包络谱峭度曲线见图12。迭代次数为107时可得到峭度最大的包络谱，故最终迭代次数为107。重构的内圈故障信号见图13，包络谱见图14。图中表明特征频率为f1=191.3 Hz，与理论值一致。对同一段原信号采用Gabor-OMP方法进行稀疏表示，重构信号见图15。该方法运行时间为462.187 s，选取原子数397。

图11 内圈故障时域信号Fig.11 Vibration signal of inner race fault

图12 不同迭代次数下重构信号的包络谱峭度(内圈故障)Fig.12 The envelop spectrum kurtosis in each step (inner race fault)

图13 内圈故障重构信号Fig.13 The reconstructed signal

图14 内圈故障重构信号包络谱Fig.14 The envelop spectrum of the reconstructed signal of inner race fault

图15 Gabor-OMP方法重构信号(内圈故障)Fig.15 The reconstructed signal by Gabor-OMP (inner race fault)

3.2.3轴承滚动体故障特征提取

轴承滚动体故障时的时域信号见图16。将原信号分割成长度为32的块，并进行K-SVD字典学习得到冗余字典。通过本文方法得到的各迭代次数下重构信号包络谱峭度曲线见图17。迭代次数为92时可得到峭度最大的包络谱，故最终迭代次数为92。重构的内圈故障信号见图18，包络谱见图19。图中表明特征频率为f2=133 Hz，与理论值一致。对同一段原信号采用Gabor-OMP方法进行稀疏表示，重构信号见图20。该方法运行时间为332.984s，选取原子数293。

图16 滚动体故障时域图Fig.16 Vibration signal of rolling element fault

图17 不同迭代次数下重构信号的包络谱峭度(滚动体故障)Fig.17 The envelop spectrum kurtosis in each step (rolling element fault)

图18 滚动体故障重构信号Fig.18 The reconstructed signal of rolling element fault

图19 滚动体故障重构信号包络谱Fig.19 The envelop spectrum of the reconstructed signal of rolling element fault

图20 Gabor-OMP方法重构信号(滚动体故障)Fig.20 The reconstructed signal by Gabor-OMP (rolling element fault)

不同故障状态下两种方法运行结果对比见表2。从表中可以看出，传统的Gabor-OMP方法需要从冗余度很高的字典中通过计算内积的方式选取原子，当原信号长度增加时，构造字典所需时间以及信号与原子之间内积运算的时间都会显著增加。相比Gabor-OMP方法，本文所提出的基于K-SVD和HBW-OOMP的方法无论是运行时间还是原子数目都有明显缩减，其原因在于K-SVD字典是一种依靠原信号的自学习训练方法，其原子结构与原信号相似度更高；另外，HBW-OOMP方法的分块思想也极大降低了进行内积运算的数据长度，从而减少了计算时间。从表示结果上看，本文方法由于选用的原子数量小，所以重构信号具有更高的稀疏度，并且对冲击信号特征的表示能力更强，更接近于原信号。

表2 两种方法运行结果对比Tab.2 Comparison of the two methods

鲁棒性是评价算法优劣的重要标准，本文通过将不同强度高斯白噪声添加到滚动体故障信号来验证算法的鲁棒性。采用峭度作为衡量信号冲击特征强弱的指标，原信号和重构信号的峭度值见图21。可以看出，在不同的信噪比下，重构信号的峭度值始终高于原信号，说明本文方法对不同强度的含噪信号均能准确地提取出故障冲击特征，具有较强的鲁棒性和适应性。

图21 不同信噪比下的信号峭度Fig.21 The kurtosis in different SNR

4 结语

本文提出了一种基于K-SVD和HBW-OOMP的轴箱轴承故障诊断方法。仿真研究表明本文方法对冲击响应信号的发生时刻以及幅值都能够准确表达，对旋转机械中的干扰信号有较强的抑制能力。将本文方法应用于轴承试验台，所得实际信号的故障特征提取结果表明，该方法能有效提取出外圈、内圈、滚动体信号的故障特征，具有较强的适应性和鲁棒性。相比Gabor-OMP方法在运行时间和表示结果上都有很大优势，具有一定的应用前景。本文方法目前以台架试验数据和仿真数据进行验证，针对实际工程中复杂工况下的故障特征提取还需做更深入的研究和优化。