Rolling element bearings localized fault diagnosis using signal differencing and median filtration

2.1. Signal differencing (differentiation)

The detection of defective bearing components -as discussed in the introduction (see Fig. 1) is based on detecting impacts (impulses) which results from the passage of a rolling element over a spalled region. Impacts are characterized by a sudden jump in acceleration and generally have high frequency content; as they tend to excite structural and sensor natural frequencies. To measure the rate of change in the acceleration response (Jerk forces) and to make these forces more prominent, the vibration signal is differentiated (signal differencing). This process will in fact provide a measure of the force change in the system (jerk) and thus has physics-based meaning.

For a digitized signal, time differencing (differentiation) simply refers to the process of taking the difference between two consecutive samples. Signal differencing has been recognized in the field of detecting hidden information. For instance, it has been utilized recently in the field of detecting hidden information with computer forensic analysis (steganalysis), e.g. [19, 20].

In ref [21], it is shown that for a time series of length $N$: $x_1,x_2,x_1,x_2,\dots ,x_N$, First differencing ($\nabla$) as can be seen from Eq. (1), leads to a series of length $N-1$ such that:

An impulsive noise or a short duration pulse tends to introduce an uncharacteristic discontinuity in a correlated signal [22]. This discontinuity becomes more detectable when the signal is differentiated and in this sense, the differentiation operation is equivalent to a decorrelation or spectral whitening process [22]. For a sinusoidal signal, as in Eq. (2):

where: $A$ – signal amplitude, $\omega$ – angular velocity, $\varphi$ – phase shift.

First differencing introduces a phase shift of $\pi /2$ and changes the amplitude by $\omega$ as can be seen from the Eq. (3):

Signal differencing can then be seen as a “high-pass filter” for the time series; as it passes and enhances the high frequency variations and attenuates low-frequencies.

