陈子旭， 朱振杰， 卢国梁

(1. 山东大学 机械工程学院，济南 250061； 2. 山东大学 高效洁净机械制造教育部重点实验室，济南 250061；3. 山东大学 机械工程国家级实验教学示范中心，济南 250061)

滚动轴承是旋转机械的关键部件，其损坏会直接影响设备的运行过程，导致设备停机、甚至引发人员伤亡。对滚动轴承的健康状态进行实时、在线诊断具有重要意义[1]。

滚动轴承频发的故障大多为局部缺陷形式的潜在故障，类型较复杂，其机械信号特点有：①强烈的非线性、非平稳特性；②常伴有强噪声干扰；②调制源弱导致早期故障特征通常很微弱。这些都给轴承故障在早期阶段的检测与识别带来很大困难[2]。轴承早期故障检测与识别的关键是从振动信号中提取故障特征信息。传统的时域方法通过提取轴承振动信号的峰度、峭度、均方根值等统计指标，对信号进行动态描述，但其容易受到噪声干扰，鲁棒性不足[3]。以傅里叶变换为基础的传统频域分析方法在处理非线性、非平稳信振动号时，常常无法准确描述信号频率随时间的动态特征，制约了其在早期故障检测中的应用[4]。为了弥补上述方法的不足，时频分析方法被相继提出，代表性的时频分析方法包括短时傅里叶变换、小波变换等。这些方法被认为可提供信号完整的时频分布，因此已被广泛应用于轴承故障诊断领域[5]。随着数字信号处理技术的不断进步，近些年来，在传统时频分析方法基础上，自适应时频分析方法得到了较大发展。Huang等[6]提出的经验模态分解(empirical modal decomposition，EMD)解决了传统时频分析存在的非自适应问题，更加适合分析多分量信号，但其也存在端点效应、模态混淆的缺点。Smith[7]提出的局部均值分解(local mean decomposition, LMD)是在EMD的基础上，对端点效应、模态混淆问题进行了改善，能获得原始信号更加精确的时频分布。尽管LMD方法在轴承故障诊断领域获得了较多应用，但其仍存在两个问题：

(1) 由LMD分解得到多尺度分量信号数据量庞大、复杂，因此仍需进一步计算以提取更低维度的数据特征。Tian等[8]提出基于乘积函数(product functions，PF)分量构建汉克尔矩阵，利用奇异值分解压缩数据尺度以获得更加稳定的特征向量值；王洪明等[9]通过计算全矢能量熵从PF分量中提取有用信息，获取得到了低维度的数据特征。

(2) 由LMD分解可以获取信号的强数据特征，而轴承早期故障特征微弱，因此需要对LMD进行增强处理以获取轴承早期微弱故障特征。最大相关峰度反褶积处理能够有效消除噪声对PF分量的影响，用于增强故障特征信息[10]；基于PF分量构建Teager能量谱能够有效识别滚动轴承故障信号中的冲击特性，从而能够提取微弱故障特征[11]；王志坚等[12]提出了一种将掩膜信号法与LMD相结合的LMD-MS(local mean decomposition-mask signal)方法，成功提取了轴承的微弱故障特征。

随着图论的不断发展，图模型作为一种高效的数学模型，已被应用于机械设备状态监测领域[13]。Lu等[14]采用图模型的方式进行数据建模，通过对模型进行距离度量实现对旋转机械的动态描述。Wang等[15]首先采用无向加权图的建模方式基于信号的周期谱进行建模，通过计算模型中位图和距离指标来描述机械设备的动态状态。Wang等[16]在之前工作的基础上进一步扩展了图建模策略，构建了时空图模型并通过比较试验验证了该方法的有效性。受以上研究的启发，本文提出了一种新的图谱域滚动轴承早期故障检测与识别方法。本方法实现了两方面改进：①建立的异常度动态指标实现了低维度轴承健康状态动态描述；②基于PF分量构建图模型进行图谱分析，实现了对LMD的增强处理。通过上述改进，本文所提方法可有效提取轴承早期故障的微弱特征。

1 方法简介

方法主要包含以下3个步骤，如图1所示。

图1 方法流程图

步骤1滚动轴承振动信号动态建模。采用LMD对滚动轴承振动信号进行分解，将分解得到的多尺度分量信号作为输入构建图模型，计算相邻图模型的相似性指标以得到能描述滚动轴承健康状态的动态指标。

