Random dynamic analysis of vertical train–bridge systems under small probability by surrogate model and subset simulation with splitting

Abstract

The response of the train–bridge system has an obvious random behavior. A high traffic density and a long maintenance period of a track will result in a substantial increase in the number of trains running on a bridge, and there is small likelihood that the maximum responses of the train and bridge happen in the total maintenance period of the track. Firstly, the coupling model of train–bridge systems is reviewed. Then, an ensemble method is presented, which can estimate the small probabilities of a dynamic system with stochastic excitations. The main idea of the ensemble method is to use the NARX (nonlinear autoregressive with exogenous input) model to replace the physical model and apply subset simulation with splitting to obtain the extreme distribution. Finally, the efficiency of the suggested method is compared with the direct Monte Carlo simulation method, and the probability exceedance of train responses under the vertical track irregularity is discussed. The results show that when the small probability of train responses under vertical track irregularity is estimated, the ensemble method can reduce both the calculation time of a single sample and the required number of samples.

1 Introduction

With the development of high-speed railways, many high-speed railway bridges have been used to replace the subgrade. Random vibration, excited by the stochastic track irregularity, is generated when a train runs on the bridge. As there is a high traffic density of trains, for example, on Beijing–Shanghai high-speed railway, where the volume of trains can reach 100 pairs per day, the railway department formulates a specific rule for monthly rail track maintenance, and preventive rail grinding is conducted at least once a year to ensure the smoothness of the rail track [1]. If we assume that the recurrence period of the maximum train responses is 1 year and 100 pairs of trains pass a bridge every day, the total number of train pairs that move on the bridge in a year is 36,500, and the exceedance probability of the maximum train response is 36,500⁻¹=2.74 × 10⁻⁵. Therefore, it is necessary to ensure the random response of the train–bridge system under small probability conditions.

For the train–bridge system, the randomness of track irregularity is an external excitation. The response of the train–bridge system is also random; thus, the dynamic response of the train–bridge system should be studied from the perspective of statistical characteristics. Based on the stochastic vibration theory and coupling dynamics of the rail track and wheel-set, the coupling vibration method [2,3,4] has been widely applied to the security analysis of high-speed trains. To analyze the reliability of moving trains on a bridge under stochastic excitations, random vibration methods, such as the pseudo-excitation method [5, 6], spectral approach [7, 8], and probability density evolution method [9] have been applied in the stochastic analysis of train–bridge systems. Additionally, the Monte Carlo simulation (MCS) method has been applied to reliability analysis of the train–bridge systems [10], and the stochastic characteristics of coupling systems have been explored in a way of the multi-sample calculation [11]. The random vibration method usually obtains statistical properties, such as the variance, power spectral density and mean value; and the reliability of the train–bridge system is evaluated by a statistical index. However, as the maximum train response may be a non-Gaussian distribution, there is deviation in the evaluation of dynamic performance under small probabilities. Though the MCS method is one of the commonly used approaches for scientific computing due to its simplicity and general applicability [12], it requires many samples to achieve high precision.

The subset simulation (SS) method [13], also a controlled MCS method, can reduce the sample number to estimate small probabilities. The principle of the SS approach is to decompose the small probability into a product of larger probabilities under the intermediate event condition. Ching et al. [14] used the concept of response trajectory splitting and proposed the method called subset simulation with splitting (SS/S). Ding and Chen [15] combined the SS/S method with the multivariate autoregressive (AR) modeling of stochastic excitations, and promoted the application of the SS/S method. The SS/S method with an improvement of an order of magnitude has higher computational efficiency than the MCS method at the same level accuracy, but it requires more samples. When the SS/S method is applied to estimate the extreme distribution of train–bridge systems, if the exceedance probability is 1.0 × 10⁻⁴, the number of samples is approximately 10⁴, which means that the calculations take a lot of time [16]. Therefore, the calculation time of one sample should be reduced.

