# Remaining Useful Life Prediction of Rolling Bearings Using PSR, JADE, and Extreme Learning Machine

Authors: Yongbin Liu; Bing He; Fang Liu; Siliang Lu; Yilei Zhao; Jiwen Zhao
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

Rolling bearings play a pivotal role in rotating machinery. The degradation assessment and remaining useful life (RUL) prediction of bearings are critical to condition-based maintenance. However, sensitive feature extraction still remains a formidable challenge. In this paper, a novel feature extraction method is introduced to obtain the sensitive features through phase space reconstitution (PSR) and joint with approximate diagonalization of Eigen-matrices (JADE). Firstly, the original features are extracted from bearing vibration signals in time and frequency domain. Secondly, the PSR is applied to embed the original features into high dimensional phase space. The between-class and within-class scatter (S
S) are calculated to evaluate the feature sensitivity through the phase point distribution of different degradation stages and then different weights are assigned to the corresponding features based on the calculated S
S. Thirdly, the JADE is employed to fuse the weighted features to obtain the advanced features which can better reflect the bearing degradation process. Finally, the advanced features are input into the extreme learning machine (ELM) to train the RUL prediction model. A set of experimental case studies are carried out to verify the effectiveness of the proposed method. The results show that the extracted advanced features can better reflect the degradation process compared to traditional features and could effectively predict the RUL of bearing.


---

## Body

## 1. Introduction

Currently, Prognostics and Health Management (PHM) is an important research area in industry, which aims at increasing availability and decreasing the downtime by predicting the residual life of machine. Rolling bearing as a fundamental component is widely used in the rotary machine. Bearing failure is one of the principal reasons for the machine damage [1]. In brief, condition monitoring of bearing is becoming more and more significant, especially in the RUL prediction, which can reduce the costs of maintenance. There are three main approaches, namely, model based, data driven, and hybrid [2–4]. Among them, the data driven methods are easier to deploy than others, which are based on the relation between degradation state and measurement data [5].One of the challenges for RUL prediction based on data driven is sensitive feature extraction from vibration signals that are always mixed by noise [6, 7]. Many modern mathematical methods are employed to overcome this challenge. These methods offer systematic, scientific, and effective tactics for extracting the features that can reflect the bearing degradation. There are some representative feature extraction methods via statistical calculation, such as time domain analysis, frequency domain analysis, and time-frequency domain analysis [8–10]. However, these original features are often redundant and with high dimensionality. Therefore, dimensionality reduction methods such as the principal component analysis (PCA) and the independent component analysis (ICA) are employed to reduce the redundancy for machine fault diagnosis [11–13]. On the other hand, the nonlinear dimensionality reduction methods are developed in recent years, such as local embedding structure (LLE) [14], local tangent space information (LTSA) [15], and local adjacency relations (LE) [16]. However, they are often with complex algorithm and need a large amount of calculation [17].In this paper, a new dimensionality reduction technique is introduced based on approximate diagonalization of Eigen-matrices (JADE) algorithm, which is originally used in blind identification. Generally, the mixed model without environmental noise can be effectively identified by JADE. However, it has a poor performance when the signals are mixed with background noise [18]. To overcome this problem, a feature weighted fusion method combining the phase space reconstitution (PSR) [19] with the between-class and within-class scatters evaluation [20] is proposed in this paper. As different kinds of signals demonstrate different structure in phase trajectory, the PSR can better grasp the nature and regularity of the time series.After extracting feature from vibration signals, machine learning algorithms can be used to establish the bearing RUL prediction model, such as the support vector machine (SVM) [21–23] and the ANN [24]. In recent years, a new machine learning algorithm called extreme learning machine (ELM) has been introduced by Huang et al. It is a novel algorithm of single hidden layer feed-forward neural networks (SLFNs), which has an extremely fast learning speed and advantageous generalization capability through randomly choosing hidden nodes and analytically determining output weights [25]. Moreover, ELM has a commendable performance in both classification and regression with a simple structure. Compared with back-propagation (BP) feed-forward network learning algorithm, the training speed of ELM is much faster while obtaining better generalization [26, 27]. In view of these advantages, the ELM is selected in this paper to establish bearing RUL prediction model.The rest of the paper is organized as follows. Section2 gives a brief introduction of PSR, JADE, and ELM. The proposed RUL prediction procedure for bearings is illustrated in Section 3. Experimental validations and comparison with other linear dimensionality reduction methods are presented in Section 4. Finally, the conclusions are drawn in Section 5.


## 2. Theory Background

### 2.1. Phase Space Reconstitution (PSR)

