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

Authors: Dongjie Zhong; Cancan Yi; Han Xiao; Houzhuang Zhang; Anding Wu
Journal: Complexity
 (2018)

---

## Abstract

Structural health monitoring and fault state identification of key components, such as rolling bearing, located in the mechanical main drive system, have a vital significance. The acquired fault signal of rolling bearing always presents the obvious nonlinear and nonstationary characteristics. Moreover, the concerned features are submerged in strong background noise. To handle this difficulty, a novel fault signal denoising scheme based on improved sparse regularization via convex optimization is proposed to extract the fault feature of rolling bearing. In this paper, the generalized minimax-concave (GMC) penalty is firstly researched to promote the sparsity of signal, which is based on traditional L1-norm and Huber function. It is designed to estimate the sparse solutions more accurately and maintain the convexity of the cost function. Then, the GMC penalty is extended to 1-D first-order total variation (TV) as nonseparability and nonconvex regularizer. Thus, a convex optimization problem, which involves a quadratic data fidelity term and a convex regularization term, is developed in this paper. To accelerate the convergence of the algorithm, it is solved by forward-backward (FB) iterative algorithm and thus the denoised signal can be obtained. In order to demonstrate its performance, the proposed method is illustrated for numerical simulation signal and applied in the feature extraction of the measured rolling bearing vibration signal.


---

## Body

## 1. 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–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–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 time-frequency 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) L1-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 L1-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 time-frequency 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 L1-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 Section2, 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 fault-bearing data from bearing test rig is analyzed in Section 4. Section 5 gives the final conclusions.


## 2. 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 seriesx∈RN can be simply expressed as x=x1,x2,…,xN. Then, the L1-norm and L2-norm is denoted as x1=∑nxn and x2=∑nxn21/2 respectively.The Huber functionsx can be defined as:
(1)sx=12x2,∣x∣≤1∣x∣−12,∣x∣≥1.The minimax-concave (MC) penaltyϕx is defined as:
(2)ϕx=∣x∣−12x2,∣x∣≤112,∣x∣≥1.Figures1(a) and 1(c) have illustrated the Huber function and MC penalty. Hence, the relationship between Huber function and MC penalty can be described as
(3)ϕx=x−sx.Figure 1
Different penalty functions.
(a)
Huber function
(b)
Scaled Huber function
(c)
MC penalty function
(d)
Scaled MC penalty functionDefining the scaled versions of the Huber function and MC penalty is convenient. Forb≠0, the scaled Huber function sbx is determined by
(4)sbx=12b2x2,∣x∣≤1b2∣x∣−12b2,∣x∣≥1b2.Similarly, the scaled MC penalty functionϕbx is given by
(5)ϕbx=x1−sbx.Figures1(b) and 1(d) demonstrate the scaled Huber function and scaled MC penalty, respectively. Let B∈RM×N and we next define the generalized Huber function SBx as
(6)SBx=infυ∈RNυ1+12Bx−υ22.Then, we propose a multivariate generalization of the MC penalty. The basic idea is to generalize (6) using the L1-norm and the generalized Huber function. Thus, we define the generalized MC (GMC) penalty function ψBx as follows:
(7)ψBx=x1−SBx.


### 2.2. The Denoising Algorithm Based on Convex Optimization

Lety∈RM be the original observed signal, A∈RM×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]:
(8)Qx=12y−Ax22+λx1.Undoubtfully, this is a typical convex optimization problem of first-order TV denoising, which is comprised by a quadratic fidelity term and an L1-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 L1-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:RN→R as follows:
(9)Fx=12y−Ax22+λψBDx,where A is an oversampled inverse discrete Fourier transform, ψB:RN→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∈Rn−1×n is defined as
(10)D=−11−11MM−11.The penaltyψB is parameterized by a matrix B and the convexity of Fx depends on B being suitably prescribed. It also should be pointed out that the choice of B will depend on A. If
(11)BTB≤1λATA,then Fx is a convex function. It is easy to satisfy the convexity condition. Given the matrix A∈RM×NN≥M, we may simply set
(12)B=γλA,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:
(13)xopt,vopt=argminx∈RNmaxv∈RNFx,v,where Fx,v=1/2y−Ax22+λDx1−λv1−γ/2ADx−v22 is the saddle function and xopt 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 Fx,vof this kind of saddle point [33]. The resulting iterative threshold algorithm uses the soft-threshold function, which is defined as
(14)softy;λ≔0y≤λy−λsignyy≥λ.The flowchart of the presented method in this paper is plotted in Figure2.Figure 2
The flowchart of the method presented in this paper.


## 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 seriesx∈RN can be simply expressed as x=x1,x2,…,xN. Then, the L1-norm and L2-norm is denoted as x1=∑nxn and x2=∑nxn21/2 respectively.The Huber functionsx can be defined as:
(1)sx=12x2,∣x∣≤1∣x∣−12,∣x∣≥1.The minimax-concave (MC) penaltyϕx is defined as:
(2)ϕx=∣x∣−12x2,∣x∣≤112,∣x∣≥1.Figures1(a) and 1(c) have illustrated the Huber function and MC penalty. Hence, the relationship between Huber function and MC penalty can be described as
(3)ϕx=x−sx.Figure 1
Different penalty functions.
(a)
Huber function
(b)
Scaled Huber function
(c)
MC penalty function
(d)
Scaled MC penalty functionDefining the scaled versions of the Huber function and MC penalty is convenient. Forb≠0, the scaled Huber function sbx is determined by
(4)sbx=12b2x2,∣x∣≤1b2∣x∣−12b2,∣x∣≥1b2.Similarly, the scaled MC penalty functionϕbx is given by
(5)ϕbx=x1−sbx.Figures1(b) and 1(d) demonstrate the scaled Huber function and scaled MC penalty, respectively. Let B∈RM×N and we next define the generalized Huber function SBx as
(6)SBx=infυ∈RNυ1+12Bx−υ22.Then, we propose a multivariate generalization of the MC penalty. The basic idea is to generalize (6) using the L1-norm and the generalized Huber function. Thus, we define the generalized MC (GMC) penalty function ψBx as follows:
(7)ψBx=x1−SBx.


