- Email: [email protected]
The spectral kurtosis: application to the vibratory surveillance and diagnostics of rotating machines
ARTICLE IN PRESS
Mechanical Systems and Signal Processing Mechanical Systems and Signal Processing 20 (2006) 308–331 www.elsevier.com/locate/jnlabr/ymssp
The spectral kurtosis: application to the vibratory surveillance and diagnostics of rotating machines Je´roˆme Antonia,, R.B. Randallb a
Laboratoire Roberval de Me´canique, UMR CNRS 6066, University of Technology of Compie`gne, 60205 Compiegne, France b School of Mechanical and Manufacturing Engineering, University of New South Wales, Sydney, Australia Received 23 March 2004; received in revised form 24 August 2004; accepted 1 September 2004 Available online 6 November 2004
Abstract In the previous paper, the authors have demonstrated the high potential of the spectral kurtosis (SK) to detect and characterise non-stationary signals. The present paper brings together these ideas and shows how the SK can be efficiently used in the vibration-based condition monitoring of rotating machines. First, and in contrast to classical kurtosis analysis, the SK provides a robust way of detecting incipient faults even in the presence of strong masking noise. Second, the SK offers an almost unique way of designing optimal filters for filtering out the mechanical signature of faults. The first property is of practical importance for monitoring purposes, whereas the second one proves very useful in diagnostics. Another originality of the paper is the introduction of the concept of kurtogram, from which optimal band-pass filters can be deduced, for instance as a prelude to envelope analysis. All the original findings presented in the paper are illustrated using actual industrial cases. r 2004 Elsevier Ltd. All rights reserved. Keywords: Spectral kurtosis; Fault detection; Blind detection; Kurtogram
Corresponding author. Tel.: +33 3 44 23 45 25; fax: +33 3 44 23 44 77.
E-mail address: [email protected] (J. Antoni). 0888-3270/$ - see front matter r 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.ymssp.2004.09.002
ARTICLE IN PRESS J. Antoni, R.B. Randall / Mechanical Systems and Signal Processing 20 (2006) 308–331
309
1. Introduction In the previous paper, the authors have studied the concept of spectral kurtosis (SK) and they have demonstrated its real potential for detecting and characterising transient signals buried in additive noise. The SK is a spectral statistic which helpfully supplements the classical power spectral density (PSD): it ideally takes zero values at those frequencies where stationary Gaussian noise only is present and high positive values at those frequencies where transients occur. This behaviour has been shown to be closely related to the local signal-to-noise ratio—as a function of frequency—and to the temporal characteristics of the transients. In spite of being particularly suited to many detection problems, the SK had rarely been used until now probably because it lacked a formal definition and a well-understood estimation procedure. The aim of Ref. [1] was to partly fill these gaps. In this paper, we apply the results established in Ref. [1] to the vibration monitoring of rotating machines and show how the SK can be efficiently used for their surveillance and diagnosis. The idea is that many faults in rotating machines produce series of impacts—at least in their early stage of development—which in turn generate transient vibration signals. This is typically the case for spalling faults, cracks or excessive clearances in critical mechanical components such as gears or rolling element bearings. The resulting vibration signals then often have an impulse-like repetitive nature due to the recurrence of the fault symptoms as the mechanical components rotate, the rate of occurrence of which, in combination with their temporal and spectral envelopes constitute the ‘‘mechanical signature’’ [2]. Many techniques have been devised to capture mechanical signatures by means of specific processing. However, a persistent difficulty is that available vibration signals are often severely corrupted by strong levels of background noise encompassing all other vibration sources in the machine under inspection. This problem can be formulated as that of detecting transient signals in strong additive noise, which is exactly the situation the SK is suited to. The paper is organised as follows. We first show in Section 2 how the SK can be used in its raw form as a simple tool for detecting incipient faults and we provide actual examples of application. The issue there is mainly concerned with machine surveillance. We then show in Section 3, how the SK can help in designing more sophisticated detection filters that extract the fault signal from the background noise, thus making it possible to capture its mechanical signature. The issue there is more concerned with machine diagnostics. Another originality of this paper is the introduction of the concept of ‘‘kurtogram’’, which we demonstrate is a particularly powerful way of using the SK as a prerequisite to envelope analysis. Again the techniques are illustrated with real-world examples. Although this paper is application oriented, it also presents a signal-processing methodology that is interesting in its own right. As a matter of fact, the kurtosis has often been employed in the signal-processing community to solve ‘‘blind’’ problems: blind identification and equalisation of systems by output kurtosis maximisation, blind separation of mixed signals by individual maximisation of the source kurtosis, etc. Here, the spectral version of the kurtosis is used to blindly identify a detection filter. The issue is only different from the previous ones in the sense that the required filter is not the exact inverse filter (i.e. which inverts the vibration signal to recover the impact force), but a filter which optimises the capacity of detection of the fault signal.
ARTICLE IN PRESS 310
J. Antoni, R.B. Randall / Mechanical Systems and Signal Processing 20 (2006) 308–331
2. Fault detection with the SK 2.1. Signal modelling As explained in the Introduction, the problem we are faced with is that of detecting a fault signal X(t) in measurement Y(t) in presence of some strong additive noise N(t) Y ðtÞ ¼ X ðtÞ þ NðtÞ:
(1)
It is assumed here that X(t) has a transient nature whereas N(t) is stationary. Indeed, many incipient faults in rotating machines produce a series of repetitive short transient forces (clearances, looseness, spalls, cracks), which in turn excite some structural resonances. Hence, a reasonably versatile model for X(t) is the generalised shot noise process [1]: X X k hðt sk Þ; (2) X ðtÞ ¼ k
where h(t) is the impulse response resulting from a single impact and where fX k gk2Z and fsk gk2Z are sequences of random variables which account for possibly random amplitudes and random occurrences of the impacts, respectively. It is assumed that for any index k, the difference skþ1 sk is a non-negative random variable whose mean specifies the average rate of repetition of the impacts. For instance, model (2) has been successfully used for describing incipient faults in rolling element bearings [3,4], where the stochasticity of the occurrences fsk g described the random slips of the rolling elements, and the amplitudes {Xk} the time-varying amplitude-modulation of
Fig. 1. Introductory example. (a) A typical vibration signal measured on a system with a faulty rolling element bearing (inner race fault). Note the high kurtosis value. (b) The same vibration signal in the case of a weak ball fault masked by surrounding noise: the kurtosis is almost zero.
ARTICLE IN PRESS J. Antoni, R.B. Randall / Mechanical Systems and Signal Processing 20 (2006) 308–331
311
the impacts. Fig. 1(a) illustrates such a signal captured on a faulty rolling element bearing. It is also trivial to adjust model (2) to represent incipient faults in gearboxes. 2.2. The SK of a fault signal The characterisation of incipient faults as just described is usually achieved by means of statistical indicators that are sensitive to the peakiness of the signal. For instance, the kurtosis and the crest factor are particularly well suited to this task. This is in contrast to the RMS value, for example, which can only characterise the intensity of the signal but does not provide information on the temporal variations of its envelope. Spectral analysis too may fail to detect the presence of a fault signal of the type (2) in particular when the random fluctuations {Xk} and fsk g on the impacts are strong, in which case there is little hope for localising the related harmonics in the frequency domain [3,4]. In contrast, the advantage of the kurtosis indicator is that (i) it takes high values in the presence of the fault signal X(t), whereas (ii) it is ideally zero when only background noise N(t) is present. Condition (i) is true when the impulse responses generated by the impacts are sufficiently well separated and when the signal-to-noise ratio is sufficiently high. Condition (ii) is true whenever the background noise can be assumed Gaussian as a result of the superposition of a large number of independent (stationary) sources. In practice, however, all these conditions are rarely satisfied. The most serious discrepancy concerns the ‘‘sufficiently high’’ signal-to-noise ratio. As a matter of fact, background noise often embodies strong vibrations from several competing sources (e.g. harmonics of rotating parts, random impulses from friction and contact forces, flow noise, etc.) which span a large frequency range and seriously mask the signal of interest. As a result, the kurtosis is unable to capture the peakiness of the fault signal and hardly departs from the zero value. This point is illustrated in the introductory example of Fig. 1, which compares an ideal case with high signal-to-noise ratio with a critical case with poor signal-to-noise ratio. In the first case the global kurtosis is 14.08 and it clearly indicates the presence of the fault, whereas in the second it is close to zero as for stationary Gaussian noise and therefore fails to detect the fault. In this situation, the kurtosis as a global indicator is not appropriate, and it is preferable to apply it locally in different frequency bands. This is exactly what the SK does. Indeed, paper [1] gave an interpretation of the SK as the kurtosis computed at the output of a perfect filter-bank, at each frequency f. Applying the theoretical results established therein, it can be shown that the SK KY(f) of signal (1) is: K Y ðf Þ ¼
K X ðf Þ ; ½1 þ rðf Þ2
(3)
where KX(f)40 is the SK of signal (2) and rðf Þ is the noise-to-signal ratio—i.e. rðf Þ ¼ S N ðf Þ=S X ðf Þ with SN(f) and SX(f) the power spectral densities of the noise and the fault signal, respectively. Hence, KY(f) will tend to KX(f) at those frequencies where the signal-to-noise is high, and will be close to zero where the masking effect of noise is strong. This explains the unique ability of the SK to scrutinise the whole frequency domain and to find those frequency bands where the fault signal can be best detected.
ARTICLE IN PRESS 312
J. Antoni, R.B. Randall / Mechanical Systems and Signal Processing 20 (2006) 308–331
2.3. Detection with the SK Based on the aforementioned observations, we now describe how to use the SK as a detection tool for machine surveillance. The idea is very simple: compute the SK and check in each frequency band for abnormally high values which may suggest the presence of a fault. The estimation of the SK was discussed in depth in paper [1], where we proposed a simple estimator based on the short-time Fourier transform (STFT). For this estimator to be valid (unbiased), we pointed out two important conditions the analysed signal should meet: C1: the non-stationarities in the signal have slow temporal evolutions as compared to the window length of the STFT, C2: the correlation length of the signal is shorter than the analysis window of the STFT. In other words, the analysed signal must be locally stationary—also coined oscillatory. Unfortunately, this is not the case for the fault signal described by model (2): the nonstationarities associated with the impulses are so rapid that the condition C1 cannot be met [1]. The consequence is that the estimated SK explicitly depends on the window length used in the STFT. To see this, let us impose that the analysis window w(n) of the STFT has duration Nw shorter than the mean spacing between two consecutive impulses; then the STFT-based SK of a P generalised point process PðtÞ ¼ k X k dðt sk Þ actually is: g K P ðf Þ ¼ 4w g4X 2; (4) pN w P P where p is the probability of occurrence of the impulse, g4w ¼ N w n jwðnÞj4 =j n jwðnÞj2 j2 is the time–bandwidth product of the square of the analysis window, and g4X ¼ EfjX k j4 g=EfjX k j2 g2 is the normalised fourth-order moment of the random variables {Xk} [1]. As expected, the estimated SK depends on Nw and therefore is biased. However, once the bias in Kp(f) has been recognised and quantified it is no longer a problem. Turning back to the signal Y(t)=h(t)P(t)+N(t) of model (1), its STFT-based SK can be approximated as g4w 1 K y ðf Þ g4X 2 : (5) pN w ½1 þ rðf Þ2 This is an approximation because it assumes that the convolution with h(t) translates into a product with the Fourier transform H(f) in the frequency domain, which can be quite crude when the analysis window of the STFT is short. Still, Eq. (5) has the advantage of being qualitatively useful and will enable the interpretation of the strategies we will propose hereafter. For instance, for a vibration signal Y(t) captured with sampling frequency fs and containing a repetitive fault with rate of occurrence fd, its SK would be approximately f s g4w 1 K Y ðf Þ ðkX þ 3Þ 2 : (6) f d Nw ½1 þ rðf Þ2 Hence the following remarks:
the SK decreases with the rate of repetition fd of the impulses, i.e. with the rotation speed of the machine,
ARTICLE IN PRESS J. Antoni, R.B. Randall / Mechanical Systems and Signal Processing 20 (2006) 308–331
313
Table 1 Number of rolling elements
12
Shaft rotation speed Ball diameter Pitch circle diameter Load angle from the radial Sampling frequency
6 Hz 7.12 mm 38.5 mm 01 48 000 Hz
the SK increases with the sampling rate fs of the signal, the SK increases with kX ; i.e. with the intensity of the fluctuations in the impulse amplitudes, the SK decreases with the window length Nw of the STFT. The last remark suggests that a short Nw will yield high SK values and consequently may be favoured. However, too short a Nw will also produce a SK with poor spectral resolution so that some details will be lost. Therefore, we propose that in practice several window lengths Nw be tried before selecting that which best emphasises the transient characteristics of the fault signal in the SK estimate; this should be paralleled with a similar procedure on the PSD estimate where the control of spectral resolution is made easier (Condition C2). Finally, the STFT-based SK should be compared with a statistical threshold that clearly indicates which of its values are significantly greater than zero and can be unambiguously assigned to a transient signal. Such a threshold can be designed from the results of paper [1] as follows. Let us suppose the measured signal contains only stationary Gaussian noise; then the STFT-based SK has an asymptotic normal distribution with zero mean and variance 4/K, where K is the number of averages. This is valid for any frequency f in the range |f – mod(fs/2)|4fs /Nw.1 Hence, a statistical threshold sa given a a level of significance is: 2 sa ¼ u1a pffiffiffiffi ; K
jf modðf s =2Þj4
fs Nw
(7)
with u1a the percentile of the normal distribution at 1 a: This means that all values beyond this level will have probability 1 a of not being produced by transients. 2.4. Examples of application 2.4.1. Application to bearing surveillance The first example concerns the analysis of rolling element bearing signals originated from a single-stage gearbox at the University of New South Wales, Sydney. The rolling element bearing characteristics are reported in Table 1. Different faults were introduced in the bearings, and acceleration signals were measured in each case. The time signals are not displayed here because they looked essentially stationary—similar to the signal in Fig. 1b—and therefore did not exhibit the presence of the faults. Their kurtosis 1
For frequencies f not in the range |f – mod(fs/2)| 4 fs /Nw, we propose to set the STFT-based SK to zero.
ARTICLE IN PRESS 314
J. Antoni, R.B. Randall / Mechanical Systems and Signal Processing 20 (2006) 308–331
Fig. 2. (a) PSD and (b) SK computed for different frequency resolution (N w ¼ 16; 32; 64; 128; 256). The analysed signal is from a rolling element bearing with an outer race fault.
values were actually close to zero in all cases. This is due to the very strong signals from the gears, which masked everything else. Fig. 2 displays the PSD and the SK in the case of an outer race fault. Different estimates are shown which correspond to window lengths N w ¼ 16; 32; 64; 128; 256 (spectral resolutions of 3000, 1500, 750, 375, 188 Hz per line). In all cases, the SK clearly reveals the presence of a transient signal in the frequency band from 7.5 kHz upwards (only values above the 1% significance level are shown). Maximum values are reached in the band [10–13] kHz where the signal-to-noise ratio is the greatest. Fig. 3 displays a similar analysis in the case of a ball fault. Here, the SK reveals the presence of a transient signal in the frequency band from 7 kHz upwards and maximum values are reached around 14 and 23 kHz. From these figures, it can be seen that the estimate based on N w ¼ 16 definitively introduces too strong a bias, whereas those for N w 464 rapidly decrease to zero. Therefore a reasonable window length for further analyses will be N w ¼ 32264; i.e. a spectral resolution of about 750–1500 Hz per line. Finally, Figs. 4 and 5 display the difference in dB-spectra corresponding to the analysed signals and compare them with the estimated SK’s. The difference in dB-spectra is the difference between the dB-spectrum under normal (non-faulty) condition and the dB-spectrum of the faulty signal. This is supposed to indicate the frequency band where the fault takes place and is widely used in vibration monitoring provided some historical data are available. These figures show that the SK
ARTICLE IN PRESS J. Antoni, R.B. Randall / Mechanical Systems and Signal Processing 20 (2006) 308–331
315
Fig. 3. (a) PSD and (b) SK computed for different frequency resolution (N w ¼ 16; 32; 64; 128; 256). The analysed signal is from a rolling element bearing with a ball fault.
and the difference in dB-spectra are roughly comparable in shape and they share more or less the same maxima. In many other experiments the authors have made the same observation—a theoretical justification to this observation will be given in Section 3. Consequently, the difference in dB-spectra validates the effectiveness we have claimed of the SK in detecting incipient faults. More than this, the SK now even appears to be a very valuable tool for replacing the difference in dB-spectra when no historical data are available on the system to be monitored. More will be said on this topic in Section 3. 2.4.2. Application to gear surveillance This second example illustrates the use of the SK in a surveillance programme. Acceleration signals were measured on a single-stage gearbox submitted to an accelerated fatigue test at the Centre d’Etudes Techniques des Industries Me´caniques (CETIM), Senlis, France. The sampling frequency was 80 kHz so as to capture the high-frequency resonances of the structure. Measurements were regularly taken (about every hour, except on week-ends where the system was stopped) during 10 days until the gears were completely damaged due to excessive spalling of their teeth. The experiment provided a vast amount of data which were used for testing the ability of the SK to detect incipient spall faults. Fig. 6 displays the SK computed for all measurements in chronological order (from number 1 to 160) in the frequency range from 0 to 40 kHz. The spectral resolution of the estimates is 78 Hz per line (in the case of surveillance, finding the optimal window length Nw is not so much a
ARTICLE IN PRESS 316
J. Antoni, R.B. Randall / Mechanical Systems and Signal Processing 20 (2006) 308–331
Fig. 4. (a) PSD, (b) difference in dB-spectra, and (c) SK in the case of the outer race fault. Spectral resolution is 750 Hz per line.
concern since in any case the STFT-based SK will steadily increase as the fault develops; the imperative condition is to keep Nw constant during the whole test). The evolution of the SK during the test exhibited interesting patterns:
the SK had temporary high values over the whole frequency range on days 1, 6, and 8. On day 1
at the very beginning of the test these were due to the system breaking in, while on days 6 and 8 they were believed to indicate the appearance of incipient spall defects on the gear teeth, which then rapidly get run in as a consequence of the high torque under which the system was operating (unfortunately the time interval between two consecutive measurements was too coarse for measuring how long these high values remained for), the SK clearly exhibited abnormal high values from 3 h before the end of the test. These coincided with the development of excessive spalling, yet very strong indication of the fault was already noticeable on day 9. the SK was more marked above 10 kHz where the sources of vibration related to the gear dynamics had no more energy, thus providing a better signal-to-noise ratio.
ARTICLE IN PRESS J. Antoni, R.B. Randall / Mechanical Systems and Signal Processing 20 (2006) 308–331
317
Fig. 5. (a) PSD, (b) difference in dB-spectra, and (c) SK in the case of the ball fault. Spectral resolution is 750 Hz per line.
Fig. 6. SK of measurements on a gearbox submitted to an accelerated fatigue test. Abnormally high values were detected on days 2, 6, 8 and 10. The test was stopped on day 12 due to excessive spall damage.
Other techniques, not shown here, were tested simultaneously with the SK. In particular, the PSD was also computed over all measurements in the same frequency range as the SK. However, it had a fairly stable behaviour during the test and could only indicate the appearance of excessive
ARTICLE IN PRESS 318
J. Antoni, R.B. Randall / Mechanical Systems and Signal Processing 20 (2006) 308–331
spalling 1 h before the test was stopped. This basic experiment had for purpose to demonstrate the high potential of the SK for surveillance; yet more expansive research is needed to find the best the best way of applying the SK to gears (influence of operating conditions, best time interval between measurements, etc.).
3. Diagnostics with the SK Although the direct analysis of a vibration signal with the SK is very valuable to detect the presence of incipient faults, it does not provide any information on the nature of the fault itself. The objective of this section is to design some ad hoc ‘‘detection filters’’ that can extract the fault signal from its noisy environment, and therefore make it possible to inspect its mechanical signature. The main originality of our approach is to use the SK to design these filters blindly without any prior information on the system. 3.1. Detection filters Many detection filters exist in the literature, most of them relying on optimising some given criterion on the filtered signal. Wide-band filters are generally required to recover a filtered signal as close as possible to the fault signal, but they are difficult to design. Besides, band-pass filters are easier to design and are well suited to recovering a filtered signal with a high signal-to-noise ratio—hence with good detectability—but this is at the expense of signal distortion. Therefore, there is a compromise to be made which depends on the chosen criterion. Typical detection filters used in vibration analysis are the Wiener filter, the matched filter, and the dB-spectrum difference. They all correspond to different criteria, which we now review. The main point is that all of them can be blindly identified from the SK. 3.1.1. The Wiener filter The Wiener filter is the optimal filter that recovers a signal X(t) buried in additive noise N(t). It is, among all possible linear filters, that filter w(t) which minimises the mean square error between X(t) and w(t)Y(t) where Y(t)=X(t)+N(t). Its principle is illustrated in Fig. 7. The Wiener filter W(f) has solution W ðf Þ ¼
1 ; 1 þ rðf Þ
(8)
where rðf Þ is a noise-to-signal ratio (ratio of the PSD SN(f) of the noise to the PSD SX(f) of the signal). However, in practice the exact solution (8) cannot be used unless rðf Þ is known a priori. This is possible when some historical data are available, but it is unfortunately not the general case. Yet, one very useful property of the SK is that it gives access to the estimation of the Wiener filter without a priori knowledge of the noise-to-signal ratio. This is based on recognising that Eq. (6) is simply proportional to the square of the Wiener filter: K Y ðf Þ ’ k W ðf Þ2 ;
jf modðf s =2Þj4
fs : Nw
(9)
ARTICLE IN PRESS J. Antoni, R.B. Randall / Mechanical Systems and Signal Processing 20 (2006) 308–331
319
Fig. 7. Principle of the Wiener filter.
Consequently, the optimal estimation of the fault signal X(t) can be recovered by merely filtering Y(t) with the square root of the SK. Intuitively, this means that frequency components where the SK is high will be enhanced in the signal whereas the filter will attenuate frequency components where the SK is near zero. While the unknown scaling factor k in Eq. (9) provides no trouble for recovering a scaled version of the fault signal X(t), it is problematic in recovering the additive noise N(t)—if ever required. This is so since the ideal filter designed for this purpose is pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi f (10) 1 W ðf Þ ¼ 1 K Y ðf Þ=k; jf modðf s =2Þj4 s ; Nw which strongly depends on k. One solution to estimate the unknown value of k is to assume that the noise has a kurtosis close to zero, so that k can be found as that value which minimises the magnitude of the kurtosis at the output of filter (10). This must be done under the constraint that pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 0p1 K Y ðf Þ=kp1 (11) in order to find a valid solution. We have experienced that this algorithm can be simply and efficiently implemented with Newton’s method of successive approximations. 3.1.2. The matched filter The previous subsection showed how to blindly identify the Wiener filter that best recovers the fault signal from background noise. In some situations one may not be concerned with recovering the exact shape of the transient signal, but would rather prefer a filtered version of it with maximum signal-to-noise ratio. The solution to this problem is obtained by applying a ‘‘matched filter’’ m(t) that maximises the ratio of the power of X(t)m(t) to the power of N(t)m(t). In the random case, the matched filter is given by the eigenvector associated with the largest eigenvalue of an autocorrelation matrix, which has for elements the inverse Fourier transform of: Mðf Þ ¼
1 : rðf Þ
(12)
ARTICLE IN PRESS 320
J. Antoni, R.B. Randall / Mechanical Systems and Signal Processing 20 (2006) 308–331
This is closely related to the Wiener filter since, from definition (9), Mðf Þ ¼
W ðf Þ : 1 W ðf Þ
This means that the matched filter can be blindly identified from the SK as pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi K Y ðf Þ=k f pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; jf modðf s =2Þj4 s ; Mðf Þ ’ Nw 1 K Y ðf Þ=k
(13)
(14)
where, again, the unknown parameter k can be determined such that it maximises the kurtosis at the output of the matched filter under constraint (11). One particularly relevant solution to (14) is when k is set such that pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi (15) 1 max K Y ðf Þ=k ¼ 51 and therefore Mðf Þ ! dðf f 0 Þ; !0
(16)
where f0 is that frequency where the maximum of KY(f) is achieved. In other words, the matched filter tends to an ideal band-pass filter (of infinite resolution) tuned on the frequency where the signal to-noise ratio is the highest.
3.1.3. The difference in dB-spectra The difference in dB-spectra is an empirical technique used in vibration monitoring; it was illustrated in the introductory example of Section 2.4. It consists in comparing the dB-spectra before and after the fault has appeared. Their difference then indicates the frequency bands where the relative change of energy induced by the fault is the greatest, thus suggesting where to bandpass the signal. The technique is usually used prior to envelope analysis. With the previous notations, the dB-spectrum difference has the expression: DdB ðf Þ ¼ 10 log10 ½S X ðf Þ þ S N ðf Þ 10 log10 ½SN ðf Þ 1 ¼ 10 log10 1 þ : rðf Þ
ð17Þ
Hence, just for the Wiener and for the matched filters, the dB-spectrum difference takes values near zero where the noise-to-signal ratio is strong, and high positive values where it is weak. Again, it can be blindly identified from the SK, since h pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffii f (18) DdB ðf Þ ’ 10 log10 1 K Y ðf Þ=k ; jf modðf s =2Þj4 s : Nw Incidentally, this proves that the difference in dB-spectra is a monotonic function of the SK and justifies the similarity we have observed between the difference in dB-spectra and the SK in the introductory example of Section 2.4.
ARTICLE IN PRESS J. Antoni, R.B. Randall / Mechanical Systems and Signal Processing 20 (2006) 308–331
321
3.1.4. Estimation issues In practice, the blind identification of any of the above detection filters requires computing the square root of the STFT-based SK. It is obvious that only the positive values of the SK should be considered in doing so, for instance those values above the statistical threshold derived in Eq. (7). It should also be mentioned that these blindly identified filters may depart somehow from their theoretical counterparts as obtained from the exact knowledge of the noise-to-signal ratio. There are at least three reasons for this: (1) the unknown scaling factor k is estimated by kurtosis optimisation, which leads to a good though not necessarily exact solution, (2) formulae (10), (14) and (18) are valid if model (5) is; non-Gaussian and/or non-stationary noise, other impulsive forces (clearance, backlash, friction, etc) are all potential sources of bias, (3) the shape of the STFT-based SK depends on the window length Nw, and so will the shapes of the estimates (10), (14) and (18). 3.2. Application to bearing diagnostics In order to illustrate the aforementioned ideas, the blindly identified Wiener filters corresponding to the signals of Example 2.4.1 are displayed in Fig. 8. They are overlaid with their theoretical counterparts given by Eq. (8), in which the noise-to-signal ratio could be estimated by taking advantage of previously measured signals under non-faulty conditions. It is seen that the blindly identified Wiener filters roughly follow the shape of the theoretical filters and evidence obvious high-pass structures. This is because the rolling element signals were strongly masked by the gearmesh harmonics that dominated the frequency range below 10 kHz. The matched filters were also computed, although not shown here; they were found to closely resemble band-pass filters centered on those frequencies where the SK is maximum. Fig. 9 displays the first signal filtered by the Wiener filter. It exhibits a series of sharp impacts whose rate of repetition coincides with the ball pass frequency on the outer race fault (BPFO). Further processing with envelope spectrum analysis has also revealed slight modulation by the cage rotation speed; that is the typical signature of a localised fault on the outer race. Fig. 10 displays the second filtered signal. As before, it exhibits a series of sharp impacts whose rate of repetition coincides with twice the ball spin frequency (BSF) and which are slightly modulated by the cage rotation speed; that is the typical signature of a ball fault. Once again, it must be emphasised that the diagnostic information, which we have been able to extract from the proposed detection filters, was by no means available from classical spectral analysis, even by using a high-frequency resolution. As previously mentioned, the reason is that bearing signals experience random fluctuations which rapidly destroy the harmonic structure of the signal when approaching high frequencies (e.g. above the hundredth harmonic of the fault [5]). 3.3. The concept of kurtogram Subsection 3.1.4 has pointed to some limitations in blindly identifying detection filters from the SK. The idea here is to propose a more robust detection filter g(t) with an imposed band-pass
ARTICLE IN PRESS 322
J. Antoni, R.B. Randall / Mechanical Systems and Signal Processing 20 (2006) 308–331
Fig. 8. Theoretical (dotted line) and estimated (full line) Wiener filters in the case of (a) the outer race fault and (b) the ball fault.
ARTICLE IN PRESS J. Antoni, R.B. Randall / Mechanical Systems and Signal Processing 20 (2006) 308–331
323
Fig. 9. (a) Raw signal and extracted fault signal with (b) the Wiener filter and (c) the matched filter in the case of the outer race fault. The extracted signal evidences a rate of impact of 29.3 Hz, which coincides with the expected frequency of an outer race fault.
Fig. 10. (a) Raw signal and extracted fault signal with (b) the Wiener filter and (c) the matched filter in the case of the ball fault. The extracted signal evidences a rate of impact of 51.6 Hz, which coincides with the expected frequency of a ball fault.
ARTICLE IN PRESS 324
J. Antoni, R.B. Randall / Mechanical Systems and Signal Processing 20 (2006) 308–331
structure with only two parameters to be blindly identified. Specifically, the objective is to find (i) the central frequency fc and (ii) the bandwidth Bf of the filter which maximises some criterion on the filtered signal. The rationale beyond this approach is the following:
a
band-pass filter has good chance of selecting the frequency band where the signal-to-noise ratio is maximum. This is reminiscent to the aforementioned matched filter. a band-pass filter has better chance to select only one impulsive source (the strongest one) in the case where several such sources are present in the signal. a band-pass filter is a necessary pre-processing stage before carrying out envelope demodulation of the signal as usually done in vibration analysis. in addition, forcing a pre-defined structure on the filter makes it independent of the assumed validity of model (5); its scope of validity is therefore wider than for the detection filters discussed in Section 2.3.
Once again, it should be emphasised that the fact that a band-pass filter modifies the waveform of the signal is no problem, in the sense that it is not the recovery of the fault signal itself that is important but that of its statistical structure—i.e. its mechanical signature. This is similar to the concept of the matched filter, as opposed to the Wiener filter. However, the matched filter was designed to maximise the signal-to-noise ratio, which did not necessarily guarantee that the impulse-like nature of the fault was optimally recovered. In order to achieve this aim, we suggest designing the band-pass filter g(t) so as to maximise the kurtosis of the envelope of the filtered signal. It is very interesting that this problem is strictly equivalent to finding the frequency f and the window length Nw which maximises the STFT-based SK over all possible choices. The justification of this statement is immediate from the filter-bank interpretation of the STFT, where f is the central frequency of a narrow-band filter of bandwidth proportional to 1/Nw [1]. In the following we shall coin the name ‘‘kurtogram’’ for the map formed by the STFT-based SK as a function of f and Nw. Hence, the optimal central frequency fc and bandwidth Bf of the band-pass filter are found as those values which jointly maximise the kurtogram. Intuitively, the mechanism of the proposed optimal band-pass filter is the following: (i) the optimal central frequency fc tunes the filter where the SK is the maximum, i.e. where the signal-to-noise ratio is the highest according to model (5), (ii) the optimal bandwidth Bf achieves the best compromise between too wide a filter G(f) which would alter the signal-to-noise ratio, and too narrow a filter G(f)—i.e., with a very long impulse response g(t)—which would alter the impulse-like nature of the filtered signal.
3.4. Example of application 3.4.1. First example: synthetic signal This first example illustrates the computation and interpretation of the kurtogram on a synthesised fault signal. The signal was modelled according to model (1)–(2) as a regular train of Dirac pulses exciting a single degree-of-freedom system with natural frequency f n ¼ 0:29 (normalised frequency) and damping ratio z ¼ 5%: This impulse-type signal was then buried in additive white noise such as to yield a signal-to-noise ratio of –10 dB. The kurtogram was
ARTICLE IN PRESS J. Antoni, R.B. Randall / Mechanical Systems and Signal Processing 20 (2006) 308–331
325
Fig. 11. Kurtogram of a synthesised signal (SK isolines as functions of f and Nw).
computed for values of Nw from 2 to 200 samples with steps of 1 sample (although such a broad range of values with such fine a resolution is not necessary in practice). It is displayed in Fig. 11 as a contour plot showing the isolines of constant kurtosis, together with 4 horizontal lines indicating the SK’s (slices of the kurtogram) one would have obtained by setting the window length to N w ¼ 16; 32, 64 and 128. The kurtogram shows a unique maximum for f ¼ 0:29 which matches the system natural frequency and for N w ¼ 34 which, for a Hanning window such as used in our STFT-based SK, corresponds to a filter bandwidth Bf 2=N w ¼ 0:059: It makes sense that the optimal frequency f* matches the system natural frequency fn since it is there that the highest signal-to-noise is achieved. It is also very interesting that the optimal bandwidth Bf* almost matches the system bandwidth 1=2z ¼ 0:056: This can actually be understood in the light of the matched filter theory, which established that in the present case (deterministic signal in white noise), the filter maximising the output signal-to-noise ratio has similar frequency gain to the system itself. This confirms our assertion that the maximum of the kurtogram yields the parameters for optimal detection. Another experiment has been conducted to confirm this fact. The previous synthetic signal was generated several times with different realisations of additive noise and then band-pass filtered for all possible values (f, Bf). The total number of times the squared envelope of the band-pass filtered signal crossed a predefined threshold is reported in Fig. 12 in the form of a probability-ofdetection surface with respect to f and Nw. Note that the maximum probability of detection is achieved for f ¼ 0:29 and N w ¼ 28 ðBf 0:071Þ; which are optimal values quite close to those previously obtained from the kurtogram. 3.4.2. Second example: application to bearing diagnostics The results presented here are a continuation of the analyses carried out in Section 2.4.1. Figs. 13 and 14 display the kurtogram computed on the rolling element bearing signals
ARTICLE IN PRESS 326
J. Antoni, R.B. Randall / Mechanical Systems and Signal Processing 20 (2006) 308–331
Fig. 12. Probability of detection surface (squared envelope detection) as a function of f and Nw.
Fig. 13. Kurtogram of a rolling element bearing signal with an outer race fault. The global maximum is achieved for f ¼ 12:5 kHz and N w ¼ 44:
corresponding to an outer race and a ball fault, respectively—see Table 1. The maximum of the kurtogram is achieved for f ¼ 12:5 kHz and N w ¼ 44 in the first case, and for f ¼ 14:0 kHz and N w ¼ 23 in the second case. Note that these results are in very good accordance with those that could be inferred from Figs. 2 and 3, where only slices of the kurtogram are shown for N w ¼ 16; 32, 64 and 128. This shows that if the computing time is a concern, then a rough estimate of the
ARTICLE IN PRESS J. Antoni, R.B. Randall / Mechanical Systems and Signal Processing 20 (2006) 308–331
327
Fig. 14. Kurtogram of a rolling element bearing with a ball fault. The global maximum is achieved for f ¼ 14:0 kHz and N w ¼ 23:
kurtogram can be obtained from a rapid evaluation of the SK’s on a logarithmic grid of Nw’s. Finally, the band-pass filters2 designed from the selected optimal parameters are displayed in Fig. 15(a) and (b) and overlaid with the SK of the analysed signals. The resulting filtered signals are displayed in Fig. 16(a) and (b) where it is noticeable that in both cases the final kurtosis of the filtered signals have significantly increased as compared to the previous results displayed in Fig. 9 and 10. Consequently, it is concluded that the mechanical signatures have been better extracted. 3.4.3. Third example: a tricky case Finally, this last example illustrates a case where the proposed band-pass filter could extract the diagnostic information where the other detection filters failed. The analysed signal is from the same test-bench as in Example 2.4.1 (Table 1, input shaft speed=10 Hz) but had been recorded with different operating conditions and a different bearing. The Wiener filter and the matched filter were blindly identified as explained in section 3.1, and revealed wide-band structures. The corresponding filtered signals did not show any spectacular increase in their kurtosis, yet they both evidenced a marked periodic modulation by the input shaft rotation. For instance, Fig. 17(b) displays the filtered signal from the matched filter. We believe this strong modulation was due to misalignment forces on the gears. Unfortunately, it completely masked a bearing fault, which was found after applying the optimal band-pass filter as obtained from the maximum of the kurtogram. The optimally band-pass filtered signal is displayed in Fig. 17(c) and it now clearly evidences the presence of an impulse-like fault whose rate of repetition was checked to coincide with an outer race fault. 2
They were obtained by using least squares FIR structures of 50 coefficients with central frequency and bandwidth as specified by f* and N nw ; respectively. An equivalent but theoretically more exact strategy is to use as filters the functions wðtÞcosð2pf tÞ; where w(t) is the N nw –long analysing window used in the STFT—based SK.
ARTICLE IN PRESS 328
J. Antoni, R.B. Randall / Mechanical Systems and Signal Processing 20 (2006) 308–331
Fig. 15. Optimal band-pass filters (thick line) as obtained from the maximum of the kurtogram, overlaid on the STFTbased SK’s estimates (dotted lines) with similar spectral resolution 1=N w : (a) Signal with an outer race fault. (b) Signal with a ball fault.
Fig. 16. Filtered signals with the optimal band-pass filters as obtained from the maximum of the kurtogram.
ARTICLE IN PRESS J. Antoni, R.B. Randall / Mechanical Systems and Signal Processing 20 (2006) 308–331
329
Fig. 17. (a) Raw signal and extracted fault signal with (b) the matched filter and (c) the optimal band-pass filter in the case of the outer race fault. Only the last extracted signal evidences a rate of impact of 49.1 Hz, which coincides with the expected frequency of an outer race fault (BPFO). Note the difference in amplitude of the different signals.
Fig. 18. Magnitude of the squared envelope spectrum. (a) Filtered signal from the matched filter. (b) Filtered signal from the optimal band-pass filter as obtained from the maximum of the kurtogram. Note the difference in amplitude of the two envelope spectra.
ARTICLE IN PRESS 330
J. Antoni, R.B. Randall / Mechanical Systems and Signal Processing 20 (2006) 308–331
For a better comparison between these results, Fig. 18 shows the spectra of the squared envelopes of the filtered signals. While the envelope from the matched filter is totally dominated by the input shaft speed harmonics, that from the optimal band-pass filter only shows the harmonics of the bearing fault.
4. Conclusion This paper has demonstrated how the SK—whose definition and properties have been introduced in a previous paper [1]—can be elegantly and very efficiently applied to the vibratory monitoring of rotating machines. The key idea is to use the high sensitivity of the SK for detecting and characterising incipient faults that produce impulse-like signals. The methodology we propose proceeds in two steps. (1) First, we suggest using the SK as a detection tool that precisely points out in which frequency band(s) the fault shows the best contrast from background noise. (2) Second, we suggest using the SK as a basis for designing ad hoc detection filters that can extract the mechanical signature of the fault. Of particular interest is the concept of kurtogram which displays the SK as a function of frequency and of spectral resolution. The maximum of the kurtogram then provides the optimal parameters from which a band-pass filter can be designed, for instance as a prelude to envelope analysis. We have shown through actual examples that step (1) finds useful applications in surveillance protocols, whereas step (2) is of evident utility for diagnostic purposes. This is particularly true when no historical data from the analysed machine are available, since then the SK offers the rather unique opportunity to find out which frequency band should be processed. Here probably lies the most important advantage of our approach as opposed to other classical techniques. This point has been extensively illustrated and validated in the case of early diagnostics of rolling element-bearing signals. Indeed, from the authors’ knowledge, this is the first instance in which a sound criterion is proposed that can automatically optimise the parameters used in envelope analysis. This result is of great importance for practical applications, and refutes the common belief that envelope analysis must rely on excessive user’s participation and expertise [6–8].
References [1] J. Antoni, R.B. Randall, The spectral kurtosis: a useful tool for characterising nonstationary signals, Mechanical Systems and Signal Processing, 2004, this issue. [2] S. Braun, Mechanical Signature Analysis: Theory and Applications, Academic Press, New York, 1986. [3] R.B. Randall, J. Antoni, S. Chobsaard, The relationship between spectral correlation and envelope analysis in the diagnostics of bearing faults and other cyclostationary machine, Mechanical Systems and Signal Processing 5 (5) (2001) 945–962. [4] J. Antoni, R.B. Randall, Differential diagnosis of gear and bearing faults, ASME Journal of Vibration and Acoustics 124 (2) (2002) 165–171. [5] D. Ho, R.B. Randall, Optimisation of bearing diagnostic techniques using simulated and actual bearing fault signals, Mechanical Systems and Signal Processing 14 (5) (2000) 763–788.
ARTICLE IN PRESS J. Antoni, R.B. Randall / Mechanical Systems and Signal Processing 20 (2006) 308–331
331
[6] P. Tse, Y.H. Peng, R. Yam, Wavelet analysis and envelope detection for rolling element bearing fault diagnosis— their effectiveness and flexibility, Transactions of the ASME, Journal of Vibration and Acoustics 123 (3) (2001) 303–310. [7] J. Altmann, J. Mathew, Multiple band-pass autoregressive demodulation for rolling-element bearing fault diagnosis, Mechanical Systems and Signal Processing 15 (5) (2001) 963–977. [8] N.G. Nikolaou, I.A. Antoniadis, Demodulation of vibration signals generated by defects in rolling element bearings using complex shifted morlet wavelets, Mechanical Systems and Signal Processing 16 (4) (2002) 677–694.
Mechanical Systems and Signal Processing Mechanical Systems and Signal Processing 20 (2006) 308–331 www.elsevier.com/locate/jnlabr/ymssp
The spectral kurtosis: application to the vibratory surveillance and diagnostics of rotating machines Je´roˆme Antonia,, R.B. Randallb a
Laboratoire Roberval de Me´canique, UMR CNRS 6066, University of Technology of Compie`gne, 60205 Compiegne, France b School of Mechanical and Manufacturing Engineering, University of New South Wales, Sydney, Australia Received 23 March 2004; received in revised form 24 August 2004; accepted 1 September 2004 Available online 6 November 2004
Abstract In the previous paper, the authors have demonstrated the high potential of the spectral kurtosis (SK) to detect and characterise non-stationary signals. The present paper brings together these ideas and shows how the SK can be efficiently used in the vibration-based condition monitoring of rotating machines. First, and in contrast to classical kurtosis analysis, the SK provides a robust way of detecting incipient faults even in the presence of strong masking noise. Second, the SK offers an almost unique way of designing optimal filters for filtering out the mechanical signature of faults. The first property is of practical importance for monitoring purposes, whereas the second one proves very useful in diagnostics. Another originality of the paper is the introduction of the concept of kurtogram, from which optimal band-pass filters can be deduced, for instance as a prelude to envelope analysis. All the original findings presented in the paper are illustrated using actual industrial cases. r 2004 Elsevier Ltd. All rights reserved. Keywords: Spectral kurtosis; Fault detection; Blind detection; Kurtogram
Corresponding author. Tel.: +33 3 44 23 45 25; fax: +33 3 44 23 44 77.
E-mail address: [email protected] (J. Antoni). 0888-3270/$ - see front matter r 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.ymssp.2004.09.002
ARTICLE IN PRESS J. Antoni, R.B. Randall / Mechanical Systems and Signal Processing 20 (2006) 308–331
309
1. Introduction In the previous paper, the authors have studied the concept of spectral kurtosis (SK) and they have demonstrated its real potential for detecting and characterising transient signals buried in additive noise. The SK is a spectral statistic which helpfully supplements the classical power spectral density (PSD): it ideally takes zero values at those frequencies where stationary Gaussian noise only is present and high positive values at those frequencies where transients occur. This behaviour has been shown to be closely related to the local signal-to-noise ratio—as a function of frequency—and to the temporal characteristics of the transients. In spite of being particularly suited to many detection problems, the SK had rarely been used until now probably because it lacked a formal definition and a well-understood estimation procedure. The aim of Ref. [1] was to partly fill these gaps. In this paper, we apply the results established in Ref. [1] to the vibration monitoring of rotating machines and show how the SK can be efficiently used for their surveillance and diagnosis. The idea is that many faults in rotating machines produce series of impacts—at least in their early stage of development—which in turn generate transient vibration signals. This is typically the case for spalling faults, cracks or excessive clearances in critical mechanical components such as gears or rolling element bearings. The resulting vibration signals then often have an impulse-like repetitive nature due to the recurrence of the fault symptoms as the mechanical components rotate, the rate of occurrence of which, in combination with their temporal and spectral envelopes constitute the ‘‘mechanical signature’’ [2]. Many techniques have been devised to capture mechanical signatures by means of specific processing. However, a persistent difficulty is that available vibration signals are often severely corrupted by strong levels of background noise encompassing all other vibration sources in the machine under inspection. This problem can be formulated as that of detecting transient signals in strong additive noise, which is exactly the situation the SK is suited to. The paper is organised as follows. We first show in Section 2 how the SK can be used in its raw form as a simple tool for detecting incipient faults and we provide actual examples of application. The issue there is mainly concerned with machine surveillance. We then show in Section 3, how the SK can help in designing more sophisticated detection filters that extract the fault signal from the background noise, thus making it possible to capture its mechanical signature. The issue there is more concerned with machine diagnostics. Another originality of this paper is the introduction of the concept of ‘‘kurtogram’’, which we demonstrate is a particularly powerful way of using the SK as a prerequisite to envelope analysis. Again the techniques are illustrated with real-world examples. Although this paper is application oriented, it also presents a signal-processing methodology that is interesting in its own right. As a matter of fact, the kurtosis has often been employed in the signal-processing community to solve ‘‘blind’’ problems: blind identification and equalisation of systems by output kurtosis maximisation, blind separation of mixed signals by individual maximisation of the source kurtosis, etc. Here, the spectral version of the kurtosis is used to blindly identify a detection filter. The issue is only different from the previous ones in the sense that the required filter is not the exact inverse filter (i.e. which inverts the vibration signal to recover the impact force), but a filter which optimises the capacity of detection of the fault signal.
ARTICLE IN PRESS 310
J. Antoni, R.B. Randall / Mechanical Systems and Signal Processing 20 (2006) 308–331
2. Fault detection with the SK 2.1. Signal modelling As explained in the Introduction, the problem we are faced with is that of detecting a fault signal X(t) in measurement Y(t) in presence of some strong additive noise N(t) Y ðtÞ ¼ X ðtÞ þ NðtÞ:
(1)
It is assumed here that X(t) has a transient nature whereas N(t) is stationary. Indeed, many incipient faults in rotating machines produce a series of repetitive short transient forces (clearances, looseness, spalls, cracks), which in turn excite some structural resonances. Hence, a reasonably versatile model for X(t) is the generalised shot noise process [1]: X X k hðt sk Þ; (2) X ðtÞ ¼ k
where h(t) is the impulse response resulting from a single impact and where fX k gk2Z and fsk gk2Z are sequences of random variables which account for possibly random amplitudes and random occurrences of the impacts, respectively. It is assumed that for any index k, the difference skþ1 sk is a non-negative random variable whose mean specifies the average rate of repetition of the impacts. For instance, model (2) has been successfully used for describing incipient faults in rolling element bearings [3,4], where the stochasticity of the occurrences fsk g described the random slips of the rolling elements, and the amplitudes {Xk} the time-varying amplitude-modulation of
Fig. 1. Introductory example. (a) A typical vibration signal measured on a system with a faulty rolling element bearing (inner race fault). Note the high kurtosis value. (b) The same vibration signal in the case of a weak ball fault masked by surrounding noise: the kurtosis is almost zero.
ARTICLE IN PRESS J. Antoni, R.B. Randall / Mechanical Systems and Signal Processing 20 (2006) 308–331
311
the impacts. Fig. 1(a) illustrates such a signal captured on a faulty rolling element bearing. It is also trivial to adjust model (2) to represent incipient faults in gearboxes. 2.2. The SK of a fault signal The characterisation of incipient faults as just described is usually achieved by means of statistical indicators that are sensitive to the peakiness of the signal. For instance, the kurtosis and the crest factor are particularly well suited to this task. This is in contrast to the RMS value, for example, which can only characterise the intensity of the signal but does not provide information on the temporal variations of its envelope. Spectral analysis too may fail to detect the presence of a fault signal of the type (2) in particular when the random fluctuations {Xk} and fsk g on the impacts are strong, in which case there is little hope for localising the related harmonics in the frequency domain [3,4]. In contrast, the advantage of the kurtosis indicator is that (i) it takes high values in the presence of the fault signal X(t), whereas (ii) it is ideally zero when only background noise N(t) is present. Condition (i) is true when the impulse responses generated by the impacts are sufficiently well separated and when the signal-to-noise ratio is sufficiently high. Condition (ii) is true whenever the background noise can be assumed Gaussian as a result of the superposition of a large number of independent (stationary) sources. In practice, however, all these conditions are rarely satisfied. The most serious discrepancy concerns the ‘‘sufficiently high’’ signal-to-noise ratio. As a matter of fact, background noise often embodies strong vibrations from several competing sources (e.g. harmonics of rotating parts, random impulses from friction and contact forces, flow noise, etc.) which span a large frequency range and seriously mask the signal of interest. As a result, the kurtosis is unable to capture the peakiness of the fault signal and hardly departs from the zero value. This point is illustrated in the introductory example of Fig. 1, which compares an ideal case with high signal-to-noise ratio with a critical case with poor signal-to-noise ratio. In the first case the global kurtosis is 14.08 and it clearly indicates the presence of the fault, whereas in the second it is close to zero as for stationary Gaussian noise and therefore fails to detect the fault. In this situation, the kurtosis as a global indicator is not appropriate, and it is preferable to apply it locally in different frequency bands. This is exactly what the SK does. Indeed, paper [1] gave an interpretation of the SK as the kurtosis computed at the output of a perfect filter-bank, at each frequency f. Applying the theoretical results established therein, it can be shown that the SK KY(f) of signal (1) is: K Y ðf Þ ¼
K X ðf Þ ; ½1 þ rðf Þ2
(3)
where KX(f)40 is the SK of signal (2) and rðf Þ is the noise-to-signal ratio—i.e. rðf Þ ¼ S N ðf Þ=S X ðf Þ with SN(f) and SX(f) the power spectral densities of the noise and the fault signal, respectively. Hence, KY(f) will tend to KX(f) at those frequencies where the signal-to-noise is high, and will be close to zero where the masking effect of noise is strong. This explains the unique ability of the SK to scrutinise the whole frequency domain and to find those frequency bands where the fault signal can be best detected.
ARTICLE IN PRESS 312
J. Antoni, R.B. Randall / Mechanical Systems and Signal Processing 20 (2006) 308–331
2.3. Detection with the SK Based on the aforementioned observations, we now describe how to use the SK as a detection tool for machine surveillance. The idea is very simple: compute the SK and check in each frequency band for abnormally high values which may suggest the presence of a fault. The estimation of the SK was discussed in depth in paper [1], where we proposed a simple estimator based on the short-time Fourier transform (STFT). For this estimator to be valid (unbiased), we pointed out two important conditions the analysed signal should meet: C1: the non-stationarities in the signal have slow temporal evolutions as compared to the window length of the STFT, C2: the correlation length of the signal is shorter than the analysis window of the STFT. In other words, the analysed signal must be locally stationary—also coined oscillatory. Unfortunately, this is not the case for the fault signal described by model (2): the nonstationarities associated with the impulses are so rapid that the condition C1 cannot be met [1]. The consequence is that the estimated SK explicitly depends on the window length used in the STFT. To see this, let us impose that the analysis window w(n) of the STFT has duration Nw shorter than the mean spacing between two consecutive impulses; then the STFT-based SK of a P generalised point process PðtÞ ¼ k X k dðt sk Þ actually is: g K P ðf Þ ¼ 4w g4X 2; (4) pN w P P where p is the probability of occurrence of the impulse, g4w ¼ N w n jwðnÞj4 =j n jwðnÞj2 j2 is the time–bandwidth product of the square of the analysis window, and g4X ¼ EfjX k j4 g=EfjX k j2 g2 is the normalised fourth-order moment of the random variables {Xk} [1]. As expected, the estimated SK depends on Nw and therefore is biased. However, once the bias in Kp(f) has been recognised and quantified it is no longer a problem. Turning back to the signal Y(t)=h(t)P(t)+N(t) of model (1), its STFT-based SK can be approximated as g4w 1 K y ðf Þ g4X 2 : (5) pN w ½1 þ rðf Þ2 This is an approximation because it assumes that the convolution with h(t) translates into a product with the Fourier transform H(f) in the frequency domain, which can be quite crude when the analysis window of the STFT is short. Still, Eq. (5) has the advantage of being qualitatively useful and will enable the interpretation of the strategies we will propose hereafter. For instance, for a vibration signal Y(t) captured with sampling frequency fs and containing a repetitive fault with rate of occurrence fd, its SK would be approximately f s g4w 1 K Y ðf Þ ðkX þ 3Þ 2 : (6) f d Nw ½1 þ rðf Þ2 Hence the following remarks:
the SK decreases with the rate of repetition fd of the impulses, i.e. with the rotation speed of the machine,
ARTICLE IN PRESS J. Antoni, R.B. Randall / Mechanical Systems and Signal Processing 20 (2006) 308–331
313
Table 1 Number of rolling elements
12
Shaft rotation speed Ball diameter Pitch circle diameter Load angle from the radial Sampling frequency
6 Hz 7.12 mm 38.5 mm 01 48 000 Hz
the SK increases with the sampling rate fs of the signal, the SK increases with kX ; i.e. with the intensity of the fluctuations in the impulse amplitudes, the SK decreases with the window length Nw of the STFT. The last remark suggests that a short Nw will yield high SK values and consequently may be favoured. However, too short a Nw will also produce a SK with poor spectral resolution so that some details will be lost. Therefore, we propose that in practice several window lengths Nw be tried before selecting that which best emphasises the transient characteristics of the fault signal in the SK estimate; this should be paralleled with a similar procedure on the PSD estimate where the control of spectral resolution is made easier (Condition C2). Finally, the STFT-based SK should be compared with a statistical threshold that clearly indicates which of its values are significantly greater than zero and can be unambiguously assigned to a transient signal. Such a threshold can be designed from the results of paper [1] as follows. Let us suppose the measured signal contains only stationary Gaussian noise; then the STFT-based SK has an asymptotic normal distribution with zero mean and variance 4/K, where K is the number of averages. This is valid for any frequency f in the range |f – mod(fs/2)|4fs /Nw.1 Hence, a statistical threshold sa given a a level of significance is: 2 sa ¼ u1a pffiffiffiffi ; K
jf modðf s =2Þj4
fs Nw
(7)
with u1a the percentile of the normal distribution at 1 a: This means that all values beyond this level will have probability 1 a of not being produced by transients. 2.4. Examples of application 2.4.1. Application to bearing surveillance The first example concerns the analysis of rolling element bearing signals originated from a single-stage gearbox at the University of New South Wales, Sydney. The rolling element bearing characteristics are reported in Table 1. Different faults were introduced in the bearings, and acceleration signals were measured in each case. The time signals are not displayed here because they looked essentially stationary—similar to the signal in Fig. 1b—and therefore did not exhibit the presence of the faults. Their kurtosis 1
For frequencies f not in the range |f – mod(fs/2)| 4 fs /Nw, we propose to set the STFT-based SK to zero.
ARTICLE IN PRESS 314
J. Antoni, R.B. Randall / Mechanical Systems and Signal Processing 20 (2006) 308–331
Fig. 2. (a) PSD and (b) SK computed for different frequency resolution (N w ¼ 16; 32; 64; 128; 256). The analysed signal is from a rolling element bearing with an outer race fault.
values were actually close to zero in all cases. This is due to the very strong signals from the gears, which masked everything else. Fig. 2 displays the PSD and the SK in the case of an outer race fault. Different estimates are shown which correspond to window lengths N w ¼ 16; 32; 64; 128; 256 (spectral resolutions of 3000, 1500, 750, 375, 188 Hz per line). In all cases, the SK clearly reveals the presence of a transient signal in the frequency band from 7.5 kHz upwards (only values above the 1% significance level are shown). Maximum values are reached in the band [10–13] kHz where the signal-to-noise ratio is the greatest. Fig. 3 displays a similar analysis in the case of a ball fault. Here, the SK reveals the presence of a transient signal in the frequency band from 7 kHz upwards and maximum values are reached around 14 and 23 kHz. From these figures, it can be seen that the estimate based on N w ¼ 16 definitively introduces too strong a bias, whereas those for N w 464 rapidly decrease to zero. Therefore a reasonable window length for further analyses will be N w ¼ 32264; i.e. a spectral resolution of about 750–1500 Hz per line. Finally, Figs. 4 and 5 display the difference in dB-spectra corresponding to the analysed signals and compare them with the estimated SK’s. The difference in dB-spectra is the difference between the dB-spectrum under normal (non-faulty) condition and the dB-spectrum of the faulty signal. This is supposed to indicate the frequency band where the fault takes place and is widely used in vibration monitoring provided some historical data are available. These figures show that the SK
ARTICLE IN PRESS J. Antoni, R.B. Randall / Mechanical Systems and Signal Processing 20 (2006) 308–331
315
Fig. 3. (a) PSD and (b) SK computed for different frequency resolution (N w ¼ 16; 32; 64; 128; 256). The analysed signal is from a rolling element bearing with a ball fault.
and the difference in dB-spectra are roughly comparable in shape and they share more or less the same maxima. In many other experiments the authors have made the same observation—a theoretical justification to this observation will be given in Section 3. Consequently, the difference in dB-spectra validates the effectiveness we have claimed of the SK in detecting incipient faults. More than this, the SK now even appears to be a very valuable tool for replacing the difference in dB-spectra when no historical data are available on the system to be monitored. More will be said on this topic in Section 3. 2.4.2. Application to gear surveillance This second example illustrates the use of the SK in a surveillance programme. Acceleration signals were measured on a single-stage gearbox submitted to an accelerated fatigue test at the Centre d’Etudes Techniques des Industries Me´caniques (CETIM), Senlis, France. The sampling frequency was 80 kHz so as to capture the high-frequency resonances of the structure. Measurements were regularly taken (about every hour, except on week-ends where the system was stopped) during 10 days until the gears were completely damaged due to excessive spalling of their teeth. The experiment provided a vast amount of data which were used for testing the ability of the SK to detect incipient spall faults. Fig. 6 displays the SK computed for all measurements in chronological order (from number 1 to 160) in the frequency range from 0 to 40 kHz. The spectral resolution of the estimates is 78 Hz per line (in the case of surveillance, finding the optimal window length Nw is not so much a
ARTICLE IN PRESS 316
J. Antoni, R.B. Randall / Mechanical Systems and Signal Processing 20 (2006) 308–331
Fig. 4. (a) PSD, (b) difference in dB-spectra, and (c) SK in the case of the outer race fault. Spectral resolution is 750 Hz per line.
concern since in any case the STFT-based SK will steadily increase as the fault develops; the imperative condition is to keep Nw constant during the whole test). The evolution of the SK during the test exhibited interesting patterns:
the SK had temporary high values over the whole frequency range on days 1, 6, and 8. On day 1
at the very beginning of the test these were due to the system breaking in, while on days 6 and 8 they were believed to indicate the appearance of incipient spall defects on the gear teeth, which then rapidly get run in as a consequence of the high torque under which the system was operating (unfortunately the time interval between two consecutive measurements was too coarse for measuring how long these high values remained for), the SK clearly exhibited abnormal high values from 3 h before the end of the test. These coincided with the development of excessive spalling, yet very strong indication of the fault was already noticeable on day 9. the SK was more marked above 10 kHz where the sources of vibration related to the gear dynamics had no more energy, thus providing a better signal-to-noise ratio.
ARTICLE IN PRESS J. Antoni, R.B. Randall / Mechanical Systems and Signal Processing 20 (2006) 308–331
317
Fig. 5. (a) PSD, (b) difference in dB-spectra, and (c) SK in the case of the ball fault. Spectral resolution is 750 Hz per line.
Fig. 6. SK of measurements on a gearbox submitted to an accelerated fatigue test. Abnormally high values were detected on days 2, 6, 8 and 10. The test was stopped on day 12 due to excessive spall damage.
Other techniques, not shown here, were tested simultaneously with the SK. In particular, the PSD was also computed over all measurements in the same frequency range as the SK. However, it had a fairly stable behaviour during the test and could only indicate the appearance of excessive
ARTICLE IN PRESS 318
J. Antoni, R.B. Randall / Mechanical Systems and Signal Processing 20 (2006) 308–331
spalling 1 h before the test was stopped. This basic experiment had for purpose to demonstrate the high potential of the SK for surveillance; yet more expansive research is needed to find the best the best way of applying the SK to gears (influence of operating conditions, best time interval between measurements, etc.).
3. Diagnostics with the SK Although the direct analysis of a vibration signal with the SK is very valuable to detect the presence of incipient faults, it does not provide any information on the nature of the fault itself. The objective of this section is to design some ad hoc ‘‘detection filters’’ that can extract the fault signal from its noisy environment, and therefore make it possible to inspect its mechanical signature. The main originality of our approach is to use the SK to design these filters blindly without any prior information on the system. 3.1. Detection filters Many detection filters exist in the literature, most of them relying on optimising some given criterion on the filtered signal. Wide-band filters are generally required to recover a filtered signal as close as possible to the fault signal, but they are difficult to design. Besides, band-pass filters are easier to design and are well suited to recovering a filtered signal with a high signal-to-noise ratio—hence with good detectability—but this is at the expense of signal distortion. Therefore, there is a compromise to be made which depends on the chosen criterion. Typical detection filters used in vibration analysis are the Wiener filter, the matched filter, and the dB-spectrum difference. They all correspond to different criteria, which we now review. The main point is that all of them can be blindly identified from the SK. 3.1.1. The Wiener filter The Wiener filter is the optimal filter that recovers a signal X(t) buried in additive noise N(t). It is, among all possible linear filters, that filter w(t) which minimises the mean square error between X(t) and w(t)Y(t) where Y(t)=X(t)+N(t). Its principle is illustrated in Fig. 7. The Wiener filter W(f) has solution W ðf Þ ¼
1 ; 1 þ rðf Þ
(8)
where rðf Þ is a noise-to-signal ratio (ratio of the PSD SN(f) of the noise to the PSD SX(f) of the signal). However, in practice the exact solution (8) cannot be used unless rðf Þ is known a priori. This is possible when some historical data are available, but it is unfortunately not the general case. Yet, one very useful property of the SK is that it gives access to the estimation of the Wiener filter without a priori knowledge of the noise-to-signal ratio. This is based on recognising that Eq. (6) is simply proportional to the square of the Wiener filter: K Y ðf Þ ’ k W ðf Þ2 ;
jf modðf s =2Þj4
fs : Nw
(9)
ARTICLE IN PRESS J. Antoni, R.B. Randall / Mechanical Systems and Signal Processing 20 (2006) 308–331
319
Fig. 7. Principle of the Wiener filter.
Consequently, the optimal estimation of the fault signal X(t) can be recovered by merely filtering Y(t) with the square root of the SK. Intuitively, this means that frequency components where the SK is high will be enhanced in the signal whereas the filter will attenuate frequency components where the SK is near zero. While the unknown scaling factor k in Eq. (9) provides no trouble for recovering a scaled version of the fault signal X(t), it is problematic in recovering the additive noise N(t)—if ever required. This is so since the ideal filter designed for this purpose is pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi f (10) 1 W ðf Þ ¼ 1 K Y ðf Þ=k; jf modðf s =2Þj4 s ; Nw which strongly depends on k. One solution to estimate the unknown value of k is to assume that the noise has a kurtosis close to zero, so that k can be found as that value which minimises the magnitude of the kurtosis at the output of filter (10). This must be done under the constraint that pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 0p1 K Y ðf Þ=kp1 (11) in order to find a valid solution. We have experienced that this algorithm can be simply and efficiently implemented with Newton’s method of successive approximations. 3.1.2. The matched filter The previous subsection showed how to blindly identify the Wiener filter that best recovers the fault signal from background noise. In some situations one may not be concerned with recovering the exact shape of the transient signal, but would rather prefer a filtered version of it with maximum signal-to-noise ratio. The solution to this problem is obtained by applying a ‘‘matched filter’’ m(t) that maximises the ratio of the power of X(t)m(t) to the power of N(t)m(t). In the random case, the matched filter is given by the eigenvector associated with the largest eigenvalue of an autocorrelation matrix, which has for elements the inverse Fourier transform of: Mðf Þ ¼
1 : rðf Þ
(12)
ARTICLE IN PRESS 320
J. Antoni, R.B. Randall / Mechanical Systems and Signal Processing 20 (2006) 308–331
This is closely related to the Wiener filter since, from definition (9), Mðf Þ ¼
W ðf Þ : 1 W ðf Þ
This means that the matched filter can be blindly identified from the SK as pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi K Y ðf Þ=k f pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; jf modðf s =2Þj4 s ; Mðf Þ ’ Nw 1 K Y ðf Þ=k
(13)
(14)
where, again, the unknown parameter k can be determined such that it maximises the kurtosis at the output of the matched filter under constraint (11). One particularly relevant solution to (14) is when k is set such that pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi (15) 1 max K Y ðf Þ=k ¼ 51 and therefore Mðf Þ ! dðf f 0 Þ; !0
(16)
where f0 is that frequency where the maximum of KY(f) is achieved. In other words, the matched filter tends to an ideal band-pass filter (of infinite resolution) tuned on the frequency where the signal to-noise ratio is the highest.
3.1.3. The difference in dB-spectra The difference in dB-spectra is an empirical technique used in vibration monitoring; it was illustrated in the introductory example of Section 2.4. It consists in comparing the dB-spectra before and after the fault has appeared. Their difference then indicates the frequency bands where the relative change of energy induced by the fault is the greatest, thus suggesting where to bandpass the signal. The technique is usually used prior to envelope analysis. With the previous notations, the dB-spectrum difference has the expression: DdB ðf Þ ¼ 10 log10 ½S X ðf Þ þ S N ðf Þ 10 log10 ½SN ðf Þ 1 ¼ 10 log10 1 þ : rðf Þ
ð17Þ
Hence, just for the Wiener and for the matched filters, the dB-spectrum difference takes values near zero where the noise-to-signal ratio is strong, and high positive values where it is weak. Again, it can be blindly identified from the SK, since h pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffii f (18) DdB ðf Þ ’ 10 log10 1 K Y ðf Þ=k ; jf modðf s =2Þj4 s : Nw Incidentally, this proves that the difference in dB-spectra is a monotonic function of the SK and justifies the similarity we have observed between the difference in dB-spectra and the SK in the introductory example of Section 2.4.
ARTICLE IN PRESS J. Antoni, R.B. Randall / Mechanical Systems and Signal Processing 20 (2006) 308–331
321
3.1.4. Estimation issues In practice, the blind identification of any of the above detection filters requires computing the square root of the STFT-based SK. It is obvious that only the positive values of the SK should be considered in doing so, for instance those values above the statistical threshold derived in Eq. (7). It should also be mentioned that these blindly identified filters may depart somehow from their theoretical counterparts as obtained from the exact knowledge of the noise-to-signal ratio. There are at least three reasons for this: (1) the unknown scaling factor k is estimated by kurtosis optimisation, which leads to a good though not necessarily exact solution, (2) formulae (10), (14) and (18) are valid if model (5) is; non-Gaussian and/or non-stationary noise, other impulsive forces (clearance, backlash, friction, etc) are all potential sources of bias, (3) the shape of the STFT-based SK depends on the window length Nw, and so will the shapes of the estimates (10), (14) and (18). 3.2. Application to bearing diagnostics In order to illustrate the aforementioned ideas, the blindly identified Wiener filters corresponding to the signals of Example 2.4.1 are displayed in Fig. 8. They are overlaid with their theoretical counterparts given by Eq. (8), in which the noise-to-signal ratio could be estimated by taking advantage of previously measured signals under non-faulty conditions. It is seen that the blindly identified Wiener filters roughly follow the shape of the theoretical filters and evidence obvious high-pass structures. This is because the rolling element signals were strongly masked by the gearmesh harmonics that dominated the frequency range below 10 kHz. The matched filters were also computed, although not shown here; they were found to closely resemble band-pass filters centered on those frequencies where the SK is maximum. Fig. 9 displays the first signal filtered by the Wiener filter. It exhibits a series of sharp impacts whose rate of repetition coincides with the ball pass frequency on the outer race fault (BPFO). Further processing with envelope spectrum analysis has also revealed slight modulation by the cage rotation speed; that is the typical signature of a localised fault on the outer race. Fig. 10 displays the second filtered signal. As before, it exhibits a series of sharp impacts whose rate of repetition coincides with twice the ball spin frequency (BSF) and which are slightly modulated by the cage rotation speed; that is the typical signature of a ball fault. Once again, it must be emphasised that the diagnostic information, which we have been able to extract from the proposed detection filters, was by no means available from classical spectral analysis, even by using a high-frequency resolution. As previously mentioned, the reason is that bearing signals experience random fluctuations which rapidly destroy the harmonic structure of the signal when approaching high frequencies (e.g. above the hundredth harmonic of the fault [5]). 3.3. The concept of kurtogram Subsection 3.1.4 has pointed to some limitations in blindly identifying detection filters from the SK. The idea here is to propose a more robust detection filter g(t) with an imposed band-pass
ARTICLE IN PRESS 322
J. Antoni, R.B. Randall / Mechanical Systems and Signal Processing 20 (2006) 308–331
Fig. 8. Theoretical (dotted line) and estimated (full line) Wiener filters in the case of (a) the outer race fault and (b) the ball fault.
ARTICLE IN PRESS J. Antoni, R.B. Randall / Mechanical Systems and Signal Processing 20 (2006) 308–331
323
Fig. 9. (a) Raw signal and extracted fault signal with (b) the Wiener filter and (c) the matched filter in the case of the outer race fault. The extracted signal evidences a rate of impact of 29.3 Hz, which coincides with the expected frequency of an outer race fault.
Fig. 10. (a) Raw signal and extracted fault signal with (b) the Wiener filter and (c) the matched filter in the case of the ball fault. The extracted signal evidences a rate of impact of 51.6 Hz, which coincides with the expected frequency of a ball fault.
ARTICLE IN PRESS 324
J. Antoni, R.B. Randall / Mechanical Systems and Signal Processing 20 (2006) 308–331
structure with only two parameters to be blindly identified. Specifically, the objective is to find (i) the central frequency fc and (ii) the bandwidth Bf of the filter which maximises some criterion on the filtered signal. The rationale beyond this approach is the following:
a
band-pass filter has good chance of selecting the frequency band where the signal-to-noise ratio is maximum. This is reminiscent to the aforementioned matched filter. a band-pass filter has better chance to select only one impulsive source (the strongest one) in the case where several such sources are present in the signal. a band-pass filter is a necessary pre-processing stage before carrying out envelope demodulation of the signal as usually done in vibration analysis. in addition, forcing a pre-defined structure on the filter makes it independent of the assumed validity of model (5); its scope of validity is therefore wider than for the detection filters discussed in Section 2.3.
Once again, it should be emphasised that the fact that a band-pass filter modifies the waveform of the signal is no problem, in the sense that it is not the recovery of the fault signal itself that is important but that of its statistical structure—i.e. its mechanical signature. This is similar to the concept of the matched filter, as opposed to the Wiener filter. However, the matched filter was designed to maximise the signal-to-noise ratio, which did not necessarily guarantee that the impulse-like nature of the fault was optimally recovered. In order to achieve this aim, we suggest designing the band-pass filter g(t) so as to maximise the kurtosis of the envelope of the filtered signal. It is very interesting that this problem is strictly equivalent to finding the frequency f and the window length Nw which maximises the STFT-based SK over all possible choices. The justification of this statement is immediate from the filter-bank interpretation of the STFT, where f is the central frequency of a narrow-band filter of bandwidth proportional to 1/Nw [1]. In the following we shall coin the name ‘‘kurtogram’’ for the map formed by the STFT-based SK as a function of f and Nw. Hence, the optimal central frequency fc and bandwidth Bf of the band-pass filter are found as those values which jointly maximise the kurtogram. Intuitively, the mechanism of the proposed optimal band-pass filter is the following: (i) the optimal central frequency fc tunes the filter where the SK is the maximum, i.e. where the signal-to-noise ratio is the highest according to model (5), (ii) the optimal bandwidth Bf achieves the best compromise between too wide a filter G(f) which would alter the signal-to-noise ratio, and too narrow a filter G(f)—i.e., with a very long impulse response g(t)—which would alter the impulse-like nature of the filtered signal.
3.4. Example of application 3.4.1. First example: synthetic signal This first example illustrates the computation and interpretation of the kurtogram on a synthesised fault signal. The signal was modelled according to model (1)–(2) as a regular train of Dirac pulses exciting a single degree-of-freedom system with natural frequency f n ¼ 0:29 (normalised frequency) and damping ratio z ¼ 5%: This impulse-type signal was then buried in additive white noise such as to yield a signal-to-noise ratio of –10 dB. The kurtogram was
ARTICLE IN PRESS J. Antoni, R.B. Randall / Mechanical Systems and Signal Processing 20 (2006) 308–331
325
Fig. 11. Kurtogram of a synthesised signal (SK isolines as functions of f and Nw).
computed for values of Nw from 2 to 200 samples with steps of 1 sample (although such a broad range of values with such fine a resolution is not necessary in practice). It is displayed in Fig. 11 as a contour plot showing the isolines of constant kurtosis, together with 4 horizontal lines indicating the SK’s (slices of the kurtogram) one would have obtained by setting the window length to N w ¼ 16; 32, 64 and 128. The kurtogram shows a unique maximum for f ¼ 0:29 which matches the system natural frequency and for N w ¼ 34 which, for a Hanning window such as used in our STFT-based SK, corresponds to a filter bandwidth Bf 2=N w ¼ 0:059: It makes sense that the optimal frequency f* matches the system natural frequency fn since it is there that the highest signal-to-noise is achieved. It is also very interesting that the optimal bandwidth Bf* almost matches the system bandwidth 1=2z ¼ 0:056: This can actually be understood in the light of the matched filter theory, which established that in the present case (deterministic signal in white noise), the filter maximising the output signal-to-noise ratio has similar frequency gain to the system itself. This confirms our assertion that the maximum of the kurtogram yields the parameters for optimal detection. Another experiment has been conducted to confirm this fact. The previous synthetic signal was generated several times with different realisations of additive noise and then band-pass filtered for all possible values (f, Bf). The total number of times the squared envelope of the band-pass filtered signal crossed a predefined threshold is reported in Fig. 12 in the form of a probability-ofdetection surface with respect to f and Nw. Note that the maximum probability of detection is achieved for f ¼ 0:29 and N w ¼ 28 ðBf 0:071Þ; which are optimal values quite close to those previously obtained from the kurtogram. 3.4.2. Second example: application to bearing diagnostics The results presented here are a continuation of the analyses carried out in Section 2.4.1. Figs. 13 and 14 display the kurtogram computed on the rolling element bearing signals
ARTICLE IN PRESS 326
J. Antoni, R.B. Randall / Mechanical Systems and Signal Processing 20 (2006) 308–331
Fig. 12. Probability of detection surface (squared envelope detection) as a function of f and Nw.
Fig. 13. Kurtogram of a rolling element bearing signal with an outer race fault. The global maximum is achieved for f ¼ 12:5 kHz and N w ¼ 44:
corresponding to an outer race and a ball fault, respectively—see Table 1. The maximum of the kurtogram is achieved for f ¼ 12:5 kHz and N w ¼ 44 in the first case, and for f ¼ 14:0 kHz and N w ¼ 23 in the second case. Note that these results are in very good accordance with those that could be inferred from Figs. 2 and 3, where only slices of the kurtogram are shown for N w ¼ 16; 32, 64 and 128. This shows that if the computing time is a concern, then a rough estimate of the
ARTICLE IN PRESS J. Antoni, R.B. Randall / Mechanical Systems and Signal Processing 20 (2006) 308–331
327
Fig. 14. Kurtogram of a rolling element bearing with a ball fault. The global maximum is achieved for f ¼ 14:0 kHz and N w ¼ 23:
kurtogram can be obtained from a rapid evaluation of the SK’s on a logarithmic grid of Nw’s. Finally, the band-pass filters2 designed from the selected optimal parameters are displayed in Fig. 15(a) and (b) and overlaid with the SK of the analysed signals. The resulting filtered signals are displayed in Fig. 16(a) and (b) where it is noticeable that in both cases the final kurtosis of the filtered signals have significantly increased as compared to the previous results displayed in Fig. 9 and 10. Consequently, it is concluded that the mechanical signatures have been better extracted. 3.4.3. Third example: a tricky case Finally, this last example illustrates a case where the proposed band-pass filter could extract the diagnostic information where the other detection filters failed. The analysed signal is from the same test-bench as in Example 2.4.1 (Table 1, input shaft speed=10 Hz) but had been recorded with different operating conditions and a different bearing. The Wiener filter and the matched filter were blindly identified as explained in section 3.1, and revealed wide-band structures. The corresponding filtered signals did not show any spectacular increase in their kurtosis, yet they both evidenced a marked periodic modulation by the input shaft rotation. For instance, Fig. 17(b) displays the filtered signal from the matched filter. We believe this strong modulation was due to misalignment forces on the gears. Unfortunately, it completely masked a bearing fault, which was found after applying the optimal band-pass filter as obtained from the maximum of the kurtogram. The optimally band-pass filtered signal is displayed in Fig. 17(c) and it now clearly evidences the presence of an impulse-like fault whose rate of repetition was checked to coincide with an outer race fault. 2
They were obtained by using least squares FIR structures of 50 coefficients with central frequency and bandwidth as specified by f* and N nw ; respectively. An equivalent but theoretically more exact strategy is to use as filters the functions wðtÞcosð2pf tÞ; where w(t) is the N nw –long analysing window used in the STFT—based SK.
ARTICLE IN PRESS 328
J. Antoni, R.B. Randall / Mechanical Systems and Signal Processing 20 (2006) 308–331
Fig. 15. Optimal band-pass filters (thick line) as obtained from the maximum of the kurtogram, overlaid on the STFTbased SK’s estimates (dotted lines) with similar spectral resolution 1=N w : (a) Signal with an outer race fault. (b) Signal with a ball fault.
Fig. 16. Filtered signals with the optimal band-pass filters as obtained from the maximum of the kurtogram.
ARTICLE IN PRESS J. Antoni, R.B. Randall / Mechanical Systems and Signal Processing 20 (2006) 308–331
329
Fig. 17. (a) Raw signal and extracted fault signal with (b) the matched filter and (c) the optimal band-pass filter in the case of the outer race fault. Only the last extracted signal evidences a rate of impact of 49.1 Hz, which coincides with the expected frequency of an outer race fault (BPFO). Note the difference in amplitude of the different signals.
Fig. 18. Magnitude of the squared envelope spectrum. (a) Filtered signal from the matched filter. (b) Filtered signal from the optimal band-pass filter as obtained from the maximum of the kurtogram. Note the difference in amplitude of the two envelope spectra.
ARTICLE IN PRESS 330
J. Antoni, R.B. Randall / Mechanical Systems and Signal Processing 20 (2006) 308–331
For a better comparison between these results, Fig. 18 shows the spectra of the squared envelopes of the filtered signals. While the envelope from the matched filter is totally dominated by the input shaft speed harmonics, that from the optimal band-pass filter only shows the harmonics of the bearing fault.
4. Conclusion This paper has demonstrated how the SK—whose definition and properties have been introduced in a previous paper [1]—can be elegantly and very efficiently applied to the vibratory monitoring of rotating machines. The key idea is to use the high sensitivity of the SK for detecting and characterising incipient faults that produce impulse-like signals. The methodology we propose proceeds in two steps. (1) First, we suggest using the SK as a detection tool that precisely points out in which frequency band(s) the fault shows the best contrast from background noise. (2) Second, we suggest using the SK as a basis for designing ad hoc detection filters that can extract the mechanical signature of the fault. Of particular interest is the concept of kurtogram which displays the SK as a function of frequency and of spectral resolution. The maximum of the kurtogram then provides the optimal parameters from which a band-pass filter can be designed, for instance as a prelude to envelope analysis. We have shown through actual examples that step (1) finds useful applications in surveillance protocols, whereas step (2) is of evident utility for diagnostic purposes. This is particularly true when no historical data from the analysed machine are available, since then the SK offers the rather unique opportunity to find out which frequency band should be processed. Here probably lies the most important advantage of our approach as opposed to other classical techniques. This point has been extensively illustrated and validated in the case of early diagnostics of rolling element-bearing signals. Indeed, from the authors’ knowledge, this is the first instance in which a sound criterion is proposed that can automatically optimise the parameters used in envelope analysis. This result is of great importance for practical applications, and refutes the common belief that envelope analysis must rely on excessive user’s participation and expertise [6–8].
References [1] J. Antoni, R.B. Randall, The spectral kurtosis: a useful tool for characterising nonstationary signals, Mechanical Systems and Signal Processing, 2004, this issue. [2] S. Braun, Mechanical Signature Analysis: Theory and Applications, Academic Press, New York, 1986. [3] R.B. Randall, J. Antoni, S. Chobsaard, The relationship between spectral correlation and envelope analysis in the diagnostics of bearing faults and other cyclostationary machine, Mechanical Systems and Signal Processing 5 (5) (2001) 945–962. [4] J. Antoni, R.B. Randall, Differential diagnosis of gear and bearing faults, ASME Journal of Vibration and Acoustics 124 (2) (2002) 165–171. [5] D. Ho, R.B. Randall, Optimisation of bearing diagnostic techniques using simulated and actual bearing fault signals, Mechanical Systems and Signal Processing 14 (5) (2000) 763–788.
ARTICLE IN PRESS J. Antoni, R.B. Randall / Mechanical Systems and Signal Processing 20 (2006) 308–331
331
[6] P. Tse, Y.H. Peng, R. Yam, Wavelet analysis and envelope detection for rolling element bearing fault diagnosis— their effectiveness and flexibility, Transactions of the ASME, Journal of Vibration and Acoustics 123 (3) (2001) 303–310. [7] J. Altmann, J. Mathew, Multiple band-pass autoregressive demodulation for rolling-element bearing fault diagnosis, Mechanical Systems and Signal Processing 15 (5) (2001) 963–977. [8] N.G. Nikolaou, I.A. Antoniadis, Demodulation of vibration signals generated by defects in rolling element bearings using complex shifted morlet wavelets, Mechanical Systems and Signal Processing 16 (4) (2002) 677–694.