Railway ballast damage detection by Markov chain Monte Carlo-based Bayesian method

Abstract

This article reports the development of a Bayesian method for assessing the damage status of railway ballast under a concrete sleeper based on vibration data of the in situ sleeper. One of the important contributions of the proposed method is to describe the variation of stiffness distribution of ballast using Lagrange polynomial, for which the order of the polynomial is decided by the Bayesian approach. The probability of various orders of polynomial conditional on a given set of measured vibration data is calculated. The order of polynomial with the highest probability is selected as the most plausible order and used for updating the ballast stiffness distribution. Due to the uncertain nature of railway ballast, the corresponding model updating problem is usually unidentifiable. To ensure the applicability of the proposed method even in unidentifiable cases, a computational efficient Markov chain Monte Carlo–based Bayesian method was employed in the proposed method for generating a set of samples in the important region of parameter space to approximate the posterior (updated) probability density function of ballast stiffness. The proposed ballast damage detection method was verified with roving hammer test data from a segment of full-scale ballasted track. The experimental verification results positively show the potential of the proposed method in ballast damage detection.

Keywords

Damage detection Bayesian approach Markov chain Monte Carlo polynomial modelling vibration railway ballast

Introduction

Vibration-based damage detection of sleeper–ballast system

Railway ballast, which retains the sleepers and rails in their required positions against forces from different directions, is a crucial component of the ballasted railway track system (see Figure 7 for a segment of ballasted track). Damage to the ballast causes uneven support of the sleepers, which increases the likelihood of track settlement and changes in the geometry of the track. These changes may result in rail buckling and derailment. Chinese Railways (CR) statistics have revealed that about 75% of the daily maintenance work on track structures is performed due to railway ballast deterioration.¹ Regular visual inspection is common practice of track maintenance adapted by industry. However, it is difficult to visually identify the ballast damage under the concrete sleeper without lifting up the sleeper. The literature reviews regarding non-model based methods of detecting ballast damage are available in Lam et al.^2,3 and are not repeated here. It is clear from the reviews that the existing methods to detect the ballast damage only focus on the ‘average’ stiffness of a given layer of substructure materials, but not the distribution of stiffness within a layer and it is difficult to locate the ballast damage. Therefore, there is a great demand of a practical method for non-destructive evaluation of ballast under a sleeper.

Theoretically, degradation of the ballast reduces the stiffness of the ballast in supporting the sleeper. When the stiffness distribution of the in situ sleeper is altered, its vibration characteristics change. This change can be reflected by variations in the modal parameters of the in situ sleeper. Therefore, it is possible to calculate the reduction in the ballast stiffness on the basis of the natural frequencies and the mode shapes of the in situ sleeper. The detection of damage by solving an inverse problem (i.e. the calculation of a set of model parameters of a system based on a set of measured responses) is not new for civil engineering structures.⁴ For example, the detection of cracks on beam and plate types of structures,^5,6 the assessment of operational condition of bridges⁷ and the structural damage detection of a benchmark through the vibration data.^8,9 Comprehensive reviews in vibration-based damage detection for civil engineering structures are available in references.^10,11 However, the extension of this idea to the detection of damage for railway ballast is new and challenging.

In 2012, Lam et al.² proposed the use of the measured modal parameters of an in situ concrete sleeper for ballast damage detection following a deterministic approach. The results of the numerical and experimental verification were very encouraging and showed that the modal parameters of in situ sleepers were sensitive to ballast damage. One of the major barriers in applying this method² in real-life applications is the high degree of uncertainty associated with the mechanical properties of the rail–sleeper–ballast system (especially for the layers of irregular ballast particles), and the uncertainties induced from the modelling error and measurement noise. To quantify this uncertainty, a probabilistic method was developed in Lam et al.³ for determining the posterior probability density function (PDF) of the ballast stiffness using the Bayesian statistical system identification framework.

Bayesian inference can update the uncertainties in the values of the model parameters^12,13 and provide a feasible identification solution for the purpose of structural health monitoring and crack growth prediction.^14,15 The original formulation of the Bayesian framework based on the measured time-domain data was developed by Beck and Katafygiotis¹⁶ in 1998. In 2000, it was extended to use measured modal data.¹⁷ Recently, Au et al.¹⁸ combined the Bayesian modal identification method and the Bayesian model updating method for the model updating of a coupled-slab system based on a set of ambient field test data.¹⁹ In this article, the Bayesian statistical system identification framework based on modal data is extended for the purpose of identification of the railway ballast damage.

Challenges in vibration-based damage detection of sleeper–ballast system

Modelling of ballast stiffness distribution

Previously, the sleeper–ballast system was modelled as a Timoshenko beam on an elastic foundation with several discrete ballast regions and the ballast stiffness in a given region was assumed to be a constant^2,3 (see Figure 1 for the discrete ballast model). Because the variation in the ballast stiffness is continuous along the sleeper in real situations, the discontinuous ‘jumps’ in the ballast stiffness at the interfaces between two discrete regions are ‘artificial’. This artificial modelling error induces the discrepancy between the model-predicted and measured modal parameters and reduces the reliability of the results. In this article, a newly developed polynomial modelling method is proposed to avoid the problem of discontinuous jumps.

Unidentifiability in model updating

Model updating is a popular method to identify the most probable set of model parameters conditional on a set of measured data under an assumed class of models.^20,21 Following the Bayesian framework, the model updating problem is considered to be locally identifiable when the posterior PDF of the uncertain parameters can be approximated by a weighted sum of multivariable Gaussian distributions.¹⁶ Depending on the quantity and quality of the measurements and the complexity of the class of models, model updating problems may not necessarily be identifiable. The definition of unidentifiable cases and the transition from identifiable to unidentifiable model updating problems are extremely complicated, and they were comprehensively studied by Katafygiotis and Lam²² in 2002. It must be pointed out that the Bayesian ballast damage detection method proposed in Lam et al.³ is applicable only for locally identifiable model updating problems. In general, a model updating problem may become unidentifiable²² when the model class has great complexity (when other factors are unchanged) and a limited measurement data is available. The objectives of this article are to overcome these limitations and develop a new and more general ballast damage detection method that is applicable even when the model updating problem is not locally identifiable, that is, handling the conditions, such as (1) model class with great complexity, (2) limited measured data and (3) high level of measurement noise. To achieve this goal, the enhanced Markov chain Monte Carlo (MCMC) simulation²³ is employed and extended in this study for model updating and model class selection of the rail–sleeper–ballast system.

