基于高维随机矩阵特征值之差的滚动轴承状态异常检测算法

朱文昌， 何雅娟， 王建波， 王 恒

(1. 南通大学 机械工程学院，江苏 南通 226019；2. 南通大学 工程训练中心，江苏 南通 226019)

滚动轴承作为旋转机械的关键部件，对其进行异常状态检测为设备预知维护提供了理论指导，避免严重故障的发生而造成的财产损失甚至人员伤亡。如何提高轴承异常状态检测结果的准确性成为研究的热点和难点[1-2]。

Leite等[3]研究了基于熵的12个特征对轴承故障检测性能的影响，使用不同的熵测度作为检测熵特征变化的指标，提出一种非参数方法检测轴承异常状态，但在轴承早期异常状态时，采集的信号成分复杂、故障信息微弱，易造成误判。刘志亮等[4]从轴承的运行状态角度考虑引入了安全域的概念，提出了一种基于核空间距离熵的安全域惩罚参数选择算法，利用支持矢量数据描述(support vector data description, SVDD)对正常样本构建安全域，结合核空间距离熵的安全域惩罚参数选择算法寻找到最优惩罚参数提高了异常检测准确率，但对高维数据处理能力不足。张西宁等[5]在格雷厄姆扫描法的基础上，引入边界软化率增大了数据点外边界的柔性，并结合射线法生成与输入样本反分布的数据集，使传统模型成为拥有精细决策边界的单类随机森林，通过输出待测数据的异常概率进行异常检测，但选用多维缩放的降维方法进行特征约简时可能会造成有用信息损失。Mao等[6]提出了一种使用半监督架构和深度特征表示的轴承早期故障在线检测的新方法, 利用自动编码器从目标轴承的正常状态数据和辅助轴承的故障状态数据中提取深度特征，引入一种安全的半监督支持向量机(safe semi-supervised support vector machine, S4VM)对轴承初期故障可有效识别，当数据样本较大时仍可实现精确检测，但对训练集数据样本要求较高，虚假特征易对S4VM造成干扰，影响最终检测精度。

随机矩阵理论主要研究随机矩阵的样本特征值、样本特征向量及它们的统计函数在一般渐近体系下的渐近行为，目前已经在频谱感知、电网配电、生物医学等领域得到了成功应用[7-9]。Zeng等[10]基于随机矩阵理论提出最大、最小特征值之比算法用于频谱感知，克服了噪声不确定性问题。刘威等[11]基于随机矩阵单环定理对电网海量数据信息进行建模，实现扰动对电网影响程度分析和影响范围评估。Korošak等[12]通过随机矩阵理论的统计特性研究了胰岛细胞的钙离子分布，进而对疾病进行检测。

为了适应设备健康监测在海量数据下得到较为精确的检测结果，本文基于随机矩阵理论，提出一种基于最大、最小特征值之差的异常检测算法对轴承状态进行检测，应用结果表明该了方法的有效性和可行性。

1 基于高维随机矩阵特征值之差的滚动轴承异常状态检测算法

1.1 高维随机矩阵构造

根据滚动轴承健康监测的特点，采用大数据分析架构对其进行处理。设轴承健康监测数据维数为P、监测时间为T、监测节点数为Q，在采样时刻ti，设备第j个节点(如测点)所监测的第n个运行状态特征量定义为监测数据子空间

ljn(ti)(i=1,2,…,T;j=1,2,…,Q;
n=1,2,…,P)

(1)

对节点j而言，监测的所有Q个特征量可以构成一个列向量，即

lj(ti)=[lj1,lj2,lj3…,ljP]T

(2)

将不同采样时刻的监测数据按照时间顺序排列，构成一个时间序列矩阵，即

L=[lj(t1),lj(t2),…,lj(ti),…],
(j=1,2,3,…,Q)

(3)

该矩阵即对滚动轴承监测所得的数据源矩阵L∈cP×Q。

为了对滚动轴承特定时间节点数据进行实时分析，采用平移时间窗对轴承不同时刻的数据信息进行锁定。设时间窗口H的规模为V×W，则在采样时刻ti，对节点j构成矩阵为

Lj(ti)=[lj(ti-W+1),lj(ti-W+2),…,lj(ti)],
(i=1,2,…,P,j=1,2,…,Q)

(4)

在采样时刻ti，窗口矩阵数据构造流程如图1所示。

图1 窗口矩阵数据构造流程图Fig.1 Flow chart of window matrix data construction

(5)

式中，θj为随机生成数，用于矩阵随机化构造。

(6)

(7)

图2 滚动轴承高维随机矩阵构建流程图Fig.2 Flowchart of constructing high-dimensional random matrix of rolling bearing

1.2 异常检测特征值指标构造

由图1、图2所示方法构建轴承ti时刻的高维监测矩阵X，由式(8)计算其协方差矩阵

