# Rolling Bearing Fault Diagnosis Based on CEEMD and Time Series Modeling

Authors: Liye Zhao; Wei Yu; Ruqiang Yan
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

Accurately identifying faults in rolling bearing systems by analyzing vibration signals, which are often nonstationary, is challenging. To address this issue, a new approach based on complementary ensemble empirical mode decomposition (CEEMD) and time series modeling is proposed in this paper. This approach seeks to identify faults appearing in a rolling bearing system using proper autoregressive (AR) model established from the nonstationary vibration signal. First, vibration signals measured from a rolling bearing test system with different defect conditions are decomposed into a set of intrinsic mode functions (IMFs) by means of the CEEMD method. Second, vibration signals are filtered with calculated filtering parameters. Third, the IMF which is closely correlated to the filtered signal is selected according to the correlation coefficient between the filtered signal and each IMF, and then the AR model of the selected IMF is established. Subsequently, the AR model parameters are considered as the input feature vectors, and the hidden Markov model (HMM) is used to identify the fault pattern of a rolling bearing. Experimental study performed on a bearing test system has shown that the presented approach can accurately identify faults in rolling bearings.


---

## Body

## 1. Introduction

Rolling element bearing failure is one of the foremost causes of failures in rotating machinery, and such failure may result in costly production loss and catastrophic accidents. Early detection and diagnosis of bearing faults while the machine is still in operation can help to avoid abnormal event progression and to reduce productivity loss [1]. Since structural defects can cause changes of the bearing dynamic characteristics as manifested in vibrations, vibration-based analysis has long been established as a commonly used technique for diagnosing bearing faults [2]. However, some nonlinear factors such as clearance, friction, and stiffness affect complexity of the vibration signals; thus it is difficult to make an accurate evaluation on the working condition of rolling bearings only through analysis in time or frequency domain as it does traditionally [3].In order to overcome limitations of the traditional techniques, autoregressive (AR) model has been successfully applied to extracting features from vibration signals for fault diagnosis in recent years [4–6]. This is because AR model is a time series analysis method whose parameters comprise important information of the system condition, and an accurate AR model can reflect the characteristics of a dynamic system [7]. For example, AR model was combined with a fuzzy classifier for fault diagnosis in vehicle transmission gear [8]. Three distinct techniques of autoregressive modeling were compared for their performance and reliability under conditions of various bearings signal lengths [9]. A diagnosis method based on the AR model and continuous HMM has also been used to monitor and diagnose the rolling bearing working conditions [10]. However, when the AR model is applied directly to the nonstationary bearing vibration signals, the analysis results are imperfect since the estimation method of the autoregression parameters of the AR model is no longer applicable. Because the vibration signal is nonstationary, whereas the AR model is suitable for stationary signal processing, it is, therefore, necessary to preprocess the vibration signals before the AR model is generated.Empirical mode decomposition (EMD) is an adaptive time-frequency signal processing method [11]. With EMD, a signal is decomposed into a series of intrinsic mode functions (IMFs) according to its own characteristics [12]. For example, a new fault feature extraction approach based on EMD method and AR model was used to process vibration signals of roller bearings [3]. However, when the EMD method is applied to the nonstationary signals containing intermittent signal components, the original signal cannot be decomposed accurately because of the problem of mode mixing [13]. To alleviate mode mixing, Wu and Huang developed ensemble empirical mode decomposition (EEMD) to improve EMD. By adding noise to the original signal and calculating the means of IMFs repeatedly, EEMD is more accurate and effective for signal decomposition [13]. Although the EEMD method has effectively resolved the mode-mixing problem, it is time consuming for implementing the large enough ensemble mean. That is to say, the algorithm efficiency will be greatly reduced. Aiming at solving this problem, the complementary ensemble EMD (CEEMD) method is proposed [14]. In this approach, the residue of added white noises can be extracted from the mixtures of data and white noises via pairs of complementary ensemble IMFs with positive and negative added white noises. The CEEMD method has the same performance as the EEMD, but the computational efficiency is greatly improved.In this paper, we combine the advantages of CEEMD and time series model and propose a new method based on CEEMD and AR model for rolling bearing fault diagnosis. The CEEMD is used as the pretreatment to filter the signal and extract the IMF which is closely correlated to the filtered signal, and then the AR model of the selected IMF is established. The AR model parameters are used as the feature vectors to a classifier, where the hidden Markov model (HMM) is used to identify the fault pattern of a rolling bearing. The rest of this paper is organized as follows. In Section2, the review of the fault diagnosis method based on AR model is presented, and the proposed method for rolling bearing fault diagnosis is discussed. The evaluations and experiments are presented in Section 3. Finally, concluding remarks are drawn in Section 4.


## 2. Theoretical Framework

### 2.1. Time Series Modeling