Chaotic time series analysis is a method for nonlinear signal processing. By extending the time series to high dimensional phase space, this method can grasp the nature of time series much better. The research of PSR on chaotic time series was started from Packard et al. [28]. Suppose that the time series is x
=
(
x
1
,
x
2
,
…
,
x
n
), the embedding dimension is chosen as d, and the time delay is selected as τ. And then a phase point is constructed as(1)
X
i
d
=
x
i
,
x
i
+
τ
,
…
,
x
i
+
d
-
2
τ
,
x
i
+
d
-
1
τ
,where n is the number of time series and i is the index of the phase point. However, the acquired signals are blended with noise, so the parameters should be determined under a reasonable rule [29]. In this paper, the autocorrelation function is employed to identify the time delay, which is defined as(2)
C
i
τ
=
1
/
n
∑
i
=
1
n
x
i
+
τ
-
x
-
x
i
-
x
-
1
/
n
∑
i
=
1
n
x
i
-
x
-
2
,where x
-
=
1
/
n
∑
i
=
1
n
x
(
i
). According to Rosenstein’s theories [30], the minimum τ that makes C
i
(
τ
) lower than 1
-
1
/
e can be considered as the most suitable time delay. And Cao’s algorithm [31] is employed to identify the embedding dimension. Suppose that the nearest phase point to X
i
(
d
) is X
n
e
a
r
i
,
d
(
d
). In accordance with the embedding theorem, if d is the suitable embedding dimension, then any two points that are close to the d dimensional reconstructed phase space will continue to be close to the d
+
1 dimensional space(3)
a
i
,
d
=
X
i
d
+
1
-
X
n
e
a
r
i
,
d
+
1
d
+
1
X
i
d
-
X
n
e
a
r
i
,
d
d
,
E
d
=
1
n
-
d
τ
∑
i
=
1
n
-
d
τ
a
i
,
d
,
E
1
d
=
E
d
+
1
E
d
.With the increase ofd, E
1
(
d
) will be gradually stabilized; therefore, the minimum d that stabilizes E
1
(
d
) is set as the embedding dimension. However, it is difficult to determine whether the curve of E
1
(
d
) is stable due to the limited length of samples. Therefore, supplementary criteria are proposed as shown in(4)
E
∗
d
=
1
n
-
d
τ
∑
i
=
1
n
-
d
τ
x
i
+
d
τ
-
x
n
e
a
r
i
,
d
+
d
τ
,
E
2
d
=
E
∗
d
+
1
E
∗
d
,where E
2
(
d
) will be changed with the increasing of d if the sample is deterministic. Besides, more details of the description are presented in [32].


### 2.2. Joint Approximate Diagonalization of Eigen-Matrices (JADE)

The JADE algorithm is widely used in blind source separation (BSS) [33]; it is based on a standard linear data model that can be expressed as follows: (5)
x
t
=
A
s
t
+
n
t
,where x
(
t
) is the source signal that mixed with noise, n
(
t
) is additive noise, and s
(
t
) is the source signal, which can be obtained by transfer matrix A. The JADE can be summarized as follows:(a)
Signal whitening is as follows:(6)
z
t
=
W
x
t
,
W
=
μ
1
-
σ
-
1
/
2
h
1
,
…
,
μ
n
-
σ
-
1
/
2
h
n
H
,
wherez
(
t
) is whitened signal, μ
1
,
μ
2
,
…
,
μ
n is the eigenvalue of sample, h
1
,
h
2
,
…
,
h
n is the corresponding eigenvector, and σ is the noise variance.
(b)
Calculate the fourth-order cumulants [34](7)
Q
z
=
C
u
m
z
i
,
z
j
,
z
k
,
z
l
∣
1
≤
i
,
j
,
k
,
l
≤
d
,
whered is the dimensionality of the vector z.
(c)
Calculate the unitary matrix(8)
U
^
=
a
r
g
m
i
n
∑
i
o
f
f
U
#
Q
z
M
i
U
,
1
≤
i
≤
n
,
where# denotes the pseudo-inverse, “off” is the square of the nondiagonal elements, ‘‘a
r
g
m
i
n’’ is the argument of a complex number, and U is a rotation matrix.
(d)
The single model can be separated as(9)
s
t
=
A
^
x
t
,
whereA
^
=
U
^
W
#.In this paper, JADE is firstly applied in dimensionality reduction and feature extraction. The single model can be obtained through Eigen-decomposition of the fourth-order cumulative matrix of the multivariate data. This not only is simpler than PCA and fast-ICA algorithm but also can improve the robustness of the results. However, the original features have different units and dimension, and it is difficult to distinguish between the contributions of each feature to the overall feature.


### 2.3. Extreme Learning Machine (ELM)

