Degradation Trend Prediction for Rotating Machinery Using Long-Range Dependence and Particle Filter Approach

Timely maintenance and accurate fault prediction of rotating machinery are essential for ensuring system availability, minimizing downtime, and contributing to sustainable production. This paper proposes a novel approach based on long-range dependence (LRD) and particle filter (PF) for degradation trend prediction of rotating machinery, taking the rolling bearing as an example. In this work, the degradation prediction is evaluated based on two health indicators time series; i.e., equivalent vibration severity (EVI) time series and kurtosis time series. Specifically, the degradation trend prediction issues here addressed have the following two distinctive features: (i) EVI time series with weak LRD property and (ii) kurtosis time series with sharp transition points (STPs) in the forecasted region. The core idea is that the parameters distribution of the LRD model can be updated recursively by the particle filter algorithm; i.e., the parameters degradation of the LRD model are restrained, and thus the prognostic results could be generated real-time, wherein the initial LRD model is designed randomly. The prediction results demonstrate that the significant improvements in prediction accuracy are obtained with the proposed method compared to some state-of-the-art approaches such as the autoregressive–moving-average (ARMA) model and the fractional order characteristic (FOC) model, etc.


Introduction
The prognostics and health management (PHM) in mechanical equipment and systems has become extremely important in many industrial applications including aero-engines, metallurgy machinery, wind turbine, etc.; owing to the fact that the PHM is able to minimize downtime and maintenance costs, maximize machinery utilization, and ensure system availability [1,2]. Prognostics is a key issue in PHM, which aims at utilizing real monitoring data to facilitate relevant degradation indicators and trends (DITs) that depict the current system health status and evaluate of remaining useful life (RUL) of a device. Therefore, reasonable prediction models with high accuracy are very crucial in the fields of industrial applications, because they should overcome different challenges and obstacles from the sampled data with external noise, variation trend, model perturbations and other underlying uncertainties, etc. [3,4].
Over the past few decades, various prognostic efforts have been proposed for mechanical condition maintenance. Current existing techniques can roughly be classified into two categories: (i) physics-based prognostics and (ii) data driven prognostics models. For a physics-based approach, the degradation model requires the knowledge of material properties, failure physical mechanisms (FPM), and the geometric structure of the mechanical equipment to estimate its failure degradation trend or remaining useful life. Generally, physics-based models with sufficient physical knowledge in a specific system may obtain an accurate prediction assessment, for example, fatigue crack growth (FCG) model [5] and fatigue spall progression life (FSPL) model [6], etc. However, it is difficult to obtain adequate modeling information regarding specific knowledge in most actual applications, especially the failure degradation process, which is too complicated. On the other hand, usually, once a physics-based model is established, the relevant model's parameters cannot be adjusted in real-time, which means that the actual degeneration in running states as reflected in the observed vibration signals will not be used. In contrast, data-driven prediction models such as the Bayesian method [7,8], machine learning techniques [9,10], support vector machine [11,12], etc., which utilize the measured historical data to establish the failure degradation model and avoid merging all kinds of physical knowledge of the system allow for an adjustment of the model parameters to capture the degradation trend of vibration signals. However, the common drawback of data-driven approaches is that they do not establish a theoretical mechanism between the physical degradation state (e.g., defect size) and the parameters obtained from the prediction models.
Today, fractional order characteristic (FOC) methods have aroused great interest for time series prediction among scholars. The FOC method emphasize long-range time correlation between the historical degradation and the predicted data and it has been demonstrated that this characteristic is difficult to restrain and control in practical applications [13,14]. In this regard, the PHM community has chosen FOC models such as long-range dependence (LRD) and fractional Brownian motion (FBM) models as the state-of-the-art approaches to address this scenario. For instance, in ref. [15], Li et al., utilized a fractionally autoregressive integrated moving average (f -ARIMA) model to investigate the degradation tendency of a rolling bearing based on the equivalent vibration intensity (EVI) health indicator, and also proved that the prediction accuracy outperforms the classical statistical time series and computational intelligence models. In ref. [16], Li et al., proposed an improved rescaled range (R/S) statistic and hybrid FBM model for the rolling bearing degradation process. In ref. [17], Song et al., developed a FBM degradation model to predict the bearing running tendency during a serious fault condition.
Although fractional order property can be observed in a wide range of mechanical operation systems and some good prediction results were obtained by the FOC models, there are two potential problems that the above mentioned references have not solved: (1) If the fractional order or LRD property of the health indicator time series (HITS) is too weak (e.g., the Hurst exponent approaching 0.5), especially the time series at the stationary operation stage, the f -ARIMA and FBM models might become unusable and the prediction accuracy reduced dramatically if the classical FOC model lacks the optimal solution. (2) Similar to the physics-based model, the model parameters of the classical LRD model such as f -ARIMA cannot be adjusted with the evolution of an actual degeneration trend of the HITS, and the prediction accuracies are also greatly decreased. (3) If some sharp transition points (STPs) are contained in the health indicator time series (HITS), especially the time series at the incipient/severe fault phase, the traditional f -ARIMA and FOC approaches treat all time series values equally, which ignore the fact that the STPs value should be preserved at a larger weight, thus limiting their effectiveness in practical application.
To overcome these issues and improve the prediction accuracy, this paper proposes a new approach for degradation prediction based on the long-range dependence (LRD) and particle filter (PF) algorithm, using EVI and kurtosis health indicators of a rolling bearing as the tested time series. The proposed method presents a prognostic tool by recursively updating the LRD model parameters with online measurements based on the particle filter (PF) algorithm. The initial parameters in the LRD model are described as probability distribution functions to incorporate the FOC property of rotating machinery degradation trend. Then, the PF is introduced to recursively estimate and update the model parameters, in which the small weighted particles will be removed in order to overcome the degeneracy phenomenon of particles, meanwhile, the updated parameters are fed into the established LRD model to predict the future health indicators step-by-step. The proposed technique is implemented to the accelerated life test of rolling bearing, which achieves a higher prognostic accuracy compared with some state-of-the-arts models such as the autoregressive-moving-average ARMA, f -ARIMA and FBM.
The remainder of this paper is organized as follows: Section 2 introduces the algorithm and theoretical derivation of the fractional Gaussian noise and LRD model. Section 3 describes the algorithm of a particle filter frame. In Section 4, the results and discussion of the proposed method are presented. Finally, conclusions and future research objectives are drawn in Section 5.

