# Phase Space Similarity as a Signature for Rolling Bearing Fault Diagnosis and Remaining Useful Life Estimation

Authors: Fang Liu; Bing He; Yongbin Liu; Siliang Lu; Yilei Zhao; Jiwen Zhao
Journal: Shock and Vibration
 (2016)

---

## Abstract

Feature extraction from vibration signal is still a challenge in the area of fault diagnosis and remaining useful life (RUL) estimation of rotary machine. In this paper, a novel feature called phase space similarity (PSS) is introduced for health condition monitoring of bearings. Firstly, the acquired signal is transformed to the phase space through the phase space reconstruction (PSR). The similar vibration always exists in the phase space due to the comparable evolution of the dynamics that are characteristic of the system state. Secondly, the normalized cross-correlation (NCC) is employed to calculate the PSS between bearing data with different states. Based on the PSS, a fault pattern recognition algorithm, a bearing fault size prediction algorithm, and a RUL estimation algorithm are introduced to analyze the experimental signal. Results have shown the effectiveness of the PSS as it can better grasp the nature and regularity of the signals.


---

## Body

## 1. Introduction

Currently, fault diagnosis is an extremely significant research field in industry; intelligent health monitoring is gradually supplemented in maintenance of the machine to ensure the stability of systems and decrease the downtime. Rolling bearing is a key component in the rotary machine, and faults occurred frequently, leading to a disastrous accident. Typical damage of bearing is located at outer race, ball, and inner race. One of the challenges for health monitoring is sensitive feature extraction from sensor signals that are often mixed by noise [1, 2]. Representative feature extraction methods are based on time-domain analysis and frequency-domain analysis [3–5]. Although the above methods cannot extract the nonlinear features from the vibration signals, the time-frequency analysis can overcome this problem, investigating a signal in both time and frequency domains [6]. In addition, empirical mode decomposition and wavelet analysis are able to express the nonstationary signals well and can also extract the sensitive features [7, 8]. Nevertheless, the extracted features are still redundant and with high dimension. Therefore, the advanced dimension reduction techniques are proposed to solve the aforementioned problem such as principal component analysis (PCA) [9] and independent component analysis (ICA) [10], and they have been employed to extract sensitive features from high-dimensional data.Chaotic time series analysis is an excellent method for nonlinear signals, which can grasp the nature and regularity of the time series [11]. And it has been verified in signal processing. In the literature [12], vibration signal is reconstructed by PSR in which the phase trajectory of different kinds of signals will demonstrate different structures. In the literature [13], the PSR is combined with time-frequency synthesis and utilized in data denoising. In these studies, the PSR is introduced for fault diagnosis and RUL estimation, which can extend time series to a high-dimensional phase space and maintain topologically equivalent. Theoretically speaking, a time series can sufficiently reconstruct an original dynamic system by PSR [14].On the one hand, the template matching is widely applied to image analysis, which can calculate the dissimilarity or similarity between the template and the image. Then the position with maximum correlation is taken to represent the instance of the template into the image [15]. Classical correlation matching methods are as below: maximum close distance (MCD), minimum absolute difference (MAD), and NCC [16]. Among them, the NCC has been widely used in correlation measure because it proved to be good at template matching [17, 18]. In these studies, NCC is utilized in sensor signal processing. Besides, the NCC is combined with PSR in order to serve for feature extraction. The new feature was called PSS. Furthermore, template matching that is based on the proposed feature is employed in the recognition of bearing fault patterns and estimation of bearing RUL.On the other hand, on the basis of these features that are extracted from the sensor signals, machine learning is widely used in machine fault diagnosis such as artificial neural networks (ANN) [19] and support vector machine (SVM) [20]. Although these strategies have been successfully utilized in fault diagnosis, the ANN still has a slow convergence and poor generalization ability and the SVM is only being used to handle a binary problem. What is more, the aforementioned methods are inefficient or even infeasible in regression problem. According to the study of the literatures [21–23], the support vector regression machine (SVR) can establish a stable nonlinear relationship between inputs and outputs and then make time series prediction. As a result, the SVR is introduced for bearing fault size prediction so as to overcome the above limitation.In this research, firstly, the autocorrelation function is used to determine the delay time. Meanwhile, Cao’s algorithm is used to determine the embedding dimension [24], and then all of data are transformed to the phase space by PSR. Thereafter, a sliding window is employed to obtain the data segment and NCC is adopted in template matching to realize the recognition of bearing fault patterns. Secondly, a SVR structure is trained by PSS, which aims to quantitatively predict the fault size. Finally, the RUL of the tested bearing is estimated by template matching that is based on PSS. Going further, three applications based on PSS are verified by case studies; the experimental results show that the proposed method can effectively realize the fault diagnosis and RUL estimation for bearing.The rest of the paper is organized as follows. Section2 gives a brief introduction of PSR, NCC, and SVR. The proposed fault pattern recognition, fault size prediction, and RUL estimation procedures for bearings are illustrated in Section 3. Experimental validations and the comparison with original signal’s similarity and root mean square (RMS) are presented in Section 4. Finally, the conclusions are drawn in Section 5.


## 2. Theory Background

### 2.1. Phase Space Reconstitution (PSR)

