Ensemble learning-based seismic vulnerability evaluation of near-fault horseshoe-shaped tunnels

Abstract

Near-fault pulse-like ground motions, characterized by high-energy velocity pulses, pose a significant threat to the seismic safety of underground structures. However, traditional seismic fragility analysis often relies on incremental dynamic analysis using scaled ground motions, which may distort the inherent distance-dependent spectral features of near-fault records. This study proposes a data-driven seismic fragility assessment framework for horseshoe-shaped tunnels utilizing Multiple Stripes Analysis (MSA) and ensemble learning architectures. A high-fidelity numerical model of the horseshoe-shaped tunnel is established, and MSA is adopted with unscaled ground motion records to preserve the inherent physical integrity of pulse-like signals. Based on this model, a large-scale dataset is generated through nonlinear dynamic simulations. Three distinct ensemble learners, Random Forest, Extreme Gradient Boosting, and Categorical Boosting (CatBoost), are developed as surrogate models to predict damage measures (DM) based on ground motion intensity measures (IMs). Results demonstrate that the CatBoost model achieves superior predictive performance, effectively capturing the complex nonlinear mapping between ground IMs and DM. Shapley Additive Explanations analysis quantifies the contributions of various IMs, revealing that velocity-based and energy-based indicators exhibit higher importance under pulse-like motions. Finally, seismic fragility curves are derived, showing good consistency with numerical results. The proposed framework provides a reliable and computationally efficient tool for the rapid seismic risk screening of tunnels in near-fault regions.

Graphical abstract

Keywords

near-fault pulse-like ground motion / ensemble learning / horseshoe-shaped tunnel / Shapley Additive Explanation / seismic fragility analysis

1 Introduction

With the development of transportation infrastructure toward greater depth, larger scale, and increasingly complex geological conditions, tunnel engineering is facing growing safety risks under strong seismic loading [1–3]. Numerous post-earthquake damage investigations indicate that, compared with far-field ground motions, near-fault ground motions tend to cause more severe damage to underground structures due to their pronounced velocity pulse effects, permanent ground displacements, and strong directional characteristics [4–6]. In particular, tunnels constructed in the vicinity of faults or crossing fault zones are highly susceptible to lining cracking, differential displacement, and even local instability during strong earthquakes, which seriously threatens the post-earthquake functional recovery of transportation systems [7–9].

Existing studies have made significant progress in revealing tunnel seismic damage mechanisms and evaluating the seismic performance of different structural forms. Observations from past earthquakes suggest that tunnels located in seismically active regions should generally be designed to withstand ground shaking induced by wave propagation or permanent ground deformation caused by earthquakes. Seismic excitation along the longitudinal axis of a tunnel is expected to induce axial deformation and longitudinal bending, while transverse excitation may result in oval deformation or combined rocking and racking deformation modes [7,8]. Damage records from previous earthquakes demonstrate that, as a typical underground structure, the seismic performance of tunnels is not only directly influenced by ground motion intensity but is also closely related to the characteristics of the seismic motion. During the 6 February 2023 Kahramanmaraş earthquake, extremely high ground motion amplitudes were recorded; however, only 10 out of 87 tunnel structures suffered limited damage [9]. In contrast, the 1999 Chi-Chi earthquake caused damage to 49 tunnels within a radius of 60 km from the epicenter, including the collapse of three tunnels [6]. Some researchers have attributed this discrepancy to the effects of near-fault velocity pulses. Wang et al. [5] conducted a statistical analysis of data from 130 near-fault seismic stations within 55 km of the Chi-Chi earthquake epicenter and reported maximum peak ground velocity and displacement values of 292 cm/s and 867 cm, respectively. The associated pronounced velocity pulse was identified as a key feature of the Chi-Chi earthquake. The waveform exhibits a distinct pulse shape with steep rise or fall, mainly resulting from forward directivity effects and slip pulse mechanisms, which input a large amount of energy into structures within a short duration. This rapid energy concentration is difficult to dissipate in a short time, a phenomenon that has been verified by model tests. Chen et al. [10] observed through physical model experiments that the transverse dynamic responses of utility tunnels and surrounding soils subjected to near-fault ground motions with velocity pulses were significantly larger than those induced by ground motions without velocity pulses.

Recently, researchers have increasingly leveraged data-driven methodologies to bridge the gap between complex geological uncertainties and predictive structural responses. Focusing on seismic vulnerability, within the framework of performance-based earthquake engineering the relationship between intensity measures and structural failure is defined, and fragility curves are derived accordingly. Peak ground acceleration and peak ground velocity are commonly adopted as seismic intensity measures (IMs). In tunnel engineering, most related studies are based on nonlinear numerical simulations, in which the exceedance probabilities of different damage states are evaluated under a large number of ground motion inputs to construct fragility curves [11,12]. With advances in computational capability and data driven methodologies, machine learning techniques have gradually been introduced into the field of earthquake engineering for structural response prediction, damage assessment, and seismic risk analysis. Compared with traditional statistical approaches, machine learning models exhibit clear advantages in handling complex nonlinear relationships and high dimensional input features [13]. Among various machine learning techniques, ensemble learning improves model generalization ability and stability by combining the predictions of multiple base learners, and has been preliminarily applied in seismic ground motion prediction and structural fragility analysis. For instance, Wen et al. [14–17] implemented ensemble algorithms to assess the seismic fragility of building clusters, achieving higher accuracy than traditional regression methods. Liu et al. [18,19] developed advanced machine learning frameworks to analyze the stability mechanisms of underground rock engineering, demonstrating the potential of ensemble algorithms in handling nonlinear geomechanical data. In the field of tunnel engineering, Wang et al. [20] proposed a machine-learning-driven approach to analyze the feature importance of seismic parameters on tunnel damage and fragility prediction. Huang et al. [21] developed a seismic fragility analysis method for rock mountain tunnels based on support vector machines. Wang et al. [22] applied extreme gradient boosting and random forest algorithms to train and predict tunnel earthquake damage states. Yang et al. [23] developed a rapid resilience assessment model for tunnels crossing strike slip faults using the Adaptive Boosting algorithm.

Despite these significant advancements, few researchers have specifically addressed the seismic performance of horseshoe-shaped tunnels, which possess inherent geometric asymmetry, under the high-energy demands of near-fault velocity pulses. Most existing studies rely on Incremental Dynamic Analysis (IDA), where researchers scale a limited set of ground motions to various intensity levels. However, as noted by previous researchers, near-fault seismic signatures exhibit strong dependence on rupture distance. Simply scaling ground motions may bias the results by distorting these distance-dependent physical features.

To address these limitations, the researchers in this study established a comprehensive evaluation framework based on Multiple Stripes Analysis (MSA) and high-fidelity numerical datasets. By utilizing 1000 unscaled near-fault records (500 pulse-like and 500 non-pulse), this study preserves the natural coupling between fault distance and structural distortion. Furthermore, this study systematically compares three distinct ensemble learning architectures, Random Forest (RF), Extreme Gradient Boosting (XGBoost) and Categorical Boosting (CatBoost) to identify the most robust learner for predicting the nonlinear damage evolution of asymmetric tunnels. Finally, the Shapley Additive Explanations method (SHAP) is employed to quantify the physical contribution of intensity measures, providing a transparent link between data-driven predictions and seismic engineering mechanisms. The results of this study contribute to enhancing the accuracy and reliability of seismic fragility assessment for near-fault tunnels and provide valuable methodological support and technical insights for the seismic design and seismic risk assessment of tunnels located in near-fault regions.

2 Methodology

2.1 Basic theory of seismic fragility analysis

2.1.1 Definition of the fragility function

Derivation of fragility curves represents one of the key objectives in performance-based earthquake engineering. The detailed formulation of fragility curves has been extensively documented in previous studies [11,12,24]. In general, seismic fragility curves are commonly described using a lognormal probability distribution function:

(1)P\left[ds⩾d{s}_i|IM\right]=\mathrm{\Phi }\left[\frac{1}{{\beta }_{\text{tot}}}\ln \left(\frac{IM}{I{M}_{d{s}_i}}\right)\right],

where P denotes the probability that a specified damage state is exceeded for a given seismic intensity level. The parameter ds represents the structural damage index, while d{s}_i denotes the defined damage state. \mathit{\Phi } is the standard normal cumulative distribution function. IM is the seismic intensity measure, which may be represented by peak ground acceleration (PGA) or peak ground velocity (PGV). I{M}_{d{s}_i} denotes the median value of the seismic intensity measure required to induce the i damage state d{s}_i, as illustrated in Fig. 1. Linear regression analysis is commonly employed to establish the relationship between engineering demand parameters and seismic intensity measures. Given the threshold value of the engineering demand parameter corresponding to each damage state, the best estimate of the median threshold value expressed in terms of the seismic intensity measure I{M}_{d{s}_i} can be obtained based on the empirical relationship between the engineering demand parameter and the intensity measure. The parameter {\beta }_{\text{tot}} represents the total lognormal standard deviation associated with each damage state, reflecting the overall variability of the seismic response. The total lognormal standard deviation {\beta }_{\text{tot}} mainly accounts for three sources of uncertainty. Assuming that these uncertainties are independent and follow lognormal distributions, the total lognormal standard deviation {\beta }_{\text{tot}} can be expressed as:

(2){\beta }_{\text{tot}}=\sqrt{{\beta }_{\text{ds}}^2+{\beta }_{\text{c}}^2+{\beta }_{\text{d}}^2},

where {\beta }_{\text{ds}} characterizes the uncertainty associated with the threshold values of damage measures or parameters used to define different limit states, for which a recommended value of 0.4 is commonly adopted for tunnel structures. The parameter {\beta }_{\text{c}} represents the variability in structural properties, such as material strength, and a recommended value of 0.3 is suggested for tunnel structures [11,25]. The parameter {\beta }_{\text{d}} accounts for an additional source of significant uncertainty related to seismic demand. This final uncertainty component reflects the dispersion of structural damage indices calculated under different seismic intensity measures. It is quantified based on the scatter obtained from regression fitting of damage data [24] and can be expressed as:

(3){\beta }_{\mathrm{d}}=\sqrt{\frac{{\displaystyle \sum _{i=1}^n{\left[\ln D{M}_i-a\ln I{M}_i-b\right]}^2}}{n-2}}.

2.1.2 Multiple stripes analysis

Nonlinear dynamic time history analyses are conducted using the MSA approach. MSA is a methodology in which nonlinear dynamic analyses are performed using different sets of ground motion records selected at multiple predefined seismic intensity levels, referred to as strips or levels [26]. In conventional IDA, each ground motion record must be repeatedly scaled and analyzed until structural collapse is observed. Such an approach may not accurately represent the characteristics of real strong earthquakes, particularly for near-fault ground motions, whose intensity and spectral characteristics are closely correlated. Scaling near-fault ground motions may therefore distort their inherent features. In contrast, MSA allows the selection of different ground motion records at each intensity level that are most consistent with the site-specific probabilistic seismic hazard analysis characteristics. This feature enables MSA to more realistically capture the randomness of ground motions when evaluating near-fault effects or earthquakes with specific spectral characteristics, thereby providing a more reliable framework for nonlinear seismic response analysis.

2.1.3 Near-fault ground motion database and intensity measures

Near-fault ground motions are generally defined as strong ground motions recorded within approximately 30 km of the seismic source. One of the most prominent characteristics of near-fault ground motions is the velocity pulse effect [4]. Due to fault rupture directivity, seismic wave energy becomes highly concentrated along the direction of rupture propagation, resulting in large-amplitude, short-duration pulse-like motions in velocity time histories. In the frequency domain, near-fault ground motions are typically dominated by low-frequency components, with their response spectra exhibiting significant amplification in the intermediate to long period range. Moreover, the presence of velocity pulses leads to a highly concentrated and impulsive pattern of seismic energy input. Extensive damage investigations and strong-motion records have demonstrated that near-fault ground motions are among the primary factors responsible for severe damage to underground structures, particularly tunnel systems.

Considering the sample size requirements for machine learning databases and in order to reduce the influence of ground motion uncertainty while ensuring the objectivity of the derived fragility curves, a total of 500 near-fault pulse-like ground motion records and 500 non-pulse ground motion records were selected from the strong-motion database of the Pacific Earthquake Engineering Research Center (PEER) [13,27]. Pulse-like ground motions are identified using the wavelet-based methodology proposed by Baker [4], as categorized in the PEER NGA-West2 database. Earthquake magnitude M_{\text{w}}, fault distance R_{\text{rup}}, and fault mechanism were adopted as key selection criteria. Specifically, the selected ground motions correspond to M_{\text{w}} ranging from 6.5 to 8.0, R_{\text{rup}} within 30 km, and reverse or strike-slip fault mechanisms. In this study, seven seismic intensity measures, IM, were considered, as listed in Table 1. These IM are capable of capturing ground motion amplitude, cumulative energy input effects, structural dynamic characteristics, and pulse-related features, thereby providing a comprehensive representation of ground motion characteristics [28].

2.1.4 Damage measure

Damage measures are structural response quantities that can be used to predict damage to both structural and nonstructural components and systems. In the seismic fragility assessment of underground structures, the rational definition of damage measures serves as the key link between external seismic IM and structural DM. This relationship is commonly established through linear or nonlinear regression between DM and seismic IM. At present, strength-based failure criteria [11,24] and deformation-based failure criteria [12,28] are most frequently adopted for tunnel structures. Strength-based indicators mainly describe the load-carrying condition of the structure in the elastic range [11]. Once the material enters the nonlinear regime, however, internal force redistribution occurs, making strength-based indicators less effective in accurately reflecting the actual damage state. In contrast, deformation-based indicators such as relative drift angle, sectional curvature, and residual strain exhibit a more direct and monotonic relationship with structural damage states. These indicators are capable of continuously characterizing the entire damage evolution process, ranging from crack initiation and plastic hinge formation to ultimate collapse. Common deformation-based damage measures for tunnels include diameter distortion ratio, offset displacement, and joint opening. For the horseshoe-shaped tunnel investigated in this study, the sectional distortion ratio is selected as the primary macroscopic damage measure to characterize global structural damage and to capture shear-dominated deformation behavior induced by near-fault ground motions [7]. The sectional distortion ratio is defined as the ratio of the maximum sectional deformation of the tunnel lining under seismic loading to the original sectional dimension, as illustrated in Fig. 2, and can be expressed as:

(4)\phi =\frac{D^{\mathrm{\prime }}-D}{D}.

The seismic design code for underground structures in China, FEMA 356 issued by the United States Federal Emergency Management Agency, the Japanese Building Standard Law, and other relevant seismic design specifications define performance levels in different ways. At present, the most commonly adopted performance levels are generally classified into four categories, namely normal operation, post-repair operation, life safety, and complete damage. By establishing a multiscale damage quantification framework, continuous structural dynamic responses can be discretized into distinct damage states, including slight damage, moderate damage, severe damage, and collapse [20,22]. This classification enables a more accurate characterization of the brittle cracking and plastic through-damage processes of tunnel linings subjected to velocity pulse excitation, and provides a unified quantitative basis for subsequent probabilistic fragility curve fitting based on ensemble learning methods. The relationship between safety-oriented damage indicator thresholds and the corresponding performance levels and fortification objectives for each damage state is presented in Table 2.

2.2 Ensemble learning models

2.2.1 Fundamentals