ELM is an efficient learning algorithm for classification and prediction in SLFNs, the structure of ELM as shown in Figure1. It can randomly generate the bias and input weights of hidden nodes. Moreover, it can immediately get the result without the puzzle of a local minimum [35].Figure 1
The structure of ELM.Suppose there areN samples (
X
i
,
t
i
), where X
i denotes the feature and t
i denotes the target value. Then the model of SLFNs with N
~ hidden nodes can be described as below:(10)
∑
i
=
1
N
~
β
i
g
i
X
j
=
∑
i
=
1
N
~
β
i
g
W
i
·
X
j
+
b
i
=
o
j
,
j
=
1
,
…
,
N
,where W
i
=
[
w
i
,
1
,
w
i
,
2
,
…
,
w
i
,
N
~
]
T is the weight vector that connects the ith hidden and the input nodes and b
i is the ith bias of hidden layer nodes, and they are randomly assigned. β
i
=
[
β
i
,
1
,
β
i
,
2
,
…
,
β
i
,
N
~
]
T is the weight vector connecting the ith hidden and the output nodes and o is the output value. To minimize the error between output o
j and target t
j, which is represented as ∑
j
=
1
N
o
j
-
t
j
=
0, there is an existing optimal value of β
i such that(11)
∑
i
=
1
N
~
β
i
g
W
i
·
X
j
+
b
i
=
t
j
,
j
=
1
,
…
,
N
.Then above N equation can be written compactly as(12)
H
β
=
T
,where(13)
H
W
,
b
,
X
=
g
W
1
·
X
1
+
b
1
⋯
g
W
N
~
·
X
1
+
b
N
~
⋮
⋯
⋮
g
W
1
·
X
N
+
b
1
⋯
g
W
N
~
·
X
N
+
b
N
~
N
×
N
~
.Thus the output weight β can be written as (14)
β
=
H
+
T
,where H
+
=
(
H
T
H
)
-
1
H
T.Although the ELM is faster than other machine learning algorithms, its accuracy is still dependent on feature extraction.


## 2.1. Phase Space Reconstitution (PSR)

Chaotic time series analysis is a method for nonlinear signal processing. By extending the time series to high dimensional phase space, this method can grasp the nature of time series much better. The research of PSR on chaotic time series was started from Packard et al. [28]. Suppose that the time series is x
=
(
x
1
,
x
2
,
…
,
x
n
), the embedding dimension is chosen as d, and the time delay is selected as τ. And then a phase point is constructed as(1)
X
i
d
=
x
i
,
x
i
+
τ
,
…
,
x
i
+
d
-
2
τ
,
x
i
+
d
-
1
τ
,where n is the number of time series and i is the index of the phase point. However, the acquired signals are blended with noise, so the parameters should be determined under a reasonable rule [29]. In this paper, the autocorrelation function is employed to identify the time delay, which is defined as(2)
C
i
τ
=
1
/
n
∑
i
=
1
n
x
i
+
τ
-
x
-
x
i
-
x
-
1
/
n
∑
i
=
1
n
x
i
-
x
-
2
,where x
-
=
1
/
n
∑
i
=
1
n
x
(
i
). According to Rosenstein’s theories [30], the minimum τ that makes C
i
(
τ
) lower than 1
-
1
/
e can be considered as the most suitable time delay. And Cao’s algorithm [31] is employed to identify the embedding dimension. Suppose that the nearest phase point to X
i
(
d
) is X
n
e
a
r
i
,
d
(
d
). In accordance with the embedding theorem, if d is the suitable embedding dimension, then any two points that are close to the d dimensional reconstructed phase space will continue to be close to the d
+
1 dimensional space(3)
a
i
,
d
=
X
i
d
+
1
-
X
n
e
a
r
i
,
d
+
1
d
+
1
X
i
d
-
X
n
e
a
r
i
,
d
d
,
E
d
=
1
n
-
d
τ
∑
i
=
1
n
-
d
τ
a
i
,
d
,
E
1
d
=
E
d
+
1
E
d
.With the increase ofd, E
1
(
d
) will be gradually stabilized; therefore, the minimum d that stabilizes E
1
(
d
) is set as the embedding dimension. However, it is difficult to determine whether the curve of E
1
(
d
) is stable due to the limited length of samples. Therefore, supplementary criteria are proposed as shown in(4)
E
∗
d
=
1
n
-
d
τ
∑
i
=
1
n
-
d
τ
x
i
+
d
τ
-
x
n
e
a
r
i
,
d
+
d
τ
,
E
2
d
=
E
∗
d
+
1
E
∗
d
,where E
2
(
d
) will be changed with the increasing of d if the sample is deterministic. Besides, more details of the description are presented in [32].


## 2.2. Joint Approximate Diagonalization of Eigen-Matrices (JADE)