In ref [23], first differencing has been recognized as the residual (pre-whitened signal) of an Autoregressive (AR) process of an order 1 (AR (1). Sawalhi and Randall [24] showed that signal differencing and in particular the 4th derivative was found to give the highest kurtosis and clarity of detecting impulses. Bozchalooi and Liang [25] gave a very detailed analysis and description to the enhancement of the signal to interference ratio (SIR) that can be obtained by the successive differentiations of the vibration signal in the case of detecting faults in rolling element bearings. Bozchalooi and Liang [25] proposed a differentiation method to enhance the fault delectability. It is reported that the iterative application of a differentiation step can enhance the relative strength of the impulsive faulty bearing signal component with respect to the vibration interferences. This has the effect of preserving the effectiveness of amplitude demodulation and hence leads to more accurate fault detection.

As a means to select the optimum differentiation order, kurtosis Eq. (4) between one differentiation process and the other can be utilized to guide the selection process; as the ultimate aim is to enhance the impulsiveness of the signal and to enhance the ability of the detection of impacts:

where:$N$ – number of samples of signal $x$, ${\mu }_x$ – mean of signal $x.$

The use of maximum correlated Kurtosis (MCK) [7] can be thought of as a better alternative to kurtosis; as it specifically enhances the impulsiveness of the fault related features. MCK of first shift can be calculated using Eq. (5), which is based on assuming a zero-mean, i.e. ${\mu }_x=$ 0. MCK requires the knowledge of the period of impact ($T$), which depends on the location of the spall. In this sense and as the algorithm is used for surveillance without prior use of the type of fault, kurtosis would be a better option. Correlated Kurtosis of first – shift:

2.2. Smoothening using a median filter

As the importance of detecting the faulty component in a bearing signal resides in analyzing the periodicity/frequency of the impulses, rather than the impulse high frequency content, it is necessary to low pass filter the squared enveloped signal (energy of the extracted fault impulses, which was found by the differencing method). If only periodic impulses exist in the squared envelope signal, then low pass filtering would be appropriate to use. However, due to the likelihood of having extraneous impulses with large magnitudes in the squared enveloped signal, a median filter is proposed here for a robust smoothening.

Median filters [26, 27] have been utilized in a number of fields such as image processing [28], speech analysis [29] and seismic signal processing [30]. A Median filter has been recognized as an effective processing tool for suppressing impulse noise, while keeping the edge information [30]. In its simplest utilization, a median filter takes $N$ samples of the input signal (usually $N$ is odd: 3, 5, 7, etc.) and reports the middle value of a sorted version of the $N$ samples. This process is a nonlinear operation which makes it impossible to get a system function (Hz) for a median filter in a signal processing sense [31].

Generally speaking, for a signal with sampling frequency $f_s$ (Hz), a median filter of length $N$ will attenuate frequency content above $f_s/N$. As the envelope spectrum is often chosen to include at least the first five to ten harmonics of the defect frequency, the upper $N$ bound of the median filter selection can be guided by the required envelope spectrum frequency limit. For instance, to include the first 10 harmonics, the length of the Median filter $N$ can be calculated using Eq. (6) as:

where: $N$: Recommended filter length (samples), $BPF$: Defect characteristic frequency defined as: BPFI, BPFO or BSF (see Fig. 1).

As impulses’ spacing is closer in the case of the inner race defect, it is more intuitive to base Eq. (6) on the inner race defect characteristic frequency (BPFI) and thus to limit the number of samples of the median filter to one tenth of the spacing between the impact impulses. This means that the analyst can select the median filter length as one of the following: (3, 5, 7, ..., one tenth of the spacing between the impulses).

Median filter can be easily implemented using the Matlab® function “medfilt1” [32]

2.3. Surveillance and diagnosis algorithm code

The following Matlab® code shows the simplicity of the proposed code. Only two parameters are needed to run the code ($D$: Differentiation order and $N$: Median filter order). Only one line of code is used in the algorithm to find the squared envelope signal. Its envelope is found using a normal fft analysis.

% Input

% $x$: time domain signal

% fs: sampling frequency

% Output

% env: Squared Enveloped signal

$D=$ 4; % Differentiation order based Kurtosis (see % section 2.1)

$N=$ 9; % Median filter order (see section 2.2, Eq. (6))

env = (medfilt1 (abs (diff ($x$, $D$)).^2, $N$)); % Only line of code to find the squared % envelope

% Envelope Spectrum

NFFT = length (env);

df = NFFT/fs;

$f=$ linspace (0, fs, NFFT);

XF = fft (env, NFFT);

plot (f (1: end/2), abs (XF (1: end/2))/NFFT);

title ('Envelope Spectrum of the median filtered signal'); xlabel ('Frequency (Hz)').

5.1. Single stage gearbox

Vibration signals collected from a single stage gearbox (spur gear, 32 teeth,1:1 ratio) with a defective (seeded fault) inner race and a defective outer race bearing are used to demonstrate the effectiveness of the presented surveillance algorithm. The gearbox under consideration was run at 10 Hz and with a torque of 50 Nm. An accelerometer was placed on the top of a defective bearing and vibration signals were collected at a sampling rate of 48 kHz. The gear mesh frequency (32X10) is 320 Hz. The estimated ball pass frequency of the inner race (BPFI) is 71.2 Hz, while that of the outer race (BPFO) is 48.9 Hz. Zoom plots representing one shaft rotation of measured vibration signals, as well as their frequency content (resolution of 5.8 Hz/line) are presented in Fig. 11. The time domain signature shows clear gearmesh events with impacts (bearing defect) of high frequency nature. Gear vibrations are noticed mainly in the frequency region up to 8 kHz in the form of discrete components, while vibrations from the bearing fault are noticed to excite resonances in the high frequency region. The interaction between gears’ vibration and the localized defective bearing vibration is additive [34] in the sense that they excite different frequency bands and can be separated. This additive nature and the excitation of high resonances by the bearing spall are the main motives behind the envelope analysis (high frequency resonance technique).

Fig. 12 shows the kurtosis values for the first one hundred orders of differentiation for both the inner race and the outer race faults. The behavior is different between the two cases, with the best results achieved around the third and fourth derivatives. In the case of the outer race fault, it is noticed that, unlike the inner race case, the kurtosis tends to decrease with the increase of differentiation order; as this seems to highlight some high frequency noise from the gears. This can be seen clearly by inspecting the differentiated signals (zoomed in around three impacts) shown in Fig. 13. The fault impacts can be clearly seen using the third/forth derivatives and no further advantage was observed when using higher order derivatives.

Fig. 11Top: raw time domain signal (0.1 seconds, which correspond to 1 shaft rotations). Bottom: power spectrum density (resolution 5.8 Hz/line). Left column: inner race, right: outer race

Fig. 12Kurtosis vs. differentiation order for the single stage gearbox: a) inner race, b) outer race

Fig. 13Left: inner race: a) raw signal, b) third differentiation, c) fiftieth differentiation. Right: outer race: d) raw signal, e) forth differentiation, f) fiftieth differentiation

In Figs. 14(a) and 14(b), the squared values of the third and fourth derivatives in the case of the inner race and outer race faults were plotted in Figs. 13(c) and 13(d) respectively. A median filter of 65 samples (approximately (48000/71×10) was used to filter signals Figs. 14(a) and 14(b) with the results shown in Figs. 14(c) and 14(d). The 65 samples filter represents roughly one tenth the spacing between two consecutive impacts. The low pass effect can be clearly seen through the removal of high frequency components and obtaining the squared envelope. Finally, the spectrum of Figs. 14(c) and 14(d) (squared envelope spectrum) was obtained and is presented in Figs. 15(a) and 15(b) respectively.

The spectra of the squared values of the differenced-median filtered signals obtained in Figs. 14(c) and 14(d) for the inner and outer races faults are plotted in Figs. 15(a) and 15(b) respectively. A clear identification of the source of the impacts is achieved through the presence of the BPFI and the BPFO. The modulation effect associated with the inner race fault is not clearly observed. There has been no established clear link (through simulations) to the effect of differentiation on signal modulation. The effect can be present due to light load and as so the signal becomes lightly modulated, Outer race fault harmonics are clear, with weak sidebands at the cage frequency. The end results of diagnosis using the proposed surveillance processing are compared to the results obtained from the semi-automated algorithm presented in section 4. The processing through the different steps are shown in Figs. 16 and 17, which identify the best band for demodulation process and provides an envelope spectrum analysis on this basis.

Fig. 14a) Squared values of signal Fig. 13(b), b) squared values of signal Fig. 13(d), c) signal Fig. 14(a) after using a median filter of order 65, d) signal Fig. 14(b) after using a median filter of order 65