The research of PSR on chaotic time series was started by Packard et al. [25]. The underlying dynamics characteristic of a system can be exhibited in the phase space. In this research, the time delay embedding is applied to characterize the nonlinear feature of time series. Suppose that the time series is x
=
(
x
1
,
x
2
,
…
,
x
n
), and 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 points in time series, i is the index of the phase point, which ranges from 1 to n
-
(
d
-
1
)
τ, the embedding dimension is chosen as d, and τ is the time delay. According to the literature [26], the embedding dimension and the time delay of the ideal time series (infinite length and noiseless) can be chosen as any value. However, the sensor signals are distributed with ambient noise. Therefore, it is difficult to determine the parameters of PSR [27].The adjacent components in the phase space will became similar when the time delay is selected as a low value. By contrast, the components will be independent in statistics with a high value of time delay. Consequently, the autocorrelation function is utilized in this research for time delay identification; it 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 [28], the minimum value of time delay τ that firstly makes C
i
(
τ
) lower than 1
-
1
/
e can be considered as the most optimal time delay, where e is the natural constant that is approximately equal to 2.718.After that, Cao’s algorithm is employed to determine the embedding dimension. In accordance with (1), the ith phase point is written as X
i
(
d
); suppose that the nearest phase point is X
near
(
i
,
d
)
(
d
). In the light of the embedding theorem, if d is qualified as embedding dimension, then any two points that are close in the d dimensional reconstructed phase space will continue to be close in the d
+
1 dimensional space:(3)
a
i
,
d
=
X
i
d
+
1
-
X
near
i
,
d
+
1
d
+
1
X
i
d
-
X
near
i
,
d
d
,
E
d
=
1
n
-
d
τ
∑
i
=
1
n
-
d
τ
a
i
,
d
,
E
1
d
=
E
d
+
1
E
d
,where a
(
i
,
d
) indicates the relative distance increment and E
(
d
) indicates the mean statistics of the a
(
i
,
d
). E
1
(
d
) will be gradually stabilized with the increase of d; therefore, the minimum d that stabilizes E
1
(
d
) is the embedding dimension. However, it is difficult to judge that E
1
(
d
) has been steadied owing to the finite length of obtained samples. Therefore, supplementary criteria are proposed as follows:(4)
E
∗
d
=
1
n
-
d
τ
∑
i
=
1
n
-
d
τ
x
i
+
d
τ
-
x
near
i
,
d
+
d
τ
,
E
2
d
=
E
∗
d
+
1
E
∗
d
.Suppose that the time series is random,E
2
(
d
) is equal to 1. In contrast, as to a deterministic sample, E
2
(
d
) will be changed with the increasing of d. The appropriate embedding dimension can be determined by observing the curve of both E
1
(
d
) and E
2
(
d
) [29].


### 2.2. Normalized Cross-Correlation (NCC)

Classical correlation matching methods are as follows: MCD, MAD, and NCC. In this research, the NCC is employed to calculate the similarity because it is invariant to a linear transformation and is less sensitive to noise [16]. Suppose that the sample phase space is Z
=
Z
1
,
Z
2
,
…
,
Z
m and the current phase space is Y
=
Y
1
,
Y
2
,
…
,
Y
m, where m is the number of phase points. The PSS is expressed as (5)
PSS
i
=
Z
-
Z
-
∗
Y
i
-
Y
-
i
′
Z
-
Z
-
2
Y
i
-
Y
-
i
2
,where Y
i
=
Y
i
,
Y
i
+
1
,
…
,
Y
i
+
l is a part of Y and it is obtained by a sliding window with length l and Z
- and Y
-
i are the mean vectors of Z and Y
i.


### 2.3. Support Vector Regression Machine (SVR)

SVR was developed for solving regression problems, which maximizes predictive accuracy and avoids overfitting simultaneously [30]. Suppose that the training sample is x
i
,
y
i, i
=
1,2
,
…
,
n, where x
i is the input features and y
i is the label of experiment condition. It can be represented by a linear function of f
(
x
)
=
w
∗
x
+
b, where w identifies the weight vector and b corresponds to the threshold coefficient. The SVR function can be more comfortable by using the ε-insensitive loss function as follows: (6)
y
i
-
w
·
x
i
-
b
≤
ε
,
w
·
x
i
+
b
-
y
i
≤
ε
,
i
=
1
,
…
,
n
.Flatness of linear function can be ensured by minimizing the1
/
2
w
2, leading to the following convex optimization problem: (7)
min
⁡
1
2
w
2
+
C
∑
i
=
1
n
ξ
i
+
ξ
i
∗
,
s.t.
y
i
-
w
·
x
i
-
b
≤
ε
+
ξ
i
,
w
·
x
i
+
b
-
y
i
≤
ε
+
ξ
i
∗
,
i
=
1
,
…
,
nwith regularization parameter C
>
0 and the slack variables ξ
i
≥
0 and ξ
i
∗
≥
0. Then we get the following convex quadratic programming problem by introducing a Lagrange function:(8)
L
=
1
2
w
2
+
C
∑
i
=
1
r
ξ
i
+
ξ
i
∗
-
∑
i
=
1
r
a
i
ε
+
ξ
i
+
y
i
-
w
·
x
i
-
b
-
∑
i
=
1
r
a
i
∗
ε
+
ξ
i
∗
+
y
i
-
w
·
x
i
-
b
-
∑
i
=
1
r
η
i
ξ
i
+
η
i
∗
ξ
i
∗
,where a
i and a
i
∗ are the Lagrange multipliers. The function can be expressed as follows:(9)
max
⁡
w
a
,
a
∗
=
-
ε
∑
i
=
1
n
y
i
a
i
∗
-
a
i
-
1
2
∑
i
,
j
=
1
n
a
i
∗
-
a
i
a
j
∗
-
a
j
x
i
·
x
j
,
s.t.
∑
i
=
1
n
a
i
-
a
i
∗
=
0
,
0
≤
a
i
,
a
i
∗
≤
C
,
i
=
1
,
…
,
n
.The regression function is(10)
f
x
=
w
·
x
+
b
=
∑
i
=
1
n
a
i
∗
-
a
i
x
i
·
x
+
b
∗
.


