# Bearing Degradation Process Prediction Based on the Support Vector Machine and Markov Model

Authors: Shaojiang Dong; Shirong Yin; Baoping Tang; Lili Chen; Tianhong Luo
Journal: Shock and Vibration
 (2014)

---

## Abstract

Predicting the degradation process of bearings before they reach the failure threshold is extremely important in industry. This paper proposed a novel method based on the support vector machine (SVM) and the Markov model to achieve this goal. Firstly, the features are extracted by time and time-frequency domain methods. However, the extracted original features are still with high dimensional and include superfluous information, and the nonlinear multifeatures fusion technique LTSA is used to merge the features and reduces the dimension. Then, based on the extracted features, the SVM model is used to predict the bearings degradation process, and the CAO method is used to determine the embedding dimension of the SVM model. After the bearing degradation process is predicted by SVM model, the Markov model is used to improve the prediction accuracy. The proposed method was validated by two bearing run-to-failure experiments, and the results proved the effectiveness of the methodology.


---

## Body

## 1. Introduction

Bearing is one of the most important components in rotating machinery. Accurate bearing degradation process prediction is the key to effective implement of condition based maintenance and can prevent unexpected failures and minimize overall maintenance costs [1, 2].To achieve effective degradation process prediction of the bearing, firstly, the features should be extracted from the collected vibration data. Then, based on the extracted features effectively prediction models should be selected [3]. Feature extraction is the process of transforming the raw vibration data collected from running equipment to relevant information of health condition. There are three types of methods to deal with the raw vibration data: time domain analysis, frequency domain analysis, and time-frequency domain analysis. The three types of methods are often chosen to extract the feature. For example, Yu [4] chose the time domain and the frequency domain transform to describe the characteristics of the vibration signals. Yan et al. [5] chose the short-time Fourier transform to extract the features. Ocak et al. [6] chose the wavelet packet transform to extract the feature of bearing wear information. Because the frequency features from FFT analysis results often tend to average out transient vibrations and thus not providing a wholesome measure of the bearing health status, in this paper, the time domain and the time-frequency domain characteristics are used to extract the original features.Although the original features can be extracted, they are still with high dimension and include superfluous information. So the original features fusion and dimensional reduction method should be used to deal with the original features so as to select the typical features. The most commonly used features fusion and dimensional reduction method is principal component analysis (PCA) [7, 8]. But the PCA is mainly used for dealing with the linear data set, while the bearing vibration features are usually suppressed by the nonlinear characteristic features, so the PCA cannot work effectively. Therefore, it is a challenging to find an effective nonlinear features fusion and dimensional reduction method. In this research a new feature extraction method local tangent space alignment (LTSA) [9] is chosen. The LTSA is an efficient manifold-learning algorithm, which can be used as a preprocessing method to transform the high dimensional data into more easily handled low dimensional data [10]; the method has been used in many fields, such as face recognition, character recognition, and image recognition [11, 12]. In this paper, the LTSA is used to achieve extracting the more sensitive features.After selecting the typical features, another challenge is how to effectively predict the bearing degradation process based on the extracted features. The existing equipment degradation process prediction methods can be roughly classified into model-based (or physics-based modals) and data-driven methods [13]. The model-based methods predict the equipment degradation process using the physical models of the components and damage propagation models based on damage mechanics [14, 15]. However, equipment dynamic response and damage propagation processes are typically very complex, and authentic physics-based models are very difficult to build [16]. Data-driven methods, also known as artificial intelligent approaches, are derived directly from routine condition monitoring data of the monitored system, which predicts the failure progression based on the learning or training process. The more prior the data is used for the training process, the more accurate is the model obtained [17]. Artificial intelligent techniques have been increasingly applied to bearing remaining life prediction recently. Lee et al. [18] presented an Elman neural network method for health condition prediction. Huang et al. [19] proposed a back-propagation network-based method for bearing degradation process prediction. However, the neural networks have the drawbacks of slow convergence; difficulty in escaping from local minima; uncertain network structure, especially when doing the bearing degradation process prediction problem with large data, and those problems will be more troublesome. The SVM [20] is most widely used recently and has a good identify and regression ability. In this paper, the SVM is used to predict the bearings degradation process.Although the SVM is effective in predicting the bearing running state, the prediction error still exists. Because any prediction methods based on the historical data for future prediction will more or less have some prediction error, it is necessary to improve the prediction results. However, the prediction error has the character of being affected by many factors, fluctuations, showing a great random, and the error points are not related. So if we want to achieve high prediction accuracy, we need to find the discipline of the prediction error and correct the error of prediction. The Markov model [21] used the state transition matrix to achieve a more precise prediction. It can be used to achieve the improvement of the prediction accuracy. So in this research, the Markov model is used to improve the prediction accuracy.The remainder of this paper is organized as follows. The methods of features extraction and the theory of dimensional reduction method LTSA are introduced in Section2. The SVM model and the Markov model for bearing degradation process prediction are described in Section 3. In Section 4, the flowchart and the procedure of this research are introduced. The case validation and actual application are presented in Section 5. Finally, the conclusions are given in Section 6.


## 2. Methods of Signal Processing and Dimensional Reduction

### 2.1. Feature Extraction

This section presents a brief discussion on original feature extraction from time domain, time-frequency domain of vibration signals. Time domain methods usually involve statistical features that are sensitive to impulsive oscillation, such as kurtosis, skewness, peak-peak (P-P), RMS, and sample variance. The 5 domain statistical features are used as original features in time domain has been used in the literatures [2, 3]:
(1)RMS=∑xi2N,
where N is the number of discrete points and xi represents the signal value at those points,
(2)Variance=s2=1(N-1)∑i=1N(xi-x-)2,
where x- is the mean value.Theith central moment for a set of data is defined as
(3)Mi=1N∑i=1N(xi-x-)i.The normalized forth moment, kurtosis, which is commonly used in bearing diagnostics, is defined as(4)Kurtosis=M4(M2)2-3.The skewness is defined as(5)Skewness=M3(M2)3/2.The peak-peak is defined as(6)peak-peak=max⁡(xi)-min⁡(xi).Empirical mode decomposition (EMD) is a powerful tool in time-frequency domain analysis. The advantage of EMD is the presentation of signals in time-frequency distribution diagrams with multiresolution, during which choosing some parameters is not needed. This property is essential in the detection of bearing faults. The EMD energy can represent the characteristic of vibration signals, and thus it is used as the input features. The (intrinsic mode function) IMF energy data sets are chosen as original features in this paper. The original features for bearing degradation process prediction based on the original features are shown in Table1.Table 1
Original features for bearing fault diagnosis and performance assessment.
Domain
Feature
Time domain
RMS, Kurtosis, Skewness, Peak-Peak, and sample variance
Time-frequency domain
EMD energy (Fe1,Fe2,Fe3,…,Fen)


### 2.2. Dimensional Reduction Based on the LTSA

Because the generated original feature sets are still with high dimension and include superfluous information, the feature extraction method LTSA is used to fuse the relevant useful features and extracts more sensitive features to work as the input of the proposed prediction model.The basic idea of LTSA is to use the tangent space of sample points to represent the geometry of the local character. Then these local manifold structures of space are lined up to construct the global coordinates. Given a data setX=[x1,x2,…,xN], xi∈Rm, a mainstream shape of d-dimension (m>d) is extracted. The LTSA feature extraction algorithm is as follows [9].(1)
Extract local information: for eachxi, i=1,2,…,N, use the Euclidean distance to determine a set xi=[xi,1,xi,2,…,xi,ki] of its neighborhood adjacent points (e.g., k the nearest neighbors).(2)
Local linear fitting: in the neighborhood of data pointsxi, a set of orthogonal basis Qi can be selected to construct the d-dimension neighborhood space of xi and the orthogonal projection of each point xi,j(j=1,2,…,N) can be calculated to the tangent space of θj(i)=QiT(xi,j-xi¯). xi¯ is the mean data for the neighborhood. The orthogonal projection in the tangent space of neighborhood data of xi is composed of local coordinate Θi=[θ(i),1,θ(i),2,…,θ(i),ki] that describes the most important information of the geometry of the xi.(3)
Global order of the local coordinates: supposing that the global coordinates ofxi converted by the Θi is Ti=[ti1,ti2,…,tiki], and then the error is
(7)Ei=Ti[I-(1k)eeT]-LiΘi,
where the I is the identity matrix; the e is the unit vector; the k is the points number of the neighborhood; the Li is the transformation matrix. In order to minimize the error, the Ti and Li should be found, and then
(8)Li=Ti(I-(1k)eeT)Θi*,Ei=Ti(I-(1k)eeT)(I-Θi*Θi),
where the Θi* is the Moor-Penrose generalized inverse of Θi. Suppose
(9)B=PWWTPT.
Let P=[P1,P2,…,PN], TPi=Ti, Pi be a selected matrix from 0-1; the T are global coordinates, and their weight matrix is
(10)W=diag⁡(W1,W2,…,WN),Wi=(I-(1k)eeT)(I-Θi*Θi).
The constraints is TTT=Id.(4)
Extract of the low-dimensional manifolds feature: since thee is the eigenvalue of matrix B, the corresponding minimum eigenvectors matrix is composed of eigenvalue. The value of section 2 to section d+1 of matrix B make of the T. T is the global coordinate mapping in the mainstream form of low-dimensional transformed from the nonlinear high-dimensional data set of X.The procedure of feature extraction can be described as follow.(1)
Use the time domain analysis methods kurtosis, skewness, peak-peak, RMS, and sample variance to extract the statistical features.(2)
Use the EMD method to decompose the collected vibration signal of each data set and get the IMF components; calculate the energy of each IMF component and get the features of the bearing in this time; then get the features of the other data sets.(3)
Use the LTSA to reduce the original features dimensions and get the main features; the extracted features are used as input of the SVM model for bearing degradation process prediction.