The JADE algorithm is widely used in blind source separation (BSS) [33]; it is based on a standard linear data model that can be expressed as follows: (5)
x
t
=
A
s
t
+
n
t
,where x
(
t
) is the source signal that mixed with noise, n
(
t
) is additive noise, and s
(
t
) is the source signal, which can be obtained by transfer matrix A. The JADE can be summarized as follows:(a)
Signal whitening is as follows:(6)
z
t
=
W
x
t
,
W
=
μ
1
-
σ
-
1
/
2
h
1
,
…
,
μ
n
-
σ
-
1
/
2
h
n
H
,
wherez
(
t
) is whitened signal, μ
1
,
μ
2
,
…
,
μ
n is the eigenvalue of sample, h
1
,
h
2
,
…
,
h
n is the corresponding eigenvector, and σ is the noise variance.
(b)
Calculate the fourth-order cumulants [34](7)
Q
z
=
C
u
m
z
i
,
z
j
,
z
k
,
z
l
∣
1
≤
i
,
j
,
k
,
l
≤
d
,
whered is the dimensionality of the vector z.
(c)
Calculate the unitary matrix(8)
U
^
=
a
r
g
m
i
n
∑
i
o
f
f
U
#
Q
z
M
i
U
,
1
≤
i
≤
n
,
where# denotes the pseudo-inverse, “off” is the square of the nondiagonal elements, ‘‘a
r
g
m
i
n’’ is the argument of a complex number, and U is a rotation matrix.
(d)
The single model can be separated as(9)
s
t
=
A
^
x
t
,
whereA
^
=
U
^
W
#.In this paper, JADE is firstly applied in dimensionality reduction and feature extraction. The single model can be obtained through Eigen-decomposition of the fourth-order cumulative matrix of the multivariate data. This not only is simpler than PCA and fast-ICA algorithm but also can improve the robustness of the results. However, the original features have different units and dimension, and it is difficult to distinguish between the contributions of each feature to the overall feature.


## 2.3. Extreme Learning Machine (ELM)

ELM is an efficient learning algorithm for classification and prediction in SLFNs, the structure of ELM as shown in Figure1. It can randomly generate the bias and input weights of hidden nodes. Moreover, it can immediately get the result without the puzzle of a local minimum [35].Figure 1
The structure of ELM.Suppose there areN samples (
X
i
,
t
i
), where X
i denotes the feature and t
i denotes the target value. Then the model of SLFNs with N
~ hidden nodes can be described as below:(10)
∑
i
=
1
N
~
β
i
g
i
X
j
=
∑
i
=
1
N
~
β
i
g
W
i
·
X
j
+
b
i
=
o
j
,
j
=
1
,
…
,
N
,where W
i
=
[
w
i
,
1
,
w
i
,
2
,
…
,
w
i
,
N
~
]
T is the weight vector that connects the ith hidden and the input nodes and b
i is the ith bias of hidden layer nodes, and they are randomly assigned. β
i
=
[
β
i
,
1
,
β
i
,
2
,
…
,
β
i
,
N
~
]
T is the weight vector connecting the ith hidden and the output nodes and o is the output value. To minimize the error between output o
j and target t
j, which is represented as ∑
j
=
1
N
o
j
-
t
j
=
0, there is an existing optimal value of β
i such that(11)
∑
i
=
1
N
~
β
i
g
W
i
·
X
j
+
b
i
=
t
j
,
j
=
1
,
…
,
N
.Then above N equation can be written compactly as(12)
H
β
=
T
,where(13)
H
W
,
b
,
X
=
g
W
1
·
X
1
+
b
1
⋯
g
W
N
~
·
X
1
+
b
N
~
⋮
⋯
⋮
g
W
1
·
X
N
+
b
1
⋯
g
W
N
~
·
X
N
+
b
N
~
N
×
N
~
.Thus the output weight β can be written as (14)
β
=
H
+
T
,where H
+
=
(
H
T
H
)
-
1
H
T.Although the ELM is faster than other machine learning algorithms, its accuracy is still dependent on feature extraction.


## 3. The Proposed Method for Bearing RUL Prediction

