CN106289774A - 一种滚动轴承故障识别与趋势预测方法 - Google Patents

Abstract

本发明公开了一种滚动轴承故障识别与趋势预测方法，属于旋转机械故障诊断领域。该方法主要包括5个步骤：获取轴承在正常及故障状态下的振动信号，作为原始信号；运用EMD将原始信号分解为有限个IMF分量；根据相关分析选择典型的IMF分量，并将这些分量进行求和，得到重组信号；运用FK算法对重组信号进行处理，自动获取用于包络分析的最佳中心频率和带宽，提取故障特征频率，从而实现故障类型识别；选择与原始信号相关系数最大的IMF分量的能量百分比作为故障程度评估指标。该方法简单有效，可大大提高信噪比，有助于准确识别故障类型，故障程度评估指标选取合理，可有效反映故障发展趋势，对轴承状态监测和故障评估有很大的实用价值。

Description

一种滚动轴承故障识别与趋势预测方法

技术领域

本发明属于旋转机械故障诊断领域，具体涉及一种基于典型IMF分量选择和快速峭度图的滚动轴承故障诊断方法。

背景技术

滚动轴承是旋转机械设备的关键零部件，同时也是设备最常见的故障来源之一，其运行状态直接影响到整台机器的安全性和可靠性。因此，实时掌握轴承运行的工作状态，了解故障的形成和发展是机械故障诊断领域最重要的课题之一。目前，滚动轴承状态监测和诊断最常用的方法是基于振动信号的分析处理，但是由于早期故障特征很微弱，信噪比低，有用的特征信号往往被淹没在背景噪声中。因此，消除噪声干扰，提高信噪比是轴承故障诊断的关键。

高阶谱(HOS)技术对高斯有色噪声具有很好的抑制能力，同时又能保持非线性系统的相位信息，因而对轴承早期故障具有较好的诊断能力。快速峭度图(FK)算法属于HOS技术的一种，它可自动为包络谱分析提供最佳中心频率和带宽。FK算法借鉴了二小波分解算法，先将原始信号经过FIR滤波器分解为若干个频段，然后计算各个频段的峭度值，自动选择峭度值最大的频段所对应的中心频率和带宽作为最佳中心频率和带宽，最后进行包络分析，从而实现滚动轴承故障识别。

经验模式分解(EMD)可以将复杂信号分解为有限个本征模态函数(IMF)，IMF分量反映了原始信号的局部特征信息。利用互相关系数筛选包含主要故障信息的IMF分量，剔除无关的或冗余的分量，实现降噪、抑制正常信息、突出故障信息的效果，从而使诊断结果更准确。

发明内容

针对故障轴承振动信号的非平稳及非高斯特性，本发明公开了一种基于典型IMF分量选择和快速峭度图的滚动轴承早期微弱故障诊断方法。

本发明包括以下步骤：

步骤一、获取滚动轴承振动加速度信号，作为原始信号；

步骤二、运用EMD将原始信号分解为有限个IMF分量；

具体进行如下处理：

步骤2.1：确定原始信号x(t)所有的局部极值点，然后用三次样条分别连接所有的局部极大值

点和局部极小值点形成上包络线、下包络线；

步骤2.2：上包络线、下包络线的均值记为m₁(t)，x(t)和m₁(t)的差值记为h₁(t)，即

x(t)-m₁(t)＝h₁(t)

理想地，如果h₁(t)是一个IMF,那么其就是x(t)的第一个分量。

步骤2.3：如果h₁(t)不是一个IMF，把h₁(t)当作新的原始信号，重复步骤2.1～步骤2.2k次，直到h_1k(t)满足IMF的两个条件，同时有如下公式

h_1(k-1)(t)-m_1k(t)＝h_1k(t)

设h_1k(t)为第一个IMF分量c₁(t)，即

c₁(t)＝h_1k(t)

步骤2.4：x(t)与c₁(t)的差值记为第一个残差，即

r₁(t)＝x(t)-c₁(t)