## 2.1. Feature Extraction

This section presents a brief discussion on original feature extraction from time domain, time-frequency domain of vibration signals. Time domain methods usually involve statistical features that are sensitive to impulsive oscillation, such as kurtosis, skewness, peak-peak (P-P), RMS, and sample variance. The 5 domain statistical features are used as original features in time domain has been used in the literatures [2, 3]:
(1)RMS=∑xi2N,
where N is the number of discrete points and xi represents the signal value at those points,
(2)Variance=s2=1(N-1)∑i=1N(xi-x-)2,
where x- is the mean value.Theith central moment for a set of data is defined as
(3)Mi=1N∑i=1N(xi-x-)i.The normalized forth moment, kurtosis, which is commonly used in bearing diagnostics, is defined as(4)Kurtosis=M4(M2)2-3.The skewness is defined as(5)Skewness=M3(M2)3/2.The peak-peak is defined as(6)peak-peak=max⁡(xi)-min⁡(xi).Empirical mode decomposition (EMD) is a powerful tool in time-frequency domain analysis. The advantage of EMD is the presentation of signals in time-frequency distribution diagrams with multiresolution, during which choosing some parameters is not needed. This property is essential in the detection of bearing faults. The EMD energy can represent the characteristic of vibration signals, and thus it is used as the input features. The (intrinsic mode function) IMF energy data sets are chosen as original features in this paper. The original features for bearing degradation process prediction based on the original features are shown in Table1.Table 1
Original features for bearing fault diagnosis and performance assessment.
Domain
Feature
Time domain
RMS, Kurtosis, Skewness, Peak-Peak, and sample variance
Time-frequency domain
EMD energy (Fe1,Fe2,Fe3,…,Fen)


## 2.2. Dimensional Reduction Based on the LTSA

Because the generated original feature sets are still with high dimension and include superfluous information, the feature extraction method LTSA is used to fuse the relevant useful features and extracts more sensitive features to work as the input of the proposed prediction model.The basic idea of LTSA is to use the tangent space of sample points to represent the geometry of the local character. Then these local manifold structures of space are lined up to construct the global coordinates. Given a data setX=[x1,x2,…,xN], xi∈Rm, a mainstream shape of d-dimension (m>d) is extracted. The LTSA feature extraction algorithm is as follows [9].(1)
Extract local information: for eachxi, i=1,2,…,N, use the Euclidean distance to determine a set xi=[xi,1,xi,2,…,xi,ki] of its neighborhood adjacent points (e.g., k the nearest neighbors).(2)
Local linear fitting: in the neighborhood of data pointsxi, a set of orthogonal basis Qi can be selected to construct the d-dimension neighborhood space of xi and the orthogonal projection of each point xi,j(j=1,2,…,N) can be calculated to the tangent space of θj(i)=QiT(xi,j-xi¯). xi¯ is the mean data for the neighborhood. The orthogonal projection in the tangent space of neighborhood data of xi is composed of local coordinate Θi=[θ(i),1,θ(i),2,…,θ(i),ki] that describes the most important information of the geometry of the xi.(3)
Global order of the local coordinates: supposing that the global coordinates ofxi converted by the Θi is Ti=[ti1,ti2,…,tiki], and then the error is
(7)Ei=Ti[I-(1k)eeT]-LiΘi,
where the I is the identity matrix; the e is the unit vector; the k is the points number of the neighborhood; the Li is the transformation matrix. In order to minimize the error, the Ti and Li should be found, and then
(8)Li=Ti(I-(1k)eeT)Θi*,Ei=Ti(I-(1k)eeT)(I-Θi*Θi),
where the Θi* is the Moor-Penrose generalized inverse of Θi. Suppose
(9)B=PWWTPT.
Let P=[P1,P2,…,PN], TPi=Ti, Pi be a selected matrix from 0-1; the T are global coordinates, and their weight matrix is
(10)W=diag⁡(W1,W2,…,WN),Wi=(I-(1k)eeT)(I-Θi*Θi).
The constraints is TTT=Id.(4)
Extract of the low-dimensional manifolds feature: since thee is the eigenvalue of matrix B, the corresponding minimum eigenvectors matrix is composed of eigenvalue. The value of section 2 to section d+1 of matrix B make of the T. T is the global coordinate mapping in the mainstream form of low-dimensional transformed from the nonlinear high-dimensional data set of X.The procedure of feature extraction can be described as follow.(1)
Use the time domain analysis methods kurtosis, skewness, peak-peak, RMS, and sample variance to extract the statistical features.(2)
Use the EMD method to decompose the collected vibration signal of each data set and get the IMF components; calculate the energy of each IMF component and get the features of the bearing in this time; then get the features of the other data sets.(3)
Use the LTSA to reduce the original features dimensions and get the main features; the extracted features are used as input of the SVM model for bearing degradation process prediction.


## 3. The SVM and Markov Model for Degradation Process Prediction

### 3.1. SVM Prediction Model

SVM is a machine learning tool that uses statistical learning theory to solve multidimensional functions. It is based on structural risk minimization principles, which overcomes the extralearning problem of ANN.The learning process of a SVM regression model is essentially a problem in quadratic programming. Given a set of data points(x1,y1)⋯(xk,yk) such that xi∈Rn as input and yi∈R as target output, the regression problem is to find a function such as
(11)y=f(x)=wϕ(x)+b,
where ϕ(x) is the high dimensional feature space, which is nonlinear mapped from the input space x, w is the weight vector, and b is the bias [22].After training, the correspondingy can be found through f(x) for the x outside the sample. The ε-support vector regression (ε-SVR) by Vapnik controls the precision of the algorithm through a specified tolerance error ε. The error of the sample is ξ, regardless of the loss, when |ξ|≤ε; else consider the loss as |ξ|-ε. First, map the sample into a high dimensional feature space by a nonlinear mapping function and convert the problem of the nonlinear function estimates into a linear regression problem in a high dimensional feature space. If we let φ(x) be the conversion function from the sample space into the high dimension feature space, then the problem of solving the parameters of f(x) is converted to solving an optimization problem (12) with the constraints in (13):(12)min⁡12∥w∥2=12(w·w)(13)00Subjecttoyi-(w·φ(x)+b)≤ε,(w·φ(x)+b)-yi≤ε,i=1,2,…,l.The feature space is one of high dimensionality and the target function is nondifferentiable. In general, the SVM regression problem is solved by establishing a Lagrange function and converting this problem to a dual optimization, that is, problem (14) with constraint of (15) in order to determine the Lagrange multipliers αi, αi*,
(14)max⁡-12∑i,j=1l(αi-αi*)(αj-αj*)(xi,xJ)-ε∑i=1l(αi+αi*)+∑i=1lyi(αi-αi*)(15)Subjectto∑i=1l(αi-αi*)=0αi,αi*∈[0,C],
where αi, αi* are Lagrange multipliers and αi,αi*≥0. αi×αi*=0. C evaluates the tradeoff between the empirical risk and the smoothness of the model.The SVM regression problem has therefore been transformed into a quadratic programming problem. The regression equation can be obtained by solving this problem. With the kernel functionK(xi,xj), the corresponding regression function is provided by
(16)f(x)=∑i=1N(αi-αi*)(αj-αj*)K(xi,xj)+b,
where the kernel function K(xi,xj) is an internal product of vectors xi and xj in feature spaces φ(xi) and φ(xj).


### 3.2. The Prediction Strategy and the Structure of the SVM Model