Autoregressive moving average (ARMA) model is the representative time series model, which can be expressed in linear difference equation form as(1)
x
t
+
φ
1
x
t
-
1
+
⋯
+
φ
n
x
t
-
n
=
a
t
+
θ
1
a
t
-
1
+
⋯
+
θ
m
a
t
-
m
,
where n and m are the parameters of the ARMA (n
,
m) model, x
t is zero mean stationary random sequence, a
t is white noise sequence, and φ
i and θ
j are model parameters to be estimated. The parameters of φ
i and θ
j are estimated by the time sequence of x
t (t
=
1,2
,
3
…), which is called the time series modeling. If φ
i
=
0, the ARMA (n
,
m) model will degrade to m order MA(m) model, and if θ
i
=
0, the ARMA (n
,
m) model will degrade as n order AR (n) model in (1). The AR model is stable and its structure is simpler than ARMA model. Therefore, the AR model will be established for characterizing the rolling bearing vibration signal, if the precision of the model is enough for expressing the system, which is expressed as
(2)
x
t
=
φ
1
x
t
-
1
+
⋯
+
φ
n
x
t
-
n
+
a
t
,
where t
=
1,2
,
…
,
N, N is the length of the time series x
t, n is the order number, and a
t
~
NID
(
0
,
σ
a
2
). The σ
a
2 is expressed as
(3)
σ
α
2
=
1
N
-
n
n
∑
t
=
n
+
1
n
(
x
t
-
∑
ι
=
1
n
φ
ι
x
t
-
1
)
2
.It is critical to determine the order number of the AR model, because the accuracy of the order not only affects the accuracy of identification of the system, but also influences the stability of the system. In order to estimate the order of the AR model correctly, FPE criterion, BIC criterion, and AIC criterion are usually used [15], and they are expressed asFPE criterion(4)
FPE
(
n
)
=
N
+
n
N
-
n
σ
a
2
,
BIC criterion(5)
AIC
(
n
)
=
N
ln
⁡
σ
a
2
+
2
n
,
AIC criterion(6)
BIC
(
n
)
=
N
ln
⁡
σ
a
2
+
n
ln
⁡
N
.After the model order is determined, the nonlinear least squares method can be used to estimate model parameters, and then the AR model with specific parameters is established.


### 2.2. Complementary Ensemble Empirical Mode Decomposition

Complementary ensemble empirical mode decomposition (CEEMD) is an improved algorithm based on empirical mode decomposition (EMD). Through EMD process, any complex time series can be decomposed into finite numbers of intrinsic mode functions (IMFs), and each IMF reflects the dynamic characteristic of the original signal. The IMF component must satisfy two conditions: (a) the number of poles and zeros is either equal to each other or differs at most by one; (b) the upper and lower envelopes must be locally symmetric about the timeline. The basic principle of EMD method is to decompose the original signalx
(
t
) into the form as shown in (7) by continuously eliminating the mean of the upper and lower envelope connected with the minimum and maximum of the signal [16]. Consider
(7)
x
t
=
∑
i
=
1
n
im
f
i
(
t
)
+
r
n
(
t
)
,
where x
t is the vibration signal, im
f
i
(
t
) is the IMF component including different frequency bands ranging from high to low, and r
n
(
t
) is the residue of the decomposition process, which is the mean trend of x
t.The EMD method is a kind of adaptive local analysis method, with each IMF highlighting the local features of the data. However, EMD decomposition results often suffer from mode mixing, which is defined as either a single IMF consisting of widely disparate scales or a signal residing in different IMF components [17]. To make it clear, a simulated signal s
(
t
) consists of a Gaussian-type impulse interference s
1
(
t
) and a cosine component with 500 Hz frequency s
2
(
t
), and a trend term s
3
(
t
) is used as an example. The equation of the simulated signal is expressed as
(8)
s
(
t
)
=
sin
(
2
π
α
t
)
e
-
(
(
t
-
t
0
)
2
/
σ
)
+
cos
⁡
(
2
π
β
t
)
+
50
t
,
where α
=
3000, β
=
500, and σ
=
10
6.The waveform of the simulated signal is shown in Figure1, and the corresponding EMD results for the signal s
(
t
) are shown in Figure 2, where the mode mixing happens.Figure 1
Signal waveforms.Figure 2
The decomposition result by EMD.To overcome the problem of mode mixing, the ensemble empirical mode decomposition (EEMD) was proposed [18], where Gaussian white noises with finite amplitude are added to the original signal during the entire decomposition process. Due to the uniform distribution statistical characteristics of the white noise, the signal with white noise becomes continuous in different time scales, and no missing scales are present. As a result, mode mixing is effectively eliminated by the EEMD process [18]. The EEMD decomposition result of signal s
(
t
) is shown in Figure 3, where the added white noise amplitude is 0.25 times the original signal standard deviation, and the number of decompositions is 200 times.Figure 3
The decomposition result by EEMD.It should be noted that, during the EEMD process, each individual trial may produce noisy results, but the effect of the added noise can be suppressed by large number of ensemble mean computations. This would be too time consuming to implement. An improved algorithm, named complementary ensemble mode decomposition (CEEMD), is suggested to improve the computation efficiency. In this algorithm, the residue of the added white noises can be extracted from the mixtures of data and white noises via pairs of complementary ensemble IMFs with positive and negative added white noises. Although this new approach yields IMF with a similar RMS noise to EEMD, it eliminates residue noise in the IMFs and overcomes the problem of mode mixing with much more efficiency [14]. The procedure on implementing CEEMD is shown below:(a)
x
1 and x
2 are constructed by adding a pair of opposite phase Gaussian white noises x
n with the same amplitude. Then x
1
=
x
+
x
n and x
2
=
x
-
x
n;
(b)
x
1 and x
2 are decomposed by EMD only a few times, and IMFx1 and IMFx2 are ensemble means of the corresponding IMF generated from each trial;
(c)
the average of corresponding component inIM
F
x
1 and IM
F
x
2 is calculated as the CEEMD decomposition results; that is,
(9)
IMF
=
(
IM
F
x
1
+
IM
F
x
2
)
2
.The flow chart of CEEMD is shown in Figure 4, where n is the decomposition trials.Figure 4
Decomposition flow chart of CEEMD.Figure5 is the decomposition result by CEEMD for the signal s
(
t
). As compared to the result shown in Figure 3, the decomposition accuracies of EEMD and CEEMD are consistent, while EEMD takes 1.62 s and CEEMD only needs 0.13 s.Figure 5
The decomposition result by CEEMD.


