# Side Wear Prediction of a Subway Outer Rail on Small Radius Curves Based on System Dynamics of Discrete Supported Track

Authors: Xianxian Yin; Linlin Kou; Xiukun Wei; Haichao Zheng
Journal: Discrete Dynamics in Nature and Society
 (2022)

---

## Abstract

When a subway train runs on a small radius curve, it causes serious side wear of the outer rail, and a large number of the outer rails of the curve are replaced frequently due to excessive side wear. This paper proposes three prediction models based on the response surface methodology (RSM), support vector machine (SVM), and the relative vector machine (RVM) to estimate the remaining useful life (RUL) of a metro outer rail at a small radius curve determined by rail side wear. In these proposed models, the directly measurable parameters, including various track geometries, axle loads, and primary suspension stiffnesses related to vehicle types, are taken into account. These parameters are proposed from the side wear mechanism and subjected to the method of global sensitivity analysis to rank them in order of importance. The rail side wear for various scenarios is calculated through a numerical simulation model of the vehicle-unballasted track and Specht friction wear. Furthermore, with the help of simulation experiments and field cases, the prediction performances of three different prediction models are analyzed in detail. The comparison indicates that the predictive model based on RVM is superior in the prediction of RUL based on rail wear. The findings presented in this paper can provide certain reference values for infrastructure managers (IMs).


---

## Body

## 1. Introduction

Rail wear at small radius curves (SRCs) has received much more attention than ever before from subway operation companies. The rail wear changes the profile of the rail head, resulting in different contact states of the wheel and rail, which has a great impact on the stability and safety of the train and the service life of each component of the vehicle-track system. A large number of outer rails at SRCs are replaced frequently due to excessive side wear. Reducing the side wear of outer rails is an urgent problem for subway operating companies that needs to be solved. Therefore, it is very important to investigate the wear mechanism and predict the remaining useful life (RUL) of railway tracks determined by rail side wear. Based on this, effective means can be proposed to reduce the wear at SCRs, which can provide certain reference values for infrastructure managers (IMs).In the last decades, some research has been done in the field of side rail wear [1–7], including numerical simulation technology of wheel-rail wear, development law and prediction method of the SRW (SRW refers to the side rail wear of outer rail in the following.) at SRCs, and also preventive maintenance strategy of rail wear. In [8], the factors influencing SRW at SRCs are obtained based on the analysis of on-the-spot track monitoring test and years of experience in rail maintenance. In [9], the influence of gauge on the side wear of outer rail with radius smaller than 650 m is studied. Based on the numerical simulation model of wheel-rail contact wear, a prediction model of wheel-rail wear is proposed in [10]. The results show that the wear of the outer rail is aggravated with the increase of the superelevation. In [11], the rail wear is divided into three stages on the basis of field investigation and analysis. In [12], the monthly side wear rate of the curve rail is calculated, in which only the influence of vehicle type to the SRW is considered. In [13], the side grinding rate curve of the curved outer rail is obtained based on the rail side wear measurements of 25 ordinary speed railway curves. Zhou Yu [14], Yan Yizhu [15], and Pan Jianjie [16], in order to study the development law of side wear of curved outer rail, tracked and measured the side wear of large volume and small radius curved rail in Shanghai, Guangzhou and Beijing Metro respectively. In [17], the wear law of high-speed railway rails is analyzed based on the Braghin wear model. In [18], the Zobory and the Braghin wear model are both introduced to study the influence of friction coefficient on wheel wear of high-speed railway. It found that when the friction coefficient is greater than 0.2, the calculated Braghin wear of high-speed railway wheels is basically unchanged. In [19], the side wear of the small radius curve of the EMU is studied on the basis of the Heuman wear index.There are also some reports on wheel-rail wear reduction. In [20], it shows that wheel-rail wear can be reduced by controlling the wheel-rail friction. The relationship between wheel-rail contact and vehicle dynamic performance is analyzed in [21]. It is proposed that the rail wear can be reduced by controlling the taper of rail tread. In [22], a design method for rail wear profiles for heavy-duty railway curves is proposed, theoretical and experimental results show that the use of this wear profile reduces SRW by 30%–40%. In [23], based on the on-the-spot investigation of the SRW on SRCs of Iranian railway, some measures to reduce rail wear are put forward, such as increasing rail hardness and improving lubrication between wheel and rail. In [24], the influence of the curve radius on the wear characteristics of rails is studied with the help of numerical calculation method and simulation test. In [25–27], it found that on the premise of ensuring driving safety, moderate underelevation is beneficial to reducing wheel and rail wear.As for wheel-rail wear prediction, a method of wheel profile estimation based on iterative wear and surface updating in the wear prediction model of friction materials is discussed in [28]. In the light of the test results of rail wear on JD-1 wheel-rail simulator, nonlinear models of running speed, attacking angle, axial load and rail wear are built up in [29]. In [30], Genetic algorithm (GA) is used to optimize the control parameters of BP neural network to predict the rail wear. In [31], the rail wear is predicted with the help of non-Hertz wheel-rail contact and energy dissipation model.To summarize [8–10], the influence of a single factor, such as the curve radius, superelevation, rail cant and gauge is analyzed based on rail maintenance experience or rail wear field tracking test results. In [12], only the rail wear laws corresponding to different types of vehicles are studied, and the quantitative prediction of SRW on different curves is not considered. Furthermore, the key causes of SRW are not analyzed by using the theory of wheel-rail contact wear on SRCs, and these parameters cannot be measured directly in practical applications. Moreover, the interaction and coupling of some factors on rail wear as well as the influence of different wear regimes on rail wear evolution are ignored, which is not comprehensive and accurate for the analysis of the factors causing rail wear. Therefore, the core work of this paper is to establish regression models for indirectly predicting the RUL of railway tracks determined by rail wear, where the wear evolution is taken into account and the input parameters are expected to be measurable in practical applications.The main contribution of this paper lies in the theoretical analysis of the key causes of the rail side wear by means of the wheel-rail contact dynamics, as well as the quantitative analysis and prediction on the development of the outer rail wear.In this paper, the RUL prediction of railway tracks at SRCs determined by side rail wear is presented by using a theoretical analysis of wheel-rail contact dynamics, a simulation model of vehicle-track-wear, and an on-site analysis of SRW. The domain factors, including various track geometry, axle load, and the primary suspension stiffness related to vehicle type are incorporated to estimate the rail side wear that can be translated to a RUL prediction of railway tracks. These parameters are first proposed from the derivation of the side wear mechanism and are ranked in order of importance with the help of the global sensitivity analysis. This paper chooses the top seven parameters for further analysis. The relationship between these parameters and rail side wear is established by using three different multivariate regression models, where different wear regimes in rail wear evolution are taken into account. Furthermore, the predicted performances of separate prediction models are analyzed for contrast. Finally, in addition to the simulation validation, validation using field measurements of the Guangzhou Metro Line 1 is also achieved.This paper is organized as follows. Section2 describes the rail side wear mechanism and the methodology of rail wear simulation, the sensitivity analysis and the concept of prediction models. Section 3 presents the results of the numerical investigations and discusses the prediction performance of different prediction models. In Section 4, validation using the field measurements from the Guangzhou Metro Line 1 is performed. Finally, conclusions are given in Section 5.


## 2. Rail Side Wear Mechanism and Methods

This section describes the rail side wear mechanism and the concepts used to derive the predictive models to estimate the RUL of railway tracks. Section2.1 presents the rail side wear mechanism and the wear calculation model. Section 2.2 focuses on the numerical simulation model for wear calculation, Section 2.3 focuses on the sensitivity analysis used to obtain the most dominant parameters for rail wear, and Section 2.4 discusses the theory of the predictive models.


### 2.1. Rail Side Wear Mechanism

When a train passes through curves, the centrifugal force on the guide wheel of the outer rail cannot be completely balanced by the lateral creep force between the wheel and rail due to many factors. The contact of the flange and the rail leads to SRW that is located at a position of 16 mm below the rail surface [11]. Figure 1 shows the wheel-rail contact forces of the leading wheelset in the pass of the curve. Here, only the wheel-rail contact forces of the leading wheelset of each bogie is introduced an example to extract the key causes of the rail side wear in the curved track from the perspective of wheel-rail contact theory. All the symbols used in this section are listed in Table S1 of the Appendix.Figure 1
Wheel-rail contact force.The guiding forceFwf provided by the outer rail flange on SRCs is as follows.(1)Fwf=ΔFcf−Fcy,where ΔFcf is the unbalanced centrifugal force, and Fcy is the lateral creep force.Considering the influence of the outer rail superelevation and the rail bottom slope, the unbalanced centrifugal forceΔFcf acting on the guide wheel of high rail is deduced as follows [32]:(2)ΔFcf=mw⋅ggt.kh⋅v2R−Δh⋅1+sL212⋅1+s212−v2R,where the coefficient kh is the function of g (in SI unit), international standard gauge (unit in mm), rail head width (unit in mm), and the ratio 3.6 when changing the unit from km/h to m/s. Other symbols are listed in Table S1 of the Appendix.As the wheel-rail creep force is saturated at small rail curves, the lateral creep forceFcy acting on wheelset is equal to the product of friction coefficient and normal force between wheel and rail [32], hence(3)Fcy=μFrw,where Frw is the wheel-rail normal contact force and μ is the friction coefficient.The unbalanced superelevation on the curves with small radius leads to the phenomenon of eccentric load on the inner and outer rails, which directly affects the creep force under the wheels on the inner and outer rails [32]. Regardless of the vertical vibration of the vehicle and the rolling vibration of each part of the vehicle that may be excited by the lateral irregularity of the track, the unbalanced load of the outer rail can be considered as the steady-state constant value. It is stated as follows [32, 33].(4)Frw=1G⋅zwcosδR+ϕw32+14mc+12mb+mwghcgt2ghgt−v2R,where G is wheel-rail contact constant and zw is the vertical displacement of wheelset. mc, mb, and mw are the vibration mass of car-body, bogie and wheelset, respectively. hc is the height from vehicle centre of gravity to rail surface.In light of equations (3)-(4), it yields(5)Fcy=μ1G⋅zwcosδR+ϕw32+14mc+12mb+mwghcgt2ghgt−v2R.Substitute equations (2), (3), (5) into (1), we obtain(6)Fwf=mwg⋅1gv2R+1gt.kh⋅v2R−Δh⋅1+sL212⋅1+s212−μ1G⋅zwcosδR+ϕw32+14mc+12mb+mwghcgt2ghgt−v2R.When the wheelset passes through rail curves, it is guided by the flange of the wheel on the outer rail. In this case, there are two or more points of contact between the guide wheel and the gauge corner of the outer rail, which results in rail wear. The Specht wear model is used to calculate the wear between wheel and rail for this case by applying different wear parameters for moderate wear and severe wear, respectively [34]. Specht wear model can be considered as a deformation of Archard model. The difference is that Specht wear model considers the multipoint contact between wheel and rail. This model introduces a jump factor on the basis of Archard model to characterize the multipoint contact condition between the wheel and the rail. When there is a multipoint surface contact between the wheel and rail, it indicates that the wheel-rail wear is serious, and the jump factor takes the empirical value of 10. Otherwise, it can be considered as a single-point contact between the wheel and rail, with slight wear, and the jump factor takes the empirical value of 1. This model assumes that the friction work generated by wheel-rail contact is linear with the volume of wear. Then, the wear volume V can be derived as follows.(7)V=kV⋅km⋅W,where W is the friction work on rail due to wheel-rail contact, and kV is a dimensionless wear coefficient related to wheel-rail material. km is the jump coefficient, and the empirical value is 1 for mild wear and 10 for severe wear [35].The wear depth is represented by the average wear depth regardless of the difference of wear depth within a contact patch, it can be derived as follows:(8)d=kV⋅km⋅FwfDHAc,where d is wear depth. D is the sliding distance. H is the hardness index, and Ac is the size of the contact area.By substituting equation (6) into (8), the wear model of outer rail is obtained as follows:(9)d=kV⋅kmH⋅DAcmwg⋅1gv2R+1gt.kh⋅v2R−Δh⋅1+sL212⋅1+s212−μ1G⋅zwcosδR+ϕw32+14mc+12mb+mwghcgt2ghgt−v2R.It is noted from equation (9) that the rail side wear depends on many parameters, such as the curve radius (R), gauge (gt), running speed (v), and load (mc, mb, mw). These parameters can be clustered into four categories. The first category includes track parameters, which consist of the curve radius (R), gauge (gt), track longitudinal gradient (sL), rail cant (s), and actual superelevation (h). The second category is wheel-rail contact parameters (Ac, zw, δR, ϕw, D), which are related to the wheel-rail profiles as well as the vehicle operation status. In addition, the vehicle parameters are considered the third category, which consists of running speed (v), axle load (mc, mb, mw), and the vehicle type characterized by the stiffness and damping of the suspension system. Finally, the fourth category includes the maintenance activities, where the grinding, lubrication, and weather conditions are all captured by the coefficient of friction (μ).As a number of parameters are required to describe the wear profiles of the wheel and rail, it is not convenient to measure directly in practical applications. In addition, the core objective of this paper is to predict rail wear with the help of directly measurable parameters. Therefore, the influence of the wheel-rail profile is not considered, while the other twelve parameters obtained from equation (9) are further explored in the sensitivity analysis presented in section 2.3. The simulation model used for rail wear calculation is introduced in the next section. With the help of the dynamic software Universal Mechanism (UM) [36], the rail wear produced by vehicles negotiating different curves is calculated for further analysis.


