Seismic fragility curves for structures using non-parametric representations

TL;DR

This paper introduces non-parametric methods, binned Monte Carlo simulation and kernel density estimation, to develop seismic fragility curves without assuming a lognormal distribution, improving accuracy in vulnerability assessments.

Contribution

The paper presents two novel non-parametric approaches for seismic fragility curve estimation, avoiding the traditional lognormal assumption and enhancing accuracy with synthetic and recorded ground motion data.

Findings

Non-parametric methods provide more accurate fragility curves than lognormal assumptions.

Curve accuracy depends on ground motion intensity measure and failure criterion.

Non-parametric approaches perform well with limited recorded ground motion data.

Abstract

Fragility curves are commonly used in civil engineering to assess the vulnerability of structures to earthquakes. The probability of failure associated with a prescribed criterion (e.g. the maximal inter-storey drift of a building exceeding a certain threshold) is represented as a function of the intensity of the earthquake ground motion (e.g. peak ground acceleration or spectral acceleration). The classical approach relies on assuming a lognormal shape of the fragility curves; it is thus parametric. In this paper, we introduce two non-parametric approaches to establish the fragility curves without employing the above assumption, namely binned Monte Carlo simulation and kernel density estimation. As an illustration, we compute the fragility curves for a three-storey steel frame using a large number of synthetic ground motions. The curves obtained with the non-parametric approaches are compared with respective curves based on the lognormal assumption. A similar comparison is presented for a case when a limited number of recorded ground motions is available. It is found that the accuracy of the lognormal curves depends on the ground motion intensity measure, the failure criterion and most importantly, on the employed method for estimating the parameters of the lognormal shape.

Figures

Tables

Table 1. Table 1: Statistics of synthetic ground motion parameters according to  [ 49 ] .

Parameter	Distribution	Support	${\mu }_X$	${\sigma }_X$
$I_a$ (s$\times$g)	Lognormal	(0, +$\infty$)	0.0468	0.164
$D_{5-95}$ (s)	Beta	[5, 45]	17.3	9.31
$t_{mid}$ (s)	Beta	[0.5, 40]	12.4	7.44
${\omega }_{mid}$/2$\pi$ (Hz)	Gamma	(0, +$\infty$)	5.87	3.11
${\omega }^{\prime }$/2$\pi$ (Hz)	Two-sided exponential	[-2, 0.5]	-0.089	0.185
${\zeta }_f$	Beta	[0.02, 1]	0.213	0.143

Table 2. Table 2: Steel frame structure - Parameters of the obtained fragility curves.

		PGA		Sa
${\delta }_o$	Approach	Median	Log-std	Median	Log-std
0.7%	MLE	$0.35g$	0.70	$0.49g$	0.36
	LR	$0.37g$	0.64	$0.44g$	0.13
	KDE	$0.36g$		$0.45g$
1.5%	MLE	$1.10g$	0.56	$1.66g$	0.31
	LR	$0.87g$	0.64	$1.47g$	0.24
	KDE	$1.08g$		$1.53g$
2.5%	MLE	$1.76g$	0.56	$2.82g$	0.37
	LR	$1.55g$	0.64	$3.29g$	0.24
	KDE	$1.82g$		$3.04g$

Table 3. Table 3: Log-standard deviation of median I ​ M 𝐼 𝑀 IM .

${\delta }_o$	Approach	$PGA$	$Sa$
0.7%	bMCS	$0.0003g$	$0.005g$
	KDE	$0.0005g$	$0.005g$
1.5%	bMCS	$0.037g$	$0.054g$
	KDE	$0.037g$	$0.050g$
2.5%	bMCS	$0.114g$	$0.090g$
	KDE	$0.120g$	$0.080g$

Table 4. Table 4: Bridge performance and respective drift-ratio threshold.

Reference	Level	Description	Damage	Drift ratio ${\delta }_o$
[69]	II	Operational	Minor	0.01
[69]	III	Life safety	Moderate	0.03
[70]	II	Operational	Minor	0.005
[70]	III	Life safety	Moderate	0.015

Equations

$$\lambda (DV)=\int \int \int P(DV|DM)\mathrm{d}P(DM|EDP)\mathrm{d}P(EDP|IM)\left|\mathrm{d}\lambda (IM)\right|,$$

$$\text{Frag}(IM;{\delta }_o)=\mathbb{P}[\Delta \ge {\delta }_o|IM],$$

$$\widehat{\text{Frag}}(IM;{\delta }_o)=\Phi \left(\frac{\ln IM-\ln \alpha }{\beta }\right),$$

$$\mathcal{L}\left(\alpha ,\beta ,\left\{I{M}_i,i=1,\dots ,N\right\}\right)=\prod _{i=1}^N{\left[\text{Frag}(I{M}_i;{\delta }_o)\right]}^{y_i}{\left[1-\text{Frag}(I{M}_i;{\delta }_o)\right]}^{1-y_i},$$

$${\left\{{\alpha }^{\ast };{\beta }^{\ast }\right\}}^{\text{T}}=\arg \max \ln \mathcal{L}\left(\alpha ,\beta ,\left\{I{M}_i,i=1,\dots ,N\right\}\right).$$

$$\ln \Delta =A\ln IM+B+\zeta Z,$$

$${\zeta }^2=\sum _{i=1}^N{e}_i^2/\left(N-2\right),$$

$$\widehat{\text{Frag}}(IM;{\delta }_o)$$

$$=\Phi \left(\frac{\ln IM-\left(\ln {\delta }_o-B\right)/A}{\zeta /A}\right).$$

$$\widetilde{{\Delta }_j}(I{M}_o)={\Delta }_j\frac{I{M}_o}{I{M}_j}.$$

$$\widehat{\text{Frag}}(I{M}_o)=\frac{N_f\left(I{M}_o\right)}{N_s\left(I{M}_o\right)},$$

$$\text{Frag}(a;{\delta }_o)=\mathbb{P}\left(\Delta \ge {\delta }_o|IM=a\right)=\underset{{\delta }_o}{\overset{+\infty }{\int }}{f}_{\Delta }(\delta |IM=a)\mathrm{d}\delta .$$

$$f_{\Delta }(\delta |IM=a)=\frac{f_{\Delta ,IM}(\delta ,a)}{f_{IM}(a)},$$

$${\hat{f}}_X\left(x\right)=\frac{1}{Nh}\sum _{i=1}^{N}K\left(\frac{x-x_i}{h}\right),$$

$${\hat{f}}_X\left(x\right)=\frac{1}{Nh}\sum _{i=1}^N\frac{1}{{\left(2\pi \right)}^{1/2}}\exp \left[-\frac{1}{2}{\left(\frac{x-x_i}{h}\right)}^2\right].$$

$${\hat{f}}_{\mathbf{X}}\left(\mathbf{x}\right)=\frac{1}{N{\left|\mathbf{H}\right|}^{1/2}}\sum _{i=1}^{N}K\left({\mathbf{H}}^{-1/2}(\mathbf{x}-{\mathbf{x}}_i)\right),$$

$${\hat{f}}_{\mathbf{X}}\left(\mathbf{x}\right)=\frac{1}{N{\left|\mathbf{H}\right|}^{1/2}}\sum _{i=1}^N\frac{1}{{\left(2\pi \right)}^{d/2}}\exp \left[-\frac{1}{2}{(\mathbf{x}-{\mathbf{x}}_i)}^{\text{T}}{\mathbf{H}}^{-1}(\mathbf{x}-{\mathbf{x}}_i)\right],$$