A surrogate model can be used to approximate a complex system and reduce the calculation time of one sample. Surrogate models such as Kriging, radial basis function, and support vector regression; and artificial neural network (ANN) have been applied in optimization, reliability analysis and wind engineering [17,18,19]. Although the surrogate models are used to predict time series [20], the current response of dynamic systems relates to historical responses and external excitation. Therefore, the nonlinear autoregressive with exogenous input (NARX) model is introduced to the surrogate model for a dynamic system. NARX model can be used to build a system model, which only relies on the input and output data, and does not need other system information [21]. Numerical examples show that the NARX surrogate model has good prediction performance, and it is widely used to predict the dynamic response of systems [22,23,24,25,26,27]. However, the training and calculation efficiency of the NARX surrogate model may be impaired by an excessive number of input sequences, such as a multivariate wind field, or lateral acceleration of the train body. Therefore, the ensemble method, which combines the SS method with the surrogate model, can significantly improve the calculation efficiency for small failure probability.

This paper is organized as follows: Section 2 reviews the coupling model of the train–bridge system. Section 3 presents the framework of ensemble methods based on NARX surrogate model and SS/S, and a two-degree car model is used to verify the ensemble method. In Sect. 4, the vertical responses of the train are discussed under small probabilities, and the accuracy and efficiency of the proposed method are proven by a direct MCS method.

2 Train–bridge systems

The numerical calculation model of the train–bridge system consists of two parts, the train subsystem and the bridge subsystem. In this section, the train model and bridge model are illustrated, respectively.

2.1 Train model

For the train subsystem, the train model is considered as a four-axle train of seven rigid bodies: one train body, two bogies and four wheel-sets. The rigid bodies are connected by springs and dampers [28]. The train model is shown in Fig. 1.

Fig. 1

Mass-spring-damping model of the train [28]

The train body and each bogie have 5 degrees of freedom, and each wheel-set has 2 degrees of freedom. Therefore, the train model has 23 degrees of freedom in total. The motion equations of each part of the train can be established by D’lembert principle [28]:

$$ \varvec{M}_{\text{v}} \ddot{\varvec{u}}_{\text{v}} + \varvec{C}_{\text{v}} \dot{\varvec{u}}_{\text{v}} + \varvec{K}_{\text{v}} \varvec{u}_{\text{v}} = \varvec{F}_{\text{v}} , $$

(1)

where M_v, C_v and K_v are the mass, damping and stiffness matrices of the trains, respectively; F_v is the vector of the wheel-rail contact force and, u_v is the vector of train displacement.

The wheel-set considers only two degrees of freedom: lateral and yaw directions; the effect of the vertical track irregularity directly acts on the bogie, which can be written as follows:

$$ {\text{Floating vibration loads}}:\;F_{{\text{t}}zj} = K_{{\text{p}}z} \sum\limits_{i = m}^{m + 1} {z_{{\text{eq}}i} } + C_{{\text{p}}z} \sum\limits_{i = m}^{m + 1} {\dot{z}_{{\text{eq}}i} } , $$

(2)

$$ {\text{Pitching vibration loads}}:\;F_{{\text{t}}\upvarphi {\text{j}}} = - l_{1} K_{{\text{p}}z} \sum\limits_{i = m}^{m + 1} {( - 1)^{i} z_{{\text{eq}}i} } - l_{1} C_{{\text{p}}z} \sum\limits_{i = m}^{m + 1} {( - 1)^{i} \dot{z}_{{\text{eq}}i}},$$

(3)

where j = 1, 2 denote the front and rear bogies and m = (j − 1) × 2 + 1; i = 1, 2, 3, 4 denote four wheel-sets; l₁ is half of the distance between two wheel-set centroids; K_pz and C_pz are the vertical stiffness and damping of the primary suspension, respectively; z_eqi and $\dot{z}_{{\text{eq}}i}$ are the equivalent irregularity of the displacement and velocity terms, respectively, which is the superposition of bridge displacement and track irregularity.

Although the train model is a 3D model, given that predicting the lateral response of the train system with a parametrical nonlinearity is difficult using the surrogate model, only the vertical track irregularity is considered and the vertical response of the train–bridge system is explored in this work.

2.2 Bridge model

The finite element method is used to build the bridge model, and the motion equation of the bridge structure is expressed as

$$ \varvec{M}_{\text{b}} \ddot{\varvec{u}}_{\text{b}} + \varvec{C}_{\text{b}} \dot{\varvec{u}}_{\text{b}} + \varvec{K}_{\text{b}} \varvec{u}_{\text{b}} = \varvec{F}_{\text{b}} , $$

