A deep learning-based remaining useful life prediction approach for bearings

TL;DR

This paper introduces a deep learning framework using CNNs and Hilbert-Huang transform for accurate, online prediction of remaining useful life in bearings, improving reliability in industrial systems.

Contribution

It presents a novel CNN-based method combined with HHT for automatic RUL prediction, reducing feature extraction efforts and enhancing adaptability across conditions.

Findings

Outperforms existing RUL prediction methods in experiments.

Successfully transfers across different bearing operating conditions.

Demonstrates high accuracy and robustness in RUL estimation.

Abstract

In industrial applications, nearly half the failures of motors are caused by the degradation of rolling element bearings (REBs). Therefore, accurately estimating the remaining useful life (RUL) for REBs are of crucial importance to ensure the reliability and safety of mechanical systems. To tackle this challenge, model-based approaches are often limited by the complexity of mathematical modeling. Conventional data-driven approaches, on the other hand, require massive efforts to extract the degradation features and construct health index. In this paper, a novel online data-driven framework is proposed to exploit the adoption of deep convolutional neural networks (CNN) in predicting the RUL of bearings. More concretely, the raw vibrations of training bearings are first processed using the Hilbert-Huang transform (HHT) and a novel nonlinear degradation indicator is constructed as the label for learning. The CNN is then employed to identify the hidden pattern between the extracted degradation indicator and the vibration of training bearings, which makes it possible to estimate the degradation of the test bearings automatically. Finally, testing bearings' RULs are predicted by using a $\epsilon$-support vector regression model. The superior performance of the proposed RUL estimation framework, compared with the state-of-the-art approaches, is demonstrated through the experimental results. The generality of the proposed CNN model is also validated by transferring to bearings undergoing different operating conditions.

Tables

Table 1. Table I: Bearing frequencies of inner race, outer race and ball [ 28 ]

Symbol	Description	Expression
$f_{inner}$	Inner race frequency	$\frac{\eta }{2}\cdot f_{\omega }\cdot \left(1+\frac{{\ell }_{ball}}{{\ell }_{pitch}}\cdot \cos \Phi \right)$
$f_{outer}$	Outer race frequency	$\frac{\eta }{2}\cdot f_{\omega }\cdot \left(1-\frac{{\ell }_{ball}}{{\ell }_{pitch}}\cdot \cos \Phi \right)$
$f_{ball}$	Ball frequency	$\frac{{\ell }_{ball}}{{\ell }_{pitch}}\cdot f_{\omega }\cdot \left(1-\frac{{\ell }_{ball}^2}{{\ell }_{pitch}^2}\cdot {\cos }^2\Phi \right)$

Table 2. Table II: Physical characteristics and operation condition

Physical parameter	Value
Number of balls of the bearing ($\eta$)	13
Ball diameter of the bearing (${\ell }_{ball}$)	3.5 mm
Pitch diameter of the bearing (${\ell }_{pitch}$)	25.6 mm
Contact angle of the bearing ($\Phi$)	$0^{\circ }$
Rotation frequency ($f_{\omega }$), bearing1	1800 r/min
Rotation frequency ($f_{\omega }$), bearing2	1600 r/min
Maximum dynamic load (F), bearing1	4000 N
Maximum dynamic load (F), bearing2	4200 N

Table 3. Table III: Parameters in the CNN model

Layer	Filters	Kernel size/Stride	Output size
Input	…	…	1x2560x1
Conv1	64	1x100/50	1x50x64
Maxpooling1	…	1x2/2	1x25x64
Conv2	64	1x2/1	1x24x64
Maxpooling2	…	1x2/2	1x12x64
Flatten	…	…	1x768
FC1	…	…	1x100
Output	…	…	1x1

Table 4. Table IV: Proposed method and tested methods for comparison

	DEI	CNN	SVR
C1	-	✓	✓
C2	✓	-	✓
Proposed method	✓	✓	✓

Table 5. Table V: Comparison of predicted results between proposed approach and other tested methods

	Bearing1_4 ($T_{failure}$ =339 s)				Bearing1_5 ($T_{failure}$ =1610 s)				Bearing1_6 ($T_{failure}$ =1460 s)
Methods	${\hat{T}}_{failure}$	$E_r\%$	ETA		${\hat{T}}_{failure}$	$E_r\%$	ETA		${\hat{T}}_{failure}$	$E_r\%$	ETA
Proposed approach	340 s	-0.29%	0.96		1500 s	6.83%	0.79		1480 s	-1.37%	0.83
C1	30 s	91.15%	0.04		820 s	49.07%	0.18		1181 s	19.11%	0.52
C2	N/A (0 s)	N/A (100%)	0.03		1140 s	29.19%	0.36		1080 s	26.02%	0.41

Table 6. Table VI: Comparison of predicted results between proposed approach and existing major approaches.

	$E_r\%$					$S_{mean}$	MAE	NRMSE
Methods	Bearing1_4	Bearing1_5	Bearing1_6	Bearing2_4	Bearing 2_6
Proposed approach	-0.29%	6.83%	-1.37%	5.75%	1.55%	0.87	46.2	0.05
FCAMN [23]	21.95%	-15.22%	-5.74%	-	-	0.35	696.7	0.13
Multi-scale CNN [24]	10.69%	-148.20%	-21.51%	-	-	0.25	1003.3	0.50
CWT-CNN [34]	20.35%	11.18%	34.93%	-1.44%	-42.64%	0.46	265.8	0.29
RNN [35]	62.07%	-22.98%	21.23%	-19.42%	-13.95%	0.17	586.0	0.57
Particle Filtering [36]	5.60%	100.00%	28.08%	8.63%	58.91%	0.42	583.8	1.29
FFT+Ratio [37]	80.00%	9.00%	-5.00%	10%	49%	0.44	252	0.32
LSTM [38]	38.69%	-99.40%	-120.07%	19.81%	17.87%	0.26	996.2	0.57
SOM [39]	-20.94%	-278.26%	19.18%	51.80%	-20.93%	0.16	1164.2	1.02

