Application of a Deep Learning Method to the Seismic Vulnerability Analysis of Cross-Fault Hydraulic Tunnels Based on MLE-IDA

Abstract

Rapidly developed deep learning methods, widely used in various fields of civil engineering, have provided an efficient option to reduce the computational costs and improve the predictive capabilities. However, it should be acknowledged that the application of deep learning methods to develop prediction models that efficiently assess the nonlinear dynamic responses of cross-fault hydraulic tunnels (CFHTs) is lacking. Thus, the objective of this study is to construct a rational artificial neural network (ANN) prediction model to generate the mass data and fragility curves of CFHTs. Firstly, an analysis of 1080 complete nonlinear dynamic time histories via incremental dynamic analysis (IDA) is conducted to obtain the mass data of the drift ratio of the CFHT. Then, the hyper-parameters of the ANN model are discussed to determine the optimal parameters based on four examined approaches to improve the prediction capacity and accuracy. Meanwhile, the traditional probabilistic seismic demand models of the predicted values obtained by the ANN model and the numerical results are compared with the statistical parameters. Eventually, the maximum likelihood estimation couping IDA method is applied to assess the seismic safety of CFHTs under different damage states. The results show that two hidden layers, ten neurons, and the ReLU activation function for the ANN model with Bayesian optimization can improve the reliability and decrease the uncertainty in evaluating the structural performance. Moreover, the amplitude of the seismology features can be used as the neurons to build the input layers of the ANN model. It is found through vulnerability analysis that the traditional seismic fragility analysis method may overestimate the earthquake resistance capacity of CFHTs compared with maximum likelihood estimation. In practical engineering, ANN methods can be regarded as an alternative approach for the seismic design and performance improvement of CFHTs.

1. Introduction

As one of the components of lifeline infrastructure systems for civil engineering, long-distance water diversion projects have been established around the globe for decades to reduce the disparity between the water supply and demand. In general, hydraulic tunnels (HTs) are the most common subterranean constructions for water conveyance infrastructure, whose primary advantage is to reduce engineering distances and circumvent geographical limitations. However, some practical engineering sites are inherently positioned in high-intensity areas due to the requirement to cross mountains and rivers, which imposes higher seismic design standards on HTs. Note that a considerable proportion of tunnels, particularly those crossing fault fracture zones, still suffer from severe damage during strong seismic events. More specifically, the seismic performance of underground structures located in the fault fracture zone is significantly poorer than that of underground structures with good geological conditions, due to the circulation dislocation and shear deformation of the fault fracture zone. For example, the severe damage to the Dagang CFHT during the Chi-Chi earthquake demonstrates that seismic activity is more common in regions with poor surrounding rock characteristics and abrupt shifts between surrounding rock masse layers [1]. By reviewing different case histories [2,3,4,5,6], it also can be inferred that the significant task for the civil engineering community is to ensure the long-term and safe operation of CFHTs in every country and society. In light of this information, it can be concluded that civil engineers and property owners must pay full attention to the seismic performance of CFHTs in regions with frequent seismic activity during seismic design.

With the widespread application of fragility analysis in civil engineering, the cloud analysis [7], incremental dynamic analysis (IDA) [8], and multiple stripe analysis [9,10] methodologies are gradually being accepted and playing irreplaceable roles in establishing robust seismic fragility curves after selecting the optimal intensity measure (IM). All three typical analytical methods require a large number of structural response results to derive the statistical seismic IM–demand measure (DM) relation [11]. In particular, several thousand heavy nonlinear time history finite element (FE) numerical simulations using a series of ground motions (GMs) are required to directly obtain the engineering demand parameters (EDPs) of underground structures. However, the long computational time and great effort required for dynamic analysis make the approach inefficient and difficult to adopt in practice [12]. In order to tackle these issues, the simple linear regression of the IM–DM relation in the logarithmic space is commonly utilized to reduce the computational costs and develop appropriate probabilistic seismic demand models (PSDMs) for underground structures, based on the work of Cornell et al. [8]. Notably, the normal distribution condition of the structural response and the linear function of IM-EDP in the logarithmic space are the primary sources of error in seismic fragility analysis [13]. Furthermore, a structure under strong GM excitation will experience nonlinear behavior, resulting in a less efficient simple linear regression. In modern civil engineering problems, another way to improve the computational efficiency is to calibrate the statistical relationships between seismic inputs and EDPs using data-driven deep learning (DL) algorithms.

