A Novel Fault Diagnosis Method for Rolling Bearing Based on Improved Sparse Regularization via Convex Optimization

Key Laboratory of Metallurgical Equipment and Control Technology, Ministry of Education, Wuhan University of Science and Technology, Wuhan 430081, China Hubei Key Laboratory of Mechanical Transmission and Manufacturing Engineering, Wuhan University of Science and Technology, Wuhan 430081, China Engineering Research Center for Metallurgical Automation and Measurement Technology of Ministry of Education, Wuhan University of Science and Technology, Wuhan 430081, China Wenzhou Special Equipment Inspection and Research Institute, Wenzhou 325007, China


Introduction
In the field of prognostics and health management (PHM) to mechanical equipment, the actual collected vibration signal contains wealthy information about operating status [1][2][3].It is worthy to note that the faulty characteristic information can always be reflected from the measured signal.Rolling bearing is the key component in the main transmission system, and its operating performance is directly related to the status of the entire plant.When some faults occur in them, it will inevitably affect the work safety and production efficiency [4,5].Therefore, analyzing the possible fault characteristics of the rolling bearings has an important practical significance [6][7][8].Modern signal processing methods provide the main technical means of rolling bearing fault diagnosis.
Generally, the acquired signal is coupled by useful information and the strong noisy component, which has the typical nonlinear and nonstationary characteristics.Hence, the main task of mechanical equipment structural health monitoring is to effectively remove the noise component and improve the signal-to-noise ratio [9,10].Traditional signal processing method is based on inner product operator, which is built on the analyzed signal and basis function.The most representative ones are short-time Fourier transform (STFT) [11] and wavelet transform (WT) [12,13].Since the size of the analyzed window is fixed, STFT lacks sufficient capacity to deal with the complex nonstationary signals.The performance of WT depends on the selection of wavelet basis function and the decomposition level.Then, local mean decomposition (LMD) algorithm is proposed as an adaptive time-frequency analysis method [14].However, it still has the problem of mode aliasing.Variational mode decomposition (VMD) [15,16] is proposed based on Wiener filtering, one-dimensional Hilbert transform, and heterodyne demodulation analysis.It is still affected by the selection of penalty parameter and the number of signal components.Presently, synchrosqueezing transform (ST) has attracted much attention due to its properties in timefrequency reassignment [17,18].Nevertheless, its performance is unsatisfied since the cross-term interference and poor scale separation.
Essentially, the signal denoising can be achieved by the calculation of a sparse approximate solution to the measured vibration signal.Thus, the novel denoising method based on convex optimization and sparsity has publicly employed in signal processing and image enhancement [19,20].It has been successfully applied in the field of mechanical fault diagnosis [21], spectral data processing, and baseline correction [22].Total variation (TV) is the main content of convex optimization algorithm, which involves a quadratic data fidelity term and a convex regularization term [23].It is developed on sparse signal models.Based on that, the first-order TV [24] and higher degree TV (HDTV) [25] regularization has been researched.However, the experimental analysis demonstrated that it may work when the signal is piecewise constant and it often produces undesirable staircase artifacts.Subsequently, wavelet total variation denoising [26] is proposed to improve the denoising performance, while the estimation of noise variance and penalty function selection is an inevitable question.Using either (1) L 1 -norm regularization and convex optimization or (2) nonconvex regularization and nonconvex optimization, the calculation of a sparse approximate solution to a linear system of equations is often performed.However, it tends to provide solutions that deviate the real values.Afterwards, a new penalty is suggested, which is a multivariate generalization of the minimax-concave (MC) penalty and is defined as generalized MC (GMC) penalty [27,28].GMC penalty simultaneously involves the generalized Huber function and the regular L 1 -norm regularization.In other words, it is an innovative nonseparability and nonconvex penalty function, which is using for achieving the sparse enhancement and signal smoothing.A nonseparable penalty has more superiority in meeting the requirement of preserving the convexity of the objective function [29].The published result fully indicates that the GMC penalty is obviously superior to the common penalty function in noise artifacts removing for timefrequency analysis [27].
For early and weak fault detection, it has the feature of latent and dynamic response as the multifactor coupling and complex transmission path [30].Since the fault is in its early state, the energy generated during the operation is low, and the signal that can be received by the sensor after attenuation is extremely weak [31].Therefore, early failure signal is easily submerged by background noise, and effective extraction of these characteristics has been a difficult problem.Inspired by the idea of jointing the nonconvex penalty and convex optimization algorithm, the traditional L 1 -norm regularization term is replaced by GMC penalty in the TV denoising scheme in this paper so as to effectively realize fault state identification.In order to verify the rationality and feasibility of the proposed method, it is used to analyze the numerical simulation and the actual fault vibration signal of rolling bearings in the bearing test rig.The result demonstrated that the proposed method has obvious advantages over traditional methods such as WT, TV, and complete ensemble empirical mode decomposition with adaptive noise (CEEMDAN) [32].
The rest of the paper is organized as follows.In Section 2, the basic ideas of GMC penalty and the researched convex optimization denoising algorithm are introduced.The simulation signal analysis is described in Section 3. The faultbearing data from bearing test rig is analyzed in Section 4. Section 5 gives the final conclusions.