Equations

$$a_{i,j}^0(t_p)=\{\begin{array}{cc}{S}_i(t_p),& j=1\\ S_i(t_p)-{\sum }_{j=1}^{n-1}{\text{IMF}}_{i,j}(t_p),& j=2,...,n\end{array}.$$

$$a_{i,j}^v(t_p)=a_{i,j}^{v-1}(t_p)-{\omega }_{i,j}^{v-1}(t_p),\forall p,$$

$${\text{IMF}}_{i,j}^A(t_p)={\text{IMF}}_{i,j}(t_p)+\text{i}{\text{IMF}}_{i,j}^H(t_p).$$

$${\text{IMF}}_{i,j}^H(t_p)=\frac{1}{\pi }{\int }_{-\infty }^{+\infty }\frac{{\text{IMF}}_{i,j}(s)}{t_p-s}ds.$$

$$\begin{array}[]{l}h_{i,j}(t_{p})=\sqrt{(\textup{IMF}_{i,j}(t_{p}))^{2}+(\textup{IMF}_{i,j}^{{}^{H}}(t_{p}))^{2}},\\ \phi_{i,j}(t_{p})=\arctan\left(\frac{\textup{IMF}^{H}_{i,j}(t_{p})}{\textup{IMF}_{i,j}(t_{p})}\right).\end{array}$$

$$f_{i,j}(t_p)=\frac{1}{2\pi }\frac{{\varphi }_{i,j}(t_p+\Delta t)-{\varphi }_{i,j}(t_p-\Delta t)}{2\Delta t}.$$

$$M_i(f_i,t_p)=\sum _{j=1}^n{h}_{i,j}(f_{i,j},t_p).$$

$$M_i(f_i)=\sum _{p=1}^P{M}_i(f_i,t_p).$$

$$L_i=\underset{f_i\in \{f_{inner},f_{outer},f_{ball}\}}{\max }{M}_i(f_i),i\in D.$$

$$L_{norm,i}=\frac{L_i-\min (\mathbf{L})}{\max (\mathbf{L})-\min (\mathbf{L})}\pm ϵ,ϵ\to 0$$

$${\mathbf{L}}_{\mathbf{norm}}=[L_{norm,1},...,L_{norm,N}].$$

$$u_m^k=W^k\ast u_e^{k-1}+B^k,$$

$$P_m^k=\underset{\gamma =1,\dots ,{\lambda }^k}{\max }{r}_{\gamma +(m-1)I_{pl}^k}^{(k-1)},$$

$$z=\frac{1}{N}\sum _{i=1}^N(L_{norm,i}-{\tilde{L}}_{norm,i}{)}^2,$$

$$\underset{w,b,{\xi }_g,{\hat{\xi }}_g}{\min }$$

$$\begin{array}{rl}f(\varphi (x_g))-L_{(g-1)s+l+1,test}& \le ϵ+{\xi }_g\\ L_{(g-1)s+l+1,test}-f(\varphi (x_g))& \le ϵ+{\hat{\xi }}_g\\ {\xi }_g\ge 0,{\hat{\xi }}_g\ge 0.& \end{array}$$

$$L(w,b,\alpha ,\hat{\alpha },\xi ,\hat{\xi },\mu ,\hat{\mu })=\frac{1}{2}{\Vert w\Vert }^2+C{\Sigma }_{g=1}^{\frac{Q-l-1}{s}+1}({\xi }_g+{\hat{\xi }}_g)$$

$$-{\Sigma }_{g=1}^{\frac{Q-l-1}{s}+1}{\mu }_g{\xi }_g-{\Sigma }_{g=1}^{\frac{Q-l-1}{s}+1}{\hat{\mu }}_g{\hat{\xi }}_g+{\Sigma }_{g=1}^{\frac{Q-l-1}{s}+1}{\alpha }_g(f(\varphi (x_g))$$

$$-L_{(g-1)s+l+1,test}-ϵ-{\xi }_g)+{\Sigma }_{g=1}^{\frac{Q-l-1}{s}+1}{\hat{\alpha }}_g(L_{(g-1)s+l+1,test}$$

$$-f(\varphi (x_g))-ϵ-{\hat{\xi }}_g).$$

$$w={\Sigma }_{g=1}^{\frac{Q-l-1}{s}+1}({\hat{\alpha }}_g-{\alpha }_g)\varphi (x_g),0={\Sigma }_{g=1}^{\frac{Q-l-1}{s}+1}({\hat{\alpha }}_g-{\alpha }_g),$$

$$C={\alpha }_g+{\mu }_g,C={\hat{\alpha }}_g+{\hat{\mu }}_g.$$

$$\begin{array}{r}\underset{\alpha ,\hat{\alpha }}{\max }\\ \end{array}$$

$$\begin{array}{rl}& {\Sigma }_{g=1}^{\frac{Q-l-1}{s}+1}({\hat{\alpha }}_g-{\alpha }_g)=0\\ & 0\le {\alpha }_g,{\hat{\alpha }}_g\le C.\end{array}$$

$$f(x)={\Sigma }_{g=1}^{\frac{Q-l-1}{s}+1}({\hat{\alpha }}_g-{\alpha }_g)\kappa (x,x_g)+b.$$

$$E_r=100\%\times \frac{T_{failure}-{\hat{T}}_{failure}}{T_{failure}}$$

$$\text{ETA}=\{\begin{array}{c}\exp (-\ln (0.5)\frac{E_r}{5})if{E}_r\le 0\\ \exp (+\ln (0.5)\frac{E_r}{20})if{E}_r>0\end{array}$$

$$S_{mean}=\frac{1}{\beta }\sum _{i=1}^{\beta }{\text{ETA}}_i$$

$$\text{MAE}=\frac{1}{\beta }\sum _{i=1}^{\beta }\mid T_{failur{e}_i}-{\hat{T}}_{failur{e}_i}\mid$$

$$\text{NRMSE}=\frac{\sqrt{\frac{1}{\beta }{\sum }_{i=1}^{\beta }(T_{failur{e}_i}-{\hat{T}}_{failur{e}_i}{)}^2}}{\frac{1}{\beta }{\sum }_{i=1}^{\beta }{\hat{T}}_{failur{e}_i}}$$

Full text

Nomenclature

CNN

Convolutional neural network.

DEI

Degradation energy indicator.

EMD

Empirical mode decomposition.

ETA

Exponential transformed accuracy.

HHT

Hilbert-huang transform.

IMF

Intrinsic mode function.

MAE

Mean average error.

NRMSE

Normalized root mean square error.

MSE

Mean square error.

REB

Rolling element bearing.

RUL

Remaining useful life.

SVR

Support vector regression.

FT or $L_{ft}$

Failure threshold.

$N$

Length of historical units of training bearing.

$Q$

Length of historical units of test bearing.

$U$

No. of predicted units of test bearing.

$\bf{S_{i}}$

Sensor measurement signal of $i$-th unit.

$\Delta t$

Sampling time.

$P$

Number of measurements in $i$-th unit.

$\tau$

Time interval between two recording phases.

$\bf{L}$

DEI of training bearing.

$\bf{L_{norm}}$

Normalized DEI of training bearing.

$\bf{L_{test}}$

Estimated DEI of test bearing.

$\bf{\widehat{L}_{U,test}}$

Predicted DEI of test bearing.

$K$

Total number of layers.

$W^{k}$

Weights in $k$-th convolutional layer.

$B^{k}$

Bias in $k$-th convolutional layer.

$\widehat{T}_{failure}$

Predicted RUL of test bearing.

$T_{failure}$

Real RUL of test bearing.

$\beta$

Number of test bearings.

$E_{r}\%$

Relative percentage error of prediction.

$S_{mean}$

Average score of prediction.

$\bf{X}$

Training set for SVR.

$\mathbb{R}$

The set of real numbers.

$\mathbb{Z}$

The set of positive intergers.

I Introduction

To constrain relative motion while reducing friction between moving parts, rolling element bearings (REBs) are one of the most widely used elements in industrial machinery. Prognostics and health management of bearings is of significance for safety, reliability and effectiveness of the mechanical systems [1, 2]. The literature [3] shows nearly half of motor failures are related to the degradation of bearings. As such, estimating the remaining useful life (RUL) (i.e., time-to-failure prognostics) of bearings has attracted a great deal of attention in recent years [4]. RUL prediction helps users monitor the condition of the bearings and provides an estimation of time left before a failure occurs. Compared with fault diagnosis, which has been well investigated over past few decades [5], the problem of RUL prediction studied in this paper is a relatively new and challenging topic due to the huge amount of uncertainties of environment and operating condition.

In general, RUL prediction approaches can be categorized into model-based and data-driven approaches. Model-based approaches aim to build a physical model to represent the degradation of the rolling bearing [6]. Li et al. [7] predicted the defect growth on a bearing unit using Paris’s law for fatigue. However, it is difficult to construct a precise physical degradation model due to the sensitivity of the model parameters and noised operating environments. This limits the practical applications of the model-based approaches. On the other hand, data-driven approaches benefit from the extensive expertise in signal processing and machine learning [8], and infer the degradation process of bearings without knowing any physics of degradation failure. The prognostic framework of the data-driven approaches mainly consists of three stages: 1) feature extraction from noisy sensory signals, which helps to build up the health indicator for the learning of system degradation behavior; 2) degradation models are trained on the training bearing using statistical or machine learning techniques; and 3) the degradation indicator of the test bearings can be estimated based on the model trained in the second stage. Then, the unknown degradation process can be predicted by applying regression techniques (i.e., $\epsilon$-SVR).

