The diagnosis of rolling bearing based on the parameters of pulse atoms and degree of cyclostationarity

2.1.1. The construction of pulse dictionary

The rotation of the bearing with localized defects generates periodic pulses [2] which induce some modes of resonances decaying exponentially, and these pulse responses reflect the status of bearings, therefore, extracting or recognizing these pulse responses is the main task in the detection of the defect. The pulse response signals could be simplified as the damped sinusoidal functions [8], and in mechanical systems this kind of signal exist widely and are usually generated by pulses, thus if the atom which is used to represent this kind of signal is termed as “pulse atom” [9, 10], it will be more concrete, clear, and easy to understand. The pulse atom is as follows [8-10]:

where $t$ is time, $T_s$ is the duration of the signal $x$. $g_{\gamma }$ is a normalized pulse atom indexed by the parameter group $\gamma$, here $\gamma =\left(d,\tau ,f,\mathrm{\varnothing }\right)$, so $g_{\gamma }=g_{\left(d,\tau ,f,\mathrm{\varnothing }\right)}$. $s$ is the normalization factor. $d$ is the damping coefficient, it denotes the decaying speed of pulse. $\tau$ is the lag of the pulse, which represents the starting point of the atom. $f$ is oscillation frequency, it indicates the damped natural frequency of the system. The inner product $\left|〈g_{\gamma },x〉\right|$ is the amplitude of the pulse, represents the energy of the atom, and $\mathrm{\varnothing }$ is the phase.

When the pulse dictionary is constructed, in the first place, appropriate ranges are chosen for each parameter according to bearing’s geometry, followed by a discretization. Let $N_d$, $N_f$, $N_{\mathrm{\varnothing }}$, $N_{\tau }$ represents the number of elements of $d$, $f$, $\mathrm{\varnothing }$, $\tau$ in dictionary, after discretization, respectively, it is easy to know that $N_{\tau }=N$, hence the total number of atoms in the dictionary is $N_d{N}_f{N}_{\mathrm{\varnothing }}N$. Apparently, the greater the value of $N_d$, $N_f$, $N_{\mathrm{\varnothing }}$ or $N_{\tau }$, the bigger the dictionary is, the sparser the signal represents. Consequently, it is necessary to choose a wide atomic dictionary for better decomposition results.

2.1.2. The problems of MP_PA implementation procedure

MP introduced by Mallat and Zhang selects atoms one by one from the dictionary [4], in each step the most suitable atom will be selected and added to the approximation model until a desired energy precision or a prespecified number of iterations is approached. In matching pursuit process of signal $s$, the atom $g_{\gamma i}$, namely the appropriate parameter group $\gamma =\left(d,\tau ,f,\mathrm{\varnothing }\right)$, which satisfies Eq. (2) is repeatedly:

where $r_i$ is the residual after $i$ iterations, and $G$ is the set of the parameter group $\gamma$.

This is a problem of solving the maximum value, in order to choose the suitable atom, the inner product of each atom and residual should be calculated and compared with each other, and it usually needs many iterations before up to the termination condition, resulting in a large amount of computation. Reference [4] has proved that, when the length of the signal is fixed, with the increase of iteration steps, the residual ${‖r_i‖}_2$ decays exponentially, therefore, $s\approx \sum _{i=0}^{\mathrm{\infty }}〈r_i,g_{ri}〉g_{ri}$. Fig. 1 is the measured vibration signal of bearing with outer race fault, Fig. 2a, b displayed the functional diagram of inner product $A=\left|〈g_r,x〉\right|$ varied with the parameter group $\gamma =\left(d,\tau ,f,\mathrm{\varnothing }\right)$, obviously, this function is a complex multimodal function.

Fig. 1The vibration signal of bearing with outer race fault

Fig. 2The diagram of the function of inner product: a) the function varies with τ and f, where ∅ and d are fixed, b) the function varies with ∅ and d, where τ and f are fixed

a)

b)

Therefore, the decomposition of vibration signals of bearings by MP_PA is solving the maximum value of a non-linear, non-differentiable, continuous and complex multimodal function, which cannot be directly solved by the techniques depending on gradient [9]. In order to meet the practical application, artificial intelligent algorithms were usually chosen to simplify the process of MP_PA, for example, genetic algorithm (GA) [6].