Theory Descriptions
2.1.Generalized Minimax-Concave (GMC) Penalty.Generally, the objective function is convex, and the optimization problem of constraint variable value in a convex set is called convex optimization problem.The penalty function is a measure of constraint violation, which makes the constraint zero when the constraint is satisfied.The one-dimension time series x ∈ R N can be simply expressed as x = x 1 , x 2 , … , x N .Then, the L 1 -norm and L 2 -norm is denoted as x 1 = ∑ n x n and x 2 = ∑ n x n 2 1/2 respectively.The Huber function s x can be defined as: Then, we propose a multivariate generalization of the MC penalty.The basic idea is to generalize (6) using the L 1 -norm and the generalized Huber function.Thus, we define the generalized MC (GMC) penalty function ψ B x as follows:

The Denoising Algorithm Based on Convex Optimization.
Let y ∈ R M be the original observed signal, A ∈ R M×N and λ > 0 is the regularization parameter.A commonly used approach to obtain an optimal sparse approximate solution is to minimize the following objective function [24]: Undoubtfully, this is a typical convex optimization problem of first-order TV denoising, which is comprised by a quadratic fidelity term and an L 1 -norm regularization term.To achieve the signal denoising and early fault feature extraction, a sparse regularization denoising algorithm based on convex optimization is proposed in this paper.Specifically speaking, the traditional L 1 -norm regularization item is replaced by the GMC penalty term and it is extended to the first-order TV denoising scheme.Thus, we redefine the objective function F R N → R as follows: where A is an oversampled inverse discrete Fourier transform, ψ B :R N → R is the generalized MC (GMC) penalty defined by (7), λ is the regularization parameter, and the symbol of Dx stands for the TV operator.The bidiagonal matrix D ∈ R n−1 ×n is defined as The penalty ψ B is parameterized by a matrix B and the convexity of F x depends on B being suitably prescribed.It also should be pointed out that the choice of B will depend on A. If then F x is a convex function.It is easy to satisfy the convexity condition.Given the matrix A ∈ R M×N N ≥ M , we may simply set  3 Complexity where M is the signal length and N is the transform length.It should be pointed out that 0 ≤ γ ≤ 1 always meets the convexity condition.Generally, we use a nominal range of 0 5 ≤ γ ≤ 0 8 so as to obtain better performance.In order to minimize the cost function by using the approximate algorithms, we rewrite it as a saddle-point problem: where 2 is the saddle function and x opt presents the denoised signal.
According to abovementioned, it is obvious that the selection of regularization parameter λ also has a significant influence in denoising performance.Commonly, we chose λ = 2 to achieve the convergence of algorithm.Therefore, the forward-backward (FB) iterative algorithm can be used to solve the problem F x, v of this kind of saddle point [33].The resulting iterative threshold algorithm uses the soft-threshold function, which is defined as The flowchart of the presented method in this paper is plotted in Figure 2.

Simulation Signal Analysis
Without loss of generality, the numerical signal analysis is used to simulate the rolling bearing fault feature identification.It is composed by frequency modulation signal,     5 Complexity harmonic signal, and strong background noise components.The numerical simulation signal can be expressed as follows:

15
where s 3 is expressed a Gaussian white noise with a variance of 1.5, the composite signal y is a typical noisy multicomponent signal and the feature frequency is, respectively, set as f 1 = 30 Hz, f 2 = 200 Hz, and f 3 = 100 Hz.
Figure 3 shows the time and frequency responses of the synthetic signal y with strong noisy signal components.Figure 3(b) suggests that only the feature frequency f 2 and f 3 can be inspected, and the baseline of frequency spectrum has generated a drift.Unfortunately, the phenomenon of frequency modulation, namely, f 2 ± f 1 , is hardly identified since the interference of noise components.Nevertheless, frequency modulation is an important tool in fault state identification.Thus, we can make a conclusion that the traditional frequency spectrum analysis method should be improved.The immediate idea is that an effective denoising algorithm is performed before frequency spectrum analysis.
Subsequently, the commonly used method wavelet denoising and total variation (TV) denoising algorithm have been employed to analyze the complex simulation signal.The wavelet denoising is performed by wavelet packet scheme.Moreover, the wavelet base function is determined as "db5," and the decomposition level is selected as 5.Then, the signal reconstruction is achieved using the     (3,0).Figure 4 has plotted the result generated by wavelet analysis in time domain and frequency domain.Judging from Figure 4, we can make a conclusion that the interesting components have been removed, and WT denoising method has failed to feature extraction.Traditional TV denoising algorithm has also been regarded as comparative analysis method in this paper, and the result is shown in Figure 5.According to Figure 5, we observe that most concerned signal components have been deleted, and its performance is still unsatisfied.Subsequently, CEEMDAN is employed to simulation signal analysis, and the result is drawn in Figure 6.Original signal is decomposed into nine intrinsic mode functions (IMF).For the parameter setting for CEEMDAN algorithm, the noise standard deviation is 0.2, and the number of realizations is 200.It can be seen from Figure 6 that only the feature frequency f 1 = 30 Hz can be found in IMF2 and IMF3.Motivated by the classical first-order total variation (TV) denoising scheme, a sparse regularization method based on convex optimization is presented in this paper.In order to verify the effectiveness of this method, the researched method based on optimization method has also been applied to it.The parameters of the proposed method are chosen as γ = 0 8 and λ = 2. Figure 7 shows the result obtained by the proposed method.It is obvious that the feature frequency f 2 and f 3 can be clearly inspected.Most importantly, the frequency modulation f 2 ± f 1 can also been determined.Comparing Figure 3(b) with Figure 7, it is obvious that the baseline in frequency domain has been corrected.The experimental results completely demonstrate that this proposed algorithm outperforms the classical wavelet denoising and TV method.
The actual noise reduction effect is directly related to the intensity of noise components.Theoretically, the smaller noise intensity will lead to the better actual denoising performance.When the intensity of Gaussian noise varies from 0.5 to 1.5, the root mean square error (RMSE) between the original noisy signal and the denoised signal is listed in Table 1.The result illustrated that the proposed method has advantage against the strong noise, and its robustness has been proved.The algorithms are run on a computer with an Intel Core i3-4160 CPU and 8.0 GB RAM.The computational costs of four signal processing methods for the simulated signal are listed in Table 2. Since the iterated operation with 200 times, the computational efficiency of CEEMDAN is low.It demonstrates that the computational complexity of the proposed method is acceptable.