Traditional forecasting methods mainly achieve single-step prediction; when those methods are used for multisteps prediction, they cannot get an overall development trend of the series. Multisteps prediction method has the ability to obtain overall information of the series which provides the possibility for long-term prediction. There are two typical alternatives to build multisteps life prediction model. One is iterated prediction and the other is direct prediction. The comparison of the two strategies can be found in a number of literatures [23]. Marcellino et al. [24] presented a large-scale empirical comparison of iterated versus direct prediction. The results show that iterated prediction typically outperforms the direct prediction. So, the iterated multisteps prediction strategy has numerous advantages and will be adopted in this paper.In order to determine the structure of the SVM, we constructed a three layers SVM prediction model. But to achieve the multisteps time series life prediction a basic problem should be suppressed. That is how many essential observations (inputs) are used for forecasting the future value (the output node number is 1), so-called embedding dimensiond. In order to suppress the problem, the CAO method [25], which is particularly efficient to determine the minimum embedding dimension through the expansion of neighbor point in the embedding space, is employed to select an appropriate embedding dimension d. Then, the SVM input node number is determined.To effectively select an appropriate embedding dimension based on the CAO method, the phase space reconstruction method should be mentioned. The fundamental theorem of phase space reconstruction is pioneered by Takens [26]. For an N-point time series X={x1,x2,…,xN}, a sequence of vectors yi in a new space can be generated as yi(d)={xi,xi+τ,…,xi+(d-1)τ}, where i=1,2,…,Nd, Nd=N-(d-1)τ is the length of the reconstructed vector yi, d is the embedding dimension of the reconstructed state space, and τ is embedding delay time. The time delay τ is chosen through the autocorrelation function [27]:
(17)C(τ)=∑i=1N-τxi′xi+τ′∑i=1N-τ(xi′)2,
where xi′=xi-x¯, x¯ is the average value of the time series. The optimal time delay τ is determined when the first minimum value of C(τ) occurs.The embedding dimensiond is chosen through CAO method, defining the quantity as follows:
(18)a(i,d)=∥yi(d+1)-yn(i,d)(d+1)∥∥yi(d)-yn(i,d)(d)∥,
where ∥·∥ is the Euclidian distance and is given by the maximum norm. yi(d) means the ith reconstructed vector and n(i,d) is an integer, so that yn(i,d)(d) is the nearest neighbor of yi(d) in the embedding dimension d. A new quantity is defined as the mean value of all a(i,d)′s:
(19)E(d)=1N-dτ∑i=1N-dτa(i,d),where  E(d) is only dependent on the dimension d and time delay τ. To investigate its variation from d to d+1, the parameter E1 is given by
(20)E1(d)=E(d+1)E(d).By increasing the value ofd, the value E1(d) is also increased and it stops increasing when the time series comes from a deterministic process. If a plateau is observed for d≥d0, then d0+1 is the minimum embedding dimension. But E1(d) has the problem of slowly increasing or has stopped changing if d is sufficiently large. CAO introduced another quantity E2(d) to overcome the problem:
(21)E2(d)=E*(d+1)E*(d),
where
(22)E*(d)=1N-dτ∑i=1N-dτ|xi+dτ-xn(i,d)+dτ|.Through CAO method, the embedding dimensiond of the SVM prediction model is chosen. The structure of the SVM model is determined.


### 3.3. SOM Clustering Method to Divide the Prediction Error

State division is the process to determine the mapping from random variable to the state space. How to obtain state division is a crux for Markov model. Traditionally, it is performed by the state division approach described as follows. LetX={x1,x2,…,xn} be the random sequence; let S={1,2,…,j,…,m} denote the state space; given a0<a1<⋯<am if random variable xi∈[aj-1,aj], where 1≤i≤n, 1≤j≤m, then the variable xi belongs to the state j, and the division of [aj-1,aj] is usually uniform divided. However, the uniform divided method depends on the people’s experience, which will affect the prediction precise. In this research, the SOM neural network [28] is used to divide the state. The SOM can be created from highly deviating, nonlinear data. After the data are input, the SOM is trained iteratively.In each training step, one sample vectorX from the input data set is chosen randomly, and the distance between it and all the weight vectors of the SOM, which are originally initialised randomly, is calculated using some distance measure. The best matching unit (BMU) is the map unit, whose weight vector is closest to X. After the BMU is identified, the weight vectors of the BMU, as well as its topological neighbors, are updated so that they are moved closer to the input vector in the input space. The vectors are updated following the learning rule:
(23)mi(t+1)=mi(t)+α(t)h(nBMU,ni,t)(X-mi(t)),
where h(nBMU,ni,t) is the neighborhood function, which is monotonically decreasing with respect to the distance between the BMU nBMU and ni in the grid, and the training time α(t) is the learning rate; a decreasing function with 0<α(t)<1.At the end of the learning process, the weight vectors are grouped into clusters depending on their distance in the input space. Unlike networks based on supervised learning, which require that target values corresponding to input vectors are known, the SOM can be used to cluster data without knowing the class membership of the input data. This character is suitable for the problem of the prediction error of the SVM model which is not clear to us, so we should classify the error without prelearning, and therefore, the function of SOM method determined is an efficient and necessary method for clustering the state. Based on the clustering results, the state is divided into some districts.


### 3.4. Markov Prediction Model to Improve the Prediction Accuracy

Consider a stochastic process{xn,n=1,2,…} that takes on a finite or countable number of possible values. Unless otherwise mentioned, this set of possible values of the process will be denoted by the set of nonnegative integers {1,2,…}. If xn=i, then the process is said to be in state i at time n. Suppose that whenever the process is in state i; there is a fixed probability pij that it will be next instating j. That is
(24)p{xn+1=j∣xn=i,xn-1=in-1,…,x1=i1}=p{xn+1=j∣xn=i}=pij
for all states i0,i1,…,in-1,i,j and all n≥0. Such a stochastic process is known as a Markov chain. Equation (24) can be interpreted as stating that for a Markov model, the conditional distribution of any state xn+1, given the past states x1,x2,…,xn-1 and the present state xn, is independent of the past states and depends only on the present state. This is called the Markovian property. The value pij represents the probability that the process will, when in state i, next make a transition into state j. Since probabilities are nonnegative and since the process must make a transition into some state, we have that
(25)pij≥0,i,j≥0;∑j=1∞pij=1,i=1,2,….If the process has a finite number of states, which means the state spaceS={1,2,…,i,j,…,s}, then the Markov chain model can be defined by the matrix of one-step transition probabilities, denoted as
(26)p=[p11p12⋯p1sp21p22⋯p2s⋮⋮⋯⋮ps1ps2⋯pss].The initial probability is computed by(27)pij=NijNi,
where Nij denotes the transition times from state i to state j and Ni denotes the number of random variables {xn,n=1,2,…,m} belonging to state i.Markov model adopts state vector and state transition matrix to deal with the prediction issue. Suppose that the state vector of momentt-1 is pt-1, the state vector of moment t is pt, and the state transition matrix is p; then the relationship is
(28)pt=pt-1p,t=1,2,…,K.Updatet from 1 to K, and then
(29)pk=p0pK,
where pK is the state vector at moment K. Equation (29) is the basic Markov prediction model, if the initial state vector and the transition matrix are given, which allows calculation of any possible future state vector.


## 3.1. SVM Prediction Model

SVM is a machine learning tool that uses statistical learning theory to solve multidimensional functions. It is based on structural risk minimization principles, which overcomes the extralearning problem of ANN.The learning process of a SVM regression model is essentially a problem in quadratic programming. Given a set of data points(x1,y1)⋯(xk,yk) such that xi∈Rn as input and yi∈R as target output, the regression problem is to find a function such as
(11)y=f(x)=wϕ(x)+b,
where ϕ(x) is the high dimensional feature space, which is nonlinear mapped from the input space x, w is the weight vector, and b is the bias [22].After training, the correspondingy can be found through f(x) for the x outside the sample. The ε-support vector regression (ε-SVR) by Vapnik controls the precision of the algorithm through a specified tolerance error ε. The error of the sample is ξ, regardless of the loss, when |ξ|≤ε; else consider the loss as |ξ|-ε. First, map the sample into a high dimensional feature space by a nonlinear mapping function and convert the problem of the nonlinear function estimates into a linear regression problem in a high dimensional feature space. If we let φ(x) be the conversion function from the sample space into the high dimension feature space, then the problem of solving the parameters of f(x) is converted to solving an optimization problem (12) with the constraints in (13):(12)min⁡12∥w∥2=12(w·w)(13)00Subjecttoyi-(w·φ(x)+b)≤ε,(w·φ(x)+b)-yi≤ε,i=1,2,…,l.The feature space is one of high dimensionality and the target function is nondifferentiable. In general, the SVM regression problem is solved by establishing a Lagrange function and converting this problem to a dual optimization, that is, problem (14) with constraint of (15) in order to determine the Lagrange multipliers αi, αi*,
(14)max⁡-12∑i,j=1l(αi-αi*)(αj-αj*)(xi,xJ)-ε∑i=1l(αi+αi*)+∑i=1lyi(αi-αi*)(15)Subjectto∑i=1l(αi-αi*)=0αi,αi*∈[0,C],
where αi, αi* are Lagrange multipliers and αi,αi*≥0. αi×αi*=0. C evaluates the tradeoff between the empirical risk and the smoothness of the model.The SVM regression problem has therefore been transformed into a quadratic programming problem. The regression equation can be obtained by solving this problem. With the kernel functionK(xi,xj), the corresponding regression function is provided by
(16)f(x)=∑i=1N(αi-αi*)(αj-αj*)K(xi,xj)+b,
where the kernel function K(xi,xj) is an internal product of vectors xi and xj in feature spaces φ(xi) and φ(xj).