步骤2.5：把c₁(t)当作原始信号，重复步骤2.1～步骤2.4n次，直到满足终止条件：当r_n(t)是一个单调函数或常数，然后获得其它IMF分量，即

$\left.\begin{array}{c}{r}_1\left(t\right)-c_2\left(t\right)=r_2\left(t\right)\\ .\\ .\\ .\\ r_{n-1}\left(t\right)-c_n\left(t\right)=r_n\left(t\right)\end{array}\right\}$

上式中r_n(t)是原始信号的中心趋势。

步骤2.6：最终，原始信号x(t)分解为n个经验模态c₁,…,c_n和一个残差r_n(t)，即

$x\left(t\right)=\underset{i=1}{\overset{n}{\Σ}}{c}_i\left(t\right)+r_n\left(t\right)$

上式中i取值1至n。

步骤三、根据相关分析选择典型的IMF分量，并将这些分量进行求和，得到重组信号；

具体进行如下处理：

步骤3.1：理论上，常用的信号s(t)由一系列真实分量c_i(t)组成，即

$s\left(t\right)=\underset{i=1}{\overset{p}{\Σ}}{c}_i\left(t\right)$

实际上，由于存在端点效应或过分解的缺陷，EMD分解存在伪分量，设常用的信号s(t)经EMD分解后得到的真分量和伪分量分别为和x_j(t)，即

$s\left(t\right)=\underset{i=1}{\overset{p}{\Σ}}\stackrel{\‾}{c_i}\left(t\right)+\underset{j=1}{\overset{q}{\Σ}}{x}_j\left(t\right)$

上式中p,q分别表示真分量和伪分量的个数。

步骤3.2：s(t)和c_i(t)之间的互相关系数的计算公式如下

$\begin{array}{c}{R}_{s,\stackrel{\‾}{c_i}}\left(\mathrm{\τ}\right)=E\[s\left(t\right)*\stackrel{\‾}{c_i}(t+\mathrm{\τ})\]=E\[\left(\underset{i=1}{\overset{p}{\Σ}}{c}_i\right(t\left)\right)*\stackrel{\‾}{c_i}(t+\mathrm{\τ})\]\\ =R_{c_i,\stackrel{\‾}{c_i}}\left(\mathrm{\τ}\right)+{\Σ}_{j=1,j\≠i}^n{R}_{c_j,\stackrel{\‾}{c_i}}\left(\mathrm{\τ}\right)\≈R_{c_i,\stackrel{\‾}{c_i}}\left(\mathrm{\τ}\right)\≈R_{c_i}\left(\mathrm{\τ}\right)\end{array}$

上式中“*”表示卷积符号，E[]表示期望运算符，τ为互相关系数的变量函数。

由于IMF分量的正交性，所以有

$\underset{j=1,j\≠i}{\overset{n}{\Σ}}{R}_{c_j,\stackrel{\‾}{c_i}}\left(\mathrm{\τ}\right)\≈0$

很显然，伪分量x_k(t)和原始信号s(t)之间的相关系数应该为零，即

$R_{s,x_k}\left(\mathrm{\τ}\right)=E\[s\left(t\right)*x_k(t+\mathrm{\τ})\]={\mathrm{\Σ}}_{j=1}^n{R}_{c_j,x_k}\left(\mathrm{\τ}\right)\≈0.$

步骤3.3：结论，原始信号与真分量之间的相关系数约等于真分量的自相关系数，而原始信号与伪分量之间的相关系数非常小，接近于零。

步骤四、运用FK算法对重组信号进行处理，提取故障特征频率，从而实现故障类型识别；具体进行如下处理：

步骤4.1：在非平稳情况下，Y(t)为激励信号x(t)的系统响应，按照Wold-Cramer进行频域分解，得到

$Y\left(t\right)={\∫}_{-\mathrm{\∞}}^{+\mathrm{\∞}}H(t,f)e^{j2\mathrm{\π}ft}dX\left(f\right)$

式中H(t,f)表示系统的时变传递函数，X(f)为x(t)的傅里叶变换。