Only accurate extraction of features from signals especially for weak fault signals, the correct diagnosis could be achieved. NPSO which has received much attention in recent years could maintain good biodiversity and reduce the probability of falling into local optimal solutions, and it performs excellently for this kind of complex multimodal function [11, 13, 14].

2.2.1. The implement of MP_PA by NPSO

The implement of MP_PA by NPSO is simple, namely, in each iteration of MP_PA, the atom which satisfies $\sup \left|〈r_j,g_{\gamma }〉\right|$ is selected by NPSO, here $j$ is the iteration number of MP, and $r_j$ is the residual signal. In NPSO, the position of $i\mathrm{t}\mathrm{h}$ particle is $X_i=\gamma =\left(d,\tau ,f,\varnothing \right)$, the fitness function is $f\left(X_i\right)=\left|〈r_j,g_{X_i}〉\right|$, after searching by particles, the position of the particle with the maximum fitness will be the parameter group $\gamma$ which we are solving, i.e., the desired atom is obtained.

On the foundation of the features of MP_PA, the optimization in this paper was conducted by FERPSO (Fitness Euclidean-distance ratio PSO, a kind of NPSO) [11] and local optimization [13] and their excellent performance has been proved in reference [13]. The following briefly describes their improvements relative to PSO, the detailed explanation for this optimization procedure can be referred to [11, 13].

FERPSO is an effective local-best niching PSO which converges quickly and adopts fitness Euclidean-distance ratio (FER) value as an evaluation index, the difference from PSO is the velocity update mode of the particles, FERPSO uses a local optimum position which is determined by the largest FER value, instead of the global optimum position in velocity update equation. The FER value of particle $j$ corresponding to particle $i$ is determined by the relationship of their personal best positions previously visited, FER value is calculated using the following equation:

where $a=\frac{‖s‖}{f\left(p_g\right)-f\left(p_w\right)}$ is a scaling factor, $‖s‖$ is the Euclidian norm of the distance between the upper and lower of the search space, $f$ is fitness function. $p_g$ is the fittest particle in the current population, whereas $p_w$ is the worst fit particle; ${pbest}_j$ and ${pbest}_i$ are the current personal best positions of particle $j$ and particle $i$, respectively. In FERPSO, the velocity update equation of each dimension for particle $i$ is modified as:

where all the symbols are the same as that in PSO except ${nbest}_i$. $t$ represents the current evolution generation count, $v_i$ is the velocity of particle $i$, $w$ is the inertia weight, $c_1$ and $c_2$ are both acceleration constants, $r_1$ and $r_2$ are two random numbers from a uniform distribution within the range of [0, 1], $X_i\left(t\right)$ is the current position of particle $i$, ${nbest}_i\left(t\right)$ is the personal best position of the particle whose FER value is the largest with particle $i$.

In order to further improve the performance of FERPSO, local optimization is added after the update of the velocity of particles, which not only enhances the fine-searching ability of the original niching PSO algorithms but also speeds up the convergence. The local optimization technique is to move ${pbest}_i$ a small step to or away the nearest personal best position of another particle (${pbest\_nearest}_i$), the process is as follows:

Step 1, find ${pbest\_nearest}_i$;

Step 2, if $f\left({pbest\_nearest}_i\right)>f\left({pbest}_i\right)$,

let $Temp={pbest}_i+c_1\bullet r_1\bullet \left({pbest\_nearest}_i-{pbest}_i\right)$;

otherwise, let $Temp={pbest}_i+c_1\bullet r_1\bullet \left({pbest}_i-{pbest\_nearest}_i\right)$.

Step 3, if $f\left(Temp\right)>f\left({pbest}_i\right)$, let ${pbest}_i=Temp$.