To extract features from raw signals, time-domain, frequency-domain, and time-frequency domain analysis are commonly adopted. Among them, the time-frequency analysis has been found to be the most efficient due to its ability to characterize transient signals over time and frequency domains [9]. Well-known time-frequency techniques for extracting bearing features include short-time Fourier transform [10], wavelets [11], Wigner-Ville distribution [12], and Hilbert-Huang transform (HHT) [13]. Implementation of the short-time Fourier transform is limited by its time-frequency resolution capability; for instance, low frequencies are difficult to identify with short windows. On the other hand, wavelets and the Winger-Ville distribution provide richer pictures than short-time Fourier transform; however, their effectiveness depend on estimating the Hurst parameter and the quality of the analyzed signal. HHT shows better computational efficiency and resolution over other time-frequency analysis. which uses the techniques of empirical mode decomposition (EMD) and Hilbert transform (HT) to decompose the original vibration signal into a number of intrinsic mode functions (IMFs) in various frequency scales. Frequency components of each IMF are related to both the sampling frequency and the signal itself, thus demonstrating that HHT is a self-adaptive signal processing technique perfectly suited to non-stationary signals. Wu et al. [14] analyzed the time-to-failure prognostics of REBs, which extracts ten statistical features using time and frequency analysis and eleven IMF features using HHT time-frequency analysis. The gear fault identification method proposed in [15] is based on the HHT, and first six IMFs are selected as inputs for SOM neural networks for fault diagnosis.

