Gear fault feature extraction and classification of singular value decomposition based on Hilbert empirical wavelet transform

2.1. Time-frequency signal decomposition based on HEWT

The HEWT is a merger of the empirical wavelet transform and Hilbert transform [26]. The EWT is used to extract adaptive modes from the vibration signal. Then, the instantaneous amplitude and frequency are performed for each mode using HT.

2.1.1. Empirical wavelet transform

The idea of EWT is to construct a set of $N$ wavelet filters (one Low pass and ($N-1$) band pass filters) capable of extracting Amplitude Modulated-Frequency Modulated components that is AM-FM components (modes) $f_k\left(t\right)$ of an input signal $f\left(t\right)$ by adaption from the processed signal [23] such as:

1

$f\left(t\right)=\sum _{k=0}^N{f}_k\left(t\right).$

For the adaptation process, the wavelet filter bank is based on the Fourier supports detected from the information contained in the processed signal spectrum by finding the local maxima then taking support boundaries ${\omega }_i$ as the middle between successive maxima [23].

As it shown in Fig. 1, the Fourier support [0, $\pi$] is partitioned into $N$ contiguous segments. Each segment is denoted as ${\mathrm{\Lambda }}_n=\left[{\omega }_{n-1},{\omega }_n\right]$, thus ${\bigcup }_{n=1}^N{\mathrm{\Lambda }}_n=\left[0,\mathrm{}\mathrm{}\pi \right]$. Centered around each ${\omega }_n$, there is transient phase $T_n$ with width $2{\tau }_n$.

The empirical wavelets are defined as band pass filters on each ${\mathrm{\Lambda }}_n$ and based on ${\mathrm{\Lambda }}_n$ a wavelet tight frame can be defined. Hence, inspired by the Meyer’s and Littlewood-Paley wavelets, a wavelet tight frame $B=\left\{{\left\{{\varphi }_n\left(t\right)\right\}}_{n=1}^{N-1},\mathrm{}\mathrm{}{\left\{{\psi }_n\left(t\right)\right\}}_{n=1}^{N-1}\right\}$ is defined. And for arbitrary $n>$ 0, their Fourier transforms, i.e. the empirical scaling function ${\widehat{\varphi }}_n\left(\omega \right)$ and the empirical wavelets ${\widehat{\psi }}_n\left(\omega \right)$ are given by the expressions of (Eqs. (2) and 3), respectively:

2