$$\text{MISE}=\mathbb{E}\left[\underset{{\mathbb{R}}^d}{\int }{\left[{\hat{f}}_{\mathbf{X}}(\mathbf{x};\mathbf{H})-f_{\mathbf{X}}(\mathbf{x})\right]}^2\mathrm{d}\mathbf{x}\right].$$

$${\hat{f}}_{IM}(a)=\frac{1}{{\left(2\pi \right)}^{1/2}N{h}_{IM}}\sum _{i=1}^N\exp \left[-\frac{1}{2}{\left(\frac{a-{IM}_i}{h_{IM}}\right)}^2\right].$$

$${\hat{f}}_{\Delta ,IM}(\delta ,a)=\frac{1}{2\pi N{\left|\mathbf{H}\right|}^{1/2}}\sum _{i=1}^N\exp \left[-\frac{1}{2}{\left(\begin{array}{c}\delta -{\Delta }_i\\ a-I{M}_i\end{array}\right)}^{\text{T}}{\mathbf{H}}^{-1}\left(\begin{array}{c}\delta -{\Delta }_i\\ a-I{M}_i\end{array}\right)\right].$$

$$\widehat{\text{Frag}}(a;{\delta }_o)=\frac{h_{IM}}{{\left(2\pi \left|\mathbf{H}\right|\right)}^{1/2}}\frac{\underset{{\delta }_o}{\overset{+\infty }{\int }}\sum _{i=1}^N\exp \left[-\frac{1}{2}{\left(\begin{array}{c}\delta -{\Delta }_i\\ a-I{M}_i\end{array}\right)}^{\text{T}}{\mathbf{H}}^{-1}\left(\begin{array}{c}\delta -{\Delta }_i\\ a-I{M}_i\end{array}\right)\right]\text{d}\delta }{\sum _{i=1}^N\exp \left[-\frac{1}{2}{\left(\frac{a-{IM}_i}{h_{IM}}\right)}^2\right]}.$$

$$\underset{y_0}{\overset{+\infty }{\int }}{f}_Y(y|X=x)\mathrm{d}y=\underset{y_0}{\overset{+\infty }{\int }}\frac{f_{X,Y}(x,y)}{f_X(x)}\mathrm{d}y=\underset{\ln y_0}{\overset{+\infty }{\int }}\frac{\frac{f_{U,V}(u,v)}{xy}}{\frac{f_U(u)}{x}}y\mathrm{d}v=\underset{\ln y_0}{\overset{+\infty }{\int }}{f}_V(v|U=u)\mathrm{d}v.$$

$$\widehat{\text{Frag}}(a;{\delta }_o)=\underset{{\delta }_o}{\overset{+\infty }{\int }}{\hat{f}}_{\Delta }(\delta |IM=a)\mathrm{d}\delta =\underset{\ln {\delta }_o}{\overset{+\infty }{\int }}{\hat{f}}_V(v|U=\ln a)\mathrm{d}v.$$

$$I_a=\frac{\pi }{2g}\underset{0}{\overset{T}{\int }}{a}^2(t)\mathrm{d}t.$$

$$I(t)=\frac{\pi }{2g}\underset{0}{\overset{t}{\int }}{a}^2(\tau )\mathrm{d}\tau ,$$

$$t_{\alpha }:I(t_{\alpha })=\alpha I_a\alpha \in [0,1].$$

$$a(t)=\frac{q(t,\mathbf{\alpha })}{{\sigma }_h(t)}\underset{0}{\overset{t}{\int }}h\left[t-\tau ,\mathbf{\lambda }\left(\tau \right)\right]\omega (\tau )\mathrm{d}\tau ,$$

$$h\left[t-\tau ,\mathbf{\lambda }(\tau )\right]=0\text{for}t<\tau$$

$$h\left[t-\tau ,\mathbf{\lambda }(\tau )\right]=\frac{{\omega }_f(\tau )}{\sqrt{1-{\zeta }_f^2(\tau )}}\exp \left[-{\zeta }_f(\tau ){\omega }_f(\tau )(t-\tau )\right]\sin \left[{\omega }_f(\tau )\sqrt{1-{\zeta }_f^2(\tau )}(t-\tau )\right]$$

$$\text{for}t\ge \tau .$$

$${\omega }_f(\tau )={\omega }_{mid}+{\omega }^{\prime }(\tau -t_{mid}),$$

$$q(t,\mathbf{\alpha })={\alpha }_1{t}^{{\alpha }_2-1}\exp (-{\alpha }_{3}t),$$

Full text

Keywords: earthquake engineering – fragility curves – lognormal assumption – non-parametric approach – kernel density estimation – epistemic uncertainty

1 Introduction

The severe socio-economic consequences of several recent earthquakes highlight the need for proper seismic risk assessment as a basis for efficient decision making on mitigation actions and disaster planning. To this end, the probabilistic performance-based earthquake engineering (PBEE) framework has been developed, which allows explicit evaluation of performance measures that serve as decision variables (DV) (e.g. monetary losses, casualties, downtime) accounting for the prevailing uncertainties (e.g. ground motion characteristics, structural properties, damage occurrence). The key steps in the PBEE framework comprise the identification of seismic hazard, the evaluation of structural response, damage analysis and eventually, consequence evaluation. In particular, the mean annual frequency of exceedance of a DV is evaluated as [1, 2, 3]:

[TABLE]

in which $P(x|y)$ is the conditional probability of $x$ given $y$, $DM$ is a damage measure typically defined according to repair costs (e.g. light, moderate or severe damage), $EDP$ is an engineering demand parameter obtained from structural analysis (e.g. force, displacement, drift ratio), $IM$ is an intensity measure characterizing the ground motion severity (e.g. peak ground acceleration, spectral acceleration) and $\lambda(IM)$ is the annual frequency of exceedance of the $IM$. Determination of the probabilistic model $P(EDP|IM)$ constitutes a major challenge in the PBEE framework since the earthquake excitation contributes the most significant part to the uncertainty in the $DV$. The present paper is concerned with this step of the analysis.

The conditional probability $P(EDP\geq\overline{edp}|IM)$, where $\overline{edp}$ denotes an acceptable demand threshold, is commonly represented graphically in the shape of the so-called demand fragility curves [4]. Thus, a demand fragility curve represents the probability that an engineering demand parameter exceeds a prescribed threshold as a function of an intensity measure of the earthquake motion. For the sake of simplicity, demand fragility curves are simply denoted fragility curves hereafter, which is also typical in the literature [5, 6]. We note however that the term fragility may also be used for $P(DM\geq\overline{dm}|IM)$ or $P(DM\geq\overline{dm}|EDP)$, i.e. the conditional probability of the damage measure exceeding a threshold $\overline{dm}$ given the ground motion intensity [7] or the engineering demand parameter [2, 3], respectively.

Originally introduced in the early 1980’s for nuclear safety evaluation [8], fragility curves are nowadays widely used for multiple purposes, e.g. loss estimation [9], assessment of collapse risk [10], design checking [11], evaluation of the effectiveness of retrofit measures [12], etc. Several novel methodological contributions to fragility analysis have been made in recent years, including the development of multi-variate fragility functions [13], the incorporation of Bayesian updating [14] and the consideration of time-dependent fragility [15]. However, the traditional fragility curves remain a popular tool in seismic risk assessment and recent literature is rich with relevant applications on various type of structures, such as irregular buildings [6], underground tunnels [16], pile-supported wharfs [17], wind turbines [18], nuclear power plant equipments [19], masonry buildings [20]. The estimation of such curves is the focus of the present paper.