There is a tendency to use DL algorithms as a valuable alternative to alleviate the computational burden in seismic safety assessment. Recently, Huang et al. [14] implemented artificial neural networks (ANNs) into the probabilistic framework for FRP-strengthened RC columns to esimtate the seismic resistance, demonstrating that the values predicted by ANNs have the potential to evaluate the structural performance. Noureldin et al. [15] proposed an ANN and genetic algorithm optimization to predict the different EDPs of low-rise buildings. Liu et al. [16] presented a fragility analysis framework using a data-driven approach, and ANNs were used to conduct the seismic resilience assessment of bridge networks under spatially correlated earthquakes. Wang et al. [13] found ANNs to be more suitable in estimating the probability of failure of nuclear power plant equipment. Li et al. [17] proposed the classical ANN model to develop non-parametric fragility curves to assess the earthquake risks of gravity dams. Yu et al. [18] recently presented an approach coupling the genetic algorithm method and ANNs to predict the seismic responses of arch dams. They found that the proposed method could significantly increase the computational efficiency and replace some dynamic analysis calculations. Dang-Vu et al. [19] provided a surrogate model that aimed to forecast the seismic responses of specific structural parts inside buildings characterized by intrinsic vertical and horizontal irregularities, leading to components with varying seismic vulnerabilities. Ghasemi and Stephens [20] presented a machine learning clustering approach to group buildings of the same type and choose key buildings to study how they respond to and are damaged by earthquakes in the area. Işık et al. [21] used the ANN method to assess the impact of different materials on seismic vulnerability. Jena et al. [22] implemented the Long Short-Term Memory model for Geospatial Information Systems to assess the earthquake vulnerability for the whole of India. Kazemi et al. [23] developed a risk assessment machine learning tool for the purpose of retrofitting and to obtain potential design strategies for RC buildings. Soleimani and Liu [24] used the ANN method for predictive PSDMs using a reputable ML approach. Urlainis and Shohet [25,26] developed the fragility parameters via a failure analysis for each damage state, using a Fault Tree Analysis and approximation of the fragility parameters in accordance with the rate of exceedance. Xu et al. [27] proposed an artificial neural network framework to simultaneously predict the nonlinear seismic responses of all buildings in a cluster subjected to multi-earthquake inputs. More recently, DL algorithms for underground structures have been mainly used for underground rock engineering design [28], underground mine planning and schedulin [29], the forecasting of rock mass and concrete lining deformation [30,31], and microseismic event recognition induced by underground excavatio [32]. There are relatively limited studies on the application of DL algorithms to predict the inelastic dynamic responses of underground structures [14]. Further studies are thus needed to utilize DL to help engineers and decision-makers in taking suitable intervention actions to prevent the possible repercussions of CFHT failure.

To appropriately address the above difficulties, this study intends to use ANNs to develop high-fidelity and high-efficiency prediction models that can replace the extensive numerical simulation in performing seismic fragility assessments of CFHTs. This manuscript is organized as follows. The methodology of the ANNs and quality measures is briefly explained in Section 2. The seismic vulnerability analysis methodology based on maximum likelihood estimation (MLE) and IDA is presented in Section 3. Then, the FE numerical model of CFHTs considering a multimedium-coupled interaction system is presented in Section 4 to estimate the structural responses. Afterward, the ANN models with different hyper-parameters are discussed in Section 5. The analysis results and discussion are presented in Section 6. Finally, the concluding remarks derived from the results are presented in Section 7.

2. A Deep Learning Method to Enhance the Deterministic Evaluation of Seismic Responses

2.1. Artificial Neural Network Methodology

The ANN methodology is utilizes deep learning models to depict the complex information processing mechanisms of the human brain and nervous system based on the knowledge of a network topology after approximating the biological NN structure and the human brain’s stimulus response mechanism [33]. Figure 1 shows a typical neuronal structure in the human body. It can be seen that a single neuron in the human body collects information transmitted by other neurons on dendrites. In the mathematical field, Most conventional ANNs use a multilayer feed-forward backpropagation network, which involves the interconnection of perceptrons with a unidirectional flow of data and computation from inputs to outputs. The objective of applying a multilayer feed-forward backpropagation network in this study is to present a high-precision regression prediction model connecting the seismic intensity (IM) and structural response. As shown in Figure 2, a classical multilayer feed-forward backpropagation network of ANNs is composed of three types of layers: (i) input layers, which receive the input data from the predicted object and are composed of neurons; (ii) hidden layers, which are regarded as a black box with mathematical functions to produce an output specific to the intended results; (iii) output layers, which are the final predicted results of the ANNs. Nodes or neurons make up these layers, representing an activation function for a specific output. Furthermore, a set of weights is used to reflect the connection between each two neuron nodes, providing priority to inputs for the task that the algorithm is seeking to learn. Given input data $x_1,x_2,\dots ,x_n$, the output data ($z_i$) of neuron $i$ in the hidden layer can be written as

$$z_i=f\left({\displaystyle \sum _{j=1}^n{w}_{ij}\cdot x_j+w_i}\right)$$

(1)

where $f\left(\cdot \right)$ indicates the activation function or transfer function; $n$ is the number of neuron nodes; k is each neuron node of the given layers; and $w_{ij}$ represents the synaptic weights. During the deep learning process, the synaptic weights $w_{ij}$ of each neuron can be written as

$$w_{ij}=w_{ij}-u\frac{\partial E}{\partial w_{ij}}, 0 < u < 1$$

(2)

where $E$ indicates the loss function, which the ANN aims to minimize. In this study, Matlab 2021a is utilized to implement the ANNs and therefore benefit from building the ANN model in this software package to predict the specifically required dataset, and the following metamodel parameters can enhance the predictive ability.

Number of layers and neurons in hidden layers: The ANNs map the mathematical relationships between the input and output databases through the connection weight relationships between each layer structure. Evidently, the regression ability of the metamodel largely depends on the complexity of the relationship between the layers of the ANNs. In general, the more hidden layers and neurons there are, the greater the nonlinear fitting ability of the ANNs. At the same time, the calculation cost required for model training calculation will increase significantly, and it will easily lead to the insufficient generalization ability of the ANN model. For the training data, an excellent fitting effect can be obtained from the ANN model. However, the ANN model cannot be well fitted to a database outside the training data, which is called the over-fitting phenomenon. In contrast, too few hidden layers and neurons fails to provide the ANNs with sufficient fitting ability, which will lead to insufficient accuracy in the model training results. Hence, the number of hidden layers and the corresponding number of neurons significantly impact the fitting performance of the ANNs.

Activation functions: The activation function is responsible for transforming the weighted sum input from the neurons of input layers into an output value that is then used as input to the following hidden layers or output layers. The existence of an activation function introduces nonlinear factors into ANNs and improves their ability to simulate nonlinear mapping relations. In the field of DL, the activation functions used mainly include the Sigmoid function, Tanh function, and ReLU function.

Optimization algorithm: The actual purpose of the ANN training process is to find the appropriate inter-layer weight coefficients so that the loss function value of the ANN model training results is as small as possible.

