基于MMFE和可拓k-medoids聚类的轴承性能退化评估

赵聪聪，刘玉梅，赵颖慧，白 杨，施继红

(1.吉林农业大学 工程技术学院,长春 130118；2.吉林大学 交通学院,长春 130022；3.一汽研发总院 智能网联开发院,长春 130011；4.一汽大众汽车有限公司 技术开发部,长春 130011)

轴承是机械系统重要的零件之一，其运行状态和使用寿命直接影响机械系统的运行可靠性[1]。轴承性能退化评估与故障诊断也一直是机械故障诊断领域的研究热点与难点。轴承实际运行状态受多种因素的影响，若能建立有效的性能退化评估模型，构建状态参数与轴承性能退化之间的映射关系，则可实现对轴承性能退化程度的实时评估。这对提高轴承的运行安全性和机械系统的工作可靠性均具有重要意义[2]。

轴承性能退化评估涉及两方面内容，即轴承状态参数的提取和性能退化评估模型的构建。目前，在轴承性能退化评估和故障诊断领域，大多以轴承垂向振动信号为研究对象，忽略了其他方向的振动信息以及不同方向振动信息之间的关联性，导致所提取的特征无法全面反应轴承的运行状态。随着非线性科学理论的发展，许多非线性动力学方法如样本熵、近似熵、模糊熵、多尺度熵等被广泛应用于机械故障诊断领域，并取得了较好效果[3]。但上述熵为单变量分析方法，不能描述不同方向振动信息之间的关联性。Ahmed等[4]将多元样本熵的概念拓展到多尺度，提出了多元多尺度样本熵(Multivariate Multi-Scale Entropy，MMSE)，以从复杂性、互预测性和长时相关性的角度来描述多通道时间序列的动态相关关系[5]。MMSE自提出以来，就在机械故障诊断领域得到了广泛应用。

为实现对轴承性能退化程度的分析，需构建定量评估模型[6]。退化评估模型包括概率模型和距离模型，前者如高斯混合模型、隐马尔可夫模型等，后者如支持向量机、支持向量描述及聚类算法等。k-medoids算法作为聚类分析方法的一种，改进了k-means算法初始中心点可能取到非样本点的缺陷，提高了k-means算法的聚类效果。Rai等[7]将EMD分解与k-medoids聚类算法相结合，利用无故障样本和失效样本建立聚类模型，实现了对轴承性能退化的有效评估。利用k-medoids算法对样本进行聚类可以得到聚类中心，但需进一步根据样本与聚类中心的距离来定量分析轴承的性能退化程度。为此，将可拓学关联函数的概念引入到轴承性能退化模型的构建，使用可拓关联度实现对轴承性能退化程度的定量评估。

基于上述分析，本文提出了一种结合多元多尺度熵、k-medoids算法和可拓学理论的轴承性能退化评估方法。以多尺度模糊熵为基础，利用轴承多向振动信号构建多元变量，提取多元多尺度模糊熵(multivariate multi-scale fuzzy entropy，MMFE)信息。对轴承正常状态样本进行k-medoids聚类得到聚类中心，根据样本与聚类中心之间的欧式距离确定可拓集合的经典域和节域，同时分析可拓关联函数的适用范围。最后根据待测样本到聚类中心的距离与轴承正常状态距离区间的关联度实现对其性能退化程度的定性定量评估，并利用轴承全寿命疲劳试验进行验证。

1 MMFE算法

1.1 多元模糊熵

为计算多元模糊熵，首先需根据Takens嵌入理论生成多元嵌入向量。

Xm(i)=[x1,i,x1,i+λ1,…,x1,i+(m1-1)λ1,x2,i,

x2,i+λ2,…,x2,i+(m2-1)λ2,…,xp,i,

xp,i+λp,…,xp,i+(mp-1)λp]=

[zi,zi+1,…,zi+m-1](i=1,2,…,N-n)

(1)