(4)

where M_b, C_b and K_b are the mass, damping and stiffness matrices of the bridge, respectively; F_b is the wheel-rail contact force, which is equal and opposite to F_v; u_b is the bridge displacement.

Commercial finite element software is commonly used to establish a finite element model of the bridge; then, the mass and stiffness matrices of the model can be derived from the commercial code. Using the modal analysis, the natural vibration characteristics of the structure will be obtained, and the damping matrix of the bridge model is obtained by the use of Rayleigh damping. This method will have some errors in the analysis of long-span suspension bridges in view of the geometry and material nonlinearity. Here, the mass, damp and stiffness matrices of the bridge are imported in MATLAB, and the separate iteration method is used to solve the train and bridge subsystems [28, 29]. The track irregularity is generated by an auto-regression (AR) model.

3 Ensemble methods based on the surrogate model and SS/S

3.1 NARX-based surrogate model of dynamic system

The NARX model, which has wide application in many fields, is expressed as follows [30, 31]:

$$ \begin{aligned} y\left( t \right) & = f(y\left( {t - 1} \right),y\left( {t - 2} \right), \ldots ,y\left( {t - n_{\text a} } \right),u\left( {t - n_{\text k} } \right), \\ & \quad u\left( {t - n_{\text k} - 1} \right), \ldots ,u\left( {t - n_{\text k} - n_{\text b} - 1} \right)) + e\left( t \right), \\ \end{aligned} $$

(5)

where $y\left( \cdot \right)$ and $u\left( \cdot \right)$ are the outputs and inputs, respectively, with different delays; $e\left( t \right)$ is the noise terms; $f\left( \cdot \right)$ denotes nonlinear functions, including linear, neural network, tree partition, wave network and sigmoid function; n_a is the number of outputs in the past; n_b is the number of past inputs, including the current time input; n_k is the starting point of input delay. To simplify the model, n_k is set as zero.

An autocorrelation analysis of the output is applied to determine n_a, and cross-correlation analysis of the input and output is applied to determine n_b. If the correlation function is not close to 0, it is believed that the past input or output of a certain delay has affected the output of the current time.

The vibration displacement and velocity of the bridge relate to the train loads. As the train–bridge system is a time-variable system, the bridge response to the moving train has some difference under various track irregularities. If the NARX model is used to predict the response of trains, the inputs should include the track irregularity that is generated by the power spectral and vibration response of the bridge. This is difficult for the NARX model because the vibration response of the bridge is unknown. In this work, the vertical vibration response of the bridge is replaced with the bridge response with all of the track irregularity set to zero. The displacement and velocity of the bridge calculated by this simplified method are similar to those calculated with consideration of the vertical track irregularity [26]. A train has four wheel-sets, so the NARX model has four external excitation inputs, and the expressions are shown in Eqs. (2) and (3). The output of the NARX model is the acceleration of the train body.

For the vertical train–bridge system, many functions can be used to predict the response of train. Function f (see Eq. (5)) is more complex, and the accuracy of the surrogate model may be higher, but complex function f will decrease the training and calculation efficiency of the surrogate model. Linear function f has a higher training efficiency than the NARX neural network model because the selection of random initial weights is avoided. Thus, it is used to fit and predict the vertical response of the train.

3.2 Subset simulation with splitting

3.2.1 Reviews of SS/S

SS/S [13] is a kind of MCS method based on conditional probability. The probability statistics in each generation sequence can be counted according to the method of MCS and is multiplied by the probability of the parent generation, to obtain the target probability level of the current generation. This method can be expressed as

$$ \Pr (Y > y_{j} ) = \mathop \prod \limits_{i = 1}^{j} p_{i} ,\quad j = 1, 2, \ldots ,m, $$

(6)

where y_j is the threshold of the jth failure region; p_i is the probability of reaching a threshold in the ith failure region.

