Incipient Fault Diagnosis of Rolling Bearings Based on Impulse-Step Impact Dictionary and Re-Weighted Minimizing Nonconvex Penalty Lq Regular Technique

: The periodical transient impulses caused by localized faults are sensitive and important characteristic information for rotating machinery fault diagnosis. However, it is very difﬁcult to accurately extract transient impulses at the incipient fault stage because the fault impulse features are rather weak and always corrupted by heavy background noise. In this paper, a new transient impulse extraction methodology is proposed based on impulse-step dictionary and re-weighted minimizing nonconvex penalty Lq regular (R-WMNPLq, q = 0.5) for the incipient fault diagnosis of rolling bearings. Prior to the sparse representation, the original vibration signal is preprocessed by the variational mode decomposition (VMD) technique. Due to the physical mechanism of periodic double impacts, including step-like and impulse-like impacts, an impulse-step impact dictionary atom could be designed to match the natural waveform structure of vibration signals. On the other hand, the traditional sparse reconstruction approaches such as orthogonal matching pursuit (OMP), L1-norm regularization treat all vibration signal values equally and thus ignore the fact that the vibration peak value may have more useful information about periodical transient impulses and should be preserved at a larger weight value. Therefore, penalty and smoothing parameters are introduced on the reconstructed model to guarantee the reasonable distribution consistence of peak vibration values. Lastly, the proposed technique is applied to accelerated lifetime testing of rolling bearings, where it achieves a more noticeable and higher diagnostic accuracy compared with OMP, L1-norm regularization and traditional spectral Kurtogram (SK) method.


Introduction
Rolling bearings are extensively used as critical elements in the transmission systems of rotating machinery, and unexpected faults may cause severe mechanical failures and great economic losses or even personal casualties.Therefore, the incipient fault diagnosis (IFD) and health management of rolling bearings are becoming more and more crucial in engineering applications [1][2][3].
Over the past decades, feature extraction techniques were applied to extract the transient characteristics from the original vibration signal, which can indicate the fault location and fault type of rolling bearings.However, fault feature extraction usually suffers from two challenges when the defects occur at the early stage: (1) at the early stage of fault degradation, fault-related components perform incompletely and incipient faults are quite different from obvious failure states, and in practical applications, it seems the normal operating stage is generally considered; (2) the useful weak fault impulse transient characteristics are often submerged in heavy-level background noise and other irrelevant components, which makes the fault diagnosis more difficult.Therefore, the major concern in bearing fault feature extraction is to determine which signal processing tools and algorithms to use to distinguish and diagnose early stage fault characteristics.Up to now, various fault diagnosis techniques have been proposed attempting to address the above challenges, such as wavelet/wavelet-packet transform [4], local mean decomposition (LMD) and its extension [5], minimum entropy deconvolution (MED) and its extension [6,7] and artificial intelligence (AI) algorithms such as artificial neural network (ANN) and fuzzy algorithm [8][9][10], Hilbert envelope spectrum [11], energy and entropy methods [12][13][14], higher order statistical techniques [15][16][17][18], to mention just a few.Unfortunately, some potential drawbacks and severe shortcomings related to the common techniques still remained unresolved.For example, an adaptive decomposition problem exists in the wavelet/wavelet-packet and LMD methods, sample training and fault severity quantitative analysis issues exist in the ANN and fuzzy algorithm, energy and entropy methods are complex and time-consuming, etc.As a consequence, if a weak fault exists, then the acquired vibration signals are rather complex, and these shortcomings may hinder the effectiveness of traditional methods.
Compared to traditional fault feature extraction approaches, sparse representation (SR) as a new signal processing method whereby a given vibration signal can be sparsely represented based on a linear combination of sparse basis or dictionary atoms.The SR has been successfully introduced on fault detection of rotating machinery.For example, Zhu and Fan [19,20] developed an optimal Laplace wavelet, tunable Q-factor wavelet transform, single-side Morlet wavelet basis combined with split variable augmented Lagrangian shrinkage algorithm (SALSA) to extract impulse components and transient features.Du [21] proposed a nuclear norm minimization that uses a weighted low-rank sparse model for bearing fault detection.Cui [22,23] introduced composite dictionary multi-atom matching and a matching pursuit algorithm based on an adaptive impulse dictionary for gear-box and bearing fault diagnosis.The defect-induced impulses redundant dictionary and matching pursuit (MP) approach was proposed by He in [24].Zhang [25] proposed a novel method called kurtosis-based weighted sparse model based on a convex optimization technique; this technique formulated the prior information into a sparse regularization problem and achieved good effect in bearing fault diagnosis.He [26] employed a local time-frequency (TF) domain sparse representation to reconstruct the native pulse waveform structure of fault transients, and proved that the proposed method was superior to traditional the MP and K-singular value decomposition (K-SVD) methods.
Although the sparse representation has achieved successful applications in fault diagnosis of rotating machinery, however, the following two situations still need to be further researched: (1) When a spalling defect or pitting corrosion is induced, a series of successive impulses will be generated during subsequent operation.However, most dictionaries and optimal wavelet-basis constructed in the previous method only use single pulse or single impact frequencies, e.g., the optimal Laplace wavelet, single-side Morlet wavelet basis, transient impulse atoms, etc. Therefore there is no guarantee that the sparse-basis construction can match the natural waveform structure of the vibration signal well.(2) In practice, due to the fluctuation of the load and speed, and the interference of the harsh working environment, some random variations will be generated between an impulse and its neighboring impulses.The traditional sparse reconstruction methods such as greedy pursuit, orthogonal matching pursuit (OMP), L1-norm regularization and iterative shrinkage algorithm ignore those time-varying physical characteristics, which leads to a lower success rate of the transient impulse reconstruction.On the other hand, the traditional sparse reconstruction approaches also treat all vibration signal values equally and thus ignore the fact that the vibration peak value may have more useful information about periodical transient impulses and should be preserved at a larger weight value.In an attempt to overcome the above shortcomings and extend the universality of SR in fault diagnosis, this paper proposes a new sparse representation algorithm called re-weighted minimizing nonconvex penalty Lq regular (R-WMNPLq) combined with an impulse-step impact dictionary for incipient bearing fault diagnosis, using a bearing accelerated life test diagnosis as a case study.Prior to the sparse representation, the original vibration signal is preprocessed by the variational mode decomposition (VMD) technique for incipient fault signal filtering.Furthermore, the impulse-step impact dictionary atoms are constructed to match the natural waveform structure of the vibration signals.Considering the time-varying physical characteristics of transient impulses and keeping a reasonable distribution consistence of peak vibration values, the R-WMNPLq (q = 0.5) technique is employed and the fault frequencies and failure location can be diagnosed accordingly.Incipient transient feature extraction results indicate that the impulse time and the period of transients can be detected more accurately and effectively in cases where previous approaches failed, which can significantly improve the performance of sparse reconstruction for extracting transient impulses from heavy noisy vibration signal.
The remainder of this paper is organized as follows: Section 2 describes the theory of the impulse-step impact dictionary.Section 3 introduces re-weighted minimizing nonconvex penalty Lq regular (q = 0.5) algorithm and the implementation details of this method.In Section 4, the diagnosis results and discussion of the proposed algorithm with other previous approaches are presented.Conclusions are made in Section 5.

