CN102866027A - 基于lmd和局域时频熵的旋转机械故障特征提取方法 - Google Patents

Abstract

本发明公开了一种基于LMD和局域时频熵的旋转机械故障特征提取方法，其技术方案的要点是，它包括如下步骤：1.利用加速度传感器对旋转机械设备进行测量，获得振动加速度信号；2.对振动加速度信号进行LMD分解，得到若干PF分量，求出各分量的瞬时幅值和瞬时频率；3.做出时频谱图，划分时频平面，计算局域时频熵；4.用局域时频熵值作为特征量，结合实验来提取故障特征。本发明实现了基于LMD的旋转机械故障诊断系统的分析过程，研究不同状态下设备的振动信号在时频分布上能量分布的差异，将局域时频熵理论引入机械故障诊断，将不同状态下的振动信号经LMD变换后再进行局域时频熵计算，并以此作为特征量来判断设备是否存在故障。

Description

基于LMD和局域时频熵的旋转机械故障特征提取方法

技术领域

本发明涉及一种机械工程领域的机械故障诊断方法，具体说本发明是一种基于LMD和局域时频熵的旋转机械故障特征提取方法。

背景技术

现今，工业生产逐步向大型化、高速化、自动化和智能化方向迈进，生产企业所使用的主要设备中，旋转类设备约占80％，这些设备能否正常运行关系到企业的巨大经济利益，如若某台设备出现故障又未能及时发现和排除，可能带来巨大的安全隐患，甚至灾难性的后果。因此，研究和应用旋转机械状态监测与故障诊断技术，对于保障生产安全、避免事故和巨额经济损失、提高设备安全管理水平，具有重要的意义。

旋转机械设备运行时，其振动信号一般是非常复杂的，各种振动因素综合起来，得到的机械系统的振动信号必然是非平稳非线性的多分量信号，特别是在出现故障的情况下，具有明显的非线性非平稳特性，其频谱随时间有较大变化，而且不同的非平稳特性也预示着不同的机械故障形式。目前，对这类非平稳信号的分析处理方法很多，常见的时频分析方法有短时傅里叶变换、Wigner分布、小波变换、经验模态分解(Empirical ModeDecomposition,简称EMD)等，但它们都有各自的局限性，如短时傅里叶变换的时频窗口大小是固定不变的；Wigner分布对多分量信号进行分析时会产生交叉项；小波变换虽具有可变的时频窗口^[9]，但和短时傅里叶变换一样是对时频面的机械格型分割，本质上它不是一种自适应的信号处理方法；EMD是一种自适应的信号处理方法，自提出后在机械故障诊断等很多领域中都得到了应用，但在理论上还存在一些问题，如EMD方法中的过包络、欠包络、模态混叠、端点效应、利用Hilbert变换计算瞬时频率产生无法解释的负频率等问题，均处在研究之中。

发明内容

本发明的目的是克服上述方法在机械故障诊断中的不足，提供一种基于LMD和局域时频熵的旋转机械故障特征提取方法。

为了实现上述目的，本发明的技术方案是：一种基于LMD和局域时频熵的旋转机械故障特征提取方法，其步骤包括：

步骤一：利用加速度传感器对旋转机械设备进行测量，获得振动加速度信号；

步骤二：对振动加速度信号进行LMD分解，得到若干PF分量，求出各分量的瞬时幅值和瞬时频率；

步骤三：做出时频谱图，划分时频平面，计算局域时频熵；

步骤四：用局域时频熵值作为特征量，结合实验来提取故障特征。

上述基于LMD和局域时频熵的旋转机械故障特征提取方法，在所述步骤二中，对振动加速度信号进行LMD分解的过程步骤如下：

1)求局部均值函数m₁₁(t)，找出加速度信号x(t)所有的局部极值点n_i，求出所有相邻的局部极值点的平均值，然后将所有相邻的平均值点用直线连接起来，并用滑动平均法进行平滑处理，得到m₁₁(t)；

2)求包络估计函数a₁₁(t)，计算相邻局部极值点的包络估计值，将所有相邻两个包络估计值a_i用直线连接，然后采用滑动平均方法进行平滑处理，得到a₁₁(t)；

3)将局部均值函数m₁₁(t)从原始信号x(t)中分离出来，得到滤掉低频信号的h₁₁(t)；

4)用h₁₁(t)除以包络估计函数a₁₁(t)以对h₁₁(t)进行解调，得到s₁₁(t)；

5)对s₁₁(t)重复步骤1)~4)，直至满足1-Δ≤a_1n(t)≤1+Δ（变动量Δ<1），则迭代过程终止，否则需要继续重复步骤1)~4)；

6)把步骤1)~4)迭代过程中产生的所有包络估计函数相乘便可以得到包络信号a₁(t)；

