马云飞，贾希胜，胡起伟，郭驰名，邢 鹏

(1.武警士官学校 军械系，杭州 310023；2.陆军工程大学石家庄校区 装备指挥与管理系，石家庄 050003)

作为装备传动系统的常见组件，轴承在各类型装备中发挥重要作用，及早发现轴承故障并进行相应维修能够有效提高装备的可用性和战备完好性[1]。因此，利用故障预测与健康管理[2]技术对轴承进行退化状态识别具有十分重要的意义。

目前，传统的退化状态识别模型建立在完备先验知识的基础上，即需要知道各退化状态的分界点，得到各个退化状态数据。但在实际中，很难准确地对轴承全寿命数据进行状态阶段划分，导致完备先验知识下的状态识别模型适用性较差[3]。而新轴承在经历短暂磨合期后，处于正常润滑阶段，以及轴承发生了较为严重的故障，处于严重摩擦阶段的情况下，监测者很容易区分出上述阶段，从而获得大量这两个阶段的状态信息。基于此，本文提出了基于正常润滑和严重摩擦两种状态信息的非完备先验知识状态评估模型。

近年来，稀疏编码理论广泛应用于降噪[4]，图像超分辨率重建[5]，图像融合[6]。Wright等[7]根据稀疏字典的判别性，提出了稀疏表示分类方法(sparse representation-based classification,SRC)。此后，学者将SRC理论应用于状态监测领域，开展了大量研究。Zhang等[8]提出基于压缩感知的轴承故障诊断模型，只需要采样和存储少量压缩观测数据，直接使用这些数据即可进行故障诊断。Tang等[9]基于稀疏表示方法和随机降维投影，对原始信号进行降维采样和SRC故障识别。Yu等[10]将时域振动信号转换到频域以减少噪音的破坏，然后利用L1范数和L2范数的组合在字典上对样本进行稀疏编码。苗中华等[11]对人工轴承数据集的振动信号进行稀疏编码特征提取，试验结果表明稀疏特征具有良好的可分性。余建波等[12]等提出基于固有时间尺度分解和稀疏编码收缩相结合的特征提取方法。王维刚等[13]将分类器训练反馈结果融入字典学习过程，将重构误差项、稀疏编码判别项和分类误差项进行整合优化。

尽管已有稀疏编码方法在故障分类领域取得了很大进展，但是：① 直接对时频域信号进行稀疏编码不够精细，不能完全反映信号本质特征；② 已有研究需要所有故障状态的完备信息进行训练，但在实际中完备信息较难获得；③ 少有文献将稀疏编码应用于装备全寿命数据评估。

针对上述问题，本文提出一种基于小波包频带的稀疏编码方法，在此基础上建立非完备信息条件下轴承退化状态评估模型。选取美国辛辛那提大学智能维修中心(intelligent maintenance systems,IMS)的轴承数据对新的稀疏编码效果进行验证，仿真结果表明新稀疏编码比传统稀疏编码能够更早地发现轴承故障。

1 稀疏编码原理

稀疏编码(sparse coding)的基本原理是将原始信号表示为一系列预先选定的字典原子线性组合。假设原始信号y∈Rm，预先训练得到冗余字典H={h1,h2,…,hn}，包含n个字典原子(冗余字典每一列称为原子)。若将y看作字典原子的线性组合，则可将其分解为

(1)

式中，x=[x1,x2,…,xn]T为稀疏系数。当原始信号y和字典H已知，求解稀疏系数x的过程称为稀疏编码。为求解原始信号的最优稀疏表示，稀疏编码通常需要解决以下两方面问题：① 字典学习，稀疏编码的前提是确定能表征任意信号的原子集合，即完备字典；② 系数求解，在给定完备字典H基础上，求解稀疏系数，常用算法有基追踪[14],正交匹配追踪[15]，基于迭代回归的凸松弛算法[16]等。