MCMC methods have been popular for sampling from a complicated probability distribution based on the construction of a Markov chain whose limit distribution is the target probability distribution.^24,25 These methods have recently been adopted for Bayesian model updating.²⁶ It has been demonstrated that samples generated by MCMC can be used to approximate the posterior PDF of uncertain parameters regardless of the identifiability of the model updating problem. In 2007, Ching and Chen²⁷ derived the transitional Markov chain Monte Carlo (TMCMC) method for Bayesian model updating and model class selection with particular consideration of complicated PDFs that have multiple modes. The TMCMC method was numerically verified with shear building data and its computational efficiency was demonstrated. To efficiently explore the parameter space, the user must appropriately determine the number of sampling levels; this issue is still a challenging research problem for the application of MCMC methods in model updating, especially for the purposes of damage detection. Recently, an enhanced MCMC method that focuses on the calculation of an appropriate number of sampling levels was developed and successfully verified through field test data from a coupled-slab system.²³ One of the main contributions of this study is to extend the formulation of the MCMC method in Lam et al.²³ for Bayesian model class selection, which forms an important element in the proposed ballast damage detection method. Following the MCMC-based Bayesian method, the railway ballast stiffness distribution under the sleeper can be estimated together with the corresponding posterior uncertainties, and the ballast damage detection results can be well exhibited using the graphical presentation, which can directly display valuable information about the uncertainties to inspectors and engineers during their visual inspection.

Proposed ballast damage detection method

The proposed method integrates the newly developed polynomial modelling method with the MCMC-based Bayesian model updating and model class selection methods. Due to the modelling error induced by the discrete modelling method,³ the polynomial modelling method is developed in this article to capture the ballast stiffness distribution using Lagrange polynomials.²⁸ In this study, the order of polynomial to be employed in modelling the ballast stiffness distribution is calculated by extending the Bayesian model class selection method in Beck and Yuen.²⁹ According to past experience, some possible model classes are relatively complex and the corresponding model updating problems are not locally identifiable. Therefore, the model class selection formulation in equation (3) of Beck and Yuen,²⁹ which works in globally identifiable cases, cannot be directly employed in this case. In the proposed method, the formulation in Beck and Yuen²⁹ was extended such that the probability of a model class can be approximated using the set of samples generated by MCMC. Once the most plausible class of models is selected, the posterior PDF calculated by the MCMC samples can be used to quantify the damage status of the ballast under the sleeper.

Polynomial modelling of ballast stiffness distribution

As verified by vibration data from a full-scale test panel, the in situ sleeper can be modelled as a Timoshenko beam located on an elastic foundation and the two rails are modelled as two independent masses on the beam.³Figure 1(a) shows an example that divides the ballast under the sleeper into six regions. Assuming the nominal value of ballast stiffness to be k_b, the ballast stiffness in a given region is equal to the non-dimensional factor of that region multiplied by the nominal value. Figure 1(b) shows an example of the updated ballast stiffness distribution. The discrete modelling method is simple, but it may not be able to capture the ‘real’ ballast stiffness distribution because the ballast stiffness in a given region is a constant, resulting in a jump discontinuity in ballast stiffness at the interface between two regions (see Figure 1(b)). This jump discontinuity does not exist in reality and is a source of a modelling error that will affect the performance of the proposed method in ballast damage detection.

To overcome the jump discontinuity problem, the ballast stiffness distribution under the sleeper should be described by a class of continuous functions. One simple and logical choice of continuous functions for this purpose is the family of polynomials with different orders.³⁰Figure 2(a) illustrates the idea by approximating the ballast stiffness distribution by a polynomial for a given order, say, the third-order polynomial (i.e. θ₁x³ + θ₂x² + θ₃x + θ₄ = 0). Because a polynomial is always a smooth function, the problem of jump discontinuity is avoided. The most straightforward way to use a polynomial as a class of models in model updating is to consider the polynomial coefficients (i.e. θ₁ to θ₄ in this illustrative example) as uncertain parameters (as in Hu and Lam³⁰). However, unlike the discrete modelling method (for which the uncertain parameters are the ballast stiffness values), the polynomial coefficients are nonparametric. The identified polynomial must be plotted before one can understand the distribution of the ballast stiffness. Furthermore, the MCMC-based model updating method calculates the posterior PDF of the polynomial coefficients, which cannot provide engineers and researchers with any direct insight into the posterior uncertainties of the ballast stiffness distribution. An additional step is required to map this information to the posterior uncertainties of the ballast stiffness distribution. To overcome this difficulty, the Lagrange polynomial is used in the polynomial modelling method to interpolate a polynomial by the ballast stiffness values at a discrete number of points under the sleeper.

The concept of polynomial modelling by Lagrange polynomial is illustrated by an example of a third-order polynomial in Figure 2(b). Four points are used to represent a third-order Lagrange polynomial. The most logical way to locate the four points under the sleeper is to distribute them uniformly (i.e. to divide the sleeper into three equal parts as shown in Figure 2(b)). The nominal value of the ballast stiffness is k_b and a non-dimensional factor is assigned to the ballast stiffness at each point. The product of the nominal value and θ₁ to θ₄ are the values of ballast stiffness at the four points. This concept can be easily extended to a general case of a polynomial of order N. In this case, the sleeper is divided into N equal parts by N + 1 points and the ballast stiffness factors at these N + 1 points are the uncertain parameters for model updating.

Assuming that the length of the sleeper is equal to l, the coordinates of these N + 1 points are (0, θ₁), (l/N, θ₂), (2l/N, θ₃),…, ((j − 1)l/N, θ_j),…, (l, θ_N_+ 1) for j = 1,…, N + 1. With this point arrangement rule, all data points are different and the Lagrange polynomial, L(x), can be formulated as a linear combination as²⁸

L (x) = \sum_{j = 1}^{N + 1} θ_{j} l_{j} (x)

(1)

where

\begin{matrix} l_{j} (x) = \underset{\binom{1 \leq m \leq N + 1}{m \neq j}}{Π} \frac{x - (m - 1) \frac{l}{N}}{(j - m) \frac{l}{N}} \\ = (\frac{x - 0}{(j - 1) \frac{l}{N}}) \dots (\frac{x - (j - 2) \frac{l}{N}}{\frac{l}{N}}) \\ (\frac{x - j \frac{l}{N}}{- \frac{l}{N}}) \dots (\frac{x - l}{(j - N - 1) \frac{l}{N}}) \end{matrix}

Finally, the $l_{j} (x)$ function can be simplified as