步骤4.2：Y(t)的四阶谱累积量C_4Y(f)定义如下，

$C_{4Y}\left(f\right)=S_{4Y}\left(f\right)-2S_{2Y}^2\left(f\right),f\≠0$

式中S(f)表示瞬时矩，用来表示复包络的能量强度，其计算公式如下，

S_2nY(f)＝E{|H(t,f)dX(f)|²ⁿ}/df

步骤4.3：SK的计算公式如下，

$K_Y\left(f\right)=\frac{C_{4Y}\left(f\right)}{S_{2Y}^2\left(f\right)}=\frac{S_{4Y}\left(f\right)}{S_{2Y}^2\left(f\right)}-2,f\≠0$

实际中，振动信号包含有用信号f(t)和噪声n(t)，即

x(t)＝f(t)+n(t)

因此，SK最终计算公式如下，

$K_x\left(f\right)=\frac{K_f\left(f\right)}{{\[1+\mathrm{\ρ}\left(f\right)\]}^2}$

式中，ρ(f)＝S_n(f)/S_f(f)表示噪声与信号的比值，S_n(f)和S_f(f)分别为n(t)和f(t)的功率谱。

步骤4.4：根据SK的快速算法，得到快速峭度图，选择最佳中心频率和带宽，然后根据包络解调提取故障特征频率，从而实现故障类型识别。

步骤五、选择与原始信号相关系数最大的IMF分量的能量百分比作为故障程度评估指标。不同状态下测得的振动信号的IMF能量百分比不同，IMF能量百分比可以用来反映轴承的不同运行状态。实验结果表明，当轴承存在故障时，IMF1与原始信号的相关度最大，能量百分比最高，可以作为故障程度评估指标。

本发明的优点与积极效果在于：

(1)双重降噪，信噪比大大提高

EMD可自适应地将复杂非平稳的振动信号分解为有限个平稳的IMF分量，这些分量的频带从高到低变化，并且反映原始信号的局部信号特征。根据相关分析，选择与原始信号相关度大的分量，即包含主要故障信息的分量，剔除相关度低的分量和伪分量，然后将选择的分量进行重组，获得重组信号，可有效提高信噪比。FK算法自动为包络分析提供最佳中心频率和带宽，克服了传统包络分析靠主观经验确定中心频率和带宽的缺陷。同时根据FK设计的滤波器处理后的包络分析更能直观清晰地反映故障特征频率，从而使诊断结果更准确。

(2)故障程度评估

本发明不仅可以对故障类型进行识别，还可以对故障程度进行评估。不同故障信号分解得到的IMF分量的能量百分比不同，选择与原始信号相关度最大的IMF分量(最能体现故障信息的分量)的能量百分比作为故障程度评估指标。结果表明，当轴承存在故障时，IMF1与原始信号相关度最大且所占能量百分比最高，可有效地反映出故障发展趋势。

附图说明

图1是本发明的滚动轴承故障诊断方法整体步骤流程图；

图2是本发明实例中case_1的包络解调图；

图3是本发明实例中case_2的包络解调图；

图4是本发明实例中case_3的包络解调图；

图5是Hilbert包络解调谱图；

图6是本发明实例中故障程度评估图。

具体实施方式

本发明是一种轴承故障诊断算法，包括故障类型识别和故障程度评估两部分，下面将结合附图，对本发明作进一步的详细说明。

以美国凯斯西储大学轴承振动数据库的部分实验数据中SKF的6205-2RS深沟球轴承为例。

本发明公开的滚动轴承故障诊断方法的整体步骤流程图如图1所示，具体步骤如下：

S1、选取振动信号

选取轴承在正常，不同程度(0.1778毫米,0.3556毫米,0.5332毫米)的单点内圈故障，外圈故障，滚动体故障状态下的振动信号作为原始信号x(t)。

S2、利用EMD将原始信号分解