Data-driven RUL prediction approaches are mainly based on statistical and machine learning techniques, such as artificial neural networks (ANN) [16], fuzzy logic systems [17], and auto-regressive (AR) models [18]. The computational cost of ANN is relatively high in terms of optimizing the weights of the model. The performances of the AR models and fuzzy logic systems require precise trend of historical observations and high-quality training data, respectively. Recently, deep learning has merged into research and industry fields, and has beaten other machine learning techniques in speech recognition and image recognition tasks [19]. Deep learning model is good at discovering high level abstractions from labeled data using a back-propagation algorithm [20]. Specifically, it learns feature representations automatically rather than designing the hand-created features by experience. As the most well-known model in deep learning, in recent years, CNN dominates the recognition and detection problems in computer vision domain, which is distinguished by three characteristics, namely local connections, shared weights, and local pooling [21, 22]. The first two characteristics show that the CNN model requires fewer parameters to detect local information than multilayer perceptron, while the last characteristic ensures shift invariance to the networks. Typically, 1-D CNN will be employed to this work to learn the latent space of input sensory time-series vibrations, which has been applied with great success to speech recognition and document reading tasks. Few attempts have been made for the prediction problem using CNN-based models [23, 24, 25, 26]. This paper exploits the adoption of CNN technique in estimating the RUL of bearings, as a prognosis problem, to learn about the nonlinear degradation behavior according to raw vibration data and an extracted label. Instead of using the CNN technique to perform the time-series prediction, the main function of the CNN model in this paper is to reveal the hidden dependencies between the vibration data and the DEI of the training bearing, which makes full use of the advantages of CNN in automatic feature extraction.

In this work, we propose a data-driven framework for predicting the RUL of REBs by applying the HHT, CNN, and $\epsilon$-Support vector regressorion ($\epsilon$-SVR). The raw vibration signals collected from sensors are processed by the HHT method and a novel time-series degradation indicator, i.e., DEI, is constructed. Subsequently, a CNN model is trained to learn the features from the input raw vibration to the DEI label on the training bearings, and used to predict the DEIs of testing bearings. Then, a $\epsilon$-SVR model is introduced so that the evolution of the degradation can be forecast till the bearing failure. The effectiveness of the proposed based framework for RUL prediction is validated on an experimental platform (i.e., PRONOSTIA). Much lower RUL prediction errors are achieved, compared with eight existing approaches in previous papers and two tested methods designed in this paper, indicating the superior performance of the proposed method.

This work makes the following contributions: 1) The proposed method successfully extracts a novel nonlinear degradation energy indicator (DEI) (see Fig. 1, compared with the linear time degradation indicator) to describe the degradation trend of the training bearing, according to the nature frequencies of bearing components; 2) The proposed CNN architecture is general and robust for similar operation conditions, it can transfer to another bearing undergoing different operating condition and obtain good prediction results, without changing CNN hyper parameters and the depth of layers; 3) The propose DEI is an integrated indicator with regards to the maximum vibration levels among different bearing components, which considers all the possible detects on the rolling element bearing. This is a more realistic indicator as the localized defects are not initially initiated in real industrial applications, meaning that all the types of defects have to be considered; and 4) CNN scales all the indicators of training bearings and test bearings into a consistent latent space. Thus, training and testing can share a same failure threshold (FT), i.e., the maximum value of the indicator for the training bearing.

The outline of this paper is as follow. Section II presents the proposed RUL prediction framework with technical details. In Section III, experimental results obtained from bearing degradation tests are carried out. Thereby, the performance of the proposed framework is validated and the results show improved accuracy in predicting the RUL compared with eight state-of-the-art approaches and two designed test methods in this paper. Section IV summarizes the paper and discusses future works.

II Degradation indicator training and RUL prediction algorithm

The overall framework for the prediction of the RUL can be decomposed into three parts. The schematic of the overall framework is shown in Fig. 2. The key challenges of this work involve: obtain the DEI to represent the degradation behaviour; establish a CNN model to map raw vibration signal to the DEI; and construct an $\epsilon$-SVR to predict the RUL. Thus, in the following subsections, the explicit expression of degradation feature extraction, CNN model, and $\epsilon$-SVR forecasting model will be derived in Section II-A, Section II-B, and Section II-C.

II-A Degradation indicator extraction

To begin with, for a training bearing, it is assumed that the raw vibration signal till the end of lifetime with $N\in\mathbb{Z}$ historical units have been acquired. The sensory vibration signal $\mathbf{S_{i}}=\{S_{i}(t_{p})\}^{P}_{p=1}$, with sampling time $\Delta t$, is measured at each historical unit for $i\in D=\left\{1,2,...,N\right\}$, where $P$ is the number of measurements recorded in each historical unit.

EMD is a self-adaptive method which is normally applied to analyze non-stationary and nonlinear signals. It decomposes the raw vibration data $\bf{S_{i}}$ into $n$ number of IMFs, illustrating the natural oscillation modes from fast to low oscillations.

