Vibration performance prediction and reliability analysis for rolling bearing

2.1. Forecasting theory

Suppose the rolling bearing vibration is given as $\mathbf{X}=\left\{e\right(i),i=\mathrm{1,2},\cdots ,T\}$, $T$ is the number of time series. The matrices $\mathbf{X}$ and $\mathbf{Y}$ can be obtained using the delay coordinate method and are shown as follows:

where $N=T-1-(m-1)\tau$, $\tau$ is delay time which is obtained by the autocorrelation function method, $m$ is embedding-dimension which is calculated by the Cao-method.

According to the Takens theory, the reconstructed phase trajectory is dynamically equivalent to the original system in the homeomorphic sense. The rolling bearing vibration signal $\mathbf{Y}\left(i\right)$ can be predicted based on the point $\mathbf{X}\left(i\right)$ in the phase space, and the mapping function $f$ is shown as follows, from which the next data point of time series can be gotten:

The own chaotic characteristics of the above prediction model, and the mapping function $f$ has a non-linear structure. Traditional statistical methods such as auto-regression, moving average, and ARIMA are not suitable for solving it. Machine learning model ELM can establish non-linear mapping, which can act as mapping function to do such a prediction work, the specific prediction principle is shown in Fig. 1.

Fig. 1Chaos prediction theory

2.2. Delay time calculation theory

The autocorrelation function method is used to calculate a delay time parameter, which is constructed based on the linear correlation degree of two motion trajectory parameters at time $t$ and $t+\tau$. The autocorrelation function $C\left(\tau \right)$ of vibration series $X=x{}_1{}^{},x_2,\cdots ,x_n$ is shown as follows:

where $\overline{x}$ is the mean value of the sample.

If the autocorrelation function shows significant attenuation trend with the growth of delay-time, the optimal delay time is equal to the time when the autocorrelation function reaches to $1-1/e$ times the initial value for the first time.

2.3. Cao method for embedding dimension

Cao Liangyue et al proposed the Cao method, which has the advantages of not-relying on the subjectivity and high calculation efficiency, let:

where, $‖•‖$ is phasor norm, $x_i(m+1)$ is the $i$th quantity of reconstructed phase space, where the embedding dimension is $m+1$, $x_{n(i,m)}$ is an integer greater than 1 and less than or equal to $N-m$, and the sequence with the shortest distance from $x_i(m+1)$ is $x_{n(i,m)}(m+1)$.

The mean value of all $a(i,m)$ is calculated using the above Eq. (5):

The embedding dimension is determined by observing the change of mean $E\left(m\right)$, along with $m$, the change ratio is shown as follows:

When the change is stable, and the value of $m$ plus 1 at this time is the required embedding dimension.

The criteria for determining the stable change of $E_1\left(m\right)$ are given below, the specific assessment process is as follows:

(1) ${\mathrm{\Delta }}_i$ is calculated as:

(2) An initial threshold $e$ is chosen according to the fluctuation of $E_1$, here, $e=\overline{{\mathrm{\Delta }}_i}$, $\overline{{\mathrm{\Delta }}_i}$ represents the average value of ${\mathrm{\Delta }}_i$, and the subscript $u$ of the first ${\mathrm{\Delta }}_i<e$ is found, the value of ${\mathrm{\Delta }}_j$ is calculated:

(3) The values of $e$, $i$ are reset as: $e=\overline{{\mathrm{\Delta }}_i}$, $j\le i\le N-1$, where $\overline{{\mathrm{\Delta }}_i}$ represents the average value of ${\mathrm{\Delta }}_i$

(4) $j\le i\le N-2$ is taken when ${\mathrm{\Delta }}_i>{\mathrm{\Delta }}_{i+1}$ and ${\mathrm{\Delta }}_{i+1}>{\mathrm{\Delta }}_{i+2}$ and ${\mathrm{\Delta }}_{i+1}<e$, embedding dimension is equal to $m=i+1$.

Through the above criteria, a relatively objective evaluation can be made when $E_1\left(m\right)$ tends to be stable, the embedding dimension calculation process is more scientific and rigorous.

2.4. ELM theory

ELM model is an algorithm proposed by Huang el al, which has its own advantages for training a single hidden layer feed-forward neural network. Its training model is shown in Fig. 2.