l_{j} (x) = \frac{(x - 0) \dots (x - (j - 2) \frac{l}{N}) (x - j \frac{l}{N}) \dots (x - l)}{{(- 1)}^{N + 1 - j} (j - 1)! (N + 1 - j)! {(\frac{l}{N})}^{N}}, for j = 1, \dots, N + 1

where θ₁, θ₂,…, θ_N +₁ are the uncertain parameters to be updated through the proposed MCMC-based Bayesian model updating method. Equation (1) is used to formulate the Lagrange polynomial for given ballast stiffness values at discrete points.

In general, a polynomial of higher order has a greater ability to capture more complicated variations (or fluctuations) in the ballast stiffness along the sleeper. Therefore, the selection of an appropriate order (or the most plausible class of models) is essential for the success of ballast damage detection.

The concept of model class in the polynomial modelling method based on Lagrange polynomial is discussed before introducing the model class selection method. Figure 3 shows a series of possible classes of models for capturing the ballast stiffness distribution under the sleeper with increasing model complexity (the simplest class at the top). The model class M ₁ (see Figure 3(a)) is a zero-order polynomial representing an undamaged situation with a uniform distribution of ballast stiffness, namely, the ballast stiffness is a constant, θ₁k_b. With damage to the ballast, the ballast stiffness distribution will not be uniform. The model class M ₂ (see Figure 3(b)) corresponds to a first-order polynomial (i.e. a straight line) with two uncertain parameters (i.e. θ₁ and θ₂), where θ₁k_b and θ₂k_b are the ballast stiffness values at the two end points of the sleeper. In Figure 3(c), the model class M ₃ corresponds to a second-order polynomial that consists of three uncertain parameters (i.e. θ₁, θ₂ and θ₃) that represent the ballast stiffness of the left end, middle and right end of the sleeper, respectively. This idea can be extended to a general case of model class M _N, which corresponds to a polynomial of order N − 1, and N points are needed to parameterize this Lagrange polynomial (see Figure 3(d)). The uncertain parameters are the N scaling factors (θ₁, θ₂,…, θ_N).

Aside from the scaling factors for the ballast stiffness, several additional uncertain parameters are considered to further reduce the level of modelling error of the class of models. A non-dimensional scaling factor θ_E is used to scale the nominal value of Young’s modulus of the concrete sleeper. This is mainly due to the prestress of the concrete sleeper and also the authors’ lack of access to the material properties of the concrete sleeper, so the nominal value is the best estimation available. The two additional uncertain parameters are the left and right rail masses, θ_L and θ_R. As a result, a total of N + 3 uncertain model parameters are considered for the model class M _N in the model updating process. They are grouped into an uncertain model parameter vector θ_N = [θ₁, θ₂,…, θ_N, θ_E, θ_L, θ_R].

The finite element method is used to model the rail–sleeper–ballast system in this study. The formulation of the element stiffness matrix of a Timoshenko beam element on an elastic foundation can be found in Lam et al.³ and Krenk.³¹ Because of the thickness of the sleeper, the effect of shear deformation cannot be neglected, especially in the vibrations of higher modes,² and therefore, the Timoshenko beam formulation is used instead of the relatively easy Euler beam formulation. The consistent mass matrix³² is used in the dynamic analysis. The element stiffness and mass matrices are used to assemble the system stiffness and mass matrices. The system’s modal parameters, such as the natural frequencies and mode shapes, can be calculated by solving the corresponding eigenvalue problem.

MCMC-based Bayesian model updating

For a given model class M _N, a set of uncertain model parameters θ_N was defined in the previous section. To simplify the formulation, the subscript N is dropped in this section. The main purpose of model updating following the Bayesian approach is to calculate the joint posterior (or updated) PDF of θ with a set of measured modal parameters D by extending the method in Vanik et al.¹⁷ Instead of considering the sample variance of the measured data as the variance of the prediction error, the effects of both the modelling error and the measurement noise are considered through the measure of percentage error, κ, which is treated as an uncertain parameter in the model updating process. By assuming that the modelling error (in percent) is the same for the natural frequencies and the mode shapes, the formulation of the posterior PDF of the uncertain parameter θ conditional on the set of measurement is given by²³

p (θ | D) = c \exp (- \frac{1}{2 κ^{2}} J (θ))

(2)

where D denotes the set of measured natural frequencies and mode shapes, and c is a normalizing constant. The goodness-of-fit function J³³ is given by

J (θ) = \sum_{a = 1}^{r} [{(\frac{{\hat{f}}_{a} - f_{a} (θ)}{{\hat{f}}_{a}})}^{2} + (1 - {| {\hat{ψ}}_{a}^{T} ψ_{a} (θ) |}^{2})]

(3)

where the subscript a denotes the mode index, r is the total number of modes to be considered in the model updating process, ${\hat{f}}_{a}$ denotes the measured natural frequency (Hz), $f_{a} (θ)$ is the calculated natural frequency whose dependence on the uncertain model parameters θ is emphasized, and ${\hat{ψ}}_{a}$ and $ψ_{a} (θ)$ denote the measured and calculated mode shapes, respectively. Interested readers are directed to Lam et al.²³ for the detailed derivation.

When the model updating problem is globally or locally identifiable, the posterior PDF in equation (2) can be approximated by a multivariable Gaussian PDF.¹⁶ However, for the purpose of model class selection, model updating for different model classes ranging from simple to complex is necessary for finding an appropriate model class such that it can well fit the measured data, without over-fitting the noise and/or error parts of the measurements. For a given set of measurements, the model updating problem may become unidentifiable when the class of models is too complex. In this situation, the posterior PDF becomes very complicated (as shown in Katafygiotis et al.³⁴) and cannot be accurately approximated by a multivariable Gaussian distribution.

A newly developed MCMC-based Bayesian model updating method was proposed by the authors in Lam et al.²³ and Beck and Au²⁶ to approximate the posterior PDF and calculate posterior uncertainties by generating samples in the important regions (i.e. the regions in the parameter space with high probability). The MCMC method in Lam et al.²³ is used in this article for model updating. The main idea of generating samples by MCMC is to perform the sampling at multiple levels with a series of bridge PDFs to ensure that the samples converge to the important regions. In the proposed method, the bridge PDFs are constructed according to equation (2) at each level in such a manner that the covered regions decrease gradually and move towards the important region of the posterior PDF.

