Weak Fault Detection of Tapered Rolling Bearing Based on Penalty Regularization Approach

Aimed at the issue of estimating the fault component from a noisy observation, a novel detection approach based on augmented Huber non-convex penalty regularization (AHNPR) is proposed. The core objectives of the proposed method are that (1) it estimates non-zero singular values (i.e., fault component) accurately and (2) it maintains the convexity of the proposed objective cost function (OCF) by restricting the parameters of the non-convex regularization. Specifically, the AHNPR model is expressed as the L1-norm minus a generalized Huber function, which avoids the underestimation weakness of the L1-norm regularization. Furthermore, the convexity of the proposed OCF is proved via the non-diagonal characteristic of the matrix BTB, meanwhile, the non-zero singular values of the OCF is solved by the forward–backward splitting (FBS) algorithm. Last, the proposed method is validated by the simulated signal and vibration signals of tapered bearing. The results demonstrate that the proposed approach can identify weak fault information from the raw vibration signal under severe background noise, that the non-convex penalty regularization can induce sparsity of the singular values more effectively than the typical convex penalty (e.g., L1-norm fused lasso optimization (LFLO) method), and that the issue of underestimating sparse coefficients can be improved.


Introduction
Condition maintenance and health management (CMHM) of rotating machines are of great significance to prevent machinery breakdown and increase system reliability in modern industrial applications [1][2][3].Vibration signals processing represents the most commonly used technique for monitoring the condition of rotating machines.However, in practice, such vibration signals are often contaminated and submerged by severe background noises and other interfering components, which results in a very low detection accuracy.Therefore, to accurately extract the useful weak information from the measured vibration signal becomes a hotspot issue in the area of mechanical fault diagnosis.
A number of state-of-the-art methodologies including wavelet/wavelet packet transform (W-WPT) [4,5], adaptive signal decomposition (ASD) [6][7][8] and high-order cyclic spectral analysis (HOCS) [9,10], have been developed for detecting weak fault information of rotating machinery under strong noise.Recently, sparse low-rank matrix regularization (SLMR) methods, that estimate the fault signal through solving optimization problems (or inverse problems), have attracted significant attention because of their ability to induce sparsity of the fault signal (or singular values) more effectively than the traditional sparse over-complete dictionaries (SODs) such as step-impulse dictionary and unit-impulse response dictionary [11][12][13][14][15].For the SLMR, in particular, we consider the problem of estimating fault signal x from the acquired signal y, i.e., y = Ax + w, where signal w represents additive noise.To estimate the fault signal x, the optimization problem can be expressed by minimizing the objective cost function (OCF) F(x), i.e., x == argmin x { 1 2 y − Ax 2 2 + λθ(x)}, λ > 0, where θ(x) is a penalty function that induces sparsity of the approximating value of x and λ is the regularization term parameter.
Methodologies for SLMR can be roughly categorized as convex and non-convex penalty regularization approaches.Among convex penalty functions, L1-norm penalty function is commonly used to regularize linear inverse problems, i.e., penalty function θ(x) = x 1 .In that case, the OCF F(x) is trivially convex and the global optimum can be obtained.Unfortunately, the L1-norm may underestimate the sparse coefficients of the underlying signal [16,17].In contrast, numerous non-convex penalties (e.g., L P pseudo-norm with p < 1 [18], reweighted L2/L1 norm [19]) and non-convex algorithms (e.g., orthogonal matching pursuit, OMP) have been proposed [20][21][22][23][24][25], which estimate non-zero components more accurately than convex regularization.However, the previous non-convex penalties were considered based on separable penalty functions and thus suffer from several issues such as sensitivity to changes with initial data and regularization term parameters, etc., and the convexity of the OCF cannot be guaranteed simultaneously.Therefore, the core objective of SLMR is to determine how to maintain convexity of the OCF when the regularizer (or penalty function) is taken to be non-convex, so that the fault component can be estimated accurately from its noisy observation.
To address such issues, a novel non-separable non-convex regularizer, namely augmented Huber penalty function, is developed in this paper, which overcomes limitations of separable non-convex regularization and ensures the strict convexity of the OCF.Specifically, the augmented penalty function is firstly defined by L1-norm and generalized Huber functions, the strict convexity and convexity condition of the proposed OCF is then derived and provided.Meanwhile, an iterative algorithm (e.g., forward-backward splitting, FBS) is presented as a solution to the proposed objective cost function.Lastly, both the simulated signal and engineering bearing data of larger reducers are employed to validate the effectiveness of the proposed approach in terms of diagnostic accuracy and signal-to-noise ratio, etc.
The layout of this work is organized as follows: Section 2 introduces the algorithm and theoretical derivation of the augmented penalty optimization; Section 3 presents the simulation evaluation of the proposed augmented Huber non-convex penalty regularization (AHNPR) method; In Section 4, the practical diagnosis results and discussion for tapered bearing using the proposed approach with other methods are presented; and Conclusions are drawn in Section 5.