步骤2早期故障检测。基于上述构建得到的动态指标，基于假设检验的方法，采用拉依达准则对轴承早期故障进行决策。

步骤3故障识别。对于检测到的故障异常，进一步采用模式识别技术进行故障类型的识别。

2 滚动轴承健康状态动态建模

滚动轴承由于复杂、恶劣的工况导致其振动信号经常受到噪声等冗余信息的干扰，并且不利于微弱故障特征的提取。LMD可有效将噪声等冗余信息同有效信息隔离，进行多尺度分析。同时，随着图论的发展，图建模方式已经被证实在提取微弱故障信息方面的重要作用。因此，针对以上方法的优势，本文基于LMD和图论方法，对滚动轴承健康状态进行动态建模。

2.1 局部均值分解

作为一种自适应分解方法，LMD能将信号分解为一系列具有物理意义的分量信号，反映信号能量在多尺度空间下的分布情况。针对监测获取的滚动轴承振动信号，记为x(t)，其经过LMD分解的步骤如下。

步骤1根据所有的局部极值点ni，求出所有局部极值平均值mi和包络估计值ai，用滑动平均法处理后得到局部均值函数m11(t)和a11(t)。

步骤3将上述迭代过程的包络函数相乘得到包络信号，即a1(t)=a11(t)·a12(t)…a1n(t)。

步骤4将纯调频信号s1n(t)与包络信号a1(t)相乘，得到第一个PF分量，即FP1(t)=s1n(t)·a1(t)。

步骤5将FP1(t)从x(t)中分离出来，得到u1(t)，将其作为振动信号重复以上计算过程得到u2(t)，直到uq(t)为单调函数时停止迭代。

步骤6滚动轴承振动信号经过上述迭代过程，被分解为一系列PF分量和一个残余分量。

(1)

本文首先对滚动轴承振动信号采用无重叠滑动窗的方式进行LMD分解，原始信号被分解为一系列PF分量，即x(t)→{FP1(t),FP2(t),…,uq(t)}。从LMD定义来看，残余分量uq(t)为一单调函数，其包含的振动信号极其微弱，具有的能量接近于零；并且振动信号分解后的前5个分量几乎包含了滚动轴承振动信号的所有频率成分，其能量和与信号的总能量近似相等[17]。因此，本文选取前5个PF分量用于图建模。

2.2 图建模

图模型通常由一系列的节点V和边L组成，即G={V,L}。通过对每对节点vi和vj之间计算距离，比如欧几里得距离、曼哈顿距离、赫凌格距离等，以表示两节点所构成的边lij的权重。图模型可以有效的表示相邻节点以及互连节点之间的连接信息，并且应用到信号的多尺度分析中。本文将分解得到的前5个PF分量作为图模型的节点序列，并采用最常用的欧氏距离度量方式，赋予图模型每条边权重进行图建模，如图2所示。本文以t1时刻为例，给出了图建模过程示意图，其建模过程分为以下几个步骤。

图2 图建模示意图

步骤1将分量FP1～FP5作为图模型的节点序列，5个分量在t1时刻的取值分别为{FP1,1,FP2,1,FP3,1,FP4,1,FP5,1}，作为图模型的节点{v1,v2,v3,v4,v5}，每对节点vi和vj之间的连接作为图模型的边lij，其中i,j∈[1,5],Z。

步骤2计算每对节点vi和vj之间的欧式距离dij，作为边lij的权重对其进行赋值。

di,j=|vi-vj|

(2)

步骤3将计算得到的欧式距离dij如式(3)所示，以邻接矩阵的形式表示，得到图模型的矩阵表示。

(3)

使用来自每个PF分量的多尺度信息在每个数据点上构建图模型。因此，每个滑动窗内部的振动信号则可以表示为一系列的图模型，即{G1,G2,…,Gl}，l为滑动窗的长度，本文将l设为5 000。构建的图模型除了可以描述时序特性，还包含信号不同尺度之间的相关信息。

在图论中，均值图是用于度量一组图模型平均水平的常见方式。对于在每个滑动窗内构建得到的图模型序列，计算其均值图用以对该段数据状态进行描述，计算公式如下

(4)

2.3 相似性度量

