基于振动信号的齿轮故障诊断方法研究硕士毕业论文

l基于EMD的齿轮故障诊断一基于EMD（经验模式分解）的振动信号特征提取1.1EMD的研究背景齿轮摩擦是齿轮箱系统中常见的非线性振动故障。当齿轮存在局部摩擦时,摩擦信号是在动静摩擦的短时间内产生的，具有频带宽的特点。在采集振动信号时混有宽频带随机噪声的噪声信号,其频带相互叠加。在振动信号中还含有能量较大的、与转速有关的背景信号,所以利用传统的滤波方法很难仅保留含有故障信息的摩擦信号而将背景信号和噪声信号滤掉。由于用小波分析方法进行故障信号特征提取存在缺点，美籍华人NordenE.Huang等人针对瞬时频率的概念进行深入研究后,创造性地提出了本征模式函数(IntrinsicModeFunction,IMF)的概念以及将任意信号分解为本征模式函数组成的新方法——基于经验模式分解(EmpiricalModeDecomposition,EMD)方法，从而赋予了瞬时频率合理的定义、物理意义和求法，初步建立了以瞬时频率为表征信号交变的基本量，以本征模式分量为时域基本信号的新的时频分析方法体系。EMD分解算法的基本思想是：对一给定信号，先获得信号极值点，通过差值获得信号包络，获得均值，与均值的差得到分解的一层信号；如此重复，获得分解结果：r，即l个IMF和一个残值r。由于EMD方法具有自适应特性，适宜于非线性、非平稳信号的分解，该方法应用于齿轮的故障诊断分析中。结果表明，该方法能够突出齿轮振动信号的故障特征，从而提高齿轮故障诊断的准确性[19]。为此，本文将EMD方法应用到齿轮的故障诊断中，利用EMD方法对齿轮信号进行分解，得到若干个IMF，然后对包含故障信息的IMF作功率谱分析，有效提取故障信号特征信息。1.2经验模式分解(EMD)方法原理经验模式分解方法的主要目的是通过对非线性、非平稳信号()xt的分解，获得一系列表征信号特征时间尺度的本征模式函数()(0,1,2,3...)ictin和一个残余分量()nrt，各个IMF是单分量的幅值或者频率调制信号。每个IMF需要满足以下几个条件信号：l(1)整个信号中零点数与极点数相等或者相差1；(2)信号上任意一点，由局部极大值点确定的包络线和由局部极小值点确定的包络线的均值均为零，即信号关于时间轴局部对称[20]。分解过程如下：（1）设信号为()xt，取其上下包络局部均值组成的系列为1()nt，令它们的差值为1()ht11()()()htxtnt(4·1)对于非线性、非平稳数据来说，一般一次处理不会形成基本模式分量，一些非对称波依然存在。因此要把1()ht看作待处理的数据继续重复上述操作。取1()ht的均值系列为2()nt，得到它们之前的差值为2()ht212()()()hthtnt(4·2)重复k次操作，就可以得到：(1)()()()kkkhthtnt(4·3)当()xt满足基本模式分量的条件后，就会获得了第一个基本模式分量1()ct1()()kctht(4·4)（2）将基本模式分量从信号中分离出来，可以得到剩余系列1()rt11()()()rtxtct(4·5)（3）把1()rt作为新的数据按上述操作进行处理，依次类推，可得到：212()()()rtrtct323()()()rtrtct……1()()()nnnrtrtct(4·6)这个处理过程在满足预先设定的停止标准后就可停止。（4）将式(4·5)和式(4·6)相加后得到1()()()ninixtctrt(4·7)这样就把一个信号()xt分解为n个基本模式函数()(0,1,2,3...)ictin和一个残余分量()nrt[19,20]，其中分解出的n个分量ic分别包含信号从高频到低频的不同频率段成分，而剩余分量nr是原始信号的中心趋势。(4·6)由于每个固有模态函数nr为平稳信号，由式(4·7)可以定义每个固有模态分量的自相关函数为l()[()()]cciiRnEckckn(4·8)1.3基于EMD的振动信号特征提取分析当齿轮发生局部摩擦故障时，设备振动信号中将出现周期性的冲击，这些周期性的冲击信号由于噪声等干扰，所以很难通过频谱分析等传统方法加以识别，因此本文提出了一种基于EMD的冲击信号提取方法。依照IMF的定义，周期性的冲击信号无法以独立的IMF形式被分解出来，但它会以类似冲击响应信号的形式出现在高频段的IMF中。假设如下仿真信号：()sin(2*50)()()xttptnt(4·9)该仿真信号是由一个50Hz的正弦信号和与之同频率的周期性冲击信号()pt构成的，并加入了一定量的均匀白噪声()nt。采样频率为2000Hz单位为微米（m）。图4·1表示了该仿真信号的波形及频谱。00.050.10.150.20.250.30.350.40.450.5-1-0.500.51该仿真信号的波形及频谱图t/sA/μm010020030040050060070080090010000246f/HzA/μm图4·1该仿真信号的波形及频谱图l图4·2表示载入的同频率的冲击信号()pt。00.050.10.150.20.250.30.350.40.450.5-0.2-0.15-0.1-0.0500.050.10.150.2t/sA/μm冲击信号图4·250Hz的冲击信号图4·3表示一定量的白噪声()nt。00.050.10.150.20.250.30.350.40.450.500.020.040.060.080.10.120.140.160.180.2白噪声t/sA/μml图4.3一定量的白噪声对这三个信号合成，时域图x(t)如下图。00.10.20.30.40.5-1.5-1-0.500.51时域合成信号t/sA/μm图4·4时域合成图对该信号进行EMD分解，按照4.1.2节所列步骤，依次进行分解得到了3个IMF1()ct，2()ct，3()ct及残余分量3()rt。图4·5表示了分解所得的3个IMF分量，由图可见第三个IMF对应于仿真信号的正弦分量，而前两个高频IMF包含了仿真信号中的周期性冲击信号和白噪声成分。观察前两个高频IMF可见，其波形出现了周期性类似于冲击响应的信号成分，这些成分就对应了原仿真信号中的周期性冲击信号。下面是EMD分解，分解得两个IMF:C1(t),C2(t),C3(t)及残余分量R3(t)。imf=emd(x);C1=imf{1};C2=imf{2};C3=imf{3};N=length(x);c=linspace(0,(N-1)*Ts,N);l00.10.20.30.40.5-0.500.5C1(t)A/μm00.10.20.30.40.5-0.500.5C2(t)A/μm00.10.20.30.40.5-202C3(t)A/μmt/s图4·5试验信号EMD分解结果图为了提取冲击信号的频率，可以对第一个IMF进行解调出来，解调后得到的结果如图4·6所示。其中上图为第一个IMF的包络，下图为包络谱。从包络谱可得，其峰值所在的频率50Hz与仿真信号冲击出现的频率一致。该仿真实验得出结论，利用EMD可以将信号中的周期性冲击信号以冲击响应信号的形式分解出来，然后再通过包络解调技术就可以确定原信号中冲击信号出现的频率。以下是C1(t)的解调过程，包络与包络谱。B1=hilbert(C1);am=abs(B1);subplot(2,1,1);plot(c,abs(B1));ylabel('A/μm');xlabel('t/s');W=0.1*fftshift(fft(am));f=linspace(-1000,1000,length(c));l00.10.20.30.40.500.20.40.60.8第一个IMF的解调结果图A/μmt/s020040060080010000510A/μmf/Hz图4·6第一个IMF的解调结果图第2节基于EMD对齿轮故障诊断的研究载入4个原始信号的原始数据如图4·7。程序实现：loadN300001j.csvs=N300001j;loadSA300001j.csvx=SA300001j;loadU03U300001j.csvy=U03U300001j;loadU05U300001j.csvz=U300001j;fs=10000;l00.020.040.060.080.10.12-0.100.1正常信号time(s)signalvalue00.020.040.060.080.10.12-0.200.2SA信号time(s)signalvalue00.020.040.060.080.10.12-0.200.2U03信号time(s)signalvalue00.020.040.060.080.10.12-0.200.2U05信号time(s)signalvalue图4·7原始信号将原始数据进行FFT变换，得到各信号的频谱，如图4·8。程序实现：S=0.01*fftshift(fft(s));X=0.01*fftshift(fft(x));Y=0.01*fftshift(fft(y));Z=0.01*fftshift(fft(z));fs=linspace(-fs,fs,length(t1));01000200030004000500060007000800090001000000.050.1正常信号频谱f(Hz)signalvalue01000200030004000500060007000800090001000000.20.4SA信号频谱f(Hz)signalvalue01000200030004000500060007000800090001000000.10.2U03信号频谱f(Hz)signalvalue01000200030004000500060007000800090001000000.20.4U05信号频谱f(Hz)signalvalue图4·8各信号的频谱图l对各信号进行一维离散小波变换，进行两层分解。程序实现：[C,L]=wavedec(x,2,'db4');%进行一维两层分解scA2=appcoef(C,L,'db4',2);%从C中提取第二层的低频系数scD2=detcoef(C,L,2);%从C中提取第二层的高频系数scD1=detcoef(C,L,1);%从C中提取第一层的高频系数00.050.1-0.0100.01正常信号高频D1time(s)signalvalue00.020.04-0.0200.02正常信号中频D2time(s)signalvalue00.020.04-0.200.2正常信号低频A2time(s)signalvalue00.050.1-0.0100.01SA信号高频D1time(s)signalvalue00.020.04-0.0200.02SA信号中频D2time(s)signalvalue00.020.04-0.500.5SA信号低频A2time(s)signalvalue00.050.1-0.0200.02U03信号高频D1time(s)signalvalue00.020.04-0.0200.02U03信号中频D2time(s)signalvalue00.020.04-0.500.5U03信号低频A2time(s)signalvalue00.050.1-0.0100.01U05信号高频D1time(s)signalvalue00.020.04-0.0200.02U05信号中频D2time(s)signalvalue00.020.04-0.500.5U05信号低频A2time(s)signalvalue图4·9各信号的高中低频对scD1进行经验模式分解，结果如图4·10。程序实现：imf=emd(scD1);C1=imf{1};C2=imf{2};C3=imf{3};N=length(scD1);c=linspace(0,(N-1)/f,N);l00.050.1-0.0100.01正常高频D1C1(t)A/μm00.050.1-505x10-3正常高频D1C2(t)A/μm00