## 2.1. Phase Space Reconstitution (PSR)

The research of PSR on chaotic time series was started by Packard et al. [25]. The underlying dynamics characteristic of a system can be exhibited in the phase space. In this research, the time delay embedding is applied to characterize the nonlinear feature of time series. Suppose that the time series is x
=
(
x
1
,
x
2
,
…
,
x
n
), and 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 points in time series, i is the index of the phase point, which ranges from 1 to n
-
(
d
-
1
)
τ, the embedding dimension is chosen as d, and τ is the time delay. According to the literature [26], the embedding dimension and the time delay of the ideal time series (infinite length and noiseless) can be chosen as any value. However, the sensor signals are distributed with ambient noise. Therefore, it is difficult to determine the parameters of PSR [27].The adjacent components in the phase space will became similar when the time delay is selected as a low value. By contrast, the components will be independent in statistics with a high value of time delay. Consequently, the autocorrelation function is utilized in this research for time delay identification; it 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 [28], the minimum value of time delay τ that firstly makes C
i
(
τ
) lower than 1
-
1
/
e can be considered as the most optimal time delay, where e is the natural constant that is approximately equal to 2.718.After that, Cao’s algorithm is employed to determine the embedding dimension. In accordance with (1), the ith phase point is written as X
i
(
d
); suppose that the nearest phase point is X
near
(
i
,
d
)
(
d
). In the light of the embedding theorem, if d is qualified as embedding dimension, then any two points that are close in the d dimensional reconstructed phase space will continue to be close in the d
+
1 dimensional space:(3)
a
i
,
d
=
X
i
d
+
1
-
X
near
i
,
d
+
1
d
+
1
X
i
d
-
X
near
i
,
d
d
,
E
d
=
1
n
-
d
τ
∑
i
=
1
n
-
d
τ
a
i
,
d
,
E
1
d
=
E
d
+
1
E
d
,where a
(
i
,
d
) indicates the relative distance increment and E
(
d
) indicates the mean statistics of the a
(
i
,
d
). E
1
(
d
) will be gradually stabilized with the increase of d; therefore, the minimum d that stabilizes E
1
(
d
) is the embedding dimension. However, it is difficult to judge that E
1
(
d
) has been steadied owing to the finite length of obtained samples. Therefore, supplementary criteria are proposed as follows:(4)
E
∗
d
=
1
n
-
d
τ
∑
i
=
1
n
-
d
τ
x
i
+
d
τ
-
x
near
i
,
d
+
d
τ
,
E
2
d
=
E
∗
d
+
1
E
∗
d
.Suppose that the time series is random,E
2
(
d
) is equal to 1. In contrast, as to a deterministic sample, E
2
(
d
) will be changed with the increasing of d. The appropriate embedding dimension can be determined by observing the curve of both E
1
(
d
) and E
2
(
d
) [29].


## 2.2. Normalized Cross-Correlation (NCC)

Classical correlation matching methods are as follows: MCD, MAD, and NCC. In this research, the NCC is employed to calculate the similarity because it is invariant to a linear transformation and is less sensitive to noise [16]. Suppose that the sample phase space is Z
=
Z
1
,
Z
2
,
…
,
Z
m and the current phase space is Y
=
Y
1
,
Y
2
,
…
,
Y
m, where m is the number of phase points. The PSS is expressed as (5)
PSS
i
=
Z
-
Z
-
∗
Y
i
-
Y
-
i
′
Z
-
Z
-
2
Y
i
-
Y
-
i
2
,where Y
i
=
Y
i
,
Y
i
+
1
,
…
,
Y
i
+
l is a part of Y and it is obtained by a sliding window with length l and Z
- and Y
-
i are the mean vectors of Z and Y
i.


## 2.3. Support Vector Regression Machine (SVR)

