Nonconvex Wavelet Thresholding Total Variation Denoising Method for Planetary Gearbox Fault Diagnosis

The vibration-based analysis is an effective technique for planetary gearbox fault diagnosis, but a difficult task is how to accurately identify fault features from noisy vibration signals. In this paper, a nonconvex wavelet thresholding total variation (WATV) denoising method is proposed for planetary gearbox fault diagnosis, which combines wavelet-domain sparsity and total variation (TV) regularization. The TV regularization algorithm is employed to modify the retained wavelet coefficients so that the occurrence of oscillations caused by wavelet thresholding is suppressed. To overcome the underestimation shortcoming of L1-norm regularization, nonconvex penalty function regularization is used to strongly promote the sparsity of estimation while guaranteeing that the global optimal solutions are obtained even though the objective function is nonconvex. Then, the split augmented Lagrangian shrinkage (SALSA) method is developed to solve the nonconvex WATV denoising problem. Two experimental studies are performed to verify the performance and effectiveness of the proposed method. Comparisons with the soft thresholding and basis pursuit denoising (BPD) methods show that the proposed method can accurately estimate the fault features from vibration signals, which means that the proposed method is an effective and promising tool for planetary gearbox fault diagnosis.


I. INTRODUCTION
Due to their small size, lightweight, large transmission ratio and strong load capacity, planetary transmission gearboxes are widely used in military aircraft, new armored vehicles, self-propelled artillery and other military equipment and civilian equipment. Due to the harsh operating conditions and intensive impact load of planetary gearboxes, gear faults such as cracks and missing teeth frequently occur [1], [2]. Failure of planetary gearbox seriously affects the safety and reliability of the equipment. If the fault is not diagnosed in time, it will cause secondary damage to the equipment and result in major economic losses and even human casualties. Thus, fault diagnosis [3] and fault trend forecasting [56] have attracted considerable attention in the past decades.
The associate editor coordinating the review of this manuscript and approving it for publication was Hao Luo .
Vibration-based analysis is an effective method for diagnosing mechanical faults [3]- [6], because the vibration signals generated by mechanical faults indicate the machinery operation condition. The information contained in the vibration signal describes not only the machine health condition but also the severity. However, in a planetary gearbox, the transmission system is more complicated than a fixed-axis gear, so the transmission paths from gear meshing points to transducers are multiple and time-varying, which may deteriorate or attenuate the vibration response of faulty components [3], [7]. In addition, the vibration signals are usually degraded by strong background noise from the working environment and other machine components. Thus, how to accurately extract the fault features from a complex noisy signal is a difficult task for planetary gearbox fault diagnosis.
Various signal processing techniques have been widely applied to vibration signal analysis and fault diagnosis during the past decades, such as time-frequency analysis [8]- [10], VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ the spectral kurtosis (SK) [11]- [13], the empirical mode decomposition (EMD) [14]- [16], the local mean decomposition (LMD) [17] and the wavelet transform (WT) [18]- [20]. The multiresolution analysis ability of the WT makes it suitable for extracting fault features from nonstationary vibration signals. Wang et al. [21] proposed a method combining the Morlet wavelet and correlation filtering to identify the impulse response parameters and the cyclic period between adjacent impulses, which is an effective way to extract features of gearbox fault diagnosis. In [22], the wavelet-based multiscale slope features were extracted from the slope of logarithmic variances calculated from the wavelet coefficients of the discrete WT, and then the extraction features were used to classify gearbox faults with high accuracy and stability.
He et al. [23] presented a wavelet packet base-selection method that selected an optimal set of time-frequency subspaces to produce discriminant features to enhance the accuracy of gearbox fault diagnosis. Wang et al. [24] constructed a new multiwavelet via an adaptive lifting scheme to extract fault features from vibration signals of a gearbox. To overcome the disadvantage that the common threshold denoising method misses important but weak features in gear fault diagnosis, Yuan et al. [25] developed a novel method combining customized multiwavelet lifting schemes with sliding window denoising. In [26], a series of sets of wavelet packet coefficients on various frequency bands were taken as an input of the deep residual networks to improve the performance of planetary gearbox fault diagnosis. Unfortunately, most of the previous studies have chosen standard and fixed wavelet basis functions independent of the given signal. In practice, the vibration signals are the dynamic integrated response of running systems, fault components, transmission paths and so on. Even the same faults generate various dynamic response signals in different transmission paths. Thus, fixed wavelet basis functions independent of the input dynamic response signals tend to reduce the accuracy of fault diagnosis [27]. On the other hand, universal thresholds are usually set in previous denoising methods, which are the same for the decomposition coefficients in the same layer. The universal thresholds often obscure the decomposition coefficients and may lead to the loss of some critical but relatively weak information in the fault feature detection, possibly reducing the accuracy of fault diagnosis [28].
Recently, a new branch of signal processing method, sparse representation has received considerable attentions in the field of mechanical fault diagnosis. The application of sparse representation in mechanical fault diagnosis was initially studied in [29]. To extract fault features from gearbox vibration signals, Li et al. [30] proposed a novel multiple enhanced sparse decomposition (MESD) method to address multiple feature extraction for gearbox compound fault vibration signals. In [31], based on the oscillatory behavior of the vibration signal, a sparsity-enabled signal decomposition method using L1-norm regularization was proposed for fault feature extraction of gearboxes. For multi-fault diagnosis of gearboxes, Zhang et al. [32] introduced a resonance-based sparse signal decomposition method with a comb filter. Yu et al. [33] proposed an improved morphological component analysis method which is a sparsity-based decomposition method, to diagnosis compound faults in a gearbox. The above sparse representation methods have satisfactory results for feature identification in gearbox fault diagnosis. In these studies, the classical L1-norm regularization is used to regularize the sparse representation problem, because the L1-norm induces the sparsity of the estimation more effectively than other convex penalties [34]. However, L1-norm solutions are not ideal for planetary gearbox fault diagnosis, because L1-norm regularization often underestimates the high-amplitude components of interest, which may still make the weak fault components obscured in noise and hard to be effectively identified [34]- [36].
To overcome the underestimation characteristic and enhance the sparsity of the estimation, many nonconvex sparse regularization methods that replace the classical L1-norm regularization with nonconvex sparsity-inducing penalties have been developed and used in mechanical fault diagnosis. To improve the decomposition accuracy for gearbox fault diagnosis, Cai et al. [35] introduced an improved sparsity-enabled signal decomposition method which used the generalized minimax-concave penalty function as a nonconvex regularizer to enhance the sparsity in the sparse approximation. In the STFT domain, Ding et al. [37] used a nonconvex penalty function to promote the sparsity of group STFT domain coefficients in an optimization problem, allowing the periodically oscillatory fault features of rotating machinery to be effectively extracted. In [36], to address the fault feature enhancement for wind turbine gearbox vibration signals, Wang et al. proposed a dual-enhanced sparse decomposition method in which the nonconvex generalized minimax-concave penalty was used to construct the sparse regularized cost function. He et al. [38] extended the overlapping group sparsity to the nonconvex regularization problem, the proposed method used a nonconvex penalty function to model the periodicity of the sparse groups, making this method suitable for feature extraction in machinery fault diagnosis.
The sparse enhancement properties of the nonconvex penalty function provide a new insight on how to accurately extract fault features. In this paper, we propose a nonconvex wavelet thresholding total variation denoising method for planetary gearbox fault diagnosis, which combines waveletdomain sparsity and total variation (WATV) regularization. The proposed method employs the TV regularization algorithm to modify the retained wavelet coefficients so that the restoration process does not degrade the information that has been considered as significant in the denoising step. Moreover, to strongly promote the sparsity of estimation, a nonconvex penalty function is employed as the regularizer of the sparse representation, and this function is guaranteed to obtain the global optimal solution even though the objective function is nonconvex. We present two experimental studies to verify the effectiveness of the proposed method in the diagnosis of localized faults in a planetary gearbox. Comparisons show that the proposed method can significantly improve the accuracy of planetary gearbox fault diagnosis compared with the results of the soft thresholding denoising and basis pursuit denoising (BPD) method.
The paper is organized as follows. Section 2 illustrates the basic preliminaries. Section 3 presents the non-convex WATV denoising model, the nonconvex penalty functions and convexity condition, and the split augmented Lagrangian shrinkage (SALSA) method is derived to solve the nonconvex WATV denoising problem. In Section 4, two experiment studies are performed to verify the effectiveness of the proposed method. Finally, conclusions are summarized in Section 5.