经验模式分解是基于这样的假设：任何信号均由不同的简单固有模态组成，即任何信号均可被分解为若干个IMF之和，每个IMF必须同时满足以下两个条件：a.在整个数据集上，极点数和过零点数相等或至多相差一个；b.在任一时间点上，信号的局部极大值所确定的上包络线和局部极小值确定的下包络线的局部均值为零。具体算法如下：

(1)确定原始信号x(t)所有的局部极值点，然后用三次样条分别连接所有的局部极大值点和局部极小值点形成上包络线、下包络线；

(2)上包络线、下包络线的均值记为m₁(t)，x(t)和m₁(t)的差值记为h₁(t)，即

x(t)-m₁(t)＝h₁(t)

理想地，如果h₁(t)满足IMF的两个条件,那么其就是x(t)的第一个分量。

(3)如果h₁(t)不是一个IMF，把h₁(t)当作新的原始信号，重复上述步骤(1)～(2)k次，直到h_1k(t)满足IMF的两个条件，同时有如下公式

h_1(k-1)(t)-m_1k(t)＝h_1k(t)

设h_1k(t)为第一个IMF分量c₁(t)，即

c₁(t)＝h_1k(t)

(4)x(t)与c₁(t)的差值为第一个残差，即

r₁(t)＝x(t)-c₁(t)

(5)把c₁(t)当作原始信号，重复(1)～(4)n次，获得其他IMF分量，即

$\left.\begin{array}{c}{r}_1\left(t\right)-c_2\left(t\right)=r_2\left(t\right)\\ .\\ .\\ .\\ r_{n-1}\left(t\right)-c_n\left(t\right)=r_n\left(t\right)\end{array}\right\}$

(6)最终，原始信号x(t)分解为n个经验模态c₁,…,c_n和一个残差r_n(t)，即

$x\left(t\right)=\underset{i=1}{\overset{n}{\Σ}}{c}_i\left(t\right)+r_n\left(t\right)$

S3、基于相关分析的IMF选择

由于存在边界效应、过分解缺陷，经EMD分解得到的IMF分量存在伪分量，这些伪分量会对诊断结果造成干扰，所以必须剔除。采用相关分析的方法，筛选包含主要故障信息的IMF分量，剔除无关的或冗余的分量，实现降噪、抑制正常信息、突出故障信息的效果，从而使诊断结果更准确。

具体原理如下：

(1)理论上，信号s(t)由一系列的真实分量c_i(t)组成，即

$s\left(t\right)=\underset{i=1}{\overset{p}{\Σ}}{c}_i\left(t\right)$

(2)实际上，EMD分解存在伪分量，设信号s(t)经EMD分解后得到的真分量和伪分量分别为和x_j(t)，即

$s\left(t\right)=\underset{i=1}{\overset{p}{\Σ}}\stackrel{\‾}{c_i}\left(t\right)+\underset{j=1}{\overset{q}{\Σ}}{x}_j\left(t\right)$

(3)s(t)和c_i(t)之间的互相关系数计算公式如下

$\begin{array}{c}{R}_{s,\stackrel{\‾}{c_i}}\left(\mathrm{\τ}\right)=E\[s\left(t\right)*\stackrel{\‾}{c_i}(t+\mathrm{\τ})\]=E\[\left(\underset{i=1}{\overset{p}{\Σ}}{c}_i\right(t\left)\right)*\stackrel{\‾}{c_i}(t+\mathrm{\τ})\]\\ =R_{c_i,\stackrel{\‾}{c_i}}\left(\mathrm{\τ}\right)+{\Σ}_{j=1,j\≠i}^n{R}_{c_j,\stackrel{\‾}{c_i}}\left(\mathrm{\τ}\right)\≈R_{c_i,\stackrel{\‾}{c_i}}\left(\mathrm{\τ}\right)\≈R_{c_i}\left(\mathrm{\τ}\right)\end{array}$

由于IMF分量的正交性，所以，有

$\underset{j=1,j\≠i}{\overset{n}{\Σ}}{R}_{c_j,\stackrel{\‾}{c_i}}\left(\mathrm{\τ}\right)\≈0$