Ensemble learning is a class of machine learning approaches that improve predictive performance and generalization capability by constructing and integrating multiple base learners. Its theoretical foundation originates from the bias-variance trade-off in statistical learning theory [29,30]. A single learner often suffers from inherent limitations in terms of bias or variance. By introducing model diversity and appropriately combining the outputs of multiple learners, ensemble learning can effectively reduce prediction variance without significantly increasing model bias, or enhance fitting capability while maintaining model stability. Consequently, ensemble learning demonstrates clear advantages in addressing problems characterized by strong nonlinearity, high parameter uncertainty, and significant data noise, and has been widely applied in structural engineering, earthquake engineering, and complex system response modeling [20,22,23,31]. According to the construction strategy of base learners and the corresponding error correction mechanisms, ensemble learning methods can generally be classified into two main categories. The first category is represented by Bagging-based algorithms, which are characterized by parallel model construction and variance reduction. These methods smooth prediction results in a statistical sense by repeatedly resampling the training data and independently training multiple models. The second category comprises Boosting-based algorithms, which are characterized by sequential model construction and bias reduction. These approaches progressively focus on samples that are difficult to fit, thereby improving overall predictive accuracy through iterative error correction. Random Forest, XGBoost, and CatBoost are representative implementations, and the basic architectures of the three algorithms are illustrated in Fig. 3.

2.2.2 Random Forest

RF is an ensemble learning method based on the Bagging paradigm [32]. Its core idea lies in constructing a large number of mutually independent decision trees and aggregating their outputs, thereby effectively reducing prediction variance. Specifically, RF first employs bootstrap sampling to generate multiple subsets from the original training dataset, with each subset used to train an individual decision tree. Owing to the differences in training samples, intrinsic diversity is naturally introduced among the base learners. During tree construction, RF further enhances diversity by randomly selecting a subset of features at each node split, rather than considering all available features. This feature subspace randomization mechanism significantly reduces the correlation among individual trees, enabling the ensemble model to preserve the low-bias property of single decision trees while substantially decreasing overall prediction variance. For regression tasks, the final output of RF is typically obtained by averaging the predictions of all decision trees. From a theoretical perspective, RF achieves a robust approximation of complex nonlinear mappings through a dual randomization mechanism that combines sample perturbation and feature perturbation. Its main advantages include relatively simple parameter tuning, robustness to outliers and noisy data, and the ability to quantitatively evaluate the contribution of input variables through feature importance measures. However, because there is no explicit error correction mechanism among the base learners, RF has limited capability to reduce model bias, which may restrict its adaptability in problems where high prediction accuracy is required.

2.2.3 XGBoost

XGBoost is a Boosting-based ensemble algorithm built upon gradient boosting decision trees. Its core principle lies in the sequential construction of weak learners, whereby each subsequent model focuses on correcting the prediction errors of the preceding one [33]. Unlike the independent and parallel modeling strategy adopted by RF the base learners in XGBoost exhibit a strict sequential dependency, and the overall model progressively approaches the true target function through iterative refinement. From a mathematical perspective, XGBoost formulates the training process as a regularized objective function minimization problem. The loss function quantifies prediction errors, while the regularization term constrains model complexity by penalizing factors such as the number of leaf nodes and the magnitude of leaf weights. By incorporating a second-order Taylor expansion of the loss function, XGBoost exploits both first-order and second-order gradient information during each iteration, thereby enhancing convergence efficiency and numerical stability. In practical implementation, XGBoost introduces several engineering optimizations, including parallelized feature splitting, automatic handling of missing values, and tree pruning strategies, which significantly improve computational efficiency and generalization performance. As a result, XGBoost demonstrates excellent predictive accuracy in problems characterized by high-dimensional feature spaces, complex nonlinear relationships, and moderate sample sizes. However, since Boosting algorithms are generally sensitive to noisy data, XGBoost typically requires careful parameter tuning and appropriate regularization to achieve optimal performance.

2.2.4 CatBoost

CatBoost is an advanced gradient boosting decision tree algorithm developed by Yandex [34,35]. Its principal advantages lie in the optimized handling of categorical features and the systematic mitigation of prediction shift. One of the key algorithmic innovations of CatBoost is the introduction of the ordered boosting strategy. By adopting a permutation scheme inspired by time series concepts, CatBoost ensures that the gradient computation for each sample depends only on preceding samples, thereby effectively suppressing model oscillation and improving generalization performance. In addition, CatBoost employs symmetric decision trees, also known as oblivious trees, as base learners, in which all nodes at the same depth share an identical splitting criterion. This tree structure acts as an implicit regularization mechanism and significantly accelerates prediction during the inference stage. At the feature processing level, CatBoost introduces an innovative target statistics encoding scheme. Categorical features are automatically transformed into numerical representations by combining prior information with ordered target statistics and moving average techniques, which effectively alleviates the curse of dimensionality associated with high-cardinality categorical variables. Compared with XGBoost, CatBoost exhibits superior robustness and strong out-of-the-box performance when dealing with datasets containing a large number of categorical features, making it particularly suitable for complex engineering problems with mixed data types.

3 Numerical modeling and dataset

3.1 Nonlinear dynamic numerical model

Based on the design parameters of a horseshoe-shaped tunnel from an actual engineering project, a nonlinear dynamic constitutive model of the tunnel–ground system was established in ABAQUS 2021 to simulate the seismic response of the tunnel, as shown in Fig. 4. The tunnel has a width of 13.7 m and a height of 11.5 m, with a lining thickness of 0.5 m. The burial depth of the tunnel is 50 m. The numerical model has a width of 170 m and a height of 180 m, and the distance from the tunnel sidewalls to the lateral model boundaries exceeds five times the tunnel width, so as to minimize boundary effects. The tunnel lining is composed of C40 concrete and is modeled using the concrete damage plasticity constitutive model [36], with an elastic modulus of 32.5 GPa and a Poisson’s ratio of 0.2. The remaining parameters are adopted from Shi et al. [37], as presented in Table 3. The surrounding rock mass is described using an elastoplastic constitutive model based on the Mohr–Coulomb strength criterion, with an elastic modulus of 900 MPa, a Poisson’s ratio of 0.32, a cohesion of 1.5 MPa, and an internal friction angle of 37°. Both the surrounding rock and the tunnel structure are discretized using solid elements. The maximum element size satisfies the following criterion [38]:

(5)\mathrm{\Delta }l\le \left(\frac{1}{8}\sim \frac{1}{10}\right){\lambda }_{min},

(6){\lambda }_{min}=\frac{v_{min}}{f_{max}},

where {\lambda }_{min} denotes the minimum wavelength of the input ground motions, v_{min} represents the minimum shear wave velocity of the surrounding rock mass, and f_{max} is the maximum frequency of the input ground motions. Rayleigh damping is adopted to construct the damping matrix of the surrounding rock, which consists of mass-proportional damping and stiffness-proportional damping, and can be expressed as:

(7)\left[C\right]=a_0\left[M\right]+a_1\left[K\right],

(8)\left[\begin{array}{c}{a}_0\\ a_1\end{array}\right]=\frac{2\zeta }{{\omega }_i+{\omega }_j}\left[\begin{array}{c}{\omega }_i\cdot {\omega }_j\\ 1\end{array}\right],

where \zeta denotes the damping ratio, which is taken as 5%. The parameter {\omega }_i represents the fundamental natural frequency of the model, while {\omega }_j corresponds to the predominant frequency of the input ground motions. This selection approach accounts for the combined effects of the structural natural frequency and the spectral characteristics of the input ground motions.

