Simultaneous faults identification of rolling element bearings and gears by combining kurtogram and independent component analysis

2.1. Introduction on envelope analysis

In theory, the envelope analysis can be explained as that the rolling elements of REBs or the meshing of gears strike a local fault. It introduces a short duration impact. Due to the impulsive feature, even the weak vibration from the incipient faults of REBs or gears can lead to some high-frequency resonances of local structures in the rotating machinery, so that bearings, bearing houses, and sensors etc. can be excited. In signal processing field, the high-frequency resonance responses can be considered as high-frequency carriers (with resonance frequencies, respectively) that are modulated by the weak vibrations. On running condition of rotating machinery, the vibration from gearboxes with local faults includes the quasi-periodic or periodic impulsive trains generated by faulty REBs or gears. If the impulsive trains included in the vibration can be separated by a band-pass filter in a suitable high-frequency band for isolating most low-frequency disturbances in the gearboxes, the envelopes of the impulsive trains can be obtained by one of popular demodulation algorithms; and then the characteristic frequencies or the periods related to the local faults can be calculated from the envelopes by envelope spectra analysis schemes.

2.2. Envelope extraction based on kurtogram

Several demodulation approaches can be employed for the envelope extraction, for instance, the Hilbert transform [2], the wavelet transform [3] and the kurtogram [4]. It is worth mentioning that the newly developed kurtogram method takes the advantage of determining the optimal demodulation frequency band adaptively for the band-pass filtering and the fast kurtogram algorithm [4], which can be calculated efficiently. Thus it will be utilized in our study. A brief introduction of the kurtogram is introduced as follows.

A higher spectral kurtosis (SK) value means stronger impulsive components contained in the corresponding frequency band. Then, SK can be used as an indicator to locate the frequency band that includes the most impulsive components and to determine the optimal demodulation band for the envelope extraction. There are various versions of the SK, the new one developed by Antoni [4] is widely recognized and applied in applications. The SK is defined as in [4] by:

where $H(t,f)$ is the time-frequency complex envelope of the band-pass filtered signal $x\left(t\right)$ at a center frequency $f$, and can be calculated by the fast kurtogram algorithm. It should be noted that from the viewpoint of a band-pass filter, the SK can be explained as the function of the center frequency $f$ and the bandwidth $∆f$. The determination of the optimal combination ($f_o$, $∆f_o$) leads to the so-called kurtogram. The fast kurtogram algorithm realizes the selection of the optimal demodulation frequency band by designing a tree-like multi-rate filter-bank structure with quasi-analytic band-pass filters, where the filter banks at different levels are composed by binary or 1/3-binary filter structures. The schematic of the binary filter structure is shown in Fig. 1.

Fig. 1Schematic of the binary frequency/frequency resolution plane of the fast kurtogram

The detailed instruction of the fast kurtogram algorithm can be found in [4]. The filtered sequence with the maximum SK value is determined as the envelope extraction result, and the corresponding combination ($f_o$, $∆f_o$) is the corresponding optimal filtering parameters. The step can be expressed by:

where $f_o$ and $∆f_o$ denote the optimal center frequency and the optimal bandwidth, respectively, $c_o\left(n\right)$ is the optimal complex envelope, the notion ‘$\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{a}\mathrm{x}$’ represents to obtain the parameters pertaining to the SK’s maximum value. Then, the optimal envelope $c_o\left(n\right)$ is determined.

3.1. Brief on ICA

As well known, the can be used to separate the observed mixtures into different independent components (ICs) according to the different statistical properties of sources. The simplest linear noiseless model of the ICA can be expressed as:

where $\mathbf{x}$ denotes an observed $m$-dimensional vector, $\mathbf{x}={\left(x_1,x_2,\dots ,x_m\right)}^T$, the element $x_i$ with $i=$1, …, $m$ denotes the $i$th observed mixture; $\mathbf{s}$ represents an $n$-dimensional random vector, $\mathbf{s}={\left(s_1,s_2,\dots ,s_m\right)}^T$, the element $s_j$ with $j=$ 1, …, $n$ represents the $j$th independent component or original source, i.e., the IC$j$; $\mathbf{A}$ is an $m$×$n$ mixing matrix, and its element $a_{ij}$ for $i=$ 1, …, $m$; $j=$ 1, …, $n$ is determined by the transfer path between the $i$th sensor and the $j$th source. According to the principles of matrix operations, if $\mathbf{x}$ and $\mathbf{A}$ can be obtained, the original sources or ICs can be obtained by:

where $\mathbf{W}$ is called the separating matrix, is the inverse of matrix $\mathbf{A}$.

3.2. Further discussion on ICA