Experimental Analyses
Rolling bearing is an important component of rotating machinery.Its main function is to support the mechanical rotary body, reduce the friction coefficient of its motion, and guarantee its accuracy.Practically, the measured bearing failure signal is more complicated than the numerical simulation signal.Hence, the proposed method is also used for experimental signal analysis.The test rig was equipped with a NICE bearing with the outer race fault shown in Figure 8.The red circle indicates the location of the outer ring fault.
The following parameters about the rolling bearing are listed in Table 3.The outer race fault conditions is expressed as follows: 300 lbs of load, input shaft rate f r = 25 Hz, and sample rate of f s = 48828 Hz.According to the theoretical calculation, the characteristic frequency of the outer ring fault is determined as f o = 80 Hz.
The time domain waveform and frequency spectrum analysis results of the collected vibration signals are plotted in Figure 9. From the Figure 9(a) of the time response about measured bearing fault signal, the impact characteristics are obvious.However, we cannot detect the characteristic frequency of the outer ring in Figure 9(b).The frequency spectrogram analysis result indicates that the components of measured faulty signal are complicated.As the interference of the strong noisy components, the fault feature frequency can hardly been identified.The ideal result is to remove unnecessary signal components by effective denoising algorithm, and to retain or amplify the related signal components.
Then, the common signal processing methods, such as wavelet analysis and TV denoising, have been applied to the real inspected signal.For the method of wavelet denoising, we perform the wavelet packet decomposition denoising scheme.It should be noted that the wavelet base function is determined as "db5," and the decomposition level is selected as 3.The coefficient of the node (2, 0) is used to signal reconstruction.Figures 10(a) and 10(b) show the results obtained by wavelet denoising and TV denoising, respectively.There is no obvious feature frequency corresponding to outer ring   10(a).Figure 10(b) suggested that the method of TV denoising has removed most of the useful signals, which is not suitable for vibration signal processing.
Similarly, CEEMDAN is also used to analyze the experimental data, and the corresponding result is plotted in Figure 11.Unfortunately, only the characteristic frequency of the outer ring fault f o and its double frequency 2f o can be inspected in IMF2, while the unwanted signal components have not been removed and it interfere with the identification of fault features.Since CEEMDAN still lacks theoretical support and exits the problem of mode aliasing, it fails to analyze the complex nonstationary vibration signal.
Finally, the proposed method based on GMC penalty function is performed to validate the method and illustrate its superiority.The concerned parameter of the proposed algorithm is chosen as γ = 0 5 and λ = 2, respectively.The computed result is drawn in Figure 12.Obviously, feature frequency of outer ring fault f o and its multiple frequencies (2f o , 3f o , 4f o , 5f o , and 6f o ) had been both identified in Figure 12.Meanwhile, the rotational frequency f r can also be determined.The abovementioned characteristics fully indicate that the fault occurs in the outer ring, which is in accordance with the actual situation.Through the comparative analysis between Figures 10 and 12, the results show that the proposed method has obvious advantages in the fault state recognition of rolling bearings under the environmental of strong noise.
Theoretically, traditional signal processing methods such as fast Fourier transform and wavelet transform are based on the idea of matching the analyzed signal to the base function.Its performance largely depends on the signal structure and the base function selection.CEEMDAN is a typical signal adaptive decomposition algorithm, while it is restricted by endpoint effect and mode mixing.Distinguished it, convex   8 Complexity optimization denoising algorithm such as the proposed method is a typical iterative optimization process for objective functions.Using the nonconvex penalty function or cost function to realize the high efficiency, it is more suitable for analyzing the complex vibration signals.

Conclusions
For rolling bearing fault diagnosis, a novel method based on the improved sparse regularization via convex optimization is proposed in this paper.The main findings of this paper include (1) based on L 1 -norm and the Huber function, a novel nonconvex penalty function has been researched, namely, GMC penalty.It is employed to achieve the sparse representation of the signal and estimate the sparse solutions more accurately.(2) The proposed method is firstly put forward by combining the GMC penalty and TV denoising scheme.It is designed to enhance the performance of noise reduction for the complex vibration signals.
(3) Compared with the traditional methods such as WT and TV denoising, the proposed method has better performance, as demonstrated by both numerical and experimental studies.

Figures 1 (
Figures 1(b) and 1(d) demonstrate the scaled Huber function and scaled MC penalty, respectively.Let B ∈ R M×N and we next define the generalized Huber function S B x as Scaled MC penalty (d) Scaled MC penalty function

Figure 2 :
Figure 2: The flowchart of the method presented in this paper.

Figure 3 :
Figure 3: The synthetic signal of simulation.

Figure 5 :
Figure 5: The result provided by the typical total variation (TV) denoising.
The result of frequency spectrum analysis by wavelet denoising

Figure 4 :
Figure 4: The result provided by wavelet denoising algorithm.

3 Figure 6 :
Figure 6: The result provided by CEEMDAN for simulation signal analysis.

Figure 7 :
Figure 7: The result provided by the proposed method.
The time response of measured bearing fault signal The frequency response of measured bearing fault signal

Figure 9 :
Figure 9: The time-frequency analysis of measured fault signal.

Figure 10 :
Figure 10: The result provided by wavelet denoising and TV denoising.

Figure 11 :Figure 12 :
Figure 11: The result provided by CEEMDAN for experimental signal analysis.

Table 1 :
The proposed algorithm in signal denoising.

Table 2 :
Computational costs of four signal processing methods for simulated signal.