# An Improved Dual-Kurtogram-BasedT2 Control Chart for Condition Monitoring and Compound Fault Diagnosis of Rolling Bearings

Authors: Zhiyuan Jiao; Wei Fan; Zhenying Xu
Journal: Shock and Vibration
 (2021)

---

## 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 deriveT2 statistics based on the covariance determinant. An alarm will be triggered when the T2 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 T2 control chart, envelope analysis, and empirical mode decomposition.


---

## Body

## 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 T2 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’sT2 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-basedT2 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. Section2 summarizes Hotelling’s T2 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 T2 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’sT2 Control Chart

### 2.1. Conventional Hotelling’sT2 Control Chart

Hotelling’sT2 control chart [22] is one of the multivariate control charts that can be used for process monitoring. The T2 control chart is calculated for μ and Σ under p variables:μ is the average value; and Σ is the covariance matrix.Each of them batches of sample data is with a sample size of n and the mean value and sample variance of the jth variable test data in the kth batch of data are(1)x¯jk=1n∑xijk,sjk2=1n−1∑i=1nxijk−x¯ij2,where xijk represents the ith observation value of the jth variable in the kth batch of data.The covariance between thejth and hth variables in the kth batch of data is(2)sjhk=1n−1∑i=1nxijk−x¯jkxihk−x¯hk.Then, the statisticsx¯jk,Sjk2, and Sjhk2 are averaged for m batches of data and p variables:(3)x¯¯j=1m∑k=1mx¯jk,j=1,2,…,p,s¯jk2=1m∑k=1msjk2,j=1,2,…,p,s¯jh=1m∑k=1msjk,j≠h,x¯¯j are elements of the vector x¯¯ and the sample covariance matrix is(4)S=s¯12s¯12⋯s¯1ps¯22⋯s¯2p⋱⋮s¯p2.TheT2 statistics can be calculated as follows:(5)Ti2=nx¯i−x¯¯′S¯−1x¯i−x¯¯.Supposing the data follow a multivariate normal distribution, the control limit of Hotelling’sT2 can be obtained using the following equation:(6)UCL=pn+1n−1n2−npFαp,n−p,where Fαp,n−p is the distribution function, n is the number of observations, p is the number of variables, and α is the false alarm rate. The process is considered in control of the T2 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. ConstructingT2 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 constructT2 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 ofT2 statistics [25]. Let T2 be a Hotelling’s statistic obtained under in-control condition. The distribution of T2 statistics can be calculated using the following kernel function:(7)f^ht=1nh∑i=1nKt−Ti2h,where K and h 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(8)Ku=12πe−1/2u2,−∞<u<∞.Therefore, the control limit ofT2 statistics based on KDE can be estimated by taking the percentile of the kernel distribution. Hence, the control limit of T2 statistics based on KDE is equal to the 1001−α-th percentile of T2distribution which can be calculated as follows:(9)CLkernel=F^ht−11−α,where F^ht 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’sT2 control chart applying to condition monitoring.


## 2.1. Conventional Hotelling’sT2 Control Chart

