Iterative generalized time–frequency reassignment for planetary gearbox fault diagnosis under nonstationary conditions

Mechanical Systems and Signal Processing ∎ (∎∎∎∎) ∎∎∎–∎∎∎

Contents lists available at ScienceDirect

Mechanical Systems and Signal Processing journal homepage: www.elsevier.com/locate/ymssp

Iterative generalized time–frequency reassignment for planetary gearbox fault diagnosis under nonstationary conditions Xiaowang Chen, Zhipeng Feng n School of Mechanical Engineering, University of Science and Technology Beijing, Beijing 100083, China

a r t i c l e i n f o

abstract

Article history: Received 29 February 2016 Accepted 15 April 2016

Planetary gearboxes are widely used in many sorts of machinery, for its large transmission ratio and high load bearing capacity in a compact structure. Their fault diagnosis relies on effective identification of fault characteristic frequencies. However, in addition to the vibration complexity caused by intricate mechanical kinematics, volatile external conditions result in time-varying running speed and/or load, and therefore nonstationary vibration signals. This usually leads to time-varying complex fault characteristics, and adds difficulty to planetary gearbox fault diagnosis. Time–frequency analysis is an effective approach to extracting the frequency components and their time variation of nonstationary signals. Nevertheless, the commonly used time–frequency analysis methods suffer from poor time–frequency resolution as well as outer and inner interferences, which hinder accurate identification of time-varying fault characteristic frequencies. Although time– frequency reassignment improves the time–frequency readability, it is essentially subject to the constraints of mono-component and symmetric time–frequency distribution about true instantaneous frequency. Hence, it is still susceptible to erroneous energy reallocation or even generates pseudo interferences, particularly for multi-component signals of highly nonlinear instantaneous frequency. In this paper, to overcome the limitations of time– frequency reassignment, we propose an improvement with fine time–frequency resolution and free from interferences for highly nonstationary multi-component signals, by exploiting the merits of iterative generalized demodulation. The signal is firstly decomposed into mono-components of constant frequency by iterative generalized demodulation. Time–frequency reassignment is then applied to each generalized demodulated mono-component, obtaining a fine time–frequency distribution. Finally, the time–frequency distribution of each signal component is restored and superposed to get the time– frequency distribution of original signal. The proposed method is validated using both numerical simulated and lab experimental planetary gearbox vibration signals. The timevarying gear fault symptoms are successfully extracted, showing effectiveness of the proposed iterative generalized time–frequency reassignment method in planetary gearbox fault diagnosis under nonstationary conditions. & 2016 Elsevier Ltd. All rights reserved.

Keywords: Planetary gearbox Fault diagnosis Nonstationary Time–frequency analysis Reassignment Generalized demodulation

n

Corresponding author. E-mail address: [email protected] (Z. Feng).

http://dx.doi.org/10.1016/j.ymssp.2016.04.023 0888-3270/& 2016 Elsevier Ltd. All rights reserved.

Please cite this article as: X. Chen, Z. Feng, Iterative generalized time–frequency reassignment for planetary gearbox fault diagnosis under nonstationary conditions, Mech. Syst. Signal Process. (2016), http://dx.doi.org/10.1016/j. ymssp.2016.04.023i

X. Chen, Z. Feng / Mechanical Systems and Signal Processing ∎ (∎∎∎∎) ∎∎∎–∎∎∎

2

1. Introduction Planetary gearboxes have extensive applications in rotorcrafts, ships, wind turbines, excavators and other industrial equipment, owing to their merits of high transmission ratio and efficiency even under heavy load and low speed conditions in a compact structure. Nevertheless, tough environment conditions (wind gust, dust, erosion, etc.) leave planetary gearboxes vulnerable to gear faults. The immediate consequence varies from poor efficiency to even catastrophic property and life losses. Hence, planetary gearbox fault diagnosis is a vital topic and has drawn increasing attentions in recent decade. In order to understand the dynamic characteristics of planetary gearboxes and enable fault diagnosis, planetary gearbox vibrations and their Fourier spectra have been investigated by many researchers. Ambarisha and Parker [1] studied the nonlinear dynamic behaviors of spur planetary gears using both a lumped parameter model and a finite element model. Mark and Hines [2] predicted the spectral behavior of planetary gearbox vibration signals based on a first-principles mathematical model. Inalpolat and Kahraman [3,4] proposed dynamic models for modulation sidebands prediction considering different gearbox configuration as well as manufacturing errors, and investigated the spread of sidebands in the vicinity of gear meshing frequency and its harmonics. Feng [5] further analyzed the modulation sources, proposed a mathematical modulation vibration signal model, and derived the equations for calculating both localized and distributed gear fault characteristic frequencies. The above researches show that planetary gearbox vibration signals feature more intricate modulation and more complex spectrum than conventional fixed-shaft gearboxes. Consequently, fault diagnosis methods for fixed-shaft gearbox often become ineffective in detecting planetary gearbox fault. Researchers have proposed some methods to extract fault feature from planetary gearbox vibration signals. Samuel and Pines [6] presented a constrained adaptive wavelet transform and CAL4 metric to diagnose the localized gear tooth fault of planetary gearboxes. Feng et al. [7,8] addressed the spectrum complexity issue via demodulation analysis based on ensemble empirical mode decomposition and energy separation, as well as local mean decomposition methods. Barszcz and Randall [9] used spectral kurtosis to detect the gear tooth crack of a wind turbine planetary gearbox. Lei [10] evaluated the health condition of a multi-stage planetary gearbox using multiclass relevance vector machine. Hong and Dhupia [11] detected and located the gear tooth defect of a planetary gearbox by integrating fast dynamic time warping and correlated kurtosis technique. Sun et al. [12] customized the optimal multi-wavelets based on improved neighboring coefficients to detect incipient gear pitting fault of a planetary gearbox in a satellite communication antennas. Most existing methods are based on a stationarity assumption, i.e. constant speed and load. However, in many real applications, for example, for wind turbine or rotorcraft transmission systems, the volatile wind power or intricate flight mission always result in nonstationary conditions. Such unsteady conditions lead to nonstationary signals, and therefore time-varying fault characteristics. Time variation of fault characteristic frequencies will undoubtedly make conventional stationarity based analysis methods ineffective. To date, reported literature on the topic of planetary gearbox fault diagnosis under nonstationary conditions have been very limited. These include [13–15] by Bartelmus, Zimroz and Bartkowiak. Based on the fact that gear fault is susceptible to external load, they proposed an indicator that reflects the linear dependence between the amplitude of meshing frequency and the operating condition, and presented a procedure for load-dependent feature processing, to monitor planetary gearbox under time-varying running conditions. These studies have enriched the literature on planetary gearbox condition monitoring. Planetary gearbox fault diagnosis essentially relies on identification of fault characteristic frequencies. Therefore, extraction of time-varying fault characteristic frequencies is a helpful alternative means for planetary gearbox fault diagnosis under nonstationary conditions. Joint time–frequency domain analysis bridges the gap between the time and frequency domains, and exhibits the time-variation of each constituent frequency component, thus providing an effective approach to nonstationary signal analysis. Among joint time–frequency domain analysis, linear and bilinear time–frequency representations are widely used [16]. Linear time–frequency representations rewrite a signal as a weighted sum of a series of bases localized in joint time–frequency domain. Typical representative methods are short time Fourier transform and wavelet transform. They require an appropriate selection of basis function adapted to the major signal structure. In addition, they are subject to Heisenberg uncertainty principle, and therefore suffer from imperfect time–frequency resolution [17]. Bilinear time–frequency representations exhibit a signal energy distribution in joint time–frequency domain. Wigner–Ville distribution is the core of bilinear time–frequency representations. It has the highest time–frequency resolution, but suffers from severe outer (cross term) interferences when analyzing multi-component signals, as well as inner interference when analyzing mono-component signals with nonlinearly time-varying instantaneous frequencies. Although improved versions, such as Cohen and affine class distributions, help to suppress these interferences, the smoothing process using kernel functions in turn deteriorates time–frequency concentration of auto-terms and may cause extra interferences [18]. In order to overcome the shortcomings of unsatisfactory time–frequency resolution and cross term interferences, Kodera et al. [19] proposed a time–frequency reassignment idea for better readability of short time Fourier transform spectrogram. Auger and Flandrin [20] further extended it to a variety of time–frequency and time–scale distributions. Since the signal energy is distributed within the time–frequency region of a kernel function, we can effectively improve the time–frequency resolution by reallocating the mean energy to the gravity center of the distribution. Reassigned time–frequency energy distribution has fine readability owing to sharp instantaneous frequency curves. However, the time–frequency reassignment is essentially based on a hypothesis of mono-component with linear instantaneous frequency. For multi-component signals, it may suffer from erroneous energy reallocation, because close signal components lead to overlap of cross terms on auto terms on time–frequency plane. Even for mono-component signals but with nonlinear instantaneous frequency, it may still Please cite this article as: X. Chen, Z. Feng, Iterative generalized time–frequency reassignment for planetary gearbox fault diagnosis under nonstationary conditions, Mech. Syst. Signal Process. (2016), http://dx.doi.org/10.1016/j. ymssp.2016.04.023i

