赵建岗， 宁 静,2， 宁云志， 陈春俊， 李艳萍

(1. 西南交通大学 机械工程学院， 成都 610031； 2. 轨道交通运维技术与装备四川省重点实验室， 成都 610031)

过去的三十年中，基于信号处理方法的特征提取算法在机械领域被广泛研究[1]，在文献[2]中，用近似熵(approximate entropy,ApEn)监测滚动轴承的健康状况。然而ApEn严重依赖于数据长度，并且其估计值总是低于处理短数据集时的预期值。为了克服这种不足，提出了样本熵(sample entropy,SampEn)[3]，由于SampEn值对数据长度不敏感并且不受信号中的噪声影响，因此它受到了很多关注。然而ApEn和SampEn仅在单一尺度上估计复杂性，当他们应用于多时间尺度时效果极差。为克服这个弱点，在文献[4]中提出了多尺度熵(multiscale entropy, MSE)。但MSE对故障诊断复杂性的估计不准确。另外用它处理长时间序列时特别耗时[5]。

近年来，引入排列熵(permutation entropy,PE)[6]来计算机械系统的复杂性[7-8]。但是与AnEn和SampEn一样，PE只能通过单一尺度分析时间序列，即可能会丢失嵌入其他尺度的许多有用的故障信息。但是某些地区机车和机床的操作环境也非常苛刻和复杂[9-10]，监控系统收集的大量信号通常是非静态和非线性信号。由于机械组件之间的相互作用和耦合效应，这些获得的信号包含多种自然振动模式，这导致单尺度特征提取方法不能理想地表达这些信号[11]。为克服这一缺点，提出了多尺度排列熵(multiscale permutation entropy,MPE)[12]来估算不同尺度上时间序列的复杂性。在文献[13]中MPE用于提取滚动轴承的故障特征，并且验证了MPE的性能优于PE。

通常在使用MPE进行特征提取之后，特征矩阵的维度将非常高。提高故障识别的准确性，需要降低高维特征矩阵的维数以获得低维主要特征，即降维[14]。主成分分析(principal component analysis,PCA)[15]是经典的降维算法，但其仅对线性和高斯数据有效，并且在处理非线性或非平稳信号时经常不能达到理想的效果。最近流形学习被用来提取嵌入在非线性高维数据中的低维内在结构，而且它已被用于机械故障诊断[16-17]。其中以线性局部切线空间排列(linear local tangent space alignment,LLTSA)[18]性能较好，本文采用LLTSA进一步提取嵌入在由MPE计算的高维特征矩阵中的低维主要特征。

丁吉等[19]用同步压缩小波变换对发电机的连杆轴承进行时频分析，为发电机相关性能和工况的保障提供了重要的估计方法。孙自强等[20]用小波分形方法计算了风机组轴承时域波形的关键维数，可以确定轴承的故障形式。向丹[21]等利用LLTSA对特征进行降维，并通过SVM对故障类别进行辨识。苏祖强等[22]通过LLTSA与加权k最近邻分类器提高了故障分类的精度。熊庆等[23]研究了多重分形趋势波动分析特征参数的敏感性和稳定性，并成功运用于滚动轴承的定量故障诊断。李永健等[24]通过改进的多尺度排列熵对列车的轴箱轴承进行了诊断研究，取得了较好的效果。

本文将MPE算法在多尺度信息处理上的优越性和LLTSA在非线性数据强适应性相结合应用于轴承的故障诊断，用两个轴承箱数据对MPE-LLTSA方法进行了测试，测试结果表明了该方法在机械模式分类和故障识别领域的应用潜力。

1 理论背景

1.1 MPE的简介

1.1.1 MPE算法

多尺度排列熵(MPE)被定义为不同尺度时间序列的一组PE值。MPE算法包括两个步骤：① 应用粗粒度过程从原始时间序列获得多尺度时间序列;② 计算每个粗粒度时间序列的PE值。这两个步骤可以简要总结如下[25]。

(1)

(2) 在MPE分析中，每个粗粒度时间序列的PE值由式(1)～(5)计算并且比例因子s的函数可以表示为

(2)

1.1.2 MPE的参数选择

在MPE算法中，需要设置四个参数，嵌入维度m，时间延迟τ，时间序列长度N和比例因子s。其中嵌入维度m决定输入的状态数量m!，并且PE值高度依赖于m的选择。如果m太小该方法将无法工作，因为不同状态太少。但是如果m太大，由于条件N≥5m!，可能会对随后的多尺寸研究造成相当大的不便[26]。Bandt等指出尺寸m应满足3≤m≤7。为了检测信号的动态变化，通常通过信息损失和计算复杂性之间的权衡关系来选择嵌入维数m。因此本文中m设为4。时间延迟τ对结果几乎没有影响，所以我们在本文中将τ设置为1。时间序列N的长度对PE值的估计有很大影响，N越大计算时间越久。而太小的N不能满足条件N≥5m!。考虑到这些约束，我们在第3章的模拟试验中将数据长度N设置为400，将比例因子s设置为7。在第4章的研究中，N设置为1 200，s设置为12。