Hotelling’sT2 control chart [22] is one of the multivariate control charts that can be used for process monitoring. The T2 control chart is calculated for μ and Σ under p variables:μ is the average value; and Σ is the covariance matrix.Each of them batches of sample data is with a sample size of n and the mean value and sample variance of the jth variable test data in the kth batch of data are(1)x¯jk=1n∑xijk,sjk2=1n−1∑i=1nxijk−x¯ij2,where xijk represents the ith observation value of the jth variable in the kth batch of data.The covariance between thejth and hth variables in the kth batch of data is(2)sjhk=1n−1∑i=1nxijk−x¯jkxihk−x¯hk.Then, the statisticsx¯jk,Sjk2, and Sjhk2 are averaged for m batches of data and p variables:(3)x¯¯j=1m∑k=1mx¯jk,j=1,2,…,p,s¯jk2=1m∑k=1msjk2,j=1,2,…,p,s¯jh=1m∑k=1msjk,j≠h,x¯¯j are elements of the vector x¯¯ and the sample covariance matrix is(4)S=s¯12s¯12⋯s¯1ps¯22⋯s¯2p⋱⋮s¯p2.TheT2 statistics can be calculated as follows:(5)Ti2=nx¯i−x¯¯′S¯−1x¯i−x¯¯.Supposing the data follow a multivariate normal distribution, the control limit of Hotelling’sT2 can be obtained using the following equation:(6)UCL=pn+1n−1n2−npFαp,n−p,where Fαp,n−p is the distribution function, n is the number of observations, p is the number of variables, and α is the false alarm rate. The process is considered in control of the T2 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. ConstructingT2 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 constructT2 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 ofT2 statistics [25]. Let T2 be a Hotelling’s statistic obtained under in-control condition. The distribution of T2 statistics can be calculated using the following kernel function:(7)f^ht=1nh∑i=1nKt−Ti2h,where K and h 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(8)Ku=12πe−1/2u2,−∞<u<∞.Therefore, the control limit ofT2 statistics based on KDE can be estimated by taking the percentile of the kernel distribution. Hence, the control limit of T2 statistics based on KDE is equal to the 1001−α-th percentile of T2distribution which can be calculated as follows:(9)CLkernel=F^ht−11−α,where F^ht 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’sT2 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 beXt, STFT can be expressed as(10)Xt,f=∫−∞+∞xτwt−τe−j2πfτdτ.Xt,f can be regarded as a complex envelope in terms of frequency and time, and SK based on STFT at frequency f can be expressed as(11)K4Xf=Xt,ft4Xt,ft22−2.Considering the presence of stationary added noise bn, the SK can be further expressed as(12)Kx+bf=Kxf1+ρf2,where Kx+bf is the SK of signal xn with stationary added noise bn and ρf 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 filterhn with a cut-off frequency of fc=1/8+ε, ε⩾0 is provided. A high-pass filter h0n and a low-pass filter h1n are constructed on basis of hn. The expression is(13)h0n=hn×ejπn/4,h1n=hn×ej3πn/4.The bandwidths of filtersh0nandh1n are 0,1/4Fs and 1/4Fs,1/2Fs, respectively, where Fs 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 k, and 2k subbands can be obtained for each layer. Let the ith subbands of the kth layer be Cki,i=0,1,2,…,2k−1. Its center frequency is fic=i+2−1×2−k−1×Fs, its bandwidth is △fk=2−k−1×Fs, and Cki’s adjacent lower subbands are(14)Ck+12i=h0n×Cki,Ck+12i+1=h1n×Cki.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, theith-level three-division filter is inserted into the ith and i+1th 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 ofh0n and h1n, the real and imaginary parts of subsignal Cki differ by 90° regardless of whether it is a high-pass or low-pass filter. By the principle of envelope demodulation, Cki is the envelope of the analytic signal, and, thus, the SK of the subsignal can be expressed as(15)KXf,△f=Cki4Cki22−2.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.Figure 1
Fast computation of the kurtogram using an arborescent filter-bank structure.Figure 2
Paving of the plane in case of 1/3-binary tree kurtogram estimator proposed by Antoni [16].


### 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 constructT2 statistics. The improved dual-kurtogram-based T2 control chart algorithm is used to monitor and provide early warning of weak faults.


## 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 beXt, STFT can be expressed as(10)Xt,f=∫−∞+∞xτwt−τe−j2πfτdτ.Xt,f can be regarded as a complex envelope in terms of frequency and time, and SK based on STFT at frequency f can be expressed as(11)K4Xf=Xt,ft4Xt,ft22−2.Considering the presence of stationary added noise bn, the SK can be further expressed as(12)Kx+bf=Kxf1+ρf2,where Kx+bf is the SK of signal xn with stationary added noise bn and ρf 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 filterhn with a cut-off frequency of fc=1/8+ε, ε⩾0 is provided. A high-pass filter h0n and a low-pass filter h1n are constructed on basis of hn. The expression is(13)h0n=hn×ejπn/4,h1n=hn×ej3πn/4.The bandwidths of filtersh0nandh1n are 0,1/4Fs and 1/4Fs,1/2Fs, respectively, where Fs 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 k, and 2k subbands can be obtained for each layer. Let the ith subbands of the kth layer be Cki,i=0,1,2,…,2k−1. Its center frequency is fic=i+2−1×2−k−1×Fs, its bandwidth is △fk=2−k−1×Fs, and Cki’s adjacent lower subbands are(14)Ck+12i=h0n×Cki,Ck+12i+1=h1n×Cki.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, theith-level three-division filter is inserted into the ith and i+1th 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 ofh0n and h1n, the real and imaginary parts of subsignal Cki differ by 90° regardless of whether it is a high-pass or low-pass filter. By the principle of envelope demodulation, Cki is the envelope of the analytic signal, and, thus, the SK of the subsignal can be expressed as(15)KXf,△f=Cki4Cki22−2.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.Figure 1
Fast computation of the kurtogram using an arborescent filter-bank structure.Figure 2
Paving of the plane in case of 1/3-binary tree kurtogram estimator proposed by Antoni [16].