Impulse-Step Impact Dictionary and Its Simulation
The periodic transient impulses of rolling bearings are mainly generated by the impact between the bearing elements, inner race and outer race.For example, when there is a pitting failure with a certain size in the outer race, the bearing element initially contacts the anterior fault edge and then exits from the lagging fault edge, as shown in Figure 1a.Thus two impacts are generated due to entry into and exit from the fault region.The first impact approaches a step-like waveform and the second impact resembles an impulse-like waveform.This phenomenon has been proven by the practical fault data from the bearing data center of Case Western Reserve University [27], as shown in Figure 1b.Based on the above analysis, in this paper, a double transient impulse dictionary atom that includes step-like and impulse-like impacts is proposed.The remainder of this paper is organized as follows: Section 2 describes the theory of the impulse-step impact dictionary.Section 3 introduces re-weighted minimizing nonconvex penalty Lq regular (q = 0.5) algorithm and the implementation details of this method.In Section 4, the diagnosis results and discussion of the proposed algorithm with other previous approaches are presented.Conclusions are made in Section 5.

Impulse-Step Impact Dictionary and Its Simulation
The periodic transient impulses of rolling bearings are mainly generated by the impact between the bearing elements, inner race and outer race.For example, when there is a pitting failure with a certain size in the outer race, the bearing element initially contacts the anterior fault edge and then exits from the lagging fault edge, as shown in Figure 1a.Thus two impacts are generated due to entry into and exit from the fault region.The first impact approaches a step-like waveform and the second impact resembles an impulse-like waveform.This phenomenon has been proven by the practical fault data from the bearing data center of Case Western Reserve University [27], as shown in Figure 1b.Based on the above analysis, in this paper, a double transient impulse dictionary atom that includes step-like and impulse-like impacts is proposed.Firstly, the period time of the bearing elements entering and then exiting from the fault region can be calculated by: where the lo is defect size on the outer race, d is the rolling element diameter, D0 is the diameter of outer race, i.e., D0 = Dp + d, Dp is the pitch diameter and fr is the shaft rotation frequency, fc is the bearing cage frequency, i.e., , α is the contact angle.As a matter of fact, when the defect size lo (mm) is smaller than the rolling element diameter d, so the rolling element cannot come Firstly, the period time of the bearing elements entering and then exiting from the fault region can be calculated by: Entropy 2017, 19, 421 where the l o is defect size on the outer race, d is the rolling element diameter, D 0 is the diameter of outer race, i.e., D 0 = D p + d, D p is the pitch diameter and f r is the shaft rotation frequency, f c is the bearing cage frequency, i.e., f c = f r 2 (1 − d D p cos α), α is the contact angle.As a matter of fact, when the defect size l o (mm) is smaller than the rolling element diameter d, so the rolling element cannot come into contact with the bottom of the pitting failure, the distance of the rolling element entering and then exiting from the fault region is half of the defect size l o (mm).Thus the period time ∆t becomes: Similarly, when there is a pitting failure with a certain size l i (mm) on the inner race, the corresponding period time ∆t i can be expressed as follows: which is the same as the period time that was derived for the same defect size on the outer race, as shown in Equation (1).We suppose the moment when the impulse-like impact response occurs is u, thus the step-like impact response occurs is u − ∆t, consequently, the single degree of freedom impulse-like impact (in the form of a decaying sinewave) can be defined as: The single step-like impact can be defined as: Therefore, the impulse-step impact impulse dictionary atom can be defined as: where f n is the system natural frequency, u the time when the impulse-like impact occurs, τ is system damping and a is the peak value ratio of impulse-like to the step-like impact.
In order to generate an impulse-step impact signal representative of that obtained from the double impact of the rolling element with the anterior and lagging fault edge, the time-domain waveforms of impulse-step impact atom, step-impulse impact atom, impulse-step impact atom, and the impulse-step impact signal without/with noise generated using Equation ( 6) are shown in Figure 2, where the simulated bearing type was NACHI 2206GK whose detailed parameters are listed in Table 1.The parameters of the impulse-step impact equation were set as follows: the system damping constant τ is 0.001, peak value ratio a is 0.3, the system natural frequency f n = 10,000 Hz, the impulse-like response happened u is 0.005, the rotor speed rotation frequency f r is 800 rpm.The time-domain waveform of the impulse-like signal with noise is shown in Figure 2e.The signal-noise ratio (SNR) is 20 dB.It can be seen that the similarities between the measured signal and the simulated signal with noise presented in Figure 1b is quite apparent.

