Long Range Dependence Prognostics for Bearing Vibration Intensity Chaotic Time Series

Abstract: According to the chaotic features and typical fractional order characteristics of the bearing vibration intensity time series, a forecasting approach based on long range dependence (LRD) is proposed. In order to reveal the internal chaotic properties, vibration intensity time series are reconstructed based on chaos theory in phase-space, the delay time is computed with C-C method and the optimal embedding dimension and saturated correlation dimension are calculated via the Grassberger–Procaccia (G-P) method, respectively, so that the chaotic characteristics of vibration intensity time series can be jointly determined by the largest Lyapunov exponent and phase plane trajectory of vibration intensity time series, meanwhile, the largest Lyapunov exponent is calculated by the Wolf method and phase plane trajectory is illustrated using Duffing-Holmes Oscillator (DHO). The Hurst exponent and long range dependence prediction method are proposed to verify the typical fractional order features and improve the prediction accuracy of bearing vibration intensity time series, respectively. Experience shows that the vibration intensity time series have chaotic properties and the LRD prediction method is better than the other prediction methods (largest Lyapunov, auto regressive moving average (ARMA) and BP neural network (BPNN) model) in prediction accuracy and prediction performance, which provides a new approach for running tendency predictions for rotating machinery and provide some guidance value to the engineering practice.


Introduction
Rolling bearings are widely used as the most common and key components in rotating machinery systems.Generally, the working environment of rolling bearings is complex, with high speeds and heavy loads, and various failures such as wear, pitting, cracking and spalling, etc. always appear.Many unexpected failures deteriorate rapidly and may cause catastrophic accidents and financial losses or even personnel casualties [1][2][3].Therefore, condition prognostics and health management (CPHM) have become a particularly crucial issue and appealing field in mechanical fault diagnosis (MFD).
Vibration intensity time series of bearing systems belongs to a special kind of highly non-stationary and nonlinear data due to bearing systems that are restricted by the complicated operating mode and conditions, hence, modeling and prognostics approaches of time series is a challenge for researchers.Vibration intensity time series are composed of various components such as Gaussian noise, impulse, variation trend, model perturbations and other uncontrolled features, etc., and fit complex nonlinear time series well, thus it has good generalization ability for different LRD time series.Finally, the suitability of this prediction model of chaotic time series for obtaining more accurate multi-step forecasts will be demonstrated later.
The remainder of the paper is organized as follows: in Section 2, the chaotic characteristics of bearing vibration intensity time series are described with phase space reconstruction, LLE [22] and Duffing-Holmes Oscillator (DHO) [23].In Section 3, the LRD modeling, Hurst parameter estimation and the prediction steps of LRD are presented.In Section 4, the forecasting results and discussion of LRD and other different previous approaches are given in detail.Section 5 presents the conclusions and a summary.

Phase Space Reconstruction
A chaos attractor can estimate whether a times series is chaotic or stochastic, thus attractor reconstruction is a key step in chaotic time series analyses.Vibration intensity time series should be expanded into a high-dimensional space to reveal its dynamic characteristics, which was proposed by Takens [24].For a given time series txpiqu , i " 1, 2, ¨¨¨, N, selecting an optimal delay time τ and an optimal embedding dimension m, the phase space can be expressed as follows: Xpiq " rxpiq, xpi `τq, ¨¨¨, xpi `pm ´1qτqs where, i " 1, 2, ¨¨¨, M, M " N ´pm ´1qτ, is the number of elements in the reconstructed phase space, so it is clear that the phase space can be reconstructed for an optimal embedding dimension m and delay time τ.

The Optimal Delay Time τ
In order to reduce the reconfiguration calculation cost and keep the nonlinear characteristics of the time series, the delay time τ can be computed based on the C-C method [25], which was obtained by statistical results and has a strong ability of resist noise.
For embedding dimension m = 2, 3, 4, 5, radius r i " i ˆ0.5ε, t is index lag, i = 1, 2, 3, 4, ε is time series standard deviation, then calculate the following three formulas: ∆Sptq " S cor ptq " ∆Sptq `|Sptq| where, Sptq and ∆Sptq can reflect the correlation of bearing vibration intensity series, and the delay time τ will correspond to the first zero crossings point times of Sptq or the first local minimum value times of ∆Sptq, and the minimum of the quantity S cor ptq, respectively.