II. PRELIMINARIES
Vibration signals observed from a faulty planetary gearbox can be modeled as where y (t) ∈ R N is the measured signal, which usually contains noise and fault signals; x (t) ∈ R N is the fault signal, which contains impulsive components with periodic characteristics caused by localized gearbox faults; and w (t) ∈ R N is white Gaussian noise. The challenging problem in the fault diagnosis of a planetary gearbox is how to accurately restore the fault signal x (t) from the noisy measured signal y (t).

A. WAVELET THRESHOLDING DENOISING
Generally, wavelet denoising contains three steps [38], [39]: forward transformation of the noisy signal to an orthogonal domain, reduction of the wavelet coefficients smaller than a given threshold, and inverse transformation of the data to the original domain. Let φ j,k (t) be an orthogonal wavelet basis function, so that the standard WT of signal y (t) can be written as where a = a j,k are the corresponding wavelet coefficients defined as Because of its simplicity, the wavelet thresholding technique has been widely used in engineering practice since the beginning of wavelet use in signal processing. One way of describing wavelet thresholding technology is to establish a set M = (j, k) ∈ K : a j,k ≥ η to record the indexes of the wavelet coefficients, then preserve all coefficients whose indexes belong to M and truncate the others to zero where τ (·) is a thresholding operator. For example, if τ is a hard thresholding operator, M is defined as the set of all coefficients whose values are larger than a given threshold. If τ is a linear thresholding operator, then M is taken as the index set of all low frequencies.