很显然，伪分量和原始信号之间的相关系数应该为零，即

$R_{s,x_k}\left(\mathrm{\τ}\right)=E\[s\left(t\right)*x_k(t+\mathrm{\τ})\]={\mathrm{\Σ}}_{j=1}^n{R}_{c_j,x_k}\left(\mathrm{\τ}\right)\≈0$

总结，原始信号与真分量之间的相关系数约等于真分量的自相关系数，而原始信号与伪分量之间的相关系数非常小，接近零。

本发明以四种转速(1797,1772,1750,1730rpm)情况下的正常信号(case_0:Normal)和三种不同程度(case_1:0.1778mm,case_2:0.3556mm,case_3:0.5332mm)的内圈故障信号为例。将信号分解16个IMF分量，然后计算各分量与原始信号之间的相关系数，表1给出了前8个分量与原始信号之间的相关系数。众所周知，当相关系数低于0.3时，则称为低相关，可忽略。根据表1，case_1的前三个分量(0.2703≈0.3)，case_2的前两个分量以及case_3的前一个分量对应的相关系数均大于或约等于0.3。所以，针对case_1、case_2、case_3这三种情形，选择的IMF分量分别为前三个、前两个、前一个，然后将对应情形的分量求和，得到重组信号。

表1原始信号与IMF分量之间的互相关系数

cases	IMF1	IMF2	IMF3	IMF4	IMF5	IMF6	IMF7	IMF8
									case_1	0.8552	0.3750	0.2703	0.1794	0.0366	0.0173	0.0035	0.0000
case_2	0.8384	0.4777	0.2291	0.1202	0.0720	0.0398	0.0143	0.0000
									case_3	0.9681	0.1536	0.1062	0.1436	0.0591	0.0182	0.0032	0.0000

1、运用FK算法对重组信号进行处理，提取故障特征频率，从而实现故障类型识别

FK算法是基于这样的假设：每类冲击存在一组最佳中心频率和带宽使得其峭度值最大。

FK是谱峭度(SK)的快速算法。SK算法具体步骤如下：

(1)在非平稳情况下，Y(t)为激励信号x(t)的系统响应，按照Wold-Cramer进行频域分解，

$Y\left(t\right)={\∫}_{-\mathrm{\∞}}^{+\mathrm{\∞}}H(t,f)e^{j2\mathrm{\π}ft}dX\left(f\right)$

式中H(t,f)表示系统的时变传递函数，X(f)为x(t)的傅里叶变换。

(2)Y(t)的四阶谱累积量C_4Y(f)定义如下，

$C_{4Y}\left(f\right)=S_{4Y}\left(f\right)-2S_{2Y}^2\left(f\right),f\≠0$

式中S(f)表示瞬时矩，用来表示复包络的能量强度，其表示如下，

$S_{2\mathrm{nY}}\left(f\right)=E\left\{\right|H(t,f)\mathrm{dX}\left(f\right){|}^{2n}\}/\mathrm{df}$

(3)SK的计算公式如下，

$K_Y\left(f\right)=\frac{C_{4Y}\left(f\right)}{S_{2Y}^2\left(f\right)}=\frac{S_{4Y}\left(f\right)}{S_{2Y}^2\left(f\right)}-2,f\≠0$

实际中，振动信号包含有用信号f(t)和噪声n(t)两部分，即

x(t)＝f(t)+n(t)

因此，SK最终计算公式如下，

$K_x\left(f\right)=\frac{K_f\left(f\right)}{{\[1+\mathrm{\ρ}\left(f\right)\]}^2}$

式中，ρ(f)＝S_n(f)/S_f(f)表示噪声与信号的比值，S_n(f)和S_f(f)分别为n(t)和f(t)的功率谱。

总结：根据K_x(f)可知，当噪声很强，即ρ(f)很大，K_x(f)约为零；反之，K_x(f)≈K_f(f)。所以，SK可有效地反映冲击信号的强弱，而故障信号多为冲击信号。