(8)

式中，σ2为矩阵X的方差。

对Bn进行特征值分解，如式(9)所示

Bn=EΛET

(9)

式中：Λ为矩阵Bn特征值对应的对角矩阵；E为正交矩阵。

滚动轴承数据信号采集中，存在未知的噪声信号，提取出Bn的特征值极值如式(10)所示

(10)

式中：λv-min与λv-max分别为振动信号协方差矩阵最小、最大值；λn为未知噪声协方差矩阵特征值。

为消除噪声干扰，提高检测指标的鲁棒性，采用最大、最小特征值之差构造异常检测指标D

D=λmax-λmin=λv-max-λv-min

(11)

结合阈值对滚动轴承异常状态进行检测，判决方法如式(12)所示

(12)

式中：Y1为检测出轴承为正常状态；Y2为检测出轴承为异常状态；γ为异常检测阈值。

1.3 异常检测阈值确定

(13)

(14)

式中，Bn为X的协方差矩阵。

且满足

(15)

(16)

将滚动轴承的实际工作状态划分为正常状态S1，异常状态S2。考虑到滚动轴承早期故障特征较为微弱，异常信号不明显，在进行异常检测时可能出现误判情况，即轴承实际已进入异常状态S2，但判断为正常状态Y1或轴承实际为正常状态S1，判断为异常状态Y2，造成误判，误警率ηw定义为

ηw=P(Y1|S2)+P(Y2|S1)

(17)

为提高检测的精度，利用误警率ηw研究检测阈值γ。由于轴承异常检测时难以给出S2状态下检测指标D的概率密度分布，无法通过P(Y1|S2)得到一个较为准确的判定阈值。因此，本文基于P(Y2|S1)情况下推导出阈值γ，如式(18)所示。

(18)

将式(13)、式(15)和式(16)代入式(18)中，得到检测阈值γ数学表达式如式(19)所示。

(19)

式中：FT-W(·)为Tracy-Widom第一分布累计函数；σ2为矩阵X的方差；N为矩阵X的列数。综上所述，基于最大、最小特征值之差的滚动轴承异常检测算法流程，如图3所示。其具体步骤如下：

图3 基于最大、最小特征值之差的滚动轴承异常检测算法流程图Fig.3 Flowchart of an algorithm for detecting abnormal state of rolling bearings based on the difference between the maximum and minimum eigenvalues

步骤1将采集到轴承各个时刻的信号按监测时间序列构造出滚动轴承数据源矩阵L；

