基于EMD包络功率谱滚动轴承故障诊断

基于EMD的包络功率谱滚动轴承故障诊断
摘要
本文简要介绍了滚动轴承的故障特征及判定方法，并基于经验模态分解EMD(Empirical Mode Decomposition)与功率谱相结合的方法。

将原始信号分解成不同尺度的固有模态函数IMF(Intrinsic Mode Function),求取IMF分量的包络,计算其包络功率谱,得到轴承的故障特征频率，进行诊断。

关键词：滚动轴承；故障诊断；EMD；包络功率谱
Abstract
This paper briefly introduces the method of judging and fault features of rolling bearing, and based on empirical mode decomposition EMD (Empirical Mode Decomposition) method combined with power spectrum. The original signal is decomposed into different scales of intrinsic mode function IMF (Intrinsic Mode Function), the envelop take IMF components, calculate the envelope power spectrum, get the fault characteristic frequency of bearing, diagnosis.
Keywords: rolling bearing; fault diagnosis; EMD; envelope power spectrum
1、概述
滚动轴承是大部分旋转机械的重要组成部分，具有效率高、摩擦阻力小、装配方便、润滑易实现等优点而应用广泛。

滚动轴承同时又是是旋转机械中较为薄弱的部分，在正常工作条件下工作一段时间总会出现各种类型的失效。

由于发现不及时引起设备停机获设备损坏，造成了生产商的巨大损失。

因此，轴承故障诊断分析显得十分必要。

滚动轴承是由内圈、外圈、滚动体和保持架四种零部件组成。

滚动轴承故障诊断分析，是根据滚动轴承运行时产生的不同的信号包含对机械状态识别与诊断非常有用的各种信息。

有效地分析、处理这些信息，建立它们和设备运行状态之间的联系，是设备故障诊断的基础。

滚动轴承在运行过程中由于各种原因产生十分复杂的振动信号。

为了确定滚动轴承的故障，通常要从振动信号的幅域、时域和频域等多角度去分析，得出准确结论需要进行多方面的验证。

本文从频域分析角度作功率谱，进而分析滚动轴承的工作情况。

相对判定标准：相对判定标准是指对轴承的同一部位定期进行振动检测，通过功率谱分析，以轴承无故障情况下的功率谱为基准，根据谱峰频率来进行判断。

2、EMD及功率谱基本原理
（1）EMD
利用EMD方法对非线性、非平稳信号的分解获得一系列表征信号特征间尺度的固有模态函数IMF(Intrinsic Mode Function)，使得各个IMF是单分量的幅值或频率调制信号.IMF要满足2个条件: ①整个数据序列的极值点与过零点的个数相等或最多相差一个；②在任意时刻，由局部极大值点形成的上包络与由局部极小值点形成的下包络的均值为零。

这2个条件实际上使得分解得到的IMF是窄带信号。

同时,EMD分解方法还建立在以下假设上: ①信号至少有2个极点,一个最大值和一个最小值；②特征时间尺度通过2个极值点之间的时间定义；③若数据缺乏极值点但有形变点，则可通过微分数据一次或几次获得极值点，然后再通过积分来获得分解结果。

对任意一个实信号x(t)进行EMD的具体步骤是:
1) 确定x(t)上的所有极大值点和极小值点，然后，将所有极大值点和所有
极小值点分别用三次样条曲线连接起来，将这两条曲线分别作为x(t) 的上下包络线。

计算出它们的平均值曲线m1 (t) ，用x(t)减去m1(t)得
h1(t)=x(t)-m1(t) （2.1.1）
如果h1(t)不满足IMF的条件，需要把h1(t)作为原信号重复上面的步骤得到 h11(t)=h1(t)-m11(t) （2.1.2）
筛选k次直到h1k(t)变为一个IMF，即
h1k(t)=h1(k-1)(t)-m1k(t) (2.1.3) 这样就从原信号中分解出了第一个IMF，称为第一阶IMF，记作
c1(t)=h1k(t) (2.1.4)
2) 从原信号中减去c1(t)得第一阶剩余信号r1(t)
r1(t) = x(t) - c1(t) (2.1.5) 由于第一阶剩余信号r1(t)还包含着更长周期的分量，因此，把r1 ( t)作为新的原信号，重复步骤1，对后面的也进行同样的筛选，这样依次分解得到r2(t)=r1(t)-c2(t)
r3(t)=r2(t)-c3(t)
…
rn(t)=rn-1(t)-cn(t) （2.1.6）
直至剩余信号rn(t)中的信息对所研究内容意义很小，或者变成一个单调函数不能再筛选出基本模式分量为止。