### 2.2. Simulation Environment

In this section, a subway vehicle-track dynamic model is built by using UM for rail wear calculation by using the vehicle and track parameters of Guangzhou metro as an example. The UM software includes different models of railway rolling stock including simplified models and 3D models, and other different software tools for the simulation of railway rolling stock dynamics are also included in this multi-body dynamics software. The vehicle (including the car-body, primary and secondary suspension systems, bogie system, etc.), track, and wheel-rail contact system in the vehicle-track coupling system in terms of “Universal mechanism” are all a subsystem, which can also be models of any complexity.The subway vehicle-track dynamic model is built referring to the parameters of Guangzhou metro vehicle and track, and the simulation calculation of rail superposition wear is carried out by combining the wheel-rail non-Hertz contact coupling dynamic model and Specht material wear model. The feasibility of the simulation model is verified by the on-site tracking monitoring data. Based on the structure of Guangzhou Metro A vehicle, a vehicle-ballasted track coupling dynamic model is established by using UM multi-body dynamics software. The vehicle model is a nonlinear model, including the car body, two bogies, four wheelsets and suspension devices, in which the car body and the bogie are regarded as rigid bodies and a total of 68 degrees of freedom are generated. The rail is considered as Euler beams supported by discrete elastic points distributed according to sleeper spacing. Detailed modelling parameters are shown in TableS2 in the Appendix. The total length of the line is 2 km. Consider that the track random irregularity has a significant influence on the nonuniform distribution rail wear along the rolling direction [37, 38], the track irregularity data of Guangzhou Metro Line 1 is used in the simulation study.The Specht wear model [34] is used to calculate the rail wear at the midpoint of the circular curve. The normal load of wheel-rail contact is calculated by multipoint contact and the normal parameters of wheel-rail contact is obtained by Kik–Piotrowski program [39]. The FASTSIM algorithm is used to calculate the tangent contact parameters of the wheel-rail contact system in dynamic calculations. According to the assumption of linear relationship between friction work and wear depth in Specht wear model, the wear depth of each point can be obtained.The calculation of rail superimposed wear is proceeded as follows. The normal wear threshold of rail is set todT, and the rail profile in the simulated model is updated for the next iteration once the maximum normal wear dimax of a certain point on the rail contact spot reaches the threshold. This paper assumes that the rail profile does not change in the same iterative calculation, that is, the rail profile when the first wheelset passes is consistent with that when the n th wheelset passes. Considering the wheel-rail contact state when the i th wheelset passing, the wear depth after the i th wheelset passes as Δi, and the cumulative wear depth Δ (Δ=∑i=1nΔi) after a train passes can be obtained. When the SRW reaches the requirement of rail replacement (Guangzhou Metro requires 13 mm), the iteration calculation is completed. The wear data and the number of cumulative vehicles passing by (NVP) are the outputs. In this simulation, the empirical value of 0.1 mm is taken as the wear threshold. The simulation model is adjusted based on the analysis results of wear measurements of Guangzhou Metro Line 1. The material wear coefficient kV calculated based on the field measurement results of side wear is 2.103e−13 m3/J.To validate if the model can be used for the further analysis, the monitoring rail wear of Yangji-Sports Centre interval of Guangzhou Metro Line 1 and the simulation results are introduced for comparative analysis. Figure2 shows the comparison of the on-the-spot measurement of SRW and the simulation results, and the Y-axis represents the SRW measured in the position of 16 mm below the rail top. The X-axis represents the number of cumulative vehicles passing by (NVP). It should be noted that 6-car formation is adopted for the train in this paper, and 4 wheelsets in each vehicle is used to calculate the rail wear in the outer curved track. Therefore, when a train passes by a curved track, a total of 24 wheelsets pass through a certain position of the outer rail in turn. In Figure 2, the radius of the measured curve is 500 m, and the lines are all underground, where the influence of environmental factors on rail wear can be captured by the friction coefficient (μ=0.4). See curve C6 in Table S3 in the Appendix for details of the line parameters and measuring points. The figure shows the monitoring data of the side rail wear of the curve from January 2014 to October 2018. During this period, no rail grinding or replacement has been performed. It shows that the measured wear has a linear growth trend with the cumulative increase in the vehicles passing by, and the numerical simulation results are in good agreement with the measurements, which verifies the feasibility of the simulation model for rail wear calculation.Figure 2
Measured and simulated wear curve.


### 2.3. Sensitivity Analysis

In light of the theoretical analysis of SRW, the influence of these parameters on wear, including the radius of the curve, axle load, friction coefficient, unbalanced superelevation, rail cant, track gradient, running speed, and stiffness of primary and secondary suspensions, cannot be ignored. To obtain the most crucial influential factors, this section introduces the global sensitivity analysis based on Sobol' method [40, 41], which is widely used at present. It studies the global influence of different parameters on the model. The range of the parameters can be extended to the whole definition domain, and different parameters can change simultaneously, reflecting the interaction between the parameters.The central idea of Sobol' method is to decompose the functionfx into 2n increasing terms, as follows:(10)fx=f0+∑i=1nfixi+∑1≤i<j≤nfi,jxi,xj+⋯+f1,2,…,nx1,x2,…,xn,where f0 is a constant term, and the integral of the other subterms to any factor is zero.The decomposition of the function is unique, and each subterm can be expressed as an integral of the functionfx. The left and right sides of equation (10) in the entire parameter domain In are squared and integrated, where In=x|0≤yi≤1;i=1,2,…,n, and(11)∫Inf2xdX−f02=∑s=1n∑i1<i2<...<isn∫fi1,i2,…,is2dxi1dxi2…dxis.The total variance is as follows:(12)D=∫f2xdX−f02.The partial variance of each order is(13)Di1,i2,…,is=∫01⋯∫01fi1,i2,…,is2xi1,xi2,…,xisdxi1dxi2…dxis.The ratio of each order's partial variance to the total variance is each order's sensitivity, and theS-order sensitivity is as follows:(14)Si1,i2,…,is=Di1,i2,…,isD.For the rail wear studied in this paper,Si is the first-order sensitivity, which shows the influence of the parameter xi on the rail wear. Si,ji≠j is a second-order sensitivity value, which shows the influence of the parameters xi and xj on the rail wear. At the same time, the coupling interaction between these two parameters is also considered. In other words, the total sensitivity not only reflects the influence of this parameter on the rail wear but also analyses the interaction with all other parameters. The larger the first-order sensitivity value of a parameter, the more significant the influence of this parameter is on the rail wear. When the difference between the first-order sensitivity value of a parameter and its total sensitivity value is large, it can be assumed that there is a coupling interaction between the parameter and other parameters [42].Before performing a global sensitivity analysis on rail wear under random parameters, the change interval and probability distribution of the parameters need to be given, and the parameters are randomly sampled within the interval. This paper builds a Sobol' sequence for random sampling calculation.The Sobol' sequence is obtained by constructing a group called the 'direction number' [43]. Assuming that mi is a positive odd number less than 2n, then(15)ai=mi2i.The direction numberai is generated by a simple polynomial fx with coefficients of only 0 and 1, as follows:(16)fx=xp+b1xp−1+⋯+bp−1x+bp.Fori>p, the recursive formula is as follows:(17)bi=b1ai−1⊕b2ai−2⊕⋯⊕bpai−p⊕ai−p2p,where ⊕ is the binary XOR. Combining the above equation with (15), we can obtain(18)mi=2b1mi−1⊕22b2mi−2⊕⋯⊕2pbpmi−p⊕mi−p.Thei th number of the Sobol' sequence is obtained as follows:(19)xi=c1a1⊕c2a2⊕c3a3⊕⋯,where c1c2c3… is a binary representation.The random number generated by the Sobol' sequence is more uniform, and there is no large blank or repetition phenomenon. Therefore, the distribution of random numbers can be described completely with fewer sampling points, which effectively improves the accuracy and efficiency of sensitivity analysis. Therefore, the Sobol' sequence method is used to randomly sample the parameters in this paper.


### 2.4. Prediction Models

Regression models can be used to find mathematical expressions of relatively complex physical processes and systems. This expression or calculation model can be regarded as a transfer function between a given input and an output [44]. Different types of regression models are available [6], such as the response surface methodology (RSM), support vector machine (SVM), relative vector machine (RVM), and artificial neural network (ANN). This study applies three different regression analysis methods, namely, RSM, SVM, and RVM, that fit our train rail wear data to predict the RUL of outer rails based on rail wear. The data points used for the regression fitting or training process are selected by the central combination design (CCD) [45, 46].


#### 2.4.1. RSM Regression Analysis

RSM can solve problems where the relationship between the objective function and the independent variable is unknown, and it is used to model and analyze the problem in which the objective function is affected by multiple variables in order to optimize the objective function [47, 48].As already mentioned in the introduction, the purpose of this study is to calculate the RUL of the rail based on wear by using a regression prediction model. Hence, for the response value of this study, the service life of the rail based on wear is selected. Here, the accumulative number of vehicles passing by (NVP) the line during the rail service is taken as the service life of the rail. Then, the rail wear undern scenarios obtained by the CCD are calculated according to the calculation process of rail wear described in Section 2.1, and the response values corresponding to that are subjected to a regression analysis to obtain a predictive model.The prediction model obtained by means of RSM has a polynomial form with a certain order. In this study, a quadratic polynomial model is chosen for fitting and is given as follows:(20)yx=a0+∑i=1naixi+aiixi2+∑j=2n∑i<jjaijxixj+ε,where a0, ai, aii, and aij are fitting coefficients and $n$ indicates the number of analysis variables.Due to the complexity of RSM polynomial fitting, this study only considers the second-order coupling cross between the influencing factors. A more complex coupling effect of various parameters on wear cannot be taken into account, which will inevitably affect the prediction performance of the RSM model. For comparative analysis, the machine learning algorithm of SVM is introduced in the next section.


#### 2.4.2. SVM Regression Analysis

