Bearing fault diagnosis based on spectrum images of vibration signals

There are many possible choices for the creation of vibration signal images. The images can be obtained in time domain, frequency domain, and time-frequency domain. In this study, we capture the FFT spectrums of vibration signals as images. The x-axis of the spectrum is frequency in Hertz, and the y-axis is the amplitude. For a given signal, the x-axis of the spectrum is determined by the sampling rate. When capturing the spectrum, its y-axis is auto-scaled. Then the parameter of the image is just the size (in pixels). With a larger image, more details of the spectrum can be depicted; while with a smaller image, some details may be lost. The flowchart of image creation in MATLAB is detailed illustrated in figure 1.

Zoom In Zoom Out Reset image size

By taking such images as the features of vibration signals, we actually make use of all information contained in the spectrum – i.e., characteristic frequencies and harmonics of bearings, the geometrical structure of the spectrum, the peak amplitude of the spectrum and so on. The images provide much more useful knowledge of bearings, and we can also take advantage of the well developed computer vision technique to realize fault diagnosis. Next, two-dimensional principal component analysis (2DPCA) is applied to process the obtained spectrum images such that the low dimensional features of the images can be obtained.

2DPCA is a feature projection method similar in principle to the conventional principal component analysis (PCA), using which we can extract the intrinsic information of images with a direct operation on the matrix. The 2DPCA projection process can be summarized as follows [23].

• Step $1:$ Suppose that there are M training image samples, and the jth training image(with $w\times h$ pixels) is denoted by an $w\times h$ matrix ${{A}_{j}},j=1,2,\ldots,M$. The average image of all the training image samples is denoted by $\overline{A}$. Then the global image scatter matrix G can be evaluated by
  $$\begin{eqnarray}&&G=\frac{1}{M}\underset{j=1}{\overset{N}{\mathop \sum}}{{({{A}_{j}}-\overline{A})}^{T}}({{A}_{j}}-\overline{A})\in {{\mathbf{R}}^{h\times h}}\end{eqnarray} \tag{ 1 }$$
  where ${{(\bullet )}^{T}}$ represents the transpose of matrix $(\bullet )$.
• Step $2:$ In order to obtain the basis vectors, it is necessary to find the eigenvectors u_k and eigenvalues ${{\lambda}_{k}}$ of the global image scatter matrix G by solving the following equation:
  $$\begin{eqnarray}&&Gu=\lambda u\end{eqnarray} \tag{ 2 }$$
  where ${{\lambda}_{k}},k=1,2,\ldots,h$, are the eigenvalues, and $u=\left[{{u}_{1}},{{u}_{2}},\ldots,{{u}_{h}}\right]$ are the corresponding eigenvectors of G. To reduce dimensions, and decrease computational expense, we next normalise and sort the eigenvectors u in decreasing order according to the corresponding eigenvalues ${{\lambda}_{k}}$. Then the $u=\left[{{u}_{\text{max}1}},{{u}_{\text{max}2}},\ldots,{{u}_{\text{max}h}}\right]$ and the corresponding eigenvalues $\lambda =\left[{{\lambda}_{\text{max}1}},{{\lambda}_{\text{max}2}},\ldots,{{\lambda}_{\text{max}h}}\right]$ can be obtained. Here ${{\lambda}_{k}}$ satisfies the following constraint:
  $$\begin{eqnarray}&&{{\lambda}_{\text{max}1}}>{{\lambda}_{\text{max}2}}>\ldots >{{\lambda}_{\text{max}h}}\end{eqnarray}$$
• Step $3:$ In order to obtain the optimal projection vectors, the first d(d  <  h) largest eigenvectors are selected to form the projection basis as
  $$\begin{eqnarray}&&U=\left[{{u}_{\text{max}1}}\quad {{u}_{\text{max}2}}\quad \cdots \quad {{u}_{\text{max}d}}\right]\end{eqnarray} \tag{ 3 }$$
• Step $4:$ Feature extraction is implemented with the projection basis obtained in the previous step. For a given image sample B, which is also the same size of $w\times h$ pixels corresponding to A_j, let
  $$\begin{eqnarray}&&{{Y}_{k}}=B{{U}_{k}}\end{eqnarray} \tag{ 4 }$$
  where, ${{U}_{k}}={{u}_{\text{max}k}}$, and $k=1,2,\ldots,d$. Then the projected feature vectors, ${{Y}_{1}},{{Y}_{2}},\ldots,{{Y}_{d}}$ can be obtained - these are called the principal components of the image sample B. At last we can obtain the eigen image of B in the form of
  $$\begin{eqnarray}&&E=\left[{{Y}_{1}},{{Y}_{2}},\ldots,{{Y}_{d}}\right]\in {{\mathbf{R}}^{w\times d}}\end{eqnarray} \tag{ 5 }$$

In order to classify the bearing faults, spectrum images of different faulty bearings must be recognized. Firstly some training image samples are processed through 2DPCA to obtain the corresponding eigen images of vibration signals with different faults. Then a nearest neighbor classification method is utilized for the classification of testing spectrum images.