Sparse Representation
Generally, estimating a sparse component x from its noisy observation y can be expressed as where w is the additive noise.It should be noted that Equation (1) belongs to a highly underdetermined equation, i.e., ill-posed or non-deterministic polynomial-time hard (N-P hard) problem [26], and there are an infinite set of solutions.Convex and non-convex optimization approaches are commonly used to estimate a sparse component from its noisy signal.A suitable optimization problem for estimating x is where λ is regularization term parameter and which controls the sparsity of the approximating value of x.The above problem is well-known as the L1-norm minimization, and the L1-norm is used as a convex function for sparsity.If the signal x is a sparse component, i.e., most of the amplitude values in x tend to zero, then the problem in Equation ( 1) can be solved by L1-norm fused lasso optimization (LFLO) algorithm, i.e., where λ 0 and λ 1 are regularization parameters.The solution of the LFLO algorithm can be given by a soft-threshold function [27], i.e., x = Soft(tvd(y, λ 1 ), λ 0 ) where tvd(•, •) is the total variation de-noising (TVD) model [28][29][30] and the soft-threshold function is defined as follows:

The Augmented Huber Non-Convex Penalty Regularization
In this section, before the elicitation of the proposed AHNPR algorithm, the following definitions and proposition for convex analysis should be reviewed first.
The Huber equation is defined through two piecewise functions [31,32], i.e., The plot of Huber function s(x) is illustrated in Figure 1.where λ is regularization term parameter and , which controls the sparsity of the approximating value of x.The above problem is well-known as the L1-norm minimization, and the L1-norm is used as a convex function for sparsity.If the signal x is a sparse component, i.e., most of the amplitude values in x tend to zero, then the problem in Equation ( 1) can be solved by L1-norm fused lasso optimization (LFLO) algorithm, i.e., where λ0 and λ1 are regularization parameters.The solution of the LFLO algorithm can be given by a soft-threshold function [27], i.e., Soft(tvd( , ), ) where tvd(, ) ⋅ ⋅ is the total variation de-noising (TVD) model [28][29][30] and the soft-threshold function is defined as follows: , Soft( , ) 0, ,

The Augmented Huber Non-Convex Penalty Regularization
In this section, before the elicitation of the proposed AHNPR algorithm, the following definitions and proposition for convex analysis should be reviewed first.
The Huber equation is defined through two piecewise functions [31,32], i.e., The plot of Huber function s(x) is illustrated in Figure 1.Let b denote a scale and the scaled Huber function (SHF) can be defined as For b = 0, i.e., s 0 (x) = 0.For b = 0, the scaled Huber function is given by For the definition of SHF, the SHF could be rewritten through a minimization form, i.e., Proof 1.According to the SHF, i.e., s b (x) = 1 b 2 s(b 2 x), b = 0, Equation ( 9) can be rewritten as Let B denote a matrix and the generalized Huber function (GHF) can be defined as In summary, if matrix B evolves into a scale, i.e., B = b, then the GHF reduces to the SHF, If B is a matrix and B T B = I is a diagonal, then the GHF is separable, and s B (x) is the sum of the SHF s b (x), i.e., Based on the definitions above, the new non-separable non-convex penalty function can be derived using the GHFs.The new non-separable non-convex penalty (NNP) function φ(x) is defined as Therefore, according to the definition of Huber function, the function φ(x) can be expressed as where s(x) is the Huber function.Similarly, the plot of NNP function φ(x) is illustrated in Figure 2. Based on the definitions above, the new non-separable non-convex penalty function can be derived using the GHFs.The new non-separable non-convex penalty (NNP) function ( ) Therefore, according to the definition of Huber function, the function ( ) x φ can be expressed as where s(x) is the Huber function.Similarly, the plot of NNP function ( ) x φ is illustrated in Figure 3. Let b denote a scale of function and the scaled NNP function can be defined as where s b (x) is SHF.For b = 0, i.e., φ 0 (x) = |x|.For b = 0, the scaled NNP function is given by Let B denote a matrix and the NNP function φ(x) will evolve into an augmented Huber function where s B (x) is the generalized Huber function (GHF).Given a matrix B, if In contrast, if B T B is a non-diagonal matrix, the AHF φ B (x) is non-separable.Thus, the proposed AHNPR method finds the component x by solving the following optimization problem: where F(x) is the objective cost function (OCF), and φ B (x) is the AHF, which is a non-convex sparsity inducing regularizer.

