A Study on Defect Identification of Planetary Gearbox under Large Speed Oscillation

Rotational speed of a reference shaft is the key information for planetary gearbox condition monitoring under nonstationary conditions. As the time-variant speed and load of planetary gearboxes result in time-variant characteristic frequencies as well as vibration magnitudes, the conventional methods tracking time-frequency ridge perform a poor robustness, especially for large speed variations. In this paper, two schemes, time-frequency ridge fusion and logarithm transformation, are proposed to track the targeted ridge curve reliably. Meanwhile, the identified ridge curve by logarithm scheme can be further refined by the timefrequency ridge fusion scheme. Hence, a procedure involving the proposed ridge estimation methods is presented to diagnose the planetary gearbox defects. Two simulation signals and a vibration signal collected from a planetary gearbox in practical engineering (provided by the conference on condition monitoring of machinery in nonstationary operations (CMMNO)) are used to verify the proposed methods. It is validated that the proposed methods can well-track the targeted ridge curve compared with two conventional methods. As a result, the characteristic frequency of each component in the planetary gearbox is clearly demonstrated and the inner race defect of one of the planet bearings is successfully discovered in the order spectrum depending on the derived expression of planet bearing fault frequency.


Introduction
Planetary gearboxes are widely applied in the drivetrains of automobiles, helicopters, wind turbines, and so forth, for the benefits of strong load-bearing capacity, compact, lightweight, and large transmission ratio [1][2][3].Due to the harsh working conditions, for example, dust, corrosion, and heavy yet unpredictable load and speed, planetary gearboxes are particularly prone to damage.Such damage can lead to a catastrophic failure of the entire mechanical system, and consequently heavy investment and productivity losses [4][5][6].Therefore, the fault diagnosis of planetary gearbox is very important for ensuring a high performance transmission.
Vibration analysis is known as an effective tool for the condition monitoring and fault diagnosis of rotary machinery [7][8][9][10].The health condition monitoring of planetary gearbox using vibration signals requires the knowledge of input or the output rotational speed, the number of stages as well as their arrangement, and the number of gear teeth.
Usually, the information about geometry parameters of gearbox is available to vibration analysis by the manufacturer's details of gearbox [11].In addition, a tachometer/encoder is installed at a reference shaft to enable the estimation of shaft speed.However, the installation of speed sensor for each machine component is not always technically feasible due to the limited space and accessibility, and not always economically viable because of the costs incurred in investment, operation, and maintenance [12].Moreover, the speed sensor measurement assumes the speed is constant between each two pulses of the key-phase signal.Therefore, it is inaccurate when the rotational speed has a large variation [13].As a conclusion in [13], several tasks of the signal processing for nonstationary signal of rotary machines are listed as follows: estimation of instantaneous rotatory speed, analysis of defect signature, and identification of rotary system parameters.It could be found that the estimation of rotational speed is the first and key task for health monitoring of gearboxes running on a nonstationary condition.