常见的字典可分为固定字典和自适应字典。其中自适应字典通过对获得的原始信号样本进行学习获得，因此稀疏表示能力更强。本文采用K奇异值分解(K-singular value decomposition，K-SVD)算法进行自适应字典学习。关于K-SVD算法的实现细节可以参考文献[17]。

2 基于小波包频带的稀疏编码方法

考虑到原始时频域信号并不能完全反映出信号本质特征，而通过小波包变换处理轴承数据，能够尽可能多的获取高低频信息，更有利于在轴承退化过程中捕捉到细微的变化。如图1所示，本文将小波包变换与稀疏编码方法相结合，改进方法主要有3项内容。

图1 基于小波包频带稀疏表示与灰色B型绝对关联度的稀疏编码方法Fig.1 The sparse coding approach based on wavelet packet sparse representation and absolute grey relational degree of B-mode

(1)小波包变换和最优频带选取

首先，对原始信号进行4层小波包变换，并根据每层结点的能量选取最优频带，最优频带能够反映最主要的信号特征。假设得到的16个小波包频带为{w1,w2,…,w16}，则可以计算每个频带的能量比重为

(2)

式中，lwp为每层频带的长度，能量比越大表明该层频带包含的信号信息越多，因此根据τ(i)选出最优频带。

(2)小波包频带稀疏分解与压缩重构

选择轴承正常状态数据作为训练样本，对每一个选定的最优频带，利用K-SVD字典学习算法分别训练一个过完备字典，得到字典集合{H1,H2,…,H}(假设ε为选择的最优频带数量)。计算新样本特征时，选用一种基于Laplace先验的多任务贝叶斯压缩感知(Laplace prior-based correlated-sparse-block Byesian compressive sensing,Lap-CBCS)算法[18]，可根据多组数据之间的相关性联合重构，重构精度较高。

具体的，对每个测试样本进行小波包变换得到频带wpi，根据设定好的观测矩阵Φ，得到观测向量

μ=Φ·wi

(3)

之后对wpi利用训练好的稀疏字典Hi进行稀疏分解。如果训练样本与测试样本结构相似，那么稀疏分解会十分有效，重构误差较小。稀疏分解过程如下

μi=Φi·wi=Φi·Hi·θi,i=1,…,

(4)

式中，θi为稀疏系数。由于μi和Φi·Hi已知，本文利用Lap-CBCS算法对多个最优频带进行联合重构。对每个测试样本，最优频带的重构误差δi可以表示如下

(5)

(3)灰色B型绝对关联度计算重构误差向量距离

由于重构误差向量是多维数据，为了对其进行降维，利用灰色B型绝对关联度计算样本与正常轴承重构误差平均值之间的距离作为稀疏编码值。灰色B型绝对关联度(absolute grey relational degree of B-mode,AGRDB)[19]方法基于灰色关联理论计算两个客观事物之间距离，目前广泛应用于工程领域。假设计算两个向量S1和S2之间的RAB，计算方法如下

(6)

其中

(7)

(8)

(9)

3 非完备信息条件下轴承退化状态评估模型

如图2，本文建立一种基于正常状态和严重摩擦状态的轴承退化评估模型，具体如下。

图2 非完备信息条件下轴承退化状态评估模型Fig.2 Flow chart of bearing degradation assessment under incomplete information conditions

(1)退化特征提取。本文采用小波包频带稀疏编码方法进行退化特征提取，提取方法见第2章。在特征提取前，应先对正常轴承信号进行小波包分解，最优频带选择，以及K-SVD字典学习。

(2)模型训练。对正常信号和严重摩擦信号样本进行小波包频带稀疏编码，并计算RAB距离，得到轴承状态边界分别为RAB，normal,RAB,failure。

(3)状态监测。对待测试信号采用相同的特征提取方法计算出RAB,sample，则当前退化状态评价指标Eunknown可以表示为

(10)