7)将包络信号a₁(t)和纯调频信号s_1n(t)相乘便可以得到原始信号的第一个PF分量PF₁(t)；

8)将第一个PF分量PF₁(t)从原始信号x(t)中分离出来，得到一个新的信号u₁(t)，将u₁(t)作为原始数据x(t)重复步骤1)~7)，循环k次，直到u_k为一个单调函数为止，则原始信号x(t)分解为k个PF分量和一个单调函数u_k之和。

上述基于LMD和局域时频熵的旋转机械故障特征提取方法，在所述步骤三中，将时频平面等分为K×L个面积相等的时频块，按照不同频段将整个时频平面划分为若干个面积相等的时频段，将每一个时频段的时频熵定义为局域时频熵，第h时频段的局域时频熵定义为：

$s_h\left(q\right)=-\underset{i=M}{\overset{N}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{L}{\mathrm{\&Sigma;}}}{q}_{i,j}\ln q_{i,j}$

其中，i＝1,…,K，j=1,…,L，q_i，j为第i×j小块平面的归一化能量。

本发明提供的一种基于LMD和局域时频熵的旋转机械故障特征提取方法的有益效果是：

1、本发明采用局部均值分解方法对选择机械振动加速度信号进行分解，自适应地将任何一个复杂的非平稳信号分解成多个瞬时频率具有物理意义的PF分量之和，其中每一个PF分量都是由一个包络信号和一个纯调频信号相乘得到，包络信号是该PF分量的瞬时幅值，而PF分量的瞬时频率则可由纯调频信号直接求出，将所有PF分量的瞬时幅值和瞬时频率组合便可以得到原始信号完整的时频分布；

2、局域时频熵给信号的时频分布的精细划分提供了一个量的度量，按照不同频段将整个时频平面划分为若干个面积相等的时频段，将每一个时频段的时频熵定义为局域时频熵。通过计算局域时频熵，可以判断出信号所代表的机械运行状态。

附图说明

图1为时频平面等分示意图；

图2为正常轴承振动加速度信号时域波形图；

图3为正常轴承局部均值分解结果图；

图4为正常轴承时频谱图；

图5为外圈故障轴承振动加速度信号时域波形图；

图6为外圈故障轴承局部均值分解结果图；

图7为外圈故障轴承时频谱图。

具体实施方式

下面结合附图与具体实施方式对本发明作进一步详细描述：

旋转机械故障诊断过程中首先需要利用加速度传感器对机械设备进行测量，获得振动加速度信号x(t)，再对振动加速度信号进行分解，提取特征值。本发明利用LMD方法对振动加速度信号进行分解。

本发明的一种基于LMD和局域时频熵的旋转机械故障特征提取方法，其具体步骤如下：

步骤一：利用加速度传感器对旋转机械设备进行测量，获得振动加速度信号；

步骤二：对振动加速度信号进行LMD分解，得到若干PF分量，包括以下步骤：

(1)求局部均值函数m₁₁(t)。找出原始信号x(t)所有的局部极值点n_i，求出所有相邻的局部极值点的平均值

$m_i=\frac{n_i+n_{i+1}}{2}---\left(1\right)$

其中，i＝1,2，...，M；M为原始信号的局部极值点的个数。然后，将所有相邻的平均值点m_i用直线连接起来，并用滑动平均法进行平滑处理，得到m₁₁(t)。

(2)求包络估计函数a₁₁(t)。包络估计值为

$a_i=\frac{|n_i-n_{i+1}|}{2}---\left(2\right)$

将所有相邻两个包络估计值a_i用直线连接，然后采用滑动平均方法进行平滑处理，得到a₁₁(t)。

(3)将局部均值函数m₁₁(t)从原始信号x(t)中分离出来，得到

h₁₁(t)=x(t)-m₁₁(t)                （3）

(4)用h₁₁(t)除以包络估计函数a₁₁(t)以对h₁₁(t)进行解调，得到

s₁₁(t)=h₁₁(t)/a₁₁(t)              （4）

理想地，s₁₁(t)是一个纯调频信号，对s₁₁(t)重复上述步骤便能得到s₁₁(t)的包络估计函数a₁₂(t)，假如a₁₂(t)不等于1，说明s₁₁(t)不是一个纯调频信号，需要重复上述迭代过程n次，直至s_1n(t)为一个纯调频信号，也即s_1n(t)的包络估计函数a_1(n+1)(t)＝1，所以有