For the completeness of this article, the MCMC algorithm for model updating is briefly summarized in the flowchart in Figure 4. In the figure, g_r denotes the number of required sampling levels. A is an algorithm parameter that controls the change rate of the bridge PDFs between two successive levels. To sample with the MCMC algorithm, in level 1, N_ss samples of the uncertain parameter θ are generated with the Metropolis–Hastings (MH) algorithm^35,36 with the bridge PDF of level 1, p₁, as the target PDF and the prior PDF as the instrumental PDF, which is used in the MH algorithm to draw samples. It should be set appropriately so the samples can be generated conveniently. In a general level g, the sampling procedure is similar to that in level 1 except that the instrumental PDF in the MH algorithm is constructed by the kernel density k_g(θ)³⁷ instead of the prior PDF. After g_r levels of sampling, the posterior PDF and marginal PDFs of the uncertain parameters can be easily obtained by the kernel density estimation³⁷ based on the samples from the final level without carrying out the high-dimensional numerical integration. The posterior marginal PDF of each component of θ (i.e. k_gr(θ(i)) in equation (4)) can be obtained conveniently by analytical integrating the posterior PDF k_gr(θ). The posterior uncertainties can then be assessed by the posterior marginal PDFs. A general form of the posterior marginal PDF of each component of θ can be expressed as²³

k_{gr} (θ (i)) = \frac{1}{N_{ss}} \sum_{j = 1}^{N_{ss}} W^{(j)} N (θ^{j} (i), σ^{(j)})

(4)

The posterior marginal PDF is the weighted sum of Gaussian PDFs $N (θ^{j} (i), σ^{(j)})$ , where W^(j) is the weighting of the jth sample. Each of these Gaussian PDFs centres at one sample with the standard deviation calculated by the generated samples. Interested readers are directed to Lam et al.²³ and Beck and Au²⁶ for the detailed formulation of the algorithm.

MCMC-based Bayesian model class selection

The order of polynomials used has a significant effect on the performance of model updating. As illustrated in Figure 3, a zero-order polynomial can only capture the ballast stiffness distribution in the undamaged situation. The higher the order, the more complicated the ballast stiffness distribution is captured. Instead of deciding the order in an arbitrary manner or following ad hoc rule-of-thumbs, the probability of a list of possible polynomial orders (i.e. model classes) conditional on a set of measurements is calculated by following the Bayesian approach.²⁹ The model class with the highest probability is considered to be the most plausible model class and is used for the purpose of ballast damage detection.

Considering the situation with N_m possible model classes, the probability of each model class conditional on the set of measurements is, p( M _j| D ) for j = 1 to N_m. This probability is proportional to the evidence of the model class, which can be expressed as²⁹

p (D | M_{j}) = \int p (D | θ_{j}, M_{j}) p (θ_{j} | M_{j}) d θ_{j}

(5)

where θ_j represents the model parameters vector in the model class M _j, p( D |θ_j, M _j) is the likelihood function, and p(θ_j| M _j) is the prior PDF. Evaluating the evidence in equation (5) is difficult because the high-dimension integration is almost impossible to calculate analytically and is too computationally demanding to calculate numerically. The asymptotic approach in Beck and Yuen²⁹ can be used to approximate the evidence when the model updating problem is locally identifiable. Because model classes with different complexity need to be considered in the selection process, it is very likely that some of the possible model classes are so complex that the model updating problem becomes unidentifiable. As a result, the formulation given in Beck and Yuen²⁹ is not applicable. In this study, the integral in equation (5) is approximated using samples generated by MCMC following the algorithm presented in section ‘MCMC-based Bayesian model updating’.

To implement this idea, equation (5) is rewritten according to the concept of importance sampling,³⁸ as

p (D | M_{j}) = \int (\frac{p (D | θ_{j}, M_{j}) p (θ_{j} | M_{j})}{p (θ_{j} | D, M_{j})}) p (θ_{j} | D, M_{j}) d θ_{j}

(6)

where the posterior PDF is introduced such that the evidence can be approximated by the following summation

p (D | M_{j}) \approx \frac{1}{N_{ss}} \sum_{h = 1}^{N_{ss}} \frac{p (D | θ_{j}^{(h)}, M_{j}) p (θ_{j}^{(h)} | M_{j})}{p (θ_{j}^{(h)} | D, M_{j})}

(7)

where the superscript h denotes the sample index. When equation (7) is directly implemented, the PDF values are evaluated at discrete samples. However, the numerical value of the goodness-of-fit function J(θ) in equation (3) is usually very large and the numerical value of the likelihood p( D |θ_j^(h), M _j) in equation (7) having a term exp(−1/(2κ²) J(θ)) becomes extremely small, resulting in a computational problem. The formulation given in equation (7) is modified to avoid this computational problem. According to the Bayesian theorem, the evidence of the model class M _j in model class selection can be expressed as²⁹

p (D | M_{j}) = \frac{p (D | θ_{j}, M_{j}) p (θ_{j} | M_{j})}{p (θ_{j} | D, M_{j})}

(8)

The right-hand side of equation (8) is part of the expression in the integral in equation (6). Since this part of the expression is independent of θ_j (as shown in equation (8)), it can be moved outside the integral. Equation (6) thus becomes

p (D | M_{j}) = (\frac{p (D | θ_{j}, M_{j}) p (θ_{j} | M_{j})}{p (θ_{j} | D, M_{j})}) \int p (θ_{j} | D, M_{j}) d θ_{j}

(9)

In equation (9), the integration of the posterior PDF over the parameter space is equal to unity by definition; therefore, equations (8) and (9) are equivalent. Taking the natural logarithm of equation (8) and making use of equation (9) give

\begin{matrix} \ln (p (D | M_{j})) = \int \ln (\frac{p (D | θ_{j}, M_{j}) p (θ_{j} | M_{j})}{p (θ_{j} | D, M_{j})}) p (θ_{j} | D, M_{j}) d θ_{j} \\ = \int [\ln (p (D | θ_{j}, M_{j})) + \ln (p (θ_{j} | M_{j}))] p (θ_{j} | D, M_{j}) d θ_{j} \\ - \int [\ln (p (θ_{j} | D, M_{j}))] p (θ_{j} | D, M_{j}) d θ_{j} \end{matrix}

(10)

Approximating equation (10) using the MCMC samples gives

\begin{matrix} \ln (p (D | M_{j})) \approx \frac{1}{N_{ss}} \sum_{h = 1}^{N_{ss}} [\ln (p (D | θ_{j}^{(h)}, M_{j})) + \ln (p (θ_{j}^{(h)} | M_{j}))] \\ - \frac{1}{N_{ss}} \sum_{h = 1}^{N_{ss}} [\ln (p (θ_{j}^{(h)} | D, M_{j}))] \end{matrix}

