基于MIGA-VMD和t-SNE的轴承故障诊断方法

王双海，米大斌，芦浩，姜文，龚思远，梁涛

(1.河北建投能源投资股份有限公司，河北石家庄 050011;2.河北建设投资集团有限责任公司，河北石家庄 050011;3.河北工业大学人工智能与数据科学学院，天津 300131)

0 前言

汽轮机是火力发电的动力机械，而轴承又是汽轮机的重要组成部件。轴承一旦发生故障，必将影响整个机械设备的运行，甚至造成经济损失以及人员伤亡。所以，在实际生产中及时诊断出轴承的运行状态至关重要。当轴承发生故障时，振动信号特性也会随之改变，若能提取出各状态的信号特征，则可进行轴承故障诊断。

文献[1]中采用经验模态分解(EMD)把原始信号分解为多个平稳的固有模态函数，并将其峭度和裕度因子作为神经网络的输入参数来识别滚动轴承故障模式。文献[2]中将EMD和奇异值分解相结合提取轴承故障特征向量，并用模糊C均值聚类来识别滚动轴承故障类型。文献[3]中针对EMD的模态混叠，采用集成经验模态分解(EEMD)来提取发动机曲轴故障特征，最后识别出不同种类的故障。文献[4]中针对EMD方法存在端点效应、模态混叠和虚假分量的问题,采用局部均值分解(LMD)和样本熵相结合来提取轴承故障特征向量，并用贝叶斯网络来识别滚动轴承的故障类型。虽然以上两种方法相对EMD效果要好，但依然无法从根本上解决模态混淆和端点效应等问题。2014年, DRAGOMIRETSKIY提出了一种新型的自适应信号处理方法——变分模态分解(Variational Mode Decomposition,VMD)方法[5]，作为一种新型的自适应信号处理方法，VMD方法引入变分模型,将信号的分解转换为约束模型最优解的寻优问题,可以避免端点效应、抑制模态混淆,并且具有很高的分解效率[6]。

排列熵作为一种非线性理论的信号分析方法，能够较好地度量时间序列的复杂性，可以将排列熵作为不同复杂度信号的特征参数。高维特征向量包含了大量故障特征信息，但是维数过高不仅会导致计算的复杂度过高的问题，还会对分类器性能造成不良影响[8]。因此，需要利用降维方法进行二次特征提取,去除冗余特征。传统的降维方法如主成分分析法(PCA)、线性判别分析法(LDA),虽然对线性结构的高维数据有很好的降维效果，但对实际获取的非线性结构数据，降维效果却是差强人意。近年来新提出的降维方法t分布随机邻近嵌入(t-SNE)是一种基于流行学习的非线性数据降维方法，具有优异的降维能力[8]。

本文作者以轴承振动信号为研究对象，提出了基于多岛遗传算法优化VMD参数、排列熵和t-SNE的特征提取方法。仿真结果表明，该方法可有效提高轴承故障诊断的正确率。

1 VMD简介

变分模态分解过程实质上是构造变分问题与求解变分问题的过程，从构造和分解两个方面介绍VMD方法。

1.1 变分问题的构造

VMD算法将分解后的IMF(Intrinsic Modal Functions)分量重新定义为一个调幅-调频信号，假设原始信号被分解为K个IMF分量，则第k个IMF分量的表达式为

uk(t)=Ak(t)cos[φk(t)]

(1)

式中:k∈{1,…,K};φk(t)是一个非递减的相位函数；Ak(t)表示包络函数。

首先，对每一个IMF分量，利用Hilbert变换求其解析信号以获得单边频谱，其频谱表达式为

(2)

其次，根据混合预估中心频率，调制频谱至对应基频带，记为

(3)

最后，计算解调信号的时间梯度L2范数的平方，估计出模态分量的带宽，并引入约束条件，构造出约束的变分模型为

(4)

式中：{uk}代表VMD分解得到的k个IMF分量；{ωk}代表IMF分量对应的中心频率；δ(t)为狄利克雷函数；⊗为卷积运算；x代表原始信号。