B. TOTAL VARIATION DENOISING
In signal denoising problems, estimation of a noise free signal x (t) from an observed noisy signal y (t) is an inverse problem as well as a classical ill-posed problem [40]. A standard and efficient way to deal with inverse problems is to define a suitable objective function F (x) consisting of a data term (consistency with the measurements) and a regularization term (based on prior information), and then to find the signal minimizing F (x) [41]. Good denoising results should provide a balance between the regularization term and the data term and can be achieved only with some form of regularization penalizing undesirable solutions. Accordingly, typical criteria have the form where y − x 2 2 is a data fidelity term used to measure the difference between y and x, φ (x) is a penalty function (or regularization term) that should be chosen as to penalize unwanted behavior in x, and λ > 0 is the regularization parameter, used to balance the trade-off between two terms.
If the signal of interest is known to be sparse, this prior information can be used to come up with a sparse regularization method. In this case, the penalty function is defined to measure the number of nonzero values in x, i.e., φ (x) = x 0 where x 0 is the l 0 -norm and defined as x 0 = N n=1 |x n | 0 . Unfortunately, with φ (x) defined as such, the regularization problem in (5) is an NP-hard problem for which the objective function F (x) cannot readily to minimize [42], [43]. Therefore, it is common to replace the l 0 -norm by the l 1 -norm as the penalty function in practical applications, because it induces sparsity most effectively and does not sacrifice the convexity of the objective function. The l 1 -norm regularization problem is given bŷ (6) is also known as BPD [41], and has been widely used in image denoising, fault diagnosis, and ultrasonic signal processing.
Another widely used nonlinear filter method based on sparse model is TV [44], [45], which operates under the assumption that the derivative of the underlying signal is sparse. It is defined as a convex optimization problem of minimizing the cost function comprising a nondifferentiable convex penalty term and a quadratic data fidelity term. The TV in the signal x ∈ R N is defined as Then, the TV denoising can be written aŝ Instead of the sparsity of the solution ofx, TV denoising reduces the sparsity of the first-order difference of the solution [46], [47]. TV denoising exploits the combination of the l 1 -norm with a derivation operator, which allows TV regularization to preserve the sharp edge information and hence offers better approximation quality for the given signals. It has been shown that TV denoising can remove noise while preserving edges, which is very important to image processing [48]. To deal with more common types of signals, TV denoising is also used in combination with other methods.

