A Prognostic Framework for Wheel Treads Integrating Parameter Correlation and Multiple Uncertainties

As crucial rotary components of high-speed trains, wheel treads in realistic operation environment usually suffer severe cyclic shocks, which damage the health status and ultimately cause safety risks. Timely and precise health prognosis based on vibration signals is an effective technology to mitigate such risks. In this work, a new parameter-related Wiener process model is proposed to capture multiple uncertainties existed in on-site prognosis of wheel treads. The proposed model establishes a quantitative relationship between degradation rate and variations, and integrates uncertainties via heterogeneity analysis of both criterions. A maximum-likelihood-based method is presented to initialize the unknown model parameters, followed by a recursive update algorithm with fully utilization of historical lifetime information. An investigation of real-world wheel tread signals demonstrates the superiority of the proposed model in accuracy improvement.


Introduction
As a fast and convenient transportation asset, high-speed trains have drawn world-wide attention in recent years. Under high-speed operating conditions, health monitoring for rotary machinery systems is a crucial technology to guarantee the safety and reliability of high speed trains [1]. With the rapid development of sensor technologies, multiple sources of operational/health data can be collected, supporting online diagnosis and prognosis [2][3][4]. For high-speed trains, engineers collect vibration signals of rotary machinery systems such as bearings, gears, axles and wheel treads via on-board monitoring device and then conduct fault recognition, location and isolation based on signal processing technologies. However, diagnosis only provides failure information and provides no guidance for preventive maintenance, which is quite crucial to mitigate malfunction risks. Therefore, it is of urgent necessary to develop advanced prognostic methodologies, which predicts real-time health information to support condition-based maintenance plan.
Remaining useful lifetime (RUL) is a main objective of health prognostic, which measures the distance between the latest monitoring time and system malfunction. Currently, RUL prediction methodologies can be substantially partitioned into two types, data-driven approaches and model-based approaches [5][6][7]. Data-driven approaches immediately utilize historical sensor data as input and outputs lifetime information. Representative data-driven technology is machine learning, including neural network [8][9][10][11], support vector machine [12,13] and artificial intelligence approaches [14][15][16]. Such technologies can adopt diverse application scenarios, but heavily rely on the quality and size of data. Model-based approaches, on the other hand, utilize statistical or physical modes to characterize the degradation paths. Combined with the measured data, model parameters can be identified and (1) A novel parameter-related Wiener process model with consideration of multiple uncertainties is proposed for accurately characterizing the real degradation process of wheel treads. (2) A recursive algorithm, which integrates KF and EM algorithm, is established to update model parameters and RUL distribution with the updating of online monitoring data. (3) An investigation of real-world wheel tread signals is used to demonstrate the superiority of the proposed prognostic framework in accuracy improvement. (4) The rest of the paper is organized as follows. Section 2 introduces the proposed RUL prediction methodology. Section 3 conducts experiment using real monitoring data of wheel treads of high-speed train. Section 4 draws the main conclusions.

Methodology
We collected vibration signals from vibration sensors mounted on the axial boxes and on-board monitoring system of a train. Root mean square (RMS) of original vibration signals is selected as the health index since it effectively describes the signal strength and reveals the impact of failure development on vibration energy. A careful review of RMS data states that there are two degradation phases, i.e., normal operation phase I and rapid degradation phase II, as shown in Figure 1. In phase I, the RMS fluctuates trivially, indicating a relatively healthy state. In phase II, the degradation speed grows abruptly until the occurrence of failure. The demarcation point between the two phases is defined as phase changing time (PCT).
Appl. Sci. 2020, 10, x FOR PEER REVIEW 3 of 16 phases, i.e., normal operation phase I and rapid degradation phase II, as shown in Figure 1. In phase I, the RMS fluctuates trivially, indicating a relatively healthy state. In phase II, the degradation speed grows abruptly until the occurrence of failure. The demarcation point between the two phases is defined as phase changing time (PCT). Precise detection of PCT is important since it influences the lifetime prediction accuracy. We introduce a 3σ-interval criterion for screening of PCT. Assume that RMS in phase I follows a normal distribution with mean μ and variance 2 σ . With both μ and σ obtained from historical data, the 3σ-interval of the extracted indicator can be constructed as [ 3 , 3 ] μ σ μ σ − + . We used this interval to monitor the state of the wheel tread by a comparison with latest RMS data. The wheel tread was treated normal if the RMS value fell within the interval, and abnormal otherwise. The abnormal state might be caused by either a fault or random noises. However, when multiple consecutive data were detected as abnormal, there was a high probability that a fault had occurred. In light of this, the following procedure was proposed for PCT detection (Algorithm 1).