经过多元嵌入重构后，得到N-n个复合延时向量Xm(i)。

max{|zi+l-1-zj+l-1|,l=1,2,…,m}

(2)

定义模糊隶属度函数。

(3)

式中，n1和r分别为模糊隶属度函数的边界梯度和宽度。

(4)

(5)

对所有i求平均值，即得到嵌入维度m时的条件概率Bm(r)。

(6)

保持其他变量不变，将多元复合延迟向量由m维嵌入到m+1维。对于一个嵌入向量M=[m1,m2,…,mp]，系统可以转换成嵌入向量为M=[m1,m2,…,mk+1,…,mp]的任意空间。由于包含p个时间序列，通过分别拓展mk+1(k=1,2,…,p)，可获得p×(N-n)个重构向量Xm+1(i)。

(7)

(8)

多元模糊熵定义如下。

(9)

1.2 多元多尺度模糊熵MMFE

MMFE考虑了多通道时间序列中各序列的自相似性和序列之间的互预测性[8]，计算过程如下。

(10)

式中，τ=1,2,…为尺度因子。

2 可拓k-mediods轴承性能评估模型

2.1 k-mediods算法

为实现上述过程，在使用k-medoids算法进行聚类时，首先在数据集A中随机选择k个点作为初始中心，即初始代表点，将剩余数据点依据就近原则分配到各类中；然后在各类中按顺序遍历选取非代表点代替初始代表点，重新计算聚类效果；通过反复迭代改进聚类质量，直至获得稳定的聚类中心[10]。

2.2 可拓距及关联函数

可拓学是我国学者蔡文教授创立的一门新兴学科，主要研究事物拓展的可能性，为解决矛盾问题和不相容问题提供了新的方法和途径[11]。可拓距作为可拓数学的基本概念，是计算可拓关联函数的基础。对于区间X0=[a,b]内的任一点x，经典数学认为x与X0的距离为0，而可拓学利用距描述x与X0的相对位置关系，将经典数学“类内即同”的思想拓展到“类内尚可分为不同层次”。设实轴上任意点x，其与区间X0的可拓距表示如下。

(11)

由式(11)知，x∉X0时，ρ(x,X0)>0；x=a或b时，ρ(x,X0)=0；x∈X0时，ρ(x,X0)<0。

可拓距的概念说明了点与区间的位置关系，而关联函数k(x)能够描述点与区间的隶属关系，且关联函数以可拓距为基础。设x0表示关联函数取最大值时的点，即最优值点；x为实域上的任意元素，区间X0=〈a,b〉，X=〈c,d〉，X0⊂X，且X0、X无公共端点，则元素x关于区间X0的关联函数与最优值点x0的位置有关。当x0在X0中点处取得时，初等关联函数k(x)表示如下。

(12)

当x0不在X0中点时，关联函数k(x)表示如下。

(13)

式(12)～(13)中，X0=〈a,b〉和X=〈c,d〉分别代表可拓集的经典域和节域；k(x)的正负和大小表明了x属于或不属于X0的程度；ρ(x,x0,X0)称为侧距。若最优值点x0∈[a,(a+b)/2]，则ρ(x,x0,X0)称为左侧距，记为ρl(x,x0,X0)；若x0∈[(a+b)/2,b]，则ρ(x,x0,X0)称为右侧距，记为ρr(x,x0,X0)。

(14)

(15)

关联度k(x)描述了事物具有某种属性的程度，能够对事物属性进行定性定量评价。由式(12)～(15)知，k(x)∈(-∞,+∞)，k(x)>0表示事物具有某种特征或满足某种属性；k(x)<0表示事物不具有该特征或不满足该属性；k(x)=0表示事物即具有该性质，又不具有该性质，属于性质变化的中间阶段。

2.3 轴承性能退化评估模型