SVM is a classification and regression analysis tool proposed by Vapnik et al. in 1995 [49] and currently is mainly used in pattern recognition and regression analysis. The greatest advantage of SVM lies in the solving of problems with a small sample and nonlinear and high-dimensional pattern recognition, which effectively overcomes the problems of large deviation, over learning, local extremum and large amounts of calculation. It can also be extended to other machine learning problems, such as function fitting [50].Given the training samplesX=x1,y1,x2,y2,…,xp,yp⊆Rn×R, where xi is the input variable, yi is the corresponding output variable, and p is the total number of sample points, according to the principle of support vector regression, the linear regression function is as follows:(21)fx=ωTϕx+b,where ω∈Rn is the weight vector of the optimal hyperplane, b is the offset threshold, ϕx is a nonlinear mapping function, and ωTϕx represents the point multiplication of two vectors.To determine the regression function coefficientsω and b, the above problem can be transformed into a constrained optimization problem as follows:(22)maxωai,ai∗=−12∑i,j=1pai−ai∗aj−aj∗Kxi,xj−ε∑i,j=1pai+ai∗+∑i=1pyiai−ai∗,∑i=1pai−∑i=1pai∗=0,0≤ai,ai∗≤C,i=1,2,…,p,where Kxi,xj=ϕxiT⋅ϕxj, Kxi,xj is the kernel function.By solving equation (22), the regression function of SVM can be obtained as(23)fx=∑i,j=1pai−ai∗Kxi,xj+b.Although SVM has been widely used in the field of prediction, its prediction results lack the ability to express uncertainty. Therefore, the RVM model with uncertainty expression is introduced in the next section for regression prediction analysis.


#### 2.4.3. RVM Regression Analysis

RVM is a classification and regression method proposed by Tipping in 2000 [51] and is one of the important research hotspots in the field of statistical learning in recent years. It has a similar function form as SVM and solves the nonlinear classification problem by introducing a kernel function and dimensionality increase.Given a data setG=xi,tii=1N, xi∈Rd represents the input vector, ti∈R represents the corresponding output target vector, and N is the number of input samples. The training model can be expressed as:(24)t=yx+ε,where y⋅ is a nonlinear function and ε is Gaussian white noise with variance σ2, that is, ε∼N0,σ2. The approximation function y can be obtained by a regression operation using the training data. The mathematical expression of RVM can be expressed as(25)t=Φω+εn,where ω=ω1,ω2,…,ωNT, ω is an N+1 dimensional vector representing the weight of RVM. Φ is the kernel function matrix, and Φ=ϕ1x,ϕ2x,…,ϕNxT, ϕix=1,Kxi,x1,…,Kxi,xN, i=1,2,…,N, and K⋅ is a kernel function.The Bayesian inference based on the Gaussian process is applied to kernel theory, and the uncorrelated sample points are removed by an autocorrelation decision theory under the structure of prior parameters. Hence, the RVM sparse probability model is obtained. When a new test sample pointxN+1 is input into the RVM, its prediction model can be expressed as(26)ptN+1|t=∫ptN+1|ω,a,σ2pω,a,σ2|tdωdadσ2,where tN+1 represents the observation target of the new test sample point xN+1.For a given new input valuex∗, the corresponding predicted output value t∗ satisfies the following Gaussian distribution:(27)pt∗|ω=∫pt∗|ω,aMP,σMP2pω|t,aMP,a,σMP2dω∼Ny∗,σ∗2,where y∗ and σ∗2 are expressed as follows:(28)y∗=μTΦX∗,(29)σ∗2=σMP2+ΦX∗T∑ΦX∗,where y∗ is the predicted output value of the RVM model with the new test data x∗, and the first term of equation (29) represents the variance of the noise estimation, and the second term represents the uncertainty of the weight estimation.To compare the prediction errors of the above prediction models, four evaluation indexes of the prediction performance are used in this paper, namely, the coefficient of determination (R2), root mean square error (RMSE), mean absolute error (MAE), mean relative error (MRE), and normalized mean square error (NMSE). The calculation formulas are as follows:(30)R2=1−∑y−y^2∑y−y¯2,RMSE=1N∑i=1Nyi−y^i2,MAE=1N∑i=1Nyi−y^i,MRE=1N∑i=1Nyi−y^iyi×100%,where yi and y^i represent the real value and predicted value, respectively, N is the total number of test samples.The closer the value ofR2 is to 1, and the closer the values of the other three indicators are to 0, indicating that the closer the prediction result is to the true value, the better the prediction performance of the regression model.


#### 2.4.4. Determination of Kernel Function and Hyperparameters

SVM and RVM models with kernel function of Gaussian Radial Basis Function (RBF) have good nonlinear processing ability. When dealing with practical engineering data, it is better to choose RBF kernel function than other kernel functions in the absence of prior knowledge [52].The hyperparameter of the regression model has a very significant impact on the accuracy of the prediction results while SVM and RVM regression models are used to predict the development trend of rail wear. In this paper, thek-fold cross validation procedure is introduced to select the penalty factor C and the parameter σ of RBF kernel function [52], so as to make the prediction result reach the optimal value. In this method, the original data is divided into k-fold, and each subset data is used as a verification set, and the remaining k−1 subset data is used as the training set, so K models can be obtained. The K models were evaluated in the validation set, and the final MSE (mean squared error) was added and averaged to get the cross validation error. The k-fold cross validation is suitable for the limited simulation data in this paper, and the evaluation results can be as close as possible to the performance of the model on the test set, so it can be used as an index for model optimization. In this paper, 10-fold cross validation is selected to optimize the model parameters.The steps to optimize model parameters are as follows:(1)
Determine the initial search interval and step size of the parameter pair (C, σ) as C∈a,b, σ∈c,d. The search steps of parameters C and σ are L0 and L1 respectively, and the k-fold cross-validation procedure is used to calculate the MSE for the selected parameters to evaluate the generalization ability of the model. According to the global minimum criterion of MSE, a preliminary optimal parameter pair is obtained as (C0, σ0). In this paper, the exponential form with the bottom of 2 is selected to perform the initial division of the parameter grid for fast search. For example, L0=1, L1=0.5, and C=2−a,2−a+1,…,2b−1,2b, σ=2−c,2−c+0.5,…,2d−0.5,2d.(2)
Expand the preliminary optimal parameter pair (C0, σ0) to obtain the quadratic grid search intervals of parameters C and σ as m,n, p,q. In the new grid, a further fine search is performed with step lengths L0′ and L1′, and the step length of the fine search satisfies: L0'<L0,L1′<L1. Using the criterion of MSE to take the global minimum, the optimal parameters (C∗,σ∗) can be obtained after the second fine grid search.The above steps are based on SVM model. In the RVM regression model, the penalty factorC is assigned automatically, so only the parameter optimization of the kernel function parameter σ is needed. The method and process are similar to it and will not be repeated.


## 2.1. Rail Side Wear Mechanism

When a train passes through curves, the centrifugal force on the guide wheel of the outer rail cannot be completely balanced by the lateral creep force between the wheel and rail due to many factors. The contact of the flange and the rail leads to SRW that is located at a position of 16 mm below the rail surface [11]. Figure 1 shows the wheel-rail contact forces of the leading wheelset in the pass of the curve. Here, only the wheel-rail contact forces of the leading wheelset of each bogie is introduced an example to extract the key causes of the rail side wear in the curved track from the perspective of wheel-rail contact theory. All the symbols used in this section are listed in Table S1 of the Appendix.Figure 1
Wheel-rail contact force.The guiding forceFwf provided by the outer rail flange on SRCs is as follows.(1)Fwf=ΔFcf−Fcy,where ΔFcf is the unbalanced centrifugal force, and Fcy is the lateral creep force.Considering the influence of the outer rail superelevation and the rail bottom slope, the unbalanced centrifugal forceΔFcf acting on the guide wheel of high rail is deduced as follows [32]:(2)ΔFcf=mw⋅ggt.kh⋅v2R−Δh⋅1+sL212⋅1+s212−v2R,where the coefficient kh is the function of g (in SI unit), international standard gauge (unit in mm), rail head width (unit in mm), and the ratio 3.6 when changing the unit from km/h to m/s. Other symbols are listed in Table S1 of the Appendix.As the wheel-rail creep force is saturated at small rail curves, the lateral creep forceFcy acting on wheelset is equal to the product of friction coefficient and normal force between wheel and rail [32], hence(3)Fcy=μFrw,where Frw is the wheel-rail normal contact force and μ is the friction coefficient.The unbalanced superelevation on the curves with small radius leads to the phenomenon of eccentric load on the inner and outer rails, which directly affects the creep force under the wheels on the inner and outer rails [32]. Regardless of the vertical vibration of the vehicle and the rolling vibration of each part of the vehicle that may be excited by the lateral irregularity of the track, the unbalanced load of the outer rail can be considered as the steady-state constant value. It is stated as follows [32, 33].(4)Frw=1G⋅zwcosδR+ϕw32+14mc+12mb+mwghcgt2ghgt−v2R,where G is wheel-rail contact constant and zw is the vertical displacement of wheelset. mc, mb, and mw are the vibration mass of car-body, bogie and wheelset, respectively. hc is the height from vehicle centre of gravity to rail surface.In light of equations (3)-(4), it yields(5)Fcy=μ1G⋅zwcosδR+ϕw32+14mc+12mb+mwghcgt2ghgt−v2R.Substitute equations (2), (3), (5) into (1), we obtain(6)Fwf=mwg⋅1gv2R+1gt.kh⋅v2R−Δh⋅1+sL212⋅1+s212−μ1G⋅zwcosδR+ϕw32+14mc+12mb+mwghcgt2ghgt−v2R.When the wheelset passes through rail curves, it is guided by the flange of the wheel on the outer rail. In this case, there are two or more points of contact between the guide wheel and the gauge corner of the outer rail, which results in rail wear. The Specht wear model is used to calculate the wear between wheel and rail for this case by applying different wear parameters for moderate wear and severe wear, respectively [34]. Specht wear model can be considered as a deformation of Archard model. The difference is that Specht wear model considers the multipoint contact between wheel and rail. This model introduces a jump factor on the basis of Archard model to characterize the multipoint contact condition between the wheel and the rail. When there is a multipoint surface contact between the wheel and rail, it indicates that the wheel-rail wear is serious, and the jump factor takes the empirical value of 10. Otherwise, it can be considered as a single-point contact between the wheel and rail, with slight wear, and the jump factor takes the empirical value of 1. This model assumes that the friction work generated by wheel-rail contact is linear with the volume of wear. Then, the wear volume V can be derived as follows.(7)V=kV⋅km⋅W,where W is the friction work on rail due to wheel-rail contact, and kV is a dimensionless wear coefficient related to wheel-rail material. km is the jump coefficient, and the empirical value is 1 for mild wear and 10 for severe wear [35].The wear depth is represented by the average wear depth regardless of the difference of wear depth within a contact patch, it can be derived as follows:(8)d=kV⋅km⋅FwfDHAc,where d is wear depth. D is the sliding distance. H is the hardness index, and Ac is the size of the contact area.By substituting equation (6) into (8), the wear model of outer rail is obtained as follows:(9)d=kV⋅kmH⋅DAcmwg⋅1gv2R+1gt.kh⋅v2R−Δh⋅1+sL212⋅1+s212−μ1G⋅zwcosδR+ϕw32+14mc+12mb+mwghcgt2ghgt−v2R.It is noted from equation (9) that the rail side wear depends on many parameters, such as the curve radius (R), gauge (gt), running speed (v), and load (mc, mb, mw). These parameters can be clustered into four categories. The first category includes track parameters, which consist of the curve radius (R), gauge (gt), track longitudinal gradient (sL), rail cant (s), and actual superelevation (h). The second category is wheel-rail contact parameters (Ac, zw, δR, ϕw, D), which are related to the wheel-rail profiles as well as the vehicle operation status. In addition, the vehicle parameters are considered the third category, which consists of running speed (v), axle load (mc, mb, mw), and the vehicle type characterized by the stiffness and damping of the suspension system. Finally, the fourth category includes the maintenance activities, where the grinding, lubrication, and weather conditions are all captured by the coefficient of friction (μ).As a number of parameters are required to describe the wear profiles of the wheel and rail, it is not convenient to measure directly in practical applications. In addition, the core objective of this paper is to predict rail wear with the help of directly measurable parameters. Therefore, the influence of the wheel-rail profile is not considered, while the other twelve parameters obtained from equation (9) are further explored in the sensitivity analysis presented in section 2.3. The simulation model used for rail wear calculation is introduced in the next section. With the help of the dynamic software Universal Mechanism (UM) [36], the rail wear produced by vehicles negotiating different curves is calculated for further analysis.


## 2.2. Simulation Environment