SVR was developed for solving regression problems, which maximizes predictive accuracy and avoids overfitting simultaneously [30]. Suppose that the training sample is x
i
,
y
i, i
=
1,2
,
…
,
n, where x
i is the input features and y
i is the label of experiment condition. It can be represented by a linear function of f
(
x
)
=
w
∗
x
+
b, where w identifies the weight vector and b corresponds to the threshold coefficient. The SVR function can be more comfortable by using the ε-insensitive loss function as follows: (6)
y
i
-
w
·
x
i
-
b
≤
ε
,
w
·
x
i
+
b
-
y
i
≤
ε
,
i
=
1
,
…
,
n
.Flatness of linear function can be ensured by minimizing the1
/
2
w
2, leading to the following convex optimization problem: (7)
min
⁡
1
2
w
2
+
C
∑
i
=
1
n
ξ
i
+
ξ
i
∗
,
s.t.
y
i
-
w
·
x
i
-
b
≤
ε
+
ξ
i
,
w
·
x
i
+
b
-
y
i
≤
ε
+
ξ
i
∗
,
i
=
1
,
…
,
nwith regularization parameter C
>
0 and the slack variables ξ
i
≥
0 and ξ
i
∗
≥
0. Then we get the following convex quadratic programming problem by introducing a Lagrange function:(8)
L
=
1
2
w
2
+
C
∑
i
=
1
r
ξ
i
+
ξ
i
∗
-
∑
i
=
1
r
a
i
ε
+
ξ
i
+
y
i
-
w
·
x
i
-
b
-
∑
i
=
1
r
a
i
∗
ε
+
ξ
i
∗
+
y
i
-
w
·
x
i
-
b
-
∑
i
=
1
r
η
i
ξ
i
+
η
i
∗
ξ
i
∗
,where a
i and a
i
∗ are the Lagrange multipliers. The function can be expressed as follows:(9)
max
⁡
w
a
,
a
∗
=
-
ε
∑
i
=
1
n
y
i
a
i
∗
-
a
i
-
1
2
∑
i
,
j
=
1
n
a
i
∗
-
a
i
a
j
∗
-
a
j
x
i
·
x
j
,
s.t.
∑
i
=
1
n
a
i
-
a
i
∗
=
0
,
0
≤
a
i
,
a
i
∗
≤
C
,
i
=
1
,
…
,
n
.The regression function is(10)
f
x
=
w
·
x
+
b
=
∑
i
=
1
n
a
i
∗
-
a
i
x
i
·
x
+
b
∗
.


## 3. Methods and Procedures

In this paper, a novel feature extraction method that is based on PSR and NCC is proposed for fault diagnosis and RUL estimation. The framework of the proposed method is shown in Figure1, where the signals come from bearing by accelerometer sensor. In the step of data reconstruction, PSR is applied to embed the vibration signals into the high-dimensional phase space. Assumed that the jth sample signal is x
j
=
(
x
1
j
,
x
2
j
,
…
,
x
n
j
), where n is the number of signal points and j ranges from 1 to M. Then the autocorrelation function is employed to determine the delay time; Cao’s algorithm is employed to determine the embedding dimension d
j in jth sample of database. To simplify, the embedding dimension is unified to d, where d
>
d
j, and the jth sample phase space is described in (11)
Z
j
=
X
i
j
,
X
i
+
1
j
,
…
,
X
n
-
d
τ
j
T
,where X
i
j
=
[
x
i
j
,
x
i
+
τ
j
,
…
,
x
i
+
(
d
-
2
)
τ
j
,
x
i
+
(
d
-
1
)
τ
j
] is the ith phase point. Suppose that y
=
(
y
1
,
y
2
,
…
,
y
n
) is the current data; in order to reduce the influence of the phase offset, the sliding window is utilized to obtain the segment from current data, where the window width is C and step length is S. The segment is written as y
k
=
(
y
i
,
y
i
+
1
,
…
,
y
i
+
C
-
1
) and the corresponding phase space is given in(12)
Y
k
=
Y
i
,
Y
i
+
1
,
…
,
Y
i
+
C
-
1
-
d
τ
T
,where Y
i
=
y
i
k
,
y
i
+
τ
k
,
…
,
y
i
+
(
d
-
2
)
τ
k
,
y
i
+
(
d
-
1
)
τ
k is ith phase point, k changes from 1 to (
n
-
C
)
/
S
+
1, and i changes from 1 to n
-
C
+
1. Thereafter, the PSS between Y
k and Z
j is calculated in accordance with (5); suppose that PSS-j is the maximum value between current data and jth sample data. Finally, the PSS is employed to health condition monitoring and it can contribute to the recognition of bearing fault patterns, prediction of fault size, and RUL estimation. The detailed description of the three applications is described as follows.Figure 1
Framework of the proposed method.(
1) Recognition of Bearing Fault Patterns. Taking all fault patterns signal as database, the maximum PSS from PSS-j will be found and the corresponding fault pattern is considered as the current fault pattern. The framework of the proposed method is shown in Figure 2.Figure 2
Framework of the fault pattern recognition.(
2) Prediction of Fault Size. As described in Figure 3, take normal vibration signals as database and bearing vibration signal with different fault size as the train data. The PSS between train data and normal signal as the training sample is to train the SVR structure, and the PSS between test data and normal signal as the testing sample is used for the trained SVR structure. Finally, the fault size of test data can be predicted by SVR.Figure 3
Framework of the fault size prediction.(
3) RUL Estimation of Rolling Bearing. The framework of the proposed method is shown in Figure 4. Take M vibration signals of the entire life cycle as database, which are acquired from bearing run-to-failure and chronological arrangement. Due to the random fluctuation of signals, the accuracy of template matching will be decreased. Therefore, the similarity metric is improved as follows: (13)
S
K
j
=
PSS
j
+
1
-
t
c
-
t
j
t
c
,where t
c is the current operating time. In this similarity metric, the factor of operating time has been considered. If S
K gets the maximum value at t
j, where the current health condition is most similar to the jth segment, then t
j is assigned to T. The RUL of the current data can be calculated as(14)
RUL
=
t
M
-
T
,where t
M denotes the failure time or lifetime that has been recorded in tested bearing sample.Figure 4
Framework of the RUL estimation.