由于可拓学关联函数能定量描述点与区间的关联程度，故将其与k-mediods算法相结合构建轴承性能退化评估模型。具体思想如下：利用k-mediods算法对轴承正常状态样本进行聚类，获得聚类中心Cnormal；通过轴承状态样本与Cnormal之间的距离确定经典域X0和节域X；对轴承任一状态样本而言，计算该样本到Cnormal的欧氏距离与经典域X0之间的关联度，根据关联度完成轴承性能退化程度的定性定量评估。具体流程如图1所示。

图1 轴承性能退化评估流程Fig.1 Assessment process of bearing performance degradation

在利用上述方法构建轴承性能退化评估模型时，应明确两点问题：①可拓集合经典域X0和节域X的确定；②关联函数的形式选择。可拓集合的确定过程如图2所示：首先利用k-mediods算法对轴承正常状态样本进行聚类，得到聚类中心Cnormal。设n1和n2是正常样本中与Cnormal之间欧式距离的最小值点和最大值点，将欧式距离记为d(Cnormal,n1)和d(Cnormal,n2)，以此作为经典域边界值的近似参考，即令a=d(Cnormal,n1)，b=d(Cnormal,n2)。为确定节域的边界值，计算全部样本到Cnormal的欧式距离，并以最近点n3和最远点n4与Cnormal之间的距离d(Cnormal,n3)和d(Cnormal,n4)作为节域的近似参考，即令c=d(Cnormal,n3)，d=d(Cnormal,n4)。本文中，样本与Cnormal之间的距离(记为dis)越趋于a，则样本越靠近聚类中心，表明其所代表的轴承运行状态属于正常的程度越大；反之，dis越趋于b，样本越远离Cnormal，其所代表的轴承运行状态属于正常的程度越弱。因此，最优点在区间X0的左端点处取得，故采用式(13)和式(14)进行关联度的计算。

图2 可拓集合参考点的选取Fig.2 Selection of reference points to extension sets

3 MMFE对轴承运行状态的表征能力分析

3.1 试验装置

为说明MMFE在表征轴承运行状态方面的优越性，现利用美国辛辛那提大学智能维护中心提供的滚动轴承全寿命疲劳试验数据进行分析[12]。试验台结构如图3所示，主轴装有四个Rexnord ZA-2115双列滚柱轴承。轴承节圆直径7.15 cm，滚动体直径0.84 cm，接触角15.17°，利用杠杆机构通过轴承座2和3向主轴施加2 721.6 kg的径向载荷。

图3 滚动轴承全寿命疲劳试验台Fig.3 Test bench for rolling bearing run-to-failure

交流电机通过皮带耦合驱动主轴，电机转速2 000 r/min。每个轴承座的水平和垂直方向装有PCB353B33压电式加速度传感器，采集8通道时间序列。轴承疲劳试验全过程共采集2 156个数据文件，采样频率20 kHz。前43个文件的采样时间间隔为5 min，剩余文件采样间隔为10 min。试验结束时，轴承3发生内圈故障，轴承4发生滚动体故障。本文以轴承3为分析对象，利用垂向和横向的振动加速度信号构建二元多尺度模糊熵。

由于所获得的样本文件较多，为便于分析，每半小时提取一次轴承振动数据作为样本，共获得690个样本。轴承在试验初期处于磨合状态，振动幅值较大，会影响分析结果的准确性。因此，本文剔除前40组样本，利用剩余的650组样本进行分析。

3.2 MMFE参数选择

由MMFE定义知，嵌入维数向量M、多维时延向量λ、时间序列长度N、模糊隶属度函数的梯度n1和宽度r都会影响MMFE的计算结果。一般嵌入维数mk=2或3、λk=1(k=1,2,…,p)时，对MMFE的影响较小。样本长度N与嵌入维数m满足N=10m～30m。由于采样频率为20 kHz，为保证频率分辨率和信息的完备性，样本长度N应大于1 000点，故嵌入维数mk=3。依据文献[13]，设定n1=2，r=0.15SD(SD为多通道数据的标准差)。现利用试验分析的方法确定样本长度N。首先利用db3小波分解方法对轴承振动信号进行降噪，而后在轴承全寿命疲劳试验的初期、中期和末期各取10组数据(初期选择第21～30组，中期为331～340组，末期为651～660组)，以代表轴承性能退化的不同阶段。分别计算样本长度N=1 000、2 000、3 000、5 000和7 000时的MMFE均值，其中尺度因子τ=20，结果如图4所示。