${\widehat{\varphi }}_n\left(\omega \right)=\left\{\begin{array}{ll}1,& \left|\omega \right|\le {\omega }_n-{\tau }_n,\\ \mathrm{c}\mathrm{o}\mathrm{s}\left[\frac{\pi }{2}\beta \left(\frac{1}{2{\tau }_n}\left(\left|\omega \right|-{\omega }_n+{\tau }_n\right)\right)\right],& {\omega }_n-{\tau }_n\le \left|\omega \right|\le {\omega }_n+{\tau }_n,\\ 0,& \text{otherwise}\text{,}\end{array}\right.$

3

${\widehat{\psi }}_n\left(\omega \right)=\left\{\begin{array}{ll}1,& {\omega }_n+{\tau }_n\le \left|\omega \right|\le {\omega }_{n+1}-{\tau }_{n+1},\\ \mathrm{c}\mathrm{o}\mathrm{s}\left[\frac{\pi }{2}\beta \left(\frac{1}{2{\tau }_{n+1}}\left(\left|\omega \right|-{\omega }_{n+1}+{\tau }_{n+1}\right)\right)\right],& {\omega }_{n+1}-{\tau }_{n+1}\le \left|\omega \right|\le {\omega }_{n+1}+{\tau }_{n+1},\\ \mathrm{s}\mathrm{i}\mathrm{n}\left[\frac{\pi }{2}\beta \left(\frac{1}{2{\tau }_n}\left(\left|\omega \right|-{\omega }_n+{\tau }_n\right)\right)\right],& {\omega }_n-{\tau }_n\le \left|\omega \right|\le {\omega }_n+{\tau }_n,\\ 0,& \text{otherwise}\text{.}\end{array}\right.$

The function $\beta \left(x\right)$ is given as follows:

4

$\beta \left(\omega \right)=\left\{\begin{array}{ll}0,& \mathrm{}x\le 0,\\ \beta \left(x\right)+\beta \left(1-x\right)=1,& \forall x\in \left[0,\mathrm{}1\right],\\ 1& x\ge 1.\end{array}\right.$

Note that the most used function that satisfies this property is:

5

$\beta \left(x\right)=x^4\left(35-84x+70x^2-20x^3\right).$

Fig. 1Fourier axis segmentation and EWT wavelets construction

Once the tight frame set of empirical wavelets is built, the definition of EWT, ${\mathcal{W}}_f^{\epsilon }\left(n,t\right)$ can be defined in the same way as the classical wavelet transform. The detail coefficients are obtained by the inner products of the input signal with the empirical wavelets [23]:

6

${\mathcal{W}}_f^{\epsilon }\left(n,t\right)=〈f,{\psi }_n〉=\int f\left(\tau \right)\overline{{\psi }_n\left(\tau -t\right)}d\tau ,$

7

$\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}={\left(\widehat{f}\left(\omega \right)\overline{{\widehat{\psi }}_n\left(\omega \right)}\right)}^{\vee }.$

The approximation coefficients are obtained by the inner product with the scaling function:

8

${\mathcal{W}}_{\mathrm{}f}^{\epsilon }\left(0,t\right)=〈f,{\varphi }_1〉\mathrm{}=\int f\left(\tau \right)\overline{{\varphi }_1\left(\tau -t\right)}d\tau .$

So, the signal is decomposed into various empirical modes $f_k$, which is given by:

9

$f_0\left(t\right)={\mathcal{W}}_f^{\mathrm{}\mathrm{}\epsilon }\left(0,t\right)\mathrm{*}{\varphi }_1\left(t\right),$

10

$f_k\left(t\right)={\mathcal{W}}_f^{\mathrm{}\epsilon }\left(k,t\right)\mathrm{*}{\psi }_k\left(t\right).$

The EWT is invertible and the signal can be reconstructed as follows:

11

$f\left(t\right)={\mathcal{W}}_{\mathrm{}f}^{\epsilon }\left(0,t\right)\mathrm{*}{\varphi }_1\left(t\right)+\sum _{n=1}^N{\mathcal{W}}_f^{\epsilon }\left(n,t\right)\mathrm{*}{\psi }_n\left(t\right).$

As the initial goal of the EWT is to get a decomposition as depicted in Eq. (1). Comparing the Eq. (1) with (11), we can deduce that each mode $f_k$ in Eq. (1) corresponds to Eqs. (9) and (10).

2.2. Singular value decomposition on the HEWT

The singular value decomposition is a matrix transformation algorithm that decomposes any given matrix $M$ ($m\times n$) into three matrices $U$, $\mathrm{\Sigma }$ and $V$ as follows:

where $U$ ($m\times m$) and $V$ ($n\times n$) are orthogonal matrices and $\mathrm{\Sigma }$ is an ($m\times n$) diagonal matrix of singular values, ${\sigma }_1\ge {\sigma }_2\cdots \ge {\sigma }_n$. The matrix $\mathrm{\Sigma }$ is represented as:

Fig. 2Feature extraction process using the proposed approach

The columns of the orthogonal matrix $U$ are called the left singular vectors and the columns of the matrix $V$ are called the right singular vectors of the matrix $M$. An important property of $U$ and $V$ is that they are mutually orthogonal [30].

As seen in Fig. 2, after decomposing the signal using EWT, the analytical signal of each extracted mode is obtained using HT (see Eq. (12)). Subsequently, the instantaneous amplitude $A_i\left(t\right)$ ($i=$1,…, $n$) is obtained using Eq. (14). The feature matrix $M$ is therefore constructed from the instantaneous amplitude $A_i\left(t\right)$, and the fault feature vector ($\sigma \left(M\right)=\left[{\sigma }_1,{\sigma }_2,\dots ,{\sigma }_n\right]$) by applying the SVD.

The use of SVD presents its own advantages. It is able to expresses the feature matrix $M$ in the form of several values (singular values), so it is endowed a dimension reduction strategy. Furthermore, the singular values have a good stability. In other words, when the feature matrix element changes, a large variance of its singular values does not occur.

2.3. Elman neural network

The Elman neural network is generally divided into four layers, including input layer, hidden layer, association layer and output layer. The basic structure of the network is seen in Fig. 3. The association between context layer and hidden layer, makes the neural networks sensitive to the history of input data. The context neurons can be treated as the memory units, so, this internal feedback improves significantly the capacity of the network to deal with dynamic information, overcoming the drawback of the feed-forward network. To this end, the Elman neural network compared with other neural networks, is robust in fault classification.

The transfer function of hidden layer is a non linear function, which is generally using a Sigmoid function. The activation function of the output layer neuron is a linear function. The mathematical model of Elman neural network is analyzed as follows:

where $y$, $x$, $u$, $x_c$ are separately representing the output vector of the network, the output vector of the hidden layer, input vector of the network and the feedback state vector. Then ${\omega }_3$, ${\omega }_2$, ${\omega }_1$ are separately representing connection weight matrix from the hidden layer to the output layer, from the input layer to the hidden layer and from the context layer to the hidden layer. $f\left(\cdot \right)$ is the transfer function of the hidden layer and $g\left(\cdot \right)$ denotes the activity function of the output neurons.

Fig. 3Structure of Elman neural network

3.1. Experimental system description

Gear vibration testing experimental apparatus is shown in Fig. 4. The rotation motion of the equipment is generated by an electric dc-motor controlled in rotational speed with a nominal speed of 3600 RPM. Torque will be transmitted to the gearbox through the coupling where several pinion fault configurations were assembled. After the reducer outputting, through gear coupling, the torque will be transferred to a magnetic powder brake capable of generating different resistive torques [31].

Fig. 4a) Experimental gearbox test rig, b) structure of the single stage gear in the gearbox