In this section, a subway vehicle-track dynamic model is built by using UM for rail wear calculation by using the vehicle and track parameters of Guangzhou metro as an example. The UM software includes different models of railway rolling stock including simplified models and 3D models, and other different software tools for the simulation of railway rolling stock dynamics are also included in this multi-body dynamics software. The vehicle (including the car-body, primary and secondary suspension systems, bogie system, etc.), track, and wheel-rail contact system in the vehicle-track coupling system in terms of “Universal mechanism” are all a subsystem, which can also be models of any complexity.The subway vehicle-track dynamic model is built referring to the parameters of Guangzhou metro vehicle and track, and the simulation calculation of rail superposition wear is carried out by combining the wheel-rail non-Hertz contact coupling dynamic model and Specht material wear model. The feasibility of the simulation model is verified by the on-site tracking monitoring data. Based on the structure of Guangzhou Metro A vehicle, a vehicle-ballasted track coupling dynamic model is established by using UM multi-body dynamics software. The vehicle model is a nonlinear model, including the car body, two bogies, four wheelsets and suspension devices, in which the car body and the bogie are regarded as rigid bodies and a total of 68 degrees of freedom are generated. The rail is considered as Euler beams supported by discrete elastic points distributed according to sleeper spacing. Detailed modelling parameters are shown in TableS2 in the Appendix. The total length of the line is 2 km. Consider that the track random irregularity has a significant influence on the nonuniform distribution rail wear along the rolling direction [37, 38], the track irregularity data of Guangzhou Metro Line 1 is used in the simulation study.The Specht wear model [34] is used to calculate the rail wear at the midpoint of the circular curve. The normal load of wheel-rail contact is calculated by multipoint contact and the normal parameters of wheel-rail contact is obtained by Kik–Piotrowski program [39]. The FASTSIM algorithm is used to calculate the tangent contact parameters of the wheel-rail contact system in dynamic calculations. According to the assumption of linear relationship between friction work and wear depth in Specht wear model, the wear depth of each point can be obtained.The calculation of rail superimposed wear is proceeded as follows. The normal wear threshold of rail is set todT, and the rail profile in the simulated model is updated for the next iteration once the maximum normal wear dimax of a certain point on the rail contact spot reaches the threshold. This paper assumes that the rail profile does not change in the same iterative calculation, that is, the rail profile when the first wheelset passes is consistent with that when the n th wheelset passes. Considering the wheel-rail contact state when the i th wheelset passing, the wear depth after the i th wheelset passes as Δi, and the cumulative wear depth Δ (Δ=∑i=1nΔi) after a train passes can be obtained. When the SRW reaches the requirement of rail replacement (Guangzhou Metro requires 13 mm), the iteration calculation is completed. The wear data and the number of cumulative vehicles passing by (NVP) are the outputs. In this simulation, the empirical value of 0.1 mm is taken as the wear threshold. The simulation model is adjusted based on the analysis results of wear measurements of Guangzhou Metro Line 1. The material wear coefficient kV calculated based on the field measurement results of side wear is 2.103e−13 m3/J.To validate if the model can be used for the further analysis, the monitoring rail wear of Yangji-Sports Centre interval of Guangzhou Metro Line 1 and the simulation results are introduced for comparative analysis. Figure2 shows the comparison of the on-the-spot measurement of SRW and the simulation results, and the Y-axis represents the SRW measured in the position of 16 mm below the rail top. The X-axis represents the number of cumulative vehicles passing by (NVP). It should be noted that 6-car formation is adopted for the train in this paper, and 4 wheelsets in each vehicle is used to calculate the rail wear in the outer curved track. Therefore, when a train passes by a curved track, a total of 24 wheelsets pass through a certain position of the outer rail in turn. In Figure 2, the radius of the measured curve is 500 m, and the lines are all underground, where the influence of environmental factors on rail wear can be captured by the friction coefficient (μ=0.4). See curve C6 in Table S3 in the Appendix for details of the line parameters and measuring points. The figure shows the monitoring data of the side rail wear of the curve from January 2014 to October 2018. During this period, no rail grinding or replacement has been performed. It shows that the measured wear has a linear growth trend with the cumulative increase in the vehicles passing by, and the numerical simulation results are in good agreement with the measurements, which verifies the feasibility of the simulation model for rail wear calculation.Figure 2
Measured and simulated wear curve.


## 2.3. Sensitivity Analysis

In light of the theoretical analysis of SRW, the influence of these parameters on wear, including the radius of the curve, axle load, friction coefficient, unbalanced superelevation, rail cant, track gradient, running speed, and stiffness of primary and secondary suspensions, cannot be ignored. To obtain the most crucial influential factors, this section introduces the global sensitivity analysis based on Sobol' method [40, 41], which is widely used at present. It studies the global influence of different parameters on the model. The range of the parameters can be extended to the whole definition domain, and different parameters can change simultaneously, reflecting the interaction between the parameters.The central idea of Sobol' method is to decompose the functionfx into 2n increasing terms, as follows:(10)fx=f0+∑i=1nfixi+∑1≤i<j≤nfi,jxi,xj+⋯+f1,2,…,nx1,x2,…,xn,where f0 is a constant term, and the integral of the other subterms to any factor is zero.The decomposition of the function is unique, and each subterm can be expressed as an integral of the functionfx. The left and right sides of equation (10) in the entire parameter domain In are squared and integrated, where In=x|0≤yi≤1;i=1,2,…,n, and(11)∫Inf2xdX−f02=∑s=1n∑i1<i2<...<isn∫fi1,i2,…,is2dxi1dxi2…dxis.The total variance is as follows:(12)D=∫f2xdX−f02.The partial variance of each order is(13)Di1,i2,…,is=∫01⋯∫01fi1,i2,…,is2xi1,xi2,…,xisdxi1dxi2…dxis.The ratio of each order's partial variance to the total variance is each order's sensitivity, and theS-order sensitivity is as follows:(14)Si1,i2,…,is=Di1,i2,…,isD.For the rail wear studied in this paper,Si is the first-order sensitivity, which shows the influence of the parameter xi on the rail wear. Si,ji≠j is a second-order sensitivity value, which shows the influence of the parameters xi and xj on the rail wear. At the same time, the coupling interaction between these two parameters is also considered. In other words, the total sensitivity not only reflects the influence of this parameter on the rail wear but also analyses the interaction with all other parameters. The larger the first-order sensitivity value of a parameter, the more significant the influence of this parameter is on the rail wear. When the difference between the first-order sensitivity value of a parameter and its total sensitivity value is large, it can be assumed that there is a coupling interaction between the parameter and other parameters [42].Before performing a global sensitivity analysis on rail wear under random parameters, the change interval and probability distribution of the parameters need to be given, and the parameters are randomly sampled within the interval. This paper builds a Sobol' sequence for random sampling calculation.The Sobol' sequence is obtained by constructing a group called the 'direction number' [43]. Assuming that mi is a positive odd number less than 2n, then(15)ai=mi2i.The direction numberai is generated by a simple polynomial fx with coefficients of only 0 and 1, as follows:(16)fx=xp+b1xp−1+⋯+bp−1x+bp.Fori>p, the recursive formula is as follows:(17)bi=b1ai−1⊕b2ai−2⊕⋯⊕bpai−p⊕ai−p2p,where ⊕ is the binary XOR. Combining the above equation with (15), we can obtain(18)mi=2b1mi−1⊕22b2mi−2⊕⋯⊕2pbpmi−p⊕mi−p.Thei th number of the Sobol' sequence is obtained as follows:(19)xi=c1a1⊕c2a2⊕c3a3⊕⋯,where c1c2c3… is a binary representation.The random number generated by the Sobol' sequence is more uniform, and there is no large blank or repetition phenomenon. Therefore, the distribution of random numbers can be described completely with fewer sampling points, which effectively improves the accuracy and efficiency of sensitivity analysis. Therefore, the Sobol' sequence method is used to randomly sample the parameters in this paper.


## 2.4. Prediction Models

Regression models can be used to find mathematical expressions of relatively complex physical processes and systems. This expression or calculation model can be regarded as a transfer function between a given input and an output [44]. Different types of regression models are available [6], such as the response surface methodology (RSM), support vector machine (SVM), relative vector machine (RVM), and artificial neural network (ANN). This study applies three different regression analysis methods, namely, RSM, SVM, and RVM, that fit our train rail wear data to predict the RUL of outer rails based on rail wear. The data points used for the regression fitting or training process are selected by the central combination design (CCD) [45, 46].


### 2.4.1. RSM Regression Analysis

RSM can solve problems where the relationship between the objective function and the independent variable is unknown, and it is used to model and analyze the problem in which the objective function is affected by multiple variables in order to optimize the objective function [47, 48].As already mentioned in the introduction, the purpose of this study is to calculate the RUL of the rail based on wear by using a regression prediction model. Hence, for the response value of this study, the service life of the rail based on wear is selected. Here, the accumulative number of vehicles passing by (NVP) the line during the rail service is taken as the service life of the rail. Then, the rail wear undern scenarios obtained by the CCD are calculated according to the calculation process of rail wear described in Section 2.1, and the response values corresponding to that are subjected to a regression analysis to obtain a predictive model.The prediction model obtained by means of RSM has a polynomial form with a certain order. In this study, a quadratic polynomial model is chosen for fitting and is given as follows:(20)yx=a0+∑i=1naixi+aiixi2+∑j=2n∑i<jjaijxixj+ε,where a0, ai, aii, and aij are fitting coefficients and $n$ indicates the number of analysis variables.Due to the complexity of RSM polynomial fitting, this study only considers the second-order coupling cross between the influencing factors. A more complex coupling effect of various parameters on wear cannot be taken into account, which will inevitably affect the prediction performance of the RSM model. For comparative analysis, the machine learning algorithm of SVM is introduced in the next section.


### 2.4.2. SVM Regression Analysis

SVM is a classification and regression analysis tool proposed by Vapnik et al. in 1995 [49] and currently is mainly used in pattern recognition and regression analysis. The greatest advantage of SVM lies in the solving of problems with a small sample and nonlinear and high-dimensional pattern recognition, which effectively overcomes the problems of large deviation, over learning, local extremum and large amounts of calculation. It can also be extended to other machine learning problems, such as function fitting [50].Given the training samplesX=x1,y1,x2,y2,…,xp,yp⊆Rn×R, where xi is the input variable, yi is the corresponding output variable, and p is the total number of sample points, according to the principle of support vector regression, the linear regression function is as follows:(21)fx=ωTϕx+b,where ω∈Rn is the weight vector of the optimal hyperplane, b is the offset threshold, ϕx is a nonlinear mapping function, and ωTϕx represents the point multiplication of two vectors.To determine the regression function coefficientsω and b, the above problem can be transformed into a constrained optimization problem as follows:(22)maxωai,ai∗=−12∑i,j=1pai−ai∗aj−aj∗Kxi,xj−ε∑i,j=1pai+ai∗+∑i=1pyiai−ai∗,∑i=1pai−∑i=1pai∗=0,0≤ai,ai∗≤C,i=1,2,…,p,where Kxi,xj=ϕxiT⋅ϕxj, Kxi,xj is the kernel function.By solving equation (22), the regression function of SVM can be obtained as(23)fx=∑i,j=1pai−ai∗Kxi,xj+b.Although SVM has been widely used in the field of prediction, its prediction results lack the ability to express uncertainty. Therefore, the RVM model with uncertainty expression is introduced in the next section for regression prediction analysis.


### 2.4.3. RVM Regression Analysis