## 4. Experiment Validation

### 4.1. Validation of Fault Patterns Recognition

The rolling bearing fault data of the Case Western Reserve University [31] are used to verify the proposed method; the vibration data are recorded on the test bench given in Figure 5. The test platform consists of four parts that are a motor, a dynamometer, a torque transducer, and a control electrometer. The resulting vibration was measured by an accelerometer being mounted to the motor housing at the fan end of the motor. Single-point faults of size 0.007, 0.014, and 0.021 inches were set on the fan end bearings by using the electric discharge machining approach. These faults were set on ball, inner race, and outer race. The motor speed was 1797 r/min, the sampling frequency is 12 kHz, and motor load is 0.Figure 5
The experimental bench.All fault locations and sizes are listed in Table1, which includes 10 patterns; every pattern contains 122426 points, and then divided into 60 samples, each sample has 512 points and is without overlap between adjacent samples; the 60 samples are divided into 30 training samples and 30 testing samples. Therefore, this is a complex classification problem of ten patterns, which requires comprehensive identifying of both the fault location and the fault size.Table 1
The description of the bearing data.
Condition
Fault size (inch)
Number of samples (train/test)
Label of pattern
Normal
0
30/30
1
Inner race
0.007
30/30
2
0.014
30/30
3
0.021
30/30
4
Outer race
0.007
30/30
5
0.014
30/30
6
0.021
30/30
7
Ball
0.007
30/30
8
0.014
30/30
9
0.021
30/30
10At the very beginning, the embedding dimensiond and the time delay τ should be determined as described in Section 2.1. Take inner-0.014-inch-fault case as an example, the curves of E
1
(
d
), E
2
(
d
), and C
i
(
τ
) are shown in Figure 6.Figure 6
(a) Curve ofE
1
(
d
) and E
2
(
d
) and (b) curve of C
i
(
τ
).
(a)
(b)It can be seen thatE
1
(
d
) tends to a stationary value when d is larger than 8 and E
2
(
d
) changes with the increasing of d. The same conclusions are drawn for other cases. Hence, the embedding dimension is set as 8. And C
i
(
τ
) is lower than 1
-
1
/
e when τ equals 2; hence, the time delay is set as 2.Thereafter, a sliding window with a window width of 64 and step length of 16 is used to intercept the segments from the current data. Take inner-0.014-inch-fault case as current data, the maximum PSS between current data and every fault pattern is calculated and the results are shown in Figure7.Figure 7
The PSS between the current data and every fault pattern.It can be found that PSS of pattern 3 is generally higher than others; therefore, the fault patterns of current data are most similar with the health condition of pattern 3. Thus, the current fault pattern is considered as inner-0.014-inch-fault pattern. Afterwards, all of testing samples are processed as above-mentioned methods that are based on PSS. In order to demonstrate the improvements, the proposed method is compared with the recognition method that is based on original signal similarity and SVM such that the classical feature RMS is selected as input. The results are shown in Table2.Table 2
Results of rolling bearing fault pattern recognition.
Fault pattern
Original signal
RMS
Phase space
Number of samples(test/true)
Accuracy
Number of samples(test/true)
Accuracy
Number of samples (test/true)
Accuracy
1
30/30
100%
30/1
3.3%
30/30
100%
2
30/30
100%
30/10
33.3%
30/30
100%
3
30/28
93.3%
30/13
43.3%
30/30
100%
4
30/29
96.7%
30/14
46.7%
30/30
100%
5
30/17
56.7%
30/25
83.3%
30/20
66.7%
6
30/16
53.3%
30/30
100%
30/18
60%
7
30/18
60%
30/26
86.7%
30/23
76.7%
8
30/26
86.7%
30/19
63.3%
30/30
100%
9
30/27
90%
30/19
63.3%
30/30
100%
10
30/30
100%
30/6
20%
30/30
100%
Total
300/251
83.7%
300/162
54%
300/271
90.3%From Table2, it can be seen that the recognition of fault patterns based on PSS has a higher accuracy. Moreover, the accuracy of the ball-fault (pattern 5 to pattern 7) recognition is lower than other fault patterns because the vibration signal of defective ball is more complicated and with strong randomization.


### 4.2. Prediction of the Fault Size

