BIC-Based Data-Driven Rail Track Deterioration Adaptive Piecewise Modeling Framework

The records of maintenance activities are required for modeling the track irregularity deterioration process. However, it is hard to guarantee the completeness and accuracy of the maintenance records. To tackle this problem, an adaptive piecewise modeling framework for the rail track deterioration process driven by historical measurement data from the comprehensive inspection train (referred to as CIT) is proposed. The identification of when maintenance activities occurred is reformulated as a model selection optimization problem based on Bayesian Information Criterion. An efficient solution algorithm utilizing adaptive thresholding and dynamic programming is proposed for solving this optimization problem. This framework’s validity and practicability are illustrated by the measurement data from the CIT inspection of the mileage section of K21 + 184 to K220 + 308 on the Nanchang-Fuzhou railway track from 2014 to 2019. The results indicate that this framework can overcome the disturbance of contaminated measurement data and accurately estimate when maintenance activities were undertaken without any historical maintenance records. What is more, the adaptive piecewise fitting model provided by this framework can describe the irregular deterioration process of corresponding rail track sections.


INTRODUCTION
Track irregularity directly impacts the running stability and safety of trains. Maintaining tracks in an acceptable condition is essential, but it consumes many physical and staff resources. In order to develop cost-effective and rational maintenance plans under limited resources, prior information about track irregularity is required. Thus, this study on predicting the deterioration of track irregularity is critical to railway operation. Many kinds of research have been carried out to forecast track irregularity. Meier-Hirmer et al., (2006) modeled the changes in standard deviation of longitudinal level within a maintenance cycle using the Gamma process. Veit and Marsching, 2010 developed an exponential function to model the behavior of track quality deterioration between two adjacent maintenance events and discussed the interrelations between deterioration rate and the initial quality. Zhu et al., (2013) applied a Gaussian random process to model track irregularities of profile and alignment and studied their power spectral densities. Considering that the evolution of track irregularity is periodic, exponential, and has multiple stages, Xu et al., (2012) employed a multi-stage linear fitting model to describe the track irregularity deterioration process between two adjacent maintenance actions. Lee et al., (2018) combined an artificial neural network (ANN) and support vector regression (SVR) to better represent the deterioration phenomena of track segments for optimizing the maintenance plans in terms of time and cost. In their experiments, at least two years of maintenance records were required to obtain a stable prediction of track deterioration. Mercier et al., (2012) conjointly utilized longitudinal and transversal leveling indicators using a bivariate Gamma process to predict track quality. Vale and Lurdes (2013) developed a stochastic model based on the Dagum distribution for longitudinal level.
Considering that maintenance activities, including tamping, grinding, and others obviously recover track irregularity and have an effect on the deterioration modes (Quiroga and Schnieder, 2010), the aforementioned studies mainly focus on the deterioration process between adjacent maintenance activities. Some studies for multiple maintenance periods have been developed under the following two main assumptions. One is that maintenance records are accessible for modeling; the other is that deterministic mathematical models can express the relationship between deterioration rates and initial qualities right after maintenance actions. Accordingly, segmenting the deterioration process of track irregularity according to maintenance activities is fundamental for exploring the deterioration rules based on historical measurement data. However, the complete and accurate records of maintenance activities are unobtainable because most of the previous records have been lost. Thus, it has become an urgent task to establish an algorithm to automatically identify when maintenance activities were carried out in the process of deterioration (referred to as maintenance-points). Each maintenance-point is tagged by the detection date, which is right after the maintenance activities.
Identification of maintenance-points in the process of track irregularity deterioration is equivalent to making inferences about unknown multiple change-points in the field of applied statistics. There are vast amounts of studies on multiple change-points analysis in different applied contexts, for example, in econometrics (Dias, 2004), in biology (Xi et al., 2011), in climatology (Reeves et al., 2007;Lu et al., 2010), and in hydrology (Perreault et al., 2000). It has also been introduced to traffic flow data for freeway incident detection. Yang et al., (2014) proposed the coupled Bayesian RPCA by extending the Bayesian robust principal component analysis (RPCA) approach for detecting unusual traffic events. The traffic events were localized based on coupling the multiple traffic data streams. Liu et al., (2008) developed an automated traffic incidents detection algorithm on the basis of the cumulative sum (CUSUM). Moreover, in order to achieve real-time defect detection of high-speed train wheels, Wang et al., (2020) utilized the Bayesian dynamic linear model (DLM) to detect change-points in strain monitoring data from high-speed train bogies. Many effective methods have been developed and verified, such as maximum likelihood, Bayes-type, cumulative sum, and others (Jandhyala et al., 2013). Among them, information criteria provides a method for multiple change-points estimation without any priori information on their locations and number (Hall et al., 2013). Bayesian Information Criterion (referred to as BIC) is popularly applied (Watanabe, 2012;Hall et al., 2015). BIC was proposed by Schwarz (1978) and is widely applied as a model selection criterion. Regarding the number of change-points as the dimension of the model, Yao (1988) applied BIC for making inferences about the change-points when the means of observations on different time periods were distinct. However, Zhang and Siegmund (2007) found that the classic BIC had poor performance when applied to irregular statistical models. Thus, Zhang proposed a modified BIC that differently penalized the model dimension components of BIC's objective function. Hannart and Naveau (2012) improved BIC for multiple change-points analysis by introducing priori information on the relative positions and amplitude of change-points and deriving a closed-form mathematical expression of the criterion based on Laplace approximation. Successes in applying BIC to other practical problems such as detecting change in acoustics have been widely reported in the literature (Chen and Gopalakrishnan, 1998;Kotti et al., 2006).
The major contribution of this paper is to propose an adaptive piecewise modeling framework that is driven by historical measurement data from CIT and enables us to describe the rail track deterioration process. This framework is capable of tolerating contaminated measurement data and automatically identifying maintenance-points in the process of deterioration. This problem is reformulated as a model selection optimization problem by taking advantage of BIC. Linear regression (referred to as LR) is applied to model each subsequence individually divided by maintenance-points. Then the objective function is derived according to the framework of BIC and is modified by incorporating an optimized weight for the model complexity component. Based on the effect of maintenance activities on deterioration rate, an efficient solution algorithm for minimizing the objective function is developed by comprehensively utilizing the adaptive thresholding and dynamic programming. The proposed framework is validated by the measurement data for the Nanchang-Fuzhou railway track through CIT collection from 2014 to 2019.
The rest of the article is organized as follows. In Modeling framework based on Bayesian Information Criterion Section, we derive an objective function based on BIC and modify it by incorporating a weight coefficient. In Solution algorithm Section, we develop a solution algorithm based on adaptive thresholding and dynamic programming. Then, we discuss the optimal value of weight coefficient. In Empirical analysis Section, the performance of the proposed framework is evaluated by practical measurement data. Finally, we summarize the research and discuss our ongoing work related to this article.