(11)

Instead of calculating the evidence, the logarithm of evidence is computed by equation (11) to avoid the computational problem. Readers are referred to Ching and Chen,²⁷ Cheung and Beck³⁹ and Muto and Beck⁴⁰ for more details about the Bayesian model class selection using stochastic simulation.

In general, a model class with greater complexity (i.e. with more uncertain model parameters) can have better fitting to the measurement, resulting in a higher likelihood value. According to equation (11), the logarithm of the evidence of a model class is not solely determined by the likelihood. It must be pointed out that the most plausible class of models cannot be obtained by the likelihood alone. In this study, the class of models to be selected is the one with the highest value of evidence among all the model classes considered. The selected model class, on one hand, is complex enough to provide a ‘good fit’ to the measurement data and, on the other hand, is simple enough that it will not ‘over fit’ the data.

MCMC-based ballast damage detection

In this section, the polynomial modelling method, MCMC-based Bayesian model updating and model class selection methods are integrated to form the proposed method of ballast damage detection, which is schematically illustrated in the flowchart in Figure 5. A set of vibration data (i.e. acceleration responses at a given degree-of-freedom under impact excitations at various degrees of freedom) is obtained from the target in situ sleeper (see process 1 in Figure 5). In the proposed method, inspectors or permanent way engineers can achieve this by a simple roving hammer test. The MODE-ID method⁴¹ is used to identify the modal parameters, such as the natural frequencies and mode shapes, from the measured time-domain responses (see process 2 in Figure 5). The set of measurements D is then ready for model updating by the MCMC method.

With different orders of polynomial, N_m model classes of the rail–sleeper–ballast system are formed (as discussed in section ‘Polynomial modelling of ballast stiffness distribution’). For each model class, the set of MCMC samples is generated on the basis of the set of measurements D by following the algorithm as shown in Figure 4 (see process 3 in Figure 5). Using the new MCMC-based model class selection formulation in equation (11), the logarithm of the evidence of all model classes can be calculated (see process 4 in Figure 5). The most plausible class of models, M _N, can be selected based on the evidence value of the model class. In this study, the model class with the maximum value of the evidence is selected as the most plausible model class. The posterior marginal PDF of the ballast stiffness for the selected model class can be calculated by equation (4) in the MCMC-based Bayesian model updating method (as discussed in section ‘MCMC-based Bayesian model updating’). Instead of focusing on the joint posterior PDF of all uncertain model parameters, the marginal posterior PDFs of the ballast stiffness parameters are important from the viewpoint of ballast damage detection. The primary outcomes from the proposed MCMC-based method are the marginal posterior PDF of various model parameters, especially the ballast stiffness value at the predefined discrete points.

According to the authors’ experience, in most cases, the posterior marginal PDFs of the ballast stiffness are similar to but not exactly the same as Gaussian distributions. The mean and standard deviation of each ballast stiffness parameter can be calculated and used to generate the ballast stiffness distribution using Lagrange polynomial.

A schematic figure is shown in Figure 6 to illustrate the idea of ballast damage detection following the proposed MCMC-based method. The solid line shows the posterior mean of the ballast stiffness distribution under the sleeper and the two dashed lines represent the upper and lower boundaries of one standard deviation from the mean. It is clear that the ballast under the left side of the sleeper is damaged, with a reduction in stiffness of approximately 40%. The upper and lower boundaries show the relatively small uncertainty of the identified ballast stiffness distribution, and hence, it reflects the great reliability of the ballast damage result.

Experimental verification using full-scale tests

To illustrate the proposed procedures and verify the proposed MCMC-based Bayesian method for detecting ballast damage, a full-scale indoor test panel consisting of eight in situ sleepers was built in the basement of a factory building in Kowloon Bay, Hong Kong (see Figure 7). The sleepers were labelled 1 to 8 from the bottom to the top of the figure. The ballasted track in the test panel was constructed according to Hong Kong MTR specifications.⁴² Two UIC60 rails (with a mass density of 60 kg/m) were fastened to the concrete sleepers. The centre-to-centre spacing of the sleepers was 700 mm and the thickness of the granite ballast layer was 350 mm. The ballast was made of crushed rock aggregate of about 50 mm (ballast of this size is considered normal under the Hong Kong MTR specifications).

Experiment description

In this study, ballast damage is defined as railway ballast degradation, which leads to a reduction in the ballast size, resulting in a decrease in the ballast’s packing level and stiffness in supporting the concrete sleeper. In the full-scale test panel, ballast damage was artificially simulated by replacing the normal-sized ballast (~50 mm) with smaller ballast (~15 mm). Sleepers 1–4 were undamaged and the artificial ballast damage was simulated under the middle of sleepers 5–8. To avoid boundary effects (the situation of the sleepers 1 and 4 is different from that of sleepers 2 and 3 (considered as undamaged), and the situation of the sleepers 5 and 8 is different from that of sleepers 6 and 7 (considered as damaged)), measured data from sleepers 1, 4, 5 and 8 were not used in model updating. In this article, data measured from sleepers 3 (undamaged) and 6 (damaged at the middle of the sleeper) were employed to demonstrate and verify the proposed method.

The equipment used in the roving hammer tests is listed in Figure 8. The target in situ sleeper was excited by an impact hammer that can be equipped with tips with different hardness (represented by different colours), as shown in Figure 8(c). A soft tip is suitable for exciting lower modes (e.g. modes 1 and 2), and a stiff tip is more appropriate for exciting higher modes (e.g. modes 3–5). One KISTELER-8776A50M3 accelerometer with sensitivity level about 100 mV/g as shown in Figure 8(a) was installed on the top surface of the concrete sleeper at a fixed location at the central line. The sensor was installed vertically to capture the system’s vertical vibration, and the impacts were applied at 11 considered locations (i.e. 11 measured degrees of freedom), as illustrated in Figure 9. The measured vibration responses from the sensor were transferred through the cables with special coating (see Figure 8(b)) and digitized using one NI-9234 module on a Compact DAQ-9178 chassis (see Figure 8(d)). For each sleeper, the entire measurement was divided into 11 parts. In each part, the excitation impact force was applied to a given measured degree of freedom, as shown in Figure 9. The measurement duration for each part was set to 20 s with a sampling frequency of 6400 Hz. During the measurement, the impact forces were applied one by one at an interval of about 1–2 s. As a result, more than 10 impulses were recorded in each part of the measurement.

