Energy weighting method and its application to fault diagnosis of rolling bearing

2.1. Time-frequency mapping

The energy weighting method was started from time-frequency analysis. As an adaptive non-parametric time–frequency analysis method, HHT offers fine time-frequency resolution and the result is free of cross-term interferences [16]. As the central part of HHT, the algorithm of empirical mode decomposition (EMD) is depicted in quotation [6]. Some individuals are decomposed from the original signal, nearly mono-component signals with ‘Hilbert-friendly’ waveforms, named as intrinsic mode function (IMF).

However, the approach also has its problems such as scale separation. To overcome the problem, a noise-assisted data analysis method named EEMD is adopted [17, 18]. The principle of EEMD is adding white noise to populate the whole-frequency space uniformly with the constituting components of different scales. The basic idea of EEMD is adding a white noise series to the original signal and decomposing the data with added white noise into IMFs. The process is repeated with different white noise series each time to obtain the means of corresponding IMFs of the decompositions $\mathbf{i}\mathbf{m}{\mathbf{f}}_i=\frac{1}{m}\sum _{k=1}^m\mathbf{i}\mathbf{m}{\mathbf{f}}_{i,k}$, $(i=\mathrm{1,2}...n;k=\mathrm{1,2}...m)$as the final result. The original signal can be decomposed into a collection of IMFs and a residue ${\mathbf{r}}_n$, it can be denoted as:

Utilize the Hilbert transform to each IMF to compute the instantaneous frequency $\mathbf{\omega }\left(t\right)$ and amplitude $\mathbf{a}\left(t\right)$.The signal can be expressed in the following form and the time-frequency map is acquired:

2.2. Multi-scale binary spectrum

The time-frequency map is a $M\times N$ matrix, $M$ is the sampling length of original signal, $N$ is frequency analysis length of HHT. For $f\in N$, time series $\mathbf{x}(t_i,f)$ stand for the time-varying energy for a particular frequency bin. If the signal contains impulse component, some energy time series may have energy transient at the instant of impulse appearance [19]. So, we put $\mathbf{B}(t_i,f)=1$ if $\left|\mathbf{x}(t_i,f)\right|=\mathrm{m}\mathrm{a}\mathrm{x}\left\{\left|\mathbf{x}(t_k,f)\right|:i-d\le k\le i+d\right\}$ i.e. energy value of the midpoint is the highest in the interval of $2d+1$ points. Otherwise, we put $\mathbf{B}(t_i,f)=0$. The binary spectrum $\mathbf{B}(t,f)$ is obtained by applying this rule for every energy time series.

Selection of value d has a significant effect on binary spectrum calculation. The process should be repeated for several times with different $d$ value to depict multi-scale binary spectrum that can reflect a variety of frequency characteristics. The first binary spectrum is corresponding to the characteristic frequency we focus on, so the interval of $2d+1$ points should cover more than two energy impulses but less than three energy impulses so that the second energy impulse can be recognized.

For the first binary spectrum, many fake energy impulses lead by noise will be existed. So, larger scopes of the interval are need in the second and other binary spectrums to reduce the affection of noise and reveal actual shock related points. The formula of interval parameter $d$ can be defined as follows:

where $f_s$ is sample rate, $f_f$ is fault frequency, $C$ is repeat count and $⌊x⌋$ is the largest even integer not greater than $x$. The larger repeat count can lead to better performance of characteristic enhancement. But it also means to be more operand and more cutting of effective results in the weights of energy merging. So, the repeat time $C$ can be 4-6 as suggestion.

As the result, $C$ binary spectrums ${\mathbf{B}}_1(t,f)$, ${\mathbf{B}}_2(t,f)$,…, ${\mathbf{B}}_C(t,f)$ are acquired.

2.3. Weights of energy analysis

In this step we calculate the weights of energy by summation of each binary spectrum in time domain:

Each binary spectrum has blanks at the head and the end, this phenomenon can be found in the weights of energy also. Blanks of weights of energy will lead to mistake of the following calculation of power spectrum. So, we must cut the blanks to prevent the mistakes. But the length of each weights of energy will be different if we delete the blanks directly because length of the blank equals to the interval parameter $d$ we choose and the parameter $d$ is different in each calculation of binary spectrum. In the calculation of the last binary spectrum, the repeat count $C$ is the maximum so the parameter $d$ is the maximum and can be written as $d_{max}$. The method of blank processing is cutting the $d_{max}$ points of the head and the end and apply this method to all weights of energy. After the process of blanks, each weights of energy ${\mathbf{W}}_c\left(t\right)$ has a length of $M-2d_{max}$ points.

Repeat the process of weights of energy calculation and blank processing to each binary spectrum. Then sum all the weights of energy to get a compound weight of energy $\widehat{\mathbf{W}}\left(t\right)$:

The power spectrum of compound weight of energy $\widehat{\mathbf{W}}\left(t\right)$ can be calculated as follows:

where ${\mathbf{F}}_W\left(\omega \right)$ is the Fourier transform of $\widehat{\mathbf{W}}\left(t\right)$.