## 3.2. The Prediction Strategy and the Structure of the SVM Model

Traditional forecasting methods mainly achieve single-step prediction; when those methods are used for multisteps prediction, they cannot get an overall development trend of the series. Multisteps prediction method has the ability to obtain overall information of the series which provides the possibility for long-term prediction. There are two typical alternatives to build multisteps life prediction model. One is iterated prediction and the other is direct prediction. The comparison of the two strategies can be found in a number of literatures [23]. Marcellino et al. [24] presented a large-scale empirical comparison of iterated versus direct prediction. The results show that iterated prediction typically outperforms the direct prediction. So, the iterated multisteps prediction strategy has numerous advantages and will be adopted in this paper.In order to determine the structure of the SVM, we constructed a three layers SVM prediction model. But to achieve the multisteps time series life prediction a basic problem should be suppressed. That is how many essential observations (inputs) are used for forecasting the future value (the output node number is 1), so-called embedding dimensiond. In order to suppress the problem, the CAO method [25], which is particularly efficient to determine the minimum embedding dimension through the expansion of neighbor point in the embedding space, is employed to select an appropriate embedding dimension d. Then, the SVM input node number is determined.To effectively select an appropriate embedding dimension based on the CAO method, the phase space reconstruction method should be mentioned. The fundamental theorem of phase space reconstruction is pioneered by Takens [26]. For an N-point time series X={x1,x2,…,xN}, a sequence of vectors yi in a new space can be generated as yi(d)={xi,xi+τ,…,xi+(d-1)τ}, where i=1,2,…,Nd, Nd=N-(d-1)τ is the length of the reconstructed vector yi, d is the embedding dimension of the reconstructed state space, and τ is embedding delay time. The time delay τ is chosen through the autocorrelation function [27]:
(17)C(τ)=∑i=1N-τxi′xi+τ′∑i=1N-τ(xi′)2,
where xi′=xi-x¯, x¯ is the average value of the time series. The optimal time delay τ is determined when the first minimum value of C(τ) occurs.The embedding dimensiond is chosen through CAO method, defining the quantity as follows:
(18)a(i,d)=∥yi(d+1)-yn(i,d)(d+1)∥∥yi(d)-yn(i,d)(d)∥,
where ∥·∥ is the Euclidian distance and is given by the maximum norm. yi(d) means the ith reconstructed vector and n(i,d) is an integer, so that yn(i,d)(d) is the nearest neighbor of yi(d) in the embedding dimension d. A new quantity is defined as the mean value of all a(i,d)′s:
(19)E(d)=1N-dτ∑i=1N-dτa(i,d),where  E(d) is only dependent on the dimension d and time delay τ. To investigate its variation from d to d+1, the parameter E1 is given by
(20)E1(d)=E(d+1)E(d).By increasing the value ofd, the value E1(d) is also increased and it stops increasing when the time series comes from a deterministic process. If a plateau is observed for d≥d0, then d0+1 is the minimum embedding dimension. But E1(d) has the problem of slowly increasing or has stopped changing if d is sufficiently large. CAO introduced another quantity E2(d) to overcome the problem:
(21)E2(d)=E*(d+1)E*(d),
where
(22)E*(d)=1N-dτ∑i=1N-dτ|xi+dτ-xn(i,d)+dτ|.Through CAO method, the embedding dimensiond of the SVM prediction model is chosen. The structure of the SVM model is determined.


## 3.3. SOM Clustering Method to Divide the Prediction Error

State division is the process to determine the mapping from random variable to the state space. How to obtain state division is a crux for Markov model. Traditionally, it is performed by the state division approach described as follows. LetX={x1,x2,…,xn} be the random sequence; let S={1,2,…,j,…,m} denote the state space; given a0<a1<⋯<am if random variable xi∈[aj-1,aj], where 1≤i≤n, 1≤j≤m, then the variable xi belongs to the state j, and the division of [aj-1,aj] is usually uniform divided. However, the uniform divided method depends on the people’s experience, which will affect the prediction precise. In this research, the SOM neural network [28] is used to divide the state. The SOM can be created from highly deviating, nonlinear data. After the data are input, the SOM is trained iteratively.In each training step, one sample vectorX from the input data set is chosen randomly, and the distance between it and all the weight vectors of the SOM, which are originally initialised randomly, is calculated using some distance measure. The best matching unit (BMU) is the map unit, whose weight vector is closest to X. After the BMU is identified, the weight vectors of the BMU, as well as its topological neighbors, are updated so that they are moved closer to the input vector in the input space. The vectors are updated following the learning rule:
(23)mi(t+1)=mi(t)+α(t)h(nBMU,ni,t)(X-mi(t)),
where h(nBMU,ni,t) is the neighborhood function, which is monotonically decreasing with respect to the distance between the BMU nBMU and ni in the grid, and the training time α(t) is the learning rate; a decreasing function with 0<α(t)<1.At the end of the learning process, the weight vectors are grouped into clusters depending on their distance in the input space. Unlike networks based on supervised learning, which require that target values corresponding to input vectors are known, the SOM can be used to cluster data without knowing the class membership of the input data. This character is suitable for the problem of the prediction error of the SVM model which is not clear to us, so we should classify the error without prelearning, and therefore, the function of SOM method determined is an efficient and necessary method for clustering the state. Based on the clustering results, the state is divided into some districts.


## 3.4. Markov Prediction Model to Improve the Prediction Accuracy

Consider a stochastic process{xn,n=1,2,…} that takes on a finite or countable number of possible values. Unless otherwise mentioned, this set of possible values of the process will be denoted by the set of nonnegative integers {1,2,…}. If xn=i, then the process is said to be in state i at time n. Suppose that whenever the process is in state i; there is a fixed probability pij that it will be next instating j. That is
(24)p{xn+1=j∣xn=i,xn-1=in-1,…,x1=i1}=p{xn+1=j∣xn=i}=pij
for all states i0,i1,…,in-1,i,j and all n≥0. Such a stochastic process is known as a Markov chain. Equation (24) can be interpreted as stating that for a Markov model, the conditional distribution of any state xn+1, given the past states x1,x2,…,xn-1 and the present state xn, is independent of the past states and depends only on the present state. This is called the Markovian property. The value pij represents the probability that the process will, when in state i, next make a transition into state j. Since probabilities are nonnegative and since the process must make a transition into some state, we have that
(25)pij≥0,i,j≥0;∑j=1∞pij=1,i=1,2,….If the process has a finite number of states, which means the state spaceS={1,2,…,i,j,…,s}, then the Markov chain model can be defined by the matrix of one-step transition probabilities, denoted as
(26)p=[p11p12⋯p1sp21p22⋯p2s⋮⋮⋯⋮ps1ps2⋯pss].The initial probability is computed by(27)pij=NijNi,
where Nij denotes the transition times from state i to state j and Ni denotes the number of random variables {xn,n=1,2,…,m} belonging to state i.Markov model adopts state vector and state transition matrix to deal with the prediction issue. Suppose that the state vector of momentt-1 is pt-1, the state vector of moment t is pt, and the state transition matrix is p; then the relationship is
(28)pt=pt-1p,t=1,2,…,K.Updatet from 1 to K, and then
(29)pk=p0pK,
where pK is the state vector at moment K. Equation (29) is the basic Markov prediction model, if the initial state vector and the transition matrix are given, which allows calculation of any possible future state vector.


## 4. Proposed Method

The flowchart of the proposed method is shown in Figure1. The method consists of four procedures sequentially: data processing and features extraction, merge of the original features, constructing-training SVM model and predicting, and Markov model for improving the prediction result. The role of each procedure is explained as follows.Figure 1
The flowchart of the proposed method.Step 1.
Data processing and features extraction. The time domain and time-frequency domain signal processing methods are used to extract the original features from the collected mass vibration data.Step 2.
Merge of the original features. The LTSA method is used to extract the typical features and reduce the dimension of the features. The extracted features are used for training the SVM model.Step 3.
Constructing the SVM model. The SVM model is constructed; the CAO method is used to determine the embedding dimension. The iterated multistep prediction method is used to forecast the future value.Step 4.
Markov model for improving the prediction result. This procedure uses the SOM method to cluster the prediction error before the Markov method; based on the state division results the Markov model is used to improve the prediction results obtained by SVM model, to get a more precise prediction.


## 5. Validation and Application

### 5.1. Validation