## 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 constructT2 statistics. The improved dual-kurtogram-based T2 control chart algorithm is used to monitor and provide early warning of weak faults.


## 4. An Improved Dual-Kurtogram-BasedT2 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 Figure3.Figure 3
Variation trend of fault characteristic frequency with the normalized value of SK and ESK.In Figure3, 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 Figure4. 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, T2 statistics are constructed by a covariance matrix and used to calculate the control limit with KDE.Figure 4
Flowchart of the algorithm.In Phase II, the algorithm in Phase I is repeated to obtain SK and ESK, and the covariance matrix was calculated to constructT2 statistics. The T2 statistics are compared with the control limit. When the T2 statistic exceeds the control limit, the subbands used to construct the T2 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 is(16)st=s1t+s2t+s3t,where s1t is the outer race fault signal, s2t is rolling element fault signal, and s3t are white noises.Among them, the expression ofs1t is(17)s1t=∑n=0N1A1e−β1t−T1cos2πf1t−T1ut−T1,and the expression of s2t is(18)s2t=1+sin500πt⋅∑n=0N2A2e−β2t−T2cos2πf2t−T2ut−T2,where Ai is the amplitude of the impact response, βi is the damping coefficient, fi is the resonance frequency, Ni is the number of pulses, T1=fBPFO−1, and T2=fBSF−1.Then,fBPFO=111.47 Hz and fBSD=51.74 Hz in this simulation, where fBPFO is outer race fault characteristics frequency and fBSD 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.Table 1
The simulation parameters.
iAiβifi148001500235003500Figure 5
(a) The simulation signal of the faulty outer race. (b) The simulation signal of faulty rolling elements. (c) The mixed-signal includes (a) (b) and Gaussian noises.
(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 increaseAi gradually and use algorithms to build the T2 statistics for monitoring the process. Thus, if the T2 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 T2 statistics begin to exceed the control limit at sample number 404 and the SNR is −10 db at this time, demonstrating that the robust T2 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.Figure 6
(a) Conventional Hotelling’sT2 control chart. The redline is UCL of 9.62. (b) An improved dual-kurtogram-based T2 control chart. The redline is UCL of 16.88.
(a)(b)In Figure7, 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 fi, 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.Figure 7
The colormap of the simulation signal via fast kurtogram.Figure 8
The colormap of the simulation signal via a fast computation of protrugram.Figure 9
(a) The original simulation signal. (b) The signal filtered by fast kurtogram. (c) The square envelope for the filtered.
(a)(b)(c)Figure 10
(a) The original simulation signal. (b) The signal filtered by fast computation of protrugram. (c) The square envelope for the filtered signal.
(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 Figure11. 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.Figure 11
The first four IMF components after EMD.Figure 12
Frequency spectrums of IMFs after EMD.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 Figure13. The fault characteristic frequency also cannot be clearly detected in Figure 13.Figure 13
Envelope spectrum analysis for the vibration signal of sample 404.Comparison studies with the EMD and Hilbert envelope spectrum analysis method show that the fault characteristics are mixed and difficult to distinguish in either Figure12 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.Figure 14
IMS test rig [30].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 Figure15(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 T2 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.Figure 15
(a) An improved dual-kurtogram-basedT2 control chart. The redline is UCL of 17.82. (b) Conventional Hotelling’s T2 control chart. The redline is UCL of 6.69.
(a)(b)In Figure16, 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 fBPFO=238.4 Hz and rolling element faults fBSF=119.3 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].Figure 16
The colormap of data via fast kurtogram.Figure 17
The colormap of data via fast computation of protrugram.Figure 18
(a) The original signal. (b) The signal filtered by fast kurtogram. (c) The square envelope for the filtered signal.
(a)(b)(c)Figure 19
(a) The original signal. (b) The signal filtered by fast computation of protrugram. (c) The square envelope for the filtered signal.
(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 Figure20. 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.Figure 20
The first four IMF components after EMD.Figure 21
Frequency spectrums of IMFs after EMD.Figure 22
Envelope spectrum analysis for data of sample 1666.


## 5.1. Simulation

To verify the effectiveness of the method, outer and rolling element faults are simulated. The signal expression is(16)st=s1t+s2t+s3t,where s1t is the outer race fault signal, s2t is rolling element fault signal, and s3t are white noises.Among them, the expression ofs1t is(17)s1t=∑n=0N1A1e−β1t−T1cos2πf1t−T1ut−T1,and the expression of s2t is(18)s2t=1+sin500πt⋅∑n=0N2A2e−β2t−T2cos2πf2t−T2ut−T2,where Ai is the amplitude of the impact response, βi is the damping coefficient, fi is the resonance frequency, Ni is the number of pulses, T1=fBPFO−1, and T2=fBSF−1.Then,fBPFO=111.47 Hz and fBSD=51.74 Hz in this simulation, where fBPFO is outer race fault characteristics frequency and fBSD 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.Table 1
The simulation parameters.
iAiβifi148001500235003500Figure 5
(a) The simulation signal of the faulty outer race. (b) The simulation signal of faulty rolling elements. (c) The mixed-signal includes (a) (b) and Gaussian noises.
(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 increaseAi gradually and use algorithms to build the T2 statistics for monitoring the process. Thus, if the T2 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 T2 statistics begin to exceed the control limit at sample number 404 and the SNR is −10 db at this time, demonstrating that the robust T2 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.Figure 6
(a) Conventional Hotelling’sT2 control chart. The redline is UCL of 9.62. (b) An improved dual-kurtogram-based T2 control chart. The redline is UCL of 16.88.
(a)(b)In Figure7, 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 fi, 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.Figure 7
The colormap of the simulation signal via fast kurtogram.Figure 8
The colormap of the simulation signal via a fast computation of protrugram.Figure 9
(a) The original simulation signal. (b) The signal filtered by fast kurtogram. (c) The square envelope for the filtered.
(a)(b)(c)Figure 10
(a) The original simulation signal. (b) The signal filtered by fast computation of protrugram. (c) The square envelope for the filtered signal.
(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 Figure11. 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.Figure 11
The first four IMF components after EMD.Figure 12
Frequency spectrums of IMFs after EMD.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 Figure13. The fault characteristic frequency also cannot be clearly detected in Figure 13.Figure 13
Envelope spectrum analysis for the vibration signal of sample 404.Comparison studies with the EMD and Hilbert envelope spectrum analysis method show that the fault characteristics are mixed and difficult to distinguish in either Figure12 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.Figure 14
IMS test rig [30].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 Figure15(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 T2 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.Figure 15
(a) An improved dual-kurtogram-basedT2 control chart. The redline is UCL of 17.82. (b) Conventional Hotelling’s T2 control chart. The redline is UCL of 6.69.
(a)(b)In Figure16, 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 fBPFO=238.4 Hz and rolling element faults fBSF=119.3 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].Figure 16
The colormap of data via fast kurtogram.Figure 17
The colormap of data via fast computation of protrugram.Figure 18
(a) The original signal. (b) The signal filtered by fast kurtogram. (c) The square envelope for the filtered signal.
(a)(b)(c)Figure 19
(a) The original signal. (b) The signal filtered by fast computation of protrugram. (c) The square envelope for the filtered signal.
(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 Figure20. 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.Figure 20
The first four IMF components after EMD.Figure 21
Frequency spectrums of IMFs after EMD.Figure 22
Envelope spectrum analysis for data of sample 1666.


## 6. Conclusion

An improved dual-kurtogram-basedT2 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 T2 control chart to enhance the sensitivity of T2 statistics. The T2 control chart adopts KDE to guarantee a precise control limit. Data are monitored through the control chart. When T2statistics 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.



---

*Source: 6649125-2021-01-13.xml*