In Fig. 2, $n$ is the number of input-variables, $h$ shows the number of hidden layer neurons, $x_i$ represents the input variable, $i=\mathrm{1,2},3,\dots ,n$, $w_{ij}$ represent the connection weight between the input variable and hidden layer neuron $j$, $w_{jo}$ represent the connection weight between output variable and hidden layer neuron $j$, $b_j$ represent the threshold of hidden layer neuron, $j=\mathrm{1,2},3,\dots ,h$, $b_o$ represents the threshold of neuron in the output layer, $f_j(•)$ represents the activation function of layered neurons, $g(•)$ represents the activation function of output layer neurons.

Fig. 2ELM network structure

Its mathematical model is:

where, $v_j$ represents the output of output layer neurons $j$.

If there are $N$ valid samples, when the output threshold is equal to 0, and the output neurons activation function is a linear activation function, then the Eq. (10) can be written as:

where, $y=\left[y\right(1),y(2),\cdots ,y(N){]}^T$ is the network output vector, $w_o=[w_{1o},w_{2o},\cdots ,w_{ho}{]}^T$ represents the output weight vector, $v$ is the output matrix of hidden neurons, input weights and threshold matrix $w$ are randomly generated:

The expression of output weight vector is:

where, $v^{*}$ represents Moore-Penrose generalized inverse of output layer matrix $v\text{,}$$y_e=\left[y_e\right(1),y_e(2),\cdots ,y_e(N){]}^T$ is desired output. If $v\in R^{N\times h}$, $N\ge h$ and $rank\left(v\right)=h$, then the Moore-Penrose generalized inverse $v^{*}$ can be expressed as:

Putting Eq. (14) into Eq. (13), one can get:

The used hidden layer activation functions are Sigmoid function, Sin function, RBF function and Hardlim function.

2.5. Kriging model

The Kriging model can be expressed approximately as a sum of one random distribution function and one polynomial, as shown in Eq. (16):

where $y\left(x\right)$ is an unknown Kriging model, $y\left(x\right)$ is known as the two order regression function of $x$, the global approximation model in the design space is provided, $\beta$ is the undetermined coefficient of regression function, its value can be estimated by the known response value, $z\left(x\right)$ is a stochastic process, it is a local deviation on the basis of global simulation, and the expectation is 0, the variance is ${\sigma }^2$, the covariance matrix can be expressed as:

where R is the correlation matrix, $R\left(x^i,x^j\right)$ represents the correlation function of any two sample points $i$, $j=\mathrm{1,2},\dots ,n$, $n$ is the number of data in the sample. $R\left(x^i,x^j\right)$ has a variety of functional forms can be selected, common correlation functions include cubic function, gauss function, linear function, spherical function, spline function and so on.

2.6. Radial basis function

Radial basis function (RBF for short) model is a kind of function which takes the Euclidean distance as the independent variable between the test point and the sample point.

RBF is one of the commonly used surrogate model, the basic form is as follows:

Basic function $\phi ={\left(\varphi \left(‖x-x_1‖\right),\dots ,\varphi \left(‖x-x_{ns}‖\right)\right)}^T\text{,}$ weight coefficient ${\beta }_r={\left(({\beta }_r{)}_1,\dots ,({\beta }_r{)}_{ns}\right)}^T$, and ${\beta }_r$ should meet the interpolation conditions:

where $y_i$ is the true value, $(f_r{)}_i$ is the prediction value, $n_s$ is the number of samples.

With:

where $\varphi$ is radial function, common radial functions include Gaussian function, Multi-quadrics function, Reciprocal Multi-quadrics function, Thin-Plate spline function.

2.7. Accuracy evaluation method

Mean Squared Error (MSE), Root Mean Squared Error (RMSE) and Mean Absolute Error (MAE) indicators are used to calculate the prediction accuracy of the model, the expressions are:

where, $y_i$ represent the true value, $\widehat{y_i}$ represents the prediction value of the machine learning model, $N$ represents the number of test samples.

3.1. Grey-bootstrap method

According to the above prediction model, the vibration state information of bearing at the next $\mathbf{U}$ step can be predicted, with the vector $\mathbf{Y}$ expressed as follows:

where, $\mathbf{Y}$ is the chaotic prediction data using the above ELM model, $y\left(u\right)$ is the $u$th data of $\mathbf{Y}$ data series, $u=\mathrm{1,2},\dots ,U$.

Using the above method, a number is randomly selected from $\mathbf{Y}$ using equal probability, the total selection number is equal to $q$, so the bootstrap sample ${\mathbf{V}}_1$ can be gotten, the above selection method is repeated $\mathbf{B}$ times, and $\mathbf{B}$ samples are obtained, the bootstrap sample is shown as follows:

where, ${\mathbf{V}}_b$ represents the $b$th bootstrap-sample, $\mathbf{B}$ represents the total bootstrap-resample number, with:

where, $g=\mathrm{1,2},\dots ,q$ and $b=\mathrm{1,2},\dots ,B$.

Based on the grey prediction model GM(1,1), supposing that the first-order accumulated generating operator for ${\mathbf{V}}_b$ is as shown below:

The grey generated model can be described as the following differential equation:

where $\mathbf{u}$ is time variable, and $c_1$ and $c_2$ are the undetermined coefficients.

Use the increment instead of differential, the above formula can be expressed as:

where, $\mathrm{\Delta }\mathbf{u}$ is equal to the unit time interval. Furthermore, the mean generated sequence vector is set as below:

Under the initial condition $y_b\left(1\right)=v_b\left(1\right)$, the least-square-solution of the grey differential equation is shown as below:

where the determined coefficients $c_1$ and $c_2$ are shown:

with:

According to the inverse-accumulated generating, the $b$th generated data are expressed as follows:

Therefore, $\mathbf{B}$ generated data for rolling bearing vibration sample can be expressed as below:

where $w_b$ is the $b$th produced data.

3.2. Poisson counting process

3.3. Modeling basic ideas

The theoretical modeling uses a variety of mathematical models such as ELM chaotic-prediction-model, grey-bootstrap method and Poisson counting process. The modeling method is shown in Fig. 3. The concrete steps are as follows:

Step 1. Based on the rolling bearing time series $\mathbf{X}$ of the vibration performance, the embedding- dimension $m$ is obtained by the Cao-method, and the delay-time $\tau$ is obtained by the autocorrelation-method, so the phase-space-reconstruction can be realized based on the chaos theory.

Step 2. 30 step vibration data are predicted using the ELM chaotic prediction model, so the small sample $\mathbf{Y}$ can be constructed with the sample number of $\mathbf{U}$.

Step 3. A statistically large sample ${\mathbf{Y}}_{\mathbf{B}}$is generated from a small sample $\mathbf{Y}$ based on the grey-bootstrap method, so the Poisson count and bearing vibration variation strength can be easily gotten.

Step 4. Under the given vibration threshold $h$, the number $\mu$ is found out beyond the optimal vibration interval $\left[-h,h\right]$ from the large sample ${\mathbf{Y}}_{\mathbf{B}}$, and then the future $\mathbf{U}$-step variation intensity $\theta$ of the bearing is acquired, then, the reliability $R\left(\theta \right)$ of the bearing for each state in the future is obtained according to the Poisson counting formula.

Fig. 3Vibration prediction and reliability analysis

4.1. Vibration test data

The data used in this paper for experimental validation were provided by the Case Western Reserve University (CWRU), Cleveland, Ohio, USA [23]. As shown in Fig. 4, the test stand consists of an electronic motor, torque transducer, dynamometer, and control electronics. The bearing to be tested supports the motor shaft. The driving end bearing is SKF6205, and the fan end bearing is SKF6203. Acceleration sensors are placed above the bearing seat of the fan and the driving end of the motor to collect the vibration acceleration signal of the bearing.

Fig. 4Vibration experiment platform of rolling bearing

Fig. 5Layout of rolling bearing simulation experiment

Vibration signal is collected by a 16-channel data logger, the power and speed are measured by a torque sensor. The test data are obtained by an acceleration sensor on the bearing at the drive end, and the rotation speed is 1796 rpm. Fig. 5 depicts a rolling bearing simulation experiment layout.

The above experimental platform is used to collect the rolling bearing vibration signal. The dynamic time data signals are collected as shown in Fig. 6 to 10. Fig. 6 to Fig. 10 respectively depicts the vibration data under normal bearing and four inner ring fault size.

Fig. 6Vibration data (normal bearing)

Fig. 7Vibration data (0.18 mm inner ring fault size)

Fig. 8Vibration data (0.36 mm inner ring fault size)

It can be seen from Figs. 6-10 that with the deepening of bearing failure, the overall vibration data of rolling bearings show an increasing trend. The vibration data of normal bearings is generally distributed between [–0.2, 0.2], and only a few data exceed this range. The bearing vibration data of the inner ring fault size of 0.18 mm is generally between [–1, 1], and data graph shows a similar saw-tooth shape. The bearing vibration data of the inner ring fault size of 0.36 mm is generally at [–1, 1]. these data near 500 series are significantly larger than other data. The bearing vibration data of the inner ring fault size of 0.54 mm is generally between [–2, 2], and the graph looks like a series of oval-shaped components. The bearing vibration data of the inner ring fault size of 0.72 mm is generally between [–4, 4], and the graph looks like a random distribution.

