Multifractional Brownian Motion and Quantum-Behaved Partial Swarm Optimization for Bearing Degradation Forecasting

Gradual degradation of the bearing vibration signal is usually studied as a nonstationary stochastic time series. Roller bearings are working at high speed in a heavy load environment so that the combination of bearing faults gradually degraded during the rotation might lead to unpredicted catastrophic accidents. )e degradation process has the property of long-range dependence (LRD), so that the fractional Brownian motion (fBm) is taken into account for a prediction model. Because of the dramatic changes in the bearing degradation process, the Hurst exponent that describes the fBm will change during the degradation process. A priori Hurst value of the conventional fBm in the prediction is fixed, thus inducing a minor accuracy of the prediction. To avoid this problem, we propose an improved prediction method. Based on the following steps, at the initial data processing, a skip-over factor is selected as the characteristics parameter of the bearing degradation process. Amultifractional Brownianmotion (mfBm) replaces the fBm for the degradation modeling. We will show that also our mfBm has the same property of long-range dependence as the fBm.Moreover, a time-varying Hurst exponentH(t) is taken to replace the constantH in fBm. Finally, we apply the quantum-behaved partial swarm optimization (QPSO) to optimize H(t) for a finite interval. Some tests and corresponding experimental results will show that our model QPSO+mfBm have a much better performance on the prediction effect than fBm.


Introduction
ROLLER bearings represent the largest application of rotating mechanisms. Since they are usually working in high speed and heavy load environment, then bearing faults might gradually lead to some catastrophic accidents [1]. erefore, it is important to study the long-range time forecasting of the bearing life to optimize the service assistance [2]. ere already exist some models for the analysis of bearing faults. e grey model (GM) was applied to forecast bearing vibration signal [3]. However, under heavy load conditions, the accuracy of forecasting for the bearing vibration signal is usually influenced by many factors, such as deterioration [4], pitting, cracking, etc. Moreover, the GM model lacks a reflection on the randomness of the data; it implies that the GM model should be enhanced. Because of the lacking information, the GM approach only considers the size of the instant forecasting value. Also the hidden Markov model (HMM) has been applied to forecast the bearing fault [5,6]. But the disadvantage of the HMM is that a state of the HMM depends only on the state itself, and this dependence is time-independent. Besides, the HMM can only predict discrete (step) catastrophe failures and cannot forecast the trend of continuous fault changes. e error back-propagation (BP) algorithm is also widely used in forecasting bearing vibration, but the selection structure of the BP algorithm is not single, its training efficiency is not high, so that overcomplicated structure may cause the overfitting phenomenon. As a result, it is impossible to achieve accurate forecasting results [7]. Compared to the artificial neural network algorithm [8], support vector machine (SVM) is based on the structural risk minimization principle, which has a better performance [9]. However, the SVM needs to set the penalty parameter that makes the prediction model more complex. Moreover, the calculation of the remaining useful life for the bearing is crucial. For this reason, some methods were developed, such as the Gaussian process (GP) model [10] and the recursive maximum likelihood estimation (RMLE) [11]. When using the RMLE to solve the unknown noise in the system caused by nonstationary operating conditions, the noise value is relatively large. e shortcoming of the GP is that the prediction results are performed on a relatively limited data set. Such studies employed the regression-based adaptive predictive models [12], the exponential model [13], and the Extended Kalman Filtering [14].
Since the bearing failure is a slowly varying process, the degradation process has a long-range time memory. e past states can affect the future degradation states. is property is called long-range dependent (LRD). ere are two kinds of methods to identify whether the process is of LRD. e autocorrelation function (ACF) is nonintegrable [15], or the power spectrum density (PSD) of the process is divergent at the origin [16].
Fractional Brownian motion (fBm) is a nonstationary stochastic process with fractal properties of self-similarity and LRD [17] and is widely applied in many forecasting problems, like, e.g., image processing [18] and network traffic analysis [19]. FBm has also some applications in degradation modeling and remaining useful life prediction [20]. Zhang et al. [21] use fBm to describe the degradation of the blast furnace wall. Qin and Lin [22] combine fBm with Delft3D, WRF model, and GIS to predict the trajectory of harmful algal blooms. Xi et al. [23] apply fBm to predict the remaining useful life of lithium batteries. Song et al. [24] use fBm for short-term power load forecasting. Gupta et al. [25] establish the relationship between DCT and discrete fBm to provide a theoretical basis for applying DCT to fBm signal reconstruction. FBm can be considered as a generalized form of the Brownian Motion (BM). Consequently, the forecasting results of the fBm model are more flexible than those obtained from the BM model. FBm is characterized by one parameter called Hurst exponent. e expression range of Hurst exponent is 0 < H < 1. Except for H � 0.5, fBm is an LRD model. Under the condition of H � 0.5, fBm degenerates to the BM, which is a model with independent increment process.
However, the conventional fBm degradation model has some restrictions. In fact, in the process of bearing degradation, the degradation process at the early stage of the fault is relatively stable, so that there will be a huge oscillation in the later stage. e dramatic change for the bearing degradation process indicates that the Hurst exponent will change at different stages of degradation, thus implying that H usually has multiple values. e Hurst value of fBm is a constant and will not change again during the current prediction phase.
In order to solve this problem, we propose an improved model called multifractional Brownian motion (mfBm).
Compared with fBm, the fixed exponent H will be replaced by a time-varying exponent H(t). Moreover the quantumbehaved particle swarm optimization (QPSO) will be proposed to optimize the fundamental parameters [26]. e main idea of QPSO is to optimize the search strategy of the particle swarm optimization (PSO) to get the global optimum solution.
In order to extract the data from bearing degradation, the skip-over factor is proposed to measure the tendency of the gradual collapse [27]. In our model, we propose an integrated approach for the mfBm and QPSO for forecasting a bearing skip-over series. e computed H for the skip-over series will be used to deduce parameters in the mfBm model. e QPSO will be used to get the global optimal H(t) value in the mfBm model; the optimal can be used to make the accurate forecasting. Finally, we adopt mfBm and mfBm-QPSO models for the forecasting experiment. Monte Carlo method is used to show that the combination of mfBm and QPSO is superior to the fBm model of bearing degradation forecasting.
is paper is organized as follows. Section 2 describes the fBm model to forecast the bearing skip-over series. e optimization of H(t) by QPSO for mfBm model is given in Section 3. Section 4 introduces the computation of the skipover factor. Some examples and experimental tests are given in Section 5, together with a discussion on our model compared with previous ones. Conclusion is given in Section 6.

The Forecasting Model of fBm
Let the fBm model be denoted by B H (t).
us, B H (t) is defined as follows [28]: where Γ(.) is the gamma function, τ is the time lag. H is the Hurst exponent. Reported methods for estimating H are variance method, absolute value method, curve fitting method, rescaled range analysis (R/S), periodic graph method, and wavelet method [28][29][30][31]. Here R/S is used to estimate the Hurst exponent. e time series of length T is divided into k subinterval of length n, where T � kn. e average value of each subinterval composes a new series Y(i). e difference between the maximum and minimum values of the new series and the standard of deviation S 2 (n) are calculated. en the ratio between the two values is as follows: 2 Complexity where X n and S 2 (n) are the mean and variance of the subinterval, respectively. e value of R/S looks like cn H as n ⟶ ∞, where c is a constant independent of n. e increment process of the fBm called fractional Gaussian noise (fGn) is stationary. e ACF of the fGn is expressed as where ε > 0 is used by smoothing fBm so that the smoothed fBm is differentiable. σ 2 is the intensity of the fGn, and it is written by en the fBm can generate a specific stochastic time series as the prediction for next power load series. e stochastic differential equations (SDEs) [32,33] of the prediction model is written by where μ is the drift parameter, σ is the diffusion coefficient. e Hurst exponent is the most significant parameter in fBm; once it is obtained, the two parameters μ and σ of the fBm model can be calculated by maximum likelihood estimation (MLE) [34,35].
Discretize (5), the increments of the series are in the form where w(t) is unit white noise. e prediction model of the fBm is rewritten as In fBm, the Hurst exponent H is a constant, and it will not change again during the current prediction phase. However, as the degradation progresses, the Hurst exponent is time-varying; thus, the series generated by (6) cannot characterize the nonlinear trend. erefore, fBm is replaced by the mfBm. For mfBm, H(t) will change at different stages of degradation, and it is a time-varying parameter. To further improve the prediction accuracy, the QPSO is used to optimize the H(t) of the mfBm.