The method to generate offspring trajectories using the SS/S method is shown in Fig. 2, where X_m is the parent response trajectory, X_c is the newly generated offspring response trajectory, y_i is the threshold of the ith failure region, and S is the splitting point, at which X_m first reaches the threshold. When the parent response trajectory first reaches the threshold value of the response failure region, the corresponding time point is used as the splitting point to regenerate the excitation trajectory; then, the response trajectory under the new excitation is generated. Combined with the response trajectory before the splitting point, an offspring response trajectory is generated.

Fig. 2

Generation of an offspring response trajectory by SS/S

3.2.2 Parameter selection

When the SS/S method is used, some factors must be determined. First, the target failure probability P_F, intermediate failure probability $\bar{p}$, and number of intermediate failure regions m are determined.

Assuming that the probability of failure and the number of samples in each intermediate failure region are equal, i.e., $p_{1} = \cdots = p_{m} = \bar{p}$, $N_{1} = \cdots = N_{m} = \bar{N}$. Then, the total number of samples is $\bar{N} + \bar{N}\left( {m - 1} \right)\left( {1 - \bar{p}} \right) = N_{T}$. The coefficient of variation (c.o.v.) of estimated value $\hat{P}_{\text F}$ is not sensitive to $\bar{p}$ when the intermediate probability is $\bar{p} \ge 0.1$ [14]. Therefore, $\bar{p} = 0.1$ is taken in this study, and the number of intermediate failure events m can be calculated according to the following formula:

$$ m = \frac{{\ln p_{\text{F}} }}{{\ln \bar{p}}}. $$

(7)

Next, the total number of samples is determined.

For the SS/S method, $\left\{ {\left[ {S^{{\left( {i,k} \right)}} } \right]:k = 1, 2 ,\ldots R_{i - 1} } \right\}$ is assumed to be the point that represents the first time that the kth offspring response trajectory reaches the threshold. $P_{i} \left( {S^{{\left( {i,k} \right)}} } \right)$ is expressed as

$$ P_{i} \left( {S^{{\left( {i,k} \right)}} } \right) \equiv P\left( {\left\{ {{\text{a }}\;{\text{trajectory}}\;{\text{hits}}\;g\left( X \right)} \right\}\left| {\left\{ {{\text{trajectories \;start\; from}} \;S^{{\left( {i,k} \right)}} } \right\}} \right.} \right). $$

(8)

According to Ref. [13], the expected value and variance of $P_{i} \left( {S^{{\left( {i,k} \right)}} } \right)$ satisfy the following equations:

$$ E_{Si} \left[ {P_{i} \left( {S^{{\left( {i,k} \right)}} } \right)} \right] = p_{i} ,\quad k = 1, 2, \ldots ,R_{i - 1} , $$

(9)

$$ 0 \le {\text{Var}}_{Si} \left[ {P_{i} \left( {S^{{\left( {i,k} \right)}} } \right)} \right] \le p_{i} \left( {1 - p_{i} } \right),\quad k = 1, 2, \ldots ,R_{i - 1} . $$

(10)

When $N_{i} \gg 1$, the variance of the estimated value $\hat{P}_{i}$ can be obtained by the following formula:

$$ {\text{Var}}\left( {\hat{P}_{i} } \right) = \frac{{p_{i} \left( {1 - p_{i} } \right)}}{{N_{i} }} + \frac{{{\text{Var}}_{Si} \left[ {P_{i} \left( {S^{{\left( {i,1} \right)}} } \right)} \right]}}{{R_{i - 1} }}. $$

(11)

The c.o.v. of the estimated value $\hat{P}_{i}$ is

$$ \delta_{i}^{2} = \frac{{{\text{Var}}\left( {\hat{P}_{i} } \right)}}{{\left\{ {E_{Si} \left[ {P_{i} \left( {S^{{\left( {i,k} \right)}} } \right)} \right]} \right\}^{2} }} = \frac{{1 - p_{i} }}{{N_{i} p_{i} }}\left( {1 + \gamma_{i} } \right), $$

(12)

where $\gamma_{i} = E\left( {\frac{{N_{1} }}{{R_{i - 1} }}} \right)\frac{{{\text{Var}}_{Si} \left[ {p_{i} \left( {S^{{\left( {i,1} \right)}} } \right)} \right]}}{{p_{i} \left[ {1 - p_{i} } \right]}}$ when $i = 1$, $\gamma_{i} = 0$.