MODELING FRAMEWORK BASED ON BAYESIAN INFORMATION CRITERION
For Chinese railways, the track quality index (TQI) is employed to quantify track irregularity. It is the sum of standard deviations of seven geometrical parameters for a 200 m-long track section (Xu et al., 2011). The standard deviation for each geometrical parameter is calculated from measurement data collected by CIT.
Among the seven geometrical parameters, track profile irregularity is particularly related to mechanized maintenance activities. Thus, the inference about maintenance-points is studied on the basis of track profile irregularity (referred to as TQI p ).
The inference of change-points based on BIC is a model selection procedure that minimizes a constrained function based on the maximum likelihood method defined by BIC (Gang and Ghosh, 2011). Accordingly, we reformulate the inference on the number and locations of maintenance-points in the deterioration process of TQI p into a model selection problem based on BIC. Denoting the set of all piecewise LR models as Ω and each model in it as Μ ∈ Ω. BIC defines the optimal fitting model from Ω as the one that minimizes Eq. 1.
wherein N is the sample size, L is the maximized likelihood of fitting model M, and K is the number of parameters to be estimated. In a certain time period, suppose that n inspections have been accomplished for a 200 m-long track section, the set of difference in days between the i th detection date and the first detection date are denoted by t (t 1 , t 2 , /, t n ) while the set of corresponding detection values of TQI p by y (y 1 , y 2 , /, y n ). And there are m maintenance-points in y (y 1 , y 2 , /, y n ). The maintenancepoints split y (y 1 , y 2 , /, y n ) into m + 1 independent subsequences. Denoting τ k (0 ≤ k ≤ m + 1) as the maintenancepoint that splits the k th and the (k + 1) st subsequences and τ 0 0, τ m+1 t n . Each subsequence is modeled by LR. Thus, the LR model M k for the k th subsequence is denoted as where β k (β k0 , β k1 ) is the corresponding parameter vector for M k and the random error term ε k is iid. Denoting the variance of ε k as σ 2 k , we obtain that ε k ∼ N(0, σ 2 k ). Based on the assumption on the distribution of ε k , for τ k−1 ≤ t i < τ k , we obtain that y i ∼ N(β k0 + β k1 t i , σ 2 k ). Denoting the value of TQI p in the i th (0 < i ≤ n) detection as y i and the probability density function of y i is expressed as The number of parameters to be estimated, including [(β 10 , β 11 ), /, (β (m+1)0 , β (m+1)1 )], (τ 1 , /, τ m ), and σ 2 is 3m. Based on Eq. 1, we obtain BIC M n(ln 2 π + 1) + ln σ 2 + 3m ln(n) (8) n(ln 2 π + 1) in Eq. 8 is fixed when the series is given. Accordingly, the objective function of BIC( M) is redefined as where ln σ 2 is the sum of squared residuals that reflects the precision of the model and 3m ln(n) is the penalty term of model complexity. We denote ζ(0 < ζ ≤ 1) as a weight coefficient for the complexity of the fitting model. The weight coefficient is determined according to a guideline which will be introduced in The optimal value of the weight coefficient Section. Thus, the object function is The optimal fitting model M for the deterioration process is defined as the one that minimizes Eq. 10. And we consider that the change-points of M are the maintenance-points which will be identified in the deterioration process.