Mathematical Problems in Engineering
A fundamental method for the rotational speed estimation is to define a suitable band around the corresponding harmonics of the reference shaft and then demodulate this band [14].However, this method can estimate the instantaneous speed only with limited speed variation [15].It is mainly restricted by the results of frequency bands overlapping when the rotational speed is varied significantly.Hence, two-step method is popularly used to estimate instantaneous speed of the reference shaft [15,16].First, the one-dimension signal is mapped into two-dimension representation by the time-frequency analysis (TFA) method.Then, tack the targeted ridge curve in the time-frequency plane.Currently, a variety of TFA methods which give insights into the complex structure of multicomponent and time-variant frequency signal have been developed to deal with nonstationary signal, including short time Fourier transform (STFT) [7], Wigner-Ville distribution (WVD) [7,10], wavelet transform (WT) [8,10], synchrosqueezing transform (SST) [9], and parameterized TFA [13,17].Considerable researches have indicated that the TFA methods could be a perfect solution for processing nonstationary signal [18][19][20][21][22][23].However all TFA methods have their own merits and deficiencies.The individual TFA method may be suitable for analyzing a specific class of signals.For WVD, it can create a timefrequency representation with high time-frequency resolution for monocomponent signal, but the cross-terms are introduced inevitably for analyzing multicomponent signal [7].WT is powerless for obtaining a precise resolution when a cluster of high frequency components existed in signal [8].The SST technology is highly dependent on the original TFA method [9].The parameterized TFA methods [13,17] show a powerful ability to achieve the high time-frequency resolution for analyzing the signal with strongly nonlinear instantaneous frequency.However, its several limitations should not be neglected, such as computational complexity and suitable matching mathematical models.Considering that the analyzed vibration data are numerous and the complicated multicomponents are buried in the signal, we utilize the simplest STFT as the TFA technique in this paper.Although the STFT could not achieve an arbitrarily high time-frequency resolution at the same time, [24,25] indicate that STFT can distinguish the adjacent frequencies from each other and is suitable for estimating the smooth ridge curve when an enough window width is given.
Subsequently, it is to identify the targeted time-frequency ridge from the time-frequency plane.As the simplest technique, the direct ridge detection algorithm, called modulus maximum method [26,27], is usually used in signal decomposition and ridge detection.The direct ridge detection algorithm is to detect the maximum magnitude point in the frequency direction for every time instant.It only considers the magnitude of each point in time-frequency plane but ignores the smoothness of the time-frequency ridge.As many peaks exist in amplitude of the time-frequency representation at each time, and their number often varies in real cases, the location of maximum amplitude may be not the actual ridge point.In these circumstances, it can be unclear which peak corresponds to which component, and which are the artifacts just induced by noise; thus it is sensitive to noise and it is not suitable to estimate instantaneous frequency of the low signal-to-noise (SNR) signal.In addition, a ridge detection algorithm for continuous wavelet transform has been proposed based on a cost function [27], which was nonsensitive to noise compared with the direct maximum method.Later, another ridge detection algorithm for STFT was proposed to detect ridge curve [24] based on the same cost function, which enhanced further the computational cost by shrinking the search region.However, its searching direction may restrict its adaptation in some applications.An improved ridge detection algorithm was proposed by minimizing the local cost functions dynamically, in which both the large amplitude in the local area and the smoothness of ridge curve were considered [25,28].In this paper, it will be shown that the improved ridge detection algorithm in [25,28] still suffers from the interferences of local bright noise and adjacent strong ridge curve when a weak ridge is extracted.Therefore, this study is mainly concentrated on the problem of robust ridge curve estimation.And then, on this basis, the order analysis is utilized to transform the nonstationary vibration signal in time-domain into stationary signal in angular domain.In this way, the smearing problem caused by the speed variation can be solved effectively and the detection of the characteristic frequencies of planetary gearbox is performed to identify the defective component in planetary gearbox.
The paper is organized as follows.The experimental setup of planetary gearbox and the motivations of this study are described in Section 2, respectively.In Section 3, two ridge extraction schemes are proposed, and the simulated signal and experimental signal are used to verify the proposed methods.The defect identification of planetary gearbox using the estimated rotational speed is presented in Section 4. A detailed discussion is given in Section 5. Finally, a conclusion is drawn in Section 6.

A Planetary Gearbox in Practical Engineering
2.1.Experimental Set-Up.A schematic of the planetary gearbox is shown in Figure 1.The planetary gearbox is employed in an actual wind turbine.The whole gearbox has three stages (one planetary and two helical parallel stages, neglecting the gear pair 8-9).Table 1 lists the gear parameters of the planetary gearbox.As described in Figure 1, the red shaft represents the input shaft and the yellow shaft stands for the high speed shaft.When the input speed of planet carrier is given at any time, the characteristic frequency of each component in the planetary gearbox can be calculated by using the expressions of planetary gearbox characteristic frequencies listed in Table 2, which are deduced from the configuration of planetary gearbox.A vibration signal collected from this planetary gearbox is provided by the conference on condition monitoring of machinery in nonstationary operations (CMMNO).The collected vibration data are sampled at 5 kHz during speed variation and the duration time is about 547.4 seconds.The vibration data are collected from an accelerometer located on the rotor side of the gearbox casing in the radial direction.It is only known that the speed of input  shaft is between 13 r/min and 15 r/min during the recording duration.Therefore, it can be deduced that the rotational speed of high speed shaft is between 1550 and 1800 r/min according to the information provided in Tables 1 and 2. The other available information about this gearbox is limited.