In order to verify the effectiveness of the proposed method, six pinions with different fault states were considered. The first one is referred as Good (G), whereas the others have several different types of faults: a Tooth Root Crack (TRC), a Chipped Tooth in Length (CTL), a Chipped Tooth in Width (CTW), a Missing Tooth (MT) and General Surface Wear (GSW) as shown in Fig. 5. Three pinions are simultaneously mounted on the input shaft of the gearbox, the engagement of each of them is done by a simple axial movement of the wheel on the output shaft (Fig. 3(b)). To record vibration signals, two accelerometers with a sensitivity 100 mV/g was mounted radially, one vertically and the other horizontally on the bearing case of the output shaft. The time sampling frequency of the accelerometer channels is 125 kHz. The cut-off frequency of the anti-aliasing filter is 27 kHz. The acquisition duration is 30 s.

Fig. 5Six pinion states: a) normal, b) tooth root crack, c) chipped tooth in length, d) chipped tooth in width, e) missing tooth, f) general surface wear

a)

b)

c)

d)

e)

f)

The accelerometer signals have been recorded for four different operating modes: under the state of loading 8 N.m: 900 RPM, under the state of loading 11 N.m: 900 RPM, under the state of loading 5 N.m: 1200 RPM, under the state of loading 8 N.m: 1200 RPM.

Fig. 6 shows the acceleration vibration signals collected from pinions with different fault modes for an operating speed 1200 RPM and under the state of loading 5 N.m.