构建得到的图模型具有较高的维度，同时，为了更好地描述滚动轴承动态特性、便于故障检测的决策，引入了相似性度量指标。本文认为LMD分解得到的各分量对于提取故障信息具有相同的贡献度，因此，采用欧式距离度量方式计算相邻图模型之间相似性，得到异常度分数序列，对滚动轴承运行过程中状态变化的可能性进行量化。计算公式如下所示

(5)

3 早期故障检测

当滚动轴承正常运转时，其振动信号分布接近正态分布，此状态下图模型结构稳定，异常度分数服从高斯分布；而当滚动轴承在某一时刻发生故障时，随之而来的是故障特征频率的出现和冲击成分的增多，其振动信号将偏离正态分布。基于此，以高斯分布为假设，采用拉依达准则对异常度分数进行假设检验(本文采用6σ准则)，实现早期故障时刻的定位。原假设和备择假设如下所示

(6)

式中，μh-1和σh-1分别为历史数据{s1,s2,…,sh-1}的平均值和标准差，计算公式如下

(7)

(8)

假设H0设定：当前时刻数据与历史数据服从同一个正态分布，即滚动轴承处于正常运行状态。备择假设H1设定：当前时刻数据出现了异于正常的状态变化，该时刻即为检测到的由正常转向故障的状态变化点。实线为异常度分数曲线，虚线为上下6σ控制限，圈表示检测到的状态变化点，如图3所示。当滚动轴承运行出现了异于正常的状态变化时，图模型结构的变化会导致异常度分数发生较大变化，异常度分数不再服从先前分布，通过判断当前时刻异常度分数值是否处于此前分布的置信区间，便可以确定当前时刻滚动轴承是否发生状态变化。

图3 早期故障检测示意图

4 故障识别

基于第3章早期故障检测，一旦检测到滚动轴承运行状态发生异常，可利用模式识别方法(如支持向量机(support vector machines, SVM)、人工神经网络、随机森林、深度学习等)实现故障类型的识别。在常见的识别方法中，K最近邻(k-nearest neighbors, KNN) 算法被广泛应用。KNN算法通过测量不同特征值之间的距离，依照多数表决的决策规则进行分类，即由待测样本的K个最近邻的训练实例中的多数类别决定待测样本的类别。本文给出了KNN算法用于分类的示意图，如图4所示。首先计算测试数据与各训练数据之间的距离，按照距离的递增关系进行排序，如果一个样本在特征空间中的K个最相似(即特征空间中最近邻的前K个点)的样本中的大多数属于某一个类别，则该样本也属于这个类别。KNN是一种非参的、惰性算法模型，其简单易用，对异常值不敏感，可用于非线性分类，同时具有较低的运算复杂度，无需大量数据进行网络训练[18]。

图4 KNN算法原理图

本文采用KNN开展故障类型识别的试验。需要特别指出的是，本文主要关注点在于提出了一种基于图的振动信号建模与故障表征方法，因此未对其他模式识别方法进行验证。实际上，本文方法也可与其他模式识别方法进行结合，以进一步提高其准确率。一旦检测到滚动轴承运行状态发生异常，则状态变化点之后的数据为待识别数据，将待测数据的均值图结合KNN算法对故障类型进行识别。值得提出的是，均值图表示为邻接矩阵的形式，在输入KNN之前需将其转换为行向量。

5 试验及讨论

5.1 试验数据

分别在XJTU-SY[19]和凯斯西储大学(Case Western Reserve University, CWRU)[20]数据集上进行了验证。其中，XJTU-SY数据集用于验证早期故障检测，该数据集是由西安交通大学设计科学研究所和基础构件研究所以及长兴晟阳科技有限公司提供的滚动轴承加速退化试验数据。该试验数据在3种工况下取得：①转速2 100 r/min，径向力12 kN；②转速2 250 r/min，径向力11 kN；③转速2 400 r/min，径向力10 kN，每个工况下有5个轴承。数据采样频率为25.6 kHz，采样间隔为1 min，每个采样持续时间为1.28 s，失效原因包括外圈磨损、内圈磨损、滚动体损坏、保持架断裂等。每一个测试轴承的详细信息，如表1所示。有关该试验系统以及试验详情见文献[19]。

表1 XJTU-SY轴承数据集信息