Fig. 9Vibration data (0.54 mm inner ring fault size)

Fig. 10Vibration data (0.72 mm inner ring fault size)

4.2. Chaotic prediction with time series

The above mentioned bearing vibration time series are identified as ${\mathbf{X}}_1$ (Normal bearing vibration data), ${\mathbf{X}}_2$(0.18 mm inner ring fault size), ${\mathbf{X}}_3$(0.36 mm inner ring fault size), ${\mathbf{X}}_4$ (0.54 mm inner ring fault size), ${\mathbf{X}}_5$ (0.72 mm inner ring fault size). The chaotic prediction method is used to predict these sequences, and the prediction step is $\mathbf{U}=$30, then the 2001-2030 original data are used to verify the accuracy and feasibility of five sequence prediction models.

Phase space parameter calculation: The auto-correlation method and Cao-method respectively have been used to acquire the delay-time and embedding-dimensions, and their results are shown in Table 1.

In order to compare the prediction accuracy of the three methods mentioned above, ELM model, Kriging model and RBF model are used to predict the next 30-steps of bearing vibration data, the comparison result between 30-steps chaotic prediction value and experimental value of above three methods are shown in Fig. 11. Among them, the red line means the experimental data, blue line means the ELM model predicted value, pink line means the Kriging model predicted value, and green line means the RBF model predicted value. The RMSE values for prediction error of three models mentioned above are shown as Fig. 12.

Fig. 11Prediction results comparison of the three methods

Fig. 12Prediction accuracy evaluation of the three methods (RMSE)

From Fig. 11 and Fig. 12, the message can be acquired that ELM model has the highest prediction accuracy compared with Kriging model and RBF model, so the ELM model is used to predict vibration signal for reliability analysis.

The determination of the phase-space parameters of the rolling bearing time series is the basis for its phase-space reconstruction, and it is required for preparing the chaotic prediction model. The comparison results between 30-step chaotic prediction value and experimental value for sequences data ${\mathbf{X}}_1$, ${\mathbf{X}}_2$, ${\mathbf{X}}_3$, ${\mathbf{X}}_4$, ${\mathbf{X}}_5$ are shown in Figs. 13-17. Among them, the red line means the experimental data, and blue line means the predicted data.

From Fig. 13, the shape of the 30-step prediction result of the sequence ${\mathbf{X}}_1$ is similar to the curve of the experimental result, the points that differ greatly from the original data are shown in step 1 and steps 19-21, where the maximum difference is only 0.0568 V.

From Fig. 14, the prediction results of the first 5 steps of the sequence ${\mathbf{X}}_2$ are almost consistent with the original data. The steps 6, 13, 20, 22, 26, 28, 29 are significantly different from the original data, and the largest step is step 22 where the value is equal to 0.4152 V, the prediction results of other steps are almost consistent with the original data.

Fig. 13Vibration prediction data (X1)

Fig. 14Vibration prediction data (X2)

Fig. 15Vibration prediction data (X3)

Fig. 16Vibration prediction data (X4)

From Fig. 15, the prediction result of the sequence ${\mathbf{X}}_3$ is relatively different from the original data, the largest step is step 10 equal to 0.2904 V, and the prediction error between steps 2, 13, 19, 26, 28, 29 and original data is relatively small.

From Fig. 16, the prediction results of the sequence ${\mathbf{X}}_4$ are relatively similar to the original data. The fitting degree between prediction result line and the experimental result line at steps 1-12 and 19-30 are higher, only the prediction results of steps 12-18 are significantly different from the original data.

From Fig. 17, the prediction results of the sequence ${\mathbf{X}}_5$ are similar to the original data, and the prediction result lines of steps 6-15 are almost identical to the original data, the change trend of the prediction lines of steps 15 to 30 is very similar to the original data, with a jagged jump between –2 and 2.

Fig. 17Vibration prediction data (X5)

In order to scientifically evaluate the prediction accuracy, the RMSE value is used to calculate the prediction accuracy, as shown in Fig. 18.

Fig. 18Prediction accuracy evaluation (RMSE)

From Fig. 18, it is known that the RMSE evaluation value is small in each bearing state. When the inner ring fault size is 0.72 m, the RMSE value is approximately equal to 0.7 at the most.

In the 5 time series prediction of ELM chaotic prediction model, the difference between the predicted value and the experimental result is very small, and the two values maintain good consistency, that indicates that the prediction model is reliable.