III. NONCONVEX TOTAL VARIATION WAVELET THRESHOLDING DENOISING METHOD A. SPARSE FAULT DIAGNOSIS MODEL
Although the TV denoising method can improve the approximation signal, it cannot be directly applied to the fault signal of a planetary gearbox in practical engineering, because the theoretical basis of the TV denoising method is based on the sparse prior information of the signal, while the fault signal of the planetary gearbox is also with nonsparse components. Recent studies have demonstrated that fault signals can be sparsely represented in an appropriate WT framework [34], [49], i.e. a = Wx and x = W T a (9) where W and W t are WT and inverse WT, respectively, and a represents the wavelet coefficients with sparse structure Under this sparse prior, the sparse theoretical model for fault signal extraction of a planetary gearbox can be obtained by combining the regularization theory shown in (6) and the TV noise reduction method shown in (8).
where λ j and β are regularization parameters that express the same meaning as λ. Then, the estimation of the fault signal x can be obtained by the inverse WT, i.e.,x = W Tâ .

B. NONCONVEX PENALTY FUNCTION
It is common to use the l 1 -norm for the penalty function φ in the sparse fault diagnosis model (shown in (10)), because the l 1 -norm can effectively promote the sparsity of the estimation while ensuring that the objective function F is a convex function [50]. Convex functions are attractive because a wealth of convex optimization theory can be used to guarantee the solution is the global optimal solution. However, the l 1 -norm is not the tightest convex envelope of sparsity, which leads to the nonzero values of the underlying signal are underestimated [51], [52], but in most cases these values constitute the signal of interest. Thus, the l 1 -norm solutions are not ideal for planetary gearbox diagnosis.
To more accurately estimate the nonzero values, nonconvex regularizers are often favored over the l 1 -norm, because nonconvex penalty functions can lead to sparser solutions than l 1 -norm. However, the use of nonconvex regularizers leads to the objective function F being nonconvex, and the signal recovery problem becoming a nonconvex optimization problem. Consequently, spurious local minima exist in which optimization algorithms may be trapped, and convergence only to a local minimum is guaranteed. In addition, solutions to the nonconvex problem are highly sensitive to perturbation of the input signal: a small perturbation in the input may lead to a large change in the output [50]- [52]. Fortunately, Nikolova [54] proposed an idea of specifying a nonconvex penalty in the formulation of convex optimization to overcome the fundamental limitation of nonconvex penalties. By restricting the nonconvex degree parameter of the nonconvex penalties to control the nonconvexity of the regularizer, the objective function can be guaranteed to be strictly convex. This idea has been applied to image processing and fault diagnosis.
In this paper, the nonconvex sparsity-inducing penalty function is used as the penalty function in (10). We use the notation φ (x; b) to denote the nonconvex penalty function, and b is a scalar parameter that controls the degree of nonconvexity of φ (x; b), with b ≥ 0. Under the assumption that such penalty function satisfies the following conditions [51], [53]: 1. φ is continuous, increasing and concave on R Several typical examples satisfying the above-listed conditions are the rational penalty function, the logarithmic penalty function and the arctangent function, and these examples are shown in Fig. 1. Obviously, for the same parameter b = 0.4, the arctangent function shows the most concavity, which means the arctangent function reduces the sparsity most strongly among the three given nonconvex functions. Thus, we use the arctangent function as the sparsity-inducing penalty function φ (x; b)

C. CONVEXITY CONDITION
In this section, we try to find a suitable condition on the parameter to ensure that the objective function is strictly convex. The objective function in (10) can be written as where Note that if F 1 and F 2 are strictly convex, then F is strictly convex. F 2 is defined as the l 1 -norm, so it is a convex function. Hence, it suffices to find the parameter b such that F 1 is strictly convex.
Combining the propositions 1 and 2 in [53], we note the following proposition.
Proposition 1: Let φ (x; b) be a nonconvex penalty that satisfies the above-listed properties, with b ≥ 0. The function G (x) defined as is strictly convex if Proof: By the Lemma A in [53], to ensure that G (x) is strictly convex, the second derivative of G (x) must be positive on By the property 5 in the part B of the part III, G 0 + = −y + λ; by symmetry G 0 − = −y − λ and λ > 0, G 0 − < G 0 + . To ensure that G (x) is positive on R/ {0}, we have the following condition: By the property 6 in the part B of the part III, we have b > 1/λ. Hence, F 1 is strictly convex if and only if Then, F 1 being a sum of strictly convex G is strictly convex. It follows that the objective function F in (10), being the sum of F 1 (strictly convex) and F 2 (convex), is strictly convex, with the condition

