#### Abstract

Piers and bearings influence each other in earthquakes, and the failure of any component will affect the whole function of bridges. Thus, it is critical to consider the correlations between multiple components in the seismic vulnerability analysis of the overall bridge system. In this paper, a method of overall seismic vulnerability analysis with nested copula model was proposed to model the correlations between multiple components. The components were combined according to their types to construct the hierarchical structure of the system. The correlation of components was modeled with copula functions, and the nested copula model of the system was eventually developed. Maximum likelihood estimation and goodness-of-fit test were used to select and optimize the copula functions. The overall seismic vulnerability of the multispan bridge system was established using the nested copula model. A four-span continuous rigid frame bridge was used to illustrate the application of the proposed method. The results show that the nested copula model accurately simulated the correlations between multiple components, which would eliminate the overestimation or underestimation by the first-order boundary method.

#### 1. Introduction

Bridges are always the weak links of traffic lines in earthquakes and more prone to be damaged under earthquakes than other structures. Accurate seismic assessments of bridges are critical to the normal operation of traffic lines [1, 2]. As an important part of performance-based earthquake engineering, seismic vulnerability analysis is widely used in seismic risk assessment of bridges because it can be used to accurately assess the impact of uncertainties in bridge structures and earthquakes from the perspective of probability [3–5]. The seismic vulnerability is defined as the failure probability of structures under specific intensity of seismic ground motions [6, 7]. Four damage states, that is, slight, moderate, serious, and complete, are usually defined according to the damage degree of structures [8]. The failure probability of each damage state is derived by theoretical or empirical methods in the seismic vulnerability analysis [4, 9]. In addition to the application in seismic performance evaluation, seismic vulnerability analysis has also been used to the seismic design optimization of bridges [10–12]. Based on the seismic vulnerability of bridges under different damage states, the design parameters of the materials and physical dimensions of components are modulated to achieve the optimal performance in the earthquakes. Together with material degradation theories, earthquake risk, and cost analysis, the seismic vulnerability analysis can be applied in the life-cycle cost-analysis of bridges with the total probability theory and convolution integral technique [13–16]. Bridge operation management, maintenance, and reinforcement can be carried out efficiently with life-cycle cost analysis when the resources are limited [17].

The accuracy of seismic vulnerability calculation is key to its application. Bridges in practice are usually multispan and consist of multiple components such as piers and bearings. These components interact with each other, and the failure of any component will affect the function of the entire bridge [18]. Compared with the study of a single component, the evaluation of the overall bridge seismic performance is more important for the normal operation of traffic lines [19]. Therefore, the research focus of seismic vulnerability is gradually shifting from a single component to the overall bridge system [20–22].

In conventional studies on the seismic vulnerability of structural systems, the first-order boundary method is most widely used due to its simplicity and ease of understanding [23, 24]. However, the first-order boundary method can only obtain the upper and lower boundaries of the system vulnerability. The difference between the upper and lower boundaries gets bigger with the increase of the total number of components in the system [25]. Thus, the system vulnerability cannot be accurately calculated with the first-order boundary method. In contract, second-order boundary method makes the difference between the upper and lower boundaries smaller [26]. But the computational complexity increased significantly with the number of components. Considering the correlation between failure modes and the impact of component damage on the system vulnerability, Monte-Carlo numerical sampling method is used to establish the vulnerability curve of the bridge system with the joint probability seismic demand model [6, 27]. However, the type of edge and joint probability distribution need to be presupposed in the establishment of the joint probability demand model, and the assumed probability distribution may not be suitable to describe the actual situation [28, 29]. In comparison with the Monte-Carlo numerical sampling method, a conditional edge product method is also applied in the calculation of the joint probability demand distribution, which has been proven simpler and more efficient in previous studies [30]. However, the conditional edge product method is based on the assumption of standard normal cumulative distribution in multidimension. Copula technology is relatively new to construct the joint probability demand model, and the types of the edge and joint probability distribution are not restricted in the technique [31]. Therefore, copula technology provides a new idea for the seismic vulnerability analysis of bridge systems.

We applied nested copula technology to separate the edge distributions of individual components from the correlation between components, thus simplifying the establishment of the joint probability distribution model. However, the dimension of the multivariate copula model dramatically increases with the number of components in the bridge system, which depends on the number of bridge spans. Furthermore, the correlations between any two of all components have to be assumed the same in the system. This may differ greatly from the actual situation. Therefore, the application of copula technology in the overall seismic vulnerability analysis of multispan bridge systems deserves to be further studied.

Based on the theory of nesting and hierarchy, we used the nested copula technique to model a multispan bridge with multicopula functions in this study. One hundred multispan bridge samples were developed by Latin Hypercube sampling method to consider the uncertainty in the bridge. The probabilistic seismic demand analyses of these multispan bridge samples were performed with a series of nonlinear time-history calculations. The correlations between the seismic demands of components were modeled and investigated with multicopula functions. The number and form of hierarchies in the nested copula model were determined depending on the type and number of components. Based on maximum likelihood estimation and goodness-of-fit test, the suitable copula functions were selected and optimized for the nested copula model. The seismic vulnerability curves of the combined components and the seismic vulnerability curve of the overall bridge system were then established with the nested copula model.

#### 2. Basic Theory

##### 2.1. Seismic Vulnerability

###### 2.1.1. Component Vulnerability

The seismic vulnerability of components is defined as the probability that the seismic demand of components under a given seismic intensity reaches or exceeds a certain limit state. Mathematically, the seismic vulnerability of components is expressed as follows [32]:where is the vulnerability of a single component; is the seismic intensity measure; *C* is the seismic capacity of components; and *D* is the seismic demand of components. The seismic capacity and the seismic demand can be approximated with lognormal distributions. Therefore, equation (1) can be further expressed as follows [33]:where is the mean of component seismic demands; is the mean of component seismic capacities; is the logarithmic standard deviation of component seismic demands; and is the logarithmic standard deviation of component seismic capacities. Previous studies verified that the can be approximatively expressed as an exponential function of IM [33]where *a* and *b* are fitting coefficients, which can be derived by probabilistic seismic demand analysis of bridges. Substituting equation (3) into equation (2), the seismic vulnerability of components is expressed as follows:

Based on the definition of normal distribution functions, equation (4) can be further expressed aswhere, represents the median of the seismic intensity, which refers to the value of the ground motion intensity corresponding to the failure probability of 50% at a specific damage state. The vulnerability of components can be conveniently represented with the median in equation (5).

###### 2.1.2. System Vulnerability

The bridge is composed of piers, bearings, and abutments. Thus, the overall seismic vulnerability of the bridge system is expressed as follows:where *P*_{Sys} is the overall seismic vulnerability of the system; ∪ represents the union; *P*_{Bi} is the vulnerability of a single bearing; *N*_{b} is the number of bearings; *P*_{Aj} is the vulnerability of a single abutment; *P*_{Pk} is the vulnerability of a single pier; and *N*_{p} is the number of piers. Based on probabilistic statistical methods for decomposition, the overall vulnerability of the system can be further expressed as follows:where *i* represents an integer from 1 to ; *P*(*X*_{i}) represents the failure probability of a single component; and is the joint failure probability of the components damaged simultaneously in the earthquake.

##### 2.2. Nested Copula Model

The joint failure probability of multiple components in equation (7) cannot be derived directly by multiplying the failure probability of each single component as in conventional methods because of the correlations between the seismic demands of multiple components. In order to model the correlations between components, the copula technique was introduced to develop the joint failure probability of multiple components. However, a single copula function cannot model the correlations between the seismic demands of multicomponents. Therefore, a nested copula model was constructed with a series of optimized copula functions through hierarchy and nesting techniques in this study.

###### 2.2.1. Definition of Copula Functions

Marginal distributions (*F*_{1}, *F*_{2},…, *F*_{d}) of random variables [*X* = (*X*_{1}, *X*_{2},…, *X*_{d})] were integrated to obtain the *d*-dimensional joint distribution function *F* of *x* = (*x*_{1}, *x*_{2},…, *x*_{d}), according to Sklar theorem [34],where *C* is the copula function, and *u*_{1}, *u*_{2},…, *u*_{d} donate the marginal distribution functions, respectively.

The copula function was used to simulate the correlations between random variables, and the arguments were the corresponding marginal distribution functions of random variables. Therefore, the joint distribution function was expressed explicitly with the marginal distributions through equation (8), thereby the computational efficiency is significantly improved. By substituting equation (8) into equation (7), the overall vulnerability of the system is expressed as follows:

With the development of copula theory, various forms of copula functions have been proposed. Elliptical copula and Archimedean copula are the two most widely used classes. In the elliptical copula class, Gaussian copula and *t*-copula are the two most popular types [35, 36]. However, they can only capture the symmetric correlations between variables. Large errors may occur if they were applied in the cases of asymmetric correlations between variables. In contract, Archimedean copula can not only describe the symmetric correlations between variables but can also describe the asymmetric correlations such as upper tail correlation or lower tail correlation between variables [37]. Therefore, the Archimedean copula has a wider range of applications.

A generator function *ψ*, a continuous and strictly monotonous decreasing function, is used in the definition of Archimedean Copula [38]. The function also satisfies *ψ* (0) = 1, *ψ* (∞) = 0. The *d*-dimensional copula is expressed as follows:where is the inverse function of the generator function, and is the *d*-dimensional unit variable. Based on different *ψ*, a series of Archimedean copula functions such as Frank, Gumbel, Clayton, and Joe Copula functions were constructed and used (Table 1).

###### 2.2.2. Definition of Nested Copula

The correlations between each two random variables are assumed the same in Archimedean Copula class, which may not correspond with the actual situation. Therefore, the techniques of nesting and hierarchy were introduced for multivariate correlations in this study. According to equation (10), Archimedean copula functions can be nested with equation (11) [37].

In *d*-dimension, the nested copula model is expressed as follows:

The hierarchical structure of equation (12) is shown in Figure 1.

In Figure 1, and are referred to the inner layer Copula functions in the hierarchical structure; is referred to the outer copula function in the hierarchical structure. Through successive nesting and layering, the correlation between multivariate variables can be accurately simulated by different types of copula functions. Finally, the joint multivariate distribution can be transformed into the joint distribution of copula functions, and the calculation difficulty of multivariate joint distribution is significantly reduced.

###### 2.2.3. Parameter Estimation

The parameters (*λ*) of copula functions in the nested copula model were obtained with maximum likelihood estimation method, a commonly used parameter estimation method. Based on the sample , the log-likelihood function is defined as follows [39]:

Through numerical optimization of equation (13), the maximum likelihood estimated value of the parameter is derived as follows:

If the number of samples *n* tends to infinity, the maximum likelihood estimation value tends to the true value, as shown as follows:

###### 2.2.4. Goodness-of-Fit Test

The suitable copula functions need to be selected after the parameter estimation with the goodness-of-fit test method to construct the hierarchical structure of nested copula model. Based on the sample, the empirical copula in the goodness-of-fit test is defined as follows:where is the empirical accumulation distribution calculated based on . The test statistics is defined as follows [40]:where is the copula function, and its parameter was estimated with Maximum likelihood estimation method. Based on null hypothesis, the sample of the test statistic was generated with equation (21).where

The *p* value of the goodness-of-fit test can be obtained by substituting equations (20) and (21) into . The copula function with the highest *p* value was then selected for the nested copula model.

#### 3. A Case Study

##### 3.1. Overview of the Bridge

A four-span concrete continuous rigid frame bridge (66 m + 2 × 120 m + 66 m) (Figure 2) was used to illustrate the application of nested copula technique in the overall seismic vulnerability assessment of multispan bridges. The main girder is a box girder with variable cross-sections. The minimum height of the girder is 2.86 m at section 2.1 (Figure 2). The maximum height of the girder is 7.26 m at section 3.2 (Figure 2). The concrete strength of the girder is 55 MPa. Double thin-walled piers with the sectional dimension of 7.25 × 1.3 m were used in the bridge. The concrete strength of the pier is 50 MPa, and the steel strength is 335 MPa. The thickness of the platforms is 5 m, and the concrete strength is 40 MPa. The pier foundation is composed of 12 piles with a diameter of 2.2 m. The piles are embedded in moderately weathered argillaceous sandstone. The concrete strength of the piles is 35 MPa. Two *U*-shaped abutments with basin-type rubber bearings are used in the bridge. The expansion joints between the girders and abutments are 0.24-m wide. The designed peak acceleration of the ground motion for the bridge is 0.10 g.

##### 3.2. Structural Modeling

The bridge model was built with the OpenSees program. The main girders were simulated with elastic beam-column elements because plastic failure rarely occurs on the main girder in earthquakes. Nonlinear beam-column elements were used to model the piers and the cross-section of piers was divided into fibers. Concrete02 model was used for the simulation of concrete. Steel02 model was selected for the simulation of steel reinforcement. The girder was rigidly connected to the top of the pier. The bearing was simulated with the BoucWen element. The top of the bearing was rigidly connected to the girder. The yield displacement, the yield force and the horizontal stiffness of the bearing was 0.004 m, 112 kN, 2.635 × 10^{3} kN/m, respectively. The abutments were modeled with zero-length elements, and the impacts of pile foundation and backfill on the abutment stiffness were considered in the simulation of abutments [41]. The longitudinal and vertical stiffnesses of the zero-length element were 1.2 × 10^{5} kN/m and 1.036 × 10^{5} kN/m, respectively. The expansion joints were simulated with collision elements in which single-slot materials were applied. The collision stiffness was taken as 1/3 of the minimum linear stiffness of the girder. The stiffness of the pile was calculated according to the bridge site condition in which the dynamic amplification effect of the soil layer was considered. The equivalent stiffnesses of the piles in six directions were derived, which are shown in Table 2. The finite-element model of the bridge is shown in Figure 3.

##### 3.3. Seismic Ground Motion Record Selection

The ground motion records used in the probabilistic seismic demand analysis of bridges were randomly selected, and the number of the selected records (>30) should meet the statistical requirement of the analysis. In this study, a total of 148 ground motion records were selected from the NGA-West2 strong earthquake database of Pacific Earthquake Engineering Research Center. This number does not only meet the statistical requirement, but also ensures the ease of calculation. Furthermore, the selected ground motion records had a wide distribution range in magnitude, epicentral distance, and duration to ensure a wide coverage of general ground motions with typical characteristics of intensity, frequency spectrum, and duration. The main properties of selected ground motion records were as follows: the soil layer shear wave velocity (*V*_{s30}), 200 m/*s* < *V*_{s30} < 540 m/s; magnitude (*M*_{w}) level, 5 < *M*_{w} < 8; 95% energy duration (*T*_{d95}), 0 *s* < *T*_{d95} < 50 s; source distance (R_{jb}), and 0 km < *R*_{jb} < 250 km. The acceleration response spectra of the selected ground motion records are depicted in logarithmic coordinates (Figure 4(a)). The magnitude, source distance, shear wave velocity, and duration of the selected ground motion records are shown in Figure 4(b).

**(a)**

**(b)**

##### 3.4. Probabilistic Seismic Demand Analysis

Bridge samples for the probabilistic seismic demand analysis were generated by Latin hypercube sampling technique. The probability distributions of the structural parameters in the bridge are summarized in Table 3 [42, 43]. In the sampling process, the 5% to 95% interval of each parameter was divided into 148 groups, and each group had the equal probability. The 0% to 5% interval and 95% to 100% interval of each group were not considered due to the influence of extreme value. A total of 148 bridge samples were generated by randomly combining the 148 groups of each parameter. Nonlinear dynamic time-history analyses were then performed for the 148 bridge samples with the 148 ground motion records.

Newton algorithm was used in the nonlinear dynamic time-history analyses, and the *γ* and *β* of Newmark integrator were set as 0.5 and 0.25, respectively. The analysis step length was the original length of the ground motion records. The structural response outputs of the time-history analyses, that is, the column curvature ductility, bearing deformations, and abutment displacements, were used as the seismic demand parameter of piers, bearings, and abutments, respectively. Probabilistic seismic demand models of bridge components were developed by the regression analysis with the results of the nonlinear dynamic time-history analyses (Figure 5).

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

**(g)**

**(h)**

**(i)**

**(j)**

Based on equation (3), mathematical expressions of the probabilistic seismic demand model of bridge components were established (Table 4).

The minimum determination coefficient (*R*) of the probabilistic seismic demand models is 0.699 (Table 4). The *p* values are all less than 0.01 (Table 4). Therefore, the probabilistic seismic demand models derived by regression analysis based on equation (3) are effective and reasonable.

##### 3.5. Nested Copula Model of the Bridge System

###### 3.5.1. Hierarchical Structure of the Nested Copula Model

The overall system of the bridge was constructed with the hierarchical structure in Figure 6(a). In earthquakes, plastic failure always takes place in the pier instead of the girder, and the falling of girders is usually caused by the large displacement of piers or bearings. Therefore, the failure of girders was not considered in this study. The piers and abutments are combined according to their location and relationship (Figure 6(a)). A nested copula model (Figure 6(b)) was developed to model the bridge system depending on the hierarchical structure (Figure 6(a)).

**(a)**

**(b)**

The nested copula model in Figure 6(b) was constructed by nesting two layers of two-dimensional copula functions. The C_{0}, C_{1}, C_{2}, C_{3}, and C_{4} were two-dimensional copula functions for the inner layer, while *C*_{k} was the copula function for the outer layer. The combined components (0#, 1#, 2#, 3#, and 4#) were constructed with corresponding components through inner layer copula functions. The overall bridge system was constructed with the combined components (0#, 1#, 2#, 3#, and 4#) through the outer layer copula function.

###### 3.5.2. Probabilistic Conversion of Seismic Demands

The definition domain of copula functions is [0, 1]. The results of nonlinear time-history analyses cannot be directly applied in the nested copula model. The probabilistic conversions of the results were performed to meet the definition field of copula functions in this study. First, the effects of ground motion intensity on the seismic demands were removed by subtracting the derived value resulted from the probabilistic seismic demand model with the nonlinear time-history analyses (Table 4). The probability distributions of the processed results were then established with the kernel density estimation. After the probabilistic conversion, the density and distribution functions of the component seismic demands are shown in Figure 7.

**(a)**

**(b)**

###### 3.5.3. Parameter Estimation and Goodness-of-Fit Test

The transformed density and distribution functions of the component seismic demands were applied in the inner layer copula functions of the nested copula model. The parameters of the two-dimensional copula functions were estimated according to equations (13) to (17) using the maximum likelihood estimation. The estimate results of the copula function parameters are listed in Table 5. Goodness-of-fit tests were then performed to select the optimal copula functions for the inner layer of the nested copula model. Based on equations (18) to (24), the *p* value of the goodness-of-fit test for each copula functions was calculated (Table 5).

According to the *p* value in Table 5, the Joe Copula is the optimal copula function for *C*_{0} and *C*_{4}, the Gumbel copula is for *C*_{1} and *C*_{3} and the Frank copula is for *C*_{2}. The density and distribution functions for the seismic demands of combined components (0#, 1#, 2#, 3#, and 4#) were then derived with the optimal copula functions. Substituting the density and distribution functions of combined components into the outer layer copula function (*C*_{k}), the parameter of *C*_{k} was calculated using the maximum likelihood estimation. The Gumbel copula function is verified as the optimal copula function for *C*_{k} with the goodness of fit test for the outer layer copula functions (Table 5).

##### 3.6. Overall Vulnerability Analysis of the Bridge System

###### 3.6.1. Limit States of Damages

Four limit states of damages, slight, moderate, serious, and complete were defined according to the bridge functional damage degree. Specifically, the four limit states of the thin-walled pier were quantified with the curvature ductility values based on the cross-section analysis of the pier. The limit values of the four damage states for the pier were 1.0, 2.0, 14.87, and 24.53, respectively. The damage states of bearings were defined with shear deformations according to the shape of the bearing. For the basin-type bearing, the corresponding shear deformations under the four limit states were taken as 100 mm, 200 mm, 250 mm, and 350 mm, respectively. The four limit states of abutments were defined with the displacements of 50 mm, 100 mm, 200 mm, and 400 mm, respectively. Coefficients of variation (COV) were used to describe the uncertainty of the limit states of piers, bearings, and abutments. The COV was set as 0.25 for the slight and moderate damage states, and 0.5 for the serious and complete damage states. The logarithmic standard deviation of component seismic capacities was derived with the COV according to equation.

###### 3.6.2. Vulnerability of Components

Combining the probabilistic demand models and the corresponding limit states of components, the seismic vulnerability curves of piers, bearings and abutments were derived with equation (4). The vulnerability curves under four damage states were depicted separately for illustration (Figure 8).

**(a)**

**(b)**

**(c)**

**(d)**

In slight damage state, the seismic vulnerabilities of piers, bearings, and abutments were generally the same, and the maximum difference among the failure probabilities of components was 27.5%. The difference between the seismic vulnerability of components became bigger with the increase of damage. In moderate, serious, and complete damage states, the maximum differences among the failure probabilities of components were 30.6%, 45.4%, and 41.7%, respectively. The failure probabilities of the bearings were smaller than the piers and abutments in slight and moderate damage states. In the serious and complete damage states, the vulnerability of piers was smaller than bearings and abutments. In general, the vulnerability of components of the same type had good consistency in all the four damage states. However, the seismic vulnerability of piers changed with the height of piers. The medians of the vulnerability of individual components under different damage states are listed in Table 6.

###### 3.6.3. Vulnerability of Combined Components

The vulnerability curves of the combined components (0#, 1#, 2#, 3#, and 4#) were derived by substituting the vulnerability of the corresponding components into the optimal Copula functions (C_{0}, C_{1}, C_{2}, C_{3} and C_{4}). (Figure 9).

**(a)**

**(b)**

**(c)**

**(d)**

The maximum differences in the seismic vulnerabilities of the five combined components were 13.54% and 14.06% in slight and moderate damage states, respectively. In serious and complete damage states, the seismic vulnerability of 0# combined component was almost the same as that of 4# combined component, and their seismic vulnerabilities was significantly greater than those of 1#, 2#, and 3# combined components. The seismic vulnerability of 3# combined component was the smallest among the five combined components. Therefore, the combined components of bearings and abutments were more vulnerable than the combined components of double thin-wall piers.

###### 3.6.4. Vulnerability of the Overall System

According to the hierarchical structure in Figure 6, the vulnerability of combined components was substituted in the outer layer copula function (*C*_{k}) to construct the overall seismic vulnerability of the bridge system with equation (9). The upper and lower boundaries of the overall vulnerability of the system were derived in order to quantify the influence of the correlations between components. If the seismic demands of components were assumed to be completely correlated, the lower boundary of the overall vulnerability of the system was calculated by equation (26):

If the seismic demands of components were assumed to be completely uncorrelated in the series system, the upper boundary of the overall vulnerability of the system was calculated by equation (27):

The overall seismic vulnerabilities of the bridge system with the Nested Copula model were compared with the aforementioned upper and lower boundaries of the overall vulnerability (Figure 10).

**(a)**

**(b)**

**(c)**

**(d)**

The median of the overall vulnerability, which is defined as the peak acceleration of ground motions (PGA) corresponding to 50% of the failure probability, was used to compare the overall vulnerability by nested copula model with the upper and lower boundaries of the overall vulnerability. Of the overall vulnerability obtained by nested copula model, the medians in the four states of slight, moderate, serious, and complete damages were 0.15 g, 0.25 g, 0.55 g, and 0.85 g, respectively. In contrast, the medians in the four damage states were 0.19 g, 0.31 g, 0.63 g, and 1.00 g, respectively, of the lower boundary of the overall vulnerability. Compared to those derived by the Nested Copula model, the increase of the medians of the overall vulnerability in the four damage states were 26.67%, 24.00%, 14.55%, and 17.65%, respectively, if the components were assumed to be completely correlated. The overall seismic vulnerability of the bridge system was underestimated with the assumption that components were completely correlated. Of the upper boundary of the overall vulnerability, the medians in the four damage states were 0.05 g, 0.09 g, 0.21 g, and 0.32 g, respectively. The decrease of the medians compared to those derived by the nested copula model in the four damage states were 66.67%, 64.00%, 61.82%, and 62.35%, respectively. The overall seismic vulnerability of the bridge system was overestimated with the assumption that components were completely uncorrelated. The comparison of the overall vulnerability obtained by nested copula model with the upper and lower boundaries of the overall vulnerability partly verified the accuracy and effectiveness of our nested copula model.

According to the overall seismic vulnerability of the bridge system obtained by the nested copula model, the failure probabilities of the continuous rigid frame bridge in the four failure states were 43.2%, 36.03%, 19.66%, and 6.00%, respectively, when the PGA was 0.10 g. The failure probabilities of the four damage states did not exceed 50% under the designed seismic intensity (*PGA* = 0.10 g), and the failure probabilities of serious and complete damage states were less than 10%.

When the PGA was 0.20 g, twice the design peak seismic acceleration, the overall failure probabilities of the four damage states were 60.85%, 42.00%, 19.16%, and 9.70%, respectively. Only the failure probability of slight damage was larger than 50%, and the failure probability of complete damage was less than 10%. When the PGA increased to 0.40 g, four times the design peak seismic acceleration, the overall failure probabilities of the slight and moderate damage states were 81.43% and 66.63%, while the probabilities of the serious and complete damages were less than 50%, that is., 40.02% and 26.27%, respectively.

In general, the overall seismic vulnerability of the continuous rigid frame bridge was at a relative low level with the designed seismic intensity. The failure probabilities of serious and complete damage states were not larger than 50% even the ground motion intensity reach the four times of the designed intensity. Therefore, the seismic performance and security of the continuous rigid frame bridge meet the needs of emergency rescue under earthquakes.

#### 4. Conclusion

In order to accurately calculate the seismic vulnerability of multispan bridges, a nested copula model was constructed with multiple Copula functions. First, the probabilistic seismic demand analysis was performed with bridge samples obtained by Latin Hypercube sampling to consider the uncertainties in structures and seismic ground motions. The density and distribution functions of the components’ seismic demands were derived, and the seismic vulnerabilities of the components were developed. The inner and outer layer copula functions in the hierarchical structure of the nested copula model were calculated by maximum likelihood estimation, and the optimal copula functions were selected by the goodness-of-fit test. By substituting the vulnerabilities of components into the nested copula model, the overall vulnerability of the bridge system was developed. A four-span continuous rigid frame bridge was used to illustrate the specific application of the proposed method. The following conclusions were obtained:(1)The nested copula model developed in this study can be directly used to simulate the correlations between multiple components. The accuracies of the overall seismic vulnerability of the bridge system were respectively improved by more than 60% and 10% compared with the assumptions that the components were completely uncorrelated and completely correlated.(2)Generally, the seismic vulnerabilities of components of the same type had good consistencies. The combined components of bearings and abutments were more vulnerable than the combined components of double thin-wall piers. That is, the continuous rigid frame bridge is more prone to failure at the places of abutments under earthquakes.(3)The overall seismic vulnerability of the continuous rigid frame bridge in the four damage states was less than 50% under the design seismic intensity (*PGA* = 0.10 g). Under the two times design peak seismic intensity (*PGA* = 0.20 g), the failure probability of only slight damage was larger than 50%, and the failure probability of complete damage was still less than 10%.

The nested copula model applies the nesting and hierarchy theory to combine multiple components into an overall system. In this study, a case of continuous rigid frame bridge was used to illustrate the application of the nested copula model method. Even so, the proposed vulnerability analysis method of overall bridge system can be applied to other bridge types with multispan for seismic vulnerability assessment.

#### Data Availability

Data, models and code that support the findings of this study are available from the corresponding author upon reasonable request.

#### Conflicts of Interest

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

#### Acknowledgments

This research was supported by the National Natural Science Foundation of China (grant no. 51808376), the China Postdoctoral Science Foundation (grant no. 2019M651076), the Technology Innovation Project of Universities in Shanxi Province (grant no. 20192L0276), Shaanxi Provincial Key R&D Program (grant no. 2020GY-096), and Shandong Provincial Department of Transportation Science and Technology Program (grant no. 2021B59). The supports are gratefully acknowledged.