In the abovementioned methods, the PSR has superior merits for characterizing nonlinear time series. The JADE focuses on feature fusion or dimensionality reduction. The ELM is an effective learning algorithm with fast training speed. A bearing RUL prediction scheme combined with the three methods is introduced in this paper.The scheme is exhibited in Figure2 and can be described as follows.Figure 2
Principle of the proposed method based on PSR and JADE.Step 1 (original feature extraction).
A total of 16 statistical features are shown in Table1 including 8 time domain features and 8 frequency domain features. F
1 to F
8 are the mean, root mean square, RMS amplitude, absolute average, skewness, waveform index, impulsion index, and kurtosis index, respectively. F
9
–
F
16 are the frequency domain feature, which indicates the degree of dispersion or concentration of the spectrum, and the change of the dominant frequency band.Table 1
List of the original extracted features.
Time domain
Frequency domain
F
1
=
1
N
∑
i
=
1
N
x
i
F
2
=
1
N
∑
i
=
1
N
x
i
2
F
9
=
1
N
∑
i
=
1
N
s
i
F
10
=
1
N
∑
j
=
1
N
s
j
-
1
N
∑
i
=
1
N
s
i
2
F
3
=
1
N
∑
i
=
1
N
x
i
2
F
4
=
1
N
∑
i
=
1
N
x
i
F
11
=
(
1
/
N
)
∑
j
=
1
N
(
s
j
-
1
/
N
∑
i
=
1
N
s
i
)
3
F
10
3
F
12
=
∑
i
=
1
N
f
i
s
i
∑
j
=
1
N
s
j
F
5
=
1
N
∑
i
=
1
N
x
i
3
F
6
=
1
/
N
∑
i
=
1
N
x
i
2
F
4
F
13
=
1
N
∑
i
=
1
N
s
i
f
i
-
F
12
2
F
14
=
∑
i
=
1
N
f
i
2
s
i
∑
j
=
1
N
s
j
F
7
=
max
⁡
x
1
/
N
∑
i
=
1
N
x
i
F
8
=
1
/
N
∑
i
=
1
N
x
i
4
1
/
N
∑
j
=
1
N
x
j
2
4
F
15
=
∑
i
=
1
N
f
i
4
s
i
∑
j
=
1
N
f
j
2
s
j
F
16
=
∑
i
=
1
N
f
i
2
s
i
∑
j
=
1
N
s
j
∑
k
=
1
N
f
k
4
s
k
WhereS
i is a spectrum for i
=
1,2
,
…
,
N, N is the number of spectrum lines, and f
i is the frequency value of the ith spectrum line.Step 2 (PSR for each feature).
The autocorrelation function is employed to select the time delay and Cao’s algorithm is used to determine the embedding dimension. All the features are reconstructed by PSR with the same parameter.Step 3 (between-class and within-class scatter (S
S) for feature weighting).
Based on the feature phase space distribution,S
S is used to describe the classification ability of features quantitatively for different degradation states. Suppose that there are C classes that are established by degradation states and the ith class vector is x
i
=
(
x
1
i
,
x
2
i
,
…
,
x
n
i
i
), where n
i is the number of samples in this class. Then S
S is defined as follows:(15)
S
S
=
t
r
a
c
e
S
b
S
w
,
S
b
=
∑
i
=
1
C
p
i
m
i
-
m
2
,
S
w
=
∑
i
=
1
C
p
i
1
n
i
∑
k
=
1
n
i
x
k
i
-
m
i
2
,
p
i
=
n
i
∑
j
=
1
C
n
j
,
m
i
=
1
n
i
∑
k
=
1
n
i
x
k
i
,
m
=
∑
i
=
1
C
p
i
m
i
,where m
i is the mean vector for samples in the ith class and m is the mean vector for all classes. The between-class scatter S
b denotes the scattered degree in different classes while the within-class scatter S
w indicates the concentration among the same class. Therefore, the classification ability can be evaluated through the S
S. The weight of each original feature is calculated as follows:(16)
w
i
=
S
S
i
∑
j
=
1
j
=
16
S
S
j
.Step 4 (JADE for feature fusion).
In this step, every feature that is extracted in Step1 can be seen as a “source signal,” and the JADE is introduced to separate the single model from the 16 source signals that have been normalized and weighted. However, the separated unidimensional vector is mixed by noise. Therefore, the exponentially weighted moving average (EWMA) is used for denoising, which is defined as follows:(17)
x
t
=
1
-
α
x
t
-
1
+
α
x
t
,where α
∈
[
0,1
] is the smoothing factor, and it generally takes between 0.05 and 0.25 [36].Step 5 (prediction model training).
Based on Steps1–4, the fused feature vector with real RUL (
X
j
,
t
j
) is input into ELM, where j is the index of the samples. The bias and input weights of hidden nodes in ELM are randomly generated. Then the model with N
~ hidden nodes can be constructed as ∑
i
=
1
N
~
β
i
g
i
(
X
j
)
=
∑
i
=
1
N
~
β
i
g
(
W
i
·
X
j
+
b
i
)
=
o
j, where w
i is the weight connecting the ith hidden and the input nodes, b
i is the bias of hidden layer nodes, which are randomly assigned, and o is the output value. β
i is the weight connecting the ith hidden and output nodes, where the best value can minimize the error between output o
j and real RUL t
j. And then the fused feature vector of testing samples is input into the trained prediction model for running bearing RUL prediction.


## 4. Experimental Results and Analysis

### 4.1. Experimental Data Set Description

In order to evaluate the effectiveness of the proposed feature fusion scheme for RUL prediction, the entire life cycle bearing data originated from the Center for Intelligent Maintenance Systems [37] is analyzed. The experimental data sets are collected from run-to-failure bearing test on the designed test platform (Figure 3). Four bearings are installed on a shaft. The rotation speed is kept at 2000 r/min by using an AC motor coupled to the shaft via rub belts. A 6000 lb radial load is applied onto the bearing and shaft by using a spring mechanism. A PCB 353B33 High Sensitivity Quartz ICP Accelerometer was installed on bearing housing. A vibration data sample of 20480 points was collected every 10 minutes by a National Instruments DAQCard-6062E data acquisition card. The sampling rate is set as 20 kHz, the vibration data are recorded on the test bench, and more details can be found in [38].Figure 3
Bearing test platform.