It can be seen that the TRC and the CTL faults (Fig. 6(b, d) respectively) do not significantly affect the vibration signal. On the other hand, a significant increase in the time signal energy is caused by the CTW, MT and GSW faults (Fig. 6(c, e, f) respectively) and as for localized defects, the presence of a repetitive shock waves at every revolution period. Similar observations were noted from the vibration signals of the other operating speed and load conditions.

Fig. 6Pinion vibration signals of gearbox with different faults: a) G, b) TRC, c) CTL, d) CTW, e) MT, f) GSW

3.2. Comparison between the proposed HEWT-SVD method with HHT-SVD, LMD-SVD and WPT in tandem with PCA for feature extraction

To verify the robustness of the proposed method, three widely used methods, HHT-SVD [32], LMD-SVD and WPT together with PCA are also implemented and their performances are compared with that of HEWT-SVD for gear feature extraction under variable conditions.

In four different operating modes, vibration signals of a normal gearbox and five kinds of pinion faults are divided into 15 groups of data. Thus, for each operating mode, we end up with 90 samples.

Since the experiment data we used are collected from test bench in laboratory setup, the vibration signals, therefore, present a low noise level, while in real word application, tremendous noise level may corrupt the vibration signals. Thus, strong background Gaussian white noises (GWN) with SNR of 0.5 dB was added to the vibration signals as seen in Fig. 7.

To begin with, we first analyze the various signals shown in Fig. 7 by the HEWT method. An example of the HEWT application for one certain fault and operating modes (gearbox with tooth root crack fault under the operating speed 1200 RPM and 5 N.m of load) is shown in Fig. 8, from which we can display the Fourier spectrum segmentation, the extracted modes and their envelopes.

It can be seen in Fig. 8(a), that the whole spectrum is adaptively divided into six regions. Consequently, six different modes are obtained in total. The analytical signal, thereafter, was computed for each mode by HT. The feature matrix can therefore be constructed for fault identification. Finally, the singular value vector which is the fault feature vector can be obtained by conducting SVD.

Fig. 7Vibration signals of gearbox with different faults corrupted with GWN of SNR = 0.5 dB: a) G, b) TRC, c) CTL, d) CTW, e) MT, f) GSW

Fig. 8HEWT results: a) Fourier spectrum segmentation, b) the extracted modes, c) envelopes of modes

Table 1, partly gives the fault feature values obtained by HEWT-SVD for an operating speed 1200 RPM and under the state of loading 5 N.m, each vector contains the first three singular values.

To make a comparison between the proposed method with HHT-SVD, LMD-SVD and WPT-PCA, they are all implemented on the same data sets.

The EMD (resp. LMD) decomposes the vibration signals into a set of IMFs (resp. product functions (PFs)). As their number is defined by the characteristic and complexity of the analyzed signal, only the first 6 IMFs (resp. PFs) were taken. To perform the HHT, the IMFs are used to get the analytical signal using HT. Next, the feature matrix was constructed and the fault feature vectors were obtained using SVD.

When using WPT-PCA, the signals are decomposed by four layers sym1 wavelet packet decomposition, and signal energy of 16 frequency bands are obtained. Subsequently, PCA is used to extract principal components from the wavelet packet energy feature for dimensionality reduction.

Fig. 9 shows the projections of the first three singular values obtained by HEWT-SVD of the normal gearbox and five kinds of pinion faults in four different operating modes in the three-dimensional plan view, the classification results for different faults of the gearbox can be seen from the figures. Obviously, for one certain fault mode and whatever the operating mode, the singular values extracted by HEWT-SVD have a high degree of coincidence. In the same context, comparing the singular value clusters of the normal gearbox with the five fault modes for different operating modes shown in Figs. 9(a-d), one can see that, even though the signals are noisy, the mode separability is satisfactory, thus a good classification and it looks better for the operating speed 1200 RPM and under the state of loading 5 N.m shown in Fig. 9(c).