In order to validate the effect of the proposed method, a validation test is proposed. The vibration signals used in this paper are provided by the Center for Intelligent Maintenance Systems (IMS), University of Cincinnati [29]. The experimental data sets are generated from bearing run-to-failure tests under constant load conditions on a specially designed test rig as shown in Figure 2. The rotation speed is kept constant at 2000 rpm. A radial load of 6000 lbs is added to the shaft and bearing by a spring mechanism. The data sampling rate is 20 kHz and the data length is 20,480 points as shown in Figure 3. It took a total of 7 days until the bearing fails. At the end, one bearing with serious wreck is used to test the proposed method as shown in Figure 4.Figure 2
Bearing test tig and sensor placement illustration [29].Figure 3
Vibration signal.The bearing components with serious wrecked after test with roller element defect and outer race defect.
(a)
(b)The time domain and time-frequency domain methods are used to deal with the collected vibration data as described in Section2, Table 1. The measurements value of kurtosis and skewness are depicted in Figures 5(a) and 5(b).The IMF1 and IMF2 energy are depicted in Figures 6(a) and 6(b).(a) Kurtosis value and (b) Skewness value of bearing vibration data.
(a)
(b)(a) IMF1 energy of bearing vibration data; (b) IMF2 energy of bearing vibration data.
(a)
(b)From the extracted features showed in Figures5 and 6 we can see the following. (1) The bearing is in normal condition during the time correlated with the first 700 points. After that time, the condition of bearing suddenly changes. It indicates that there are some faults occurring in this bearing. (2) Different features reflect the bearing running state in different shapes. For example, the kurtosis and the IMF1 indict the bearing running state with an upward trend, while the skewness indicts the bearing running state with a downward trend. The kurtosis, skewness, and the IMF1, indicting the bearing running state, have a large fluctuation after the 700-point, while the IMF2, indicting the bearing running state, have a large fluctuation till after the 850-point. So using these features is not appropriate to reflect the bearing running condition. So it is very necessary to extract a sensitive feature by the original features to appropriately reflect the bearing running condition.We performed feature extraction by means of LTSA to extract a sensitive feature and reduce the dimensionality of calculated features. After LTSA is used (in this article the parameters of the neighborhood factork equals 8, the embedding dimension d equals 1), the bearing running state features dataset is got. The first main projected vector is chosen as the input of the SVM model. The result is shown in Figure 7. In comparison with the LTSA method, we also extracted the features through the PCA method, and the first main principal component is chosen. The result is shown in Figure 8.Figure 7
The first main projected vector of test bearing based on LTSA method.Figure 8
The first principal component of test bearing based on PCA method.From Figures7 and 8 we can see that the LTSA method can extract an effective feature dataset, which is sensitive to the changes of bearing running state, while the extracted feature is based on PCA method with a bad effect, before the 700 point; we even cannot see the fluctuation of the bearing running state, and the trend convert is also not obvious, from which we cannot know the bearing running state effectively. This result indicates that information extracted by LTSA could be more effective than that extracted by PCA.After extracting the typical features, the CAO method is used to determine the embedding dimension of the SVM model, based on the theorem of phase space reconstruction. We first choose the delay timeτ through the autocorrelation function.The optimal time delayτ is determined when the first minimum value of C(τ) occurs.Based on the extracted features dataset, the delay timeτ is set to 3 for the projected vector values through the autocorrelation function as shown in Figure 9.Figure 9
Selection of the delay time by the autocorrelation function value of LTSA1.Then the embedding dimension is selected by the CAO method. The result is shown in Figure10; the optimal embedding dimension d for the projected vector is chosen as 10.Figure 10
Selection of the embedding dimension by the CAO method of LTSA1.Based on the selected optimal embedding dimension, the SVM model is used to achieve the multisteps prediction. The RBF kernel function is used:(30)K(xi,xj)=exp⁡(-12σ2∥x-xi∥2).In this research the regularity parameterC is set to 90.3, the kernel function parameter σ2 is set to 20, and the ε is set to 0.001. The parameters C and ε are selected by Particle Swarm Optimization algorithm (PSO) [30]. The popular size of the PSO is set to 100, the interaction number of the PSO method is set to 20, and the fitness function of the PSO method is set to choose the parameters which make the SVM model fitting error in the training process the smallest. The error goal is set to 0.05, the dimension of the PSO is set to 2, and the inertia weight is set to w=0.5, c1=c2=1.2. Based on the selected typical features, the features dataset is used to train SVM model and the input features number of SVM is 9 determined by embedding dimension. Then, the trained SVM model is used to predict the bearing running state. Before the 700th points, the bearing is working in a normal state, so the 701–900 points are used to train the SVM model and the following 85 points are employed for testing. In order to evaluate the predicting performance, the root-mean square error (RMSE) is utilized as follows:
(31)RMSE=∑i=1N(yi-y^i)2N,
where N represents the total number of data points in the test set, yi is actual value in training set or test set, and y^i represents the predicted value of the model. The actual value and the predicted result are shown in Figure 11.Figure 11
Prediction result based on the SVM model.From Figure11 we can see that the trend of the bearing running state can be predicted by SVM model. From the prediction, we can get a general understanding of the bearing running state in the future, but the predicted result is not accurate, especially the stage of 70–85 point, through calculation. The RMSE of the actual and the predicted result is 0.0469, so the predict result is not satisfied and the prediction error of SVM model (the predicted results subtract the actual data) is shown in Figure 12.Figure 12
The prediction error of the SVM model.Then the SOM method is used to divide the error into some districts, the iteration number if SOM is 2000; the structure of the state classification matrix is[3×1]. The results of the state division by SOM are shown in Table 2.Table 2
The result of the state division and the prediction results.
Prediction number
Prediction error
State
Prediction number
Prediction error
State
Prediction number
Prediction error
State
1
−0.001622
2
31
−0.014505
2
61
−0.009996
2
2
−0.003811
2
32
−0.017872
2
62
−0.015532
2
3
−0.003841
2
33
−0.022864
2
63
−0.041807
1
4
−0.007212
2
34
−0.031673
1
64
−0.022312
2
5
−0.008162
2
35
−0.033233
1
65
−0.034537
1
6
−0.02057
2
36
−0.05774
1
66
−0.060772
1
7
−0.023467
1
37
−0.090116
1
67
−0.10752
1
8
−0.026518
1
38
−0.04199
1
68
−0.096034
1
9
−0.02786
1
39
−0.034686
1
69
−0.12397
1
10
−0.027196
1
40
−0.032578
1
70
−0.13711
1
11
−0.038575
1
41
−0.010655
2
71
−0.13712
1
12
−0.025033
1
42
−0.013265
2
72
−0.03195
1
13
−0.018556
2
43
−0.006294
2
73
−0.02958
1
14
−0.019786
2
44
−0.007407
2
74
−0.19645
1
15
−0.01983
2
45
−0.021696
2
75
−0.18364
1
16
−0.020994
2
46
−0.017219
2
76
0.24355
3
17
−0.032989
1
47
−0.037728
1
77
0.43612
3
18
−0.019086
2
48
−0.011678
2
78
0.20092
3
19
−0.013994
2
49
−0.001478
2
79
−0.05844
1
20
−0.014597
2
50
−0.003989
2
80
0.76262
3
21
−0.013641
2
51
−0.00038
2
81
0.14341
3
22
−0.014209
2
52
−0.00578
2
82
0.28256
3
23
−0.013578
2
53
−0.025056
1
83
0.1326
3
24
−0.022126
2
54
−0.045813
1
84
0.3338
3
25
−0.019049
2
55
−0.01996
2
85
0.18377
3
26
−0.013825
2
56
−0.05286
1
27
−0.012682
2
57
−0.019165
2
28
−0.014012
2
58
−0.039591
1
29
−0.012089
2
59
−0.015265
2
30
−0.01567
2
60
−0.058697
1It can be seen from Table2 that when there is downward trend or the point value is less than 0, the state is set to 1 and when there is upward trend, the state is set to 3. According to the classification of the SOM model, the Markov state is divided into the following districts: [-0.10752,-0.00038], [-0.00038,0.1326], and [0.1326,0.43612]. The Markov model is used to improve the prediction error. For example, the 78 point, the Markov state transition matrix from 1–77 point is
(32)p=[0.6970.27270.03030.23260.76740001].Chose three points 77, 76, and 75 which are recent to the 78 point and set the transfer step as 1, 2, 3; the state prediction results based on the Markov model are shown in Table3.Table 3
Table of the probability status.
Point
Status probability
Transfer step
State 1
State 2
State 3
77
3
1
0
0
1
76
3
2
0
0
1
75
1
3
0.4757
0.4562
0.0681
Total
0.4757
0.4562
2.0681From Table3 we can see that the accumulated value of stage 3 is the largest, so the stage of 78 point is set to 3 and the result is the same with the SOM clustering method.In this research, according to the stage of the prediction of the Markov model, the correct value is calculated byx~(0)(k)=x^(0)(k)-ε, where ε is the median value of the divided stage area and x^(0)(k) is the value predicted though SVM model.For the 78 point, the corrected value is0.57298-(0.1326+0.43612)/2=0.2956, where the actual value at this point is 0.37206. Then other points have also been corrected though this method. The corrected results of the point 70–85 are shown in Table 4.Table 4
The corrected results of the points 70–85 based on Markov model.
Point
SVM model prediction error
Relative error (%)
SVM-Markov model prediction error
Relative error (%)
70
−0.13711
84.37%
−0.1025
63.07%
71
−0.13712
79.91%
−0.1235
71.97%
72
−0.03195
16.92%
−0.0121
6.41%
73
−0.02958
12.85%
0.0965
41.92%
74
−0.19645
84.44%
−0.09512
40.88%
75
−0.18364
84.21%
−0.10236
46.94%
76
0.24355
93.60%
0.01043
4.01%
77
0.43612
79.88%
0.03542
6.49%
78
0.20092
54.00%
−0.07646
20.55%
79
−0.05844
23.60%
−0.01574
6.36%
80
0.76262
344.98%
0.13262
59.99%
81
0.14341
26.80%
0.03562
6.66%
82
0.28256
114.19%
0.20302
82.05%
83
0.1326
49.59%
−0.00251
0.94%
84
0.3338
2061%
0.03291
203.20%
85
0.18377
1132.15%
0.06548
403.40%From Table4 we can see that the Markov model makes the results more precise, which validate the necessary to use the Markov model to improve the effect of the proposed method. The RMSE of the actual and the predicted result is 0.0091, so the prediction accuracy improved significantly. However, because those points are so far away compared to the predicted point of 1–69, the results still have some error.In order to compare the predict effect, the most usually used prediction model BP neural networks is used to predict the bearing running state based on the selected features. The learning rate of the neural network and its momentum coefficient are 0.01; the weights are initialized to uniformly distribute random values between −0.1 and 0.1; the iteration number is 2000; the training error is 0.001; the input number is 9; the hidden number is 15; the output node number is 1. The prediction results are shown in Figure13.Figure 13
Prediction result based on the BPNN model.From Figure13 we can see that the prediction results based on the BPNN model is not working effectively. There are some peaks while in the same position the actual status is not obvious. The RMSE of the predicted result is 0.0932, so the prediction results of the traditional BPNN model is not more effective than the SVM model. In addition, the prediction method based on the BPNN has the problem of prediction results which are unstable; when the same data is used to train and predict, the results are different and even the neural network may fall into the local optimum as shown in Figures 14 and 15.Figure 14
The training process curve of the BPNN method.Figure 15
The actual data and the result predicted by BPNN fall into the local optimum.With the data, the proposed method has also been compared with other methods that had been proposed in relative research.(1) The principal signal features extracted by PCA are utilized by HMM to predict the bearing running state [31]. (2) The time domain and frequency domain features have been directly used as the input of the prediction model, and the result has been predicted by Neural Network algorithm [1]. (3) The original features have been extracted by PCA as the input of the SVM prediction model [32]. (4) The proposed method in this research. The RMSE of the different methods predicted results is shown in Table 5.Table 5
The RMSE results of different prediction methods.
Method
The features extracted by PCA are utilized by HMM prediction model
The features directly used as the input of the Neural Network prediction model
The features extracted by PCA as the input of the SVM prediction model
The proposed method
RMSE
0.0698
0.1023
0.0521
0.0091From Table5, we can see that the RMSE of different prediction methods is very different. The prediction method that the original features have been directly used as the input of the NN model works the worst. This is the reason why the original features are still with high dimension and include superfluous information, which is not appropriate for state prediction; in addition, the NN prediction model has the drawbacks of slow convergence and difficulty in escaping from local minima. The prediction method based on the HMM model works not more effective than the method based on the SVM model; that is the reason why the HMM is not appropriate for long time forecast. The proposed method works the best; this is because the LTSA features extraction method can effectively extract the typical features and reduce the dimension and the SVM-Markov model can predict the state more precisely than the SVM and Markov model only. So through the comparison we can get that the proposed method is very effective in bearing running state prediction.