### 4.2. Original Feature Extraction

In this paper, 100 groups of training samples are selected from the run-to-failure bearing data, and the data samples adjacent to the training data samples are regarded as testing samples. The length of every group’s data is set to 1680 points. Take outer-race-fault case as an example; the acquired run-to-failure bearing data is(
X
1
,
X
2
,
…
,
X
984
), where X
i is the ith acquired data sample and i ranges from 1 to 984. The training samples can be selected as X
90
,
X
99
,
…
,
X
981 and the testing samples can be selected as X
91
,
X
100
,
…
,
X
982, where X
j
=
x
1
,
x
2
,
…
,
x
1680 and j is the index of samples and ranges from 1 to 100.The 16 original features (as shown in Table1) are extracted from the training samples and illustrated in Figures 4–6. However, these features are not robust in characterizing the degradation process. Take F
12 as an example; it performs well in inner-race-fault case and roller-fault case but performs badly in outer-race-fault case. Therefore, it is very important to extract a robust and favorable feature from the original features, which can better express the bearing degradation. In addition, from these extracted features, it can probably be found that the time ratio of each operation stage (including normal, slight degradation, and serious degradation) in the outer-race-fault case, inner-race-fault case, and roller-fault case is about 0.5 : 0.4 : 0.1, 0.6 : 0.3 : 0.1, and 0.5 : 0.2 : 0.3, respectively.Figure 4
The extracted 16 features of outer-race-fault case.Figure 5
The extracted 16 features of inner-race-fault case.Figure 6
The extracted 16 features of roller-fault case.


### 4.3. Phase Space Reconstitution

As the high dimensional phase space can grasp the nature of time series much better, the PSR is then applied to reconstruct the original feature vector to the phase space before the following feature contribution evaluation usingS
S.In this paper, the autocorrelation function and Cao’s algorithm are employed to identify the time delay and embedding dimension. TakeF
1 of outer-race-fault case as an example; the results are shown in Figure 7. Figure 7(a) shows that C
i
(
τ
) changes with the increase of delay time, and the minimum τ that makes C
i
(
τ
) lower than 1
-
1
/
e could be considered as the most optimal time delay. From Figure 7(b), it can be seen that E
1
(
d
) keeps stable when d is larger than 10, and E
2
(
d
) randomly fluctuates with d. Consequently, the time delay is set to 1 and embedding dimension is set to 10. Then other features are processed as the abovementioned methods, and the phase points distribution of different operational stages is shown in Figures 8–10 (two dimensions are chosen to reflect the phase point distribution for better observation experience).Figure 7
(a) Curve ofC
i
(
τ
) and (b) curve of E
1
(
d
) and E
2
(
d
).
(a)
(b)Figure 8
The phase points distribution of outer-race-fault case in two dimensions.Figure 9
The phase points distribution of inner-race-fault case in two dimensions.Figure 10
The phase points distribution of roller-fault case in two dimensions.


### 4.4. Feature Fusion Based on JADE and SS

S
S is then employed to evaluate the contributions of the feature. As shown in Figure 11, the abscissa is the feature number which has been described in Table 1. The ordinate is the value of S
S which is the evaluation index of feature. The greater the S
S, the better the identifiable capacity of bearing degradation states. Take inner-race-fault case as an example; it can be seen that F
1, F
5, F
7, and F
8 cannot effectively characterize or distinguish different degradation stages, and F
3 provides excellent preference for distinguishing of different degradation states. The weights of every feature can be calculated as w
i
=
S
S
i
/
∑
j
=
1
j
=
16
S
S
j. Besides, F
12 performs well in roller-fault case but performs badly in outer-race-fault case. Therefore, it is very important to extract a robust feature that could be adapted to different cases.Figure 11
TheS
S of 16 features in (a) outer-race-fault case, (b) inner-race-fault case, and (c) roller-fault case.
(a)
(b)
(c)The JADE is then applied to fuse the weighted features. In order to improve the sensitivity and reliability of the fused feature in the bearing slight degradation, the EWMA is employed for denoising according to (17), where the smoothing factor is selected as 0.2. After all of the steps, the new feature vector can be extracted and shown in Figure 12.Figure 12
The extracted features of bearing degradation in three cases.


### 4.5. Prediction Performance