The results obtained by HHT-SVD for the operating speed 900 RPM and under the state of loading 8 and 11 N.m respectively shown in Fig. 10(a and b) and those obtained by LMD-SVD for the operating speed 1200 RPM and under the state of loading 5 respectively shown in Fig. 11(c) are almost comparable to the results obtained by HEWT-SVD shown in Fig. 9. A noticeable separability between faults and a good coincidence of the singular value vectors. However, under the remaining operating modes, the classification effect is worst between normal gear and gear with TRC and between gears with CTW and CTL as shown in Figs. 10(c, d) and Figs. 11(a, d) and between normal gear and gear with CTW as shown in Fig. 11(b) from which, the clustering gap is very small or even inexistent that may lead to the misclassification.

Fig. 9Scatter plots of the first three singular values obtained by HEWT-SVD for different faults: a) load 8 N.m, speed 900 RPM, b) Load 11 N.m, speed 900 RPM, c) Load 5 N.m, speed 1200 RPM, d) Load 8 N.m, speed 1200 RPM

Fig. 10Scatter plots of the first three singular values obtained by HHT-SVD for different faults: a) load 8 N.m, speed 900 RPM, b) load 11 N.m, speed 900 RPM, c) load 5 N.m, speed 1200 RPM, d) load 8 N.m, speed 1200 RPM

Fig. 11Scatter plots of the first three singular values obtained by LMD-SVD for different faults: a) load 8 N.m, speed 900 RPM, b) load 11 N.m, speed 900 RPM, c) load 5 N.m, speed 1200 RPM, d) load 8 N.m, speed 1200 RPM

Fig. 12Scatter plots of the first three principal components obtained by WPT-PCA for different faults: a) load 8 N.m, speed 900 RPM, b) load 11 N.m, speed 900 RPM, c) load 5 N.m, speed 1200 RPM, d) load 8 N.m, speed 1200 RPM

On the other hand, comparison between Figs. 12(a, c, d), Figs. 10(a, c, d) and Figs. s11(a, c, d) shows that: utilizing the PCA in the process of fault feature extraction for wavelet packet frequency band energy, the effect is significantly better than the results obtained by HHT-SVD and LMD-SVD. Nevertheless, there exist major fluctuations for all states, the clustering intervals between the faults states are relatively small. Moreover, the effect of classification is worst between gears with TRC and CTW as shown in Fig. 12(b).

Known from the above analysis, using HEWT with SVD for fault feature extraction under different operating modes is much better than HHT-SVD, LMD-SVD and WPT in tandem with PCA.

Since the instantaneous amplitude matrices of the vibration signals of different fault in the pinions are different, adaptively extracted by the HEWT, thus, they can be used as the feature of diagnostic fault. As the goal is to extract more efficient and more accurate fault features for fault recognition and classification, dimensionality reduction was done by the SVD, thus fault can be accurately identified.

It was obvious in the scatter plots that the various fault singles have a large difference and no overlapping portions in Figs. 9(a-d), from which we suggest that the various faults can be distinguished well.

3.3. State classification based on Elman neural network

Here, the fault feature vectors obtained by HEWT–SVD are used to identify and classify the gear states by using the Elman neural network.

The data sets are divided into training dataset with 168 groups (7 groups of data for each fault status) while the remaining data set (192 groups) are used to test the effectiveness and the accuracy of the classifier. State classification results are partly shown in Table 2, from which we can see that even under variable operating modes and even the signals are noisy the actual outputs of the Elman neural network are consistent with the target outputs. Therefore, combining HEWT-SVD with Elman neural network is effective for gear fault diagnosis under variable operating modes.

In order to compare the classification effect, another classification technique Back Propagation (BP) neural network is applied in classification. Comparison results are shown in Table 3 in which 5 examples are given to calculate the mean value of classification accuracy. As shown in the table, the Elman neural network has an advantage over BP in classification accuracy. Actually, both classifiers have classification accuracies higher than 0.9, because the obtained feature vectors have good separability. This gives confirmatory evidence about the effectiveness of the proposed feature extraction method.