The computational procedure consists of three main stages, namely initial geostatic stress equilibrium, tunnel excavation with lining installation, and seismic dynamic time history analysis. During the geostatic stress equilibrium stage, the bottom boundary of the model is fixed, while the lateral boundaries are constrained in the normal direction to establish the initial in situ stress field. After tunnel excavation and lining support are completed, a second stress equilibrium analysis is performed to eliminate any residual imbalance and to prevent its influence on the subsequent dynamic analysis. During the dynamic time history analysis, the static boundaries are replaced by viscoelastic artificial boundaries to suppress seismic wave reflections. Equivalent reaction forces are applied at the model boundaries to maintain equilibrium with the in-situ stresses [39,40]. The seismic loads are introduced by transforming the ground motion records into equivalent nodal forces at the model boundaries using the wave input method, which can be expressed as:

(9)F_{\mathrm{b}}=(K_{\mathrm{b}}{u}_{\mathrm{b}}^{\mathrm{f}\mathrm{f}}+C_{\mathrm{b}}{\dot{u}}_{\mathrm{b}}^{\mathrm{f}\mathrm{f}}+{\sigma }_{\mathrm{b}}^{\mathrm{f}\mathrm{f}}n)A_{\mathrm{b}},

where K_{\text{b}} denotes the spring stiffness of the viscoelastic boundary, and C_{\text{b}} represents the corresponding damping coefficient. The terms u_{\text{b}}^{\text{ff}} and {\dot{u}}_{\text{b}}^{\text{ff}} denote the displacement and velocity of the free field, respectively, while {\sigma }_{\text{b}}^{\text{ff}} is the free-field seismic stress. The vector n represents the outward normal direction of the boundary, and A_{\text{b}} is the control area associated with the boundary node.

3.2 Characteristics of near-fault ground motions

The analysis of characteristic indices of near-fault ground motions is critical for elucidating structural damage mechanisms. Near-fault ground motions exhibit distinctive features in both the time and frequency domains, and their characterization provides important reference information for evaluating the predictive accuracy of numerical and data-driven models. Figure 5 illustrates the distributions of key characteristics for the selected 500 near-fault pulse-like ground motions and 500 non-pulse ground motions. In addition to the seismic intensity measures listed in Table 2, these characteristics include earthquake magnitude, fault distance, and pulse period. A comparison between pulse-like and non-pulse ground motions reveals that the mean values of magnitude, PGA, PGV, PGD, PGV/PGA, Arias intensity, CAV, and S_{\text{a}} are generally higher for pulse-like records. In particular, the energy-related intensity measures, namely CAV and Arias intensity, are markedly larger for pulse-like ground motions, indicating their substantially higher seismic energy content. R_{\text{rup}} shows a clear negative correlation with pulse characteristics: pulse-like ground motions are predominantly distributed within 10 km of the fault, whereas non-pulse ground motions are mainly concentrated in the range of 10–30 km. This trend reflects the attenuation of pulse effects with increasing R_{\text{rup}}. The correlation between M_{\text{w}} and the occurrence of velocity pulses is relatively weak. The modal magnitude of pulse-like ground motions is higher than that of non-pulse motions, suggesting that stronger earthquakes are more likely to generate pulse-type ground motions. Loh et al. [41] proposed that a PGV/PGA ratio greater than 0.2 can be used as an indicator for identifying near-fault pulse-like ground motions, as this ratio directly affects the boundary period between acceleration-sensitive and velocity-sensitive regions of the response spectrum. As shown in Fig. 5, pulse-like records are densely clustered in the higher range above 0.2, confirming their pronounced long-period characteristics and strong spectral energy concentration. The pulse periods of pulse-like ground motions span a wide range from 0.7 to 13 s, which further highlights their long-period characteristic. Overall, the selected records cover a broad range of frequencies, amplitudes, and pulse periods, enabling a comprehensive representation of tunnel seismic responses under near-fault ground motions.

Figure 6 illustrates the correlation heat maps among characteristic indices of near-fault pulse-like and non-pulse ground motions. The correlation coefficients quantitatively reveal the interdependence and potential redundancy among different parameters, providing valuable guidance for input feature selection in subsequent seismic fragility assessments based on ensemble learning methods. From the perspective of source parameters, R_{\text{rup}} exhibits a moderate negative correlation with PGA, PGV, and S_{\text{a}}, which is consistent with the classical attenuation law of ground motions with increasing distance. This correlation is more pronounced for pulse-like ground motions, indicating a stronger dependence of pulse characteristics on R_{\text{rup}}. Notably, the influence of T_{\text{p}} on ground motion intensity is particularly evident in the correlation analysis of pulse-like records, showing moderate to strong correlations with PGA, PGD, and the PGV/PGA ratio. This observation suggests that near-fault pulse effects predominantly govern the low- to intermediate-frequency components of ground motions, which are more effectively captured by displacement-related intensity measures. Clear differences are observed in the correlations among peak IM for pulse-like and non-pulse ground motions. Strong correlations are found between PGA and PGV, as well as between PGV and PGD, reflecting the inherent kinematic relationships between acceleration and velocity, and between velocity and displacement. In contrast, the relatively weak correlation between PGA and PGD for pulse-like ground motions further indicates a fundamental distinction between high-frequency acceleration components and long-period, large-displacement effects. The PGV/PGA ratio shows significant correlations with PGD and pulse period, while exhibiting a negative correlation with PGA. This confirms that PGV/PGA effectively characterizes the dominance of velocity pulses relative to high-frequency acceleration components, making it a robust indicator for identifying near-fault pulse-like ground motions and assessing their damage potential. Regarding energy-based intensity measures, Arias intensity and CAV display an extremely strong correlation for pulse-like ground motions. Both measures also exhibit high correlations with PGA and PGV, highlighting their ability to integrate the combined effects of ground-motion amplitude and duration. In addition, spectral acceleration S_{\text{a}} shows strong correlations with PGA, PGV, Arias intensity, and CAV, indicating that S_{\text{a}} effectively captures both amplitude characteristics and spectral properties of ground motions. The correlation analysis reveals that pulse-like ground motions are dominated by velocity- and energy-related parameters, whereas non-pulse-like motions exhibit strong coupling among traditional IM, highlighting the necessity of differentiated IM selection for tunnel seismic fragility assessment.

3.3 Relationship between intensity measure and damage measure

The influence of near-fault velocity pulses on the dataset is first examined from a seismic fragility perspective, which provides essential insight for interpreting the subsequent machine learning results. Under the loading conditions and surrounding rock properties considered in this study, PGA, PGV, and PGD are separately adopted as IMs to investigate the distribution of observed damage responses. Figure 7 illustrates the logarithmic relationship between the tunnel DM and PGA for pulse-like and non-pulse ground motions. In this figure, the horizontal axis represents lnPGA and the vertical axis denotes lnDM, together with the corresponding linear regression fits and coefficients of determination. From an overall perspective, both pulse-like and non-pulse ground motions exhibit a clear linear increasing trend of lnDM with lnPGA, indicating that PGA exerts a pronounced controlling effect on tunnel damage evolution. However, distinct differences are observed between the two motion types in terms of regression slope, data dispersion, and damage evolution characteristics, as summarized below. First, the coefficient of determination for pulse-like ground motions reaches 0.9336, demonstrating an excellent fitting quality with data points tightly clustered around the regression line. In contrast, the corresponding value for non-pulse ground motions is 0.8417, which is noticeably lower and reflects a higher degree of scatter. This indicates that tunnel damage exhibits a more stable and consistent relationship with PGA under pulse-like ground motions. Second, the regression slope for pulse-like ground motions is close to unity, suggesting a high sensitivity of the damage measure to variations in PGA. Even a modest increase in seismic intensity can result in a pronounced escalation of structural damage. This behavior is consistent with the amplification of structural response induced by velocity pulses and long-period effects inherent in near-fault pulse-like ground motions. By comparison, the regression slope for non-pulse ground motions is smaller, implying a more gradual increase in tunnel damage with increasing PGA. This contrast highlights the critical role of pulse effects in triggering severe tunnel damage. Third, in terms of damage state distribution, pulse-like ground motions are more likely to drive the tunnel into moderate, severe, or even collapse damage states at the same PGA level, whereas non-pulse ground motions tend to concentrate damage within the slight to moderate range. This observation suggests that conventional IM dominated by PGA may be adequate for non-pulse ground motions but may underestimate the actual damage potential under pulse-like motions. Overall, the results indicate that pulse-like ground motions not only intensify tunnel damage but also strengthen the correlation between damage and PGA, leading to a more deterministic and steeper damage evolution trend. These findings quantitatively demonstrate the significant amplification effect of near-fault pulse characteristics on tunnel seismic response and damage mechanisms. They further underline the necessity of distinguishing between pulse-like and non-pulse ground motions in near-fault seismic fragility analyses and of incorporating intensity measures capable of capturing low-frequency and pulse characteristics to improve the rationality and reliability of seismic risk assessments.

