CN110991422A - 基于多元时移多尺度排列熵的滚动轴承故障诊断方法 - Google Patents

Abstract

本发明公开一种基于多元时移多尺度排列熵的滚动轴承故障诊断方法，属于故障诊断技术领域。该方法首先采集待诊断物体的原始故障振动信号，然后提取原始故障振动信号的多元时移多尺度排列熵值，采用拉普拉斯分值法对多元时移多尺度排列熵进行降维，将降维后的故障特征样本分为多个训练样本和测试样本，然后采用多个训练样本对基于蝙蝠算法优化的支持向量机的多故障特征分类器进行训练，采用已训练完成的多故障特征分类器对测试样本进行分类，最后根据分类结果识别物体的故障类型以及程度。本发明在处理传感器采集信号的多通道信号中具有较高的创新性，在故障识别过程中具有较高的识别度。

Description

基于多元时移多尺度排列熵的滚动轴承故障诊断方法

技术领域:

本发明属于故障诊断技术领域，特别涉及基于多元时移多尺度排列熵的滚动轴承故障诊断方法。

背景技术：

设备故障诊断技术一直都是故障诊断与维护领域的研究重点，由于复杂的机械系统在运转过程中往往表现为非线性行为，振动信号往往呈现非线性和非平稳特性。传统的线性分析方法很难提取出隐藏在振动信号中的障特征信息。因此，非线性分析方法成为故障诊断技术领域相关学者和技术人员不断探索的热点。近年来，近似熵，样本熵，模糊熵，散布熵，排列熵(Permutation entropy,PE)等已经被广泛地应用到机械故障诊断领域，并取得了非常良好的故障诊断效果。

排列熵是最近提出的一种用于检测时间序列随机性和混乱性的方法，但排列熵局限于时间序列单一尺度下的分析，其他尺度上隐藏的丰富信息被忽略。在此种情况下，多尺度排列熵(Multi-scale Permutation Entropy,MPE)应运而生，多尺度排列熵被定义为不同尺度下的排列熵，能够衡量时间序列在不同尺度下的复杂性。但是，多尺度排列熵计算过程中存在不成熟的粗粒化过程，粗粒化过程中，随着尺度因子的增大，粗粒化得到的时间序列的长度会不断减少，进而造成得到的排列熵值的偏差不断增大，最终导致诊断结果的稳定性和精确性降低。另外，多尺度排列熵只能处理单通道的振动信号，对于多通道的振动信号束手无策。

发明内容：

本发明针对现有多尺度排列熵方法的不足，提供一种基于多元时移多尺度排列熵(Multivariate Time-Shifting Multi-scale Permutation Entropy,MTSMPE)滚动轴承的滚动轴承故障诊断方法。

首先，时移多尺度化时间序列能够克服多尺度排列熵中排列熵随着尺度因子的增大而偏差较大的问题；其次，多元思想的引入，能够处理传感器测取的多个通道的振动信息。最后，采用拉普拉斯分值法(Laplace Score,LS)对振动信号的特征进行降维，获取的特征集输入到蝙蝠算法优化的支持向量机(Bat Aalgorithm optimization Support VectorMachine,BA-SVM)中，对故障信息的类型和程度进行判断，避免了人工分类的繁琐和计算耗时等问题。本发明方法在特征提取过程中能够更好地提取故障振动信号的非线性故障特征，同时，能够在故障识别过程中也具有较高的故障识别度。

本发明提供的一种基于多元时移多尺度排列熵的滚动轴承故障诊断方法包括以下步骤：

(1)采集待诊断物体的原始故障振动信号；

(2)提取原始故障振动信号的多元时移多尺度排列熵；

(3)采用拉普拉斯分值法对多元时移多尺度排列熵进行降维，得到降维后的故障特征样本；

(4)将降维后的故障特征样本分为多个训练样本和测试样本；

(5)采用多个训练样本对基于蝙蝠算法优化的支持向量机)的多故障特征分类器进行训练；

(6)采用已训练完成的多故障特征分类器对测试样本进行分类；

(7)根据分类结果识别物体的工作状态和故障类型以及程度。

进一步地，步骤(2)中所测取原始故障信息的多元时移多尺度排列熵值的过程包括：

(2-1)对获取的原始故障振动信号进行多元时移粗粒化；