### 2.3. Fault Diagnosis Based on CEEMD and Time Series Model

Based on CEEMD and time series model, a hybrid fault diagnosis approach can be designed. The hybrid approach combines the advantages of CEEMD method in the nonstationary signal decomposition with the ability of time series modeling in feature extraction. The flow chart of the developed approach is shown in Figure6.Figure 6
The flow chart of the proposed method.The main steps are as follows.Step  1. The rolling bearing vibration signal is sampled and then decomposed by CEEMD with the process shown in Figure 4.Step  2. The product of energy density and average period of the IMFs which is a constant value according to [19] is calculated using (10) and parameter R
P
j is calculated using (11). Then the signal is filtered by comparing the parameter R
P
j and the given threshold value; that is to say, when R
P
j
⩾
1, the previous j
-
1 IMFs with the trend term need to be removed as noise and to rebuild the residual IMFs as filtered signal [19, 20]:
(10)
P
j
=
E
j
×
T
j
,
(11)
R
P
j
=
|
P
j
-
(
(
1
/
(
j
-
1
)
)
∑
i
=
1
j
-
1
P
j
)
(
1
/
(
j
-
1
)
)
∑
i
=
1
j
-
1
P
j
|
,
(
j
≥
2
)
,
where E
j
=
(
1
/
N
)
∑
i
=
1
N
[
A
j
(
i
)
]
2 is the energy density of the j
th IMF, T
j
=
2
N
/
O
j is the average period of the j
th IMF, N is the length of each IMF, A
j is the amplitude of the j
th IMF, and O
j is the total number of extreme points of j
th IMF.Step  3. Equation (12) is used to calculate the correlation coefficient between the filtered signal and each IMF, and the IMF which is closely correlated to the filtered signal is selected for AR modeling [21]:
(12)
ρ
x
y
=
∑
k
=
1
N
x
(
k
)
y
(
k
)
[
∑
k
=
1
N
x
(
k
)
2
∑
k
=
1
N
y
(
k
)
2
]
1
/
2
.Step  4. The least square method is used to estimate the parameters vectors of the AR model established in Step3, and the parameters vectors are considered as the model feature vector.Step  5. After scalar quantization by index calculation formula of Lloyds algorithm in (13) [22], the feature vector is used to train the HMM of each bearing working condition:
(13)
indx
(
x
)
=
{
1
x
≤
partition
(
i
)
i
+
1
partition
(
i
)
<
x
≤
partition
(
i
+
1
)
N
partition
(
N
-
1
)
<
x
,
where N is the length of the codebook vector, partition (i) is the partition vector with the length of N
-
1, and x is the feature vector for scalar quantization.Step  6. A test vibration signal can then be acquired for diagnosis, and the model feature vector is first extracted. After scalar quantization, the feature vector is put into the well-trained HMMs, and the corresponding HMM which has the maximum probability is regarded as the classification result [23].


## 2.1. Time Series Modeling

Autoregressive moving average (ARMA) model is the representative time series model, which can be expressed in linear difference equation form as(1)
x
t
+
φ
1
x
t
-
1
+
⋯
+
φ
n
x
t
-
n
=
a
t
+
θ
1
a
t
-
1
+
⋯
+
θ
m
a
t
-
m
,
where n and m are the parameters of the ARMA (n
,
m) model, x
t is zero mean stationary random sequence, a
t is white noise sequence, and φ
i and θ
j are model parameters to be estimated. The parameters of φ
i and θ
j are estimated by the time sequence of x
t (t
=
1,2
,
3
…), which is called the time series modeling. If φ
i
=
0, the ARMA (n
,
m) model will degrade to m order MA(m) model, and if θ
i
=
0, the ARMA (n
,
m) model will degrade as n order AR (n) model in (1). The AR model is stable and its structure is simpler than ARMA model. Therefore, the AR model will be established for characterizing the rolling bearing vibration signal, if the precision of the model is enough for expressing the system, which is expressed as
(2)
x
t
=
φ
1
x
t
-
1
+
⋯
+
φ
n
x
t
-
n
+
a
t
,
where t
=
1,2
,
…
,
N, N is the length of the time series x
t, n is the order number, and a
t
~
NID
(
0
,
σ
a
2
). The σ
a
2 is expressed as
(3)
σ
α
2
=
1
N
-
n
n
∑
t
=
n
+
1
n
(
x
t
-
∑
ι
=
1
n
φ
ι
x
t
-
1
)
2
.It is critical to determine the order number of the AR model, because the accuracy of the order not only affects the accuracy of identification of the system, but also influences the stability of the system. In order to estimate the order of the AR model correctly, FPE criterion, BIC criterion, and AIC criterion are usually used [15], and they are expressed asFPE criterion(4)
FPE
(
n
)
=
N
+
n
N
-
n
σ
a
2
,
BIC criterion(5)
AIC
(
n
)
=
N
ln
⁡
σ
a
2
+
2
n
,
AIC criterion(6)
BIC
(
n
)
=
N
ln
⁡
σ
a
2
+
n
ln
⁡
N
.After the model order is determined, the nonlinear least squares method can be used to estimate model parameters, and then the AR model with specific parameters is established.


## 2.2. Complementary Ensemble Empirical Mode Decomposition

Complementary ensemble empirical mode decomposition (CEEMD) is an improved algorithm based on empirical mode decomposition (EMD). Through EMD process, any complex time series can be decomposed into finite numbers of intrinsic mode functions (IMFs), and each IMF reflects the dynamic characteristic of the original signal. The IMF component must satisfy two conditions: (a) the number of poles and zeros is either equal to each other or differs at most by one; (b) the upper and lower envelopes must be locally symmetric about the timeline. The basic principle of EMD method is to decompose the original signalx
(
t
) into the form as shown in (7) by continuously eliminating the mean of the upper and lower envelope connected with the minimum and maximum of the signal [16]. Consider
(7)
x
t
=
∑
i
=
1
n
im
f
i
(
t
)
+
r
n
(
t
)
,
where x
t is the vibration signal, im
f
i
(
t
) is the IMF component including different frequency bands ranging from high to low, and r
n
(
t
) is the residue of the decomposition process, which is the mean trend of x
t.The EMD method is a kind of adaptive local analysis method, with each IMF highlighting the local features of the data. However, EMD decomposition results often suffer from mode mixing, which is defined as either a single IMF consisting of widely disparate scales or a signal residing in different IMF components [17]. To make it clear, a simulated signal s
(
t
) consists of a Gaussian-type impulse interference s
1
(
t
) and a cosine component with 500 Hz frequency s
2
(
t
), and a trend term s
3
(
t
) is used as an example. The equation of the simulated signal is expressed as
(8)
s
(
t
)
=
sin
(
2
π
α
t
)
e
-
(
(
t
-
t
0
)
2
/
σ
)
+
cos
⁡
(
2
π
β
t
)
+
50
t
,
where α
=
3000, β
=
500, and σ
=
10
6.The waveform of the simulated signal is shown in Figure1, and the corresponding EMD results for the signal s
(
t
) are shown in Figure 2, where the mode mixing happens.Figure 1
Signal waveforms.Figure 2
The decomposition result by EMD.To overcome the problem of mode mixing, the ensemble empirical mode decomposition (EEMD) was proposed [18], where Gaussian white noises with finite amplitude are added to the original signal during the entire decomposition process. Due to the uniform distribution statistical characteristics of the white noise, the signal with white noise becomes continuous in different time scales, and no missing scales are present. As a result, mode mixing is effectively eliminated by the EEMD process [18]. The EEMD decomposition result of signal s
(
t
) is shown in Figure 3, where the added white noise amplitude is 0.25 times the original signal standard deviation, and the number of decompositions is 200 times.Figure 3
The decomposition result by EEMD.It should be noted that, during the EEMD process, each individual trial may produce noisy results, but the effect of the added noise can be suppressed by large number of ensemble mean computations. This would be too time consuming to implement. An improved algorithm, named complementary ensemble mode decomposition (CEEMD), is suggested to improve the computation efficiency. In this algorithm, the residue of the added white noises can be extracted from the mixtures of data and white noises via pairs of complementary ensemble IMFs with positive and negative added white noises. Although this new approach yields IMF with a similar RMS noise to EEMD, it eliminates residue noise in the IMFs and overcomes the problem of mode mixing with much more efficiency [14]. The procedure on implementing CEEMD is shown below:(a)
x
1 and x
2 are constructed by adding a pair of opposite phase Gaussian white noises x
n with the same amplitude. Then x
1
=
x
+
x
n and x
2
=
x
-
x
n;
(b)
x
1 and x
2 are decomposed by EMD only a few times, and IMFx1 and IMFx2 are ensemble means of the corresponding IMF generated from each trial;
(c)
the average of corresponding component inIM
F
x
1 and IM
F
x
2 is calculated as the CEEMD decomposition results; that is,
(9)
IMF
=
(
IM
F
x
1
+
IM
F
x
2
)
2
.The flow chart of CEEMD is shown in Figure 4, where n is the decomposition trials.Figure 4
Decomposition flow chart of CEEMD.Figure5 is the decomposition result by CEEMD for the signal s
(
t
). As compared to the result shown in Figure 3, the decomposition accuracies of EEMD and CEEMD are consistent, while EEMD takes 1.62 s and CEEMD only needs 0.13 s.Figure 5
The decomposition result by CEEMD.


## 2.3. Fault Diagnosis Based on CEEMD and Time Series Model

Based on CEEMD and time series model, a hybrid fault diagnosis approach can be designed. The hybrid approach combines the advantages of CEEMD method in the nonstationary signal decomposition with the ability of time series modeling in feature extraction. The flow chart of the developed approach is shown in Figure6.Figure 6
The flow chart of the proposed method.The main steps are as follows.Step  1. The rolling bearing vibration signal is sampled and then decomposed by CEEMD with the process shown in Figure 4.Step  2. The product of energy density and average period of the IMFs which is a constant value according to [19] is calculated using (10) and parameter R
P
j is calculated using (11). Then the signal is filtered by comparing the parameter R
P
j and the given threshold value; that is to say, when R
P
j
⩾
1, the previous j
-
1 IMFs with the trend term need to be removed as noise and to rebuild the residual IMFs as filtered signal [19, 20]:
(10)
P
j
=
E
j
×
T
j
,
(11)
R
P
j
=
|
P
j
-
(
(
1
/
(
j
-
1
)
)
∑
i
=
1
j
-
1
P
j
)
(
1
/
(
j
-
1
)
)
∑
i
=
1
j
-
1
P
j
|
,
(
j
≥
2
)
,
where E
j
=
(
1
/
N
)
∑
i
=
1
N
[
A
j
(
i
)
]
2 is the energy density of the j
th IMF, T
j
=
2
N
/
O
j is the average period of the j
th IMF, N is the length of each IMF, A
j is the amplitude of the j
th IMF, and O
j is the total number of extreme points of j
th IMF.Step  3. Equation (12) is used to calculate the correlation coefficient between the filtered signal and each IMF, and the IMF which is closely correlated to the filtered signal is selected for AR modeling [21]:
(12)
ρ
x
y
=
∑
k
=
1
N
x
(
k
)
y
(
k
)
[
∑
k
=
1
N
x
(
k
)
2
∑
k
=
1
N
y
(
k
)
2
]
1
/
2
.Step  4. The least square method is used to estimate the parameters vectors of the AR model established in Step3, and the parameters vectors are considered as the model feature vector.Step  5. After scalar quantization by index calculation formula of Lloyds algorithm in (13) [22], the feature vector is used to train the HMM of each bearing working condition:
(13)
indx
(
x
)
=
{
1
x
≤
partition
(
i
)
i
+
1
partition
(
i
)
<
x
≤
partition
(
i
+
1
)
N
partition
(
N
-
1
)
<
x
,
where N is the length of the codebook vector, partition (i) is the partition vector with the length of N
-
1, and x is the feature vector for scalar quantization.Step  6. A test vibration signal can then be acquired for diagnosis, and the model feature vector is first extracted. After scalar quantization, the feature vector is put into the well-trained HMMs, and the corresponding HMM which has the maximum probability is regarded as the classification result [23].


## 3. Evaluation of the Method Based on CEEMD and AR Model

### 3.1. Evaluation Using Simulated Signals

To demonstrate the validity of the method proposed in this study, three signalsx
1
(
t
), x
2
(
t
), and x
3
(
t
) are simulated as shown in Figure 7. The signal x
1
(
t
) consists of a Gaussian-type impulse interference, a cosine component with 10 Hz frequency, a trend term, and white noise. The signal x
2
(
t
) consists of a Gaussian-type impulse interference, a square wave with 65% duty ratio, a trend term, and white noise. The signal x
3
(
t
) consists of a Gaussian-type impulse interference, a sawtooth wave with 15 Hz frequency, a trend term, and white noise.Figure 7
Signal waveforms ofx
1
(
t
), x
2
(
t
), and x
3
(
t
).Figure8 shows the results of the CEEMD of signals x
1
(
t
), x
2
(
t
), and x
3
(
t
). Correlation coefficients between filtered signal and each IMF are illustrated in Table 1.Table 1
Correlation coefficients between filtered signal and each IMF.
Signal
Correlation coefficient
IMF1
IMF2
IMF3
IMF4
IMF5
IMF6
IMF7
IMF8
x
1
(
t
)
−0.0031
−0.0009
0.0371
0.4096
0.9668
0.2428
0.1273
−0.0448
x
2
(
t
)
0.0051
0.0004
0.0435
0.2111
0.4695
0.8887
0.7214
−0.0201
x
3
(
t
)
−0.0234
−0.0154
0.0286
0.5900
0.8953
0.1649
0.1887
−0.0210The decomposition results by CEEMD.
(a)
x
1
(
t
)
(b)
x
2
(
t
)
(c)
x
3
(
t
)It can be seen in Table1 that the IMF which is closely correlated to the filtered signal is IMF5 for both signal x
1
(
t
) and signal x
3
(
t
) and IMF6 for signal x
2
(
t
). They are used to construct the AR models, and the corresponding feature vectors are estimated as shown in Table 2. After scalar quantization, the feature vectors are used to train the HMM for signal classification.Table 2
Model parameter estimation results.
Signal
Model parameter
φ
1
φ
2
φ
3
φ
4
φ
5
φ
6
x
1
(
t
)
4.7183
−9.1103
9.2034
−5.1408
1.5207
−0.1914
x
2
(
t
)
4.8894
−9.8269
10.3945
−6.1194
1.9153
−0.2531
x
3
(
t
)
4.8718
−9.9955
11.1616
−7.2529
2.6430
−0.4282A total of 90 feature vectors were collected from three groups of signals using the proposed approach. One-third of the feature vectors in each condition were used for training the classifier and others were used for testing. The results of the signal classification are listed in Table3.Table 3
Signal classification results.
Signal type
Test sample
Classification results
Classification rate [%]
Overall classification rate [%]
x
1
(
t
)
x
2
(
t
)
x
3
(
t
)
x
1
(
t
)
20
19
1
0
95
96.7
x
2
(
t
)
20
0
19
1
95
x
3
(
t
)
20
0
0
20
100Results in Table3 indicate that the presented method based on CEEMD and time series modeling can effectively identify different signals, and the overall classification rate is 96.7%. For the purpose of comparison, the signal classification rates use the method based on time series modeling only, and the method based on EMD and time series modeling is also calculated. 88.3% and 93.3% classification rates are obtained, respectively. It is obvious that efficiency of the signal classification method proposed in this paper is improved to a certain extent.


### 3.2. Evaluation Using Experimental Data

In order to illustrate the practicability and effectiveness of the proposed method, a bearing fault data set from the electrical engineering laboratory of Case Western Reserve University is analyzed [24]. The data set is acquired from the test stand shown in Figure 9, where it consists of a 2 hp motor, a torque transducer, a dynamometer, and control electronics. The test bearings support the motor shaft which is the deep groove ball bearings with the type of 6205-2RS JEMSKF. Vibration data was collected at 12,000 samples per second using accelerometers, which are attached to the housing with magnetic bases. The motor load level was controlled by the fan in the right side of Figure 9.Figure 9
Bearing test stand.Figure10 illustrates representative waveforms of the sample vibration signals measured from the test bearings under four initial conditions: (a) signal from a healthy bearing, (b) signal from a bearing with inner ring defect, (c) signal from a bearing with rolling element defect, and (d) signal from a bearing with outer ring defect. These signals were measured under 0 hp motor load with the motor speed of 1797 rpm. The decomposed IMFs of these signals are shown in Figure 11.Vibration signal waveforms of different conditions.
(a)
(b)
(c)
(d)The decomposition results by CEEMD under different conditions.
(a)
No defect
(b)
Inner ring defect
(c)
Rolling element defect
(d)
Outer ring defectCorrelation coefficients calculated between the filtered signal and each IMF are shown in Table4.Table 4
Correlation coefficients between filtered signals and each IMF.
Signal
Correlation coefficient
IMF1
IMF2
IMF3
IMF4
IMF5
IMF6
IMF7
IMF8
(a)
0.4135
0.7538
0.4381
0.4880
0.4356
0.1792
0.0971
−0.0056
(b)
0.8794
0.4275
0.2583
0.1337
0.0421
0.0285
−0.0009
−0.0074
(c)
0.9509
0.2180
0.2325
0.1337
0.0821
0.0350
−0.0017
0.0009
(d)
0.9878
0.1267
0.0636
0.0509
0.0136
0.0060
−0.0008
−0.0068The IMF which is closely correlated to the filtered signal is IMF2 for signal (a) and IMF1 for signals (b), (c), and (d), respectively. These IMFs are used for AR model construction. The model order estimation curves of the four conditions based on the principle of FPE criterion are shown in Figure12. We can see that when the model order is 6, each model's residual tends to be stable. Therefore the model order is selected as 6, and the results of parameters estimation are listed in Table 5.Table 5
Model parameter estimation results.
Signal
Model parameter
φ
1
φ
2
φ
3
φ
4
φ
5
φ
6
(a)
3.1280
−4.7797
4.2245
−2.1489
0.4241
0.0356
(b)
0.2084
−1.3585
0.5142
−0.6356
0.3471
−0.0422
(c)
0.1335
−1.6472
0.3941
−0.8473
0.2142
−0.1011
(d)
−0.1172
−1.2159
0.1178
−0.1283
0.1467
0.2533The model order estimation curves.
(a)
No defect
(b)
Inner ring defect
(c)
Rolling element defect
(d)
Outer ring defectThe parameters in Table5 were quantified by Lloyds algorithm in (12) as feature vectors for training the HMMs of different conditions. The results of quantization are revealed in Figure 13.The results of quantization.
(a)
No defect
(b)
Inner ring defect
(c)
Rolling element defect
(d)
Outer ring defectA total of 160 feature vectors were collected from the four conditions, half of the feature vectors were used for training the classifier and others for signal classification, and the classification results are listed in Table6. Out of 80 test feature vectors, just two cases were not correctly classified, and the overall classification rate is 97.5%.Table 6
Fault diagnosis using CEEMD and time series model.
Fault type
Test sample
Classification results
Classification rate [%]
Overall classification rate [%]
No defect
Inner ring defect
Rolling element defect
Outer ring defect
No defect
20
20
0
0
0
100
97.5
Inner ring defect
20
0
19
1
0
95
Rolling element defect
20
0
1
19
0
95
Outer ring defect
20
0
0
0
20
100For comparison, Tables7 and 8 list classification results based on time series modeling using measured signal directly and based on EMD and time series model method. From the comparison results, the proposed method is efficient for rolling bearing fault diagnosis, and the overall classification rate of the proposed method is higher to a certain extent than the other two methods mentioned above.Table 7
Fault diagnosis using time series model only.
Fault type
Test sample
Classification results
Classification rate [%]
Overall classification rate [%]
No defect
Inner ring defect
Rolling element defect
Outer ring defect
No defect
20
19
1
0
0
95
90.0
Inner ring defect
20
1
17
2
0
85
Rolling element defect
20
0
2
17
1
85
Outer ring defect
20
0
0
1
19
95Table 8
Fault diagnosis using EMD and time series model.
Fault type
Test sample
Classification results
Classification rate [%]
Overall classification rate [%]
No defect
Inner ring defect
Rolling element defect
Outer ring defect
No defect
20
19
1
0
0
95
93.75
Inner ring defect
20
0
18
2
0
90
Rolling element defect
20
0
1
19
0
95
Outer ring defect
20
0
0
1
19
95


## 3.1. Evaluation Using Simulated Signals

To demonstrate the validity of the method proposed in this study, three signalsx
1
(
t
), x
2
(
t
), and x
3
(
t
) are simulated as shown in Figure 7. The signal x
1
(
t
) consists of a Gaussian-type impulse interference, a cosine component with 10 Hz frequency, a trend term, and white noise. The signal x
2
(
t
) consists of a Gaussian-type impulse interference, a square wave with 65% duty ratio, a trend term, and white noise. The signal x
3
(
t
) consists of a Gaussian-type impulse interference, a sawtooth wave with 15 Hz frequency, a trend term, and white noise.Figure 7
Signal waveforms ofx
1
(
t
), x
2
(
t
), and x
3
(
t
).Figure8 shows the results of the CEEMD of signals x
1
(
t
), x
2
(
t
), and x
3
(
t
). Correlation coefficients between filtered signal and each IMF are illustrated in Table 1.Table 1
Correlation coefficients between filtered signal and each IMF.
Signal
Correlation coefficient
IMF1
IMF2
IMF3
IMF4
IMF5
IMF6
IMF7
IMF8
x
1
(
t
)
−0.0031
−0.0009
0.0371
0.4096
0.9668
0.2428
0.1273
−0.0448
x
2
(
t
)
0.0051
0.0004
0.0435
0.2111
0.4695
0.8887
0.7214
−0.0201
x
3
(
t
)
−0.0234
−0.0154
0.0286
0.5900
0.8953
0.1649
0.1887
−0.0210The decomposition results by CEEMD.
(a)
x
1
(
t
)
(b)
x
2
(
t
)
(c)
x
3
(
t
)It can be seen in Table1 that the IMF which is closely correlated to the filtered signal is IMF5 for both signal x
1
(
t
) and signal x
3
(
t
) and IMF6 for signal x
2
(
t
). They are used to construct the AR models, and the corresponding feature vectors are estimated as shown in Table 2. After scalar quantization, the feature vectors are used to train the HMM for signal classification.Table 2
Model parameter estimation results.
Signal
Model parameter
φ
1
φ
2
φ
3
φ
4
φ
5
φ
6
x
1
(
t
)
4.7183
−9.1103
9.2034
−5.1408
1.5207
−0.1914
x
2
(
t
)
4.8894
−9.8269
10.3945
−6.1194
1.9153
−0.2531
x
3
(
t
)
4.8718
−9.9955
11.1616
−7.2529
2.6430
−0.4282A total of 90 feature vectors were collected from three groups of signals using the proposed approach. One-third of the feature vectors in each condition were used for training the classifier and others were used for testing. The results of the signal classification are listed in Table3.Table 3
Signal classification results.
Signal type
Test sample
Classification results
Classification rate [%]
Overall classification rate [%]
x
1
(
t
)
x
2
(
t
)
x
3
(
t
)
x
1
(
t
)
20
19
1
0
95
96.7
x
2
(
t
)
20
0
19
1
95
x
3
(
t
)
20
0
0
20
100Results in Table3 indicate that the presented method based on CEEMD and time series modeling can effectively identify different signals, and the overall classification rate is 96.7%. For the purpose of comparison, the signal classification rates use the method based on time series modeling only, and the method based on EMD and time series modeling is also calculated. 88.3% and 93.3% classification rates are obtained, respectively. It is obvious that efficiency of the signal classification method proposed in this paper is improved to a certain extent.


## 3.2. Evaluation Using Experimental Data

In order to illustrate the practicability and effectiveness of the proposed method, a bearing fault data set from the electrical engineering laboratory of Case Western Reserve University is analyzed [24]. The data set is acquired from the test stand shown in Figure 9, where it consists of a 2 hp motor, a torque transducer, a dynamometer, and control electronics. The test bearings support the motor shaft which is the deep groove ball bearings with the type of 6205-2RS JEMSKF. Vibration data was collected at 12,000 samples per second using accelerometers, which are attached to the housing with magnetic bases. The motor load level was controlled by the fan in the right side of Figure 9.Figure 9
Bearing test stand.Figure10 illustrates representative waveforms of the sample vibration signals measured from the test bearings under four initial conditions: (a) signal from a healthy bearing, (b) signal from a bearing with inner ring defect, (c) signal from a bearing with rolling element defect, and (d) signal from a bearing with outer ring defect. These signals were measured under 0 hp motor load with the motor speed of 1797 rpm. The decomposed IMFs of these signals are shown in Figure 11.Vibration signal waveforms of different conditions.
(a)
(b)
(c)
(d)The decomposition results by CEEMD under different conditions.
(a)
No defect
(b)
Inner ring defect
(c)
Rolling element defect
(d)
Outer ring defectCorrelation coefficients calculated between the filtered signal and each IMF are shown in Table4.Table 4
Correlation coefficients between filtered signals and each IMF.
Signal
Correlation coefficient
IMF1
IMF2
IMF3
IMF4
IMF5
IMF6
IMF7
IMF8
(a)
0.4135
0.7538
0.4381
0.4880
0.4356
0.1792
0.0971
−0.0056
(b)
0.8794
0.4275
0.2583
0.1337
0.0421
0.0285
−0.0009
−0.0074
(c)
0.9509
0.2180
0.2325
0.1337
0.0821
0.0350
−0.0017
0.0009
(d)
0.9878
0.1267
0.0636
0.0509
0.0136
0.0060
−0.0008
−0.0068The IMF which is closely correlated to the filtered signal is IMF2 for signal (a) and IMF1 for signals (b), (c), and (d), respectively. These IMFs are used for AR model construction. The model order estimation curves of the four conditions based on the principle of FPE criterion are shown in Figure12. We can see that when the model order is 6, each model's residual tends to be stable. Therefore the model order is selected as 6, and the results of parameters estimation are listed in Table 5.Table 5
Model parameter estimation results.
Signal
Model parameter
φ
1
φ
2
φ
3
φ
4
φ
5
φ
6
(a)
3.1280
−4.7797
4.2245
−2.1489
0.4241
0.0356
(b)
0.2084
−1.3585
0.5142
−0.6356
0.3471
−0.0422
(c)
0.1335
−1.6472
0.3941
−0.8473
0.2142
−0.1011
(d)
−0.1172
−1.2159
0.1178
−0.1283
0.1467
0.2533The model order estimation curves.
(a)
No defect
(b)
Inner ring defect
(c)
Rolling element defect
(d)
Outer ring defectThe parameters in Table5 were quantified by Lloyds algorithm in (12) as feature vectors for training the HMMs of different conditions. The results of quantization are revealed in Figure 13.The results of quantization.
(a)
No defect
(b)
Inner ring defect
(c)
Rolling element defect
(d)
Outer ring defectA total of 160 feature vectors were collected from the four conditions, half of the feature vectors were used for training the classifier and others for signal classification, and the classification results are listed in Table6. Out of 80 test feature vectors, just two cases were not correctly classified, and the overall classification rate is 97.5%.Table 6
Fault diagnosis using CEEMD and time series model.
Fault type
Test sample
Classification results
Classification rate [%]
Overall classification rate [%]
No defect
Inner ring defect
Rolling element defect
Outer ring defect
No defect
20
20
0
0
0
100
97.5
Inner ring defect
20
0
19
1
0
95
Rolling element defect
20
0
1
19
0
95
Outer ring defect
20
0
0
0
20
100For comparison, Tables7 and 8 list classification results based on time series modeling using measured signal directly and based on EMD and time series model method. From the comparison results, the proposed method is efficient for rolling bearing fault diagnosis, and the overall classification rate of the proposed method is higher to a certain extent than the other two methods mentioned above.Table 7
Fault diagnosis using time series model only.
Fault type
Test sample
Classification results
Classification rate [%]
Overall classification rate [%]
No defect
Inner ring defect
Rolling element defect
Outer ring defect
No defect
20
19
1
0
0
95
90.0
Inner ring defect
20
1
17
2
0
85
Rolling element defect
20
0
2
17
1
85
Outer ring defect
20
0
0
1
19
95Table 8
Fault diagnosis using EMD and time series model.
Fault type
Test sample
Classification results
Classification rate [%]
Overall classification rate [%]
No defect
Inner ring defect
Rolling element defect
Outer ring defect
No defect
20
19
1
0
0
95
93.75
Inner ring defect
20
0
18
2
0
90
Rolling element defect
20
0
1
19
0
95
Outer ring defect
20
0
0
1
19
95


## 4. Conclusions

Aiming at diagnosing rolling bearing faults, a hybrid approach based on CEEMD and time series modeling is proposed in this paper. The CEEMD method can decompose the nonstationary signal into a series of IMFs with low computation. AR model is an effective approach to extract the fault feature of the vibration signals and the fault pattern can be identified directly by the extracted fault features without establishing the mathematical model and studying the fault mechanism of the system. In this paper, the CEEMD method is used as a pretreatment, which can increase the accuracy of the AR model for the measured signal, and the AR model of the IMF which is closely correlated to the filtered signal is established to extract the fault feature parameters. Comparing to the EMD-AR approach and the direct modeling approach where raw signals are directly used as input for AR modeling, a higher classification rate was shown to be achieved by using the new approach (e.g., 96.7% for simulated signals and 97.5% for experimental data). Meanwhile we anticipate that the proposed method can also be used for incipient fault diagnosis in rolling bearing, where further experiments are needed to verify the accuracy. Since the approach presented in this study is generic in nature, it can be readily adapted to a broad range of applications for machine fault diagnosis.



---

*Source: 101867-2014-07-07.xml*