2.2. Prediction of Structural Responses by the Deep Learning Method

The main objective of the applied ANNs is to establish a reliabile mathematical regression model related to the structural responses and earthquake IMs to improve the computational efficiency of the numerical model and decrease the uncertainty of the selected GMs and dynamic interaction systems. Before constructing the ANN predicton models, thousands of finite element (FE) simulations are conducted to develop the primary training database. Meanwhile, it is also necessary to select different earthquake IMs to reflect the various seismology characteristics and reveal the effect of different hyper-parameters. Figure 3 illustrates the operational flow and the main steps are as follows.

(I)

A database reflecting the inelastic dynamic response analysis of the FE simulations is generated using the CFHT–surrounding mass–fluid dynamic interaction system. After the FE simulations, the structural responses’ information time histories can be thus extracted.

(II)

The feature selection of the seismology characteristics of GM records is determined, such as the peak ground acceleration (PGA), peak ground velocity (PGV), and peak ground displacement (PGD), which are significant in the construction of the metamodel. Notably, one disadvantage is that the more feature elements that are selected, the more calculations will be required, because this approach requires a large number of repeated training steps of the metamodel.

(III)

The extracted structural response information and the seismology features of the GM records are merged together to predict the structural responses. During the process, the ratio of the training database size and testing database size is set as 7:3 to ensure the accuracy of the metamodel.

(IV)

The number of neurons and layers in the hidden layers, the activation function, and the optimization algorithm are studied in depth to determine the optimizable hyper-parameters of the ANNs, which can enhance the prediction accuracy.

(V)

Evaluations of the prediction accuracy and efficiency are performed using a range of quality measures, as well as discussions of the limitations and benefits of ANNs.

2.3. Quality Measures

Any artificial database derived from metamodels needs to be examined with respect to a real FE simulation database. To determine the validity of the artificial database, the predictive capacity is evaluated using the mean square error (MSE), root mean square error (RMSE), mean absolute error (MAE), and coefficient of correlation (R²). These metrics reflect the gap that exists between the calculated values from the FE method and the predicted value. The global error of RMSE, MAE, and R can describe the difference between the expected structural response predicted by the ML-PCA and the numerical database from the following relations:

$$\mathit{MSE}=\frac{1}{n_f}{\displaystyle \sum _{i=1}^{n_f}{\left(v_i-{\widehat{v}}_i\right)}^2}$$

(3)

$$\mathit{RMSE}=\sqrt{\frac{1}{n_f}{\displaystyle \sum _{i=1}^{n_f}{\left(v_i-{\widehat{v}}_i\right)}^2}}$$

(4)

$$\mathit{MAE}=\frac{{\displaystyle \sum _{i=1}^{n_f}\max }\left|v_i-{\widehat{v}}_i\right|}{n_f}$$

(5)

$$R^2=1-\frac{{\displaystyle {\sum }_{n_f}{\left(v_i-{\widehat{v}}_i\right)}^2}}{n_f{\sigma }^2}$$

(6)

where $n_f$ is the total number of EDP points in the database; $v_i$ and ${\widehat{v}}_i$ represent the structural response values and predict values, respectively. Note that the MSE can be obtained by squaring the RMSE, and the formula is not given for brevity.

3. Seismic Vulnerability Analysis Coupling Maximum Likelihood Estimation and Incremental Dynamic Analysis

Seismic fragility analysis is a powerful statistical analytic tool to determine a structure’s seismic capability, which reflects the exceedance probability of the different DSs of structures under different earthquake intensities. At the same time, the vulnerability curve can be used to improve the seismic capacity of the structure and optimize the emergency rescue plan through transparent and scientific quantification. The vulnerability curves can be described with the following lognormal cumulative distribution function [9,10,34]:

$$p_f\left(\left.\mathit{EDP}\ge \mathit{edp}\right|\mathit{IM}=\mathit{im}\right)=1-\mathsf{\Phi }\left[\frac{\mathit{ln}(\mathit{edp})-\mathit{ln}({\lambda }_{\mathit{EDP}\left|\mathit{IM}\right.})}{\sqrt{{\beta }_{\mathit{EDP}\left|IM\right.}^2+{\beta }_C^2+{\beta }_{DS}^2}}\right]$$

(7)

where $\mathsf{\Phi }\left(\cdot \right)$ is the cumulative distribution function; ${\beta }_{\mathit{EDP}\left|\mathit{IM}\right.}^{}$ represents the standard deviation; ${\beta }_C^{}$ is the aleatory uncertainty related to the capacity of structures [35]; ${\beta }_{DS}^{}$ is the epistemic uncertainty related to the damage state (DS) of structures; ${\lambda }_{\mathit{EDP}\left|\mathit{IM}\right.}$ represents the mean value at different DSs of structures.

Generally, the classical PSDM related to the DM and earthquake IM in the log-linear function is calculated as

$$ln({\lambda }_{\mathit{EDP}\left|\mathit{IM}\right.})=\mathit{lna}+\mathit{bln}(\mathit{IM})+{\epsilon }_i$$

(8)

where $a$ and $b$ indicate the simple linear regression values, respectively; ${\epsilon }_i$ represents the residual error.

For Equation (8), it is widely accepted that the IM and EDP have a log-normal relationship, and the trend models can use the rational expression of the randomness and the uncertainty of the GM records and underground structures. However, not all IM–EDP relations can be determined by a log-normal linear regression model, which means that Equation (8) lack’s universality in seismic performance assessment. To overcome this limitation, MLE is introduced to conduct the vulnerability analysis of CFHTs. MLE is a method of calculating statistical parameters from a collection of structural responses with the same distribution, and it may be used to quantify the discrepancy between the statistical trend model and the EDP in the ANN-based PSDM. To identify the $p_f$, the highest probability of observing the collapse data that originated from a statistics-based trend model based on the MLE-IDA approach is utilized. The statistical parameters of the log-normal vulnerability function belonging to the MLE-IDA method can be estimated as follows [9,36,37]:

$$\mathit{Likelihood}=\mathsf{\Phi }\left(\frac{\ln (I{M}_i)-In(\theta )}{\beta }\right)$$

(9)

where $I{M}_i$ indicates the $i\mathrm{th}$ earthquake IM value; $\theta$ and $\beta$ are the median and the logarithmic standard deviation. Under the assumption that the $i\mathrm{th}$ earthquake IM value is isolated, the product of the individual likelihoods can reflect the likelihood of all databases:

$$\mathit{Likelihood}=\left({\displaystyle \prod _{i=1}^m\mathsf{\Phi }\left(\frac{\ln \left(I{M}_i\right)-\ln \theta }{\beta }\right)}\right){\left(1-\mathsf{\Phi }\left(\frac{\ln \left(I{M}_{\max }\right)-\ln \theta }{\beta }\right)\right)}^{n-m}$$

(10)

where $n$ and $m$ represent the predicted total number of GM records and the collapse number of structural responses, respectively; $I{M}_{\max }$ indicates the pre-defined maximum IM values. Mathematically, Equation (10) can be equivalent to the following formula:

$$\left\{\widehat{\theta },\widehat{\beta }\right\}=\underset{\theta ,\beta }{{\displaystyle \arg \max }}{\displaystyle \sum _{i=1}^m\left[\ln \left(\mathsf{\Phi }\left(\frac{\ln (I{M}_i)-\ln \theta }{\beta }\right)\right)+\left(n-m\right)\ln \left(1-\mathsf{\Phi }\left(\frac{\ln (I{M}_{\max })-\theta }{\beta }\right)\right)\right]}$$

(11)

With the primary goal of producing vulnerability curves for CFHTs based on the predicted values, a detailed framework description is provided in Figure 4. Simultaneously, each step is listed as follows:

(a)

Select the various seismology features of GM records from the PEER database and obtain the corresponding equivalent load time history from the velocity records and displacement records;

(b)

Establish the FE model of CFHTs considering the fluid–structure interaction effect and artificial boundary conditions;

(c)

Perform the IDA analysis on the CFHTs’ underground motion and excitation records;

(d)

Construct a training metamodel of ANNs for the CFHTs based on the numerical results;

(e)

Determine the control parameters and subsequently develop the probability damage curves according to the MLE-IDA.

4. Numerical Method

4.1. Numerical Model

The numerical model of the typical geometric features of the HT crossing fault fracture zone is built based on the FE code “ABAQUS”, which is frequently utilized to perform a series of FE dynamic analyses of underground structures, as displayed in Figure 5. As can be seen, the scale of the numerical model is $290 \mathrm{m}\times 150 \mathrm{m}\times 127 \mathrm{m}$ (length × width × height), and the angle and width of the fault fracture zone are, respectively, ${60}^{\circ }$ and 10 m. The external width of the HT is 15.30 m, and the thickness of the concrete lining and the exterior height is 1.00 m and 17.15 m, separately. Eight-node reduced integration elements C3D8R in the element library are adopted to model the footwall, hanging wall, fault fracture zone, and structures. The corresponding total number of elements and nodes is 61,440 and 69,986, respectively. Moreover, the eight-node fluid element in the dynamic analysis is used to model the fluid–structure interaction dynamics, which is regarded as a linearly elastic, inviscid, irrotational, and compressible medium. Simultaneously, the maximum C3D8R size is set to 7.0 m, according to the work of Kuhlemeyer and Lysmer [38], which complies with the transmitted wavelengths and may consequently reduce the element mesh size’s influence in the dynamic analysis.

Considering the linear/nonlinear dynamic behavior and plastic deformation of rock, the Drucker–Prager model is introduced to express the mechanism of failure at which the state of the surrounding rock mass reaches an inelastic state [39]. The material properties for the hanging wall and footwall are as follows: mass density = 2700 kg/m³, elastic modulus = 5.0 GPa, Poisson’s ratio = 0.29, friction angle = 35°, cohesion = 0.6 MPa. The material properties for the fault fracture zone are as follows: mass density = 2000 kg/m³, elastic modulus = 0.3 GPa, Poisson’s ratio = 0.33, friction angle = 24.2°, cohesion = 0.1 MPa. Furthermore, the concrete plastic damage model supplemented in the ABAQUS material library is used to describe the damage and cracking of the concrete lining induced by the stiffness reduction laws, as shown in Figure 6. The material properties for the concrete lining are as follows: mass density = 2450 kg/m³, elastic modulus = 28 GPa, Poisson’s ratio = 0.167, compressive yield stress = 16.7 MPa, tensile yield stress = 1.78 MPa.

The applied loads during the static analysis of the FE model are as follows: (1) self-weight load; (2) internal hydrostatic pressure; (3) external hydrostatic pressure. Note that the initial stress state of the rock without the tunnel is balanced. Subsequently, the CFHT considering the stress release of the rock is constructed with 40% of the equivalent nodal force, as shown in Figure 7. During Step 1 and Step 2 of the static analysis process, it can be seen that the bottom truncated boundary condition is fixed in the three directions. The lateral truncated boundary conditions are constant on the normal vector. After the entire tunnel is built, the initial stress state of the system and the truncated boundary node reaction forces are extracted for dynamic analysis. In addition to the loads imposed by the static analysis, hydrodynamic pressure and seismic loads are considered in the dynamic analysis. The dynamic analysis is “restarted” (only the displacement is reset to zero, without affecting the extracted starting stress) after the first static one to simplify the viewing of the deformed shape. Additionally, the viscous-spring artificial boundary (VSAB) is activated to absorb the outward-traveling wave and is used to reflect the recovery ability of the medium outside the FE model (see Figure 7) [41].