1.2 线性局部切线空间排列

线性局部切线空间排列(LLTSA)[27]是一种经典的流形学习算法，它可以提取高维数据集的固有几何信息，即挖掘隐藏在高维观察空间中嵌入的低维流形，从而不会丢失数据中的重要信息[28-29]。

假设一个高维数据集X={x∈(R)D,i=1,2,…,N}，它是从属于(R)d的较低维特征空间的基础流形中采样的，其中D和d(d

Y=ATXHN

(3)

式中:HN=I-eeT/N表示中心矩阵;I表示单位矩阵;e表示所有元素等于1的列向量。 LLTSA的基本假设是每个样本附近的局部结构信息可以用其局部切线空间表示。因此为构建从高维输入空间到低维特征空间的投影矩阵，可以在全局低维特征空间中重新排列所有样本的局部切线空间。

在LLTSA算法中，根据样本之间的欧几里德距离构造采样点xi(i=1,2,…,N)的邻域xi=[xi1,xi2,…,xik]，其中k表示最近样本的数量。然后用局部变换矩阵Qi来将邻域xi映射到局部低维切线空间。因此样本点xi附近的局部结构可以近似表示为

(4)

获得局部结构后，为得到X的全局低维表示Y，将所有样本的局部切向空间重新排列在全局低维特征空间中，选择矩阵为S=[S1,S2,…,SN]其中Si(i=1,2,…,N)是0-1的选择向量，因此Yi=YSi，其中Yi=[yi1,yi2,…,yik]是Xj的全局低维表示。此步骤的目标函数可以转换为以下最小问题

(5)

(6)

式中,B=SWWTST。用ATXHNBHNXTA=Id确定唯一的Y。最后式(4)的解转化为广义特征值问题为

ATXHNBHNXTα=λXHNXTα

(7)

可以由对应于式(7)的特征向量λ1≤λ2≤…≤λd的d个特征向量α1,α2,…,αd得到转换矩阵A，即A=α1,α2,…,αd。

2 基于MPE-LLTSA的特征提取方法

故障诊断中，对于经常出现在机械领域中的非线性信号，特征选择是一个难题。EEMD因为其自适应性而被广泛用于处理这种非平稳信号。流形学习在非线性信号中具有更强的降维能力[30]，因此本文引入了流形学习算法LLTSA。本文将该特征提取方法缩写为EEMD-LLTSA, EEMD-LLTSA的技术路线及其识别方法如图1(a)所示。

(a) EEMD-LLTSA

然而在复杂的机械设备中，监测系统收集的大量非线性信号通常包含多种自然振动模式，这导致单尺度特征提取通常无法表征这些信号。因此引入多尺度方法MPE来改善传统基于单尺度的非线性动态参数的性能。首先通过使用MPE计算信号x(t)以获得一组熵值{HMN1,HMN2,…,HMNs},其中s是比例因子。然后构建具有高维度的样本矩阵V，用LLSTA挖掘高维样本特征矩阵V中的固有结构以适应低维特征矩阵。最后用LSSVM建立诊断模型进行状态分类。本文将该特征提取方法缩写为MPE-LLTSA，MPE-LLTSA特征提取方法及其条件识别的技术路线如图1(b)所示。详细步骤描述如下：

步骤1参数选择：时间序列x(i),k=1,2,…,N，相应参数有嵌入维数m，时间延迟τ，时间序列N的长度和比例因子s。

步骤2MPE算法计算信号的多尺度熵值为

(8)

步骤3对所有样本进行步骤1和2，构建高维样本特征矩阵V为

(9)

式中,p代表样本数量。

步骤4通过流形学习降低维数：使用LLTSA降低样本特征矩阵V的维数以获得低维特征向量，如下面的公式所示

(10)

式中:fn表示第n维特征;d是嵌入维度。

步骤5训练和识别：将获得的低维流形特征分为训练样本和测试样本，并使用LSSVM实现状态识别。

3 试验探究

3.1 模拟试验

3.1.1 信号模拟

非线性信号通常具有调幅[31]的特性，这经常出现在机械领域。根据振动信号构造两个调幅信号来验证MPE-LLTSA特征提取方法的性能。

w(t)=am·[1+mcos(2πfmt)]·sin(2πfct)

(11)

式中，am=1，m=8,fm=0.3 Hz，fm=0.3 Hz。对于第一信号w1(t)，fc,1=3，对于第二信号w2(t)，fc,2=3.1 Hz，采样频率为20 Hz。

下面用随机生成的混合矩阵模拟三个复杂信号。

(12)

通过上述矩阵M对w1(t)和w2(t)执行矩阵运算

(13)

添加三个随机高斯噪声n1(t),n2(t) 和n3(t)以分别生成新信号a(t),b(t)和c(t)，如图2所示。

(a)

(14)

3.1.2 特征评估

如图2所示，箭头处的那三个模拟信号的幅度是不同的，但差异不是那么明显。根据它们的波形无法准确区分这三个信号。为了通过频率特征区分三个信号，将FFT方法引入上述三个模拟信号，分别生成其频域图，如图3所示。可以看出，通过FFT获得的这三个频域图之间的差异并不明显。通过FFT，我们同样无法准确地区分这些有调制特性的非平稳信号。最后引入Hilbert-Huang变换(HHT)来分析这三类信号的时频特性，谱图的结果如图4所示。它表明这些简单的时频分析方法无法准确识别出三个信号。

(a)

(a)

上述时频分析方法无法准确识别这三种信号。为了解决这个问题，引入了MPE-LLTSA特征提取方法。为每种类型的信号选择60个样本，总共180个。每个样本的数据长度N设置为400，时间长度为20秒。如1.1节所示，在MPE算法中，分别设置三个参数，嵌入维数m，时间延迟τ和比例因子s。在本次研究中，我们设定m=

4,τ=1和s=7。如1.2节所示，在算法LLTSA中，根据文献[32]，两个参数邻居数k和内在维度d分别设置为7和3。a(t),b(t)和c(t) 的每个样本的MPE结果显示在图5中。通过MPE计算所有样本以重建180×12高维特征矩阵。然后通过使用LLTSA方法获得180×3的低维特征矩阵。降维后的结果如图6(a)所示。可以看出这三个模拟信号可以通过MPE-LLTSA算法完全分离，得到的不同信号的三维特征高度集中，聚类效果极佳。为了证明LLTSA方法的优越性，还引入了PCA方法。聚类效果如图6(b)所示。结果表明，MPE-PCA特征提取方法 还可以准确识别三个模拟信号，聚类效果也非常好。它表明MPE特性可以有效地克服信号中的噪声影响，适用于非线性信号处理。

图5 三个模拟信号的MPE特征

(a) MPE-LLTSA

为了证明MPE-LLTSA优于EEMD-LLTSA方法，引入了机械中常用的能量和SampEn功能进行仿真分析。根据文献[32]在EEMD算法中集合数和信噪比(SNR)两个参数分别设置为100和0.3。首先所有样本都通过EEMD进行分解，α(t),b(t)和c(t)的每个样本的结果如图7所示。由于EEMD是主成分分析方法，因此信号的主要信息包含在前几个IMF中。此外所有样本获得的IMF数量须大于6。因此为了确保EEMD-LLTSA方法的后续处理，提取前六个IMF的能量和SampEn特征以分别构造180×6高维特征矩阵。然后通过LLTSA算法减少特征矩阵以生成180×3的低维特征矩阵(参见图8)。从图7和图8可知，EEMD-energy-LLTSA特征提取方法可以分离三类信号，但聚类效果明显差于MPE-LLTSA和MPE-PCA。EEMD- SampEn-LLTSA方法甚至无法分离这三个模拟信号，因为在特征提取过程中，单尺度特征提取方法仅考虑信号的一个单尺度信息。在某种程度上，它忽略了其他维度中包含的故障信息。这种单尺度特征提取方法经常无法在具有严重故障信息耦合的机械设备实现期望的结果。

图7 模拟信号α(t),b(t)和c(t)的EEMD结果

(a) EEMD-energy-LLTSA

3.1.3 聚类评价

为了准确地评估聚类能力，引入了三个评价指标，即类内散射Sw，类间散射Sb和比率Sw/Sb以描述聚类效应。对于特征向量[f1，f2,…,fd]，d是嵌入维度，参数Sw和Sb的定义为

(15)

(16)

通过计算类内离散度Sw，类间离散度Sb和比率Sw/Sb，可以获得3.1.2节中实现的特征的准确评估，如表1所示。结果表明，MPE-LLTSA方法具有良好的分类和聚类能力，优于MPE-PCA和EEMD-LLTSA方法。

表1 对三种信号的不同特征进行聚类评估

3.2 机械故障诊断试验结果

从滚动轴承系统测量的信号非常复杂。为了进一步验证MPE-LLTSA方法的有效性，在这一节中将使用Case Western Reserve University提供的滚动轴承的试验数据[33]。轴承试验系统及其示意图如图9所示。它由2马力的马达(左)，扭矩传感器(中心)，测功机(右)和控制电子设备(未示出)组成。试验轴承是6205-2RS JEM SKF深沟球轴承，局部单点损坏分别使用电火花加工机设置在外滚道，内滚道和滚动体的位置。本文中的振动信号取自安装在驱动端的加速度计，采样频率为12 kHz。

图9 滚动轴承试验系统草图

3.2.1 四类情况

这种情况下，分析包括健康，内在缺陷(IRD)，滚动元件缺陷(RED)和外在缺陷(ORD)在内的四类情况。本节中的振动信号从损坏0.178 cm的轴承中收集。数据长度N设置为1 200。训练和测试样本的数量在表2中给出。引入MPE-LLTSA来识别这四个信号。如1.1.3节所述，在MPE方法中我们设置m=4，τ=1和s=12。四个信号的波形如图10(a)及其MPE结果如图10(b)所示。对所有样品通过MPE构建高维样本特征矩阵之后用LLTSA降维得到低维特征。在LLTSA方法中，我们设置k=7和d=3。计算结果如图11(a)所示，表明MPE-LLTSA可以完全分离这四个信号。为了验证我们提出的MPE-LLTSA方法的优越性，还计计算了MPE-PCA和EEMD-LLTSA结果，如图11(b)，图11(c)和图11(d)所示。为得到精确的评估，还计算了类内离散度Sw，类间离散度Sb和Sw/Sb，如表3中所示。最后引入LSSVM方法识别四种状态，正确的识别率如表3所示。从图11和表3可知MPE-LLTSA方法可以区分这四种状态，且其聚类效果优于MPE-PCA和EEMD-LLTSA。

表2 不同检测的试验数据的详细描述

(a)

(a) MPE-LLTSA

表3 不同检测下不同特征的聚类评价

3.2.2 五种状态

在表4所示的情况下，分析了滚动体的五种不同的滚动轴承故障(健康的，0.178 cm，0.356 cm，0.533 cm和1.02 cm)，详细信息显示在表4中。数据长度N选择为1 200，训练和测试样本的数量分别为50和50。在MPE方法中，我们设置m=4，τ=1和s=12。在LLTSA方法中，我们设置k=7和d=3。五个信号的波形图如图12(a)及其MPE结果如图12(b)所示。MPE-LLTSA的计算结果如图13(a)所示。为了做比较，在图13(b)，图13(c)和图13(d)中列出了MPE-PCA和EEMD-LLTSA的结果。为了评估更精确，还在表5中计算了三个参数Sw，Sb和Sw/Sb。最后引入LSSVM方法识别五种状态，并在表5中列出正确的识别率。从图13和表5，我们可以看到MPE-LLTSA可以完全分离出五类故障，分类结果和聚类效果优于MPE-PCA和EEMD-LLTSA。测试结果表明，MPE-LLTSA方法可用于轴承缺陷严重性评估。

表4 五个检测严重程度的试验数据的详细描述

(a)

(a) MPE-LLTSA

3.3 讨 论

我们提出的MPE-LLTSA方法的优点主要是由于以下两个原因：

(1) 能量和SampEn等单尺度特征提取方法只能分析单个尺度的时间序列，导致隐藏在其他尺度中的有用信息丢失。然而非线性或非平稳信号通常包含多种自然振荡模式，并且故障信息通常是复杂的并且总是隐藏在多尺度时间序列中。MPE方法可以估算不同尺度的时间序列的复杂度，可以有效地克服信号中的噪声影响。

(2) 在模式识别领域，由EEMD，MPE等计算出的特征矩阵的维数通常很高，难以直接识别这些高维特征。为了提高识别精度，有必要进一步提取嵌入在高维特征空间中的低维主要特征。LLTSA具有良好的减小高维非线性数据维度的能力。

4 结 论

本文提出了一种用于机械故障诊断的MPE-LLTSA特征提取方法。由于机械设备结构复杂，操作环境恶劣，监测系统采集的大量非线性信号通常包含多种自然振动信号，单尺度特征提取方法往往无法表征这些信号。为了解决这个问题，引入MPE来提取非线性信号的多尺度特征并构建高维特征矩阵。随后用LLTSA挖掘嵌入的内在结构并进行低维特征提取。最后用LSSVM进行机械故障诊断。同时还引入了一些特征提取方法，包括MPE-PCA和EEMD-LLTSA进行比较。仿真和实际案例的测试结果表明，我们提出的MPE-LLTSA方法具有良好的分类和聚类能力，与其他方法的比较结果证明了我们提出的MPE-LLTSA方法在机械故障诊断领域的价值和应用潜力。