Figure 8 illustrates the logarithmic relationship between the DM and PGV. Under both pulse-like and non-pulse ground motions, lnDM exhibits a clear linear correlation with lnPGV, indicating that PGV is a reasonable indicator for characterizing the seismic damage evolution of tunnel structures. Nevertheless, pronounced differences are observed between the two motion types in terms of regression slope, data dispersion, and damage evolution characteristics, as detailed below. First, for pulse-like ground motions, the coefficient of determination is 0.7076. Although the data points generally align with the fitted regression line, the scatter is relatively large, reflecting the inherent uncertainty associated with velocity pulses and their amplification effects on structural response. In contrast, non-pulse ground motions yield a higher coefficient of determination of 0.8046, with data points more tightly clustered around the regression line. This indicates that, for non-pulse motions, PGV provides a more stable and consistent explanation of tunnel damage, and the governing mechanism of velocity effects on damage is comparatively simpler. Second, in terms of regression slope, the values for pulse-like and non-pulse ground motions are 0.9062 and 0.8570, respectively. The steeper slope associated with pulse-like motions implies that tunnel damage is more sensitive to changes in PGV when pronounced velocity pulses are present. In such cases, an identical increase in PGV may trigger substantially more severe structural damage. Last, From the perspective of damage-state distribution, pulse-like ground motions are more likely to drive the structure into severe damage states within the same PGV range. If near-fault pulse effects are neglected and PGV is used alone as the IM, the damage risk of tunnels subjected to pulse-like motions may be underestimated.

Figure 9 presents the relationship between the DM and PGD under pulse-like and non-pulse ground motions, together with the corresponding linear regression fits and coefficients of determination. Overall, a positive correlation between PGD and DM is observed for both motion types, indicating that tunnel damage tends to increase with displacement amplitude. However, compared with PGA and PGV, the correlation between PGD and DM is markedly weaker, as evidenced by the substantial data scatter and the limited quality of the regression fits. The coefficients of determination for pulse-like and non-pulse motions are only 0.1296 and 0.3684, respectively, with a large proportion of data points deviating from the fitted lines. Furthermore, the boundaries between different damage states in the PGD domain are poorly defined, particularly for pulse-like motions, where significant overlap among damage levels is observed. This demonstrates that, although displacement-based indicators have clear physical relevance for describing permanent deformation and large-scale structural response, they are unable to adequately capture the rapid energy input and velocity-pulse-dominated mechanisms that govern damage evolution under near-fault ground motions. Consequently, PGD alone is insufficient for effective discrimination of tunnel damage in near-fault conditions.

In summary, it can be concluded that PGA and PGV exhibit relatively stable and consistent correlations with tunnel damage across different ground-motion types, whereas PGD shows limited explanatory capability. Therefore, under near-fault ground motions, particularly in the presence of pronounced pulse effects, PGA and PGV are more suitable as key seismic intensity measures for tunnel fragility analysis.

3.4 Fragility curves

Figure 10 presents the seismic fragility curves of the tunnel under pulse-like and non-pulse ground motions for different damage states, including no damage, slight damage, moderate damage, severe damage, and collapse. The horizontal axis denotes PGA, while the vertical axis represents the probability of exceeding a given damage state. A comparison of the fragility curves under the two types of ground motions provides direct insight into the influence of near-fault pulse effects on tunnel seismic performance. Under pulse-like ground motions, the fragility curves for all damage states exhibit a noticeably steeper ascending trend. As PGA increases, the probabilities of exceeding the no-damage and slight-damage states rapidly approach unity even at relatively low PGA levels, indicating that pulse-like motions are highly effective in triggering initial tunnel damage. For moderate, severe, and collapse damage states, the corresponding fragility curves are significantly elevated within the PGA range of 0.4g–0.8g, demonstrating that, at the same PGA level, pulse-like ground motions are more likely to induce higher levels of structural damage. This behavior reflects the pronounced deformation demands imposed by the large-energy velocity pulses and long-period components inherent in near-fault pulse-like ground motions. In contrast, the fragility curves under non-pulse ground motions are consistently shifted downward, particularly for moderate and higher damage states. When PGA is below 0.4g, tunnel damage is predominantly limited to no damage or slight damage, and the probability of exceeding moderate or more severe damage states increases only gradually. Even at higher PGA levels, the exceedance probabilities associated with severe damage and collapse remain relatively low, indicating that non-pulse ground motions exert a comparatively milder destructive effect on tunnel structures.

A direct comparison between the two motion types reveals that pulse-like ground motions substantially increase the likelihood of reaching higher damage states at the same PGA level, resulting in a pronounced leftward shift of the fragility curves toward lower PGA values. This discrepancy clearly indicates that relying solely on PGA as an intensity measure, without distinguishing between pulse-like and non-pulse ground motions, may lead to an underestimation of tunnel seismic risk under near-fault conditions. Therefore, it is essential to explicitly account for ground motion type in tunnel seismic fragility assessments and to incorporate intensity measures or correction approaches capable of capturing pulse effects, in order to enhance the reliability and robustness of the evaluation results.

4 Results of ensemble learning models

4.1 Training of ensemble learning models

Seismic IMs, together with M_{\text{w}}, R_{\text{rup}}, and T_{\text{p}}, are integrated into a unified dataset for training the ensemble learning models. Prior to model training, the seismic response dataset of the horseshoe-shaped tunnel is subjected to data cleaning and normalization. The dataset is then randomly divided into training and testing subsets with a ratio of 7:3. To fully exploit the information contained in the limited number of ground motion records, a five-fold cross-validation strategy is adopted during the training process. Specifically, the training set is evenly partitioned into five subsets, and each subset is alternately used as a validation set to evaluate the generalization performance of the model under the current hyperparameter configuration. During each training iteration, the base learners iteratively fit the complex nonlinear relationships between near-fault seismic intensity measures and tunnel damage indicators, such as lining strain and deformation ratio. Model biases are continuously corrected by minimizing the loss function through gradient-based updates, enabling the ensemble models to accurately capture the response characteristics of structurally vulnerable regions in horseshoe-shaped tunnels subjected to dynamic loading [13].

Model performance is highly sensitive to hyperparameter selection. Given the strong interactions among key parameters such as learning rate, tree depth, and the number of leaf nodes in ensemble learning algorithms, conventional grid search is not adopted in this study. Instead, Bayesian optimization is employed due to its superior computational efficiency. This approach constructs a Gaussian process surrogate model to approximate the relationship between hyperparameter combinations and model performance, and uses an acquisition function to balance global exploration and local exploitation. As a result, the optimal hyperparameter set that minimizes the mean squared error is identified within a limited number of iterations. Considering the high variability inherent in near-fault seismic data, particular emphasis is placed on tuning the L1 and L2 regularization coefficients to enhance model robustness when switching between pulse-like and non-pulse ground motions. The optimized models are finally validated through independent testing, demonstrating their strong reliability and predictive capability in the seismic fragility assessment of horseshoe-shaped tunnels. The following is a summary analysis of the parameters of various machine learning models as well as the best parameter scheme (Tables 4–6).