After the c.o.v. of estimated value $\hat{P}_{i}$ is obtained, the c.o.v. of estimated value $\hat{P}_{\text{F}}$ can be obtained [13]:

$$ \delta^{2} = \mathop \sum \limits_{i = 1}^{m} \delta_{i}^{2} \approx \frac{{1 - p_{1} }}{{N_{1} p_{1} }} + \mathop \sum \limits_{i = 2}^{m} \frac{{1 - p_{i} }}{{N_{i} p_{i} }}\left( {1 + \gamma_{i} } \right). $$

(13)

Equation (13) can be simplified as

$$ \delta \approx \sqrt {\frac{{\left( {1 - p_{1} } \right)\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{\gamma } }}{{N_{1} p_{1} }}} , $$

(14)

where $\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{\gamma } = \mathop \sum \nolimits_{i = 1}^{m} \left( {1 + \gamma_{i} } \right)$ when $i = 1$, and $\gamma_{i} = 0$.

If $\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{\gamma } \ge 0$, the lower bounds of sample number of the SS/S method can be written as

$$ N_{1} \ge \frac{{m\left( {1 - p_{1} } \right)}}{{\delta^{2} p_{1} }}. $$

(15)

Equation (15) can be used to estimate the minimum number of samples. Generally, a small δ is selected to obtain a good precision, such as δ = 0.1.

3.2.3 Generation of stochastic excitation

In the framework of SS/S, Ding and Chen [15] proposed the AR model to simulate the stochastic process. The track irregularity is a univariate stochastic process, and the p-order AR model of track irregularity can be expressed as

$$ W\left( t \right) = \mathop \sum \limits_{k = 1}^{p} A\left( k \right)W\left( {t - k\Delta t} \right) + Ln\left( t \right), $$

(16)

where p is the order of the AR model, which is generally taken as 4; A(k) is the autoregressive coefficient; W(t) is the value of track irregularity at the current time; Δt is the length of the time step; n is Gaussian white noise with the expected value of 0 and variance of 1; L is a coefficient.

Equation (16) is multiplied by W(t − jΔt) on the left and right sides. Taking the mathematical expectations, we obtain the following results:

$$ R_{\text{w}} \left( {j{{\Delta }}t} \right) = \mathop \sum \limits_{k = 1}^{p} R_{\text{w}} \left[ {\left( {j - k} \right){{\Delta }}t} \right]A\left( k \right), $$

(17)

$$ R_{0} = R_{\text{w}} \left( 0 \right) - \mathop \sum \limits_{k = 1}^{p} A\left( k \right)R_{\text{w}} \left( {k{{\Delta }}t} \right), $$

(18)

where R_w(jΔt) is the autocorrelation function of track irregularity.

The power spectral density function and autocorrelation function of the stochastic process satisfy Wiener–Khinchin formula:

$$ R_{\text{w}} \left( {j{{\Delta }}t} \right) = \mathop \int \limits_{0}^{\infty } S_{\text{w}} \left( f \right)\cos (2{\uppi}{f} \cdot j{{\Delta }}t){\text{d}}f, $$

(19)

where S_w(f) is the track irregularity power spectrum function. Therefore, the German low interference spectrum is used to generate the track irregularity.

The autocorrelation function of the track irregularity time history can be calculated according to Eq. (19), when j = 0, 1, 2,…, p. The autoregressive coefficient A(k) is obtained by substituting the autocorrelation function into Eq. (17). Then R₀ is obtained by substituting A(k) and the autocorrelation function into Eq. (18), and L is determined as follows:

$$ L = \sqrt {R_{0} } . $$

(20)

For the track irregularity time history, the length of time step Δt is obtained by substituting A(k) and L into Eq. (16).

Taking the German high-speed track spectrum as an example, the generation of offspring samples in the SS/S method is simulated, and the power spectrum and time history of offspring samples are shown in Fig. 3. As Fig. 3 shows, the simulated spectrum is consistent with the target spectrum, and the simulated sample is validated. The offspring sample has an identical trajectory to the parent sample before the splitting point, and the trajectory only changes after the splitting point, which satisfies the requirements of the SS/S method.

Fig. 3

Simulation of track irregularity: a power spectrum and b track irregularity of offspring by splitting