SOLUTION ALGORITHM
Since the number of maintenance-points is unknown, a large amount of computation is needed for attaining the optimal fitting model based on Eq. 10. In order to reduce computation load and to make the algorithm practical, we propose an efficient solution algorithm based on the characteristics of maintenance-points.

The Different Characteristics of Maintenance-Points and Contaminated Measurement Data
This paper is targeted to automatically identify the maintenancepoints in the deterioration process of TQI p for exploring the deterioration rules of track irregularity. However, outliers in the deterioration process caused by contaminated measurement data might interfere with the identification of maintenance-points. Each outlier is tagged by the corresponding detection date. Maintenance-points and outliers are characterized as follows.
The deterioration process of TQI p of a 200 m-long track section on the Nanchang-Fuzhou railway from 2014 to 2019 is shown in Figure 1. As shown in Figure 1, the value of TQI p drops obviously after maintenance activities. Denoting the first order difference of y (y 1 , y 2 , /, y n ) as d (d 1 , d 2 , /, d n ) where d 1 0 and d i y i − y i−1 . d i is much greater/smaller than the neighboring values if maintenance activities were carried out at t i . The outliers caused by contaminated measurement data display the same characteristics. The difference between the maintenance-points and outliers is their different impact on the current deterioration process. The maintenance-points terminate the current deterioration cycle, reduce the value of TQI p to a specified scope, and start a new deterioration cycle. Outliers show significant deviations from the current deterioration process but have no impact on the current deterioration rate.