### 5.2. Application

After validating the effectiveness of the proposed method, the method has been used to the actual application. The test rig is shown in Figure16.Figure 16
Bearing test tig.The bearings are hosted on the shaft and the shaft is driven by AC motor. The rotation speed is kept at 1000 rpm and a radial load of 3 kg is added to the bearing. The data sampling rate is 25600 Hz and the data length is 102400 points collected on the date of 2011.11.25 as shown in Figure17. Every 2 hours, the vibration data are collected for one time. The collected data from 2011.11.25 to 2011.12.17 are analyzed after running for 1 year.Figure 17
Vibration signal.The time domain and time-frequency domain methods are used to deal with the collected vibration data as described in Section2, Table 1. Then the features are normalized through x⌢=(x-xmin⁡)/(xmax⁡-xmin⁡) and processed into the interval [0,1]. The LTSA is used to reduce the dimensionality of calculated features and the result is shown in Figure 18.Figure 18
The first main projected vector of test bearing based on LTSA method.From Figure18, we can see that the bearing running state has a fluctuation and upward trend. Especially at 150 points, there is a sudden change of trend, which reflects the bearing’s working status change at this moment.After extracting the typical features, the CAO method is used to determine the embedding dimension of the SVM model. The delay timeτ is set as 2 for the projected vector values though the autocorrelation function as shown in Figure 19.Figure 19
Selection of the delay time by the autocorrelation function value of LTSA1.The embedding dimension is selected by the CAO method. The result is shown in Figure20 and the optimal embedding dimension d for the projected vector is chosen as 10.Figure 20
Selection of the embedding dimension by the CAO method of LTSA1.Based on the selected optimal embedding dimension, the SVM model is used to achieve the prediction. The regularity parameterC is set as 909.5 and the ε is set to 0.01 selected by PSO method. Based on the selected typical features, the features dataset is used to train SVM model, the 1–130 points are used to train the SVM model, and the following 20 points are employed for testing. The prediction results are shown in Figure 21.Figure 21
Prediction result based on the SVM model.From Figure21 we can see that the trend of the bearing running state can be predicted by SVM model; from the prediction, we can get a general upward trend similar to the actual status. In addition, the sudden change of points 8 to 12 (near the 150 points in original signal as mentioned in Figure 17) is also showed out. However, the results are still not precise. In order to improve the prediction effect, the SOM method is used to divided the prediction error into some districts, the iteration number if SOM is 1000; the structure of the state classification matrix is [3×1]. The results of the state division by SOM are shown in Table 6.Table 6
The result of the state division and the prediction results.
Point
Prediction error
State
Point
Prediction error
State
1
0
3
11
0.053935
3
2
−0.063821
1
12
−0.0052358
1
3
0.002906
3
13
−0.058422
1
4
0.0074271
3
14
−0.094039
1
5
−0.021835
2
15
−0.032173
2
6
0.024138
3
16
−0.034682
2
7
0.0059543
3
17
−0.014101
2
8
0.016244
3
18
−0.13365
1
9
0.022316
3
19
−0.06811
1
10
0.016692
3
20
−0.070675
1Based on the classification of the SOM model, the Markov state is divided into the following districts[-0.13365,-0.034682], [-0.034682,0], and [0,0.053935]; the Markov model is used to improve the prediction error. The corrected results of the points 1–20 are shown in Table 7.Table 7
The corrected results of the points 1−20 based on Markov model.
Point
SVM model prediction error
Relative error (%)
SVM-Markov model prediction error
Relative error (%)
1
0
0
0
0
2
−0.0103821
579.21%
−0.006361
178.31%
3
0.002906
8.5739%
0.000603
1.7791%
4
0.0074271
20.311%
−0.001021
2.7922%
5
−0.021835
274.63%
−0.009362
117.75%
6
0.024138
30.429%
0.013254
16.708%
7
0.0059543
12.666%
0.002615
5.5628%
8
0.016244
22.459%
0.010027
13.863%
9
0.022316
33.565%
0.009565
14.386%
10
0.016692
32.417%
0.009473
18.397%
11
0.053935
73.825%
0.021563
29.515%
12
−0.0052358
34.947%
−0.002355
15.719%
13
−0.058422
107.06%
−0.019347
35.455%
14
−0.094039
122.79%
−0.014236
18.589%
15
−0.032173
263.1%
−0.009878
80.779%
16
−0.034682
159.3%
−0.008969
41.197%
17
−0.004101
355.3%
−0.001892
199.5%
18
−0.13365
117.18%
−0.050841
44.576%
19
−0.06811
120.92%
−0.028323
50.282%
20
−0.070675
135.04%
−0.014267
27.152%From Table7, we can see that the Markov model make the results more precise.Through the validation and actual application result we can see that the proposed method can predict the future status of the bearing, which is necessary for us to make some plan and do maintenance to reduce the risk of unnecessary accident.


## 5.1. Validation