X. Chen, Z. Feng / Mechanical Systems and Signal Processing ∎ (∎∎∎∎) ∎∎∎–∎∎∎

3

be susceptible to wrong energy reallocation, because nonlinear instantaneous frequency results in extra inner interferences, thus leading to an erroneous gravity center of the energy distribution. This means that time–frequency reassignment is effective only for mono-component signals with linear instantaneous frequency. For an arbitrary complex signal, if we could decompose it into some constituent mono-components and in the meantime convert their associated instantaneous frequency into linear one, then the merits of time–frequency reassignment could be well exploited. Recently, Olhede and Walden [21] proposed a flexible signal demodulation method termed generalized demodulation that can map any curvilinear instantaneous frequency trajectory into a constant frequency, i.e. a horizontal line parallel to the time axis whereas vertical to the frequency axis on the time–frequency plane. This method provides an approach to linearizing arbitrary nonlinear instantaneous frequency, and meanwhile it makes the target component separable by bandpass filtering. However, the original method applies generalized demodulation only once, thus being unable to decompose arbitrary complex multi-component signals into mono-components of constant frequency. The instantaneous frequency trajectories of constituent frequency components are usually far different from each other. Even though the target component is mapped into a constant frequency after a single application of generalized demodulation, the instantaneous frequency curves of rest components may still be nonlinear and even overlap with others in the frequency domain, thus being not ready for reassignment or even inseparable by conventional bandpass filters. In order to separate an arbitrary complex signal into constituent mono-components and map their instantaneous frequencies into constant ones, thus being ready for time–frequency reassignment, iterative generalized demodulation is necessary [22]. The generalized demodulation method is applied iteratively to the signal, each time with a demodulation function specially designed for a target component. By iterating the generalized demodulation and bandpass filtering, all the time-varying components of interest can be separated as mono-components of constant frequency one by one, thus being ready for best exploiting the merits of time–frequency reassignment. This paper dedicates to integrate the iterative generalized demodulation with conventional time–frequency reassignment methods, termed as iterative generalized time–frequency reassignment, thus overcoming the shortcomings of erroneous energy reallocation and pseudo interferences existing with conventional time–frequency reassignment, and generating a quality time–frequency distribution of higher resolution and free of both cross and inner interferences. We then use this method to pinpoint the time-varying characteristic frequencies from nonstationary vibration signals, thus addressing the issue of planetary gearbox fault diagnosis under nonstationary conditions. Hereafter, this paper is structured as follows. In Section 2, we first evaluate the performance of traditional and iterative generalized time–frequency reassignment via analysis of a synthetic signal, and then elaborate time–frequency reassignment algorithm as well as iterative generalized demodulation. In Section 3, a numerical simulated planetary gearbox vibration signal is generated to illustrate the principle and demonstrate the performance of iterative generalized time–frequency reassignment. In Section 4, the lab experimental vibration signals of a planetary gearbox encompassing both localized and distributed gear fault cases are analyzed to further validate the effectiveness of proposed method in planetary gearbox fault diagnosis under nonstationary conditions. Finally, conclusions are drawn in Section 5.

2. Iterative generalized time–frequency reassignment 2.1. Motivation Time–frequency analysis relies on an accurate instantaneous frequency estimation of each true signal component [16]. Conventional time–frequency analysis methods, including both linear and bilinear time–frequency representations, suffer from either poor time–frequency resolution or both outer and inner interferences, thus may fail to pinpoint the instantaneous frequency finely. In order to improve the time–frequency readability, reassignment method has been proposed [19]. This method is actually based on an assumption of symmetric energy distribution about the true instantaneous frequency. While it is effective in enhancing time–frequency resolution for a simple case, i.e. mono-component signals of linear instantaneous frequency, time–frequency reassignment is not generally valid for arbitrary complex signal, especially for multi-component signals of nonlinear instantaneous frequency. For multi-component signals, if the instantaneous frequency trajectories are close to each other, the time–frequency tilling cell may not be fine enough to pinpoint the true instantaneous frequency from each other, due to the outer interferences. Moreover, the existence of outer interferences may lead to an overall energy distribution asymmetric about the true instantaneous frequency. Even for mono-component signals, if the instantaneous frequency is curvilinear, inner interferences will be resulted, and the energy distribution is no longer symmetric about the true instantaneous frequency. In such cases, the gravity center of the time–frequency distribution does not coincide with but deviate from the true instantaneous frequency. These factors will undoubtedly hinder the time–frequency reassignment from reallocating signal energy distribution to the theoretical location, i.e. along the true instantaneous frequency. To illustrate the above issues existing with the time–frequency reassignment method, a synthetic signal s (t ) composed of two time-varying frequency components f1 (t ) and f2 (t ) is generated

s (t ) = cos [2π

∫ f1 (t ) dt ] +

cos [2π

∫ f2 (t ) dt ],

(1)

Please cite this article as: X. Chen, Z. Feng, Iterative generalized time–frequency reassignment for planetary gearbox fault diagnosis under nonstationary conditions, Mech. Syst. Signal Process. (2016), http://dx.doi.org/10.1016/j. ymssp.2016.04.023i

X. Chen, Z. Feng / Mechanical Systems and Signal Processing ∎ (∎∎∎∎) ∎∎∎–∎∎∎

4

where

⎧ 100t + 260, 0 ≤ t < 0.73 f1 (t ) = ⎨ , ⎩ 333, 0.73 ≤ t < 1

(2)

⎧ 100t + 240, 0 ≤ t < 0.25 ⎪ f2 (t ) = ⎨ 2240t 2 − 2240t + 660 − 48.3 sin (22πt + 0.45π ), 0.25 ≤ t < 0.73. ⎪ 666t − 242.98, 0.73 ≤ t ≤ 1 ⎩