Motivations of This Study.
The time-domain waveform and the frequency spectrum of the raw vibration signal collected from planetary gearbox are shown in Figure 2, respectively.It can be seen from the time-domain waveform that its amplitude performs the intense fluctuation.It is indicated that the energy of the component buried in the raw signal exists in the remarkable change with increasing of time.As shown in Figure 2(b), the regular spectrum lines cannot be clearly detected in the Fourier spectrum.Thus, it is powerless for condition monitoring of planetary gearbox using the conventional technique.Due to the deficiency of the Fourier transform, it is necessary to find a more effective method to analyze nonstationary signals.As described in Figure 3, the procedure is usually employed for condition monitoring of planetary gearbox under nonstationary condition.For the second step, the STFT is utilized to analyze the raw vibration signal which has been discussed in the introduction section, where a hamming window of 2 15 samples is used to obtain a high frequency resolution, considering the input rotational speed between 13 r/min and 15 r/min.As presented in Figure 4, the time-frequency representation obtained by STFT can reveal more information than the Fourier transform.A lot of time-frequency ridge curves occur in the time-frequency plane and oscillate with the change of time.Several ridge curves containing specific properties are pointed out in the time-frequency plane.The ridge curve indicated by line 1 is attributed to the high energy time-frequency ridge curve.The ridge curve indicated by line 3 and line 4 belonging to the inherent frequency components could be caused by the structure resonance and interference of the electrical noise.According to the information of the high speed shaft, the meshing frequency of high speed shaft locates roughly between 750 Hz and 870 Hz due to the 29 teeth of gear  7 .Figure 4(b) is a close-up view of Figure 4(a).It can be seen from Figure 4(b) that the ridge curve of line 2 is between 750 Hz and 870 Hz roughly.Therefore, line 2 is likely the meshing frequency curve of high speed shaft.For the third step, it is firstly to determine the targeted ridge curve and then extract this ridge curve.Generally, the characteristic frequencies of the components in planetary gearbox modulated by the strong energy sidebands [2] are unsuitable as the targeted ridge curve.In addition, the high energy ridge curve out of sync with the rotational speed of the reference shaft is also not suitable as the targeted ridge curve, for example, line 1, which will be proved in Section 3. As shown in Figure 4(b), the modulation sidebands around the meshing frequency of high speed shaft are almost invisible because the gear on the high speed shaft is far away from the center of planetary gear train.Hence, the ridge curve indicated by line 2 is selected as the targeted ridge curve.Next, as discussed in Introduction, the method proposed in [25,28] is firstly attempted to identify these ridge curves from the time-frequency plane.
The method proposed in [25,28] is simply described as follows.
(1) Determine the initial searching point (  , f  ) in the time-frequency plane (TF(  ,   ),  = 1, 2, . . ., ,  = 1, 2, . . ., ), where  and  correspond to the number of discrete time and discrete frequency, respectively.In this case, it is defined by the local maximum in the time-frequency plane.As shown in Figure 5(a), the white rectangles near the corresponding targeted ridge curve represent the local regions where the point with the maximum is selected as initial searching point.
(2) Search the targeted ridge curve based on the cost function  , for the other  − 1 time instants defined as Table 2: Characteristic frequencies of the components in planetary gearbox.