$\left\{\begin{array}{c}{h}_{11}\left(t\right)=x\left(t\right)-m_{11}\left(t\right)\\ h_{12}\left(t\right)=s_{11}\left(t\right)-m_{12}\left(t\right)\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ h_{1n}\left(t\right)=s_{1(n-1)}\left(t\right)-m_{1n}\left(t\right)\end{array}\right.---\left(5\right)$

式中，

$\left\{\begin{array}{c}{s}_{11}\left(t\right)=h_{11}\left(t\right)/a_{11}\left(t\right)\\ s_{12}\left(t\right)=h_{12}\left(t\right)/a_{12}\left(t\right)\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ s_{1n}\left(t\right)=h_{1n}\left(t\right)/a_{1n}\left(t\right)\end{array}\right.---\left(6\right)$

迭代终止的条件为

$\underset{n\&RightArrow;\&infin;}{\lim }{a}_{1n}\left(t\right)=1---\left(7\right)$

实际应用中，在不影响分解效果的前提下，为了减少迭代次数，降低运算时间，可以设定一个变动量Δ，用

1-Δ≤a_1n(t)≤1+Δ                    （8）

作为迭代终止的条件。

(5)把迭代过程中产生的所有包络估计函数相乘便可以得到包络信号(瞬时幅值函数)

$a_1\left(t\right)=a_{11}\left(t\right)a_{12}\left(t\right)\&CenterDot;\&CenterDot;\&CenterDot;a_{1n}\left(t\right)=\underset{q=1}{\overset{n}{\mathrm{\&Pi;}}}{a}_{1q}\left(t\right)---\left(9\right)$

(6)将包络信号a₁(t)和纯调频信号s_1n(t)相乘便可以得到原始信号的第一个PF分量

PF₁(t)＝a₁(t)s_1n(t)            （10）

它包含了原始信号中最高的频率成分，是一个单分量的调幅—调频信号，其瞬时幅值就是包络信号a₁(t)，其瞬时频率f₁(t)则可由纯调频信号s_1n(t)求出，即

$f_1\left(t\right)=\frac{1}{2\mathrm{\&pi;}}\frac{d\left[\arccos \left(s_{1n}\left(t\right)\right)\right]}{\mathrm{dt}}---\left(11\right)$

(7)将第一个PF分量PF₁(t)从原始信号x(t)中分离出来，得到一个新的信号u₁(t)，将u₁(t)作为原始数据重复以上步骤，循环k次，直到u_k为一个单调函数为止。

$\left\{\begin{array}{c}{u}_1\left(t\right)=x\left(t\right)-{\mathrm{PF}}_1\left(t\right)\\ u_2\left(t\right)=u_1\left(t\right)-{\mathrm{PF}}_2\left(t\right)\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ u_k\left(t\right)=u_{k-1}\left(t\right)-{\mathrm{PF}}_k\left(t\right)\end{array}\right.---\left(12\right)$

原始信号x(t)能够被所有的PF分量和u_k重构，即

$x\left(t\right)=\underset{p=1}{\overset{k}{\mathrm{\&Sigma;}}}{\mathrm{PF}}_P\left(t\right)+u_k\left(t\right)---\left(13\right)$

说明LMD分解没有造成原信号信息的丢失。至此，将所有PF分量的瞬时幅值和瞬时频率组合便可以得到原始信号x(t)完整的时频分布。

步骤三：做出时频谱图，划分时频平面，计算局域时频熵；

时频熵给不同状态下信号的时频分布提供了一个量的度量，通过计算不同信号的时频熵，可以判断出信号所代表的机械运行状态。但是，实际振动信号的频率成分比较多，得到的时频分布比较复杂，如果对比整体的时频熵值可能不足以细致地反映系统真实的运行状态，基于此，提出了局域时频熵的概念。将时频平面等分为K×L个面积相等的时频块，如图1所示，每块内的能量为W_i，j(i＝1,…,K，j＝1,…,L)，整个时频平面的能量为A，对每块进行能量归一化，得到q_i，j＝W_i，j/A，于是就有

$\underset{i=1}{\overset{K}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{L}{\mathrm{\&Sigma;}}}{q}_{i,j}=1---\left(14\right)$

符合计算信息熵的初始归一化条件。然后，按照不同频段将整个时频平面划分为若干个面积相等的时频段，将每一个时频段的时频熵定义为局域时频熵，仿照信息熵的计算公式，第h时频段的局域时频熵定义为：

$s_h\left(q\right)=-\underset{i=M}{\overset{N}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{L}{\mathrm{\&Sigma;}}}{q}_{i,j}\ln q_{i,j}---\left(15\right)$

计算能量集中时频段内的局域时频熵s₁(q),s₂(q),…,s_k(q)，而忽略能量分布较稀疏的时频段。这样，一方面提高时频分辨率，细致地反映了各时频区域能量分布的差别；另一方面可以减小计算量，提高运算时间。

步骤四：用局域时频熵值作为特征量，结合实验来提取故障特征。

设定轴承故障模拟试验台的转子转速为1725r/min，模拟轴承外圈单点电蚀故障，故障尺寸为0.007inches，采样频率为2048Hz，轴承外圈故障特征频率为103.3Hz，转轴基频为28.75Hz。采集一组正常轴承振动加速度信号和一组轴承外圈故障振动加速度信号，分别如图2、图5所示，可以看出故障信号振幅增大，且出现周期性冲击成分。

本发明中采用局部均值分解方法对正常轴承振动加速度信号进行分解，得到4个PF分量和1个余量，如图3所示，做出其时频谱，如图4所示，振动信号的能量在整个时频平面内比较分散。将采集到的故障振动加速度信号进行如上处理，得到5个PF分量和1个余量，如图6所示，做出其时频谱，如图7所示，振动信号的能量集中在故障特征频率（103.3Hz）及其二倍频（206.6Hz）、转轴基频的20倍频（575Hz）频段处。

分别将图6、图7所示时频平面等分为262144个时频块，在能量比较集中的103.3Hz、206.6Hz及575Hz附近划分时频段，分别命名为时频段1~3，计算局域时频熵值，如表1所示。正常轴承的时频熵值比故障状态下振动信号的时频熵值大，与时频分析中能量分布相对应，这与最大熵原理是一致的；进一步，相比于正常状态，轴承的故障特征频率（时频段1）及其二倍频（时频段2）附近的局域时频熵值明显减小，可以判断轴承外圈发生了故障。对于时频段3，该局域时频熵值（为5.3418）占总时频熵值的82.73％，这是轴承由于本身结构特点在正常运转时产生振动和噪声。说明在强背景噪声的干扰下，本发明提出的方法是有效的。而且，故障信号局域时频熵值之和（为6.3769）占总时频熵值的98.77％，几乎没有能量的损失。

基于此，通过对比局域时频熵值的大小，可以判断轴承的工作状态。

表1

Claims (3)

1.一种基于LMD和局域时频熵的旋转机械故障特征提取方法，其特征在于：其步骤包括：

步骤一：利用加速度传感器对旋转机械设备进行测量，获得振动加速度信号；

步骤二：对振动加速度信号进行LMD分解，得到若干PF分量，求出各分量的瞬时幅值和瞬时频率；

步骤三：做出时频谱图，划分时频平面，计算局域时频熵；

步骤四：用局域时频熵值作为特征量，结合实验来提取故障特征。

2.根据权利要求1所述的一种基于LMD和局域时频熵的旋转机械故障特征提取方法，其特征在于：在所述步骤二中，对振动加速度信号进行LMD分解的过程步骤如下：

1)求局部均值函数m₁₁(t)，找出加速度信号x(t)所有的局部极值点n_i，求出所有相邻的局部极值点的平均值，然后将所有相邻的平均值点用直线连接起来，并用滑动平均法进行平滑处理，得到m₁₁(t)；

2)求包络估计函数a₁₁(t)，计算相邻局部极值点的包络估计值，将所有相邻两个包络估计值a_i用直线连接，然后采用滑动平均方法进行平滑处理，得到a₁₁(t)；

3)将局部均值函数m₁₁(t)从原始信号x(t)中分离出来，得到滤掉低频信号的h₁₁(t)；