(2-2)计算同一尺度因子τ下生成的τ个符号序列的概率；

(2-3)对同一尺度因子下的所有符号概率求平均，通过香农熵的定义得到原始故障振动信号的多元时移多尺度排列熵值；

(2-4)对所有的尺度因子重复步骤(2-2)到(2-3)的操作，得到振动信号在所有尺度因子下的多元时移多尺度排列熵值。

进一步地，步骤(3)所述采用拉普拉斯分值法对多元时移多尺度排列熵进行特征降维的步骤包括：

(3-1)根据多元时移多尺度排列熵的最大尺度因子τ，构建一个含有τ个样本点的近邻图Q，判断样本点i与样本点j是否连通；

(3-2)若样本点i与样本点j不连通，令S_ij＝0；若样本点i与样本点j连通，则令

S_ij＝exp(-||x_i-x_j||²/c)

式中，c为常数，第i个样本点对应x_i，x_i是特征值元素，S_ij为加权矩阵S的元素；

(3-3)定义

f_r＝[f_r1,f_r2,…,f_rn]^T

其中，f_ri为第i个样本点的第r个特征值(i＝1,2,…,n)；T表示转置，D表示矩阵SI的对角矩阵，矩阵L为近邻图Q的拉普拉斯矩阵；

对各个特征值进行去均值化处理得到：

表示去均值后的特征值，I^T和f_r ^T分别表示I和f_r的转置；

(3-4)计算第r个特征值的拉普拉斯分值L_r：

其中，[N/τ]为第r个特征值的方差，

表示

的转置。

进一步地，步骤(3)中拉普拉斯分值进行从小到大排序，选取分值较小的前15个多元时移多尺度排列熵值作为故障特征值。

进一步地，步骤(5)中所述蝙蝠算法优化的支持向量机用于对故障特征样本中各样本的工作状态和故障类型进行分类，并分别根据已经训练完成的多故障特征分类器中的每单一蝙蝠算法优化的支持向量机的输出O(y)是否是+1进行判断；具体判断步骤包括：

(5-1)若输出是O(y)＝+1，则停止输入到下一个支持向量机，输出该测试样本集的分类；

(5-2)若输出是O(y)＝-1，则将该测试样本输入到下一个支持向量机，直到输出结果为+1时，输出测试样本的分类。

进一步地，步骤(2-1)中所述多元时移粗粒化过程包括:

(6-1)对于给定的尺度因子τ和时间序列

经过多元时移的处理，能够得到新的时间序列：

其中，τ是尺度因子，M是振动信号的通道数，k(1≤k≤τ)和β(β＝τ)是正整数，分别表示时间序列的起点和间隔时间。

是个四舍五入的整数，表示上边界个数。

(6-2)计算每个多元时移粗粒化序列

的排列熵值，同一尺度因子下熵值求平均得到多元时移多尺度排列熵：

其中，m为嵌入维度，t代表延迟时间。

进一步地，步骤(2-2)及步骤(2-3)具体步骤如下：

(7-1)假设对于给定原始时间数据长度为N的时移粗粒化序列Z：Z₁,Z₂,...Z_N，对其进行相空间重构，得到嵌入维度矩阵，其矩阵表达式包括：

其中，K＝N-m+1，K代表矩阵的行数，且每一行定为一模式向量Z(r)，故K也代表模式向量的个数；

(7-2)将每一个模式向量Z(r)的m个数据按照升序重新排列，即

z(r)＝{z(i+(j₁-1)t)≤z(i+(j₂-1)t)≤…≤z(i+(j_m-1)t)}

如果存在z(i+(j₁-1)t)＝z(i+(j₂-1)t)，按j值的大小来进行排序，即当j_k1<j_k2，有z(i+(j₁-1)t)≤z(i+(j₂-1)t),所以任何一个模式向量Z(r)都可以得到一组符号序列：S(g)＝{j₁,j₂,…,j_m}

其中，g＝1,2,…,k,K≤m！，m个不同的符号{j₁,j₂,…,j_m}共有m！种不同的排列，对应地，共有m！种不同的符号排列，s(g)是m！种符号序列中的一种；

(7-3)计算每一种符号序列出现的概率

此时，时间序列Z的排列熵按照香农熵的形式定义为：