4.2 Model prediction performance

To systematically evaluate the performance of ensemble learning regression models in predicting the seismic fragility of near-fault tunnels, multiple statistical metrics are employed to quantitatively assess both prediction accuracy and generalization capability. By comparing model predictions with numerical simulation results, the fitting performance of each model across different levels of seismic intensity is comprehensively examined. The primary evaluation metrics include the mean absolute error (MAE), mean squared error (MSE), root mean squared error (RMSE), and the coefficient of determination (R²). Among these metrics, the R² is adopted as the principal indicator to quantify the proportion of variance in the target variable that can be explained by the model. An R² value approaching unity indicates a strong predictive capability and an excellent agreement between predicted and simulated tunnel damage responses. In addition, RMSE provides an intuitive measure of the average deviation between predicted and actual values, which has clear physical significance and practical relevance in engineering applications.

(10)RMSE=\sqrt{\frac{1}{n}\sum _1^n{(y_i-{\hat{y}}_i)}^2},

(11)R^2=1-\frac{{\displaystyle \sum _1^n{(y_i-{\hat{y}}_i)}^2}}{{\displaystyle \sum _1^n{(y_i-\overline{y})}^2}},

where y_i denotes the reference values obtained from numerical simulations, {\hat{y}}_i represents the corresponding predictions generated by the ensemble learning models, \overline{y} is the mean of the observed sample values, and n is the total number of samples.

Figure 11 presents a comparative visualization of the regression prediction performance of RF, XGBoost, and CatBoost, for datasets corresponding to pulse-like and non-pulse ground motions. In the figure, the horizontal axis represents the reference DM obtained from numerical simulations, while the vertical axis denotes the values predicted by the ensemble learning models. The red dashed line indicates the ideal one-to-one prediction, the shaded gray band represents a ±15% error margin, and the R^2 together with the RMSE are reported as quantitative performance metrics. For pulse-like ground motions, all three models successfully capture the overall linear relationship between predicted and reference values. The corresponding R² values are 0.9211, 0.9419, and 0.9303, with RMSE values of 0.4029, 0.3458, and 0.3786 for RF, XGBoost, and CatBoost, respectively. Most data points are distributed close to the ideal prediction line, and the majority fall within the ±15% error bounds. Among the three models, XGBoost exhibits the strongest fitting capability for tunnel damage prediction under pulse-like ground motions. Accordingly, the prediction accuracy of the ensemble models under pulse-like conditions follows the order XGBoost > CatBoost > RF. For non-pulse ground motions, the three models also demonstrate a strong linear correlation between predicted and reference values. The corresponding R² values are 0.9182, 0.9484, and 0.9581. In this case, the CatBoost model performs best, achieving both the highest R² and the lowest RMSE, which indicates superior generalization capability and prediction stability under non-pulse ground motions. The prediction accuracy ranking under non-pulse conditions is CatBoost > XGBoost > RF. From the scatter distribution, the prediction points associated with non-pulse ground motions are more tightly clustered, with the majority concentrated near the ideal fitting line at small distortion levels. This suggests that severe tunnel damage is less likely under non-pulse ground motions and that the mapping between seismic IMs and structural damage is more stable, which is more favorable for data-driven model training and generalization.

A comprehensive comparison of the two datasets reveals that all models achieve higher prediction accuracy and lower error levels under non-pulse ground motions, with RMSE values consistently below 0.3. In contrast, RMSE values for pulse-like ground motions exceed 0.3 for all models. This observation reflects the higher seismic intensity and increased response complexity associated with pulse-like ground motions, which substantially amplify tunnel damage and impose greater challenges for regression modeling. Moreover, the performance discrepancies among the models across the two motion types highlight their differing capabilities in capturing nonlinear features and complex variable interactions.

Overall, all three ensemble learning models demonstrate strong predictive capability for tunnel seismic damage indicators, although subtle differences in performance are observed. Specifically, CatBoost exhibits the best overall performance under non-pulse ground motions, whereas XGBoost shows superior fitting ability under pulse-like ground motions. These results confirm the effectiveness of ensemble learning approaches for tunnel seismic fragility assessment and underscore the necessity of distinguishing between pulse-like and non-pulse ground motions in data-driven modeling and analysis.

In Subsection 3.4, fragility curves for pulse-like and non-pulse ground motions were constructed based on the results of 500 numerical simulations for each case. In the present Subsection, fragility curves are generated using 100 data points predicted by three ensemble learning models, namely RF, XGBoost, and CatBoost, as shown in Fig. 12. From the perspective of damage states, the prediction accuracy is highest at the slight damage level, where the predicted points almost coincide with the reference curves. At the severe damage and collapse levels, some dispersion is observed due to the reduced number of samples and increased response uncertainty. Nevertheless, the predicted fragility curves still preserve the correct monotonic increasing trend, without pronounced overestimation or underestimation. This suggests that the models retain an acceptable level of predictive reliability even at high damage states.

Among the three models, CatBoost and XGBoost exhibit predictions that are generally closer to the reference fragility curves, with more concentrated scatter distributions. The RF model shows slightly larger dispersion at certain damage levels, although its overall trend remains consistent with the reference curves. Taken together, these results indicate that ensemble learning models can predict structural fragility curves with high accuracy and successfully capture the nonlinear evolution characteristics associated with different ground motion types and damage states. This confirms the effectiveness and engineering applicability of machine learning methods for seismic fragility assessment of tunnel structures.

4.3 Shapley additive explanations

SHAP is a model interpretability approach grounded in game theory, specifically based on Shapley values, and is designed to quantitatively characterize the marginal contribution of each input feature to the model prediction [15,42]. The core concept of SHAP is to interpret the prediction process of a machine learning model as a cooperative game, in which each feature is treated as a player and the model output is regarded as the total payoff of the game. By evaluating the marginal gain associated with different combinations of features, SHAP assigns a unique contribution value to each feature, thereby providing a fair and additive explanation of the model output. In engineering and scientific applications, SHAP analysis facilitates the interpretation of data-driven predictions by revealing the underlying physical or mechanical mechanisms embedded in the model. By examining the distribution patterns and variation trends of SHAP values, key influencing factors and their threshold effects can be systematically identified. This information provides valuable guidance for model refinement, parameter sensitivity analysis, and engineering decision-making. Consequently, SHAP has become a widely adopted and powerful tool for model interpretability, particularly well-suited for mechanism exploration and result validation in high-dimensional and strongly nonlinear problems.

Figure 13 presents the SHAP-based feature importance results of different ensemble learning models under pulse-like and non-pulse ground motion conditions. The horizontal axis represents the mean absolute SHAP value, which quantifies the global contribution of each seismic intensity measure to the prediction of tunnel seismic damage. A larger SHAP value indicates a more pronounced influence of the corresponding feature on the model output. Under pulse-like ground motions, the three models exhibit a highly consistent hierarchy of dominant features. PGA emerges as the most influential parameter, with mean SHAP values significantly exceeding those of the other features. This finding indicates that strong transient inertial effects remain the primary controlling factor governing structural response under pulse-like excitation. In addition, PGV, Arias intensity, and S_{\text{a}} are consistently ranked among the most important secondary features, reflecting the pronounced low-frequency velocity pulses and concentrated energy input characteristic of pulse-like ground motions, which play a critical role in the accumulation of structural damage. For non-pulse ground motions, the distribution of feature importance becomes more uniform. Although all three models still identify PGA as the most important predictor, its relative dominance is noticeably reduced compared with the pulse-like case, suggesting that structural response under non-pulse excitation is not solely governed by peak inertial forces. The SHAP values associated with Arias intensity, CAV, M_{\text{w}}, and R_{\text{rup}} are generally lower, indicating that energy concentration and near-fault effects exert a comparatively limited influence on damage prediction for non-pulse ground motions.