To evaluate the effectiveness of the proposed feature for fault size prediction, the SVR is adopted in this section. The PSS based on two-class mode is proposed to predict the fault size. At first, take normal bearing vibration signal as one class and bearing vibration signal with different fault sizes as another class. The results of PSS are shown in Figure8.Figure 8
The PSS between two fault sizes.
(a)
(b)
(c)From Figure8, the classes 1–4 are represented as fault sizes 0, 0.007, 0.014, and 0.021 inches, respectively. The height of the bar is the PSS between corresponding fault size signal and normal signal. It can be seen that the PSS gradually decreases with the increasing of the fault size. After that, consider the PSS of classes 1, 2, and 4 as the training sample to train the SVR structure and the PSS of class 3 as the testing sample. In this study, the kernel function of SVR is set as RBF kernel K
(
x
i
,
x
)
=
exp
⁡
(
-
x
i
-
x
2
/
2
σ
2
), where x and x
i are represented as feature vectors, x
i
-
x
2 can be recognized as the squared Euclidean distance between the two feature vectors, and σ is a free parameter. Finally, the fault size of the testing sample will be regression-analyzed and predicted by SVR. The results of prediction are presented in Figure 9.Figure 9
The prediction of fault size based on SVR.
(a)
Inner race
(b)
Ball
(c)
Outer raceThe maximum error is utilized to determine the accuracy of the fault size prediction and is defined asmax
⁡
(
f
predicted
-
f
actual
); the error is 0.043, 0.049, and 0.019 inches, respectively.The values of RMS of three bearing fault cases are shown in Figure10. It can be seen that the values of RMS are nonmonotonic with growth of fault size, so it was not suitable in fault size prediction.Figure 10
The RMS of three bearing fault cases.
(a)
(b)
(c)


### 4.3. RUL Estimation of Rolling Bearing

To evaluate the effectiveness of the proposed feature for RUL estimation, the entire life cycle data of bearing is originated from the Center for Intelligent Maintenance Systems [32]. The experimental data sets are generated from run-to-failure bearing, which is tested on the designed test platform. 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 at 20 kHz, the vibration data are recorded on the test bench given in Figure 11, and more details can be found in [33].Figure 11
Bearing test platform and sensor placement illustration.
(a)
(b)Taking outer-race-fault case as an example, the curves ofE
1
(
d
), E
2
(
d
), and C
i
(
τ
) are shown in Figure 12. From Figure 12(a), it can be seen that E
1
(
d
) keeps stable when d is larger than 8, and E
2
(
d
) randomly fluctuates with d. Figure 12(b) shows 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. Consequently, the time delay is set to 2 and embedding dimension is set to 10. Then, the entire life data is reconstructed by PSR, and the similarity between the two adjacent data is calculated by NCC and shown in Figure 13.Figure 12
(a) Curve ofE
1
(
d
) and E
2
(
d
) and (b) curve of C
i
(
τ
).
(a)
(b)Figure 13
Bearing performance assessment based on PSS.There are three health states: normal operation stage, slight degradation stage, and serious degradation stage. It can be seen that the value of PSS becomes higher in the serious degradation stage. The reason is that the vibration becomes more regular when the failure occurs. PSS can better reveal the potential links between adjacent signals, which have been proved in the front part of this research. At the same time, PSS begins to increase in the region of 500-phase point, which can be considered as the threshold for the beginning of the RUL estimation. After that, the original data are divided into reference data and test data; reference data are viewed as the database; test data are considered as the current data to estimate the RUL. The results are shown in Figures14–16.Figure 14
RUL estimation for the outer-race-fault case based on PSS.Figure 15
RUL estimation for the inner-race-fault case based on PSS.Figure 16
RUL estimation for the ball-fault case based on PSS.It can be seen that the estimated RUL is close to the actual RUL. Furthermore, the comparisons with the original signal similarity and RMS are conducted to prove the superiority of the proposed method in RUL estimation. The results of estimation based on original signal similarity are shown in Figures17–19. Three criteria defined in the following equation are applied to evaluate the accuracy of the RUL estimation: (15)
F
1
=
1
N
∑
i
=
1
N
x
i
-
y
i
,
F
2
=
1
N
∑
i
=
1
N
x
i
-
y
i
2
,
F
3
=
1
N
∑
i
=
1
N
x
i
-
y
i
2
,where x
i is the estimated RUL, y
i is the actual RUL, and N is the length of the data points. The estimation results of the last 15 h are shown in Table 3.Table 3
Comparison of RUL estimation.
F
1
F
2
F
3
Outer race
RMS
6.0578
8.5576
5.1652
Original signal
1.5323
2.0537
1.0665
Phase space
1.3226
1.9092
0.8540
Inner race
RMS
10.8441
12.5821
8.2846
Original signal
4.1250
5.4658
3.0780
Phase space
3.6250
4.8606
2.7941
Ball
RMS
8.5964
10.2489
7.0659
Original signal
3.1875
4.1155
2.4702
Phase space
3
3.7583
2.3632Figure 17
RUL estimation for the outer-race-fault case based on original signal similarity.Figure 18
RUL estimation for the inner-race-fault case based on original signal similarity.Figure 19
RUL estimation for the ball-fault case based on original signal similarity.As for the findings in Table3, the proposed method has a higher accuracy, which indicates that the PSR can better grasp the nature and regularity of signal. It is effective in the RUL estimation of bearing, especially in the serious degradation stage.


## 4.1. Validation of Fault Patterns Recognition