RVM is a classification and regression method proposed by Tipping in 2000 [51] and is one of the important research hotspots in the field of statistical learning in recent years. It has a similar function form as SVM and solves the nonlinear classification problem by introducing a kernel function and dimensionality increase.Given a data setG=xi,tii=1N, xi∈Rd represents the input vector, ti∈R represents the corresponding output target vector, and N is the number of input samples. The training model can be expressed as:(24)t=yx+ε,where y⋅ is a nonlinear function and ε is Gaussian white noise with variance σ2, that is, ε∼N0,σ2. The approximation function y can be obtained by a regression operation using the training data. The mathematical expression of RVM can be expressed as(25)t=Φω+εn,where ω=ω1,ω2,…,ωNT, ω is an N+1 dimensional vector representing the weight of RVM. Φ is the kernel function matrix, and Φ=ϕ1x,ϕ2x,…,ϕNxT, ϕix=1,Kxi,x1,…,Kxi,xN, i=1,2,…,N, and K⋅ is a kernel function.The Bayesian inference based on the Gaussian process is applied to kernel theory, and the uncorrelated sample points are removed by an autocorrelation decision theory under the structure of prior parameters. Hence, the RVM sparse probability model is obtained. When a new test sample pointxN+1 is input into the RVM, its prediction model can be expressed as(26)ptN+1|t=∫ptN+1|ω,a,σ2pω,a,σ2|tdωdadσ2,where tN+1 represents the observation target of the new test sample point xN+1.For a given new input valuex∗, the corresponding predicted output value t∗ satisfies the following Gaussian distribution:(27)pt∗|ω=∫pt∗|ω,aMP,σMP2pω|t,aMP,a,σMP2dω∼Ny∗,σ∗2,where y∗ and σ∗2 are expressed as follows:(28)y∗=μTΦX∗,(29)σ∗2=σMP2+ΦX∗T∑ΦX∗,where y∗ is the predicted output value of the RVM model with the new test data x∗, and the first term of equation (29) represents the variance of the noise estimation, and the second term represents the uncertainty of the weight estimation.To compare the prediction errors of the above prediction models, four evaluation indexes of the prediction performance are used in this paper, namely, the coefficient of determination (R2), root mean square error (RMSE), mean absolute error (MAE), mean relative error (MRE), and normalized mean square error (NMSE). The calculation formulas are as follows:(30)R2=1−∑y−y^2∑y−y¯2,RMSE=1N∑i=1Nyi−y^i2,MAE=1N∑i=1Nyi−y^i,MRE=1N∑i=1Nyi−y^iyi×100%,where yi and y^i represent the real value and predicted value, respectively, N is the total number of test samples.The closer the value ofR2 is to 1, and the closer the values of the other three indicators are to 0, indicating that the closer the prediction result is to the true value, the better the prediction performance of the regression model.


### 2.4.4. Determination of Kernel Function and Hyperparameters

SVM and RVM models with kernel function of Gaussian Radial Basis Function (RBF) have good nonlinear processing ability. When dealing with practical engineering data, it is better to choose RBF kernel function than other kernel functions in the absence of prior knowledge [52].The hyperparameter of the regression model has a very significant impact on the accuracy of the prediction results while SVM and RVM regression models are used to predict the development trend of rail wear. In this paper, thek-fold cross validation procedure is introduced to select the penalty factor C and the parameter σ of RBF kernel function [52], so as to make the prediction result reach the optimal value. In this method, the original data is divided into k-fold, and each subset data is used as a verification set, and the remaining k−1 subset data is used as the training set, so K models can be obtained. The K models were evaluated in the validation set, and the final MSE (mean squared error) was added and averaged to get the cross validation error. The k-fold cross validation is suitable for the limited simulation data in this paper, and the evaluation results can be as close as possible to the performance of the model on the test set, so it can be used as an index for model optimization. In this paper, 10-fold cross validation is selected to optimize the model parameters.The steps to optimize model parameters are as follows:(1)
Determine the initial search interval and step size of the parameter pair (C, σ) as C∈a,b, σ∈c,d. The search steps of parameters C and σ are L0 and L1 respectively, and the k-fold cross-validation procedure is used to calculate the MSE for the selected parameters to evaluate the generalization ability of the model. According to the global minimum criterion of MSE, a preliminary optimal parameter pair is obtained as (C0, σ0). In this paper, the exponential form with the bottom of 2 is selected to perform the initial division of the parameter grid for fast search. For example, L0=1, L1=0.5, and C=2−a,2−a+1,…,2b−1,2b, σ=2−c,2−c+0.5,…,2d−0.5,2d.(2)
Expand the preliminary optimal parameter pair (C0, σ0) to obtain the quadratic grid search intervals of parameters C and σ as m,n, p,q. In the new grid, a further fine search is performed with step lengths L0′ and L1′, and the step length of the fine search satisfies: L0'<L0,L1′<L1. Using the criterion of MSE to take the global minimum, the optimal parameters (C∗,σ∗) can be obtained after the second fine grid search.The above steps are based on SVM model. In the RVM regression model, the penalty factorC is assigned automatically, so only the parameter optimization of the kernel function parameter σ is needed. The method and process are similar to it and will not be repeated.


## 2.4.1. RSM Regression Analysis

RSM can solve problems where the relationship between the objective function and the independent variable is unknown, and it is used to model and analyze the problem in which the objective function is affected by multiple variables in order to optimize the objective function [47, 48].As already mentioned in the introduction, the purpose of this study is to calculate the RUL of the rail based on wear by using a regression prediction model. Hence, for the response value of this study, the service life of the rail based on wear is selected. Here, the accumulative number of vehicles passing by (NVP) the line during the rail service is taken as the service life of the rail. Then, the rail wear undern scenarios obtained by the CCD are calculated according to the calculation process of rail wear described in Section 2.1, and the response values corresponding to that are subjected to a regression analysis to obtain a predictive model.The prediction model obtained by means of RSM has a polynomial form with a certain order. In this study, a quadratic polynomial model is chosen for fitting and is given as follows:(20)yx=a0+∑i=1naixi+aiixi2+∑j=2n∑i<jjaijxixj+ε,where a0, ai, aii, and aij are fitting coefficients and $n$ indicates the number of analysis variables.Due to the complexity of RSM polynomial fitting, this study only considers the second-order coupling cross between the influencing factors. A more complex coupling effect of various parameters on wear cannot be taken into account, which will inevitably affect the prediction performance of the RSM model. For comparative analysis, the machine learning algorithm of SVM is introduced in the next section.


## 2.4.2. SVM Regression Analysis

SVM is a classification and regression analysis tool proposed by Vapnik et al. in 1995 [49] and currently is mainly used in pattern recognition and regression analysis. The greatest advantage of SVM lies in the solving of problems with a small sample and nonlinear and high-dimensional pattern recognition, which effectively overcomes the problems of large deviation, over learning, local extremum and large amounts of calculation. It can also be extended to other machine learning problems, such as function fitting [50].Given the training samplesX=x1,y1,x2,y2,…,xp,yp⊆Rn×R, where xi is the input variable, yi is the corresponding output variable, and p is the total number of sample points, according to the principle of support vector regression, the linear regression function is as follows:(21)fx=ωTϕx+b,where ω∈Rn is the weight vector of the optimal hyperplane, b is the offset threshold, ϕx is a nonlinear mapping function, and ωTϕx represents the point multiplication of two vectors.To determine the regression function coefficientsω and b, the above problem can be transformed into a constrained optimization problem as follows:(22)maxωai,ai∗=−12∑i,j=1pai−ai∗aj−aj∗Kxi,xj−ε∑i,j=1pai+ai∗+∑i=1pyiai−ai∗,∑i=1pai−∑i=1pai∗=0,0≤ai,ai∗≤C,i=1,2,…,p,where Kxi,xj=ϕxiT⋅ϕxj, Kxi,xj is the kernel function.By solving equation (22), the regression function of SVM can be obtained as(23)fx=∑i,j=1pai−ai∗Kxi,xj+b.Although SVM has been widely used in the field of prediction, its prediction results lack the ability to express uncertainty. Therefore, the RVM model with uncertainty expression is introduced in the next section for regression prediction analysis.


## 2.4.3. RVM Regression Analysis

RVM is a classification and regression method proposed by Tipping in 2000 [51] and is one of the important research hotspots in the field of statistical learning in recent years. It has a similar function form as SVM and solves the nonlinear classification problem by introducing a kernel function and dimensionality increase.Given a data setG=xi,tii=1N, xi∈Rd represents the input vector, ti∈R represents the corresponding output target vector, and N is the number of input samples. The training model can be expressed as:(24)t=yx+ε,where y⋅ is a nonlinear function and ε is Gaussian white noise with variance σ2, that is, ε∼N0,σ2. The approximation function y can be obtained by a regression operation using the training data. The mathematical expression of RVM can be expressed as(25)t=Φω+εn,where ω=ω1,ω2,…,ωNT, ω is an N+1 dimensional vector representing the weight of RVM. Φ is the kernel function matrix, and Φ=ϕ1x,ϕ2x,…,ϕNxT, ϕix=1,Kxi,x1,…,Kxi,xN, i=1,2,…,N, and K⋅ is a kernel function.The Bayesian inference based on the Gaussian process is applied to kernel theory, and the uncorrelated sample points are removed by an autocorrelation decision theory under the structure of prior parameters. Hence, the RVM sparse probability model is obtained. When a new test sample pointxN+1 is input into the RVM, its prediction model can be expressed as(26)ptN+1|t=∫ptN+1|ω,a,σ2pω,a,σ2|tdωdadσ2,where tN+1 represents the observation target of the new test sample point xN+1.For a given new input valuex∗, the corresponding predicted output value t∗ satisfies the following Gaussian distribution:(27)pt∗|ω=∫pt∗|ω,aMP,σMP2pω|t,aMP,a,σMP2dω∼Ny∗,σ∗2,where y∗ and σ∗2 are expressed as follows:(28)y∗=μTΦX∗,(29)σ∗2=σMP2+ΦX∗T∑ΦX∗,where y∗ is the predicted output value of the RVM model with the new test data x∗, and the first term of equation (29) represents the variance of the noise estimation, and the second term represents the uncertainty of the weight estimation.To compare the prediction errors of the above prediction models, four evaluation indexes of the prediction performance are used in this paper, namely, the coefficient of determination (R2), root mean square error (RMSE), mean absolute error (MAE), mean relative error (MRE), and normalized mean square error (NMSE). The calculation formulas are as follows:(30)R2=1−∑y−y^2∑y−y¯2,RMSE=1N∑i=1Nyi−y^i2,MAE=1N∑i=1Nyi−y^i,MRE=1N∑i=1Nyi−y^iyi×100%,where yi and y^i represent the real value and predicted value, respectively, N is the total number of test samples.The closer the value ofR2 is to 1, and the closer the values of the other three indicators are to 0, indicating that the closer the prediction result is to the true value, the better the prediction performance of the regression model.


## 2.4.4. Determination of Kernel Function and Hyperparameters

SVM and RVM models with kernel function of Gaussian Radial Basis Function (RBF) have good nonlinear processing ability. When dealing with practical engineering data, it is better to choose RBF kernel function than other kernel functions in the absence of prior knowledge [52].The hyperparameter of the regression model has a very significant impact on the accuracy of the prediction results while SVM and RVM regression models are used to predict the development trend of rail wear. In this paper, thek-fold cross validation procedure is introduced to select the penalty factor C and the parameter σ of RBF kernel function [52], so as to make the prediction result reach the optimal value. In this method, the original data is divided into k-fold, and each subset data is used as a verification set, and the remaining k−1 subset data is used as the training set, so K models can be obtained. The K models were evaluated in the validation set, and the final MSE (mean squared error) was added and averaged to get the cross validation error. The k-fold cross validation is suitable for the limited simulation data in this paper, and the evaluation results can be as close as possible to the performance of the model on the test set, so it can be used as an index for model optimization. In this paper, 10-fold cross validation is selected to optimize the model parameters.The steps to optimize model parameters are as follows:(1)
Determine the initial search interval and step size of the parameter pair (C, σ) as C∈a,b, σ∈c,d. The search steps of parameters C and σ are L0 and L1 respectively, and the k-fold cross-validation procedure is used to calculate the MSE for the selected parameters to evaluate the generalization ability of the model. According to the global minimum criterion of MSE, a preliminary optimal parameter pair is obtained as (C0, σ0). In this paper, the exponential form with the bottom of 2 is selected to perform the initial division of the parameter grid for fast search. For example, L0=1, L1=0.5, and C=2−a,2−a+1,…,2b−1,2b, σ=2−c,2−c+0.5,…,2d−0.5,2d.(2)
Expand the preliminary optimal parameter pair (C0, σ0) to obtain the quadratic grid search intervals of parameters C and σ as m,n, p,q. In the new grid, a further fine search is performed with step lengths L0′ and L1′, and the step length of the fine search satisfies: L0'<L0,L1′<L1. Using the criterion of MSE to take the global minimum, the optimal parameters (C∗,σ∗) can be obtained after the second fine grid search.The above steps are based on SVM model. In the RVM regression model, the penalty factorC is assigned automatically, so only the parameter optimization of the kernel function parameter σ is needed. The method and process are similar to it and will not be repeated.