CWRU数据集用于验证故障类型识别，CWRU数据集来自凯斯西储大学轴承数据中心。数据集包含正常数据、采样率12 kHz的驱动端故障数据、采样率48 kHz的驱动端故障数据、风扇端故障数据。本文采用48 kHz驱动端故障数据进行试验，试验数据在4种工况下获得：①负载0，转速1 797 r/min；②负载745.69 W，转速1 772 r/min；③负载1 491.39 W，转速1 750 r/min；④负载2 237.07 W，转速1 730 r/min。每种工况下包含滚动轴承4种工作状态：①正常(Normal)；②外圈磨损(OF)；③内圈磨损(IF)；④滚动体磨损(RF)，其中，3种故障类型均是通过放电加工产生，每种故障类型包括3种不同的破坏尺寸，如表2所示。数据集共包括10种不同类型的数据。有关该数据集更为详细的介绍见文献[20]。

表2 CWRU 驱动端故障数据类型

5.2 故障检测试验

XJTU-SY数据集既包含缓慢失效的故障信号，又包含突然失效的故障信号。本文方法在故障检测阶段的两个例子，如图5、图6所示。图5(a)、图6(a)为XJTU-SY数据集中数据Bearing1_3和Bearing3_3，分别对应这两种不同退化形式的信号。图5(b)、图6(b)为本文检测结果，其中圈为检测得到的状态变化点。若故障发生后，方法发出预警，认为方法成功检测到状态变化，并且预警时刻与真正开始发生故障时刻的延迟越小则检测效果越好；若方法发出预警，但轴承此时实际并未发生故障，则认为方法此时为误报。对于突然失效的故障信号，其时域波形发生较大突变，通过肉眼便能确定故障开始时刻[21]；而对于缓慢失效的故障信号，通常对应轴承退化的多个阶段，很难确定故障开始时刻[22]。因此针对缓慢失效故障信号，对检测到的状态变化点所在滑动窗数据时域波形和频谱分析，并与正常数据进行比较，以确定方法检测的准确性。本文以Bearing1_3数据为例进行说明，该轴承在试验时发生了轴承外圈故障，正常数据和检测到的故障数据的时域波形和频谱图，如图7、图8所示。从时域来看，检测到的故障数据与正常数据相比波形上有较大差距，出现了较多故障发生时典型的调幅和冲击成分。该数据集轴承一倍转频理论值为34.38 Hz，从频谱图来看，正常数据中主要包含35.85 Hz，102.40 Hz，291.90 Hz，327.70 Hz频率成分，分别对应轴承转动的基频和倍频；检测到的故障数据频谱中，除了一倍转频成分35.85 Hz外，还出现了107.50 Hz及其二倍频215.10 Hz、三倍频322.60 Hz等倍频成分，该频率成分接近外圈故障特征频率理论值107.91 Hz，由此可以获知本文方法成功做出状态变化预警。

图5 Bearing1_3检测结果

图6 Bearing3_3检测结果

图7 正常数据波形及频谱

图8 故障数据波形及频谱

为了进一步验证方法的优势，本文将所提方法与广泛采用的均方根值(root mean square, RMS)、峭度(Kurtosis)、偏度(Skewness)、均方误差(mean square error, MSE)、标准差(standard deviation, SD)、方差(variation, Var)、特征能量(signature energy, SE)、特征能量比(signature energy ratio, SER)共8种指标进行了比较。将每种方法首次检测到的异常警报作为轴承早期故障，利用查准率、查全率和综合测度(F1_score)3个指标对每个方法进行综合评价，计算公式如下

(9)

(10)

(11)

所有方法的评价结果，如表3所示。其中，本文方法同SE、SER、MSE、Var指标均达到了查全率100%，查准率100%，综合测度1.0，优于RMS、Kurtosis、SD、Skewness指标。进一步地，本文在表4给出了本文方法同SE、SER、MSE、Var 4种方法预警时刻的对比。以表4中Bearing1_3为例，本文方法在第386个滑动窗处发出预警，早于其他方法预警时刻。从表4可知，除Bearing1_2和Bearing3_2，本文方法发出故障预警的时刻均早于其他方法，说明本文方法能够更有效提取滚动轴承振动信号中的早期微弱故障信息。

表3 试验对比结果

表4 不同方法的预警时刻