Modal identification

A set of measured time-domain data is shown in Figure 10. The time-domain responses were transformed to the frequency domain via fast Fourier transformation. The modal parameters, such as the natural frequencies and mode shapes, were identified by MODE-ID.⁴² Due to the space limitation, detail procedures for obtaining a set of measured natural frequencies and mode shapes are not presented, and interested readers are redirect to Lam et al.^2,3 The first five modes were identified with high accuracy and were used in ballast damage detection; the natural frequencies of the first five modes for both undamaged and damaged cases were summarized in Table 1. It is clear from the table that the ballast damage in general reduces the natural frequencies of the in situ sleeper. The first two modes are very sensitive to the ballast damage (with over 30% reduction in natural frequency), while mode 3 can also be considered as sensitive (with about 10% reduction in natural frequency). In this case study, modes 4 and 5 are not sensitive to the ballast damage.

Table 1.

Natural frequencies for the undamaged (sleeper 3) and damaged (sleeper 6) cases (Hz).

	Mode 1	Mode 2	Mode 3	Mode 4	Mode 5
Undamaged	90.96	98.93	174.40	415.78	829.12
Damaged	59.05	68.61	156.08	408.56	786.57
Difference (%)	35.08	30.65	10.50	1.74	5.13

Lagrange polynomials were used in the polynomial modelling method. Following the MCMC-based model updating method, experimental verification was carried out on both the undamaged and damaged cases in the following sections. According to Lam et al.^2,3 and Kumaran et al.,⁴³ the nominal values of Young’s modulus and the density of the concrete sleeper were set to 38 × 10⁹ N/m² and 2200 kg/m³, respectively, and the nominal values of the ballast stiffness and rail mass were set to 120 × 10⁷ N/m² and 42 kg, respectively.

Undamaged case

The proposed MCMC–based method of ballast damage detection was first tested using the measurements from the undamaged system (i.e. sleeper 3). The class of zero-order polynomials (i.e. model class M ₁) was used and there were four uncertain parameters (θ_E, θ_L, θ_R and θ₁). The MCMC method was used for generating samples for this model class. Figure 11 shows the MCMC samples in the final sampling level (i.e. level 20). In each level, 2000 samples were generated, and they are presented as dots in various two-dimensional projections in Figure 11. The distribution of the samples shows the important range in the parameter space. It is clear from the figure that the important region for θ_E is from 1 to 1.6, θ_L and θ_R are in the range from 0 to 0.6 and the scaling factor of ballast stiffness θ₁ is from 0.8 to 1.4. It is clear that there is only one important region in the parameter space of interest.

The evidence of this model class was calculated by equation (11) using the MCMC samples. The value of the logarithm of the evidence is shown in the first row of Table 2. Next, a more complicated model class was tested. Similar procedures were carried out for model class M ₂ (i.e. the first-order polynomial). Because of space limitations, the generated samples are not shown. The value of the logarithm of the evidences of the model class M ₂ is shown in the second row of Table 2. Since the logarithm of the evidence of M ₁ was greater than that of M ₂, the model class M ₁ (i.e. zero-order Lagrange polynomial) was selected in this case.

Table 2.

Evidences for different classes of models in the undamaged and damaged cases.

	Class of models	Logarithm of evidence
Undamaged	M ₁	−11.34
	M ₂	−12.88
Damaged	M ₁	−13.50
	M ₂	−11.74
	M ₃	−11.44
	M ₄	−12.44

Once the most plausible model class ( M ₁) was selected, its MCMC samples were used to calculate the posterior marginal PDFs of all uncertain parameters (four in this case); they are presented in Figure 12. The shape of the posterior marginal PDFs was similar to a Gaussian distribution. The tails for both θ_E and θ₁ approached zero; however, the left tails for the PDFs of θ_L and θ_R were cut at zero. By fitting the marginal posterior PDFs with Gaussian distributions, the estimated means and standard deviations of uncertain model parameters were calculated and are summarized in the first two rows of Table 3. From Table 3, the optimal factor for ballast stiffness was close to unity (1.06) showing that the nominal value provides a good estimation. The corresponding standard deviation was 0.14, which shows that the uncertainty associated with the updated ballast stiffness in this case was not great. Because the class of zero-order polynomials was selected, there was only one uncertain ballast stiffness parameter, implying that the ballast stiffness distribution was uniform. The ballast damage detection results are summarized in Figure 13. The solid line represents the posterior mean of the ballast stiffness distribution and the two dashed lines represent the upper and lower boundaries at one standard deviation from the mean. Obviously, the flat distribution of the ballast stiffness implies that there was no ballast damage in this case and the damage detection results from the proposed method were correct.

Table 3.

Optimal parameters and the corresponding standard deviations in the undamaged and damaged cases.

	${\hat{θ}}_{1}$	${\hat{θ}}_{2}$	${\hat{θ}}_{3}$	${\hat{θ}}_{E}$	${\hat{θ}}_{L}$	${\hat{θ}}_{R}$
Undamaged case	1.06	–	–	1.30	0.10	0.07
Standard deviation	0.14	–	–	0.15	0.19	0.17
Damaged case	1.56	0.37	2.12	1.33	0.06	1.25
Standard deviation	0.31	0.17	0.61	0.12	0.17	0.91

Other uncertain model parameters identified in Table 3 are also discussed. The updated Young’s modulus of the sleeper was 30% greater than the nominal value. Since the effect of prestressed force in the sleeper was not considered in the modelling, the nominal value underestimated Young’s modulus of the sleeper. The updated rail masses on the left and right sides were both smaller than the nominal value due to the uneven levelling of the neighbouring sleepers.³

Damaged case: ballast damage at the middle under the sleeper

The measurements from the damaged system (i.e. sleeper 6) were used to test the proposed method. The model classes M ₁ to M ₄ were first considered; the logarithms of the evidences are summarized in rows 3–6 of Table 2. Based on the results in Table 2 for the first four models, it can be seen that model M ₃ is selected as the most probable model due to the highest evidence.

Based on the set of MCMC samples generated for model class M ₃, the posterior marginal PDFs of the uncertain parameters were calculated and are presented in Figure 14. The shapes of the PDFs were similar to a Gaussian distribution and the left tails for the PDFs of θ_L and θ_R were cut. The estimated means and standard deviations for all of the uncertain parameters were calculated and are summarized in rows 3 and 4 of Table 3.