由图4可知，随着样本长度的增加，各尺度下MMFE波动减小，趋于平稳。这说明样本长度越大，MMFE对轴承运行状态的表征能力越稳定。但应注意，样本长度增加，计算量随之增加。当样本长度N=1 000时，滚动轴承任一运行状态下振动信号的二元MMFE波动较大，此时MMFE对轴承运行状态的表征能力较差。轴承二元振动信号各尺度下的MMFE在N=5 000和7 000时相差不大，几近相等。当N=3 000时，各尺度下MMFE波动较小，且与N=5 000或7 000时的MMFE计算结果相差较小，特别是在疲劳试验末期。综合上述分析，本文取样本长度N=3 000，得到滚动轴承在不同运行状态下的MMFE变化情况，如图5所示。

(a) 滚动轴承初期振动信号的MMFE

图5 N=3 000时MMFE变化Fig.5 The change curve of MMFE with N=3 000

由图5知，在轴承性能退化的不同阶段，各尺度下的MMFE不同，说明MMFE可作为表征轴承运行状态的特征参数。随着疲劳试验的进行，轴承从正常状态经一系列性能退化过程到最终失效，振动信号的随机性减弱，确定性增强，信号构成由复杂逐渐变得简单，故MMFE逐渐减小，与图5结果相符。

3.3 MMFE与MFE性能对比

为说明轴承垂向振动信号与横向振动信号之间的关联性，表明MMFE比多尺度模糊熵(MFE)能更有效地表征轴承运行状态，现分别对前文滚动轴承三种状态下的30组垂向和横向样本进行分析，计算每种状态下10组样本的MFE均值。参数设置同MMFE算法，数据点数N=3 000，结果如图6所示。

(a) 滚动轴承横向振动信号MFE

由图6可以看出，单向振动信号的MFE对轴承运行状态的区分能力弱于MMFE，尤其在小尺度和高尺度因子的情况下。此外，对比图5和图6可知，MFE难以区分轴承的初期和中期状态，而MMFE在小尺度因子下能有效区分这两种状态，说明MMFE比MFE能更有效地识别轴承早期性能退化。为进一步说明轴承横向与垂向振动信息之间的相关性，现以滚动轴承运行中期的10组样本为例，分析两向振动信号MMFE和单向振动信号MFE的均值，结果如图7所示。

图7 MMFE与MFE特征对比Fig.7 Comparison of MMFE and MFE

根据相关性理论知，两组数据之间的相关性越高，复杂度越低，熵值越小；反之亦然。由于MMFE均值在整个尺度因子上小于MFE，表明MMFE考虑了轴承两向振动信息之间的相关性。相比于轴承单向振动信号的MFE特征，各尺度下MMFE变化更为平稳，且在较高尺度下MMFE趋于平稳。因此，利用轴承两向振动信号所构建的MMFE特征能更准确的描述轴承运行状态。

4 轴承全寿命疲劳试验分析

4.1 最优尺度MMFE选取

在利用20个尺度下的MMFE构建高维特征向量进行聚类时，会导致计算量过大，聚类精度降低。此外，由图5知，不同尺度下的MMFE所包含的信息量不同：1≤τ≤9时，MMFE能较好地区分滚动轴承的不同运行状态；10≤τ≤20时，MMFE的区分能力减弱。因此，需选取最优尺度下的MMFE来构建轴承运行状态特征向量。主成分分析(principal component analysis，PCA)具有去除冗余信息、降低特征向量维数、提高识别率、保留主要特征的特点[14]。因此，本文利用PCA方法对20维MMFE特征向量进行降维。前10个尺度的MMFE累计贡献率见表1。