For the $i$-th unit, the $j$-th mode ($j=1,\ldots,n$) of IMF, IMF${}_{i,j}(t_{p})=a^{v}_{i,j}(t_{p})$ for $p=1,\ldots,P$, is calculated iteratively associated with the iteration number $v=0,1,\ldots$.

First let $v=0$ and initialize $a^{0}_{i,j}(t_{p})$ for $p=1,\ldots,P$ by:

[TABLE]

Define IMF${}_{i,j}(t_{p})=a^{0}_{i,j}(t_{p})$ for $p=1,\ldots,P$ only if the IMF meets the following two conditions:

(I)

The IMF should have one or zero difference between the extrema number and the number of zero crossing,

(II)

along the time axis, the average value of upper and lower bound of the IMF should be zero everywhere.

Otherwise, update the IMF${}_{i,j}(t_{p})=a^{v}_{i,j}(t_{p})$ for $p=1,\ldots,P$ through the following iteration procedure with $v=1,2,\ldots$ until the IMF satisfies both aforementioned conditions (I) and (II):

[TABLE]

where $\omega^{v-1}_{i,j}(t_{1}),\ldots,\omega^{v-1}_{i,j}(t_{P})$ are the mean values of the upper envelope and the lower envelope of $a^{v-1}_{i,j}(t_{1}),\ldots,a^{v-1}_{i,j}(t_{P})$.

Once obtaining all IMFs, the analytical form of each IMF can be written as:

[TABLE]

where i is the imaginary part of the $\textup{IMF}^{A}_{i,j}(t_{p})$. $\textup{IMF}^{H}_{i,j}(t_{p})$ is the Hilbert transformation of $\textup{IMF}_{i,j}(t_{p})$ by convolution with function $\frac{1}{\pi t}$, given as:

[TABLE]

By this means, we can calculate the instantaneous amplitude $h_{i,j}(t_{p})$ and phase $\phi_{i,j}(t_{p})$

[TABLE]

Accordingly, it is easy to derive the instantaneous frequency

[TABLE]

for $p=1,\ldots,P$.

Then, the Hilbert Spectrum of $\bf{S_{i}}$ is obtained by

[TABLE]

The marginal Hilbert spectrum (MHS) $M_{i}(f_{i})$ can be written as:

[TABLE]

Nature frequencies of bearing components depend on the geometry of the bearing and its rotation speed. Expression of these frequencies are given in Table I, in where $\eta$ is the number of balls, $f_{\omega}$ is the rotation frequency, $\Phi$ is the contact angle, and $\ell_{ball}$ and $\ell_{pitch}$ are the ball diameter and pitch diameter, respectively.

With the bearing frequencies of different components, the value of the DEI, $L_{i}$, at historical unit $i$ is defined as the maximum value of the MHS by substituting the $f_{inner}$, $f_{outer}$ and $f_{ball}$ into $M_{i}(f_{i})$, given that

[TABLE]

The extracted DEI \bf{L}$$=[L_{1},...,L_{N}] is nromalized before training:

[TABLE]

for $i=1,\ldots,N$, and $\epsilon$ is an infinitesimal that used to avoid the value of label equal to 0 or 1. Thus, the normalized DEI is

[TABLE]

II-B DEI pattern learning

In this work, layers with repeated components are stacked in a CNN architecture, including convolutional layers, pooling layers, fully connected layers, and a regression layer [21].

Convolutional layer contains organized patches in convolutional layers, each patch is calculated by composing the features of the previous layer through a filter bank with the following equation:

[TABLE]

where $u^{k}_{m}\in\mathbb{R}$ denotes the output of the $m$-th unit in the $k$-th layer. $u^{(k-1)}_{e}\in\mathbb{R}^{1\times o^{k}}$ is the input data the $e$-th sub-vector in the previous layer $k-1$, where $o^{k}$ is kernel size in layer $k$. $W^{k}\in\mathbb{R}^{1\times o^{k}}$ and $B^{k}\in\mathbb{R}$ denote the connecting weights and bias in the $k$-th layer, respectively. ‘$\ast$’ means the convolution operation. It is noted that when $k=1$, $u^{(k-1)}_{e}$ is a sub-vector of the raw vibration data $\bf{S_{i}}$. Here we define all neurons in each layer is ${\bf{u}}^{k}=[u^{k}_{1},\ldots,u^{k}_{m},\ldots,u^{k}_{G^{k}}]$ for $k\in D_{k}=$ {1,2,$\ldots$,$K$}, where $G^{k}\in\mathbb{Z}$ is the number of neurons in the $k$-th layer and $K$ is the total number of layers. For convolutional layer, $G^{k}=(G^{k-1}-o^{k})/I^{k}_{cv}+1$ and $I^{k}_{cv}\in\mathbb{Z}$ is the stride in convolutional layer.

Activation function is introduced after convolutional layer. Among various activation functions, Rectified linear unit (ReLU) ${r}^{k}_{m}=\max(0,{u}^{k}_{m})$ is chosen as the nonlinear activation function to prevent the issue of vanishing gradient which may significantly increase the training time or even lead to the non-convergence.

Pooling layer is then used as a nonlinear down-sampling layer to extract the maximum feature values in each patches of the input data. Its function is to save computation time and downsize the number of parameters of the model as well as control overfitting. More specifically, pooling transforms small windows into single values by maxing or averaging. Consequently, the features extracted within the small window are similar and therefore illustrating the shift invariance property of CNN. Max-pooling layer is selected in this work as it is an algorithmic choice to ensure the generalization of neural networks [29], which is given by:

