Incipient Fault Feature Extraction for Rotating Machinery Based on Improved AR-Minimum Entropy Deconvolution Combined with Variational Mode Decomposition Approach

Qing Li 1,* , Xia Ji 1 and Steven Y. Liang 1,2 1 College of Mechanical Engineering, Donghua University, Shanghai 201620, China; jixia0206@163.com 2 George W. Woodruff School of Mechanical Engineering, Georgia Institute of Technology, Atlanta, GA 30332-0405, USA; steven.liang@me.gatech.edu * Correspondence: suesliqing@163.com or 1159155@mail.dhu.edu.cn; Tel.: +86-021-6779-2583; Fax: +86-021-6779-2439


Introduction
Rolling bearings and gearboxes are widely used in rotating mechanical drive systems in manufacturing.Rotating machinery always endure heavy loads under complicated and adverse operating conditions, and as mechanical systems become increasingly complex and coupled, some slight faults such as pitting, spalling and breakage fault" may lead to a grave danger or industrial catastrophe [1][2][3].Therefore, the high-efficiency and accurate fault diagnosis and condition-based maintenance technique are important in both reducing operating and maintenance costs and improving plant safety.
Recently, the vibration signal processing techniques for fault characteristic extraction have attracted considerable attention and also have been proved to be effective in fault monitoring of roller their center frequencies can be estimated along with the decomposition online.Inspired by the VMD technique, this research utilizes VMD to detect the multiple signatures caused by gearbox multiple-faults.To this end, this paper proposes a novel hybrid method based on improved AR-MED and variational mode decomposition for single and multiple weak faults diagnosis of rotating machinery.Results of experimental and engineering application indicate that the periodic transient impulses can be detected more accurately and effectively in cases where previous approaches failed, which shows a significant improvement in early incipient fault detection of rotating machinery.
The remainder of the paper is organized as follows: in Section 2, the proposed improved AR-MED method is given.Section 3 gives some fundamental definitions, computation algorithm of VMD method and an example discussion.Section 4 validates the effectiveness of the proposed approach by simulated fault diagnosis of rotating machinery.The experimental validation and engineering application validation are presented in Section 5. Finally, Section 6 concludes this paper.

Improved Autoregressive Minimum Entropy Deconvolution (Improved-AR-MED)
When a pitting corrosion or spalling defect on one bearing surface strikes another contact surface, a series of successive impulses and structure resonance sequences will be generated during consecutive operation, regardless of the localized failure on the inner race, outer race or rolling elements.Generally, the fault vibration signal Sig(n) collected from the mechanical equipment can be described as a convolution between the response function of the systems transmission path h(n) and the combined structure resonance signals d(n) (deterministic signals) for the system dynamics, a fault impulse function x(n) and the background noise e(n) (illustrated in Equation ( 1)): where symbol * represents convolution.The objective of many residual signal methods, including the autoregressive (AR) model-based method, is to remove the deterministic parts from the original vibration signal.

Preprocessing: AR Model-Based Method
For a stationary random sequence x n {n = 1, 2, . . ., N}, the p-th order AR(p) model is defined as follows [26]: where, ϕ i is the AR model coefficients and ε n is a residual error.The calculation process of residual signal ε n obtained from the AR model-based is as follows: Step 1: AR model order p.The optimal order of AR model can be calculated by the Akaike information criterion (AIC) [27].
Step 3: Residual error ε n .Perform step1 ahead prediction on the potentially faulty vibration data and calculate the residual prediction error.
As a consequence, by subtracting the predicted parts from the original vibration signal, the residual signal can be obtained including the impulses and stationary background noise.Therefore, the fault impulses parts will be left in the residual ε n of the AR model-based equation.

Improved AR Minimum Entropy Deconvolution
Minimum entropy deconvolution (MED) was originally proposed by Wiggins in [30] and has shown its effectiveness to enhance the impulse characteristic from a mixture vibration fault signal.
Entropy 2017, 19, 317 4 of 26 The fundamental purpose of the MED is to find an inverse filter that counteracts the influence of transmission path and to search an optimum set of filter coefficients that recover the output signal which is closer to the original impulses with the maximum value of kurtosis.
Figure 1 illustrates the basic deconvolution process performed by the AR-MED filtering.This process eliminates the structure resonance signal (deterministic signal).The residual ε n obtained from AR model-based can be used as input signal to the MED deconvolution process.The AR-MED deconvolution filter g(n) generates output signals ∧ y(n) which represents the fault impulses as close as possible to the original input x(n).
Entropy 2017, 19, 317 4 of 26 transmission path and to search an optimum set of filter coefficients that recover the output signal which is closer to the original impulses with the maximum value of kurtosis.Figure 1 illustrates the basic deconvolution process performed by the AR-MED filtering.This process eliminates the structure resonance signal (deterministic signal).The residual n ε obtained from AR model-based can be used as input signal to the MED deconvolution process.The AR-MED deconvolution filter g(n) generates output signals ( ) y n ∧ which represents the fault impulses as close as possible to the original input x(n).