To test the performance of MP_PA which is implemented by FERPSO with local optimization, a simulated vibration signal of a bearing with defect was decomposed. In FERPSO, the number of the particles in the population was 30; the max iteration number was 300. The signal $x$ is a simulated vibration signal of a bearing with localized defects; it consists of a pulse sequence $x\text{'}$ (as Eq. (5)) and Gaussian noise with signal to noise ratio (SNR) of 1dB. These signals are showed in Fig. 3a, it shows that the left pulse of signal $x$ is almost completely submerged in noise. After iterating three times of MP_PA using FERPSO, the results are illustrated in Fig. 3b, a comparison is also made by PSO and GA respectively in Fig. 3c, d. The comparison shows that MP_PA using FERPSO extracted the pulse components very well, and the SNR is as high as to 19 dB between the reconstructed signal $x$ and signal $x\text{'}$, and the SNR are 9 dB and 8 dB by PSO and GA, respectively. What’s more, we also test the three methods when the signal $x\text{'}$ was added Gaussian noise with different SNR, and the results, which represent the average values from ten-time calculations, are listed as in Table 1, they furthermore indicate that FERPSO performs better when solving this problem.

where $g_{\left(d,\tau ,f,\varnothing \right)}$ is a normalized pulse atom, the sampling frequency $fs=$12000 Hz.

Fig. 3The original signal and its decompositions: a) the simulated signal x' and signal x; reconstruction and residual by MP_PA simplified by different methods: b) FERPSO, c) PSO, d) GA

a)

b)

c)

d)

2.2.2. The simplification of MP_PA by SJ

When applying MP_PA to decompose a signal directly without the help of artificial intelligence algorithms, the selection of one suitable atom needs the computation complexity as follow:

where $c_N$ represents the computation of calculating one time the inner product of two signals whose length are both $N$. With the same accuracy, if the signal’s length grows to $2N$, and the computation complexity will be:

Obviously, the amount of computation of $c_{2N}$ is double to that of $c_N$, and the amount of computation of $C_2$ is four times to that of $C_1$. If the signal whose length is $2N$ is divided into two signals whose length are both $N$, because the total number of the selected atoms are the same, and then the computation complexity would reduce to one half. It means that decomposing a signal after segmentation will reduce the amount of computation, and the more the segments it has the less computation that needs. For example, if a signal is divided into eight segments with the same length, the computation complex will reduce to 1/8 of the origin. The experimental results in reference [10] show that, even if applying genetic algorithm, such kind of artificial intelligence algorithms to implement MP_PA, almost the same ratio of computation complexity will also be reduced due to segmentation.

However, the fixed length of segmentation of original signal may divide the pulse components which locates at the endpoint of segments into two parts, this may cause error even wrong in process of decomposition (namely reconstruction), so an appropriate segmentation is necessary [17]. There are several segmentation techniques which have been introduced and successfully applied in different areas, such as an improved algebraic segmentation technique [17], the Kramers-Moyal coefficient method [18], statistical methods [19] and so on. In this paper, a simple and reasonable technique is employed, as shown in Fig. 4, let the segments overlap one by one, and make the overlap contain the pulse of the former segment, in joint process the overlap in the former segment will be deleted. The core of this technique is that the pulse which locates at the endpoint of a segment will be decomposed in next segment, and the beginning of one segment was selected adaptively according to the former one. The process is as follows: first, a proper fixed length of segments should be determined; here we let the duration of a segment be about one rotation of the bearing. Firstly, segment 1 is selected and decomposed (namely reconstructed), and then from the end of its reconstruction to forward, we find the point whose amplitude is less than a certain prespecified threshold which is near zero, then delete the part behind the selected point in its reconstruction, and make the point correspond to the selected point in original signal be the beginning of segment 2. Afterwards we decompose next segment one by one, until the whole signal is reconstructed, at last, all the reconstructions of segments are joined together one by one.

Fig. 4The process of SJ, the original signal is the same as Fig. 1

2.3.1. The extraction of characteristic defect frequency based on DCS

It has been proved by Antoni [15], owing to the slight random slip between the rolling element and the cage, the pulse response induced by localized fault is not the same as the process of faulty gears which is of a purely periodic process, but close to second-order cyclostationary process. In the spectral correlation density of a failure bearing, the characteristic defect frequency displays as a prominent peak.

The methods based on cyclostationarity could eliminate the interference of random noises and improve the SNR, they have received more and more attentions in the application of the bearing fault diagnosis [15, 16, 20, 21]. As an important second-order cyclic statistic, DCS is defined on the foundation of cyclic autocorrelation function (CAF) or cyclic spectrum density (CSD) as follows:

where $\left|R^a\left(\tau \right)\right|$ is the CAF, $\alpha$ is cyclic frequency, $\tau$ is delayed time, $S^{\alpha }\left(f\right)$ is the CSD. It is known that DCS indicates the ratio of the energy of cyclic frequency $\alpha$ to the energy of stationary part in the signal; DCS is a physical quantity which is able to completely describe the cyclostationarity. DCS is a function with cyclic frequency $\alpha$ as the single variable, thus it is simpler than CAF or CSD which is a function of two variables. In theory, only if the signal is cyclostationary at the cyclic frequency $\alpha$, it will be that $DCS\left(\alpha \right)>0$, otherwise $DCS\left(\alpha \right)=0$, in actual fact, measured signal is complicated so that DCS is greater than 0 at many cyclic frequencies, in order to recognize the cyclic frequency clearly, the DCS is normalized in this paper.

After MP_PA, the conventional methods reconstructing signal accord to $s\approx \sum _{i=0}^{\infty }〈r_i,g_{i\left(d,\tau ,f,\varnothing \right)}〉g_{i\left(d,\tau ,f,\varnothing \right)}$, but in this way it will also take the high resonance frequencies which are not help to extract characteristic defect frequency into the reconstruction. If we reconstruct signal in accordance with Eq. (9), namely each pulse atom is simplified to a pulse which occurs at the delayed time $\tau$ and the pulse’s amplitude is the same as the amplitude of atom, it should be noted that the reconstruction in this way is termed as pulse sequence, because the reconstruction will filter the resonance components and only preserves the low frequency components, the DCS of pulse sequence will become more simple and easy to identify. What’s more, because the DCS is obtained on the foundation of matching pursuit and DCS, it is termed as MP_DCS in this paper.

Taking into account the calculation error in actual fact, if we directly use the amplitude of the characteristic defect frequency as the feature value, it will inevitably be inconsistent with the practice. Therefore, the weighted sum of the amplitude of the characteristic defect frequency and its neighbors is taken as the feature value, as follows:

where $WS$ is weighted sum, $f$ is characteristic defect frequency, $A$ is amplitude, $a_i$ is weighted factor, and $1>a_1>a_2>a_3>\dots$, $∆f$ is frequency interval which equals the frequency resolution of DCS, and $n$ represents the number of weighted points.

2.3.2. The extraction of additional multi features

Besides the feature of characteristic defect frequency in the vibration signal of bearing with localized defects, there are other features, for example, the pulse energy will be much more than that of normal bearings; the signal will be dominated by high frequency due to some modes of resonance; what’s more, compared with normal bearings, the damping coefficients of pulses in fault ones will also vary. The parameters of pulse atom have explicit physical significances and almost contain all the information of pulses. If the features are extracted directly from these parameters, not only the process is more direct and convenient, but also the multi fault features could be extracted from multi dimensions, therefore, the diagnosis and detection would be more reliable and accurate. Consequently, besides characteristic defect frequency, the pulse components’ mean energy, frequency center and mean damping coefficients were all taken as features for better diagnosis, the equations based on parameters of pulse atoms are as follows:

where $E$ is the mean energy of pulse atoms, $A$ is the amplitude of atom, $n$ is the number of atoms.

where $FC$ is frequency center, $f$ is resonance frequency of the pulse atom.

where $MD$ is mean damping coefficient, $d$ is the damping coefficient of pulse atom.

2.3.3. Fault classification based on SVM

SVM is a machine learning algorithm based on statistical learning theory and structure risk minimization [22], it converts the optimization problem to a convex quadratic programming problem, and overcomes the problems of neural network such as its structure is complex, it is easy to fall into local minima and over learning, and has weak generalization capacity. SVM also has excellent learning ability and adaptability for small samples, and has been widely used in fault diagnosis in recent years [9, 23].

SVM multi-class classifier is generally constructed by multiple SVM two-classifiers combining together [24], for instance, one-against-all algorithm, and one-against-one algorithm. For the pattern recognition of bearing failure, to determine whether the bearing is failure or not is the most important issue, the second is the judgment which kind of failure. Therefore, the one-against-all algorithm is selected in this paper.