表1 MMFE特征主成分分析Tab.1 Principal component analysis of MMFE

由表1知，前7个尺度下的MMFE累计贡献率大于90%，涵盖了绝大部分的有用信息。由图4和图5知，轴承不同运行状态的MMFE在第8个尺度取得最大值。故为保证有用信息的完备性，利用前8个尺度下的MMFE构建滚动轴承的运行状态特征向量。

4.2 可拓k-medoids聚类分析

在利用k-medoids算法获得表征轴承正常运行状态的聚类中心Cnormal时，所采用的样本全部来自于轴承正常运行状态，故样本类别k=1。轴承正常状态样本数量Num的多少会影响Cnormal的准确性。样本量过少，无法覆盖轴承正常运行状态的可能范围，Cnormal与实际值存在较大偏差。由于无法确定轴承性能退化的初始时刻，若样本量过大，部分样本可能来源于轴承早期性退化或故障状态，导致样本分布变得分散，同样会使Cnormal偏离实际值。本文采用不同的样本数量进行分析。令样本量Num=20～440，每次分析时使样本量增加20，计算不同样本数量下各样本点与Cnormal之间欧式距离的平均值，记为Dis，结果如图8所示。

由图8知，当样本量Num=20时，平均距离Dis取得最小值。但由于样本量过小，所获得的Cnormal为局部中心，不能代表轴承全部正常状态样本的聚类中心。随着Num的增加，样本分布先分散后集中，Dis先增大后减小。当Num=180时，Dis取得极小值。随着样本量的进一步增加，样本的聚类性减弱，平均距离再次增大，表明样本分布变得分散。故利用k-dedoids聚类方法确定Cnormal时，样本数Num=180，由此得到Cnormal=[0.371 6,0.415 4,0.478 2,0.496 2,0.541 0,0.572 7,0.559 6,0.615 2]。

图8 样本与Cnormal之间的平均距离Fig.8 The average distance between the samples and Cnormal

为表明Cnormal与轴承全部运行状态下聚类中心Call的位置区别，现利用k-dedoids方法对代表轴承全部状态的650组样本进行聚类。此时不需要划分轴承的性能退化阶段，样本代表了轴承运行中的任一状态，故样本类别k=1。聚类结果如图9所示，聚类中心Call=[0.350 0,0.400 5,0.473 8,0.485 4,0.542 7,0.566 3,0.546 8,0.603 5]。

(a) 尺度因子τ=1 ～3

在获得聚类中心Cnormal后，结合图2利用样本点与Cnormal之间的欧式距离来确定经典域X0和节域X。经典域X0代表了轴承正常状态样本与Cnormal之间欧式距离的可能范围，故计算轴承180组正常状态样本与Cnormal之间的欧式距离，并以其最小值和最大值作为X0的下限和上限。节域代表了轴承任一状态样本与Cnormal之间欧氏距离的可能范围。同理，根据轴承全部的650组样本与Cnormal之间的欧式距离来确定节域X。通过计算，180组正常状态样本与Cnormal之间欧式距离的最小值和最大值分别为0.030 5和0.090 6，650组样本与Cnormal之间欧式距离的最小值和最大值分别为0.029 1和0.162 4。因此，a=0.030 5，b=0.090 6，c=0.029 1，d=0.162 4，即X0=〈0.030 5,0.090 6〉，X=〈0.029 1,0.162 4〉。

在确定经典域X0和节域X后，进一步根据式(13)和式(14)计算待测样本到Cnormal的欧氏距离与轴承正常状态距离区间的关联度，由此获得轴承性能退化曲线。由于轴承两向振动信号的获取不可避免的受环境影响，其二元MMFE变化不平稳，导致关联度的计算存在波动，如图10中实线所示。为减小关联度波动对轴承性能退化分析的影响，现利用四次多项式对关联度曲线进行拟合，获得关联度的变化趋势。拟合关系见式(16)，拟合结果如图10中虚线所示。