D. SALSA ALGORITHM FOR NONCONVEX WATV
With the convexity condition (21), we can reliably obtain via convex optimization the global minimum of (10) as long as the parameter b is chosen to satisfy (21). In this paper, the split augmented Lagrangian shrinkage (SALSA) [55] method is used to solve the nonconvex WATV denoising problem. The good convergence of SALSA has been proven in practice [42], [55]. With variable splitting, problem (10) can be transformed into the constrained optimization problem min a,c where Using the augmented Lagrangian method to represent the problem (24), we have where d is a multiplier vector to constraint a = c. The solution of (25) can be obtained by iteratively minimizing with respect to each variable alternately, as proven in [55]. Each iteration step of the SALSA algorithm is given by Substituting f 1 (a) and f 2 (c) into (26) and (27) respectively, the explicit form of the SALSA algorithm can be obtained as follows: Table 1 summarizes the whole algorithm for solving problem (10), which is guaranteed to converge to the unique VOLUME 8, 2020 global minimizer. Both minimization problems, namely, iteratively thresholding and TV denoising, can be solved exactly. By running SALSA until the stopping criterion (the relative variation of the objection function falls below 10 −5 ) is satisfied, the global optimal solution of the nonconvex WATV denoising problem can be found, as shown in Fig. 2.   Fig. 3. The test rig is composed mainly of a driving motor, transmission box, clutch, planetary gearbox, loading motor, hydraulic station, and speed and torque instrument and test system. The driving motor transmits power to the planetary gearbox through the gearbox and clutch, and then the power is transferred to the load generator through the gearbox. The driving motor provides power for the planetary gearbox, and the speed of the motor is controlled by a speed controller, which accommodates adjustment in the range of 0-1500 r/min. The load is provided by the loading motor connected to the output shaft of the planetary gearbox, and the load torque can be adjusted by controlling the platform in the range of 0-900 N•m. The hydraulic station is responsible for providing lubricant oil pressure and shift pressure for the planetary gearbox.

B. EXPERIMENTAL SCHEME
In a planetary gearbox, sun gear teeth are the most vulnerable part to fail, since their multiplicity of meshes with the planet gear increases the possibility of damage to the sun gear. In this paper, faults of broken teeth and crack on the sun gear in the K3 row are considered because broken teeth and cracks are common fault types in gearboxes. In practice, cracks usually begin at the maximum stress point of the gear teeth and then develops along the normal line of the root curve. Thus, the cracks are processed along the normal line of the tooth's root curve by machine tools in our experiments. In addition, the surface wear fault is simulated by cutting a part of the tooth along the axis direction of the tooth. Pictures of the damaged sun gears are shown in Fig. 4. It is well known that vibration-based analysis is an effective method for diagnosing mechanical system faults, the vibra-  Table 2.

C. EXPERIMENTAL RESULTS AND DISCUSSION
In this section, we apply the proposed nonconvex WATV denoising method to diagnosis the faults of the planetary gearbox as introduced in the part B. To further demonstrate the effectiveness of the proposed method, the analysis results are compared with those of the soft thresholding denoising method and the BPD method. Compared with the vibration signals of the vertical direction, the signals of the horizontal direction are more sensitive to the damage. Because the vertical vibrations are constrained by the gravity of the test rig and the basement, their vibration amplitudes are not as large as the horizontal ones. Therefore, the vibration signals of the horizontal direction are considered in this paper. Moreover, the faults are located on the sun gear in the K3 row, so the vibration signal of the measuring point 3 closest to the fault source is selected for analysis.
In this work, the wavelet threshold is where N l is the length of the level l detail signal of wavelet analysis, σ l is the noise standard deviation of level l detail signal. And the regularization parameter is suggested using where j is the wavelet scale, η is the relative weight control parameter with a nominal value of η = 0.95.

