Abstract
Condition monitoring and compound fault diagnosis are crucial key points for ensuring the normal operation of rotating machinery. A novel method for condition monitoring and compound fault diagnosis based on the dual-kurtogram algorithm and multivariate statistical process control is established in this study. The core idea of this method is the capability of the dual-kurtogram in extracting subbands. Vibration data under normal conditions are decomposed by the dual-kurtogram into two subbands. Then, the spectral kurtosis (SK) of Subband I and the envelope spectral kurtosis (ESK) of Subband II are formulated to construct a control limit based on kernel density estimation. Similarly, vibration data that need to be monitored are constructed into two subbands by the dual-kurtogram. The SK of Subband I and the ESK of Subband II are calculated to derive statistics based on the covariance determinant. An alarm will be triggered when the statistics exceed the control limit and suitable subbands for square envelope analysis are adopted to obtain the characteristic frequency. Simulation and experimental data are used to verify the feasibility of the proposed method. Results confirm that the proposed method can effectively perform condition monitoring and fault diagnosis. Furthermore, comparison studies show that the proposed method outperforms the traditional control chart, envelope analysis, and empirical mode decomposition.
1. Introduction
Rolling bearings are important components of rotating machinery. However, many failures of these bearing-host systems often occurred by bearing faults, and, thus, bearing condition monitoring is crucial for guaranteeing the safe operation of systems. However, actual vibration signals, particularly weak fault vibration signals, are frequently contaminated by different degrees of noise [1]. Therefore, effective resonance demodulation for extracting fault information with narrow bands is imperative for condition monitoring and diagnosis. Simultaneously, compound faults often occur under actual working conditions, so the diagnosis of compound faults is indispensable. Overall, the feature extraction of weak faults, real-time monitoring to provide a warning, and compound fault diagnosis are important for ensuring the normal operation of rotating machinery.
The health of a machine is highly dependent on the strength and reliability of its bearings, and different defects in bearings may arise during their manufacture. Hence, detecting such defects is highly necessary for the condition monitoring and quality inspection of bearings [2]. Popular statistical parameters in the time and frequency domains are frequently used to monitor the condition of rolling element bearings [3–5]. Therefore, the use of statistical process control to establish a control chart with these statistics is highly effective. However, single variable monitoring is far from sufficient under complex working conditions. Thus, conducting condition monitoring to provide early warning is indispensable to prevent heavy losses from defects, and this task can be accomplished through multivariate statistical process control (MSPC). MSPC is currently used in various fields, but it is rarely applied to the condition monitoring of rotating machinery. Hadian and Rahimifard [6] utilized multivariate statistical control charts and a complex process capability index to monitor project duration and costs, providing the manager with reliable information between cost and actual performance. Prastyo et al. [7] used multivariate control charts to monitor the water quality of a factory, effectively improving the high false alarm rate of the previous water quality monitoring technique. He et al. [8] applied grayscale images and MSPC to a manufacturing workshop to effectively monitor production. Some researchers have applied MSPC to rotating machinery. Jung [9] proposed a wavelet-based nonparametric multivariate control chart that solved the problem of parameter distribution assumptions; however, the bootstrap method he used in his study may introduce estimated deviations. Henneberg et al. [10] considered his study as the first implementation of control charts on the condition monitoring of thruster gears. Given the excellent performance of MSPC in other fields, the current authors are keen to introduce it into rotating machinery.
However, the vibration signals of rotating machinery are frequently nonstationary and submerged in considerable noises. To address these problems, several issues need to be considered when constructing the control charts: (1) theoretical distribution of the signals with noises and (2) constructing suitable outliers submerged by strong noises. To estimate the appropriate probability density function and overcome the high false alarm rate, the kernel density estimation (KDE) can be used to provide an accurate control limit for the charts. Some studies have adjusted the control limit of Hotelling’s chart to improve its performance in monitoring data with unknown distribution by using nonparametric techniques [11, 12]. The adaptability of KDE techniques can reduce the false alarm rate of a multivariate control chart. On the other hand, the fault signal will be overwhelmed by a large number of noises and other interference components, resulting in inconspicuous fault characteristics and overdue fault prediction. We need to extract the resonance narrowband for enhancing the feature, and the method introduced in the next section can enhance the fault feature, so the details of the method are provided in the next section.
Diagnosing faults is imperative after condition monitoring. Envelope analysis [13] can effectively extract the characteristic frequency of faulty bearing signals, but it is susceptible to interference components and noises. To address this issue, narrowband analysis is typically used to improve the performance of envelope analysis. A kurtogram [14, 15] exhibits good performance in narrowband signal analysis. The spectral kurtosis (SK) is a function of the short-time Fourier transform (STFT) window and it is used to locate the sensitive fault frequency band of vibration signals. However, a kurtogram is difficult to apply in engineering because of its low computational efficiency. To improve the efficiency of the kurtogram, Antoni proposed an algorithm called fast kurtogram [16]. Calculation efficiency is considerably improved by using a series of finite impulse response filter banks (FIR) based on a 1/3-binary tree structure. Fast kurtogram changes the selection of the optimal demodulation frequency band on the basis of the researcher’s experience. Therefore, the introduction of a fast kurtogram strongly supports the selection of the optimal frequency band. Considering the performance of fast kurtogram, a series of improved optimal demodulation frequency band algorithms based on fast kurtogram have been proposed thereafter. Lei et al. [17] suggested using the wavelet packet transform instead of STFT or FIR filter to enhance the decomposition effect because the former can process nonstationary signals more effectively than STFT. Wang and Liang [18] proposed an improved algorithm that uses window superposition technology to find the optimal filter. This method can adaptively estimate the center frequency and bandwidth. Chen et al. [19] designed an improved kurtogram using redundant second-generation wavelet packet transform and correlated kurtosis; this kurtogram can effectively extract fault features from the local frequency band and prevent frequency aliasing. In this method, the criterion parameter is determined by calculating the frequency band of the narrowband time-domain signal. When the criterion parameter is used to determine the optimal frequency band, the filter is designed on basis of SK. Although SK yields valid results under certain conditions, it also fails in some cases, for example, in the presence of a relatively strong, non-Gaussian noise with high peaks or a relatively high repetition rate of fault impulses. To solve this problem, Barszcz and Jablonsk [20] proposed protrugram, in contrast with a fast kurtogram, protrugram calculates the kurtosis of the envelope spectrum amplitude to replace SK, effectively locating the best resonance frequency band for high-impact faults. From the preceding discussion, fast kurtogram and protrugram can be combined. In this manner, a narrowband with increased fault information is established, improving the diagnosis of faults.
Motivated by the aforementioned issues, an improved dual-kurtogram-based control chart is proposed to provide early and effective warning and diagnosis of compound faults. In contrast with the conventional kurtogram, the dual-kurtogram algorithm can simultaneously diagnose high and low-density impact faults [21]. Simultaneously, dual-kurtogram with MSPC can provide an early warning of fault. Therefore, the dual-kurtogram can achieve the purposes of condition monitoring and fault diagnosing rotating machinery.
The remainder of this paper is structured as follows. Section 2 summarizes Hotelling’s control chart and the KDE algorithm. The concept of the dual-kurtogram algorithm is interpreted in Section 3. The combined dual-kurtogram and multivariate Hotelling’s control chart, that is, for the proposed method, is presented in Section 4. Section 5 validates the proposed method through simulations and experiments. The conclusions are drawn and future work directions are provided in Section 6.
2. Multivariate Hotelling’s Control Chart
2.1. Conventional Hotelling’s Control Chart
Hotelling’s control chart [22] is one of the multivariate control charts that can be used for process monitoring. The control chart is calculated for and under variables: is the average value; and is the covariance matrix.
Each of the batches of sample data is with a sample size of and the mean value and sample variance of the th variable test data in the th batch of data arewhere represents the th observation value of the th variable in the th batch of data.
The covariance between the th and th variables in the th batch of data is
Then, the statistics , and are averaged for batches of data and variables: are elements of the vector and the sample covariance matrix is
The statistics can be calculated as follows:
Supposing the data follow a multivariate normal distribution, the control limit of Hotelling’s can be obtained using the following equation:where is the distribution function, is the number of observations, is the number of variables, and is the false alarm rate. The process is considered in control of the statistics in equation (5) lower than the upper control limit.
However, the vibration signal of a bearing is submerged by strong noises and does not follow a normal distribution, Thus, we need to perform enhancing features to avoid strong noises interference and use KDE to calculate the control limit to ensure the accuracy of the control limit.
2.2. Constructing Statistics with Dual-Kurtogram
The faulty signal is submerged by strong noises and other interference components, which will lead the fault characteristics not obvious and the fault indication not in time. Considering the above circumstance, the kurtogram algorithm is often used for enhancing features to make the fault characteristics more obvious. Therefore, an improved dual-kurtogram algorithm mentioned is used to enhance the fault characteristics, and the values of SK and ESK in the subbands extracted by the dual-kurtogram are, respectively, constructed into sequences to construct statistics with equation (5).
2.3. KDE-Based Control Limit
The probability density function of an unknown random variable can be estimated by applying the KDE method. This nonparametric technique was first introduced by Rosenblatt [23] and Parzen [24] and, thus, it is sometimes called the Rosenblatt–Parzen kernel density estimator, which is an extension of the histogram estimator.
KDE has been used to estimate the distribution of statistics [25]. Let be a Hotelling’s statistic obtained under in-control condition. The distribution of statistics can be calculated using the following kernel function:where and define the kernel function and the estimated smoothing parameter, respectively. One of the most used kernels is the Gaussian kernel, which is also adopted in the current analysis. The Gaussian kernel is defined as
Therefore, the control limit of statistics based on KDE can be estimated by taking the percentile of the kernel distribution. Hence, the control limit of statistics based on KDE is equal to the -th percentile of distribution which can be calculated as follows:where is the empirical distribution function and is the significance level.
From the proceeding introduction, signals are filtered by the dual-kurtogram algorithm for extracting SK and ESK to construct a Hotelling’s control chart applying to condition monitoring.
3. Dual-Kurtogram Algorithm
Fast kurtogram [16] is an effective tool for the resonance demodulation of nonstationary signals. It can select a frequency band with more bearing fault information through filtering, enhancing bearing characteristics, and reducing noise interference. The core idea is to design filter banks with different hierarchical decomposition bandwidths to obtain the SK of the filtered signals and identify the fault resonance frequency band with the most abundant fault information. This method is widely used in the fault diagnosis of rotating machinery.
3.1. Fast Kurtogram
In the case of nonstationary signals, SK is a highly effective indicator. Antoni [15] designed a kurtogram based on STFT to filter the resonance frequency band of rolling bearings, realizing early fault diagnosis of rolling bearings.
Letting the signal be , STFT can be expressed as can be regarded as a complex envelope in terms of frequency and time, and SK based on STFT at frequency can be expressed asConsidering the presence of stationary added noise , the SK can be further expressed aswhere is the SK of signal with stationary added noise and is the signal-to-noise ratio (SNR).
However, a kurtogram based on STFT is required to calculate SK under nearly all combinations of window length and center frequency, necessitating a considerable amount of calculation. Antoni [16] proposed using a series of FIR filter banks instead of STFT to improve calculation efficiency. Such improvement is called fast kurtogram, and its basic steps are as follows: Step 1: a low-pass filter with a cut-off frequency of , is provided. A high-pass filter and a low-pass filter are constructed on basis of . The expression is The bandwidths of filters are and , respectively, where is sample frequency. Adopting the tree structure shown in Figure 1, the signal is repeatedly decomposed from top to bottom using the aforementioned high- and low-pass filters. Meanwhile, the number of decomposition layers is , and subbands can be obtained for each layer. Let the ith subbands of the kth layer be . Its center frequency is , its bandwidth is , and ’s adjacent lower subbands are Step 2: by constructing a three-division filter bank to make the frequency band perfect, the principle is similar to that of a two-division filter bank, which converts only the two divisions of each layer into three divisions. Thereafter, the th-level three-division filter is inserted into the th and th two-division filter to produce a one-to-one correspondence. The resulting filter bank structure becomes a 1/3-binary tree structure, as shown in Figure 2. Step 3: from the formulas of and , the real and imaginary parts of subsignal differ by 90° regardless of whether it is a high-pass or low-pass filter. By the principle of envelope demodulation, is the envelope of the analytic signal, and, thus, the SK of the subsignal can be expressed as The SK of all the subsignals is calculated, and the subband used for demodulation can be obtained by selecting the center frequency and bandwidth corresponding to the maximum SK.
3.2. Fast Computation of Protrugram
After the proposal of a fast kurtogram, improvements to this algorithm have been proposed. Although the quality of the demodulation frequency band can be improved to a certain extent, this band is susceptible to the impact of impulse noise irrelevant to bearing faults. To overcome this defect, Barszcz and Jablonski [20] proposed the use of ESK instead of SK to overcome the problem in which the signal is in the presence of a relatively strong non-Gaussian noise with high peaks or a relatively high repetition rate of fault impulses.
Protrugram and fast kurtogram exhibit two major differences: (1) the metric is converted from SK to ESK and (2) the variable frequency band is converted into a fixed narrowband. Considering the two differences, we use ESK instead of SK and apply it directly to fast kurtogram instead of the fixed narrowband in protrugram. Accordingly, the highest ESK of subbands can be calculated directly using the algorithm.
3.3. Dual-Kurtogram Algorithm
Fast kurtogram has been widely used in previous studies. However, it can hardly extract the fault frequency band in the case of compound faults, and this subband is frequently found in only one fault. When a fast kurtogram is used in the study of compound fault characteristics, it can effectively locate the low-density impact fault subband and the high-impact density fault subband will be submerged [26]. To compensate for this defect, a protrugram in 1/3-binary tree form [27] is used in this algorithm. This algorithm can locate the high-impact density fault subband in the experiment. Then, a dual-kurtogram algorithm is constructed for extracting two subbands and detecting compound faults effectively.
The dual-kurtogram algorithm can locate the best resonance frequency bands with low- and high-impact density faults. Considering its advantages and those of the dual-kurtogram algorithm, combining them and the extracted subband will produce the largest SK and ESK. Then, the largest SK and the largest ESK should be used to construct statistics. The improved dual-kurtogram-based control chart algorithm is used to monitor and provide early warning of weak faults.
4. An Improved Dual-Kurtogram-Based Control Chart
Bearing vibration signals can effectively reflect the health of bearings, and they are widely used in fault detection. These signals are obtained during the whole operation of a bearing. Many statistical parameters or characteristics of vibration signals can be used to evaluate the health of bearings qualitatively and quantitatively.
In this proposed method, spectral kurtosis (SK) and envelope spectral kurtosis (ESK) are utilized as two time-frequency domain indicators. To examine the sensitivity of the two indicators to the characteristic frequency of faults, different signals simulating bearing outer race fault with different characteristic frequencies are conducted. The fault characteristic frequency ranges from 20 Hz to 120 Hz, the step size is 0.1 Hz, other parameters are with the same values except for the fault characteristic frequency, the sampling frequency is set as 12000 Hz, and the time length is 0.3 s. The analysis result is shown in Figure 3.
In Figure 3, the horizontal axis represents the fault characteristic frequency of the stimulation signal, the vertical axis represents the normalized value, and the normalized value is calculated by dividing the actual value (the value of SK and ESK) by the maximum of these values. As shown in Figure 3, the spectral kurtosis shows a downward trend with the decrease of fault characteristic frequency, and the envelope spectral kurtosis shows an upward trend with the decrease of fault characteristic frequency. Simultaneously, the variation trend of fault characteristic frequency with the normalized value of SK and ESK represents the density of its impact. Through the above analysis, it can be concluded that the spectral kurtosis is more sensitive to low-density impact. The envelope spectral kurtosis is more sensitive to high-density impact. Considering the proposed dual-kurtogram method, Subband I can be calculated with the highest spectral kurtosis value, and Subband II can be calculated with the highest envelope spectral kurtosis value which can be extracted respectively. Therefore, in practical application, the fast kurtogram can recognize low-density impact while the fast computation of protrugram can identify high-density impact. Given that the dual-kurtogram algorithm focuses on SK and ESK, SK and ESK are ideal indicators for measuring the health of bearings in this study. MSPC exerts a good effect on the real-time monitoring of multiple indicators. Therefore, fast kurtogram and fast computation of protrugram are combined to find a reasonable fusion rule for monitoring faults under low- and high-impact density faults in rolling bearings, realizing the monitoring and diagnosis of compound faults in rolling bearings.
The flowchart of the proposed fault monitoring and diagnosis method involving two phases is shown in Figure 4. Data are gathered and analyzed to construct the control chart in Phase I. Normal bearing vibration signals are applied to the dual-kurtogram algorithm, and two band-pass filters are designed to obtain the subbands. The SK value is obtained from Subband I, and the ESK is obtained from Subband II. The covariance matrix was calculated by two values of each subband. Then, statistics are constructed by a covariance matrix and used to calculate the control limit with KDE.
In Phase II, the algorithm in Phase I is repeated to obtain SK and ESK, and the covariance matrix was calculated to construct statistics. The statistics are compared with the control limit. When the statistic exceeds the control limit, the subbands used to construct the statistic are extracted, and envelope analysis is performed on these subbands. The fault characteristic frequencies corresponding to the two subbands are obtained. In this way, the fault location of the rolling bearings can be obtained.
5. Experimental Verification
To verify the introduced method, simulations and experimental tests are conducted in our study. The results positively validate the proposed approach. A simulation of a compound fault is introduced in the following subsections.
5.1. Simulation
To verify the effectiveness of the method, outer and rolling element faults are simulated. The signal expression iswhere is the outer race fault signal, is rolling element fault signal, and are white noises.
Among them, the expression of isand the expression of iswhere is the amplitude of the impact response, is the damping coefficient, is the resonance frequency, is the number of pulses, , and .
Then, Hz and Hz in this simulation, where is outer race fault characteristics frequency and is rolling element fault characteristics fault. The relevant parameters are listed in Table 1, and the result of the simulation signal is presented in Figure 5.
(a)
(b)
(c)
With the basic model, the entire simulation process will begin as follows: first, a normal bearing vibration signal is constructed through simulation, and a control chart is developed using the algorithm described in the flowchart. Then, considering the long-term work of a bearing, our simulation designs a step to simulate the increased degree for faults. We increase gradually and use algorithms to build the statistics for monitoring the process. Thus, if the statistics exceed the control limit, then an alarm can be issued, resulting in the analysis of Subband I and Subband II through the envelope spectrum.
During the monitoring process in the actual simulation, each sample includes 4 subsamples and the subsamples contain 4096 data points. The calculation of the control limits is 9.62 through equation (8) and the calculation of the control limit is 16.88 through kernel estimation (which is equation (9)). As shown in Figure 6(b), the statistics begin to exceed the control limit at sample number 404 and the SNR is −10 db at this time, demonstrating that the robust control chart can indicate early faults. In accordance with the early warning criteria of the control chart, the small diagram can be seen in Figure 6(b). If four out of five consecutive points exceed the control limit, then an early warning should be issued [28]. Thus, an alarm will be issued, and fault diagnosis will be performed on the subband that calculated the dual-kurtogram algorithm at this period. For comparison, conventional control charts are also used for monitoring, as shown in Figure 6(a), and the sensitivity of conventional control charts to abnormal signals is found to be completely lower than that of the proposed method.
(a)
(b)
In Figure 7, the center frequency is 1500 Hz and the bandwidth is 1000 Hz; in Figure 8, the center frequency is 3500 Hz and the bandwidth is 1000 Hz in each subband. These values are consistent with our simulated resonance frequencies , which are 1500 Hz and 3500 Hz. Simultaneously, square envelope analysis should be performed on the subbands, and the characteristic fault frequencies, which are 110 Hz and 49 Hz, can be extracted as shown in Figures 9 and 10. These findings are also consistent with our simulation results.
(a)
(b)
(c)
(a)
(b)
(c)
For comparison, EMD, a traditional fault diagnosis method that can decompose the complicated signal function into a number of intrinsic mode functions (IMFs) is applied to the simulation signal. The decomposed first six IMFs are shown in Figure 11. The fast Fourier transform is adopted to the intrinsic mode functions (IMFs) to obtain the frequency spectrums, which are illustrated in Figure 12. The fault characteristic frequency cannot be clearly detected in Figure 12.
Envelope spectrum analysis draws the envelope of the signal and identifies the fault frequency from the envelope spectrum. The result of applying the Hilbert envelope spectrum analysis method is shown in Figure 13. The fault characteristic frequency also cannot be clearly detected in Figure 13.
Comparison studies with the EMD and Hilbert envelope spectrum analysis method show that the fault characteristics are mixed and difficult to distinguish in either Figure 12 or 13, also indicating a better performance of our proposed method.
5.2. Bearing’s Run-to-Failure Data Set
In this study, bearings’ run-to-failure data are used to test the effectiveness of the proposed algorithm in condition monitoring and fault diagnosis. Spherical roller bearing run-to-failure life data are obtained from a test conducted for 35 days [29]. All the bearings are force-lubricated. Rexnord ZA-2115 double row bearings are installed on the shaft as shown in Figure 14. The rotation speed is kept constant at 2000 RPM via an AC motor coupled to the shaft by rubber belts. A radial load of 2722 KG is applied to the shaft and bearing by a spring mechanism. PCB 353B33 accelerometers are used to acquire the bearing’s vibration signals. Vibration signals are recorded by a computer-based NI-DAQ at a sampling rate of 20 kHz for 1 s each time and repeated every 10 min.
An inner race defect was discovered in test bearing 3, and a roller element defect and outer race defect were found in test bearing 4. In this study, the data of bearing 4 was chosen to verify the feasibility of the proposed algorithm.
In Figure 15(a), each sample includes 20480 data points. The calculation of the control limit is 6.69 through equation (6) and the calculation of the control limit is 17.82 through kernel estimation, after sample 1666, five continuous points are above the control limit; thus, the alarm should be issued. Then, for the subband of sample 1666, the characteristic fault frequency can be calculated using the square envelope. For comparison, a conventional Hotelling’s control chart is shown in Figure 15(b). The data do not follow a normal distribution, and, thus, the control limit of the traditional method is estimated incorrectly and exceeded in a large number of normal working stages from sample 200 to sample 500.
(a)
(b)
In Figure 16, the center frequency is 8125 Hz and bandwidth is 1250 Hz; in Figure 17, the center frequency is 1458.3 Hz and bandwidth is 416.6 Hz in each subband. At the same time, the square envelope analysis should be performed on the subbands. From the results presented in Figures 18 and 19 there are outer race faults Hz and rolling element faults Hz are found in this dataset, and severe faults occur at sample number 1666 at approximately 30 days. These results are similar to the conclusion obtained by Qiu [30].
(a)
(b)
(c)
(a)
(b)
(c)
For comparison, the EMD and envelop spectrum analysis are applied to the data of sample 1666. The first four IMFs after EMD are shown in Figure 20. The frequency spectrums of IMFs are illustrated in Figure 21. The result of applying the envelope spectrum is shown in Figure 22. Comparison studies with the EMD and Hilbert envelope spectrum analysis method show that the fault characteristics are mixed and difficult to distinguish in either Figure 21 or 22, also indicating a better performance of our proposed method. In summary, an improved dual-kurtogram-based control chart can smoothly carry out an early warning of faults and can effectively diagnose compound faults under actual circumstances.
6. Conclusion
An improved dual-kurtogram-based control chart is proposed for the condition monitoring and compound fault diagnosis of rolling bearings in this study. The dual-kurtogram algorithm is used to construct Subband I and Subband II, which include high- and low-impact density faults. The SK extracted by Subband I and the ESK extracted by Subband II is used to construct a dual-kurtogram-based control chart to enhance the sensitivity of statistics. The control chart adopts KDE to guarantee a precise control limit. Data are monitored through the control chart. When statistics exceed the control limit, an alarm is issued, and square envelope analysis is performed on Subband I and Subband II to obtain the characteristic fault frequency. Although this method achieves effective condition monitoring and fault diagnosis, it still exhibits several shortcomings. The following situation will be the focus of our next research. The operating state of equipment is more complicated when experiencing practical problems, and static condition monitoring cannot meet actual requirements. Therefore, an adaptive control chart for dynamic monitoring should be established. Some fault subbands are submerged by the dual-kurtogram algorithm to a certain extent, and, thus, we must research further on the use of this algorithm to obtain hidden fault subbands and achieve complete fault diagnosis.
Data Availability
All the data used in the study are available. Previously reported data were used to support this study and are available. These prior studies are cited at relevant places within the text as references.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This research was supported by the National Natural Science Foundation of China (Grant no. 51905218) and the Senior Talent Foundation of Jiangsu University (Grant no. 19JDG040).