Based on model class M ₃, the ballast stiffness at the left end, middle and right end of the sleeper was calculated through the scaling factors θ₁, θ₂ and θ₃, respectively. The ballast stiffness factors at the left and right ends were 1.56 and 2.12, respectively, and that at the middle was 0.37. Before the ballast stiffness distribution was calculated, the results already showed that the ballast at the middle of the sleeper was damaged, and that the ballast stiffness reduction at the middle was about 63% of the ballast stiffness. The ballast damage detection results of the proposed method are summarized in Figure 15, which shows the mean (solid line) and the upper and lower boundaries (dashed lines at one standard deviation from the mean) of the ballast stiffness distribution. It is very clear from the figure that the ballast damage was located at the middle part of the sleeper and the uncertainty of the identified ballast distribution was reasonable, which indicates that the proposed method successfully identified the simulated ballast damage.

To get an idea of the accuracy of the match between the measured and model-predicted (after model updating) modal parameters, the ‘most probable’ model (i.e. the mean values estimated in Table 3) of M ₃ was used to calculate the natural frequencies and mode shapes of the system. Figure 16 shows the measured and model-predicted natural frequencies and mode shapes of the first five modes. It is very clear from the figure that the natural frequencies were very similar and the mode shapes were matched for all five modes except for the first mode. It is believed that the mode shape discrepancy in mode 1 is mainly due to modelling error (there may be ‘some sort of’ nonlinear properties in stiffness of the ballast which will be presented as modelling error in the proposed method). The amount of information that can be extracted from the measured data is limited, and therefore, the uncertainties associated with the damage detection results will be very high if the order of polynomial is too high. The Bayesian model class selection method adopted in the proposed method kept a balance between the likelihood (its value will be reduced for an increased polynomial order) and polynomial order. The percentage differences between the measured and model-predicted natural frequencies for all five modes are summarized in Table 4. It is clear that all differences were lower than 3% indicating that the good performance of the identified model class M ₃ for given measurement data.

Table 4.

Measured and model-predicted (after model updating) natural frequencies for the damaged case (Hz).

	Mode 1	Mode 2	Mode 3	Mode 4	Mode 5
Updated	60.45	67.09	157.85	418.08	763.07
Measured	59.05	68.61	156.08	408.56	786.57
Difference (%)	2.37	−2.22	1.13	2.33	−2.99

Conclusion

A new polynomial modelling method has been developed to determine the ballast stiffness distribution under a sleeper using Lagrange polynomials to prevent artificial jumps in the ballast stiffness. The newly developed modelling method has been integrated with the MCMC-based Bayesian model updating and model class selection methods to form a practical method for detecting ballast damage.

To ensure that the method is applicable for both identifiable and unidentifiable cases, MCMC has been used to generate samples for approximating the probability of a model class conditional on a set of measurements and for approximating the posterior PDF of uncertain model parameters. This article has derived the new formulations for using MCMC samples in Bayesian model class selection.

The use of Lagrange polynomials in modelling the variation in ballast stiffness under a sleeper not only reduces the level of modelling error, but also provides physical meaning to the uncertain model parameters in the model updating process. It should be noted that it is difficult to judge the ballast damage situation from the posterior PDF of uncertain parameters if the polynomial coefficients are directly considered as the uncertain parameters in the model updating process.

There were few investigations in the literature regarding the application of MCMC-based Bayesian methods for the detection of damage in structural systems, and this article is the first to deal with the detection of damage to railway ballast. Experimental verification using measured data from the full-scale test panel has clearly proved the feasibility and practical value of the proposed method in detecting ballast damage. It provides valuable information for inspectors and engineers about the condition of the ballast under concrete sleepers during their visual inspection.

If the selected model class shows a uniform distribution of ballast stiffness under a sleeper and the identified value of the ballast stiffness is close to the nominal value, the target in situ sleeper is treated as undamaged. Otherwise, the ballast at certain locations under the sleeper may be damaged. The Bayesian model class selection is the key component in the proposed method for the selection of the most plausible model class conditional on a given set of measurements.

Using the Bayesian model updating method, the ballast stiffness distribution under the sleeper can be estimated together with the corresponding posterior uncertainties. It turns out that Bayesian model updating has high potential of application in railway ballast damage detection.

The damage simulated in this study was relatively large at the development stage of the proposed method to ensure the damage-induced change in modal parameters was not contaminated by measurement noise and modelling error. Additional experiments are needed in the future to study the effects of ballast damage extents by simulating ballast damage with various ballast sizes (e.g. 20–30 mm).

Roving hammer test data has been used in this study, and it can provide supplementary information to the inspector during visual inspection. In the final stage of this approach, the inspector can graphically see the ballast stiffness distribution and quantify the condition of the rail–sleeper–ballast system.

The frequency domain method assumes that the system behaviour is linear. However, the ballast only behaves linearly under small-amplitude vibrations. For large-amplitude vibrations, it is believed that the nonlinear effects cannot be neglected and the proposed method can only provide an approximate solution. Research is ongoing to perform Bayesian model updating in the time domain with the consideration of the nonlinear effects in the vibration. In addition, it is in general believed that time-domain data are more reliable because they are free from error, which may be introduced by conversion of the time-domain data to other domains, such as the frequency and modal domains.

It is clear from the literature that the natural frequencies of a structural system depend on the temperature and humidity.⁴⁴ In the development stage, the proposed method has been verified using measured data from the indoor test panel, which allows the environmental factors, such as temperature and humidity, to be well controlled. Since the proposed method directly estimates the ballast distribution without a set of baseline data (i.e. a set of measured modal parameters of the system in the undamaged status), the problem of temperature difference in the two measurements (one for the baseline and one for the possibly damaged system) is not valid for the proposed method. A long-term monitoring program for recording the modal parameters and the corresponding temperatures of several in situ sleepers is in progress for studying the effects of environmental factors on the modal parameters of in situ sleepers.

Footnotes

Acknowledgements

The authors would like to thank Dr Man-Tat Wong and Hong Kong MTRC for their help in the development of the full-scale test panel in Kowloon Bay, Hong Kong. The help from Dr Hua-Yi Peng and Mr Alabi Stephen in the experiment is highly appreciated.

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship and/or publication of this article.

Funding

The work described in this article was fully supported by a grant from the Research Grants Council of the Hong Kong Special Administrative Region, China (Project No. 9041889 (CityU 115413)).

References

Zeng

. Dynamic tests for heavy haul railway tracks. J China Railw Soc 1988; 10: 66–77.

Lam

Wong

Yang