Convexity Condition
In this section, the convexity of the proposed OCF will be proved to consider how to set the parameter τ and the relationship between matrix A and matrix B.
If y ∈ R M , A ∈ R M×N is a matrix, the regularization term parameter λ > 0 and the proposed OCF where A T is the complex conjugate transpose of matrix A, and B T is the complex conjugate transpose of matrix B. The last term max v∈R g(x, v) is convex due to the point-wise maximum of a set of convex functions [31].Hence, the proposed OCF For convexity condition B T B ≤ 1 λ A T A, it can be rewritten as B T B = (τ/λ)A T A, in which the parameter τ satisfies τ ≤ 1.Hence, given a matrix A, the matrix B could be set as B = √ τ/λA, 0 ≤ τ ≤ 1.The parameter τ controls the non-convexity of the AHF φ B (x).If τ = 0, then B = 0 and function φ B (x) reduces to the L1-norm.If τ = 1, then B T B = (τ/λ)A T A and the function φ B (x) are maximally non-convex.
Here, the nominal range of parameter τ is discussed by a simulated discrete signal.Given a simulated systematic signal of where the amplitudes of M 0 and M 1 are 1 and 2, respectively, the signal length is N = 100 with the frequencies f 1 = 0.25 and f 2 = 0.4.For the de-noising experiment, the Gaussian noise with σ = 1.5 is added to acquire the simulated systematic signal.The systematic signal f (t) is approximately sparse in the frequency domain, as illustrated in Figure 3, which can be represented by f = Ax, where A is an over-sampled inverse discrete Fourier transform and x is a sparse vector of Fourier coefficients.Specifically, matrix A is defined as with M = 100 and N = 256.Assuming that matrix A T is the complex conjugate transpose of matrix A and it can be calculated that B T B is a diagonal matrix.
Algorithms 2018, 11, x FOR PEER REVIEW 8 of 21 where the amplitudes of M0 and M1 are 1 and 2, respectively, the signal length is N = 100 with the frequencies f1 = 0.25 and f2 = 0.4.For the de-noising experiment, the Gaussian noise with σ = 1.5 is added to acquire the simulated systematic signal.The systematic signal f(t) is approximately sparse in the frequency domain, as illustrated in Figure 5, which can be represented by f = Ax, where A is an over-sampled inverse discrete Fourier transform and x is a sparse vector of Fourier coefficients.Specifically, matrix A is defined as ) exp( ( 2) ) − , with M = 100 and N = 256.Assuming that matrix A T is the complex conjugate transpose of matrix A and it can be calculated that B T B is a diagonal matrix.In order to obtain a nominal range of parameter 0 1 τ ≤ ≤ , we vary the parameter τ from 0.05 to 0.95 (by increments of 0.05) and calculate each root-mean-square error (RMSE) between the actual value 1 ( ) f t and its estimated value 1 ( ) f t ∧ , and also we calculate the optimized coefficients (or Fourier coefficients amplitude) at frequencies f1 = 0.25 and f2 = 0.4 for each parameter τ.The RMSE values and optimized coefficients for each parameter τ are shown in Figure 6.From Figure 6a, it can be seen that the larger value of parameter τ would get a lower value of RMSE, and the optimal parameters τ are concentrated in the 0.8 0.9 τ ≤ ≤ range.As a matter of fact, the smaller value of parameters τ would make the estimated signal noisier and less sparse.From Figure 6b, it is found that the optimized coefficients at f1 = 0.25 increase and then slightly decrease with increasing parameter τ; meanwhile, the optimized coefficients at f1 = 0.4 increase as parameter τ increases.From the point of view that the significant coefficients are not underestimated and get a clearer signal 1 ( ) f t , the optimal parameters τ are also concentrated in the 0.8 0.9 τ ≤ ≤ range.Therefore, in practice the nominal range of parameter τ is set to 0.8 0.9 τ ≤ ≤ .For the selection of the regularization parameter λ, we have where σ is the standard deviation (SD) of the additive noise, 0 β is the constant so as to maximize the signal-to-noise ratio (SNR).Here, parameters β0 and γ are typically set to be constant values, i.e., β0 = [0.5, 1], γ = [7.5,8].In practice, the standard deviation (SD) of the background noise in Equation ( 22) can be computed using the fault signal and healthy data under the same operating environment.Moreover, when the healthy data are not available or are unknown, the standard deviation (SD) of background noise can still be estimated by the following formula: In order to obtain a nominal range of parameter 0 ≤ τ ≤ 1, we vary the parameter τ from 0.05 to 0.95 (by increments of 0.05) and calculate each root-mean-square error (RMSE) between the actual value f 1 (t) and its estimated value f1 (t), and also we calculate the optimized coefficients (or Fourier coefficients amplitude) at frequencies f 1 = 0.25 and f 2 = 0.4 for each parameter τ.The RMSE values and optimized coefficients for each parameter τ are shown in Figure 4. From Figure 4a, it can be seen that the larger value of parameter τ would get a lower value of RMSE, and the optimal parameters τ are concentrated in the 0.8 ≤ τ ≤ 0.9 range.As a matter of fact, the smaller value of parameters τ would make the estimated signal noisier and less sparse.From Figure 4b, it is found that the optimized coefficients at f 1 = 0.25 increase and then slightly decrease with increasing parameter τ; meanwhile, the optimized coefficients at f 1 = 0.4 increase as parameter τ increases.From the point of view that the significant coefficients are not underestimated and get a clearer signal f 1 (t), the optimal parameters τ are also concentrated in the 0.8 ≤ τ ≤ 0.9 range.Therefore, in practice the nominal range of parameter τ is set to 0.8 ≤ τ ≤ 0.9.
For the selection of the regularization parameter λ, we have where σ is the standard deviation (SD) of the additive noise, β 0 is the constant so as to maximize the signal-to-noise ratio (SNR).Here, parameters β 0 and γ are typically set to be constant values, i.e., 5,8].In practice, the standard deviation (SD) of the background noise in Equation ( 22) can be computed using the fault signal and healthy data under the same operating environment.Moreover, when the healthy data are not available or are unknown, the standard deviation (SD) of background noise can still be estimated by the following formula: which is a traditional estimator of the noise level used for wavelet de-noising [33], where MAD(y) represents the median absolute deviation (MAD) of signal y, i.e., As an example, if we take the constant parameters to be β 0 = 0.7 and γ = 7.5, with the regularization parameter λ = 7.5 × 0.7 × 1.5 = 7.875, the solutions of LFLO and proposed AHNPR method with parameter τ = 0.8 are shown in Figure 5. Comparing the results of the LFLO and proposed AHNPR methods, the proposed AHNPR solution is sparser than the results of LFLO in frequency domain (see Figure 5c) and the LFLO underestimates significant coefficients of the underlying signal.
As an example, if we take the constant parameters to be β0 = 0.7 and γ = 7.5, with the regularization parameter λ = 7.5 × 0.7 × 1.5 = 7.875, the solutions of LFLO and proposed AHNPR method with parameter τ = 0.8 are shown in Figure 7. Comparing the results of the LFLO and proposed AHNPR methods, the proposed AHNPR solution is sparser than the results of LFLO in frequency domain (see Figure 7c) and the LFLO underestimates significant coefficients of the underlying signal.
As an example, if we take the constant parameters to be β0 = 0.7 and γ = 7.5, with the regularization parameter λ = 7.5 × 0.7 × 1.5 = 7.875, the solutions of LFLO and proposed AHNPR method with parameter τ = 0.8 are shown in Figure 7. Comparing the results of the LFLO and proposed AHNPR methods, the proposed AHNPR solution is sparser than the results of LFLO in frequency domain (see Figure 7c) and the LFLO underestimates significant coefficients of the underlying signal.Therefore, in practice, based on the above analysis, the optimal parameter τ is set to be 0.8 ≤ τ ≤ 0.9, the constant parameters β 0 and γ are set to be β 0 = 0.7 and γ = 7.5, and the regularization parameter λ is set to be λ = γ × β 0 × σ.
Finally, in order to minimize the F(x), and assuming B T B = (τ/λ)A T A with parameter 0.8 ≤ τ ≤ 0.9, the optimization problem could be transformed into saddle-point problem, i.e., The solution of the saddle-point problem can be obtained using the forward-backward splitting (FBS) algorithm, [34], as listed in Algorithm 1.
Algorithm 1 Iterative algorithm for the proposed AHNPR method Initialization:

