经验模式分解方法

1 周期平稳信号处理 经验模式分解 2 Research Center of Condition Monitoring & Fault Diagnosis ( RCCMFD ), Southeast University，025-83792220 周期平稳信号处理 周期平稳信号分析 •现实世界中许多物理现象由于其本身具有周期性，因而所产 生的随机数据的概率模型具有周期变化的时变参数，如机器齿 轮、皮带、链条、轴、轴承及活塞的运动，地球公转和自转造 成的季节变化周期性变化等。 •数学上把均值和自相关函数呈现周期性或接近为周期性的信 号称为周期性平稳随机信号。 •齿轮箱、滚动轴承等的故障信号一般为非平稳信号，传统的 基于平稳随机信号的分析方法就无能为力。而时频分析方法， 如小波分析、Wigner-Ville分布等虽然有助于这些非平稳信号 的分析，但由于难以定量等原因，在诊断方面的应用还存在许 多问题。 3 Research Center of Condition Monitoring & Fault Diagnosis ( RCCMFD ), Southeast University，025-83792220 周期平稳信号处理 •但齿轮、滚动轴承等具有周期运动的特点，当某个齿或滚动 体有故障时，其故障信号也呈现周期性，即只有当故障齿或滚 动体受力时才有故障信号出现。 •其振动信号的统计参数具有周期变化的时变特性，其信号具 有周期性平稳的特征，即旋转机械出现故障后，其均值、自相 关函数、甚至高阶矩函数都呈现出周期性。 •这类故障信号的信噪比较低，影响故障诊断的准确性和故障 定位的精确性。 4 Research Center of Condition Monitoring & Fault Diagnosis ( RCCMFD ), Southeast University，025-83792220 周期平稳信号处理 •周期性平稳随机信号的研究可以追溯到20世纪50年代，但 是，由于缺少有效的分析工具，特别是数学处理手段，所得到 的研究成果寥寥无几。 •70年代初研究集中在周期性平稳随机信号的分解上，将周期 性平稳信号分解为一组平稳随机信号的表达式，这些平稳随机 信号之间是相关的，之后研究的重点又放在信号的自相关函数 和功率谱上，由于这两者是周期的或几乎周期的，因而可以将 其用傅里叶级数来表示。 •随后Willian A.Garder又提出了谱相关理论和谱冗余概念，并 应用于各种实际问题的解决中，在通讯、雷达、模式识别、水 文以及气象等领域取得了一些重要的研究成果。 •近期被引入机械故障诊断 5 Research Center of Condition Monitoring & Fault Diagnosis ( RCCMFD ), Southeast University，025-83792220 周期平稳信号处理 •当平稳随机信号具有遍历性时，计算随机信号的数字特征可 以用时间平均代替统计平均。用样本平均的方法进行均值、相 关函数、功率谱分析等时频域处理。 •周期平稳信号具有非平稳性质，不能使用时间平均，唯有信 号统计量变化的周期性是可以利用的重要信息。 •以均值的估计为例，虽然无法直接使用时间平均来估计信号 的均值，但如果已知周期平稳信号的周期T0，就可以对过程以 T0为周期进行采样，这样的采样值在周期平稳的前提下满足遍 历性，仍然可以用样本平均来估计其均值（一阶矩函数）： å -= ¥® + + == N NnN xx nTtxN tMt )( 12 1lim)( )( 01,m 式中的t 时间t +T0以代替，其均值保持不变 6 Research Center of Condition Monitoring & Fault Diagnosis ( RCCMFD ), Southeast University，025-83792220 周期平稳信号处理 •对均值作周期函数的傅里叶级数展开 )( 1,1, tj m xx eMtM aaå ¥ -¥= = 令2πm / T0 =α，傅里叶系数可以写为 ò- -= 2 2 1, 0 1, 0 0 )(1 T T tj xx dtetMT M aa )(1lim )( )12( 1lim )( )12( 1lim)( 2 2 2 2 0 2 2 0 0 1, 00 00 0 0 dtetx T dtetx TN dtenTtx TN tM T T tj T tjTNT TNTN tj N Nn T TNx ò ò å ò - - ¥® -+ --¥® - -= -¥® = + = + + = a a aa 时变均值的α频率分量 称为循环均值，α为一阶循环频率)(1, tM x a 7 Research Center of Condition Monitoring & Fault Diagnosis ( RCCMFD ), Southeast University，025-83792220 周期平稳信号处理 )()(lim),( 0 * 0 tt -++= å -= ¥® nTtxnTtxtR N NnN x )()(),( 0 2 åå ¥ -¥= ¥ -¥= == m tj x m mt T j xx eReRtR aa p a ttt 傅里叶级数 ),(1)( 2/ 2/ 0 0 0 dtetR T R tj T T xx aa tt - -ò= 傅里叶系数（循环自相关函数） )()(1lim)( 2/ 2/ dtetxtx T R T T tj Tx ò- -* ¥® -= aa tt 相关函数 )2/()2/(1lim)( 2/ 2/ dtetxtx T R T T tj Tx ò- -* ¥® -+= aa ttt对称形式 8 Research Center of Condition Monitoring & Fault Diagnosis ( RCCMFD ), Southeast University，025-83792220 周期平稳信号处理 ≠0的频率α称为二阶循环频率)(taxR 一个循环平稳信号的循环频率α可能有多个，其中，零循环频 率对应信号的平稳部分，只有非零的循环频率才刻画信号的 循环平稳特性。 三阶循环矩函数可用同样的方法定义得 dtetxtxtx T M tj T TTx aa tttt - -¥® ++= ò )()()( 1lim),( 2 *2/ 2/ 1213, ≠0的频率称为三阶循环频率),( 213, ttaxM 9 Research Center of Condition Monitoring & Fault Diagnosis ( RCCMFD ), Southeast University，025-83792220 周期平稳信号处理 例 一实验轴承为单列向心球轴承，外圈固定，其内圈有一点 蚀，转轴频率为30Hz，故障频率为181 Hz，采样频率为 13605 Hz，采集到的为振动加速度信号，采样点数为2048个。 其时域波形、幅值谱图、循环自相关函数谱图分别如图所示。 0 0 . 0 2 0 .0 4 0 . 0 6 0 .0 8 0 . 1 0 .1 2 0 . 1 4 0 .1 6 -2 0 0 0 -1 5 0 0 -1 0 0 0 -5 0 0 0 5 0 0 1 0 0 0 A m pl itu de u m 2 / s t s e c o n d 实验数据内圈故障时域波形 振动信号呈现出周 期脉冲现象（内圈 故障），脉冲强弱 呈现不规则的变化 10 0 1 0 0 0 2 0 0 0 3 0 0 0 4 0 0 0 5 0 0 0 6 0 0 0 7 0 0 0 - 1 - 0 . 5 0 0 . 5 1 1 . 5 x 1 0 4 A m pl itu de 0 .1 *u m 2 / s F r e q u e n c y H z 0 2 0 0 4 0 0 6 0 0 8 0 0 1 0 0 0 1 0 0 0 2 0 0 0 3 0 0 0 4 0 0 0 5 0 0 0 6 0 0 0 7 0 0 0 8 0 0 0 9 0 0 0 A m pl itu de 0 .1 *u m 2 / s F re q u e n c y H z 1 8 1 H z 3 0 H z 3 6 2 H z 5 4 3 H z 实验数据内圈故障幅值谱图 实验数据内圈故障循环自相关谱图 在5740 Hz处，幅值比 较突出，此频率为轴 承系统受到冲击以后 的自振频率，故障频 率辨别起来非常困难 100 Hz频率以下，30 Hz的频率很突出，à 转轴的频率 181 Hz及其倍频处， 幅值也很突出à内圈 的故障频率 从其边频也可以得到 转轴的频率 11 The empirical mode decomposition （EMD） for nonlinear and non- stationary time series analysis 经验模式分解方法 12 Research Center of Condition Monitoring & Fault Diagnosis ( RCCMFD ), Southeast University Main question How can we analyze non-stationary and nonlinear time series? v Stationary time series: ensemble means are time-shift invariant §Weakly stationary: true for 1st and 2nd order means ï ï þ ïï ý ü = ¥< = -stst t t rr xE constxE , 2 )( )( • A non-stationary example is a transient signal like a delta function 经验模式分解 13 Research Center of Condition Monitoring & Fault Diagnosis ( RCCMFD ), Southeast University vNonlinear time series: generated by an underlying dynamical system obeying nonlinear equations • Often linear approximations can be made -> complex time series can be decomposed into a superposition of simple solutions (e.g. sinusoidal waves) What do we want to gain from the analysis? • Determine characteristic time / frequency scales for the energy 经验模式分解 14 Research Center of Condition Monitoring & Fault Diagnosis ( RCCMFD ), Southeast University Fault diagnosis is filled with complex time series • Rub, oil oscillation, gear box vibration, etc. Motivation 经验模式分解 15 Research Center of Condition Monitoring & Fault Diagnosis ( RCCMFD ), Southeast University Existing methods Dominant approach: decomposition of the time series into component basis functions satisfying a) Completeness of the basis b) Orthogonality of the basis Examples: •Fourier methods •Wavelet analysis •Principal component analysis •etc. . . . 经验模式分解 16 Research Center of Condition Monitoring & Fault Diagnosis ( RCCMFD ), Southeast University Fourier methods Method •Decompose time series, x(t), into global sinusoidal components of fixed amplitude, aj ,(e.g. FFT) -> complete & orthogonal basis Interpretation •The spectral amplitudes, aj, yield the energy contributed by a sinusoid at frequency j that spans the whole time series 经验模式分解 17 Research Center of Condition Monitoring & Fault Diagnosis ( RCCMFD ), Southeast University Most useful when . . . . •The underlying process is linear so that the superposition of sinusoidal solutions makes physical sense •The time series is stationary, since the aj are constant ü If not, the spectral energy spreads, often requiring carefully-phased, global (possibly harmonic) sinusoidal components to reconstruct the non-uniform time series üThe Fourier spectrum looses track of time location for events -> not a local description 经验模式分解 18 Research Center of Condition Monitoring & Fault Diagnosis ( RCCMFD ), Southeast University If piece-wise stationarity can be assumed •Construct spectrogram by sliding a window across the timeseries and performing Fourier analysis 经验模式分解 19 Research Center of Condition Monitoring & Fault Diagnosis ( RCCMFD ), Southeast University But even if weak stationarity holds, Fourier methods may not give an efficient (adaptive) representation •Single decaying sinusoid is compactly described as: tetx t p322 01.0 cos)( -= 经验模式分解 20 Research Center of Condition Monitoring & Fault Diagnosis ( RCCMFD ), Southeast University Wavelet analysis Method •Decompose time series, x(t), into local, time-dilated and time-translated wavelet components, y -> complete (not necessarily orthogonal) basis dt a bttx a baW )()(1),( * -= ò ¥ ¥- y Interpretation •W represents the energy in x of temporal scale a at t=b 经验模式分解 21 Research Center of Condition Monitoring & Fault Diagnosis ( RCCMFD ), Southeast University Attractive because. . . . •Local, although higher frequencies are more localized •Uniform temporal resolution for all frequency scales But, resolution is limited by the basic wavelet •Useful for characterizing gradual frequency changes But. . . . •Non-adaptive -- same basic wavelet is used for all data 经验模式分解 22 Research Center of Condition Monitoring & Fault Diagnosis ( RCCMFD ), Southeast University Principal component analysis Method •Decompose time series, z(t), into eigenbases of the covariance matrix -> complete, orthogonal basis å =×=×= n kkkjkkjkk ffCffxftatxz 1 ,),()(),( ld Interpretation •Modes fj often interpreted as “directions” of independent variations, but they may not be physical modes Attractive because. . . . •Adaptive -- derived from the data 经验模式分解 23 Research Center of Condition Monitoring & Fault Diagnosis ( RCCMFD ), Southeast University But. . . . •Distribution of eigenvalues do not yield characteristic time or frequency scales •Eigenmodes themselves need not be linear or stationary (and therefore are still not easily analyzed by spectral methods) 经验模式分解 24 Research Center of Condition Monitoring & Fault Diagnosis ( RCCMFD ), Southeast University Instantaneous vs. global frequency Global frequency = average frequency derived by weighting by the Fourier power spectrum, |S(ù)|2 ò= wwww dS 2)( Instantaneous frequency of a signal, X(t) •Form an analytic signal, Z(t), from X(t) Take the Hilbert transform of X(t) ò ¥ ¥- =+= - = )('' ' )()()()(,)(1)( tietatiytxtzdt tt txPty q p 经验模式分解 25 Research Center of Condition Monitoring & Fault Diagnosis ( RCCMFD ), Southeast University Define instantaneous frequency as the rate of phase change dt td )(q w = üPhysical interpretation: y(t) is the best local fit (since the 1/t emphasizes local properties) of a trig function to x(t) üPolar coordinate description for z(t) allows the phase q(t) to be defined •z(t) has the same positive frequency spectrum as x(t), but zero negative frequency spectrum 经验模式分解 26 Research Center of Condition Monitoring & Fault Diagnosis ( RCCMFD ), Southeast University Goals of the new method 1) Decompose time series into superposition of components with well defined instantaneous frequency －> Intrinsic Mode Functions (IMF) s Components should (approximately) obey earlier requirements of completeness, orthogonality, locality and adaptiveness. s To get an IMF, need to • LOCALLY eliminate riding waves • LOCALLY eliminate asymmetries (defined by envelope of extrema) 经验模式分解 27 Research Center of Condition Monitoring & Fault Diagnosis ( RCCMFD ), Southeast University 2）Construct the Hilbert spectrum of each IMF, representing it in the amplitude – instantaneous frequency - time plane Properties of IMF’s An intrinsic mode function (IMF) is a function that satisfies two conditions: (1)in the whole data set, the number of extrema and the number of zero crossings must either equal or differ at most by one; and (2) at any point, the mean value of the envelope defined by the local maxima and the envelope defined by the local minima is zero. 经验模式分解 28 Research Center of Condition Monitoring & Fault Diagnosis ( RCCMFD ), Southeast University 1) corresponds loosely to finding “narrow-band” signals, or eliminating “riding-waves” 2) ensures that the instantaneous frequency will not have fluctuations arising from an asymmetric wave forms The name `intrinsic mode function' is adopted because it represents the oscillation mode imbedded in the data. With this definition, the IMF in each cycle, defined by the zero crossings, involves only one mode of oscillation, no complex riding waves are allowed. With this definition, an IMF is not restricted to a narrow band signal, and it can be both amplitude and frequency modulated. In fact, it can be non-stationary. 经验模式分解 29 Research Center of Condition Monitoring & Fault Diagnosis ( RCCMFD ), Southeast University Finding IMFs: Emperical Mode Decomposition Most data are not naturally IMF’s Unfortunately, most of the data are not IMFs. At any given time, the data may involve more than one oscillatory mode; that is why the simple Hilbert transform cannot provide the full description of the frequency content for the general data as reported by Long et al. (1995). What are we looking for? The essence of the method is to identify the intrinsic oscillatory modes by their characteristic time scales in the data empirically, and then decompose the data accordingly. 经验模式分解 30 Research Center of Condition Monitoring & Fault Diagnosis ( RCCMFD ), Southeast University Assumptions 1) At least two extrema -- one max, one min 2) Characteristic time scale defined by the time between extrema 3) If no extrema, differentiation will reveal extrema (integrate at end) 经验模式分解 31 Research Center of Condition Monitoring & Fault Diagnosis ( RCCMFD ), Southeast University The sifting process Procedure •Identify the extrema (both maxima and minima) of the data, x(t) •Generate the envelope by connecting maxima points with a cubic spline, and minima points with a cubic spline •Determine the LOCAL mean, m1, by averaging the envelope •Since IMF should have zero local mean, subtract out the mean from the data 11)( hmtx =- 经验模式分解 32 Research Center of Condition Monitoring & Fault Diagnosis ( RCCMFD ), Southeast University •h1 is probably not an IMF; repeat as necessary until it is -> End up with an IMF As described above, the process is indeed like sifting: to separate the finest local mode from the data first based only on the characteristic time scale. The sifting process, however, has two effects: (a) to eliminate riding waves; and (b) to smooth uneven amplitudes. 经验模式分解 33 Figure 3. Illustration of the sifting processes: (a) the original data; (b) the data in thin solid line, with the upper and lower envelopes in dot-dashed lines and the mean in thick solid line; (c) the difference between the data and m1. This is still not an IMF, for there are negative local maxima and positive minima suggesting riding waves. 34 Figure 4. Illustration of the effects of repeated sifting process: (a) after one more sifting of the result in figure 3c, the result is still asymmetric and still not a IMF; (b) after three siftings, the result is much improved, but more sifting needed to eliminate the asymmetry. The final IMF is shown in figure 2 after nine siftings. 35 Research Center of Condition Monitoring & Fault Diagnosis ( RCCMFD ), Southeast University Finding all the IMFs Procedure •Once the first set of “siftings” results in an IMF, define This 1st component contains the finest temporal scale in the signal •Generate the residue, r1, of the data by subtracting out c1 •The residue now contains information about longer periods -> resift to find additional components •Form the superposition of the components to reconstruct the data Linear superposition of modes Efficient/adaptive (based on data) representation khc 11 = 11)( rctx =- nnn rcrrcr =-=- -1221 ,, K 经验模式分解 36 Example of a turbulence data set Decomposition Total of 9 components •Efficient representation (relative to Fourier) of a turbulence data set •Never have oscillations of the same scale at the same time in different modes 37 38 Research Center of Condition Monitoring & Fault Diagnosis ( RCCMFD ), Southeast University Figure. The resulting empirical mode decomposition components from the wind data: (a) the original data and the components c1-c4; (b) the components c5-c9. Notice the last component, c9, is not an IMF; it is the trend. This example illustrates the advantage of adopting the successive extrema for defining the time scale; it also illustrates the difficulties of dealing with non-stationary data: even a meaningful mean is impossible to define, but for EMD this difficulty is eliminated. Figures 6a; b summarize all the IMF obtained from this repeated sifting processes. We have a total of nine components. Comparing this with the traditional Fourier expansion, one can immediately see the e_ciency of the EMD: the expansion of a turbulence data set with only nine terms. From the result, one can see a general separation of the data into locally non- overlapping time scale components. In some components, such as c1 and c3, the signals are intermittent, then the neighbouring components might contain oscillations of the same scale, but signals of the same time scale would never occur at the same locations in two different IMF components. 经验模式分解 39 Research Center of Condition Monitoring & Fault Diagnosis ( RCCMFD ), Southeast University •artificially generated multichromatic wave frequency range from 0.714Hz - 1.66Hz •calculated by second order irregular wave theory for multichromatic waves 经验模式分解 40 Research Center of Condition Monitoring & Fault Diagnosis ( RCCMFD ), Southeast University 经验模式分解 41 Research Center of Condition Monitoring & Fault Diagnosis ( RCCMFD ), Southeast University 图e-5是从轴承座上实测的外圈具有故障的6311型滚动球轴承 的振动信号时域波形，实验时轴的转频为25Hz，采样频率为 4096Hz，采样点数为1024，经计算滚动轴承外围故障特征频 率为76Hz，为了突出故障信息特征，用EMD对滚动轴承振动 信号进行分解，得到第一个满足IMF条什的分量c1，如图e-6 所示。 经验模式分解 42 Research Center of Condition Monitoring & Fault Diagnosis ( RCCMFD ), Southeast University 从原始信号的第一个IMF分量c1的波形可以看出明显的冲击特征，在一个 采样周朋内大约有19个脉冲，因此可以求得每两个脉冲之间的时间间隔为 T=0.25/19=0. 0132s，其频率为f=1/T=76Hz，正好是滚动轴承的外圈故障 特征频率，与实际情况相符。而用小波分析方法，对相同的原始信号进行 5层小波分解后得到32个结点，选择第5层的第11个结点重构得到的时域波 形如图e-7所示。可以看出有明显的脉冲，但是在一个采样周期内大约有 16个脉冲，经计算脉冲产生的频率为f=16/0.25=64Hz，与滚动轴承的外圈 故障特征频率有一点差别。 经验模式分解 43 Research Center of Condition Monitoring & Fault Diagnosis ( RCCMFD ), Southeast University 在齿轮故障诊断中的应用 由于齿轮故障的振动信息在传递过程中所经环节较多，这样使测得的振 动信号含有大量的噪声，而且高频振动信号在传递过程中大多丢失，增加 了齿轮故障诊断的难度．图e-8是一测得的断齿齿轮振动加速度信号的时 域波形，齿轮的齿数为37，转频为7Hz，采样频率为1024Hz，从图中可以 看出调幅信号的一些特征，但是要判断齿轮的缺陷还需作进一步的分析。 采用EMD对齿轮振动信号进行分解，得到第一个满足IMF条件的分量c1如 图e-8所示。 经验模式分解 44 Research Center of Condition Monitoring & Fault Diagnosis ( RCCMFD ), Southeast University 从断齿齿轮振动信号的第一个IMF分量c1的波形可以看出明 显的冲击特征，在一个采样周期内大约有7个脉冲，因此脉 冲产生的频率为7Hz，正好与齿轮的转频相等。 经验模式分解