1) CASE 1
The testing result of case 1 is shown in Fig. 6 the time-domain waveform of the vibration signal with a duration of 2 s is illustrated in Fig. 6 (a). The fault signature cannot be identified because the useful periodic pulses are buried in the strong background noise. The frequency spectrum and the magnified spectrum of the signal are shown in Fig. 6 (b) and (d), respectively. The rotating frequency of the motor f 0 and the meshing frequency of the fixed-shaft gear f m , as well as their multiples are predominant. Moreover, there are many sidebands equally spaced around the meshing frequency f m and its multiples, and the sideband spacing is the rotating frequency f 0 . The main frequency components in the Hilbert envelope spectrum (shown in Fig. 6 (c)) are the rotating frequency f 0 and its double frequency. Meanwhile, the meshing frequency f m and its sidebands are also prominent in Fig. 6 (c). However, the fault features of the planetary gearbox cannot be identified from Fig. 6, because the useful information is buried in the strong background noise. The meshing frequency of the main pump f m1 and the return pump f m2 cannot be identified in the spectrum. Because the measuring point 3 is the farthest point from the hydraulic station, so the vibration signals cannot be effectively collected by the vibration sensor of measuring point 3. Even the vibration signal can be collected, considering the complexity of the transmission path, those vibration characteristics are relatively weak in the collected signal, which cannot be identified. The meshing frequency of the K3 row f mK3 is also not identified because of the weak vibration characteristics. To extract the useful fault information, the proposed method is utilized to process the vibration signals. The results are shown in Fig. 7. From the time domain waveform of Fig. 7 (a), many irrelevant interference components and noises have been removed and strong periodic impulses are clearly revealed. From the frequency and the magnified spectrum of Fig. 7 (b) and (d), respectively, it is obvious that except for the meshing frequency f m and its multiples, as well as the sidebands with an equal interval of the rotating frequency f 0 , the gear fault characteristic frequency f s3 and its double frequency are successfully extracted in the spectrums. The fault characteristic frequency f s3 is also clearly revealed in the Hilbert envelope spectrum, as shown in Fig. 7 (c). Thus, the proposed nonconvex WATV method successfully detects the crack fault of the planetary gearbox. More specifically, the fault features of the planetary gearbox are successfully extracted utilizing the proposed method.
For comparison, we also use the soft threshold denoising method to process the same vibration signal, and the results are shown in Fig. 8. Although periodic pulses exist in the time-domain waveform as shown in Fig. 8 (a), they are not as obvious as those in the time-domain waveform as shown in Fig. 7 (a), which indicates that many noise and irrelevant interference components still exist in the analysis results. In both frequency spectrum and the Hilbert envelope spectrum, only the meshing frequency f m , rotating frequency f 0 and their multiples can be obtained; the fault characteristic frequency used to monitor the health status of the gearbox cannot be observed from Fig. 8 (b), (c), or (d).  For further comparison, the BPD denoising method is sequentially applied to the same vibration signals, with the same parameters as those used in the proposed method. The results are shown in Fig. 9. From the time domain waveform of Fig. 9 (a), the estimated periodic transient components are not as accurate as those estimated by the proposed method, because the amplitudes of the transients are underestimated. Moreover, the analysis results are not sparse enough and there are still many interference components, which effect the fault diagnosis. Thus, no useful fault characteristic frequency can be extracted from the spectrums, as shown in Fig. 9 (b), (c), and (d).
According to the mathematical principle, the analysis results of the soft threshold and the BPD method should be the same. However, compared Fig. 8 and Fig. 9, we find that the results are not the same. This is mainly because the soft threshold method is employed to deal with the decomposition coefficient after wavelet decomposition, and then reconstruct the signal using the estimated wavelet coefficients, so as to realize the noise reduction of the measured signal. The BPD method is directly to analyze the original measured signal. In addition, the parameter settings of the two methods are different. In soft threshold denoising, the threshold is set by equation (31), while in BPD algorithm, the regularization parameter is determined by equation (32). Therefore, the results of the two methods are different.