. A feasibility study on railway ballast damage detection utilizing measured vibration of in situ concrete sleeper. Eng Struct 2012; 45: 284–298.

Lam

Wong

. The Bayesian methodology for the detection of railway ballast damage under a concrete sleeper. Eng Struct 2014; 81: 289–301.

Friswell

Peng

JET

Wilson

DAL

. Using vibration data and statistical measures to locate damage in structures. Modal Anal 1994; 9(4): 239–254.

Lam

Leung

AYT

. Multicrack detection on semirigidly connected beams utilizing dynamic data. J Eng Mech: ASCE 2008; 134(1): 90–99.

Yin

Lam

Chow

. A Bayesian probabilistic approach for crack characterization in plate structures. Comput-Aided Civ Inf 2010; 25: 375–386.

Salawu

Williams

. Bridge assessment using forced-vibration testing. J Struct Eng: ASCE 1995; 121(2): 161–173.

Ching

Beck

. New Bayesian model updating algorithm applied to a structural health monitoring benchmark. Struct Health Monit 2004; 3(4): 313–332.

Johnson

Lam

Katafygiotis

et al . Phase I IASC-ASCE structural health monitoring benchmark problem using simulated data. J Eng Mech: ASCE 2004; 130(1): 3–15.

10.

Fan

Qiao

. Vibration-based damage identification methods: a review and comparative study. Struct Health Monit 2011; 10(1): 83–111.

11.

Salawu

. Detection of structural damage through changes in frequency: a review. Eng Struct 1997; 19(9): 718–723.

12.

Jeffreys

. Theory of probability. 3rd ed. Oxford: Clarendon Press, 1961.

13.

Box

GEP

Tiao

. Bayesian inference in statistical analysis. Reading, MA: Addison-Wesley, 1973.

14.

Yuen

. Bayesian methods for structural dynamics and civil engineering. New York: John Wiley & Sons, 2010.

15.

Choi

Kim

. Identification of correlated damage parameters under noise and bias using Bayesian inference. Struct Health Monit 2011; 11(3): 293–303.

16.

Beck

Katafygiotis

. Updating models and their uncertainties. I: Bayesian statistical framework. J Eng Mech: ASCE 1998; 124(4): 455–461.

17.

Vanik

Beck

. Bayesian probabilistic approach to structural health monitoring. J Eng Mech: ASCE 2000; 126(7): 738–745.

18.

Zhang

et al . Full-scale dynamic testing and modal identification of a coupled floor slab system. Eng Struct 2012; 37: 167–178.

19.

Lam

Peng

. Development of a practical algorithm for Bayesian model updating of a coupled slab system utilizing field test data. Eng Struct 2014; 79: 182–194.

20.

Friswell

Mottershead

. Finite element model updating in structural dynamics. Dordrecht: Springer Science and Business Media, 1995.

21.

Mares

Friswell

Mottershead

. Model updating using robust estimation. Mech Syst Signal Pr 2002; 16(1): 169–183.

22.

Katafygiotis

Lam

. Tangential-projection algorithm for manifold representation in unidentifiable model updating problems. Earthq Eng Struct D 2002; 31(4): 791–812.

23.

Lam

Yang

. Bayesian model updating of a coupled-slab system using field test data utilizing an enhanced Markov chain Monte Carlo simulation algorithm. Eng Struct 2015; 102(1): 144–155.

24.

Geyer

Thompson

. Annealing Markov chain Monte Carlo with applications to ancestral inference. J Am Stat Assoc 1995; 90: 909–920.

25.

Tierney

. Markov chains for exploring posterior distributions. Ann Stat 1994; 22: 1701–1762.

26.

Beck

. Bayesian updating of structural models and reliability using Markov chain Monte Carlo simulation. J Eng Mech: ASCE 2002; 128(4): 380–391.

27.

Ching

Chen

. Transitional Markov chain Monte Carlo method for Bayesian model updating, model class selection, and model averaging. J Eng Mech: ASCE 2007; 133(7): 816–832.

28.

Burden

Faires

. Numerical analysis. Boston, MA: Brooks/Cole, 2001.

29.

Beck

Yuen

. Model selection using response measurements: Bayesian probabilistic approach. J Eng Mech: ASCE 2004; 130(2): 192–203.

30.

Lam

. Identification of the degradation of railway ballast under a concrete sleeper. In: Proceedings of the 2013 world congress on advances in structural engineering and mechanics (ASEM-13), Jeju, Korea, 8–12 September 2013.

31.

Krenk

. Mechanics and analysis of beams, columns and cables – a modern introduction to the classic theories. 2nd ed. Berlin; Heidelberg: Springer Science and Business Media, 2001.

32.

Przemieniecki

. Theory of matrix structural analysis. Mineola, NY: Courier Corporation, 1985.

33.

Papadimitriou

Ntotsios

Giagopoulos

et al . Variability of updated finite element models and their predictions consistent with vibration measurements. Struct Control Hlth 2012; 19(5): 630–654.

34.

Katafygiotis

Lam

Papadimitriou

. Treatment of unidentifiability in structural model updating. Adv Struct Eng 2000; 3(1): 19–39.

35.

Hastings

. Monte Carlo sampling methods using Markov chains and their applications. Biometrika 1970; 57(1): 97–109.

36.

Metropolis

Rosenbluth

et al . Equation of state calculations by fast computing machines. J Chem Phys 1953; 21(6): 1087–1091.

37.

Beck

. A new adaptive importance sampling scheme for reliability calculations. Struct Saf 1999; 21: 135–158.

38.

Robert

Casella

. Monte Carlo statistical methods. New York: Springer Science and Business Media, 2013.

39.

Cheung

Beck

. Calculation of posterior probabilities for Bayesian model class assessment and averaging from posterior samples based on dynamic system data. Comput-Aided Civ Inf 2010; 25(5): 304–321.

40.

Muto

Beck

. Bayesian updating and model class selection for hysteretic structural models using stochastic simulation. J Vib Control 2008; 14(1–2): 7–34.

41.

Beck

. Determining models of structures from earthquake records. PhD Thesis, California Institute of Technology, Pasadena, CA, 1978.

42.

MTR. HK MTR materials and workmanship specification for track work. Hong Kong: MTR Corporation Limited, 2009.

43.

Kumaran

Menon

Nair

. Dynamic studies of railtrack sleepers in a track structure system. J Sound Vib 2003; 268: 485–501.

44.

Xia

Hao

Zanardo

et al . Long term vibration monitoring of an RC slab: temperature and humidity effect. Eng Struct 2006; 28(3): 441–452.