当P_g＝1/m！时，H_p(m)达到最大值ln(m！)，因此，通过ln(m！)将排列熵H_p(m)进行标准化处理，即H_p＝H_p(m)/ln(m！)。

进一步地，蝙蝠算法优化支持向量机步骤包括：

(8-1)设置蝙蝠种群，种群大小为sol、维度为D、脉冲速度为R，脉冲频率为F，最大迭代次数为N；

(8-2)蝙蝠个体初始化。利用基本公式产生蝙蝠个体，构成初始蝙蝠种群；

(8-3)计算每个蝙蝠的适应度值并确定最优蝙蝠的位置。将所有蝙蝠个体作为支持向量机的参数在训练集上训练，并在测试集上测试，得到的数值作为对应蝙蝠的适应度值并且输出，根据返回的适应度值的大小，找到最优的蝙蝠个体；

(8-4)生成新蝙蝠个体。根据传统蝙蝠算法基本公式分别更新脉冲频率与飞行速度，生成新的蝙蝠个体；

(8-5)更新蝙蝠个体。产生一个随机数，如果随机数大于当前脉冲速度，则对当前群体中最优蝙蝠个体的邻域进行随机扰动，获得一个新的蝙蝠个体，而后用新蝙蝠个体替换当前蝙蝠个体，当前蝙蝠的平面位置坐标定义为(C_best,g_best)，输出最终的分类模型；；

(8-6)将测试样本输入到最优的蝙蝠算法优化支持向量机的模型，输出物体故障类型和程度。

本发明从故障物体振动信号中提取多元时移多尺度排列熵。多元时移多尺度排列熵是一种衡量时间序列非线性动力学行为和随机性突变的有效分析方法，能够有效地提取蕴藏在振动信号的更丰富、更全面的深层故障特征信息，在特征提取的过程中有较高的创新性。另外，基于多元时移多尺度排列熵的特征提取方法与基于蝙蝠算法优化的支持向量机相结合，在故障识别过程中能够得到较高的故障识别率。本发明具有以下显著技术特点：

(1)本发明所提的多元时移多尺度排列熵，主要利用了“多元+时移多尺度”的思想，克服了传统多尺度排列熵中粗粒化时间序列的缺陷。时移多尺度化的过程，使得到的时移粗粒化序列对数据长度N的依赖性大大降低(几乎不受影响)，多元主要解决是传感器测量得到的多通道信号。

(2)本发明所提蝙蝠算法优化的支持向量机是一种智能算法，能够较大程度地解决人工操作的繁琐和耗时问题，并且能够得到更高的识别率。

(3)本发明将基于多元时移多尺度排列熵的特征提取算法，基于拉普拉斯分值法的特征维度选择方法与基于蝙蝠算法优化的支持向量机的故障识别算法相结合，系统性地提出了一种新的故障诊断方法。

附图说明：

图1为本发明基于多元时移多尺度排列熵的滚动轴承故障诊断方法的流程图；

图2为本发明从振动信号中提取多元时移多尺度排列熵的流程图；

图3为本发明中蝙蝠算法优化支持向量机的工作原理流程图；

图4(a)为某大学滚动轴承故障数据的三通道时域波形图；

图4(b)为某大学滚动轴承故障数据的X通道时域波形图；

图4(c)为某大学滚动轴承故障数据的Y通道时域波形图；

图4(d)为某大学滚动轴承故障数据的Z通道时域波形图。

图5为13种状态下的振动信号的多元多尺度排列熵均值标准差曲线图；

图6为13种状态下的振动信号的多元时移多尺度排列熵均值标准差曲线图；

图7为基于多元移动多尺度排列熵，拉普拉斯分值与蝙蝠算法优化的支持向量机的滚动轴承故障方法识别率。

图8为对于LS前15个特征输入到BA-SVM的识别显示出13类状态样本的分类图。

具体实施方式：

本发明方法包括如下步骤：

(1)采集待诊断物体的原始故障振动信号；

(2)提取原始故障振动信号的多元时移多尺度排列熵；

(3)采用拉普拉斯分值法对多元时移多尺度排列熵进行降维，得到降维后的故障特征样本；

(4)将降维后的故障特征样本分为多个训练样本和测试样本；

(5)采用多个训练样本对基于蝙蝠算法优化的支持向量机)的多故障特征分类器进行训练；