Fractional Gaussian Noise and LRD Model
Let X(t) be a random fractional Gaussian noise (FGN) series, then its auto-correlation function (ACF) denoted by r(τ; ε, H) is given by [18], where 0 < H < 1 is the Hurst exponent, and σ 2 = (Hπ) −1 Γ(1 − 2H) cos(Hπ) is the intensity of FGN [19,20]. Parameter ε > 0 is utilized by regularizing fractional Brownian motion (FBM) so that the regularized FBM is differentiable (see ref. [18]). Without loss of generality, let ε = 1, the ACF form of the FGN in the discrete case is expressed by, Generally, the random fractional Gaussian noise (FGN) series correspond to three subclasses of the time series, and the subclasses can be defined by the specific Hurst exponent H: (a) LRD process with 0.5 < H < 1; (b) short-range dependence process (SRD) with 0 < H < 0.5; and (c) anti-persistence process (the time series displays a random walk behavior) with H = 0.5. Thus, the LRD time series is special case of the FGN.
For a stochastic process x(t), its ACF is denoted by r xx (τ) = E[x(t)x(t + τ)], and τ denotes the time lag. When x(t) exhibits LRD characteristic, its ACF is non-integrable, that is From the perspective of asymptotic approximation, equivalently, two almost equivalent definitions with ACF and two-sided power spectral density (PSD) S x (τ) are as follow, where c 1 and c 2 are constants. Note that the parameter β is the exponent of the LRD. As a matter of fact, the parameter β can be computed by Hurst H, i.e., β = 2 − 2H or H = 1 − β/2. The f -ARIMA as a typical fractional order characteristic (FOC) model, which includes all three cases is the fractional ARIMA (p, d, q) process, i.e., f -ARIMA (p, d, q), the f -ARIMA (p, d, q) formulation is described by [21,22], where 0 < d < 0.5, ε t is a Gaussian white noise, B is the backshift operator, i.e., Bx t = x t−1 , Bx 2 t = x t−2 , · · · . The function Φ(B) and Θ(B) are polynomials of the backshift operator B, Since d denotes fractional value, the operator ∆ d = (1 − B) d is defined by binomial series, we have, where Γ is the Gamma function. It can be seen from the Equations (5) and (6) that d is a parameter to measure LRD feature, which can be calculated by the Hurst exponent; i.e., d = H − 0.5. Usually, the larger the Hurst exponent value, the stronger the long-range persistence. Besides, for the Hurst parameter H, several estimators can be used to estimate the value of the H, including the absolute value method (AVM), rescaled range (R/S) analysis method, Abry-Veitch, Whittle estimator, periodogram method, etc. The detailed formulas or estimation methodologies are presented in [23,24]. In this paper, we will employ the R/S method to investigate the LRD characteristic of a bearing degradation time series.