Numerical Simulation
A simulation signal is utilized to investigate the effectiveness of the proposed AHNPR approach.Considering the periodical impulse features that are generated by the localized fault in rotating machinery, the simulation signal is constructed as where a = 0.1 is a damping coefficient, the structural frequency of system is f n = 2000 Hz, signal length is N = 2048, and the sampling frequency is f s = 20 KHz.In addition, the additive Gaussian noise with σ = 0.4 is added to acquire the simulated synthetic signal.Figure 6a,c depict the obtained periodical impulses with localized faults and the simulated synthetic signal.Figure 6b,d  Therefore, in practice, based on the above analysis, the optimal parameter τ is set to be 0.8 0.9 τ ≤ ≤ , the constant parameters β0 and γ are set to be β0 = 0.7 and γ = 7.5, and the regularization parameter λ is set to be Finally, in order to minimize the F(x), and assuming

B B
A A with parameter 0.8 ≤ τ ≤ 0.9, the optimization problem could be transformed into saddle-point problem, i.e., ( , ) arg min max ( , ) The solution of the saddle-point problem can be obtained using the forward-backward splitting (FBS) algorithm, [34], as listed in Algorithm 1.

Numerical Simulation
A simulation signal is utilized to investigate the effectiveness of the proposed AHNPR approach.Considering the periodical impulse features that are generated by the localized fault in rotating machinery, the simulation signal is constructed as where a = 0.1 is a damping coefficient, the structural frequency of system is fn = 2000 Hz, signal length is N = 2048, and the sampling frequency is fs = 20 KHz.In addition, the additive Gaussian noise with σ = 0.4 is added to acquire the simulated synthetic signal.Figure 8a,c depict the obtained periodical impulses with localized faults and the simulated synthetic signal.Following this, the proposed AHNPR method is applied to process the simulated synthetic signal.Specifically, matrix A is defined as , ( 1) exp( ( 2) ) − .Assuming that matrix A T is the complex conjugate transpose of matrix A, it can be calculated that A T A is a diagonal matrix, the parameter τ is set to τ = 0.8, and thus matrix B is = B A τ λ .In addition, taking the constant parameters to be β0 = 0.7 and γ = 7.5, the parameter λ could be calculated by Equation ( 22), i.e., λ = 7.5 × 0.7 × 0.4 = 2.1.Since B T B is a nondiagonal matrix, the augmented Huber penalty function The de-noising results of the LFLO method and proposed AHNPR method are shown in Figure 9a,c, respectively.Figure 9b,d  Following this, the proposed AHNPR method is applied to process the simulated synthetic signal.Specifically, matrix A is defined as Assuming that matrix A T is the complex conjugate transpose of matrix A, it can be calculated that A T A is a diagonal matrix, the parameter τ is set to τ = 0.8, and thus matrix B is B = √ τ/λA.In addition, taking the constant parameters to be β 0 = 0.7 and γ = 7.5, the parameter λ could be calculated by Equation ( 22), i.e., λ = 7.5 × 0.7 × 0.4 = 2.1.Since B T B is a non-diagonal matrix, the augmented Huber penalty function φ B (x) is non-separable.
The de-noising results of the LFLO method and proposed AHNPR method are shown in Figure 7a,c, respectively.Figure 7b,d are the corresponding Morlet wavelet time-frequency diagrams of Figure 7a,c, respectively.It can be noted that the periodic characteristics in Figure 7c are more obvious than in Figure 7a.The sparse Fourier coefficients of LFLO and the proposed AHNPR method are shown in Figure 8a,b.The non-convex penalty function method aims to induce sparsity in sparse coefficients more effectively than the LFLO method.Moreover, the LFLO method underestimates the significant (large-amplitude) sparse coefficients of the de-noised signal.
coefficients more effectively than the LFLO method.Moreover, the LFLO method underestimates the significant (large-amplitude) sparse coefficients of the de-noised signal.