(6)采用已训练完成的多故障特征分类器对测试样本进行分类；

(7)根据分类结果识别物体的工作状态和故障类型以及程度。

本发明在特征提取中有较高的创新性，在故障识别过程中具有较高的识别率。

由于多尺度排列熵在进行多尺度粗粒化过程中，对于不同尺度因子下的粗粒化序列，只采用了与起始点有关的序列，而忽略了该尺度因子下的其他序列。在不同的尺度因子下，由于粗粒化序列长度变短对熵值影响较大，其计算熵值存在必然的误差。

鉴于以上原因，为了克服多尺度排列熵值随着尺度因子增大而偏差较大的问题，本实施例创新性地应用所提出的多元时移多尺度排列熵作为原始故障振动信号的特征值，有效地抑制传统的粗粒化过程中时间序列变短对熵值曲线的影响，以振动信号X为例，若假设f_s和T₀分别为振动信号的采样频率和采样时长，令N＝f_s×T₀，则多通道的振动信号可以写成时间序列的形式:

K是振动信号的

通道数，即元数；N是时间序列的长度，采用如下的步骤对其进行多元时移化。

(2-1)对获取的原始故障振动信号进行多元时移粗粒化；

(2-2)计算同一尺度因子τ下生成的τ个符号序列的概率；

(2-3)对同一尺度因子下的所有符号概率求平均，通过香农熵的定义得到原始故障振动信号的多元时移多尺度排列熵值；

(2-4)对所有的尺度因子重复步骤(2-2)到(2-3)的操作，得到振动信号在所有尺度因子下的多元时移多尺度排列熵值。

步骤(2-1)中所述多元时移粗粒化过程包括:

(6-1)对于给定的尺度因子τ和时间序列

经过多元时移的处理，可以得到新的时间序列：

其中，τ是尺度因子，M是振动信号的通道数，k(1≤k≤τ)和β(β＝τ)是正整数，分别表示时间序列的起点和间隔时间。

是个四舍五入的整数，表示上边界个数。

(6-2)计算多元时移粗粒化序列

的排列熵，同一尺度因子下熵值求平均。

其中，m为嵌入维度，t代表延迟时间。

为了说明多元时移多尺度排列熵的优越性，本实施例以滚动轴承作为故障物体说明方法的有效性，并对不同工况的滚动轴承故障信号进行分析。实验采用某大学振动监测诊断中心的滚动轴承试验测试数据，轴承全寿命试验台为ID-25/30型，ID-25/30型轴承全寿命试验台包括驱动装置、支撑装置、加载装置及缓冲装置，测试轴承型号为SKF 6206-2Z深沟球轴承，通过电火花加工技术在滚动轴承上布置单点故障，信号采集使用DASP软件采集，采样频率为10240HZ，其实验数据描述如表1所示，每种状态取50组数据，每组数据长度为4096个数据点，13种状态下的三通道(X,Y,Z)原始信号时域波形图如图4(a)-图4(d)所示。

将上述13种故障滚动轴承的振动信号分别进行多元时移多尺度排列熵和多元多尺度排列熵(Multivariate Multi-scale Permutation entropy,MMPE)分析，计算出熵值后分别用LS进行特征排序，MTSMPE与MMPE经过LS排序的特征排序情况如表2与表3所示，特征排序后选择前15个特征绘制其均值标准差曲线如图5和图6所示。

表1某大学振动监测诊断中心滚动轴承试验测试数据描述

表2 LS排序后MMPE特征排序

表3 LS排序后MTSMPE特征排序

分析图5和图6可知，首先，分析排序后MMPE与MTSMPE的稳定性，利用MATLAB中errobar函数绘制的均值方差图明显的显示出MTSMPE要比MMPE要稳定很多，尤其表现在正常轴承，外圈故障轴承。其次，分析在同一特征下各个信号的MMPE，MTSMPE的分布密集性，MTSMPE的分布相对于MMPE的分布没有明显的优越性，但是值得注意的是，在较差的尺度特征上，MTSMPE相对于MMPE熵值的区分度越来越明显了，这就是说明在较差特征上以MTSMPE作为轴承故障特征表征的特征识别要优越于以MMPE作为轴承故障特征表征的特征识别。