## 3. Results and Discussion

The evolution process of the rail profile due to side wear is described in the first part of this section. It is important to complete this process before proceeding with regression analysis as the wear profile evolution due to SRW is related to whether the prediction model needs to be analyzed in stages and needs to be accurately described. Thereafter, the results for the sensitivity analysis and the comparative analysis of the prediction models are presented.


### 3.1. The Wear Profile Evolution due to SRW

The side wear of the outer rail at SRCs is calculated by using a dynamic simulation model, and the evolution of the worn rail profile during its whole life stage is calculated as shown in Figure3, where the curve radius is 500 m, and other parameters of the curve are introduced in Table S3 under curve C6 in the Appendix.Figure 3
The wear profile evolution due to SRW.
(a)(b)(c)(d)Figure3 shows that the evolution of SRW is mainly divided into 3 stages, and this is in accordance with the results of [51].Stage I. This stage occurs during the period from when the new rail starts its service to the time when wear occurs at the position of 16 mm below the top of the track, as shown in Figure3(b). In this stage, the wear concentrates on the gauge corner and gradually evolves to the area of the rail indicated in Figure 3(c) as the NVP increases. Wear has not yet occurred at 16 mm below the rail top in this stage, so mild wear at this time cannot be detected according to the definition of side wear.Stage II. In this stage, SRW occurs and increases linearly with increasing NVP. In addition, the wear position continues to expand to the top of the rail as well as the lower part of the rail side, as shown in Figure3(c). Generally, catastrophic wear occurs when the SRW reaches 6 mm–8 mm, indicating the end of the severe wear stage.Stage III. As the wear position continues to evolve downward, the SRW obviously intensifies. Then, the wear position gradually remains below the rail top by 26 mm–27 mm and does not continue to extend downward. At this stage, the rail wear is catastrophic and has the appearance of plastic deformation of the rail, as shown in Figure3(d).


### 3.2. Sensitivity Analysis of SRW

In this section, the global sensitivity analysis described in section2.3 is performed to evaluate the importance of each influencing factor by calculating the fluctuation of the response value caused by the change in a single factor or a combination of parameters. Before that, the concept of the average wear rate (SWR) of rail side wear rwi is introduced as the response value of the sensitivity analysis, and it indicates the increment of SRW caused by the NVP at a certain time i. The specific definition of SWR is shown as follows:(31)rwi=ΔdΔN,where, Δd stands for the increment in SRW in a certain period of time i. ΔN denotes the NVP on the line in a corresponding period.In light of the theoretical analysis of SRW, the occurrence of SRW at SRCs is the result of many factors interacting together that are closely related to the track parameters, the vehicle parameters, the operating conditions and the wheel-rail contact state. Due to specific modelling challenges, not all parameters shown in equation (9) are explored in this study. This paper concerns the effect of the available track parameters and vehicle parameters on the rail service life based on wear, with the aim of evaluating the rail wear state in a timely manner and predicting the remaining life of the rail. The parameters to be further explored in the sensitivity analysis and their respective intervals are presented in Table 1, and these parameter ranges are selected based on field experience. Furthermore, parameter fluctuations related to lubrication, pollution and weather conditions are all reflected by the friction coefficient.Table 1
Distribution range of influence parameters.
ParameterPhysical meaningUnitLower limitUpper limitvSpeedkm/h4580sRail cant‰2530sLLongitudinal slope of line‰-2020ΔhUnbalanced superelevationmm-4040gtGaugemm14331441RRadius of circular curvem200900μFriction coefficient—0.20.6laAxle loadkg875013400k1xPrimary longitudinal stiffnessN/m52000007800000k1yPrimary transverse stiffnessN/m52000007800000k2xSecondary longitudinal stiffnessN/m20320003048000k2ySecondary transverse stiffnessN/m20320003048000Following the strategy of Sobol' global sensitivity analysis, the Sobol' sequence is determined to randomly sample the given parameters within the interval and to calculate the response values in this scenario by means of the simulation model. Since the core problem of this paper is to identify the rail wear state by using the performance state parameters that can be measured online, the local effects caused by the wheel-rail profile are neglected. CHN60 is chosen as the rail profile, and s1002 is chosen as the wheel profile. The responses for each scenario obtained by the Sobol' sequence are calculated, and the calculation results of the sensitivity indexes are shown in Figure4.Figure 4
Global sensitivity analysis results.The friction coefficient has the largest influence on rail wear. The sensitivity of vehicle suspension parameters (k1x, k1y, k2x, k2y) is relatively high among the 12 factors. In addition, the sensitivity values of the curve radius, gauge and axle load are slightly less than that of the primary suspension parameter. The first-order sensitivity values of unbalanced super-elevation and rail cant are 3.122 and 3.883, respectively, and have a significant effect on rail wear. Among the above factors, the sensitivity value of the track gradient is lower than that of the running speed (0.7514) and has the least impact on rail wear. Furthermore, the total sensitivity values of the other influencing parameters are significantly higher than the first-order sensitivity values except for that of the track gradient. This analysis result is consistent with the theoretical analysis results. For the wear of the curved track, there exists a certain interaction and coupling on the rail wear among the track parameters, vehicle structure parameters and operating conditions parameters.The friction coefficient directly affects the contact roughness between the wheel and the rail, which in turn significantly affects the friction and wear between the wheel-rail materials. As for the stiffness of the suspension system, due to the stiffness of the suspension material affects the structural flexibility of the vehicle when cornering, it indirectly affects the wear and tear of materials. The curve radius plays a decisive role in the contact state of the wheel and rail, which determines the function degree of wheel flange guidance when the vehicle is cornering, so it is also the key factor affecting the rail wear. Here, only the influence of the longitudinal stiffness of the primary suspension and the curve radius on the rail wear is taken as an example to illustrate. The calculation results show that with the increase of the longitudinal stiffness of the primary suspension system, the average wear rate of rail side wear increases linearly. Within the range of variation, the average wear rate of rail side wear increases by 0.08 mm/104 train. When the stiffness reaches 6825000 N/m, the average wear rate of rail side wear slows down obviously, and the increment is only 0.003 mm/104 train. Therefore, the spring stiffness can be appropriately reduced to reduce rail wear on the premise that the rotation coefficient of suspension system is not more than 0.1, which meets the requirements of national standards.The radius of circular curve, as the main cause of side wear of curve outer rail, has a significant effect on the average wear rate of side wear. When the curve radius is increased from 200 m to 900 m, the average wear rate on the outer rail side of the curve is reduced from 0.5396 mm/104 train to 0.0727 mm/104 train. It can be seen that with the increase of the radius of the curve, the rail-side grinding becomes slower and slower, with the radius of 400 m as the limit, and the average wear rate of the rail-side mill has a significantly different rate of change. On the small radius curve with the radius between 200 m and 400 m, with the increase of the curve radius, the average wear rate of rail side wear decreases rapidly from 0.5396 mm/104 train to 0.2244 mm/104 train, and the average wear rate of rail decreases by 29.2% with the increase of curve radius of 100 m. And on the curve with the radius between 400 m and 900 m, along with the increase of the curve radius, the average wear rate of rail side wear is reduced more slowly, from 0.2244 mm/104 train to 0.0727 mm/104 train, and the average wear rate with the increase of curve radius of 100 m is only reduced by 7.5%. Therefore, it is considered that the rail side wear is relatively slow on the curve with radius greater than 400 m. While on the curve with small radius less than 400 m, the rail side wear increases rapidly with the decrease of the curve radius.In summary, the first-order sensitivity value of each parameter enables the selection of the most important influencing factors, and its comparison with the total sensitivity value reveals in more detail how the coupling relationship between these parameters affects the responses. After the global sensitivity analyses, the seven most important parameters are selected to be applied in the following derivation of the regression prediction model in the next section: friction coefficientμ, primary longitudinal stiffness k1x, radius of circular curve R, track gauge gt, axle load la, unbalanced superelevation Δh, and cant s.


### 3.3. Prediction Modelling Results

This section presents the purpose of the response values used for the regression prediction process and the comparative analysis of applying different regression models to prediction.


#### 3.3.1. Determination of Response Values

With the aim of predicting the RUL of a rail based on side wear, the service life of rail is introduced as the response variable of the regression prediction model. It is defined as the cumulative number of vehicles passing by during rail service. In light of the results from Section3.1, there are three wear regimes of rail side wear, namely, mild, severe and catastrophic. The corresponding service lives for these three stages are T1, T2, T3. Due to different types of wear mechanisms, it is not appropriate to capture different wear states with a single predictive model. Hence, it is necessary to develop separate regression prediction models corresponding to the various wear regimes. As the SRW can be alleviated in a timely manner and the service life of the curve can be extended if maintenance measures are carried out before the SRW evolves into stage III, prediction models of rail service life in middle and sever wear are proposed by using the abovementioned seven parameters that can be measured online. The response value T1 for CCD scenarios in the mild wear regime is designed to be the cumulative vehicles passing by (NVP, unit in 104 trains) the line from the time the rail is put into use until side wear occurs. Similarly, the NVP during the severe wear regime is proposed as the response value T2 for CCD scenarios in wear stage II.An important step before proceeding with the regression prediction is the determination of the data points and the simulation scenarios. Considering that the two prediction methods of SVM and RVM have no specific requirements for the settings of the sample points, this article follows the solution method of the RSM prediction model to determine the simulation scenarios [47]. First, the parameters for regression analysis are normalized. The seven key factors extracted from the sensitivity analysis are standardized according to the upper and lower limits of the parameters set in Table 1.Define the normalized variables of the original parametersμ, R, k1x, gt, la, s, and Δh as xi, i=1,2,3,4,5,6,7, it is(32)xi=Xi−x¯iui−x¯i,where, Xi is the original variable of the normalized variable xi, x¯i=1/2ui−li, ui and li are the upper and lower limits of xi, respectively.After that, the simulation scenarios are designed by the central combination design method (CCD) [45]. This method reduces the total number of observations to be evaluated while maintaining a certain accuracy of the second-order response surface model. Following this method, the number of observations per regression model is designed as 152 for the two wear regimes of mild and severe. The predictive results are described in the next section.


#### 3.3.2. Comparative Analysis of RUL Prediction