Fragility curves are typically classified into four categories according to the data sources, namely analytical, empirical, judgment-based or hybrid fragility curves [21]. Analytical fragility curves are derived from data obtained by analyses of structural models. Empirical fragility curves are based on the observation of earthquake-induced damage reported in post-earthquake surveys. Judgment-based curves are estimated by expert panels specialized in the field of earthquake engineering. Hybrid curves are obtained by combining data from different sources. Each of the aforementioned categories has its own advantages and limitations. In this paper, analytical fragility curves based on data collected from numerical structural analyses are of interest.

The typical approach to compute analytical fragility curves presumes that the curves have the shape of a lognormal cumulative distribution function [22, 23]. This approach is therefore considered parametric. The parameters of the lognormal distribution are determined either by maximum likelihood estimation [22, 24, 13] or by fitting a linear probabilistic seismic demand model in the log-scale [5, 25, 26, 27]. The assumption of lognormal fragility curves is almost unanimous in the literature due to the computational convenience as well as due to the ease of combining such curves with other elements of the seismic probabilistic risk assessment framework. However, the validity of such assumption remains questionable (see also [28]).

In this paper, we present two non-parametric approaches for establishing the fragility curves, namely binned Monte Carlo simulation (bMCS) and kernel density estimation (KDE). The main advantage of bMCS over existing techniques also based on Monte Carlo simulation is that it avoids the bias induced by scaling ground motions to predefined intensity levels. In the KDE approach, we introduce a statistical methodology for fragility estimation, which opens new paths for estimating multi-dimensional fragility functions as well. The proposed methods are subsequently used to investigate the validity of the lognormal assumption in a case study, where we develop fragility curves for different thresholds of the maximum drift ratio of a three-story steel frame subject to synthetic ground motions. The comparison between KDE-based and lognormal fragility curves is also shown for a concrete bridge column subject to recorded motions using results from an earlier study by the authors [29]. The proposed methodology can be applied in a straightforward manner to other types of structures or classes of structures or using different failure criteria.

The paper is organized as follows: in Section 2, the different approaches for establishing the fragility curves, namely the lognormal and the proposed bMCS and KDE approaches, are presented. In Section 3, the method recently developed by Rezaeian and Der Kiureghian [30] for generating synthetic ground motions, which is employed in the following numerical investigations, is briefly recalled. The case studies are presented in Sections 4 and 5 and the results are discussed in Section 6. The paper concludes with a summary of the main findings and perspectives on future research.

2 Computation of fragility curves

Fragility curves represent the probability of failure of a system, associated with a specified criterion, for a given intensity measure ($IM$) of the earthquake motion. Failure herein represents the exceedence of a prescribed demand limit. A commonly used demand parameter in earthquake engineering is the maximal drift ratio $\Delta$, i.e. the maximal relative horizontal displacement normalized by the corresponding height [6]. Thus, the fragility function is cast as follows:

[TABLE]

in which $\text{Frag}(IM;\,\delta_{o})$ denotes the fragility at the given $IM$ for a threshold $\delta_{o}$ of $\Delta$. In order to establish the fragility curves, a number $N$ of transient finite element analyses of the structure under consideration are used to provide paired values $\left\{\left(IM_{i},\Delta_{i}\right),i=1,\,\dots\,,N\right\}$.

2.1 Lognormal approach

The classical approach for establishing fragility curves consists in assuming a lognormal shape for the curves described in Eq. (2). Two techniques are typically used to estimate the parameters of the lognormal fragility curves, namely maximum likelihood estimation and linear regression. These are presented below.

2.1.1 Maximum likelihood estimation

One assumes that the fragility curves can be written in the following general form:

[TABLE]

where $\Phi(\cdot)$ denotes the standard Gaussian cumulative distribution function (CDF), $\alpha$ is the “median” and $\beta$ is the “log-standard deviation” of the lognormal curve. Shinozuka et al. [22] proposed the use of maximum likelihood estimation to determine these parameters as follows: One denotes by $\omega$ the event that the demand threshold $\delta_{o}$ is reached or exceeded and assumes that $Y(\omega)$ is a random variable with a Bernoulli distribution. In particular, $Y$ takes the value 1 with probability $\text{Frag}(\cdot;\,\delta_{o})$ and the value 0 with probability $1-\text{Frag}(\cdot;\,\delta_{o})$. Considering a set of $i=1,\,\dots\,,N$ ground motions, the likelihood function reads:

[TABLE]

where $IM_{i}$ is the intensity measure of the $i^{th}$ seismic motion and $y_{i}$ represents a realization of the Bernoulli random variable $Y$. The latter takes the value 1 or 0 depending on whether the structure under the $i^{th}$ ground motion sustains the demand threshold $\delta_{o}$ or not. The parameters $(\alpha,\,\beta)$ are obtained by maximizing the likelihood function. In practice, a straightforward optimization algorithm is applied on the log-likelihood function:

[TABLE]

2.1.2 Linear regression

One first assumes a probabilistic seismic demand model, which relates a structural response quantity of interest (herein drift ratio) to an intensity measure of the earthquake motion. Specifically, the demand $\Delta$ is assumed to follow a lognormal distribution of which the log-mean value is a linear function of $\ln IM$, leading to:

[TABLE]

where $Z\sim{\mathcal{N}}(0,1)$ is a standard normal variable. Parameters $A$ and $B$ are determined by means of ordinary least squares estimation in a log-log scale. Parameter $\zeta$ is obtained by:

[TABLE]

where $e_{i}$ the residual between the actual value $\ln\Delta$ and the value predicted by the linear model: $e_{i}=\ln\Delta_{i}-A\ln\left(IM_{i}\right)-B$. Then, Eq. (2) rewrites:

[TABLE]

A comparison to Eq. (3) shows that the median and log-standard deviation of the lognormal fragility curve in Eq. (8) are $\alpha=\exp\left[\left(\ln\delta_{o}-B\right)/A\right]$ and $\beta=\zeta/A$, respectively. This approach to fragility estimation is widely employed in the literature, see e.g. [23, 31, 32, 33] among others.

The two methods described in this section are parametric because they impose the shape of the fragility curves (Eq. (3) and Eq. (8)), which is that of a lognormal CDF when considered as a function of $IM$. We note that by using the linear-regression approach, one accepts two additional assumptions, namely the linear function for the log-mean value of $\Delta$ and the constant dispersion (or homoscedasticity) of the residuals independently of the $IM$ level. Effects of these assumptions have been investigated by Karamlou and Bocchini [28]. In the sequel, we propose two non-parametric approaches to compute fragility curves without relying on the lognormality assumption.

2.2 Binned Monte Carlo simulation

Having at hand a large sample set $\left\{\left(IM_{j},\Delta_{j}\right),\,j=1,\,\dots\,,N\right\}$, it is possible to use binned Monte Carlo simulation (bMCS) to compute the fragility curves, as described next. Let us consider a given abscissa $IM_{o}$. Within a small bin surrounding $IM_{o}$, say $\left[IM_{o}-h,IM_{o}+h\right]$ one assumes that the maximal drift $\Delta$ is linearly related to the $IM$. This assumption is exact in the case of linear structures, but would only be an approximation in the nonlinear case. Therefore, the maximal drift $\Delta_{j}$, which is related to $IM_{j}\in\left[IM_{o}-h,IM_{o}+h\right]$, is converted into the drift $\widetilde{\Delta_{j}}(IM_{o})$, which is related to the $j^{th}$ input signal scaled to have an intensity measure equal to $IM_{o}$:

[TABLE]