Optimal Embedding Dimension and Saturated Correlation Dimension
The saturated correlation dimension can measure the complexity of the reconstruction attractor in phase space.The saturated correlation dimension D and optimal embedding dimension m are calculated using the Grassberger-Procaccia (G-P) algorithm [26,27].This algorithm uses the phase space reconstruction of time series in Equation (1), which can be estimated by embedding dimension m through the principle that the correlation dimension increases with the increases of embedding dimension and reaching a saturated state.Given a small values of m 0 , the correlation integral Cpr, mq is given by: Hpr ´||Xpiq ´Xpjq||q (5) where, r is the radius of the sphere centered on Xpiq or Xpjq, ||Xpiq ´Xpjq|| is calculated as Euclidean norm, N is the length of time series, H(x) is the Heaviside function, defined as follows: There is a linear relationship between attractor dimension and cumulative distribution function, namely, Dpmq " | lnCpr, mq lnr |, thus the asymptotic curve was established by lnC(r,m)-lnr, then the correlation dimension D and embedding dimension m can be obtained by the formula m ě 2D `1.

Largest Lyapunov Exponent (LLE)
Lyapunov exponents generally measure the exponential rates of divergence or convergence of nearby trajectories in phase space, which are calculated to characterize chaotic time series.Given a time series, which is embedded in K-dimensional phase space, and has K Lyapunov exponents (λ 1 , λ 2 , ¨¨¨, λ k ), if the largest LLE is positive, λ ą 0, it means that the time series is chaotic.Similarly, the phase plane trajectory is a chaotic state, which can be obtained by the Duffing-Holmes Oscillator [23].In this paper, the LLE is calculated by using Wolf's algorithm [22], which is given as follows: given a initial point of reconstruct phase space X(t 0 ), locate a nearest neighbor point to the initial point X 0 (t 0 ), and denote the distance between those two points L 0 .At a later time t 1 , if the distance outside of the parameters specified ε, that is L 1 0 " |Xpt 1 q ´X0 pt 1 q| ą ε, then this procedure is repeated by selecting the new neighbor points at the replacement points on the attractor until the fiducial trajectory has traversed the entire data, the Lyapunov exponents λ is given by: where, L 1 i and L i are the Euclidean distances calculated between the two nearest neighboring points on the different trajectories, M is the total number of iteration steps.Therefore m and τ must be selected for the LLE mentioned above first.In addition, the maximum predictable steps (MPS) term is the reciprocal of the largest Lyapunov index, namely, t MPS " 1{λ.

Theory of Long Range Dependence
Long range dependence (LRD) of a process can be defined by an asymptotic power-law decrease of its autocorrelation function (ACF) and power spectral density (PSD) [28].Let X " px t : t " 1, 2, 3...q be a stochastic process with the ACF r xx pτq " Erxptqxpt `τqs, then X is called SRD series if r xx pτq is integrable [16,29], that is ş 8 0 r xx pτqdτ ă 8, on the other side, X is LRD if r xx pτq is nonintegrable, that is ş 8 0 r xx pτqdτ " 8.Moreover, the autocorrelation function follows asymptotically: where c 1 > 0 is a constant and 0 < β < 1. Equivalently, its two-sided PSD follows asymptotically: Entropy 2016, 18, 23 where, c 2 > 0 is a constant.The f -ARIMA (p, d, q) [30,31] process arises as a special case of LRD when the differencing parameter d in the range (´0.5,0.5).Let X " px t : t " 1, 2, 3...q represent the time series of vibration intensity.Then, LRD for the series can be formulated as: where ε t is white Gaussian noise (WGN), d " p´0.5, 0.5q, B the backshift operator, defined by BX t " X t´1 , BX 2 t " X t´2 , ¨¨¨, ΦpBq and ΘpBq are polynomial functions of the backshift operator B given by: ΘpBq " and ∆ " p1 ´Bq is fractional differencing operator defined by: where, , and Γ is the Gamma function.
Generally, we judge whether a stochastic time series is long-range dependent or not, from the perspective of the definition directly, the discrete integral of the time series is difficult and there might exists trend errors in the result; from another perspective, references [32,33] proposed that LRD exhibits self similarity that is characterized by the parameter H or Hurst parameter H P p0.5, 1q, thus the degree of self-similarity is expressed by the Hurst parameter.Clearly, fractional parameter d in f -ARIMA (p, d, q) is the indicator of LRD, d can be calculated by the Hurst exponent H.This parameter represents the speed of decay of a process' autocorrelation function.If 0.5 < H < 1, then the series is consistent with LRD and the closer it is to 1, the stronger the long memory or self-similarity; Otherwise H = 0.5 and 0 < H < 0.5 indicates that the time series is an anti-persistent process (or standard Brownian motion) [32] and SRD, respectively.The Hurst exponent H can be estimated by several methods, in this paper, R/S method is introduced to calculate the parameter H [33].