图10 轴承性能退化曲线Fig.10 The curve of bearing performance degradation

f(k)=-0.064 91×k4-0.194 9×k3-

0.113 5×k2-0.002 5×k+0.551 5

(16)

式中，k为关联度的计算值。

结合关联度定义和图10关联度拟合曲线知，从试验开始到第440个样本时，关联度大于0.5且其拟合曲线变化平稳，表明待测样本与Cnormal之间的欧式距离较小，待测样本较靠近Cnormal，其所代表的轴承运行状态良好。第440个样本到第540个样本时，关联度拟合曲线出现下降趋势，表明待测样本与Cnormal之间的欧式距离增大，样本逐渐远离Cnormal，轴承开始出现早期性能退化。由于关联度大于0，轴承仍处于正常运行状态。试验进行到第540个样本时，关联度为0，说明待测样本与Cnormal之间的欧式距离达到经典域的上限，即达到正常样本与Cnormal之间欧式距离的最大值。第540个样本后，拟合曲线的关联度值小于0，说明待测样本与Cnormal之间的欧式距离已超出其正常运行范围，轴承性能出现恶化。

4.3 时域特征对比分析

为说明本文所提方法对轴承性能退化程度的表征能力，提取轴承3两向振动信号的典型有量纲和无量纲时域特征，结果如图11所示。

(a) 峭度指标

由图11知，随着疲劳试验的进行，轴承故障程度逐渐加重，时域指标随之发生变化。峭度在第540个样本处发生明显波动，说明轴承发生故障，与图10所示结果一致，但峭度无法有效检测轴承的早期性能退化。平均值在轴承疲劳试验中变化较为平缓，仅在试验结束时才出现明显波动，说明平均值不能有效表征轴承的性能退化过程。波形因子的变化趋势与轴承性能退化趋势相一致，但也仅在试验末期才出现较大波动，无法检测轴承的早期性能退化。

综上分析，以上时域指标不能有效检测轴承的早期性能退化，且无法定量表述轴承的性能退化程度。相比而言，本文所提方法可以描述轴承的性能退化趋势，并能够定量表示其性能退化程度。由图11还可看出，上述时域指标在垂向和横向存在明显的相关关系，说明轴承两向振动信号之间存在较大的相关性，进行两向振动信号的综合分析能获得更准确的特征参数。

4.4 轴承早期性能退化分析

由前文分析知，轴承在第440个样本处开始出现早期性能退化，在第540个样本处发生内圈故障。通过计算得到轴承内圈故障的特征频率理论值为296.9 Hz，利用频谱分析法对以上两个样本的垂向振动信号进行功率谱分析，结果如图12所示。

(a) 第440个样本的降噪信号

由图12可知，第440个样本的特征频率为271 Hz，表明轴承开始出现早期性能退化。第540个样本的功率谱中出现了与轴承内圈故障频率相近的291 Hz及其倍频成分，说明轴承发生内圈故障，与本文分析结果相符。

5 结 论

(1) 将MMFE引入到轴承运行状态的特征提取。利用轴承垂向和横向振动信号构建了二元变量，并提取二元多尺度模糊熵。与轴承单向振动信号的MFE特征相比，二元MMFE能更有效的表征轴承运行状态。

(2)将可拓学关联函数与k-medoids聚类算法相结合构建了轴承性能退化评估模型。在对可拓集合和关联函数形式进行分析的基础上，利用待测样本到Cnormal的欧式距离与轴承正常状态距离区间之间的关联度实现了轴承性能退化的定量评估。与常用的时域指标相比，本文所提方法能有效检测轴承的早期性能退化，并能实现对轴承性能退化程度的定量表述。