3.3 Analytical procedure

The time history of the train response is split at threshold point S (Fig. 2) in the SS/S method, and the offspring track irregularity is generated by the AR model. The NARX model can be considered as a surrogate model of the time series, and it is applied to approximate the time history of train responses. Therefore, it is easy to create a link between the SS/S method and the NARX surrogate model if the surrogate model has sufficient accuracy.

The ensemble method, which is composed of the NARX surrogate model and SS/S method, is called the NARX-SS/S method. In the ensemble method, the vertical response of the train in Eq. (1) is replaced by the NARX model, which is expected to reduce the calculation time of each sample, and the SS/S method is used to reduce the total number of samples.

There are two important objectives in the NARX-SS/S method: obtain the NARX model of the train bridge system and use the NARX model in the framework of the SS/S method. Therefore, the analysis process of NARX-SS/S is summarized as follows:

• Step 1 The time history of the bridge response at each wheel-set position is obtained from Eqs. (2)–(4) when all track irregularities are set as zero, which is an approximate input to replace the dynamic response of the bridge.

• Step 2 Equations (1) and (4) are solved under a set of arbitrary random track irregularities with some initial steps set as zero to ensure identical initial condition for all cases.

• Step 3 A correlation analysis is conducted to determine the initial time delays n_a and n_b, and the accurate time delays are obtained by trial calculation. The sample is trained by the NARX model, and the surrogate model of the vertical response of the train is obtained.

• Step 4 The track irregularity is generated by the AR model, N₁ samples of first level are performed by the NARX surrogate model, and the maximum sample in one level is calculated by Eqs. (1) and (4). If the maximum absolute error between the target value and the prediction value exceeds 5%, the input and output for the maximum value are used to retrain the NARX surrogate model, and N₁ samples of the first level are recalculated by a new NARX surrogate model [15].

• Step 5 The parent samples of the first level, which correspond to the intermediate failure probability, are searched, the splitting point is determined, and the conditional track irregularity is generated by the AR model. The offspring samples are generated by the NARX surrogate model. Similarly, the maximum absolute error is determined by Step 4.

• Step 6 Step 5 is repeated until the target failure probability is reached, and the SS/S method is used to calculate the probability of exceedance of the train response based on Eq. (6).

3.4 Numerical verification

To quickly verify the ensemble NARX-SS/S method, a two-degree-of-freedom model is adopted. The model diagram is shown in Fig. 4.

Fig. 4

Quarter train model

The equation of the quarter train model can be written as follows:

$$ \left[ {\begin{array}{*{20}c} {m_{1} } & 0 \\ 0 & {m_{2} } \\ \end{array} } \right] \cdot \left\{ {\begin{array}{*{20}c} {\ddot{x}_{1} } \\ {\ddot{x}_{2} } \\ \end{array} } \right\} + \left[ {\begin{array}{*{20}c} {c_{1} } & { - c_{1} } \\ { - c_{1} } & {c_{1} + c_{2} } \\ \end{array} } \right] \cdot \left\{ {\begin{array}{*{20}c} {\dot{x}_{1} } \\ {\dot{x}_{2} } \\ \end{array} } \right\} + \left[ {\begin{array}{*{20}c} {k_{1} } & { - k_{1} } \\ { - k_{1} } & {k_{1} + k_{2} } \\ \end{array} } \right] \cdot \left\{ {\begin{array}{*{20}c} {x_{1} } \\ {x_{2} } \\ \end{array} } \right\} = \left[ {\begin{array}{*{20}c} 0 \\ {k_{2} z + c_{2} \dot{z}} \\ \end{array} } \right], $$

(21)

where m₁ = 24 t; m₂ = 3.2 t; k₁ = 800 kN/m; k₂ = 2080 kN/m; c₁ = 66 kN/(m/s); c₂ = 60 kN/(m/s); $\dot{z} = z^{\prime}V$; V is the train speed; z is the track irregularity with a step of 0.2 m, which is generated by the AR model, and $z^{\prime}$ is the derivative of $z$ with respect to displacement.

Equation (21) is solved by the Newmark-β method, with a train speed of 200 km/h, a step length of 8000, a time step of 0.0036 s and a space length of 1600 m. The German high-speed track spectrum is used to reconstruct the track irregularity.

