Classification of machinery vibration signals based on group sparse representation

3.1. Group sparse representation model

Suppose a test sample $\mathbf{y}\in R^m$ from class $k$ is expressed as a linear combination of training signal samples in the following formula:

where:

is a training dataset as the concatenation of $n$ training samples of all $K$ object classes, and:

is a combination coefficient vector. The GSR emphasizes that $\mathbf{\alpha }$ should have nonzero elements, corresponding to class $k$, and zero or close to zero elements elsewhere (i.e. ${\mathbf{\alpha }}_k\ne 0$, and ${\mathbf{\alpha }}_i=0$ for $i=1,\dots ,k-1,k+1,\dots ,K$). This “group sparse” problem is equivalent to the following optimization model:

where the mixed norm ${‖\mathbf{\alpha }‖}_{\mathrm{2,0}}={\sum }_{i=1}^{K}I$$({‖{\mathbf{\alpha }}_i‖}_2>0)$, here:

Obviously, solving the problem Eq. (5) is also NP-hard. A convex relaxation optimization substitutes the problem Eq. (5) as:

where ${‖\mathbf{\alpha }‖}_{\mathrm{2,1}}={\sum }_{i=1}^K{‖{\mathbf{\alpha }}_i‖}_2$. From model Eq. (6), it can be seen that GSR seeks a group sparse coefficient vector by imposing L1-norm constraint between the classes for achieving the sparsity among classes and L2-norm constraint within each class for ensuring the homogenous samples can be sufficiently exploited in the representation of the test sample [13].

3.2. Coefficients solving algorithm

The regularized least squares model gotten from Eq. (6) is expressed as:

where $\gamma \in [0,\mathrm{}1]$ is the regularization parameter. Problem Eq. (7) is a group lasso model [24] and can be solved by several existing techniques, such as Block-coordinate Descent (BCD) [15], spectral projected gradient method (SPGL1) [16] and Sparse Learning with Euclidean projection (SLEP) [17], etc. For its computational rapidity and global convergence property, we choose SLEP as a basic tool for solving problem Eq. (7). The algorithm contains the following steps [18].

Firstly, $loss\left(\mathbf{\alpha }\right)=\frac{1}{2}{‖\mathbf{y}-\mathbf{A}\mathbf{\alpha }‖}_2^2$ is denoted as a smooth convex loss function and the following model is constructed for approximating the objective function $f(\bullet )$:

where the first-order Taylor expansion at the point ${\mathbf{\alpha }}_0$ for the smooth loss function $loss(\bullet )$ is applied and the regularization term $\frac{L}{2}{‖\mathbf{\alpha }-{\mathbf{\alpha }}_0‖}_2^2$ prevents $\mathbf{\alpha }$ walking far away from ${\mathbf{\alpha }}_0$.

Secondly, the Armijo Goldstein line search schemes and accelerated gradient descent are used for search the approximate solution. Specifically, suppose ${\mathbf{\alpha }}_i$ is the approximate solution of the $i$th iteration, the search point ${\mathbf{s}}_i$ can be gotten as:

where ${\lambda }_i$ is a tunable coefficient. Then the approximate solution of the $(i+1)$th iteration ${\mathbf{\alpha }}_{i+1}$ is computed as the minimizer of ${{\rm M}}_{L_i,s_i}\left(\mathbf{\alpha }\right)$, i.e.:

where $L_i$ is determined by the Armijo Goldstein line search rule and keeps ${\mathbf{\alpha }}_{i+1}$ in the neighborhood of ${\mathbf{s}}_i$. Problem Eq. (10) can be solved as ${\mathbf{\alpha }}_{i+1}={\pi }_{12}({\mathbf{s}}_i-loss\text{'}({\mathbf{s}}_i)/L_i,\mathrm{\gamma }/L_i)$, here ${\pi }_{12}(\bullet )$ is the L1L2-norm Euclidean projection problem expressed as:

Problem Eq. (11) can decouple into a set of $K$ independent L2-norm ones, i.e.:

where ${\mathbf{\alpha }}_k$ and ${\mathbf{v}}_k$ denote the sub-vector of $\mathbf{\alpha }$ and $\mathbf{v}$ corresponding to the $k$th group $(k=\text{1,}\text{2,}\text{…}\text{,}K)$, respectively. Obviously, the solution of problem Eq. (12) is:

Thirdly, update parameter ${\lambda }_{i+1}$ as:

where $t_i=(1+\sqrt{1+4t_{i-1}^2})/2$. The three steps are repeated lots of times until the objective function $f(\bullet )$ converges to a minimum value. To better showing the solving algorithm of GSR, we summarize it in Table 1

5.1. Simulation analysis

This subsection mainly demonstrates that the proposed method has strong noise immunity and its classification performance is almost unaffected by the sampling starting point. The classification performance includes three factors, i.e. reconstruction errors (RE) (Eq. (19)), SCI (Eq. (20)) [8] and accuracy rate (AR) (Eq. (21)):

Seven types of signals, including a normal signal and six fault signals, are simulated as seven conditions of a rolling element bearing, i.e. no fault, two outer-race faults with different defect sizes, two ball faults with different defect sizes and two inner-race faults with different defect sizes, respectively. The normal signal is set as a white Gaussian noise (WGN) signal with unit standard deviation. For the six fault signals, their model is expressed as:

where $f_A$ denotes amplitude modulation frequency, $\beta$ is attenuation factor, $i$ is impulse number, $f_n$ and $f_0$ are natural frequency and fault feature frequency, respectively. These parameters are set as in Table 2. The seven signals with sampling frequency of 2048 Hz are plotted in Fig. 2.

Each of the seven signals is mixed with 50 different WGN keeping SNR in the range of –2 dB to 2 dB such that 350 (50×7) training samples are obtained. The fault 1 signal is chosen as a test object. At first, it is mixed with WGN to get a noisy signal with SNR of 0 dB. After implementation of the proposed method, the sparse coefficients and RE of the noisy signal are obtained (in Fig. 3), both of which confirm that the test sample is of fault 1 category. Then, the WGN intensity is enlarged to reduce SNR to –2 dB and the proposed method is performed as above.

Although the sparse coefficients spread to three classes and the minimal RE magnifies a lot (in Fig. 4), the test sample still can be classified correctly according to the minimal RE, which shows the strong noise immunity of the proposed method. At last, the fault 1 signal with SNR of 0 dB is circularly shifted right 50 points and tested as above. From Fig. 5, it can be seen that the classification performance is almost unaffected by the sampling starting point.

Fig. 2The waveform of seven simulation signals, which simulate seven conditions of a rolling element bearing (from above to below: no fault, two outer-race faults with different defect sizes, two ball faults with different defect sizes and two inner-race faults with different defect sizes)

Fig. 3The sparse coefficients and reconstruction errors (RE) of the fault 1 test sample with SNR of 0 dB, obtained by the proposed method

Fig. 4The sparse coefficients and reconstruction errors (RE) of the fault 1 test sample with SNR of –2 dB, obtained by the proposed method

Fig. 5The sparse coefficients and reconstruction errors (RE) of the fault 1 test sample with SNR of 0 dB and circularly shifting right 50 points, obtained by the proposed method

Most tests are conducted. 50 signals of each category with SNR in the range of –2 dB to 2 dB are generated as the test samples. So far, the “simulation dataset”, containing 350 training samples and 350 test ones, has been constructed, which will be used again in discussion section. The proposed method is executed to classify the 350 test samples. The results (summarized in Table 5) preliminarily demonstrate the classification performances of the proposed method.

5.2. Classification results of bearing vibration data

To investigate the effectiveness of the proposed method for vibration signal classification, bearing vibration signals are considered in this subsection.

The vibration dataset from Case Western Reserve University Bearing Data Center [19] has been analyzed by lots of researchers [20-21] and considered as a benchmark. In our study, this dataset is adopted as well. The experimental platform consists of a motor, a dynamometer, a torque transducer, and control electronics. Single point faults of size 0.007, 0.014, 0.021 and 0.028 in. were set on the drive-end bearings (Type 6205-2RS JEM SKF) at the location of outer raceway, inner raceway and rolling element (ball), respectively. The vibration data were measured by using an accelerometer being attached to the motor housing with the sampling frequency of 12 kHz.

In this experiment, twelve fault-types of vibration data are chosen to construct training and test datasets, including a normal, three types of outer race fault, four types of inner race fault, and four types of ball fault. Each data is split with an overlapping length of 512-point into lots of segments, whose length are set to 2048-point. From these segments, we randomly select 50 ones as training samples and other 25 ones as testing samples. Totally, 600 (12×50) samples and 300 (12×25) samples with 12 classes of bearing data are used as training set and test one, respectively. The descriptions of the “bearing dataset” are shown in Table 3.

Without loss of generality, we randomly select a test sample from each class and classify them by the proposed method. Their sparse coefficients and RE are plotted in Fig. 6 and Fig. 7, respectively, from which we can see that almost all sparse coefficients of each test sample appear in their own corresponding class and the minimal RE of each test sample just corresponds to its class with value significantly smaller than the other RE. Then, all the 300 test samples are classified by the proposed method. The results (summarized in Table 5), including AR, average SCI and average minimal reconstruction error (AMRE), show that our method can correctly category all the test samples with a high average SCI and relatively low AMRE.

5.3. Classification results of gearbox vibration data

The proposed method is further applied in gearbox fault diagnosis. The experimental platform consists of a motor, a drive shaft seat, a gearbox and a damper, etc. The vibration signals are acquired by an acceleration sensor placed in the output shaft bearing seat. A normal situation and five fault-ones, including three single faults, i.e. tooth-broken and point-corrosion of large gear, wear-out of small gear, and two combination faults, i.e. broken-wear and point-wear faults, are considered. Rotating speed is set as 1500 r/min. The horizontal vibration signals are collected with a sampling frequency of 5120 Hz and the sampling time is about 20 s in each situation.

Fig. 6The sparse coefficients of test samples from each class in bearing dataset

Fig. 7The reconstruction errors of test samples from each class in bearing dataset

Like the experiment of bearing data classification, in the preprocessing stage, each of the six gearbox vibration signals is split with an overlapping length of 128-point into lots of segments, whose length are set to 2048-point. From these segments, we randomly select 50 ones as training samples and other 20 ones as testing samples. Totally, there are 300 (6×50) training samples and 120 (6×20) test samples. The descriptions of the “gearbox dataset” are shown in Table 4.

Without loss of generality, we randomly select a test sample from each class and classify them by the proposed method. Their sparse coefficients and RE are plotted in Fig. 8 and Fig. 9, respectively. From Fig. 8, it can be seen that most of sparse coefficients of each test sample appear in their own corresponding class except for “PCL” and “PCL-WOS” samples. The reason for this phenomenon is that the feature frequency of point-corrosion fault is just equal to rotating frequency and the instantaneous value of vibration signal will almost do not change compared with normal conditions when the fault size is little. The minimal RE of each test sample in Fig. 9 indicates the proposed method can correctly classify it.

The classification results of the 120 test samples are displayed in Table 5. The accuracy rate has reached to 97.5 % (117/120), which verifies that our method can be applied in the gearbox fault diagnosis. Moreover, the high average SCI (0.9723) demonstrates again that group sparse representation based classification method can represent a test vibration signal using only signals from the same fault-type.

6.1. The influence of parameters

In our proposed method, two important parameters should be considered. One is the total points of DFT $N$, and the other is the regularization parameter of GSR $\gamma$. In all of the tests above, including simulation analysis and two experiments, $N$ is set equally to the length of samples (i.e. $N=m=\text{2048}$) and $\gamma$ is set as 0.5.

Fig. 10 shows the relationship between classification accuracy rate and the total points of DFT. It can be seen that the accuracy rate is improved significantly with the increase of the total points of DFT. Theoretically, the more DFT points, the more spectral information we can get. Therefore, we can improve the accuracy rate by increasing DFT points. However, group sparse representation will need more computational time to complete when the number of DFT point increases. For machinery fault diagnosis, vibration signal is analyzed in real-time and diagnostic conclusion should be determined in the crucial early hours. This requires the computational time of a diagnosis method as shorter as possible. As far as the proposed method, it needs to balance the classification accuracy rate and computational time to determine DFT points. In general, we can set DFT points as the half of sample length if accuracy rate is in range of acceptance.

Fig. 8The sparse coefficients of test samples from each class in gearbox dataset

Fig. 9The reconstruction errors of test samples from each class in gearbox dataset

The relationship between classification accuracy rate and the regularization parameter $\gamma$ is visualized in Fig. 11, which demonstrates that the accuracy rate almost keeps unchanged with a high value when the parameter $\gamma$ varies in [0.0001, 0.8], while the accuracy rate will drop sharply if $\gamma$ approaches to 0 or 1. In other words, when group sparse representation are utilized in classification of machinery vibration signals, the influence of the regularization parameter $\gamma$ is very limited when it is set in a reasonable range.

Besides $N$ and $\gamma$, the termination tolerance $\epsilon$ and sample length $m$ have some effect on the classification results as well. In all of the tests above, $\epsilon$ and $m$ is set respectively as 0.001 and 2048. From Table 1, it can be seen that $\mathrm{\epsilon }$ directly determines the runtime of coefficient solving algorithm and the reconstruction error. For classification-based diagnosis method, the reconstruction precision is not as important as for other signal approximation problems. Therefore, we set the termination condition of coefficient solving algorithm as the difference of objective function of two adjacent iterations less than $\mathrm{\epsilon }$to reduce runtime. The parameter $m$ has similar effect with $N$ on classification accuracy rate and computational time. In practical application, if the training samples of each class are enough, a time-series vibration signal should no longer be split into a number of segments unless its length is too long. In this case, $N$ should be regulated rather than $m$.

Fig. 10The accuracy rate of the three datasets with the change of DFT points N

Fig. 11The accuracy rate of the three datasets with the change of regularization parameter γ

6.2. Comparison of the proposed method with some other methods

In this subsection, the proposed method is compared with SRC and SVM to demonstrate its superiority.

The L1-norm based SRC [8] is utilized to classify the test samples of three datasets for comparison. It is known there have been many effective methods for solving L1-minimization problems. Here, l1ls [20], SLEP [17] and GPSR [23] are adopted respectively to solve the problem. From the results in Table 5, it is seen that their classification performances are not as good as that of the proposed method. Specifically, AR and average SCI obtained by the three methods are respectively lower than that of proposed method roughly 10 and 50 percent, while average minimal RE is higher 20 percent or so.

For further comparison, the three datasets are tested by SVM method. SVM is a pattern recognition classification algorithm indicating favorable generalization performance based on statistical learning theory [24]. Fault classifications are achieved with these methods as characteristic parameters (e.g., root mean square value, kurtosis, mean square frequency, etc.) from time and frequency domains are extracted and adjusted based on the actual situation [5]. Modeled on reference [25], vibration signals are translated into different frequency bands by wavelet package transform (WPT); then, the optimal features are selected based on the distance evaluation technique from the statistical characteristics of raw signals and wavelet package coefficients, and the energy characteristics of decomposition frequency band, finally, the optimal features are input the SVM ensemble with SVM toolbox [26] to identify the fault type. The results recorded in Table 5 show that the AR obtained by the SVM ensemble method is improved significantly compared with that of SRC, yet it is still lower than that of the proposed method.

In machinery fault diagnosis, classification efficiency is another important factor. We investigate the computational time (CT) of these methods on the three datasets and tabulate the results in Table 5. Matlab 2012b is run with a computer of CPU: 3.6 GHz, RAM: 4G. From the observations, the proposed method is as fast as SRC solved by SLEP and GPSR and significantly faster than SVM ensemble and SRC solved by l1l s.