In order to validate the effect of the proposed method, a validation test is proposed. The vibration signals used in this paper are provided by the Center for Intelligent Maintenance Systems (IMS), University of Cincinnati [29]. The experimental data sets are generated from bearing run-to-failure tests under constant load conditions on a specially designed test rig as shown in Figure 2. The rotation speed is kept constant at 2000 rpm. A radial load of 6000 lbs is added to the shaft and bearing by a spring mechanism. The data sampling rate is 20 kHz and the data length is 20,480 points as shown in Figure 3. It took a total of 7 days until the bearing fails. At the end, one bearing with serious wreck is used to test the proposed method as shown in Figure 4.Figure 2
Bearing test tig and sensor placement illustration [29].Figure 3
Vibration signal.The bearing components with serious wrecked after test with roller element defect and outer race defect.
(a)
(b)The time domain and time-frequency domain methods are used to deal with the collected vibration data as described in Section2, Table 1. The measurements value of kurtosis and skewness are depicted in Figures 5(a) and 5(b).The IMF1 and IMF2 energy are depicted in Figures 6(a) and 6(b).(a) Kurtosis value and (b) Skewness value of bearing vibration data.
(a)
(b)(a) IMF1 energy of bearing vibration data; (b) IMF2 energy of bearing vibration data.
(a)
(b)From the extracted features showed in Figures5 and 6 we can see the following. (1) The bearing is in normal condition during the time correlated with the first 700 points. After that time, the condition of bearing suddenly changes. It indicates that there are some faults occurring in this bearing. (2) Different features reflect the bearing running state in different shapes. For example, the kurtosis and the IMF1 indict the bearing running state with an upward trend, while the skewness indicts the bearing running state with a downward trend. The kurtosis, skewness, and the IMF1, indicting the bearing running state, have a large fluctuation after the 700-point, while the IMF2, indicting the bearing running state, have a large fluctuation till after the 850-point. So using these features is not appropriate to reflect the bearing running condition. So it is very necessary to extract a sensitive feature by the original features to appropriately reflect the bearing running condition.We performed feature extraction by means of LTSA to extract a sensitive feature and reduce the dimensionality of calculated features. After LTSA is used (in this article the parameters of the neighborhood factork equals 8, the embedding dimension d equals 1), the bearing running state features dataset is got. The first main projected vector is chosen as the input of the SVM model. The result is shown in Figure 7. In comparison with the LTSA method, we also extracted the features through the PCA method, and the first main principal component is chosen. The result is shown in Figure 8.Figure 7
The first main projected vector of test bearing based on LTSA method.Figure 8
The first principal component of test bearing based on PCA method.From Figures7 and 8 we can see that the LTSA method can extract an effective feature dataset, which is sensitive to the changes of bearing running state, while the extracted feature is based on PCA method with a bad effect, before the 700 point; we even cannot see the fluctuation of the bearing running state, and the trend convert is also not obvious, from which we cannot know the bearing running state effectively. This result indicates that information extracted by LTSA could be more effective than that extracted by PCA.After extracting the typical features, the CAO method is used to determine the embedding dimension of the SVM model, based on the theorem of phase space reconstruction. We first choose the delay timeτ through the autocorrelation function.The optimal time delayτ is determined when the first minimum value of C(τ) occurs.Based on the extracted features dataset, the delay timeτ is set to 3 for the projected vector values through the autocorrelation function as shown in Figure 9.Figure 9
Selection of the delay time by the autocorrelation function value of LTSA1.Then the embedding dimension is selected by the CAO method. The result is shown in Figure10; the optimal embedding dimension d for the projected vector is chosen as 10.Figure 10
Selection of the embedding dimension by the CAO method of LTSA1.Based on the selected optimal embedding dimension, the SVM model is used to achieve the multisteps prediction. The RBF kernel function is used:(30)K(xi,xj)=exp⁡(-12σ2∥x-xi∥2).In this research the regularity parameterC is set to 90.3, the kernel function parameter σ2 is set to 20, and the ε is set to 0.001. The parameters C and ε are selected by Particle Swarm Optimization algorithm (PSO) [30]. The popular size of the PSO is set to 100, the interaction number of the PSO method is set to 20, and the fitness function of the PSO method is set to choose the parameters which make the SVM model fitting error in the training process the smallest. The error goal is set to 0.05, the dimension of the PSO is set to 2, and the inertia weight is set to w=0.5, c1=c2=1.2. Based on the selected typical features, the features dataset is used to train SVM model and the input features number of SVM is 9 determined by embedding dimension. Then, the trained SVM model is used to predict the bearing running state. Before the 700th points, the bearing is working in a normal state, so the 701–900 points are used to train the SVM model and the following 85 points are employed for testing. In order to evaluate the predicting performance, the root-mean square error (RMSE) is utilized as follows:
(31)RMSE=∑i=1N(yi-y^i)2N,
where N represents the total number of data points in the test set, yi is actual value in training set or test set, and y^i represents the predicted value of the model. The actual value and the predicted result are shown in Figure 11.Figure 11
Prediction result based on the SVM model.From Figure11 we can see that the trend of the bearing running state can be predicted by SVM model. From the prediction, we can get a general understanding of the bearing running state in the future, but the predicted result is not accurate, especially the stage of 70–85 point, through calculation. The RMSE of the actual and the predicted result is 0.0469, so the predict result is not satisfied and the prediction error of SVM model (the predicted results subtract the actual data) is shown in Figure 12.Figure 12
The prediction error of the SVM model.Then the SOM method is used to divide the error into some districts, the iteration number if SOM is 2000; the structure of the state classification matrix is[3×1]. The results of the state division by SOM are shown in Table 2.Table 2
The result of the state division and the prediction results.
Prediction number
Prediction error
State
Prediction number
Prediction error
State
Prediction number
Prediction error
State
1
−0.001622
2
31
−0.014505
2
61
−0.009996
2
2
−0.003811
2
32
−0.017872
2
62
−0.015532
2
3
−0.003841
2
33
−0.022864
2
63
−0.041807
1
4
−0.007212
2
34
−0.031673
1
64
−0.022312
2
5
−0.008162
2
35
−0.033233
1
65
−0.034537
1
6
−0.02057
2
36
−0.05774
1
66
−0.060772
1
7
−0.023467
1
37
−0.090116
1
67
−0.10752
1
8
−0.026518
1
38
−0.04199
1
68
−0.096034
1
9
−0.02786
1
39
−0.034686
1
69
−0.12397
1
10
−0.027196
1
40
−0.032578
1
70
−0.13711
1
11
−0.038575
1
41
−0.010655
2
71
−0.13712
1
12
−0.025033
1
42
−0.013265
2
72
−0.03195
1
13
−0.018556
2
43
−0.006294
2
73
−0.02958
1
14
−0.019786
2
44
−0.007407
2
74
−0.19645
1
15
−0.01983
2
45
−0.021696
2
75
−0.18364
1
16
−0.020994
2
46
−0.017219
2
76
0.24355
3
17
−0.032989
1
47
−0.037728
1
77
0.43612
3
18
−0.019086
2
48
−0.011678
2
78
0.20092
3
19
−0.013994
2
49
−0.001478
2
79
−0.05844
1
20
−0.014597
2
50
−0.003989
2
80
0.76262
3
21
−0.013641
2
51
−0.00038
2
81
0.14341
3
22
−0.014209
2
52
−0.00578
2
82
0.28256
3
23
−0.013578
2
53
−0.025056
1
83
0.1326
3
24
−0.022126
2
54
−0.045813
1
84
0.3338
3
25
−0.019049
2
55
−0.01996
2
85
0.18377
3
26
−0.013825
2
56
−0.05286
1
27
−0.012682
2
57
−0.019165
2
28
−0.014012
2
58
−0.039591
1
29
−0.012089
2
59
−0.015265
2
30
−0.01567
2
60
−0.058697
1It can be seen from Table2 that when there is downward trend or the point value is less than 0, the state is set to 1 and when there is upward trend, the state is set to 3. According to the classification of the SOM model, the Markov state is divided into the following districts: [-0.10752,-0.00038], [-0.00038,0.1326], and [0.1326,0.43612]. The Markov model is used to improve the prediction error. For example, the 78 point, the Markov state transition matrix from 1–77 point is
(32)p=[0.6970.27270.03030.23260.76740001].Chose three points 77, 76, and 75 which are recent to the 78 point and set the transfer step as 1, 2, 3; the state prediction results based on the Markov model are shown in Table3.Table 3
Table of the probability status.
Point
Status probability
Transfer step
State 1
State 2
State 3
77
3
1
0
0
1
76
3
2
0
0
1
75
1
3
0.4757
0.4562
0.0681
Total
0.4757
0.4562
2.0681From Table3 we can see that the accumulated value of stage 3 is the largest, so the stage of 78 point is set to 3 and the result is the same with the SOM clustering method.In this research, according to the stage of the prediction of the Markov model, the correct value is calculated byx~(0)(k)=x^(0)(k)-ε, where ε is the median value of the divided stage area and x^(0)(k) is the value predicted though SVM model.For the 78 point, the corrected value is0.57298-(0.1326+0.43612)/2=0.2956, where the actual value at this point is 0.37206. Then other points have also been corrected though this method. The corrected results of the point 70–85 are shown in Table 4.Table 4
The corrected results of the points 70–85 based on Markov model.
Point
SVM model prediction error
Relative error (%)
SVM-Markov model prediction error
Relative error (%)
70
−0.13711
84.37%
−0.1025
63.07%
71
−0.13712
79.91%
−0.1235
71.97%
72
−0.03195
16.92%
−0.0121
6.41%
73
−0.02958
12.85%
0.0965
41.92%
74
−0.19645
84.44%
−0.09512
40.88%
75
−0.18364
84.21%
−0.10236
46.94%
76
0.24355
93.60%
0.01043
4.01%
77
0.43612
79.88%
0.03542
6.49%
78
0.20092
54.00%
−0.07646
20.55%
79
−0.05844
23.60%
−0.01574
6.36%
80
0.76262
344.98%
0.13262
59.99%
81
0.14341
26.80%
0.03562
6.66%
82
0.28256
114.19%
0.20302
82.05%
83
0.1326
49.59%
−0.00251
0.94%
84
0.3338
2061%
0.03291
203.20%
85
0.18377
1132.15%
0.06548
403.40%From Table4 we can see that the Markov model makes the results more precise, which validate the necessary to use the Markov model to improve the effect of the proposed method. The RMSE of the actual and the predicted result is 0.0091, so the prediction accuracy improved significantly. However, because those points are so far away compared to the predicted point of 1–69, the results still have some error.In order to compare the predict effect, the most usually used prediction model BP neural networks is used to predict the bearing running state based on the selected features. The learning rate of the neural network and its momentum coefficient are 0.01; the weights are initialized to uniformly distribute random values between −0.1 and 0.1; the iteration number is 2000; the training error is 0.001; the input number is 9; the hidden number is 15; the output node number is 1. The prediction results are shown in Figure13.Figure 13
Prediction result based on the BPNN model.From Figure13 we can see that the prediction results based on the BPNN model is not working effectively. There are some peaks while in the same position the actual status is not obvious. The RMSE of the predicted result is 0.0932, so the prediction results of the traditional BPNN model is not more effective than the SVM model. In addition, the prediction method based on the BPNN has the problem of prediction results which are unstable; when the same data is used to train and predict, the results are different and even the neural network may fall into the local optimum as shown in Figures 14 and 15.Figure 14
The training process curve of the BPNN method.Figure 15
The actual data and the result predicted by BPNN fall into the local optimum.With the data, the proposed method has also been compared with other methods that had been proposed in relative research.(1) The principal signal features extracted by PCA are utilized by HMM to predict the bearing running state [31]. (2) The time domain and frequency domain features have been directly used as the input of the prediction model, and the result has been predicted by Neural Network algorithm [1]. (3) The original features have been extracted by PCA as the input of the SVM prediction model [32]. (4) The proposed method in this research. The RMSE of the different methods predicted results is shown in Table 5.Table 5
The RMSE results of different prediction methods.
Method
The features extracted by PCA are utilized by HMM prediction model
The features directly used as the input of the Neural Network prediction model
The features extracted by PCA as the input of the SVM prediction model
The proposed method
RMSE
0.0698
0.1023
0.0521
0.0091From Table5, we can see that the RMSE of different prediction methods is very different. The prediction method that the original features have been directly used as the input of the NN model works the worst. This is the reason why the original features are still with high dimension and include superfluous information, which is not appropriate for state prediction; in addition, the NN prediction model has the drawbacks of slow convergence and difficulty in escaping from local minima. The prediction method based on the HMM model works not more effective than the method based on the SVM model; that is the reason why the HMM is not appropriate for long time forecast. The proposed method works the best; this is because the LTSA features extraction method can effectively extract the typical features and reduce the dimension and the SVM-Markov model can predict the state more precisely than the SVM and Markov model only. So through the comparison we can get that the proposed method is very effective in bearing running state prediction.