The target probability of exceedance P_F is 10⁻⁴, and the basic setups of the splitting method are m = 4 and p_i= 0.1 (i = 1,2, …, m). If the c.o.v. of δ is set to 0.1, 3600 samples are required for each intermediate failure region by the lower bounds of the SS/S method [13], i.e. N₁≥ 3600. In this work, 4000 samples are used for each intermediate failure region. Thus, 4000 + 3600 + 3600 + 3600 = 14,800 samples are required in total.

When we look for the splitting point of the parent sample, the splitting point may not be in the time step (Fig. 5), and there are two methods to address this issue. One method is to first obtain the values of z(t − Δt + Δs) time steps from the splitting point by interpolation, take them as the initial values, and use the AR model to generate the samples; then, the values at the original spacing point are obtained by interpolation. Another method is to directly take the value of one step after the split point (as shown in Fig. 5, taking z(t) as the splitting point). The second method can reduce the rejected ratio in the generation of offspring samples; so it has been used to determine the splitting point in the course of verifying the ensemble model. Additionally, the train with four wheel-sets has multiple inputs, and the track irregularity of the first wheel-set is used to determine the splitting point.

Fig. 5

Split point of the responses and excitations

Figure 6 shows the exceedance probability of the absolute maximum accelerations of m₁, which is obtained by the direct MCS method and NARX-SS/S method. From Fig. 6a, a theoretical curve of probability of exceedance based on Eq. (21) is obtained by 1 million MCS samples, and the estimated c.o.v. of the MCS is 0.1. The two methods keep consistent with the target probabilities of exceedance of 10⁻⁴. Although the c.o.v. of the NARX-SS/S method is 0.188, which is 1.88 times higher than that of the MCS method when the probability of exceedance is 10⁻⁴, the c.o.v. of the two methods is relatively small, and the final result still fits well. Figure 6b shows that under different failure probabilities, when the MCS method uses as many samples as the NARX-SS/S method, the c.o.v. values of the two methods show that the NARX-SS/S method increase in efficiency with the decrease in exceedance probability. In other words, the number required is significantly less than that of the direct MCS method.

Fig. 6

Comparison of MCS and NARX-SS/S methods: a exceedance probability of train body and b comparison of c.o.v. for different estimators

4 Case study

Since the track irregularity sample is simulated by the track power spectrum, the sample length will affect the coincidence between the simulated and targeted spectra. A longer sample size corresponds to higher accuracy. To ensure that the track irregularity sample and the responses of the train are the stationary process, the bridge model consists of 50 pre-stressed concrete simple-supported beams, each of which has a span of 32 m, and the total length of the bridge is 1600 m. The main beam section is 12 m wide and 3.05 m high. The heights of all piers are set to 10 m, the first 3 spans of the main beam are rigid, and the track irregularity is set as zero. Suppose that these main beams do not deform and initial conditions are identical for every sample. The bridge model is simulated by a spatial beam element, the damping ratio of the bridge structure is 2%, and the cross-section of the bridge is shown in Fig. 7. The train model of CRH3 is introduced in Sect. 2.1. The track irregularity is simulated by the German high-speed track spectrum, and the space step is 0.2 m.

Fig. 7

Bridge cross-section (mm)

Evaluation indices such as the acceleration of the train body and wheel load reduction rate can be used to evaluate the safety and comfort of the train. The threshold values of the vertical acceleration and wheel load reduction rate in the specifications are 1.0 m/s² and 0.6, respectively. The thresholds of evaluation indices do not have a definite value of exceedance probability, but the target probability of exceedance has an important effect on the dynamic performance evaluation of the train. The vibration of the train is mainly caused by track irregularities without considering other factors, such as train failures and bridge damages. Generally, when the maintenance period of the track is longer, the recurrence period is smaller and the computing time of reliability analysis is more. Therefore, the target exceedance probability of the train evaluation indices can be determined by the maintenance period of the track.