Experimental Verification
To demonstrate the validity of the proposed approach in engineering applications, the tapered rolling bearing with a faulty outer race is implemented.Experimental vibration data were collected from several accelerometers instrumented on a bearing end bracket from a large reducer, as shown in Figure 11.The geometrical parameters of the tested tapered bearing are listed in Table 1, and the fault frequency of the bearing outer race is 118.8Hz.The sampling frequency is 51.2 KHz (right side of input shaft) and 5.12 KHz (left side of input shaft), the rotation frequency is 16.67 Hz and sampling length is 6.4 s.The tapered rolling bearing (FAG-32212-A) fault was located on the right side of the input shaft.In this experiment, the vibration data collected from the right side of input shaft (horizontal direction, i.e., non-drive end) were used in the first experiment.To create a more challenging condition, the vibration data collected from the left side of the input shaft (horizontal direction, drive end) were also analyzed in the second experiment.

Experimental Verification
To demonstrate the validity of the proposed approach in engineering applications, the tapered rolling bearing with a faulty outer race is implemented.Experimental vibration data were collected from several accelerometers instrumented on a bearing end bracket from a large reducer, as shown in Figure 9.The geometrical parameters of the tested tapered bearing are listed in Table 1, and the fault frequency of the bearing outer race is 118.8Hz.The sampling frequency is 51.2 KHz (right side of input shaft) and 5.12 KHz (left side of input shaft), the rotation frequency is 16.67 Hz and sampling length is 6.4 s.The tapered rolling bearing (FAG-32212-A) fault was located on the right side of the input shaft.In this experiment, the vibration data collected from the right side of input shaft (horizontal direction, i.e., non-drive end) were used in the first experiment.To create a more challenging condition, the vibration data collected from the left side of the input shaft (horizontal direction, drive end) were also analyzed in the second experiment.The first experiment had two tasks, i.e., to estimate the weak fault signal from the original vibration signal and to calculate the sparsity of the sparse coefficient.The raw vibration signal (right side of input shaft, 4096 sampling points were selected, i.e., 0.08 s) and its envelop spectrum are displayed in Figure 12a,b, respectively.As shown in Figure 12b, although the spectrum peak at 237.5 Hz, which is consistent with the twice outer-race fault frequency, can be detected; however, the spectrum peak is masked by heavy background noise and frequency features are not clear enough to detect practical fault.The first experiment had two tasks, i.e., to estimate the weak fault signal from the original vibration signal and to calculate the sparsity of the sparse coefficient.The raw vibration signal (right side of input shaft, 4096 sampling points were selected, i.e., 0.08 s) and its envelop spectrum are displayed in Figure 10a,b, respectively.As shown in Figure 10b, although the spectrum peak at 237.5 Hz, which is consistent with the twice outer-race fault frequency, can be detected; however, the spectrum peak is masked by heavy background noise and frequency features are not clear enough to detect practical fault.The first experiment had two tasks, i.e., to estimate the weak fault signal from the original vibration signal and to calculate the sparsity of the sparse coefficient.The raw vibration signal (right side of input shaft, 4096 sampling points were selected, i.e., 0.08 s) and its envelop spectrum are displayed in Figure 12a,b, respectively.As shown in Figure 12b, although the spectrum peak at 237.5 Hz, which is consistent with the twice outer-race fault frequency, can be detected; however, the spectrum peak is masked by heavy background noise and frequency features are not clear enough to detect practical fault.The proposed AHNPR approach and the LFLO method are respectively employed to process the measured vibration signal and the obtained results are displayed in Figure 11.The parameter τ was set to τ = 0.8 and matrix A is also defined as It can be calculated that B T B is a non-diagonal matrix, thus the augmented Huber function (AHF) φ B (x) is non-separable.Here, it should be noted that the healthy data are not available or are unknown, according to Equations ( 23) and (24).Taking the constant parameters to be β 0 = 0.7 and γ = 7.5, the parameters of the proposed method can be obtained as follows: the standard deviation σ = 4.234, and regularization term parameters λ = 7.5 × 0.7 × 4.234 = 22.23.Comparing results with the time-domain waveform and envelope spectrum, one can observe improvements at two levels: the first is the significant enhancement in the contribution of periodical characteristic.The periodical characteristic with LFLO is weaker than that related to the proposed method.The second improvement is that the amplitudes of fault frequency and its harmonic are worthy of further enhancement.Note that some other unrelated components in Figure 11d are greatly reduced compared with the envelope spectra of the raw signal and the de-noised signal generated by the LFLO method.23) and (24).Taking the constant parameters to be β0 = 0.7 and γ = 7.5, the parameters of the proposed method can be obtained as follows: the standard deviation σ = 4.234, and regularization term parameters λ = 7.5 × 0.7 × 4.234 = 22.23.Comparing results with the time-domain waveform and envelope spectrum, one can observe improvements at two levels: The first is the significant enhancement in the contribution of periodical characteristic.The periodical characteristic with LFLO is weaker than that related to the proposed method.The second improvement is that the amplitudes of fault frequency and its harmonic are worthy of further enhancement.Note that some other unrelated components in Figure 13d are greatly reduced compared with the envelope spectra of the raw signal and the de-noised signal generated by the LFLO method.In addition, the distribution of sparse coefficients generated by LFLO and the proposed AHNPR method are given in Figure 14a,b, respectively.One can see that the sparse coefficients are underestimated by the LFLO method along with the sparse coefficients associated with additive noise that reduce sparsity.In comparison with the LFLO result, the sparse coefficients of the proposed AHNPR method are enhanced and amplified.In addition, the distribution of sparse coefficients generated by LFLO and the proposed AHNPR method are given in Figure 12a,b, respectively.One can see that the sparse coefficients are underestimated by the LFLO method along with the sparse coefficients associated with additive noise that reduce sparsity.In comparison with the LFLO result, the sparse coefficients of the proposed AHNPR method are enhanced and amplified.In addition, the distribution of sparse coefficients generated by LFLO and the proposed AHNPR method are given in Figure 14a,b, respectively.One can see that the sparse coefficients are underestimated by the LFLO method along with the sparse coefficients associated with additive noise that reduce sparsity.In comparison with the LFLO result, the sparse coefficients of the proposed AHNPR method are enhanced and amplified.As to the second experiment, Figure 13a,b present the raw vibration signal on the drive end (i.e., the left side of the input shaft, 4096 sampling points were selected, i.e., 0.8 s) and its envelope spectrum.In this case, the outer race fault frequency cannot be observed.This may be caused by the fact that the outer race fault is quite weak and the transfer path of the accelerometers located on the left side are longer than in the previous experiment.As such, the proposed AHNPR algorithm is applied to separate the fault feature from other unrelated components and the results are displayed in Figure 13c,d.The parameters of the proposed method can be obtained as follows: the standard deviation σ = 4.44, regularization term parameter λ = 7.5 × 0.7 × 4.4 = 23.31, and parameter τ is set to τ = 0.8.
The outer race fault frequency and its harmonics are detected in the corresponding envelope spectrum.In addition, interestingly, the sidebands are distributed in both fault frequency and its harmonics, after calculating the difference between the fault frequencies and sidebands.It should be noted that the difference values are consistent with the rotational frequency (16.67 Hz).Therefore, the outer race fault of bearings can be more clearly diagnosed.The above result demonstrates that the proposed algorithm can detect the bearing fault and eliminate the background noise interference, regardless of whether the accelerometers are located on the right or left side of the input shaft.As to the second experiment, Figure 15a,b present the raw vibration signal on the drive end (i.e., the left side of the input shaft, 4096 sampling points were selected, i.e., 0.8 s) and its envelope spectrum.In this case, the outer race fault frequency cannot be observed.This may be caused by the fact that the outer race fault is quite weak and the transfer path of the accelerometers located on the left side are longer than in the previous experiment.As such, the proposed AHNPR algorithm is applied to separate the fault feature from other unrelated components and the results are displayed in Figure 15c,d.The parameters of the proposed method can be obtained as follows: the standard deviation σ = 4.44, regularization term parameter λ = 7.5 × 0.7 × 4.4 = 23.31, and parameter τ is set to τ = 0.8.
The outer race fault frequency and its harmonics are detected in the corresponding envelope spectrum.In addition, interestingly, the sidebands are distributed in both fault frequency and its harmonics, after calculating the difference between the fault frequencies and sidebands.It should be noted that the difference values are consistent with the rotational frequency (16.67 Hz).Therefore, the outer race fault of bearings can be more clearly diagnosed.The above result demonstrates that the proposed algorithm can detect the bearing fault and eliminate the background noise interference, regardless of whether the accelerometers are located on the right or left side of the input shaft.

