- Email: [email protected]
SINGULARITY ANALYSIS USING CONTINUOUS WAVELET TRANSFORM FOR BEARING FAULT DIAGNOSIS
Mechanical Systems and Signal Processing (2002) 16(6), 1025–1041 doi:10.1006/mssp.2002.1474, available online at http://www.idealibrary.com on
SINGULARITYANALYSIS USING CONTINUOUS WAVELET TRANSFORM FOR BEARING FAULT DIAGNOSIS Q. Sun Department of Mechanical and Manufacturing Engineering, University of Calgary, 2500 University Dr., N.W., Calgary, Alberta, Canada T2N 1N4. E-mail: [email protected]
and Y. Tang Faculty of Mechanical Engineering, Beijing Science and Technology University, University Road 23, Beijing 100083, People’s Republic of China. E-mail: [email protected] (Received 11 November 2001, accepted 1 February 2002) In this paper, wavelet transform is applied to detect abrupt changes in the vibration signals obtained from operating bearings being monitored. In particular, singularity analysis across all scales of the continuous wavelet transform is performed to identify the location (in time) of defect-induced bursts in the vibration signals. Through modifying the intensity of the wavelet transform modulus maxima, defect-related vibration signature is highlighted and can be easily associated with the bearing defect characteristic frequencies for diagnosis. Due to the fact that vibration characteristics of faulty bearings are complex and defect-related vibration signature is normally buried in the wideband noise and high frequency structural resonance, simple signal processing cannot be used to detect bearing fault. We show, through experimental results, that the proposed method has the ability to discriminate noise from the signal significantly and is robust to bearing operating conditions, such as load and speed, and severity of the bearing damage. These properties are desirable for automatic detection of machine faults. # 2002 Published by Elsevier Science Ltd.
1. INTRODUCTION
Bearing defect diagnosis is a typical problem in machine and process fault detection and diagnosis [1]. Once defects occur on the working surface of bearing elements, impact is generated when mating elements encounter the defects. Sharp transient response accompanied by damped oscillation can often be observed in the measured vibration signals. The periodic transient behaviour in the vibration signal contains rich information about bearing health. On the other hand, however, defect-induced impulses are often buried in the wideband noise accompanied by high frequency structural resonance. This imposes a great deal of difficulty in identifying defect-induced impulses from the measured vibration signals. Depending on many other factors, such as machinery operating conditions and severity of bearing defects, the measured vibration signals are highly nondeterministic and non stationary. A major concern in bearing defect diagnostics is to capture the defect-induced impulses in the vibration signal and associate them with bearing kinematics with known operational conditions. 0888–3270/02/+$35.00/0
# 2002 Published by Elsevier Science Ltd.
1026
Q. SUN AND Y. TANG
Common techniques used for bearing defect detection include time and frequency domain analyses [2]. Statistical information of the time domain signal can be used as trend parameters [3–8]. They can provide information such as the energy level of the vibration signals and the shape of the amplitude probability distributions. High energy level of the vibration signal measured by the RMS values may indicate severely damaged components, while spikiness of the signal measured by higher order moments of the probability distribution function indicates incipient defects [4]. Calculations of the relative energy level of bearing vibrations improve the robustness of trending through time domain parameters to varying operational conditions [5]. However, effectiveness of applying the time domain trending parameters depends on bearing failure mode and its operational environment. On-set defects affect only a narrow frequency band. The overall peak amplitude will not be affected until the amplitude in the frequency band becomes the largest overall component [8], which implies poor sensitivity of the time domain analysis to incipient defects. Furthermore, time domain statistical parameters could be very sensitive to the environmental noise and therefore unreliable in detecting bearing defects. Other than focussing on the time series of the signal, spectrum analysis [9] provides power spectrum in the frequency domain. Major frequency components at which the power spectrum peaks can be identified and compared with the bearing characteristic frequencies to determine whether the bearing is experiencing a particular fault. The disadvantage of spectrum comparison, however, is that peaks in the power spectrum are usually very sensitive to fluctuations of the shaft speed. To overcome this problem, patterns of the power spectrum peaks as opposed to isolated frequency peaks are generated by cepstrum analysis [10] for trending. In recent years, envelope analysis [11, 12] has become another popular method used for bearing defect diagnostics. In particular, envelope analysis has been recognised as being an effective approach for detecting bearing incipient failure in high noise environment. Early version of the envelope analysis is through applying a band pass filter to the analog vibration signal centred at a structural resonance frequency. Rectification is then applied followed by a smoothing circuit to recover the enveloped signal. It has become more popular these days that digital signal processing techniques such as those through the Hilbert transform are applied instead. The Hilbert transform allows convenient construction of both the amplitude and phase demodulations of a signal. Envelope analysis provides a mechanism for extracting periodic excitation or demodulated defectinduced structural resonance. Sensitivity of the envelope analysis to incipient defects is improved using band pass filtering. However, one of the difficulties with envelope analysis is in determining the best frequency band which encompasses the frequency of structural resonance. Time frequency domain analysis, such as short time Fourier transform (STFT) is generally used for processing transient and non stationary signals [13]. However, STFT suffers from the fixed time frequency resolution. Wavelet transforms overcome the drawback of the STFT in that it provides a good compromise between location and frequency resolution [14]. A bigger window is used for low frequency estimates while a smaller window is used for high frequency estimates. It allows more accurate local descriptions of localised signal characteristics. Local analysis with wavelet transform can be performed by multiresolution formalism [15], which is designed to represent signals in finer and finer details. Special aspects of the signal such as trends, breakdown points, discontinuities in high derivatives and selfsimilarity are revealed on suitable scales. In addition, wavelet transform possesses the ability of filtering the polynomial behaviour to some predefined degree, which can be better illustrated by singularity analysis [16–21]. By detecting the modulus maxima in
SINGULARITY ANALYSIS
1027
the wavelet transform, singularity analysis provides an effective means in extracting sharp transient characteristics from the signal. In this light, we explore the application of singularity analysis for bearing defect diagnosis. The structure of the paper is organised as follows. In Section 2, we focus on the relevant aspects of the wavelet transformation with an emphasis on singularity analysis to characterise scale-free behaviour through the Lipschitz exponent. In Section 3, signal reconstruction from wavelet transform modulus maxima and feature extraction from the reconstructed signal will be discussed. Details of implementation with the application to bearing defect diagnosis and experimental results are presented in Section 4. Section 5 closes the paper with conclusions.
2. CONTINUOUS WAVELET TRANSFORM AND SINGULARITY ANALYSIS
Wavelet transform is best known for its ability of analysing the local behaviour of functions. In signal processing for pattern recognition, it is the irregularity not the smoothness in the signal that provides the most interesting and discriminating nature. Fourier transform, known as the global transformation assumes infinite energy in the signal and cannot be used to identify signals with short duration. short time Fourier transform uses fixed moving windows and therefore cannot achieve desirable resolutions in both time and frequency simultaneously. Wavelet transform overcomes the problem by using basis functions representing small waves with compact support, called wavelets. In particular, all the basis functions are obtained from their mother wavelet or kernel jðxÞ by scaling, that is, js ¼ ð1=sÞjðx=sÞ; where s 2 Rþ is called the scale factor which ‘adapts’ the width of the wavelet kernel to the required microscopic resolution and thus changing its frequency contents. For any signal f ðxÞ 2 L2 ðRÞ, the wavelet transform of the signal is performed through the convolution of the signal with the scaled and translated kernel: Z þ1 Wf ðs; xÞ ¼ f ðxÞ js ðxÞ ¼ f ðuÞjs ðx uÞ du ð1Þ 1
From the above equation, note that the wavelet transform Wf ðs; xÞ is a function of the scale s and the spatial position x. For one-dimensional vibration signals, x represents time. The plane defined by (s, x) is called the scale-space plane. As s decreases, we are looking into finer and finer details (with higher frequency contents) of the signal at a particular location x. A standard way of representing the continuous wavelet transform is to use a two- or three-dimensional plot called scalogram which plots the modulus of the wavelet transform jWf ðs; xÞj as functions of location x over a range of scales s. Figure 1(b) shows a two-dimensional scalogram of the vibration signal [Fig. 1(a)] obtained from a freight car bearing with seeded defects on the bearing outer raceway. The gray scale from black to white represents the increased magnitude of jWf ðs; xÞj. 2.1. FUNCTION REGULARITY AND THE LIPSCHITZ EXPONENT In mathematics, the local regularity of a function f ðxÞ at a particular point x0 is often measured by the Lipschitz exponent: j f ðxÞ f ðx0 Þj4Ajx x0 ja
ð2Þ
where A > 0; ðx; x0 Þ 2 a; b½2 ; 04a41: According to the above equation, the function is called uniformly Lipschitz a over the interval a; b½. The larger the value of a, the smoother the function f ðxÞ. Lipschitz 1 corresponds to a continuously differentiable function at point x0 ; Lipschitz a (05a51) means that function f ðxÞ is continuous at x0 but the derivative of the function at that point is discontinuous; Lipschitz 0 indicates the function
1028
Q. SUN AND Y. TANG
Figure 1. Scalogram and wavelet transform modulus.
being discontinuous but bounded in the neighbourhood of x0 and therefore the function is singular at x0 [16]. The above concept can be extended for values a > 1. For this case, we are interested in the order n up to which the function f ðxÞ is continuously differentiable at x0 . If there exists a polynomial pn ðxÞ of degree n and n5a4n þ 1, such that the following relation holds: j f ðxÞ pn ðxÞj4jx x0 ja
ð3Þ
the function is said Lipschitz a at x0 . The polynomial pn ðxÞ is often associated with the Taylor’s expansion of f ðxÞ around x0 but equation (3) holds even if such an expansion does not exist. If a positive Lipschitz a (a > 0) provides information on the degree of differentiability of a function, one would naturally expect that a negative Lipschitz a may reveal properties of singularity of a function f ðxÞ at a localised interval. Mallat and Hwang [16] defined the negative Lipschitz a for tempered distribution of finite order through its primitive and proved that if f ðxÞ is Lipschitz a then its primitive is Lipschitz a+1 and vice versa. This is extremely useful in signal processing dealing with noise. For instance, white noise is a distribution that is almost everywhere singular with a ¼ 12 e for any e > 0 [16]. A Dirac is uniformly Lipschitz-1 in the neighbourhood of 0 because its primitive is Lipschitz 0. A classical tool for measuring the global regularity of a function f ðxÞ is to look at the asymptotic decay of its Fourier transform f#ðoÞ [16]. To follow the same analogy of global vs local properties between the Fourier and wavelet transforms, the latter is sought to characterise the local regularity of the function at x0 . To illustrate this idea, let us assume
SINGULARITY ANALYSIS
1029
that function f ðxÞ can be characterised by the Lipschitz exponent at x0 , and can be locally described as f ðxÞ ¼ a0 þ a1 ðx x0 Þ þ þ an ðx x0 Þn þ Ajx x0 ja : Let us also assume a wavelet jðxÞ with at least n þ 1 vanishing moments, that is Z þ1 xk jðxÞ dx ¼ 0 8k; 04k5n þ 1:
ð4Þ
ð5Þ
1
Using equations (1), (4) and (5), the wavelet transform of f using jðxÞ at the fixed location x0 becomes Z þ1 Aju x0 ja js ðx0 uÞ du Wf ðs; x0 Þ ¼ 1 Z þ1 a ð6Þ ¼ As jnja js ðnÞdn / sa : 1
Therefore, proportionality between the wavelet transform and the Lipschitz exponent a exists. In fact, exact theorems have been developed [16] which provides a necessary and sufficient condition on the asymptotic decay of jWf ðs; xÞj when the scale s approaches zero: jWf ðs; xÞj4Ae sa
ð7Þ
where e is a positive constant and f ðxÞ 2 L2 ðRÞ is uniformly Lipschitz a over the intervals ðx; x0 Þ 2 a þ e; b e½2 . It is important to note that the above condition (7) is localised on intervals. At fine scales, the condition characterises the local regularity behaviour of the signal in the neighbourhood of x0 . 2.2. WAVELET TRANSFORM MODULUS MAXIMA FOR SINGULARITY DETECTION Condition (7) can be used for different applications with different purposes. It can be used to estimate the Lipschitz exponent a for determining the type of regularities of the signal at localised intervals. For machinery fault diagnosis, we are concerned with capturing certain types of singularities related to machinery conditions. During bearing operations, impact is generated if a defect exists which can be described by an impulsive force of short duration. According to Newton’s second law of motion, this results in a discontinuous acceleration response due to the impact. If the vibration signal is obtained from accelerometers mounted on the bearing housing, we expect to observe discontinuity directly in the signal, that is, a ¼ 0. If velocity pickups are used to obtain the vibration signal, singularities caused by the defect-induced impact exist at the derivative level of the velocity profile, that is, 05a51. Similarly, if proximity probes are used in signal acquisition, we are then interested in catching the locations of discontinuity of the second derivatives of the signal: 15a52. From the previous discussion, one may be able to detect local singularities of a function by estimating the ratio of decay of jWf ðs; xÞj across all scales. It is in fact not practical since it imposes measure of decay of jWf ðs; xÞj in the two-dimensional scale-space (s, x) in the neighbourhood of x0 which requires intensive computations. One solution was suggested to only measure the decay of jWf ðs; x0 Þj at fine scales. It is, however, not reliable due to the simplification of a two-dimensional measure into a one-dimensional measure [16]. A more practical and popular way of applying condition (7) has been developed in the past by various researchers [16, 17]. To introduce the idea of detecting singularities through the modulus maxima first proposed by Mallat and Hwang [16], let us assume
1030
Q. SUN AND Y. TANG
that f ðxÞ is Lipschitz 0 at some localised intervals. We use a smoothing function yðxÞ which can be viewed as the impulse response of a low pass filter. The scaled version can be defined as 1 x ys ¼ y : ð8Þ s s Discontinuity of f ðxÞ at x can be defined as local sharp variation point of f ðxÞ smoothed by ys ðxÞ. In other words, a singular point in f ðxÞ corresponds to the inflection point in f ys . We can define the wavelet as the derivative of the smoothing function jðxÞ ¼ The wavelet transform of f ðxÞ is given by
dyðxÞ : dx
dys d Wf ðs; xÞ ¼ f ðxÞ js ðxÞ ¼ f ðxÞ s ðxÞ ¼ s ðf ys ÞðxÞ: dx dx
ð9Þ
ð10Þ
From equation (10), the wavelet transform is proportional to the derivative of f smoothed by ys . The local maxima of jWf ðs; xÞj correspond to the inflection point with sharp variation on f ys . As well, values of jWf ðs; xÞj at local maxima increase as scale s increases. This is an important observation as it can be used to discriminate modulus maxima created by irrelevant noise from those by the signal. To implement the above idea, we trace the ridge of the wavelet transform modulus. If jWf ðs; xÞj has, besides x0 , N maxima at xj ; 14j4N, in the localised interval including x0 , we can connect these points at different scale s and generate the lines to form ridges. From equation (7), we note that these ridgelines will converge to the singularity point x0 as s ! 0. Figure 1(c) shows the ridgelines corresponding to the wavelet transform modulus as shown in Fig. 1(b). One could observe some long ridgelines covering the whole scale range and almost periodically, two of them converge to a single point at fine scales. Ridgelines that do not propagate to large scales are due to noise or high frequency structural resonance. It is these long ridgelines covering the whole scale space and their periodic patterns that reveal the relevant information about the health of machine components. In terms of choosing the wavelet kernel, a wavelet c with one vanishing moment is preferred to detect singularities with a ¼ 0. As such, singularities at the derivatives of the signal will not be detected which allows us to focus on the defect-related abrupt changes in the signal. A wavelet with two vanishing moments will be ideal for processing velocity signals and a wavelet with three vanishing moments should be employed for displacement signals.
3. SIGNAL RECONSTRUCTION
The property of exponential decay of the wavelet transform modulus maxima at the singular point as mentioned in the previous section is very useful for diagnosis applications. Mallat and Hwang showed that maxima caused by noise especially white noise with Gaussian distribution of certain variance do not propagate to large scales [16]. In this work, we noticed that values of the wavelet transform modulus intensify significantly as the scale increases. This provides the motivation for signal denoising through thresholding the wavelet transform maxima and recover the signal from the remaining information to enhance the interesting nature of the signal. Mathematical analysis of this inverse problem is quite involved and no proof has yet been found that such reconstruction is possible and stable. However, in this paper,
SINGULARITY ANALYSIS
1031
we wish to show, through numerical experiments, that it is indeed quite useful for bearing defect diagnosis. 3.1. SIGNAL RECOVERY FROM WAVELET TRANSFORM MAXIMA Although general and complete theorem does not exist, Kicey and Lennard showed that wavelet transform modulus maxima define a complete and stable signal representation [18], if the signal has a band-limited Fourier transform and if the Fourier transform of the wavelet has a complete support. For general dyadic wavelets, Meyer [19] and Berman and Baras [20] proved that exact reconstruction is not possible for general dyadic wavelets. However, signals with the same wavelet maxima differ from each other only slightly, which explains the success of numerical reconstruction. Suppose that c is a dyadic wavelet, which means that there exist A; B > 0 such that þ1 X # j 2 8o 2 R f0g; A4 ð11Þ cð2 oÞ 4B j¼1
j
the dyadic wavelet transform fWf ð2 ; xÞgj2Z is a complete and stable representation. Singularities create sequences of maxima that converge towards the corresponding location at fine scales, and the Lipschitz regularity is calculated from the amplitude of the wavelet transform modulus maxima. At each scale 2j , if the local maxima of Wf ð2j ; xÞ is located at fxj;p g, we are interested in the magnitude of Z þ1 x x 1 j;p j Wf ð2 ; xj;p Þ ¼ f cj;p ¼ pffiffiffiffi f ðxÞc dx: ð12Þ 2j 2j 1 Signal recovery is achieved by finding another function f* such that W f*ð2j ; xj;p Þ ¼ f* cj;p ¼ f cj;p
ð13Þ
and whose maxima are located at xj;p . Mallat and Zhong introduced the alternative orthogonal projection algorithm to obtain f* [16]. Cetin and Ansari proposed an iterative algorithm using projections onto convex sets with guaranteed convergence [22]. Several other algorithms have been proposed along the same line [23, 24]. 3.2. FEATURE EXTRACTION THROUGH THRESHOLDING In an effort to recover the signal from the wavelet transform maxima in order to highlight the singularities pertinent to the original signal, one could also enhance the strength of these singularities by modifying the amplitude of the maxima. In specific applications such as bearing defect diagnosis, prior knowledge about bearing kinematics and dynamics helps discriminating our attention to different signal singularities. We may even remove some singularities by suppressing the corresponding maxima if they are obviously unrelated to bearing defect characteristics. In this work, we observed that the singularities in the acceleration signal caused by periodic impact present increased the intensity of the wavelet transform modulus maxima as the analysing scale increases. In contrast, magnitude of the maxima corresponding to noise remains as smaller values and does not propagate to larger scales. Furthermore, singularities caused by high frequency structural resonance also do not propagate to large scales by the very nature of scale s being inversely proportional to the frequency. Therefore, as long as the wavelet transform is performed for large enough scale space, amplitude of the modulus maxima will be dominated by singularities due to low frequency periodic impulses. This observation
1032
Q. SUN AND Y. TANG
motivated us to propose the following thresholding scheme for signal denoising and recovery: Step 1: At each scale 2j , all wavelet transform coefficients that are not maxima of wavelet transform modulus are set to zero, that is ( Wf ð2j ; xp Þ > Wf ð2j ; xÞ8x=xp if f c j;p ð14Þ Wf ð2j ; xÞ ¼ 0 otherwise: j Stepj 2: At each scale 2 , remove all wavelet transform maxima whose amplitude Wf ð2 ; xÞ is below a preset threshold Tj : ( f cj;p if Wf ð2j ; xp Þ5Tj j ð15Þ Wf ð2 ; xÞ ¼ 0 otherwise:
Step 3: Signal is reconstructed through the inverse wavelet transform based on the modified wavelet transform modulus on all scales. The first step ensures the reconstruction algorithm to zoom in singular points in the signal. The second step aims at cancelling small signal variations. Small signal variations disappear because the corresponding maxima are removed, but sharp signal transitions corresponding to large maxima are not affected. We found it effective in this work that the threshold Tj can be determined corresponding to 90% of the signal energy. By adopting the two constraint conditions of equations (14) and (15), denoising can be achieved. After signal is reconstructed based on the modified wavelet transform maxima, only sharp transitional portion of the original signal is retained in the reconstructed signal.
4. APPLICATIONS TO BEARING DEFECT DIAGNOSIS
We apply the methodology of detecting singularities through wavelet transform maxima to bearing defect diagnosis. In particular, due to the availability of the experimental data, we focus on double row tapered rolling element bearings used in freight cars as shown in Fig. 2. Vibration signals were provided by the Association of American Railroads (AAR) with a bearing test right set up at the Transportation Technology Center (TTC) in Colorado for the development of efficient and reliable bearing fault diagnosis systems. We wish to show our analysis through two typical examples, that is, bearings with defects located on the outer raceways and on the inner raceways. Various experiments were performed to reflect bearing conditions such as different types of defects and their severity, as well as varying operating conditions such as axial loads on the bearings and their rotation speeds. Our objective is to extract bearing condition signature embedded in the vibration signals through detecting the abrupt change or impulsive nature in the signal relevant to bearing characteristic frequencies. We start with choosing the smoothing function with which the wavelet kernel can be derived. We use the Gaussian function with standard deviation s ¼ 1, that is 1 2 yðxÞ ¼ pffiffiffiffiffiffi ex =2 : ð16Þ 2p The derivative of the smoothing function yields the wavelet kernel: x 2 cðxÞ ¼ pffiffiffiffiffiffi ex =2 : 2p
ð17Þ
It can be verified that the above wavelet has just one vanishing moment. This is desirable since we want to detect singularities in the measured signal representing the acceleration profile. Such wavelet kernel has also been shown to be very effective for image edge
1033
SINGULARITY ANALYSIS
Seal wear ring
Roller
Spacer ring
End cap Inner race or cone
Outer ring or cup
Backing ring
Figure 2. Schematic diagram of railway freight car bearing.
Figure 3. Bearing with single cup spall.
detection in computer vision [17]. As numerical experiments, we wish to demonstrate the method through two typical cases: bearing with defects on the outer and inner raceways. 4.1. CASE 1: DIAGNOSIS OF DEFECTS ON THE OUTER RACEWAY We use the experimental data representing signals collected from bearings with ‘seeded’ defects on the outer raceway or cup. Figure 3 shows the outer raceway of type F bearing with a single spall representing an incipient defect. Outer race ring is fixed stationary on the bearing housing while inner race is shrink fit with the wheel shaft. Vibration signals are
1034
Q. SUN AND Y. TANG
collected from accelerometers mounted on the bearing housing. Tests were carried out with bearings under different operating conditions, such as wheel travelling speed (30–80 mph) and axial load representing full and empty carloads. When defect such as fatigue spall develops on the cup, it is likely to happen inside the load zone of the bearing [25]. According to bearing kinematics and dynamics, impact occurs each time when a roller in the bearing encounters the spall. Damped impulse response is expected from the vibration signal at the frequency of ball-passing-outerraceway fbpo [2]. However, impulses caused by the cyclic impact are obscure in the time series of the signal as shown in Fig. 4. These impulses are buried in the wideband noise consisting of high frequency machinery resonance and environmental noise. Using singularity analysis as introduced above, we first perform dyadic wavelet transform on the acceleration signal. Maxima locations can be identified and thresholded by Tj at scale 2j for j ¼ 0; 1; :::; 6 as introduced in Section 3.2 for denoising. Through the inverse wavelet transform of the thresholded modulus, the reconstructed signal after denoising is shown in Fig. 5. Distinct evenly spaced impulse clusters can be observed from the reconstructed signal. For this set of data, the wheel speed is 30 mph and the data were digitised with a sampling frequency of 13 kHz. Therefore, from Fig. 5, the frequency of the impulse clusters is identified as 47.8 Hz.Theoretical estimation for ball-passing-outer-raceway frequency from the bearing kinematics gives fbpo ¼ 48:2 Hz: It should be noted that the theoretical estimation of bearing characteristic frequencies is based on an idealised model with idealised operating patterns and conditions. In our case, several reasons should be counted for the discrepancy between the real and estimated bearing characteristic frequencies. First, the tested bearings consist of tapered rollers instead of balls as used in the theoretical model [2]. This will alter the kinematics and therefore the characteristic frequencies. Second, theoretical calculation assumes no slip which is almost impossible in bearing operations. Finally, manufacturing error and incorrect design prevent bearing geometry to be perfect as assumed. Consequently, theoretical estimations of bearing characteristic frequencies should be used only as guidelines instead of benchmarks for defect diagnosis. In this case, 0.8% error between the estimated and measured frequencies for ball-passing-outer-raceway places a strong indication that fatigue spalls had been
Figure 4. Vibration signal of bearing with single cup spall at a speed of 30 mph.
SINGULARITY ANALYSIS
1035
Figure 5. Reconstructed signal from Fig. 4.
developed on the bearing cup. In addition, it is worthwhile to stress from Fig. 5 that only defect-induced impulse clusters are retained in the reconstructed signal. It indicates the effectiveness of the proposed algorithm in cancelling out the high frequency structural resonance and environmental noise. The phenomenon of an impulse cluster instead of a single impulse emerging in the reconstructed signal can be explained with bearing dynamics. Due to localised bearing defects such as spalling, as shown in Fig. 3, sharp impact is generated during the operation, which further excites structural resonance. Under the effect of structural damping, the resulting vibration by the defect-induced impact consists of a leading impulse followed by exponentially decaying sinusoidal oscillations with high frequencies. Singularities exist in the acceleration profile of the vibration with larger intensity at the location of leading impulses. By focussing on detecting the discontinuity in the acceleration profile, we are able to identify the time locations of the leading impulses and therefore spatial locations of the defect. To analyse the severity of the bearing damage, the damped oscillation patterns provide more detailed information. The same bearing was tested when the wheel traveles at 80 mph. Although the impulsive nature can be observed from the time domain acceleration profile, no definite periodicity can be identified from the original signal as shown in Fig. 6. Compared with the original signal, impulse clusters can be observed as evenly distributed on the reconstructed signal as shown in Fig. 7. The frequency of these impulse clusters is estimated as 127.4 Hz. The theoretical value of fbpo for a wheel speed of 80 mph is 128.7 Hz. Again, these two frequency values are close enough, with 1% difference, to indicate spalling on the cup. We now show another test on a type E bearing with multiple spalls on the outer race as shown in Fig. 8. Types E and F differ from each other on geometric dimensions. With seeded multiple spalls, we intend to study the effect of more severely damaged bearings as compared to bearing with a single spall. Figure 9 shows the acceleration profile of the bearing at a wheel travelling speed of 70 mph. Samples are recorded at the same sampling frequency (13 kHz). The reconstructed signal from the wavelet transform modulus maxima again reliably reveals impulse clusters spaced regularly at the time intervals corresponding to the ball-passing-outer-race period as shown in Fig. 10. With multiple spalls, however,
1036
Q. SUN AND Y. TANG
Figure 6. Vibration signal of bearing with single cup spall at a speed of 80 mph.
Figure 7. Reconstructed signal from Fig. 6.
the measured response is a sum of impact-induced impulse response from every single spalls. The phase difference between the two adjacent impulses either causes cancellation or reinforcement of the bearing defect frequencies. This is shown from Fig. 10 by the very different intensity of impulse response each time when a roller strikes the spalling area. As damage to the bearing becomes more severe, the periodicity of these impulses is expected to decrease. In addition, the measured characteristic frequency may be different from the theoretical estimations due to the aforementioned reasons. The reconstructed signal in this case shows impulse clusters that appear with a frequency of 124.1 Hz while theoretical estimation gives a frequency of 118.3 Hz. This shows a 5% error on the characteristics frequency estimation at a wheel travelling speed of 70 mph. To further improve the
SINGULARITY ANALYSIS
1037
Figure 8. Bearing with multiple cup spalls.
Figure 9. Vibration signal of bearing with multiple cup spalls.
reliability of the method, pattern recognition as developed earlier by one of the authors maybe applied to process the reconstructed signals [5, 25]. 4.2. CASE 2: DIAGNOSIS OF DEFECTS ON THE INNER RACEWAY Figure 11 shows a picture of the tested bearing with multiple defects on the bearing inner raceway or cone. When defects occur on the bearing inner raceway, they rotate relative to the outer raceway in shaft rotation speed. Due to bearing clearance, impact is only generated when a cone defect rotates into the load zone [25] and is encountered by a passing roller. The characteristic frequency in this case is the ball-passing-inner-race frequency fbpi [2]. The predominant feature of the vibration signal as a response to defectinduced impact consists of harmonics with inner race defect characteristic frequencies modulated with the inner race or shaft rotating frequencies fi . Similar to defects on the
1038
Q. SUN AND Y. TANG
Figure 10. Reconstructed signal from Fig. 9.
Figure 11. Bearing with multiple cone spalls.
cup, multiple cone defects can be treated as the sum of many localised defects located inside the bearing loading zone at different phase angles around the shaft centre. Consequently, bearing defect frequencies will be reinforced or cancelled depending on their relative phase difference. Similar characteristics can also be found in the vibration signals of bearing with roller defects. The only difference is that the defective roller revolves relative to the outer raceway in cage frequency fc and the defects contact the cone and cup alternately. The defect characteristic frequency is thus twice of the roller-spinfrequency: 2fbsf . The approach is similar to that treating the cone defects. This involves identifying the defect characteristic frequencies, fbpi for the cone defect and 2fbsf for the broken roller, modulated with the inner race rotation (for cone defects) or cage rotation
SINGULARITY ANALYSIS
1039
Figure 12. Vibration signal of bearing with multiple cone spalls.
(for broken rollers). Therefore, for simplicity of illustration, we focus on treating the cone defect. We consider a type F bearing with multiple seeded cone spalls as shown in Fig. 11. The vibration signal illustrated in Fig. 12 was obtained with a sampling frequency of 2.6 kHz. The measurement was carried out at a wheel travelling speed of 40 mph and under an axial load of 33 k lbs representing full carloads. Repetitive peaks corresponding to rollers passing inner race frequency and modulation with shaft rotation can be observed from the time-domain acceleration profile. However, precise location of time and intensity of the impact cannot be measured precisely because of noise resulted from the high frequency structural resonance. In contrast, through singularity detection and thresholding denoising, defect-induced impact is clearly identified from the reconstructed signal as shown in Fig. 13. It shows the effectiveness of the de-noising property of the proposed algorithm. As such, frequencies of the dominant harmonics can be easily calculated. In our case, the frequency of impulses is 84.3 Hz, while the frequency of the impulse clusters is 5.9 Hz. Compared with values from theoretical estimations, the ball-passing-inner-race frequency is fbpi=79.3 Hz and the shaft rotation frequency is 6.2 Hz, generating a relative error of 6.3 and 4.8%, respectively. Once again considering the error factors, test results are close enough to the theoretical estimation to indicate cone spalls. It is ascertained now that the reconstructed signal shown in Fig. 13 is the characterised signal pattern of the bearing with multiple inner race defects. Comparing Fig. 12 with Fig. 13, the proposed algorithm shows a great promise in highlighting the defect-induced impact in the vibration signals for bearing fault diagnosis.
5. CONCLUSIONS
In this paper, we applied singularity analysis through the continuous wavelet transform to bearing defect diagnosis. We showed that through identifying the ridgelines in the wavelet transform modulus that propagate from fine to large scales and converge to singular points at fine scales, we can effectively capture the time location of impact in the
1040
Q. SUN AND Y. TANG
Figure 13. Reconstructed signal from Fig. 12.
bearing vibration signals. To increase the signal-to-noise ratio, we proposed to threshold the wavelet transform modulus maxima to retain 90% of the signal energy and reconstruct the signal from the remaining portion of the modulus to highlight singularities caused by defect-induced impact. As such, features can be extracted and associated with the bearing characteristic frequencies for fault classifications. The proposed method was tested by processing vibration signals of bearings with different seeded defects. We demonstrated through these tests that the proposed method is very effective in isolating the defect-related signature from the original signals and robust to bearing operating conditions. Further development would be to incorporate the proposed algorithm into bearing diagnosis systems. Simple processing of the denoised signal, such as frequency spectrum calculation can be performed to extract trending parameters, such as the dominant frequencies for diagnosis. The trending parameters can further be processed through pattern recognition for fault classification [5, 25]. It offers a great potential for the development of an automatic machinery fault diagnosis system.
ACKNOWLEDGEMENT
The authors would like to thank the National Research Council of Canada and the Association of American Railroads for providing the testing data.
REFERENCES 1. B. K. N. Rao 1996 Handbook of Condition Monitoring. Oxford: Elsevier Advanced Technology. 2. I. Howard 1994 Defense Science and Technology Organization Technical Report DSTO-RR0013, Australia. A review of rolling element bearing vibration}detection, diagnosis and prognosis 3. D. Dyer and R. M. Stewart 1978 ASME Journal of Mechanical Design 100, 229–235. Detection of rolling element bearing damage by statistical vibration analysis. 4. C. Q. Li and C. J. Pickering 1982 Condition Monitoring and Diagnostic Technology 2, 81–84. Robustness and sensitivity of non-dimensional amplitude parameters for diagnosis of fatigue spalling.
SINGULARITY ANALYSIS
1041
5. Q. Sun, F. Xi and G. Krishnappa 1998 Proceedings of the Canadian Society of Mechanical Engineering Forum, 423–429, Toronto. Signature analysis of rolling element bearing defects. 6. C. J. Li and J. Ma 1992 Proceedings of the 46th Meeting of the Mechanical Failures Prevention Group, 53–46, Virginia. Bearing localized defect detection through wavelet decomposition of vibrations. April 7–9. 7. H. R. Martin, F. Ismail and A. Sakuta 1992 Mechanical Systems and Signal Processing 6, 317–327. Algorithms for statistical moment evaluation for machine health monitoring. 8. R. B. Randall 1985 Maintenance Management International 5, 161–167. Computer aided vibration spectrum trend analysis for condition monitoring. 9. S. Braun 1986 Mechanical Signature Analysis. New York: Academic Press. 10. G. T. Zheng and W. J. Wang 2001 ASME Journal of Vibration and Acoustics 123, 222–229. A new cepstral analysis procedure of recovering excitations for transient components of vibration signals and applications to rotating machinery condition monitoring. 11. M. S. Darlow, R. H. Badgley and G. W. Hogg 1974 US Army Air Mobility Research and Development Laboratory Technical Report-74–77, October. Application of high frequency resonance techniques for bearing diagnostics in helicopter gearboxes. 12. P. W. Tse, Y. H. Peng and R. Yam 2001 ASME Journal of Vibration and Acoustics 123, 303–310. Wavelet analysis and envelope detection for rolling element bearing fault diagnosis}their effectiveness and flexibilities. 13. J. B. Allen and L. R. Rabiner 1977 Proceedings of the IEEE 65, 1558–1564. A unified approach to short time Fourier analysis and synthesis. 14. C. S. Burrus, R. A. Gopinath and H. T. Guo 1998 Introduction to Wavelets and Wavelet Transforms. Englewood Cliffs, NJ: Prentice-Hall. 15. S. G. Mallat 1989 IEEE Transactions on Pattern Recognition and Machine Intelligence 11, 674– 693. A theory for multiresolution signal decomposition: the Wavelet Representation. 16. S. G. Mallat and W. L. Hwang 1992 IEEE Transactions on Information Theory 38, 617–643. Singularity detection and processing with wavelets. 17. S. G. Mallat and S. Zhong 1992 IEEE Transactions on Pattern Recognition and Machine Intelligence 14, 710–732. Characterization of signals from multiscale edges. 18. C. J. Kicey and C. J. Lennard 1997 Fourier Analysis and Applications 3, 63–82. Unique reconstruction of band-limited signals by a Mallat–Zhong wavelet transform algorithm. 19. Y. Meyer 1993 Wavelets: Algorithms and Applications. Philadelphia: SIAM. [Translated and revised by R. D. Ryan]. 20. Z. Berman and J. S. Baras 1993 IEEE Transactions on Signal Processing 41, 3216–3231. Properties of the multiscale maxima and zero-crossings representations. 21. F. J. Herrmann 2001 Applied and Computational Harmonic Analysis 11, 64–88. Singularity characterization by monoscale analysis: application to seismic imaging. 22. A. E. Cetin and R. Ansari 1994 IEEE Transactions on Signal Processing 42, 194–196. Signal recovery from wavelet transform maxima. 23. R. Carmona 1995 Wavelets and Statistics, pp. 96–108. Berlin: Springer-Verlag. Extrema reconstruction and spline smoothing: variations on an algorithm of Mallat and Zhong. 24. Z. Cvetkovic and M. Vetterli 1995 IEEE Transactions on Signal Processing 43, 681–693. Discrete-time wavelet extrema representation: design and consistent reconstruction. 25. P. Chen, Q. Sun and F. Xi 2000 The IASTED International Conference on Applied Simulation and Modelling, Banff, Canada, 40–44, July 24–26. Feature extraction using segmentation analysis for bearing defect diagnosis.
SINGULARITYANALYSIS USING CONTINUOUS WAVELET TRANSFORM FOR BEARING FAULT DIAGNOSIS Q. Sun Department of Mechanical and Manufacturing Engineering, University of Calgary, 2500 University Dr., N.W., Calgary, Alberta, Canada T2N 1N4. E-mail: [email protected]
and Y. Tang Faculty of Mechanical Engineering, Beijing Science and Technology University, University Road 23, Beijing 100083, People’s Republic of China. E-mail: [email protected] (Received 11 November 2001, accepted 1 February 2002) In this paper, wavelet transform is applied to detect abrupt changes in the vibration signals obtained from operating bearings being monitored. In particular, singularity analysis across all scales of the continuous wavelet transform is performed to identify the location (in time) of defect-induced bursts in the vibration signals. Through modifying the intensity of the wavelet transform modulus maxima, defect-related vibration signature is highlighted and can be easily associated with the bearing defect characteristic frequencies for diagnosis. Due to the fact that vibration characteristics of faulty bearings are complex and defect-related vibration signature is normally buried in the wideband noise and high frequency structural resonance, simple signal processing cannot be used to detect bearing fault. We show, through experimental results, that the proposed method has the ability to discriminate noise from the signal significantly and is robust to bearing operating conditions, such as load and speed, and severity of the bearing damage. These properties are desirable for automatic detection of machine faults. # 2002 Published by Elsevier Science Ltd.
1. INTRODUCTION
Bearing defect diagnosis is a typical problem in machine and process fault detection and diagnosis [1]. Once defects occur on the working surface of bearing elements, impact is generated when mating elements encounter the defects. Sharp transient response accompanied by damped oscillation can often be observed in the measured vibration signals. The periodic transient behaviour in the vibration signal contains rich information about bearing health. On the other hand, however, defect-induced impulses are often buried in the wideband noise accompanied by high frequency structural resonance. This imposes a great deal of difficulty in identifying defect-induced impulses from the measured vibration signals. Depending on many other factors, such as machinery operating conditions and severity of bearing defects, the measured vibration signals are highly nondeterministic and non stationary. A major concern in bearing defect diagnostics is to capture the defect-induced impulses in the vibration signal and associate them with bearing kinematics with known operational conditions. 0888–3270/02/+$35.00/0
# 2002 Published by Elsevier Science Ltd.
1026
Q. SUN AND Y. TANG
Common techniques used for bearing defect detection include time and frequency domain analyses [2]. Statistical information of the time domain signal can be used as trend parameters [3–8]. They can provide information such as the energy level of the vibration signals and the shape of the amplitude probability distributions. High energy level of the vibration signal measured by the RMS values may indicate severely damaged components, while spikiness of the signal measured by higher order moments of the probability distribution function indicates incipient defects [4]. Calculations of the relative energy level of bearing vibrations improve the robustness of trending through time domain parameters to varying operational conditions [5]. However, effectiveness of applying the time domain trending parameters depends on bearing failure mode and its operational environment. On-set defects affect only a narrow frequency band. The overall peak amplitude will not be affected until the amplitude in the frequency band becomes the largest overall component [8], which implies poor sensitivity of the time domain analysis to incipient defects. Furthermore, time domain statistical parameters could be very sensitive to the environmental noise and therefore unreliable in detecting bearing defects. Other than focussing on the time series of the signal, spectrum analysis [9] provides power spectrum in the frequency domain. Major frequency components at which the power spectrum peaks can be identified and compared with the bearing characteristic frequencies to determine whether the bearing is experiencing a particular fault. The disadvantage of spectrum comparison, however, is that peaks in the power spectrum are usually very sensitive to fluctuations of the shaft speed. To overcome this problem, patterns of the power spectrum peaks as opposed to isolated frequency peaks are generated by cepstrum analysis [10] for trending. In recent years, envelope analysis [11, 12] has become another popular method used for bearing defect diagnostics. In particular, envelope analysis has been recognised as being an effective approach for detecting bearing incipient failure in high noise environment. Early version of the envelope analysis is through applying a band pass filter to the analog vibration signal centred at a structural resonance frequency. Rectification is then applied followed by a smoothing circuit to recover the enveloped signal. It has become more popular these days that digital signal processing techniques such as those through the Hilbert transform are applied instead. The Hilbert transform allows convenient construction of both the amplitude and phase demodulations of a signal. Envelope analysis provides a mechanism for extracting periodic excitation or demodulated defectinduced structural resonance. Sensitivity of the envelope analysis to incipient defects is improved using band pass filtering. However, one of the difficulties with envelope analysis is in determining the best frequency band which encompasses the frequency of structural resonance. Time frequency domain analysis, such as short time Fourier transform (STFT) is generally used for processing transient and non stationary signals [13]. However, STFT suffers from the fixed time frequency resolution. Wavelet transforms overcome the drawback of the STFT in that it provides a good compromise between location and frequency resolution [14]. A bigger window is used for low frequency estimates while a smaller window is used for high frequency estimates. It allows more accurate local descriptions of localised signal characteristics. Local analysis with wavelet transform can be performed by multiresolution formalism [15], which is designed to represent signals in finer and finer details. Special aspects of the signal such as trends, breakdown points, discontinuities in high derivatives and selfsimilarity are revealed on suitable scales. In addition, wavelet transform possesses the ability of filtering the polynomial behaviour to some predefined degree, which can be better illustrated by singularity analysis [16–21]. By detecting the modulus maxima in
SINGULARITY ANALYSIS
1027
the wavelet transform, singularity analysis provides an effective means in extracting sharp transient characteristics from the signal. In this light, we explore the application of singularity analysis for bearing defect diagnosis. The structure of the paper is organised as follows. In Section 2, we focus on the relevant aspects of the wavelet transformation with an emphasis on singularity analysis to characterise scale-free behaviour through the Lipschitz exponent. In Section 3, signal reconstruction from wavelet transform modulus maxima and feature extraction from the reconstructed signal will be discussed. Details of implementation with the application to bearing defect diagnosis and experimental results are presented in Section 4. Section 5 closes the paper with conclusions.
2. CONTINUOUS WAVELET TRANSFORM AND SINGULARITY ANALYSIS
Wavelet transform is best known for its ability of analysing the local behaviour of functions. In signal processing for pattern recognition, it is the irregularity not the smoothness in the signal that provides the most interesting and discriminating nature. Fourier transform, known as the global transformation assumes infinite energy in the signal and cannot be used to identify signals with short duration. short time Fourier transform uses fixed moving windows and therefore cannot achieve desirable resolutions in both time and frequency simultaneously. Wavelet transform overcomes the problem by using basis functions representing small waves with compact support, called wavelets. In particular, all the basis functions are obtained from their mother wavelet or kernel jðxÞ by scaling, that is, js ¼ ð1=sÞjðx=sÞ; where s 2 Rþ is called the scale factor which ‘adapts’ the width of the wavelet kernel to the required microscopic resolution and thus changing its frequency contents. For any signal f ðxÞ 2 L2 ðRÞ, the wavelet transform of the signal is performed through the convolution of the signal with the scaled and translated kernel: Z þ1 Wf ðs; xÞ ¼ f ðxÞ js ðxÞ ¼ f ðuÞjs ðx uÞ du ð1Þ 1
From the above equation, note that the wavelet transform Wf ðs; xÞ is a function of the scale s and the spatial position x. For one-dimensional vibration signals, x represents time. The plane defined by (s, x) is called the scale-space plane. As s decreases, we are looking into finer and finer details (with higher frequency contents) of the signal at a particular location x. A standard way of representing the continuous wavelet transform is to use a two- or three-dimensional plot called scalogram which plots the modulus of the wavelet transform jWf ðs; xÞj as functions of location x over a range of scales s. Figure 1(b) shows a two-dimensional scalogram of the vibration signal [Fig. 1(a)] obtained from a freight car bearing with seeded defects on the bearing outer raceway. The gray scale from black to white represents the increased magnitude of jWf ðs; xÞj. 2.1. FUNCTION REGULARITY AND THE LIPSCHITZ EXPONENT In mathematics, the local regularity of a function f ðxÞ at a particular point x0 is often measured by the Lipschitz exponent: j f ðxÞ f ðx0 Þj4Ajx x0 ja
ð2Þ
where A > 0; ðx; x0 Þ 2 a; b½2 ; 04a41: According to the above equation, the function is called uniformly Lipschitz a over the interval a; b½. The larger the value of a, the smoother the function f ðxÞ. Lipschitz 1 corresponds to a continuously differentiable function at point x0 ; Lipschitz a (05a51) means that function f ðxÞ is continuous at x0 but the derivative of the function at that point is discontinuous; Lipschitz 0 indicates the function
1028
Q. SUN AND Y. TANG
Figure 1. Scalogram and wavelet transform modulus.
being discontinuous but bounded in the neighbourhood of x0 and therefore the function is singular at x0 [16]. The above concept can be extended for values a > 1. For this case, we are interested in the order n up to which the function f ðxÞ is continuously differentiable at x0 . If there exists a polynomial pn ðxÞ of degree n and n5a4n þ 1, such that the following relation holds: j f ðxÞ pn ðxÞj4jx x0 ja
ð3Þ
the function is said Lipschitz a at x0 . The polynomial pn ðxÞ is often associated with the Taylor’s expansion of f ðxÞ around x0 but equation (3) holds even if such an expansion does not exist. If a positive Lipschitz a (a > 0) provides information on the degree of differentiability of a function, one would naturally expect that a negative Lipschitz a may reveal properties of singularity of a function f ðxÞ at a localised interval. Mallat and Hwang [16] defined the negative Lipschitz a for tempered distribution of finite order through its primitive and proved that if f ðxÞ is Lipschitz a then its primitive is Lipschitz a+1 and vice versa. This is extremely useful in signal processing dealing with noise. For instance, white noise is a distribution that is almost everywhere singular with a ¼ 12 e for any e > 0 [16]. A Dirac is uniformly Lipschitz-1 in the neighbourhood of 0 because its primitive is Lipschitz 0. A classical tool for measuring the global regularity of a function f ðxÞ is to look at the asymptotic decay of its Fourier transform f#ðoÞ [16]. To follow the same analogy of global vs local properties between the Fourier and wavelet transforms, the latter is sought to characterise the local regularity of the function at x0 . To illustrate this idea, let us assume
SINGULARITY ANALYSIS
1029
that function f ðxÞ can be characterised by the Lipschitz exponent at x0 , and can be locally described as f ðxÞ ¼ a0 þ a1 ðx x0 Þ þ þ an ðx x0 Þn þ Ajx x0 ja : Let us also assume a wavelet jðxÞ with at least n þ 1 vanishing moments, that is Z þ1 xk jðxÞ dx ¼ 0 8k; 04k5n þ 1:
ð4Þ
ð5Þ
1
Using equations (1), (4) and (5), the wavelet transform of f using jðxÞ at the fixed location x0 becomes Z þ1 Aju x0 ja js ðx0 uÞ du Wf ðs; x0 Þ ¼ 1 Z þ1 a ð6Þ ¼ As jnja js ðnÞdn / sa : 1
Therefore, proportionality between the wavelet transform and the Lipschitz exponent a exists. In fact, exact theorems have been developed [16] which provides a necessary and sufficient condition on the asymptotic decay of jWf ðs; xÞj when the scale s approaches zero: jWf ðs; xÞj4Ae sa
ð7Þ
where e is a positive constant and f ðxÞ 2 L2 ðRÞ is uniformly Lipschitz a over the intervals ðx; x0 Þ 2 a þ e; b e½2 . It is important to note that the above condition (7) is localised on intervals. At fine scales, the condition characterises the local regularity behaviour of the signal in the neighbourhood of x0 . 2.2. WAVELET TRANSFORM MODULUS MAXIMA FOR SINGULARITY DETECTION Condition (7) can be used for different applications with different purposes. It can be used to estimate the Lipschitz exponent a for determining the type of regularities of the signal at localised intervals. For machinery fault diagnosis, we are concerned with capturing certain types of singularities related to machinery conditions. During bearing operations, impact is generated if a defect exists which can be described by an impulsive force of short duration. According to Newton’s second law of motion, this results in a discontinuous acceleration response due to the impact. If the vibration signal is obtained from accelerometers mounted on the bearing housing, we expect to observe discontinuity directly in the signal, that is, a ¼ 0. If velocity pickups are used to obtain the vibration signal, singularities caused by the defect-induced impact exist at the derivative level of the velocity profile, that is, 05a51. Similarly, if proximity probes are used in signal acquisition, we are then interested in catching the locations of discontinuity of the second derivatives of the signal: 15a52. From the previous discussion, one may be able to detect local singularities of a function by estimating the ratio of decay of jWf ðs; xÞj across all scales. It is in fact not practical since it imposes measure of decay of jWf ðs; xÞj in the two-dimensional scale-space (s, x) in the neighbourhood of x0 which requires intensive computations. One solution was suggested to only measure the decay of jWf ðs; x0 Þj at fine scales. It is, however, not reliable due to the simplification of a two-dimensional measure into a one-dimensional measure [16]. A more practical and popular way of applying condition (7) has been developed in the past by various researchers [16, 17]. To introduce the idea of detecting singularities through the modulus maxima first proposed by Mallat and Hwang [16], let us assume
1030
Q. SUN AND Y. TANG
that f ðxÞ is Lipschitz 0 at some localised intervals. We use a smoothing function yðxÞ which can be viewed as the impulse response of a low pass filter. The scaled version can be defined as 1 x ys ¼ y : ð8Þ s s Discontinuity of f ðxÞ at x can be defined as local sharp variation point of f ðxÞ smoothed by ys ðxÞ. In other words, a singular point in f ðxÞ corresponds to the inflection point in f ys . We can define the wavelet as the derivative of the smoothing function jðxÞ ¼ The wavelet transform of f ðxÞ is given by
dyðxÞ : dx
dys d Wf ðs; xÞ ¼ f ðxÞ js ðxÞ ¼ f ðxÞ s ðxÞ ¼ s ðf ys ÞðxÞ: dx dx
ð9Þ
ð10Þ
From equation (10), the wavelet transform is proportional to the derivative of f smoothed by ys . The local maxima of jWf ðs; xÞj correspond to the inflection point with sharp variation on f ys . As well, values of jWf ðs; xÞj at local maxima increase as scale s increases. This is an important observation as it can be used to discriminate modulus maxima created by irrelevant noise from those by the signal. To implement the above idea, we trace the ridge of the wavelet transform modulus. If jWf ðs; xÞj has, besides x0 , N maxima at xj ; 14j4N, in the localised interval including x0 , we can connect these points at different scale s and generate the lines to form ridges. From equation (7), we note that these ridgelines will converge to the singularity point x0 as s ! 0. Figure 1(c) shows the ridgelines corresponding to the wavelet transform modulus as shown in Fig. 1(b). One could observe some long ridgelines covering the whole scale range and almost periodically, two of them converge to a single point at fine scales. Ridgelines that do not propagate to large scales are due to noise or high frequency structural resonance. It is these long ridgelines covering the whole scale space and their periodic patterns that reveal the relevant information about the health of machine components. In terms of choosing the wavelet kernel, a wavelet c with one vanishing moment is preferred to detect singularities with a ¼ 0. As such, singularities at the derivatives of the signal will not be detected which allows us to focus on the defect-related abrupt changes in the signal. A wavelet with two vanishing moments will be ideal for processing velocity signals and a wavelet with three vanishing moments should be employed for displacement signals.
3. SIGNAL RECONSTRUCTION
The property of exponential decay of the wavelet transform modulus maxima at the singular point as mentioned in the previous section is very useful for diagnosis applications. Mallat and Hwang showed that maxima caused by noise especially white noise with Gaussian distribution of certain variance do not propagate to large scales [16]. In this work, we noticed that values of the wavelet transform modulus intensify significantly as the scale increases. This provides the motivation for signal denoising through thresholding the wavelet transform maxima and recover the signal from the remaining information to enhance the interesting nature of the signal. Mathematical analysis of this inverse problem is quite involved and no proof has yet been found that such reconstruction is possible and stable. However, in this paper,
SINGULARITY ANALYSIS
1031
we wish to show, through numerical experiments, that it is indeed quite useful for bearing defect diagnosis. 3.1. SIGNAL RECOVERY FROM WAVELET TRANSFORM MAXIMA Although general and complete theorem does not exist, Kicey and Lennard showed that wavelet transform modulus maxima define a complete and stable signal representation [18], if the signal has a band-limited Fourier transform and if the Fourier transform of the wavelet has a complete support. For general dyadic wavelets, Meyer [19] and Berman and Baras [20] proved that exact reconstruction is not possible for general dyadic wavelets. However, signals with the same wavelet maxima differ from each other only slightly, which explains the success of numerical reconstruction. Suppose that c is a dyadic wavelet, which means that there exist A; B > 0 such that þ1 X # j 2 8o 2 R f0g; A4 ð11Þ cð2 oÞ 4B j¼1
j
the dyadic wavelet transform fWf ð2 ; xÞgj2Z is a complete and stable representation. Singularities create sequences of maxima that converge towards the corresponding location at fine scales, and the Lipschitz regularity is calculated from the amplitude of the wavelet transform modulus maxima. At each scale 2j , if the local maxima of Wf ð2j ; xÞ is located at fxj;p g, we are interested in the magnitude of Z þ1 x x 1 j;p j Wf ð2 ; xj;p Þ ¼ f cj;p ¼ pffiffiffiffi f ðxÞc dx: ð12Þ 2j 2j 1 Signal recovery is achieved by finding another function f* such that W f*ð2j ; xj;p Þ ¼ f* cj;p ¼ f cj;p
ð13Þ
and whose maxima are located at xj;p . Mallat and Zhong introduced the alternative orthogonal projection algorithm to obtain f* [16]. Cetin and Ansari proposed an iterative algorithm using projections onto convex sets with guaranteed convergence [22]. Several other algorithms have been proposed along the same line [23, 24]. 3.2. FEATURE EXTRACTION THROUGH THRESHOLDING In an effort to recover the signal from the wavelet transform maxima in order to highlight the singularities pertinent to the original signal, one could also enhance the strength of these singularities by modifying the amplitude of the maxima. In specific applications such as bearing defect diagnosis, prior knowledge about bearing kinematics and dynamics helps discriminating our attention to different signal singularities. We may even remove some singularities by suppressing the corresponding maxima if they are obviously unrelated to bearing defect characteristics. In this work, we observed that the singularities in the acceleration signal caused by periodic impact present increased the intensity of the wavelet transform modulus maxima as the analysing scale increases. In contrast, magnitude of the maxima corresponding to noise remains as smaller values and does not propagate to larger scales. Furthermore, singularities caused by high frequency structural resonance also do not propagate to large scales by the very nature of scale s being inversely proportional to the frequency. Therefore, as long as the wavelet transform is performed for large enough scale space, amplitude of the modulus maxima will be dominated by singularities due to low frequency periodic impulses. This observation
1032
Q. SUN AND Y. TANG
motivated us to propose the following thresholding scheme for signal denoising and recovery: Step 1: At each scale 2j , all wavelet transform coefficients that are not maxima of wavelet transform modulus are set to zero, that is ( Wf ð2j ; xp Þ > Wf ð2j ; xÞ8x=xp if f c j;p ð14Þ Wf ð2j ; xÞ ¼ 0 otherwise: j Stepj 2: At each scale 2 , remove all wavelet transform maxima whose amplitude Wf ð2 ; xÞ is below a preset threshold Tj : ( f cj;p if Wf ð2j ; xp Þ5Tj j ð15Þ Wf ð2 ; xÞ ¼ 0 otherwise:
Step 3: Signal is reconstructed through the inverse wavelet transform based on the modified wavelet transform modulus on all scales. The first step ensures the reconstruction algorithm to zoom in singular points in the signal. The second step aims at cancelling small signal variations. Small signal variations disappear because the corresponding maxima are removed, but sharp signal transitions corresponding to large maxima are not affected. We found it effective in this work that the threshold Tj can be determined corresponding to 90% of the signal energy. By adopting the two constraint conditions of equations (14) and (15), denoising can be achieved. After signal is reconstructed based on the modified wavelet transform maxima, only sharp transitional portion of the original signal is retained in the reconstructed signal.
4. APPLICATIONS TO BEARING DEFECT DIAGNOSIS
We apply the methodology of detecting singularities through wavelet transform maxima to bearing defect diagnosis. In particular, due to the availability of the experimental data, we focus on double row tapered rolling element bearings used in freight cars as shown in Fig. 2. Vibration signals were provided by the Association of American Railroads (AAR) with a bearing test right set up at the Transportation Technology Center (TTC) in Colorado for the development of efficient and reliable bearing fault diagnosis systems. We wish to show our analysis through two typical examples, that is, bearings with defects located on the outer raceways and on the inner raceways. Various experiments were performed to reflect bearing conditions such as different types of defects and their severity, as well as varying operating conditions such as axial loads on the bearings and their rotation speeds. Our objective is to extract bearing condition signature embedded in the vibration signals through detecting the abrupt change or impulsive nature in the signal relevant to bearing characteristic frequencies. We start with choosing the smoothing function with which the wavelet kernel can be derived. We use the Gaussian function with standard deviation s ¼ 1, that is 1 2 yðxÞ ¼ pffiffiffiffiffiffi ex =2 : ð16Þ 2p The derivative of the smoothing function yields the wavelet kernel: x 2 cðxÞ ¼ pffiffiffiffiffiffi ex =2 : 2p
ð17Þ
It can be verified that the above wavelet has just one vanishing moment. This is desirable since we want to detect singularities in the measured signal representing the acceleration profile. Such wavelet kernel has also been shown to be very effective for image edge
1033
SINGULARITY ANALYSIS
Seal wear ring
Roller
Spacer ring
End cap Inner race or cone
Outer ring or cup
Backing ring
Figure 2. Schematic diagram of railway freight car bearing.
Figure 3. Bearing with single cup spall.
detection in computer vision [17]. As numerical experiments, we wish to demonstrate the method through two typical cases: bearing with defects on the outer and inner raceways. 4.1. CASE 1: DIAGNOSIS OF DEFECTS ON THE OUTER RACEWAY We use the experimental data representing signals collected from bearings with ‘seeded’ defects on the outer raceway or cup. Figure 3 shows the outer raceway of type F bearing with a single spall representing an incipient defect. Outer race ring is fixed stationary on the bearing housing while inner race is shrink fit with the wheel shaft. Vibration signals are
1034
Q. SUN AND Y. TANG
collected from accelerometers mounted on the bearing housing. Tests were carried out with bearings under different operating conditions, such as wheel travelling speed (30–80 mph) and axial load representing full and empty carloads. When defect such as fatigue spall develops on the cup, it is likely to happen inside the load zone of the bearing [25]. According to bearing kinematics and dynamics, impact occurs each time when a roller in the bearing encounters the spall. Damped impulse response is expected from the vibration signal at the frequency of ball-passing-outerraceway fbpo [2]. However, impulses caused by the cyclic impact are obscure in the time series of the signal as shown in Fig. 4. These impulses are buried in the wideband noise consisting of high frequency machinery resonance and environmental noise. Using singularity analysis as introduced above, we first perform dyadic wavelet transform on the acceleration signal. Maxima locations can be identified and thresholded by Tj at scale 2j for j ¼ 0; 1; :::; 6 as introduced in Section 3.2 for denoising. Through the inverse wavelet transform of the thresholded modulus, the reconstructed signal after denoising is shown in Fig. 5. Distinct evenly spaced impulse clusters can be observed from the reconstructed signal. For this set of data, the wheel speed is 30 mph and the data were digitised with a sampling frequency of 13 kHz. Therefore, from Fig. 5, the frequency of the impulse clusters is identified as 47.8 Hz.Theoretical estimation for ball-passing-outer-raceway frequency from the bearing kinematics gives fbpo ¼ 48:2 Hz: It should be noted that the theoretical estimation of bearing characteristic frequencies is based on an idealised model with idealised operating patterns and conditions. In our case, several reasons should be counted for the discrepancy between the real and estimated bearing characteristic frequencies. First, the tested bearings consist of tapered rollers instead of balls as used in the theoretical model [2]. This will alter the kinematics and therefore the characteristic frequencies. Second, theoretical calculation assumes no slip which is almost impossible in bearing operations. Finally, manufacturing error and incorrect design prevent bearing geometry to be perfect as assumed. Consequently, theoretical estimations of bearing characteristic frequencies should be used only as guidelines instead of benchmarks for defect diagnosis. In this case, 0.8% error between the estimated and measured frequencies for ball-passing-outer-raceway places a strong indication that fatigue spalls had been
Figure 4. Vibration signal of bearing with single cup spall at a speed of 30 mph.
SINGULARITY ANALYSIS
1035
Figure 5. Reconstructed signal from Fig. 4.
developed on the bearing cup. In addition, it is worthwhile to stress from Fig. 5 that only defect-induced impulse clusters are retained in the reconstructed signal. It indicates the effectiveness of the proposed algorithm in cancelling out the high frequency structural resonance and environmental noise. The phenomenon of an impulse cluster instead of a single impulse emerging in the reconstructed signal can be explained with bearing dynamics. Due to localised bearing defects such as spalling, as shown in Fig. 3, sharp impact is generated during the operation, which further excites structural resonance. Under the effect of structural damping, the resulting vibration by the defect-induced impact consists of a leading impulse followed by exponentially decaying sinusoidal oscillations with high frequencies. Singularities exist in the acceleration profile of the vibration with larger intensity at the location of leading impulses. By focussing on detecting the discontinuity in the acceleration profile, we are able to identify the time locations of the leading impulses and therefore spatial locations of the defect. To analyse the severity of the bearing damage, the damped oscillation patterns provide more detailed information. The same bearing was tested when the wheel traveles at 80 mph. Although the impulsive nature can be observed from the time domain acceleration profile, no definite periodicity can be identified from the original signal as shown in Fig. 6. Compared with the original signal, impulse clusters can be observed as evenly distributed on the reconstructed signal as shown in Fig. 7. The frequency of these impulse clusters is estimated as 127.4 Hz. The theoretical value of fbpo for a wheel speed of 80 mph is 128.7 Hz. Again, these two frequency values are close enough, with 1% difference, to indicate spalling on the cup. We now show another test on a type E bearing with multiple spalls on the outer race as shown in Fig. 8. Types E and F differ from each other on geometric dimensions. With seeded multiple spalls, we intend to study the effect of more severely damaged bearings as compared to bearing with a single spall. Figure 9 shows the acceleration profile of the bearing at a wheel travelling speed of 70 mph. Samples are recorded at the same sampling frequency (13 kHz). The reconstructed signal from the wavelet transform modulus maxima again reliably reveals impulse clusters spaced regularly at the time intervals corresponding to the ball-passing-outer-race period as shown in Fig. 10. With multiple spalls, however,
1036
Q. SUN AND Y. TANG
Figure 6. Vibration signal of bearing with single cup spall at a speed of 80 mph.
Figure 7. Reconstructed signal from Fig. 6.
the measured response is a sum of impact-induced impulse response from every single spalls. The phase difference between the two adjacent impulses either causes cancellation or reinforcement of the bearing defect frequencies. This is shown from Fig. 10 by the very different intensity of impulse response each time when a roller strikes the spalling area. As damage to the bearing becomes more severe, the periodicity of these impulses is expected to decrease. In addition, the measured characteristic frequency may be different from the theoretical estimations due to the aforementioned reasons. The reconstructed signal in this case shows impulse clusters that appear with a frequency of 124.1 Hz while theoretical estimation gives a frequency of 118.3 Hz. This shows a 5% error on the characteristics frequency estimation at a wheel travelling speed of 70 mph. To further improve the
SINGULARITY ANALYSIS
1037
Figure 8. Bearing with multiple cup spalls.
Figure 9. Vibration signal of bearing with multiple cup spalls.
reliability of the method, pattern recognition as developed earlier by one of the authors maybe applied to process the reconstructed signals [5, 25]. 4.2. CASE 2: DIAGNOSIS OF DEFECTS ON THE INNER RACEWAY Figure 11 shows a picture of the tested bearing with multiple defects on the bearing inner raceway or cone. When defects occur on the bearing inner raceway, they rotate relative to the outer raceway in shaft rotation speed. Due to bearing clearance, impact is only generated when a cone defect rotates into the load zone [25] and is encountered by a passing roller. The characteristic frequency in this case is the ball-passing-inner-race frequency fbpi [2]. The predominant feature of the vibration signal as a response to defectinduced impact consists of harmonics with inner race defect characteristic frequencies modulated with the inner race or shaft rotating frequencies fi . Similar to defects on the
1038
Q. SUN AND Y. TANG
Figure 10. Reconstructed signal from Fig. 9.
Figure 11. Bearing with multiple cone spalls.
cup, multiple cone defects can be treated as the sum of many localised defects located inside the bearing loading zone at different phase angles around the shaft centre. Consequently, bearing defect frequencies will be reinforced or cancelled depending on their relative phase difference. Similar characteristics can also be found in the vibration signals of bearing with roller defects. The only difference is that the defective roller revolves relative to the outer raceway in cage frequency fc and the defects contact the cone and cup alternately. The defect characteristic frequency is thus twice of the roller-spinfrequency: 2fbsf . The approach is similar to that treating the cone defects. This involves identifying the defect characteristic frequencies, fbpi for the cone defect and 2fbsf for the broken roller, modulated with the inner race rotation (for cone defects) or cage rotation
SINGULARITY ANALYSIS
1039
Figure 12. Vibration signal of bearing with multiple cone spalls.
(for broken rollers). Therefore, for simplicity of illustration, we focus on treating the cone defect. We consider a type F bearing with multiple seeded cone spalls as shown in Fig. 11. The vibration signal illustrated in Fig. 12 was obtained with a sampling frequency of 2.6 kHz. The measurement was carried out at a wheel travelling speed of 40 mph and under an axial load of 33 k lbs representing full carloads. Repetitive peaks corresponding to rollers passing inner race frequency and modulation with shaft rotation can be observed from the time-domain acceleration profile. However, precise location of time and intensity of the impact cannot be measured precisely because of noise resulted from the high frequency structural resonance. In contrast, through singularity detection and thresholding denoising, defect-induced impact is clearly identified from the reconstructed signal as shown in Fig. 13. It shows the effectiveness of the de-noising property of the proposed algorithm. As such, frequencies of the dominant harmonics can be easily calculated. In our case, the frequency of impulses is 84.3 Hz, while the frequency of the impulse clusters is 5.9 Hz. Compared with values from theoretical estimations, the ball-passing-inner-race frequency is fbpi=79.3 Hz and the shaft rotation frequency is 6.2 Hz, generating a relative error of 6.3 and 4.8%, respectively. Once again considering the error factors, test results are close enough to the theoretical estimation to indicate cone spalls. It is ascertained now that the reconstructed signal shown in Fig. 13 is the characterised signal pattern of the bearing with multiple inner race defects. Comparing Fig. 12 with Fig. 13, the proposed algorithm shows a great promise in highlighting the defect-induced impact in the vibration signals for bearing fault diagnosis.
5. CONCLUSIONS
In this paper, we applied singularity analysis through the continuous wavelet transform to bearing defect diagnosis. We showed that through identifying the ridgelines in the wavelet transform modulus that propagate from fine to large scales and converge to singular points at fine scales, we can effectively capture the time location of impact in the
1040
Q. SUN AND Y. TANG
Figure 13. Reconstructed signal from Fig. 12.
bearing vibration signals. To increase the signal-to-noise ratio, we proposed to threshold the wavelet transform modulus maxima to retain 90% of the signal energy and reconstruct the signal from the remaining portion of the modulus to highlight singularities caused by defect-induced impact. As such, features can be extracted and associated with the bearing characteristic frequencies for fault classifications. The proposed method was tested by processing vibration signals of bearings with different seeded defects. We demonstrated through these tests that the proposed method is very effective in isolating the defect-related signature from the original signals and robust to bearing operating conditions. Further development would be to incorporate the proposed algorithm into bearing diagnosis systems. Simple processing of the denoised signal, such as frequency spectrum calculation can be performed to extract trending parameters, such as the dominant frequencies for diagnosis. The trending parameters can further be processed through pattern recognition for fault classification [5, 25]. It offers a great potential for the development of an automatic machinery fault diagnosis system.
ACKNOWLEDGEMENT
The authors would like to thank the National Research Council of Canada and the Association of American Railroads for providing the testing data.
REFERENCES 1. B. K. N. Rao 1996 Handbook of Condition Monitoring. Oxford: Elsevier Advanced Technology. 2. I. Howard 1994 Defense Science and Technology Organization Technical Report DSTO-RR0013, Australia. A review of rolling element bearing vibration}detection, diagnosis and prognosis 3. D. Dyer and R. M. Stewart 1978 ASME Journal of Mechanical Design 100, 229–235. Detection of rolling element bearing damage by statistical vibration analysis. 4. C. Q. Li and C. J. Pickering 1982 Condition Monitoring and Diagnostic Technology 2, 81–84. Robustness and sensitivity of non-dimensional amplitude parameters for diagnosis of fatigue spalling.
SINGULARITY ANALYSIS
1041
5. Q. Sun, F. Xi and G. Krishnappa 1998 Proceedings of the Canadian Society of Mechanical Engineering Forum, 423–429, Toronto. Signature analysis of rolling element bearing defects. 6. C. J. Li and J. Ma 1992 Proceedings of the 46th Meeting of the Mechanical Failures Prevention Group, 53–46, Virginia. Bearing localized defect detection through wavelet decomposition of vibrations. April 7–9. 7. H. R. Martin, F. Ismail and A. Sakuta 1992 Mechanical Systems and Signal Processing 6, 317–327. Algorithms for statistical moment evaluation for machine health monitoring. 8. R. B. Randall 1985 Maintenance Management International 5, 161–167. Computer aided vibration spectrum trend analysis for condition monitoring. 9. S. Braun 1986 Mechanical Signature Analysis. New York: Academic Press. 10. G. T. Zheng and W. J. Wang 2001 ASME Journal of Vibration and Acoustics 123, 222–229. A new cepstral analysis procedure of recovering excitations for transient components of vibration signals and applications to rotating machinery condition monitoring. 11. M. S. Darlow, R. H. Badgley and G. W. Hogg 1974 US Army Air Mobility Research and Development Laboratory Technical Report-74–77, October. Application of high frequency resonance techniques for bearing diagnostics in helicopter gearboxes. 12. P. W. Tse, Y. H. Peng and R. Yam 2001 ASME Journal of Vibration and Acoustics 123, 303–310. Wavelet analysis and envelope detection for rolling element bearing fault diagnosis}their effectiveness and flexibilities. 13. J. B. Allen and L. R. Rabiner 1977 Proceedings of the IEEE 65, 1558–1564. A unified approach to short time Fourier analysis and synthesis. 14. C. S. Burrus, R. A. Gopinath and H. T. Guo 1998 Introduction to Wavelets and Wavelet Transforms. Englewood Cliffs, NJ: Prentice-Hall. 15. S. G. Mallat 1989 IEEE Transactions on Pattern Recognition and Machine Intelligence 11, 674– 693. A theory for multiresolution signal decomposition: the Wavelet Representation. 16. S. G. Mallat and W. L. Hwang 1992 IEEE Transactions on Information Theory 38, 617–643. Singularity detection and processing with wavelets. 17. S. G. Mallat and S. Zhong 1992 IEEE Transactions on Pattern Recognition and Machine Intelligence 14, 710–732. Characterization of signals from multiscale edges. 18. C. J. Kicey and C. J. Lennard 1997 Fourier Analysis and Applications 3, 63–82. Unique reconstruction of band-limited signals by a Mallat–Zhong wavelet transform algorithm. 19. Y. Meyer 1993 Wavelets: Algorithms and Applications. Philadelphia: SIAM. [Translated and revised by R. D. Ryan]. 20. Z. Berman and J. S. Baras 1993 IEEE Transactions on Signal Processing 41, 3216–3231. Properties of the multiscale maxima and zero-crossings representations. 21. F. J. Herrmann 2001 Applied and Computational Harmonic Analysis 11, 64–88. Singularity characterization by monoscale analysis: application to seismic imaging. 22. A. E. Cetin and R. Ansari 1994 IEEE Transactions on Signal Processing 42, 194–196. Signal recovery from wavelet transform maxima. 23. R. Carmona 1995 Wavelets and Statistics, pp. 96–108. Berlin: Springer-Verlag. Extrema reconstruction and spline smoothing: variations on an algorithm of Mallat and Zhong. 24. Z. Cvetkovic and M. Vetterli 1995 IEEE Transactions on Signal Processing 43, 681–693. Discrete-time wavelet extrema representation: design and consistent reconstruction. 25. P. Chen, Q. Sun and F. Xi 2000 The IASTED International Conference on Applied Simulation and Modelling, Banff, Canada, 40–44, July 24–26. Feature extraction using segmentation analysis for bearing defect diagnosis.