( ) ( ) ( ) ( ) ( ) n y n x n h n e n
x n e n d n h n = + + * The objective of the MED deconvolution is to find an optimal inverse filter g(n) to recover the periodic impulse signal, that is: where K represents the length of the inverse filter g(n).
By taking the partial derivative with respect to l for Equation (4), we have: During the process of periodic impulse recovery x(n), the objective of inverse filter g(n) is to recover the main characteristics information and to maximize the kurtosis of x(n).The kurtosis is taken as the normalized 4-th order moment given by: Optimum filter settings occur when the derivative of the Equation ( 6) is zero such as: According to Equation (5), by using ( ) and thus the Equation ( 7) can be The objective of the MED deconvolution is to find an optimal inverse filter g(n) to recover the periodic impulse signal, that is: where K represents the length of the inverse filter g(n).
By taking the partial derivative with respect to l for Equation (4), we have: During the process of periodic impulse recovery x(n), the objective of inverse filter g(n) is to recover the main characteristics information and to maximize the kurtosis of x(n).The kurtosis is taken as the normalized 4-th order moment given by: Entropy 2017, 19, 317 5 of 26 Optimum filter settings occur when the derivative of the Equation ( 6) is zero such as: According to Equation (5), by using ∂g(n)/∂g(l) = y(n − l) and thus the Equation ( 7) can be rewritten as: Equation ( 8) can be written in matrix form: The vector f is obtained from the cross correlation of input and output signal ε(n) and ∧ y(n) of the inverse filter g(n), A is the Toeplitz autocorrelation matrix of the ε(n) with dimension (L × L) and g is the optimum filter coefficient.Here, f (l) can be expressed as: where, a = Finally, the inverse filter g can be obtained as: The iterative steps of the MED are tabulated in below: (1) Assume initial filter g(n) (0) as a centered impulse, i.e., g(n ( (4) Solve for new filter coefficients, i.e., g (i) = A −1 f (i) .
(5) Compute the error criterion, Err = g (i) − g (i−1) 2 2 , if Err ≤ tolerance, the iteration process is finished, or Err > tolerance, then i = i + 1 and repeat the process from step 2.
However, there are a series of limitations (such as optimal solution problem) with existing AR-MED in rotating machinery fault diagnosis, the topmost being that the AR-MED approach only deconvolves a single impulse or a finite sequence of impulses, while in the rotating machinery fault diagnosis, the purpose is to solve a transient impulse sequence to be better suited for fault detection.In the next section, an improved algorithm (namely, multiple impulses AR-MED approach) aims at overcoming the defect of a single impulse or a finite sequence of impulses problem of traditional AR-MED is given.

Definition.
Multiple points D-Norm, see reference [32]: Entropy 2017, 19, 317 6 of 26 where x = g × y and The multiple impulses D-Norm (MID) maximization problem can be expressed as: where parameter t is a constant vector that represents the location of multiple impulses.For example, if t = [0100010100] T , it means that the location of output impulses at t = 2, t = 6 and t = 8 will be deconvolved.Since g is the maximum of t T |x k |/ x , if and only if, g is an extremal point of t T |x k |/ x , the equivalent of Equation ( 13) can be established.
Hence, the maxima and minima are found by differentiating with respect to the filter g and equating to 0, thus we have: x k y k+L−l , therefore, the Equation ( 13) can be written as: Let , for all l = 1, 2, ..., L, and converting the Equation (15) to matrix form, then: where According to Equation (4), i.e., x = Y 0 T g, Equation ( 16) can be expressed as: If g is a solution of Equation (17), then any multiple is also a solution, namely, g = c ∧ g, we have: Entropy 2017, 19, 317 7 of 26 Therefore multiples of are non-trivial solutions for g, we have the solution as follows: Finally, for all k = 1, 2, ..., N − L, the array of possible solutions g, namely, G = [G 1 , G 2 , ..., G N−L ], which can be calculated as G = (Y 0 Y 0 T) −1 Y 0 t.Furthermore, the array of outputs X = [X 1 , X 2 , ..., X N−L ] can be obtained, i.e., X In the improved AR minimum entropy deconvolution method, the multiple transient impulses can be represented by the output solution X k .

Review of VMD
The purpose of the variational mode decomposition (VMD) is to decompose a multi-component original signal into k discrete number of band-limited Intrinsic Mode Functions (IMFs) such that each IMF is considered to compact around their respective center frequencies and its bandwidth can be estimated using the H 1 Gaussian smoothness of the demodulated signal [25].Unlike EMD and LMD methods, the number of IMF components to be involved in the decomposition is chosen prior to the decomposition processing.Here, the IMF component u k (t) can be expressed as follows: where,A k (t) is the instantaneous amplitude, the phase φ k (t) is a non-decreasing function and the ω k (t) is instantaneous phase, i.e., ω k (t) = φ k (t) = dφ k (t)/dt.In addition, contrary to the EMD and LMD methods, the IMF component generated by VMD is starting from the coarsest level, that is to say, the low IMFs represents fast oscillations (i.e., high frequency modes), and high order IMFs represents slow oscillations (i.e., low frequency modes).The bandwidth of each IMF component u k (t) can be evaluated with following steps: Step 1: For each IMF component u k (t), the real signal is converted to the analytic signal by means of the Hilbert transform to obtain a unilateral frequency spectrum as follows: Step 2: For each IMF component u k (t), the frequency spectrum of u k (t) is translated to baseband by using modulation property to the respective estimated center frequencies.That is: Step 3: The bandwidth of each IMF component u k (t) is estimated through the H 1 Gaussian smoothness of the demodulated signal, i.e., the squared L 2 -norm of the gradient.So the resulting is written as a constrained variational problem as follows: where, f is the original signal, {u k } represents the K-th IMF model with centre frequency {ω k }, δ(t) is the Dirac function, k is the number of IMF models, t is time script and symbol * denotes convolution, respectively.To address the variational problem above, a quadratic penalty α and Lagrangian multipliers λ(t) are introduced, thus the augmented Lagrangian L is given by: Then, Equation ( 24) can be solved by the alternate direction method of multipliers method (ADMM), the detailed ADMM algorithms of which are given in [33][34][35].Finally, all the IMF models obtained from this solution are written as: The iterative steps of the VMD are as follows: ( (3) For all ω ≥ 0, update ∧ u k is as follows: (4) Update ∧ ω k is as follows: (5) Update λ is as follows: (6) Iterate steps (2) to (5) until the residual satisfies a given stopping criterion.That is: Entropy 2017, 19, 317 9 of 26 In the end, the k IMF modes are obtained by the inverse Fast Fourier transform (IFFT) of the signal

Example Discussion
In order to compare the decomposition performance of the EEMD, LMD and VMD methods, the IMFs, PFs and Hilbert spectrum plots are discussed.Here, a multi-component AM-FM signal x(t) is given as follows: where sampling point N = 1000, t = 1/1000:1/1000:1, sampling frequency Fs = 1024 Hz.
Note that the number of modes k is adaptively obtained and unchangeable in EEMD and LMD methods, however, the number of IMFs is required to predetermine for VMD method, thus the number of IMFs k in VMD was set to be 3 for the further decomposition.The time domain signal waveforms of x(t), x 1 (t), x 2 (t) and x 3 (t) are shown in Figure 2. Examples of IMFs and Product functions (PF) modes obtained respectively by the EEMD, LMD and VMD are illustrated in Figure 3a,c,e.The Hilbert spectrum plots of the mono-components derived from the three decomposition methods are shown in Figure 3b,d,f, respectively.It can be found from Figure 3a-d that the IMFs and PFs derived from EEMD and LMD method are all anamorphic with serve scale-mixing problem.While in the decomposition of VMD method (see Figure 3e,f), it can be clearly found that three IMFs components, i.e., IMF1(t), IMF2(t) and IMF3(t) are corresponding to the x 1 (t), x 2 (t) and x 3 (t) of the original AM-FM signal x(t), respectively.Moreover, the three horizontal lines (instantaneous frequency are 4 Hz, 14 HZ and 30 Hz) in Figure 3f has a narrower bandwidth than that of the EEMD and LMD methods.Therefore, the comparison results illustrate that the IMF components can be accurately extracted from the mixed signal by using the VMD method, which is much closer approximation to the real components of x(t).

Example Discussion
In order to compare the decomposition performance of the EEMD, LMD and VMD methods, the IMFs, PFs and Hilbert spectrum plots are discussed.Here, a multi-component AM-FM signal x (t) is given as follows: where sampling point N = 1000, t = 1/1000:1/1000:1, sampling frequency Fs = 1024 Hz.
Note that the number of modes k is adaptively obtained and unchangeable in EEMD and LMD methods, however, the number of IMFs is required to predetermine for VMD method, thus the number of IMFs k in VMD was set to be 3 for the further decomposition.The time domain signal waveforms of x(t), x1(t), x2(t) and x3(t) are shown in Figure 2. Examples of IMFs and Product functions (PF) modes obtained respectively by the EEMD, LMD and VMD are illustrated in Figure 3a,c,e.The Hilbert spectrum plots of the mono-components derived from the three decomposition methods are shown in Figure 3b,d,f, respectively.It can be found from Figure 3a-d that the IMFs and PFs derived from EEMD and LMD method are all anamorphic with serve scale-mixing problem.While in the decomposition of VMD method (see Figure 3e,f), it can be clearly found that three IMFs components, i.e., IMF1(t), IMF2(t) and IMF3(t) are corresponding to the x1(t), x2(t) and x3(t) of the original AM-FM signal x(t), respectively.Moreover, the three horizontal lines (instantaneous frequency are 4 Hz, 14 HZ and 30 Hz) in Figure 3f has a narrower bandwidth than that of the EEMD and LMD methods.Therefore, the comparison results illustrate that the IMF components can be accurately extracted from the mixed signal by using the VMD method, which is much closer approximation to the real components of x(t).

Simulation Analysis for Weak Fault Extraction of Rotating Machinery
In order to verify the effectiveness of the above algorithms, in this section, a bearing simulation fault signal can be described by the following mathematical model: in which A 0 = 1 denotes the amplitude of modulating signal, a = 0.1 is a damping coefficient, sampling point N = 8192, t = 1/2000:1/2000:0.4096,sampling frequency Fs = 20,000 Hz, the frequencies of the above signal are as follows: f 1 = 20, f 2 = 35, and f 3 = 2000 represents the natural frequency of excited structure.Moreover, a stationary white noise x 3 (t) was added to obtain the simulated signals with SNR = −6dB.
Figure 4a depicts a periodic impulse signal and Figure 4b depicts the raw simulation fault signal.It is clear to see that the periodic impulse would be easily lost or submerged within the heavy background noise.Figure 4c depicts the envelope spectrum of the raw simulation fault signal.It should be noted that the fault characteristic frequency f = 100 Hz, which convolves the resonance frequency band that cannot be detected.

Simulation Analysis for Weak Fault Extraction of Rotating Machinery
In order to verify the effectiveness of the above algorithms, in this section, a bearing simulation fault signal can be described by the following mathematical model: Figure 4a depicts a periodic impulse signal and Figure 4b depicts the raw simulation fault signal.It is clear to see that the periodic impulse would be easily lost or submerged within the heavy background noise.Figure 4c depicts the envelope spectrum of the raw simulation fault signal.It should be noted that the fault characteristic frequency f = 100 Hz, which convolves the resonance frequency band that cannot be detected.Hence, in order to simplify the appearance of a signal and obtain the periodic impulses from the noisy mixed signal, the AR-MED filter method was introduced directly.Figure 5a shows the AR-MED filtered signal waveform and Figure 5b   Hence, in order to simplify the appearance of a signal and obtain the periodic impulses from the noisy mixed signal, the AR-MED filter method was introduced directly.Figure 5a shows the AR-MED filtered signal waveform and Figure 5b illustrates the envelope spectrum of AR-MED filtered signal, respectively.It can be seen in Figure 5b that the fault characteristic frequency is also lost in the background noise.This phenomenon can be explained as follows: (1) AR-MED is optimizing the norm Kurtosis which prefers a solution of a single impulse or deterministic impulse.For the bearing fault signals, a series of periodic impulses features maybe not suit for this approach.(2) There are still a heavy noise remained in the output signal that generated by AR-MED filter, which means the impulsive convolved in the suppressed resonance frequency bands cannot be detected by AR-MED.
Figure 5c,d respectively illustrate the filtered signal waveform and its envelope spectrum with improved AR minimum entropy deconvolution method (here, the length of filter-size L = 1000).Comparing the filtering result and envelope spectrum obtained from improved AR minimum entropy deconvolution, it can be observed that multiple strong peaks were deconvolved and the fault characteristic frequency and its harmonics can be clearly observed.However, under highly noisy background the spectrum magnitude in Figure 5d is weak, and there are many interference components exist on both left and right of fault characteristic frequencies.
Entropy 2017, 19, 317 12 of 26 signal, respectively.It can be seen in Figure 5b that the fault characteristic frequency is also lost in the background noise.This phenomenon can be explained as follows: (1) AR-MED is optimizing the norm Kurtosis which prefers a solution of a single impulse or deterministic impulse.For the bearing fault signals, a series of periodic impulses features maybe not suit for this approach.(2) There are still a heavy noise remained in the output signal that generated by AR-MED filter, which means the impulsive convolved in the suppressed resonance frequency bands cannot be detected by AR-MED.
Figure 5c,d respectively illustrate the filtered signal waveform and its envelope spectrum with improved AR minimum entropy deconvolution method (here, the length of filter-size L = 1000).Comparing the filtering result and envelope spectrum obtained from improved AR minimum entropy deconvolution, it can be observed that multiple strong peaks were deconvolved and the fault characteristic frequency and its harmonics can be clearly observed.However, under highly noisy background the spectrum magnitude in Figure 5d is weak, and there are many interference components exist on both left and right of fault characteristic frequencies.Consequently, the VMD approach was employed to decompose the filtered signal generated by improved AR-MED method.Since the VMD requires a predetermined number of sub-signals K, for instance, if the number of IMFs is M, then the parameter M used to determine the number of Entropy 2017, 19, 317 13 of 26 decomposition result, therefore, the resulting central frequency corresponding to different parameter K from 2 to 7 are tabulated in Table 1.As shown in Table 1, it can be seen that the K varying from 2 to 6 meets the condition that all center frequencies are better separated by the uniform distribution, while K = 7 shows that the central frequency 1699 Hz is very close to the 2100 Hz, and the natural frequency of excited structure (2000 Hz from Equation ( 31)) disappeared, thus the best result is given by VMD with K = 6.The decomposed results using VMD and their magnitude spectrum have been shown in Figure 6a,b.It is observed that all the signatures are successfully identified and clearly shown in 1st to 6-th IMF model.Here, according to the criterion of maximum kurtosis, i.e., the maximum kurtosis criterion component contains the main fault information of the mixed signal.The kurtosis is calculated to select the best fault sub-band for the envelope analysis, and the kurtosis of each IMF mode is as follows: K IMF1 = 2.2802, K IMF2 = 6.7916,K IMF3 = 3.2967, K IMF4 = 3.0748, K IMF5 = 3.0472 and K IMF6 = 3.1231.Hence, the IMF2 component was elected to extract the fault characteristic frequency.
Consequently, the VMD approach was employed to decompose the filtered signal generated by improved AR-MED method.Since the VMD requires a predetermined number of sub-signals K, for instance, if the number of IMFs is M, then the parameter M used to determine the number of decomposition result, therefore, the resulting central frequency corresponding to different parameter K from 2 to 7 are tabulated in Table 1.As shown in Table 1, it can be seen that the K varying from 2 to 6 meets the condition that all center frequencies are better separated by the uniform distribution, while K = 7 shows that the central frequency 1699 Hz is very close to the 2100 Hz, and the natural frequency of excited structure (2000 Hz from Equation ( 31)) disappeared, thus the best result is given by VMD with K = 6.The decomposed results using VMD and their magnitude spectrum have been shown in Figure 6a,b.It is observed that all the signatures are successfully identified and clearly shown in 1st to 6-th IMF model.Here, according to the criterion of maximum kurtosis, i.e., the maximum kurtosis criterion component contains the main fault information of the mixed signal.The kurtosis is calculated to select the best fault sub-band for the envelope analysis, and the kurtosis of each IMF mode is as follows: KIMF1 = 2.2802, KIMF2 = 6.7916,KIMF3 = 3.2967, KIMF4 = 3.0748, KIMF5 = 3.0472 and KIMF6 = 3.1231.Hence, the IMF2 component was elected to extract the fault characteristic frequency.Figure 6c shows the envelop spectrum of IMF2 component.It was observed that the fault characteristic frequency f = 100 Hz and its harmonics (i.e., 2f, 3f, 4f, 5f, 6f and 7f ) are extracted by using VMD method.Meanwhile, the amplitudes of spectral peak increase dramatically, thus making the extraction of fault frequencies easy and accurate.In summary, the experimental results fully illustrate that the proposed approach was able not only to detect the fault frequencies in a weak fault signal but also to remove noise well.

Case 1-Bearing Fault Experimental Setup
The bearing fault vibration signal was generated by the Intelligent Maintenance Systems (IMS) and University of Cincinnati (Cincinnatti, OH, USA) [36].Four ZA-2115 double row bearings (Rexnord, Milwaukee, WI, USA) were installed on the shaft as shown in Figure 7.The rotation speed is kept constant at 2000 rpm and driven by an AC motor and coupled by rub belts.Vibration data was collected every 10 min by a NI 6062E-DAQ-Card (National Instruments, Austin, TX, USA), the sampling rate set at 20 kHz and the data length is 20,480 points.This accelerated life test was carried out successively for eight days (from 10:32:39 on 12 February 2004 to 06:22:39 on 19 February 2004, from normal to severe faults, i.e., 9840 min) until the magnetic plug exceeds a certain level and causes an electrical switch to close.The geometric parameters and ball pass frequency outer (BPFO) race of the tested bearing are presented in Table 2.
Entropy 2017, 19, 317 14 of 26 Figure 6c shows the envelop spectrum of IMF2 component.It was observed that the fault characteristic frequency f = 100 Hz and its harmonics (i.e., 2f, 3f, 4f, 5f, 6f and 7f) are extracted by using VMD method.Meanwhile, the amplitudes of spectral peak increase dramatically, thus making the extraction of fault frequencies easy and accurate.In summary, the experimental results fully illustrate that the proposed approach was able not only to detect the fault frequencies in a weak fault signal but also to remove noise well.

Case 1-Bearing Fault Experimental Setup
The bearing fault vibration signal was generated by the Intelligent Maintenance Systems (IMS) and University of Cincinnati (Cincinnatti, OH, USA) [36].Four ZA-2115 double row bearings (Rexnord, Milwaukee, WI, USA) were installed on the shaft as shown in Figure 7.The rotation speed is kept constant at 2000 rpm and driven by an AC motor and coupled by rub belts.Vibration data was collected every 10 min by a NI 6062E-DAQ-Card (National Instruments, Austin, TX, USA), the sampling rate set at 20 kHz and the data length is 20,480 points.This accelerated life test was carried out successively for eight days (from 10:32:39 on 12 February 2004 to 06:22:39 on 19 February 2004, from normal to severe faults, i.e., 9840 min) until the magnetic plug exceeds a certain level and causes an electrical switch to close.The geometric parameters and ball pass frequency outer (BPFO) race of the tested bearing are presented in Table 2.

Case 1-Results and Discussion of Bearing Single Fault
This accelerated life test was carried for 8 days, the equivalent vibration intensity (EVI) entire life-cycle curve and kurtosis entire life-cycle curve of rolling bearings 1 are respectively shown in Figure 8a,b.The detailed algorithm of EVI is given in our previous work [37].Figure 8b reveals that the evolution of kurtosis can be divided into three stages, i.e., normal operation, incipient fault and severe fault stage.It also can be seen in Figure 8b that the failure scope of bearing

Case 1-Results and Discussion of Bearing Single Fault
This accelerated life test was carried for 8 days, the equivalent vibration intensity (EVI) entire life-cycle curve and kurtosis entire life-cycle curve of rolling bearings 1 are respectively shown in Figure 8a,b.The detailed algorithm of EVI is given in our previous work [37].Figure 8b reveals that the evolution of kurtosis can be divided into three stages, i.e., normal operation, incipient fault and severe fault stage.It also can be seen in Figure 8b that the failure scope of bearing 1 is from the 6470th minute (incipient failure point) to the 9840th minute (severe damage point).During the first five days of operation, no underlying trend was observed.After the test had been carried out for 6490 min (at the end of the 5th day), the EVI and kurtosis of the vibration measurements increased significantly in the serious damage stage.
It also can be seen in Figure 8a,b that the range of failure probability is from incipient failure time to severe failure time, namely 6210th minute-8700th minute.Considering the early weak fault maybe submersed in heavy background noise, four sets of data (3rd day, 4th day, 5th day and 6th day, with the corresponding operation times being 2980, 4420, 5860 and 7300 min, respectively) from the experimental system were analyzed.
During the first five days of operation, no underlying trend was observed.After the test had been carried out for 6490 min (at the end of the 5th day), the EVI and kurtosis of the vibration measurements increased significantly in the serious damage stage.It also can be seen in Figure 8a,b that the range of failure probability is from incipient failure time to severe failure time, namely 6210th minute-8700th minute.Considering the early weak fault maybe submersed in heavy background noise, four sets of data (3rd day, 4th day, 5th day and 6th day, with the corresponding operation times being 2980, 4420, 5860 and 7300 min, respectively) from the experimental system were analyzed.A summary of the vibration acceleration signals of the 3rd day discussed in this paper is given in Figure 9.The signal in the top row Figure 9a represents the raw signal without any processing.It is seen that no impacts can be located and the kurtosis value is 3.5192 (it is approximately a normal distribution), it is easy to conclude that the bearing is running stably.In the second row Figure 9b, the improved AR-MED filtered signal is presented.It can be seen that the amplitude is decreased substantially by the improved AR-MED method.The third row Figure 9c shows the IMF components of improved AR-MED denoised signal decomposed by using the VMD method, which can be help distinguish the periodic impulse signals from the mixed noisy signals.As result, the ball pass frequency outer (BPFO) signal is completely contaminated by the heavy background noise which can be confirmed from the envelope spectrum as presented in Figure 9d, here, the kurtosis value of the second IMF model is the largest.A similar analysis result can be implemented in Figure 10, where the ball pass frequency outer (BPFO) also cannot be found in the corresponding envelope spectrum.A summary of the vibration acceleration signals of the 3rd day discussed in this paper is given in Figure 9.The signal in the top row Figure 9a represents the raw signal without any processing.It is seen that no impacts can be located and the kurtosis value is 3.5192 (it is approximately a normal distribution), it is easy to conclude that the bearing is running stably.In the second row Figure 9b, the improved AR-MED filtered signal is presented.It can be seen that the amplitude is decreased substantially by the improved AR-MED method.The third row Figure 9c shows the IMF components of improved AR-MED denoised signal decomposed by using the VMD method, which can be help distinguish the periodic impulse signals from the mixed noisy signals.As result, the ball pass frequency outer (BPFO) signal is completely contaminated by the heavy background noise which can be confirmed from the envelope spectrum as presented in Figure 9d, here, the kurtosis value of the second IMF model is the largest.A similar analysis result can be implemented in Figure 10, where the ball pass frequency outer (BPFO) also cannot be found in the corresponding envelope spectrum.Figure 11a-d, respectively depict the raw time domain signal of the 4th day (5860 min) along with its improved AR-MED filtered signal, the IMF components and the corresponding envelope spectrum.In this case, the envelop spectrum shown in Figure 11d indicates that the characteristic frequency at 230.7 Hz and its harmonic (230.7 × 2 Hz) are extracted (it is roughly equal to the BPFO frequency), which means that there is a fault in the outer race, whereas the low-frequency noise frequencies located around it are still significant.
Entropy 2017, 19, 317 17 of 26 Figure 11a-d, respectively depict the raw time domain signal of the 4th day (5860 min) along with its improved AR-MED filtered signal, the IMF components and the corresponding envelope spectrum.In this case, the envelop spectrum shown in Figure 11d indicates that the characteristic frequency at 230.7 Hz and its harmonic (230.7 × 2 Hz) are extracted (it is roughly equal to the BPFO frequency), which means that there is a fault in the outer race, whereas the low-frequency noise frequencies located around it are still significant.Figure 12d depicts the corresponding envelope spectrum with peaks at 230.7 Hz and its harmonics (2f0, 3f0), which shows that the BPFO peak is highlighted while the noise interference can be ignored.By comparing Figures 11d and 12d, as the failure progresses, the amplitudes of the BPFO peaks increase, thus making the feature monitoring easier and more accurate.
Moreover, disassembling the roller bearing and checking the running performance, the engineers found that an outer race with localized wear defect was discovered in test bearing 1 as shown in Figure 13.Indeed, it should be noted that 5860th minute is earlier than the incipient fault time 6490th minute observed in the kurtosis curves (see Figure 8b), on this occasion, the impulse features are buried in the strong background noise and are difficult to identify by the traditional kurtosis curve.Figure 12d depicts the corresponding envelope spectrum with peaks at 230.7 Hz and its harmonics (2f 0 , 3f 0 ), which shows that the BPFO peak is highlighted while the noise interference can be ignored.By comparing Figures 11d and 12d, as the failure progresses, the amplitudes of the BPFO peaks increase, thus making the feature monitoring easier and more accurate.
Moreover, disassembling the roller bearing and checking the running performance, the engineers found that an outer race with localized wear defect was discovered in test bearing 1 as shown in Figure 13.Indeed, it should be noted that 5860th minute is earlier than the incipient fault time 6490th minute observed in the kurtosis curves (see Figure 8b), on this occasion, the impulse features are buried in the strong background noise and are difficult to identify by the traditional kurtosis curve.Furthermore, in order to confirm the advantages of our proposed method, the same vibration signal of 4th day (5860 min) is processed via the common wavelet and LMD methods.Figure 14a and 15a show the wavelet sub-band and PF components, in which the raw vibration signals were divided into eight components, respectively.As illustrated in Figures 14a and 15a, where the high-frequency models represent fast oscillations and the low-frequency models represent low oscillations in Wavelet and LMD.Generally, the useful fault information was contained in high-frequency models, based on the criterion of maximum kurtosis, the sub-band component a7 and component PF2 were applied for fault feature extraction, respectively.
The kurtogram of the a7 sub-band component is displayed in Figure 14b, and the kurtogram of the PF2 component is displayed in Figure 15b, from which the optimal demodulation frequency band, namely (0 Hz-4896 Hz) and (5104 Hz-7396 Hz) can be detected.Thus a band-pass filter is designed to extract the potential features from the a7 sub-band component and PF2 component.Lastly, the envelope spectrum is applied to the filtered signal, the corresponding envelop spectrums are illustrated in Figures 14c and 15c, respectively.Furthermore, in order to confirm the advantages of our proposed method, the same vibration signal of 4th day (5860 min) is processed via the common wavelet and LMD methods.Figure 14a and 15a show the wavelet sub-band and PF components, in which the raw vibration signals were divided into eight components, respectively.As illustrated in Figures 14a and 15a, where the high-frequency models represent fast oscillations and the low-frequency models represent low oscillations in Wavelet and LMD.Generally, the useful fault information was contained in high-frequency models, based on the criterion of maximum kurtosis, the sub-band component a7 and component PF2 were applied for fault feature extraction, respectively.
The kurtogram of the a7 sub-band component is displayed in Figure 14b, and the kurtogram of the PF2 component is displayed in Figure 15b, from which the optimal demodulation frequency band, namely (0 Hz-4896 Hz) and (5104 Hz-7396 Hz) can be detected.Thus a band-pass filter is designed to extract the potential features from the a7 sub-band component and PF2 component.Lastly, the envelope spectrum is applied to the filtered signal, the corresponding envelop spectrums are illustrated in Figures 14c and 15c, respectively.Furthermore, in order to confirm the advantages of our proposed method, the same vibration signal of 4th day (5860 min) is processed via the common wavelet and LMD methods.Figures 14a and  15a show the wavelet sub-band and PF components, in which the raw vibration signals were divided into eight components, respectively.As illustrated in Figures 14a and 15a, where the high-frequency models represent fast oscillations and the low-frequency models represent low oscillations in Wavelet and LMD.Generally, the useful fault information was contained in high-frequency models, based on the criterion of maximum kurtosis, the sub-band component a7 and component PF2 were applied for fault feature extraction, respectively.
The kurtogram of the a7 sub-band component is displayed in Figure 14b, and the kurtogram of the PF2 component is displayed in Figure 15b, from which the optimal demodulation frequency band, namely (0 Hz-4896 Hz) and (5104 Hz-7396 Hz) can be detected.Thus a band-pass filter is designed to extract the potential features from the a7 sub-band component and PF2 component.Lastly, the envelope spectrum is applied to the filtered signal, the corresponding envelop spectrums are illustrated in Figures 14c and 15c      As can be seen, the fault characteristic frequency 236.4 Hz cannot be detected by the wavelet and LMD methods in the envelope spectrum and it is also hard to distinguish the fault location from the incipient fault signal.On the whole, the proposed technique presents good performance in  As can be seen, the fault characteristic frequency 236.4 Hz cannot be detected by the wavelet and LMD methods in the envelope spectrum and it is also hard to distinguish the fault location from the incipient fault signal.On the whole, the proposed technique presents good performance in enhancing weak periodic signals corrupted by heavy background noise, thus re-confirming the effectiveness and efficiency in terms of the weak fault characteristic extraction.

Case 2-Gearbox Fault Experimental Setup
The gearbox fault vibration signals were collected by Qian-Peng Diagnosis Engineering Co., Ltd.(Zhenjiang, China).The rotation speed of drive gear shaft was kept constant at 1470 rpm.Vibration data was collected in the driving end by an accelerometer ADA16-8/2 (LPCI, Zhenjiang, China), the sampling frequency is 5120 Hz and the data length is 53,248 points.This multiple-fault in the gearbox couples a localized pitting in the bull gear and a localized wear fault in the pinion gear.Gear material is S45C steel, the module of bull and pinion gear is 2, the tooth number of pinion gear and bull gear is Z 1 = 55 and Z 2 = 75, respectively, i.e., gear number ratio is Z 1 :Z 2 = 11:15, the pressure angle is 20

Case 2-Results and Discussion of Gearbox Multiple-Fault
Figure 17a illustrates the time domain signal of gear pair without any preprocessing.Figure 17b shows envelope spectrum of the raw signal.From Figure 17b, it has no distinctive fault frequency, which was buried in the strong background noise.In order to accurately extract multiple fault impulses from the raw signal, the proposed improved AR-MED method was used to the raw signal, and the result of improved AR-MED denoising signal is given in Figure 18a, where the filtered signal clearly shows multiple impulses that are near the harmonics of the defect frequencies.Moreover, from Figure 18a, we also can see that the interval between the two periodic impulses is approximately equal to 0.0562s, which is roughly equal to the pitting fault frequency (17.96Hz) of the tested gear.Then the VMD method is utilized to analyze the filtered signal.Figure 18b,c show the spectrum of IMF modes extracted by VMD, where it is quite clear that the center frequencies coincide, a phenomenon (mode duplication) that emerged in model number K = 20, which demonstrate that model number K = 19 is better than K = 20.In order to accurately extract multiple fault impulses from the raw signal, the proposed improved AR-MED method was used to the raw signal, and the result of improved AR-MED denoising signal is given in Figure 18a, where the filtered signal clearly shows multiple impulses that are near the harmonics of the defect frequencies.Moreover, from Figure 18a, we also can see that the interval between the two periodic impulses is approximately equal to 0.0562 s, which is roughly equal to the pitting fault frequency (17.96Hz) of the tested gear.Then the VMD method is utilized to analyze the filtered signal.Figure 18b,c show the spectrum of IMF modes extracted by VMD, where it is quite clear that the center frequencies coincide, a phenomenon (mode duplication) that emerged in model number K = 20, which demonstrate that model number K = 19 is better than K = 20.In order to accurately extract multiple fault impulses from the raw signal, the proposed improved AR-MED method was used to the raw signal, and the result of improved AR-MED denoising signal is given in Figure 18a, where the filtered signal clearly shows multiple impulses that are near the harmonics of the defect frequencies.Moreover, from Figure 18a, we also can see that the interval between the two periodic impulses is approximately equal to 0.0562s, which is roughly equal to the pitting fault frequency (17.96Hz) of the tested gear.Then the VMD method is utilized to analyze the filtered signal.Figure 18b,c show the spectrum of IMF modes extracted by VMD, where it is quite clear that the center frequencies coincide, a phenomenon (mode duplication) that emerged in model number K = 20, which demonstrate that model number K = 19 is better than K = 20.which leads to a lot of inconvenience in the engineering application, so an adaptive calculation algorithm should be established in this regard in future works.

Figure 2 .
Figure 2. The time domain waveform of the multi-component AM-FM signal x(t).

Figure 2 .Figure 3 .
Figure 2. The time domain waveform of the multi-component AM-FM signal x(t).

Figure 3 .
Figure 3. Decomposition results and time-frequency distribution of the multivariate signal with noise using EEMD, LMD and VMD methods.(a) IMF components of x(t) decomposed by using EEMD; (b) Hilbert spectrum distribution of the EEMD results; (c) PF components of x(t) decomposed by using LMD; (d) Hilbert spectrum distribution of the LMD results; (e) IMF components of x(t) decomposed by using VMD; (f) Hilbert spectrum distribution of the VMD results.

Figure 4 .
Figure 4. Simulated signal and its envelope spectrum.(a) Periodic impulse signal; (b) Raw fault signal with heavy noise; (c) Envelope spectrum of the raw fault signal.

Figure 4 .
Figure 4. Simulated signal and its envelope spectrum.(a) Periodic impulse signal; (b) Raw fault signal with heavy noise; (c) Envelope spectrum of the raw fault signal.

Figure 5 .
Figure 5. Filtered signal and its envelope spectrum with AR-MED and improved AR-MED method.(a) Filtered signal with AR-MED method; (b) Envelope spectrum of AR-MED filtered signal; (c) Filtered signal with improved AR-MED method; (d) Envelope spectrum of improved AR-MED filtered signal.

2 Figure 5 .
Figure 5. Filtered signal and its envelope spectrum with AR-MED and improved AR-MED method.(a) Filtered signal with AR-MED method; (b) Envelope spectrum of AR-MED filtered signal; (c) Filtered signal with improved AR-MED method; (d) Envelope spectrum of improved AR-MED filtered signal.

Figure 6 .Figure 6 .
Figure 6.Decompose results using the VMD and spectrum analysis.(a) IMF components of improved AR-MED de-noised signal decomposed by VMD; (b) The amplitude spectrum of each IMF component; (c) Envelope spectrum of IMF2.

Figure 7 .
Figure 7. Bearing test rig; (a) Ball bearing fault test rig; (b) The schematic diagram of accelerated life test.

Figure 7 .
Figure 7. Bearing test rig; (a) Ball bearing fault test rig; (b) The schematic diagram of accelerated life test.

Figure 8 .
Figure 8.Time feature (a) equivalent vibration intensity and (b) kurtosis of bearing 1 for the whole life cycle.

Figure 8 .
Figure 8.Time feature (a) equivalent vibration intensity and (b) kurtosis of bearing 1 for the whole life cycle.

Figure 9 .Figure 9 .Figure 10 .
Figure 9. Results using the proposed method of method the 3rd day's data (2980 min); (a) The original bearing vibration signal; (b) The improved AR-MED filtered signal; (c) The IMF components of improved AR-MED filtered signal generated by VMD method; (d) The envelope spectrum of IMF2.

Figure 10 .
Figure 10.Results using the proposed results of the 4th day's data (4420 min); (a) The original bearing vibration signal; (b) The improved AR-MED filtered signal; (c) The IMF components of improved AR-MED filtered signal generated by VMD method; (d) The envelope spectrum of IMF2.

Figure 11 .
Figure 11.Results using the proposed method method on the 5th day's data (5860 min); (a) The original bearing vibration signal; (b) The improved AR-MED filtered signal; (c) The IMF components of improved AR-MED filtered signal generated by VMD method; (d) The envelope spectrum of IMF2.

Figure 11 .
Figure 11.Results the proposed method method on the 5th day's data (5860 min); (a) The original bearing vibration signal; (b) The improved AR-MED filtered signal; (c) The IMF components of improved AR-MED filtered signal generated by VMD method; (d) The envelope spectrum of IMF2.

Figure 12 .
Figure 12. Results using the proposed method on the 6th day's data (7300 min); (a) The original bearing vibration signal; (b) The improved AR-MED filtered signal; (c) The IMF components of improved AR-MED filtered signal generated by VMD method; (d) The envelope spectrum of IMF2.

Figure 13 .
Figure 13.Photo of outer race defect in bearing 1 after the test.

Figure 12 .Figure 12 .
Figure 12. Results using the proposed method on the 6th day's data (7300 min); (a) The original bearing vibration signal; (b) The improved AR-MED filtered signal; (c) The IMF components of improved AR-MED filtered signal generated by VMD method; (d) The envelope spectrum of IMF2.

Figure 13 .
Figure 13.Photo of outer race defect in bearing 1 after the test.

Figure 13 .
Figure 13.Photo of outer race defect in bearing 1 after the test.

Figure 14 .
Figure 14.Results using Wavelet method.(a) Sub-band components of original signal decomposed by using Wavelet; (b) Kurtogram of a7 sub-band component; (c) Fourier transform magnitude of the squared envelope.

Figure 14 .
Figure 14.Results using Wavelet method.(a) Sub-band components of original signal decomposed by using Wavelet; (b) Kurtogram of a7 sub-band component; (c) Fourier transform magnitude of the squared envelope.

Figure 15 .
Figure 15.Results using LMD method.(a) PF components of original signal decomposed by using LMD; (b) Kurtogram of PF2 component; (c) Fourier transform magnitude of the squared envelope.

Figure 15 .
Figure 15.Results using LMD method.(a) PF components of original signal decomposed by using LMD; (b) Kurtogram of PF2 component; (c) Fourier transform magnitude of the squared envelope.
• and face width of teeth is 20 mm.The overall instrument configuration and the real equipment layout figure of the gearbox test rig are shown in Figure 16a,b, respectively.The pitting fault frequency f pitting of the bull gear is f pitting = 17.96Hz and the wear fault frequency of pinion gear is f wear = 33.41Hz.Entropy 2017, 19, 317 21 of 26 enhancing weak periodic signals corrupted by heavy background noise, thus re-confirming the effectiveness and efficiency in terms of the weak fault characteristic extraction.5.3.Case 2-Gearbox Fault Experimental SetupThe gearbox fault vibration signals were collected by Qian-Peng Diagnosis Engineering Co., Ltd.(Zhenjiang, China).The rotation speed of drive gear shaft was kept constant at 1470 rpm.Vibration data was collected in the driving end by an accelerometer ADA16-8/2 (LPCI, Zhenjiang, China), the sampling frequency is 5120 Hz and the data length is 53,248 points.This multiple-fault in the gearbox couples a localized pitting in the bull gear and a localized wear fault in the pinion gear.Gear material is S45C steel, the module of bull and pinion gear is 2, the tooth number of pinion gear and bull gear is Z1 = 55 and Z2 = 75, respectively, i.e., gear number ratio is Z1:Z2 = 11:15, the pressure angle is 20° and face width of teeth is 20 mm.The overall instrument configuration and the real equipment layout figure of the gearbox test rig are shown in Figure16a,b, respectively.The pitting fault frequency fpitting of the bull gear is fpitting = 17.96Hz and the wear fault frequency of pinion gear is fwear = 33.41Hz.

Figure 16 .
Figure 16.Gearbox test rig; (a) The engineering drawing of gearbox test rig; (b) The equipment layout figure of real gearbox test rig.

Figure 16 .
Figure 16.Gearbox test rig; (a) The engineering drawing of gearbox test rig; (b) The equipment layout figure of real gearbox test rig.

Entropy 2017, 19 , 317 22 of 26 5. 4 .
Figure17aillustrates the time domain signal of gear pair without any preprocessing.Figure17bshows envelope spectrum of the raw signal.From Figure17b, it has no distinctive fault frequency, which was buried in the strong background noise.

Figure 17 .
Figure 17.Raw vibration signal and its envelope spectrum.(a) Raw vibration signal; (b) Envelope spectrum of the raw fault signals.

Figure 17 .
Figure 17.Raw vibration signal and its envelope spectrum.(a) Raw vibration signal; (b) Envelope spectrum of the raw fault signals.

Entropy 2017, 19 , 317 22 of 26 5. 4 .
Figure17aillustrates the time domain signal of gear pair without any preprocessing.Figure17bshows envelope spectrum of the raw signal.From Figure17b, it has no distinctive fault frequency, which was buried in the strong background noise.

Figure 17 .
Figure 17.Raw vibration signal and its envelope spectrum.(a) Raw vibration signal; (b) Envelope spectrum of the raw fault signals.
illustrates the envelope spectrum of AR-MED filtered

Table 1 .
Central frequency corresponding to different K values.

Table 1 .
Central frequency corresponding to different K values.

Table 2 .
Geometric parameters and expected characteristic frequencies of the tested bearing.

Table 2 .
Geometric parameters and expected characteristic frequencies of the tested bearing.