From the perspective of ensemble learning algorithms, the overall trends in feature importance rankings are highly consistent across the three models, demonstrating the robustness of the SHAP-based interpretations. The XGBoost model exhibits the strongest dependence on PGA, whereas RF and CatBoost allocate feature importance more evenly among secondary indicators such as PGV and S_{\text{a}}. Cross-validation of SHAP results across multiple ensemble models enhances the reliability of the conclusions and further highlights the necessity of employing multiple ensemble learning approaches in the seismic fragility assessment of tunnel structures.

Figure 14 presents SHAP summary plots illustrating the seismic feature importance derived from different ensemble learning models under pulse-like and non-pulse ground motion conditions. Compared with Fig. 13, the summary plots provide a more detailed visualization of the SHAP value associated with each individual data point, thereby indicating the positive or negative contribution of each feature to the predicted tunnel damage measure. The color scale ranges from blue to red, representing low to high feature values, respectively. For the pulse-like ground motion dataset, PGA, PGV, and S_{\text{a}} consistently rank as the most influential features across all three models, indicating that seismic IM remain the dominant factors governing tunnel seismic response and damage evolution under near-fault pulse-like motions. PGA exhibits the widest SHAP value distribution, with high feature values predominantly corresponding to positive SHAP values, demonstrating that strong PGA significantly increases the probability of tunnel damage. PGV and S_{\text{a}} follow closely, highlighting the importance of velocity pulses and their interaction with the structural dynamic characteristics of the tunnel. The T_{\text{p}} and the PGV/PGA ratio show notable contributions in the XGBoost and CatBoost models, suggesting that these algorithms are particularly sensitive to spectral characteristics and energy input features of ground motions. Arias intensity and CAV, which represent cumulative energy effects, exhibit SHAP values concentrated mainly in the positive range, confirming that energy-based intensity measures play a non-negligible role in promoting tunnel damage. For the non-pulse ground motion dataset, PGA and PGV remain the most important features; however, their SHAP value ranges are significantly narrower than those observed under pulse-like conditions. This indicates that, in the absence of pronounced velocity pulses, the destructive effects of ground motions on tunnels are relatively milder. The importance rankings of M_{\text{w}} and R_{\text{rup}} increase noticeably in the XGBoost and CatBoost models, suggesting that, under non-pulse excitation, the influences of conventional intensity measures and source parameters become more balanced. Overall, the results demonstrate that, under pulse-like ground motions, velocity- and energy-related indicators such as PGV, S_{\text{a}}, Arias intensity, and T_{\text{p}} gain substantially greater importance, confirming that near-fault pulse effects constitute a key controlling mechanism that amplifies seismic damage in horseshoe-shaped tunnels.

With respect to model-to-model comparisons, XGBoost exhibits the widest SHAP value distribution under pulse-like ground motions, indicating a higher sensitivity to high-intensity seismic scenarios. CatBoost and RF produce a more balanced feature ranking, assigning non-negligible importance to parameters beyond PGA, S_a, and PGV, such as M_w and R_{\text{rup}}. Under non-pulse ground motions, CatBoost yields the most compact SHAP value distribution, reflecting greater prediction stability. XGBoost remains sensitive to high feature values, although with reduced variability compared to the pulse-like case, while RF exhibits the smallest dispersion of feature contributions. These findings demonstrate that adopting a multi-model comparative framework is essential for enhancing the reliability of tunnel fragility assessment. The SHAP-based interpretation results further provide a sound theoretical basis for ground motion selection and parameter control in the seismic design of tunnels located in near-fault regions.

5 Discussion

The ensemble learning framework developed in this study not only achieves superior statistical predictive accuracy but also elucidates the complex nonlinear seismic response mechanisms inherent to horseshoe-shaped tunnels under near-fault ground motions. Methodologically, the MSA framework offers a distinct advantage over IDA. While IDA relies on scaling ground motions, which disrupts the intrinsic fault-distance dependency and spectral signatures, MSA utilizes unscaled records to preserve the physical coupling between source parameters and site response. This ensures that the derived fragility curves reflect the genuine seismological characteristics of near-fault regions, providing a more rigorous link between machine learning predictions and structural mechanics.

Practically, this framework enables real-time seismic damage screening when integrated with Structural Health Monitoring systems, allowing for rapid post-earthquake assessment based on IMs. Furthermore, its computational efficiency allows designers to perform large-scale probabilistic risk scanning to optimize liner support in high-risk zones.

However, certain limitations remain. The high predictive accuracy primarily indicates the model’s fidelity in emulating finite element simulations; thus, external validation using post-earthquake reconnaissance data is necessary. Additionally, the model’s generalizability to diverse geological conditions or different tunnel scales should be further extended via transfer learning. Finally, incorporating uncertainties in geotechnical parameters will be essential for a more comprehensive probabilistic risk assessment in future studies.

6 Conclusions

This study evaluates the seismic fragility of near-fault horseshoe-shaped tunnels using MSA and three ensemble learning algorithms. A total of 500 near-fault pulse-like ground motions and 500 non-pulse ground motions were selected, and a dataset linking seismic IMs with tunnel DM was constructed through nonlinear numerical simulations to train the machine learning models. Model interpretability was further investigated using SHAP analysis. The main conclusions can be summarized as follows.

1) Near-fault pulse-like ground motions induce significantly more severe damage to tunnel structures than non-pulse ground motions. When PGA is adopted as the intensity measure, the probabilities of moderate and severe damage under pulse-like ground motions are markedly higher than those under non-pulse ground motions.

2) The three ensemble learning models, RF, XGBoost, and CatBoost, demonstrate good capability in predicting tunnel structural damage and are able to accurately reconstruct seismic fragility curves. The CatBoost model achieves superior predictive performance, effectively capturing the complex nonlinear mapping between ground motion IMs and tunnel DM.

3) SHAP-based interpretability analysis reveals that PGA and PGV are the most influential indicators governing tunnel damage states. Under pulse-like ground motions, Arias intensity, and S_{\text{a}} are considered next-level important features, reflecting the dominant effects of low-frequency velocity pulses and concentrated energy input on cumulative structural damage. In contrast, feature importance under non-pulse ground motions is more evenly distributed among the considered parameters.

References

[1] Zhang D L , Sun Z Y , Fang Q . Scientific problems and research proposals for Sichuan–Xizang railway tunnel construction. Underground Space, 2022, 7(3): 419–439

[2] Zhu H H , Yu H T , Han F Q , Wei Y B , Yuan Y. . Seismic resilience design principles and key issues for tunnels crossing active faults.. China Journal of Highway and Transport, 2023, 36(11): 193–204

[3] Zhu H H , Yan J X , Liang W H . Challenges and development prospects of ultra-long and ultra-deep mountain tunnels. Engineering, 2019, 5(3): 384–392

[4] Baker J W . Quantitative classification of near-fault ground motions using wavelet analysis. Bulletin of the Seismological Society of America, 2007, 97(5): 1486–1501

[5] Wang G Q , Zhou X Y , Zhang P Z , Igel H . Characteristics of amplitude and duration for near fault strong ground motion from the 1999 Chi-Chi, Taiwan Earthquake. Soil Dynamics and Earthquake Engineering, 2002, 22(1): 73–96

[6] Wang T T , Kwok O L A , Jeng F S . Seismic response of tunnels revealed in two decades following the 1999 Chi-Chi earthquake (Mw 7. 6) in Taiwan: A review. Engineering Geology, 2021, 287: 106090