Review of Sparse Representation
The basic idea of sparse representation is that a vibration signal can be represented as a linear superposition of a few sparse atoms with residual component [28][29][30][31][32]. Denoting a vibration signal p y R ∈ , the approximating process can be represented as: where x is the approximating signal, n represents residual component, , n , and α k an sparse coefficients of the vibration signal y.As demonstrated in Equation ( 7), there are two issues to be solved: (1) Designing a redundant dictionary D. The first important issue is how to construct a redundancy dictionary D that suitable for the transient behavior of fault impulse components.(2) Recovering sparse coefficients α .Another important issue is how to design an optimization algorithm to calculate the sparse coefficients of vibration fault signal.

Review of Sparse Representation
The basic idea of sparse representation is that a vibration signal can be represented as a linear superposition of a few sparse atoms with residual component [28][29][30][31][32]. Denoting a vibration signal y ∈ R p , the approximating process can be represented as: where x is the approximating signal, n represents residual component, D ∈ R p×n , p << n called redundant dictionary, which consists of n sparse atoms and α k an sparse coefficients of the vibration signal y.As demonstrated in Equation ( 7), there are two issues to be solved: (1) Designing a redundant dictionary D. The first important issue is how to construct a redundancy dictionary D that suitable for the transient behavior of fault impulse components.(2) Recovering sparse coefficients α.Another important issue is how to design an optimization algorithm to calculate the sparse coefficients of vibration fault signal.