In this step, the fused feature vector and the corresponding real RUL are input into the ELM to construct the RUL prediction model. The data samples adjacent to the training data samples are used to test the model. Results are shown in Figure13. It can be seen that the predicted RUL is coinciding fairly well with the real RUL, which indicates that the feature fusion method and constructed ELM model can achieve an ideal result. However, the RUL prediction accuracy is relatively high in the first stage as the vibration is relatively stable. The slight fault stage is always complex and sharp fluctuation also happens at this stage. The vibration at the serious fault stage is more regular, which can lead to higher prediction accuracy.Figure 13
Predicted RUL: (a) outer-race-fault case, (b) inner-race-fault case, and (c) roller-fault case.
(a)
(b)
(c)The PCA and ICA as the representative linear dimensionality reduction methods are applied on the same data to compare with the proposed method. The root mean square error (RMSE) is selected to evaluate the RUL prediction accuracy which can be defined as(18)
R
M
S
E
=
∑
i
=
1
n
t
i
-
y
i
2
,
i
=
1,2
,
…
,
n
,where t
i is the real RUL and y
i is the predicted RUL.Considering the problem of repeatability, the ELM is applied on the same training data 10 times, and the average of the RMSE is selected to evaluate the prediction accuracy. The results of comparison are shown in Tables2–4. It can be seen that the proposed feature extraction is better than the traditional two methods, which proves the merits of the proposed method for fusing original features.Table 2
Predicted accuracy of outer-race-fault case.
Training
Testing
Number of training data
RMSE
Number of testing data
RMSE
PCA
100
12.549
100
13.526
ICA
100
12.605
100
14.108
JADE
100
8.872
100
10.018Table 3
Predicted accuracy of inner-race-fault case.
Training
Testing
Number of training data
RMSE
Number of testing data
RMSE
PCA
100
16.073
100
17.548
ICA
100
16.165
100
17.103
JADE
100
14.552
100
14.846Table 4
Predicted accuracy of roller-fault case.
Training
Testing
Number of training data
RMSE
Number of testing data
RMSE
PCA
100
9.668
100
9.124
ICA
100
9.650
100
9.077
JADE
100
9.105
100
8.227


## 4.1. Experimental Data Set Description

In order to evaluate the effectiveness of the proposed feature fusion scheme for RUL prediction, the entire life cycle bearing data originated from the Center for Intelligent Maintenance Systems [37] is analyzed. The experimental data sets are collected from run-to-failure bearing test on the designed test platform (Figure 3). Four bearings are installed on a shaft. The rotation speed is kept at 2000 r/min by using an AC motor coupled to the shaft via rub belts. A 6000 lb radial load is applied onto the bearing and shaft by using a spring mechanism. A PCB 353B33 High Sensitivity Quartz ICP Accelerometer was installed on bearing housing. A vibration data sample of 20480 points was collected every 10 minutes by a National Instruments DAQCard-6062E data acquisition card. The sampling rate is set as 20 kHz, the vibration data are recorded on the test bench, and more details can be found in [38].Figure 3
Bearing test platform.


## 4.2. Original Feature Extraction

In this paper, 100 groups of training samples are selected from the run-to-failure bearing data, and the data samples adjacent to the training data samples are regarded as testing samples. The length of every group’s data is set to 1680 points. Take outer-race-fault case as an example; the acquired run-to-failure bearing data is(
X
1
,
X
2
,
…
,
X
984
), where X
i is the ith acquired data sample and i ranges from 1 to 984. The training samples can be selected as X
90
,
X
99
,
…
,
X
981 and the testing samples can be selected as X
91
,
X
100
,
…
,
X
982, where X
j
=
x
1
,
x
2
,
…
,
x
1680 and j is the index of samples and ranges from 1 to 100.The 16 original features (as shown in Table1) are extracted from the training samples and illustrated in Figures 4–6. However, these features are not robust in characterizing the degradation process. Take F
12 as an example; it performs well in inner-race-fault case and roller-fault case but performs badly in outer-race-fault case. Therefore, it is very important to extract a robust and favorable feature from the original features, which can better express the bearing degradation. In addition, from these extracted features, it can probably be found that the time ratio of each operation stage (including normal, slight degradation, and serious degradation) in the outer-race-fault case, inner-race-fault case, and roller-fault case is about 0.5 : 0.4 : 0.1, 0.6 : 0.3 : 0.1, and 0.5 : 0.2 : 0.3, respectively.Figure 4
The extracted 16 features of outer-race-fault case.Figure 5
The extracted 16 features of inner-race-fault case.Figure 6
The extracted 16 features of roller-fault case.


## 4.3. Phase Space Reconstitution

As the high dimensional phase space can grasp the nature of time series much better, the PSR is then applied to reconstruct the original feature vector to the phase space before the following feature contribution evaluation usingS
S.In this paper, the autocorrelation function and Cao’s algorithm are employed to identify the time delay and embedding dimension. TakeF
1 of outer-race-fault case as an example; the results are shown in Figure 7. Figure 7(a) shows that C
i
(
τ
) changes with the increase of delay time, and the minimum τ that makes C
i
(
τ
) lower than 1
-
1
/
e could be considered as the most optimal time delay. From Figure 7(b), it can be seen that E
1
(
d
) keeps stable when d is larger than 10, and E
2
(
d
) randomly fluctuates with d. Consequently, the time delay is set to 1 and embedding dimension is set to 10. Then other features are processed as the abovementioned methods, and the phase points distribution of different operational stages is shown in Figures 8–10 (two dimensions are chosen to reflect the phase point distribution for better observation experience).Figure 7
(a) Curve ofC
i
(
τ
) and (b) curve of E
1
(
d
) and E
2
(
d
).
(a)
(b)Figure 8
The phase points distribution of outer-race-fault case in two dimensions.Figure 9
The phase points distribution of inner-race-fault case in two dimensions.Figure 10
The phase points distribution of roller-fault case in two dimensions.