步骤2利用平移时间窗锁定数据源矩阵各个时间段轴承数据，将锁定的数据矩阵通过分段、随机化、矩阵扩增及维度重构等方法构造出行列之比c(0

步骤3提取高维随机矩阵最大、最小特征值构造出特征值之差指标D，利用随机矩阵理论中的M-P及Tracy-Widom第一分布计算检测阈值γ，结合检测指标对轴承运行状态进行判断。

2 应用研究

2.1 数据来源

试验数据来自辛辛那提大学智能维护系统(intelligent maintenance system，IMS)中心的滚动轴承全寿命试验。试验采用4个Rexnord ZA-2115双列滚动轴承安装在同一主轴上，轴承径向施加2 722 kg恒定载荷，电机转速2 000 r/min。试验结束时轴承1外圈发生故障并失效，此过程共经历164 h，每10 min对轴承进行一次数据采集，采样频率20 kHz，得到采集数据文件984组。本文对轴承1数据进行分析研究。

对采集到的轴承1数据通过式(1)～式(3)可得到监测数据源矩阵L∈c20 480×984，设平移时间窗口H的规模为20 000×2，对轴承不同时刻的数据信息进行锁定；经过4次模拟矩阵的构建，并通过矩阵分段、随机化、扩增及维度重构，最终得到ti时刻高维随机特征矩阵X∈c400×500(行列之比c=0.8)，按照2.2节、2.3节所述方法求得轴承1不同时刻的检测指标D与相应的检测阈值γ，利用式(12)对轴承状态异常检测及判决。

2.2 误警率对异常状态检测的影响

由式(19)可知，检测阈值γ与误警率ηw有关，不同误警率会影响轴承异常检测结果。取不同误警率ηw对轴承1异常点进行检测，结果如图4～图6所示。

图4 当ηw=0.10时轴承1早期异常检测结果Fig.4 Early abnormal detection results of bearing 1 when ηw=0.10

图5 当ηw=0.05时轴承1早期异常检测结果Fig.5 Early abnormal detection results of bearing 1 when ηw=0.05

图6 当ηw=0.01时轴承1早期异常检测结果Fig.6 Early abnormal detection results of bearing 1 when ηw=0.01

由图4～图6可知，当ηw=0.10时特征值指标与阈值在文件号532出现交点，可判断早期异常的发生。文件号532前特征值指标位于阈值下方，进入异常状态后特征值指标整体位于阈值上方。但在交点附近，出现误判现象，即通过阈值判定轴承已进入异常状态Y2，但仍存在少数特征值指标D小于阈值γ。当ηw=0.01时，此时阈值曲线可较好地将轴承1正常与异常状态进行划分，阈值与特征值指标的“混叠”现象基本消除，提高了滚动轴承早期异常检测的准确率，但检测的早期异常点时间滞后。

为进一步研究不同误警率ηw取值时阈值γ与特征值指标D的关系，分别取ηw为0.01，0.05，0.10，对动态阈值曲线与特征值指标曲线之间的绝对误差值之和进行分析，如图7所示。随着误警率ηw减小，检测阈值与特征值指标之间的绝对误差值在逐渐减小，阈值曲线从整体上更逼近于检测指标，提高了检测的精度。因此，可根据采集到的轴承数据复杂程度及实际生产对异常检测准确率进行综合考虑，选择合适的误警率对检测阈值进行调节，以达到最佳的检测效果。

图7 不同ηw下检测阈值与检测指标绝对误差值之和Fig.7 The sum of the detection threshold and the absolute error value of the detection index under different ηw

2.3 滚动轴承异常状态检测结果分析

由2.2节可知，当ηw=0.01时，轴承1的检测效果最佳，检测结果见图6。轴承1检测指标与检测阈值在文件号543时相交，在交点之后轴承1特征值指标上升趋势明显，阈值曲线位于特征值指标曲线下方，判定文件号543为轴承1早期异常点，与文献[16]检测出的早期异常点一致，与峭度指标(如图8所示)检测结果(文件号649)相比，本文所提算法可提前17.67 h检测出轴承异常状态。

图8 基于峭度指标法检测轴承异常状态检测结果Fig.8 Detection results of bearing abnormal state based on kurtosis index method

文献[17]提出了采用随机矩阵最大、最小特征值之比的检测指标，并给出了检测阈值t′，如式(20)所示

(20)

式中：M和N分别为矩阵的行和列；c为矩阵行、列之比；FT-W(·)为Tracy-Widom第一分布累计分布函数。

取误警率ηw=0.01，c=0.8，基于最大、最小特征值之比的特征指标与阈值对轴承1进行异常状态检测，并与基于特征值之差的检测结果进行比较，如图10所示。采用特征值之比在文件号598时检测出轴承异常状态，与本文所提出算法相比滞后约9.2 h，说明最大、最小特征值之差算法对轴承早期异常更为敏感。

图9 基于不同特征值指标的轴承异常状态检测结果比较Fig.9 Comparison of bearing abnormal state detection results based on different eigenvalue indexes

最大、最小特征值之比指标在轴承早期异常点后波动剧烈，单调性较差，这是由于轴承进入退化阶段，提取的最小与最大特征值之间差异逐步增大，λmax/λmin将此差异放大了，造成相邻时刻指标变化加剧，降低了指标的单调性。而采用特征值之差指标对轴承1运行状态进行表征时，数据曲线波动小、稳定性高，可有效消除噪声带来的干扰，降低了轴承在退化阶段数据差异对指标单调性的影响，进一步提高检测指标鲁棒性。

此外，由式(20)可知，最大、最小特征值算法给出的异常阈值主要由矩阵规模c和误警率ηw确定，在两者确定的条件下，异常状态检测阈值为常数。但是，滚动轴承健康监测是动态过程，不同的时间段采集的数据不同，在保证准确率的情况下，检测阈值应随着轴承监测信号变化而改变。本文利用平移时间窗对不同时刻的数据信息进行锁定，结合当前时刻数据方差σ2与规模c给出了动态阈值的计算公式，当轴承处于正常状态时，采集到的数据分布稳定，方差σ2变化较小，阈值接近于常数。当滚动轴承进入异常状态时，随着退化程度的加剧，不同采样时刻数据间差异性越来越大，方差σ2变化较大，阈值随之增加，其变化趋势为随特征指标变化的动态曲线，可实现对轴承异常状态的有效检测。

3 结 论

(1) 利用平移时间窗对滚动轴承不同时刻数据进行锁定，并通过分段、随机化、扩增和维度重构等方法进行高维随机特征矩阵的构造。

(2) 基于随机矩阵理论提出一种最大、最小特征值之差的滚动轴承异常状态检测算法，利用特征值之差构建出的检测指标可降低未知噪声带来的干扰，提高指标的表征能力，结合误警率所推导出动态检测阈值，可实现对早期异常状态的检测。

(3) 利用最大、最小特征值之差算法对美国IMS轴承1进行异常检测，结果表明，与最大、最小特征值之比算法比较，不仅可及早检测出滚动轴承早期异常的发生，且检测效果更优，更符合设备实际退化过程。