Suppose that the ith projection feature matrix of the M training image samples is ${{E}_{i}}=[Y_{1}^{(i)},Y_{2}^{(i)},\ldots,Y_{d}^{(i)}]$, where $i=1,2,\ldots,M$, and that of the jth testing image is ${{T}_{j}}={{E}_{j}}=[Y_{1}^{(\,j)},Y_{2}^{(\,j)},\ldots,Y_{d}^{(\,j)}]$. Here we apply Euclidean distance [24] to measure the distance between E_i and T_j as follows

$$\begin{eqnarray}&&{{d}_{i}}\left({{E}_{i}},{{T}_{j}}\right)={{L}_{p=2}}\left({{E}_{i}},{{T}_{j}}\right)=\underset{r=1}{\overset{d}{\mathop \sum}}\Vert Y_{r}^{(i)}-Y_{r}^{\left(\,j\right)}{{\Vert}_{2}}\end{eqnarray} \tag{ 6 }$$

where $\Vert Y_{r}^{(i)}-Y_{r}^{\left(\,j\right)}{{\Vert}_{2}}$ denotes the Euclidean distance between $Y_{r}^{(i)}$ and $Y_{r}^{\left(\,j\right)}$, and $Y_{r}^{(i)}$, $Y_{r}^{\left(\,j\right)}$ are the rth vector of E_i and T_j.

Suppose $L={{s}_{1}},{{s}_{2}},\ldots,{{s}_{N}},\left(N\leqslant M\right)$ is the category label set of the M training samples. In order to classify the jth testing image, it is necessary to find the subscript η, which satisfies

$$\begin{eqnarray}&&{{d}_{\eta}}=\min ~\left({{d}_{i}}\right)\end{eqnarray} \tag{ 7 }$$

Then if the ηth training image is assigned as ${{s}_{k}}\left({{s}_{k}}\in L\right)$, the jth testing image is classified as s_k.

The main procedure of fault diagnosis process based on spectrum images is summarized as follows.

• Data Acquisition: The spectrum images of vibration signals are firstly created through FFT as described in section 2.1. The image database can be constructed with these spectrum images.
• Eigen-images Extraction: Once the image database is set up, the eigen images can be obtained through 2DPCA, as detailed in section 2.2.
• Fault Classification: A testing spectrum image can be classified by comparing with training spectrum images using the method given in section 2.3.

The flow chart of bearing fault diagnosis based on spectrum images is shown in figure 2.

Zoom In Zoom Out Reset image size

The tests are firstly demonstrated with PCA. According to the feature dimensionality reduction criterion in [26], the so-called contribution of selected components with 20%, 40%, 60%, 80%, 90% and 100% are firstly performed to determine the reduced dimension. Thereafter, dimension reduction with 90% is designated in our research, as this was the case in which the average classification rate was the highest. The experimental results are shown in tables 2 and 3.

From table 2 we can see that the classification rates are mostly larger than $90\%$, but the classification rates of several tests are relatively low. For example, we only obtain a classification rate of about 75% when using images under Load1 as training data and images under Load3 as testing data. The actual output with n  =  3 in this case is shown as figure 7. We can see intuitively that the IF (target output is: 2) images are classified incorrectly. Some testing images of IF are categorized into BF (target output is: 3), and the others are categorized into OF (target output is: 4), incorrectly. Through observing and analyzing the spectrum images of IF (Load3) and BF (Load1) carefully, it is plausible that the spectrum images of IF (Load3), BF (Load3) and BF (Load1) looked very similar, which resulted in the low classification rate.

Zoom In Zoom Out Reset image size

Similar tests are carried out on the spectrum images with fault size being 0.021, and the classification results are presented in table 3. We can observe that the overall classification rate is higher than the classification rate with fault size 0.014. It is worth mentioning that an acceptable classification rate can be achieved by using only a single training image (n  =  1).

Having tested the performance of the 2DPCA based method with different dimension reduction, we determined d  =  10 defined in formula (3) as the best selection. The diagnostic results with 2DPCA for fault size being 0.014 and 0.021 are given in tables 4 and 5. Taking Load0 as training samples, the classification rate could reach $100\%$ with n  =  5. When the sampling number of each training class is equal to or greater than 5, most of the test cases could achieve a considerable classification rate in excess of 90%. Unfortunately two cases (marked with $\star$) are still around 75%, and the reason for this is similar to that discussed earlier.

Table 4. The classification rate based on 2DPCA with fault size being 0.014.

# of test	na	Testing data
		Test1(%)	Test2(%)	Test3(%)	Test4(%)	T(s)b
1	1	Load0(99.97)	Load1(99.85)	Load2(99.41)	Load3(94.50)	51
	3	Load0(100.00)	Load1(100.00)	Load2(100.00)	Load3(99.95)	56
	5	Load0(100.00)	Load1(100.00)	Load2(100.00)	Load3(100.00)	63
	10	Load0(100.00)	Load1(100.00)	Load2(100.00)	Load3(100.00)	81