This procedure is illustrated in Figure 1. The fragility curve at $IM_{o}$ is then obtained by a crude Monte Carlo estimator:

[TABLE]

where $N_{f}\left(IM_{o}\right)$ is the number of points in the bin such that $\widetilde{\Delta_{j}}(IM_{o})\geqslant\delta_{o}$ and $N_{s}(IM_{o})$ is the total number of points that fall into the bin $\left[IM_{o}-h,IM_{o}+h\right]$.

We note that the bMCS approach bears similarities to the stripe analysis introduced by Shome et al. [34]. However, when using stripe analysis, one scales all ground motions to the intensity level of interest. As a result, certain signals are scaled with factors that are considerably larger or smaller than unity, which may lead to gross approximations of the corresponding responses [35, 36, 37]. The reader is referred to [37] for some illustrations of the effects of the scale factor on the introduced bias, with the latter represented by the ratio of the median nonlinear response of the considered system subject to the scaled motions to the respective median response of the system subject to natural motions with all motions characterized by the same IM level. In general, the bias ratio tends to become larger with increasing deviation of the scale factor from unity. On the other hand, the scaling in binned MCS is confined in the vicinity of the intensity level $IM_{o}$, where the vicinity is defined by the bin width $2h$ chosen so that the scale factors are close to unity. Accordingly, the bias due to ground motion scaling is negligible in bMCS.

Following the above discussion, it should be noted that bias from scaling can be avoided by a proper selection of ground motions. For instance, Shome et al. [34] showed that the scaling of motions that correspond to a narrow interval of earthquake magnitudes and source-to-site distances does not introduce bias into the nonlinear response estimates. Furthermore, Luco and Bazzurro [35] showed that the bias can be reduced by selecting records that have appropriate response spectrum shapes. According to Bazzurro et al. [38] and Vamvatsikos and Cornell [39], the existence of scale-induced bias also depends on several other factors, such as the structural characteristics and the considered intensity and damage measures. The topic of ground motion scaling is complex and falls outside the scope of this paper. We underline that by using the bMCS approach, we avoid introducing bias in the results independently of the ground motion characteristics or other factors. In the following case studies, the resulting fragility curves serve as reference for assessing the accuracy of the various considered techniques for fragility estimation.

2.3 Kernel density estimation

The fragility function defined in Eq. (2) may be reformulated using the conditional probability density function (PDF) $f_{\Delta|IM}$ as follows:

[TABLE]

By definition, this conditional PDF is given as:

[TABLE]

where $f_{\Delta,IM}(\cdot)$ is the joint distribution of the vector $(\Delta,\,IM)$ and $f_{IM}(\cdot)$ is the marginal distribution of the $IM$. If these quantities were known, the fragility function in Eq. (11) would be obtained by a mere integration. In this section, we propose to determine the joint and marginal PDFs from a sample set $\left\{\left(IM_{i},\Delta_{i}\right)\right.$, $\left.i=1,\,\dots\,,N\right\}$ by means of kernel density estimation (KDE).

For a single random variable $X$ for which a sample set $\left\{x_{1},\,\dots\,,x_{N}\right\}$ is available, the kernel density estimate of the PDF reads [40]:

[TABLE]

where $h$ is the bandwidth parameter and $K(\cdot)$ is the kernel function which integrates to one. Classical kernel functions are the Epanechnikov, uniform, normal and triangular functions. The choice of the kernel is known not to affect strongly the quality of the estimate provided the sample set is large enough [40]. In case a standard normal PDF is adopted for the kernel, the kernel density estimate rewrites:

[TABLE]

In contrast, the choice of the bandwidth $h$ is crucial since an inappropriate value of $h$ can lead to an oversmoothed or undersmoothed PDF estimate [41].

Kernel density estimation may be extended to a random vector $\boldsymbol{X}\in{\mathbb{R}}^{d}$ given an i.i.d sample $\left\{\boldsymbol{x}_{1},\,\dots\,,\boldsymbol{x}_{N}\right\}$ [40]:

[TABLE]

where $\boldsymbol{H}$ is a symmetric positive definite bandwidth matrix with determinant denoted by $\left|\boldsymbol{H}\right|$. When a multivariate standard normal kernel is adopted, the joint distribution estimate becomes:

[TABLE]

where $(\cdot)^{\textsf{T}}$ denotes the transposition. For multivariate problems (i.e. $\boldsymbol{X}\in{\mathbb{R}}^{d}$), the bandwidth matrix typically belongs to one of the following classes: spherical, ellipsoidal and full matrix, which respectively contain 1, $d$ and $d(d+1)/2$ independent unknown parameters. The matrix $\boldsymbol{H}$ can be computed by means of plug-in or cross-validation estimators. Both estimators aim at minimizing the asymptotic mean integrated squared error (MISE):

[TABLE]

However, the two approaches differ in the formulation of the numerical approximation of MISE. For further details, the reader is referred to Duong [41]. In the most general case when the correlations between the random variables are not known, the full matrix should be used. In this case, the smoothed cross-validation estimator is the most reliable among the cross-validation methods [42].

Eq. (14) is used to estimate the marginal PDF of the $IM$, namely $\hat{f}_{IM}(a)$, from a sample $\left\{IM_{i},\,i=1,\,\dots\,,N\right\}$:

[TABLE]

Eq. (16) is used to estimate the joint PDF $\hat{f}_{\Delta,IM}(\delta,a)$ from the data pairs $\left\{(IM_{i},\,\Delta_{i}),\,i=1,\,\dots\,,N\right\}$:

[TABLE]

The conditional PDF $f_{\Delta}(\delta|IM=a)$ is eventually estimated by plugging the estimations of the numerator and denominator in Eq. (12). The proposed estimator of the fragility function eventually reads:

[TABLE]

The choice of the bandwidth parameter $h$ and the bandwidth matrix $\boldsymbol{H}$ plays a crucial role in the estimation of fragility curves, as seen in Eq. (20). In the above formulation, the same bandwidth is considered for the whole range of the $IM$ values. However, there are typically few observations available corresponding to the upper tail of the distribution of the $IM$. This is due to the fact that the annual frequency of seismic motions with $IM$ values in the respective range (e.g. $PGA$ exceeding $1g$) is low (see e.g. [43]). This is also the case when synthetic ground motions are used, since these are generated consistently with statistical features of recorded motions. Preliminary investigations have shown that by applying the KDE method on the data in the original scale, the fragility curves for the higher demand thresholds tend to be unstable in their upper tails [44]. To reduce effects from the scarcity of observations at large $IM$ values, we propose the use of KDE in the logarithmic scale, as described next.

Let us consider two random variables $X$, $Y$ with positive supports, and their logarithmic transformations $U=\ln X$ and $V=\ln Y$. One has:

[TABLE]

Accordingly, by substituting $X=IM$ and $Y=\Delta$, the fragility function in Eq. (11) can be obtained in terms of $U=\ln IM$ and $V=\ln\Delta$ as:

[TABLE]

The use of a constant bandwidth in the logarithmic scale is equivalent to the use of a varying bandwidth in the original scale, with larger bandwidths corresponding to larger values of $IM$. The resulting fragility curves are smoother than those obtained by applying KDE with the data in the original scale.

2.4 Epistemic uncertainty of fragility curves

It is of major importance in fragility analysis to investigate the variability in the estimated curves arising due to epistemic uncertainty. This is because a fragility curve is always computed based on a limited amount of data, i.e. a limited number of ground motions and related structural analyses. Large epistemic uncertainties may affect significantly the total variability of the seismic risk assessment outcomes. Characterizing and propagating epistemic uncertainties in seismic loss estimation has therefore attracted attention from several researchers [2, 45, 46].