The rolling bearing fault data of the Case Western Reserve University [31] are used to verify the proposed method; the vibration data are recorded on the test bench given in Figure 5. The test platform consists of four parts that are a motor, a dynamometer, a torque transducer, and a control electrometer. The resulting vibration was measured by an accelerometer being mounted to the motor housing at the fan end of the motor. Single-point faults of size 0.007, 0.014, and 0.021 inches were set on the fan end bearings by using the electric discharge machining approach. These faults were set on ball, inner race, and outer race. The motor speed was 1797 r/min, the sampling frequency is 12 kHz, and motor load is 0.Figure 5
The experimental bench.All fault locations and sizes are listed in Table1, which includes 10 patterns; every pattern contains 122426 points, and then divided into 60 samples, each sample has 512 points and is without overlap between adjacent samples; the 60 samples are divided into 30 training samples and 30 testing samples. Therefore, this is a complex classification problem of ten patterns, which requires comprehensive identifying of both the fault location and the fault size.Table 1
The description of the bearing data.
Condition
Fault size (inch)
Number of samples (train/test)
Label of pattern
Normal
0
30/30
1
Inner race
0.007
30/30
2
0.014
30/30
3
0.021
30/30
4
Outer race
0.007
30/30
5
0.014
30/30
6
0.021
30/30
7
Ball
0.007
30/30
8
0.014
30/30
9
0.021
30/30
10At the very beginning, the embedding dimensiond and the time delay τ should be determined as described in Section 2.1. Take inner-0.014-inch-fault case as an example, the curves of E
1
(
d
), E
2
(
d
), and C
i
(
τ
) are shown in Figure 6.Figure 6
(a) Curve ofE
1
(
d
) and E
2
(
d
) and (b) curve of C
i
(
τ
).
(a)
(b)It can be seen thatE
1
(
d
) tends to a stationary value when d is larger than 8 and E
2
(
d
) changes with the increasing of d. The same conclusions are drawn for other cases. Hence, the embedding dimension is set as 8. And C
i
(
τ
) is lower than 1
-
1
/
e when τ equals 2; hence, the time delay is set as 2.Thereafter, a sliding window with a window width of 64 and step length of 16 is used to intercept the segments from the current data. Take inner-0.014-inch-fault case as current data, the maximum PSS between current data and every fault pattern is calculated and the results are shown in Figure7.Figure 7
The PSS between the current data and every fault pattern.It can be found that PSS of pattern 3 is generally higher than others; therefore, the fault patterns of current data are most similar with the health condition of pattern 3. Thus, the current fault pattern is considered as inner-0.014-inch-fault pattern. Afterwards, all of testing samples are processed as above-mentioned methods that are based on PSS. In order to demonstrate the improvements, the proposed method is compared with the recognition method that is based on original signal similarity and SVM such that the classical feature RMS is selected as input. The results are shown in Table2.Table 2
Results of rolling bearing fault pattern recognition.
Fault pattern
Original signal
RMS
Phase space
Number of samples(test/true)
Accuracy
Number of samples(test/true)
Accuracy
Number of samples (test/true)
Accuracy
1
30/30
100%
30/1
3.3%
30/30
100%
2
30/30
100%
30/10
33.3%
30/30
100%
3
30/28
93.3%
30/13
43.3%
30/30
100%
4
30/29
96.7%
30/14
46.7%
30/30
100%
5
30/17
56.7%
30/25
83.3%
30/20
66.7%
6
30/16
53.3%
30/30
100%
30/18
60%
7
30/18
60%
30/26
86.7%
30/23
76.7%
8
30/26
86.7%
30/19
63.3%
30/30
100%
9
30/27
90%
30/19
63.3%
30/30
100%
10
30/30
100%
30/6
20%
30/30
100%
Total
300/251
83.7%
300/162
54%
300/271
90.3%From Table2, it can be seen that the recognition of fault patterns based on PSS has a higher accuracy. Moreover, the accuracy of the ball-fault (pattern 5 to pattern 7) recognition is lower than other fault patterns because the vibration signal of defective ball is more complicated and with strong randomization.


## 4.2. Prediction of the Fault Size

To evaluate the effectiveness of the proposed feature for fault size prediction, the SVR is adopted in this section. The PSS based on two-class mode is proposed to predict the fault size. At first, take normal bearing vibration signal as one class and bearing vibration signal with different fault sizes as another class. The results of PSS are shown in Figure8.Figure 8
The PSS between two fault sizes.
(a)
(b)
(c)From Figure8, the classes 1–4 are represented as fault sizes 0, 0.007, 0.014, and 0.021 inches, respectively. The height of the bar is the PSS between corresponding fault size signal and normal signal. It can be seen that the PSS gradually decreases with the increasing of the fault size. After that, consider the PSS of classes 1, 2, and 4 as the training sample to train the SVR structure and the PSS of class 3 as the testing sample. In this study, the kernel function of SVR is set as RBF kernel K
(
x
i
,
x
)
=
exp
⁡
(
-
x
i
-
x
2
/
2
σ
2
), where x and x
i are represented as feature vectors, x
i
-
x
2 can be recognized as the squared Euclidean distance between the two feature vectors, and σ is a free parameter. Finally, the fault size of the testing sample will be regression-analyzed and predicted by SVR. The results of prediction are presented in Figure 9.Figure 9
The prediction of fault size based on SVR.
(a)
Inner race
(b)
Ball
(c)
Outer raceThe maximum error is utilized to determine the accuracy of the fault size prediction and is defined asmax
⁡
(
f
predicted
-
f
actual
); the error is 0.043, 0.049, and 0.019 inches, respectively.The values of RMS of three bearing fault cases are shown in Figure10. It can be seen that the values of RMS are nonmonotonic with growth of fault size, so it was not suitable in fault size prediction.Figure 10
The RMS of three bearing fault cases.
(a)
(b)
(c)