## 2.2. The Denoising Algorithm Based on Convex Optimization

Lety∈RM be the original observed signal, A∈RM×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]:
(8)Qx=12y−Ax22+λx1.Undoubtfully, this is a typical convex optimization problem of first-order TV denoising, which is comprised by a quadratic fidelity term and an L1-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 L1-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:RN→R as follows:
(9)Fx=12y−Ax22+λψBDx,where A is an oversampled inverse discrete Fourier transform, ψB:RN→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∈Rn−1×n is defined as
(10)D=−11−11MM−11.The penaltyψB is parameterized by a matrix B and the convexity of Fx depends on B being suitably prescribed. It also should be pointed out that the choice of B will depend on A. If
(11)BTB≤1λATA,then Fx is a convex function. It is easy to satisfy the convexity condition. Given the matrix A∈RM×NN≥M, we may simply set
(12)B=γλA,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:
(13)xopt,vopt=argminx∈RNmaxv∈RNFx,v,where Fx,v=1/2y−Ax22+λDx1−λv1−γ/2ADx−v22 is the saddle function and xopt 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 Fx,vof this kind of saddle point [33]. The resulting iterative threshold algorithm uses the soft-threshold function, which is defined as
(14)softy;λ≔0y≤λy−λsignyy≥λ.The flowchart of the presented method in this paper is plotted in Figure2.Figure 2
The flowchart of the method presented in this paper.


## 3. 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, harmonic signal, and strong background noise components. The numerical simulation signal can be expressed as follows:(15)s1=0.5cos2πf2t+sin2πf1t,s2=0.3sin2πf3t,s3=1.5⋅randn1,1024,y=s1+s2+s3,where s3 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 f1=30 Hz, f2=200 Hz, and f3=100 Hz.Figure3 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 f2 and f3 can be inspected, and the baseline of frequency spectrum has generated a drift. Unfortunately, the phenomenon of frequency modulation, namely, f2±f1, 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.Figure 3
The synthetic signal of simulation.
(a)
The waveform in time response
(b)
The frequency spectrum analysisSubsequently, 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 coefficient of the node (3, 0). Figure4 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.Figure 4
The result provided by wavelet denoising algorithm.
(a)
The time-domain waveform by wavelet denoising
(b)
The result of frequency spectrum analysis by wavelet denoisingFigure 5
The result provided by the typical total variation (TV) denoising.
(a)
The waveform in time response by TV method
(b)
The frequency spectrum analysis by TV methodSubsequently, CEEMDAN is employed to simulation signal analysis, and the result is drawn in Figure6. 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 f1=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 f2 and f3 can be clearly inspected. Most importantly, the frequency modulation f2±f1 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.Figure 6
The result provided by CEEMDAN for simulation signal analysis.Figure 7
The result provided by the proposed 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 Table1. 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.Table 1
The proposed algorithm in signal denoising.
Noise variance
0.5
0.9
1.3
1.5
RMSE
0.1522
0.1713
0.1854
0.2026Table 2
Computational costs of four signal processing methods for simulated signal.
Method
Wavelet denoising
TV denoising
CEEMDAN
Proposed method
Computational costs
0.52
0.46
52
0.49


## 4. 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 Figure8. 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 fr=25 Hz, and sample rate of fs=48828 Hz. According to the theoretical calculation, the characteristic frequency of the outer ring fault is determined as fo=80 Hz.Figure 8
A photo of outer race fault.Table 3
Rolling element bearing parameters.
Roller diameter/mm
Pitch diameter/mm
Number of elements
Contact angle
0.235
1.245
8
0The time domain waveform and frequency spectrum analysis results of the collected vibration signals are plotted in Figure9. 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.Figure 9
The time-frequency analysis of measured fault signal.
(a)
The time response of measured bearing fault signal
(b)
The frequency response of measured bearing fault signalThen, 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. Figures10(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 fault in Figure 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.Figure 10
The result provided by wavelet denoising and TV denoising.
(a)
The result computed by wavelet denoising
(b)
The result processed by TV denoisingSimilarly, CEEMDAN is also used to analyze the experimental data, and the corresponding result is plotted in Figure11. Unfortunately, only the characteristic frequency of the outer ring fault fo and its double frequency 2fo 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.Figure 11
The result provided by CEEMDAN for experimental signal analysis.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 fo and its multiple frequencies (2fo, 3fo, 4fo, 5fo, and 6fo) had been both identified in Figure 12. Meanwhile, the rotational frequency fr 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.Figure 12
The result provided by the proposed method.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 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.


## 5. 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 L1-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.



---

*Source: 2169364-2018-07-03.xml*