1.2 变分问题的求解

引入二次惩罚因子α和拉格朗日乘法系数λ(t)，将约束变分问题变为非约束变分问题。其中二次惩罚因子可在高斯噪声存在的情况下保证信号的重构精度，拉格朗日系数使得约束条件保持严格性，增广的拉格朗日函数如下：

L({uk},{ωk},λ)=〈λ(t),x(t)-

(5)

然后用交替方向乘子方法(ADMM)对上式进行极值求解得到IMF分量的频域表达式如下所示：

(6)

(7)

(8)

通过以上分析,VMD的算法流程如下：

(2)根据式(6)和式(7)更新uk和ωk；

(3)根据式(8)更新λ；

2 MIGA优化VMD参数

在VMD分解过程中，分解的IMF分量个数K和惩罚参数α对分解结果有较大影响。若K值过大，分解结果中会产生虚假分量，误导对实验结果的判断；若K值较小，则会产生模态混叠现象。若惩罚参数α越小，分解后各个IMF分量的带宽就越小；相反，若惩罚参数α越大，分解后各个IMF分量的带宽就越大。除K、α以外，其他参数对分解效果影响较小，设置为经验值，即tau=0,init=1,DC=0,ε=1×10-7。因此，在对信号进行VMD分解前，需要选取适当的K和α。

在VMD参数的优化过程中，采用粒子群优化算法对VMD的参数进行优化[9]。但是粒子群算法容易陷入局部最优，搜索路径也比较复杂。“遗传算法”的概念是由密西根大学的HOLLAND于1975年首次提出的，模仿自然界中“优胜劣汰，适者生存”的法则进行优化过程[7]。在参数优化过程中，采用遗传算法对VMD的参数进行优化。它通过染色体选择、交叉、变异等迭代操作得到参数的优化解。该算法具有隐式并行搜索的特点，在求解过程中具有很强的全局搜索能力，但是遗传算法容易过早收敛和陷入局部最优。针对遗传算法经常出现的过早收敛、产生局部最优解的问题，有学者提出了多岛遗传算法(Multi-Island Genetic Algorithm, MIGA)。多岛遗传算法通过反复恰当地使用遗传算法的算子和选择原则，从亲代到子代，从子代到孙代，从孙代到重孙代，不停地繁衍，使得种群对环境的适应性不断升高[10]。本文作者采用多岛遗传算法进行参数优化，获取最优的参数组合。多岛遗传算法寻优的具体步骤如下：

步骤1，初始化种群P0。

步骤2,将初始种群P0划分为多个“岛”，即子种群。

步骤3,计算各个“岛”上每个个体的适应度值。

步骤4,在每个岛上进行标准遗传算法中的选择、交叉和变异等操作。

步骤5,若满足迁移条件，则从一个岛迁移到另一个岛。

步骤6,若当前迭代次数n小于最大迭代次数MaxGen，则继续进行步骤3,否则进行步骤7。

步骤7,输出种群中适应度值最优的个体作为最优解。

多岛遗传算法的参数设计如下：初始种群划分为10个岛，每个岛的种群数设为10，优化10代，则计算个体数为1 000个。变异概率设为0.01，交叉概率设为1，迁移概率设为0.01，精英个体数量设为1。

为了验证MIGA的收敛性和优化性能，测试了不同的算法，包括PSO算法、GA算法、MIGA算法。通过以下函数作为适应度函数比较它们的优化性能

(9)

对于PSO算法，最大迭代次数设为50，种群规模为100。对于遗传算法，最大迭代次数设为50，种群规模为100，交叉概率为1，变异概率为0.01。对于多岛遗传算法，初始种群划分为10个岛，每个岛的种群数设为10，最大迭代次数为50，交叉概率为1，变异概率为0.01，迁移概率为0.01，精英个体数量为1。优化过程的迭代曲线如图1所示。可以看出：对于等式(9)中的适应度函数，PSO算法的收敛速度优于GA算法，但是PSO算法和GA算法都陷入了局部最优解；文中所用的MIGA算法不仅收敛速度优于PSO算法和GA算法，并且获得的最优值更接近理论最优值。结果表明：MIGA算法具有较快的收敛速度、较强的全局优化能力。