Conclusions
This paper proposes a novel weak fault detection approach for tapered rolling bearing based on the augmented Huber nonconvex penalty regularization (AHNPR) approach.The purpose of this paper is to overcome the following two problems: (1) The traditional convex regularizer (e.g., L1norm) may underestimate sparse coefficients of the underlying signal; and (2) the convexity of the objective function cannot be guaranteed using traditional separable non-convex penalties regularizers.
In this work, we proposed an AHNPR for solving the low-rank matrix approximation (LRMA) problem (or inverse problem).Compared to the LFLO algorithm, where the singular values would be calculated by a soft-threshold function, the proposed method does not underestimate singular values.The AHNPR is known to estimate singular values more accurately.In addition, an efficient algorithm using the forward-backward splitting (FBS) technique to solve the proposed OCF is derived.Finally, the effectiveness of the proposed method is demonstrated by the simulated vibration signal practical rolling bearing under severe additive background noise.
A possible future direction of research involves the use of AHFs, which are focused on multifault diagnosis under variable speed or variable harsh running environments.

Figure 1 .Figure 1 .
Figure 1.The Huber function s(x).Let b denote a scale and the scaled Huber function (SHF) can be defined as

Algorithms 2018, 11 ,Figure 2 .
Figure 2. The generalized Huber function (GHF) with different matrices B. (a) The B T B is a diagonal matrix; (b) the B T B is a non-diagonal matrix.