Re-Weighted Minimizing Nonconvex Penalty Lq Regular and Its Simulation Experiment
It should be noted that y = Dα + n is a highly underdetermined equation [33], for which there is an infinite set of solutions.In the present method, using an optimization approach, the problem of signal reconstruction by sparse representation under residual error constraints can be calculated by: where c is a threshold of the residual error.Moreover, the prior knowledge of the original signal is usually utilized to regularize the solution under residual error constraint is expressed as where λ is regularization weight and ζ(x) regularization term.From the perspective of Bayesian estimation, the Dα − y 2 2 and ζ(x) can be viewed as the likelihood part and prior knowledge part, respectively.Therefore, the ζ(x) prior knowledge part plays a significant role in signal reconstruction based on sparse representation.
Usually, a fault vibration signal can also be divided into two types: the first is the periodic transient impulse containing step-like impact and impulse-like impact, and the second is the smoothing regions between the impulse and its neighbor impulse, as shown in Figure 2e.For the first one, the physical structure and fault type determine the similarity between two impulses, and the influence of external noise is relatively weak.However, in smoothing regions, due to the fluctuation of the varying load and speed, and the interference of the harsh working environment, the influence of external noise will play a critical role in signal reconstruction.If the noise level is strong, the information of noise in smoothing regions is regarded as structural information in sparse coefficients.Meanwhile, the classical optimization and regularization approaches also treat all vibration signal values equally and thus ignore a fact that the vibration peak value may have more useful information of periodical transient impulses, and cannot remove the false structural information contained in the sparse coefficients, and the traditional methods may cause instability and obvious artifacts in the reconstructed signal.
To overcome this issue, a new re-weighted minimizing nonconvex penalty Lq (0 < q ≤ 1) regular (R-WMNPLq) method is introduced, that is: where q is regular operator, ε(ε > 0) is smoothing parameter, λ(λ > 0) is penalty parameter and b is measurement vector.It should be mentioned that the smoothing parameter plays a critical role in signal reconstruction in terms of smoothing regions.The detailed update procedure is shown in Algorithm 1.
Update α by formula ε k+1 = min ε k , β • r(α (k+1) ) s+1 , where r(α) represents the rearrangement of absolute values of r(α (k+1) ) in the decreasing order, and r(α) s+1 is the (s + 1) th component value of r(α).Note that, if ε k+1 = 0, choose α (k+1) to be an approximation of sparse solution and stop this iteration.11: Let k = k + 1, and return to step 4 to continue.12: Output: Sparse coefficients α; 13: End For the R-WMNPLq method, the following theorem summarizes the results for 0 < q ≤ 1, thus we have the following theorem which can prove the above proposed algorithm: Theorem 1. Error estimation theorem.Suppose that x o is a sparse signal with sparsity level s which satisfies Dx o = b.Without loss of generality, here the sparse coefficient α is substituted by vector x.The smooth parameter ε k → ε * with k → ∞ .Matrix D satisfies the restricted isometry property (RIP) [29,30,32] of order 2 s with δ 2s < 1, when ε * > 0, the sequence {x (k) } has at least one convergent subsequence.Suppose that the limit ε k = ε * is a local optimal solution for Equation (10), we have: where δ s (x ε * ) 2 is the approximate error of x ε * , which satisfies δ s (x ε * ) 2 = inf y 2,0 ≤s x ε * − y 2 .For the special case, when ε * = 0, there must exist a convergent subsequence converging to point x o , it satisfies, where C 1 , C 2 and C 3 are independent positive constants.To prove the theorem 1, the following two lemmas are required.
where C 4 is an independent positive constant.
Herein, combining the above inequalities in Lemmas 1 and 2, Theorem 1 can be proved ultimately.For simplicity, the detailed proof process was derived and presented in the Appendix A. In the next section, the choice of regular operator q will be discussed in detail via a simulation experiment.
For the choice of q (0 < q ≤ 1), we assume q varying among a data scope {0.1, 0.5, 0.7, 1}.Firstly, the dictionary D was generated by a rand-function rand (64, 256), and the test signal has t non-zeros narrow-pulse that subject to the standard Gaussian distribution (SGD), the locations of non-zeros were generated randomly, and the number t varying among {8, 10, 12, . . ., 32}.The penalty parameter λ = 10 −6 is small enough which satisfies Dx = Dx 0 .Taking the R-WMNPLq algorithm iterative 1000 times, if the recovery error (RE) satisfy RE = x r − x 0 2 / x 0 2 ≤ 10 −3 , the algorithm iteration is stopped, where x r stands for non-zeros narrow-pulse.Figure 3a shows the random signal with 32 non-zero pulses and Figure 3b shows the recovery success rate (RSR) curves with different regular operator q.From Figure 3b, it can be seen that q = 0.1, q = 0.5 performed better than q = 0.7 and much better than q = 1.Moreover, the RSR curve with q = 0.5 is slightly higher than q = 0.1.Therefore, in this paper, regular operator q = 0.5 was chosen as the optimal operator.(  , , ) ε λ is monotonically decreasing sequence, hence: for all k ≥ 1 and 1 ≤ i ≤ n, there exists a positive number β which satisfies ( ) , hence: Let , and thus Lemma 2 is proved conclusively.□ Herein, combining the above inequalities in Lemmas 1 and 2, Theorem 1 can be proved ultimately.For simplicity, the detailed proof process was derived and presented in the Appendix section.In the next section, the choice of regular operator q will be discussed in detail via a simulation experiment.For the choice of q (0 < q ≤ 1), we assume q varying among a data scope {0.1, 0.5, 0.7, 1}.Firstly, the dictionary D was generated by a rand-function rand (64, 256), and the test signal has t non-zeros narrow-pulse that subject to the standard Gaussian distribution (SGD), the locations of non-zeros were generated randomly, and the number t varying among {8, 10, 12, …, 32}.The penalty parameter , the algorithm iteration is stopped, where r x stands for non-zeros narrow-pulse.Figure 3a shows the random signal with 32 non-zero pulses and Figure 3b shows the recovery success rate (RSR) curves with different regular operator q.
From Figure 3b, it can be seen that q = 0.1, q = 0.5 performed better than q = 0.7 and much better than q = 1.Moreover, the RSR curve with q = 0.5 is slightly higher than q = 0.1.Therefore, in this paper, regular operator q = 0.5 was chosen as the optimal operator.

Experimental Evaluation
The experimental setup of our roller bearing accelerated life test is shown in