(3)

Fig. 1(a) shows the true theoretical instantaneous frequency curves of the synthetic signal. It embodies three typical challenging issues existing with time–frequency reassignment: (1) more than one close instantaneous frequency trajectories within a time–frequency tiling cell, see region A in Fig. 1(a); (2) crossing between instantaneous frequency trajectories, see 500

500 B 400 A

Frequency [Hz]

Frequency [Hz]

400

300

200

100

300 200

100 C

0

0

0.2

0.4

0.6

0.8

0

1

0

0.2

500

500

400

400

300 200

0

0.2

0.4

0.6

0.8

1

0.6

0.8

1

0.8

1

200

0

1

0

0.2

0.4

Time [s]

500

500

400

400

Frequency [Hz]

Frequency [Hz]

0.8

300

Time [s]

300 200

300 200 100

100 0

0.6

100

100

0

0.4

Time [s]

Frequency [Hz]

Frequency [Hz]

Time [s]

0

0.2

0.4

0.6

Time [s]

0.8

1

0

0

0.2

0.4

0.6

Time [s]

Fig. 1. Synthetic signal analysis. (a) True instantaneous frequency. (b) Short time Fourier transform spectrogram. (c) Reassigned short time Fourier transform spectrogram. (d) Iterative generalized reassigned short time Fourier transform spectrogram. (e) Continuous wavelet transform scalogram. (f) Reassigned continuous wavelet transform scalogram. (g) Iterative generalized reassigned continuous wavelet transform scalogram. (h) Wigner–Ville distribution. (i) Reassigned smoothed pseudo Wigner–Ville distribution. (j) Iterative generalized reassigned smoothed pseudo Wigner–Ville distribution.

Please cite this article as: X. Chen, Z. Feng, Iterative generalized time–frequency reassignment for planetary gearbox fault diagnosis under nonstationary conditions, Mech. Syst. Signal Process. (2016), http://dx.doi.org/10.1016/j. ymssp.2016.04.023i

500

500

400

400

Frequency [Hz]

Frequency [Hz]

X. Chen, Z. Feng / Mechanical Systems and Signal Processing ∎ (∎∎∎∎) ∎∎∎–∎∎∎

300 200

100

0

300 200

100

0

0.2

0.6

0.4

0.8

0

1

0

0.2

0.4

0.6

0.8

1

0.6

0.8

1

Time [s]

500

500

400

400

Frequency [Hz]

Frequency [Hz]

Time [s]

300

200 100

0

5

300 200

100

0

0.2

0.4

0.6

0.8

1

0

Time [s]

0

0.2

0.4

Time [s] Fig. 1. (continued)

region B; (3) curvilinear instantaneous frequency trajectory, see region C. These cases are quite common in real world vibration signals, especially for planetary gearbox vibration signals under nonstationary conditions. We use this synthetic signal to illustrate the limitations of time–frequency reassignment and the principle of our proposed improvement. Short time Fourier transform spectrogram, continuous wavelet transform scalogram and smoothed pseudo Wigner–Ville distribution are three typical time–frequency representations. Fig. 1(b), (c), (e), (f), (h) and (i) show their time–frequency distribution and corresponding reassigned version of the synthetic signal respectively. The original time–frequency distributions roughly extract the time-varying frequency components and their time evolution feature. Their reassignment versions improve the time–frequency resolution and suppress partially outer interferences. But this is only limited to the upper instantaneous frequency within the interval [0.4, 0.6] s. The reassigned time–frequency representations approach to the true upper instantaneous frequency, because the two components are far away from each other and the upper one changes linearly with time. However, for the specific time–frequency structure within the three regions A, B and C, both the original time–frequency representations and their reassignment versions fail to reveal the true time-varying behaviors. For region A, due to the poor time–frequency resolution, the closely spaced instantaneous frequency trajectories are not figured out by short time Fourier transform and continuous wavelet transform, see Fig. 1(b) and (e). In the meantime, some crossings exist in between the two instantaneous frequency trajectories, because of outer interferences. So does the reassigned short time Fourier transform spectrogram and continuous wavelet transform scalogram, as shown in Fig. 1(c) and (f). In the Wigner–Ville distribution, significant outer interferences appear in between the two instantaneous frequency trajectories, and still exist even after reassignment, see Fig. 1(h) and (i). Within region B, the true instantaneous frequency trajectories intersect. This leads to severe outer interferences in the vicinity of the crossing, see Fig. 1(b), (e) and (h). The outer interferences are not completely removed even after time– frequency reassignment, thus resulting in distortion of the true instantaneous frequency curves, see Fig. 1(c), (f) and (i). Within region C, the two instantaneous frequencies are relatively further away from each other. Hence the outer interferences are not so strong as that within regions A or B. However, the lower instantaneous frequency curve is nonlinear, leading to inner interferences, see Fig. 1(b), (e) and (h). Such inner interference further results in deviation of reassigned time–frequency distribution from the true instantaneous frequency curve, see (c), (f) and (i). Real world nonstationary signals are usually multi-component of arbitrary curvilinear instantaneous frequency. For multi-component signals, outer interferences inevitably exist in between neighboring signal components on the time– Please cite this article as: X. Chen, Z. Feng, Iterative generalized time–frequency reassignment for planetary gearbox fault diagnosis under nonstationary conditions, Mech. Syst. Signal Process. (2016), http://dx.doi.org/10.1016/j. ymssp.2016.04.023i

X. Chen, Z. Feng / Mechanical Systems and Signal Processing ∎ (∎∎∎∎) ∎∎∎–∎∎∎

6

Fig. 2. Inner interference effect.

frequency plane. When the signal instantaneous frequency is nonlinear, inner interferences also arise at the inner side of an instantaneous frequency curve on the time–frequency plane, even for a mono-component signal. See the interval [0, 0.75] s in Fig. 2 for an example. The synthetic mono-component signal in Fig. 2 is generated according to Eqs. (4) and (5) as below

s (t ) = cos [500πt + 2π

∫ f (t ) dt ],

(4)

where

⎧ ⎛ 40 ⎞ ⎪ 100 sin ⎜ πt ⎟, 0 ≤ t < 0.75 ⎝ 7 ⎠ f (t ) = ⎨ . ⎪ 100, 0.75 ≤ t < 1 ⎩

(5)

The inner interference will lead to estimation error in the gravity center of time–frequency distribution, thus resulting in deviation of the gravity center from the true instantaneous frequency curve. Fortunately, mono-components of linear instantaneous frequency do not suffer from such unwanted effects, because they are free of inner interference. In this case, the gravity center of time–frequency distribution matches well with the true instantaneous frequency curve, see the interval [0.75, 1] s in Fig. 2. The above simulation analysis shows that time–frequency reassignment is effective only for monocomponent of linear instantaneous frequency, but is ineffective for both multi-component signals and even mono-component signals of nonlinear instantaneous frequency. Therefore, how to further extend time–frequency reassignment to arbitrary complex signals remains a challenging topic. Based on the above analyses, mono-component and symmetric energy distribution about the true instantaneous frequency curve on the time–frequency plane are necessary premises for the best exploitation of time–frequency reassignment. Mono-component signals of linear instantaneous frequency satisfy such a requirement. Therefore, how to decompose complex multi-component signal into mono-components and map their instantaneous frequencies into linear ones are the key to success in extending time–frequency reassignment to arbitrary complex signals. Iterative generalized demodulation is capable of decomposing complex multi-component signals into mono-components with constant frequencies, thus providing an approach to making signals meet the requirement of mono-components with a symmetric time–frequency distribution. 2.2. Iterative generalized demodulation Generalized demodulation can map an arbitrary time-variant frequency component into a constant frequency one [21]. For an arbitrary mono-component signal x (t ) = A (t ) cos [2πφ (t )], if we multiply it with a demodulation function exp [ − j2πυ (t )], then we have s (t ) = x (t ) exp [ − j2πυ (t )], where φ (t ) and υ (t ) are instantaneous phase functions. The Fourier transform of the signal s (t ) ∞