## 4.4. Feature Fusion Based on JADE and SS

S
S is then employed to evaluate the contributions of the feature. As shown in Figure 11, the abscissa is the feature number which has been described in Table 1. The ordinate is the value of S
S which is the evaluation index of feature. The greater the S
S, the better the identifiable capacity of bearing degradation states. Take inner-race-fault case as an example; it can be seen that F
1, F
5, F
7, and F
8 cannot effectively characterize or distinguish different degradation stages, and F
3 provides excellent preference for distinguishing of different degradation states. The weights of every feature can be calculated as w
i
=
S
S
i
/
∑
j
=
1
j
=
16
S
S
j. Besides, F
12 performs well in roller-fault case but performs badly in outer-race-fault case. Therefore, it is very important to extract a robust feature that could be adapted to different cases.Figure 11
TheS
S of 16 features in (a) outer-race-fault case, (b) inner-race-fault case, and (c) roller-fault case.
(a)
(b)
(c)The JADE is then applied to fuse the weighted features. In order to improve the sensitivity and reliability of the fused feature in the bearing slight degradation, the EWMA is employed for denoising according to (17), where the smoothing factor is selected as 0.2. After all of the steps, the new feature vector can be extracted and shown in Figure 12.Figure 12
The extracted features of bearing degradation in three cases.


## 4.5. Prediction Performance

In this step, the fused feature vector and the corresponding real RUL are input into the ELM to construct the RUL prediction model. The data samples adjacent to the training data samples are used to test the model. Results are shown in Figure13. It can be seen that the predicted RUL is coinciding fairly well with the real RUL, which indicates that the feature fusion method and constructed ELM model can achieve an ideal result. However, the RUL prediction accuracy is relatively high in the first stage as the vibration is relatively stable. The slight fault stage is always complex and sharp fluctuation also happens at this stage. The vibration at the serious fault stage is more regular, which can lead to higher prediction accuracy.Figure 13
Predicted RUL: (a) outer-race-fault case, (b) inner-race-fault case, and (c) roller-fault case.
(a)
(b)
(c)The PCA and ICA as the representative linear dimensionality reduction methods are applied on the same data to compare with the proposed method. The root mean square error (RMSE) is selected to evaluate the RUL prediction accuracy which can be defined as(18)
R
M
S
E
=
∑
i
=
1
n
t
i
-
y
i
2
,
i
=
1,2
,
…
,
n
,where t
i is the real RUL and y
i is the predicted RUL.Considering the problem of repeatability, the ELM is applied on the same training data 10 times, and the average of the RMSE is selected to evaluate the prediction accuracy. The results of comparison are shown in Tables2–4. It can be seen that the proposed feature extraction is better than the traditional two methods, which proves the merits of the proposed method for fusing original features.Table 2
Predicted accuracy of outer-race-fault case.
Training
Testing
Number of training data
RMSE
Number of testing data
RMSE
PCA
100
12.549
100
13.526
ICA
100
12.605
100
14.108
JADE
100
8.872
100
10.018Table 3
Predicted accuracy of inner-race-fault case.
Training
Testing
Number of training data
RMSE
Number of testing data
RMSE
PCA
100
16.073
100
17.548
ICA
100
16.165
100
17.103
JADE
100
14.552
100
14.846Table 4
Predicted accuracy of roller-fault case.
Training
Testing
Number of training data
RMSE
Number of testing data
RMSE
PCA
100
9.668
100
9.124
ICA
100
9.650
100
9.077
JADE
100
9.105
100
8.227


## 5. Conclusions

A novel feature extraction method that integrates the PSR,S
S, and JADE is introduced in this paper. The contribution of the original features is considered according to S
S that derived from phase point distribution after PSR. The JADE is employed to fuse the features. Based on the extracted advanced features, the ELM is employed to construct the RUL prediction model.A set of experimental case studies are presented to verify the effectiveness of the proposed method. The results indicate that the proposed method is appropriate for RUL prediction and even performs better than the traditional methods. The introduced JADE joint PSR andS
S feature extraction method is a multifeature fusion method that can inhibit the nonsensitive features while enhancing the superior features. It is also suitable for other machine RUL predictions, such as cutting tools and gears.



---

*Source: 8623530-2016-04-18.xml*