## 5.2. Application

After validating the effectiveness of the proposed method, the method has been used to the actual application. The test rig is shown in Figure16.Figure 16
Bearing test tig.The bearings are hosted on the shaft and the shaft is driven by AC motor. The rotation speed is kept at 1000 rpm and a radial load of 3 kg is added to the bearing. The data sampling rate is 25600 Hz and the data length is 102400 points collected on the date of 2011.11.25 as shown in Figure17. Every 2 hours, the vibration data are collected for one time. The collected data from 2011.11.25 to 2011.12.17 are analyzed after running for 1 year.Figure 17
Vibration signal.The time domain and time-frequency domain methods are used to deal with the collected vibration data as described in Section2, Table 1. Then the features are normalized through x⌢=(x-xmin⁡)/(xmax⁡-xmin⁡) and processed into the interval [0,1]. The LTSA is used to reduce the dimensionality of calculated features and the result is shown in Figure 18.Figure 18
The first main projected vector of test bearing based on LTSA method.From Figure18, we can see that the bearing running state has a fluctuation and upward trend. Especially at 150 points, there is a sudden change of trend, which reflects the bearing’s working status change at this moment.After extracting the typical features, the CAO method is used to determine the embedding dimension of the SVM model. The delay timeτ is set as 2 for the projected vector values though the autocorrelation function as shown in Figure 19.Figure 19
Selection of the delay time by the autocorrelation function value of LTSA1.The embedding dimension is selected by the CAO method. The result is shown in Figure20 and the optimal embedding dimension d for the projected vector is chosen as 10.Figure 20
Selection of the embedding dimension by the CAO method of LTSA1.Based on the selected optimal embedding dimension, the SVM model is used to achieve the prediction. The regularity parameterC is set as 909.5 and the ε is set to 0.01 selected by PSO method. Based on the selected typical features, the features dataset is used to train SVM model, the 1–130 points are used to train the SVM model, and the following 20 points are employed for testing. The prediction results are shown in Figure 21.Figure 21
Prediction result based on the SVM model.From Figure21 we can see that the trend of the bearing running state can be predicted by SVM model; from the prediction, we can get a general upward trend similar to the actual status. In addition, the sudden change of points 8 to 12 (near the 150 points in original signal as mentioned in Figure 17) is also showed out. However, the results are still not precise. In order to improve the prediction effect, the SOM method is used to divided the prediction error into some districts, the iteration number if SOM is 1000; the structure of the state classification matrix is [3×1]. The results of the state division by SOM are shown in Table 6.Table 6
The result of the state division and the prediction results.
Point
Prediction error
State
Point
Prediction error
State
1
0
3
11
0.053935
3
2
−0.063821
1
12
−0.0052358
1
3
0.002906
3
13
−0.058422
1
4
0.0074271
3
14
−0.094039
1
5
−0.021835
2
15
−0.032173
2
6
0.024138
3
16
−0.034682
2
7
0.0059543
3
17
−0.014101
2
8
0.016244
3
18
−0.13365
1
9
0.022316
3
19
−0.06811
1
10
0.016692
3
20
−0.070675
1Based on the classification of the SOM model, the Markov state is divided into the following districts[-0.13365,-0.034682], [-0.034682,0], and [0,0.053935]; the Markov model is used to improve the prediction error. The corrected results of the points 1–20 are shown in Table 7.Table 7
The corrected results of the points 1−20 based on Markov model.
Point
SVM model prediction error
Relative error (%)
SVM-Markov model prediction error
Relative error (%)
1
0
0
0
0
2
−0.0103821
579.21%
−0.006361
178.31%
3
0.002906
8.5739%
0.000603
1.7791%
4
0.0074271
20.311%
−0.001021
2.7922%
5
−0.021835
274.63%
−0.009362
117.75%
6
0.024138
30.429%
0.013254
16.708%
7
0.0059543
12.666%
0.002615
5.5628%
8
0.016244
22.459%
0.010027
13.863%
9
0.022316
33.565%
0.009565
14.386%
10
0.016692
32.417%
0.009473
18.397%
11
0.053935
73.825%
0.021563
29.515%
12
−0.0052358
34.947%
−0.002355
15.719%
13
−0.058422
107.06%
−0.019347
35.455%
14
−0.094039
122.79%
−0.014236
18.589%
15
−0.032173
263.1%
−0.009878
80.779%
16
−0.034682
159.3%
−0.008969
41.197%
17
−0.004101
355.3%
−0.001892
199.5%
18
−0.13365
117.18%
−0.050841
44.576%
19
−0.06811
120.92%
−0.028323
50.282%
20
−0.070675
135.04%
−0.014267
27.152%From Table7, we can see that the Markov model make the results more precise.Through the validation and actual application result we can see that the proposed method can predict the future status of the bearing, which is necessary for us to make some plan and do maintenance to reduce the risk of unnecessary accident.


## 6. Conclusions

(1)
The time domain and time-frequency domain methods are used to extract the original features from the mass vibration data, and in order to reduce the original features dimension and the superfluous information of the original features, the multifeatures fusion technique LTSA is used to fusion the original features and reduce the dimension.
(2)
Use the proposed SVM model to achieve bearing running state prediction. The proposed approach is validated by real-world vibration signals. The results show that the proposed methodology is of high accuracy, which is effective for the bearing running state prediction.
(3)
This research gives an example of combined approaches for the bearing running state prediction. Through analysis and validation we can get that the proposed method takes good use of the advantages of each part and achieve a high recognition accuracy and efficiency.
(4)
As the redundancy increases, the complexity of computation increases as well. This is one of the main shortcomings of the proposed method, which will be explored in the future.



---

*Source: 717465-2014-03-05.xml*