由式(10)可知，Eunknown可以用来描述测试信号与正常轴承信号的差异程度：当Eunknown→0时，表示测试信号的状态接近于正常；当Eunknown→1时，表示测试信号的状态接近于严重摩擦。此外，可以预先设定报警的门限值Q：当Eunknown>Q时，则表明轴承退化严重，发布预警信息并制定维修计划；若Eunknown≤Q，则继续监测下一个样本。

4 公开数据验证

采用辛辛那提大学智能维修中心全寿命轴承数据进行方法验证。该试验装置包含4个(Rexnord)ZA-2115双排滚子轴承，在每个轴承的外壳上安装了PCB 353 V33型号的高灵敏度石英ICP加速度计。总共进行3次试验，数据描述如表1所示。轴的转速为2 000 r/min，采样频率为20 kHz，数据采集设备采用NI公司DAQ 6062E采集卡。

表1 IMS轴承数据集参数Tab.1 The parameters of IMS bearing data sets

(1)小波包分解与最优频带选择

选取3组试验中的故障轴承数据进行分析，首先对其进行4层小波包分解，得到对应的频谱和能量分布如图3所示。按照能量从大到小排序，选择其中最靠前的8个小波包频带作为最优频带得到表2所示结果。

表2 IMS轴承数据集小波包分解最优频带Tab.2 The best frequency bands of IMS bearing data sets after wavelet packet decomposition

(2)退化特征提取

对选取的最优频带分别训练一个稀疏字典，由于全寿命数据前一部分数据认为是正常状态，选取每组数据前200～前500个样本作为训练样本。为获取过完备字典，要求字典的规格为m×n(m>n)，本文选取的稀疏字典规格如表3所示，其中线性组合数量为稀疏表示选取的原子数量。为保证收敛，K-SVD学习的迭代次数设置为50。

表3 稀疏字典规格Tab.3 Sparse dictionary specification

训练完成后，对全寿命数据每个数据样本分别按照第2章所示方法进行特征提取，得到小波包分解8个频带的重构误差。如图4所示，以第二组数据集全寿命数据为例，采用模糊C算法对获得的重构误差进行聚类。

图4 小波包频带重构误差特征聚类效果Fig.4 The clustering effect of reconstruction errors for wavelet packet frequency bands

(a)第一组数据

首先将特征归一化，并设置类别数为3，将所有特征值聚成3堆，其中类别1为正常状态，类别2为中间状态，而类别3为严重摩擦状态。为保证聚类既考虑特征值又考虑时间顺序因素，除了获得的8通道小波包特征，还增加了时间因子序列{0/984,1/984,…,983/984}作为第9个特征。由聚类结果可知，重构误差特征在没有先验知识的情况下，可以将全部样本按照退化时间顺序聚为3类。除了中间状态和严重摩擦状态临界处有部分样本混淆外，大致反映了轴承退化规律，由此可以验证小波包频带重构误差具有判别性。

(3)状态监测与预警

得到稀疏编码特征后，按照式(6)计算出轴承状态边界分别为RAB，normal,RAB，failure。其中，计算RAB，normal时选用正常状态数据样本，计算RAB，failure时选用严重故障数据样本。样本区间可根据全寿命退化数据的变化趋势确定(例如轴承退化最后阶段发生突变)。如图5～图7，数据集1和数据集2很容易根据曲线变化确定出正常运行和严重摩擦状态，从而得到训练样本，样本区间如表4所示。但观察数据集3发现其退化趋势不明显，无法区分出正常运行和严重摩擦状态边界，该数据集不适合用本文方法进行状态识别。

表4 训练样本区间Tab.4 The training sample interval

图5 数据集1训练样本选取Fig.5 The training samples selection of data set No.1

图6 数据集2训练样本选取Fig.6 The training samples selection of data set No.2

图7 数据集3退化特征序列Fig.7 The degradation feature sequence of data set No.3

