CN107991097A - 一种基于多尺度符号动力学熵的轴承故障诊断方法 - Google Patents

Abstract

本发明提供一种基于多尺度符号动力学熵的轴承故障诊断方法，结合了多尺度、符号动力学与熵的概念，符号动力学熵是对时间序列复杂度的一种定量度量，也是时间序列一个重要的非线性特征，熵值会随着系统状态的变化而变化。此外，利用了多尺度熵的概念，计算不同尺度上的符号动力学熵，从不同的时间尺度上衡量信号的复杂度。将多尺度符号动力学熵应用到滚动轴承的故障诊断中，利用其提取滚动轴承的故障信息，同时结合分类器，准确识别出滚动轴承3种不同故障类型。

Description

一种基于多尺度符号动力学熵的轴承故障诊断方法

技术领域

本发明涉及数字信号处理领域。

背景技术

奥地利物理学家玻尔兹曼把熵与相容于每一个宏观态的微观状态数联系起来，为熵作出了微观的解释，即著名的玻尔兹曼公式。玻尔兹曼关系式不仅说明了微观状态数的物理意义，还给出了熵函数的统计解释(微观意义)。从这一式中可以看出，熵的大小由这一状态对应的微观态数目的多少来决定，熵的增加意味着系统包含的微观态的数目的增加。代表摘大小的微观状态数可以用粒子在空间分布的混乱度(无序度)的概念来解释，粒子越集中，数密度越大，即有序度越高则熵值越小，反之熵就越大。也就是说，微观态的多少反映了系统混乱度(无序度)的大小。因此，玻尔兹曼关系式表明，熵增加的过程就是系统混乱度(无序度)增大的过程，较小的熵就意味着系统的混乱度小，较大的摘则说明系统的混乱度较大。实际上，这一关系式揭示了熵的本质，即:熵代表了一个系统的混乱程度。玻尔兹曼从统计意义上对熵的微观解释使其物理意义更加明确，也为熵概念的泛化奠定了基础。近几十年来，人们对熵的应用早己超出热力学、物理学的范畴，不断泛化熵的概念和理论，并应用于信息论、数学、非线性动力学以及社会科学等各个领域中。

香农(C.E.Shamion)第一次将熵的概念引入到信息论当中，证明熵与信息内容的不确定程度有等价关系，并把这种不确定程度称为信息熵。信息熵为信息的量化度量建立了一个统一的计量方法，奠定了现代信息论的基础。一个系统越是混乱，信息熵就越大，也就是不确定性越大，包含的信息量越多，反之，信息熵就越小，不确定性越小，包含的信息越少。

滚动轴承是旋转机械中最常用且最重要的组件之一。由于滚动轴承在运行中的故障率较高，使得滚动轴承的故障诊断得到广泛重视。目前，常见的滚动轴承故障诊断方法主要包括振动信号分析，声学信号分析，轴承温度监测和铁谱分析法等。由于振动信号具有易于采集和分析的优点，基于振动信号的故障诊断方法应用最为广泛。

当滚动轴承发生故障时，测得振动信号呈现非平稳的特征。传统时频域的特征提取方法往往不能有效提取隐藏的故障特征，导致故障识别率降低。由于滚动轴承在不同位置发生故障时，故障冲击引起的系统响应不同，从而使得振动信号的混乱程度有所不同。基于熵值理论的分析方法可以不经过信号的分解或变换，直接度量信号的复杂度，完成滚动轴承不同故障位置的识别。目前，常用的熵值包括：样本熵(Sample entropy，SE)和排列熵(Permutation entropy，PE)。然而SE对于大数据的计算效率较低，而PE虽然能够提高计算效率，但是对噪声过于敏感，不能够准确提取轴承故障信息。因此，需要寻找和研究新的非线性动力学指标。

发明内容

本发明解决的技术问题是：本发明为了解决机械系统中信号出现非平稳非线性导致故障特征不明显的问题，提出了基于多尺度符号动力学和熵的信号特征提取方法。

本发明的技术方案是：一种基于多尺度符号动力学熵的轴承故障诊断方法，包括以下步骤：