快速峭度图算法采用二小波算法，首先将信号分解为一系列的频带，然后计算这些频带的峭度值；峭度值最大对应的频带的中心频率和带宽即为最佳参数，然后根据最佳参数进行滤波、包络解调，得到故障特征频率，从而实现故障类型识别。

本发明中轴承在不同转速情况下的内圈故障特征频率如表2所示，运用FK处理重组信号和正常信号，结果如表3，图2，图3和图4所示。其中K_max表示最大峭度值，B_w表示最佳带宽，f_c表示最佳中心频率。根据图2(c)、图3(c)和图4(c)，其对应的内圈故障特征频率(154.9Hz,155.9Hz,155.6Hz)及其倍频均可被清楚地提取出来，即该方法可准确无误地识别出内圈故障类型。

为了体现本发明在故障类型识别方面的优越性，还采用了Hilbert包络解调方法作为对比。采用Hilbert包络解调方法处理1730rpm下的case_2内圈故障信号，结果如图5所示。根据图5，虽然故障特征频率(155.9Hz)可以被提取出来，但是包络谱中还存在一些幅值非常明显且无法解释的频率成分，如，69.5Hz和242.2Hz，这些成分会严重干扰诊断结果，甚至可能导致错误的结果。相比之下，本发明图图2(c)、图3(c)和图4(c)中提取出的故障特征频率相对其他频率成分来说非常明显。

表2不同转速下的轴承内圈故障频率

转速/rpm	1730	1750	1772	1797
					频率/Hz	156.1	157.9	159.8	162.1

表3不同内圈故障重组信号的FK参数值

cases	case_0	case_1	case_2	case_3
					Kmax	0.9	3.6	17.5	9.3
Bw(Hz)	62.5	750	500	6000
					fc(Hz)	4968.75	4125	3750	3000

S4、选择故障程度评估指标

实验结果表明，IMF能量百分比会随着轴承状态的改变而发生变化，其表达式如下

$P_i=E_i/{\mathrm{\Σ}}_{i=1}^n{E}_i$

式中E_i，P_i分别表示第i个IMF分量的能量及其百分比。

实验结果表明，当轴承存在故障时，IMF1与原始信号的相关度最大，所占能量百分比最高，所以本发明选择IMF1的能量百分比作为故障程度评估指标。

本发明计算了四种转速条件下的四种不同故障情形的IMF1能量百分比，结果如表4和图6所示。表4和图6表明：随着故障程度的越来越严重，IMF1的能量百分比越来越大，且正常情形与故障情形之间的变化趋势相对明显，不同程度的故障之间变化趋势相对平缓。不同转速下的总体变化趋势一致，说明该指标具有较好的鲁棒性。

表4不同情形的IMF1能量百分比

cases	1730(rpm)	1750(rpm)	1772(rpm)	1797(rpm)
					case_0	0.2675	0.1585	0.1580	0.2938
case_1	0.8071	0.7949	0.7993	0.7907
					case_2	0.8431	0.9210	0.9348	0.8763
case_3	0.9570	0.9656	0.9602	0.9380

通过以上对滚动轴承故障类型识别及故障程度评估方法的详细描述，可见本发明的基于典型IMF分量选择和快速峭度图的滚动轴承故障诊断方法具有明显的优势：

1、双重降噪，信噪比大大提高

选择与原始信号相关度大的分量，即包含主要故障信息的分量，剔除相关度低的分量和伪分量，缩小了信号处理范围，减少了计算量，同时可有效提高信噪比。FK算法自动为包络分析提供最佳中心频率和带宽，克服了传统包络分析靠主观经验确定中心频率和带宽的缺陷。同时根据FK设计的滤波器处理后的包络分析更能直观清晰地反映故障特征频率，从而使诊断结果更准确。

2、故障程度评估

本发明不仅可以对故障类型进行识别，还可以对故障程度进行评估。选择与原始信号相关度最大的IMF分量(最能体现故障信息的分量)的能量百分比作为故障程度评估指标。结果表明，当轴承存在故障时，IMF1与原始信号相关度最大且所占能量百分比最高，可有效地反映出故障发展趋势。