Hurst with R/S Method
Given a stochastic sample series X " px t : t " 1, 2, 3...q, with the average xpnq and the variance S 2 pnq, then R/S method can be expressed by: where, w k " px 1 `x2 `¨¨¨`x k q ´kxpnq, k " 1, 2, 3 ¨¨¨, n.
If the sample series satisfy LRD process, then: If the sample series consistent with SRD process, then: Ep Rpnq Spnq q " Cn 0.5 , n Ñ 8 For the formulas above, C > 0 is a constant, and the Hurst exponent H can be obtained by the least squares fitting algorithm for the R/S curve (Log(n)-logRS curve) in the logarithmic coordinates.

LRD Prediction Steps
The vibration intensity can reflect the overall running trend of the equipment, thus the bearing vibration intensity time series before the predicted point are selected as forecast samples, and the prediction steps of LRD algorithm are given below: Step 1: By means of abnormal data deletion and zero mean normalization for original time series samples in preprocessing operation in order to improve the prediction precision.
Step 2: The Hurst exponent of vibration intensity series was estimated by the R/S method.
Step 3: Then the differencing parameter d can be calculated by the formula H " d `0.5.Then the fractional differencing can be defined as follows: Generally, for fractional differencing, it is the key of the whole prediction process.
Step 4: Process Y " tY t , t " 1, 2, 3 ¨¨¨, Nu is considered as an ARMA model, and then the model parameters and orders can be determined by the Akaike information criterion (AIC), and the concrete operation process can be stated as follows: first, we determine the AR model parameters rφ 1, φ 2, ¨¨¨, φ p s with the Yule-Walker equation: where,rpkq is autocorrelation function (ACF).Then, process Y " tY t , t " 1, 2, 3 ¨¨¨, Nu is filtered by the Finite Impulse Response (FIR) filter ZpBq " 1 `p ř i"1 φ i B ´k, and treating process Z " tZ t , t " 1, 2, 3 ¨¨¨, Nu as a MA(q) model.Finally, the m orders AR model can be established with Z " tZ t , t " 1, 2, 3 ¨¨¨, Nu, where m ąą q, and using this model for linear prediction, that is equivalent to a AR(q), and the rθ 1, θ 2, ¨¨¨, θ q s can be generated by the Yule-Walker equation again.
Step 5: Using ARMA model to predict future time series.
Step 6: The prediction values of finally actual vibration intensity can be acquired by anti-fractional difference from Step 5, and the differential value is -d.The specific forecast process is as shown in Figure 1.