步骤一：轴承的振动加速度信号对应为一个时间序列Y{y(i'),i'＝1,2,…,N}，其中i'表示时间序列对应某一时刻的数值，N表示时间序列的长度；该序列中的数值，按照尺度因子τ进行粗粒化分割，形成若干粗粒向量

其中τ为正整数(建议τ取值范围:1-20)；若干粗粒向量形成新的子时间序列,可将其进一步表述为X{x(i),i＝1,2,…,N₀}。其中i表示子时间序列在i时刻所对应

的数值，N₀表示分割后子时间序列的长度；

步骤二：对步骤一得到的每个子时间序列，分别求取符号动力学熵。

符号动力学熵的求解过程包括以下子步骤：

子步骤一：采用等概率区间划分准则，将子时间序列的幅值域划分成ε (3≤ε≤12)个区间，将ε称为符号数。用符号σ代替时间序列中元素的数值，从而形成符号序列Z{z(k),k＝1,2,…,N₀}，z(k)表示原始时间序列符号化后第k个数值所对应的符号σ；

子步骤二：人为给定的嵌入维数m和时间延迟λ，将符号序列

Z{z(k),k＝1,2,…,N₀}分割成一系列的嵌入向量，

其中

j＝1,2,...,N₀-(m-1)λ，用q^ε,m,λ表示子向量出现的模式种类，同时统计状态模式q^ε,m,λ在所有嵌入向量中出现的概率P(q^ε,m,λ)，并构建状态模式矩阵ε^m表示为状态模式数；

子步骤三：符号序列视为由一系列连续的状态模式组成的，相邻两个

状态模式可视为由上一个状态模式迁移到下一个状态模式；当已经观

测到状态模式q^ε,m,λ时，之后出现符号为σ的概率为状态迁移概率，

表示为P(σ|q^ε,m,λ)，计算公式如下：

P(σ|q^ε,m,λ)＝P{z(j+mλ)＝σ|j:j≤N-mλ,type(Z^ε,m,λ)＝q^ε,m,λ}

在状态迁移概率的基础上，构建ε^m×ε的状态迁移矩阵：

子步骤四：根据信息理论中香农熵的定义，定义符号动力学熵SDE等于

状态模式概率熵和状态迁移概率熵之和，计算公式如下：

式中，m 是嵌入维数，λ是时间延迟，ε是符号数；

子步骤五：符号动力学熵SDE(X,m,λ,ε)取得最大值为ln(ε^m+1)；将SDE 进行归一化：

SDE(X,m,λ,ε)＝SDE(X,m,λ,ε)/ln(ε^m+1)

步骤三：根据每个时间子序列的符号动力学熵，求取多尺度符号动力学熵 MSDE，公式如下：

步骤四：根据多尺度符号动力学熵MSDE，采用最小冗余最大相关方法(mRMR)对特征进行特征排序，选出最优特征构建成新的特征向量，输入到最小二乘支持向量机(LSSVM)中识别滚动轴承不同的故障类型。

发明效果

本发明的技术效果在于：本专利提出了一种新的复杂度计算方法—符号动力学熵(Symbol dynamic entropy,SDE)。SDE以符号动力学滤波为基础，将时间序列符号化，并构建状态模式矩阵和状态迁移矩阵，从而极大地保留了时间序列的状态信息。鉴于 SDE是从单一尺度上度量时间序列的复杂性，在SDE的基础上，再引入多尺度分析的概念，提出了多尺度符号动力学熵(Multi-scale symbolic dynamic entropy，MSDE)。 MSDE能够从不同尺度下衡量时间序列的动力学特性，极大地丰富了SDE的含义。

结合了多尺度、符号动力学与熵的概念，符号动力学熵是对时间序列复杂度的一种定量度量，也是时间序列一个重要的非线性特征，熵值会随着系统状态的变化而变化。此外，利用了多尺度熵的概念，计算不同尺度上的符号动力学熵熵，从不同的时间尺度上衡量信号的复杂度。把其应用到旋转机械的故障诊断中，利用熵值大小作为故障特征识别不同的故障类型及故障严重程度。

附图说明

附图1是本发明的方法流程图

附图2为粗粒化过程的示意图

附图3为实验一的轴承振动信号不同状态下的时域波形及频谱：(a)时域波形；(b)频谱

具体实施方式

参见图1，基于多尺度符号动力学和熵的信号特征提取方法，包括以下步骤：

步骤1、对于给定的时间序列Y{y(i'),i'＝1,2,…,N}，其中i′表示时间序列对应某一时刻的数值,按照粗粒化分割的方法，将原始的时间序列转化为多尺度的时间序列 X{x(i),i＝1,2,…,N₀}，其中N表示时间序列的长度，τ为尺度因子(建议τ取值范围:1-20；

步骤2、将得到的多尺度的时间序列X{x(i),i＝1,2,…,N₀},其中i表示子时间序列在i时刻所对应的数值，N₀表示分割后子时间序列的长度,转换成为符号序列(符号化)，即将步骤1中的时间序列用符号σ表示，从而形成符号序列Z{z(k),k＝1,2,…,N₀}，z(k)表示原始时间序列符号化后第k个数值所对应的符号σ；

步骤3、根据符号序列Z{z(k),i＝1,2,…,N₀}，确定嵌入维数m和时间延迟λ并构建一系列的嵌入向量其中j＝1,2,...,N₀-(m-1)λ；

步骤4、统计步骤3中嵌入向量可能出现模式种类q^ε,m,λ，并计算所有嵌入向量中每种模式出现的概率P(q^ε,m,λ)；

步骤5、统计在步骤4中每种模式q^ε,m,λ出现时，下一个符号σ出现的条件概率 P(σ|q^ε,m,λ)；

步骤6、计算步骤5中计算得到的每种模式概率P(q^ε,m,λ)和状态迁移概率P(σ|q^ε,m,λ)。

步骤7、根据信息理论中香农熵的定义，定义符号动力学熵SDE等于状态模式概率熵和状态迁移概率熵之和。

步骤8、对于步骤7中计算得到的符号动力学熵值SDE，进行归一化，得到归一化符号动力学熵值SDE(X,m,λ,ε)。

步骤9、将符号动力学熵扩展到多尺度领域，得到多尺度符号动力学熵 MSDE(y,τ,m,λ,ε)。

1.符号动力学熵

步骤1：将时间序列X{x(i),i＝1,2,…,N₀}，N₀表示为时间序列的长度转换为符号序列(简称：符号化)。

将时间序列的幅值域划分成ε个区间，将ε称为符号数。从而使得时间序列中的每一个元素根据其数值大小对应且仅对应了唯一的区间，用符号σ代替时间序列中元素的数值，从而形成符号序列Z{z(k),k＝1,2,…,N₀}。

步骤2：根据符号序列生成子向量，计算状态模式概率并构建状态模式矩阵。

根据嵌入维数m和时间延迟λ，可以将符号序列Z{z(k),k＝1,2,…,N₀}分割成一系列的子向量j＝1,2,...,N₀-(m-1)λ。每个子向量的符号排列模式都是唯一的，用q^ε,m,λ表示子向量出现的模式种类，同时统计状态模式q^ε,m,λ在所有子向量中出现的概率P(q^ε,m,λ)，并构建状态模式矩阵ε^m表示为状态模式数。

步骤3：计算状态迁移概率并构建状态迁移矩阵。

符号序列可视为由一系列连续的状态模式组成的，相邻两个状态模式可视为由上一个状态模式迁移到下一个状态模式。当已经观测到状态模式q^ε,m,λ时，之后出现的符号为σ的概率即为状态迁移概率，将其表示为P(σ|q^ε,m,λ)。

P(σ|q^ε,m,λ)＝P{z(j+mλ)＝σ|j:j≤N-mλ,type(Z^ε,m,λ)＝q^ε,m,λ} (1)

式中，ε是符号数，ε^m是状态模式数。

在状态迁移概率的基础上，构建ε^m×ε的状态迁移矩阵，状态迁移矩阵的形式为：

步骤4：根据信息理论中香农熵的定义，定义符号动力学熵SDE等于状态模式概率熵和状态迁移概率熵之和。

式中，m是嵌入维数，λ是时间延迟，ε是符号数。

步骤5：符号动力学熵SDE(X,m,λ,ε)取得最大值为ln(ε^m+1)。将SDE进行归一化，如式(4)所示：

SDE(X,m,λ,ε)＝SDE(X,m,λ,ε)/ln(ε^m+1) (4)

2.多尺度符号动力学熵

(1)长度为N的原始序列Y{y(i'),i'＝1,2,…,N}，对于给定嵌入维数m,时间延迟λ和符号数ε，进行粗粒化分割，从而形成新的粗粒向量如图1所示。

式中，τ为尺度因子(建议τ取值范围:1-20)。

(2)计算尺度时间序列的符号动力学熵SDE：

3.总体步骤及实现措施

基于MSDE和mRMR轴承不同故障位置特征提取方法的步骤如下：

(1)首先采用多尺度动力学熵提取轴承不同故障类型的特征，得到20个不同尺度下的符号动力学熵值。在本文中，设置嵌入维数m＝3、时间延迟λ＝1、符号数ε＝12和尺度τ＝20；

(2)采用mRMR对特征进行特征排序；

(3)选出最优的4个特征构建成新的特征向量；

(4)将得到的新特征向量输入到最小二乘支持向量机(LSSVM)中识别滚动轴承不同的故障类型。

基于MSDE、mRMR和LSSVM的滚动轴承磨损故障位置特征提取方法流程

为了更好地说明SDE的定义，下面给出一个SDE的计算实例。给定一个随机生成的时间序列X＝{-10,1,6,12,16,1,-14,-7,-10,24,-6,8,-1,9,-7,-14,-14,5,-1,-1}，假定符号数ε＝2，嵌入维数m＝2，时间延迟λ＝1，计算其符号动力学熵SDE。

首先，将时间序列符号化，生成符号序列。将时间序列的幅值域分割为两个区间，分别用符号α和β表示。采用等概率区间划分准则，将分割向量取值为-1，符号化后结果为：Z＝{α,β,β,β,β,β,α,α,α,β,α,β,β,β,α,α,α,β,β,β}。

其次，构建状态模式并计算状态模式概率。当符号数ε＝2，嵌入维数m＝2时，由符号序列可以生成N-(m-1)λ＝19个子向量，可能出现的状态模式有4种：

统计每种状态模式出现的概率，并构建状态模式矩阵：

然后，根据式(2-21)计算状态迁移矩阵：

最后，根据状态模式概率和状态迁移概率计算SDE值并归一化：

SDE(X,m,λ,ε)＝SDE(X,m,λ,ε)/ln(ε^m+1)＝0.9388 (2-24)

本专利效果可以通过以下实验进一步说明：

实验条件：本次实验总共选取不同故障尺寸和不同载荷下的三次故障模拟实验来验证本章所提方法的有效性。每一次实验选取同种故障尺寸和载荷下的轴承四种不同健康状态：正常，内圈故障，外圈故障，滚动体故障。

表1是本次故障模拟实验的实验条件数据。图1为轴承振动信号不同状态下的时域波形及频谱。当轴承出现不同故障类型，能够出现比较明显的故障冲击，然而仅仅从时域波形和频域上很区分出轴承具体的故障类型，因此有必要对其进行特征提取，来进一步确定轴承故障发生的具体位置。

表1.实验一数据列表

首先，采用MSDE对每一次实验的数据进行处理，提取故障特征。然后，采用 mRMR对特征进行优选，选取最优的4个特征输入到LSSVM中做模式识别。最后采用最小二乘支持向量机(LSSVM)进行模式识别。识别结果如表2所示。可以发现， MSDE能够有效地提取轴承振动信号的故障特征，结合LSSVM实现轴承不同故障位置的准确诊断。

表2.MSDE和LSSVM的方法识别率

Claims (1)

1.一种基于多尺度符号动力学熵的轴承故障诊断方法，其特征在于，包括以下步骤：

步骤一：轴承的振动加速度信号对应为一个时间序列Y{y(i'),i'＝1,2,…,N}，其中i'表示时间序列对应某一时刻的数值，N表示时间序列的长度；该序列中的数值，按照尺度因子τ进行粗粒化分割，形成若干粗粒向量

<mrow> <msubsup> <mi>X</mi> <mi>j</mi> <mi>&amp;tau;</mi> </msubsup> <mo>=</mo> <mfrac> <mn>1</mn> <mi>&amp;tau;</mi> </mfrac> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mrow> <mo>(</mo> <mi>j</mi> <mo>-</mo> <mn>1</mn> <mo>)</mo> </mrow> <mi>&amp;tau;</mi> <mo>+</mo> <mn>1</mn> </mrow> <mrow> <mi>j</mi> <mi>&amp;tau;</mi> </mrow> </munderover> <msub> <mi>y</mi> <mi>i</mi> </msub> <mo>,</mo> <mn>1</mn> <mo>&amp;le;</mo> <mi>j</mi> <mo>&amp;le;</mo> <mfrac> <mi>N</mi> <mi>&amp;tau;</mi> </mfrac> </mrow>

其中τ为正整数(建议τ取值范围:1-20)；若干粗粒向量形成新的子时间序列,可将其进一步表述为X{x(i),i＝1,2,…,N₀}。其中i表示子时间序列在i时刻所对应的数值，N₀表示分割后子时间序列的长度；

步骤二：对步骤一得到的每个子时间序列，分别求取符号动力学熵。

符号动力学熵的求解过程包括以下子步骤：

子步骤一：采用等概率区间划分准则，将子时间序列的幅值域划分成ε(3≤ε≤12)个区间，将ε称为符号数。用符号σ代替时间序列中元素的数值，从而形成符号序列Z{z(k),k＝1,2,…,N₀}，z(k)表示原始时间序列符号化后第k个数值所对应的符号σ；

子步骤二：人为给定的嵌入维数m和时间延迟λ，将符号序列

Z{z(k),k＝1,2,…,N₀}分割成一系列的嵌入向量，

<mrow> <msubsup> <mi>Z</mi> <mi>j</mi> <mrow> <mi>m</mi> <mo>,</mo> <mi>&amp;lambda;</mi> </mrow> </msubsup> <mo>{</mo> <mi>z</mi> <mrow> <mo>(</mo> <mi>j</mi> <mo>)</mo> </mrow> <mo>,</mo> <mi>z</mi> <mrow> <mo>(</mo> <mi>j</mi> <mo>+</mo> <mi>&amp;lambda;</mi> <mo>)</mo> </mrow> <mo>,</mo> <mo>...</mo> <mo>,</mo> <mi>z</mi> <mrow> <mo>(</mo> <mi>j</mi> <mo>+</mo> <mo>(</mo> <mrow> <mi>m</mi> <mo>-</mo> <mn>1</mn> </mrow> <mo>)</mo> <mi>&amp;lambda;</mi> <mo>)</mo> </mrow> <mo>}</mo> </mrow>

其中j＝1,2,...,N₀-(m-1)λ，用q^ε,m,λ表示子向量出现的模式种类，同时统计状态模式q^ε,m,λ在所有嵌入向量中出现的概率P(q^ε,m,λ)，并构建状态模式矩阵ε^m表示为状态模式数；

子步骤三：符号序列视为由一系列连续的状态模式组成的，相邻两个

状态模式可视为由上一个状态模式迁移到下一个状态模式；当已经观

测到状态模式q^ε,m,λ时，之后出现符号为σ的概率为状态迁移概率，

表示为P(σ|q^ε,m,λ)，计算公式如下：

P(σ|q^ε,m,λ)＝P{z(j+mλ)＝σ|j:j≤N-mλ,type(Z^ε,m,λ)＝q^ε,m,λ}

在状态迁移概率的基础上，构建ε^m×ε的状态迁移矩阵：

子步骤四：根据信息理论中香农熵的定义，定义符号动力学熵SDE等于状态模式概率熵和状态迁移概率熵之和，计算公式如下：

式中，m是嵌入维数，λ是时间延迟，ε是符号数；

子步骤五：符号动力学熵SDE(X,m,λ,ε)取得最大值为ln(ε^m+1)；将SDE进行归一化：

SDE(X,m,λ,ε)＝SDE(X,m,λ,ε)/ln(ε^m+1)

步骤三：根据每个时间子序列的符号动力学熵，求取多尺度符号动力学熵MSDE，公式如下：

<mrow> <mi>M</mi> <mi>S</mi> <mi>D</mi> <mi>E</mi> <mrow> <mo>(</mo> <mi>y</mi> <mo>,</mo> <mi>&amp;tau;</mi> <mo>,</mo> <mi>m</mi> <mo>,</mo> <mi>&amp;lambda;</mi> <mo>,</mo> <mi>&amp;epsiv;</mi> <mo>)</mo> </mrow> <mo>=</mo> <mi>S</mi> <mi>D</mi> <mi>E</mi> <mrow> <mo>(</mo> <msubsup> <mi>X</mi> <mi>j</mi> <mi>&amp;tau;</mi> </msubsup> <mo>,</mo> <mi>m</mi> <mo>,</mo> <mi>&amp;lambda;</mi> <mo>,</mo> <mi>&amp;epsiv;</mi> <mo>)</mo> </mrow> </mrow>

步骤四：根据多尺度符号动力学熵MSDE，采用最小冗余最大相关方法(mRMR)对特征进行特征排序，选出最优特征构建成新的特征向量，输入到最小二乘支持向量机(LSSVM)中识别滚动轴承不同的故障类型。