Fig. 15a) Spectrum of signal Fig. 14(c) (inner race fault), b) spectrum of signal Fig. 14(d) (outer race fault)

The wavelet kurtogram highlighting the excited bands of the preprocessed signals is shown in Fig. 16(a) and 16(b) for the inner and outer race respectively. Filtered signals are presented in Figs. 16(c) and 16(d). Note the correspondence between the kurtosis of the differenced signal (Fig. 12) and the kurtosis of the filtered signals after pre-processing and SK filtration.

Fig. 16a) Wavelet kurtogram of AR-MED filtered signal of the inner race fault, b) inner race AR-MED-SK filtered signal, c) wavelet kurtogram of AR-MED filtered signal of the outer race fault, d) outer race AR-MED-SK filtered signal

Finally, the squared envelop spectrum extracted using the semi-automated algorithm is presented in Fig. 17 and can be compared with Fig. 15 to illustrate the effectiveness. Comparisons indicate the success of the proposed algorithm in detecting and diagnosing the localized faults. The modulation of the inner race by the shaft speed is not prominent and aligns with the squared envelope result obtained in Fig. 15; thus, indicating that the load is not heavy enough to cause this modulation.

Fig. 17a) Squared envelope spectrum for the inner race faulty signal obtained using the semi-automated process using the filtered signal of Fig. 16(b), b) squared envelope spectrum for the outer race faulty signal obtained using the semi-automated process using the filtered signal of Fig. 16(d)

5.2. Inner race fault in a planetary bearing of a helicopter gearbox

The second example used to demonstrate the practicality of the proposed approach is an interesting case of a faulty planetary bearing (both inner race and roller faults) from a 3-planet gearbox. This case represents a more challenging and complicated gearbox system compared to the one presented in section 5.1. The signal analyzed was captured at a rate of 51200 Hz. The faulty planetary bearing had a severe inner race spalling and defective rolling elements. The BPFI was estimated at 117.5 Hz, BSF estimated at 37 Hz while the cage frequency was estimated at 9.81 Hz. The analysis of this case using the semi-automated approach was presented in reference [20].

Fig. 18 presents the time domain signal and the power spectrum density of the vibration signal captured from one of the accelerometers placed on the ring. The presences of the gearmesh frequencies as well as their harmonics/sidebands are noticed over the full frequency bandwidth up to 22 kHz.

Fig. 18a) Time domain signal (just over three revolutions), b) power spectrum density (1.25 Hz/line)

The time domain signal presented in Fig. 18 was differenced and the kurtosis values were plotted for each order. The result, as presented in Fig. 19, shows that the 9th differentiation order gives the highest kurtosis.

A zoom in, showing the first 30.000 samples (just over three rotations of the output speed), of the raw signal and its 9th differentiation order is plotted in Fig. 20. This reveals some of the hidden impacts which are attributed to the faulty bearing

Fig. 19Kurtosis vs. differentiation order for the planetary gearbox signal (maximum Kurtosis at order 9)

The effect of squaring and applying a median filter on the 9th differentiated signal of an order of 45 (one tenth the spacing between two inner race impacts) are shown in Figs. 21(a) and 21(b) respectively. Note the effect of the median filter on enhancing the impacts and making them protrude more from the background noise. The spectrum of the signal obtained in Fig. 21(b) is plotted in Fig. 22 (two zoom-in regions (a) from 0 to 100 Hz, and (b) from 90 Hz to 500 Hz). The cage frequency and the BPFI can be clearly seen (with clear four harmonics), thus provides a clear diagnosis of the fault.

The wavelet kurtogram After the removal of all shaft harmonics, signal pre-whitening and MED application is plotted in Fig. 23 [33]. This shows that the best band to enhance impulsiveness and extract fault components is present in the vicinity of the 20 kHz frequency. The result obtained from signal differencing and median filtration in Fig. 22 is very well comparable to the result of the AR-MED-SK, which is plotted in Fig. 24 [33].

Fig. 20Planetary gearbox signal: a) raw signal, b) ninth differentiation

Fig. 21a) Squared values of signal Fig. 20(b), b) signal Fig. 21(a) after using a median filter of order 45

Fig. 22a) Spectrum of signal Fig. 21(b) (0-100 Hz), b) spectrum of signal Fig. 21(b) (90-500 Hz)

Fig. 23Wavelet kurtogram of AR-MED filtered 3-planet signal [19]

Fig. 24SK automated envelope analysis: a) squared envelope spectrum (0-100 Hz), b) squared envelope spectrum (100-500 Hz) [33]