Algorithm 1.
PCT detection procedure (1) Let 0 w = , find the first time 0 τ when the extracted RMS value exceeds 3σ-interval; (2) Let 1 w w = + , find the time w τ when the consecutive Step (2); else, output w τ as the PCT value.
According to the detected PCT, we could divide the whole degradation process of the wheel tread into two phases and the RMS values in phase II were further used for RUL prediction.

RUL Prediction Model
For a specific unit, let ( ) v t denote the degradation rate of the extracted RMS. We constructed where λ represents the basic degradation rate, ( )  Precise detection of PCT is important since it influences the lifetime prediction accuracy. We introduce a 3σ-interval criterion for screening of PCT. Assume that RMS in phase I follows a normal distribution with mean µ and variance σ 2 . With both µ and σ obtained from historical data, the 3σ-interval of the extracted indicator can be constructed as [µ − 3σ, µ + 3σ]. We used this interval to monitor the state of the wheel tread by a comparison with latest RMS data. The wheel tread was treated normal if the RMS value fell within the interval, and abnormal otherwise. The abnormal state might be caused by either a fault or random noises. However, when multiple consecutive data were detected as abnormal, there was a high probability that a fault had occurred. In light of this, the following procedure was proposed for PCT detection (Algorithm 1).
According to the detected PCT, we could divide the whole degradation process of the wheel tread into two phases and the RMS values in phase II were further used for RUL prediction.

RUL Prediction Model
For a specific unit, let v(t) denote the degradation rate of the extracted RMS. We constructed where λ represents the basic degradation rate, ϕ(t; θ) is a function of time with parameter θ, m is a constant and B(t) is the standard Brownian motion. Below are some interpretations of Equation (1). (3) Standard Brownian motion B(t) represents temporal uncertainty of the degradation process.
Accordingly, the constructed degradation rate v(t) follows a normal distribution, i.e., v(t) ∼ N(λ · ϕ(t; θ), mλ). Note that λ · ϕ(t; θ) denotes the deterministic term related to material performance and failure mode, and √ mλ · dB(t) ∼ N(0, mλ) denotes the stochastic term which describes random errors caused by inhomogeneity of material. The fluctuation caused by inhomogeneity of material is proportional to the material itself.
Taking integration on both sides of (1), we could obtain the degradation path as where x 0 describes initial value. Without loss of generality, we assumed x 0 = 0. Since measurement errors are inevitable in condition monitoring, the observation function is formulated as where ε describes measurement uncertainty, which follows a normal distribution with zero mean and variance σ 2 m . A notable feature of the above-mentioned degradation model is that it integrates multiple uncertainties extensively existed in wheel tread, which can better reflects the actual health condition. In addition, the model relates degradation variation (described by the diffusion parameter in Equation (2), i.e., √ mλ) to degradation rate (described by the drift parameter in Equation (2), i.e., λ) quantitatively. This is consistent with actual observation that a severer degradation provides an acute vibration variation [28].

Parameter Initialization and Updating Algorithm
RUL prediction starts with offline parameter estimation. A total of c monitoring time points starting from the PCT, i.e., t 1 , t 2 , · · · , t c , were selected as the training set, as shown in Figure 2. Maximum-likelihood was chosen as the estimation algorithm. (2) The function of time ( ) ; t ϕ θ is determined by specific failure mode; (3) Standard Brownian motion ( ) B t represents temporal uncertainty of the degradation process.
Accordingly, the constructed degradation rate ( ) v t follows a normal distribution, i.e., Taking integration on both sides of (1), we could obtain the degradation path as where 0 x describes initial value. Without loss of generality, we assumed 0 0 x = . Since measurement errors are inevitable in condition monitoring, the observation function is formulated as where ε describes measurement uncertainty, which follows a normal distribution with zero mean and variance 2 m σ .
A notable feature of the above-mentioned degradation model is that it integrates multiple uncertainties extensively existed in wheel tread, which can better reflects the actual health condition. In addition, the model relates degradation variation (described by the diffusion parameter in Equation (2), i.e., mλ ) to degradation rate (described by the drift parameter in Equation (2), i.e., λ ) quantitatively. This is consistent with actual observation that a severer degradation provides an acute vibration variation [28].