The theoretical approach to determine the variability of an estimator relies on repeating the estimation with an ensemble of different random samples. However, this approach is not feasible in earthquake engineering because of the high computational cost. In this context, the bootstrap resampling technique is deemed appropriate [2]. Given a set of observations ${\mathcal{X}}=\left(\boldsymbol{X}_{1},\,\dots\,,\boldsymbol{X}_{n}\right)$ of $\boldsymbol{X}$ following an unknown probability distribution, the bootstrap method allows estimation of the statistics of a random variable that depends on $\boldsymbol{X}$ in terms of the observed data ${\mathcal{X}}$ and their empirical distribution [47].

To estimate statistics of the fragility curves with the bootstrap method, we first draw $M$ independent random samples with replacement from the original data set $\left\{\left(IM_{i},\Delta_{i}\right),i=1,\,\dots\,,N\right\}$. These represent the so-called bootstrap samples. Each bootstrap sample has the same size $N$ as the original sample, but the observations are different: in a particular sample, some of the original observations may appear multiple times while others may be missing. Next, we compute the fragility curves for each bootstrap sample using the approaches in Sections 2.1, 2.2 and 2.3. Finally, we perform statistical analysis of the so-obtained $M$ bootstrap curves. In the subsequent example illustration, the above procedure is employed to evaluate the median and 95% confidence intervals of the estimated fragility curves and also, to assess the variability of the $IM$ value corresponding to a $50\%$ probability of failure.

3 Synthetic ground motions

3.1 Properties of recorded ground motions

Let us consider a recorded earthquake accelerogram $a(t)$, $t\in[0,T]$ where $T$ is the total duration of the motion. The peak ground acceleration is $PGA=\max\limits_{\begin{subarray}{c}t\in[0,T]\end{subarray}}\,\lvert a(t)\lvert$. The Arias intensity $I_{a}$ is defined as:

[TABLE]

Defining the cumulative square acceleration as:

[TABLE]

one determines the time instant $t_{\alpha}$ by:

[TABLE]

In addition to the Arias intensity, important properties of the accelerogram in the time domain include the effective duration, defined as $D_{5-95}=t_{95\%}-t_{5\%}$, and the instant $t_{mid}$ at the middle of the strong-shaking phase [49]. Based on the investigation of a set of recorded ground motions, Rezaeian and Der Kiureghian [49] proposed that $t_{mid}$ is taken as the time when $45\%$ of the Arias intensity is reached i.e. $t_{mid}\equiv t_{45\%}$.

Other important properties of the accelerogram are related to its frequency content. Analyses of recorded ground motions indicate that the time evolution of the predominant frequency of an accelerogram can be represented by a linear model, whereas its bandwidth can be considered constant [30]. Rezaeian and Der Kiureghian [30] describe the evolution of the predominant frequency in terms of its value $\omega_{mid}$ at the time instant $t_{mid}$ and the slope of the evolution $\omega^{\prime}$. The same authors describe the bandwidth in terns of the bandwidth parameter $\zeta$. A procedure for estimating the parameters $\omega_{mid}$, $\omega^{\prime}$ and $\zeta$ for a given accelerogram is presented in [30], whereas a simplified version is proposed in [49].

Next, we describe a method for simulating synthetic accelerograms in terms of the set of parameters $\left(I_{a},D_{5-95},t_{mid},\omega_{mid},\omega^{\prime},\zeta_{f}\right)$; this method will be used to generate the seismic motions in a subsequent case study.

3.2 Simulation of synthetic ground motions

The use of synthetic ground motions has been attracting an increasing interest from the earthquake engineering community. This practice overcomes the limitations posed by the small number of records typically available for a design scenario and avoids the need to scale the motions. Use of synthetic ground motions allows one to investigate the structural response for a large number of motions, which is nowadays feasible with the available computer resources (see e.g. [48]).

Different stochastic ground motion models can be found in the literature, which can be classified in three types [49]: record-based parameterized models that are fit to recorded motions, source-based models that consider the physics of the source mechanism and wave travel-path, and hybrid models that combine elements from both source- and record-based models. Vetter and Taflanidis [50] compared the source-based model by Boore [51] with the record-based model by Rezaeian and Der Kiureghian [49] with respect to the estimated seismic risks. It was found that the latter leads to higher estimated risks for low-magnitude events, but the risks are quantified in a consistent manner exhibiting correlation with the hazard characteristics. This model is employed in the present study to generate a large suite of synthetic ground motions that are used to obtain pairs of the ground motion intensity measure and the associated structural response, $(IM,\Delta)$, in order to conduct fragility analysis. The approach, originally proposed in [30], is summarized below.

The seismic acceleration $a(t)$ is represented as a non-stationary process. In particular, the non-stationarity is separated into two components, namely a spectral and a temporal one, by means of a modulated filtered Gaussian white noise:

[TABLE]

in which $q(t,\boldsymbol{\alpha})$ is the deterministic non-negative modulating function, the integral is the non-stationary response of a linear filter subject to a Gaussian white-noise excitation and $\sigma_{h}{(t)}$ is the standard deviation of the response process. The Gaussian white-noise process denoted by $\omega(\tau)$ will pass through a filter $h\left[t-\tau,\boldsymbol{\lambda}(\tau)\right]$, which is selected as the pseudo-acceleration response of a single-degree-of-freedom (SDOF) linear oscillator:

[TABLE]

In the above equation, $\boldsymbol{\lambda}(\tau)=\left(\omega_{f}(\tau),\zeta_{f}(\tau)\right)$ is the vector of time-varying parameters of the filter $h$, with $\omega_{f}(\tau)$ and $\zeta_{f}(\tau)$ respectively denoting the filter’s natural frequency and damping ratio at instant $\tau$. Note that $\omega_{f}(\tau)$ corresponds to the evolving predominant frequency of the ground motion represented, while $\zeta_{f}(\tau)$ corresponds to the bandwidth parameter of the motion. As noted in Section 3.1, $\zeta_{f}(\tau)$ may be taken as a constant $\left(\zeta_{f}(\tau)\equiv\zeta\right)$, while $\omega_{f}(\tau)$ is approximated as a linear function:

[TABLE]

where $t_{mid}$, $\omega_{mid}$ and $\omega^{\prime}$ are as defined in Section 3.1. After being normalized by the standard deviation $\sigma_{h}{(t)}$, the integral in Eq. (26) becomes a unit-variance process with time-varying frequency and constant bandwidth. The non-stationarity in intensity is then captured by the modulating function $q(t,\boldsymbol{\alpha})$, which determines the shape, intensity and duration $T$ of the signal. This is typically described by a Gamma-like function [49]:

[TABLE]

where $\boldsymbol{\alpha}=\left\{\alpha_{1},\alpha_{2},\alpha_{3}\right\}$ is directly related to the energy content of the signal through the quantities $I_{a}$, $D_{5-95}$ and $t_{mid}$ defined in Section 3.1 (see [49] for details).

For computational purposes, the acceleration in Eq. (26) can be discretized as follows:

[TABLE]

where the standard normal random variable $U_{i}$ represents an impulse at instant $t_{i}=i\times\dfrac{T}{n}\,,\,i=1,\,\dots\,,n$, ($T$ is the total duration) and $s_{i}(t,\boldsymbol{\lambda}(t_{i}))$ is given by:

[TABLE]