Candidate Breakpoints Identified by Adaptive Thresholding Method
The maintenance-points and outliers in the deterioration process are collectively referred to as "candidate breakpoints". Distinguishing the maintenance-points from outliers within candidate breakpoints will greatly reduce computation load. Accordingly, we develop a method for identifying candidate breakpoints in the deterioration process based on the aforementioned characteristics of maintenance-points and outliers. Constant thresholding is not feasible since track  irregularity recovers at different degrees after maintenance among track sections. What is more, outliers cannot be within a predetermined range. Adaptive thresholding provides a solution to this problem (Breier and Branišová, 2015;Wang, 2015). On the basis of adaptive thresholding, we develop a method combining the autoregressive model (referred to as AR) to identify candidate breakpoints in the deterioration process. Candidate breakpoints are localized by applying this method to the first order difference d (d 1 , d 2 , /, d n ) of TQI p . The values of (d 1 , d 2 , /, d n ) are dynamically stable within a small range if there is no candidate breakpoint, while the similarity in the distribution of (d 1 , d 2 , /, d n ) is destroyed if there is a candidate breakpoint. Thus, we define a sliding window, and the value at the current moment is predicated based on the historical values which are selected into the sliding window. AR is applied to predicate the value at the current moment based on the historical values. The values of the upper threshold and lower threshold are adjusted via the predicated value and variance of historical values in the sliding window. The difference between predicated value and actual value at the current moment decides whether there is a candidate breakpoint or not. We consider t i as a candidate breakpoint if d i exceeds the preset thresholds. The method is divided into five steps as follows.
Step one: calculate the first order difference d (d 1 , d 2 , /, d n ) of TQI p .
Step two: denote the sliding window as w (w 1 , w 2 , /, w l ) where l is the window size and the absolute value of every element in d as d abs (|d 1 |, |d 2 |, /, |d n |). q 90 is defined as the 90%-quantile of d abs . Starting with the first element in d abs , d i is added into w if |d i | ≤ q 90 while t i is tagged as a candidate breakpoint if |d i | > q 90 , then the remaining elements are recursively checked in sequence until the sliding window is full.
Step three: fit w with AR(p) through the Yule-Walker method (Brockwell et al., 1987), where p is the order of the AR model. Denote the predicated value at the current moment by AR(p) asd i .
Step four: according to the Pauta criterion (3σ criterion), the proportion of outliers in a series is less than 0.3% under the constraint of 3σ (Li et al., 2016). We denote the upper threshold as T upper and the lower threshold as T lower , then where Sd is the standard deviation of historical values in w.
Step five: if d i ∉ [T lower , T upper ], t i is tagged as a candidate breakpoint, and w is not changed. Otherwise, d i is added into w while the earliest element in w is removed. Return to Step three until all of the elements in d have been detected.
The candidate breakpoints identified by the aforementioned method are denoted by τ c (τ 1 , τ 2 , /, τm) wherem is the number of candidate breakpoints. The composite diagram of a typical realization is shown in Figure 2, in which dots represent the first order difference of TQI p , the shaded part represents the limitation range of thresholds, and red crosses represent the identified candidate breakpoints. Moreover, the pseudocode of the adaptive thresholding method is provided in Figure 3.

Dynamic Programming for Finding Optimal Fitting Model
Dynamic programming is a multi-stage optimization method and is applicable to various practical problems (Bellman and Dreyfus, 1962). We now consider a method based on the principle of dynamic programming to find an optimal fitting model that achieves the minimum of Eq. 10. Suppose that r(1 ≤ r ≤m) breakpoints are selected from all the candidate breakpoints τ c (τ 1 , τ 2 , /, τm) and then the series y (y 1 , y 2 , /, y n ) is divided c(τ k−1 , τ k ) is defined as the sum of squared residuals for the subsequence, which is constrained in (τ k−1 , τ k ). We obtain For 0 ≤ j ≤ r, defining b r−j (τ j ) as the minimum of S(τ 1 , τ 2 , /, τ r ) on the basis that the first j breakpoints are confirmed, we obtain b r−j τ j Min τj+1,/,τr Considering the tamping will not be operated for a rail track section twice in a month, we set the constraint that τ j+1 − τ j > 30 and suppose that τ 0 t 0 , τ r+1 t n . Searching the optimal combination of candidate breakpoints is equivalent to solving the following recursive problem: To sum up, the recurrence formulas are From the previous, it is concluded that Min[BIC( M r)] b r (τ 0 ). What is more, Min[BIC( M r + 1)] can be calculated based on the intermediate results for computing b r (τ 0 ). Denoting f r (τ z ) as the optimal results for the subsequence which begins at τ z under the circumstance that the number of selected candidate breakpoints is r, then the recurrence formulas of f r (τ z ) are Searching the optimal results under the assumption that there are r maintenance-points in the deterioration process of TQI p is equivalent to calculating f r (τ 0 ). The iteration is terminated if the aforementioned constraint cannot be satisfied. For 1 ≤ r ≤m, the set of optimal results with a different number of selected candidate breakpoints is denoted by f (τ 0 ) [f 1 (τ 0 ), f 2 (τ 0 ), /fm(τ 0 )] and the optimal piecewise fitting model is the one that achieves Min[f (τ 0 )]. The change-points of this model which are selected from the set of candidate breakpoints, are the maintenance-points, while others are outliers.