Polynomial fitting and training are performed on all simulation scenes and response values in different wear regimes, and an RSM-based model of quadratic form as given in equation (20) is fitted. Figure 5 shows the prediction results for the three regression prediction models.Figure 5
Comparison of predicted and observed values from numerical simulation. (a) Wear regime I. (b) Wear regime II.
(a)(b)Figure5(a) shows that the prediction results (the rail service life in wear regime I) of the three prediction models are close to the simulation values and show good prediction stability over time. By contrast, the prediction results of rail service life for wear regime II in Figure 5(b) show that the predicted values based on RSM and SVM deviate significantly from the simulated values over time, confirming the lack of long-term trend prediction accuracy. Furthermore, these three prediction models based on RSM, SVM, and RVM can effectively predict the rail service life due to side wear, which can preliminarily prove the effectiveness of the three methods proposed in this paper. However, it can also be seen from the figure that the prediction results based on the RSM and SVM models deviate from the true value in a relatively large range. In comparison, the prediction results based on the RVM model are closer to the true value with less fluctuation. This shows that the prediction method based on the RVM model has less uncertainty and the prediction effect is more stable.To compare the prediction performance of the above three different prediction models in more detail, Table2 gives a numerical comparison of the prediction performance indicators.Table 2
Prediction performance indexes of different prediction models based on simulation data.
ModelNRSMSVsRVsR2RMSEMAEMRE(a) Wear regime IRSM152——0.90110.067640.58210.8314%SVM—133—0.93040.035720.30315.7663%RVM——130.98860.013430.11322.1096%(b) Wear regime IIRSM152——0.91420.08300.065112.953%SVM—135—0.94910.06190.048110.1997%RVM——140.98820.02720.01993.3516%Due to space limitations, Table2 (a) is only used as an example for explanation. The RMSE value of RVM is 0.1343, which is significantly lower than that of RSM and SVM, indicating that the prediction accuracy based on the RVM prediction model is higher. The same conclusion can be obtained by comparing the other two prediction error evaluation indexes. The coefficients of determination R2 of these models are 0.9011, 0.9304, and 0.9886, respectively. Although the R2 of SVM and RVM are slightly larger, those of the three models are close to 1, which shows that the fitting degree of these prediction models is high, and the RVM model has the highest fitting degree. Additionally, compared with the number of RSM training samples of 152, the number of support vectors (SVs) of the multivariate SVM is 133, accounting for almost 87.5% of the total number of training samples (152). The number of correlation vectors (RVs) of the multivariate RVM is only 13, accounting for only 8.6% of the total number of training samples, which shows that RVM has stronger sparsity in the application of regression prediction.As explained before, the RVM prediction model has the ability to express uncertainty. The upper and lower bounds of the prediction error constitute the 95% confidence interval for the prediction value. For the prediction of rail service life in wear regime I, the number of simulation values falling into the RVM prediction confidence interval is 140, accounting for 92.1% of the total number of training samples. In comparison, 137 simulation values of rail service life in wear regime II fall into the RVM prediction confidence interval, accounting for 90.13% of the total training samples. This verifies the prediction accuracy of the multivariable RVM prediction model.In summary, all the above conclusions are obtained by simulation tests. In the next section, the rail wear measured from the Guangzhou Metro is introduced to verify these conclusions.


## 3.1. The Wear Profile Evolution due to SRW

The side wear of the outer rail at SRCs is calculated by using a dynamic simulation model, and the evolution of the worn rail profile during its whole life stage is calculated as shown in Figure3, where the curve radius is 500 m, and other parameters of the curve are introduced in Table S3 under curve C6 in the Appendix.Figure 3
The wear profile evolution due to SRW.
(a)(b)(c)(d)Figure3 shows that the evolution of SRW is mainly divided into 3 stages, and this is in accordance with the results of [51].Stage I. This stage occurs during the period from when the new rail starts its service to the time when wear occurs at the position of 16 mm below the top of the track, as shown in Figure3(b). In this stage, the wear concentrates on the gauge corner and gradually evolves to the area of the rail indicated in Figure 3(c) as the NVP increases. Wear has not yet occurred at 16 mm below the rail top in this stage, so mild wear at this time cannot be detected according to the definition of side wear.Stage II. In this stage, SRW occurs and increases linearly with increasing NVP. In addition, the wear position continues to expand to the top of the rail as well as the lower part of the rail side, as shown in Figure3(c). Generally, catastrophic wear occurs when the SRW reaches 6 mm–8 mm, indicating the end of the severe wear stage.Stage III. As the wear position continues to evolve downward, the SRW obviously intensifies. Then, the wear position gradually remains below the rail top by 26 mm–27 mm and does not continue to extend downward. At this stage, the rail wear is catastrophic and has the appearance of plastic deformation of the rail, as shown in Figure3(d).


## 3.2. Sensitivity Analysis of SRW

In this section, the global sensitivity analysis described in section2.3 is performed to evaluate the importance of each influencing factor by calculating the fluctuation of the response value caused by the change in a single factor or a combination of parameters. Before that, the concept of the average wear rate (SWR) of rail side wear rwi is introduced as the response value of the sensitivity analysis, and it indicates the increment of SRW caused by the NVP at a certain time i. The specific definition of SWR is shown as follows:(31)rwi=ΔdΔN,where, Δd stands for the increment in SRW in a certain period of time i. ΔN denotes the NVP on the line in a corresponding period.In light of the theoretical analysis of SRW, the occurrence of SRW at SRCs is the result of many factors interacting together that are closely related to the track parameters, the vehicle parameters, the operating conditions and the wheel-rail contact state. Due to specific modelling challenges, not all parameters shown in equation (9) are explored in this study. This paper concerns the effect of the available track parameters and vehicle parameters on the rail service life based on wear, with the aim of evaluating the rail wear state in a timely manner and predicting the remaining life of the rail. The parameters to be further explored in the sensitivity analysis and their respective intervals are presented in Table 1, and these parameter ranges are selected based on field experience. Furthermore, parameter fluctuations related to lubrication, pollution and weather conditions are all reflected by the friction coefficient.Table 1
Distribution range of influence parameters.
ParameterPhysical meaningUnitLower limitUpper limitvSpeedkm/h4580sRail cant‰2530sLLongitudinal slope of line‰-2020ΔhUnbalanced superelevationmm-4040gtGaugemm14331441RRadius of circular curvem200900μFriction coefficient—0.20.6laAxle loadkg875013400k1xPrimary longitudinal stiffnessN/m52000007800000k1yPrimary transverse stiffnessN/m52000007800000k2xSecondary longitudinal stiffnessN/m20320003048000k2ySecondary transverse stiffnessN/m20320003048000Following the strategy of Sobol' global sensitivity analysis, the Sobol' sequence is determined to randomly sample the given parameters within the interval and to calculate the response values in this scenario by means of the simulation model. Since the core problem of this paper is to identify the rail wear state by using the performance state parameters that can be measured online, the local effects caused by the wheel-rail profile are neglected. CHN60 is chosen as the rail profile, and s1002 is chosen as the wheel profile. The responses for each scenario obtained by the Sobol' sequence are calculated, and the calculation results of the sensitivity indexes are shown in Figure4.Figure 4
Global sensitivity analysis results.The friction coefficient has the largest influence on rail wear. The sensitivity of vehicle suspension parameters (k1x, k1y, k2x, k2y) is relatively high among the 12 factors. In addition, the sensitivity values of the curve radius, gauge and axle load are slightly less than that of the primary suspension parameter. The first-order sensitivity values of unbalanced super-elevation and rail cant are 3.122 and 3.883, respectively, and have a significant effect on rail wear. Among the above factors, the sensitivity value of the track gradient is lower than that of the running speed (0.7514) and has the least impact on rail wear. Furthermore, the total sensitivity values of the other influencing parameters are significantly higher than the first-order sensitivity values except for that of the track gradient. This analysis result is consistent with the theoretical analysis results. For the wear of the curved track, there exists a certain interaction and coupling on the rail wear among the track parameters, vehicle structure parameters and operating conditions parameters.The friction coefficient directly affects the contact roughness between the wheel and the rail, which in turn significantly affects the friction and wear between the wheel-rail materials. As for the stiffness of the suspension system, due to the stiffness of the suspension material affects the structural flexibility of the vehicle when cornering, it indirectly affects the wear and tear of materials. The curve radius plays a decisive role in the contact state of the wheel and rail, which determines the function degree of wheel flange guidance when the vehicle is cornering, so it is also the key factor affecting the rail wear. Here, only the influence of the longitudinal stiffness of the primary suspension and the curve radius on the rail wear is taken as an example to illustrate. The calculation results show that with the increase of the longitudinal stiffness of the primary suspension system, the average wear rate of rail side wear increases linearly. Within the range of variation, the average wear rate of rail side wear increases by 0.08 mm/104 train. When the stiffness reaches 6825000 N/m, the average wear rate of rail side wear slows down obviously, and the increment is only 0.003 mm/104 train. Therefore, the spring stiffness can be appropriately reduced to reduce rail wear on the premise that the rotation coefficient of suspension system is not more than 0.1, which meets the requirements of national standards.The radius of circular curve, as the main cause of side wear of curve outer rail, has a significant effect on the average wear rate of side wear. When the curve radius is increased from 200 m to 900 m, the average wear rate on the outer rail side of the curve is reduced from 0.5396 mm/104 train to 0.0727 mm/104 train. It can be seen that with the increase of the radius of the curve, the rail-side grinding becomes slower and slower, with the radius of 400 m as the limit, and the average wear rate of the rail-side mill has a significantly different rate of change. On the small radius curve with the radius between 200 m and 400 m, with the increase of the curve radius, the average wear rate of rail side wear decreases rapidly from 0.5396 mm/104 train to 0.2244 mm/104 train, and the average wear rate of rail decreases by 29.2% with the increase of curve radius of 100 m. And on the curve with the radius between 400 m and 900 m, along with the increase of the curve radius, the average wear rate of rail side wear is reduced more slowly, from 0.2244 mm/104 train to 0.0727 mm/104 train, and the average wear rate with the increase of curve radius of 100 m is only reduced by 7.5%. Therefore, it is considered that the rail side wear is relatively slow on the curve with radius greater than 400 m. While on the curve with small radius less than 400 m, the rail side wear increases rapidly with the decrease of the curve radius.In summary, the first-order sensitivity value of each parameter enables the selection of the most important influencing factors, and its comparison with the total sensitivity value reveals in more detail how the coupling relationship between these parameters affects the responses. After the global sensitivity analyses, the seven most important parameters are selected to be applied in the following derivation of the regression prediction model in the next section: friction coefficientμ, primary longitudinal stiffness k1x, radius of circular curve R, track gauge gt, axle load la, unbalanced superelevation Δh, and cant s.


## 3.3. Prediction Modelling Results

This section presents the purpose of the response values used for the regression prediction process and the comparative analysis of applying different regression models to prediction.


### 3.3.1. Determination of Response Values

With the aim of predicting the RUL of a rail based on side wear, the service life of rail is introduced as the response variable of the regression prediction model. It is defined as the cumulative number of vehicles passing by during rail service. In light of the results from Section3.1, there are three wear regimes of rail side wear, namely, mild, severe and catastrophic. The corresponding service lives for these three stages are T1, T2, T3. Due to different types of wear mechanisms, it is not appropriate to capture different wear states with a single predictive model. Hence, it is necessary to develop separate regression prediction models corresponding to the various wear regimes. As the SRW can be alleviated in a timely manner and the service life of the curve can be extended if maintenance measures are carried out before the SRW evolves into stage III, prediction models of rail service life in middle and sever wear are proposed by using the abovementioned seven parameters that can be measured online. The response value T1 for CCD scenarios in the mild wear regime is designed to be the cumulative vehicles passing by (NVP, unit in 104 trains) the line from the time the rail is put into use until side wear occurs. Similarly, the NVP during the severe wear regime is proposed as the response value T2 for CCD scenarios in wear stage II.An important step before proceeding with the regression prediction is the determination of the data points and the simulation scenarios. Considering that the two prediction methods of SVM and RVM have no specific requirements for the settings of the sample points, this article follows the solution method of the RSM prediction model to determine the simulation scenarios [47]. First, the parameters for regression analysis are normalized. The seven key factors extracted from the sensitivity analysis are standardized according to the upper and lower limits of the parameters set in Table 1.Define the normalized variables of the original parametersμ, R, k1x, gt, la, s, and Δh as xi, i=1,2,3,4,5,6,7, it is(32)xi=Xi−x¯iui−x¯i,where, Xi is the original variable of the normalized variable xi, x¯i=1/2ui−li, ui and li are the upper and lower limits of xi, respectively.After that, the simulation scenarios are designed by the central combination design method (CCD) [45]. This method reduces the total number of observations to be evaluated while maintaining a certain accuracy of the second-order response surface model. Following this method, the number of observations per regression model is designed as 152 for the two wear regimes of mild and severe. The predictive results are described in the next section.


### 3.3.2. Comparative Analysis of RUL Prediction

