A knowledge-based prognostics framework for railway track geometry degradation

This paper proposes a paradigm shift to the problem of infrastructure asset management modelling by focusing towards forecasting the future condition of the assets instead of using empirical modelling approaches based on historical data. The proposed prognostics methodology is general but, in this paper, it is applied to the particular problem of railway track geometry deterioration due to its important implications in the safety and the maintenance costs of the overall infrastructure. As a key contribution, a knowledge-based prognostics approach is developed by fusing on-line data for track settlement with a physics-based model for track degradation within a filtering-based prognostics algorithm. The suitability of the proposed methodology is demonstrated and dis- cussed in a case study using published data taken from a laboratory simulation of railway track settlement under cyclic loads, carried out at the University of Nottingham (UK). The results show that the proposed methodology is able to provide accurate predictions of the remaining useful life of the system after a model training period of about 10% of the process lifespan.


Introduction
In most developed countries, the continuous ageing and the growing demand of use of critical infrastructures calls for advanced Prognostics and Health Management (PHM) concepts for optimal infrastructure asset management [1]. In particular for the railway infrastructure, the continuous deterioration of the track due to traffic loading represents the main ageing factor requiring periodic interventions to restore the track to an acceptable geometry [2,3]. These maintenance interventions not only represent a significant part of the railway operation expenses, but also imply temporary line closures and a reduction of the effective network capacity, so they need to be planned months in advance. It is in this context of anticipated decision making where the benefits of PHM can be fully exploited for maintenance cost optimization and downtime reduction.
As evident from the literature, railway track degradation and maintenance modelling to date has a strong empirical retrospective character, mainly grounded on data-based (stochastic or phenomenological) models with limited prospective capability. A systematic review of these models is proposed by Dahlberg [4] for track settlement modelling, and more recently by Soleimanmeigouni et al. [5] and Higgins and Liu [6] focusing also in track maintenance modelling. The prediction accuracy of those models strongly depends on the quality and quantity of the available historic data, and thus they are prone to misjudgements especially under medium-to-long term future scenarios. To overcome this limitation, some physics-based models have been proposed to predict the progressive mechanical degradation of the track under cyclic loadings from first geomechanical principles, like the cyclic densification model by Suiker and de Borst [7], and the elastoplastic models by Indraratna et al. [8] and more recently by Sun et al. [9]. The referred physics-based models are transparent to human understanding and they can adapt to different material and loading conditions without much training, which is essential to confer the required early predictive capability. However, they are unable to account for the uncertainty in the predictions since they are based on deterministic input-output relationships. In this context, a number of authors have adopted statistical methods to quantify the modelling uncertainty in track geometry degradation, such as regression analysis [10], stochastic processes [11][12][13], Markov models [14,15], and Bayesian analysis [16], to cite but a few. Notwithstanding, as stated before, the accuracy of these methods relies almost entirely on historical data, which bounds their prospective capability and the reliability of their predictions.
In this paper, a paradigm shift for track geometry degradation modelling and maintenance is proposed. Instead of making maintenance decisions based on either a retrospective data-based modelling of the track, or, alternatively, using a purely deterministic physicsbased approach, a knowledge-based prognostics framework is proposed herein. This approach fuses information from physics-based models and available data about track degradation within a Bayesian learning paradigm to sequentially reduce the initial modelling uncertainty [17] as long as new data are collected, so as to obtain increasingly accurate forecasts of the future condition of the track. These forecasts enable the testing of various load and utilisation scenarios, and more importantly, allow us to answer the question of "when will failure occur" with quantified uncertainty. The proposed methodology enables informed, anticipated, and risk-based decisions about optimal railway track asset management, and can be extended to other railway assets to enable decision-making at a system infrastructural level. To the authors' best knowledge, the existing literature on track degradation prognostics is currently limited to Mishra et al. [18], focusing on railway switches based on a phenomenological model for track geometry degradation.
For the problem of physical modelling of the track deterioration, an elasto-plastic model for track settlement is developed based on first principles of the plasticity of soils [19] and Critical State Soil Mechanics (CSSM) [20]. This model provides the temporal evolution of the elementary volumetric and deviatoric plastic strains of the track based on some geomechanical inputs. Apart from the traffic loads, the inputs required by the model reduce to the critical state parameters [20], which can be obtained from routine triaxial tests, and two uncertain parameters controlling the cyclic hardening of the granular layers of the track, which are sequentially updated using condition monitoring data. For the sake of Bayesian data assimilation, the proposed model is embedded within a hidden Markov model [21] where the deviatoric and volumetric plastic strains constitute the latent states, and the vertical permanent axial strain is the observed state. Next, a sequential Bayesian updating framework is used to update both the latent states and the uncertain model parameters as long as new data become available. The resulting equations for sequential updating are analytically intractable, hence a Sequential Monte Carlo approach based on particle filters [22] is adopted whereby the probability density functions (PDFs) are approximated by weighted samples or particles, which represent random realisations of the updated track degradation states. Finally, by extrapolating those particles into the future, a probability-based estimation of the time (load cycles) to reach a predefined functional limit for track degradation is subsequently obtained. This time defines what in the Prognostics and Health Management community is commonly known as the remaining useful life (RUL) [1] of the system.
As a case study, the proposed methodology is tested against experimental data taken from Aursudkij et al. [23] about permanent axial strain in a ballasted railway track, carried out at the Nottingham Railway Test Facility [24]. The results show that the proposed prognostics methodology is able to accurately anticipate the future states of deformation of the track after a training period of about 10% of the total length of the process, which corresponds to the time required by the model to assimilate the data. This prognostics performance is compared with the performance obtained for the same dataset using an empirical logarithmic model for track settlement. The assessment is carried out using a newly developed prognostics metric named as the relative prospective evidence, which measures the relative prediction accuracy of two different model classes in relation to the amount of information that the models need to extract from the data to perform their predictions. The results show that the physics-based model proposed in this paper is the one with the larger and earlier anticipation capacity, needing less support from condition monitoring data to make accurate prognostics estimates.
The remainder of the paper is organised as follows. Section 2 presents the physical fundamentals behind the proposed model for track deterioration. Section 3 overviews the mathematical basis and computational aspects of Bayesian state estimation and prognostics, which have been both specialised for their application to the railway track degradation problem. In Section 4, the proposed framework is applied to railway track settlement data to serve as a case study. A discussion about the suitability of model-based prognostics for railway track maintenance planning is provided in Section 5. Finally, Section 6 provides concluding remarks.