The process of weights of energy analysis has been repeated for times. To balance the performance of characteristic enhancement and computation quantity, the repeat time can be 4-6. Each weights of energy can reflect the distribution of shock energy in different frequency scale by variation of parameter $d$. Analyze the compound weight of energy which synthesizes multi-scale details of frequency characteristics with power spectrum. Then the characteristic frequency can be preliminary observed.

4.1. Wheel set experiment

The wheel set test stand is shown in Fig. 9(a). The constitution includes the frame, a wheel set (a healthy bearing and a damaged bearing included), wheel set lifting device, driving device (electric motors and rubber wheels), brake, hydraulic system and control system. Two bearings are installed on the wheel set by interference fit. Lifting device with groove on the top lift the wheel set and embed it into the groove of frame cantilever. After wheel set fixation, electric motor with rubber wheel rub the side face of wheel to drive it. Rotating speed of wheel is 469 rpm. Model of bearings applied in this experiment is 352226x2-2z, its parameters are listed in Table 1. Vibration signals are acquired by NI PXI DAQ device with acceleration sensors. The sensors shown in Fig. 9(b) are installed with magnet. We take the signal acquired by the vertical fixed sensor for signal processing because it can reflect the fault more directly.

Fig. 9Scene of wheel set test

a) Wheel set test stand

b) Installation of the sensors

4.2. Bearing outer-race defect identification

Fault bearing applied in the experiment is shown in Fig. 10. A defect on outer-race can be seen in the photo. The defect is caused by abrasion in daily use, so it has an irregular shape. Characteristic frequency of the outer-race defect is 67.3 Hz. Vibration acceleration signal is acquired and presented in Fig. 11. Sampling rate of data acquisition is 5120 Hz. Sample time is 2 s.

With direct FFT analysis, FFT amplitude spectrum of the vibration signal is shown in Fig. 12. It’s difficult to find characteristic frequencies caused by the outer-race defect. For further analysis, the envelope spectrum of the vibration signal is shown in Fig. 13, only rotating frequency and resonances can be identified also. The peak at the fault frequency about 67 Hz is similar with another at 75 Hz. So, it’s unable to justify the outer-race defect.

Fig. 10Outer-race fault bearing

Fig. 11Vibration acceleration signal of outer-race fault bearing test

Fig. 12FFT amplitude spectrum of outer-race fault bearing vibration signal

Fig. 13Envelope spectrum of outer-race fault bearing vibration signal

The vibration signal is processed with EEMD. Set the noise level as 0.2 and ensemble number as 100. After the EEMD preprocessor-based HHT a time-frequency map is obtained and shown in Fig. 14.

In this situation, parameters of Eq. (3) is $f_s=$5120, $f_f=$67.3 and the progress of multi-scale binary spectrum calculation can be repeated for 6 times. So, the constant $d$ can be selected as $d_1=$38, $d_2=$76, $d_3=$114, $d_4=$152, $d_5=$190, $d_6=$228 corresponding to $C_1=$1, $C_2=$2, $C_3=$3, $C_4=$4, $C_5=$5, $C_6=$6. The results of multi-scale binary spectrum calculation are shown in Fig. 15.

Fig. 14Time-frequency map of outer-race fault bearing vibration signal

Fig. 15Binary spectrum of outer-race fault bearing vibration signal

a)${\mathbf{B}}_1\left(t,f\right)$, $d_1=$38

b)${\mathbf{B}}_2\left(t,f\right)$, $d_2=$76

c)${\mathbf{B}}_3\left(t,f\right)$, $d_3=$114

d)${\mathbf{B}}_4\left(t,f\right)$, $d_4=$152

e)${\mathbf{B}}_5\left(t,f\right)$, $d_5=$190

f)${\mathbf{B}}_6\left(t,f\right)$, $d_6=$228

Weights of energy are calculated and depicted in Fig. 16. Here the blanks of each weights of energy are treated with 228 points cutting at the head and the end. Then the weights of energy are emerged and the compound weight of energy ${\widehat{W}}_r\left(t\right)$ is shown in Fig. 17. The power spectrum of ${\widehat{W}}_r\left(t\right)$ is acquired as Fig. 18.

Fig. 16Weights of energy with blank processing

a)${\mathbf{W}}_1\left(t\right)$

b)${\mathbf{W}}_2\left(t\right)$

c)${\mathbf{W}}_3\left(t\right)$

d)${\mathbf{W}}_4\left(t\right)$

e)${\mathbf{W}}_5\left(t\right)$

f)${\mathbf{W}}_6\left(t\right)$

Fig. 17Compound weight of energy W^r(t)

Fig. 18Power spectrum of W^r(t)

For this experiment, impulses generated by bearing defect have been overwhelmed by noise of drive system. So, the fault characteristic is difficult to be extracted. According to the analysis of energy weighting method, fault frequency and the resonances have been enhanced by multi-scale binary spectrum and converted to time domain by weights of energy. The compound weight of energy synthesized results of each analysis scales. Its power spectrum can reveal the fault characteristic clearly.