[TABLE]

where $\lambda^{k}\in\mathbb{Z}$ is the pooling size, and $I^{k}_{pl}$ is the stride in max pooling layer.

Fully connected layer and regression layer, like a classic ANN network, take the results of the convolution and max-pooling processes and use them to generate a predicted label. Since we use a normalized DEI $\mathbf{L_{norm}}\in\mathbb{R}^{1\times N}$ as the label for learning, the sigmoid function $\text{sigm}({\bf{u}}^{K-1})$ with an output value between (0,1) is applied to the last layer for normalized output. Hereby, mean square error (MSE) function is used to compute the loss with the expression:

[TABLE]

in where the proposed CNN model is minimizing the loss function $z$ between ground label DEI $\mathbf{L_{norm}}$ and predicted label $\tilde{\bf{L}}_{\bf{norm}}$. Algorithm 1 outlines the proposed CNN modeling procudure.

II-C RUL prediction

With the obtained CNN model, for a new test bearing with $Q\in\mathbb{Z}$ historical units, the estimated DEI ${\bf{L_{test}}}=[L_{1,test},\ldots,L_{Q,test}]\in\mathbb{R}^{1\times Q}$ can be automatically generated by the trained CNN with the new vibration signal $\bf{S_{i}}$, where $i\in{D_{test}}=\left\{1,\ldots,Q\right\}$. Then, to predict the RUL $\widehat{T}_{failure}\in\mathbb{R}^{+}$, a $\epsilon$-SVR forecasting model [31] is formalized to predict the upcoming degradation $\widehat{L}_{Q+1,test}$, $\widehat{L}_{Q+2,test}$, $\ldots$ based on the estimated DEI $\bf{L_{test}}$ by sliding window method. The forecasting process contains three sub-steps:

Extract training features from the estimated DEI $\bf{L_{test}}$ over a sliding window. The schematic of this step is illustrated in Fig. 4. The estimated DEI $\bf{L_{test}}$ is decomposed into overlapping windows associated with sampling window size $l$ and sliding size $s$. $x_{g}=(\mu_{g},\sigma_{g}^{2})$ for $g\in\{1,\ldots,\frac{Q-l-1}{s}+1\}$ represents a training feature for $\epsilon$-SVR, where $\mu$ denotes the mean value and $\sigma^{2}$ denotes the variance of each sampling window. Thus, the training set for $\epsilon$-SVR ${\bf{X}}=[(x_{1},L_{l+1,test}),\ldots,(x_{g},L_{(g-1)s+l+1,test}),\ldots,(x_{\frac{Q-l-1}{s}+1},$ $L_{Q,test})]$ is obtained, where $L_{(g-1)s+l+1,test}$ corresponds to the next value in $\bf{L_{test}}$ of the $g$-th sampling window; 2. 2.

$\epsilon$-SVR modeling is described in Algorithm 2, while at the application level, two parameters (distance limit $\epsilon\in\mathbb{R}$ and penalty parameter $C\in\mathbb{R}$) can be set manually when training the prediction model. A radial basis function (RBF) is necessary when we intend to train a nonlinear model; 3. 3.

The SVR model $f(x)$ learned in Algorithm 2 is then used to predict the RUL by sliding window method (with the same $l$ and $s$ in step 1). Since the test bearing undergoing same operating condition as the training bearing, it is reasonable to define the FTs (denoted as $L_{ft}$) of the test bearing equals to the last feature of the DEI of the training bearing, such that $L_{ft}$=$L_{N}$. Hence, the first prediction can be calculated as $\widehat{L}_{Q+1,test}=f(x_{\frac{Q-l}{s}+1})$, and the predicted DEI ${\bf{\widehat{L}_{U,test}}}=[\widehat{L}_{Q+1,test},\ldots,\widehat{L}_{Q+U,test}]$ can be obtained by shifting the sampling window, with $\widehat{L}_{Q+U-1,test}\leq L_{ft}\leq\widehat{L}_{Q+U,test}$. This will consequently lead to $\widehat{T}_{failure}=U\times\tau$, where $\tau$ is the time interval between two recording phases.

III Experiments

III-A Data description

The validation of the proposed RUL prediction framework is conducted on an experimentation platform named PRONOSTIA (see Fig. 5). This platform is built as a combination of three parts: a rotating part, a loading part, and a measurement part. The rotation of the test bearing is driven by the low speed shaft whose rotating torque is transmitted from an AC motor. A radial force generated by this loading part is applied on the external ring of the testing bearing. Since this external radial force exceeds the bearing’s allowable dynamic load (4000N), the degradation behavior is accelerated so that we can observe its degradation process within a relatively short time in few hours. During experiment tests, two high-frequency accelerometers (Type DYTRAN 3035B) are placed orthogonally on the external race of the test bearing to acquire the horizontal and vertical vibrations respectively. The accelerometer bandwidth is 0.5$\sim$10 kHz ($\pm$ 5%), and its natural/resonant frequency is 45 kHz. In this work, we extract our degradation labels by using the horizontal vibrations.

Bearing1 is chosen to validate the proposed framework. More specifically, the training set bearing1_2 is used for extracting the DEI and training the CNN model. Test sets bearing1_4, bearing1_5, and bearing1_6 are then used for estimating their DEIs and predicting the RULs by applying $\epsilon$-SVR forecasting model. Results of the Bearing2 under different rotational frequency and external dynamic load are also provided and compared. The geometry parameters and the operation conditions of the bearing1 and bearing2 are listed in Table II. Sampling frequency of the vibration sensor is 25.6 kHz. 0.1 s accelerometer vibration signals are recorded at a fixed time interval $\tau=10$ s. Therefore, each recording phase contains $p=2560$ measurements. More detailed description of the data set, bearings, and sensors can refer to the data description in [32].

