Nested Copula Model for Overall Seismic Vulnerability Analysis of Multispan Bridges

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.


Introduction
Bridges are always the weak links of tra c 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 tra c 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][4][5]. e seismic vulnerability is de ned as the failure probability of structures under speci c intensity of seismic ground motions [6,7]. Four damage states, that is, slight, moderate, serious, and complete, are usually de ned according to the damage degree of structures [8]. e 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][11][12]. Based on the seismic vulnerability of bridges under di erent 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][14][15][16]. Bridge operation management, maintenance, and reinforcement can be carried out e ciently with life-cycle cost analysis when the resources are limited [17]. e 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.
ese 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]. erefore, the research focus of seismic vulnerability is gradually shifting from a single component to the overall bridge system [20][21][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.
e difference between the upper and lower boundaries gets bigger with the increase of the total number of components in the system [25]. us, the system vulnerability cannot be accurately calculated with the firstorder 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]. erefore, 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. is may differ greatly from the actual situation. erefore, 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. e probabilistic seismic demand analyses of these multispan bridge samples were performed with a series of nonlinear time-history calculations. e correlations between the seismic demands of components were modeled and investigated with multicopula functions. e 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. e 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.

Seismic Vulnerability
e 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 P Single is the vulnerability of a single component; IM is the seismic intensity measure; C is the seismic capacity of components; and D is the seismic demand of components. e seismic capacity and the seismic demand can be approximated with lognormal distributions. erefore, equation (1) can be further expressed as follows [33]: where D is the mean of component seismic demands; C is the mean of component seismic capacities; β D|IM is the logarithmic standard deviation of component seismic demands; and β C is the logarithmic standard deviation of component seismic capacities. Previous studies verified that the D 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 as 2 Shock and Vibration where, (C/a) − b 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. e vulnerability of components can be conveniently represented with the median in equation (5).

System
Vulnerability. e bridge is composed of piers, bearings, and abutments. us, 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 N b + N p + 2; P(X i ) represents the failure probability of a single component; and P(X 1 X 2 . . . X N b +N p +2 ) is the joint failure probability of the N b + N p + 2 components damaged simultaneously in the earthquake.

Nested Copula
Model. e 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. erefore, a nested copula model was constructed with a series of optimized copula functions through hierarchy and nesting techniques in this study.

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. e copula function was used to simulate the correlations between random variables, and the arguments were the corresponding marginal distribution functions of random variables. erefore, 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]. erefore, 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]. e function also satisfies ψ (0) � 1, ψ (∞) � 0. e d-dimensional copula is expressed as follows: where ψ − 1 is the inverse function of the generator function, and I d 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).

Definition of Nested Copula.
e correlations between each two random variables are assumed the same in Archimedean Copula class, which may not correspond with the actual situation. erefore, 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].

Shock and Vibration
In d-dimension, the nested copula model is expressed as follows: e hierarchical structure of equation (12) is shown in Figure 1.
In Figure 1, C(·; ψ 1 ) and C(·; ψ S ) are referred to the inner layer Copula functions in the hierarchical structure; C(·; ψ 0 ) is referred to the outer copula function in the hierarchical structure. rough 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.

Parameter Estimation.
e 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 u i , i ∈ 1, . . . , n { }, the log-likelihood function is defined as follows [39]: rough 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:

Goodness-of-Fit
Test. e 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 F(X i ) is the empirical accumulation distribution calculated based on X 1 , . . . , X n . e test statistics is defined as follows [40]: where C θn (u) 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).   Shock and Vibration where e p value of the goodness-of-fit test can be obtained by substituting equations (20) and (21) . e copula function with the highest p value was then selected for the nested copula model.

Structural
Modeling. e bridge model was built with the OpenSees program. e 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. e girder was rigidly connected to the top of the pier. e bearing was simulated with the BoucWen element. e top of the bearing was rigidly connected to the girder. e 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. e 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]. e longitudinal and vertical stiffnesses of the zero-length element were 1.2 × 10 5 kN/m and 1.036 × 10 5 kN/m, respectively. e expansion joints were simulated with collision elements in which single-slot materials were applied. e collision stiffness was taken as 1/3 of the minimum linear stiffness of the girder. e stiffness of the pile was calculated according to the bridge site condition in which the dynamic amplification effect of the soil layer was considered. e equivalent stiffnesses of the piles in six directions were derived, which are shown in Table 2. e finite-element model of the bridge is shown in Figure 3.

Seismic Ground Motion Record Selection.
e 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. is 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  Shock and Vibration ground motions with typical characteristics of intensity, frequency spectrum, and duration. e 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. e acceleration response spectra of the selected ground motion records are depicted in logarithmic coordinates (Figure 4(a)). e magnitude, source distance, shear wave velocity, and duration of the selected ground motion records are shown in Figure 4(b).

Probabilistic Seismic Demand Analysis.
Bridge samples for the probabilistic seismic demand analysis were generated by Latin hypercube sampling technique. e 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. e 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 timehistory 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 c and β of Newmark integrator were set as 0.5 and 0.25, respectively. e analysis step length was the original length of the ground motion records. e 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).

Shock and Vibration
Based on equation (3), mathematical expressions of the probabilistic seismic demand model of bridge components were established (Table 4). e minimum determination coefficient (R) of the probabilistic seismic demand models is 0.699 (Table 4). e p values are all less than 0.01 (Table 4). erefore, the probabilistic seismic demand models derived by regression analysis based on equation (3) are effective and reasonable.

Hierarchical Structure of the Nested Copula Model.
e 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. erefore, the failure of   Shock and Vibration girders was not considered in this study.
e 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)). e nested copula model in Figure 6(b) was constructed by nesting two layers of two-dimensional copula functions. e 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. e combined components (0#, 1#, 2#, 3#, and 4#) were constructed with corresponding components through inner layer copula functions. e overall bridge system was constructed with the combined components (0#, 1#, 2#, 3#, and 4#) through the outer layer copula function.

Probabilistic Conversion of Seismic Demands.
e definition domain of copula functions is [0, 1]. e results of nonlinear time-history analyses cannot be directly applied in the nested copula model. e 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). e 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.

Parameter Estimation and Goodness-of-Fit Test.
e transformed density and distribution functions of the component seismic demands were applied in the inner layer copula functions of the nested copula model. e parameters of the two-dimensional copula functions were estimated according to equations (13) to (17) using the maximum likelihood estimation. e 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 goodnessof-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    (a) C k (C 0 , C 1 , C 2 , C 3 , C 4 ) C 0 (uz 0 , ut 0 ) C 4 (uz 4 , ut 4 ) C 2 (u 2-1 , u 2-2 ) C 3 (u 3-1 , u 3-2 ) C 1 (u 1-1 , u 1-2 ) uz 0 ut 0 uz 4 ut 4 u 1-1 u 1-2 u 2-1 u 2-2 u 3-1 u 3-2 (b) Figure 6: (a) Hierarchical structure of the bridge system and (b) the nested copula model of the bridge system. Shock and Vibration 9 for C 1 and C 3 and the Frank copula is for C 2 . e 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. e 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).

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. e limit values of the four damage states for the pier were 1.0, 2.0, 14.87, and 24.53, respectively. e 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. e 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. e COV was set as 0.25 for the slight and moderate damage states, and 0.5 for the serious and complete damage states. e logarithmic standard deviation of component seismic capacities β C was derived with the COV according to equation.

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). e vulnerability curves under four damage states were depicted separately for illustration (Figure 8).
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%. e 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. e 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.
e medians of the vulnerability of individual components under different damage states are listed in Table 6.

Vulnerability of Combined Components.
e 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). e 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.
e seismic vulnerability of 3# combined component was the smallest among the five combined components. erefore, 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). e 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):   Table 5 are the optimal functions we have chosen and used in the paper.
e 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). e 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. e 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. e 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. e overall seismic vulnerability of the bridge system was overestimated with the assumption that components were completely uncorrelated. e 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. e 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. e 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. erefore, the seismic performance and security of the continuous rigid frame bridge meet the needs of emergency rescue under earthquakes.

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. e density and distribution functions of the components' seismic demands were derived, and the seismic vulnerabilities of the components were developed. e 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 goodnessof-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. e following conclusions were obtained: (1) e nested copula model developed in this study can be directly used to simulate the correlations between multiple components. e 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. e combined components of bearings and abutments were more vulnerable than the combined components of double thin-wall piers.
at is, the continuous rigid frame bridge is more prone to failure at the places of abutments under earthquakes.
(3) e 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%. e 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
e authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.