LRD Prediction Steps
The vibration intensity can reflect the overall running trend of the equipment, thus the bearing vibration intensity time series before the predicted point are selected as forecast samples, and the prediction steps of LRD algorithm are given below: Step 1: By means of abnormal data deletion and zero mean normalization for original time series samples in preprocessing operation in order to improve the prediction precision.
Step 2: The Hurst exponent of vibration intensity series was estimated by the R/S method.
Step 3: Then the differencing parameter d can be calculated by the formula 0.5 . Then the fractional differencing can be defined as follows: Generally, for fractional differencing, it is the key of the whole prediction process.
Step 4: is considered as an ARMA model, and then the model parameters and orders can be determined by the Akaike information criterion (AIC), and the concrete operation process can be stated as follows: first, we determine the AR model parameters with the Yule-Walker equation: where, ( ) r k is autocorrelation function (ACF).
Then, process , and treating process , where m q  , and using this model for linear prediction, that is equivalent to a AR(q), and the 1, 2, [ , ] q     can be generated by the Yule-Walker equation again.
Step 5: Using ARMA model to predict future time series.
Step 6: The prediction values of finally actual vibration intensity can be acquired by antifractional difference from Step 5, and the differential value is -d.The specific forecast process is as shown in Figure 1.In this study, vibration intensity was selected as forecast samples instead of the original acceleration signals or another indicator because the vibration fault acceleration signal has potential cyclical and non-stationary characteristics, whereas the vibration intensity time series is highly nonlinear, non-stationary and aperiodic.Therefore the vibration intensity time series is more suitable for the LRD method.The vibration intensity is defined as the root mean square (RMS) value of the vibration velocity ranging from 10 Hz to 1000 Hz, which is integrated and an effective characteristic to reflect the equipment operational condition.Vibration intensity value is the vibration level of the vibration velocity, and it includes the main parameters and features of the vibration signals.Moreover, it reflects the overall operating state of the mechanical equipment, hence the value of vibration intensity was selected as the sensitive factor for the bearing status [36].Given a discrete vibration velocity signal ( ) v n with N points, the vibration intensity can be expressed as: Next, vibration acceleration signals that are captured from each location and direction must be changed into velocity signals.Velocity signals can be obtained from accelerometer signals using numerical integration, however, the velocity signals always contain a trend of errors.Reference [37] proved that the errors trend can be eliminated by a polynomial fitting of the extreme value, thus more precise velocity signals can be acquired from this algorithm, and the vibration intensity can be calculated with Equation (17), thus the equivalent vibration intensity (EVI) is defined as: where, Vx, Vy, Vz are the vibration intensity direction of axis X, Y, Z, and the unit is mm/s; Nx, Ny, Nz are measuring points of axis X, Y, Z, respectively.Therefore, 2156 sets of vibration velocity points of bearing 3 are captured using Equation ( 18), Figure 3 shows the relationship between the whole accelerated life data and experimental data in this paper.According to equipment vibration international standard ISO-2372 (Class I), the vibration intensity across the fault critical line 7.1 in the In this study, vibration intensity was selected as forecast samples instead of the original acceleration signals or another indicator because the vibration fault acceleration signal has potential cyclical and non-stationary characteristics, whereas the vibration intensity time series is highly nonlinear, non-stationary and aperiodic.Therefore the vibration intensity time series is more suitable for the LRD method.The vibration intensity is defined as the root mean square (RMS) value of the vibration velocity ranging from 10 Hz to 1000 Hz, which is integrated and an effective characteristic to reflect the equipment operational condition.Vibration intensity value is the vibration level of the vibration velocity, and it includes the main parameters and features of the vibration signals.Moreover, it reflects the overall operating state of the mechanical equipment, hence the value of vibration intensity was selected as the sensitive factor for the bearing status [36].Given a discrete vibration velocity signal vpnq with N points, the vibration intensity can be expressed as: Next, vibration acceleration signals that are captured from each location and direction must be changed into velocity signals.Velocity signals can be obtained from accelerometer signals using numerical integration, however, the velocity signals always contain a trend of errors.Reference [37] proved that the errors trend can be eliminated by a polynomial fitting of the extreme value, thus more precise velocity signals can be acquired from this algorithm, and the vibration intensity can be calculated with Equation (17), thus the equivalent vibration intensity (EVI) is defined as: where, V x , V y , V z are the vibration intensity direction of axis X, Y, Z, and the unit is mm/s; N x , N y , N z are measuring points of axis X, Y, Z, respectively.Therefore, 2156 sets of vibration velocity points of bearing 3 are captured using Equation ( 18), Figure 3 shows the relationship between the whole accelerated life data and experimental data in this paper.According to equipment vibration international standard ISO-2372 (Class I), the vibration intensity across the fault critical line 7.1 in the last few points (see Figure 3a), and this sudden fault was caused by an artificially increased load in the laboratory.However, in a practical engineering application, incipient bearing faults occur based on the slowly changing characteristics over the whole life.Thus the vibration intensity time series in Figure 3b were selected to investigate the prediction efficiency of the LRD model, which could have some guidance value for engineering practice.
Entropy 2015, 17, 1-15 on the slowly changing characteristics over the whole life.Thus the vibration intensity time series in Figure 3b were selected to investigate the prediction efficiency of the LRD model, which could have some guidance value for engineering practice.