Parameter Initialization and Updating Algorithm
RUL prediction starts with offline parameter estimation. A total of c monitoring time points starting from the PCT, i.e., 1 2 , , , c t t t  , were selected as the training set, as shown in Figure 2.
Maximum-likelihood was chosen as the estimation algorithm.  Denote RMS values at these time points as y c = [y 1 , · · · , y c ] , where [·] represents the vector transposition, and y i = y(t i ), i = 1, 2, · · · , c. Based on (3), y i can be expressed as with , · · · , 1} c×c and MN represents the multivariate normal distribution.
Proof of Proposition 1 can be found in Appendix A. Denote unknown model parameters as For simplification, a parameter transformation is carried out for Ω, i.e., and the log-likelihood function is transformed into The first partial derivatives of ln l Θ| y c with respect to µ λ and σ 2 λ are given as Let ∂ ln l Θ| y c /∂µ λ = ∂ ln l Θ| y c /∂σ 2 λ = 0, then the maximum likelihood estimations (MLEs) of µ λ and σ 2 λ based on the monitoring data up to t c , denoted as µ λ,0 and σ 2 λ,0 , are derived as Substitute Equation (10) into Equation (8), we obtained a simplified version of the log-likelihood function with model parameters θ, m and σ 2 m . In this way, the MLEs of θ, m and σ 2 m , denoted as θ 0 , m 0 and σ 2 m,0 , can be obtained using some numerical optimization algorithms. Finally, we could obtain µ λ,0 and σ 2 λ,0 by substituting θ 0 , m 0 and σ 2 m,0 into Equation (10).
After offline training of model parameters based on the extracted RMS values up to t c , we could further carry out online RUL prediction when a new monitoring data is available. Based on the proposed model, the RUL of the unit at t i , denoted as r i , is defined as where d represents the failure threshold, According to the definition, the conditional probability density function (PDF) of r i can be derived as where Equation (12) states that the conditional PDF of r i uses the latest monitoring data instead of the entire historical data. This is conflict with the degradation evolution pattern introduced in Equation (11). To incorporate the history of observations while maintain the nice property of Wiener process model, we considered an updating procedure for parameter λ by using a linear random walk model [32], i.e., where η ∼ N(0, κ). This model is chosen because of its simplicity and practical rationality. According to the property of random walk, the hidden drift at time point t i , λ i is a time-dependent variable dependent on λ i−1 and adjusted by a noise term. Accordingly, the state space model for observation y i can be constructed as where Since the existence of hidden states in both state transition function Equation (13) and observation function Equation (14), we could not utilize maximum likelihood estimation method to update the model parameters and system states directly. To this end, we developed a recursive algorithm by integrating KF and EM for parameter update and estimation. The proposed algorithm consists of two main steps. The KF-step focuses on the estimate and update of the system state, and the EM-step aims to estimate other unknown and non-updating parameters in the observation equation.
Recall that the initial mean and variance of drift is µ λ,0 and σ 2 λ,0 , i.e., λ 0 ∼ N µ λ,0 , σ 2 λ,0 . In addition, the hidden state λ i can be estimated once the observed monitoring data up to t i , i.e., y i = y 1 , · · · , y i , are available. Denote the mean and variance of λ i by E i|i = E λ i |y i and Var i|i = Var λ i |y i respectively.
To obtain E i|i and Var i|i , we need to calculate the conditional PDF of λ i given y i , denoted by p λ i |y i .
Recursion solution of p λ i |y i can be estimated from p λ i−1 |y i−1 using the well-adopted Bayesian rule, i.e., Combing Equation (13) to Equation (15), we could obtain that p λ i |y i follows Gaussian distribution with mean E i|i and variance Var i|i , i.e., where both E i|i and Var i|i can be estimated by using the following KF-step (Algorithm 2).
(4) Estimate the filter gain: Following the KF-step, the historical monitoring data can be captured via recursively updating E i|i and Var i|i , based on which the conditional PDF of λ i , i.e., p λ i |y i , is derived.
On the basis of the estimated p λ i |y i , we need to update the unknown parameters including µ λ , , then the log-likelihood function of y i can be given as where λ i = {λ 0 , λ 1 , · · · , λ i−1 } denotes the history of the updated drift parameters, and p y i , λ i Θ denotes the joint PDF of the observed monitoring data y i . According to Equations (13) and (14), When ignoring the constant item, the log-likelihood function can be formulated by Note that Equation (18) involves hidden variables λ i = {λ 0 , λ 1 , · · · , λ i−1 }, so estimating unknown model parameters by maximizing this function is intractable. To solve the problem, we considered EM algorithm that provides a possible framework for estimating the unknown parameters involving hidden variables. A basic assumption of EM algorithm is that the hidden variables can be estimated by the monitoring data. Our problem meets the assumption since p λ i |y i was obtained by Equation (16). A detail of the EM-step for our problem is given as follows (Algorithm 3).
To calculate the conditional expectation E ln λ j . In this paper, these conditional expectations are estimated based on the properties of variance-covariance and Rauch-Tung-Striebel (RTS) smoothing algorithm [32]. Finally, the PDF of the updated RUL at the current monitoring time point t i can be obtained by using the law of total probability, i.e., where the model parameters have been derived above.
where k is the number of iterations andˆ Θ denotes the estimated parameters in the k-th step conditional on y i .