Quantum-Behaved Particle Swarm for the Hurst Exponent Optimization
e QPSO introduces quantum computing into the PSO, so that each particle in the space has the quantum behavior. Since each particle of the PSO moves in the search space with a velocity vector to find its local optimum and the global optimum, the disadvantage of the PSO is that it will be trapped in its optimal local solution and cannot reach the optimal global one. To avoid this problem, the particles in QPSO have a probability to appear anywhere in space; thus, QPSO has the ability to avoid falling into a locally optimal solution by searching a globally optimal solution. It has also fewer parameters to be adjusted. e optimization computation of H(t) by QPSO will be exhibited as follows.
Clerc and Kennedy [36] have proved that each particle must converge to the local attractor p i,j to provide the convergence for the PSO algorithm. ey defined it as Comparing with the convergence of the PSO, the QPSO defines one-dimensional delta potential well for each dimension at the point p i (t), and each particle has a quantum behavior, which can be described by the Schrödinger [37]. e distribution function of the particle position is given as where L i,j (t) is the standard deviation of the distribution. e position of each particle can be obtained by [38] where u ∼ U(0, 1), m(t) is defined as the average of a global mean best position: where N is the population size, and p i is the local best position of ith particle. e values of L i,j (t) and the position H i,j (t + 1) are evaluated as follows: where β < 1.781, which is the only parameter that can guarantee the convergence of the particle. Continuing with the previous example, we use the previous 30 point to forecast next 30 point.  Table 1. Initialize the current position, individual optical positions, and the population's global optimal position. en each particle is used to find the global optimal solution.
(2) Initialize the local optimal position P i,j (0) for the current position H i,j (0). e mean optimal position Complexity m j (t) is computed by (11). Select a suitable valueÎ 2 and evaluate the objective function value f(H i,j (t)). (3) e local attractor p i,j (t) is computed by (8). If rand(0, 1) > 0.5, then else (4) e optimal position p i,j (t) of each particle H i,j (t) is updated by P g,j � arg min P i (P i ). (5) e global optimal position is given by H gbest � arg min P g (P g ). (6) End of iteration.
During the process of updating the particle's position, it is necessary to have a comparison between the fitness value of the current individual position of each particle and the fitness value of the global optimal position. e loop iteration is used to find the global optimal solution. Once the iteration number or minimum error is reached in QPSO, the optimal approximation power curve for the H can be found.
is effective value of H parameter will be used by the mfBm model for the next data series forecasting.

Computing the Skip-Over Series
e dimensionless parameter skip-over factor is essentially an improvement of variance. Variance is an indicator related to energy. Energy-related indices will be disturbed by drastic operating load changes. e skip-over factor reduces the sensitivity to energy; thus, the gradual degradation trend can be reflected better. e skip-over factor is the average of the sum of squares of the difference between the data and the mean of the minimum and the measure of the difference between the source data and the expected value of the minimum value. Based on the whole-period sampling, the vibration waveforms of bearings often show skip-over characteristic, such as the failure of foundation loosening. With the faults deteriorating, the skip-over factor also changes. e essence of the skip-over factor is amplitude modulation of a waveform. e computation of skip-over factor is as follows: (1) Firstly, the bearing vibration signals are sampled and normalized: where k is a scale factor. (2) e vibration signals of the bearing after normalization are divided into x 11 , x 12 , . . . , x 1m ; . . . ; x n1 , x n2 , . . . , x nm . (3) Let m � 1024 [27], among set x 11 , x 12 , . . . , x 1m , extract minimum value x min � x 1,p , 1 < p < m; then get the average value for all m sets: (4) Calculated variance