2	1	Load0(85.78)	Load1(100.00)	Load2(100.00)	Load3${{(75.45)}^{\star}}$	50
	3	Load0(88.67)	Load1(100.00)	Load2(100.00)	Load3${{(75.00)}^{\star}}$	56
	5	Load0(91.17)	Load1(100.00)	Load2(100.00)	Load3${{(76.46)}^{\star}}$	63
	10	Load0(92.56)	Load1(100.00)	Load2(100.00)	Load3${{(78.26)}^{\star}}$	81
3	1	Load0(87.79)	Load1(100.00)	Load2(100.00)	Load3${{(77.69)}^{\star}}$	51
	3	Load0(91.56)	Load1(99.96)	Load2(100.00)	Load3${{(80.74)}^{\star}}$	58
	5	Load0(92.24)	Load1(100.00)	Load2(100.00)	Load3${{(80.66)}^{\star}}$	65
	10	Load0(93.45)	Load1(100.00)	Load2(100.00)	Load3${{(83.13)}^{\star}}$	84
4	1	Load0(99.25)	Load1(86.17)	Load2(99.51)	Load3(100.00)	52
	3	Load0(99.55)	Load1(85.94)	Load2(99.47)	Load3(100.00)	59
	5	Load0(99.65)	Load1(84.63)	Load2(99.26)	Load3(100.00)	66
	10	Load0(99.89)	Load1(83.64)	Load2(99.53)	Load3(100.00)	85

It is worth mentioning that the time consumption of 2DPCA is considerably smaller than PCA, especially when n is larger. A detailed comparison is given in table 6. It is clear that 2DPCA method is more efficient than PCA method when using FFT spectrum images for fault diagnosis of bearings.

Table 6. The time consumption diversity of Load0 as training with fault size being 0.014.

Training data	Testing data	n	${{T}_{\text{pca}}}$(s)	${{T}_{2\text{dpca}}}$(s)	$\bigtriangleup T$(s)a
Load0	Load0	1	68	51	17
		3	82	58	24
		5	96	64	32
		10	144	82	62
	Load1	1	71	52	19
		3	85	59	26
		5	95	65	30
		10	145	83	62
	Load2	1	65	51	14
		3	80	57	23
		5	99	64	35
		10	143	82	61
	Load3	1	66	51	15
		3	80	56	24
		5	95	63	32
		10	145	81	64

^a$\bigtriangleup T={{T}_{\text{pca}}}-{{T}_{2\text{dpca}}}$, that is: the time consumption difference between PCA and 2DPCA.

The spectrum image of a given vibration signal is constructed based on FFT. In fact the spectrum is just a vector of amplitudes. Nevertheless it is not easy to mine useful knowledge directly from the vector. The spectrum image is a different view of the vector and provides a new way to extract information. We also carry out similar tests by taking FFT spectrum amplitudes as the features, where PCA and the same minimum distance method are used for fault classification. Similarly, after testing the performance of different dimension reduction with PCA and selecting the best case, the contribution of selected components with 90% is also designated in this section.

The results are shown in tables 7 and 8. Obviously the classification performance using FFT amplitude as features is inferior to that using the spectrum images, especially when the testing load condition is different from the training load condition.

Some other remarks are as follows.

• (1)  
  By comparing tables 2 and 3 with tables 4 and 5, it can be found that bearing fault diagnosis based on 2DPCA can achieve better performance than PCA in most cases.
• (2)  
  Time consumption of processing based on 2DPCA is much less than PCA, especially when the number of training samples per class is large.
• (3)  
  In general, a larger n could obtain a higher classification rate (see tables 2–5)). When using the spectrum image as the feature, an acceptable classification rate can still be achieved with only one single training image.

In order to illustrate the potential application of the proposed method in bearing fault diagnosis, a comparative study between the present work and published literature is presented in table 9. Adopting the same faulty bearing data collected from the Case Western Reserve University [25], most of the previous works considered only single load condition, where the training load condition and the testing load condition are the same. Only a few works evaluated the vibration data of multiple loading conditions. As shown in table 9, bearing fault diagnosis based on SLLEP under load 3hp condition, was carried out to classify bearings with fault size being 0.021in using minimum-distance classifier in [27]. In [28], with FDF as feature, SVMs and fractal dimension were employed to diagnose the bearings with fault size being 0.014in and 0.021in. In [29], feature extraction based on LMD analysis method and MSE was put forward to perform fault diagnosis of bearings under load 3hp condition. In [30], bearings under load 2hp condition, were subjected to fault diagnosis based on multiscale permutation entropy (MPE) and improved support vector machine based binary tree (ISVM-BT). Moreover taking all load conditions into account, the improved distance evaluation technique and ANFIS were also employed to diagnose bearings with seven healthy cases in [31].

However the classification rate of the proposed method can achieve 100% in cases where the training and testing load condition are the same, as shown in tables 4 and 5. And the classification rate is still high when the training and testing load condition are different.