[7] Tsinidis G , de Silva F , Anastasopoulos I , Bilotta E , Bobet A , Hashash Y M A , He C , Kampas G , Knappett J , Madabhushi G . et al. Seismic behaviour of tunnels: From experiments to analysis. Tunnelling and Underground Space Technology, 2020, 99: 103334

[8] Hashash Y M A , Hook J J , Schmidt B , I-Chiang Yao J . Seismic design and analysis of underground structures. Tunnelling and Underground Space Technology, 2001, 16(4): 247–293

[9] Apostolaki S , Karahan S , Riga E , Tsinidis G , Gokceoglu C , Pitilakis K . Seismic performance of tunnels and verification of available seismic risk models for the 2023 Kahramanmaraş earthquakes. Tunnelling and Underground Space Technology, 2025, 156: 106185

[10] Chen H J , El Naggar M H , Chu J , Li X J , He Q M , Wang L H , Liu X W , Zhou L Y . Transverse response of utility tunnel under near fault ground motions: Multi-shake table array tests. Soil Dynamics and Earthquake Engineering, 2023, 174: 108135

[11] Argyroudis S A , Pitilakis K D . Seismic fragility curves of shallow tunnels in alluvial deposits. Soil Dynamics and Earthquake Engineering, 2012, 35: 1–12

[12] Andreotti G , Lai C G . Use of fragility curves to assess the seismic vulnerability in the risk analysis of mountain tunnels. Tunnelling and Underground Space Technology, 2019, 91: 103008

[13] Nguyen H D , Lee Y J , LaFave J M , Shin M . Seismic fragility analysis of steel moment frames using machine learning models. Engineering Applications of Artificial Intelligence, 2023, 126: 106976

[14] Wen H J , Hu J W , Xiong F G , Zhang C , Song C H , Zhou X Z . A random forest model for seismic-damage buildings identification based on UAV images coupled with RFE and object-oriented methods. Natural Hazards, 2023, 119(3): 1751–1769

[15] Wen H J , Liu B , Di M R , Li J Y , Zhou X Z . A SHAP-enhanced XGBoost model for interpretable prediction of coseismic landslides. Advances in Space Research, 2024, 74(8): 3826–3854

[16] Wen H J , Wu J N , Zhang C , Zhou X Z , Liao M Y , Xu J H . Hybrid optimized RF model of seismic resilience of buildings in mountainous region based on hyperparameter tuning and SMOTE. Journal of Building Engineering, 2023, 71: 106488

[17] Wen H J , Zhou X Z , Zhang C , Liao M Y , Xiao J F . Different-classification-scheme-based machine learning model of building seismic resilience assessment in a mountainous region. Remote Sensing, 2023, 15(9): 2226

[18] Liu B , Wen H J , Di M R , Huang J H , Liao M Y , Yu J Y , Xiang Y T . Mapping and interpretability of aftershock hazards using hybrid machine learning algorithms. Journal of Rock Mechanics and Geotechnical Engineering, 2025, 17(8): 4908–4932

[19] Liu B , Wen H J , Di M R , Liao M Y , Huang J H , Qian L , Chu Y B . Robustness study of seismic hazard assessment models based on multi-source data. Gondwana Research, 2025, 144: 87–108

[20] Wang Q , Geng P , Wang L J , He D W , Shen H M . Machine learning-driven feature importance appraisal of seismic parameters on tunnel damage and seismic fragility prediction. Engineering Applications of Artificial Intelligence, 2024, 137: 109101

[21] Huang G , Qiu W G , Zhang J R . Modelling seismic fragility of a rock mountain tunnel based on support vector machine. Soil Dynamics and Earthquake Engineering, 2017, 102: 160–171

[22] Wang L J , Geng P , Chen J B , Wang T Q . Machine learning-based fragility analysis of tunnel structure under different impulsive seismic actions. Tunnelling and Underground Space Technology, 2023, 133: 104953

[23] Yang L J , Xin C L , Wang Z , Yu X Y , Hajirasouliha I , Feng W K . Seismic resilience analysis of high-speed railway tunnels across fault zones using ensemble learning approach. Underground Space, 2025, 25: 99–131

[24] Qiu W G , Huang G , Zhou H C , Xu W P . Seismic vulnerability analysis of rock mountain tunnel. International Journal of Geomechanics, 2018, 18(3): 04018002

[25] Argyroudis S A , Mitoulis S Α , Winter M G , Kaynia A M . Fragility of transport assets exposed to multiple hazards: State-of-the-art review toward infrastructural resilience. Reliability Engineering & System Safety, 2019, 191: 106567

[26] Baker J W . Efficient analytical fragility function fitting using dynamic structural analysis. Earthquake Spectra, 2015, 31(1): 579–599

[27] Nguyen D D , Park D , Shamsher S , Nguyen V Q , Lee T H . Seismic vulnerability assessment of rectangular cut-and-cover subway tunnels. Tunnelling and Underground Space Technology, 2019, 86: 247–261

[28] Zhang X Y , Pan Y Y , Su L , Ou E T , Liu H , Liu C , Cui J . Seismic fragility assessment of shield tunnels in liquefiable soil-rock strata using fuzzy method for IM optimization. Tunnelling and Underground Space Technology, 2024, 152: 105957

[29] Dietterich T G. Ensemble methods in machine learning. In: The 1st International Workshop on Multiple Classifier Systems. Cagliari: Springer, 2000, 1–15

[30] Nordhausen K . Ensemble methods: Foundations and algorithms by Zhi-Hua Zhou. International Statistical Review, 2013, 81(3): 470

[31] Wang J D , Huang L C , Liang Y , Dong J , Sun W . Transfer learning enhanced physics-informed neural network for buried pipeline deformation analysis under permanent ground deformation. Computers and Geotechnics, 2026, 189: 107630

[32] Ho T K . The random subspace method for constructing decision forests. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1998, 20(8): 832–844

[33] Chen T QGuestrin C. XGBoost: A scalable tree boosting system. In: Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining. San Francisco: ACM, 2016, 785–794

[34] Prokhorenkova LGusev GVorobev ADorogush A VGulin A. CatBoost: Unbiased boosting with categorical features. In: Proceedings of the 32nd International Conference on Neural Information Processing Systems. Montréal: Curran Associates Inc., 2018, 6639–6649

[35] Dorogush A V, Ershov V, Gulin A. CatBoost: Gradient boosting with categorical features support. 2018, arXiv: 1810.11363

[36] Lee J , Fenves G L . Plastic-damage model for cyclic loading of concrete structures. Journal of Engineering Mechanics, 1998, 124(8): 892–900

[37] Shi C , Tao L J , Ding P , Wang Z G , Jia Z B . Study on seismic response characteristics and failure mechanism of giant-span flat cavern. Tunnelling and Underground Space Technology, 2023, 140: 105328

[38] Lysmer J , Kuhlemeyer R L . Finite dynamic model for infinite media. Journal of the Engineering Mechanics Division, 1969, 95(4): 859–877

[39] Du X L , Zhao M , Wang J T. . Stress artificial boundary in FEA for near-field wave problem.. Chinese Journal of Theoretical and Applied Mechanics, 2006, 38(1): 49–56

[40] Liu J B , Gu Y , Li B , Wang Y . An efficient method for the dynamic interaction of open structure-foundation systems. Frontiers of Architecture and Civil Engineering in China, 2007, 1(3): 340–345

[41] Loh C H , Wan S , Liao W I . Effects of hysteretic model on seismic demands: Consideration of near-fault ground motions. The Structural Design of Tall Buildings, 2002, 11(3): 155–169

[42] Lundberg S MLee S I. A unified approach to interpreting model predictions. In: Proceedings of the 31st International Conference on Neural Information Processing Systems. Long Beach: Curran Associates Inc., 2017, 4768–4777

Rights & permissions

The Author(s). This article is published with open access at link.springer.com and journal.hep.com.cn