Abstract
Order tracking has become one of the most effective methods for fault detection of rotating machinery under the time-varying shaft speed conditions. The transient phase estimation is very important for order tracking, especially when the tachometer installation is not convenient. The transient phase is usually obtained by integrating the instantaneous frequency (IF), so the IF estimation has attracted a great deal of concerns. This article describes a new IF estimation method based on the joint application of the synchrosqueezing transform (SST) and the multiscale chirplet path pursuit (MSCPP) method. The SST method as its high frequency resolution merits is used to estimate the frequency parameters for the parameter settings of the MSCPP method, that will resolve the high computation problem of the MSCPP method to a certain degree, so as to extensively use the high accuracy of the MSCPP method in IF estimation. The order spectrum based on the estimated IF can provide the demodulation information for the bearing fault diagnosis. The performance of the proposed method has been validated by both simulation and experimental data.
1. Introduction
Bearings are widely used in various rotating machinery. Their failure could cause machine breakdowns, production loss, and even catastrophic accidents. Thus, bearing fault detection has been a very important topic for industry and received a lot of academic attentions [1, 2]. Various FFT-based demodulation methods, such as Hilbert transform (HT) and Teager energy operator (TEO), have been commonly used for bearing fault detection [3]. However, under the time-varying rotating speed cases, the speed variation may cause “smearing” of the discrete frequencies in the frequency domain, so that the fault-related frequencies may no longer manifest as discrete frequency lines, and thus the fault features cannot be easily detected. That is to say, the traditional FFT-based demodulation analysis is no longer effective for bearing fault detection in the time-varying speed cases. However, the nonstationary vibration signals gathered under the time-varying speed conditions, especially in startup or shutdown process of the rotating machinery, which often contains more useful information that reflects certain fault symptoms [4, 5]. Fortunately, through the equal angle-interval resampling, order tracking technique can convert the nonstationary signals (in time domain) to stationary signals (in angle domain), that enables the above FFT-based demodulation methods to be revalidated for fault detection in time-varying speed cases [6]. It is therefore that order tracking has become one of the most commonly used methods for fault detection of rotating machinery under the time-varying shaft speed cases.
For order tracking technique, the equal angle-interval resampling should resample the nonstationary vibration signals at a rate proportional to the shaft rotating speed, which means that the shaft rotating speed should be known beforehand. The straightforward way to obtain the shaft rotating speed is to install speed measuring devices on the rotating machinery, such as tachometers and encoders. Nevertheless, due to the design reason or cost concerns, these speed measuring devices are not always allowed to be installed. As a consequence, it is necessary to extract the instantaneous rotating speed (IRS) information from the vibration signals. So far, there have been many IRS estimation related researches reported in the literature. Gao et al. put forward a method of order tracking based on instantaneous frequency (IF) estimation [6]. In detail, the IF is estimated from vibration signal via time-frequency spectrum peak searching algorithm (TFSPSA), and thus the IRS is estimated by means of polynomial fitting of the extracted IF; then, numerical integration and interpolation arithmetic are utilized to implement the equal angle-interval resampling. Most of the proposed order tracking methods are based on the abovementioned idea and framework; just the IF extraction approaches are different. For example, Li et al. extracted the IF based on the joint use of empirical mode decomposition (EMD) and Hilbert transform [7]. Ding et al. estimated the IF based on energy centrobaric correction method [8]. In [5, 9], the multiscale chirplets and sparse signal decomposition are joint-introduced into the estimation of instantaneous gear meshing frequency under the speed fluctuation state. Likewise with [6], polynomial fitting, numerical integration, and interpolation arithmetic are in turn used to accomplish the order spectrum analysis (OSA) of the nonstationary vibration signal in [5, 10].
Besides the order tracking technique, time-frequency analysis (TFA), as it helps reveal insight into the complex structure of the nonstationary vibration signals [11], is widely used in fault detection of rotating machinery under time-varying speed working condition independently. Kim et al. carried out a comparative research on damage detection in speed-up and coastdown process of grinding spindle-typed rotor-bearing system by use of short time Fourier transform (STFT), Winger-Ville distribution (WVD), and discrete wavelet transform (DWT) [12]. Meltzer and Ivanov studied the application of Cohen-class time-frequency approaches to fault detection in gear drives with nonstationary rotational speed [13]. Dong and Chen used Wigner–Ville spectrum based on cyclic spectral density (CSWVS) for extracting bearing fault patterns, when the bearing signals are under influences of random noise and gear vibrations [14]. Cai and Hu put forward a gear fault feature extraction method based on Sobel operator and WHT [15]. A review about the recent advances in time-frequency analysis methods for machinery fault diagnosis has been studied and published in 2013 [16]. In recent 5 years, synchrosqueezing transform (SST) [17] as a time-frequency analysis technique has received extensive attention. Li et al. successively adopted SST for fault diagnosis of gearbox under nonstationary conditions [18–20].
The most commonly used TFA method is the linear TFA which is expressed as the inner product of signal with a dictionary of time-frequency atoms, such as STFT and DWT. However, the TF resolution of this class of TFA methods is limited by the Heisenberg uncertainty principle; that is, the best time localization and highest frequency resolution cannot be achieved simultaneously [19]. In addition, they lack adaptability to match the complicated components since the basis in either the STFT or WT is fixed [21]. Another category of TFA method is bilinear TFA, which includes WVD and its modified variants. Though WVD has desirable theoretical characteristics, it is subjected to negative effect of cross-terms due to multicomponents [18]. Various modified variants have attempted to exploit the best TF resolution of WVD and meanwhile suppress the negative effect of cross-terms [22, 23]. However, they will compromise with TF resolution and autointegrity. In addition to linear and bilinear TFA, Hilbert–Huang transform (HHT) is another well-known TFA method. However, HHT is susceptible to singularity signal and may produce pseudointrinsic mode functions, thus interfering the correct time-frequency distribution [24]. Whether One TFA method is effective or not is due to its capability of concentrating the energy of TFR at or around the true IF in the TF plane. For this, Daubechies et al. proposed a wavelet based time-frequency enhancement method, that is, synchrosqueezing transform (SST) [17]. This technique eliminates the blur in frequency (scale) dimension via squeezing the TFR along the frequency (scale) axis, thus acquiring a better frequency resolution. However, the SST method only improves the frequency resolution; there still exists diffusion in time domain, and as a wavelet based approach for TFA enhancement, it is still limited as stated by the Heisenberg uncertainty principle. In recent years, to overcome the aforementioned diffusion problem in time domain, Dong et al. have successively put forward the improved methods [14–16], and as stated these improved SST methods have been used in fault diagnosis of gearbox under time variable speed conditions.
The TFA of vibration signals can provide the time-frequency structure of signal components. However, even ignoring the time-frequency ambiguity caused by the theoretical defects of various TFA methods, there is no TFA method that can visually reveal the multiples relationship of the frequencies containing vibration signals, and the amplitude ratio of the signal components also cannot be visually observed; these are not conducive to quantitative and qualitative identification of bearing faults. In view of the description above, comparing with the TFA methods, order spectrum can display the fault symptoms more intuitively. Therefore, this paper applies the order spectrum analysis for bearing fault detection.
As mentioned before, without any speed measuring devices installed, the IRS (or the IF) information should be accurately extracted from vibration signals for the equal angle-interval resampling of the order spectrum analysis. The TFSPSA has been widely used to extract the IF from signals. However, vibration signal always contains multicomponents, and the peaks of its time-frequency spectrum do not always appear in the same signal component, and the peaks always correspond to the local maximum value, not the global maximum. In other words, the IF extraction methods based on the TFSPSA are not adapted to all types of vibration signals.
The multiscale chirplets path pursuit (MSCPP) method was first put forward to detect nonstationary signal from very noisy data by Candès et al. [25], which has been extended to process the multicomponent nonstationary signals by Feng et al. [11]. It uses the linear frequency of the chirplets to adaptively approximate the IF curve piece by piece. To some extent, the MSCPP method can be called a flexible and adaptive method for the IF estimation. First of all, thanks to the multiscale dividing time interval characteristic, the MSCPP method can flexibly select shorter chirplet atoms whenever the IF curve is complex and shift to select longer chirplet templates if the IF becomes simpler [5]. Moreover, the MSCPP method does not require any prior knowledge about the phase and amplitude of the analyzed signal [5, 10]. The MSCPP method has a high IF estimation accuracy which has been verified in [10]. Nevertheless, the MSCPP method needs to search the best matching chirplets in the whole time-frequency plane, that will incur a huge computational cost and cannot meet the real-time requirement of the engineering application. Therefore, the SST method as it has higher frequency resolution is introduced to rough estimation of the IF information, which helps to select the appropriate parameters for the MSCPP method and thus avoids the blindness of global search, that will improve the computational efficiency. Consequently, the MSCPP method and the SST method are combined to accurately estimate the IF for the OSA of the bearing vibration signals collected under time-varying speed cases.
The paper hereafter is organized as follows. Section 2 presents the proposed method, including brief overviews of MSCPP and SST and detailed introduction to the proposed OSA method. And simulation and experiment studies on the application of the proposed method for bearing fault detection are, respectively, reported in Sections 3 and 4. The conclusion is drawn in Section 5.
2. The Presentation of the Proposed Methodology
The proposed method involves three basic steps. First, the frequency parameters are extracted via the SST method, which offers precondition for the second step. Second, the IF is estimated based on the MSCPP method. Like the first step, the second step provides precondition for the final step. Finally, the order spectrum analysis is completed based on the estimated IF. The overall structure of the proposed method is displayed in Figure 1. And the algorithms of the SST, MSCPP, and OSA methods are successively detailed in the following.
2.1. The Synchrosqueezing Transform (SST) Method [21]
The SST method is a reallocation method for the time-frequency spectrum, whose aim is to enhance the TFR. In order to understand the idea of the SST method, let us consider such signals with the general form , where represents noise and each component is an oscillating function with smoothly time-varying frequency and amplitude. For such signals , the SST approach can be performed as the following three steps.
The first step starts from the Continuous Wavelet Transform (CWT) of the signal defined as , where is the scale and is the time offset and is the chosen wavelet which is concentrated on the positive frequency axis: for ( is the Fourier transform of wavelet function ). For the purely harmonic signal , .
The second step is to calculate the FM-demodulated frequency . A candidate IF estimation for the signal is deduced by the differential operation as [21]. For the purely harmonic signal , we have .
The final step is reassigning energy in the time-scale plane to the time-frequency plane according to the map . Define the divisions with and the frequency divisions with . The SST of signal is defined only at the centers of the frequency bins , by the summing different contributions as .
2.2. The Multiscale Chirplets Path Pursuit (MSCPP) Method [5, 10]
The main idea of the MSCPP method is to use the chirplet atoms whose instantaneous frequency is linear to adaptively approximate the nonstationary signals, piece by piece. These nonstationary signals can be expressed as , where is amplitude, is initial phase, is the noise, and is the continuous and derivable instantaneous phase of the signal . The chirplet atoms adopted in the MSCPP method is from the chirplet dictionary which can be expressed as a family of functions with the following form: where is the dyadic time interval with the binary scale division defined as , is the total sampling time, is the number of sampling points, is the time-scale coefficient, and , ; then, . Moreover, and are the slope and offset coefficients, respectively, which depend on the scale coefficient . From (1), it is easy to think that the IF of a chirplet atom is linear and equals . According to the sampling theorem, should be less than ( is the sampling rate). is a rectangular window function, which is when and 0 when . The amplitude coefficient is the normalization factor that makes .
For each time interval , the chirplet atom which has the highest correlation to signal is selected from the chirplet dictionary; this process is accomplished by computing the maximum projection coefficient in each time interval; that is, ( is the inner product operator). Actually, includes the amplitude and initial phase information of the signal component in the time interval . If we denote as the signal component expressed by in the time interval , then can be represented as follows:
The MSCPP algorithm can be divided as the following two steps.
Step 1. Select a chirplet that is closest to signal in each time interval and then obtain the corresponding signal component .
Step 2. Use the best path algorithm [22] to form a matching signal whose time span is equal to that of the analytic signal and whose energy reaches the maximum among all the possible paths. For better understanding, this step can be mathematically expressed as follows.
Denote to be a set in which every element represents a possible connecting path; that is, the combined time interval corresponding to each connecting path should cover the whole analysis time period and with no overlap. The formed signal should satisfy , where signal is connected by the elements from and the connecting path is .
From the foregoing discussion, for a multicomponent signal, such as a bearing vibration signal, the signal component formed by the MSCPP method is the component with the largest energy in the original vibration signal. In the meantime, the IF of this signal component can be estimated by concatenating the piecewise linear frequencies of the chirplets.
2.3. The Order Spectrum Analysis (OSA)
The OSA method mainly contains the following three steps.
Step 1 (polynomial fitting of the extracted IF). Suppose , , are the discrete coordinates of the IF estimated by the MSCPP method. is the number of sampling points in time domain. Then, a polynomial function, such as a simple quadric function , is used as a curve fitting function. And based on Least Squares Method (LSM) [6], the fitting coefficients , , and can be calculated.
Step 2 (computing the resampling times). Supposing , resampling times, , is the number of sampling points in frequency domain; that is, equals to FFT length of order analysis. (, the max order range) is constant angular increment. Then, the numerical integration equation is established. Thus, the resampling times can be calculated with the equation given as .
Step 3 (resampling in angular domain). The resampling times are used to do the angle interpolation to the original sampling data , and then the equiangular resampling signal is gotten. At last, the FFT about the envelope of resampled signal is performed to get the order envelope demodulation spectrum.
3. Simulation Validation
In the case of a failure, the bearing vibration signals always manifest as repeated decaying pulse FM signals. And the modulation frequencies (rotational frequency or fault characteristic frequency) are changing with the time-varying rotational speed. Therefore, to validate the proposed approach, a decay pulse signal is constructed to imitate the impulse resulting from faults. And the rotational frequency and the fault characteristic frequency are, respectively, expressed as and . The decay pulse signal and the modulation frequencies and are, respectively, shown in Figures 2 and 3. The vibration signal mainly contains two FM components, namely, the fault characteristic frequency modulation component and the rotational frequency modulation component. Furthermore, when the failure occurs, the magnitude of the former is always higher than that of the latter, so in this study, the amplitude ratio is set as 1.5 : 1. Suppose the sampling frequency is 1024 Hz and the number of sampling points is 8192. Then, the waveform of the simulated vibration signal without noise is displayed in Figure 4. With the MATLAB function awgn.m, the noise is added to the simulated vibration signal with ; then, the waveform of the polluted signal is exhibited in Figure 5.
(a)
(b)
The time-frequency spectrum of the noisy vibration signal is obtained with the SST method and the result is shown in Figure 5(a). From Figure 5, we can see that there exists a curve with relative deeper color in the red dotted box. In order to discern the time-frequency curve more clearly, the contour plot is applied and exhibited in Figure 5(b). See that in Figure 5(b) the curve in the same red dotted box is more clearer than that in Figure 5(a), and the frequency range of the curve is about [85 Hz, 240 Hz].
Based on Figure 5(b), in the following MSCPP method, the frequency search range is reduced from to , and the maximum frequency difference between two connected chirplets in the joint point is set as 10 Hz. Consider the sampling point of 8192; the scale coefficient is equal to 13. Combined with Figure 5, the variation of the curve in the dotted box is gentler, so there is no need to set the scale coefficient from 0 to 13. That is to say, there is no need to excessively divide the time domain, and that will greatly increase the amount of calculation and thus cause the waste of time. In this study, the scale coefficient is set from 8 to 11. Using the MSCPP method with the parameter settings above, the signal component with the largest energy is formed by piecewise connecting the chirplet atoms, and the IF curve is simultaneously formed, respectively, (see Figures 6 and 7).
Based on the extracted IF (in Figure 7), the quadric polynomial fitting, numerical integration, and cubic spline interpolation (see details in Section 2.3) are performed in turn to obtain the equiangular resampling signal. Setting the max order range as 100, the envelope signal of the equiangular resampling signal and its corresponding order envelope demodulation spectrum are, respectively, shown in Figures 8 and 9. Although the equiangular resampling signal is stationary, influenced by the modulation and noise, the fault cannot be easily identified from Figure 8. However, the fault characteristic frequency can be seen clearly in Figure 9. To further observe the order spectrum in Figure 9, the order-axis display range is narrowed to and , respectively, shown in Figures 10(a) and 10(b). In Figure 10(a), except for the fault characteristic frequency , the order envelope demodulation spectrum is dominated by the rotational frequency and its multiples. Furthermore, the fault characteristic frequency and its second and third harmonics are highlighted in Figure 10(b). The abovementioned phenomenon is consistent with the simulated vibration signal. Researchers can recognize the faults according to the order spectrum.
(a)
(b)
4. Experiment Tests
The effectiveness of the proposed method is further investigated using the experimental data of faulty bearings. The experimental setup is shown in Figure 11. The rotor is well balanced, to respectively make a groove with 0.15 mm width and 0.13 mm depth on the outer race and the inner race of two normal bearings via laser cutting technology. The parameters of faulty bearing are given in Table 1. In Table 1, BPFO (ball pass frequency, outer race) and BPFI (ball pass frequency, outer race) represent the fault characteristic frequency of the outer race and the inner race, respectively, and is the shaft rotational frequency. The sampling frequency is set as 1024 Hz; the sampling time is 4 s.
4.1. Experiment 1 (Detection of the Inner Race Fault)
The vibration acceleration signal and the corresponding rotation speed signal of the inner race bearing are shown in Figures 12(a) and 12(b), respectively. The time-frequency spectrum of the vibration signal obtained by the SST method is exhibited in Figure 13(a). From Figure 13(a), we can see that the frequency range of interest is from 0 to 100 Hz. Therefore, the time-frequency spectrum in the frequency range is displayed in Figure 13(b). In Figure 13(b), the peaks of the time-frequency spectrum do not appear in the same signal component. There are some peaks existing in Curve.1, and there are also some peaks existing in other curves, such as the ellipses in Curve.2. Therefore, the IF extraction methods based on the time-frequency spectrum peak searching algorithm [4] are invalid in this study.
(a)
(b)
(a)
(b)
According to Figure 13(b), in the following MSCPP method, the frequency search range is reduced to , and the maximum frequency difference between two connected chirplets in the joint point is set as 5 Hz, and the scale coefficient is set from 8 to 11. Based on such parameters settings, the IF of the signal component with largest energy is extracted and presented in Figure 14. According to the rotation speed signal in Figure 12(b), the extracted IF in Figure 14 is approximated to the fault characteristic frequency; this result coincides well with the actual. In actual, once fault occurs; the signal component with the largest energy should be the component corresponding to the fault characteristic frequency. However, since the energy of the signal component corresponding to the fault characteristic frequency is the global maximum, and not the local maximum, the fault characteristic frequency curve in Figure 13 is not obvious.
Based on the extracted IF (in Figure 14), the vibration signal in Figure 12(a) is resampled and transformed to the stationary signal in angular domain (for the detailed algorithm, see Section 2.3). Figure 15 shows the envelope signal of the resampled signal. In the resampling process above, the max order is set as 50. The order envelope demodulation spectrum is exhibited in Figure 16. To further determine the fault, the order-axis display range in Figure 17 is reduced to . Figure 17 shows that the order envelope demodulation spectrum is dominated by the shaft rotational frequency , the fault characteristic frequency BPFI, and their respective multiples.
4.2. Experiment 2 (Detection of the Outer Race Fault)
The waveform of the outer race vibration acceleration signal and the corresponding shaft speed signal are, respectively, plotted in Figures 18(a) and 18(b). And the time-frequency spectrum obtained by means of the SST method is exhibited in Figure 19(a). In Figure 19(a), except for a few peaks of the time-frequency spectrum, most peaks are located in Curve.1. Combining the shaft speed signal shown in Figure 18(b) with the parameters information , the Curve.1 in Figure 19(a) may be the fault characteristic frequency curve BPFO. For further research, the MSCPP method is applied to extract the fault characteristic frequency curve. According to the time-frequency spectrum in Figure 19(a), the frequency search range is set as [30 Hz, 80 Hz], and the maximum frequency difference between two connected chirplets in the joint point is set as 10 Hz, and the scale coefficient is set from 8 to 11. The fault characteristic frequency curve is extracted from the vibration signal and shown with blue line in Figure 19(b). From Figure 19(b), we can see that the extraction frequency curve is very close to the true fault characteristic frequency curve. In the following, using the algorithm in Section 2.3, the vibration signal is resampled and the order envelope demodulation spectrum is obtained and shown in Figure 20. From Figure 20, we can see the demodulation frequencies are mainly the fault characteristic frequency BPFO and its multiples. So the outer race fault is easily recognized.
(a)
(b)
(a)
(b)
5. Conclusions
In this paper, the SST method and the MSCPP algorithm have been jointly applied to extract the IF information for the order spectrum analysis of the bearing nonstationary vibration signals And then the modulation information can be accurately captured from the envelope demodulation spectrum, which provides critical information for bearing fault diagnosis.
The main advantage of the proposed method developed by combined use of the SST method and the MSCPP algorithm includes the following: (a) comparing with the time-frequency spectrum (can be obtained by SST) peak searching algorithm, the proposed method has a more extensive application range in the IF estimation, because the MSCPP method as an algorithm based on the global search is effective in the IF estimation, even when the local maximum phenomenon exists in the time-frequency spectrum, which always makes the time-frequency spectrum peak searching algorithm useless; (b) based on the time-frequency spectrum obtained by the SST method, the frequency parameter settings of the MSCPP method can be more reasonable, which avoids the time waste of the blind global computation.
Nevertheless, there are some deficiencies inevitably existing in the proposed method, including the following: (a) the proposed method is valid only when the shaft speed curve is smooth. That is because the resampling process is achieved via polynomial interpolations, which makes the proposed method helpless when the shaft speed varies with the sinusoidal or the jump curves. (b) The proposed method is suitable for the shaft speed varies in a relative small range; otherwise, the high computation will prevent the timely detection of faults.
Competing Interests
The authors declare that they have no competing interests.
Acknowledgments
This work was supported by Natural Science Foundation of China (Grant nos. 51605405, 51605406), Science & Technology Innovation Project of Fujian Province, China (Grant no. 2016H2003), and Natural Science Foundation of Fujian Province, China (Grant nos. 2014J05065, 2014H0049, and 2015J01672). Their support is gratefully acknowledged.