As a summary, the considered seismic motion generation model consists of the three temporal parameters $\left(\alpha_{1},\alpha_{2},\alpha_{3}\right)$, which are related to $\left(I_{a},D_{5-95},t_{mid}\right)$, the three spectral parameters $\left(\omega_{mid},\omega^{\prime},\zeta_{f}\right)$ and the standard Gaussian random vector $\boldsymbol{U}$ of size $n$. Rezaeian and Der Kiureghian [49] proposed a methodology for determining the temporal and spectral parameters according to earthquake and site characteristics, i.e. the type of faulting of the earthquake (strike-slip fault or reverse fault), the closest distance from the recording site to the ruptured area and the shear-wave velocity of the top 30 m of the site soil. For the sake of simplicity, in this paper these parameters are directly generated from the statistical models given in [49], which are obtained from analysis of a large set of recorded ground motions.

4 Steel frame structure subject to synthetic ground motions

4.1 Problem setup

We determine the fragility curves for the three-storey three-span steel frame shown in Figure 2. The dimensions of the structure are: storey-height $H=3$ m, span-length $L=5$ m. The vertical load consists of dead load (weight of frame elements and supported floors) and live load (in accordance with Eurocode 1 [52]) resulting in a total distributed load on the beams $q=20$ kN/m. In the preliminary design stage, the standard European I beams with designation IPE 300 A and IPE 330 O are chosen respectively for the beams and columns. The steel material has a nonlinear isotropic hardening behavior following the uniaxial Giuffre-Menegotto-Pinto steel model as implemented in the finite element software OpenSees [53]. Ellingwood and Kinali [5] have shown that uncertainty in the properties of the steel material has a negligible effect on seismic fragility curves. Therefore, the mean material properties are used in the subsequent fragility analysis: $E_{0}=210,000$ MPa for the Young’s modulus (initial elastic tangent in the stress-strain curve), $f_{y}=264$ MPa for the yield strength [54, 55] and $b=0.01$ for the strain hardening ratio (ratio of post-yield to initial tangent in the stress-strain curve). Figure 2 depicts the hysteretic behavior of the steel material at a specified section for an example ground motion. The structural components are modelled with nonlinear force-based beam-column elements characterized by distributed plasticity along their lengths, while use of fiber sections allows modelling the plasticity over the element cross-sections [56]. The connections between structural elements are modeled with rigid nodes. The first two natural periods of the building obtained by modal analysis are $T_{1}=0.61$ s and $T_{2}=0.181$ s, corresponding to natural frequencies $f_{1}=1.64$ Hz and $f_{2}=5.53$ Hz. Rayleigh damping is considered with the damping ratio of the first two modes set equal to 2%.

The structure is subject to seismic motions represented by synthetic acceleration time histories at the ground level. Each time history is modelled in terms of six randomized parameters $\left(\alpha_{1},\alpha_{2},\alpha_{3},\omega_{mid},\omega^{\prime},\zeta_{f}\right)$ directly related to the parameters in Table 1 and a Gaussian input vector $\boldsymbol{U}$ as described in Section 3. The statistics of the parameters in Table 1 are taken from [49]; in the latter study, the authors derived the listed distributions and associated parameters by analyzing a set of recorded ground motions corresponding to strong strike-slip and reserve earthquakes with moment magnitudes in the range 6-8 and rupture distances in the range of 10-100 km. The reader is referred to [49] for viewing the correlations between these parameters. The duration of each time history is computed from the corresponding set of parameters $\left(\alpha_{1},\alpha_{2},\alpha_{3}\right)$ and is used to determine the size of the Gaussian vector $\boldsymbol{U}$. Two example synthetic acceleration time histories are shown in Figure 3. Transient dynamic analyses of the frame are carried out for a total of $N=20,000$ synthetic motions using the finite element software OpenSees.

Numerous types of $IM$ can be used to describe the earthquake severity, see e.g. [57]. Peak ground acceleration ($PGA$) is a convenient measure that is straightforward to obtain from a given time history and has been traditionally used in attenuation relationships and design codes. However, structural responses may exhibit large dispersions for a certain $PGA$, since they are also highly dependent on other features of earthquake motions, e.g. the frequency content and duration of the strong motion phase. Structure-specific $IM$s, such as the spectral acceleration $Sa$ and the pseudo spectral acceleration $Psa$, tend to be better correlated with structural responses [57, 58]. In the following, we compute fragility curves considering both $PGA$ and $Sa$ as $IM$s. $Sa$ represents $Sa(T_{1})$ i.e. the spectral acceleration for a single-degree-of-freedom system with period equal to the fundamental period $T_{1}$ of the frame and viscous damping ratio equal to 2%.

The engineering demand parameter commonly considered in fragility analysis of steel buildings is the maximal inter-storey drift ratio, i.e. the maximal difference of horizontal displacements between consecutive storeys normalized by the storey height (see e.g. [5, 59, 60]). Accordingly, we herein develop fragility curves for three different thresholds of the maximal inter-storey drift ratio over the frame. To gain insight into structural performance, we consider the thresholds 0.7%, 1.5% and 2.5%, which are associated with different damage states in seismic codes. In particular, the thresholds 0.7% and 2.5% are recommended in [61] to respectively characterize light and moderate damage for steel frames, while the threshold 1.5% corresponds to the damage limitation requirement for buildings with ductile non-structural elements according to Eurocode 8 [62]. These descriptions only serve as rough damage indicators, since the relationship between drift limit and damage in the PBEE framework is probabilistic.

4.2 Fragility curves

As described in Section 2, the lognormal approach relies on assuming that the fragility curves have the shape of a lognormal CDF and estimating the parameters of this CDF. Using the maximum likelihood estimation (MLE) approach, the observed failures for each drift threshold are modeled as outcomes of a Bernoulli experiment and the parameters ($\alpha$, $\beta$) of the fragility curves are determined by maximizing the respective likelihood function. Using the linear regression (LR) technique, the parameters of the lognormal curves are derived by fitting a linear model to the paired data $(\ln IM,\,\ln\Delta)$. Figure 4 depicts the paired data $(\ln PGA,\,\ln\Delta)$ and $(\ln Sa,\,\ln\Delta)$ together with the fitted models based on linear regression. It can be seen that a single linear model is not appropriate for the cloud of points $(\ln Sa,\,\ln\Delta)$ and thus, bilinear regression is used in this case (see also [63, 64, 65] for use of a similar model). The break point in the bilinear model ($Sa=0.45\,g$) is determined according to the method presented in [66] using the R package segmented. When $PGA$ is used as $IM$, the coefficient of determination of the fitted linear model is $R^{2}=0.663$; when $Sa$ is used as $IM$, it is $R_{1}^{2}=0.978$ and $R_{2}^{2}=0.785$ for the first and second part of the bilinear model, respectively. Note that use of $Sa$ as $IM$ leads to a smaller dispersion, i.e. a smaller $\zeta$ in Eq. (6), as compared to $PGA$; this is expected since $Sa$ is a structure-specific $IM$. In the bMCS method, the bin width $h$ is set equal to $0.25\,IM_{o}$. The resulting scale factors vary in the range $[0.75,\,1.25]$ corresponding to a bias ratio approximately equal to unity [37]. The KDE approach requires estimation of the bandwidth parameter and the bandwidth matrix. Using the cross-validation estimation implemented in R [67], these are determined as $h=0.133$, $\boldsymbol{H}=\begin{bmatrix}0.031&0.024;~{}0.024&0.027\end{bmatrix}$ when $PGA$ is used as $IM$, and $h=0.155$, $\boldsymbol{H}=\begin{bmatrix}0.023&0.023;~{}0.023&0.024\end{bmatrix}$ when $Sa$ is used as $IM$.