Characteristics frequency Expression
Meshing frequency of gear pairs 1-2 and 2-3 Meshing frequency of gear pair 4-5 Meshing frequency of gear pair 6-7 is the rotational frequency of the red input shaft, that is, the rotational frequency of planet carrier.Step 2 Step 4 Step 3 Step 5 Step 1 where is the weight factor to adjust the proportion between the amplitude and smoothness of the ridge curve and f  is the frequency minimizing the cost function  , : Generally, all frequencies at   are taken as the candidates, but the computational load is heavy and the large amplitude of adjacent ridge curve could cause the failure of the ridge curve extraction.Thus, the searching range is dynamically restricted to a small frequency band FB  : where   is the half bandwidth.In this case,   is set as 20 Hz which is a tradeoff between the computational time and stability because the space between line 1 and line 2 approximates 100 Hz.Moreover, the frequency band FB  is dynamically modified with the change of f  due to the essential property of the ridge extraction method in [25,28].
As shown in Figure 5(a), the high energy ridge curve indicated by line 1 could be well-extracted by the ridge extraction algorithm in [25,28].For the targeted ridge curve indicated by line 2, we make an attempt to use two different initial searching points to track it.However, some remarkable deviations are produced when the initial points are selected from the local regions in left and right hands of the time-frequency plane.In addition, some close-up views of Figure 5(a) are shown in Figure 5(b).It can be seen that some small deviations exist between the two extracted ridge curves using the different initial searching points.
As discussed above, it is almost impossible to achieve a satisfactory result for the targeted ridge curve based on the current time-frequency ridge detection scheme.Consequently, it is necessary to develop alternative methods to address this problem.Intuitively, it can be concluded that two solutions could be utilized to settle the problem from the above analysis.One way to obtain a closer result of the targeted ridge curve is to combine and fuse the useful information buried in the ridge curves extracted by using different initial points.The data fusion technique [29] is the promising tool to accomplish this objective.Another way is to design a more rational cost function for the current time-frequency ridge detection scheme.As displayed in Figures 6(a

Algorithm of Time-Frequency Ridge Fusion.
As aforementioned, the ridge curves obtained by the method in [25,28] using different initial searching points perform their own deficiencies and merits, and the fusion technique can be used to identify the useful information in the candidate ridge curves.Hence, we propose a novel scheme named as the time-frequency ridge fusion to achieve a more accurate ridge curve than that provided by a sole ridge curve.The basic idea of this scheme is to retain the useful features in the extracted ridge curves based on the smoothness property of real ridge curve and singularity of illusive ridge curve.These messages can be uncovered by the differential result of the multiple ridge curves corresponding to different initial searching point.It should be noted that the proposed scheme is different from the time-frequency data fusion technique in [21,30].The algorithm in [21,30] belongs to Step 2 in Figure 3. Figure 7 gives an illustration of the time-frequency ridge fusion scheme.The detailed procedure of the proposed fusion scheme can be described as follows.
First.Give several initial searching points for a targeted ridge curve.These initial searching points are defined by the maximum in the local regions which are usually located in both ends of the analyzed signal.It is worth noting that these regions should be close enough to the targeted ridge curve.In addition, the multiple local regions can be distributed roughly evenly in the time-frequency plane when more messages are needed to represent a whole targeted ridge curve under the worse measurement condition.
Second.Search the time-frequency ridge curves f1 , f2 , . . ., f ( represents the number of initial searching points) based on the cost function scheme in [25,28] using the initial searching points given in the above step.Usually, these extracted ridge curves are not identical:  Then, it is known that all of these time-frequency ridge curves can be regarded as a complete set : And the real targeted time-frequency ridge curve f is identified as set : Third.Design a fusion engine to define set  based on the smoothness property of real ridge curve and singularity of illusive ridge curve.In particular, data fusion engine is designed as follows.
(1) Calculate the differential result between the ridge curves f and f(+1) obtained by using the adjacent initial searching points (  ,   ) and ( +1 ,  +1 ): (2) Define the fusion region   according to the differential curve Δ f as where    and    represent the left and right endpoints of fusion region   defined as where  ∈ [   ,   +1 ] and    and   +1 stand for the serial number of   and  +1 , respectively. is a threshold (1∼2 times of frequency resolution).There are  − 1 fusion regions in the whole time-domain.
(3) Search the specific fusion segments in the fusion region   according to the differential curve Δ f .Generally, the segments are considered as the demanded fusion locations where the ridge curves f and f(+1) are not in good agreement with each other.Therefore, the range between a pair of rising and falling edges in the differential curve Δ f is a small fusion segment defined as Δ f ( : ). and  are calculated as where Δ represents the difference between the length of ridge curves f ( : ) and f(+1) ( : ) in a small fusion segment as follows: where Finally.Output the fusion ridge curve f as the approximately targeted ridge curve.
Remark.The regions [1,    ] and [  +1 , ] at both ends of the analyzed signal cannot be processed in this procedure.But they are usually some very small pieces and close to the initial searching points.Therefore, f1 (1 :    ) and f (  +1 : ) could be directly regarded as f .

Simulation Evaluation.
In this subsection, a simulated signal [19] is employed to examine the performance of the proposed fusion scheme as follows: where   () represents the rotational component to simulate the speed variation given as As shown in Figure 8(a), the time-frequency representation of the simulated signal at a SNR −8 dB is obtained by STFT where a hamming window of 5000 samples is employed to acquire an enough frequency resolution.The time-frequency representation uncovers the three components buried in the analyzed signal.Then, the targeted ridge curve is tracked by the cost function scheme in [25,28] using two initial searching points.The parameters used in the cost function are set the same as the contents in Section 2.2.However, the targeted ridge curve cannot be represented by both of the extracted ridge curves.Thus, the fusion scheme is employed to obtain the targeted ridge curve.As described in Figure 8(b), the useful information in the fusion process is represented in the differential result of the two ridge curves.Moreover, some small segments which should be further improved are clearly shown in the differential curve.As presented in Figure 8(c), the fusion result is well consistent with the theoretical ridge curve.To compare with the proposed method, the results obtained by the cost function scheme in [25,28] and modulus maximum method in [26,27] are shown in Figure 8(d), respectively.It can be seen that the ridge curve obtained by modulus maximum method exists in the large deviation in the early stage.The case of divergency also occurs in the ridge curve extracted by the cost function scheme at a later stage.And some small deviations exist in the middle place of the ridge curve extracted by the cost function scheme.All of these results verify that the proposed method can obtain an approving result approaching to the theoretical ridge curve.
To simulate the worse measurement condition, the simulated signal at SNR −10 dB is utilized to verify the proposed method.Figure 9(a) demonstrates the time-frequency representation of the simulated signal where the parameters are set as the above analysis.Then, two initial searching points are defined by the maximum in the local regions.However, the  ridge curves obtained by the cost function scheme in [25,28] cannot represent the whole targeted ridge curve shown in Figure 9(b).And then, an additional initial searching point is supplemented in the middle place of the analyzed signal.As presented in Figure 9(b), the three ridge curves using different initial searching points represent seemingly the targeted ridge curve.Figure 9(c) gives the differential curves of them.The fusion result is provided in Figure 9(d).It can be seen that the targeted ridge curve can be well-performed by the fusion scheme at an extremely low SNR.

Experimental Evaluation.
In this subsection, we validate the proposed fusion scheme using the planetary gearbox vibration data given in Section 2.1.As introduced in Section 2.2, two ridge curves obtained by the cost function scheme using the initial searching points at the left and right ends of the signal cannot represent the whole targeted ridge curve.It can be seen from Figure 10 that the two fail regions are clearly described in the differential curve and the enlarged view of its medium region displays a lot of fusion segments.Figure 11 shows the result obtained by the proposed fusion scheme.The fused ridge curve can track the variation of the targeted ridge curve by visual observation.According to the obtained ridge curve, the rotational speed of the high speed shaft is calculated due to the 29 teeth of gear  7 .As shown in Figure 12, the rotational speed of the high speed shaft obtained by the fusion scheme is between 1550 and 1800 r/min at most operating time which is in good agreement with the discussion in Section 2.1.
Furthermore, the results of the rotational speed of high speed shaft obtained by modulus maximum method [26,27] and cost function scheme in [25,28] are employed to compare with the proposed method.As demonstrated in Figure 12, the rotational speed extracted by the modulus maximum method has a little dance in the early stage and produces a big break in the place of near 100 seconds.The rotational speed extracted by the cost function scheme in [25,28] presents the large jumps at both ends of the signal and the small fluctuations arise in some places.On the contrary, the fusion scheme gives a smooth rotational speed curve in the whole time-domain.The above analysis results verify that our proposed method delivers a good robustness for identifying the weak and intense fluctuant ridge curve.The accuracy of the extracted time-frequency ridge will be further demonstrated later.

Algorithm of the Logarithm Scheme.
We have summarized that the logarithm transformation of time-frequency representation could produce a better demonstration than the original one in Section 2.2.Hence, the cost function  , in [25,28] is redesigned as All other parameters used in the logarithm scheme are the same as the contents of cost function scheme introduced in Section 2.2.

Simulation Evaluation.
The simulated signals presented in Section 3.1.2are employed to verify the logarithm scheme in this subsection.As shown in Figure 13(a), the ridge curve for the simulated signal at SNR −8 dB extracted by the logarithm scheme can track the targeted ridge curve roughly and perform a more robust result than the cost function scheme in [25,28] (result of cost function scheme was given in Section 3.1.2).However, some small deviations still exist in the extracted ridge curves.Then, the result can be further refined by the fusion scheme shown in Figure 13(b).
Similarly, the simulated signal at SNR −10 dB is processed by the logarithm scheme.As described in Figure 14(a), the ridge curves extracted by the logarithm scheme only using two initial searching points can represent the complete targeted ridge curve, but three initial searching points for the cost function scheme in [25,28] are needed to obtain a whole targeted ridge curve as given in Section 3.1.2.Although the extracted ridge curves also produce some deviations, the fusion scheme can be used to further improve the result as shown in Figure 14(b).

Experimental Evaluation.
In this subsection, the logarithm scheme is applied to analyze the experimental signal given in Section 2.1.As demonstrated in Figure 15(a), the targeted ridge curve is extracted by the logarithm scheme and the evident failure region does not arise in the extracted ridge curve.But they still exist in some local deviations as described in Figure 15(b).Then, the fusion scheme proposed in Section 3.1 is used to further process the extracted ridge curves in Figure 15(a).As shown in Figure 16, the rotational speed of the high speed shaft obtained by the logarithm scheme (containing the fusion scheme) is also between 1550 and 1800 r/min at most of the time.In addition, the high energy curve indicated by line 1 can be also well-tracked.

Further Validation of the Proposed Methods.
To remove the effect of speed fluctuations, order analysis is usually performed to convert the equal time sampling to equal circumferential angle sampling [12,16].It could be expected that the spectrum lines in order spectrum corresponding to the meshing frequency of each component in the planetary gearbox will be evidently observed when the targeted ridge curve is accurately identified by the proposed methods.For this reason, the order analysis is employed to validate the accuracy of the extracted rotational speed for the experimental signal in Sections 3.1 and 3.2.
As shown in Figure 16, the order maps are achieved by order analysis using the rotational speed of high speed shaft.The results obtained by the fusion scheme and logarithm scheme are consistent with each other.It can be seen that a large amount of straight lines are synchronous to the targeted ridge curve in the low order region.In addition, the ratio between the ridge curves indicated by line 1 and line 2 is given in Figure 18.The high energy ridge curve indicated by line 1 is not synchronous to line 2.This phenomenon is also clearly shown in Figure 17.As the time-frequency representations shown in Figure 17 are similar, we only show one type of the obtained order spectra to save space in the following study.Figure 19 shows the order spectrum of the planetary gearbox.Comparing with the original Fourier spectrum in Figure 2(b), we can detect a lot of peaks in Figure 19.According to the information in Table 2, the relationship between the rotational frequency of high speed shaft and the characteristic frequencies of all components in the planetary gearbox is given in Table 4. Then the order corresponding to each component can be calculated using the gear parameters in Table 1.Next, we will give the physical meaning of these peaks in the order spectrum.
As presented in Figure 19, the meshing frequency  45 of gear pair 4-5 and its harmonics are clearly shown as discrete lines.In addition, the meshing frequency  67 of gear pair 6-7 is also clearly presented in Figure 19.A zoom-in in the region between order 0 and order 6.5, as demonstrated in Figure 20, reveals the meshing frequency  12 of planet gear and sun gear, equally for planet gear and gear ring, and its harmonics markedly.Moreover, the modulated sidebands (planet carrier rotational frequency   ) around the meshing frequency  12 of planet gear and sun gear are also clearly described in Figure 21.As the above discussion, the one-toone correspondence relationship between all these evident discrete order lines and characteristic frequencies of the components in planetary gearbox is uncovered.It is verified that the reference speed signal extracted by our proposed methods can well remove speed fluctuations from the raw vibration signal.And its accuracy is good enough to analyze to the vibration signal of planetary gearbox in practical engineering.

Defective Identification of Planetary Gearbox
Most failures of gearbox could initiate in bearings and later advance to gears and other components [31].Additionally, the faulty rolling bearing existing in the planetary gearbox has been informed by CMMNO.As for the defective bearing in planetary gearbox, it is either located on the fixed-axis shaft or planet carrier shaft.If the faulty bearing is located on the fixed-axis shaft, it is easy to calculate the defective frequency as [32].However, the bearings supporting the planet gears, that is, planet bearings, exhibit a high failure rate and most gearbox failures initiate in them [33].Moreover, the spectral structure of defective planet bearing vibration signal is very complex due to periodically time-variant working condition, as well as the effect of vibration transfer path in planetary gearbox.Thus, it is necessary to establish the vibration signal model of the defective planet bearing according to the above discussion.To our best knowledge, the vibration response of a planet bearing with localized defects (spalls or pits) has been studied in only a very few papers [31,33,34].A lumped parameter model of a planetary drivetrain containing a planet bearing with localized defects was built in [33], which did not consider the dominant meshing effect of gear pairs.In addition, a fault diagnosis method for planet bearing using an internal vibration sensor was proposed in [34] where an accelerometer was mounted internally on the planet carrier to address the issues of variable transmission path.In fact, the faulty planet bearing is regarded as the fixed-axis bearing in [34].In this section, the explicit equations for calculating the characteristic frequencies of defective planet bearing are derived and the vibration signature of planet bearing presented in Fourier spectrum is also discussed.Later, the defective component in planetary gearbox is detected in the order spectrum.

Vibration Signal Model of Defective Planet
Bearing.In the case of the defective planet bearing, a variety of sidebands which are caused by the amplitude and frequency modulations appear in the frequency spectrum.In the subsection, how to calculate characteristic frequency of planet bearing and the sideband behavior in the frequency spectrum will be presented and discussed.A clear and explicit expression of the movement relationship between the planet bearing and planet carrier shaft should be firstly given.As the outer race of planet bearing is usually fixed to the planet gear, the movement direction between the planet gear and planet carrier is an inverse relationship; thus the relative rotational frequency of a planet gear to the planet carrier is given as where  2 is the absolute rotational frequency of planet gear and can be obtained by the fixed planet carrier method as Then, the equations of characteristic frequencies of faulty planet bearing could be derived from the movement relationship in (16).As given in (18), the fault characteristic  frequencies of planet bearing are different from the fixed-axis shaft bearing: where  represents the coefficient of the local defective frequency listed in Table 3;   ,   , and   are ball pass frequency, outer race (BPFO); ball pass frequency, inner race (BPFI); and ball spin frequency (BSF), respectively.In addition, the vibration signature of a defective planet bearing is different from a fixed-axis bearing because of the complex and time-varying vibration transmission path.The reasons include the following: (1) the variation in vibration transmission path is caused by the rotation of planet carrier [35,36], (2) the variation in impulse magnitude is caused by the rotation of an outer race defect relative to the load zone, and (3) the variation between the impulse direction of an outer race defect and planet gear meshing direction is caused by the rotation of planet gear [33].Therefore, the Fourier spectrum of vibration signal for a faulty planet bearing is very complex even if the vibration data are collected from the planetary gearbox under the stationary condition.As for BPFI, since the defect does not move relative to the load zone, the spectral structure of inner race defect of planet bearing can be described by  ⋅   +  ⋅   , where  and  are integers, respectively, due to the variation in vibration transmission path caused by the planet carrier rotation.As for BPFO, the variations in vibration transmission path, impulse magnitude, and force direction between impulse and planet gear meshing exist.Therefore, the spectral structure of outer race defect of planet bearing can be represented by  ⋅   +  ⋅   +  ⋅    2 , where , , and  are integers, respectively.As for the BSF, they also exist in the variation in vibration transmission path caused by the rotation of planet carrier.In addition, the faulty characteristic frequency of element is modulated by the self-rotating frequency of bearing cage   which can be calculated the same way as (18).Thus, the spectral structure of ball defect of planet bearing can be shown by  ⋅   +  ⋅   +  ⋅   , where , , and  are integers, respectively.

Result of Defective Identification.
As described in [33], most failures of gearbox initiate in bearings and the planet bearings are considered as one of the most critical components with very high failure rate.Meanwhile, the meshing frequency of planetary gear and sun gear and its harmonics are dominantly modulated by the planet carrier passing frequency 3  (3 planet gears) shown in Figure 22.And a close-up view in the region between order 0 and order 0.8, as demonstrated in Figure 23, also displays the planet carrier passing frequency and its harmonics.These phenomena in the order spectrum indicate that the three planet gears are serious asymmetry.Therefore, the most suspected defective bearing is one of the three planet bearings according to the information shown in Figures 22 and 23.Furthermore, as described in Figure 24, two groups of the characteristic frequencies 3  ±   and 4  ±   are clearly displayed in the order spectrum between order 0.5 and order 0.8.These characteristic frequencies are in good agreement with the spectral structure of inner defect of planet bearing discussed in Section 4.1.According to above analysis, it can be judged that a local defect exists in the inner race of one of the planet bearings.

Discussions
The reference rotational speed estimation is the first step and also a vital procedure for realizing fault diagnosis of planetary gearbox under the varying speed condition.However, the installation of speed sensor for each and every machine component is not always technically feasible due to the different limited factors.Therefore, the time-frequency ridge fusion scheme and logarithm scheme were proposed to extract the rotational speed of reference shaft, respectively.And the result obtained by logarithm scheme needs to be further enhanced by the fusion scheme.The analysis results verified that the proposed methods possessed the capacity to overcome the drawback of conventional ridge detection algorithms.However, it is worth noting that the fusion scheme has still much space to be improved in the future.Currently, considering the smoothness of the targeted ridge curve, we used the shortest route principle in the fusion scheme.Some other criterions, such as averaging the candidate ridge curves Mathematical Problems in Engineering in the fusion region and using the gradient information of the extracted ridge curve, are worth investigation in the future.In addition, the number of initial searching points can be determined by specific requirement of engineering practices.Especially for the worse measurement condition, more initial searching points need to be used.Moreover, the definition of the targeted ridge curve is also critical for realizing the fault diagnosis of planetary gearbox.In fact, the ridge curve with the high energy may be not synchronous to the reference shaft and some synchronous ridge curves in the low frequency region are usually modulated by the strong sidebands.Thus, the ridge curves under these situations are unsuitable as the targeted ridge curve.In this study, it is verified that the ridge curve corresponding to the gear meshing frequency of high speed shaft is suitably regarded as the targeted ridge curve.As a result, the research in this paper has the following main contributions and advantages.(1) A time-frequency ridge fusion scheme for the rotational speed estimation is proposed to remedy the deficiency of the conventional ridge detection algorithm.(2) A logarithm scheme for timefrequency ridge extraction is proposed from the view of strengthening the weak ridge and reducing the amplitude discrepancy.Meanwhile, the extracted ridge curves by the logarithm scheme can be further refined by the fusion scheme.
(3) Considering that the movement formation of planet bearing is different from that of the fixed-axis bearing, the explicit equations of planet bearing defective frequencies are deduced and the spectral structures of these defects are discussed.( 4) One implication of this study is that the analyzed vibration data are collected from the planetary gearbox in practical situation, but many researches are focused only on the case of seeded fault for planetary gearbox.

Conclusions
A study on the defect identification of planetary gearbox under the intensely nonstationary condition has been introduced in this paper.To estimate the instantaneous speed of the reference shaft reliably, two novel schemes named as the time-frequency ridge fusion scheme and logarithm scheme are proposed.Then, the simulated signals and experimental signal are employed to evaluate the proposed methods, respectively.It is validated that the proposed methods can well-identify the targeted ridge curve compared with the conventional methods.Moreover, the accuracy of the estimated rotational speed by the proposed methods is further validated by the order analysis.Afterwards, the order spectrum is utilized to detect the fault in planetary gearbox.As a consequence, the inner race defect of one of the planet bearings is confirmed by detecting the order spectrum based on the detailed derivation for the expression of planet bearing faulty frequency.For future research, it is planned for us to design the more reliable and rational fusion criterions, that is, the advanced fusion engine in the fusion scheme for extracting the ridge curve in the worse working environments.

Figure 3 :
Figure 3: Flow chart for condition monitoring of planetary gearbox under nonstationary condition.

Figure 4 :
Figure 4: (a) Time-frequency representation obtained by STFT and (b) close-up view of (a).
) and 6(b), the logarithm transformation of time-frequency representation by STFT produces seemingly a better demonstration than that in Figures4(a) and 4(b).It is because the logarithm transformation can weaken the disparity of amplitude.In the next section, we will design the more robust ridge extraction methods based on the views of fusion and logarithm transformation.