获得待测试样本数据后，可计算出该样本稀疏编码特征，并得到RAB,sample。根据式(10)可求出最终的退化状态指标，并与阈值Q进行比较。本文设置两级预警，其中黄色预警线对应阈值Q=0.07，表明轴承只出现了轻微异常，处于早期故障状态；红色预警线对应阈值Q=0.7，表明轴承已经发生了严重故障。如图8和图9所示，采用两种方法计算数据集2的稀疏编码特征：① 直接进行K-SVD字典学习与压缩重构，即传统稀疏编码;② 采用小波包频带稀疏编码，获得对应重构误差与AGRDB稀疏编码。对两种方法得到的退化特征，分别画出对应的黄色预警线和红色预警线。根据局部放大图可以发现，利用传统稀疏编码方法得到的两级预警点分别为580，959，而利用本文小波包频带稀疏编码得到的两级预警点分别为549，935，本文方法两级预警临界点均早于传统方法，有利于轴承严重摩擦故障与早期故障的预警。此外，在图8(a)和图9(a)中700～800出现了一个较大波动，与其他特征对比，发现是数据集本身工况变化等原因造成的。。

(a)退化状态监测

如果根据实际需要，两级预警线阈值Q发生改变，根据两种方法以及小波包能量退化曲线分别得到的3种预警点如表5和表6所示(当Q≥0.5时，设置红色预警线；当Q<0.5时，设置黄色预警线)。发现小波包频带稀疏编码的绝大部分预警临界点小于传统稀疏编码特征和小波包能量特征。同理，对数据集1计算不同阈值Q下的预警临界点得到表7。观察两组数据集结果，发现当阈值Q较小时，小波包频带稀疏编码的预警临界点可能会大于传统稀疏编码，这主要是因为此时轴承处于正常状态，退化趋势不明显，小波包稀疏编码无法反应出退化特征微小变化。此外，本文在原始数据基础上增加高斯随机白噪声，噪声等级按照信噪比计算

表5 数据集2红色预警线临界点对比Tab.5 The comparison of red warning curve critical points for data set No.2

表6 数据集2黄色预警线临界点对比Tab.6 The comparison of yellow warning curve critical points for data set No.2

表7 数据集1预警线临界点对比Tab.7 The comparison of warning curve critical points for data set No.1

(13)

式中，u和u′分别为原始信号和带噪音信号。如图10，得到6种噪音等级下的稀疏编码特征趋势，重点观察两级预警线区域的特征变化，发现噪音等级对稀疏编码特征的影响不大，且不同特征曲线的整体趋势基本一致，表明本文所提特征具有一定的抗噪性。

(a)退化状态监测

图10 不同噪音下小波包频带稀疏编码特征趋势Fig.10 The trend of sparse coding on wavelet packet subbands under different noises

此外，对比传统稀疏编码和小波包频带稀疏编码特征计算需要的时间复杂度。编码特征计算包括两部分：稀疏字典训练(K-SVD算法)和重构误差求解(Lap-CBCS算法)。由表8可知，尽管小波包频带稀疏编码进行字典学习的时间较短，但其重构时间较长。这主要是因为本文所提特征需要对小波包分解后的每个频带分别进行重构，尽管每个频带长度小于原始信号长度，但把若干个频带重构时间加和后大于对原始信号进行重构的时间。

表8 传统稀疏编码与小波包频带稀疏编码时间复杂度对比Tab.8 The time complexity comparison of traditional sparse coding and sparse coding on wavelet packet frequency bands

5 结 论

本文提出一种基于小波包频带稀疏分解和灰色B型绝对关联度的编码方法，并构建了非完备信息条件下轴承退化状态评估模型，根据稀疏编码计算退化特征评价指标，通过与故障阈值比较进行预警发布。利用IMS中心的轴承公开数据进行仿真实验，结果表明新稀疏编码特征相较于传统稀疏编码特征能够更早的发现故障，有利于及时对轴承进行预防性维修，提高轴承寿命。尽管新稀疏编码特征的计算复杂度要高于传统方法，但考虑到数据传输到状态监测后台可以利用的计算资源较多，且在某些特定场景需要尽早发现故障征兆，因此新方法在实际中存在一定的应用价值。