Complexity
Define D x as the skip-over factor. Figure 1 shows the prediction flowchart of the QPSO + mfBm.

Experiment
e actual data was generated from accelerated life test of roller bearing Center for Intelligent Maintenance Systems (IMS), University of Cincinnati. Four Rexnord ZA-2115 double-row bearings were installed on the shaft as shown in Figure 2. e shaft is driven by an AC motor and coupled by rub belts. Accelerometers were installed on the bearing housing. e rotation speed of the shaft was 2000 rpm, the sampling rate was set at 20 kHz, and the data length was 20,480 points. Bearing 1 of vibration signals with an outer race failure is selected. Figure 3 shows the full-lifecycle tendency of bearing within 7 days of the skip-over series. e skip-over series has relative stability and sensitivity to initial fault. e initial fault appears about 7000 min. Here the ACF method is used to calculate the Hurst exponent of the series. Denote the measured ACF of the skip-over series by r(k), R(k) be the modeled ACF. For simplicity, we use the normalized ACF. Equation (7) is rewritten as

Complexity
Let M � E[(R − r) 2 ] be the mean square errors. As the M is minimized, the corresponding H value is the Hurst exponent. Figure 4 shows the fitting ACF curve, the calculated Hurst exponent is H � 0.849. e high fit of the measured curve to the model curve also indicates that the skip-over series is in accordance with the fBm.
e H of this skip-over series is 0.5 < H < 1. We focus on fBm, PSO + mfBm, and QPSO + mfBm models. e value of each 30 points is selected to calculate one H, the model parameters μ, and σ of the fBm change with H. en, we can deduce X t+1 ⟵X t by (6), the fBm can generate the approximation series. e parameter estimation results are shown in Table 2. In this prediction model, the estimation of parameters is affected by two aspects: the observed value and the time interval. e feature extraction process is to integrate multiple points of the original vibration signal into one point, that is, skip-over factor. After feature extraction, the values of skip-over factor are determined, and the time interval is changed. When the time interval is small, the calculated parameter values will increase. e larger the time interval, the smaller the corresponding parameter values. Here, the time interval is set as Δt � 0.1, resulting in a large calculated value. Figure 5 shows the prediction by three models. In the fBm prediction, the average relative error of fBm is 4.0%. In PSO + mfBm prediction, the average relative error of mfBm is 3.2%. In QPSO + mfBm prediction, the average relative error of the QPSO + mfBm can be reduced to 2.5%. Obviously, QPSO + mfBm has performed better on prediction accuracy.
In order to obtain the best prediction accuracy, different combinations of observation data length and prediction step    Table 3 and the average relative error is shown in Table 4.
Obviously, 30 points predict next 30 points has minimum average error in statistic meaning. Figure 7 shows prediction results by three models from 600th to 981th points. Figure 8 shows a box-scatter plot describing the relative error distribution of the forecasting skip-over series. e mean relative error of the QPSO + mfBm is lower than the other two models. Meanwhile, the relative error distribution of the QPSO + mfBm is more concentrated than the other two models; it means the QPSO + mfBm has lower relative error and stable prediction performance.

Conclusion
In this paper, a hybrid framework (mfBm and QPSO) for predicting a bearing skip-over series is proposed. We have used the ACF method to estimate the Hurst value of the entire degradation process. e fitting results show that the degradation process of the bearing has the LRD property; moreover, the fBm-based model gives the more reasonable description of the process. However, we have shown that the fBM-based model can be significantly improved by adopting the QPSO + mfBm. By using this improved method, we have shown that the predicted value is more consistent with the actual value, and the average relative error is significantly lower. In terms of relative error distribution, the QPSO + mfBm is more efficient. e experimental results have also shown that the improvement to the LRD prediction method for fBm is effective.

Data Availability
e data used to support the findings of this study are available at (https://ti.arc.nasa.gov/tech/dash/groups/pcoe/ prognostic-data-repository/#bearing).

Conflicts of Interest
e authors declare that they have no conflicts of interest.