图1 优化过程的迭代曲线

用多岛遗传算法进行参数优化时，需选定一个适应度函数。文献[11]中提出了包络熵Ep的概念，故障信号经VMD分解后，IMF分量的包络熵值Ep可反映出该分量的稀疏特性。如果分解后IMF分量中包含的噪声较多，从而掩盖故障冲击特征，则该IMF分量的稀疏性较弱，包络熵较大；相反，若IMF分量中包含较多的规律性故障冲击，则该IMF分量具有较强的稀疏性，包络熵较小。在K与α的影响下，K个分量就会有K个包络熵值，选择K个包络熵值中最小的一个作为局部最小包络熵minEp，该最小熵值对应的分量有着丰富的故障特征信息[6]。将局部最小熵作为多岛遗传算法中的适应度函数，整个搜索过程就是要找到全局最小包络熵以及对应的最佳分量所在的最优参数组合[K,α]。所以，构建如下的适应度函数：

minF=minEp

(10)

式中包络熵Ep的计算公式为

(11)

式中:j=1,2,3,…,N;a(j)是x(j)经Hilbert解调得到的包络信号；ej是对a(j)进行归一化得到的；Ep是根据信息熵的计算规则得到的。

为了验证局部包络熵作为适应度函数的有效性，构造了一个滚动轴承的模拟信号。其表达式如下：

(12)

式中:c(t)表示加噪声的轴承模拟信号;y0为位移常数;fn为轴承的固有频率;ξ为阻尼系数;n(t)为趋于真实噪声的高斯白噪声;t为采样时间;采样频率fs=20 kHz，采样点数量N=4 096。这里，设置y0=3，fn=3 000 Hz，ξ=0.09。

图2展示了在无噪声、噪声强度为1、噪声强度为2和噪声强度为3的情况下模拟信号的波形。通过4幅图的比较，可以发现信号中的周期性脉冲完全被强背景噪声所淹没；噪声强度与包络熵呈正相关：随着噪声的增强(减弱)，周期脉冲变得模糊(显著)，包络熵变大(变小)。这也证明了信号稀疏度越大、包络熵越小的规律。因此，包络熵最小的模态分量可以确定为最佳分量。

图2 不同噪声强度下模拟信号的时域波形

3 相关系数准则

相关系数是反映变量之间相关关系密切程度的统计指标[16]。文中采用皮尔逊相关系数分析各IMF分量与原始信号的相关性程度。设样本X和样本Y，则两者的相关系数为

(13)

式中：r为相关系数；Cov(X,Y)为样本X和样本Y的协方差；D(X)为样本X的方差；D(Y)为样本Y的方差。

相关系数r的值越大，代表两个样本的相关性越高。通过计算IMF分量与原始信号的相关系数来选取最能代表轴承振动信号的部分分量。

4 排列熵

4.1 排列熵原理

排列熵是描述时间序列复杂性和不规律程度的一种参数，具有计算简单、抗噪能力强等优点[14]。排列熵的算法简介如下。

对于长度为n的时间序列{[x(i),i=1,2,…,n]}，参考相空间延迟坐标法对其进行重构，每个采样点取其连续的m个样本点，得到点x(i)的m维重构向量[12]：

X(i)={x(i),x(i+τ),…,x[i+(m-1)τ]}

(14)

其中:m是嵌入维数；τ是延迟时间；i=1,2,…,n-(m-1)τ,即X(i)向量的个数为n-(m-1)τ。将X(i)中的元素按从大到小的顺序排列，即

x[i+(j1-1)τ]≤x[i+(j2-1)τ]≤…≤

x[i+(jm-1)τ]

(15)

式中:j1、j2,…,jm为X(i)中各元素的位置。

如果X(i)中出现大小相等的两个元素，如x[i+(j1-1)τ]=x[i+(j2-1)τ]，则按照j的大小来排序，即x[i+(j1-1)τ]≤x[i+(j2-1)τ]。

所以任一向量X(i)一定有一种符号序列：

S(l)=(j1,j2,…,jm)

(16)

其中:l=1,2,…,k且k≤m!。根据排列组合定理，m个元素最多有m!种排列方式，S(l)只是其中的一种。若S(l)这种排列方式出现的概率为Pl，则序列{x(i),i=1,2,…,n}的排列熵可定义为

(17)

若H(m,τ)越小说明时间序列越规则，若H(m,τ)越大说明时间序列无序程度越高，信号复杂度越大。

4.2 参数选择

排列熵的计算值大小主要受3个参数的影响，分别是信号序列长度n、嵌入维度m、延迟时间τ。作者依据文献[12]分别设置嵌入维度m=6，延迟时间τ=1。

4.3 几组典型信号的排列熵值

为了证明排列熵具有表征信号特征的能力，下面计算6组具有代表性信号的排列熵值。

x1为长度为3 000、服从均匀分布白噪声信号。

x2(t)=0.5t2+0.9t-0.3

x3(t)=sin(2π·900t)

x4(t)=sin(2π·30t)

x5(t)=sin(2π·30t)·sin(2π·90t)

x6(t)=0.5πt+5

经计算，x1白噪声的排列熵值最大为6.456 2；其次x5高频正弦信号排列熵值为3.021 4；x3、x4非高频正弦信号、调幅信号为1.099 2，1.671 8；x2、x6时域的曲线和直线信号排列熵值最小都为0。由此分析，随着信号复杂度越大，排列熵值越大，排列熵可以作为信号分类的特征。

5 t-SNE流行学习算法

t分布随机邻近嵌入(t-SNE)算法是一种深度学习的非线性流行学习算法，对高维非线性数据集具有优异的降维效果[13]。t-SNE算法具体步骤[17]如下:

(1)设高维样本X={x1,x2,…,xn}，任意两个样本数据点xi和xj相似的条件概率为

(18)

式中：δi为以xi为中心点的高斯分布方差，可根据用户指定的困惑度Perp经二分搜索确定，困惑度定义为

(19)

(2)计算高维样本数据的联合概率密度：

(20)

(3)初始化低维数据解y0={y1,y2,…,yn}。

(4)计算低维数据相似度：

(21)

(5)使用KL散度度量高维数据分布P和低维数据分布Q的相似度：

(22)

(6)计算梯度

(23)

(7)获取低维数据：

(24)

式中：η为学习率;α(t)为优化参数;t代表第t次迭代。

(8)重复步骤(4)—(7),直到满足迭代次数T。经T次迭代得到低维数据y(T)={y1,y2,…,yn}。

6 基于MIGA-VMD和t-SNE的滚动轴承故障诊断流程

基于上述理论，基于MIGA-VMD分解、排列熵和t-SNE的轴承故障诊断模型如图3所示。具体步骤如下：

图3 故障诊断流程

(1)获取各状态的振动信号。

(2)用MIGA算法搜索VMD方法在各状态下的最佳参数组合[K0,α0],采用经参数优化的VMD方法分解各状态下的振动信号得到K0个IMF分量。

(3)计算各IMF分量与原始信号的相关性系数。

(4)计算与原始信号相关性较高的6个IMF分量的排列熵，组成特征向量。

(5)利用t-SNE方法降维得到三维特征向量。

(6)将训练集的特征向量输入到SVM进行训练，得到SVM分类模型。

(7)将测试集样本的特征向量输入到训练好的SVM模型，实现故障诊断。

7 仿真实验

7.1 实验数据描述

采用凯斯西储大学电气工程实验室的数据进行仿真。选取轴承的负载为0，转速为1 797 r/min,损伤直径为0.177 8 mm,采样频率为12 kHz的正常状态、内圈故障、滚珠故障、外圈故障振动信号。每种状态各取前10 s的数据，样本长度设为2 000，即每种状态各有60组样本，共240组样本，其中160组样本作为训练样本，80组样本作为测试样本。