III-B Degradation indicator extraction

The DEI of the bearing1_2 is extracted by substituting its outer ring frequency $f_{outer}$= 168 Hz, inner ring frequency $f_{inner}$ = 221 Hz, and ball frequency $f_{ball}$ = 215.4 Hz into Eq. (9). The evolution of the final extracted DEI $\bf{L}$ of the bearing1_2 is showed in Fig. 6(a). We also show the time evolutions of intermediate features $M_{i}(f_{inner})$, $M_{i}(f_{outer})$, and $M_{i}(f_{ball})$ of Eq. (9) in Fig. 7. It can be observed that the magnitude value of each time point in Fig. 6(a) is the maximum of $M_{i}(f_{inner})$, $M_{i}(f_{outer})$, and $M_{i}(f_{ball})$ at that time.

III-C Degradation indicator estimation

We use the vibration signal in the horizontal direction of the bearing1_2 as the input of the CNN, and the normalized DEI is used as the label which contains historical units $N=871$. The normalized DEI, $\bf{L_{norm}}$, is shown in Fig. 6(b). A less complex architecture of the CNN model is designed to improve the robustness of the network. As shown in Fig. 3, our finalized CNN model consists of $K=6$ layers: two convolutional layers (Conv1 and Conv2), two max-pooling layers (Maxpooling1 and Maxpooling2), one fully connected layer (FC1), and one regression layer for output. Before model training, Adam is set to be the optimizer with a small value 0.00001 as it guarantees a quick loss convergence compared with a larger or smaller learning rate for the CNN training process. The activation function in the output layer is the Sigmoid function, while ReLU function is used in the previous layers. Details of parameters in the proposed CNN model are concluded in Table III. The convolutional window sizes (kernel sizes) of convolutional layers are set to a large value 100 and a small value 2, respectively. The kernel size of Conv1 is relevantly large in order to extract more features from the raw vibration signal for more impressive power, meanwhile, small kernel size is selected for Conv2 to prevent overfitting. Hyper-parameters are obtained after 1000 iterations of training.

The estimated DEIs $\bf{L_{test}}$ as the output of the CNN model are shown in Fig. 8. In Fig. 8(a), estimated DEI of the training bearing1_2 shows the similar time evolution as the DEI label in Fig. 6(b). The final estimated DEI value of the bearing1_2, $L_{ft}=0.9756$, is defined as the failure threshold for the test bearings in Fig. 8(b)-(d).

III-D RUL prediction

As presented in Section III-C, the estimated DEIs have obtained from the trained CNN model. However, the DEIs of the test sets shown in Fig. 8(b)-(d) do not reach their fault limit, which need a regression algorithm to predict the estimated DEI till the end-of-life of each test bearing. A $\epsilon$-SVR method is proposed to predict the upcoming degradation process of the test bearings. We conduct a training on the predicted DEI of the bearing1_2. The sampling window size $l$ and the moving size $s$ in Fig. 4 is set to 50 and 1, respectively. The kernel function used in the prediction case is an RBF and penalty parameter C of the error term is chosen as 5.09 by grid search and cross validation using the method of GridSearchCV [33]. We estimate the DEI after 1000 steps based on the existing DEI, and using the maximum value of predicted DEI of bearing1_2 (i.e., $L_{tf}=0.9756$) to limit the termination time of the test sets.

Fig. 9(b)-(d) show the predicted RULs $\widehat{T}_{failure}$ of the test bearings till the failures occur. Red lines represent the predicted evolution of the bearings’ degradation behavior using the $\epsilon$-SVR method. The RULs are calculated as the difference between the final time when DEI reaching the failure threshold and the time of the last known point of the test bearings. For bearing1_4, the predicted RUL is 340s, while 1500s and 1480s are the predicted RULs for bearing1_5 and bearing1_6, respectively.

III-E Comparison and discussion

To assess the accuracy of the proposed method and compare to other existing approaches, two metrics are commonly adopted: 1) The relative percentage error ($Er\%$) which is given by Eq. (15); and 2) The exponential transformed accuracy (ETA) proposed in IEEE PHM 2012 [32]. ETA is an assessment index to distinguish the seriousness of the underestimate and overestimate of RUL prediction. It is clearly that underestimate (i.e., early warning) is preferred than overestimation (i.e., warning after damage) to prevent more severe damage of the bearing. The formulas are expressed in Eq. (16).

[TABLE]

[TABLE]

where where $T_{failure}$ is the real RUL for the test bearing. A higher \mid$$Er\%$$\mid means a worse RUL prediction result. On the other hand, ETA value varies from 0 to 1, and a higher score means a better RUL prediction result. In this work, $Er$ is the common choice of most previous literature, thus it will be used for comparison.

In addition to the two metrics that evaluate the prediction performance for a specific bearing, three more assessment metrics, namely average score $S_{mean}$, mean average error MAE, and the normalized root mean square error NRMSE, to make a comprehensive comparison of different methods, which are with the expression forms

[TABLE]

[TABLE]

[TABLE]

where $\beta$ is the number of test bearings.

To verify the benefits of the DEI and CNN techniques on RUL prediction, here we also develop other two tested methods for comparison purpose (see Table IV). The proposed method and the tested methods are denoted and explained as follows:

C1: CNN and $\epsilon$-SVR: This tested method uses a conventional linear time degradation label for CNN training rather than the nonlinear DEI. By this means, we can illustrate the impact of the DEI on RUL prediction. 2. 2.

C2: DEI and $\epsilon$-SVR: Without training the CNN model for feature extraction, DEI in this tested method is extracted manually and the $\epsilon$-SVR is followed for the RUL prediction. Note that computing a DEI for a new bearing requires high computational power and longer time. In the meantime, the FT of each test bearing has to be pre-defined artificially, which increase uncertainties of the RUL prediction affected by different working conditions. By this means, we can illustrate the impact of CNN modeling on the prediction of the final RUL. 3. 3.

Proposed method: DEI-based CNN and $\epsilon$-SVR: This is the proposed framework which integrates DEI extraction, CNN, and $\epsilon$-SVR into one framework.

The predicted numerical errors of the test bearings with the proposed approach and the tested methods are listed in Table V. Our approach achieves $Er\%$ of -0.29%, 6.83%, and -1.37% for bearing1_4, bearing1_5, and bearing1_6, respectively, which are much more smaller than the C1 and C2 methods. C1 uses a linear time degradation label for the training of the CNN model. The results show more than 19% prediction errors for test bearings and even 91.15% prediction error is obtained from bearing1_4, indicating that time degradation label is less effective than the DEI for the CNN training process. C2 is the method extracting the degradation indicator of the test bearings and define the FTs manually. As testing bearing1_4, 1_5, and 1_6 operate under same working conditions of bearing1_2, we employ the maximum and minimum value of bearing1_2 to normalize the extracted DEIs of C2 method. With same working condition and same normalization parameters, FT of the proposed method could be reasonablely used in C2 as well. To evaluate the impact of CNN modeling on the estimation of the final RUL, for the test bearings, the estimated DEIs extracted using HHT and calculated by trained CNN model are compared in Fig. 10. Without the CNN modeling procedure, one of the main drawback of the C2 method is that it requires long time to calculate ($\sim$2 s of each sampling period). This limits the practical application of this method in industry. Moreover, Fig. 10(a) shows that due to the uncertainties and huge mount of noise, the DEI extracted by C2 method has already exceeded the FT before the exact failure time, resulting in a 100% $E_{r}$. Similarly, in Fig. 10(c), the DEI extracted of C2 method is much noised than it of proposed method. At 16470 s, the magnitude of DEI:C2 is almost close to the FT, this phenomenon might lead to a waste of sources due to much underestimated of RUL. The comparison results in Table V demonstrate the benefits of using DEI and CNN in estimating RULs of REBs.

To demonstrate the generality of our proposed RUL estimation framework, other test bearings (i.e, bearing2_4 and bearing2_6) that operates undergoing a different external load and rotational speed are analyzed. Bearing2_2 is used for training, and the trained CNN model is obtained without changing any hyper-parameters and the architecture of model in Table III. We just fine-tune the $\epsilon$-SVR forecasting model by changing penalty parameter C from 5.09 to 7.09, resulting in 5.75% and 1.55% $E_{r}$ for Bearing2_4 and bearing2_6, respectively. Good RUL prediction results of bearings with different operating conditions indicate the repeatability and robustness of our proposed method, with respect to the hyper parameters and the architecture of the CNN.

To further verify the proposed approach, the predicted numerical errors of RULs generated by the proposed method and eight published methods are compared and listed in Table VI. Other published approaches include a recurrent neural network method based health indicator [35], the method proposed by the winner of the IEEE PHM 2012 prognostic [37], a convolutional long-short-term memory network (LSTM) method [38], and a self-organization Mapping (SOM) method, etc.. In addition, recent CNN-based approaches, including the frozen convolution and activated memory network (FCAMN) [23], multi-scale CNN [24], and continuous wavelet transform CNN (CWT-CNN) [34] are also compared. The results of the comparison shown in Table VI confirm that our approach significantly outperforms the referenced methods with the highest average score $S_{mean}=0.87$, and the smallest MAE = 46.2 and NRMSE = 0.05. In particular, a -0.29% $E_{r}$ of bearing1_4 is achieved, owing to a 1s absolute time error. This result benefits from a good nonlinear degradation indicator extracted using the HHT method. In addition, the CNN is a powerful tool for discovering the hidden pattern of the extracted degradation indicator and the underlying bearing system, further increasing the accuracy of the predicted RUL.

It can be concluded from the experiment results that the proposed data-driven RUL estimation approach has much smaller prediction errors, compared with both the tested methods in this work and the other published methods in previous studies.

IV Conclusion

In this paper, a data-driven framework for RUL prediction of rolling element bearing is presented using the HHT method, a CNN model, and an $\epsilon$-SVR forecasting model. A nonlinear degradation indicator DEI is first extracted from the raw vibration signals using the HHT method, which is defined as the label for the training. A CNN model is trained to discover the hidden pattern between the extracted DEI and the raw vibration data of the training bearing. In this way, predicted DEIs are automatically obtained when applying the trained CNN model to the test bearings. Finally, the RULs of the testing bearings are obtained using an $\epsilon$-SVR forecasting model. An experimentation platform that allows to observe the accelerated degradation process of bearings is employed to validate the proposed framework. The proposed framework achieves much smaller prediction errors for RUL predictions than previous published approaches.

Future work includes the application of the proposed framework to a wider range of case studies on experimental data in other applications [40], and the investigation of other potential degradation labels to achieve even higher accuracy in estimating RUL.