在上述实验数据中，故障振动的信号共有13种，每种故障振动信号有50组，共计650组样本，现随机从每种故障振动信号的50组样本中选取20组样本作为训练样本，剩下的30组样本作为测试样本，共计260组训练样本与390组测试样本；分别将训练集与测试集的前1到前15个特征输入到BA-SVM中，得到基于MMPE与MTSMPE的故障识别率如图7所示,LS降维后的前15个特征特征输入到BA-SVM的故障识别率如图8所示。

从图7中可以看到，基于MTSMPE，LS与BA-SVM的故障识别率曲线一直在基于MMPE，LS与BA-SVM的故障识别率曲线之上，这说明充分了基于MTSMPE,LS与BA-SVM的故障诊断方法具有较高识别率，在故障识别方法领域较优越的性能。即便在大部分特征上识别率差别仅在1％之间，但是随着样本数据库的增大，这1％的误差就会避免很多样本分类的错误，所以不影响基于MTSMPE,LS与BA-SVM的故障诊断方法在此领域的优越性。图8中对于LS前15个特征输入到BA-SVM的识别进行详细作图，显示出13类状态样本都被正确的分类。

Claims (8)

1.基于多元时移多尺度排列熵的滚动轴承故障诊断方法，其特征在于该方法包括如下具体步骤：

(1)采集待诊断物体的原始故障振动信号；

(2)提取原始故障振动信号的多元时移多尺度排列熵；

(3)采用拉普拉斯分值法对多元时移多尺度排列熵进行降维，得到降维后的故障特征样本；

(4)将降维后的故障特征样本分为多个训练样本和测试样本；

(5)采用多个训练样本对基于蝙蝠算法优化的支持向量机)的多故障特征分类器进行训练；

(6)采用已训练完成的多故障特征分类器对测试样本进行分类；

(7)根据分类结果识别物体的工作状态和故障类型以及程度。

2.根据权利要求1所述的一种基于多元时移多尺度排列熵的滚动轴承故障诊断方法，其特征在于步骤(2)中所测取原始故障信息的多元时移多尺度排列熵值的过程包括：

(2-1)对获取的原始故障振动信号进行多元时移粗粒化；

(2-2)计算同一尺度因子τ下生成的τ个符号序列的概率；

(2-3)对同一尺度因子下的所有符号概率求平均，通过香农熵的定义得到原始故障振动信号的多元时移多尺度排列熵值；

(2-4)对所有的尺度因子重复步骤(2-2)到(2-3)的操作，得到振动信号在所有尺度因子下的多元时移多尺度排列熵值。

3.根据权利要求1所述的一种基于多元时移多尺度排列熵的滚动轴承故障诊断方法，其特征在于步骤(3)采用所述拉普拉斯分值法对多元时移多尺度排列熵进行特征降维的步骤包括：

(3-1)根据多元时移多尺度排列熵的最大尺度因子n，构建一个含有n个样本点的近邻图Q，判断样本点i与样本点j是否连通；

(3-2)若样本点i与样本点j不连通，令S_ij＝0；若样本点i与样本点j连通，则令

S_ij＝exp(-||x_i-x_j||²/c)

式中，c为常数，第i个样本点对应x_i，x_i是特征值元素，S_ij为加权矩阵S的元素；

(3-3)定义

f_r＝[f_r1,f_r2,…,f_rn]^T

其中，f_ri为第i个样本点的第r个特征值(i＝1,2,…,n)；T表示转置，D表示矩阵SI的对角矩阵，矩阵L为近邻图Q的拉普拉斯矩阵；

对各个特征值进行去均值化处理得到：

表示去均值后的特征值，I^T和f_r ^T分别表示I和f_r的转置；

(3-4)计算第r个特征值的拉普拉斯分值L_r：

其中，[N/τ]为第r个特征值的方差，

表示

的转置。

4.根据权利要求3所述的一种基于多元时移多尺度排列熵的滚动轴承故障诊断方法，其特征在于所述步骤(3)中拉普拉斯分值进行从小到大排序，选取分值较小的前15个多元时移多尺度排列熵值作为故障特征值。

5.根据权利要求1所述的一种基于多元时移多尺度排列熵的滚动轴承故障诊断方法，其特征在于步骤(5)中所述蝙蝠算法优化的支持向量机用于对故障特征样本中各样本的工作状态和故障类型进行分类，并分别根据已经训练完成的多故障特征分类器中的每单一蝙蝠算法优化的支持向量机的输出O(y)是否是+1进行判断；具体判断步骤包括：