Case Study
So far, some types of train have installed on-board real-time monitoring system, based on which online vibration signals of wheel treads can be monitored. A complete set of on-board real-time monitoring system is shown in Figure 3. We could find that it includes shock vibration sensor, data pre-processor, diagnostic instrument, display terminal, etc. When the train is running, the system can output original vibration data to data analysis center via the wireless local area network (WLAN).
Appl. Sci. 2020, 10, x FOR PEER REVIEW 9 of 16 amplitudes at some monitoring points were missing and this could be reflected by the blank intervals in Figure 5.   To demonstrate the applicability of the proposed RUL framework, vibration signals of wheel treads of a high-speed train collected from the system were adopted as the data set. Figure 4 shows a schematic for the wheel tread with sensor. The marked point 11 and point 16 are installation positions of the sensors used for monitoring the wheel treads. Figure 5 presents the original vibration signals of two wheel treads. It is easy to find that the monitoring system recorded the running mileage of wheel treads and corresponding amplitude (dB) at each monitoring point. The reference value for measuring dB values was 10 −5 mm/s. Note that the recorded amplitude reflected the normalized strength of vibration shocks by eliminating the effects of running speed and axle diameter of the train. The recorded running mileage started from 402,590 km and ended at 447,870 km and 450,000 km for the first and the second wheel tread, respectively. In addition, because of human factors, the amplitudes at some monitoring points were missing and this could be reflected by the blank intervals in Figure 5.   To understand the failure processes of two wheel treads more clearly, failure mode analysis were carried out based on physical technologies. Results show that the degradation of the first wheel tread was caused by out-of-round and the second one was caused by stripping. Engineering experience shows that out-of-round and stripping are two most common failure modes of the wheel tread of train. For an illustrative purpose, an illustration of the two failure modes was provided, as shown in Figure 6a,b respectively. In particular, out-of-round is always caused by uneven wear of wheel treads. During the running state, because of various factors such as possible scratches in traction and braking,    To understand the failure processes of two wheel treads more clearly, failure mode analysis were carried out based on physical technologies. Results show that the degradation of the first wheel tread was caused by out-of-round and the second one was caused by stripping. Engineering experience shows that out-of-round and stripping are two most common failure modes of the wheel tread of train. For an illustrative purpose, an illustration of the two failure modes was provided, as shown in Figure 6a,b respectively. In particular, out-of-round is always caused by uneven wear of wheel treads. During the running state, because of various factors such as possible scratches in traction and braking, To understand the failure processes of two wheel treads more clearly, failure mode analysis were carried out based on physical technologies. Results show that the degradation of the first wheel tread was caused by out-of-round and the second one was caused by stripping. Engineering experience shows that out-of-round and stripping are two most common failure modes of the wheel tread of train. For an illustrative purpose, an illustration of the two failure modes was provided, as shown in Figure 6a,b respectively. In particular, out-of-round is always caused by uneven wear of wheel treads. During the running state, because of various factors such as possible scratches in traction and braking, foreign body bruising and so on, the wheel rolling circle loses its roundness and presents a polygon phenomenon. Vibration and noise caused by wheel polygon not only affect the comfort of the ride but also increase the vertical contact force of the wheel and rail, resulting in accelerated damage of the track and the components of train and shortening the service life. Different from this, stripping is always caused by fatigue stress, which is inevitable in all of rotary machinery systems. Tiny defects in the material will be continuously amplified under the effect of fatigue stress until part of the material is stripped from the entire wheel tread. Around the occurrence of stripping, degradation velocity of tread varies significantly because the balance of rotary system has been broken. This phenomenon will be verified in the following. Note that when the degradation caused by either out-of-round or stripping reaches a certain extent, appropriate and seasonable maintenance measures must be taken. Otherwise, the train may happen catastrophic failure during its high-speed running state. For the wheel treads provided in our case, turning maintenance was carried out at the ending points of their monitoring processes, which means that the wheel treads have been considered to be failed in practical applications. material is stripped from the entire wheel tread. Around the occurrence of stripping, degradation velocity of tread varies significantly because the balance of rotary system has been broken. This phenomenon will be verified in the following. Note that when the degradation caused by either outof-round or stripping reaches a certain extent, appropriate and seasonable maintenance measures must be taken. Otherwise, the train may happen catastrophic failure during its high-speed running state. For the wheel treads provided in our case, turning maintenance was carried out at the ending points of their monitoring processes, which means that the wheel treads have been considered to be failed in practical applications. According to the characteristics of the vibration signals, the mileage, i.e., the X-axis of Figure 5, was selected as a generalized time, which would be used in our proposed model. Meanwhile, the RMS of original vibration signals was extracted as the degradation indicator. To obtain the RMS, the recorded mileages were separated into equally distributed sliding intervals, each of which was 200 km in length. In this way, the first RMS value was calculated based on the amplitudes between 0 and 200 km, and the second RMS value was calculated based on the amplitudes between 100 and 300 km, and so on. Figure 7 shows the extracted RMS values of the two wheel treads. There were in total 389 extracted data points for the first wheel tread and 456 data points for the second wheel tread. We found that the extracted RMS values were relatively stable with some fluctuations at the beginning, but present obvious increasing trends within the following monitoring process. Therefore, we detected the PCT of each wheel tread and divided the degradation path into normal operation phase and rapid degradation phase. Results are also presented in Figure 7. It is observed that the RMS was sensitive to the incipient fault of wheel treads. Except this, we found that the second wheel tread performed a two-phase degradation behavior within its rapid degradation phase. Combined with its failure mode, we though that this was reasonable as the health states of the wheel tread could be different from each other before and after the occurrence of stripping. In this view, the changing point of the two-state degradation should be the time of stripping, and when the balance of rotary wheel tread was broken, more severe vibrations accelerated its degradation. For convenience, the changing point was denoted as stripping t in the following.
To carry out RUL prediction, we need to determine the failure threshold of wheel treads. As mentioned above, the two wheel treads were considered to be failed in practical applications. Thus, we selected the mean of the final RMS values of the two wheel treads as the threshold, i.e., 35 d = . According to the characteristics of the vibration signals, the mileage, i.e., the X-axis of Figure 5, was selected as a generalized time, which would be used in our proposed model. Meanwhile, the RMS of original vibration signals was extracted as the degradation indicator. To obtain the RMS, the recorded mileages were separated into equally distributed sliding intervals, each of which was 200 km in length. In this way, the first RMS value was calculated based on the amplitudes between 0 and 200 km, and the second RMS value was calculated based on the amplitudes between 100 and 300 km, and so on. Figure 7 shows the extracted RMS values of the two wheel treads. There were in total 389 extracted data points for the first wheel tread and 456 data points for the second wheel tread. We found that the extracted RMS values were relatively stable with some fluctuations at the beginning, but present obvious increasing trends within the following monitoring process. Therefore, we detected the PCT of each wheel tread and divided the degradation path into normal operation phase and rapid degradation phase. Results are also presented in Figure 7. It is observed that the RMS was sensitive to the incipient fault of wheel treads. Except this, we found that the second wheel tread performed a two-phase degradation behavior within its rapid degradation phase. Combined with its failure mode, we though that this was reasonable as the health states of the wheel tread could be different from each other before and after the occurrence of stripping. In this view, the changing point of the two-state degradation should be the time of stripping, and when the balance of rotary wheel tread was broken, more severe vibrations accelerated its degradation. For convenience, the changing point was denoted as t stripping in the following.
To carry out RUL prediction, we need to determine the failure threshold of wheel treads. As mentioned above, the two wheel treads were considered to be failed in practical applications. Thus, we selected the mean of the final RMS values of the two wheel treads as the threshold, i.e., d = 35. In addition, ϕ(t; θ) = θt θ−1 was selected for the 1-st wheel tread and ϕ(t; θ) = θ was selected for the 2-nd wheel tread because the extracted RMS values always increase with the running mileage in a power-law form (or a linear form) for the wheel treads whose failure mode is out-of-round (or stripping) according to our previous engineering experience. Under these assumptions, we illustrated our proposed prediction model and corresponding parameter updating algorithm based on RMS values of the rapid degradation phase.
2-nd wheel tread because the extracted RMS values always increase with the running mileage in a power-law form (or a linear form) for the wheel treads whose failure mode is out-of-round (or stripping) according to our previous engineering experience. Under these assumptions, we illustrated our proposed prediction model and corresponding parameter updating algorithm based on RMS values of the rapid degradation phase. To verify the rationality of the proposed model, we separated the degradation paths of the two wheel treads into 10 equal intervals respectively and calculated degradation rate and degradation diffusion of each interval. According to the results, a correlation test between degradation rate and degradation diffusion was carried out. Figure 8 shows the test results, which mean that there really existed a linear relationship between the drift parameters and the square of diffusion parameters, as mentioned in our proposed model.
Then, based on the parameter initialization method, we initialized the unknown parameters by using the first twenty RMS values from each changing point including two PCTs and a stripping point. Results are shown in Table 1. Next, the model parameters were updated continuously based on KF and EM algorithm with the increasing of RMS values. The parameter updating procedure is presented in Figure 9. It is observed that all parameters were not estimated accurately at the beginning of the predicting process. However, when enough RMS values were put into parameter estimation, these parameters gradually converged to fixed values. Note that, for the second wheel tread whose failure mode was stripping, the parameter updating procedure included two phases. This was consistent with the above-mentioned explanations.
To identify the effectiveness of the proposed model, we predicted the degradation state and the RUL distribution at different monitoring points of the two wheel treads. Figure 10 compares the predicted degradation path with the original degradation path. We can observe that the prediction results matched the original degradation path well except for the first twenty points, which means that our proposed model had a better goodness-of-fit with vibration monitoring data of the wheel treads. Figure 11 presents the PDFs of the RUL of the wheel treads at different monitoring points calculated by our proposed model. We found that our proposed model could cover actual RULs with an acceptable precision. To verify the rationality of the proposed model, we separated the degradation paths of the two wheel treads into 10 equal intervals respectively and calculated degradation rate and degradation diffusion of each interval. According to the results, a correlation test between degradation rate and degradation diffusion was carried out. Figure 8 shows the test results, which mean that there really existed a linear relationship between the drift parameters and the square of diffusion parameters, as mentioned in our proposed model.  Then, based on the parameter initialization method, we initialized the unknown parameters by using the first twenty RMS values from each changing point including two PCTs and a stripping point. Results are shown in Table 1. Next, the model parameters were updated continuously based on KF and EM algorithm with the increasing of RMS values. The parameter updating procedure is presented in Figure 9. It is observed that all parameters were not estimated accurately at the beginning of the predicting process. However, when enough RMS values were put into parameter estimation, these parameters gradually converged to fixed values. Note that, for the second wheel tread whose failure mode was stripping, the parameter updating procedure included two phases. This was consistent with the above-mentioned explanations.     To identify the effectiveness of the proposed model, we predicted the degradation state and the RUL distribution at different monitoring points of the two wheel treads. Figure 10 compares the predicted degradation path with the original degradation path. We can observe that the prediction results matched the original degradation path well except for the first twenty points, which means that our proposed model had a better goodness-of-fit with vibration monitoring data of the wheel treads. Figure 11 presents the PDFs of the RUL of the wheel treads at different monitoring points calculated by our proposed model. We found that our proposed model could cover actual RULs with an acceptable precision.    For comparison, we further compared the expected RULs derived by our model (defined as Model 1) and the following model (defined as Model 2), i.e., which does not consider the quantitative relationship between degradation rate and degradation variation, with the actual RULs at various monitoring points. Table 2 shows the relative errors (REs) of the predicted results. Noticeably, our proposed model incurred smaller REs compared with the model in Equation (20). Therefore, the results could demonstrate that our proposed model could help to improve the accuracy of degradation path tracking and RUL estimation. In the meanwhile, it should be mentioned that the REs at the first two monitoring points of the 2-nd wheel tread were a little big, which means that the predicted results were not quite satisfactory. This is reasonable because the two points were particularly close to the stripping point of the wheel tread and there was not much degradation information available for RUL prediction. With the increasing of degradation data, we could see that the REs of the predicted results became smaller and smaller, especially for our proposed model. For comparison, we further compared the expected RULs derived by our model (defined as Model 1) and the following model (defined as Model 2), i.e., which does not consider the quantitative relationship between degradation rate and degradation variation, with the actual RULs at various monitoring points. Table 2 shows the relative errors (REs) of the predicted results. Noticeably, our proposed model incurred smaller REs compared with the model in Equation (20). Therefore, the results could demonstrate that our proposed model could help to improve the accuracy of degradation path tracking and RUL estimation. In the meanwhile, it should be mentioned that the REs at the first two monitoring points of the 2-nd wheel tread were a little big, which means that the predicted results were not quite satisfactory. This is reasonable because the two points were particularly close to the stripping point of the wheel tread and there was not much degradation information available for RUL prediction. With the increasing of degradation data, we could see that the REs of the predicted results became smaller and smaller, especially for our proposed model.