在实际应用中，计算时间同样对早期故障检测有着重要影响[23]。本文以“Bearing1_3”为测试数据，在MATLAB R2017b，i5-5200U CPU 2.2 GHz，RAM 4.00 GB操作环境下，对不同方法计算时间进行测试。测试数据前386个滑动窗内数据每个滑动窗内的平均计算时间，如表5所示。几种方法前386个滑动窗数据总的计算时间，如图9所示。相比4种对比方法，本文方法计算时间较高，并且随着数据量增加升高明显。但是，值得提出的是，本文方法在第386个滑动窗处已成功发出预警，而其他方法在较晚时刻，甚至是完全失效阶段才发出预警。相比其他方法，本文方法对早期故障阶段的信息更为敏感，能在故障发生后的较早阶段发出故障预警，并且对于缓慢失效过程尤为明显。因此，本文通过LMD和图论方法提取信号多尺度信息，相比其他方法能够更有效提取早期微弱故障特征。在实际应用中，在设备故障的早期阶段发出预警可帮助操作人员及时排查设备故障原因、避免损失进一步扩大，为停机操作和制定故障维修计划争取更多的时间。作为未来工作方向之一，我们将继续改进所提方法，例如采用并行计算等方式提高方法计算效率，使其更加适用于实际应用。

表5 平均计算时间

图9 测试数据计算时间

5.3 故障识别试验

首先利用t-SNE方法将故障数据映射到低维空间进行三维可视化。不同形状的散点表示不同类型的轴承数据，如图10所示。如IF7为轴承内圈故障，故障尺寸为0.177 8 mm，其他具体图例含义同表2。由图10可知，相同类型的数据具有较好的聚集性，并且不同类型数据之间具有较好的分离性。这说明所建图模型可以有效的对不同类数据进行描述，便于分类器进行分类实现故障诊断。

图10 故障数据在特征空间的三维可视化

本文将构建的图模型的50%用于训练，50%用于测试。在数据集给出的4种负载条件下，对这10种类型数据利用KNN算法进行分类识别。本文方法在4种负载条件下，分别能达到90%、100%、100%、100%的准确率，为了进一步说明方法的优越性，本文同现有的方法进行了对比，本文方法和对比方法的试验结果，如表6所示。从表6可知，对于在负载745.69 W、1 491.39 W、2 237.07 W条件下，本文方法的表现为最优。在0条件下，杨云博等[24-26]准确率优于本文方法。杨云博等利用一种改进的人工鱼群算法对SVM进行优化，通过提取到的振动信号时域特征进行故障诊断；王太勇等提出了一种基于DenseNet的卷积核dropout(KD)智能故障诊断模型KD-DenseNet。其中，杨云博等只考虑了正常、外圈磨损、内圈磨损、滚动体磨损4种故障类型，并且并未将故障破坏尺寸考虑在内；王太勇等试验中的故障类型只有7种，没有考虑故障破坏尺寸为0.533 4 mm时的故障类型。本文认为故障类型数目的不同是导致这两种方法准确率均高于所提方法的原因。并且，需要指出，在实际应用中，轴承的故障类型多种多样，不同的工况会导致不同的故障破坏尺寸，本文方法将CWRU数据集中的10种故障类型全部考虑在内，分类更加精确，具有更好的现实意义。Zeng等提出了一种基于柔性凸壳的最大余量分类方法，并将其应用于滚动轴承的故障诊断，该方法将83%的数据用于训练，17%的数据用于测试，在负载0时，其准确率达到99.60%，并且故障种类数为10。相比该方法，本文方法使用了更少的训练数据，具有更高的实际应用价值。在实际应用过程中，轴承故障类型多种多样，轴承故障类数据又通常难以大量获取，因此，使用较少的故障数据完成故障诊断任务更加符合实际应用需求。因此，相比其他方法，本文所提方法具有较大优势。

表6 故障识别结果

6 结 论

针对滚动轴承健康状态监测问题，本文提出了一种新的图谱域滚动轴承早期故障检测与识别方法。通过基于PF分量构建图模型进行图谱分析，能够有效获取滚动轴承早期故障特征，通过计算相邻图模型间的相似性指标对轴承健康状态进行低维度动态描述，实现了滚动轴承早期故障检测。利用LMD对滚动轴承振动信号进行多尺度分析，通过图建模的方式对得到的多尺度数据增强，能够很好地用于故障类型描述，达到滚动轴承故障识别的目标。在后续工作中，将研究提高方法的运算效率，使其更好的应用于实际场合中。