Figure 3 .Figure 2 .
Figure 3.The non-separable non-convex penalty function.Let b denote a scale of function and the scaled NNP function can be defined as ( ) ( ) b b x x s x φ

Figure 5 .
Figure 5.The noise-free signal, simulated systematic signal and its frequency spectrum: (a) Noise-free signal; (b) simulated systematic signal; (c) frequency spectrum of simulated systematic signal.

Figure 3 .
Figure 3.The noise-free signal, simulated systematic signal and its frequency spectrum: (a) Noise-free signal; (b) simulated systematic signal; (c) frequency spectrum of simulated systematic signal.

Figure 9 .Figure 7 .Figure 10 .
Figure 9.The solutions of the LFLO and the proposed AHNPR method.(a) The result of the LFLO method; (b) Morlet wavelet time-frequency diagram of the result generated by LFLO; (c) the result of the proposed AHNPR method; (d) Morlet wavelet time-frequency diagram of the result generated by the proposed AHNPR method.

Figure 8 .
Figure 8.The results of sparse Fourier coefficients using the LFLO and the proposed AHNPR method.(a) The result of the LFLO method; (b) the result of the proposed AHNPR method.

Figure 11 .
Figure 11.Experimental set up for the tapered bearing test.

Figure 12 .
Figure 12.Original vibration signal and its envelope spectrum.(a) Original vibration signal; (b) the envelope spectrum of the original vibration signal.