Chaotic Characteristic Analysis of Bearing Vibration Intensity Time Series
In this section, the chaotic characteristic can be evaluated based on the proposed 200 sets of bearing vibration intensity time series.

Chaotic Characteristic Analysis of Bearing Vibration Intensity Time Series
In this section, the chaotic characteristic can be evaluated based on the proposed 200 sets of bearing vibration intensity time series.Sptq and ∆Sptq curves of vibration intensity time series with C-C method as shown in Figure 4a,b.The result of Figure 4a pointed out that the first zero crossings point times of Sptq is t = 2, similarily, the first local minimum value times of ∆Sptq is t = 2 in Figure 4b, therefore the optimal delay time of vibration intensity time series is τ = 2.The curves between lnC(r,m) and lnr according to the G-P method are shown in Figure 5, and Figure 6 presents the relation between the embedding dimension m and the saturated correlation dimension D. Figures 5 and 6 prove that the curve region reaches a saturated state gradually as the embedding dimension m increases.Figure 6 shows that the convergence value of the correlation dimension is 6.5336, hence the embedding dimension of the vibration intensity time series is m = 14 The curves between lnC(r,m) and lnr according to the G-P method are shown in Figures 5 and 6 presents the relation between the embedding dimension m and the saturated correlation dimension D. Figures 5 and 6 prove that the curve region reaches a saturated state gradually as the embedding dimension m increases.Figure 6 shows that the convergence value of the correlation dimension is 6.5336, hence the embedding dimension of the vibration intensity time series is m = 14 according to the formula m ě 2D `1, which means that vibration intensity time series can be completely revealed in 14-dimensional space.
Entropy 2015, 17, 1-15 according to the formula  , which means that vibration intensity time series can be completely revealed in 14-dimensional space.Further, the largest Lyapunov exponent of the vibration intensity time series is calculated using the Wolf method, and the result is 0.0316

 
, that is 0   , which means that the vibration intensity time series is chaotic.In addition, the maximum predictable steps 1 31.64 , therefore, the maximum prediction time is 31 (point number or steps), and the prediction step is creditable under the condition of 31 t  .Figure 7 is the phase plane trajectory of vibration intensity time series with the Duffing-Holmes Oscillator (DHO), which also indicates that the vibration intensity time series belongs to a chaotic sequence. , which means that vibration intensity time series can be completely revealed in 14-dimensional space.Further, the largest Lyapunov exponent of the vibration intensity time series is calculated using the Wolf method, and the result is 0.0316

 
, that is 0   , which means that the vibration intensity time series is chaotic.In addition, the maximum predictable steps 1 31.64 , therefore, the maximum prediction time is 31 (point number or steps), and the prediction step is creditable under the condition of 31 t  .Figure 7 is the phase plane trajectory of vibration intensity time series with the Duffing-Holmes Oscillator (DHO), which also indicates that the vibration intensity time series Further, the largest Lyapunov exponent of the vibration intensity time series is calculated using the Wolf method, and the result is λ " 0.0316, that is λ ą 0, which means that the vibration intensity time series is chaotic.In addition, the maximum predictable steps t mps " 1{λ " 31.64,therefore, the maximum prediction time is 31 (point number or steps), and the prediction step is creditable under the condition of t ď 31.Figure 7 is the phase plane trajectory of vibration intensity time series with the Duffing-Holmes Oscillator (DHO), which also indicates that the vibration intensity time series belongs to a chaotic sequence.Further, the largest Lyapunov exponent of the vibration intensity time series is calculated using the Wolf method, and the result is 0.0316   , that is 0   , which means that the vibration intensity time series is chaotic.In addition, the maximum predictable steps , therefore, the maximum prediction time is 31 (point number or steps), and the prediction step is creditable under the condition of 31 t  .Figure 7 is the phase plane trajectory of vibration intensity time series with the Duffing-Holmes Oscillator (DHO), which also indicates that the vibration intensity time series belongs to a chaotic sequence.

Prediction Analysis of LRD
The Hurst exponent H is obtained by the R/S method as shown in Figure 8.The resulting H = 0.5933 > 0.5 shows that the vibration intensity chaotic time series is consistent with the LRD property, equivalently, the ACF r xx pτq is nonintegrable, thus the LRD model can be applied to chaotic time series to predict the future running vibration intensity points.Considering the prediction steps obtained above is 31, thus the LRD is applied to predict 30 points vibration intensity ahead, that is 300 min.

Prediction Analysis of LRD
The Hurst exponent H is obtained by the R/S method as shown in Figure 8.The resulting H = 0.5933 > 0.5 shows that the vibration intensity chaotic time series is consistent with the LRD property, equivalently, the ACF ( ) xx r  is nonintegrable, thus the LRD model can be applied to chaotic time series to predict the future running vibration intensity points.Considering the prediction steps obtained above is 31, thus the LRD is applied to predict 30 points vibration intensity ahead, that is 300 min.The predicted values of the chaotic time series generated by the LRD model, largest Lyapunov model (L-Lyap), ARMA model and BP neural network (BPNN) are illustrated in Figure 9. From Figure 9, it can be observed that the models are able to reasonably track the vibration intensity variation trend.The quantitative prediction results, prediction absolute error (AE) and relative error (RE) of four models are shown in Appendix.The forecasting errors that differ between actual and predicted values with the proposed models are shown in Figure 10.It is clear that the predicted results using the LRD model have better a match with the actual vibration intensity than the predicted result using other prediction methods, and the LRD method can offer more accurate forecasting result than other methods.From Figure 9, it can be observed that the models are able to reasonably track the vibration intensity variation trend.The quantitative prediction results, prediction absolute error (AE) and relative error (RE) of four models are shown in Appendix.The forecasting errors that differ between actual and predicted values with the proposed models are shown in Figure 10.It is clear that the predicted results using the LRD model have better a match with the actual vibration intensity than the predicted result using other prediction methods, and the LRD method can offer more accurate forecasting result than other methods.The predicted values of the chaotic time series generated by the LRD model, largest Lyapunov model (L-Lyap), ARMA model and BP neural network (BPNN) are illustrated in Figure 9. From Figure 9, it can be observed that the models are able to reasonably track the vibration intensity variation trend.The quantitative prediction results, prediction absolute error (AE) and relative error (RE) of four models are shown in Appendix.The forecasting errors that differ between actual and predicted values with the proposed models are shown in Figure 10.It is clear that the predicted results using the LRD model have better a match with the actual vibration intensity than the predicted result using other prediction methods, and the LRD method can offer more accurate forecasting result than other methods.Furthermore, in order to determine quantitatively the better model, the following error criteria are employed for a comprehensive model comparison as shown in Table 1.

Predict Error Criterions Computational Formula
Mean Absolute Error The four models discussed in the paper are applied to vibration intensity chaotic time series, and the performance comparison results are shown in Table 2.It can be observed in all the cases that Mean Absolute Error (MAE), Root-Mean-Square Error (RMSE), Normalized Mean Square Error (NMSE), Maximum of Absolute Error (Max-AE) obtained with the LRD model are considerably lower than those calculated with the other models, which clear shows that the LRD model for chaotic time series has high accuracy, therefore theLRD prediction model provides a new way for predicting the running tendency or remaining useful lifetime of rotating machinery in practical engineering.Furthermore, in order to determine quantitatively the better model, the following error criteria are employed for a comprehensive model comparison as shown in Table 1.The four models discussed in the paper are applied to vibration intensity chaotic time series, and the performance comparison results are shown in Table 2.It can be observed in all the cases that Mean Absolute Error (MAE), Root-Mean-Square Error (RMSE), Normalized Mean Square Error (NMSE), Maximum of Absolute Error (Max-AE) obtained with the LRD model are considerably lower than those calculated with the other models, which clear shows that the LRD model for chaotic time series has high accuracy, therefore theLRD prediction model provides a new way for predicting the running tendency or remaining useful lifetime of rotating machinery in practical engineering.

Conclusions
In this paper, a chaotic time series prediction approach for bearing vibration intensity based on LRD method is proposed.Results obtained from this research are as follows: 1.
The vibration intensity time series have chaotic properties, and the chaotic characteristics of vibration intensity time series can be jointly determined by the largest Lyapunov exponent and phase plane trajectory.

2.
Experimental results and analysis demonstrate the unique ability of the f -ARIMA model for chaotic time series.The long range dependence prediction method is powerful for problems characterized by chaotic time series, small samples and nonlinearity.The LRD prediction method is better than other prediction methods in prediction accuracy and prediction performance, thus providing a new approach for the running tendency prediction of rotating machinery and some guidance value for engineering practice.In addition, the Hurst exponent can be calculated by the time series self-feature, instead of a subjective integer difference.

3.
High prediction accuracy is the key and critical to CPHM for improving reliability and reducing overall maintenance costs.The significant economic impact of prediction accuracy or prediction errors in rotating machinery status prediction on the mechanical equipment market and the systematic chaotic and self-similarity properties of different rotating machinery prediction models would be meaningful topics for further research.

Figure 1 .
Figure 1.Forecasting process of the long range dependence.

Figure 1 .
Figure 1.Forecasting process of the long range dependence.
The vibration signal was generated by the NSF I/UCR Center for Intelligent Maintenance Systems (Cincinnati, OH, USA; IMS-www.imscenter.net) and the University of Cincinnati[34,35].Four Rexnord ZA-2115 double row bearings were installed on the shaft as shown in Figure2.The shaft is driven by an AC motor and coupled by rub belts.Accelerometers were installed on bearing housing.Four thermocouples were attached to the outer race of each bearing to record bearing temperature for the purpose of monitoring the lubrication.Vibration data was collected every 10 min by a NI 6062E DAQ Card.The sampling rate set at 20 kHz and the data length is 20,480 points.Inner race fault vibration signal of bearing 3 was produced by a National Instruments LabVIEW program from 22 October 2003 12:06:24 to 25 November 2003 23:39:56 (from normal to severe fault), that is to say, the overall test data had 2156 sets.Entropy 2015, 17, 1-15 Four Rexnord ZA-2115 double row bearings were installed on the shaft as shown in Figure 2. The shaft is driven by an AC motor and coupled by rub belts.Accelerometers were installed on bearing housing.Four thermocouples were attached to the outer race of each bearing to record bearing temperature for the purpose of monitoring the lubrication.Vibration data was collected every 10 min by a NI 6062E DAQ Card.The sampling rate set at 20 kHz and the data length is 20,480 points.Inner race fault vibration signal of bearing 3 was produced by a National Instruments LabVIEW program from 22 October 2003 12:06:24 to 25 November 2003 23:39:56 (from normal to severe fault), that is to say, the overall test data had 2156 sets.

Figure 3 .
Figure 3.The relationship between the whole accelerated life data and experiment data (vibration intensity time series).

Figure 3 .
Figure 3.The relationship between the whole accelerated life data and experiment data (vibration intensity time series).

8
In this section, the chaotic characteristic can be evaluated based on the proposed 200 sets of bearing vibration intensity time series.( ) S t and ( ) S t  curves of vibration intensity time series with C-C method as shown in Figure4a,b.The result of Figure4apointed out that the first zero crossings point times of ( ) S t is t = 2, similarily, the first local minimum value times of ( ) S t  is t = 2 in Figure4b, therefore the optimal delay time of vibration intensity time series is τ = 2.

Figure 4
Figure 4. ( ) S t and ( ) S t  curves for time delay.(a) ( ) S t curve for time delay; (b) ( ) S t  curve

Figure 6 .
Figure 6.Saturated correlation dimension method for the embedding dimension.

Figure 6 .
Figure 6.Saturated correlation dimension method for the embedding dimension.

Figure 6 .
Figure 6.Saturated correlation dimension method for the embedding dimension.

Figure 6 .
Figure 6.Saturated correlation dimension method for the embedding dimension.

Figure 7 .
Figure 7. Phase plane trajectory with Duffing-Holmes Oscillator (x is the signals produced by the D-H system, x' is the derivative of x).

Figure 7 .
Figure 7. Phase plane trajectory with Duffing-Holmes Oscillator (x is the signals produced by the D-H system, x' is the derivative of x).

Figure 9 .
Figure 9.Comparison of the forecasting for chaotic time series using the proposed models.

Figure 10 .
Figure 10.The forecasting errors distribution of vibration intensity.

Table 2 .
Performance comparison for vibration intensity chaotic time series.