In vibration analysis, the observed vector $\mathbf{x}$ can be obtained directly from a multi-channel synchronous acquisition by sensors, and the data from each channel represents an element (mixture) of $\mathbf{x}$. Generally, the matrix $\mathbf{W}$ is not obvious and the solution for $\mathbf{s}$ may not exist, but the ICA is able to estimate the separating transfer matrix $\mathbf{W}$ based on the statistical properties of the element $s_j$ in the vector $\mathbf{s}$. Once $\mathbf{W}$ is obtained, the vector $\mathbf{s}$ can be calculated by Eq. (4). It is worth pointing out that there are two key constraints on utilizing the separability of the popular algorithms. Firstly, the dimension of independent components is required to be not more than the number of the observed vectors viz. $m\le n$ in the standard ICA algorithms. Secondly, the independent component must have non-Gaussian distributions (at most one of the sources $s_j$ is Gaussian). Please refer to [5] for more details about the ICA.

4.1. Difficulty on gearbox faults identification at early stage

As aforementioned, if several REBs and gears faults exist simultaneously in a gearbox, then one of these faults corresponds to a stronger vibration level. The results from conventional envelope analysis often miss the characteristic frequencies of the weaker impulsive trains from incipient faults of REBs or gears. To address this issue, some schemes have been proposed such as the discrete/random separation (DRS) [6] and the cepstral editing procedure (CEP) [7]. In this paper, a new approach is proposed and introduced in next subsection.

4.2. Envelope separation scheme

Envelope analysis is a powerful tool for obtaining the characteristic frequencies of weak impulse trains generated by incipient faults of REBs and gears in gearboxes, and the is a novel method for its separation ability based on the observed mixtures. Intuitively, the combination of envelope and the ICA should contribute to obtain the characteristic frequencies from a gearbox where multiple faulty sources exist simultaneously. Inspired by this idea, this study has been investigated. The main steps of the proposed scheme are listed as follows.

1) Multi-channel data acquisition. The observed vibrations are synchronously acquired from a gearbox by the sensors (three piezoelectric acceleration sensors are employed in our study).

2) Envelopes extraction. The kurtogram scheme is utilized for the envelopes extraction. The envelope from each channel is extracted, respectively. It is worth mentioning that after this step, the strong noise beyond the demodulation frequency band is isolated and the incipient faults related impulsive components are enhanced.

3) Envelope ICs separation. As aforementioned, in the conditions where more than one fault of rolling element bearings (REBs) and gears exist simultaneously in a gearbox, the conventional envelope analysis may be invalid. To address this issue, the ICA is employed to separate the envelope mixtures according to the natural properties of the sources. After this operation, the faults related envelope mixtures from different sources in the same gearbox can be separated into different envelope ICs, which are corresponding to different sources, such as the faulty REB, the faulty gear and disturbances etc.

4) Squared envelope spectrum analysis. The squared spectrum analysis is applied to each envelope IC, and the corresponding characteristic frequencies can be exposed clearly for the diagnosis of the gearbox.

5.1. Simulation

Firstly, the vibration sources in the gearbox are simulated, which include the impulse trains generated by a faulty REB and a faulty gear respectively, the shaft related low-frequency harmonics, the gear mesh vibration and a Gauss noise.

A single impulsive response for simulating the vibration generated by a REB with a single local fault is given by:

where $r\left(t\right)$ is the impulsive response, $f_r=$8,000 Hz denotes a resonant frequency of the system, the constants 0.05 and 500 are the amplitude and the damping coefficient of the response respectively. Then, we repeat the $r\left(t\right)$ 92 times per second (92 Hz) to simulate the impulse train $R\left(t\right)$ generated by the faulty REB. In the same way, an impulsive response for simulating a gear with a single local fault is generated, which is defined by:

where $g\left(t\right)$ is the impulsive response, the constant 0.5 defines the initial phase of the response respectively. Then, by repeating the $g\left(t\right)$ 13 times per second viz. setting the shaft rotating frequency $f_1=$13 Hz, one can simulate the periodic vibration $G\left(t\right)$ generated by the faulty gear. The shaft related harmonic vibration is given by:

where $H\left(t\right)$ denotes the harmonic vibration. The gear mesh vibration is given as in [8] by:

where $M\left(t\right)$ denotes the gear mesh vibration, $f_m$ is the gear mesh frequency, $A_n$, ${\phi }_n$, $a_n\left(t\right)$ and ${\varphi }_n\left(t\right)$ are the amplitude, the initial phase, the amplitude modulation function and the phase modulation function of the $n$th harmonic of the mesh vibration respectively. In the study, we assume that $f_m=$ 413 Hz and the mesh vibration includes 3 harmonics, $N=$ 3, and let:

The Gauss white noise $n\left(t\right)$ is simulated with power spectral density (PSD) –5 dB. Then, the sources in the gearbox are simulated, $\mathbf{s}={\left\{R\left(t\right),G\left(t\right),H\left(t\right),M\left(t\right),n\left(t\right)\right\}}^T$, which include five different sources. By setting the sampling rate as 40,960 Hz, the waveforms of the five simulated different sources are shown in Fig. 2 respectively.

Fig. 2Waveforms of simulation sources

Secondly, Eq. (3) is used to generate the observed mixtures by setting the values of the elements in the matrix $\mathbf{A}$ randomly. In this study, three simulated vibrations are generated by:

