A Novel Fault Diagnosis Scheme for Rolling Bearing Based on Convex Optimization in Synchroextracting Chirplet Transform

Synchroextracting transform (SET) developed from synchrosqueezing transform (SST) is a novel time-frequency (TF) analysis method. Its concentrated TF spectrum is obtained by applying a synchroextracting operator into TF transformation co-efficients on the TF plane. For this class of post-processing TF analysis methods, the main research focuses on the accurate estimation of instantaneous frequency (IF). However, the performance of TF analysis is greatly affected by the strong frequency modulation (FM) signal. In particular, the actual measured mechanical vibration signals always contain strong background noise, which decreases the resolution of TF representation, resulting in an inaccurate ridge extraction. To solve this problem, an improved penalty function based on the convex optimization scheme is firstly introduced for signal denoising. Based on the superiority of the linear chirplet transform (LCT) in dealing with modulated signals, the synchroextracting chirplet transform (SECT) is employed to sharpen the TF representation after the convex optimization denoising operation. To verify the effectiveness of the proposed method, the numerical simulation signals and the measured fault signals of rolling bearing are carried out, respectively. The results demonstrate that the proposed method leads to a better solution in rolling bearing fault feature extraction.


Introduction
As one of the most important components in rotating machines, the stable and regular running of rolling bearing can ensure the reliability of whole mechanical equipment. Therefore, it is of great significance to realize the structural health monitoring of rolling bearing [1,2]. Due to the complexity of the mechanical transmission chain, the measured vibration signals of rolling bearing always present strong non-stationary characteristics [3,4]. However, the traditional frequency-domain analysis method based on a fast Fourier transform (FFT) can only be utilized to process stationary signals. Therefore, the joint time-frequency analysis approach, which can reveal the time-varying properties of a signal, has been greatly studied by researchers [5][6][7].
Commonly, the collected vibration signal is a typical non-stationary time series in engineering practice. More characteristic information can be obtained if the signal is transformed to the time-frequency (TF) domain. Undoubtedly, the TF transform is an important technique for representations can be identified by the SECT method, and the proposed method in this paper is called the improved SECT (ISECT). In order to verify the performance of the proposed method ISECT, a numerical simulation signal analysis and a feature extraction to measure the rolling bearing vibration signal with outer ring fault have been carried out, respectively.
The rest of this paper is arranged as follows: In Section 2, we give a theoretical introduction to convex optimization based on the GMC penalty and the synchroextracting chirplet transform (SECT). In Section 3, the frequency modulation signal with added Gaussian white noise is presented to verify the effectiveness of the proposed ISECT method. In addition, the results acquired by STFT, SST, and SET are employed to make a comparison. In Section 4, the measured rolling bearing fault signals from the experiment rig and the industrial site are used to further demonstrate the capacity of ISECT. Conclusions are given in Section 5.

Convex Optimization Based on GMC Penalty
To obtain a more accurate approximation solution of the raw signal, a non-convex penalty function called the generalized mini-max concave (GMC) penalty is employed to substitute the traditional regularization term l 1 -norm in a convex optimization problem.
Firstly, we define the sparse approximate solution to the original signal. The specific solution process can be interpreted by minimizing the least squares cost function: where y ∈ R M is the raw signal, x 1 denotes the l 1 -norm operator to x, A ∈ R M×N is an over-sampled inverse discrete short-time Fourier transform operator, and λ > 0 corresponds to the regularization parameter. We implement the STFT as a normalized tight frame, i.e., AA H = I. Theoretically, l 0 norm is difficult to be optimized (NP problem), and l 1 norm is the optimal convex approximation of l 0 norm. It also should be pointed out that l 1 norm is easier to be optimized than l 0 norm. However, this may underestimate the true solution. Some studies have shown that a non-convex penalty function as a regularization term can maintain the convexity of the cost function. In this paper, we introduce the non-convex penalty function GMC as the regularization [25], which is defined as a generalization of the Huber function and the mini-max concave penalty function.
The Huber function is firstly defined as: Then, the mini-max concave penalty function φ(x) can be expressed as: According to Equations (2) and (3), we obtain the mini-max concave penalty function φ(x) as follows: where s(x) is the Huber function, and |x| represents the absolute value of x. Subsequently, a scaling version of the Huber function s b is proposed for the multi-variate generalization of φ(x): where b 0. From Equation (5), the corresponding scaled MC penalty is described as: Similar to Equation (4), the scaling penalty function φ b (x) can be denoted as: where s b is the scaling MC penalty function. Then, let B ∈ R M×N , and we can thus define the generalized Huber function S B : By generalizing the l 1 -norm and the Huber function, the generalized mini-max concave penalty function can be expressed as: where S B is the generalized Huber function and x 1 represents the l 1 -norm of x.
Ultimately, Equation (1) can be rewritten as: where ψ B is the generalized mini-max concave penalty function, which is parameterized by matrix B.
Additionally, the selection of matrix B significantly depends on A. In order to maintain the convexity of the above objective function defined in Equation (10), the convexity condition is presented as follows: When the inverse short-time Fourier transform matrix A is given, Equation (11) can be illustrated as B = γ/λ A, 0 ≤ γ ≤ 1. It is obvious that the penalty function is equivalent to the l 1 -norm when we let γ = 0. In the experiment, the optimal range of γ can be chosen as 0.5 ≤ γ ≤ 0.8.
To address the problem of minimization in Equation (10), we need to transform it into a typical saddle-point problem: where is the saddle function, and x opt is the optimized result. Furthermore, we can solve the minimization problem by the forward-backward (FB) algorithm [27].