It is assumed that 55 pairs of trains pass a bridge every day and the maintenance period of the track is half a year. The number of train pairs that move on the bridge in half a year is 55 × 180 = 9900, and the exceedance probability is 9900⁻¹≈ 1 × 10⁻⁴. The target probability of exceedance P_F for the vertical acceleration of the train body is set as 10⁻⁴; the basic setting of the SS/S method is m = 4 and p_i= 0.1 (i = 1,2,…, m). If δ = 0.1, then, N₁ ≥ 3600. In this work, the sample number of each intermediate failure area is set as 4000. Thus, 4000 + 3600 + 3600 + 3600 = 14,800 samples are required in total.

Since directly solving the train–bridge coupling equation has very low efficiency, the MCS method in this study calculates only 30,000 samples using Eqs. (1) and (4). Figure 8 presents the exceedance probability of the absolute maximum accelerations of the first vehicle body, which is obtained by the MCS method and NARX-SS/S method. In Fig. 8a, when the MCS method uses only 30,000 samples, if the target probability of exceedance is less than 10⁻³, the MCS method and NARX-SS/S method have consistent results. This is predictable because when the MCS method uses only 30,000 samples, and the target probability of exceedance is 10⁻³, its estimated c.o.v. is 0.183, and the c.o.v. of NARX-SS/S method is 0.181; i.e., the MCS method has almost an identical c.o.v. to the NARX-SS/S method. When the target probability of exceedance is 10⁻⁴, the c.o.v. values of the estimated value of the MCS and NARX-SS/S methods are 0.577 and 0.217, respectively; i.e., the MCS method has a far greater c.o.v. than the NARX-SS/S method. Therefore, when the target probability of exceedance is within 10⁻³–10⁻⁴, the two methods have different results. The difference is due to the insufficient number of samples in the MCS method, and the predicted value of the surrogate model, which is selected from the maximum sample of the SS/S method, is consistent with the target value obtained from Eqs. (1) and (4) by directly solving the train–bridge coupling equation (Fig. 9). When there are sufficient samples (e.g., 10⁶ samples) in the MCS method, the results of the two methods should be identical to Fig. 6a, and show consistency. The vertical acceleration of the vehicle body is less than the threshold with an acceleration of 1.0 m/s in the target probability of exceedance P_F = 10⁻⁴ because the train running speed is lower than the design vehicle speed, the simply supported beam bridge has a large stiffness, and a good track irregularity spectrum is used in this work.

Fig. 8

Comparison of MCS and NARX-SS/S methods: a exceedance probability of the train body acceleration and b comparison of the c.o.v. for different estimators

Fig. 9

Comparison of the predicted value of the surrogate model and the target value: a time history of acceleration and b errors

Figure 8b shows that the c.o.v. of the MCS method takes as many samples as the NARX-SS/S method under different failure probabilities. Figure 8b shows the same conclusion as Fig. 6b. When the target probability of exceedance is sufficiently large (e.g. 10⁻²), the efficiency difference between the two methods is not significant. With the decrease in probability of exceedance, the NARX-SS/S method has higher efficiency; i.e., it requires fewer samples than the MCS method.

With the MCS method, the calculation time of Eqs. (1) and (4) is around 320 s, while it only takes approximately 0.9 s to calculate a sample of the surrogate model, and the train time of a surrogate model is less than 5 min (a PC with 2.6 GHz Intel Xeon CPU and 64-bit OS was used for all numerical calculations), which is significantly faster than directly solving Eqs. (1) and (4). The number of samples in the SS/S method is significantly less than that in the MCS method. Therefore, the NARX-SS/S method can be used to analyze the probability distribution, and the computational efficiency can be greatly improved. The ensemble methods can be used to acquire the small probability characteristic of the dynamic system under the stochastic excitation. The accuracy of the ensemble method is mainly determined by the precision of the surrogate model.

5 Conclusions

In this work, an ensemble method based on the NARX surrogate model and SS/S method is presented, and the correctness and efficiency of the proposed method are verified by the MCS method. The ensemble method is applied to acquire the small probability distribution of the train response under vertical track irregularity. The results show that the ensemble method NARX-SS/S has the advantages of the SS/S method and surrogate model and its computational efficiency is obviously improved. In other words, the calculation time of a single sample is decreased, and the required number of samples is greatly reduced when the small probability of train responses under vertical track irregularity is estimated.