2) CASE 2
The waveform of the vibration signal for case 2 and its corresponding spectrums are shown in Fig. 10. As seen from  Fig. 10 (b), (c), and (d), the main frequency components in the frequency spectrum and the Hilbert envelope spectrum are the meshing frequency and the rotating frequency, as well as their multiples. Neither the fault characteristic frequency nor its multiples can be identified in Fig. 10. Then, the proposed nonconvex WATV denoising method is used to process the vibration signal in Fig. 10 (a), and the results are shown in Fig. 11. In the time-domain waveform in Fig. 11 (a), the transients in the analysis results are apparent, and the period of the transients is obvious. From the spectrums of Fig. 11 (b), (c), and (d), it is noted that the dominant frequencies are not only the meshing frequency f m and the rotating frequency f 0 , but also their multiples and sidebands; the fault characteristic frequency f s3 and its multiple also clearly appear in the Hilbert envelope spectrum (shown in Fig. 11 (c)) and the magnified spectrum (shown in Fig. 11 (d)). Based on the above-mentioned analysis, we can conclude that there is a localized fault on the sun gear.
For further comparison, both the soft thresholding denoising method and the BPD method are employed to analyze the vibration signals, and the results are shown in Fig. 12 and Fig. 13, respectively. Obviously, there are still many noise and interference components in the timedomain waveform of the results of the two methods, which results in the periodic characteristics of the estimated transient components in Fig. 12 (a) and Fig. 13 (a) being less prominent than those in Fig. 11 (a). In the spectrums shown in Fig. 12, in addition to the main frequency components  as shown in Fig. 10, there is also a harmonic component of the fault characteristic frequency 2f s3 in the Hilbert envelope spectrum shown in Fig. 12 (c). Despite the existence of the harmonic frequency, other evidence is necessary for the definite diagnosis of a sun gear fault. However, no other useful features related to the sun gear can be extracted from the analysis results. Similarly, in the results of the BPD method, only the second-harmonic component of the fault characteristic frequency 2f s3 exists in the Hilbert envelope spectrum shown in Fig. 12 (c), and again, no more useful features related to the fault can be extracted. Compared with the results of the proposed method shown in Fig. 11, neither the soft thresholding denoising method nor the BPD method provide satisfactory performance in the fault diagnosis of the planetary gearbox.

V. CONCLUSION
This paper proposes a nonconvex WATV denoising method for the purpose of detecting faults in a planetary gearbox. The propose method formulates wavelet thresholding and TV as a unified problem. We use a TV minimization algorithm to reconstruct the retained wavelet coefficients, so the pseudo-Gibbs oscillations phenomena caused by pure wavelet thresholding are removed. In order to strongly induce the wavelet sparsity, we consider a modification of the TV optimization problem where the L1 norm regularizer is replaced by a nonconvex penalty function. However, the use of a nonconvex regularizer converts the denoising problem into a nonconvex optimization problem. Thus, the conditions on the nonconvex degree parameter of the nonconvex penalties are addressed to ensure that the total objective function is strictly convex. Then, the convex optimization solver SALSA is introduced to find the global optimal solution.
A planetary gearbox test rig is established and two modes of faults are simulated. The vibration signals obtained under two different motor speeds are used to verify the effectiveness of the proposed method in the diagnosis of localized faults in a sun gear. Comparisons to the soft thresholding denoising and BPD methods demonstrate that the proposed method can better preserve feature components of interest and can significantly improve the estimation accuracy, thus providing improved fault detection accuracy. The proposed nonconvex WATV denoising method is concluded to be an effective tool for fault diagnosis of a planetary gearbox. FUZHOU FENG received the B.Eng. degree and the M.Eng. degree in vehicle application engineering from The Academy of Armored Forces Engineering, Beijing, China, and the Ph.D. degree in mechanical engineering from Tsinghua University, Beijing. Since 2005, he has been a Professor with the Department of Vehicle Engineering, Army Academy of Armored Forces. His current research interests include reliability assessment, fault prognostics, and health management.