Synchroextracting Chirplet Transform
The chirplet transform (CT) is generalized from the wavelet transform (WT) and has many excellent properties of the short-time Fourier transform (STFT), which is defined as where f (t) is the analyzed signal, g(t) is a real and even window function in the Schwartz class, and β is the chirplet ratio (CR) parameter.
Based on the CT model, a new rotating parameter α is introduced: where T s corresponds to the sampling time, and F s is the sampling frequency. We assume the parameter α can be described as: Based on the former description, we can rewrite the expression of a chirplet transform as: Theoretically, CT has a satisfied TF representation around its IF with concentrated energy.
Moreover, the amplitude C g f (t, ω, β) can reach the maximum among all values. Thus, the best optimization solution of α from the amplitude of C g f (t, ω, β) can be expressed as: Then, the matching CT C g f (t, ω, β) can be calculated as where C g f (t, ω, α * ) > λ and β * = F s 2T s tan(α * ). Therefore, a novel TF representation called the synchroextracting chirplet transform (SECT) [20] is proposed as: It should be pointed out that the δ(ω − ∼ ω(t, w)) is a synchroextracting operator, and it can be expressed as: To achieve the high-resolution representations of a multi-component signal, the proposed method ISECT has combined the advantage of a convex optimization scheme based on the GMC penalty and the TF analysis method of SECT. The specific process of the approach proposed in this paper is shown in Figure 1.

Simulated Signal Analysis
To verify the effectiveness of the proposed method, a simulated signal analysis is performed firstly. Commonly, mechanical equipment vibration signals consist of amplitude modulation (AM) and frequency modulation (FM) signals under strong background noise [28]. Thus, the simulated signal is defined as: x 2 = cos(2π(10t + 5 cos(t))) (23) where the simulated signal x is composed of three parts: x 1 , x 2 , and n. Obviously, x 1 is a typical AM-FM signal, and x 2 is a FM signal. The characteristic frequency is defined as f 1 = 15 Hz, f 2 = 20 Hz, and f 3 = 1 Hz. It should be noted that the ideal instantaneous frequency of x 1 and x 2 can be calculated as IF 1 = 20 − 2 sin(2πt) and IF 2 = 10 − 5 sin(t), respectively. The symbol n represents Gaussian noise components with SNR = 5 dB. Figure 2 shows the time-domain waveform of signals x 1 and x 2 . The multi-component numerical simulation signal is displayed in Figure 3. It is obvious that the traditional signal processing methods such as FFT cannot be effective in identifying time-varying features such as IF 1 and IF 2 . The main reason is that FFT lacks the ability to analyze non-stationary and aperiodic signals. Subsequently, mainstream time-frequency analysis methods, such as STFT, SST, SET, and SECT, are used to analyze the simulated signal. The TF representations results are shown in Figure 4.   According to Figure 4, a traditional time-frequency analysis method like STFT has a poor ability in achieving concentrated TF representation. Compared with STFT, popular methods such as SST, SET, and SECT have obvious superiority in a TF analysis. To some extent, the ridge corresponding to the AM-FM signal x 1 and the FM signal x 2 can be inspected. Theoretically, based on the chirplet transform and the synchroextracting operator, SECT has better analysis results. SET and SECT are based on the mathematical foundations of SST. When the instantaneous frequencies of the interested components are squeezed and rearranged, the noise components are still distributed on the time-frequency plane. Thus, the calculated performance of the above-mentioned methods should be improved due to the existing noise components. The most effective solution is to first reduce the noisy components and then optimize the time-frequency representations. Figure 5 is the ideal time-frequency curve corresponding to x 1 and x 2 . According to the previous analysis, we can make a conclusion that the noise components have a negative effect in sharpening the TF ridges. Thus, the GMC denoising approach is conducted firstly, and then the method of SECT is used to achieve high-resolution TF representation. The result calculated by the proposed method is shown in Figure 6. It should be noted that the extracted ridge form of the optimized TF plane is shown in Figure 6b. Comparing Figure 5 with Figure 6b, it can be found that the proposed method has an obvious advantage in identifying the time-varying signal feature. Commonly, Rayleigh entropy is regarded as a significant criterion in reflecting the time-frequency aggregation. The smaller value of Rayleigh entropy always indicates a better performance of TF representation. The Rayleigh entropy calculated by different methods is listed in Table 1. The results demonstrate that the proposed method ISECT has a better performance in denoising and in sharpening TF ridges.