3.1.1. The analysis of normal bearing

Fig. 5a is the vibration signal of normal bearing, we cannot determine which state the bearing is from Fig. 5a. Fig. 5b is the pulse sequence correspondingly, due to the high frequency interferences eliminated, only the pulses were reserved, as a result the pulse sequence is very sparse and easy to be calculated.

Fig. 5c is MP_DCS, the rotation frequency and its double could be seen clearly, and the amplitudes of the other frequencies are all very small, we could easily determine that the bearing is normal. Fig. 6a is WES corresponding to Fig. 5a, d1 to d4 are the envelope spectrum of the first to fourth levels of wavelet high frequency coefficients, a4 is the envelope spectrum of wavelet low frequency coefficient, the rotation frequency and its double could be seen in d2, d3 with many interference and including a peak of characteristic defect frequency of rolling element failure, it is difficult to determine the bearing condition and may cause confusion in a routine diagnostic survey. Fig. 6b is EMDES correspondingly, the rotation frequency or its double could be clearly identified only in IMF1, IMF2, however, there exists many interferential components resulting in difficult to diagnose.

Fig. 5The vibration signal of normal bearing and its decomposition: a) vibration signal of normal bearing, b) the pulse sequence, c) the DCS of pulse sequence

a)

b)

c)

Fig. 6The decomposition corresponding to Fig. 5(a): a) WES, b) EMDES

a)

b)

3.1.2. The analysis of outer race failure bearing

Fig. 7a is the vibration signal of outer race failure, Fig. 7b, Fig. 8a, b are MP_DCS, WES and EMDES respectively of the fault signal. The outer characteristic defect frequency and its harmonics (up to eighteenth harmonic) are all clear and with few noises, these features indicate clearly that the bearing carries an outer race defect. In contrast, although the outer characteristic defect frequency and its harmonics could be recognized in WES or EMDES, they provide less explicit information about the condition of bearing, in addition, the random noise in the spectrum also blurs the featured components.

Fig. 7The vibration signal of the bearing with outer race failure and its MP_DCS: a) vibration signal of outer race, b) MP_DCS

a)

b)

Fig. 8The decomposition corresponding to Fig. 7(a): a) WES, b) EMDES

a)

b)

3.1.3. The analysis of inner race failure bearing

The inner race fault vibration signal is shown in Fig. 9a. Fig. 9b, and Fig. 10a, b are its MP_DCS, WES, and EMDES respectively. In Fig. 9b, the frequency of the peak with the largest amplitude is 158.2 Hz, is approximately the inner race characteristic defect frequency, hence we could determine that the bearing was failure in inner race. For WES, nevertheless there exist peaks of 158.2 Hz in d3 and d4, the amplitude of other frequencies, for example, 58.6 Hz (double of rotation frequency), are relatively high, thus it is a challenge to determine the bearing condition, a similar situation also occurs in EMDES.

Fig. 9The vibration signal of the bearing with inner race failure and its MP_DCS: a) vibration signal of outer race, b) MP_DCS

a)

b)

Fig. 10The decomposition corresponding to Fig. 9(a): a) WES, b) EMDES

a)

b)

3.1.4. The analysis of ball failure bearing

The ball fault vibration signal is shown in Fig. 11(a). Fig. 11(b), and Fig. 12(a), (b) are its MP_DCS, WES, EMDES respectively. In Fig. 11(b), from the peak with the frequency of 137.7 Hz, we can speculate that the bearing was ball failure. For WES, despite there exists ball characteristic defect frequency in d3, but the peak is not prominent enough, and it has been almost submerged in interferential components. For EMDES, it displays as a peak of ball characteristic defect frequency, as well as the inner race characteristic defect frequency, on the other hand it produces more irrelevant components which do not contribute to fault features and can further lead to inaccurate diagnosis.

Fig. 11The vibration signal of the bearing with ball failure and its MP_DCS: a) vibration signal of outer race, b) MP_DCS

a)

b)

Fig. 12The decomposition corresponding to Fig. 11(a): a) WES, b) EMDES

a)

b)