4.2. Interaction System and Earthquake Records

During the GM excitation process, it is required to express the relative squeezing, reclosing, sliding, and dislocation between the surrounding rock and concrete lining. To reflect the complex dynamic contact mechanics, the Lagrange multiplier method and penalty function are adopted to establish the rock–structure interface relationship. At the same time, the interface related to the concrete lining and surrounding rock mass is also defined to simulate the normal compressive stress mechanical transferring characteristics and the shear direction contact relationship. Furthermore, a nonlinear shear-displacement-dependent wakening friction law equation can be used to decribe the ultimate shearing stress of the contact surface:

$$\tau ={\tau }_{\mathit{resid}}+({\tau }_{\mathit{peak}}-{\tau }_{\mathit{resid}}){\left(\frac{\alpha }{\alpha +u_s^p}\right)}^{\beta }=c_{\mathit{resid}}+{\sigma }_n\tan {\phi }_{\mathit{resid}}+(c_{\mathit{peak}}+{\sigma }_n\tan {\phi }_{\mathit{peak}}-c_{\mathit{resid}}-{\sigma }_n\tan {\phi }_{\mathit{resid}}){\left(\frac{\alpha }{\alpha +u_s^p}\right)}^{\beta }$$

(12)

where ${\tau }_{\mathit{resid}}$ and ${\tau }_{\mathit{peak}}$ are, respectively, the residual shear strength and peak shear strength; $c_{\mathit{resid}}$ and $c_{\mathit{peak}}$ represent the residual cohesion and peak cohesion, respectively; ${\phi }_{\mathit{resid}}$ and ${\phi }_{\mathit{peak}}$ indicate the residual and peak friction angle, respectively; ${\sigma }_n$ is the normal stress; $\alpha$ and $\beta$ are the control parameters of exponential degradation; and $u_s^p$ indicates the inelastic shear displacement. The parameters for the hanging wall and footwall are as follows: shear stiffness = 4.5 GPa/m, normal stiffness = 4.0 GPa/m, peak cohesion = 1.0 MPa, residual cohesion = 0.4, peak friction angle = 41.5°, residual friction angle = 34.2°, $\alpha$ = 0.005, $\beta$ = 1.2. For the fault fracture zone and footwall/hanging wall, the Hook and Coulomb friction law equation is used to express the normal and tangential contact relationship with the friction coefficient = 0.4. As for the fluid–structure dynamic interaction, the stress–strain relationship of the contact system can be derived as follows:

$$\left\{\begin{array}{c}\kappa \\ {\kappa }_x\\ {\kappa }_y\\ {\kappa }_z\end{array}\right\}=\left[\begin{array}{cccc}{ℂ}_{11}& & & \\ & {ℂ}_{22}& & \\ & & {ℂ}_{33}& \\ & & & {ℂ}_{44}\end{array}\right]\left\{\begin{array}{c}{\epsilon }_v\\ w_x\\ w_y\\ w_z\end{array}\right\}$$

(13)

where $\kappa$ is the hydraulic pressure; ${ℂ}_{11}$ and ${\epsilon }_v$ represent the bulk modulus and volumetric strain of fluid elements; ${\kappa }_i$ indicates rotational stresses ($i=x,y,z$); ${ℂ}_{22}$, ${ℂ}_{33}$, ${ℂ}_{44}$ are constraint parameters; and $w_i$ is rotation ($i=x,y,z$).

To express the randomness and uncertainty of seismic events, a set of 90 GM records is selected from the PEER center of the United States, including the pulse-type near-fault, near-fault without velocity pulse, and far-field, as suggested and used in the work of Sun et al. [42]. The distribution of the magnitude and fault distance of the selected earthquakes is displayed in Figure 8. It can be seen from Figure 8 that the earthquake magnitude and fault distance of the selected GM records cover a wide range. Figure 9 shows the spectral acceleration of the equivalent damping ratio of 5% for single selected earthquake records. It is worth noting that the predominant source of aleatory uncertainty in the dynamic analysis system and seismic performance assessment is the unpredictability of the GM from the seismology features. Sun et al. [40] provide a thorough list of potential earthquake IMs for CFHTs, comparing up to 20 parameters based on correlation, efficiency, proficiency, and sufficiency. The final calibrated earthquake IMs are used as the predictive variables of the ANN model to reduce the aleatory uncertainty, as shown in

Table 1. The IMs considered in this investigation.

Description	Symbol	Definition	Unit	Reference
Peak Ground Acceleration	PGA	$\mathit{PGA}=\max \left|{\ddot{u}}_g(t)\right|$	g	-
Peak Ground Velocity	PGV	$\mathit{PGV}=\max \left|{{\displaystyle \dot{u}}}_g(t)\right|$	m/s	-
Acceleration Root-Mean-Square	Arms	$a_{\mathrm{RMS}}=\sqrt{1/t_{\mathrm{tot}}{\int }_0^{t_{\mathrm{tot}}}{\left[{\ddot{u}}_g(t)\right]}^{2}dt}$	g	[43]
Velocity RMS	Vrms	$v_{\mathrm{RMS}}=\sqrt{1/t_{\mathrm{tot}}{\displaystyle {\int }_0^{t_{\mathrm{tot}}}{\left[{\dot{u}}_g(t)\right]}^{2}dt}}$	m	[44]
Arias Intensity	IA	${{\displaystyle I}}_a=\frac{\pi }{2g}{\displaystyle \underset{0}{\overset{t_{\mathrm{tot}}}{\int }}{\left[{\ddot{u}}_g(t)\right]}^2}dt$	m/s	[45]
Acceleration Spectrum Intensity	ASI	$\mathit{ASI}={\displaystyle {\int }_{0.1}^{0.5}{S}_a(\xi =0.05,T)dT}$	m*s	[46]
Velocity Spectrum Intensity	VSI	$\mathit{VSI}={\displaystyle {\int }_{0.1}^{0.5}{S}_v(\xi =0.05,T)dT}$	m	[46]

. Meanwhile, the 90 GM records are adjusted from 0.1 g to 1.2 g to generate an appropriate structural response.

On the other hand, the dynamic boundary effect induced by seismic waves may be mitigated by utilizing a virtual artificial FE boundary, which involves extracting the calculation region from the semi-infinite elastic medium. In the recent literature, the VSAB is frequently employed as it is capable of imitating the elastic recovery capabilities of the medium and absorbing the dispersed wave during dynamic analysis. More specifically, the VSAB can connect the FE model’s truncated boundary through three continuous springs and dampers in different directions $i=x,y,z$, as demonstrated in the work of Liu et al. [47]. Then, the nodal boundary forces in the direction $i=x,y,z$ of strong GMs are applied at the FE boundary, which can be expressed as

$$f_{li}=K_{li}{d}_{li}^f+C_{li}{v}_{li}^f+A_l{\sigma }_{li}^f$$

(14)

where $K_{li}$ and $C_{li}$ indicate the stiffness coefficients and damping coefficients, respectively; $A_l$ is the average area of all elements around truncated boundary node l; $d_{li}^f$, $v_{li}^f$, and ${\sigma }_{li}^f$ indicate the time history of the displacement, velocity, and nodal stress of the FE node in the free-field seismic wave, respectively.

4.3. Seismic Performance Assessment Index

Prior to the seismic performance assessment, it is essential to define the DSs of CFHTs with different degrees of cumulative damage. In particular, the structural damage index (SDI) can be used to identify the deformation characteristics of tunnels induced by failure mechanisms. A very useful SDI for the identification of the failure mechanism is the relative displacement between the crowns of the arch section and the inverted arch section related to the equivalent cross-sectional diameter of the tunnel, named the drift ratio, based on the work of Andreotti and Lai [48]. The SDI is expressed by three DSs: “no damage”, “slight/moderate damage”, and “extensive damage”, as listed in

Table 2. Definitions of the DS, damage level, and damage measures [48].

Damage State	Damage Level	Damage Measures
		Drift Ratio	Log Std
DS0	No damage	-	-
DS1	Slight/moderate damage	${\left(\mathsf{\delta }/{2\mathrm{R}}_{\mathrm{eq}}\right)}_{{\mathrm{DS}}_1}=0.16\%$	${\mathsf{\beta }}_{\mathrm{DS}1}=0.22$
DS2	Extensive damage	${\left(\mathsf{\delta }/{2\mathrm{R}}_{\mathrm{eq}}\right)}_{{\mathrm{DS}}_2}=0.18\%$	${\mathsf{\beta }}_{\mathrm{DS}2}=0.24$

.

5. Prediction Accuracy of Different Hyper-Parameter Configurations for the Proposed Trained ANN Model

5.1. Predictive Capability Evaluation of the Number of Layers and Neurons of Hidden Layers

Initially, different numbers of layers and neurons in the hidden layers of the trained ANN model are consisted to choose the optimal number of layers and neurons for the PSDM. Under a constant number of 10 neurons and the ReLu activation function, one hidden layer, two hidden layers, and three hidden layers are utilized to reveal the effect of the number of layers on the predicted database. Meanwhile, the PGV of the earthquake IMs, namely the correlation, efficiency, practicality, and proficiency scalar-valued IMs used for deep underground structures, is used for the input layer. Figure 10 displays the MSE, RMSE, MAE, and R² of the ANN based on incorporating a different number of layers. As shown in Figure 10a, the most significant correlation coefficient of 0.81 belongs to the case with two layers, which means that the predicted value strongly correlates with the numerical simulation results. In addition, the MSE, RMSE, and MAE for the case of two layers are lower than or equal to the cases with one layer and three layers, while similar trends can be found in the testing database (see Figure 10b). For the case with three layers, it can be found from Figure 10b that the prediction of the ANN model compared to the testing database appears to be overfitted, which means that the number of hidden layers needs to be reasonably determined.

Figure 11 presents the values of MSE, RMSE, MAE, and R² for the training database and testing database with different numbers of neurons for the ANN models. It is observed from Figure 11a that the coefficient of RMSE exhibits the lowest value for 10 neurons among all the examined ANN models. As the number of neurons increases, the coefficient of RMSE gradually stabilizes at 0.36. Compared to the examined ANN models, the result for 5 neurons, 10 neurons, and 15 neurons of RMSE, MSE, and MAE is the same, while their corresponding coefficients are 0.35, 0.13, and 0.25, respectively. Furthermore, the use of 5 neurons and 10 neurons in the ANN models may be regarded as more efficient than the use of other numbers of neurons according to the coefficients of R², RMSE, MSE, and MAE. From a comprehensive perspective, it is suggested that 10 neurons can be utilized to predict the structural EDP, which results in the highest prediction accuracy.

5.2. Predictive Capability Evaluation of Activation Functions

Activation functions can compensate for the shortcomings of the ANN model in learning and understanding very complex and nonlinear functions of underground structures. As described in Section 2, an appropriate activation function can accelerate the calculation speed and convergence speed of the ANN model, reduce the dependence between model parameters, and effectively alleviate the overfitting phenomenon. The linear function, Sigmoid function, Tanh function, and ReLU function are studied and the corresponding coefficients of the MSE, RMSE, MAE, and R² are illustrated in Figure 12. It can be concluded that the predicted values of the Sigmoid function, Tanh function, and ReLU function in the ANN models are correlated with the training database because they have a more significant correlation coefficient than the linear function. In addition, the smallest MAE, MSE, and RMSE of 0.24, 0.12, and 0.35 are attributed to the ReLU function, which means that this function is the most efficient activation function for the training database. For the testing database, the MSE of the ReLU function has the lowest dispersion with regard to the other activation functions, and the ranking of the ANN models is ReLU > Tanh > Sigmoid > linear function, as displayed in Figure 12b. Compared with the linear function, the correlation coefficients (R² = 0.78) of the Tanh activation function and Sigmoid activation function are obviously more significant than 0.76. In addition, the error of the PSDM model established by traditional linear analysis may be large, which will distort the vulnerability assessment results for the CFHTs.

5.3. Predictive Capability Evaluation of Optimization Algorithm

A reasonable optimization function can significantly improve the predictive ability and accuracy of the neural network algorithm model while reducing the value of the MSE. Figure 13 compares the coefficients of MSE, RMSE, MAE, and R² in the training database and testing database using the Bayesian optimization algorithm, random search optimization algorithm, and grid search optimization algorithm and without an optimization algorithm. As expected, the largest and smallest coefficients occur with the Bayesian optimization algorithm. The corresponding values are 0.11, 0.34, 0.23, and 0.81 for the training database, as shown in Figure 13a. Conversely, the grid search optimization algorithm is not suitable for the construction of ANN models. The main reason is that the predicted value is not closely related to the structural EDP and the error is significantly larger than for the other optimization algorithms. Moreover, the applicability of the grid optimization algorithm is lower than when not using an optimization algorithm, both in the training database and testing database (see Figure 13). Thus, not all optimization algorithms can improve the predictive ability of ANN models, especially for a highly nonlinear database. Based on the present results, in the next section, the hyper-parameters selected for the ANN models are as follows: 10 neuros, two hidden layers, the ReLU activation function, and the Bayesian optimization algorithm. These are applied to predict the structural response.

6. Analysis Results and Discussion

6.1. Investigation of the Effect of Seismology Features on Response Uncertainties

Selecting the seismology features of the input layers is a critical step in building ANN models, mainly because the seismology features are closely related to the structural responses. To reduce the error associated with the ANN model, six different types of seismology characteristics are studied to select the appropriate neurons for the input layers.

Table 3. Comparison of RMSE, MSE, MAE, and R2 for the training and testing datasets under different seismology features.

Case	Seismology Features	Database	RMSE	MSE	MAE	R2
Case 1	PGV	Train	0.34	0.11	0.23	0.81
		Test	0.34	0.11	0.24	0.81
Case 2	PGV and Vrms	Train	0.32	0.10	0.22	0.84
		Test	0.33	0.11	0.24	0.82
Case 3	Peak values (PGA, PGV, PGD)	Train	0.29	0.09	0.21	0.86
		Test	0.30	0.09	0.22	0.86
Case 4	Amplitude	Train	0.24	0.06	0.17	0.91
		Test	0.22	0.05	0.16	0.92
Case 5	Integral	Train	0.47	0.22	0.37	0.66
		Test	0.46	0.21	0.36	0.65
Case 6	Frequency content	Train	0.37	0.13	0.27	0.80
		Test	0.33	0.11	0.24	0.79

reports the values of RMSE, MSE, MAE, and R² of the ANN models for the training and testing datasets under different seismology features. The trends of the RMSE, MSE, MAE, and R² in Table 3 are significantly different among the examined cases. The smallest RMSE, MSE, and MAE of Case 4 for the training database and testing database compared to other examined cases indicate that the discreteness of the predicted database for the amplitude type is much higher than for seismology features such as PGA, PGV, and PGD. It is also found that the values of R2 for Cases 1–3 and Case 5 in Table 3 are smaller than those calculated for Case 4. Moreover, it can be found from Table 3 that there is a lower correlation coefficient regarding the integral and frequency content of seismology features compared with the other examined cases. Simultaneously, the coefficients of RMSE, MSE, and MAE for Case 5 and Case 6 are more significant than those for the other examined cases, which means that the predicted database of the corresponding ANN models is not suitable for the subsequent vulnerability analysis. The results also show that each set of information tends to compensate for its limitations, limiting the application of the seismic performance assessment.

6.2. The Feasibility of the Predicted Database by the Proposed ANN Model

The 1080 numerical results calculated by FE simulation are compared with the predicted database of the ANN model in the logarithmic space, as shown in Figure 14. With regard to the correlation from the simple linear regression, the values predicted by the ANN model are strongly associated with the calculated values of the FE dynamic analysis of the CFHT, as shown in Figure 14. In general, the higher the correlation coefficient between the predicted database and the simulation database, the better the ability of the ANN model to predict the structural response. Meanwhile, the statistical distribution characteristics of the structural response for the ANN model and simulation results are studied, as shown in Figure 15. In this figure, the purple rectangle represents the relative frequency of SDIs, and the orange curve is the final result of the Gaussian function fitting. It can be found that the distribution of the predicted EDP of the CFHT under the ANN model has a slight difference from the target database. The difference in the corresponding mean value and the standard deviation is 0.01, which means that the proposed ANN model can predict the nonlinear behavior of CFHTs during low-strength and high-strength seismic events.

To further demonstrate the application of the ANN model in seismic performance assessment, the PSDM of the traditional seismic vulnerability analysis is utilized to derive the statistical seismic IM–demand measure relation, as shown in Figure 16. It can be seen from Figure 16a that the correlation coefficient of the PSDM from the FE simulation data is 0.86, showing a better correlation between the targeted database and PGV. Mathematically, the higher the correlation coefficient between the PGV and EDP, the better the ability of the PGV to explain the EDP and to decrease the uncertainty induced by seismic events. For the predicted database, the correlation coefficients (R² = 0.90) of the PSDM between the PGV and EDP are significantly larger than that of the PGV and FE simulation EDP. In other words, the ANN model has sufficient potential and reliability to bridge the gap between the PGV and EDPs for structural performance assessment. Moreover, the efficiency (${\beta }_{\mathit{EDP}\left|\mathit{IM}\right.}^{}$) of the PSDM related to the PGV and EDP is displayed in Figure 16b. With the consideration of the ANN model and FE simulation, it can be found that the PSDM developed by the predicted database presents better efficiency, ${\beta }_{\mathit{EDP}\left|\mathit{IM}\right.}^{}=0.34$, than the target database, ${\beta }_{\mathit{EDP}\left|\mathit{IM}\right.}^{}=0.39$. Obviously, the ANN model has the potential to reduce the computational costs and develop appropriate PSDMs for the CFHT.

6.3. Seismic Fragility Analysis Using the Probabilistic ANN Model

A comprehensive fragility curve of CFHTs produced by the probabilistic ANN model is presented according to the MLE-IDA methodology and traditional fragility analysis, as displayed in Figure 17. Under slight/moderate damage, the exceedance probability of the MLE-IDA fragility curve of the CFHT is equal to 98.66% at $\mathrm{PGV}=1.5 \mathrm{m}/\mathrm{s}$, which exceeds the value of the traidiotnal fragility curve by approximately 40.48%. Furthermore, it can be found that only when the PGV exceeds 0.35 m/s does the probability of damage to the MLE-IDA exceed that of the traditional fragility curve. However, the difference between the MLE-IDA and the traditional approach is very small. Specifically, the traditional approach to evaluating the damage probability in the slight/moderate DS of CFHTs is generally poorer than the MLE-IDA method. Similarly, the exceedance probability of the CFHT based on the MLE-IDA method is also significantly larger than that of the traditional method when the seismic intensity exceeds 0.41 m/s. For example, the exceedance probability of the CFHT under the MLE-IDA and traditional method is approximately 31.12% and 84.74%, for $\mathrm{PGV}=1.2 \mathrm{m}/\mathrm{s}$, respectively. As the PGV increases, the size of the difference between the MLE-IDA and traditional method also increases gradually. It is worth emphasizing that the above phenomenon induced by the PSDM of underground structures tends to overestimate the EDP under extremely seismic excitation [49].

7. Final Remarks

The work presented herein explored the application of a deep learning methodology to predict the EDPs and generate fragility curves. For this purpose, the regression algorithms of artificial neural networks (ANNs) were investigated by tests to identify the optimal hyper-parameters to predict the structural response. The main seismology characteristics were also considered in different cases, and multivariate earthquake IM features can be integrated into one seismic IM feature through dimensionality reduction. Then, a three-dimensional (3D) numerical cross-fault hydraulic tunnel (CFHT) considering the nonlinear interaction system was chosen as a typical structure to conduct the dynamic time history analysis. Given the input training database regarding the integrated IM and the structural response of the finite element model, the proposed ANN method helps to produce different prediction models. Subsequently, the statistical parameters of the structural responses in the prediction models were compared to test the ANN methodology, while the ANN-based fragility curves of IDA were constructed via the maximum likelihood estimation (MLE) method. From the attained results, the following conclusions can be obtained.

(1)

Compared with the hyper-parameters in terms of the mean square error, root mean square error, mean absolute error, and correlation based on the training database and testing database, it was found that two hidden layers and ten neurons are most suitable parameters for the ANN model. Meanwhile, it was demonstrated that the Bayesian optimization algorithm has higher predictive power for structural responses than other examined optimization algorithms. Furthermore, the direct use of the ANN model is not recommended to predict the structural response for the CFHT since different hyper-parameters have various abilities to change the predicted capacity of the model.

(2)

The predicted database of the integrated peak-value-type IMs of the input layers is more reliable than the scalar-valued IMs according to the examined methods. For the studied case, it was also confirmed that the amplitude type of seismology characteristics predicted by the ANN model has the highest prediction capacity compared with the other examined seismology features. On the other hand, the integral type of seismology characteristics is not suitable for the input layers of the ANN model to predict the structural response.

(3)

Only a slight difference was observed from the PSDM between the predicted EDPs and the calculated value of the CFHT, indicating that it is a reliable tool for the estimation of the seismic performance of CFHTs. Meanwhile, the ANN model is most suitable for the fragility analysis of the CFHTs because it can improve the computational efficiency and overcome the normally distributed assumption. For the seismic fragility curves, the traditional fragility analysis method may underestimate the probability of damage to CFHTs compared with MLE-IDA.

In regard to this study, some limitations should be noted. For instance, the nonlinear dynamic plastic model was assumed to be the Drucker and Prager model, which may lead to dangerously biased conclusions for the tunnel under dynamic conditions. In addition, only the incoming angles of the vertical SV waves were utilized to describe the input of the earthquake load. Moreover, there is still a great deal of work that needs to be done in order to improve the proposed strategy. Different seismology features of earthquake records and multiple evaluation indicators of the structural response should be employed in order to more efficiently express the probability of damage of structures. The corrosion mechanism of the concrete lining induced by complex engineering environments should be further discussed in future studies. Meanwhile, the use of different metamodels, power or nonlinear regression models, statistical approaches, and probabilistic seismic hazard analysis can be considered to improve the prediction accuracy of the deep learning model.