The Optimal Value of the Weight Coefficient
The value of the weight coefficient ζ in Eq. 10 has a significant impact on the accuracy and reliability of the results identified by the aforementioned framework. Relying on the historical measurement data from 2014 to 2019 for the nearly 200 kmlong track sections of the Nanchang-Fuzhou rail line (as shown in Figure 4), we obtain the optimal value of ζ which enables the identified maintenance-points to almost correspond with the actual ones. The measurement data are acquired from CIT, which inspects the railways twice a month on average in China.
In particular, the preprocessing and transforming of measurement data, which include mileage correction, historical waveform data alignment, and the TQI calculation of each geometric parameter, have been completed relying on the system developed by our team (Xu et al., 2015). Thus, we consider that the data are complete and reliable. Meanwhile, we obtain a set of actual maintenance-points of each track section via manual analysis. They were considered as correctly identified maintenance-points if they were also included in the set of the actual ones.
To assess the performance of the proposed framework with different values of ζ, we employ the precision (referred to as PRC) and recall rates (referred to as RCL) given by: where T est denotes the number of correctly identified maintenancepoints from candidate breakpoints , F est is the number of erroneous ones, and M man is the number of actual maintenance-points from manual analysis. F 1 defined by Eq. 22 is also considered. The higher the value of F 1 , better performance is obtained.
To find the optimal value of ζ, the maintenance-points of each rail track section are estimated by the proposed framework, whose ζ gradually increases by 0.05. By contrast, using the estimated maintenance-points to the actual ones, we obtain the PRC and RCL of each section. For each ζ, the proportion of sections whose PRC 100% and sections whose RCL 100% are counted separately, while the results are tabulated in Table 1. From Table 1, it is concluded that RCL is getting smaller while PRC is getting larger with the increase of ζ. Based on F 1 , the optimal value of ζ is determined as 0.65.

EMPIRICAL ANALYSIS
In this section, 33 track sections of 200 m in length are further analyzed to demonstrate the performance of the presented framework with ζ 0.65. Then, the calculation procedure of this framework is displayed through two track sections in detail.

Performance Analysis
PRC and RCL for each track section are calculated based on the comparison between the estimated maintenance-points and the actual ones. To evaluate the accuracy of the presented framework with the interference of contaminated measurement data, we have investigated the outliers due to the contaminated measurement data of each track section. N outlier is denoted as the number of outliers. The results of the metrics are tabulated in Table 2. Based on that, we obtain the distributions of PRC and RCL, which are tabulated in Table 3.
From Table 3, we obtain that the proposed framework owns a high PRC and RCL for most rail track sections. It indicates that this framework is capable of overcoming the disturbance of contaminated measurement data and accurately distinguishing the maintenance-points from outliers within candidate breakpoints. Meanwhile, the RCL of a few sections are not satisfactory. Through analyzing further, we find that it has resulted from the fact that not all maintenance-points are included in the set of candidate breakpoints.