Conclusions
This paper constructed a parameter-related Wiener process model with a recursive algorithm for RUL prediction of wheel treads. A distinguished superiority of our model was the consideration of parameter correlation and multiple uncertainties, which better fits the evolution pattern of degradation path. A recursive algorithm based on KF and EM technologies was developed to realize joint estimation and update of lifetime parameters. Experiments on real-world vibration signals of wheel treads demonstrated the superiority of the proposed methodology in prediction accuracy.
A possible extension is to develop an early warning system based on RUL prediction results. According to the safety requirements, we could design different warning thresholds so that precise time points for making maintenance decisions could be given. In other words, when the degradation value of the wheel tread surpasses a given percentage of the failure threshold, we expected that early warning would occur. This technology could contribute to safety-critical predictive maintenance plans [33].

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

Appendix A
Proof of Proposition 1. According to the statistical properties of Brownian motion, we have E(B(t i )) = 0 and D(B(t i )) = t i , where E(·) and D(·) denote the mean and variance, respectively. In this way, the mean of the health indicators at t i is clearly derived as and then E y c = [E(y 1 ), · · · , E(y c )] = µ λ t c .
For the covariance between y i and y j , we calculate it by using Cov y i , y j = Cov λ where each item can be given as Cov √ mλB(t i ), Denote the results in matrix and we derive D y c = σ λ 2 t c t c + mµ λ Q c σ m 2 E c . (A7)