Figure 5 :
Figure 5: (a) Ridge curves extracted by the method in [25, 28] and (b) local failure regions of the targeted ridge curve.

Figure 6 :Figure 7 :
Figure 6: (a) Logarithm transformation of the time-frequency representation and (b) close-up view of (a).

)
= 10 for simulating the medium-high frequency components, without loss of generality   = 0, and  10 = 0.18,  11 = 0.1, and  12 = 0.18.The second component is regarded as the targeted ridge curve with the weak energy.A Gaussian white noise () is added to mimic the background noise and interferential component.The sampling frequency is 10000 Hz and sampling time is 10 seconds.

Figure 8 :Figure 9 :
Figure 8: Simulated signal at SNR −8 dB: (a) time-frequency representation, (b) differential curve, (c) fusion result, and (d) comparison between the proposed method and the conventional methods.

Figure 10 :
Figure 10: Differential result of the time-frequency ridge curves using two different initial searching points.

Figure 11 :Figure 12 :
Figure 11: Result obtained by the fusion scheme.

Figure 13 :Figure 14 :
Figure 13: Simulated signal at SNR −8 dB: (a) results obtained by logarithm scheme and (b) fusion result of the ridge curves in (a).

Figure 15 :Figure 16 :
Figure 15: Ridge curve extraction by logarithm scheme: (a) overall ridge curve in time-frequency plane and (b) local failure regions of the targeted ridge curve.

Figure 17 :Ratio between line 2 and line 1 Figure 18 :
Figure 17: Order map (a) obtained by fusion scheme and (b) obtained by logarithm scheme (containing the fusion process).

Figure 19 :
Figure 19: Order spectrum obtained by order analysis using the rotational speed of high speed shaft.

Figure 22 :Figure 23 :
Figure 22: Meshing frequencies of gear pair 1-2 and their demodulated sidebands by planet carrier passing frequency.

Table 3 :
Coefficient of the characteristic frequency of defective planet bearing.

Table 4 :
Relationship between the rotational frequency   of high speed shaft and the characteristic frequencies of all other components.