Particle Filter Frame
The PFF is a sequential series simulation-based algorithm based on a recursive Bayesian filter (RBF), which utilizes the Monte Carlo (MC) algorithm to constantly adjust the weight and position of the samples (called particles) and revise the original experience conditional distribution according to the adjusted particle information [25][26][27]. Actually, it uses the particles and its weights as a discrete random metric to approximate and accordingly update the related probability distribution (RPD).

Bayesian Filter Algorithm
To illustrate the PFF algorithm, the theoretical background of general Bayesian filter theory should be firstly introduced. For a non-linear stochastic system, we use the state equation and measurement equation to describe the dynamic state-space formulation, i.e., where x k and y k represent the state variable value and data measurement at k, respectively. Symbols n k−1 and u k denote the process noise representing uncertainty and measurement update noises, and f k is state transition function. Symbol h k is measurement function, which denotes the nonlinear relationship between the dynamic system states and the noisy measurements. Within a Bayesian framework, the recursive process (or posterior distribution estimation) is basically divided into two parts.
(1) Suppose the prior probability distribution p(x k−1 |y 1:k−1 ) is known, the number of samples is N from the posterior distribution according to Equation (7). Considering the first-order Markov process, the approximation of the posterior distribution can be given by, (2) The posterior distribution function estimation can be updated by the Bayes formula as follows, p(x k |y 1:k ) = p(y k |x k )p(x k |y 1:k−1 ) p(y k |y 1:k−1 ) where p(y k |x k ) is likelihood function at time k, and p(y k |y 1:k−1 ) is the normalizing constant which can be calculated as follow, p(x k |y 1:k−1 ) = p(x k |y 1:k−1 )p(y k |x 1:k−1 )dx k (10)