For the two types of $IM$ and the three drift limits considered, Table 2 lists the medians and log-standard deviations of the lognormal curves obtained with both the MLE and LR approaches. The median determines the position where the curve attains the value 0.5, whereas the log-standard deviation is a measure of the steepness of the curve. Note that the MLE approach yields a distinct log-standard deviation for each drift threshold, whereas a single log-standard deviation for all drift thresholds is obtained with the LR approach. The medians of the KDE-based curves are also computed and are listed in Table 2 for comparison. The KDE-based medians, which serve as the reference values, may be overestimated or underestimated by the lognormal approach, depending on the method used to estimate the parameters, the considered $IM$ and the drift threshold; the absolute deviations tend to be larger for larger drift thresholds.

For the case when $PGA$ is considered as $IM$, Figure 5 (left) shows the fragility curves obtained with the MLE- and LR-based lognormal approaches and the bMCS- and KDE-based non-parametric approaches. One first observes a remarkable consistency between the curves obtained with the two non-parametric approaches despite the distinct differences in the underlying algorithms. This validates the accuracy of the proposed methods. For the lower threshold ($\delta_{o}=0.7\%$), both parametric curves are in good agreement with the non-parametric ones. For the two higher thresholds, the LR-based lognormal curves exhibit significant deviations from the non-parametric ones leading to an overestimation of the failure probabilities. Note that for $\delta_{o}=1.5\%$ and $\delta_{o}=2.5\%$, the median $PGA$ (leading to 50% probability of exceedance) is respectively underestimated by 19% and 15% when the LR aproach is used (see Table 2). In contrast, the MLE-based lognormal curves are in a fair agreement with their non-parametric counterparts with the largest discrepancies observed for the highest threshold $\delta_{o}=2.5\%$.

Figure 5 (right) shows the resulting fragility curves when $Sa$ is considered as $IM$. The non-parametric curves based on bMCS and KDE remain consistent independently of the drift threshold. For $\delta_{o}=0.7\%$, the fragility curves are steep, which is due to the strong correlation between $Sa$ and $\Delta$ when the structure behaves linearly. For this threshold, the LR-based curve is closer to the non-parametric curves than the MLE-based one. For the two larger thresholds, the MLE-based curves are fairly accurate, whereas the LR-based curves exhibit significant deviations from their non-parametric counterparts. In particular, the LR-based curves overestimate the failure probabilities for $\delta_{o}=1.5\%$ and underestimate the failure probabilities for $\delta_{o}=2.5\%$. Note that for $\delta_{o}=1.5\%$, the median $Sa$ is underestimated by 4%, whereas for $\delta_{o}=1.5\%$, the median $Sa$ is overestimated by 8% when the LR aproach is used (see Table 2).

Summarizing the above results, the MLE-based lognormal approach yields fragility curves that are overall close to the non-parametric ones; however, it smooths out some details of the curves that can be obtained with the non-parametric approaches. On the contrary, the LR-based lognormal curves can be highly inaccurate. As noted in Section 2.1.2, the LR approach assumes that the residuals of the fitted model in the log-scale (Eq. (6)) follow a normal distribution with a constant standard deviation independently of the $IM$ level. Figure 6 shows histograms of $\ln\Delta$ at two example levels of $PGA$ and $Sa$ together with the fitted normal distributions according to Eq. (6). The responses $\Delta$ at each $IM$ level are obtained consistently with the bMCS approach. Obviously, the assumption of a normal distribution is not valid, which is more pronounced when $Sa$ is used as $IM$. This explains the inaccuracy of the LR-based fragility curves for both types of $IM$, despite the relatively high coefficients of determination of the fitted models in the case of $Sa$.

4.3 Estimation of epistemic uncertainty by bootstrap resampling

In the following, we use the bootstrap resampling technique (see Section 2.4) to investigate the epistemic uncertainty in the fragility curves estimated with the proposed non-parametric approaches.

We examine the stability of the estimated curves by comparing those with the bootstrap medians, and the variability in the estimation by computing bootstrap confidence intervals. For the two considered $IM$s and the three drift thresholds of interest, Figure 7 shows the median bMCS- and KDE-based fragility curves and the 95% confidence intervals obtained by bootstrap resampling with 100 replications together with the respective estimated curves (also shown in Figure 5). Figure 7 clearly shows that both the bMCS-based and the KDE-based median fragility curves obtained with the bootstrap method do not differ from the curves estimated with the original set of observations. This shows the stability of the proposed approaches. For a specified $IM$ and drift limit, the confidence intervals of the bMCS- and KDE-based curves have similar widths. The interval widths tend to increase with increasing drift limit and increasing $IM$ value.

In order to quantify the effects of epistemic uncertainty, one can estimate the variability of the median $IM$, i.e. the $IM$ value leading to 50% probability of exceedance. Assuming that the median $IM$ ($PGA$ or $Sa$) follows a lognormal distribution [68], the median $IM$ is determined for each bootstrap curve and the log-standard deviation of the distribution of the median is computed. Table 3 lists the log-standard deviations of the median $IM$ values for the same cases as in Figure 7. These results demonstrate that epistemic uncertainty is increasing with increasing threshold $\delta_{o}$. In all cases, the log-standard deviations are relatively small indicating a low level of epistemic uncertainty, which is due to the large number of transient analyses ($N=20,000)$ considered in this study. Although use of such large sets of ground motions is not typical in practice, it is useful for the refined analysis presented here.

5 Concrete column subject to recorded ground motions

To demonstrate the comparison between the lognormal and the non-parametric approaches for the case when fragility curves are based on recorded ground motions, we herein briefly summarize a case study by the authors originally presented in [29].

In this study, we estimate the fragility of a reinforced concrete column with a uniform circular cross-section, representing a pier of a typical California highway overpass bridge [63] (see Figure 8). The column is modelled in the finite element code OpenSees as a fiberized nonlinear beam-column element. For details on the modelling of the concrete material and the steel reinforcement, the reader is referred to [29]. The loading-unloading behavior and the pushover curve of the column are shown in Figure 8. Three-dimensional time-history analyses of the bridge column are conducted for $N=531$ earthquake records (each comprising three orthogonal component accelerograms). These records are obtained from the PEER strong motion database and cover a wide range of source-to-site distances and earthquake moment magnitudes [63]. The developed fragility curves represent the probability of the maximal drift ratio $\Delta$ in the transverse direction exceeding specified thresholds $\delta_{o}$ as a function of the peak ground acceleration $PGA$ or the pseudo-spectral acceleration $Psa$ corresponding to the first transverse mode ($T_{1}=0.535~{}s$). The considered drift ratio thresholds, shown in Table 4, are recommended for the operational and life safety levels by two different sources [69, 70].

The fragility curves are established with the MLE-based and LR-based lognormal approaches and the KDE-based non-parametric approach. Due to the relatively small number of data, the bMCS method is not considered herein. Figure 9 depicts the clouds of points $(IM_{i},\,\Delta_{i})$ in the logarithmic scale for the two $IM$s together with the linear fitted models. The coefficients of determination of the latter are $R^{2}=0.729$ for the case of $PGA$ and $R^{2}=0.963$ for the case of $Psa$. Note that for small values of $Psa$ ($Psa<$0.2 g) a linear function provides a perfect fit, which is due to the fact that in this range of $Psa$, the column behavior can be represented by a linear single-degree-of-freedom model. Further details on the parameters of the different fragility functions can be found in [29].