∫−∞ x (t ) exp { − j2π [υ (t ) + ft ]} dt = S [f + υ ̇(t )].

(6)

If exp { − j2π [υ (t ) + ft ] is taken as the Fourier basis function, then the generalized Fourier transform of original signal x (t ) is defined as

SG (f ) =

∞

∫−∞ x (t ) exp { − j2π [υ (t ) + ft ]} dt.

(7)

Therefore, the inverse generalized Fourier transform is derived as Please cite this article as: X. Chen, Z. Feng, Iterative generalized time–frequency reassignment for planetary gearbox fault diagnosis under nonstationary conditions, Mech. Syst. Signal Process. (2016), http://dx.doi.org/10.1016/j. ymssp.2016.04.023i

X. Chen, Z. Feng / Mechanical Systems and Signal Processing ∎ (∎∎∎∎) ∎∎∎–∎∎∎

x (t ) =

∞

7

∞

∫−∞ SG (f ) exp {j2π [υ (t ) + ft ]} df = exp [j2πυ (t )] ∫−∞ SG (f ) exp (j2πft ) df .

(8)

If SG (f ) = δ (f − f0 ), i.e. the time-variant frequency φ̇ (t ) of the signal is converted into a constant frequency f0 , then x (t ) = exp {j2π [υ (t ) + f0 t ]}. Hence comes a simple yet effective way to map a time-variant frequency into a constant one, as long as the time-variant phase function υ (t ) satisfies

φ̇ (t ) = υ ̇ (t ) + f0 ,

(9)

where the instantaneous frequency φ̇ (t ) can be estimated via conventional time–frequency analysis. After the mono-component of interest is mapped into a constant frequency f0 , it can be separated from other constituent components via a bandpass filter centered at f0 . Real world time-varying multi-component signals can be modeled as N

x (t ) =

N

∑ xi (t ) = ∑ ai (t ) cos [2πϕi (t )]. i=1

i=1

(10)

In most cases, the instantaneous frequency of each mono-component has a changing rate different from each other at any instant. For such signals, after generalized demodulation, only one target mono-component can be separated and mapped into a constant frequency, while other constituent components are still left time-varying in the residual signal. Therefore, in order to decompose all the constituent components into mono-components of constant frequency, it is necessary to iteratively apply the generalized demodulation, each time with a specially designed demodulation function. Hence comes the iterative generalized demodulation [22]. Its procedure is summarized briefly as below: (1) For any signal x (t ), estimate the constituent frequency components by conventional time–frequency analysis, such as short time Fourier transform. (2) Estimate the instantaneous frequency fi (t ) of each constituent component xi (t ) by data fitting to its time–frequency ridge. (3) Establish the corresponding generalized demodulation phase function υi (t ) = ∫ [fi (t ) − f0 ] dt , and further construct the demodulation function exp [ − j2πυi (t )], for each target constituent component xi (t ). with the demodulation function exp [ − j2πυi (t )], obtaining (4) Frequency modulate the signal x (t ) y (t ) = x (t ) exp [ − j2πυi (t )], thus mapping the instantaneous frequency of the ith constituent component xi (t ) into a constant frequency f0 . (5) Extract the ith generalized demodulated component yi (t ) using a bandpass filter, with a center frequency of f0 , and a bandwidth equal to the minimum frequency spacing between the instantaneous frequency of the ith component and upper/lower neighboring ones. (6) Repeat steps (2) through (5), until all target frequency components are extracted. Taking the real part of each separated component yi (t ), now they are mono-components of constant frequency and ready for time–frequency reassignment.

2.3. Time–frequency reassignment Both linear and bilinear time–frequency representations suffer from poor time–frequency resolution and outer and/or inner interferences (particularly for Cohen and Affine class distributions of multi-component signals). Time–frequency reassignment aims to improve the time–frequency resolution and reduce the interferences, thus enhancing the readability of time–frequency representations [19,20]. For any signal x (t ), Cohen class distributions are essentially a time–frequency smoothed Wigner–Ville distribution WVDx (t , f ) with a kernel function ϕ (t , f )

TFR x (t , f ) =

+∞

+∞

∫−∞ ∫−∞

ϕ (t − t′, f − f ′) WVDx (t′, f ′) dt′df ′.

(11)

They can be considered as a sum of weighted distributions around the time–frequency location (t , f ). Actually, TFRx (t , f ) is the mean energy in the vicinity of (t , f ). Although the averaging operation can effectively suppress the outer (cross term) interference, it deteriorates the time–frequency concentration of useful auto terms. In order to improve the time–frequency resolution, it is necessary to properly reallocate the mean energy TFRx (t , f ) at location (t , f ) to the gravity center of the distribution. Then the new coordinates become +∞

+∞

−∞

−∞

∫ ∫−∞ τϕ (τ, υ) WVDx (t − τ, f − υ) dτdυ t^ (t , f ) = t − −∞ , +∞ +∞ ∫ ∫ ϕ (τ, υ) WVDx (t − τ, f − υ) dτdυ

(12)

and Please cite this article as: X. Chen, Z. Feng, Iterative generalized time–frequency reassignment for planetary gearbox fault diagnosis under nonstationary conditions, Mech. Syst. Signal Process. (2016), http://dx.doi.org/10.1016/j. ymssp.2016.04.023i

X. Chen, Z. Feng / Mechanical Systems and Signal Processing ∎ (∎∎∎∎) ∎∎∎–∎∎∎

8 +∞

+∞

−∞

−∞

∫ ∫−∞ υϕ (τ, υ) WVDx (t − τ, f − υ) dτdυ ^ f (t , f ) = f − −∞ , +∞ +∞ ∫ ∫ ϕ (τ, υ) WVDx (t − τ, f − υ) dτdυ

(13)

respectively. Summing up the reassigned distributions at new locations, yields the reassigned Cohen class time–frequency distribution

RTFR x (t′, f ′) =

+∞

+∞

∫−∞ ∫−∞

^ TFR x (t , f ) δ [t′ − t^ (t , f )] δ [f ′ − f (t , f )] dt df ,

(14)

where δ (⋅) is the Dirac delta function. Affine class distributions are a time–scale smoothed Wigner–Ville distribution using a kernel function ϕ (c , b)

TSR x (t , a) =

+∞

+∞

∫−∞ ∫−∞

⎛ t − t′ ⎞ ϕ⎜ , af ′⎟ WVDx (t′, f ′) dt′df ′. ⎝ a ⎠

(15)

In the neighborhood of a time–scale location (t , a), the coordinates of gravity center of the affine class distribution are +∞

+∞

−∞

−∞

τ

∫ ∫−∞ τϕ ( a , υ0 − aυ) WVDx (t − τ, f − υ) dτdυ t^ (t , a) = t − −∞ , +∞ +∞ ∫ ∫ ϕ ( aτ , υ0 − aυ) WVDx (t − τ, f − υ) dτdυ

(16)

and

a^ (t , a) =

υ0 ∫

+∞

−∞ +∞

+∞

∫−∞ ϕ ( aτ , υ0 − aυ ) WVDx (t − τ, f − υ) dτdυ +∞

∫−∞ ∫−∞ υϕ ( aτ , υ0 − aυ) WVDx (t − τ, f − υ) dτdυ

. (17)

Then the reassigned affine class time–scale distribution

RTSR x (t′, a′) =

+∞

+∞

∫−∞ ∫−∞

TSR x (t , a) δ [t′ − t^ (t , a)] δ [a′ − a^ (t , a)] dt da.

(18)

For an affine class time–scale distribution, it can be converted into a time–frequency distribution according to the relationship between scale and frequency f = fNyquist /a , where the Nyquist frequency fNyquist equals half of the sampling frequency. A quality reassignment of both Cohen and affine class distributions should have high time–frequency resolution and be free of both outer and inner interferences. According to Eqs. (12–14) and (16–18), the reassignment of both Cohen and affine class distributions highly relies on the coordinates of gravity center. As such, it requires the calculated gravity center of the energy distribution at any time–frequency location coincide exactly with the true instantaneous frequency curve. This implicitly imposes a strict constraint on the signal, i.e. mono-component with a linear instantaneous frequency. For such signals, their time–frequency distribution is symmetric about the true instantaneous frequency. Thus the calculated gravity center coincides with the coordinates of true instantaneous frequency. Otherwise, the reallocated energy distribution will deviate from the true time–frequency location, because: multi-component will result in outer interference, and nonlinear

Fig. 3. Analysis procedure.

Please cite this article as: X. Chen, Z. Feng, Iterative generalized time–frequency reassignment for planetary gearbox fault diagnosis under nonstationary conditions, Mech. Syst. Signal Process. (2016), http://dx.doi.org/10.1016/j. ymssp.2016.04.023i

X. Chen, Z. Feng / Mechanical Systems and Signal Processing ∎ (∎∎∎∎) ∎∎∎–∎∎∎

9

instantaneous frequency will lead to inner interference even for mono-component. For such signals, the outer interferences will mislead identification of true signal components and their instantaneous frequencies, if they are close on the time– frequency plane. Meanwhile, the existence of inner interference will result in an energy distribution asymmetric about the true instantaneous frequency, thus making the calculated gravity center deviate from the true instantaneous frequency. In summary, the outer and inner interferences due to multi-component and nonlinear instantaneous frequency severely disturb accurate calculation of the gravity center, and mislead estimation of the true instantaneous frequency, thus hindering a quality time–frequency reassignment. However, real world signals are often multi-component, and have complex nonlinear instantaneous frequencies. To extend time–frequency reassignment to arbitrary complex signals, we improve it by exploiting the merit of iterative generalized demodulation in decomposing arbitrary complex signal into mono-components of constant frequency. 2.4. Iterative generalized time–frequency reassignment In order to make time–frequency reassignment effective to analyze arbitrary complex signals, iterative generalized demodulation is used to decompose any signal into mono-components of constant frequency. Through such a preprocessing, the requirement of mono-component and linear instantaneous frequency will be met, and the merits of time–frequency reassignment can be best exploited, yielding a quality time–frequency distribution of fine time–frequency resolution and outer and inner interferences free nature. Fig. 3 shows the flowchart of this proposed reassigned time–frequency distribution improved via iterative generalized demodulation. The analysis procedure is detailed as below. (1) Decompose a signal x (t ) into mono-components yi (t ) of a constant frequency f0 , each time using a specifically designed demodulation function according to an estimated instantaneous frequency fi (t ) of the target component based on conventional time–frequency analysis, following the procedure introduced in Section 2.2. (2) Calculate the time–frequency energy distribution TFRyi (t , f ) of each mono-component yi (t ) using time–frequency reassignment method, such as the reassigned version of short time Fourier transform spectrogram, continuous wavelet transform scalogram, and smoothed pseudo Wigner–Ville distribution. (3) Given the estimated instantaneous frequency fi (t ) of original component and the time–frequency representation magnitude TFRyj (t , f ) of associated generalized demodulated component yi (t ), the time–frequency representation of the ith original component xi (t ) can be obtained as TFRxi (t , f ) = TFRyi [t , f + f0 − fi (t )]. (4) Repeat steps (2) through (3), until the time–frequency representation of all constituent components are calculated. Then the time–frequency representation of the original signal x (t ) can be obtained by superposing the time–frequency representation of all components TFRx (t , f ) = ∑i TFRxi (t , f ). The proposed iterative generalized time–frequency reassignment method has not only fine time–frequency resolution, but is also free of both outer and inner interferences. With the aid of iterative generalized demodulation as a preprocessing tool, each constituent component of a complex signal is generalized demodulated into constant frequency and separated into a mono-component. Such a mono-component of constant frequency has neither outer nor inner interferences on the time–frequency plane, so its time–frequency distribution is symmetric about its instantaneous frequency. It exactly meets the requirement by time–frequency reassignment, i.e. the gravity center of energy distribution coincides with the true instantaneous frequency. Therefore, the merits of time–frequency reassignment can be best exploited. The time–frequency distribution of original signal is a linear superposition of the reassigned time–frequency representation of each constituent mono-component, so it inherits the merits from the time–frequency reassignment of constant frequency mono-components, i.e. fine time–frequency concentration as well as both outer and inner interferences free nature. These merits enable the proposed iterative generalized time–frequency reassignment to be a potentially effective approach for analysis of nonstationary complex multi-component signals. We recall the illustration shown by Fig. 1. The iterative generalized reassigned version of typical time–frequency representations, including short time Fourier transform spectrogram, continuous wavelet transform scalogram, and smoothed pseudo Wigner–Ville distribution, all have a higher time–frequency resolution, and are more free of both outer and inner interferences, than their respective original and reassigned versions. However the constituent components behave, either closely spaced within region A, cross over each other within region B, or nonlinearly fluctuating instantaneous frequency within region C, they all exhibit clearly a time–frequency structure exactly consistent with the true one. This demonstrates that the proposed method improves considerably the time–frequency readability.

3. Simulation evaluation In this section, we use a numerical simulated signal to test the performance of proposed iterative generalized time– frequency reassignment method in analyzing nonstationary planetary gearbox vibration signals. The numerical simulated signal is generated according to the planetary gearbox vibration model in [5] and in a sun gear fault case. Without loss of generality, only the fundamental frequencies are considered, i.e., the gear meshing frequency fm (t ), sun gear rotating frequency f s(r ) (t ) and sun gear fault characteristic frequency fs (t ). Then the simulated signal is modeled as Please cite this article as: X. Chen, Z. Feng, Iterative generalized time–frequency reassignment for planetary gearbox fault diagnosis under nonstationary conditions, Mech. Syst. Signal Process. (2016), http://dx.doi.org/10.1016/j. ymssp.2016.04.023i

X. Chen, Z. Feng / Mechanical Systems and Signal Processing ∎ (∎∎∎∎) ∎∎∎–∎∎∎

6

0.35

4

0.3 0.25

2

Amplitude

Amplitude

10

0 -2

0.1

-4 -6 0

0.2 0.15

0.05

0.2

0.4

0.6

0.8

1

Time [s]

0 0

200

400

600

800

1000

Frequency [Hz]

Frequency [Hz]

500

400

300

200

100 0

0.2

0.4

0.2

0.4

Time [s]

0.6

0.8

1

0.6

0.8

1

0.6

0.8

1

Frequency [Hz]

500

400

300

200

100 0

Time [s]

Frequency [Hz]

500

400

300

200

100 0

0.2

0.4

Time [s] Fig. 4. Simulated signal analysis. (a) Waveform. (b) Fourier spectrum. (c) Reassigned short time Fourier transform spectrogram. (d) Iterative generalized reassigned short time Fourier transform spectrogram. (e) Reassigned continuous wavelet transform scalogram. (f) Iterative generalized reassigned continuous wavelet transform scalogram. (g) Reassigned smoothed pseudo Wigner–Ville distribution. (h) Iterative generalized reassigned smoothed pseudo Wigner–Ville distribution.

Please cite this article as: X. Chen, Z. Feng, Iterative generalized time–frequency reassignment for planetary gearbox fault diagnosis under nonstationary conditions, Mech. Syst. Signal Process. (2016), http://dx.doi.org/10.1016/j. ymssp.2016.04.023i

X. Chen, Z. Feng / Mechanical Systems and Signal Processing ∎ (∎∎∎∎) ∎∎∎–∎∎∎

2

1

3

4

6

5

11

7

Fig. 5. Wind turbine drivetrain test rig: (1) motor, (2) tachometer, (3) fixed-shaft gearbox, (4) planetary gearbox stage 1, (5) planetary gearbox stage 2, (6) accelerometers, (7) brake.

x (t ) =

⎡ ⎤ ⎡ − cos ⎣⎢ 2π ∫ f (t ) dt ⎦⎥ } { 1 + A cos ⎣⎢ 2π ∫ f (t ) dt + ϕ} cos {2π ∫ f {1   (r ) s

AM by

sun rotation

m (t ) dt

s

AM by

sun

fault

⎡ ⎤ + B sin ⎣⎢ 2π fs (t ) dt + φ⎦⎥ + θ} + n (t ), 

∫

FM by

(19)

sun fault

where A = 0.8 and B = 0.2 are the AM and FM magnitude respectively, the initial phrases ϕ = φ = θ = 0, and a Gaussian white noise n (t ) at a SNR (signal to noise ratio) of 3 dB is added to mimic background noise. In this simulation, the frequency parameters are set the same as those of the lab planetary gearbox in Section 4.1, see Table 2 for details of gear configuration. To simulate unsteady states, a quadric speed up process for an example, the sun gear rotating frequency is set to fs(r ) = − 12.9t 2 + 20t + 12. Then according to the gear configuration, the gear meshing frequency and the sun gear fault characteristic frequency can be calculated as fm (t ) = (50/3) fs(r ) (t ) and fs (t ) = (10/3) fs(r ) (t ) respectively. Given the above information, a synthetic signal can be generated according to Eq. (17). Fig. 4(a) and (b) show the synthetic signal waveform and its Fourier spectrum respectively. It is difficult to discern the constituent components, because of time-variation of the frequency components and their overlapping in the frequency domain. Time–frequency analysis is a better approach to analysis of such a nonstationary signal. According to the reassigned short time Fourier spectrogram shown in Fig. 4(c), the synthetic signal is composed of nine characteristic frequency components. They are the meshing frequency fm (t ), and its sum or difference combinations with the sun gear rotating frequency f s(r ) (t ) and sun gear fault frequency fs (t ), i.e., fm (t ) ± fs(r ) (t ), fm (t ) ± fs (t ), fm (t ) ± fs (t ) ± fs(r ) (t ). By polynomial data fitting to the time–frequency ridges of nine frequency components, their instantaneous frequencies fi (t ) can be estimated. Integrating fi (t ) over time t , yields the corresponding instantaneous phase 2πυi (t ). Then the generalized demodulation function exp [ − j2πυi (t )] can be constructed for each instantaneous frequency fi (t ). Following the proposed iterative generalized time–frequency reassignment procedure, we obtain the iterative generalized reassigned version of short time Fourier spectrogram, continuous wavelet transform scalogram and smoothed pseudo Wigner–Ville distribution, as shown in Fig. 4(d), (f) and (h) respectively. They all have a fine time–frequency resolution and are free of both outer and inner interferences. The good time–frequency readability enables better extraction of the timevarying sidebands including both the frequency value and sideband spacing at any time instant. All the nine instantaneous frequency curves are clearly shown, exactly matching the nine theoretical components fm (t ) + fFM (t ), fm (t ) ± fs(r ) (t ) + fFM (t ), fm (t ) ± fs (t ) + fFM (t ) and fm (t ) ± fs (t ) ± fs(r ) (t ) + fFM (t ), where fFM (t ) = B cos [2π ∫ fs (t ) dt ] fs (t ) stands for the oscillating frequency due to gear fault FM effect. The small instantaneous frequency ripple around each component fFM (t ) and its variation with the time-varying speed are also clearly exhibited, as displayed in the close-up view at top right corner of Fig. 4(d), (f) and (h). On the other hand, the reassigned time–frequency representations of short time Fourier spectrogram, continuous Table 1 Fixed-shaft gearbox configuration. Gear

Input Intermediate Output

Number of gear teeth Stage 1

Stage 2

32 80 –

– 40 72

Please cite this article as: X. Chen, Z. Feng, Iterative generalized time–frequency reassignment for planetary gearbox fault diagnosis under nonstationary conditions, Mech. Syst. Signal Process. (2016), http://dx.doi.org/10.1016/j. ymssp.2016.04.023i

X. Chen, Z. Feng / Mechanical Systems and Signal Processing ∎ (∎∎∎∎) ∎∎∎–∎∎∎

12

Table 2 Planetary gearbox configuration. Gear

Number of gear teeth

Ring Planet Sun

Stage 1

Stage 2

100 40(4) 20

100 36(4) 28

Note: The number of planet gears in the parenthesis.

Fig. 6. Sun gear damage. (a) stage 1 sun gear wear. (b) stage 2 sun gear chipping.

Table 3 Planetary gearbox characteristic frequencies. Frequency

Stage 1

Stage 2

Meshing fm(t) (r) Sun rotating fs (t) Planet carrier rotating fc(t) Sun fault fs(t) Planet fault fp(t) Ring fault fr(t)

(100/27) fd(t) (2/9) fd(t) (1/27) fd(t) (20/27) fd(t) (5/54) fd(t) (4/27) fd(t)

(175/216) fd(t) (1/27) fd(t) (7/864) fd(t) (175/1512) fd(t) (175/7776) fd(t) (7/216) fd(t)

wavelet transform scalogram and smoothed pseudo Wigner–Ville distribution, as shown by Fig. 4(c), (e) and (g), roughly reveal the trend of nine frequency components, but mistake the small instantaneous frequency ripple fFM (t ) as irregular texture in between neighboring frequency components, because of imperfect time–frequency readability. The above analysis clearly demonstrates the effectiveness of the proposed method in identifying the complex time-varying signal components. This shows its potential in extracting the time-varying fault characteristic frequencies of planetary gearbox from vibration signals under nonstationary conditions.

4. Experimental evaluation 4.1. Experiment setting Experiments are carried out on a wind turbine drivetrain test rig in University of Ottawa lab, to further validate the proposed iterative generalized time–frequency reassignment analysis method. The test rig encompasses both a two stage Please cite this article as: X. Chen, Z. Feng, Iterative generalized time–frequency reassignment for planetary gearbox fault diagnosis under nonstationary conditions, Mech. Syst. Signal Process. (2016), http://dx.doi.org/10.1016/j. ymssp.2016.04.023i

X. Chen, Z. Feng / Mechanical Systems and Signal Processing ∎ (∎∎∎∎) ∎∎∎–∎∎∎

8

0.06 0.05

4

Amplitude [m/s 2]

Amplitude [m/s 2]

6

2 0 -2 -4

0.04 0.03 0.02 0.01

-6 -8 0

13

1

2

3

4

0 0

5

50

100

150

Time [s]

200

250

300

350

400

Frequency [Hz]

400

65

fm1+ (7/4) fs1

Frequency [Hz]

Speed [Hz]

60

55

200

100

50

45 0

300

1

2

Time [s]

3

4

fm1- fs1

0 0

5

1

2

Time [s]

3

4

5

Frequency [Hz]

80

fd

60

40 (r) f s2

20

0 0

1

2

Time [s]

3

4

5

Fig. 7. Baseline signal. (a) Waveform. (b) Fourier spectrum. (c) Drive motor speed. (d) Iterative generalized reassigned time–frequency representation 0– 400 Hz. (e) Iterative generalized reassigned time–frequency representation 0–80 Hz.

planetary gearbox and a two stage fixed-shaft gearbox, as shown in Fig. 5. The power flow path is: drive motor-fixed shaft gearbox-planetary gearbox-electromagnetic powder brake. Gear configuration of the fixed-shaft gearbox and the planetary gearbox are listed in Tables 1 and 2 respectively. Accelerometers are mounted on the top of planetary gearbox casing to measure vibration signals at a sampling frequency of 20,000 Hz. The electromagnetic powder brake applies a 16.3 N∙m external load to the output shaft connecting the planet carrier of planetary gearbox stage 2. To simulate the nonstationary conditions, the drive motor is set to speed down from 60 Hz to 50 Hz in 5 s. The corresponding wind turbine blade speed reduces continuously from 0.486 Hz to 0.405 Hz, covering the rated speed 0.417 Hz of typical wind turbines in real world. In the meantime, the drive motor speed is recorded by an encoder at a sampling frequency of 20 Hz. To simulate typical localized as well as distributed gear faults, two types of gear damage are introduced to the sun gears of planetary gearbox stage 1 and stage 2 respectively, while other gears are healthy. Fig. 6(a) and (b) show the sun gear fault. Please cite this article as: X. Chen, Z. Feng, Iterative generalized time–frequency reassignment for planetary gearbox fault diagnosis under nonstationary conditions, Mech. Syst. Signal Process. (2016), http://dx.doi.org/10.1016/j. ymssp.2016.04.023i

X. Chen, Z. Feng / Mechanical Systems and Signal Processing ∎ (∎∎∎∎) ∎∎∎–∎∎∎

14

0.06

8

0.05

4

Amplitude [m/s 2 ]

Amplitude [m/s 2]

6

2 0 -2 -4

0.03 0.02 0.01

-6 -8 0

0.04

1

2

3

4

0 0

5

50

100

65

400

60

300

200

250

fm1+ (7/4) fs1

55

300

350

400

fm1+ 3 fs1

200

100

50

(r) fm1-4 f s1

45 0

150

Frequency [Hz]

Frequency [Hz]

Speed [Hz]

Time [s]

1

2

3

Time [s]

4

5

0 0

1

fm1- fs1

fm1

2

3

4

5

Time [s]

Fig. 8. Stage 1 sun gear wear signal. (a) Waveform. (b) Fourier spectrum. (c) Drive motor speed. (d) Iterative generalized reassigned time–frequency representation.

For stage 1 sun gear, wear is introduced to every gear tooth. For stage 2 sun gear, chipping is created to one gear tooth. Three sets of experiment are carried out, i.e. baseline when all gear are healthy, wear and chipping. Given the drive motor speed and gear configuration, the characteristic frequency at any time instant can be calculated accordingly [5]. For ease of characteristic frequency identification, we list the calculation equation of each characteristic frequency in terms of the drive motor speed fd (t ), in Table 3. At the instant t = 0, the drive motor speed fd (t ) ¼60 Hz, the (r ) characteristic frequencies are calculated as fm1 (t ) ¼222.2 Hz, f s1 (t ) ¼13.3 Hz, and fs1 (t ) ¼44.4 Hz for planetary gearbox stage (r ) 1, and fm2 (t ) ¼ 48.6 Hz, f s2 (t ) ¼ 2.2 Hz, and fs1 (t ) ¼6.9 Hz for planetary gearbox stage 2. 4.2. Baseline signal analysis In baseline case, all gears are healthy. Fig. 7(a), (b) and (c) show the signal waveform, its Fourier spectrum, and the corresponding drive motor speed. According to Table 3, the maximum meshing frequency is 222.2 Hz for planetary gearbox stage 1, and 48.6 Hz for stage 2. Since gear faults are usually carried and manifested by the sidebands around the meshing frequency and its harmonics, we focus on 0–400 Hz and 0–80 Hz, which cover 1.5 times the meshing frequency of stage 1 and stage 2 respectively, to find the gear fault feature of stage 1 and 2. Fig. 7(d) and (e) present the iterative generalized reassigned time–frequency representation (of short time Fourier transform spectrogram, since all the iterative generalized reassigned time–frequency representations have equal performance) of the baseline signal within the frequency band of 0–400 Hz and 0–80 Hz respectively. In Fig. 7(d), targeting for stage 1, two dominant frequency components fm1 − fs1 and fm1 + (7/4) fs1 are present as sidebands centered around the stage 1 meshing frequency fm1. In Fig. 7(e), focusing on stage 2, the present frequency components correspond to the drive motor rotating frequency fd and the stage 2 sun gear rotating frequency fs2(r ). The presence of these frequency components is reasonable, because manufacturing errors and minor defects are inevitable. 4.3. Detection of sun gear wear Fig. 8(a), (b) and (c) present the sun gear wear signal waveform, its Fourier spectrum and corresponding drive motor speed. We concentrate on the frequency band of 0–400 Hz to examine time-varying sidebands around the stage 1 meshing Please cite this article as: X. Chen, Z. Feng, Iterative generalized time–frequency reassignment for planetary gearbox fault diagnosis under nonstationary conditions, Mech. Syst. Signal Process. (2016), http://dx.doi.org/10.1016/j. ymssp.2016.04.023i

X. Chen, Z. Feng / Mechanical Systems and Signal Processing ∎ (∎∎∎∎) ∎∎∎–∎∎∎

15

0.35

1

Amplitude [m/s 2]

Amplitude [m/s 2]

0.3 0.5

0

-0.5

0.25 0.2 0.15 0.1 0.05

-1 0

1

2

3

4

0 0

5

10

20

80

60

60

Frequency [Hz]

Speed [Hz]

65

55

40

50

60

70

80

fm2+ 2 fs2 40

fd

(r) fm2- 3[ f s2 + fs2]

20

50

45 0

30

Frequency [Hz]

Time [s]

1

2

3

Time [s]

4

5

0 0

1

2

3

4

5

Time [s]

Fig. 9. Stage 2 sun gear chipping signal. (a) Waveform. (b) Fourier spectrum. (c) Drive motor speed. (d) Iterative generalized reassigned time–frequency representation.

frequency, because the wear is on the sun gear of stage 1. Fig. 8(d) shows the iterative generalized reassigned time–frequency representation. In addition to the two components fm1 − fs1 and fm1 + (7/4) fs1 present in the baseline signal as shown in Fig. 7(d), more sidebands are clearly identified, including the stage 1 meshing frequency fm1, and its sum and difference combinations with the sun gear rotating frequency and sun gear fault frequency of stage 1, such as fm1 − 4fs1(r ) and fm1 + 3fs1. These extra frequency components all relate to the stage 1 sun gear fault, thus indicating fault existence on the stage 1 sun gear. 4.4. Detection of sun gear chipping Fig. 9 displays the analysis result of the stage 2 sun gear chipping signal. In this case, we focus on the frequency band of 0–80 Hz, which covers 1.5 times the maximum stage 2 meshing frequency (48.6 Hz). In the iterative generalized reassigned time–frequency representation, Fig. 9(d), the drive motor rotating frequency component fd is present, the same as in the baseline signal, see Fig. 7(d). Besides, two new frequency components emerge. They are the stage 2 meshing frequency plus or minus the stage 2 sun gear rotating frequency and fault characteristic frequency harmonics or their combinations, present (r ) as fm2 − 3 [f s2 + fs2 ] and fm2 + 2fs2 . These new frequencies link to the stage 2 sun gear fault, implying fault occurrence on the stage 2 sun gear. Although planetary gearbox vibration signals are very complex, and become even more intricate due to time-varying sidebands under nonstationary conditions, the time-varying sun gear fault characteristic frequency components are successfully extracted, thanking to the good readability of iterative generalized reassigned time–frequency representation. All the findings are consistent with the real experimental setting. This demonstrates the effectiveness of the proposed method in detecting planetary gearbox fault under nonstationary conditions.

5. Conclusions A quality time–frequency analysis method is necessary to extract planetary gearbox fault features under unsteady conditions, since nonstationary planetary gearbox vibration signals feature complex time-varying sidebands. To this end, we Please cite this article as: X. Chen, Z. Feng, Iterative generalized time–frequency reassignment for planetary gearbox fault diagnosis under nonstationary conditions, Mech. Syst. Signal Process. (2016), http://dx.doi.org/10.1016/j. ymssp.2016.04.023i

16

X. Chen, Z. Feng / Mechanical Systems and Signal Processing ∎ (∎∎∎∎) ∎∎∎–∎∎∎

have proposed an iterative generalized time–frequency reassignment method, by exploiting the uniqueness of iterative generalized demodulation to decompose nonstationary multi-component signals into mono-components of constant frequency, thus meeting the requirement of mono-component with linear instantaneous frequency by time–frequency reassignment and further utilizing its merit in improving time–frequency readability. The proposed method has fine time– frequency resolution, and is free of both outer and inner interferences, thus being suited for complex nonstationary signal analysis. This has been validated via both numerical simulated and lab experimental signals of planetary gearboxes under nonstationary conditions. The time-varying fault frequency components of both localized and distributed gear fault are successfully identified.

Acknowledgments This work is supported by National Natural Science Foundation of China (11272047, 51475038), Program for New Century Excellent Talents in University (NCET-12-0775), and Natural Sciences and Engineering Research Council of Canada (RGPIN 121433-2011).

References [1] V.K. Ambarisha, R.G. Parker, Nonlinear dynamics of planetary gears using analytical and finite element models, J. Sound Vib. 302 (2007) 577–595. [2] W.D. Mark, J.A. Hines, Stationary transducer response to planetary-gear vibration excitation with non-uniform planet loading, Mech. Syst. Signal Process. 23 (2009) 1366–1381. [3] M. Inalpolat, A. Kahraman, A theoretical and experimental investigation of modulation sidebands of planetary gear sets, J. Sound Vib. 323 (2009) 677–696. [4] M. Inalpolat, A. Kahraman, A dynamic model to predict modulation sidebands of a planetary gear set having manufacturing errors, J. Sound Vib. 329 (2010) 371–393. [5] Z. Feng, M.J. Zuo, Vibration signal models for fault diagnosis of planetary gearboxes, J. Sound Vib. 331 (2012) 4919–4939. [6] P.D. Samuel, D.J. Pines, Constrained adaptive lifting and the CAL4 metric for helicopter transmission diagnostics, J. Sound Vib. 319 (2009) 698–718. [7] Z. Feng, M. Liang, Y. Zhang, S. Hou, Fault diagnosis for wind turbine planetary gearboxes via demodulation analysis based on ensemble empirical mode decomposition and energy separation, Renew. Energy 47 (2012) 112–126. [8] Z. Feng, M.J. Zuo, J. Qu, et al., Joint amplitude and frequency demodulation analysis based on local mean decomposition for fault diagnosis of planetary gearboxes, Mech. Syst. Signal Process. 40 (2013) 56–75. [9] T. Barszcza, R.B. Randall, Application of spectral kurtosis for detection of a tooth crack in the planetary gear of a wind turbine, Mech. Syst. Signal Process. 23 (2009) 1352–1365. [10] Y. Lei, Z. Liu, et al., Health condition identification of multi-stage planetary gearboxes using a mRVM-based method, Mech. Syst. Signal Process. 60–61 (2015) 289–300. [11] L. Hong, J.S. Dhupia, A time domain approach to diagnose gearbox fault based on measured vibration signals, J. Sound Vib. 333 (2014) 2164–2180. [12] H. Sun, Y. Zi, Z. He, et al., Customized multiwavelets for planetary gearbox fault detection based on vibration sensor signals, Sensors 13 (2013) 1183–1209. [13] W. Bartelmus, R. Zimroz, Vibration condition monitoring of planetary gearbox under varying external load, Mech. Syst. Signal Process. 23 (2009) 246–257. [14] W. Bartelmus, R. Zimroz, A new feature for monitoring the condition of gearboxes in non-stationary operation conditions, Mech. Syst. Signal Process. 23 (2009) 1528–1534. [15] R. Zimroz, A. Bartkowiak, Two simple multivariate procedures for monitoring planetary gearboxes in non-stationary operating conditions, Mech. Syst. Signal Process. 38 (2013) 237–247. [16] Z. Feng, M. Liang, F. Chu, Recent advances in time–frequency analysis methods for machinery fault diagnosis: A review with application examples, Mech. Syst. Signal Process. 38 (2013) 165–205. [17] F. Hlawatsch, G.F. Boudreaux-Bartels, Linear and quadratic time–frequency signal representations, IEEE Signal Process. Mag. 9 (1992) 21–67. [18] L. Cohen, Time–frequency distributions: a review, Proc. IEEE 77 (1989) 941–981. [19] K. Kodera, R. Gendrin, C. de Villedary, Analysis of time-varying signals with small BT values, IEEE Trans. Acoust. Speech Signal Process. 26 (1978) 64–76. [20] F. Auger, F. Flandrin, Improving the readability of time-frequency and time-scale representations by the reassignment method, IEEE Trans. Signal Process. 43 (1995) 1068–1089. [21] S. Olhede, A.T. Walden, A generalized demodulation approach to time-frequency projections for multi component signals, Proc. R. Soc. A–Math. Phys. Eng. Sci. 461 (2005) 2159–2179. [22] Z. Feng, F. Chu, M.J. Zuo, Time–frequency analysis of time-varying modulated signals based on improved energy separation by iterative generalized demodulation, J. Sound Vib. 330 (2011) 1225–1243.

Please cite this article as: X. Chen, Z. Feng, Iterative generalized time–frequency reassignment for planetary gearbox fault diagnosis under nonstationary conditions, Mech. Syst. Signal Process. (2016), http://dx.doi.org/10.1016/j. ymssp.2016.04.023i