Particle Filter Frame Algorithm
The specific PFF algorithm is described as follows: (a) Initialization.
Sampling N particles x i 0 N i=1 from the prior probability distribution function, and the initialization weight is ω i 0 = 1/N accordingly. (b) Sequential importance sampling (SIS). (1) Resample independently N times from the above discrete distribution. (2) The prior weights are used to update the new weights, as follow, where p(y k x i k ) is the likelihood function.
(3) Normalize the importance weights, i.e., k is the normalized weight of the i-th particle at time k, and N the number of particles. (c) Resampling.
After calculating the normalized weights, resampling is used to avoid the particle degeneracy phenomenon of the PFF algorithm by means of removing the small weighted particles and preserving the promising particles with the larger weight.
The threshold of resampling is defined as effective sample then generating a new set of particles x i, * k N i=1 by resampling N-times from p(x k |y 1: , and set the weights of all the new particles to 1/N accordingly. Finally, for the sampled particles with the associated weights, the state vector ∧ x k can be estimated with Figure 1 illustrates the filter evolution procedures of particles in the recursive processes, where the dash line represents the optimal solution. Algorithms 2018, 11, x FOR PEER REVIEW 6 of 18

LRD Predicting Operator Driven by Particle Filter Frame
The flow chart of the degradation trend prognostics for the rolling bearings based on the proposed method is illustrated in Figure 2. The main procedures of the degradation trend prognostics can be described as follows: Step 1. Particle initialization.
(a) The f-ARIMA degradation models are established by the tested EVI time series or kurtosis time series. In this stage, the initial f-ARIMA model order and coefficients distribution are determined, i.e., (b) By drifting all of the coefficients set randomly, that is, is N, and their weights are initialized as 0 1 Step 2. For k = 1, 2, …, L (L is step length): (a) For k = k − 1, predict the time series output When the new output time series point yk arrives, update each particle's weight: (c) Normalize the particles weight, i.e., (d) Resampling is performed to remove the small weight particles and generate promising new particles with the larger weight, thus the new particles set is Step 3. The state vector k x ∧ can be estimated as

LRD Predicting Operator Driven by Particle Filter Frame
The flow chart of the degradation trend prognostics for the rolling bearings based on the proposed method is illustrated in Figure 2. The main procedures of the degradation trend prognostics can be described as follows: Step 1. Particle initialization. (a) The f -ARIMA degradation models are established by the tested EVI time series or kurtosis time series. In this stage, the initial f -ARIMA model order and coefficients distribution are determined, i.e., , where p and q are model orders, d 0 is initial fractional value, and φ(i) and θ(i) are model coefficients.
By drifting all of the coefficients set randomly, that is, When the new output time series point y k arrives, update each particle's weight: (c) Normalize the particles weight, i.e., Resampling is performed to remove the small weight particles and generate promising new particles with the larger weight, thus the new particles set is { ∧ x i k }.
Step 3. The state vector ∧ x k can be estimated as Step 4. Train the f -ARIMA model using the update state vector ∧ x k to get the prediction value of ∧ y k+1 .
Step 5. Algorithm stops when k is equal to step length L, otherwise, go to step (2).
Algorithms 2018, 11, x FOR PEER REVIEW 7 of 18 Step 5. Algorithm stops when k is equal to step length L, otherwise, go to step (2).

Experimental Setup
The experimental data were acquired by the Center for Intelligent Maintenance Systems (IMS) [28,29]. The experimental setup and its schematic diagram are displayed in Figure 3a,b, respectively. The vibration acceleration signals were collected every 10 min while input shaft speed was 2000 rpm (33.33 Hz) and the sampling frequency was 20 kHz. Two kinds of accelerated life test experiments were investigated: (1) the first test was carried out successively for 34 days (about 21,560 min), and (2) the second test was carried out successively for 8 days (about 9840 min). In this experiment, the severe wear failure in the inner race of bearing 3 and severe wear failure in the outer race of bearing 1 were respectively detected after dismantling. More details about the experimental description and experiment parameters can be found at the following website: http://www.imscenter.net/. The health indicators time series; i.e., the EVI degeneration curve of bearing 3, and the Kurtosis degeneration curve of bearing 1-are illustrated in Figure 4a,b, respectively. The detailed calculation process of health indicator EVI can be found in ref. [15]. The kurtosis waveform presented in Figure 4b can be divided into running status based on vibration amplitude; i.e., (1) normal operation, (2)

Experimental Setup
The experimental data were acquired by the Center for Intelligent Maintenance Systems (IMS) [28,29]. The experimental setup and its schematic diagram are displayed in Figure 3a,b, respectively. The vibration acceleration signals were collected every 10 min while input shaft speed was 2000 rpm (33.33 Hz) and the sampling frequency was 20 kHz. Two kinds of accelerated life test experiments were investigated: (1) the first test was carried out successively for 34 days (about 21,560 min), and (2) the second test was carried out successively for 8 days (about 9840 min). In this experiment, the severe wear failure in the inner race of bearing 3 and severe wear failure in the outer race of bearing 1 were respectively detected after dismantling. More details about the experimental description and experiment parameters can be found at the following website: http://www.imscenter.net/.
Algorithms 2018, 11, x FOR PEER REVIEW 7 of 18 Step 5. Algorithm stops when k is equal to step length L, otherwise, go to step (2).

Experimental Setup
The experimental data were acquired by the Center for Intelligent Maintenance Systems (IMS) [28,29]. The experimental setup and its schematic diagram are displayed in Figure 3a,b, respectively. The vibration acceleration signals were collected every 10 min while input shaft speed was 2000 rpm (33.33 Hz) and the sampling frequency was 20 kHz. Two kinds of accelerated life test experiments were investigated: (1) the first test was carried out successively for 34 days (about 21,560 min), and (2) the second test was carried out successively for 8 days (about 9840 min). In this experiment, the severe wear failure in the inner race of bearing 3 and severe wear failure in the outer race of bearing 1 were respectively detected after dismantling. More details about the experimental description and experiment parameters can be found at the following website: http://www.imscenter.net/. The health indicators time series; i.e., the EVI degeneration curve of bearing 3, and the Kurtosis degeneration curve of bearing 1-are illustrated in Figure 4a,b, respectively. The detailed calculation process of health indicator EVI can be found in ref. [15]. The kurtosis waveform presented in Figure 4b can be divided into running status based on vibration amplitude; i.e., (1) normal operation, (2) The health indicators time series; i.e., the EVI degeneration curve of bearing 3, and the Kurtosis degeneration curve of bearing 1-are illustrated in Figure 4a,b, respectively. The detailed calculation process of health indicator EVI can be found in ref. [15]. The kurtosis waveform presented in Figure 4b can be divided into running status based on vibration amplitude; i.e., (1) normal operation, (2)   To demonstrate the validity and superiority of the proposed approach, the EVI and kurtosis degeneration curves are presented as test cases. As a first test case, the EVI time series with weak LRD property has been measured. The second test case is a kurtosis time series with some STPs. Meanwhile, some state-of-the-art prediction methods including the ARMA and FOC methods (f-ARIMA model [15] and FBM model [16]) are also employed for comparison. Figure 5a shows the EVI curve of bearing 3, from the 801st-point to 1200th-point. The Hurst exponent is calculated via the R/S approach, as shown in Figure 5b. For the time series X(t), the R/S statistic can be expressed by,

Case 1: The EVI Time Series with Weak LRD Property
where S(n) is the standard deviation (SD); i.e.,  To demonstrate the validity and superiority of the proposed approach, the EVI and kurtosis degeneration curves are presented as test cases. As a first test case, the EVI time series with weak LRD property has been measured. The second test case is a kurtosis time series with some STPs. Meanwhile, some state-of-the-art prediction methods including the ARMA and FOC methods (f -ARIMA model [15] and FBM model [16]) are also employed for comparison. Figure 5a shows the EVI curve of bearing 3, from the 801st-point to 1200th-point. The Hurst exponent is calculated via the R/S approach, as shown in Figure 5b. For the time series X(t), the R/S statistic can be expressed by,

Case 1: The EVI Time Series with Weak LRD Property
where S(n) is the standard deviation (SD); i.e., S(n) = 1 , and R(n) is the range in each segment. The Hurst H can be generated from the fitted straight line of log(R(n)/S(n)) vs. log(S(n)). The resulting H = 0.5455 is close to 0.5, thus the fractional differencing parameter d = H − 0.5 = 0.0455 > 0 indicates that the EVI time series is to conform to the long-range dependence property. However, the LRD feature is very weak, that is to say, the performance in the future degradation trend has equal probability to be below and above the performance in the current degradation trend. Usually, the prediction accuracy generated by the FOC models may be decreased greatly in this case [13,14,16]. Further, the EVI prognostic result with the proposed method will be described and the results are also compared with the classical ARMA, f -ARIMA, and FBM models.
Algorithms 2018, 11, x FOR PEER REVIEW 9 of 18 0.5 = 0.0455 > 0 indicates that the EVI time series is to conform to the long-range dependence property. However, the LRD feature is very weak, that is to say, the performance in the future degradation trend has equal probability to be below and above the performance in the current degradation trend. Usually, the prediction accuracy generated by the FOC models may be decreased greatly in this case [13,14,16]. Further, the EVI prognostic result with the proposed method will be described and the results are also compared with the classical ARMA, f-ARIMA, and FBM models. Initially, the LRD model (e.g., f-ARIMA) is defined as ( )(1 ) ( ) , the particular modeling process can refer to [15,30,31], the AR coefficients ϕ(1), ϕ(2), ϕ(3), …, ϕ(p), MA coefficients θ(1), θ(2), θ(3), …, θ(q), wherein the parameters p and q can be calculated by the Akaike information criterion (AIC) [32,33]. The AR coefficients, MA coefficients, and the fractional differencing parameter d are modeled as probability distributions following a certain distribution (called particles). The initial distribution of model parameters can determine the performance of sequential importance sampling (SIS) of PFF. In the update stage, the unknown parameters combined with the state transition probability p(xk|xk−1) can be estimated recursively. Once the update measurements stop, the posterior distribution function p(xk+l|xk) can be calculated to predict a one-step ahead EVI time series based on the latest updated parameters distribution. Furthermore, multi-step (here the step length is L = 40) ahead prognostics can be obtained by performing the prediction process in Bayesian inference. In this experiment, the proposed model is applied to predict 40 points vibration intensity ahead; i.e., 400 min in advance. The EVI prediction results and errors generated by the proposed method and the ARMA and FOC methods are shown in Figure 6, using the 400 EVI points (from the 801st point to 1200th point) as the historical-based data. The results indicate that the ARMA, f-ARIMA and FBM models cannot generally track the EVI fluctuation trend well under the time series with weak LRD property. As shown in Figure 6d,e, the predicted EVI time series fit the actual time series better than the ARMA, f-ARIMA, and FBM models. It can be seen from Figure 6d,e that, compared with the ARMA and FOC methods, the EVI predictions by the proposed method are closer to the actual values, the prediction results having smaller error values. The specific parameters for quantitative assessment can be found in Table A1 in Appendix A.  Initially, the LRD model (e.g., f -ARIMA) is defined as Φ(B)(1 − B) d x t = Θ(B)ε t , the particular modeling process can refer to [15,30,31], the AR coefficients φ(1), φ(2), φ(3), . . . , φ(p), MA coefficients θ(1), θ(2), θ(3), . . . , θ(q), wherein the parameters p and q can be calculated by the Akaike information criterion (AIC) [32,33]. The AR coefficients, MA coefficients, and the fractional differencing parameter d are modeled as probability distributions following a certain distribution (called particles). The initial distribution of model parameters can determine the performance of sequential importance sampling (SIS) of PFF. In the update stage, the unknown parameters combined with the state transition probability p(x k |x k−1 ) can be estimated recursively. Once the update measurements stop, the posterior distribution function p(x k+l |x k ) can be calculated to predict a one-step ahead EVI time series based on the latest updated parameters distribution. Furthermore, multi-step (here the step length is L = 40) ahead prognostics can be obtained by performing the prediction process in Bayesian inference. In this experiment, the proposed model is applied to predict 40 points vibration intensity ahead; i.e., 400 min in advance.
The EVI prediction results and errors generated by the proposed method and the ARMA and FOC methods are shown in Figure 6, using the 400 EVI points (from the 801st point to 1200th point) as the historical-based data. The results indicate that the ARMA, f -ARIMA and FBM models cannot generally track the EVI fluctuation trend well under the time series with weak LRD property. As shown in Figure 6d,e, the predicted EVI time series fit the actual time series better than the ARMA, f -ARIMA, and FBM models. It can be seen from Figure 6d,e that, compared with the ARMA and FOC methods, the EVI predictions by the proposed method are closer to the actual values, the prediction results having smaller error values. The specific parameters for quantitative assessment can be found in Table A1 in Appendix A.    Figure 7a shows the kurtosis curve bearing 1 at incipient fault stage, from the 546th point to 945th point. The Hurst exponent is also computed via the R/S approach, as shown in Figure 7b. The resulting H = 0.9640, which is much larger than 0.5, and d = H − 0.5 = 0.4640 indicates that the kurtosis time series have a strong LRD property. However, it should be noted that some STPs are contained in this kurtosis time series (here is historical region, which means the region is formed by the points 546th to point 945th), which illustrates the uncertainties of the predicted data. in this kurtosis time series (here is historical region, which means the region is formed by the points 546th to point 945th), which illustrates the uncertainties of the predicted data. In this experiment, 400 kurtosis points (from the 546th point to 945th point) are used as the historical-based data to predict the 40 points; i.e., 400 mins in advance. The prognostic results and prognostic errors with the proposed method and other methods are displayed in Figure 8. It is observed that the prediction values by the f-ARIMA and FBM models gradually deviate from the actual values, especially in the last 15 points (some STPs are contained in the predicted region), which shows that the FOC models have failed when the STPs are contained in the predicted region. By comparison, the proposed approach performs better in multi-step prediction than the existent ARMA model and the FOC models. The specific parameters for quantitative assessment can be found in Table A2 in Appendix A. Overall, the above experimental results verify that the improvement of the proposed approach in terms of the weak-LRD property and STPs time series, which is able to provide a more accurate prediction for the rolling bearing degeneration. In this experiment, 400 kurtosis points (from the 546th point to 945th point) are used as the historical-based data to predict the 40 points; i.e., 400 mins in advance. The prognostic results and prognostic errors with the proposed method and other methods are displayed in Figure 8. It is observed that the prediction values by the f -ARIMA and FBM models gradually deviate from the actual values, especially in the last 15 points (some STPs are contained in the predicted region), which shows that the FOC models have failed when the STPs are contained in the predicted region. By comparison, the proposed approach performs better in multi-step prediction than the existent ARMA model and the FOC models. The specific parameters for quantitative assessment can be found in Table A2 in Appendix A. Overall, the above experimental results verify that the improvement of the proposed approach in terms of the weak-LRD property and STPs time series, which is able to provide a more accurate prediction for the rolling bearing degeneration.

Case 2: The Kurtosis Time Series with STPs
Algorithms 2018, 11, x FOR PEER REVIEW 11 of 18 in this kurtosis time series (here is historical region, which means the region is formed by the points 546th to point 945th), which illustrates the uncertainties of the predicted data. In this experiment, 400 kurtosis points (from the 546th point to 945th point) are used as the historical-based data to predict the 40 points; i.e., 400 mins in advance. The prognostic results and prognostic errors with the proposed method and other methods are displayed in Figure 8. It is observed that the prediction values by the f-ARIMA and FBM models gradually deviate from the actual values, especially in the last 15 points (some STPs are contained in the predicted region), which shows that the FOC models have failed when the STPs are contained in the predicted region. By comparison, the proposed approach performs better in multi-step prediction than the existent ARMA model and the FOC models. The specific parameters for quantitative assessment can be found in Table A2 in Appendix A. Overall, the above experimental results verify that the improvement of the proposed approach in terms of the weak-LRD property and STPs time series, which is able to provide a more accurate prediction for the rolling bearing degeneration.

Conclusions
This paper proposes a novel prognosis approach for the degradation trend of rotating machinery based on a long-range dependence and particle filter algorithm. The main objective is to address the problem of prognostics in FOC models while health indicators time series possess a weak LRD property and the STPs are contained in the predicted range. The two main contributions of this work can be summarized as: (1) For engineering applications, a new condition prognostic method is introduced as an effective tool for the long lifetime prediction and assessment of mechanical equipment in the PHM filed, the mechanical running degradation pattern can be predicted well in real-time and the most vital components of the mechanical equipment can be repaired and replaced prior to actual catastrophic failure. (2) For theoretical analysis, the particle filter frame (PFF) is embedded in the FOC model, an update scheme for online model parameters has been introduced to adapt the initial state model based on the PFF algorithm. Due to the adaptability of the parameters, the performance of our proposed model is remarkably efficient than classical time-series ARMA model and the FOC models, even the initial health status is unknown and the failure degradation behavior is timevarying. That is to say, the main advantage of the proposed approach is that the prognostic model has an ability to generate the reliable probabilistic results despite the uncertainties of the

Conclusions
This paper proposes a novel prognosis approach for the degradation trend of rotating machinery based on a long-range dependence and particle filter algorithm. The main objective is to address the problem of prognostics in FOC models while health indicators time series possess a weak LRD property and the STPs are contained in the predicted range. The two main contributions of this work can be summarized as: (1) For engineering applications, a new condition prognostic method is introduced as an effective tool for the long lifetime prediction and assessment of mechanical equipment in the PHM filed, the mechanical running degradation pattern can be predicted well in real-time and the most vital components of the mechanical equipment can be repaired and replaced prior to actual catastrophic failure. (2) For theoretical analysis, the particle filter frame (PFF) is embedded in the FOC model, an update scheme for online model parameters has been introduced to adapt the initial state model based on the PFF algorithm. Due to the adaptability of the parameters, the performance of our proposed model is remarkably efficient than classical time-series ARMA model and the FOC models, even the initial health status is unknown and the failure degradation behavior is time-varying.
That is to say, the main advantage of the proposed approach is that the prognostic model has an ability to generate the reliable probabilistic results despite the uncertainties of the initial model parameters.
As we all known, in the case of offline, the vibration data has been collected for further analysis (e.g., trend prediction or remaining useful life, i.e., RUL), and many adaptive prediction methods (such as Wavelet neural network (WNN) and ARMA-recursive least squares, i.e., ARMA-RLS, etc.) with parameters adjustment strategy can be performed only under the known data. Thus, the forecast errors are adjusted based on the predicted data and original data. If the original data is known before the prediction the adaptive prediction methods can work well, but in most cases, the original data is unknown, especially when the equipment is working under normal conditions, so the parameters adjustment strategy does not make sense. Based on this principle, if the original data is unknown, the proposed model is actually able to correctly identify when the bearing lifetime will enter the incipient fault region based on the data from the normal operational region. If not, the proposed model may lose its function.
Future research will investigate the suitability of the proposed model under harsh working conditions or different operational conditions, and the prediction for the remaining useful life (RUL) of the bearing will be investigated. As well, the proposed model might be expanded to other complex industrial applications such as diesel engine failure degradation and aero-engine fault monitor, etc.

Conflicts of Interest:
The authors declare no conflict of interest.