Experimental Evaluation
The experimental setup of our roller bearing accelerated life test is shown in Figure 4. Bearing accelerated vibration signals were generated by an Intelligent Maintenance System (IMS) [35,36].The sampling rate is 20 kHz.The authors analyzed the vibration acceleration data from bearing 1 that the accelerated life test was carried out successively for 8 days (from 12 February 2004 10:32:39 to 19 February 2004 06:22:39, from normal to severe fault, i.e., 9840 min).Meanwhile, four Rexnord ZA-2115 bearings (pitch diameter is 71.501 mm, roller diameter 8.4074 mm, roller number 16 and contact angle 15.17 deg) were detected using acceleration sensors and thermocouples.Therefore, by calculating, the fault characteristic frequency of bearing 1# outer race is 236.4Hz. contact angle 15.17 deg) were detected using acceleration sensors and thermocouples.Therefore, by calculating, the fault characteristic frequency of bearing 1# outer race is 236.4Hz.In order to avoid propagation of damages to the whole experimental platform and for security reasons, the accelerated life test was stopped when the vibration amplitude of the vibration raw signal surpassed 5 m/s −2 .Figure 5a shows the vibration raw signal of the whole life-cycle of bearing 1. Figure 5b shows the Kurtosis curve over the whole life-cycle of bearing 1 and indicates that there is a long time with normal operation in whole life-cycle, but the period of fault from incipient stage to severe stage is relatively short.As shown in Figure 5, there is an obvious transient feature at point 647 in the incipient fault.However, due to the interference of harsh working environment and background noises, the engineer cannot sure whether the fault is happened before point 647 or not.Hence, to verify the effectiveness of the proposed method for bearing incipient fault diagnosis, the experimental data at point 535 was chosen which has no obviously wave phenomenon in whole life-cycle.In order to avoid propagation of damages to the whole experimental platform and for security reasons, the accelerated life test was stopped when the vibration amplitude of the vibration raw signal surpassed 5 m/s −2 .Figure 5a shows the vibration raw signal of the whole life-cycle of bearing 1. Figure 5b shows the Kurtosis curve over the whole life-cycle of bearing 1 and indicates that there is a long time with normal operation in whole life-cycle, but the period of fault from incipient stage to severe stage is relatively short.As shown in Figure 5, there is an obvious transient feature at point 647 in the incipient fault.However, due to the interference of harsh working environment and background noises, the engineer cannot sure whether the fault is happened before point 647 or not.Hence, to verify the effectiveness of the proposed method for bearing incipient fault diagnosis, the experimental data at point 535 was chosen which has no obviously wave phenomenon in whole life-cycle.contact angle 15.17 deg) were detected using acceleration sensors and thermocouples.Therefore, by calculating, the fault characteristic frequency of bearing 1# outer race is 236.4Hz.In order to avoid propagation of damages to the whole experimental platform and for security reasons, the accelerated life test was stopped when the vibration amplitude of the vibration raw signal surpassed 5 m/s −2 .Figure 5a shows the vibration raw signal of the whole life-cycle of bearing 1. Figure 5b shows the Kurtosis curve over the whole life-cycle of bearing 1 and indicates that there is a long time with normal operation in whole life-cycle, but the period of fault from incipient stage to severe stage is relatively short.As shown in Figure 5, there is an obvious transient feature at point 647 in the incipient fault.However, due to the interference of harsh working environment and background noises, the engineer cannot sure whether the fault is happened before point 647 or not.Hence, to verify the effectiveness of the proposed method for bearing incipient fault diagnosis, the experimental data at point 535 was chosen which has no obviously wave phenomenon in whole life-cycle.Figure 6a-d show the original signal at point 535 (1024 sampling points were chosen, i.e., about 0.05 s), the time-frequency distribution a with short-time Fourier transform (STFT), the amplitude spectrum and Hilbert envelope spectrum of original vibration signal, respectively.From Figure 6a, the periodical impulses are submerged in heavy noise and fault type cannot be determined yet.From Figure 6d, although the spectrum peak at 240 Hz and its harmonic frequencies which consists with the outer-race fault frequency can be detected without de-nosing, however, the spectrum peak masked by heavy background noise and features are not be evident enough to detect fault.Considering the complexity of bearing vibration signals with different frequencies and the repetition of fault patterns, the variational mode decomposition (VMD) [37][38][39] method was used to preprocess the original signal.We can obviously observe that the amplitude spectrum structures were displayed more clearly when the modal number is 20, however, the modal duplication phenomenon starts to appear when the modal number reaches 21, which demonstrates that the model number K = 20 is better than K = 21.Figure 8 shows 20 intrinsic mode functions (IMF) models of original vibration signal decomposed by VMD method, which can be help to distinguish the periodic impulse from the mixed noisy signal.Here, the kurtosis of the 19th IMF is the maximum, which means the transient impulses feature may be contained in 19th IMF model based on criterion of maximum kurtosis [39].
Entropy 2017, 19, x 12 of 20 Figure 7a,b show the amplitude spectrum of each IMF mode with K = 20 and K = 21, respectively.We can obviously observe that the amplitude spectrum structures were displayed more clearly when the modal number is 20, however, the modal duplication phenomenon starts to appear when the modal number reaches 21, which demonstrates that the model number K = 20 is better than K = 21.Figure 8 shows 20 intrinsic mode functions (IMF) models of original vibration signal decomposed by VMD method, which can be help to distinguish the periodic impulse from the mixed noisy signal.Here, the of the 19th IMF is the maximum, which means the transient impulses feature may be contained in 19th IMF model based on criterion of maximum kurtosis [39].2) IMF( 3) IMF( 4) IMF( 5) IMF( 6) IMF( 7) IMF( 8) IMF( 9) IMF( 10) IMF( 11) IMF( 12) IMF( 13) IMF( 14) IMF( 15) IMF( 16) IMF( 17) IMF( 18) IMF( 19) IMF( 20 4) IMF( 5) IMF( 6) IMF( 7) IMF( 8) IMF( 9) IMF( 10) IMF( 11) IMF( 12) IMF( 13) IMF( 14) IMF( 15) IMF( 16) IMF( 17) IMF( 18) IMF( 19) IMF( 20 Furthermore, the proposed R-WMNPLq (q = 0.5) algorithm was employed on the 19th IMF model and its related parameters are illustrated in Table 3.The reconstructed periodic impulses signal, the time-frequency distribution with short-time Fourier transform (STFT) and the Hilbert envelope spectrum of the reconstructed periodic impulses signal are depicted in Figure 9a-c respectively.It can be observed that the proposed R-WMNPLq (q = 0.5) algorithm combined with the impulse-step impact dictionary not only extracts the transient impulse components clearly but also the noise components in reconstructed vibration signal have been removed evidently, and the signalto-noise ratio (SNR) in Figure 9a is −27.7460db.Compared with the original vibration signal as shown in Figure 6, the time-frequency distribution combines the information in time and frequency domains, it can be easily seen from Figure 9a,b that there are no interference noises among the extracted transient impulses.The transient impulses time-frequency distribution is more clearly, which effectively reveal the fault feature from the incipient vibration signal.The Hilbert envelop spectrum is shown in Figure 9c.As can be seen, the characteristic frequency 240 Hz (close to the theoretical fault frequency of outer race 236.4 Hz) and its harmonic frequencies (3fo, 4fo, 5fo, 6fo, 7fo and 8fo) are clearly detected, therefore, the proposed method is exactly suitable for incipient fault bearing signal.
Table 3. Parameters of the proposed R-WMNPLq (q = 0.5) method.Furthermore, the proposed R-WMNPLq (q = 0.5) algorithm was employed on the 19th IMF model and its related parameters are illustrated in Table 2.The reconstructed periodic impulses signal, the time-frequency distribution with short-time Fourier transform (STFT) and the Hilbert envelope spectrum of the reconstructed periodic impulses signal are depicted in Figure 9a-c respectively.It can be observed that the proposed R-WMNPLq (q = 0.5) algorithm combined with the impulse-step impact dictionary not only extracts the transient impulse components clearly but also the noise components in reconstructed vibration signal have been removed evidently, and the signal-to-noise ratio (SNR) in Figure 9a is −27.7460db.Compared with the original vibration signal as shown in Figure 6, the time-frequency distribution combines the information in time and frequency domains, it can be easily seen from Figure 9a,b that there are no interference noises among the extracted transient impulses.The transient impulses time-frequency distribution is more clearly, which effectively reveal the fault feature from the incipient vibration signal.The Hilbert envelop spectrum is shown in Figure 9c.As can be seen, the characteristic frequency 240 Hz (close to the theoretical fault frequency of outer race 236.4 Hz) and its harmonic frequencies (3f o , 4f o and 5f o ) are clearly detected, therefore, the proposed method is exactly suitable for incipient fault bearing signal.

Regular
Table 2. Parameters of the proposed R-WMNPLq (q = 0.5) method.A considerable amount of literature has been published on the application of orthogonal matching pursuit (OMP) and L1-norm regularization algorithms in mechanical fault diagnosis [40][41][42].To further validate the superiority of the proposed method, the OMP and L1-norm regularization techniques were sequentially used on the 19th IMF model vibration signal.The running iteration time is set as 50.Figures 10a,c,e and 10b,d,f are the results of the OMP and L1-norm regularization methods, respectively.The signal-to-noise ratio (SNR) in Figure 10a,b is −29.1315db and −35.4638 db, respectively.As shown in Figure 10a, strong noises still remained in the reconstructed impulsive signal.Besides, by comparing Figure 6a with Figure 10a, it should be noted that the conventional L1norm technique removes too much energy of the original vibration signal to effectively reduce the noises but also shrinks the fault feature frequency.Compared with the results of Hilbert envelope spectrum shown in Figure 10e,f, it can be seen that the OMP and L1-norm regularization technique do not have a satisfactory performance in incipient fault diagnosis of rolling bearings.A considerable amount of literature has been published on the application of orthogonal matching pursuit (OMP) and L1-norm regularization algorithms in mechanical fault diagnosis [40][41][42].To further validate the superiority of the proposed method, the OMP and L1-norm regularization techniques were sequentially used on the 19th IMF model vibration signal.The running iteration time is set as 50. Figure 10a,c,e and Figure 10b,d,f are the results of the OMP and L1-norm regularization methods, respectively.The signal-to-noise ratio (SNR) in Figure 10a,b is −29.1315db and −35.4638 db, respectively.As shown in Figure 10a, strong noises still remained in the reconstructed impulsive signal.Besides, by comparing Figure 6a with Figure 10a, it should be noted that the conventional L1-norm technique removes too much energy of the original vibration signal to effectively reduce the noises but also shrinks the fault feature frequency.Compared with the results of Hilbert envelope spectrum shown in Figure 10e,f, it can be seen that the OMP and L1-norm regularization technique do not have a satisfactory performance in incipient fault diagnosis of rolling bearings.In addition to the sparse representation method methods, the superiority and effectiveness requires further validation with a traditional approach.Therefore, the same signal, namely the 19th IMF model, was also processed by the spectral Kurtogram (SK) method [43,44].The Kurtogram of the 19th IMF model is displayed in Figure 11a, from which the optimal demodulation frequency band, namely 5333-10,000 Hz, can be detected.Thus a band-pass filter was designed to extract the potential features from the 19th IMF model, then the envelope spectrum was applied to the filtered signal, the corresponding envelope spectrum is shown in Figure 11b.As can be seen, no explicit fault characteristic frequencies are observed and it is also hard to distinguish the fault location from the incipient fault signal.In addition to the sparse representation method methods, the superiority and effectiveness requires further validation with a traditional approach.Therefore, the same signal, namely the 19th IMF model, was also processed by the spectral Kurtogram (SK) method [43,44].The Kurtogram of the 19th IMF model is displayed in Figure 11a, from which the optimal demodulation frequency band, namely 5333-10,000 Hz, can be detected.Thus a band-pass filter was designed to extract the potential features from the 19th IMF model, then the envelope spectrum was applied to the filtered signal, the corresponding envelope spectrum is shown in Figure 11b.As can be seen, no explicit fault characteristic frequencies are observed and it is also hard to distinguish the fault location from the incipient fault signal.

Conclusions
This work originated from a study on the sparse representation approach and incipient fault diagnosis of rolling bearings.Although a lot of works have achieved successful application in fault diagnosis of rotating machinery based upon sparse representation methods such as greedy pursuit, orthogonal matching pursuit (OMP), L1-norm regularization, the developed approaches are not satisfactory for reconstructing periodic transient impulses and identifying the physical structure information of periodic impulses, especially when the fault is in an incipient stage.
This paper proposes a novel feature extraction method for incipient bearing fault diagnosis combining re-weighted minimizing nonconvex penalty Lq (R-WMNPLq, q = 0.5) regular and impulse-step impact dictionary.The proposed method provides a new point of view for periodic instantaneous impulse reconstruction by introducing a penalty parameter, smoothing parameter and regular operator on a sparse representation model to guarantee the reasonable distribution consistence of peak vibration values.On the other hand, the original physical structure information was formed by impulse-step impact dictionary atoms.Effectiveness in the extraction of transient impulse is verified by an accelerated life test.The experimental analysis shows that the proposed method can achieve good performance in reducing noises and extracting fault characteristic from raw vibration signals in comparison with the matching pursuit method (OMP), L1-norm regularization and spectral Kurtogram (SK) method, especially for vibration signals with heavy background noises, and it is well suited for on-line practical applications.
However, the proposed methodology is only applicable for accelerated life test of rolling bearings under constant operating conditions, and variable conditions such as variable speed, torque and variable harsh working environments should be considered in the future which may help generalizing the proposed method.Moreover, the proposed sparse representation method can be also improved for the detection of multiple faults concurrence in our future work.

Conclusions
This work originated from a study on the sparse representation approach and incipient fault diagnosis of rolling bearings.Although a lot of works have achieved successful application in fault diagnosis of rotating machinery based upon sparse representation methods such as greedy pursuit, orthogonal matching pursuit (OMP), L1-norm regularization, the developed approaches are not satisfactory for reconstructing periodic transient impulses and identifying the physical structure information of periodic impulses, especially when the fault is in an incipient stage.
This paper proposes a novel feature extraction method for incipient bearing fault diagnosis combining re-weighted minimizing nonconvex penalty Lq (R-WMNPLq, q = 0.5) regular and impulse-step impact dictionary.The proposed method provides a new point of view for periodic instantaneous impulse reconstruction by introducing a penalty parameter, smoothing parameter and regular operator on a sparse representation model to guarantee the reasonable distribution consistence of peak vibration values.On the other hand, the original physical structure information was formed by impulse-step impact dictionary atoms.Effectiveness in the extraction of transient impulse is verified by an accelerated life test.The experimental analysis shows that the proposed method can achieve good performance in reducing noises and extracting fault characteristic from raw vibration signals in comparison with the matching pursuit method (OMP), L1-norm regularization and spectral Kurtogram (SK) method, especially for vibration signals with heavy background noises, and it is well suited for on-line practical applications.
However, the proposed methodology is only applicable for accelerated life test of rolling bearings under constant operating conditions, and variable conditions such as variable speed, torque and variable harsh working environments should be considered in the future which may help generalizing the proposed method.Moreover, the proposed sparse representation method can be also improved for the detection of multiple faults concurrence in our future work.

Figure 1 .
Figure 1.The time-domain waveform of the fault signal for a single pitting failure.(a) The physical model; (b) time-domain waveform of the fault signal.

Figure 1 .
Figure 1.The time-domain waveform of the fault signal for a single pitting failure.(a) The physical model; (b) Time-domain waveform of the fault signal.

Figure 2 .
Figure 2. The time-domain waveform of (a) impulse-like impact atom; (b) step-impulse impact atom; (c) impulse-step impact atom; (d) impulse-step impact signal without noise and; (e) impulse-step impact signal with a SNR of 20 dB.

Figure 2 .
Figure 2. The time-domain waveform of (a) impulse-like impact atom; (b) step-impulse impact atom; (c) impulse-step impact atom; (d) impulse-step impact signal without noise and; (e) impulse-step impact signal with a SNR of 20 dB.

Figure 3 .
Figure 3.Comparison results of recoverability with different q.(a) Random signal with 32 non-zero pulses; (b) Comparison results of RSR with different q.

=
Dx Dx .Taking the R-WMNPLq algorithm iterative 1000 times, if the recovery error (RE) satisfy

Figure 4 .Figure 3 .
Figure 3.Comparison results of recoverability with different q.(a) Random signal with 32 non-zero pulses; (b) Comparison results of RSR with different q.

Figure 5 .
Figure 5.The vibration raw signal and the Kurtosis curve of the whole life-cycle of bearing 1.(a) The vibration raw signal of the whole life-cycle of bearing 1; (b) The Kurtosis curve of the wholse life-cycle of bearing 1.

Figure 5 .Figure 5 .
Figure 5.The vibration raw signal and the Kurtosis curve of the whole life-cycle of bearing 1.(a) The vibration raw signal of the whole life-cycle of bearing 1; (b) The Kurtosis curve of the wholse life-cycle of bearing 1.

Figure 6 .
Figure 6.Original vibration signal and its time-frequency analysis.(a) Original vibration signal; (b) Time-frequency distribution of original vibration signal; (c) Amplitude spectrum of the original vibration signal; (d) Hilbert envelope spectrum of original vibration signal.

Figure
Figure 7a,b show the amplitude spectrum of each IMF mode with K = 20 and K = 21, respectively.We can obviously observe that the amplitude spectrum structures were displayed more clearly when the modal number is 20, however, the modal duplication phenomenon starts to appear when the modal number reaches 21, which demonstrates that the model number K = 20 is better than K = 21.Figure8shows 20 intrinsic mode functions (IMF) models of original vibration signal decomposed by VMD method, which can be help to distinguish the periodic impulse from the mixed noisy signal.Here, the kurtosis of the 19th IMF is the maximum, which means the transient impulses feature may be contained in 19th IMF model based on criterion of maximum kurtosis[39].

Figure 7 .
Figure 7.The comparison of amplitude spectrum of the IMF modes.(a) The amplitude spectrum of IMF modes with K = 20 and α = 2000; (b) The amplitude spectrum of IMF modes with K = 21 and α = 2000.

Figure 7 .Figure 8 .
Figure 7.The comparison of amplitude spectrum of the IMF modes.(a) The amplitude spectrum of IMF modes with K = 20 and α = 2000; (b) The amplitude spectrum of IMF modes with K = 21 and α = 2000.

Figure 9 .
Figure 9.The identified results using the proposed method.(a) The reconstructed signal; (b) Timefrequency distribution of the reconstructed signal; (c) Hilbert envelope spectrum of the reconstructed signal.

Figure 9 .
Figure 9.The identified results using the proposed method.(a) The reconstructed signal; (b) Time-frequency distribution of the reconstructed signal; (c) Hilbert envelope spectrum of the reconstructed signal.

Figure 10 .
Figure10.The identified results using the proposed method.(a) The reconstructed signal using the OMP method; (b) Time-frequency distribution of the reconstructed signal using the OMP method; (c) Hilbert envelope spectrum of the reconstructed signal using the OMP method; (d) The reconstructed signal using the L1-Norm regularization method; (e) Time-frequency distribution of the reconstructed signal using the L1-Norm regularization method; (f) Hilbert envelope spectrum of the reconstructed signal using the L1-Norm regularization method.

2 Figure 10 .
Figure10.The identified results using the proposed method.(a) The reconstructed signal using the OMP method; (b) Time-frequency distribution of the reconstructed signal using the OMP method; (c) Hilbert envelope spectrum of the reconstructed signal using the OMP method; (d) The reconstructed signal using the L1-Norm regularization method; (e) Time-frequency distribution of the reconstructed signal using the L1-Norm regularization method; (f) Hilbert envelope spectrum of the reconstructed signal using the L1-Norm regularization method.

Figure 11 .
Figure 11.Diagnosis result using the spectral kurtogram method.(a) Kurtogram of 19th IMF model component; (b) The Hilbert envelope spectrum of the band-pass filtered signal.

Figure 11 .
Figure 11.Diagnosis result using the spectral kurtogram method.(a) Kurtogram of 19th IMF model component; (b) The Hilbert envelope spectrum of the band-pass filtered signal.
Entropy 2017, 19, 421 3 of 20 incipient bearing fault diagnosis, using a bearing accelerated life test diagnosis as a case study.Prior to the sparse representation, the original vibration signal is preprocessed by the variational mode decomposition (VMD) technique for incipient fault signal filtering.Furthermore, the impulse-step impact dictionary atoms are constructed to match the natural waveform structure of the vibration signals.Considering the time-varying physical characteristics of transient impulses and keeping a reasonable distribution consistence of peak vibration values, the R-WMNPLq (q = 0.5) technique is employed and the fault frequencies and failure location can be diagnosed accordingly.Incipient transient feature extraction results indicate that the impulse time and the period of transients can be detected more accurately and effectively in cases where previous approaches failed, which can significantly improve the performance of sparse reconstruction for extracting transient impulses from heavy noisy vibration signal.

Table 1 .
The parameters of the bearing NACHI 2206GK.