Case 1: Data Analysis of a Fault Test Rig
Through the dynamic analysis, it can be seen that the natural frequency decreases and the characteristic frequency increases during the fault deterioration of rolling bearing. Therefore, the natural frequency of a vibration signal can reflect the severity of the rolling bearing fault. To prove the effectiveness of the proposed method, an experimental data analysis of a fault test rig is performed. To facilitate the rolling bearing analysis, the data analysis is used for the fault feature frequency identification. The experimental dataset was provided by the Machinery Failure Prevention Technology (MFPT) Society [29]. The test rig with a NICE bearing is performed to gather acceleration data for baseline conditions at 300 lbs of load. The structure parameters of the rolling bearing is described in Table 2. The data acquisition parameters are 25 Hz and 48,828 Hz for the input shaft rate and sample rate, respectively. The detailed calculation process can be found in [30]; the failure frequency of outer ring is f 0 = 80 Hz. Figure 7 shows the outer ring fault.  Figure 8b. Thus, the advanced signal processing method should be considered. Then, the envelope spectrum analysis is performed, and the corresponding result is plotted in Figure 9. Unfortunately, the fault characteristic frequency still cannot be inspected. TF representations of the envelope signal are plotted in Figure 10. Because of time-frequency ambiguity resulting from noise, we can only find the tendency of the outer ring fault characteristic frequency.   The phenomenon of noise interference is the main problem to be solved in a TF analysis. One of the pre-processing means is to realize the signal denoising of complex multi-component signals. In this paper, the GMC denoising based on a convex optimization scheme is used to process the experimental data firstly. Subsequently, the SECT method is applied to the denoised signal, and the TF analysis result is shown in Figure 11. In order to clearly identify the features, a multiple ridge extraction to the SECT method is performed, and the result is plotted in Figure 12. According to Figure 12, we can clearly inspect the fault characteristic frequency of the outer ring and its harmonics. The diagnosis is consistent with the facts, and the superiority of the proposed ISECT method is demonstrated.

Case 2: Data Analysis of a Large Desulfurization Fan
The large desulfurization fan is a key equipment in a steel plant, and its main function is to remove dust and control environmental pollution. The structure diagram of the large desulfurization fan is shown in Figure 13. Once the equipment breaks down in the production process, the dust inside the converter cannot be removed, which not only pollutes the surrounding environment but also causes huge potential security problems. In order to ensure the normal operation of the large desulfurization fan, status monitoring and a diagnosis analysis on each part of the fan should be performed. The input shaft speed of the hydraulic coupler is 850 r/min, and the output shaft speed is 740 r/min. The specific model of the rolling bearings on both sides of the fan is 22344CA (NSK, Kunshan City, China) , and the number of rolls is 13. Bearing operation parameters, data acquisition parameters, and bearing failure frequency are shown in Table 3. It should be noted that the outer ring fault feature frequency is calculated as f 0 = 64 Hz. The measured rolling bearing vibration signal of the large desulfurization fan is drawn in Figure 14. We still cannot identify the outer ring fault feature frequency in the frequency spectrum. Moreover, the noise interference components are easy to be inspected. Then, the envelope spectrum analysis was conducted, and the corresponding result is plotted in Figure 15b. Nevertheless, the fault characteristic components of the signal are masked by noise. Thus, the analysis results provided by the FFT and the envelope spectrum are unsatisfied. TF representations of the envelope signal of the large desulfurization fan are shown in Figure 16. We can hardly identify the phenomenon of the outer ring fault features from the fuzzy time-frequency plane.   Subsequently, the denoising operation is executed by the GMC penalty based on the convex optimization framework, which can improve the signal to noise ratio. Then, SST and SET, a novel TF analysis method, is applied to the vibration signal obtained by the proposed GMC denoising. The results provided by the SST and SET methods are plotted in Figure 17, respectively. However, the ridge curves related to the fault characteristic components are not obvious. Thus, the resolution of the TF representations of the fan vibration signal needs to be further improved. Based on the GMC denoising scheme, SECT is put forward to achieve a concentrated TF plane. The denoised signal is processed by the SECT method, and the result is shown in Figure 18. The multiple ridge identification is plotted in Figure 19. This suggests that the interested signal component of the outer ring fault characteristic has a sharpened ridge, which is an obvious distinction from other extraneous components. This shows that the outer ring fault feature frequency and its multiplication can be easily extracted-that is to say, the effectiveness of the proposed method has been proven.

Conclusions
High-resolution time-frequency representations of multi-component signals are a key part in mechanical fault diagnosis. The noise is inevitable in the actual environment, resulting in large disturbances and unsatisfactory results. To overcome this shortcoming, the convex optimization with the GMC penalty function combined with SECT is introduced in this paper. The main contributions of this paper are summarized as follows: (1) The convex optimization algorithm with the GMC penalty function is used for signal denoising, which can enhance the readability of time-frequency representation; (2) the SECT algorithm is employed to achieve a concentrated time-frequency plane, which aims to accurately extract the ridge curves related to fault characteristic components; and (3) the proposed method ISECT can effectively extract the rolling bearing fault feature. Both the analysis results of the numerical simulation signal and the measured bearing failure data demonstrate that the proposed method has improved the time-frequency representations of multi-component signals. We recommend that the proposed approach can be applied to intelligent fault diagnosis.