4)用h₁₁(t)除以包络估计函数a₁₁(t)以对h₁₁(t)进行解调，得到s₁₁(t)；

5)对s₁₁(t)重复步骤1)~4)，直至满足1-Δ≤a_1n(t)≤1+Δ（变动量Δ<1），则迭代过程终止，否则需要继续重复步骤1)~4)；

6)把步骤1)~4)迭代过程中产生的所有包络估计函数相乘便可以得到包络信号a₁(t)；

7)将包络信号a₁(t)和纯调频信号s_1n(t)相乘便可以得到原始信号的第一个PF分量PF₁(t)；

8)将第一个PF分量PF₁(t)从原始信号x(t)中分离出来，得到一个新的信号u₁(t)，将u₁(t)作为原始数据x(t)重复步骤1)~7)，循环k次，直到u_k为一个单调函数为止，则原始信号x(t)分解为k个PF分量和一个单调函数u_k之和。

3.根据权利要求1所述的一种基于LMD和局域时频熵的旋转机械故障特征提取方法，其特征在于：在所述步骤三中，将时频平面等分为K×L个面积相等的时频块，按照不同频段将整个时频平面划分为若干个面积相等的时频段，将每一个时频段的时频熵定义为局域时频熵，第h时频段的局域时频熵定义为：

$s_h\left(q\right)=-\underset{i=M}{\overset{N}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{L}{\mathrm{\&Sigma;}}}{q}_{i,j}\ln q_{i,j}$

其中，i＝1,…,K，j=1,…,L，q_i，j为第i×j小块平面的归一化能量。

Images