Figure 10 depicts the obtained fragility curves for the two types of $IM$ and the four drift-ratio thresholds of interest. When $PGA$ is used as $IM$, the curves obtained with the two lognormal approaches are close to each other for all thresholds, but exhibit deviations from the non-parametric curves, which tend to be larger for higher PGA levels and larger drift limits. When $Psa$ is used as $IM$, the MLE-based curves are in a fair agreement with the KDE-based ones; however, the former smooth out some details that can be obtained with the non-parametric approach. In contrast, the LR-based curves are inaccurate for all but the smaller drift threshold. Overall, the LR-based curves exhibit larger deviations from the non-parametric ones for $Psa$ than for $PGA$ as $IM$, although the $R^{2}$ coefficient of the linear fit is higher for $Psa$. This can be explained by the fact that the assumption of homoscedastic errors, inherent in Eq. (6), is not valid for the specific data set $(Psa,\Delta)$, as one can observe in Figure 9.

6 Discussion

Using the non-parametric fragility curves as reference, the accuracy of the lognormal curves is found to depend on the method used to estimate the parameters of the underlying CDF, the considered $IM$ and the drift threshold of interest. In most cases, the MLE-based curves are fairly close to the non-parametric ones, whereas the LR-based curves exhibit significant deviations. The lognormal curves tend to deviate more from the non-parametric ones for larger drift limits. Considering both case studies, the MLE-based curves are more accurate for a structure-specific $IM$ ($Sa$, $Psa$) than for $PGA$. Different $IM$s have been recommended in the literature for structures of different type, size and material [58, 71]. Accordingly, the accuracy of the lognormal fragility curves may depend on those factors as well. Possible dependence of the accuracy of the lognormal curves on the considered response quantity needs to be investigated as well.

As noted in Section 2.2, the bMCS approach bears similarities with the so-called stripe analysis [34, 38, 72]. A comparison between the stripe and cloud analyses, where the latter corresponds to the LR-based lognormal approach [59, 73, 74, 75], was carried out by Celik and Ellingwood [74]. In the mentioned study, concrete structures were subject to 40 synthetic ground motions. Differences in the response statistics obtained with the two methods for three $IM$ levels were found insignificant and hence, use of the cloud analysis was justified. In contrast, Baker [73] showed that cloud analysis can significantly underestimate the mean annual rate of exceeding a large maximum interstory drift. Karamlou and Bocchini [28] recently conducted large-scale simulations on bridge structures in order to investigate the underlying assumptions of the cloud analysis. Their results showed that, in general, the conditional distribution of a demand parameter for a given $IM$ level is not lognormal. In addition, it was found that the assumptions of a linear function for the probabilistic seismic demand model in the log-scale (power function in the normal scale) and of constant dispersion of the respective errors can lead to significant errors in fragility analysis. These findings are consistent with our results shown in Figure 6. The limitations of the LR-based approach have also been mentioned by Jalayer et al. [75].

Based on the results of our case studies and the above discussion, we recommend the use of the MLE approach if fragility curves are developed in a parametric manner. The superiority of the MLE over the LR approach relies on the fact that the former avoids the assumptions of the linear model and the homoscedasticity of the errors that are inherent in the latter. However, when a detailed description of the fragility function is important, a non-parametric approach should be used. The bMCS method requires a large number of data, which can be typically obtained by use of synthetic motions; note that due to the current computer capacities and the use of distributed computing, large-scale simulations are becoming increasingly popular among both researchers and practitioners. On the other hand, the KDE approach can be employed even with a limited number of recorded motions at hand, as shown in our second case study. We again emphasize that the two non-parametric approaches lead to almost identical curves in the case when they could be applied independently with the same (large) dataset.

7 Conclusions

Seismic demand fragility evaluation is one of the basic elements in the framework of performance-based earthquake engineering (PBEE). At present, the classical lognormal approach is widely used to establish such fragility curves mainly due to the fact that the lognormality assumption makes seismic risk analysis more tractable. The approach consists in assigning the shape of a lognormal cumulative distribution function to the fragility curves. However, the validity of this assumption remains an open question.

In this paper, we introduce two non-parametric approaches in order to examine the validity of the classical lognormal approach, namely the binned Monte Carlo simulation and the kernel density estimation. The former computes the crude Monte Carlo estimators for small subsets of ground motions with similar values of a selected intensity measure, while the latter estimates the conditional probability density function of the structural response given the ground motion intensity measures using the kernel density estimation technique. The proposed approaches can be used to compute fragility curves when the actual shape of these curves is not known as well as to validate or calibrate parametric fragility curves. Herein, the two non-parametric approaches are confronted to the classical lognormal approach in two case studies, considering synthetic and recorded ground motions.

In the case studies, the fragility curves are established for various drift thresholds and different types of the ground motion intensity measure, namely the peak ground acceleration ($PGA$), and the structure-specific spectral acceleration ($Sa$) and pseudo-spectral acceleration ($Psa$). The two non-parametric curves are always consistent, which proves the validity of the proposed techniques. Accordingly, the non-parametric curves are used as reference to assess the accuracy of the lognormal curves. The parameters of the latter are estimated with two approaches, namely by maximum likelihood estimation and by assuming a linear probabilistic seismic demand model in the log-scale. The maximum likelihood estimation approach is found to approximate fairly well the reference curves in most cases, especially when a structure-specific intensity measure is used; however, it smooths out some details that can be obtained with the non-parametric approaches. In contrast, the assumption of a linear demand model in the log-scale is found overall inaccurate. When integrated in the PBEE framework, inaccuracy in fragility estimation may induce errors in the probabilistic consequence estimates that serve as decision variables for risk mitigation actions. The bootstrap resampling technique is employed to assess effects of epistemic uncertainty in the non-parametric fragility curves. Results from bootstrap analysis validate the stability of the fragility estimates with the proposed non-parametric methods.

Recently, fragility surfaces have emerged as an innovative way to represent the vulnerability of a system [13]; these represent the failure probability conditional on two intensity measures of the earthquake motions. The computation of these surfaces is not straightforward and requires a large computational effort. The present study opens new paths for establishing the fragility surfaces: similarly to the case of fragility curves, one can use kernel density estimation to obtain fragility surfaces that are free of the lognormality assumption and consistent with the surfaces obtained by Monte Carlo simulation.

We note that the computational cost of the two proposed approaches is significant when they are based on large Monte Carlo samples. In order to reduce this cost, alternative approaches may be envisaged. Polynomial chaos (PC) expansions [76, 77] appear as a promising tool. Based on a smaller sample set (typically a few hundreds of finite element runs), PC expansion provides a polynomial approximation that surrogates the structural response. The feasibility of post-processing PC expansions in order to compute fragility curves has been shown in [78, 79] in the case a linear structural behavior is assumed. The extension to nonlinear behavior is currently in progress.

We underline that the proposed non-parametric approaches are essentially applicable to other probabilistic models in the PBEE framework, relating decision variables with structural damage and structural damage with structural response. Once all the non-parametric probabilistic models are available, they can be incorporated in the PBEE framework by means of numerical integration. Then a full seismic risk assessment may be conducted by avoiding potential inaccuracies introduced from simplifying parametric assumptions at any step of the analysis. Optimal high-fidelity computational methods for incorporating non-parametric fragility curves in the PBEE framework will be investigated in the future.

Acknowledgements

The authors are thankful to the two anonymous reviewers for various valuable comments that helped improve the quality of the manuscript. Discussions with Dr. Sanaz Rezaeian, who provided clarifications on the stochastic ground motion model used in this study, are also acknowledged.