The waveforms of the observed mixtures are shown in Fig. 3. Then one of the observed mixtures ($x_1\left(t\right)$ is selected in the study) is used to calculate the kurtogram by the fast kurtogram algorithm [4] for obtaining the optimal demodulation band in envelope analysis according to Eq. (2). The kurtogram is shown in Fig. 4.

In Fig. 4, the optimal band ($f_o$,$∆f_o$) determined by the kurtogram is (8,000 Hz, 213.33 Hz). The center frequency ($f_o=$8,000 Hz) for the band-pass filtering in envelope analysis is consistent with the setting resonant frequency ($f_r=$8000 Hz) in Eq. (5) and Eq. (6). Consequently, the envelopes of $x_1\left(t\right)$, $x_2\left(t\right)$ and $x_3\left(t\right)$ are extracted respectively, and the squared envelope spectra are shown in Fig. 5. It is clear that the setting characteristic frequency 13 Hz (12.96 Hz in Fig. 5) of the simulated faulty gear and its harmonics are all exposed. However the characteristic frequency of the simulated faulty REB is missed.

Fig. 3Waveforms of simulated observed mixtures

Fig. 4Kurtogram of x1(t), (SKmax=1.2 at level 6.5, ∆fo= 213.33, fo= 8000 Hz)

Fig. 5Squared envelope spectra of simulated observed mixtures

By utilizing the proposed method, the three envelopes are separated into three envelope ICs by the fastICA algorithm [5] at first, then the squared envelope spectra of envelope ICs are calculated. The results are shown in Fig. 6.

The characteristic frequency of the simulated faulty gear 13 Hz (12.97 Hz in Fig. 6) and the characteristic frequency of the simulated faulty REB 92 Hz (91.97 Hz in Fig. 6) are clearly identified in the spectra corresponding different envelope ICs simultaneously.

Fig. 6Squared envelope spectra of envelope ICs

5.2. Experimental test on a gearbox

The picture of the gearbox test rig is shown in Fig. 7. A faulty REB with outer race fault and a faulty gear with a broken tooth have been machined artificially and mounted on the input shaft end. The REB type is N1007 (the number of rolling elements $Z=$ 15, the contact angle $\alpha =$ 0, the diameter of the rolling element $d=$ 6.75 mm, and the pitch circle diameter $D=$ 48.5 mm). The rotating speed of the input shaft is about 600 rev/min (the rotating frequency of the input shaft $f_{r1}=$10 Hz). The faulty REB is mounted on the input shaft.

Theoretically, the ball pass frequency of outer race (BPFO) of the faulty REB can be calculated as in [9] by:

Fig. 7Gearbox test rig

The broken tooth gear is mounted on the output shaft, which is running at the rotating speed about 400 rev/min (the rotating frequency of the output shaft $f_{r2}=$6.67 Hz). A NI USB9234 (4-channel synchronous sampling) is utilized for data acquisition with the sampling rate 20 kHz; three piezoelectric accelerometers (140609R#, sensitivity 12.6 pC/g, 140608R#, 12.9 pC/g and 140607R#, 12.5 pC/g) are mounted on the bearing house of the input shaft (see Fig. 7) and a charge amplifier is employed.

To evaluate the proposed scheme, comparisons of the conventional envelope analysis scheme and the proposed scheme have been conducted and introduced as follows. The kurtogram of vibration from the channel 1 is shown in Fig. 8, where $S{K}_{max}=$3 is located at level 5, and the optimal combination ($f_o$,$∆f_o$) are $f_o=$4843.75 Hz and $∆f_o=$312.5 Hz.

Fig. 8Kurtogram of vibration from channel 1, (SKmax= 3 at level 5, ∆fo= 312.5, fo= 4843.75 Hz)

By performing the conventional envelope analysis without the ICA, the obtained squared envelope spectra are shown in Fig. 9, and it can be found that only the characteristic frequency harmonics corresponding to the faulty gear can be clearly found, where 13.43 Hz is about the 2× component of the output shaft speed. However, the characteristic frequencies related to the faulty REB are hardly to identify.

Fig. 9Squared envelope spectra of vibrations form channel 1, channel 2 and channel 3 by conventional scheme

By applying the proposed scheme, the obtained squared envelope spectra are shown in Fig. 10. The characteristic frequencies of the faulty gear can be found in the spectrum of the envelope IC₁, and the characteristic frequencies of the faulty REB, 65.64 Hz component and its harmonics can be clearly identified in the spectrum of the envelope IC₂ simultaneously. It is worth mentioning that the rotating speed of the input shaft is set on the control panel, and the rotating speed of the output shaft is calculated by the gear ratio, so that the actual rotating speed of the input shaft and the output shaft are slightly different to the setting value. This may explain why the characteristic frequencies of the faulty gear and the characteristic frequencies of the faulty REB are slightly different to the ideal values.

Fig. 10Squared envelope spectra of envelope IC1, envelope IC2 and envelope IC3 by proposed scheme