## 4.3. RUL Estimation of Rolling Bearing

To evaluate the effectiveness of the proposed feature for RUL estimation, the entire life cycle data of bearing is originated from the Center for Intelligent Maintenance Systems [32]. The experimental data sets are generated from run-to-failure bearing, which is tested on the designed test platform. 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 at 20 kHz, the vibration data are recorded on the test bench given in Figure 11, and more details can be found in [33].Figure 11
Bearing test platform and sensor placement illustration.
(a)
(b)Taking outer-race-fault case as an example, the curves ofE
1
(
d
), E
2
(
d
), and C
i
(
τ
) are shown in Figure 12. From Figure 12(a), it can be seen that E
1
(
d
) keeps stable when d is larger than 8, and E
2
(
d
) randomly fluctuates with d. Figure 12(b) shows 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. Consequently, the time delay is set to 2 and embedding dimension is set to 10. Then, the entire life data is reconstructed by PSR, and the similarity between the two adjacent data is calculated by NCC and shown in Figure 13.Figure 12
(a) Curve ofE
1
(
d
) and E
2
(
d
) and (b) curve of C
i
(
τ
).
(a)
(b)Figure 13
Bearing performance assessment based on PSS.There are three health states: normal operation stage, slight degradation stage, and serious degradation stage. It can be seen that the value of PSS becomes higher in the serious degradation stage. The reason is that the vibration becomes more regular when the failure occurs. PSS can better reveal the potential links between adjacent signals, which have been proved in the front part of this research. At the same time, PSS begins to increase in the region of 500-phase point, which can be considered as the threshold for the beginning of the RUL estimation. After that, the original data are divided into reference data and test data; reference data are viewed as the database; test data are considered as the current data to estimate the RUL. The results are shown in Figures14–16.Figure 14
RUL estimation for the outer-race-fault case based on PSS.Figure 15
RUL estimation for the inner-race-fault case based on PSS.Figure 16
RUL estimation for the ball-fault case based on PSS.It can be seen that the estimated RUL is close to the actual RUL. Furthermore, the comparisons with the original signal similarity and RMS are conducted to prove the superiority of the proposed method in RUL estimation. The results of estimation based on original signal similarity are shown in Figures17–19. Three criteria defined in the following equation are applied to evaluate the accuracy of the RUL estimation: (15)
F
1
=
1
N
∑
i
=
1
N
x
i
-
y
i
,
F
2
=
1
N
∑
i
=
1
N
x
i
-
y
i
2
,
F
3
=
1
N
∑
i
=
1
N
x
i
-
y
i
2
,where x
i is the estimated RUL, y
i is the actual RUL, and N is the length of the data points. The estimation results of the last 15 h are shown in Table 3.Table 3
Comparison of RUL estimation.
F
1
F
2
F
3
Outer race
RMS
6.0578
8.5576
5.1652
Original signal
1.5323
2.0537
1.0665
Phase space
1.3226
1.9092
0.8540
Inner race
RMS
10.8441
12.5821
8.2846
Original signal
4.1250
5.4658
3.0780
Phase space
3.6250
4.8606
2.7941
Ball
RMS
8.5964
10.2489
7.0659
Original signal
3.1875
4.1155
2.4702
Phase space
3
3.7583
2.3632Figure 17
RUL estimation for the outer-race-fault case based on original signal similarity.Figure 18
RUL estimation for the inner-race-fault case based on original signal similarity.Figure 19
RUL estimation for the ball-fault case based on original signal similarity.As for the findings in Table3, the proposed method has a higher accuracy, which indicates that the PSR can better grasp the nature and regularity of signal. It is effective in the RUL estimation of bearing, especially in the serious degradation stage.


## 5. Conclusions

In this paper, a novel method based on PSS is proposed to recognize the fault patterns, predict the fault size, and estimate the RUL of rolling bearing. The PSR is employed to transform the vibration signal into high-dimensional phase space. The NCC is employed to calculate the similarity between two health conditions, which is not affected by the scaling and shifting of the signals. The proposed method has two important characteristics. Firstly, the nonlinear vibration signal of bearing is represented by phase space, which can better grasp the nature and regularity of the signals. Secondly, the proposed feature of PSS is utilized in template matching and as a powerful feature for data-based bearing fault size prediction. To verify the effectiveness of the proposed method, the PSS is evaluated by the actual data acquired from bearing and compared with the original signal similarity. The experimental results show that the proposed method has a satisfying performance. According to the analysis, the PSS is also been suited to other machine faults diagnosis and RUL estimation, such as gears and cutting tools.Moreover, owning to the limited samples, the SVR structure cannot be fully trained to obtain a better prediction result. Part of the validation is not very comprehensive, which is the deficiency of this paper. Therefore, more experiments will be carried out in the future. At the same time, the technical fields in the proposed PSS can be further improved and strengthened.



---

*Source: 5341970-2016-04-26.xml*