Physics-based model for track settlement
A constitutive model for ballast plastic deformation under cycling loading is developed here following the general theory of plasticity, CSSM, and Pender's postulates about plastic deformation of overconsolidated soils [25]. The model predicts the evolution of the plastic vertical strain p 1 of a representative track section, which is the variable of interest here, whereby the settlement of the section can be straightforwardly derived. The reader is referred to [25] for a detailed overview of the foundations of the proposed model, but for the sake of clarity, the key formulation and definitions are reproduced here with a uniform notation.

Physical fundamentals
Following the theory of plasticity, a granular medium under a stress state represented by the stress tensor σ ij , will experience an increment in the plastic strains given by Hill's flow rule [26], as where d ij p is the corresponding plastic strain increment tensor, df is the differential of function f, which denotes the yield criterion (or yield surface), g is the plastic potential, and h is the hardening function. These elements are briefly explained in Definitions 1 to 3, respectively. If the stress invariants p and q are adopted instead of the stress tensor σ ij , Eq. (1) can be rewritten as: with v p and s p being the corresponding plastic volumetric and distortional strain invariants. The reader is referred to Appendix A for further insight about stress and strains invariants. Having v p and s p as known, the vertical plastic strain of the track section can be straightforwardly computed as: as shown in Eq. (A.4b) in Appendix A. It follows that functions f, g and h need to be determined so as to define a closed-form expression for the constitutive model in Eq. (2). [19]). The yield criterion defines the limit of the elastic behaviour and the onset of plastic deformation under any possible combination of stresses. Mathematically:

Definition 1 (Yield criterion
Assuming Pender's hypothesis that overconsolidated soils experience plastic strain only when there is a change in the stress ratio ( = q p / ) [25], the yield criterion can be assumed to be given by: where η i defines a particular constant stress ratio yield locus i (refer to Fig. 2 in [25]). By differentiating Eq. (4), = df dq dp, i and substituting = + dq dp pd i obtained from the definition of η, the required expression for df in Eq. (2) can be shown to be given by: Definition 2 (Plastic potential [19]). The plastic potential is given by the generic expression = g p q ( , ) 0 and relates the incremental plastic strains with the stress state following the plastic flow rule in Eq. (1).
An analytical expression for g can be obtained by noting that the plastic strains are normal to the plastic potential = g p q ( , ) 0, thus implying that [19] where M is the slope of the critical state line in the p q space [28], and η* is a function of the stress ratio given by = p p * / , cs with p cs being defined in Eq. (13) below. The stress-dilatancy criterion in Eq. (7) differs from that proposed by Rowe [27] in that it provides positive values for plastic volumetric strain increments (i.e., contraction), as long as η* does not exceed the critical state value M. By substituting Eq. (7) into (6) and solving the resulting differential equation, an expression for the plastic potential function is obtained, as: where C is an integration constant. [19]). As a work hardening material, ballast experiences a dynamical change of its plastic response during yielding. This change is accomplished by modifying the yield surface as plastic flow occurs. The hardening function (h in Eq. (2)) accounts for such effect.

