- Email: [email protected]
APPLICATION OF WAVELETS TO GEARBOX VIBRATION SIGNALS FOR FAULT DETECTION
Journal of Sound and Vibration (1996) 192(5), 927–939
APPLICATION OF WAVELETS TO GEARBOX VIBRATION SIGNALS FOR FAULT DETECTION W. J. W Department of Mechanical and Manufacturing Engineering, De Montfort University, Leicester, England
P. D. MF Department of Engineering Science, University of Oxford, Oxford OX1 3PJ, England (Received 19 May 1994, and in final form 11 September 1995) The wavelet transform is used to represent all possible types of transients in vibration signals generated by faults in a gearbox. It is shown that the transform provides a powerful tool for condition monitoring and fault diagnosis. The vibration signal from a helicopter gearbox is used to demonstrate the application of the suggested wavelet by a simple computer algorithm. The major advantage of the wavelet transform for analyzing the signal is that it possesses multi-resolutions for localizing short-time components so that all possible types of gear faults can be displayed by a single time–scale distribution resulting from the transform. 7 1996 Academic Press Limited
1. INTRODUCTION
To determine the condition of an inaccessible gear in an operating machine, the vibration signal of the machine can be measured at some convenient location on the casing, together with a signal which is synchronous with the rotation of a shaft or gear in the machine. The time domain average of the vibration of an individual gear can then be extracted from the original vibration signal, representing in the time domain the vibration of that gear alone over one revolution [1, 2]. Although the time domain average of the vibration of a gear contains all the information describing the gear condition, clear symptoms of any damage may not be directly visible in the time domain average, particularly in the early stages. Further processing of the time domain average may be necessary to detect and diagnose the first signs of gear damage, and a wide variety of different techniques have been explored over the years. Most techniques have sought to represent the gear vibration signal in either the time domain or the frequency domain, with each approach having its advantages and disadvantages. However, a more recent trend has been toward representation in the time–frequency or time–scale domain, seeking to combine the advantages of both approaches. A time–frequency distribution describes simultaneously when a signal component occurs and how its frequency spectrum develops with time [3–5]. One important time–frequency representation, the windowed Fourier transform, gives a constant resolution in the time and frequency domains because the window width in the time domain is fixed. Reducing 927 0022–460X/96/200927 + 13 $18.00/0
7 1996 Academic Press Limited
928
. . . .
the window width will increase the time resolution, but will reduce the frequency resolution, because under the uncertainty principle, a higher time and frequency resolution cannot be achieved simultaneously [5]. This may not be a problem when only one resolution is required to display the components of interest. However, if resolutions are required to be different to suit a variety of signal components of different duration, for example in fault diagnosis where suitable resolutions may be not known in advance, this fixed-window time–frequency distribution does not meet the requirement of displaying all possible sizes of components. Thus the windowed Fourier transform is really only suited to the analysis of signals where all of the patterns appear at approximately the same size [6]. In the time–frequency representation, the fixed resolution represents only one of all the possible sizes. The features of other sizes in the signal cannot be expressed with reasonably high resolutions. For example, in the existing time–frequency representation known as the Gabor spectrogram [5, 7], the recommended width of the window for the detection of gear tooth damage is related to the meshing period of one tooth. This fixed width may enable representation of certain sizes of local damage on a tooth, but may not be suitable for other faults, such as the long duration of distributed faults, or the very short duration of early cracks. The recently developed wavelet transform uses a changing width function so that not one but a series of resolutions can be achieved in a single time–scale display. Unlike the Fourier transform, in which a signal is decomposed on to a sinusoid function basis, the wavelet transform uses a more general function as the basis. This produces, on the one hand, more comprehensive transform results but, on the other hand, a variety of possible explanations. The selection of a wavelet basis thus becomes a major issue in the application of the wavelet transform. With the rapid theoretical progress being made in the wavelet transform over the past ten years, its applications have rapidly penetrated areas of one-dimensional and two-dimensional signal analysis in many fields of science and engineering. It is suggested here that if the property of the time–scale localization of the wavelet transform is applied to the analysis of local phenomena in vibration signals, such as the short-time variations caused by faults, it may be possible to give time and scale information for all of the possible sizes of features so as to find the faults and assess their nature from a single three-dimensional representation [8, 9]. The vibration signal picked up from an operating machine has long been used to diagnose the condition of an inaccessible moving component inside the machine. Of course, this vibration is assumed to be associated with the component, and normally the component is one of the vibration generators. Therefore, any damage, such as a fatigue crack, will cause a relatively short-time transient in the vibration. In the early stages, this change is small, so that a high detection sensitivity is required. In gear fault diagnosis, although the information about the gear condition can be conveyed by the vibration signal transmitted through the gearbox casing, and the time domain average [1, 2] produces a high signal-to-noise ratio for the gear in question, clear symptoms may not be directly observed from the average and more information about its nature is very hard to determine directly from the averaged time history. Time–frequency distributions have already been applied to detect faults, and have been proved very useful [5]. However, in more general cases, many faults with different durations are possible in a mechanical component. Faults in a gear, such as wear and spalling, may occur with many different sizes. Since the wavelet transform [6] uses a series of sizes of windows to compare with all sections of the signal, it is therefore possible to display them all simultaneously.
929
2. WAVELET TRANSFORM
To overcome the limitation of the fixed resolution of the windowed Fourier transform in the time and frequency domains, Grossmann and Morlet [10] decomposed the signal x(t) into a family of functions which are the translation and dilation of a unique-valued function c(t), to give the wavelet transform which is defined by WTx (t, s) =
g
a
x(t)zsc(s(t − t)) dt,
(1)
−a
where s is a scaling factor which produces dilation, t is the time or some other spatial co-ordinate, and c(t) is called a wavelet. The corresponding wavelet family is generated by translation of the wavelet in the time domain and dilation in the scale domain, given by (zsc(s(t − t)))(t,s) $ R2 ,
(2)
where R denotes the set of real numbers. Any function in L2(R), the vector space of measurable square integrable one-dimensional functions x(t), can be characterized by its decomposition of the wavelet family of equation (2). A wavelet transform can be interpreted as a decomposition of a signal into a set of frequency channels. A process of dilation in a wavelet family is illustrated in Figure 1. As the scaling factor increases, from bottom to top, the width of the wavelet envelope decreases. In Figure 2 are shown, from left to right, the corresponding Fourier spectra for the wavelet family. As the scaling factor increases, both the centre frequency and the width of the frequency band increase. By using all of the scales, the full frequency range can be covered, from a low frequency to the required cut-off frequency.
Figure 1. Wavelets in the time domain. (a) Real part; (b) imaginary part.
. . . .
930
Figure 2. Wavelets in the frequency domain; spectra at different scales.
The reconstruction of the signal x(t) from the wavelet transform is x(t) =
1 Cc
g g a
−a
a
WTx (t, s)zsc(s(t − t)) ds dt,
(3)
0
where Cc =
g
a
0
=C( f )=2 df f
(4)
and C( f ) is the Fourier transform of c(t). In order to reconstruct x(t), Cc must satisfy Cc Q +a.
(5)
An important particular case of the discrete wavelet transform is that some wavelets c(t) exist such that the family z2 jc(2 j(t − 2−jn))( j,n) $ Z2 ,
(6) 2
where Z denotes the set of integers, is an orthonormal basis of L (R) [8, 11, 12]. The orthonormal wavelet transform is defined as WTx (n, j) =
g
a
x(t)z2 jc(2 j(t − 2−jn)) dt.
(7)
−a
For orthogonal wavelets, the reconstruction of signal x(t) is then x(t) = s s WTx (n, j)z2 jc(2 j (t − 2−jn)).
(8)
j$Z n$Z
Orthogonal wavelets, such as the well-known Daubechies wavelet series, offer faster algorithms and no redundancy in the decomposition [8, 12]. The original signal can be reconstructed by a minimum number of wavelet terms. That is, the original signal is decomposed into a compact wavelet series, which is advantageous in signal reconstruction. However, due to its special shape and limited number of scales, a single wavelet amplitude
931
map has not enough scales to describe all the details, both large and small, of the signal. Since N = 2M samples normally provide M + 1 different widths of wavelets, all scales of transients in the analyzed signal will be scaled only into the nearest size of wavelet. In gear damage detection, a small change in the vibration signal may ‘‘miss the chance’’ of matching the wavelet of the same size during the translation and dilation, as the locations of the orthogonal wavelets are limited, and there is no intermediate scale and location provided, so that no wavelet matches well the features of the signal. Another problem of using the orthogonal wavelets is the time translation variant. That is, the same transients at different time may be shown as different patterns. Therefore, using the non-orthogonal wavelet transform is preferred for the time–scale representation of short-time features, avoiding the mentioned problems. The wavelet transform means that the function x(t) is characterized by its decomposition of a wavelet family with a series of different frequency bandwidths. This can be explained by working out the bandwidth of the wavelet family. The Fourier transform of the wavelet family zsc(st) is Cs ( f ) = (1/zs)C( f/s).
(9)
Let f0 be the centre of the passband of C( f ). Then
g
a
g
a
( f − f0 )=C( f )=2 df = 0,
0
so that f0 =
>g
f =C( f )=2 df
0
a
=C( f )=2 df.
(10)
0
For the wavelet family, one finds, by substituting equation (9) into (10), that the centre frequency of the passband is
g
a
f =Cs ( f )=2 df
0
>g
a
=Cs ( f )=2 df = sf0 .
(11)
0
The root mean square bandwidth of the wavelet around f0 is given by sf , where sf2 =
g
a
( f − f0 )2 =C( f )=2 df.
(12)
0
With C( f ) replaced by Cs ( f ) and f0 by sf0 in the above equation, the root mean square bandwidth is seen to be ssf . That is, the bandwidth of the wavelet is proportional to the scaling factor s. Thus the wavelet transform produces a series of frequency channels with increasing bandwidth as the scale s changes.
3. DAMAGE DETECTION
Early damage to a gear tooth usually causes a variation in the associated vibration signal over a short time, initially less than one tooth meshing period, taking the form of modulated or unmodulated oscillation. In later stages, the duration of the abnormal variation becomes longer, lasting more than one tooth meshing period. Other distributed faults, such as eccentricity and wear, may cover the most part of the whole revolution of the gear.
. . . .
932
To suit these features, the wavelet c(t) can be chosen to be a Gaussian-enveloped oscillation, given by c(t) = c exp(−s 2t 2 − i2pf0 t),
(13)
where c, s and f0 are positive parameters, s determines the width of the wavelet and hence the width of the frequency band, and f0 is the frequency of oscillation which is the centre frequency of the band. The parameter c can usually be unity. It is to be expected that by changing the scaling factor in its corresponding wavelet family, this function will be suitable for the modelling of all sizes of oscillating transients in a gear vibration signal. The Fourier transform of c(t) is given by C( f ) = (czp/s) exp(−(p/s)2( f + f0 )2 ).
(14)
It can be seen that the frequency of oscillation f0 appears at the centre of the frequency band of the wavelet. From the definition of equation (1), the wavelet transform is WTx (t, s) =
g
a
x(t)czs exp(−s 2s 2(t − t)2 − i2pf0 s(t − t)) dt.
(15)
−a
For transients, the wavelet c(t) can also be chosen from other functions, such as a symmetrical Gaussian-enveloped oscillating function, given by c(t) = c exp(−s 2t 2 ) cos (2pf0 t)
(16)
C( f ) = (czp/s)[exp(−(p/s)2( f + f0 )2 ) + exp(−(p/s)2( f − f0 )2 )].
(17)
with a Fourier transform of
The wavelet family czs exp(−s 2s 2(t − t)2 − i2pf0 s(t − t)) is illustrated in Figure 1, which models the oscillating transients with all oscillating frequencies as the scale changes. At the same time, the width of the wavelet is reduced as the scale becomes larger. The Fourier spectra for all scales of wavelets are shown in Figure 2. When a signal is decomposed on to one of these wavelet bases, large values will be obtained if there is a transient with the same shape as the wavelet. The wavelet family of equation (13) is well suited to the detection of the transients in the signal, such as the transients caused in a vibration signal by a cracked gear tooth.
4. COMPUTER IMPLEMENTATION
For gear damage diagnosis, the function x(t) can be either the time domain average of the vibration of a gear or the corresponding analytic function, both of which are periodic with the gear rotation. The selected wavelet of equation (16) satisfies c(−t) = c (t), so that, by using the convolution theorem, equation (15) can be written as WTx (t, s) = (1/zs)F−1 (X( f ) C ( f/s)),
(18)
where F represents the Fourier transform, and C ( f ) represents F(c (t)). For the Gaussian enveloped oscillating wavelet in equation (13), the wavelet transform can be written as WTx (t, s) = (c/zss)F−1 (X( f ) exp(−(p/s)2(( f/s) − f0 )2 )).
(19)
Thus the calculation of the wavelet transform can be implemented by an inverse Fourier transform, taking advantage of the efficient fast Fourier transform (FFT) algorithm.
933
The Fourier transform of the wavelet family can be written as (c/zss) exp(−(p/ss)2( f − sf0 )2 ),
(20)
which implies that the centre frequency of the wavelet is sf0 , proportional to s. Since the half-power bandwidths in the time and frequency domains are tB = 1·18/s and fB = 0·375s respectively, then at scaling factor s, they are therefore 1·18/ss and 0·375ss respectively. The narrower the wavelet is in the time domain, the wider it is in the frequency domain. This increases the frequency overlap between the wavelets of adjacent high scales, which brings about redundancy. Consider a new discrete scaling factor s = 2kl,
k = 0, 1, 2, . . . , M − 1,
(21) M
where M is an integer. For N samples of signal, M can be defined by 2 = N. The centre frequency and half-power bandwidth of the kth wavelet become f0 (k) = 2klf0 ,
fB (k) = 0·375 × 2kls.
(22, 23)
If the cut-off frequency is fC , the analysis frequency range then will be [( f0 − 0·187s), fC ]. At the scaling factor k = M − 1, the cut-off frequency is fC = 2(M − 1)lf0 + 0·187 × 2(M − 1)ls.
(24)
Let adjacent frequency bands overlap at half-power points, so that 2klf0 + 0·187 × 2kls = 2(k + 1)lf0 − 0·187 × 2(k + 1)ls,
k = 0, 1, 2, . . . , M − 1,
which can be simplified to f0 + 0·187s = 2lf0 − 0·187 × 2ls.
(25)
In Figure 3, the top and bottom curves show upper and lower frequencies, and the middle curve is the centre frequency. From equation (25) f0 = 0·187{(2l + 1)/(2l − 1)}s.
Figure 3. Frequency bands of wavelets at different scaling factors.
(26)
. . . .
934
From equations (22), (23) and (26), the lower and upper frequencies of each band or channel are, respectively, fL (k) = f0 (k) − 12 fB (k) = 0·375
2kl s, 2 −1
fU (k) = f0 (k) + 12 fB (k) = 0·375
l
2(k + 1)l s. 2l − 1
Therefore the scaling factor s = 2kl brings about l (for example 12 or 13 ) octave frequency bands, so that fU (k) = 2lfL (k).
(27)
For a given cut-off frequency, there are four parameters, f0 , s, l and M, to be adjusted. Two independent equations, (24) and (25), can be used to determine two parameters given the values of the others. In total, there are 4 × 3/(2 × 1) = 6 combinations of situation, but probably only two of these are likely. Firstly, M and l can be set to find f0 and s. By considering equations (24) and (26), f0 and s can be found as f0 = {(2l + 1)/2Ml + 1}fC ,
s = 5·35{(2l − 1)/2Ml + 1}fc .
(28, 29)
Secondly, l and f0 can be set to find s and M. Parameters f0 and s can be found as s = 5·35
(2l − 1) f, (2l + 1) 0
M=
0
1
f 2l + 1 1·44 ln C l + 1 + 1. l f0 2
(30, 31)
Upon assuming that N samples exactly cover one period of the signal, the corresponding discrete form of equation (18) is given by the circular convolution WTx (n, k) =
N−1
1 z2
kl
s X(r)C (2−klr) exp(2pirn/N),
r=0
n = 0, 1, 2, . . . , N − 1;
k = 0, 1, 2, . . . , M − 1,
(32)
where X(r) and C (2−klr) take the discrete values of the Fourier transforms of x(t) and c(t) respectively. A typical time–scale distribution of the wavelet transform will require M FFTs of N samples. Usually, M is small, and therefore a complete map of the wavelet transform can be completed quickly. For the calculation of the discrete form of the FFT, aliasing-free conditions must be satisfied for a given cut-off frequency. Equation (32) can be implemented easily on a personal computer by a short program, in which a standard FFT algorithm can be called. The resulting distribution can be displayed in contour form as a two-dimensional chart with a typical size of 1024 × 10 pixels, representing the rotation angle of the gear and wavelet transforms with different scaling factors s.
5. EXAMPLES
The time domain averages described in this section were obtained from a helicopter main rotor gearbox undergoing a fatigue test in a full-scale back-to-back gearbox test facility [13]. During sustained testing at very high load, a fatigue crack began at the root of one tooth in the input spiral bevel pinion and subsequently propagated along the length of the tooth. In Figure 4, the top left curve shows the time domain average of the vibration of the gear obtained before the fatigue crack began to propagate. The gear has 22 teeth, but in the time domain average the fourth harmonic at 88 orders is dominant, either due to a
935
Figure 4. Time domain averages and residuals for helicopter gear.
high frequency response in the vibration propagation path, or else the particular gear tooth geometry. By removing the fundamental and harmonics of the tooth meshing frequency, the residual signal shown in the bottom left of Figure 4 is obtained, representing the difference between the actual vibration signal and that caused by a hypothetical average tooth [14]. At this stage, there is no indication of tooth damage. Figure 5 is the time–scale distribution of the wavelet transform for the residual signal, based upon the Gaussian-enveloped oscillation wavelet family. Absolute values of the
Figure 5. Time–scale distribution of residual for undamaged gear (1/2 octave passband).
936
. . . . c
Figure 6. Time–scale distribution of residual for damaged gear (1/2 octave passband).
wavelet transform are presented. Due to the variations between the teeth caused by normal manufacturing errors, extensive patterns appear at all scales. All sizes of variation in the signal are fully displayed. There is no dominant peak in the distribution, showing that there is no damage. At the low scale of 2, three large ‘‘hills’’ spaced at intervals of one third of a revolution are apparent. The cause of these peaks is not certain, but they could originate from a coupling external to the gearbox. Small sizes of variation can be found in the higher scale area. However, the pattern is relatively uniform and at a low level, confirming that the gear is in good condition. At low values of s, at which large sizes of variations in the signal should be apparent, the distribution is also relatively uniform through the whole rotation of the gear. The top right curve of Figure 4 shows the time domain average for the vibration of the same gear obtained shortly after the fatigue crack had started to propagate. Note that there is still no indication of damage which can be observed as the time domain average is dominated by the strong vibration at the harmonics of the tooth meshing frequency. After removal of the fundamental and harmonics, the residual signal of the bottom right of Figure 4 shows a large burst of vibration around sample 750, near 270 degrees of rotation, caused by the fatigue crack, followed by a smaller burst slightly later at sample 950, about 345 degrees. Figure 6 represents the time–scale distributions of the residual signal obtained for the wavelet transform. In the distribution, when using absolute values, clear peaks near sample 750 show details of the components. At low scales, no significant changes appear. In Figure 7 is shown the distribution when using 20 scales, which gives finer details along the scale co-ordinate, which may sometimes be necessary but at the cost of increased computing time.
937
Figure 7. Time–scale distribution of residual for damaged gear (1/8 octave passband).
6. DISCUSSION
There are strong similarities between the time–scale distributions shown in Figures 5–7 and the time–frequency distribution of the Gabor spectrogram [5, 7], but for the important difference that in the wavelet transform the resolution becomes higher as the scale s becomes larger. Low resolution is used at low scales to show large sizes of faults and high resolution is used at high scales to display small sizes of faults. The examples given demonstrate that the time–scale distribution of the wavelet transform can successfully represent, in a single display, both small and large sizes of variations in the vibration of a gear, such as the small variation from one tooth to the next due to normal manufacturing errors, or the large scale changes due to tooth damage, eccentricity or other faults. This may be contrasted with the time–frequency distribution, in which the width of the window remains constant across the distribution, so that a selected width may only be suitable for one particular type of fault under study. Some progress has already been made toward automated interpretation of time–frequency distributions [15]. The time–scale distribution may need to be inspected and interpreted by a sophisticated system since, firstly, each scale value corresponds to a bandpass filter with pre-defined boundaries and a centre frequency and, secondly, a different wavelet family gives a completely different distribution in the time–scale domain. The choice of types can be varied, depending on the objective of the signal processing. The features to be identified in the time–scale distribution are different from those in the time–frequency distribution, and new rules to suit a particular time–scale distribution need to be worked out for automatic interpretation. The basis of selecting the wavelet family to perform the transform is to find a function to compare the local components of interest. Therefore, the transform gives a high level if at a certain scale and time there is a similarity between the wavelet and the analyzed signal, otherwise it gives a low level. The change of scale can be linear or octave,
938
. . . .
and the number of scales can also be chosen according to the requirement of frequency resolution.
7. CONCLUSIONS
This paper has outlined the definition of the wavelet transform and then demonstrated how it can be applied to the analysis of the vibration signals produced by a gear in a helicopter gearbox in order to represent gear condition and detect faults. It has been shown that the time–scale distribution of the wavelet transform can be employed to analyze the local features of vibration signals. Unlike the time–frequency distribution, which incorporates a constant time and frequency resolution, the wavelet transform can display simultaneously both the large and small sizes in a signal, enabling the detection of both distributed and local faults. From the possible choices of wavelets, it is suggested that the Gaussian-enveloped oscillating wavelet is well-suited to detecting various sizes of gear faults. By adopting a series of different scales and shifting along the time axis via all locations, the wavelet can be arranged to compare all sections of the analyzed signal so that all faults, large size and small size, can be represented in a single display. In the implementation, effective wavelet bandwidths and speed of dilation have been recommended, to ensure that the number of scales just covers the frequency band of interest and that redundancy of computation is minimized.
ACKNOWLEDGMENT
The authors wish to thank the Mechanical Research Division of Westland Helicopters Limited for supplying the Wessex gearbox data.
REFERENCES 1. S. G. B 1986 Mechanical Signature Analysis. Orlando, Florida: Academic Press. See pp. 35–66. 2. P. D. MF 1987 Mechanical Systems and Signal Processing 1(1), 83–95. A revised model for the extraction of periodic waveforms by time domain averaging. 3. T. A. C. M. C and W. F. G. M¨ 1980 Philips Journal of Research 35, 372–389. The Wigner distribution—a tool for time–frequency signal analysis, part 3: relations with other time–frequency signal transformations. 4. L. C 1989 Proceedings of the IEEE 77(7), 941–981. Time–frequency distribution—a review. 5. W. J. W and P. D. McF 1993 Mechanical Systems and Signal Processing 7(3), 193–203. Early detection of gear failure by vibration analysis, I: calculation of the time–frequency distribution. 6. S. G. M 1989 IEEE Transactions on Acoustics Speech and Signal Processing 37(12), 2091–2110. Multifrequency channel decompositions of images and wavelet models. 7. D. G 1946 Journal of the Institution of Electrical Engineers, London 93, 429–457. Theory of communication. 8. D. E. N 1993 An Introduction to Random Vibration, Spectral and Wavelet Analysis. London: Longman Scientific & Technical. 9. W. J. W and P. D. MF 1993 Transactions of the American Society of Mechanical Engineers, Structural Dynamics and Vibration, 16th Annual Energy-Sources Technology Conference and Exhibition, Houston, 31 January–3 February, PD-Vol. 52, 15–20. Application of the wavelet transform to gearbox vibration analysis. 10. G. G and G. M 1984 SIAM Journal of Mathematics 15, 723–736. Decomposition of Hardy functions into square integrable wavelets of constant shape. 11. C. K. C 1992 An Introduction to Wavelets, Volume 1. Orlando, Florida: Academic Press. 12. I. D 1992 Ten Lectures on Wavelets. Pennsylvania: SIAM.
939
13. M. J. S 1981 Westland Helicopters Ltd, Mechanical Research Report MRR-20019. Wessex main rotor gearbox (MOD 1762) fatigue test—vibration monitoring—trial no. 3 (final). 14. P. D. MF 1987 Mechanical Systems and Signal Processing 1(2), 173–183. Examination of a technique for the early detection of failure in gears by signal processing of the time domain average of the meshing vibration. 15. W. J. W and P. D. MF 1983 Mechanical Systems and Signal Processing 7(3), 205–215. Early detection of gear failure by vibration analysis, II: interpretation of the time–frequency distribution using image processing techniques.
APPLICATION OF WAVELETS TO GEARBOX VIBRATION SIGNALS FOR FAULT DETECTION W. J. W Department of Mechanical and Manufacturing Engineering, De Montfort University, Leicester, England
P. D. MF Department of Engineering Science, University of Oxford, Oxford OX1 3PJ, England (Received 19 May 1994, and in final form 11 September 1995) The wavelet transform is used to represent all possible types of transients in vibration signals generated by faults in a gearbox. It is shown that the transform provides a powerful tool for condition monitoring and fault diagnosis. The vibration signal from a helicopter gearbox is used to demonstrate the application of the suggested wavelet by a simple computer algorithm. The major advantage of the wavelet transform for analyzing the signal is that it possesses multi-resolutions for localizing short-time components so that all possible types of gear faults can be displayed by a single time–scale distribution resulting from the transform. 7 1996 Academic Press Limited
1. INTRODUCTION
To determine the condition of an inaccessible gear in an operating machine, the vibration signal of the machine can be measured at some convenient location on the casing, together with a signal which is synchronous with the rotation of a shaft or gear in the machine. The time domain average of the vibration of an individual gear can then be extracted from the original vibration signal, representing in the time domain the vibration of that gear alone over one revolution [1, 2]. Although the time domain average of the vibration of a gear contains all the information describing the gear condition, clear symptoms of any damage may not be directly visible in the time domain average, particularly in the early stages. Further processing of the time domain average may be necessary to detect and diagnose the first signs of gear damage, and a wide variety of different techniques have been explored over the years. Most techniques have sought to represent the gear vibration signal in either the time domain or the frequency domain, with each approach having its advantages and disadvantages. However, a more recent trend has been toward representation in the time–frequency or time–scale domain, seeking to combine the advantages of both approaches. A time–frequency distribution describes simultaneously when a signal component occurs and how its frequency spectrum develops with time [3–5]. One important time–frequency representation, the windowed Fourier transform, gives a constant resolution in the time and frequency domains because the window width in the time domain is fixed. Reducing 927 0022–460X/96/200927 + 13 $18.00/0
7 1996 Academic Press Limited
928
. . . .
the window width will increase the time resolution, but will reduce the frequency resolution, because under the uncertainty principle, a higher time and frequency resolution cannot be achieved simultaneously [5]. This may not be a problem when only one resolution is required to display the components of interest. However, if resolutions are required to be different to suit a variety of signal components of different duration, for example in fault diagnosis where suitable resolutions may be not known in advance, this fixed-window time–frequency distribution does not meet the requirement of displaying all possible sizes of components. Thus the windowed Fourier transform is really only suited to the analysis of signals where all of the patterns appear at approximately the same size [6]. In the time–frequency representation, the fixed resolution represents only one of all the possible sizes. The features of other sizes in the signal cannot be expressed with reasonably high resolutions. For example, in the existing time–frequency representation known as the Gabor spectrogram [5, 7], the recommended width of the window for the detection of gear tooth damage is related to the meshing period of one tooth. This fixed width may enable representation of certain sizes of local damage on a tooth, but may not be suitable for other faults, such as the long duration of distributed faults, or the very short duration of early cracks. The recently developed wavelet transform uses a changing width function so that not one but a series of resolutions can be achieved in a single time–scale display. Unlike the Fourier transform, in which a signal is decomposed on to a sinusoid function basis, the wavelet transform uses a more general function as the basis. This produces, on the one hand, more comprehensive transform results but, on the other hand, a variety of possible explanations. The selection of a wavelet basis thus becomes a major issue in the application of the wavelet transform. With the rapid theoretical progress being made in the wavelet transform over the past ten years, its applications have rapidly penetrated areas of one-dimensional and two-dimensional signal analysis in many fields of science and engineering. It is suggested here that if the property of the time–scale localization of the wavelet transform is applied to the analysis of local phenomena in vibration signals, such as the short-time variations caused by faults, it may be possible to give time and scale information for all of the possible sizes of features so as to find the faults and assess their nature from a single three-dimensional representation [8, 9]. The vibration signal picked up from an operating machine has long been used to diagnose the condition of an inaccessible moving component inside the machine. Of course, this vibration is assumed to be associated with the component, and normally the component is one of the vibration generators. Therefore, any damage, such as a fatigue crack, will cause a relatively short-time transient in the vibration. In the early stages, this change is small, so that a high detection sensitivity is required. In gear fault diagnosis, although the information about the gear condition can be conveyed by the vibration signal transmitted through the gearbox casing, and the time domain average [1, 2] produces a high signal-to-noise ratio for the gear in question, clear symptoms may not be directly observed from the average and more information about its nature is very hard to determine directly from the averaged time history. Time–frequency distributions have already been applied to detect faults, and have been proved very useful [5]. However, in more general cases, many faults with different durations are possible in a mechanical component. Faults in a gear, such as wear and spalling, may occur with many different sizes. Since the wavelet transform [6] uses a series of sizes of windows to compare with all sections of the signal, it is therefore possible to display them all simultaneously.
929
2. WAVELET TRANSFORM
To overcome the limitation of the fixed resolution of the windowed Fourier transform in the time and frequency domains, Grossmann and Morlet [10] decomposed the signal x(t) into a family of functions which are the translation and dilation of a unique-valued function c(t), to give the wavelet transform which is defined by WTx (t, s) =
g
a
x(t)zsc(s(t − t)) dt,
(1)
−a
where s is a scaling factor which produces dilation, t is the time or some other spatial co-ordinate, and c(t) is called a wavelet. The corresponding wavelet family is generated by translation of the wavelet in the time domain and dilation in the scale domain, given by (zsc(s(t − t)))(t,s) $ R2 ,
(2)
where R denotes the set of real numbers. Any function in L2(R), the vector space of measurable square integrable one-dimensional functions x(t), can be characterized by its decomposition of the wavelet family of equation (2). A wavelet transform can be interpreted as a decomposition of a signal into a set of frequency channels. A process of dilation in a wavelet family is illustrated in Figure 1. As the scaling factor increases, from bottom to top, the width of the wavelet envelope decreases. In Figure 2 are shown, from left to right, the corresponding Fourier spectra for the wavelet family. As the scaling factor increases, both the centre frequency and the width of the frequency band increase. By using all of the scales, the full frequency range can be covered, from a low frequency to the required cut-off frequency.
Figure 1. Wavelets in the time domain. (a) Real part; (b) imaginary part.
. . . .
930
Figure 2. Wavelets in the frequency domain; spectra at different scales.
The reconstruction of the signal x(t) from the wavelet transform is x(t) =
1 Cc
g g a
−a
a
WTx (t, s)zsc(s(t − t)) ds dt,
(3)
0
where Cc =
g
a
0
=C( f )=2 df f
(4)
and C( f ) is the Fourier transform of c(t). In order to reconstruct x(t), Cc must satisfy Cc Q +a.
(5)
An important particular case of the discrete wavelet transform is that some wavelets c(t) exist such that the family z2 jc(2 j(t − 2−jn))( j,n) $ Z2 ,
(6) 2
where Z denotes the set of integers, is an orthonormal basis of L (R) [8, 11, 12]. The orthonormal wavelet transform is defined as WTx (n, j) =
g
a
x(t)z2 jc(2 j(t − 2−jn)) dt.
(7)
−a
For orthogonal wavelets, the reconstruction of signal x(t) is then x(t) = s s WTx (n, j)z2 jc(2 j (t − 2−jn)).
(8)
j$Z n$Z
Orthogonal wavelets, such as the well-known Daubechies wavelet series, offer faster algorithms and no redundancy in the decomposition [8, 12]. The original signal can be reconstructed by a minimum number of wavelet terms. That is, the original signal is decomposed into a compact wavelet series, which is advantageous in signal reconstruction. However, due to its special shape and limited number of scales, a single wavelet amplitude
931
map has not enough scales to describe all the details, both large and small, of the signal. Since N = 2M samples normally provide M + 1 different widths of wavelets, all scales of transients in the analyzed signal will be scaled only into the nearest size of wavelet. In gear damage detection, a small change in the vibration signal may ‘‘miss the chance’’ of matching the wavelet of the same size during the translation and dilation, as the locations of the orthogonal wavelets are limited, and there is no intermediate scale and location provided, so that no wavelet matches well the features of the signal. Another problem of using the orthogonal wavelets is the time translation variant. That is, the same transients at different time may be shown as different patterns. Therefore, using the non-orthogonal wavelet transform is preferred for the time–scale representation of short-time features, avoiding the mentioned problems. The wavelet transform means that the function x(t) is characterized by its decomposition of a wavelet family with a series of different frequency bandwidths. This can be explained by working out the bandwidth of the wavelet family. The Fourier transform of the wavelet family zsc(st) is Cs ( f ) = (1/zs)C( f/s).
(9)
Let f0 be the centre of the passband of C( f ). Then
g
a
g
a
( f − f0 )=C( f )=2 df = 0,
0
so that f0 =
>g
f =C( f )=2 df
0
a
=C( f )=2 df.
(10)
0
For the wavelet family, one finds, by substituting equation (9) into (10), that the centre frequency of the passband is
g
a
f =Cs ( f )=2 df
0
>g
a
=Cs ( f )=2 df = sf0 .
(11)
0
The root mean square bandwidth of the wavelet around f0 is given by sf , where sf2 =
g
a
( f − f0 )2 =C( f )=2 df.
(12)
0
With C( f ) replaced by Cs ( f ) and f0 by sf0 in the above equation, the root mean square bandwidth is seen to be ssf . That is, the bandwidth of the wavelet is proportional to the scaling factor s. Thus the wavelet transform produces a series of frequency channels with increasing bandwidth as the scale s changes.
3. DAMAGE DETECTION
Early damage to a gear tooth usually causes a variation in the associated vibration signal over a short time, initially less than one tooth meshing period, taking the form of modulated or unmodulated oscillation. In later stages, the duration of the abnormal variation becomes longer, lasting more than one tooth meshing period. Other distributed faults, such as eccentricity and wear, may cover the most part of the whole revolution of the gear.
. . . .
932
To suit these features, the wavelet c(t) can be chosen to be a Gaussian-enveloped oscillation, given by c(t) = c exp(−s 2t 2 − i2pf0 t),
(13)
where c, s and f0 are positive parameters, s determines the width of the wavelet and hence the width of the frequency band, and f0 is the frequency of oscillation which is the centre frequency of the band. The parameter c can usually be unity. It is to be expected that by changing the scaling factor in its corresponding wavelet family, this function will be suitable for the modelling of all sizes of oscillating transients in a gear vibration signal. The Fourier transform of c(t) is given by C( f ) = (czp/s) exp(−(p/s)2( f + f0 )2 ).
(14)
It can be seen that the frequency of oscillation f0 appears at the centre of the frequency band of the wavelet. From the definition of equation (1), the wavelet transform is WTx (t, s) =
g
a
x(t)czs exp(−s 2s 2(t − t)2 − i2pf0 s(t − t)) dt.
(15)
−a
For transients, the wavelet c(t) can also be chosen from other functions, such as a symmetrical Gaussian-enveloped oscillating function, given by c(t) = c exp(−s 2t 2 ) cos (2pf0 t)
(16)
C( f ) = (czp/s)[exp(−(p/s)2( f + f0 )2 ) + exp(−(p/s)2( f − f0 )2 )].
(17)
with a Fourier transform of
The wavelet family czs exp(−s 2s 2(t − t)2 − i2pf0 s(t − t)) is illustrated in Figure 1, which models the oscillating transients with all oscillating frequencies as the scale changes. At the same time, the width of the wavelet is reduced as the scale becomes larger. The Fourier spectra for all scales of wavelets are shown in Figure 2. When a signal is decomposed on to one of these wavelet bases, large values will be obtained if there is a transient with the same shape as the wavelet. The wavelet family of equation (13) is well suited to the detection of the transients in the signal, such as the transients caused in a vibration signal by a cracked gear tooth.
4. COMPUTER IMPLEMENTATION
For gear damage diagnosis, the function x(t) can be either the time domain average of the vibration of a gear or the corresponding analytic function, both of which are periodic with the gear rotation. The selected wavelet of equation (16) satisfies c(−t) = c (t), so that, by using the convolution theorem, equation (15) can be written as WTx (t, s) = (1/zs)F−1 (X( f ) C ( f/s)),
(18)
where F represents the Fourier transform, and C ( f ) represents F(c (t)). For the Gaussian enveloped oscillating wavelet in equation (13), the wavelet transform can be written as WTx (t, s) = (c/zss)F−1 (X( f ) exp(−(p/s)2(( f/s) − f0 )2 )).
(19)
Thus the calculation of the wavelet transform can be implemented by an inverse Fourier transform, taking advantage of the efficient fast Fourier transform (FFT) algorithm.
933
The Fourier transform of the wavelet family can be written as (c/zss) exp(−(p/ss)2( f − sf0 )2 ),
(20)
which implies that the centre frequency of the wavelet is sf0 , proportional to s. Since the half-power bandwidths in the time and frequency domains are tB = 1·18/s and fB = 0·375s respectively, then at scaling factor s, they are therefore 1·18/ss and 0·375ss respectively. The narrower the wavelet is in the time domain, the wider it is in the frequency domain. This increases the frequency overlap between the wavelets of adjacent high scales, which brings about redundancy. Consider a new discrete scaling factor s = 2kl,
k = 0, 1, 2, . . . , M − 1,
(21) M
where M is an integer. For N samples of signal, M can be defined by 2 = N. The centre frequency and half-power bandwidth of the kth wavelet become f0 (k) = 2klf0 ,
fB (k) = 0·375 × 2kls.
(22, 23)
If the cut-off frequency is fC , the analysis frequency range then will be [( f0 − 0·187s), fC ]. At the scaling factor k = M − 1, the cut-off frequency is fC = 2(M − 1)lf0 + 0·187 × 2(M − 1)ls.
(24)
Let adjacent frequency bands overlap at half-power points, so that 2klf0 + 0·187 × 2kls = 2(k + 1)lf0 − 0·187 × 2(k + 1)ls,
k = 0, 1, 2, . . . , M − 1,
which can be simplified to f0 + 0·187s = 2lf0 − 0·187 × 2ls.
(25)
In Figure 3, the top and bottom curves show upper and lower frequencies, and the middle curve is the centre frequency. From equation (25) f0 = 0·187{(2l + 1)/(2l − 1)}s.
Figure 3. Frequency bands of wavelets at different scaling factors.
(26)
. . . .
934
From equations (22), (23) and (26), the lower and upper frequencies of each band or channel are, respectively, fL (k) = f0 (k) − 12 fB (k) = 0·375
2kl s, 2 −1
fU (k) = f0 (k) + 12 fB (k) = 0·375
l
2(k + 1)l s. 2l − 1
Therefore the scaling factor s = 2kl brings about l (for example 12 or 13 ) octave frequency bands, so that fU (k) = 2lfL (k).
(27)
For a given cut-off frequency, there are four parameters, f0 , s, l and M, to be adjusted. Two independent equations, (24) and (25), can be used to determine two parameters given the values of the others. In total, there are 4 × 3/(2 × 1) = 6 combinations of situation, but probably only two of these are likely. Firstly, M and l can be set to find f0 and s. By considering equations (24) and (26), f0 and s can be found as f0 = {(2l + 1)/2Ml + 1}fC ,
s = 5·35{(2l − 1)/2Ml + 1}fc .
(28, 29)
Secondly, l and f0 can be set to find s and M. Parameters f0 and s can be found as s = 5·35
(2l − 1) f, (2l + 1) 0
M=
0
1
f 2l + 1 1·44 ln C l + 1 + 1. l f0 2
(30, 31)
Upon assuming that N samples exactly cover one period of the signal, the corresponding discrete form of equation (18) is given by the circular convolution WTx (n, k) =
N−1
1 z2
kl
s X(r)C (2−klr) exp(2pirn/N),
r=0
n = 0, 1, 2, . . . , N − 1;
k = 0, 1, 2, . . . , M − 1,
(32)
where X(r) and C (2−klr) take the discrete values of the Fourier transforms of x(t) and c(t) respectively. A typical time–scale distribution of the wavelet transform will require M FFTs of N samples. Usually, M is small, and therefore a complete map of the wavelet transform can be completed quickly. For the calculation of the discrete form of the FFT, aliasing-free conditions must be satisfied for a given cut-off frequency. Equation (32) can be implemented easily on a personal computer by a short program, in which a standard FFT algorithm can be called. The resulting distribution can be displayed in contour form as a two-dimensional chart with a typical size of 1024 × 10 pixels, representing the rotation angle of the gear and wavelet transforms with different scaling factors s.
5. EXAMPLES
The time domain averages described in this section were obtained from a helicopter main rotor gearbox undergoing a fatigue test in a full-scale back-to-back gearbox test facility [13]. During sustained testing at very high load, a fatigue crack began at the root of one tooth in the input spiral bevel pinion and subsequently propagated along the length of the tooth. In Figure 4, the top left curve shows the time domain average of the vibration of the gear obtained before the fatigue crack began to propagate. The gear has 22 teeth, but in the time domain average the fourth harmonic at 88 orders is dominant, either due to a
935
Figure 4. Time domain averages and residuals for helicopter gear.
high frequency response in the vibration propagation path, or else the particular gear tooth geometry. By removing the fundamental and harmonics of the tooth meshing frequency, the residual signal shown in the bottom left of Figure 4 is obtained, representing the difference between the actual vibration signal and that caused by a hypothetical average tooth [14]. At this stage, there is no indication of tooth damage. Figure 5 is the time–scale distribution of the wavelet transform for the residual signal, based upon the Gaussian-enveloped oscillation wavelet family. Absolute values of the
Figure 5. Time–scale distribution of residual for undamaged gear (1/2 octave passband).
936
. . . . c
Figure 6. Time–scale distribution of residual for damaged gear (1/2 octave passband).
wavelet transform are presented. Due to the variations between the teeth caused by normal manufacturing errors, extensive patterns appear at all scales. All sizes of variation in the signal are fully displayed. There is no dominant peak in the distribution, showing that there is no damage. At the low scale of 2, three large ‘‘hills’’ spaced at intervals of one third of a revolution are apparent. The cause of these peaks is not certain, but they could originate from a coupling external to the gearbox. Small sizes of variation can be found in the higher scale area. However, the pattern is relatively uniform and at a low level, confirming that the gear is in good condition. At low values of s, at which large sizes of variations in the signal should be apparent, the distribution is also relatively uniform through the whole rotation of the gear. The top right curve of Figure 4 shows the time domain average for the vibration of the same gear obtained shortly after the fatigue crack had started to propagate. Note that there is still no indication of damage which can be observed as the time domain average is dominated by the strong vibration at the harmonics of the tooth meshing frequency. After removal of the fundamental and harmonics, the residual signal of the bottom right of Figure 4 shows a large burst of vibration around sample 750, near 270 degrees of rotation, caused by the fatigue crack, followed by a smaller burst slightly later at sample 950, about 345 degrees. Figure 6 represents the time–scale distributions of the residual signal obtained for the wavelet transform. In the distribution, when using absolute values, clear peaks near sample 750 show details of the components. At low scales, no significant changes appear. In Figure 7 is shown the distribution when using 20 scales, which gives finer details along the scale co-ordinate, which may sometimes be necessary but at the cost of increased computing time.
937
Figure 7. Time–scale distribution of residual for damaged gear (1/8 octave passband).
6. DISCUSSION
There are strong similarities between the time–scale distributions shown in Figures 5–7 and the time–frequency distribution of the Gabor spectrogram [5, 7], but for the important difference that in the wavelet transform the resolution becomes higher as the scale s becomes larger. Low resolution is used at low scales to show large sizes of faults and high resolution is used at high scales to display small sizes of faults. The examples given demonstrate that the time–scale distribution of the wavelet transform can successfully represent, in a single display, both small and large sizes of variations in the vibration of a gear, such as the small variation from one tooth to the next due to normal manufacturing errors, or the large scale changes due to tooth damage, eccentricity or other faults. This may be contrasted with the time–frequency distribution, in which the width of the window remains constant across the distribution, so that a selected width may only be suitable for one particular type of fault under study. Some progress has already been made toward automated interpretation of time–frequency distributions [15]. The time–scale distribution may need to be inspected and interpreted by a sophisticated system since, firstly, each scale value corresponds to a bandpass filter with pre-defined boundaries and a centre frequency and, secondly, a different wavelet family gives a completely different distribution in the time–scale domain. The choice of types can be varied, depending on the objective of the signal processing. The features to be identified in the time–scale distribution are different from those in the time–frequency distribution, and new rules to suit a particular time–scale distribution need to be worked out for automatic interpretation. The basis of selecting the wavelet family to perform the transform is to find a function to compare the local components of interest. Therefore, the transform gives a high level if at a certain scale and time there is a similarity between the wavelet and the analyzed signal, otherwise it gives a low level. The change of scale can be linear or octave,
938
. . . .
and the number of scales can also be chosen according to the requirement of frequency resolution.
7. CONCLUSIONS
This paper has outlined the definition of the wavelet transform and then demonstrated how it can be applied to the analysis of the vibration signals produced by a gear in a helicopter gearbox in order to represent gear condition and detect faults. It has been shown that the time–scale distribution of the wavelet transform can be employed to analyze the local features of vibration signals. Unlike the time–frequency distribution, which incorporates a constant time and frequency resolution, the wavelet transform can display simultaneously both the large and small sizes in a signal, enabling the detection of both distributed and local faults. From the possible choices of wavelets, it is suggested that the Gaussian-enveloped oscillating wavelet is well-suited to detecting various sizes of gear faults. By adopting a series of different scales and shifting along the time axis via all locations, the wavelet can be arranged to compare all sections of the analyzed signal so that all faults, large size and small size, can be represented in a single display. In the implementation, effective wavelet bandwidths and speed of dilation have been recommended, to ensure that the number of scales just covers the frequency band of interest and that redundancy of computation is minimized.
ACKNOWLEDGMENT
The authors wish to thank the Mechanical Research Division of Westland Helicopters Limited for supplying the Wessex gearbox data.
REFERENCES 1. S. G. B 1986 Mechanical Signature Analysis. Orlando, Florida: Academic Press. See pp. 35–66. 2. P. D. MF 1987 Mechanical Systems and Signal Processing 1(1), 83–95. A revised model for the extraction of periodic waveforms by time domain averaging. 3. T. A. C. M. C and W. F. G. M¨ 1980 Philips Journal of Research 35, 372–389. The Wigner distribution—a tool for time–frequency signal analysis, part 3: relations with other time–frequency signal transformations. 4. L. C 1989 Proceedings of the IEEE 77(7), 941–981. Time–frequency distribution—a review. 5. W. J. W and P. D. McF 1993 Mechanical Systems and Signal Processing 7(3), 193–203. Early detection of gear failure by vibration analysis, I: calculation of the time–frequency distribution. 6. S. G. M 1989 IEEE Transactions on Acoustics Speech and Signal Processing 37(12), 2091–2110. Multifrequency channel decompositions of images and wavelet models. 7. D. G 1946 Journal of the Institution of Electrical Engineers, London 93, 429–457. Theory of communication. 8. D. E. N 1993 An Introduction to Random Vibration, Spectral and Wavelet Analysis. London: Longman Scientific & Technical. 9. W. J. W and P. D. MF 1993 Transactions of the American Society of Mechanical Engineers, Structural Dynamics and Vibration, 16th Annual Energy-Sources Technology Conference and Exhibition, Houston, 31 January–3 February, PD-Vol. 52, 15–20. Application of the wavelet transform to gearbox vibration analysis. 10. G. G and G. M 1984 SIAM Journal of Mathematics 15, 723–736. Decomposition of Hardy functions into square integrable wavelets of constant shape. 11. C. K. C 1992 An Introduction to Wavelets, Volume 1. Orlando, Florida: Academic Press. 12. I. D 1992 Ten Lectures on Wavelets. Pennsylvania: SIAM.
939
13. M. J. S 1981 Westland Helicopters Ltd, Mechanical Research Report MRR-20019. Wessex main rotor gearbox (MOD 1762) fatigue test—vibration monitoring—trial no. 3 (final). 14. P. D. MF 1987 Mechanical Systems and Signal Processing 1(2), 173–183. Examination of a technique for the early detection of failure in gears by signal processing of the time domain average of the meshing vibration. 15. W. J. W and P. D. MF 1983 Mechanical Systems and Signal Processing 7(3), 205–215. Early detection of gear failure by vibration analysis, II: interpretation of the time–frequency distribution using image processing techniques.