Claims (4)

1.一种滚动轴承故障识别与趋势预测方法，其特征在于：该方法为一种基于典型IMF分量选择和快速峭度图的滚动轴承早期微弱故障诊断方法；

该方法包括以下步骤：

步骤一、获取滚动轴承振动加速度信号，作为原始信号；

步骤二、运用EMD将原始信号分解为有限个IMF分量；

步骤三、根据相关分析选择典型的IMF分量，并将这些分量进行求和，得到重组信号；

步骤四、运用FK算法对重组信号进行处理，提取故障特征频率，从而实现故障类型识别；

步骤五、选择与原始信号相关系数最大的IMF分量的能量百分比作为故障程度评估指标；不同状态下测得的振动信号的IMF能量百分比不同，IMF能量百分比可以用来反映轴承的不同运行状态。

2.根据权利要求1所述的一种滚动轴承故障识别与趋势预测方法，其特征在于：

步骤二、运用EMD将原始信号分解为有限个IMF分量；具体进行如下处理，

步骤2.1：确定原始信号x(t)所有的局部极值点，然后用三次样条分别连接所有的局部极大值点和局部极小值点形成上包络线、下包络线；

步骤2.2：上包络线、下包络线的均值记为m₁(t)，x(t)和m₁(t)的差值记为h₁(t)，即

x(t)-m₁(t)＝h₁(t)

理想地，如果h₁(t)是一个IMF,那么其就是x(t)的第一个分量；

步骤2.3：如果h₁(t)不是一个IMF，把h₁(t)当作新的原始信号，重复步骤2.1～步骤2.2k次，直到h_1k(t)满足IMF的两个条件，同时有如下公式

h_1(k-1)(t)-m_1k(t)＝h_1k(t)

设h_1k(t)为第一个IMF分量c₁(t)，即

c₁(t)＝h_1k(t)

步骤2.4：x(t)与c₁(t)的差值记为第一个残差，即

r₁(t)＝x(t)-c₁(t)

步骤2.5：把c₁(t)当作原始信号，重复步骤2.1～步骤2.4n次，直到满足终止条件：当r_n(t)是一个单调函数或常数，然后获得其它IMF分量，即

$\left.\begin{array}{c}{r}_1\left(t\right)-c_2\left(t\right)=r_2\left(t\right)\\ \begin{array}{c}.\\ .\\ .\end{array}\\ r_{n-1}\left(t\right)-c_n\left(t\right)=r_n\left(t\right)\end{array}\right\}$

上式中r_n(t)是原始信号的中心趋势；

步骤2.6：最终，原始信号x(t)分解为n个经验模态c₁,…,c_n和一个残差r_n(t)，即

$x\left(t\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{c}_i\left(t\right)+r_n\left(t\right)$

上式中i取值1至n。

3.根据权利要求1所述的一种滚动轴承故障识别与趋势预测方法，其特征在于：

步骤三、根据相关分析选择典型的IMF分量，并将这些分量进行求和，得到重组信号；

具体进行如下处理：

步骤3.1：理论上，常用的信号s(t)由一系列真实分量c_i(t)组成，即

$s\left(t\right)=\underset{i=1}{\overset{p}{\Σ}}{c}_i\left(t\right)$

实际上，由于存在端点效应或过分解的缺陷，EMD分解存在伪分量，设常用的信号s(t)经EMD分解后得到的真分量和伪分量分别为和x_j(t)，即

$s\left(t\right)=\underset{i=1}{\overset{p}{\Σ}}\stackrel{\‾}{c_i}\left(t\right)+\underset{j=1}{\overset{q}{\Σ}}{x}_j\left(t\right)$

上式中p,q分别表示真分量和伪分量的个数；

步骤3.2：s(t)和c_i(t)之间的互相关系数的计算公式如下

$\begin{array}{c}{R}_{s,\stackrel{\‾}{c_i}}\left(\mathrm{\τ}\right)=E\[s\left(t\right)*\stackrel{\‾}{c_i}\left(t+\mathrm{\τ}\right)\]=E\[\left(\underset{i=1}{\overset{p}{\mathrm{\Σ}}}{c}_i\left(t\right)\right)*\stackrel{\‾}{c_i}\left(t+\mathrm{\τ}\right)\]\\ =R_{c_i,\stackrel{\‾}{c_i}}\left(\mathrm{\τ}\right)+{\mathrm{\Σ}}_{j=1,j\≠i}^n{R}_{c_j,\stackrel{\‾}{c_i}}\left(\mathrm{\τ}\right)\≈R_{c_i,\stackrel{\‾}{c_i}}\left(\mathrm{\τ}\right)\≈R_{c_i}\left(\mathrm{\τ}\right)\end{array}$

上式中“*”表示卷积符号，E[]表示期望运算符，τ为互相关系数的变量函数；

由于IMF分量的正交性，所以有

$\underset{j=1,j\≠i}{\overset{n}{\Σ}}{R}_{c_j,\stackrel{\‾}{c_i}}\left(\mathrm{\τ}\right)\≈0$

很显然，伪分量x_k(t)和原始信号s(t)之间的相关系数应该为零，即

$R_{s,x_k}\left(\mathrm{\τ}\right)=E\[s\left(t\right)*x_k(t+\mathrm{\τ})\]={\mathrm{\Σ}}_{j=1}^n{R}_{c_j,x_k}\left(\mathrm{\τ}\right)\≈0;$

步骤3.3：结论，原始信号与真分量之间的相关系数约等于真分量的自相关系数，而原始信号与伪分量之间的相关系数非常小，接近于零。

4.根据权利要求1所述的一种滚动轴承故障识别与趋势预测方法，其特征在于：

步骤四、运用FK算法对重组信号进行处理，提取故障特征频率，从而实现故障类型识别；

具体进行如下处理：

步骤4.1：在非平稳情况下，Y(t)为激励信号x(t)的系统响应，按照Wold-Cramer进行频域分解，得到

$Y\left(t\right)={\∫}_{-\mathrm{\∞}}^{+\mathrm{\∞}}H(t,f)e^{j2\mathrm{\π}ft}dX\left(f\right)$

式中H(t,f)表示系统的时变传递函数，X(f)为x(t)的傅里叶变换；

步骤4.2：Y(t)的四阶谱累积量C_4Y(f)定义如下，

$C_{4Y}\left(f\right)=S_{4Y}\left(f\right)-2S_{2Y}^2\left(f\right),f\≠0$

式中S(f)表示瞬时矩，用来表示复包络的能量强度，其计算公式如下，

S_2nY(f)＝E{|H(t,f)dX(f)|²ⁿ}/df

步骤4.3：SK的计算公式如下，

$K_Y\left(f\right)=\frac{C_{4Y}\left(f\right)}{S_{2Y}^2\left(f\right)}=\frac{S_{4Y}\left(f\right)}{S_{2Y}^2\left(f\right)}-2,f\≠0$

实际中，振动信号包含有用信号f(t)和噪声n(t)，即

x(t)＝f(t)+n(t)

因此，SK最终计算公式如下，

$K_x\left(f\right)=\frac{K_f\left(f\right)}{{\[1+\mathrm{\ρ}\left(f\right)\]}^2}$

式中，ρ(f)＝S_n(f)/S_f(f)表示噪声与信号的比值，S_n(f)和S_f(f)分别为n(t)和f(t)的功率谱；

步骤4.4：根据SK的快速算法，得到快速峭度图，选择最佳中心频率和带宽，然后根据包络解调提取故障特征频率，从而实现故障类型识别；

步骤五、选择与原始信号相关系数最大的IMF分量的能量百分比作为故障程度评估指标；

不同状态下测得的振动信号的IMF能量百分比不同，IMF能量百分比可以用来反映轴承的不同运行状态。