至此，信号x ( t)已被分解成n个基本模式分量ci ( t)和一个剩余信号rn(t)。

这样,由式(2.1.5)和式(2.1.6)得到:
（2.1.7）
进一步,各个IMF分量可通过Hilbert变换进行包络解调。

但是,由于所分析信号的有限长度、信号的两端点不能确定是极点，那么，在进行三次样条插值的时候，必然使得信号的上下包络在信号的两端附近严重扭曲，即产生端点效应，可用极值点对称延拓法来处理该问题。

（2）双谱分析
双谱属于高阶谱,与功率谱相比具有如下特点: ①功率谱是实数,不包含相位信息，而高阶谱是复数，因而保留了相位信息；②能抑制噪声的影响；③保
留了系统的非线性信息。

因此，用高阶谱分析振动信号更容易获得特征信息。

本文采用双谱估计的直接法进行计算, 即将观测数据分段，利用FFT计算数据段的离散傅里叶变换，进而估计各阶频域矩，利用累积量谱与矩谱之间的关系
求得双谱估计。

当轴承发生故障时, 采样信号的特征为受干扰的冲击调制信号,即
(10)
式中,ωi 为调制源(包括轴承故障特征频率及其谐波频率)；
ω0 为载波频率；
b为任意常数。

因此故障轴承振动信号解调后的信号包含了故障特征频率的一簇谐波,且相位是互相关联的,即存在二次相位耦合现象。

若设ωF 为轴承的故障特征频率,则双谱的(ωF ,ωF )处必然出现相位耦合现象,从而双谱在(ωF ,ωF )处会有明显的谱峰。

3、诊断结果
根据以上分析, 可得双谱估计为
当ω =ωF时,必然出现明显的峰值,将峰值对应的频率与理论计算的轴承的故障特征频率相比较, 就可以得出正确的结论,同时减小计算量,增强频谱图的可视性。

本文中使用Matlab编程实现,首先对原始振动信号进行零均值化处理,并实施EMD分解,对分解后IMF分量进行包络解调,最后利用双谱分析提取轴承的故障特征。

图1、2、3即为故障轴承利用上述方法得到的。

图 1 正常/外圈故障/内圈故障的信号图
图2 正常/外圈故障/内圈故障IMF图
图3 正常/外圈故障/内圈故障的包络图
4、结论
通过对比正常轴承的功率谱，外圈故障和内圈故障的功率谱有着明显不同，表明存在着异常损伤，能够判定轴承故障。

同时，外圈故障和内圈故障功率谱之间也存在着明显的差异，与所给数据故障类型吻合。

由此可知，将EMD与双谱相结合,提取轴承的故障特征频率，能准确地区别正常及故障滚动轴承，该方法适用于滚动轴承的诊断。

参考文献
[1]张键.机械故障诊断技术.机械工业出版社机械工业出版社.2008
[2]李成超.基于经验模态分解的轴承故障诊断方法.大连大学
[3]任玥.基于Hilbert_Huang变换的滚动轴承智能诊断方法研究.西南交通大学.2007
附：
编程程序
正常轴承（内圈故障/外圈故障）：
%采样频率
fs=12000;
load 97.mat;%正常，内圈故障209.mat，外圈故障130.mat
xdata=X097_DE_time(1:1024);
xdata=(xdata-mean(xdata))/std(xdata,1);
%时域波形
figure(1); N=1024;
plot(1:N,xdata,'k-');
xlabel('时间 t/n');
ylabel('电压 V/v'); %db10小波进行4层分解 %一维小波分解
[c,l]=wavedec(xdata,4,'db10'); %重构第1~4层细节信号
d4=wrcoef('d',c,l,'db10',4);
d3=wrcoef('d',c,l,'db10',3);
d2=wrcoef('d',c,l,'db10',2);
d1=wrcoef('d',c,l,'db10',1); %显示细节信号
figure(2);
subplot(4,1,1);
plot(d4,'k-','LineWidth',2); ylabel('d4'); subplot(4,1,2); plot(d4,'k-','LineWidth',2); ylabel('d3'); subplot(4,1,3); plot(d4,'k-','LineWidth',2); ylabel('d2'); subplot(4,1,4); plot(d4,'k-','LineWidth',2); ylabel('d1');
xlabel('时间 t/s'); %第1层细节信号的包络谱
y=hilbert(d1);
ydata=abs(y);
y=y-mean(y);
nfft=10240;
p=abs(fft(ydata,nfft));
figure(3);
plot((0:nfft/2-l)/nfft*fs,p(l:nfft/2),'k-');
xlabel('频率 f/Hz'); ylabel('功率谱 P/W');。