4.3. Forecast results generated by grey-bootstrap method

Now, the prediction results of the next 30 steps of each sequence are processed by grey-bootstrap method to simulate the large-scale generated data of the vibration performance for each state of the bearing in the next 30 steps. The generated data is doing Poisson count under the given threshold value, and then find the performance reliability of the bearing in the next 30 steps. In the grey-bootstrap generation, setting the sampling number $q=$30 and repeated executions $\mathbf{B}=$10000, A large amount of generated-data ${\mathbf{Y}}_{\mathbf{B}}$ in the next 30 steps for sequence ${\mathbf{X}}_1$, ${\mathbf{X}}_2$, ${\mathbf{X}}_3$, ${\mathbf{X}}_4$, ${\mathbf{X}}_5$ is shown in Figs. 19-23.

Fig. 19Generated data of series X1 in future 30 steps

Fig. 20Generated data of series X2in future 30 steps

Fig. 21Generated data of series X3in future 30 steps

From the Figs. 19-23, the 10000 generated data simulated by grey-bootstrap-method has obvious randomness, and the corresponding value for sequences ${\mathbf{X}}_1$, ${\mathbf{X}}_2$, ${\mathbf{X}}_3$, ${\mathbf{X}}_4$ and ${\mathbf{X}}_5$ are concentrated in [–0.02, 0.06], [–0.2, 0.2], [–0.15, 0.3], [–0.4, 0.6] and [–0.5, 1]. The five sets of data show a non-linear increasing trend, which demonstrates that when the degree of failure deepens, the bearing vibration signal value gradually increases, and the inherent change law of the bearing vibration is closely related to the service life of bearing. Then, the evolution mechanism of performance reliability is determined and influenced.

Fig. 22Generated data of series X4in future 30 steps

Fig. 23Generated data of series X5 in future 30 steps

4.4. Future state reliability assessment

The bearing vibration threshold $h=$0.1 is set according to the generated data ${\mathbf{Y}}_{\mathbf{B}}$ in Figs. 19-23, and the number of generated data is to be found for each sequence exceeding the vibration threshold $h$, that is, the number $\mu$ that falls outside the optimal vibration interval [–0.1, 0.1] for the 10000 generated data of each sequence which can be calculated, so the variation intensity $\theta$ can also be obtained from the Eq. (39), then the reliability dynamic prediction result of the future 30 steps for each sequence can be calculated from Eq. (41). The results are shown in Table 2.

As it can be seen from Table 2, the data generated in the next 30 steps fall in the best vibration range for sequence ${\mathbf{X}}_1$, the variation intensity is 0, the reliability reaches 100 % that indicates that the bearing running state is relatively stable, no trace of bad behavior variation has occurred, and the state of maintaining the best vibration performance is very good.

In the sequence ${\mathbf{X}}_2$, the number of generated data in the next 30 steps that falls within the prescribed optimal vibration interval is 1768, the variation intensity is small and equal to 0.1768, and the reliability reaches 83.78 %, indicating that the bearing has a certain variation.

For the sequence ${\mathbf{X}}_3$, the number of generated data in the next 30 steps that falls within the prescribed optimal vibration interval is 3414, the variation intensity is relatively small and equal to 0.341, and the reliability is 71.11 %.

For the sequence ${\mathbf{X}}_4$, the number of generated data in the next 30 steps that falls within the specified optimal vibration interval is 5526, the variation intensity is more and equal to 0.553, and the reliability is 57.52 %.

For the sequence ${\mathbf{X}}_5$, the number of generated data in the next 30 steps that falls within the prescribed optimal vibration interval is 7432, the variation intensity is large and equal to 0.743, and the reliability is 47.56 %.

In order to describe the change of reliability intuitively, the change curve of the reliability parameter is shown in Fig. 24.

Fig. 24Change of reliability parameter R

As shown in Fig. 24, as the degree of bearing failure deepens, the reliability gradually decreases from 100 % with zero-failure to 47.56 % when the inner ring failure size is 0.7112 mm.

As shown above, in the process of chaos prediction of bearing vibration time series, the prediction method of chaos prediction model based on ELM is accurate and reliable that meets the general prediction requirements of engineering practice. Combined with the grey-bootstrap method, the prediction value of each state sequence data in the next 30 steps is sampled and processed, and a large number of generated signals participated in bearing performance degradation is accurately simulated. The variation intensity obtained by counting process can effectively describe the variation degree of bearing, and can reveal the influence mechanism on the bearing state on the variation process of performance reliability in operation.