Polynomial fitting and training are performed on all simulation scenes and response values in different wear regimes, and an RSM-based model of quadratic form as given in equation (20) is fitted. Figure 5 shows the prediction results for the three regression prediction models.Figure 5
Comparison of predicted and observed values from numerical simulation. (a) Wear regime I. (b) Wear regime II.
(a)(b)Figure5(a) shows that the prediction results (the rail service life in wear regime I) of the three prediction models are close to the simulation values and show good prediction stability over time. By contrast, the prediction results of rail service life for wear regime II in Figure 5(b) show that the predicted values based on RSM and SVM deviate significantly from the simulated values over time, confirming the lack of long-term trend prediction accuracy. Furthermore, these three prediction models based on RSM, SVM, and RVM can effectively predict the rail service life due to side wear, which can preliminarily prove the effectiveness of the three methods proposed in this paper. However, it can also be seen from the figure that the prediction results based on the RSM and SVM models deviate from the true value in a relatively large range. In comparison, the prediction results based on the RVM model are closer to the true value with less fluctuation. This shows that the prediction method based on the RVM model has less uncertainty and the prediction effect is more stable.To compare the prediction performance of the above three different prediction models in more detail, Table2 gives a numerical comparison of the prediction performance indicators.Table 2
Prediction performance indexes of different prediction models based on simulation data.
ModelNRSMSVsRVsR2RMSEMAEMRE(a) Wear regime IRSM152——0.90110.067640.58210.8314%SVM—133—0.93040.035720.30315.7663%RVM——130.98860.013430.11322.1096%(b) Wear regime IIRSM152——0.91420.08300.065112.953%SVM—135—0.94910.06190.048110.1997%RVM——140.98820.02720.01993.3516%Due to space limitations, Table2 (a) is only used as an example for explanation. The RMSE value of RVM is 0.1343, which is significantly lower than that of RSM and SVM, indicating that the prediction accuracy based on the RVM prediction model is higher. The same conclusion can be obtained by comparing the other two prediction error evaluation indexes. The coefficients of determination R2 of these models are 0.9011, 0.9304, and 0.9886, respectively. Although the R2 of SVM and RVM are slightly larger, those of the three models are close to 1, which shows that the fitting degree of these prediction models is high, and the RVM model has the highest fitting degree. Additionally, compared with the number of RSM training samples of 152, the number of support vectors (SVs) of the multivariate SVM is 133, accounting for almost 87.5% of the total number of training samples (152). The number of correlation vectors (RVs) of the multivariate RVM is only 13, accounting for only 8.6% of the total number of training samples, which shows that RVM has stronger sparsity in the application of regression prediction.As explained before, the RVM prediction model has the ability to express uncertainty. The upper and lower bounds of the prediction error constitute the 95% confidence interval for the prediction value. For the prediction of rail service life in wear regime I, the number of simulation values falling into the RVM prediction confidence interval is 140, accounting for 92.1% of the total number of training samples. In comparison, 137 simulation values of rail service life in wear regime II fall into the RVM prediction confidence interval, accounting for 90.13% of the total training samples. This verifies the prediction accuracy of the multivariable RVM prediction model.In summary, all the above conclusions are obtained by simulation tests. In the next section, the rail wear measured from the Guangzhou Metro is introduced to verify these conclusions.


## 3.3.1. Determination of Response Values

With the aim of predicting the RUL of a rail based on side wear, the service life of rail is introduced as the response variable of the regression prediction model. It is defined as the cumulative number of vehicles passing by during rail service. In light of the results from Section3.1, there are three wear regimes of rail side wear, namely, mild, severe and catastrophic. The corresponding service lives for these three stages are T1, T2, T3. Due to different types of wear mechanisms, it is not appropriate to capture different wear states with a single predictive model. Hence, it is necessary to develop separate regression prediction models corresponding to the various wear regimes. As the SRW can be alleviated in a timely manner and the service life of the curve can be extended if maintenance measures are carried out before the SRW evolves into stage III, prediction models of rail service life in middle and sever wear are proposed by using the abovementioned seven parameters that can be measured online. The response value T1 for CCD scenarios in the mild wear regime is designed to be the cumulative vehicles passing by (NVP, unit in 104 trains) the line from the time the rail is put into use until side wear occurs. Similarly, the NVP during the severe wear regime is proposed as the response value T2 for CCD scenarios in wear stage II.An important step before proceeding with the regression prediction is the determination of the data points and the simulation scenarios. Considering that the two prediction methods of SVM and RVM have no specific requirements for the settings of the sample points, this article follows the solution method of the RSM prediction model to determine the simulation scenarios [47]. First, the parameters for regression analysis are normalized. The seven key factors extracted from the sensitivity analysis are standardized according to the upper and lower limits of the parameters set in Table 1.Define the normalized variables of the original parametersμ, R, k1x, gt, la, s, and Δh as xi, i=1,2,3,4,5,6,7, it is(32)xi=Xi−x¯iui−x¯i,where, Xi is the original variable of the normalized variable xi, x¯i=1/2ui−li, ui and li are the upper and lower limits of xi, respectively.After that, the simulation scenarios are designed by the central combination design method (CCD) [45]. This method reduces the total number of observations to be evaluated while maintaining a certain accuracy of the second-order response surface model. Following this method, the number of observations per regression model is designed as 152 for the two wear regimes of mild and severe. The predictive results are described in the next section.


## 3.3.2. Comparative Analysis of RUL Prediction

Polynomial fitting and training are performed on all simulation scenes and response values in different wear regimes, and an RSM-based model of quadratic form as given in equation (20) is fitted. Figure 5 shows the prediction results for the three regression prediction models.Figure 5
Comparison of predicted and observed values from numerical simulation. (a) Wear regime I. (b) Wear regime II.
(a)(b)Figure5(a) shows that the prediction results (the rail service life in wear regime I) of the three prediction models are close to the simulation values and show good prediction stability over time. By contrast, the prediction results of rail service life for wear regime II in Figure 5(b) show that the predicted values based on RSM and SVM deviate significantly from the simulated values over time, confirming the lack of long-term trend prediction accuracy. Furthermore, these three prediction models based on RSM, SVM, and RVM can effectively predict the rail service life due to side wear, which can preliminarily prove the effectiveness of the three methods proposed in this paper. However, it can also be seen from the figure that the prediction results based on the RSM and SVM models deviate from the true value in a relatively large range. In comparison, the prediction results based on the RVM model are closer to the true value with less fluctuation. This shows that the prediction method based on the RVM model has less uncertainty and the prediction effect is more stable.To compare the prediction performance of the above three different prediction models in more detail, Table2 gives a numerical comparison of the prediction performance indicators.Table 2
Prediction performance indexes of different prediction models based on simulation data.
ModelNRSMSVsRVsR2RMSEMAEMRE(a) Wear regime IRSM152——0.90110.067640.58210.8314%SVM—133—0.93040.035720.30315.7663%RVM——130.98860.013430.11322.1096%(b) Wear regime IIRSM152——0.91420.08300.065112.953%SVM—135—0.94910.06190.048110.1997%RVM——140.98820.02720.01993.3516%Due to space limitations, Table2 (a) is only used as an example for explanation. The RMSE value of RVM is 0.1343, which is significantly lower than that of RSM and SVM, indicating that the prediction accuracy based on the RVM prediction model is higher. The same conclusion can be obtained by comparing the other two prediction error evaluation indexes. The coefficients of determination R2 of these models are 0.9011, 0.9304, and 0.9886, respectively. Although the R2 of SVM and RVM are slightly larger, those of the three models are close to 1, which shows that the fitting degree of these prediction models is high, and the RVM model has the highest fitting degree. Additionally, compared with the number of RSM training samples of 152, the number of support vectors (SVs) of the multivariate SVM is 133, accounting for almost 87.5% of the total number of training samples (152). The number of correlation vectors (RVs) of the multivariate RVM is only 13, accounting for only 8.6% of the total number of training samples, which shows that RVM has stronger sparsity in the application of regression prediction.As explained before, the RVM prediction model has the ability to express uncertainty. The upper and lower bounds of the prediction error constitute the 95% confidence interval for the prediction value. For the prediction of rail service life in wear regime I, the number of simulation values falling into the RVM prediction confidence interval is 140, accounting for 92.1% of the total number of training samples. In comparison, 137 simulation values of rail service life in wear regime II fall into the RVM prediction confidence interval, accounting for 90.13% of the total training samples. This verifies the prediction accuracy of the multivariable RVM prediction model.In summary, all the above conclusions are obtained by simulation tests. In the next section, the rail wear measured from the Guangzhou Metro is introduced to verify these conclusions.


## 4. Field Analysis

In this paper, the proposed prediction models are validated with rail wear measurements of the Guangzhou Metro Line 1 since January 2014. The response values used in this paper are first calculated on the basis of the prediction models with the specified input parameters. After that, a comparative analysis is performed between the response values obtained from the rail wear measurements and those calculated values. If the difference between the verification data and the fitted values of these prediction models is small enough, the obtained predictive model is acceptable.Figure6 shows the development of SRW measured at the midpoints of circular curves. It shows that the development laws of SRW on the different curves are similar. SRW does not occur on a new rail immediately. The wear of the gauge corner is aggravated with the increase in the NVP, while the SRW occurs gradually. With the increase in NVP, the SRW of different curves increases linearly at different rates. Then, the SWR of each curve is aggravated suddenly when the wear reaches 6 mm–8 mm. After that, it shows a linear increase with the NVP again after the wear reaches 9 mm–10 mm. This is consistent with the development process of the SRW predicted above. In addition, the time when SRW occurs in different curves is also different. In the steady wear stage (stage II), although the side wear of different curves increases linearly with the NVP, the growth rate is completely different. This shows that the SWR in each stage is affected by different track parameters and operating conditions.Figure 6
Development process of SRW in different railway curves.These evaluation indexes of prediction performance are calculated for validation purposes, as shown in Table3. The coefficients of determination (R2) of the predicted points for wear regimes I and II are very close to 1, proving the high fitting degree of these prediction models. The RMSE values of these predictive points for rail wear regime I are equal to 0.04214, 0.03938, and 0.0206, respectively, which proves their high prediction accuracy.Table 3
Prediction performance indexes of different prediction models based on field measurements.
ModelWear regime IWear regime IIR2RMSEMAEMRER2RMSEMAEMRERSM0.97850.042140.34449.716%0.97040.05070.048317.5733%SVM0.98130.039380.32468.0695%0.97550.04620.03658.9506%RVM0.99490.02060.1171.9827%0.98560.03540.0275.7753%Table4 gives the confidence interval obtained from the RVM prediction. It shows that the on-site observation values of rail service life for wear regime I all fall within the confidence interval obtained by the RVM prediction. For wear regime II, all the response values observed fall into the confidence interval given by the RVM model except for the set of data displayed in bold. From that, the better prediction performance of the RVM model is validated. From these given plots and values, it can be concluded that these predictive models proposed in this paper provide results with considerable accuracy. Among them, the prediction performance of the RVM model is the best, and it maintains a higher prediction accuracy while greatly reducing the computational effort.Table 4
Confidence interval obtained from RVM prediction.
CurvesT1104trainsT2104trainsMeasurementsRVM fitted valuesLower limit of RVM fittedUpper limit of RVM fittedMeasurementsRVM fitted valuesLower limit of RVM fittedUpper limit of RVM fitted10.620.618460.604300.65285.752645.4125.3355.84020.70.69100.67330.74087.843147.25446.78848.187332.64562.61742.47162.734117.8518319.29117.31821.38442.9442.832.613373.031221.54419.89919.30522.47254.87284.88114.43595.277033.003332.15228.73933.26365.16885.1805.0755.25219.9534421.70020.82622.12578.16647.9147.6418.61227.95929.25527.89132.71788.6729.1838.2959.58629.835929.55927.37231.809


## 5. Conclusions

Three regression prediction models based on RSM, SVM, and RVM proposed in this paper are demonstrated to enable RUL prediction of railway tracks due to rail wear. In these proposed models, seven measurable parameters, including track parameters and axle load, as well as the primary suspension stiffness related to the vehicle type, are taken into account. These parameters are proposed from the side wear mechanism and are the top seven important factors obtained by means of the global sensitivity analysis. The relationship between these parameters and the rail side wear is established by using three different regression models, where the rail side wear for various scenarios is calculated through a numerical simulation model of the vehicle-unballasted track and Specht friction wear. These predictive models provide results with coefficients of determination higher than 0.9, proving the considerable accuracy of these models. Among them, the RVM model maintains a high prediction accuracy while greatly reducing the computational effort. Finally, along with numerical validation, as presented in this work, validation through field measurements is also performed.



---

*Source: 7037655-2022-04-12.xml*