(5-1)若输出是O(y)＝+1，则停止输入到下一个支持向量机，输出该测试样本集的分类；

(5-2)若输出是O(y)＝-1，则将该测试样本输入到下一个支持向量机，直到输出结果为+1时，输出测试样本的分类。

6.根据权利要求2所述的一种基于多元时移多尺度排列熵滚动轴承故障诊断方法，其特征在于：步骤(2-1)中所述多元时移粗粒化过程包括:

(6-1)对于给定的尺度因子τ和时间序列

经过多元时移的处理，能够得到新的时间序列：

其中，τ是尺度因子，M是振动信号的通道数，k(1≤k≤τ)和β(β＝τ)是正整数，分别表示时间序列的起点和间隔时间。

是个四舍五入的整数，表示上边界个数；

(6-2)计算每个多元时移粗粒化序列

的排列熵值，同一尺度因子下熵值求平均得到多元时移多尺度排列熵：

其中，m为嵌入维度，t代表延迟时间。

7.根据权利要求2所述的一种基于多元时移多尺度排列熵的滚动轴承故障诊断方法，其特征在于所述步骤(2-2)及步骤(2-3)具体步骤如下：

(7-1)假设对于给定原始时间数据长度为N的时移粗粒化序列Z：Z₁,Z₂,...Z_N，对其进行相空间重构，得到嵌入维度矩阵，其矩阵表达式包括：

其中，K＝N-m+1，K代表矩阵的行数，且每一行定为一模式向量Z(r)，故K也代表模式向量的个数；

(7-2)将每一个模式向量Z(r)的m个数据按照升序重新排列，即

z(r)＝{z(i+(j₁-1)t)≤z(i+(j₂-1)t)≤…≤z(i+(j_m-1)t)}

如果存在z(i+(j₁-1)t)＝z(i+(j₂-1)t)，按j值的大小来进行排序，即当j_k1<j_k2，有

z(i+(j₁-1)t)≤z(i+(j₂-1)t),所以任何一个模式向量Z(r)都可以得到一组符号序列：

S(g)＝{j₁,j₂,…,j_m}

其中，g＝1,2,…,k,K≤m！，m个不同的符号{j₁,j₂,…,j_m}共有m！种不同的排列，对应地，共有m！种不同的符号排列，s(g)是m！种符号序列中一种；

(7-3)计算每一种符号序列出现的概率

此时，时间序列Z的排列熵按照香农熵的形式定义为：

当P_g＝1/m！时，H_p(m)达到最大值ln(m！)，因此，通过ln(m！)将排列熵H_p(m)进行标准化处理，即H_p＝H_p(m)/ln(m！)。

8.根据权利要求1所述的一种基于多元时移多尺度排列熵滚动轴承故障诊断方法，其特征在于所述蝙蝠算法优化支持向量机的构建包括以下步骤：

(8-1)设置蝙蝠种群，种群大小为sol、维度为D、脉冲速度为R，脉冲频率为F，最大迭代次数为N；

(8-2)蝙蝠个体初始化，利用基本公式产生蝙蝠个体，构成初始蝙蝠种群；

(8-3)计算每个蝙蝠的适应度值并确定最优蝙蝠的位置。将所有蝙蝠个体作为支持向量机的参数在训练集上训练，并在测试集上测试，得到的数值作为对应蝙蝠的适应度值并且输出，根据返回的适应度值的大小，找到最优的蝙蝠个体；

(8-4)生成新蝙蝠个体，根据传统蝙蝠算法基本公式分别更新脉冲频率与飞行速度，生成新的蝙蝠个体；

(8-5)更新蝙蝠个体，产生一个随机数，如果随机数大于当前脉冲速度，则对当前群体中最优蝙蝠个体的邻域进行随机扰动，获得一个新的蝙蝠个体，而后用新蝙蝠个体替换当前蝙蝠个体，当前蝙蝠的平面位置坐标定义为(C_best,g_best)，输出最终的分类模型；

(8-6)将测试样本输入到最优的蝙蝠算法优化支持向量机的模型，输出物体故障类型和程度。

Images