Definition 3 (Hardening function
An expression for h can be obtained following the premise that no volume change takes place during undrained deformation [20]. Mathematically: where d v e is the elastic volumetric strain increment given by [20] = + d dp  (11) In Eq. (11), Pender's assumption of parabolic undrained stress paths in the p q space [25] is adopted to relate the stress derivatives dp and dη. This assumption is expressed as: where η 0 is the initial stress ratio (corresponding to the minimum load within the cycle), p 0 is the initial value of the mean stress invariant p, and p cs is the value of p at the point on the critical state line corresponding to the current voids ratio e, obtained as [28,29]: = p e exp cs cs (13) with Γ and λ cs being material parameters. These parameters, together with M and κ (see nomenclature Table A.3), are the critical state parameters of the model and can be evaluated from a series of drained triaxial compression tests [30]. Having obtained the differential dp using Eq. (12), and substituting it into Eq. (11), the hardening function becomes after some algebraic manipulation 1 : It should be emphasised that the last equation is nominally valid for the first load cycle. For subsequent loading, h should account for other complex phenomena observed in experiments under cyclic loading conditions, such as the Bauschinger effect, the effect of the stress ratio and loading history, etc. [19]. A semi-empirical approach is adopted here to account for such processes following Indraratna and Salim [8], which essentially consists of multiplying the hardening function in Eq. (14) by a semi-empirical factor ϕ given by where α and β are empirical fitting parameters, n is the accumulated number of cycles, ⟨ · ⟩ are the Macauley brackets, and p e is given by the expression [8]: with p min and p max being the minimum and maximum mean stress invariant within load cycle n, respectively. The corresponding deviatoric components q q , min max and q e are obtained using Eq. (A.4a), defined in Appendix A.

Constitutive equations of proposed model
Substituting the expressions for f, g and h obtained above into Eq. (2), and after some algebraic manipulation, the following expressions for the evolution of the plastic volumetric and deviatoric strains are obtained: where · | n denotes a variable evaluated at load cycle n, and = q p / max max max . Note that for the nth loading cycle, the initial condition for both deviatoric and volumetric plastic strains is given by the corresponding plastic strains at the previous cycle n 1. The initial conditions ( ) v p 0 0 and ( ) s p 0 0 are assumed as known. Observe also that the expressions in Eq. (17) are non-linear differential equations of the plastic strains, thus they need to be integrated numerically by replacing the derivatives by their finite difference approximation (e.g., using explicit Euler method). An algorithmic description of the numerical method proposed to integrate the constitutive Eqs. (17a) and (17b) is presented in Algorithm 1.

Track-degradation prognostics
Prognostics is concerned with predicting the future state of health of engineering systems or components given the current degree of wear or damage, and, based on that, estimating the remaining time in which the system is expected to continue performing its intended function within desired specifications [1]. The predictions are carried out using no additional evidence but a sequence of measurements up to the time of prediction, and the available models for system behaviour. A physicsbased filtering-based prognostics framework is proposed herein as depicted in Fig. 1. The definitions in Section 3.1 are used to provide the 1 In the normal range of stresses, ballast remains in states denser than the critical state [29], thus, the sign of the hardening function in Eq. (14) is reversed. Set min and max deviatoric stresses: q min , Set min and max mean stresses: p min , Set initial stress state q 0 , Set initial voids ratio: e i=1 ← e| n 6: Compute stress state: p i ,

15:
Obtain total volumetric strain: end for 18: Step-by-step description of model for track settlement.
required mathematical background about prognostics in the context of the track degradation problem investigated here.

Fundamentals about model-based prognostics
Definition 4 (Hidden Markov model [21]    The future states prediction are represented using shaded areas for the 5% 95% and 25% 75% probability bands.
to be conditionally independent given the states x n .
A discrete-time state-space representation is typically adopted to describe the hidden Markov model, as follows [1]: where x is the state transition equation, represented by a possibly nonlinear function of the state variable x n along with a set u n u of input parameters to the system (loadings, geomechanical inputs, etc.). The function × × g : n n n n n u x y denotes the observation equation and provides the expected observation or measurement y n given the hidden state x n . Both functions f and g may depend on a set of uncertain model parameters n (e.g., fitting parameters). The terms n n x and w n n y in Eq. (18) are stochastic variables that represent the modelling error and the measurement error, respectively. Supported by the Principle of Maximum Information Entropy [31], ν n and w n can be conservatively modelled as zero-mean Gaussian distributions [1]. Thus, the state space model defined in Eq. (18) can be probabilistically rewritten as 2 : where × n n x x and × w n n y y are the covariance matrices of the model error and the measurement error, respectively. Note in Eq. (19) and in what follows that the conditioning on the model inputs u n is dropped for simpler notation.
In this work, the (latent) state of the system at loading cycle n is assumed to be described by the plastic deviatoric and volumetric strains, i.e., for the measurement error, with σ ν,s and σ ν,v being respectively the standard deviations of the deviatoric and volumetric components of the model error, and σ w the standard deviation of the measurement error. With these settings, the deterministic physics-based model proposed in Section 2 is stochastically embedded within a hidden Markov model by defining the functions f and g in Eqs. (19a) and (19b), as: where functions f s and f v in Eq. (20a) are given by In the last equation, s p n and v p n are the increments in the plastic deviatoric and volumetric strains within loading cycle n, respectively. These strain increments can be obtained as: where subscripts (s, v) denote the deviatoric and volumetric components, respectively, and N S is the number of discretising steps of the loading ramp at cycle n (the unloading ramp is regarded as elastic). The reader is referred to Algorithm 1 for details about the numerical computation of the plastic deviatoric and volumetric strains for load cycle n. for those measurements, which is given by Bayes' Theorem as [1]:

Definition 5 (Sequential state estimation
Observe from Eq. (24) that model parameters θ n are virtually timeevolving although they are essentially not dependent on time, which is a key problem that typically arises when sequentially updating the state variable as an augmented state [32]. A common solution is to add a white-noise perturbation to the set of updated parameters at time n 1 before evolving to the next predicted state at time n, i.e., = + , This induces a Markoviantype artificial dynamics to the model parameters whereby the required PDF ( ) n n 1 is prescribed [32,33], as follows: where is a covariance matrix specified here as n n n j n n ,1 , , , being the variance of the random walk of the jth component of the parameter vector θ at time or cycle n. Note that such artificial evolution imposes a loss of information in θ over time (e.g., increasingly larger spread in π(z 0:n |y 1:n )) since additional uncertainties are artificially added to the model parameters. Several methods have been proposed in the literature to overcome this drawback, with the most popular being those that impose a shrinkage over n as long as new data are available [32]. An efficient method of this class is adopted in this work [34], which consists in sequentially modifying the variances = j n , 1, , , 2 n j , as follows: where RMAD(θ n, j ) is the relative median absolute deviation of π(θ n, j |y 1:n ), RMAD* j is the target RMAD for π(θ n, j |y 1:n ), and P* [0, 1] j is a scaling constant that controls the speed of convergence to RMAD* j . Both RMAD* j and P* j need to be tuned by the modeller. A comprehensive discussion about the optimal choice for P* j and RMAD* j is found in [34].
Definition 6 (Predicted system state). Given the most up-to-date assessment of the system at time or cycle n, π(z n |y 1:n ), a probability-based prediction of the future states of degradation of the system ℓ-steps forward in time can be obtained by Total Probability Theorem as [1,35]: where ℓ > 1 is the number of prediction steps in absence of new data, and z z ( ) t t 1 is the state transition equation defined in Eq. (24), which, for t > n, encloses the future behaviour of the system.
To numerically solve this multi-dimensional integral, an approximation can be readily obtained by conditional sampling, using recursively the one-step transition equation defined in (24) Definition 7 (End of Life (EOL) [1]). The EOL is defined as the earliest time or loading cycle t > n when the predicted event z [ ] t takes place, where z t ∼ π(z t |y 1:n ) and is the subset of states where the system behaviour becomes unacceptable.
In the context of the proposed track degradation problem, represents the set of values for the vertical plastic strain exceeding the threshold , which delimits the boundary between the failure region and the safe region = .
Observe that since the predicted states z t are uncertain (recall Eq. (27)), the EOL will be an uncertain variable distributed according to the PDF π(EOL n |y 1:n where k is a normalising constant and z ( ) t is an indicator function that assigns the unity if z ( ), t and makes zero the rest.

Definition 8 (Remaining
Useful Life (RUL) [1]). The RUL is the remaining time in which the system is foreseen to perform within the region of acceptable performance, where = , the non-empty subset of authorised states of the system. A prediction of the RUL from the time of prediction n can be straightforwardly obtained from EOL n as = n RUL EOL n n . Definition 9 (Prospective evidence for a model). Given a degradation process n observed up to current time or load cycle n, and a particular model class which idealises such a process (e.g., the one proposed in Section 2), the prospective evidence for model class is the probability of the subsequent degradation process, denoted as +, n to be predicted by . The prospective evidence can be obtained by Total Probability Theorem as: Note that for this definition, data for the future degradation process + are assumed to be known, which holds true for offline and pseudoonline prognostics analyses, but it could be a hard assumption in a purely online prognostics scenario. In that case, we assume that + can be available through experiments run in similar conditions or by the existence of repeated tests.

Algorithms for prognostics
The methodology for track degradation prognostics explained above involves the evaluation of multidimensional integrals (recall Eqs. (23) and (27), respectively) which are analytically intractable, except for some of especial linear cases using Gaussian uncertainties [22]. Approximating methods like Sequential Monte Carlo [22] are required to overcome those integrals by numerical approximation. In prognostics, these approximating methods are typically integrated within algorithms which (1) sequentially estimate the joint PDF of states and model parameters as long as new data are collected, referred to as particle filtering algorithms, and (2) extrapolate the updated state estimations into the future in absence of new observations until the EOL is reached. Depending upon implementation both algorithms may not be separable, however, for the sake of clarity, they will be presented here separately in Sections 3.2.1 and 3.2.2, respectively.

Sequential state estimation algorithm
With particle filters, the PDF of the most up to date state of degradation of the track π(z 0:n |y 1:n ) can be approximated through a set of N weighted particles or samples where δ is the Dirac delta and . The value of the particle can be obtained by sequential importance sampling [22] each time a new data point arrives, as explained in the pseudocode given as Algorithm 2. Note that a systematic resampling step is proposed in Algorithm 2 (see line 15) to avoid the well-known weight degeneracy problem associated to the sequential importance sampling method [36]. If necessary, a control step on this degeneracy by using the effective sample size (ESS) [36] may be incorporated before the resampling step. The interested reader is referred to Orchard and coworkers [35] for further details about sequential state estimation for model-based prognostics.

Particle-filtering based prognostics
Using the particle-filtering approach presented above, the goal is to obtain an up-to-date estimate of the EOL of the track as the PDF of the earliest time when the predicted degradation states reach the failure region . To this end, a particle filter approximation of the predictive Eq. (27) is needed, which can be shown to be obtained as [1]: Observe that the last equation cannot be solved analytically, however it can be sampled by drawing one conditional sample sequence = , {updated particles at time n} Algorithm 2. SIS particle filter with parameter adaptation.
being an indicator function that assigns the unity if + z ( ), n and makes zero the rest. An algorithmic description for filtering-based prognostics is provided below as Algorithm 3.

Prognostics metrics
A prognostic performance evaluation is required to assess the goodness of the EOL/RUL estimates so as to avoid subsequent maintenance decisions based on poor predictions, which may increase system risk. Two main attributes are typically used to evaluate the performance of a prognostic estimation, namely [37]: (a) correctness, which is related to the prediction accuracy when compared with observed (future) outcomes; and (b) confidence, which deals with the uncertainty in the EOL/RUL predictions. Based on these attributes, a prognostic error measure where denotes the mathematical expectation, std is the standard deviation, and EOL* represents a benchmark EOL, which is considered as an uncertain variable with known statistics (e.g., EOL* might be represented by a Gaussian PDF with known mean and standard deviation). Observe from Eq. (35) that n is contributed by the relative discrepancy between the expected values of EOL n and EOL* (first term), which accounts for the correctness of the prediction, and the relative difference between the spreads in the PDFs describing the EOL n and EOL* (second term), which accounts for the confidence of the prediction.

Remark (Non symmetry).
In situations where the computed PDF of EOL n lacks of smoothness or symmetry, n is preferably formulated using the median (Md), as a measure of location, and the interquartile range (IQR), as a measure of spread, as follows: In addition to quantifying the prognostics performance for a particular model class or a particular algorithm, one may be interested in assessing and ranking two different model classes, e.g., and , according to their relative prognostics performance. To this end, a novel relative prognostics metric is proposed herein based on the ratio between the prospective evidence for each candidate model class (refer to Definition 9), as follows: where RPE n stands for the relative prospective evidence at time or load cycle n, which provides a measure of the relative anticipation capacity of the two model classes at different time instants during the process, such that RPE n ≶1 when + + ( , ) ( , ). The reason for proposing a relative performance metric ratio based on evidence is because the computation of the evidence has been shown to automatically enforce a quantitative expression of the Principle of Model Parsimony or Ockham's razor [38,39], so the extremes of overfitting (at the cost of too much information extracted from data), or under-fitting of the data become naturally penalised. In other words, if the future degradation process + n is predicted equally well by the two models, then the "simpler", i.e., the one which needs less support from the data, would be favoured by the proposed RPE metric. More insight about the proposed prognostics metric will be given below in the context of a case study.
Remark. When dealing with online prognostics, data about the future degradation process + n might be uncertain and known only through a PDF ). In this case, the RPE metric can be obtained as a mathematical expectation, as The multidimensional integral in Eq. (38) can be numerically evaluated by the Monte Carlo method as with + ( ) being T sample degradation trajectories drawn from strain in a ballasted railway track taken from the literature [23]. The test, as reported by Aursudkij et al. [23], was conducted on the Railway Test Facility (RTF) of the University of Nottingham (UK) [24]. [KPa] (σ 3 represents the mean pressure due to geostatic stresses in the railway track section). The settlement of the central sleeper was measured using LVDT sensors at a set of non-regularly scheduled load cycles, whereby the permanent axial strain data were obtained by dividing the measured settlements by the ballast section height prior to test (0.3 [m]). A summary of the dataset used for this case study is provided in Table 1. The reader is referred to Aursudkij et al. [23] for further information about the experimental setup, and to Brown et al. [24] for a detailed description of the Nottingham RTF.
Following the methodology for sequential state estimation presented in Section 3.2, probability-based predictions about the evolution of the permanent axial strain are obtained. The results, which are presented in Fig. 2 for several time instants (load cycles) during the loading process, show a good agreement between the track settlement as measured by the LVDT sensors and that estimated by the particle filter algorithm. For these calculations, the voids ratio of the intact ballast was set to = e 0. 8,0 and the values of the critical state parameters shown in Table 2 were adopted to evaluate the underlying constitutive geomechanical model. These critical state parameters have the desirable property of being independent of the state of the material (e.g., voids ratio and stress condition), and, for this case study, they have been obtained through initial fitting tests with the data. In the engineering practice, they can be obtained from routine triaxial tests. Here, the model error term ν defining the state transition equation is assumed to be the same for both the volumetric and deviatoric components, therefore = diag ( , ) in Eq. (19a), where σ ν is an uncertain parameter that is sequentially updated with the data, as shown further below. . These parameters resulted to be the most sensitive model parameters after a preliminary sensitivity analysis. A uniform PDF is conservatively adopted for the prior PDF of the model parameters π(θ j ), = j 1, 2, 3, as specified in Table 2 (see fourth column), as a way of representing our prior state of ignorance about the initial values of such uncertain parameters.
At each prediction step, which corresponds to the loading cycles when settlement data are collected, Algorithm 2 is run by using = N 5000 particles to obtain an updated estimation of both the permanent axial strain and the uncertain model parameters. For this algorithm, the scaling variables = RMAD RMAD * 0.1· j j 0, and = P* 0.9 j (recall Eq. (26)) are adopted to control the artificial dynamics of the model parameters = = { } , j j 1,2,3 where RMAD 0, j is the relative median absolute deviation of the marginal priors π(θ j ). The diagonal elements j 0, of the covariance matrix 0 in Eq. (25) are appropriately selected through initial test runs and set to 5% of the 5th-95th inter-percentile range of the marginal priors π(θ j ). To reveal the sequential uncertainty reduction in the updated model parameters, the updated mean estimation of the jth component of θ, as well as its 25% 75%, 5% 95% uncertainty bands, are plotted against the loading cycles in Fig. 3 for = j 1, 2, 3. Next, based on the most up-to-date predictions of the permanent axial strain of the track, Algorithm 3 is used for EOL and RUL calculation. To this end, the time-ahead predictions of the track settlement, which are depicted in Fig. 2 using shaded areas to illustrate the prediction uncertainty, are extrapolated forward in time until the permanent axial strain crosses the threshold between the useful domain and the failure domain , whereby the PDF of EOL is obtained as a first-passage time distribution problem. For this example, the threshold is set to . To speed up the EOL calculation, which requires the simulation of = N 5000 particle trajectories thousands (or even millions) of loading cycles ahead until reaching the failure domain, the "natural" load cycles are converted to duty cycles by dividing the load cycles by a fixed factor > DC 1 [40]. Here, a duty cycle is defined as an appropriate fixed amount of load cycles where a sensible amount of settlement is accumulated. After some pilot tests, = DC 50 was revealed as a suitable value for this example. Fig. 4a shows the sequence of estimated PDFs of EOL for the different load cycles when track settlement data are available. By comparing between consecutive PDFs in Fig. 4a, one can observe that the mean EOL prediction gradually tends to converge to the reference mean value for EOL* (blue marker in Fig. 4a), and the prediction uncertainty (i.e., the spread of the PDFs) tends to decrease, as a consequence of the sequential updating of the model as long as new data are collected. The gradual increase in the prediction accuracy is also revealed in Fig. 4b, where the mean RUL estimations are plotted against time of prediction using an accuracy cone type of representation [41]. Two cones of accuracy at 10% and 20% of the reference RUL, denoted here as RUL* and obtained as = n RUL* EOL* , are included to visually interpret the prognostics performance in regards to prediction accuracy. Observe that the expected RUL prediction is appreciably inaccurate for the initial 5000 load cycles (values out of accuracy cones), which suggests that a training stage is required by the model to learn from the data. This learning period is also manifested in Fig. 4c, where the proposed prognostics error n in Eq. (36) is plotted against loading cycles assuming a Gaussian PDF centred at 75,000 and 10% c.o.v. for EOL*.
Finally, the prognostics performance of the track geometry degradation model proposed in Section 2 is compared with the performance of a well-known class of phenomenological models for track settlement [30,42,43], given by the logarithmic model = + A B n log , p 1 where A and B are fitting parameters. For comparison purposes, this model is referred to as model class , 0 while the one proposed in Section 2 is denoted as model class 1 . A discrete-time state-space representation of model class 0 is obtained and subsequently embedded within the SIS particle filtering algorithm (refer to Algorithm 2), following the methodology in Section 3.1. respectively. The same uncertainty structure (i.e., modelling and measurement errors) that is assumed for model 1 is also assumed for 0 . The assessment of the two model classes is carried out by computing the RPE n metric proposed in Section 3.3. The results are shown in Fig. 5 for the load cycles where measurements are available. In view of Fig. 5, the physicsbased model class 1 proposed in this paper is revealed to be the one with the larger and earlier anticipation capability (RPE n > 1), since it provides larger RPE n values in the earlier stages of the process. More insight and discussion about prognostics performance is provided in the next section.

On the case study results
The proposed methodology for knowledge-based prognostics for railway track degradation has been exemplified using the case study presented in the previous section based on data from the Nottingham Railway Test Facility. As a first interpretation of the results in Fig. 2, the proposed methodology is able to accurately anticipate the evolution of the plastic vertical strain of the track with quantified uncertainty after an initial learning stage using limited data. The length of such a learning stage is revealed in Fig. 3, where the updated mean of the uncertain model parameters is shown to reach an asymptotic behaviour after the first 5000 load cycles (about 7% of the reference EOL*), which denotes that after this learning period the model has extracted enough information from the data to make accurate predictions. This is in agreement with the results shown in Fig. 4b for the RUL estimations, and also with those in Fig. 4c, where = n 5000 is revealed as the time when the prognostics error measure n starts a slower and gradual decrease; i.e., the predictions become increasingly accurate for cycles n > 5000, therefore making a decision based on the predicted EOL n is recommended from that time. Note also in Fig. 4c that the overall prognostics error n is mostly contributed by the uncertainty in the predictions (second term in Eq. (35)), while the error term accounting for the average prognostics accuracy becomes almost negligible after the aforementioned learning period. This early correctness of the prognostics results is a consequence of having adopted a physics-based model for track deterioration instead of a data-based approach. At this point, it is worth mentioning that the learning period required by the model, and consequently, the actual anticipation capacity of the prognostics algorithm, will certainly depend on the early availability and quality of the condition monitoring data, but also on the prior uncertainty of model parameters and their relative "importance" in relation to the overall model response [44]. This model parameter importance is expected to be lower for physics-based prognostics formulations, like the one proposed here, than for purely data-based approaches, where the model response becomes fully dominated by the Fig. 3. Trace of the mean values of model parameters = { , , } against load cycles. Shaded areas represent the 25% 75% (darker color) and 5% 95% probability bands, respectively.
(uncertain) value of the parameters [45]. This point is reinforced by the values RPE n > 1 in Fig. 5, which reveals that the proposed physicsbased model performs better predictions than the data-based logarithmic model used as a benchmark (denoted as 0 ), needing less support from the data to train the fitting parameters. This tendency is more accentuated in the earlier stages of the degradation process due to the initial lack of condition monitoring data. As the process evolves towards the end, both models have extracted enough information from the data to allow them to perform similarly.
In terms of the prediction uncertainty, observe from Figs. 2 and 4 a that the spread in the future predictions, which is initially high due to the high prior uncertainty in the model parameters (especially in the hardening parameters α and β), tends to decrease as the model is sequentially updated with new data, which is a common and desirable feature in any prognostics problem. However, it is also observed that the reduction of this prediction uncertainty does not take place gradually along the full extent of the process but concentrated at the initial stage, approximately coinciding with the learning period of the model where most of the data are available. In this respect, it should be pointed out that the observed uncertainty reduction is a particular consequence of the data collection pattern in this case study, i.e., collected over a set of non-regularly scheduled load cycles mostly concentrated at the beginning of the process, and not a general feature of the proposed physics-based prognostics framework.
Besides the modelling uncertainties discussed above, there is an additional source of uncertainty coming from the prognostics algorithm itself, which is related to the asymptotic approach of the predicted plastic axial strain to the established threshold (as observed in Fig. 2). This asymptotic behaviour makes that many sample predictions need very long simulations to reach the threshold, which leads to poor EOL estimations unless the amount of particles (N in Algorithm 2) is sufficiently large. A suitable solution to overcome this issue is by the adoption of dedicated prognostics algorithms based on high-efficiency sampling techniques, like the one recently proposed by the authors in [46] based on SubSet simulation, or the one developed by Yan et al. [47] based on adaptive Lebesgue sampling, among others. This task constitutes a desirable future improvement of this work.

On the extensibility to track maintenance
As stated in Section 1, the degradation of the railway track geometry results in a key asset management problem for the railway industry with important implications in safety and cost. The potential and also the limitations of the proposed prognostics framework in positively contributing to optimal track asset management is discussed in this section.
First, the authors remark that an apparent limitation of the proposed approach in the context of railway track asset management is that a physical variable such as the plastic strain (or settlement) of the track is adopted instead of more common management variables like the track geometry index [48] or the standard deviation (SD) of the settlement [3]. However, the aforementioned management variables can be expressed as functions of the settlement of the track [49] and these can be included within the hidden Markov model in Eq. (18), whereby the proposed filtering-based prognostics methodology is subsequently derived.
Next, in regards to the applicability of the proposed methodology to a real life scenario, it should be pointed out that the physical model described in Section 2 provides the evolution of the permanent settlement of the ballast as a pseudo-plane strain problem. This implies that the railway track is assumed to (approximately) have homogeneous material and geometrical characteristics along the longitudinal direction of loading. For a real world railway track, such assumption would lead to the need for a segmentation of the track into different homogenous sections where the proposed prognostics framework can be applied. Those sections would constitute predictable components of the track asset for a specific railway line, and other assets related to the track maintenance such as the drainage, switches and crossings might be taken into account to define the length of those sections. In this context, it should be emphasised that the proposed methodology is general and as such, it can be extended to not only other railway assets, but also to other engineering systems and infrastructures. This can be achieved simply by considering the appropriate degradation models and condition monitoring data, and embedding them within a hidden Markov model, as explained in Section 3.1, where the subsequent prognostics methodology is derived. From this standpoint, state-of-theart methodologies for infrastructure asset management like Petri nets [50,51] can be adopted to perform modelling activities such as optimal maintenance planning and life-cycle cost analyses at system level based on prognostics, which constitutes an immediate next step of this research.
Finally, to illustrate the potential of the proposed methodology for improved maintenance decisions and optimum asset availability, a conceptual scheme has been provided in Fig. 6. In this figure, a degradation process, which is represented using a black solid line, is being monitored from = n 0 up to present time n. At this time, a probabilistic prediction of the future evolution of the process is performed whereby a PDF of the EOL is obtained given a predefined threshold that delimit the end of the safe region and the onset of the failure region . Based on this prediction, a decision is made on how much time a corrective maintenance activity can be delayed based on the risk of reaching the failure region during the delay time. To illustrate this decision-making process, two foreseen intervention scenarios are represented in blue and red considering two sample delay times, τ and τ′ respectively, with τ′ > τ. Observe that if the intervention is delayed until + n (scenario 2), a longer lifetime extension is obtained (red PDF) in relation to that obtained for the shorter delay scenario (blue PDF), but it is at the cost of assuming a higher probability of failure (red shaded area in the EOL PDF). Conversely, the shorter delay time scenario (in blue) is more conservative in terms of probability of failure, but it leads to a shorter lifetime extension. This simple example reveals a relevant trade-off between lifetime extension and reliability which is typically faced by infrastructure owners and operators based on experience, which can be rigorously assessed and quantified with the proposed prognostics methodology for infrastructure asset management.

Conclusions
A knowledge-based prognostics methodology for railway track asset management has been developed in this paper. The proposed methodology is general but in this paper it was applied to the railway track due to its impact on the safety and the maintenance cost of the overall infrastructure. A geomechanical model for ballast settlement has been provided based on first principles and postulates about the yielding of granular materials. This model is embedded within a particle filtering algorithm for sequential state estimation and prognostics. The algorithm obtains probabilistic estimations of the remaining useful life (RUL) of the railway track as the PDF of the time to reach a threshold for tolerable vertical axial strain. A case study has been presented using experimental data for track settlement from the literature to illustrate and discuss the proposed methodology. The results show that the proposed prognostics framework is able to provide reasonably accurate predictions of the RUL after a short model training period. More research effort is needed to exploit the potential of prognostics in the context of optimal track asset management at whole-life whole-system level by integrating prognostics within state-of-art methodologies for infrastructure asset management. Fig. 6. Illustrative scheme on prognostics-based maintenance. The blue and red dashed lines represent the predicted system behaviour under two alternative intervention scenarios. Observe that delaying an intervention leads to a larger extension of EOL, although it is at the cost of an increasing risk of failure. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Assessment of Novel Material Railway Drainage Systems" along with RSSB who jointly and equally support it, and also the Lloyd's Register Foundation which partially provides support to this work. RSSB is a rail industry body. Through research, analysis and insight, RSSB supports its members and stakeholders to deliver a safer, more efficient and sustainable rail system. The Lloyd's Register Foundation is a charitable foundation in the UK helping to protect life and property by supporting engineering-related education, public engagement and the application of research. The authors would also like to thank Aursudkij, B., McDowell, G. R., & Collop, A. C. [23] for the experimental dataset.

Appendix A. Modelling details
For granular materials like ballast and suballast under a three-dimensional stress state given by the stress tensor σ ij , = i j , 1, 2, 3, the following relationships are used to obtain the stress invariants p and q: where ϵ ij is the (i, j)-th element of the strain tensor. Under the assumption of axisymmetric stress state ( = The reader is referred to the nomenclature Table A.3 for further information of the geomechanical parameters involved in this work.

Appendix B. Prospective evidence estimation
The prospective evidence  (B.4a)