7.2 故障特征提取

首先，利用MIGA算法搜索VMD算法的最优参数组合。以内圈故障信号为例，VMD参数寻优过程的迭代曲线如图4所示。经过8次迭代搜索到了全局的最优解，此时的包络熵值最小，为4.238。该全局最优解对应的参数α=2 999、K=6，将其引入到VMD算法的参数设置中。图5、图6为对内圈故障信号进行VMD处理得到的6个模态分量对应的信号波形和频谱。

图4 寻优过程的迭代曲线

图5 内圈故障信号波形

对4种状态信号分别进行VMD参数优化，对应的最佳参数组合如表1所示。

根据表1中的参数组合[K0,α0]设置VMD方法的K、α，利用参数优化的VMD方法分解样本得到多个IMF分量。

图7为4种状态下各分量与原始信号的相关程度。

图7 4种状态下各IMF分量与原始信号的相关程度

选取与原始信号相关性较高的6个IMF分量并计算其排列熵做为特征向量。

7.3 故障特征简约

利用t-SNE方法对六维的特征向量R240×6进行降维得到低维特征向量R240×3，为了验证t-SNE的降维效果，同时采用PCA、LDA、t-SNE 3种方法进行降维处理，降维结果如图8所示。

图8 3种方法降维结果

从图8可以看出：在3种降维方法中，经t-SNE降维后4种类型的数据特征聚集性最好，特征区分最为明显。因此，本文作者采用t-SNE方法进行降维，得到三维的特征向量。

7.4 故障分类

将160组训练样本的特征向量输入到SVM分类器进行训练，得到用于故障诊断的SVM模型。通过训练好的SVM模型对4种状态下的80组测试样本分类，测试样本的分类结果如图9所示。在SVM的训练和测试中，采用数字标签代表轴承的4种运行状态：正常状态(标签1)，内圈故障状态(标签2)，滚珠故障状态(标签3)，外圈故障状态(标签4)。

图9 测试样本的分类结果

由图9可以看出：采用基于MIGA-VMD和排列熵、t-SNE的特征提取方法，轴承状态的诊断精度可达到100%。

7.5 方法对比

传统的VMD参数设置方法参考文献[15]的中心频率选择法，通过观察不同模态数K下对应的中心频率变化，来确定K的取值。若中心频率变化较大，说明出现了欠分解，则增大K值；若中心频率变化较小，说明出现了过分解，则减小K值。

为证明所提方法的有效性，选用同样的样本数据，分别采用EMD、EEMD、LMD、传统VMD(K=6，α=2 000)、MIGA-VMD 5种方法进行信号分解，同样采用排列熵构建特征向量，并用t-SNE降维处理得到低维特征向量来训练SVM模型，用训练好的SVM模型对4种状态的测试样本分类。5种方法得到的低维特征分布与SVM分类准确率如图10、表2所示。

由图10、表2可以看出:基于MIGA-VMD和排列熵、t-SNE的特征提取方法，诊断正确率优于其他4种方法，可以从干扰信号中提取出微小故障特征，实现对不同故障正确诊断识别。

图10 低维特征分布

表2 5种方法的SVM分类准确率

8 结论

提出基于MIGA-VMD和排列熵、t-SNE的特征提取方法,借助参数优化的VMD方法分解振动信号,计算分解后与原始信号相关性较高的部分IMF分量的排列熵作为故障特征，利用t-SNE方法进行降维得到低维特征向量并将其输入SVM实现轴承的故障诊断。经MIGA算法优化参数的VMD方法与EMD、EEMD、LMD、传统VMD等方法相比,提取的故障特征更为丰富,有利于对轴承状态的诊断识别。t-SNE非线性降维方法与传统线性降维方法PCA、LDA相比，在轴承诊断中具有更佳的降维效果。仿真实验表明,采用此MIGA-VMD和排列熵、t-SNE的特征提取方法,轴承的故障诊断精度有明显的提高。