Sensitivity Analysis
A large range of thresholds in Candidate Breakpoints Identified by the Adaptive Thresholding Method Section leads to an excessive number of candidate breakpoints, requiring more time to obtain the minimum of Eq. 10. However, the actual maintenance-points might be left out if the range of thresholds is too small. Accordingly, the sensitivity of this framework to the range of thresholds is discussed in this section. According to Eqs 11, 12, we find that the range of thresholds is significantly affected by the times of Sd (the standard deviation of historical values selected into the sliding window), which is denoted by F Sd . To access the sensitivity, we apply this framework to some of the sections in Table 2 for each F Sd ∈ {2.5, 3, 4}. The other values of parameters in the implementation of this framework are kept consistent with those in Performance Analysis Section. The comparison results among different values of F Sd are tabulated in Table 4. The computation time in seconds and the number of candidate breakpoints are given in column 3, 4, while the estimated maintenance-points of each section with different values of F Sd are given in column 5. From Table 4, we obtain that compared with F Sd 3, the computation time is increased by 145% on average for F Sd 2.5. What is more, the estimated maintenance-points are the same between F Sd 3 and F Sd 2.5. Although it requires less computation time when F Sd 4, some of the actual maintenance-points are left out as the estimation results of the sections starting at K51.367, K84.601, and K112.791 indicate. Thus, it is reasonable that F Sd 3 in Candidate Breakpoints Identified by the Adaptive Thresholding Method Section. Section One: K117 + 137-K117 + 337 This section is on a tangent track. The candidate breakpoints identified by the adaptive thresholding method are shown in Figure 5A. The value of Min[BIC( M r)] for a different number of selected candidate breakpoints are shown in Figure 5B and BIC( M) obtains the minimum when the number is 2. The identified maintenancepoints are 2016-05-26 and 2018-04-13. The estimated maintenance-points and piecewise fitting model are shown in Figure 5C. The maintenance-points estimated by the proposed framework are exactly as the actual ones. The deterioration process is divided into three subprocesses. Although there are lots of outliers caused by contaminated measurement data in the first subprocess, the maintenance-points are accurately identified.
Section Two: K162 + 900-K163 + 100 This section is on a curved track. The candidate breakpoints identified by the adaptive thresholding method are shown in Figure 6A, and the value of Min[BIC( M r)] for this section are shown in Figure 6B. The number of selected candidate breakpoints corresponding to the optimal result is 3. The estimated maintenance-points and piecewise fitting model are shown in Figure 6C. The identified maintenance-points are 2016-3-14, 2018-01-12, and 2018-06-26. The outliers caused by contaminated measurement data are mainly located in the first subprocess. The slope of the fitting model for each subprocess reflects its deterioration rate. It is obvious that the deterioration rates are different before and after maintenance activities. Thus, we believe that the deterioration rates might be affected by maintenance activities.

CONCLUSION
In this paper, a rail track deterioration modeling framework driven by historical measurement data from CIT is proposed.
The modeling framework requires no historical maintenance records and does not assume the quality of track measurement data. The proposed framework formulates the identification of maintenance activities with a model selection optimization problem, based on a modified Bayesian Information Criterion by incorporating an optimized weight for the model complexity component into the objective function. An efficient solution algorithm utilizing adaptive thresholding and dynamic programming is proposed for the model selection problem, taking the characteristics of the effect of maintenance on track deterioration trend.
The proposed track deterioration modeling framework is applied to the historical measurement data from 2014 to 2019 for the nearly 200 km-long track sections of the Nanchang-Fuzhou rail line. Based on that application, the optimal value of the weight coefficient which is incorporated for the model complexity is discussed in The Optimal Value of the Weight Coefficient Section. Moreover, the assessment indicators are calculated based on 33 200 m-long track sections. As the assessment indicators indicate, the proposed framework is capable of accurately identifying the maintenance-points and creating an adaptive piecewise model of the deterioration process.
However, for a few track sections, the estimated maintenancepoints are less than the actual ones, which resulted from the fact that not all maintenance-points were included in the set of candidate breakpoints. Therefore, one of the emphases for the next step will be on improving the algorithm to ensure that the set of candidate breakpoints contains all of the maintenancepoints.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation