# Seismic Fragility Analysis of the Reinforced Concrete Continuous Bridge Piers Based on Machine Learning and Symbolic Regression Fusion Algorithms

Authors: Hanbo Zhu; Changqing Miao
Journal: Shock and Vibration
 (2021)

---

## Abstract

In the fragility analysis, researchers mostly chose and constructed seismic intensity measures (IMs) according to past experience and personal preference, resulting in large dispersion between the sample of engineering demand parameter (EDP) and the regression function with IM as the independent variable. This problem needs to be solved urgently. Firstly, the existing 46 types of ground motion intensity measures were taken as a candidate set, and the composite intensity measures (IMs) based on machine learning methods were selected and constructed. Secondly, the modified Park–Ang damage index was taken as EDP, and the symbolic regression method was used to fit the functional relationship between the composite intensity measures (CIMs) and EDP. Finally, the probabilistic seismic demand analysis (PSDA) and seismic fragility analysis were performed by the cloud-stripe method. Taking the pier of a three-span continuous reinforced concrete hollow slab bridge as an example, a nonlinear finite element model was established for vulnerability analysis. And the composite IM was compared with the linear composite IM constructed by Kiani, Lu Dagang, and Liu Tingting. The functions of them were compared. The analysis results indicated that the standard deviation of the composite IM fragility curve proposed in this paper is 60% to 70% smaller than the other composite indicators which verified the efficiency, practicality, proficiency, and sufficiency of the proposed machine learning and symbolic regression fusion algorithms in constructing composite IMs.


---

## Body

## 1. Introduction

The earthquake contains many pieces of random information due to different conditions, such as source, propagation path, and site. To meet the requirements of effectiveness, benefit, adequacy, and robustness [1] in structural fragility analysis, the seismic intensity measures (IMs) were proposed based on three seismic aspects of amplitude, spectrum, and duration [2].The most basic IMs of seismic amplitude are the maximum peak ground motions, which include the peak ground acceleration (PGA), peak ground velocity (PGV), and peak ground displacement (PGD), denoted asaPG, vPG, and dPG. Nuttli et al. [3] took the third peak value in the seismic acceleration and velocity time history as the ground motion intensity measure (IM), denoted by the sustained maximum acceleration (SMA) and the sustained maximum velocity (SMV), respectively. The first and third peaks of the seismic acceleration-time history after filtering the high-frequency noise were used as the ground motion intensity index, which is called the effective design acceleration (EDA) [4].Based on a simplified single-point elastic or elastoplastic model and considering the basic natural vibration period of the structure, the spectrum IMs were proposed. In order to include the dynamic characteristics of structures with different fundamental frequencies, Housner [5] proposed the velocity spectrum IM SIvζ of the structure period, where ζ is the damping ratio. Mackie [6] extended the velocity spectrum intensity to the acceleration spectrum intensity SIvζ. For large-span and high-rise structures, the damage of the structure during ground motion was the result of the superposition of multiple modes. Shome [7] and Brozovicˇ [8] found that the contribution of higher-order modal effects could not be ignored during the damage process. Cordova [9], Lu [10], Tsantaki [11], Adam [12], and so on used the multimode geometric mean response acceleration as the earthquake IMs, respectively, denoted as S¯aC, S¯aL, S¯aT, and S¯aA. Rathje [13] performed component fusion on the Fourier spectrum of the ground motion time history curves and used the average period Tm as the ground motion IMs.The seismic duration IMs refer to the indexes related to the duration of the ground motion. These IMs have little effect on the peak deformation response of the structure, but they will cause a huge difference in cumulative damage [14–17]. Trifunac [18] took the duration t5to95 of cumulative energy consumption from 5% to 95% as a seismic IM. Foschaar [19] took the duration t5to75 of cumulative energy consumption from 5% to 75% as a seismic IM. Malhotra et al. [20] used the maximum response amplitude equivalent cycles Ncy as a ground motion IM.The abovementioned three types of ground motion IMs based on the amplitude, spectrum, and duration have relatively low correlations, while the correlations of IMs in each type are relatively high. According to this characteristic, many scholars have proposed the multifactor seismic IMs considering the three factors through empirical means. Arias [21] took the integral value Ia of energy and square of ground motion acceleration as the ground motion IM. Ayala et al. [22] introduced spectrum parameters into the modified Arias intensity measure Pd. Nau et al. [23] integrated the square function of the acceleration time history over the total duration to obtain the squared cumulative velocity index vsq. Cosenza et al. [24] normalized squared cumulative acceleration ID by dimensionless processing the squared cumulative acceleration asq with aPG and vPG. Kramer [25] proposed the integral of the absolute value of the acceleration time history to obtain the absolute velocity index vCA. Based on the construction principle similar to the Arias index, Vanmarcke [26] proposed average acceleration, velocity, and displacement IMs, which were denoted as aRMS, vRMS, and dRMS. Park et al. [27] constructed composite IM Ic based on aRMS and t5to95. Riddell [28] used ground motion amplitude parameters aPG, vPG, and dPG and duration parameters t5to75 to combine displacement IM ICd, velocity IM ICv, and acceleration IM ICa. In addition, the earthquake IMs also include the average period Tm and frequency fm of the response spectrum, the maximum yield limit fyK of the response spectrum, and the maximum ductility coefficient μK of the elastic single-degree-of-freedom system. While taking into account the effectiveness and benefit of multifactor IMs, their adequacy is constantly improved. While taking into account the effectiveness and usefulness of multifactor IMs, the adequacy is constantly improved. However, the consideration process of multi-indicator factors contains more subjective factors. For the characteristics of different structures, there is no universal process for the selection, fitting, and construction of a unified IM. The robustness of the combined IMs cannot be guaranteed.To further improve the efficiency, practicality, proficiency, and sufficiency of ground motion IMs [29, 30], scholars have proposed vector-type parameter IMs and linear composite IMs. These multiple linear combination IMs are constructed through experience and have achieved good results in specific cases, but there are also many defects. Kiani et al. [31] used vPG and vCA as earthquake IMs and at the same time considered the influence of the two duration parameters of t5to95 and t5to75. It had a good effect on the fragility analysis of specific structures. However, the multicollinearity of t5to95 and t5to75 was not excluded. Cheng et al. [32] used two-dimensional IM vectors [Samax,e, aRMS] and [Samax, e, vCA] for probabilistic seismic demand analysis (PSDA) and obtained better fitting results than single IM. Wang Xiaoping [33] used the residual sum of squares (RSS) and goodness of fit to prove that the two-dimensional vector IMs parameter Samax,e,ICa was relatively effective. However, Lv Dagang and Wang Xiaoping used the cloud-stripe method to reduce the dimensionality of parameter indicator of the two-dimensional vector IM1,IM2, which believed that the variability was caused by IM2 and eliminated the influence of IM1. The method was only suitable for considering a single factor, which could not consider multiple factors. Using the partial least squares (PLS) method, Liu Tingting [29, 34] proved that the log-linear combination of the six IMs aPG,asq,vPG,vsq,Pd,SIvζ=0.05 has the best effect. However, this method could not rule out the strong nonlinear relationship among the selected IMs.In order to take into account efficiency, practicality, proficiency, and sufficiency of ground motion IMs, based on machine learning and symbolic regression fusion algorithms, composite earthquake IMs and fragility analysis of reinforced concrete bridge piers were proposed in this paper. Firstly, the ground motions were initially selected according to the site conditions, and then the greedy algorithm was used to reselect ground motions based on the designed acceleration target spectrum. Secondly, the finite element model was established, and the random dynamic analysis was carried out. Thirdly, the modified Park–Ang index was used as the pier damage index, and the damage grade was divided. Fourthly, based on the sample space of IMs, a variety of machine learning methods are used to select a specific number of IMs. The complex nonlinear relationship between IMs and engineering demand parameter (EDP) was fitted and optimized by symbolic regression, and then the probabilistic seismic demand model was obtained. Finally, the fragility curves of piers were established using the modified cloud map-strip method, and the seismic fragility analysis was carried out.


## 2. Finite Element Modeling and Seismic Wave Selection

### 2.1. Finite Element Modeling

The Wuyi Avenue Bridge was a three-span reinforced concrete continuous beam bridge. The bridge length was 65.04 m. The bridge deck adopted 3 × 20m prestressed concrete hollow slab beams. Each span of it was composed of 9 prestressed concrete hollow slabs. The hollow slabs were prefabricated with C50 concrete. The wet joints and the bridge deck were, respectively, cast with C50 and C40 concrete on-site. In the substructure, the bridge piers were column piers with a height of 6.3 m and a circular cross-sectional diameter of 1.1 m. The bridge abutments were ribbed abutments with a height of 6.3 m. The foundation was bored piles with a height of about 27.2 m and a circular cross-sectional diameter of 1.2 m. C30 concrete was used for the substructure. Grade 270 prestressing strands, HRB335, and HRB400 steel bars were mainly used for all components of the bridge. GJZ200 × 250 × 42 bearings were adopted at piers 1 and 2, and GJZF4200 × 250 × 44 at abutments.OpenSees software was employed to build the finite element model of this bridge. ConcreteCM and Steel4 were used to simulate the nonlinear material constitutive model of concrete and steel bar, respectively. Elastic_Beam_Column element was used to simulate the hollow slab. The impact element was used to simulate the interaction between the main beam and abutment anticollision block. The Elastomeric Bearing Plasticity element was used to simulate the bearing. Hinge Beam_Column element was used to simulate reinforced concrete piers. Fiber section was used to simulate the plastic hinge part. Rigid body characteristics of main abutment nodes were connected by the Rigid_Link element. The PySimple1 model was used to simulate abutment-soil interaction. Link element was used to simulate pile-soil interaction. PySimple1, TzSimple1, and QzSimple1 were used to simulate horizontal resistance of pile side, vertical friction force of pile body, and vertical resistance of pile bottom, respectively. The others were simulated by Elastic_Beam_Column element without considering nonlinearity. The bridge layout, typical section diagram, and modeling method of the Wuyi Avenue Bridge are shown in Figure1.Figure 1
Analysis model of the reinforced concrete continuous girder bridge; (a) reinforced concrete column; (b) rubber bearing; (c) reinforced concrete pile; (d) elastic section beam.


### 2.2. Selection of Ground Motions

The Wuyi Avenue Bridge was a class B bridge with two-stage fortification, which met the requirements of no damage under minor earthquakes and no collapse under severe earthquakes. Minor earthquakes and severe earthquakes refer to the earthquake intensity with the exceeding probability of approximately 63% and 3%, respectively, in the region within 50 years. It was 6-degree fortification intensity. The soil of the site belonged to medium-hard soil, with a shear wave velocity of 250–500 (m/s) and an overburden layer thickness of about 8 m. Therefore, it belonged to the class II site. The PGA was between 0.05 g and 0.8 g. According to the above conditions, 516 ground motions were selected from 23014 seismic records in the database [34].The design acceleration response spectrum of the highway bridge seismic code was(1)ST=Smax0.6TT0+0.4T≤T0SmaxT0<T≤TgSmaxTgTTg<T≤10,where T is the period; T0 is the maximum period of the linear ascending segment of the reaction spectrum, with a value of 0.1 s; Tg is the characteristic period, with a value of 0.35s [36, 37]; and Smax is the maximum value of the designed acceleration response spectrum, as shown in the following formula:(2)Smax=2.5CiCsCdA.In Formula (2), Ci was the seismic importance coefficient, and the importance coefficient of two-level fortification was 0.5 and 1.7, respectively; the site coefficient Cs was 1.0; the damping ratio of the structure was 0.05, the damping adjustment coefficient Cd was 1.0, and the peak acceleration A of horizontal basic ground motion was 0.05 g.The seismic acceleration records were modulated in the range of 0.2–5 times of amplitudes [37]. The horizontal acceleration records were combined according to the square root sum of squares (SRSS) method. The greedy optimization algorithm [38] was used to select 100 combined ground motion records, and the mean square error was 0.72%. The acceleration response spectrum, average acceleration response spectrum, and target acceleration response spectrum of each seismic wave are shown in Figure 2. The acceleration response spectrum of 200 ground motions without combined amplitude modulation also is shown in Figure 2.Figure 2
Seismic acceleration response spectrum: (a) the amplitude modulated spectrum by the maximum SRSS; (b) unmodulated spectrum.
(a)(b)


### 2.3. Engineering Demand Parameters and Damage Grade Classification of Piers

Li et al. [39] defined the damage grade of piers according to the curvature ductility ratio of the pier section as the damage index, which did not consider the low-cycle fatigue effect and used the damage value calculated by the modified Park–Ang damage model as the damage index, denoted as D [40]. The criterion of damage states proposed by Stong et al. [41] and Guo et al. [42] is shown in Table 1.Table 1
Criterion of damage states for reinforced concrete continuous bridge piers.
Damage states (D)No damage or minor local damageRepairableIrreparableCollapse or complete destructionStong et al. [41]D<0.110.11≤D<0.40.4≤D<0.770.77≤D<1Guo et al. [42]D<0.080.08≤D<0.30.4≤D<0.680.68≤D<1Compared with the damage state criterion prosed by Strong, the damage state criterion considering the three-dimensional extension of the criterion through experiments proposed by Guo et al. could better reflect the real damage state of structures, which was suitable for reinforced concrete piers in this article.


## 2.1. Finite Element Modeling

The Wuyi Avenue Bridge was a three-span reinforced concrete continuous beam bridge. The bridge length was 65.04 m. The bridge deck adopted 3 × 20m prestressed concrete hollow slab beams. Each span of it was composed of 9 prestressed concrete hollow slabs. The hollow slabs were prefabricated with C50 concrete. The wet joints and the bridge deck were, respectively, cast with C50 and C40 concrete on-site. In the substructure, the bridge piers were column piers with a height of 6.3 m and a circular cross-sectional diameter of 1.1 m. The bridge abutments were ribbed abutments with a height of 6.3 m. The foundation was bored piles with a height of about 27.2 m and a circular cross-sectional diameter of 1.2 m. C30 concrete was used for the substructure. Grade 270 prestressing strands, HRB335, and HRB400 steel bars were mainly used for all components of the bridge. GJZ200 × 250 × 42 bearings were adopted at piers 1 and 2, and GJZF4200 × 250 × 44 at abutments.OpenSees software was employed to build the finite element model of this bridge. ConcreteCM and Steel4 were used to simulate the nonlinear material constitutive model of concrete and steel bar, respectively. Elastic_Beam_Column element was used to simulate the hollow slab. The impact element was used to simulate the interaction between the main beam and abutment anticollision block. The Elastomeric Bearing Plasticity element was used to simulate the bearing. Hinge Beam_Column element was used to simulate reinforced concrete piers. Fiber section was used to simulate the plastic hinge part. Rigid body characteristics of main abutment nodes were connected by the Rigid_Link element. The PySimple1 model was used to simulate abutment-soil interaction. Link element was used to simulate pile-soil interaction. PySimple1, TzSimple1, and QzSimple1 were used to simulate horizontal resistance of pile side, vertical friction force of pile body, and vertical resistance of pile bottom, respectively. The others were simulated by Elastic_Beam_Column element without considering nonlinearity. The bridge layout, typical section diagram, and modeling method of the Wuyi Avenue Bridge are shown in Figure1.Figure 1
Analysis model of the reinforced concrete continuous girder bridge; (a) reinforced concrete column; (b) rubber bearing; (c) reinforced concrete pile; (d) elastic section beam.


## 2.2. Selection of Ground Motions

The Wuyi Avenue Bridge was a class B bridge with two-stage fortification, which met the requirements of no damage under minor earthquakes and no collapse under severe earthquakes. Minor earthquakes and severe earthquakes refer to the earthquake intensity with the exceeding probability of approximately 63% and 3%, respectively, in the region within 50 years. It was 6-degree fortification intensity. The soil of the site belonged to medium-hard soil, with a shear wave velocity of 250–500 (m/s) and an overburden layer thickness of about 8 m. Therefore, it belonged to the class II site. The PGA was between 0.05 g and 0.8 g. According to the above conditions, 516 ground motions were selected from 23014 seismic records in the database [34].The design acceleration response spectrum of the highway bridge seismic code was(1)ST=Smax0.6TT0+0.4T≤T0SmaxT0<T≤TgSmaxTgTTg<T≤10,where T is the period; T0 is the maximum period of the linear ascending segment of the reaction spectrum, with a value of 0.1 s; Tg is the characteristic period, with a value of 0.35s [36, 37]; and Smax is the maximum value of the designed acceleration response spectrum, as shown in the following formula:(2)Smax=2.5CiCsCdA.In Formula (2), Ci was the seismic importance coefficient, and the importance coefficient of two-level fortification was 0.5 and 1.7, respectively; the site coefficient Cs was 1.0; the damping ratio of the structure was 0.05, the damping adjustment coefficient Cd was 1.0, and the peak acceleration A of horizontal basic ground motion was 0.05 g.The seismic acceleration records were modulated in the range of 0.2–5 times of amplitudes [37]. The horizontal acceleration records were combined according to the square root sum of squares (SRSS) method. The greedy optimization algorithm [38] was used to select 100 combined ground motion records, and the mean square error was 0.72%. The acceleration response spectrum, average acceleration response spectrum, and target acceleration response spectrum of each seismic wave are shown in Figure 2. The acceleration response spectrum of 200 ground motions without combined amplitude modulation also is shown in Figure 2.Figure 2
Seismic acceleration response spectrum: (a) the amplitude modulated spectrum by the maximum SRSS; (b) unmodulated spectrum.
(a)(b)


## 2.3. Engineering Demand Parameters and Damage Grade Classification of Piers

Li et al. [39] defined the damage grade of piers according to the curvature ductility ratio of the pier section as the damage index, which did not consider the low-cycle fatigue effect and used the damage value calculated by the modified Park–Ang damage model as the damage index, denoted as D [40]. The criterion of damage states proposed by Stong et al. [41] and Guo et al. [42] is shown in Table 1.Table 1
Criterion of damage states for reinforced concrete continuous bridge piers.
Damage states (D)No damage or minor local damageRepairableIrreparableCollapse or complete destructionStong et al. [41]D<0.110.11≤D<0.40.4≤D<0.770.77≤D<1Guo et al. [42]D<0.080.08≤D<0.30.4≤D<0.680.68≤D<1Compared with the damage state criterion prosed by Strong, the damage state criterion considering the three-dimensional extension of the criterion through experiments proposed by Guo et al. could better reflect the real damage state of structures, which was suitable for reinforced concrete piers in this article.


## 3. The Calculation Process of Seismic Fragility

### 3.1. Earthquake IMs

According to the amplitude, spectrum, duration, and combination of the earthquake IMs, 46 different IMs are listed in Table2. In addition to the IMs mentioned above, there were the following IMs: IMs of acceleration, velocity, and displacement at the peak of elastic seismic response spectrum were Samax,e, Svmax,e, and Sdmax,e, respectively. IMs of peaks of elastic pseudo acceleration spectrum and pseudo velocity spectrum were Samax,ep and Svmax,ep, respectively. IMs of acceleration, velocity, and displacement at the peak of the elastoplastic acceleration response spectrum were Samax,p, Svmax,p, and Sdmax,p, respectively. IMs of peaks of elastoplastic pseudo acceleration spectrum and pseudo velocity spectrum were Samax,pp and Svmax,pp, respectively. IMs of period and frequency of reaction spectrum were Tm and fm, respectively. IMs of the maximum values of the plastic yield limit and ductility coefficient of the reaction spectrum of the elastic single-degree-of-freedom system were fyK and μK, respectively.Table 2
Classification of IMs.
ClassificationAmplitudeSpectrumDurationCompositesIMsaPGvPGdPGSMASMVEDA1EDA3Samax,eSvmax,eSdmax,eSamax,epSvmax,epSamax,pSvmax,pSdmax,pSamax,ppSvmax,ppSIvζSIaζS¯a.CS¯a.LS¯a.TS¯a.ATmfyKfmμKt5to95t5to75NcyEcumIaeIaPdasqvsqIDIFvCAaRMSvRMSdRMSIcICdICvICaTotal: 46720316


### 3.2. Selection of Earthquake IMs Based on Machine Learning Method

46 different IMs of 100 seismic records in each of the two horizontal directions were calculated. After normalization, a ridge regression algorithm was used to initially filter out irrelevant IMs, and the least squares regression of the dependent variable and the independent variable was fitted by the regularization coefficientβ(k) as follows [43]:(3)βk=XTX+kI−1XTy,where X is a 100×46 full-rank matrix, the row represents the number of samples, and the type of IMs is tabulated; y is the vector of EDP observations; I is the unit vector; k is the regularization coefficient, and the range of interpolation value is from 0.1 to 1.Ask increased, if β(k) continued to tend to 0, it indicated that this parameter had a low correlation with other parameters and was representative. The parameter was retained in the primary election. The result of parameter selection was not changed by ridge regression analysis. However, the calculation efficiency was improved by ridge regression analysis.According to the correlation between dependent variables, the communication with local agent (CLA) clustering method [44], K-means clustering method [45], and K-medoids clustering method [46] were used for characteristic parameter clustering. Among them, K-means clustering and CLA clustering methods were suitable for the sample clustering analysis with double precision data and small data set size; the K-medoids clustering method was suitable for cluster analysis with complex characteristic parameter types (including strings and numbers) and large sample size. In this case, the number of samples was small, and K-means and CLA clustering methods were adopted to cluster the parameters into 3 groups.For the three groups of parameters, the univariate correlation detection algorithm was used to calculate the correlation between the sample of characteristic parameters and the observed values within the group. The univariate correlation detection algorithms were the Pearson correlation coefficient (PCC) algorithm [47], Spearman’s correlation coefficient (SCC) algorithm [48], Kendall’s correlation coefficient (KCC) algorithm [49], distance correlation coefficient (DCC) algorithm [50], maximal information coefficient (MIC) algorithm [51], and so on. The characteristic comparison of the above correlation algorithms is shown in Table 3. The three IMs with the highest correlation with the dependent variable in each category were selected as the final IMs.Table 3
The characteristic comparison of the correlation algorithms.
Correlation algorithmApplication scopeComputational complexityRobustnessCorrelation criterionPCCLinearLowLow0⟶±1, increasedSCCLinear and monotone nonlinearLowMedium0⟶±1, increasedKCCLinear and monotone nonlinearLowMedium0⟶±1, increasedDCCLinear and nonlinearHighHigh0⟶+∞, increasedMICLinear and nonlinearMediumHigh0⟶+∞, increased


### 3.3. PSDA Based on Symbolic Regression Method

The nonlinear relationship between IMs and EDP was fitted and optimized by symbolic regression [52, 53]. An initial set of functions was defined that contains any combination of basic mathematical operators, trigonometric functions, and exponential logarithmic operators. Mean absolute error (EMA) was used as the evaluation index. The greater the index value, the worse the fitness. In order to give consideration to calculation efficiency and stability, EMA>0.6, complexity less than 30, and a number of parameters less than or equal to 3 were taken as the elimination criteria of the genetic algorithm. Population size and maximum genetic algebraic limits were not set, and EMA <0.01 was only taken as the final genetic target.PSDA and seismic fragility were analyzed using the cloud-stripe method. To make up for the shortage of sample size of incremental dynamic analysis (IDA) and PSDA, the PSDA model based on cloud map was changed to the probabilistic intensity model based on strip approach using the modified cloud-stripe method. It was assumed that the probabilistic intensity model conformed to lognormal distribution [54]. The fragility function (cumulative distribution function) of components under specific combination IMs was(4)Pd≥Dk|IMs=x=ϕInd−InSDβDk|IMs,where d represents the component response (EDP) of a component under random seismic loading; Dk represents the damage state; SD represents the EDP of a component; and βDk|IMs represents the logarithmic standard deviation of seismic demand.The logarithmic functionSD of multiple IMs InIM1,InIM2,...InIMj,...InIMm and established by symbolic regression was(5)InSD=fInIM1,InIM2,…InIMj,…InIMm.Formula (5) could be abbreviated as(6)InSD=fIMs.In formula (6), m represents the number of basic IMs in multiple IMs; was the logarithmic combination of InIM1,InIM2,...InIMj,...InIMm.According to formula (6), random sample di of component response was extracted between IMsk−1 and IMsk+1. The logarithmic standard deviation βDk|IMs of damage state Dk was(7)βDk|IMs=∑Indi≥Dk−fIMs2N−n.On this basis, the minimum and maximum values within the range of response bands of demand parameters in the samples were calculated and denoted asdi,min and di,max, respectively. The confidence factor ωi of random sample di was constructed according to the logarithmic distance of di,min⟶Dk or Dk⟶di,max, and ωi is as follows:(8)ωi=1−InDk/dimaxInDk/di,min,Indi,max/Dk2.The confidence factorωi was taken as the weight to correct formula (7). Then, the following could be obtained:(9)βDk|IMs=∑ωiIndi≥Dk−fIMs2∑ωi1−n/N,where n is the number of parameters of the fitting function and N is the number of samples. By comparing formulas (9) and (7), it could be seen that the closer to Dk, the higher the reliability of the sample value in calculating the fragility function at the state.


## 3.1. Earthquake IMs

According to the amplitude, spectrum, duration, and combination of the earthquake IMs, 46 different IMs are listed in Table2. In addition to the IMs mentioned above, there were the following IMs: IMs of acceleration, velocity, and displacement at the peak of elastic seismic response spectrum were Samax,e, Svmax,e, and Sdmax,e, respectively. IMs of peaks of elastic pseudo acceleration spectrum and pseudo velocity spectrum were Samax,ep and Svmax,ep, respectively. IMs of acceleration, velocity, and displacement at the peak of the elastoplastic acceleration response spectrum were Samax,p, Svmax,p, and Sdmax,p, respectively. IMs of peaks of elastoplastic pseudo acceleration spectrum and pseudo velocity spectrum were Samax,pp and Svmax,pp, respectively. IMs of period and frequency of reaction spectrum were Tm and fm, respectively. IMs of the maximum values of the plastic yield limit and ductility coefficient of the reaction spectrum of the elastic single-degree-of-freedom system were fyK and μK, respectively.Table 2
Classification of IMs.
ClassificationAmplitudeSpectrumDurationCompositesIMsaPGvPGdPGSMASMVEDA1EDA3Samax,eSvmax,eSdmax,eSamax,epSvmax,epSamax,pSvmax,pSdmax,pSamax,ppSvmax,ppSIvζSIaζS¯a.CS¯a.LS¯a.TS¯a.ATmfyKfmμKt5to95t5to75NcyEcumIaeIaPdasqvsqIDIFvCAaRMSvRMSdRMSIcICdICvICaTotal: 46720316


## 3.2. Selection of Earthquake IMs Based on Machine Learning Method

46 different IMs of 100 seismic records in each of the two horizontal directions were calculated. After normalization, a ridge regression algorithm was used to initially filter out irrelevant IMs, and the least squares regression of the dependent variable and the independent variable was fitted by the regularization coefficientβ(k) as follows [43]:(3)βk=XTX+kI−1XTy,where X is a 100×46 full-rank matrix, the row represents the number of samples, and the type of IMs is tabulated; y is the vector of EDP observations; I is the unit vector; k is the regularization coefficient, and the range of interpolation value is from 0.1 to 1.Ask increased, if β(k) continued to tend to 0, it indicated that this parameter had a low correlation with other parameters and was representative. The parameter was retained in the primary election. The result of parameter selection was not changed by ridge regression analysis. However, the calculation efficiency was improved by ridge regression analysis.According to the correlation between dependent variables, the communication with local agent (CLA) clustering method [44], K-means clustering method [45], and K-medoids clustering method [46] were used for characteristic parameter clustering. Among them, K-means clustering and CLA clustering methods were suitable for the sample clustering analysis with double precision data and small data set size; the K-medoids clustering method was suitable for cluster analysis with complex characteristic parameter types (including strings and numbers) and large sample size. In this case, the number of samples was small, and K-means and CLA clustering methods were adopted to cluster the parameters into 3 groups.For the three groups of parameters, the univariate correlation detection algorithm was used to calculate the correlation between the sample of characteristic parameters and the observed values within the group. The univariate correlation detection algorithms were the Pearson correlation coefficient (PCC) algorithm [47], Spearman’s correlation coefficient (SCC) algorithm [48], Kendall’s correlation coefficient (KCC) algorithm [49], distance correlation coefficient (DCC) algorithm [50], maximal information coefficient (MIC) algorithm [51], and so on. The characteristic comparison of the above correlation algorithms is shown in Table 3. The three IMs with the highest correlation with the dependent variable in each category were selected as the final IMs.Table 3
The characteristic comparison of the correlation algorithms.
Correlation algorithmApplication scopeComputational complexityRobustnessCorrelation criterionPCCLinearLowLow0⟶±1, increasedSCCLinear and monotone nonlinearLowMedium0⟶±1, increasedKCCLinear and monotone nonlinearLowMedium0⟶±1, increasedDCCLinear and nonlinearHighHigh0⟶+∞, increasedMICLinear and nonlinearMediumHigh0⟶+∞, increased


## 3.3. PSDA Based on Symbolic Regression Method

The nonlinear relationship between IMs and EDP was fitted and optimized by symbolic regression [52, 53]. An initial set of functions was defined that contains any combination of basic mathematical operators, trigonometric functions, and exponential logarithmic operators. Mean absolute error (EMA) was used as the evaluation index. The greater the index value, the worse the fitness. In order to give consideration to calculation efficiency and stability, EMA>0.6, complexity less than 30, and a number of parameters less than or equal to 3 were taken as the elimination criteria of the genetic algorithm. Population size and maximum genetic algebraic limits were not set, and EMA <0.01 was only taken as the final genetic target.PSDA and seismic fragility were analyzed using the cloud-stripe method. To make up for the shortage of sample size of incremental dynamic analysis (IDA) and PSDA, the PSDA model based on cloud map was changed to the probabilistic intensity model based on strip approach using the modified cloud-stripe method. It was assumed that the probabilistic intensity model conformed to lognormal distribution [54]. The fragility function (cumulative distribution function) of components under specific combination IMs was(4)Pd≥Dk|IMs=x=ϕInd−InSDβDk|IMs,where d represents the component response (EDP) of a component under random seismic loading; Dk represents the damage state; SD represents the EDP of a component; and βDk|IMs represents the logarithmic standard deviation of seismic demand.The logarithmic functionSD of multiple IMs InIM1,InIM2,...InIMj,...InIMm and established by symbolic regression was(5)InSD=fInIM1,InIM2,…InIMj,…InIMm.Formula (5) could be abbreviated as(6)InSD=fIMs.In formula (6), m represents the number of basic IMs in multiple IMs; was the logarithmic combination of InIM1,InIM2,...InIMj,...InIMm.According to formula (6), random sample di of component response was extracted between IMsk−1 and IMsk+1. The logarithmic standard deviation βDk|IMs of damage state Dk was(7)βDk|IMs=∑Indi≥Dk−fIMs2N−n.On this basis, the minimum and maximum values within the range of response bands of demand parameters in the samples were calculated and denoted asdi,min and di,max, respectively. The confidence factor ωi of random sample di was constructed according to the logarithmic distance of di,min⟶Dk or Dk⟶di,max, and ωi is as follows:(8)ωi=1−InDk/dimaxInDk/di,min,Indi,max/Dk2.The confidence factorωi was taken as the weight to correct formula (7). Then, the following could be obtained:(9)βDk|IMs=∑ωiIndi≥Dk−fIMs2∑ωi1−n/N,where n is the number of parameters of the fitting function and N is the number of samples. By comparing formulas (9) and (7), it could be seen that the closer to Dk, the higher the reliability of the sample value in calculating the fragility function at the state.


## 4. Fragility Analysis of the Bridge Piers

### 4.1. Structure and Fragility Analysis of Composite Seismic IMs

Considering the coupling effects of bending, shear, and bond-slip of reinforced concrete columns, a finite element model was established and the top force-displacement curve under monotone loading was simulated, as shown in Figure3.Figure 3
The top force-displacement curve under monotone loading.Then, the parameters of the modified Park–Ang model were calculated. The ultimate displacement ductility ratioμu,mo was 13.76 mm. If the finite element calculation conditions were not allowed, the curves and empirical value of μu,mo could be calculated by referring to [55, 56], but the accuracy of it would be affected.The loading path in the random calculation was extracted. The damage of the piers was evaluated by the modified Park–Ang model. According to the criterion of damage states in Table1, the damage grade of the bridge piers is divided according to the limit states in Table 4.Table 4
Limit states of bridge pier.
Limit statesSlightModerateExtensiveCollapseD0.080.30.681To analyze the correlation between 46 different IMs and damage indexD, irrelevant IMs were preliminarily excluded by the least square regression coefficient curve calculated using a ridge regression algorithm. After grouping cross-validation and normalization, the ridge regression curves of ground motion IMs are shown in Figure 4. k was the regularization coefficient, and its value ranges from 0 to 1. However, in this example, when the value of k ranges from 0 to 0.04, the ridge regression curve fluctuates significantly, which is convenient to distinguish its fluctuation law.Figure 4
Ground motion IMs ridge regression curve.In Figure4, the dotted line did not show a monotonous increase or decrease trend, nor did it converge to 0, indicating that D had a weak correlation with parameters such as aPG, vPG, Ecum, pd, vsq, vCA, vRMS, ICv, Sdmax,e, SIa, SaC, SaL, SaT, Samax,pp, Samax,p, uk, Fm, and Ncy. Parameters with weak correlation with D were removed. Others were clustered into 3 groups by the K-means clustering method and CLA clustering method for the remaining IMs. The K-means clustering method had high randomness in the clustering trial calculation. The clustering result using the K-means clustering method was unstable, which was due to the strong dependence of the method on initial point selection and poor clustering effect on high-dimensional data. However, the CLA clustering method had high accuracy and stability, and the clustering results are shown in Table 5.Table 5
Seismic IMs grouping based on the CLA clustering method.
Groups123IMsdPGSMVasqIaIFaRMSdRMSIcICdSvmax,eSvmax,epSIvζS¯a.ASvmax,pSdmax,pSvmax,ppTmt5to95t5to75IDSamax,eSamax,epICaSMAEDA1EDA3μKComparing the characteristics of the five univariate correlation algorithms in Table3, the correlation ranking of IMs and D in the group was carried out. Considering that the calculation principles of the PCC, SCC, and KCC algorithms were similar and the sorting results were basically the same, the PCC and SCC algorithms were excluded. The sequence numbers obtained by the remaining three algorithms (KCC, DCC, and MIC) were added and reordered according to the size, and the relevance ranking of the comprehensive algorithm was obtained. Taking the first group of IMs as an example, the calculation results of correlation ranking are shown in Table 6.Table 6
Relevance ranking of IMs andD in the first group.
Correlation algorithm IMsKCCDCCMICTotalComprehensive algorithmdPG157141SMV4110217asq71327542Ia61224489IF5211215aRMS12426504dRMS339188Ic8824486ICd2610203Svmax,e1111367512Svmax,ep1514459114SIvζ137295416S¯a.A1415448910Svmax,p99296013Sdmax,p1616458711Svmax,pp1010326615Tm17175110217It could be seen from Table6 that the correlation between dPG and D was the highest in the above algorithms, so dPG was used to represent the final parameter of the first group. Using the same method for each group, IMs with less correlation between dPG, ICa, and SMA were selected which were higher correlation with D.Random samples were extracted. Then, the fragility can be directly estimated from the capacity of each damage state (see Table4) as well as the probabilistic seismic demand model (PSDM) parameters obtained from the regression analysis in Figure 5. These parameters are utilized to generate the fragility curves of the piers using formula (9). The logarithmic function of D and composite IMs for bridge column pier in the transverse direction were obtained as formula (11). Fragility curves in the transverse direction are established as shown in Figure 6.(10)InD=21⋅InICa⋅exp−59.024⋅In2ICa−0.320⋅In8.045e−74.483+IndPG−IndPG.Figure 5
The logarithmic function ofD (EDP) and composite IMs for column pier in the transverse direction.Figure 6
Seismic vulnerability curve of the bridge column pier in the transverse direction.Similarly, the seismic IMs along the bridge weredPG, fyK, and SMA. The PSDM parameters of column pier in the longitudinal direction are obtained from the regression analysis as shown in Figure 7. The logarithmic function of D (EDP) and composite IMs in the longitudinal direction was obtained as formula (12). Fragility curves in the longitudinal direction are established as shown in Figure 8.(11)InD=1.101+0.409⋅IndPG−0.106⋅InfyK⋅InSMA−0.105⋅In2dPG⋅exp0.204⋅IndPG−595.Figure 7
The logarithmic function ofD (EDP) and composite IMs for column pier in the longitudinal direction.Figure 8
Seismic vulnerability curve of the bridge column pier in the longitudinal direction.


### 4.2. Comparison of Fragility Analysis Results of Various IMs

At present, the research on the combination of IMs has not formed a system. In recent years, the log-linear composite IMs have been proposed by Kiani [30], Lv Dagang [32], and Liu Tingting [29], which were selected as the comparison group. The optimal single IMs could also be selected by the comprehensive algorithm of KCC, DCC, and MIC. Based on the random ground motion response, the linear combination coefficients of the four IMs were refitted. The basic IMs and logarithmic fitting functions of D in the direction of the longitudinal and transverse bridge are shown in Table 7. A comparison of EDP fitting degree evaluation indexes of single IM and various multivariate IMs is shown in Table 8.Table 7
The EDP fitting function of multivariate IMs and single IMs.
IMsDirectionBasic IMThe EDP fitting functionSingleTransversedPGInD=6.325+1.907⋅IndPGLongitudinaldPGInD=5.244+1.617⋅IndPGKiani [30]TransversevPG, t5to95, t5to75, vCAInD=9.128+3.107⋅InvPG+0.992⋅Int5to75−0.925⋅Int5to95−1.1⋅InvCALongitudinalvPG, t5to95, t5to75, vCAInD=6.515+2.295⋅InvPG+0.906⋅Int5to75−0.891⋅Int5to95−0.544⋅InvCALv Dagang [32]TransverseSamax,e, vCAInD=2.023⋅InSamax,e−0.184⋅Inv−CA3.022LongitudinalSamax,e, vCAInD=1.209⋅InSamax,e+0.305⋅Inv−CA2.759Liu Tingting [29]TransverseaPG, asq, vPG, vsq, Pd, SIvInD=8.712+1.916⋅InvPG+0.783⋅InSIv+0.275⋅Invsq−0.162⋅InPd+0.142⋅InaPG−0.880⋅InasqLongitudinalaPG, asq, vPG, vsq, Pd, SIvInD=6.622+0.908⋅Invsq+0.458⋅InvPG+0.343⋅InSIv+0.286⋅InaPG−0.129⋅InPd−0.693⋅InasqTable 8
The comparison of fitness evaluation index of single IM and various multivariate IMs.
DirectionIMsGoodness of fitCorrelation coefficientMSEMAEFitting complexityNumber of parametersTransverseMachine learning and symbolic regression fusion algorithms0.9840.9930.1000.234305Single0.7010.8371.8801.07022Kiani0.7110.8431.8161.15655Lv Dagang0.5640.7512.7361.20033Liu Tingting0.7150.8461.7901.07177LongitudinalMachine learning and symbolic regression fusion algorithms0.9970.9980.0120.084305Single0.7960.8920.7550.73122Kiani0.7630.8740.8750.78855Lv Dagang0.6240.7911.3870.89133Liu Tingting0.7930.8910.7650.69777MSE is the mean square error; MAE is the mean absolute error.In Table8, the goodness of fit and correlation coefficient of the combination index based on the machine learning algorithm and the fitting function of the modified Park–Ang damage index of piers both exceed 98%, much higher than the other indexes which were 56%–90%. The mean square errors and absolute errors were lower than 25%, much lower than 73%–188% of other indexes. The formula complexity was also high, but the basic condition of 30 was guaranteed. It could be seen that the machine learning composite IMs were superior to other single and multiple composite IMs in describing the response of piers in structures.The logarithmic standard deviation of fragility function in different damage states was calculated according to formulas (4)∼(9). The logarithmic standard deviation of bridge piers fragility curve in different IMs is shown in Table 9.Table 9
βDk|IMs of column pier vulnerability curve.
DirectionIMsDamage statesDS1DS2DS3DS4TransverseMachine learning and symbolic regression fusion algorithms0.3850.3750.3720.371Single1.2571.2461.2421.240Kiani1.3621.3551.3531.352Lv Dagang1.4551.4311.4231.420Liu Tingting1.2711.2651.2621.261LongitudinalMachine learning and symbolic regression fusion algorithms0.1120.1060.1050.105Single0.8670.8340.8270.825Kiani1.1451.1261.1221.121Lv Dagang1.1411.0661.0521.047Liu Tingting0.9190.8820.8750.873Note. DS1, DS2, DS3, and DS4 represent four types of damage states: slight damage, moderate damage, severe damage, and complete damage.The linear fitting coefficients of Kiani, Lv Dagang, and Liu Tingting composite IMs in the corresponding reference cases were 0.970, 0.642, and 0.903, and the logarithmic standard deviation was 0.271, 0.5, and 0.3, respectively. The standard deviations of Kiani, Lv Dagang, and Liu Tingting composite IMs were about 1.35, 1.44, and 0.88 in this case. It indicated the robustness of the composite IMs was even lower than that of the single IM because of the bias in the selection of the seismic IMs.The standard deviation of fragility curves of the composite IMs based on the fusion algorithm of machine learning and symbolic regression was 60%∼70% smaller than that proposed by Kiani, Lv Dagang, and Liu Tingting, indicating that the machine learning method was more effective, beneficial, adequate, and robust in selecting the index fitting function.


## 4.1. Structure and Fragility Analysis of Composite Seismic IMs

Considering the coupling effects of bending, shear, and bond-slip of reinforced concrete columns, a finite element model was established and the top force-displacement curve under monotone loading was simulated, as shown in Figure3.Figure 3
The top force-displacement curve under monotone loading.Then, the parameters of the modified Park–Ang model were calculated. The ultimate displacement ductility ratioμu,mo was 13.76 mm. If the finite element calculation conditions were not allowed, the curves and empirical value of μu,mo could be calculated by referring to [55, 56], but the accuracy of it would be affected.The loading path in the random calculation was extracted. The damage of the piers was evaluated by the modified Park–Ang model. According to the criterion of damage states in Table1, the damage grade of the bridge piers is divided according to the limit states in Table 4.Table 4
Limit states of bridge pier.
Limit statesSlightModerateExtensiveCollapseD0.080.30.681To analyze the correlation between 46 different IMs and damage indexD, irrelevant IMs were preliminarily excluded by the least square regression coefficient curve calculated using a ridge regression algorithm. After grouping cross-validation and normalization, the ridge regression curves of ground motion IMs are shown in Figure 4. k was the regularization coefficient, and its value ranges from 0 to 1. However, in this example, when the value of k ranges from 0 to 0.04, the ridge regression curve fluctuates significantly, which is convenient to distinguish its fluctuation law.Figure 4
Ground motion IMs ridge regression curve.In Figure4, the dotted line did not show a monotonous increase or decrease trend, nor did it converge to 0, indicating that D had a weak correlation with parameters such as aPG, vPG, Ecum, pd, vsq, vCA, vRMS, ICv, Sdmax,e, SIa, SaC, SaL, SaT, Samax,pp, Samax,p, uk, Fm, and Ncy. Parameters with weak correlation with D were removed. Others were clustered into 3 groups by the K-means clustering method and CLA clustering method for the remaining IMs. The K-means clustering method had high randomness in the clustering trial calculation. The clustering result using the K-means clustering method was unstable, which was due to the strong dependence of the method on initial point selection and poor clustering effect on high-dimensional data. However, the CLA clustering method had high accuracy and stability, and the clustering results are shown in Table 5.Table 5
Seismic IMs grouping based on the CLA clustering method.
Groups123IMsdPGSMVasqIaIFaRMSdRMSIcICdSvmax,eSvmax,epSIvζS¯a.ASvmax,pSdmax,pSvmax,ppTmt5to95t5to75IDSamax,eSamax,epICaSMAEDA1EDA3μKComparing the characteristics of the five univariate correlation algorithms in Table3, the correlation ranking of IMs and D in the group was carried out. Considering that the calculation principles of the PCC, SCC, and KCC algorithms were similar and the sorting results were basically the same, the PCC and SCC algorithms were excluded. The sequence numbers obtained by the remaining three algorithms (KCC, DCC, and MIC) were added and reordered according to the size, and the relevance ranking of the comprehensive algorithm was obtained. Taking the first group of IMs as an example, the calculation results of correlation ranking are shown in Table 6.Table 6
Relevance ranking of IMs andD in the first group.
Correlation algorithm IMsKCCDCCMICTotalComprehensive algorithmdPG157141SMV4110217asq71327542Ia61224489IF5211215aRMS12426504dRMS339188Ic8824486ICd2610203Svmax,e1111367512Svmax,ep1514459114SIvζ137295416S¯a.A1415448910Svmax,p99296013Sdmax,p1616458711Svmax,pp1010326615Tm17175110217It could be seen from Table6 that the correlation between dPG and D was the highest in the above algorithms, so dPG was used to represent the final parameter of the first group. Using the same method for each group, IMs with less correlation between dPG, ICa, and SMA were selected which were higher correlation with D.Random samples were extracted. Then, the fragility can be directly estimated from the capacity of each damage state (see Table4) as well as the probabilistic seismic demand model (PSDM) parameters obtained from the regression analysis in Figure 5. These parameters are utilized to generate the fragility curves of the piers using formula (9). The logarithmic function of D and composite IMs for bridge column pier in the transverse direction were obtained as formula (11). Fragility curves in the transverse direction are established as shown in Figure 6.(10)InD=21⋅InICa⋅exp−59.024⋅In2ICa−0.320⋅In8.045e−74.483+IndPG−IndPG.Figure 5
The logarithmic function ofD (EDP) and composite IMs for column pier in the transverse direction.Figure 6
Seismic vulnerability curve of the bridge column pier in the transverse direction.Similarly, the seismic IMs along the bridge weredPG, fyK, and SMA. The PSDM parameters of column pier in the longitudinal direction are obtained from the regression analysis as shown in Figure 7. The logarithmic function of D (EDP) and composite IMs in the longitudinal direction was obtained as formula (12). Fragility curves in the longitudinal direction are established as shown in Figure 8.(11)InD=1.101+0.409⋅IndPG−0.106⋅InfyK⋅InSMA−0.105⋅In2dPG⋅exp0.204⋅IndPG−595.Figure 7
The logarithmic function ofD (EDP) and composite IMs for column pier in the longitudinal direction.Figure 8
Seismic vulnerability curve of the bridge column pier in the longitudinal direction.


## 4.2. Comparison of Fragility Analysis Results of Various IMs

At present, the research on the combination of IMs has not formed a system. In recent years, the log-linear composite IMs have been proposed by Kiani [30], Lv Dagang [32], and Liu Tingting [29], which were selected as the comparison group. The optimal single IMs could also be selected by the comprehensive algorithm of KCC, DCC, and MIC. Based on the random ground motion response, the linear combination coefficients of the four IMs were refitted. The basic IMs and logarithmic fitting functions of D in the direction of the longitudinal and transverse bridge are shown in Table 7. A comparison of EDP fitting degree evaluation indexes of single IM and various multivariate IMs is shown in Table 8.Table 7
The EDP fitting function of multivariate IMs and single IMs.
IMsDirectionBasic IMThe EDP fitting functionSingleTransversedPGInD=6.325+1.907⋅IndPGLongitudinaldPGInD=5.244+1.617⋅IndPGKiani [30]TransversevPG, t5to95, t5to75, vCAInD=9.128+3.107⋅InvPG+0.992⋅Int5to75−0.925⋅Int5to95−1.1⋅InvCALongitudinalvPG, t5to95, t5to75, vCAInD=6.515+2.295⋅InvPG+0.906⋅Int5to75−0.891⋅Int5to95−0.544⋅InvCALv Dagang [32]TransverseSamax,e, vCAInD=2.023⋅InSamax,e−0.184⋅Inv−CA3.022LongitudinalSamax,e, vCAInD=1.209⋅InSamax,e+0.305⋅Inv−CA2.759Liu Tingting [29]TransverseaPG, asq, vPG, vsq, Pd, SIvInD=8.712+1.916⋅InvPG+0.783⋅InSIv+0.275⋅Invsq−0.162⋅InPd+0.142⋅InaPG−0.880⋅InasqLongitudinalaPG, asq, vPG, vsq, Pd, SIvInD=6.622+0.908⋅Invsq+0.458⋅InvPG+0.343⋅InSIv+0.286⋅InaPG−0.129⋅InPd−0.693⋅InasqTable 8
The comparison of fitness evaluation index of single IM and various multivariate IMs.
DirectionIMsGoodness of fitCorrelation coefficientMSEMAEFitting complexityNumber of parametersTransverseMachine learning and symbolic regression fusion algorithms0.9840.9930.1000.234305Single0.7010.8371.8801.07022Kiani0.7110.8431.8161.15655Lv Dagang0.5640.7512.7361.20033Liu Tingting0.7150.8461.7901.07177LongitudinalMachine learning and symbolic regression fusion algorithms0.9970.9980.0120.084305Single0.7960.8920.7550.73122Kiani0.7630.8740.8750.78855Lv Dagang0.6240.7911.3870.89133Liu Tingting0.7930.8910.7650.69777MSE is the mean square error; MAE is the mean absolute error.In Table8, the goodness of fit and correlation coefficient of the combination index based on the machine learning algorithm and the fitting function of the modified Park–Ang damage index of piers both exceed 98%, much higher than the other indexes which were 56%–90%. The mean square errors and absolute errors were lower than 25%, much lower than 73%–188% of other indexes. The formula complexity was also high, but the basic condition of 30 was guaranteed. It could be seen that the machine learning composite IMs were superior to other single and multiple composite IMs in describing the response of piers in structures.The logarithmic standard deviation of fragility function in different damage states was calculated according to formulas (4)∼(9). The logarithmic standard deviation of bridge piers fragility curve in different IMs is shown in Table 9.Table 9
βDk|IMs of column pier vulnerability curve.
DirectionIMsDamage statesDS1DS2DS3DS4TransverseMachine learning and symbolic regression fusion algorithms0.3850.3750.3720.371Single1.2571.2461.2421.240Kiani1.3621.3551.3531.352Lv Dagang1.4551.4311.4231.420Liu Tingting1.2711.2651.2621.261LongitudinalMachine learning and symbolic regression fusion algorithms0.1120.1060.1050.105Single0.8670.8340.8270.825Kiani1.1451.1261.1221.121Lv Dagang1.1411.0661.0521.047Liu Tingting0.9190.8820.8750.873Note. DS1, DS2, DS3, and DS4 represent four types of damage states: slight damage, moderate damage, severe damage, and complete damage.The linear fitting coefficients of Kiani, Lv Dagang, and Liu Tingting composite IMs in the corresponding reference cases were 0.970, 0.642, and 0.903, and the logarithmic standard deviation was 0.271, 0.5, and 0.3, respectively. The standard deviations of Kiani, Lv Dagang, and Liu Tingting composite IMs were about 1.35, 1.44, and 0.88 in this case. It indicated the robustness of the composite IMs was even lower than that of the single IM because of the bias in the selection of the seismic IMs.The standard deviation of fragility curves of the composite IMs based on the fusion algorithm of machine learning and symbolic regression was 60%∼70% smaller than that proposed by Kiani, Lv Dagang, and Liu Tingting, indicating that the machine learning method was more effective, beneficial, adequate, and robust in selecting the index fitting function.


## 5. Conclusion

Based on machine learning and symbolic regression fusion algorithms, structure and fragility analysis of composite earthquake IMs for reinforced concrete bridge piers was presented in this article. The main conclusions were as follows:(1)
The existing seismic parameters were sorted out and 46 common earthquake IMs were given, which were used to describe seismic wave characteristics.(2)
Based on the ridge regression algorithm, CLA clustering algorithm, comprehensive algorithm including KCC, DCC, and MIC, and symbolic regression method, the grouping and selection of seismic IMs were complied by MATLAB software. The functional relationship betweenD and composite IMs was finally obtained.(3)
A probabilistic intensity model based on the modified cloud-stripe method was proposed. This method revised random samples under various damage states by introducing confidence factors, which made up for the inaccurate calculation of the fragility function of the original PSDA in the case of insufficient samples.(4)
A pier column of a three-span continuous reinforced concrete hollow slab bridge was taken as an example, and the composite IMs were obtained based on machine learning and symbolic regression fusion algorithms. Compared with other composite IMs proposed by Kiani, Lv Dagang, and Liu Tingting, it could be concluded that the goodness of fit and correlation coefficient of the IMs and EDP in this article was over 98%, and the mean square error and absolute error were less than 25%, indicating that the IMs proposed in this article could better meet the requirements of seismic fragility evaluation.(5)
At present, this method cannot calculate the component composite IMs index under multiobjective EDPs. In the future, new fitting and optimization methods need to be used.



---

*Source: 8969389-2021-12-06.xml*