Figure 9 .
Figure 9. Experimental set up for the tapered bearing test.

Algorithms 2018 , 21 Figure 11 .
Figure 11.Experimental set up for the tapered bearing test.

Figure 12 .Figure 10 .
Figure 12.Original vibration signal and its envelope spectrum.(a) Original vibration signal; (b) the envelope spectrum of the original vibration signal.

Figure 13 .
Figure 13.The results of LFLO and the proposed AHNPR method.(a) The time-domain waveform generated by LFLO; (b) envelope spectrum of the results based on the LFLO method; (c) the timedomain waveform generated by the proposed AHNPR method; (d) envelope spectrum of the results based on the proposed AHNPR method.

Figure 14 .Figure 11 .
Figure 14.The optimized coefficients of the de-noised signal generated by LFLO and the proposed AHNPR method.(a) The result using the LFLO method; (b) the result using the proposed AHNPR method.

Figure 13 .
Figure 13.The results of LFLO and the proposed AHNPR method.(a) The time-domain waveform generated by LFLO; (b) envelope spectrum of the results based on the LFLO method; (c) the timedomain waveform generated by the proposed AHNPR method; (d) envelope spectrum of the results based on the proposed AHNPR method.

Figure 14 .Figure 12 .
Figure 14.The optimized coefficients of the de-noised signal generated by LFLO and the proposed AHNPR method.(a) The result using the LFLO method; (b) the result using the proposed AHNPR method.

Figure 16 .
Figure 16.The diagnosis results for the right and and left side of the input shaft generated by the asymmetric convex penalty regularization (ACPR) method.(a) The sparse component of the right side of the input shaft, demodulated by the ACPR method; (b) the sparse component of the left side of the input shaft, demodulated by the ACPR method; (c) the envelope spectrum of the sparse component of the right side of the input shaft; (d) envelope spectrum of the sparse component of the left side of the input shaft.
are the corresponding Morlet wavelet time-frequency diagrams of periodical impulses and simulated synthetic signals, respectively.From Figure6c,d, it is seen that the periodic impulses are completely buried in heavy background noise.

Table 1 .
The geometrical parameters of the tested tapered rolling bearing.

Table 1 .
The geometrical parameters of the tested tapered rolling bearing.
Algorithms 2018, 11, x FOR PEER REVIEW 15 of 21The proposed AHNPR approach and the LFLO method are respectively employed to process the measured vibration signal and the obtained results are displayed in Figure13.The parameter τ is non-separable.Here, it should be noted that the healthy data are not available or are unknown, according to Equations (