基于贝叶斯分层模型的液化侧移稳健的易损性分析方法

摘要

提出了一种基于贝叶斯分层模型的液化侧移稳健的易损性分析方法。采用贝叶斯分层模型量化不同增量动力分析（IDA）曲线间的差异，结合抽样方法预测潜在侧移的分布，建立液化侧移稳健的易损性曲线和超越概率曲线。以一处实际发生过液化侧移的场地为例，展示了稳健的易损性曲线及超越概率曲线的建立方法，并与相关方法进行比较。结果表明，所提出的方法可以较好地模拟IDA曲线的分布，较为准确地量化易损性曲线和超越概率曲线的不确定性。

关键词

液化侧移; 贝叶斯分层模型; 稳健的易损性分析; 基于性能的抗震设计; 增量动力分析（IDA）

地震引起的土体液化是城市地区常见的地震次生灾害，其表现包括砂沸、地表沉降、侧向变形甚至液化流滑。土体侧移和液化滑流通常发生在水平微倾场地或附近具有自由面的场地，易在城市地区造成大规模破坏和损失。液化侧移导致的巨大破坏引发了国际学界对土体液化研究的重视，液化侧移的风险评估问题已成为地震和岩土工程领域的研究热点^［1］。

液化侧移的易损性分析以侧移大小为状态指标，侧移超过目标限值即达到失效状态。增量动力分析（IDA）^［2］是目前使用最广泛的建立失效状态与地震动参数概率关系的方法^［3-5］。然而，目前传统增量动力分析法在建立液化侧移易损性曲线的过程中存在一些问题：由于数值分析方法计算量较大，因此通常只能采用有限数量的模拟样本，而基于少数样本建立的易损性曲线可能具有较高的不确定性。目前传统增量动力分析法无法量化有限数量样本产生的估计不确定性，更无法准确量化易损性曲线和超越概率曲线的不确定性。此外，传统增量动力分析法不能保证预测结果在邻近侧移的连续性，不同侧移的易损性曲线可能发生交错的现象^［6］，也不能保证符合概率逻辑的小位移超越概率必须大于等于大位移超越概率，导致无法获得连贯统一的侧移超越概率曲线。因此，基于增量动力分析法建立考虑多种不确定性且在邻近位移连续稳定的液化侧移易损性曲线值得进一步研究。

对于传统增量动力分析法，Vamvatsikos等^［2］分别拟合单个IDA曲线，再对拟合参数进行统计，将参数的统计值代表IDA曲线簇的形状。然而，传统增量动力分析法生成的数据通常以地震波为单位成组排列，不同地震波生成的IDA曲线通常具有较大差异^［6-8］，存在明显的组间变异性。此外，同一条地震波缩放到不同强度时产生的响应也不完全服从特定模式，可能存在诸如“非单调”“回折”“平台阶段”等现象^［9］。传统增量动力分析法难以量化IDA曲线的组间和组内变异性，无法考虑样本数量对组间和组内变异性估计的影响，而贝叶斯分层模型非常适合处理这类成组排列且存在组间变异性的数据。贝叶斯分层模型通过分层架构区分组间变异性与组内变异性，同时通过贝叶斯方法考虑样本数量与质量对参数估计不确定性的影响。

针对采用传统增量动力分析法进行液化侧移易损性分析遇到的问题，本文提出了一种基于贝叶斯分层模型的液化侧移稳健的易损性分析方法。首先，介绍了液化侧移的易损性分析与传统增量动力分析法，并结合具体案例展示传统增量动力分析法应用在液化侧移易损性分析中的不足；然后，介绍了稳健的易损性分析的概念与方法，包括Jalayer等^［10］提出的直接回归法和本文提出的贝叶斯分层模型法；最后，将2种方法建立的液化侧移稳健的易损性曲线和超越概率曲线进行了比较。

1 液化侧移的易损性分析与传统增量动力分析法

1.1　液化侧移的易损性分析

液化侧移的易损性分析是基于性能的液化侧移分析的重要组成，而基于性能的液化侧移分析源于Cornell^［11］提出的地震风险评估方法，目前属于基于性能的抗震设计方法（PBEE）的一个分支。基于性能的抗震设计方法是一种评估工程设施在不同来源、强度地震动下性能的概率方法^［12-13］。太平洋地震工程研究中心（PEER）^［14］提出的基于性能的抗震设计框架可以表示为

$$\lambda (S)=\iiint G(S|D)\mathrm{d}G(D|d)\mathrm{d}G(d|I)\mathrm{d}\lambda (I)$$

（1）

式中：G（d|I）为给定地震动强度I的情况下工程需求参数d超过限值的概率函数；D为损伤度量；S为决策变量；λ（I）为地震动强度I的年超越概率，由地震动风险模型计算得到。

对于基于性能的液化侧移分析，液化侧移d_x即为工程需求参数d。液化侧移达到限值D_x的年平均概率λ（d_x=D_x）即为液化侧移的超越概率λ（d_x），可以表示为：

$$\lambda (d_x)=\lambda (d_x=D_x)=\int G(d_x|I)\mathrm{d}\lambda (I)=$$

$$\int P(d_x\ge D_x|I)\mathrm{d}\lambda (I)$$

（2）

式中：$P(d_x\ge D_x|I)$为给定地震动强度I下模型响应超过限值的概率。G（d_x|I）又可以写成P（d_x≥D_x|I），即液化侧移易损性分析的目标函数，表征给定地震动强度I下侧移超过目标限值的概率。

1.2　传统增量动力分析法

增量动力分析法是最常用的建立易损性模型的方法。一般增量动力分析法计算流程为：选择一组与目标场地具有相似性质的地震波数据集，通过缩放将所有地震波调整到目标地震动强度水平，然后将地震波输入数值模型并进行非线性动力响应分析，每条地震波的数据汇总后建立IDA曲线。增量动力分析法采用多条地震波，并考虑了地震波的变异性。此外，土体参数的不确定性也会对模型响应产生较大影响^［15］。为此，Vamvatsikos^［8］提出了基于拉丁超立方抽样^［16］的增量动力分析法，该方法可以同时考虑地震波和土体参数不确定性。

传统建立易损性曲线的增量动力分析法分为2种^［6］（下文简称方法1和方法2）。方法1，假定模型响应d_x在达到限值D_x时，地震动强度I服从对数正态分布。根据限值对应的地震动强度I统计数据点，计算临界地震动强度I的中位值η和对数标准差β（见文献［17］），再根据下式计算给定地震动强度I下模型响应超过限值的概率：

$$P(d_x\ge D_x|I)=P(I\ge I_D)=\Phi \left(\frac{\mathrm{l}\mathrm{n}I-\mathrm{l}\mathrm{n}\eta }{\beta }\right)$$

（3）

式中：Ф为标准正态分布的累积分布函数。方法2，直接计算给定地震动强度I下d_x统计数据点超过限值的概率：

$$P(d_x\ge D_x|I)=\frac{n}{N}$$

（4）

式中：n为d_x统计数据点中超过限值的样本数量；N为d_x统计数据点的总数量。然后，根据式（4）计算不同地震动强度I下的概率样本点，采用对数正态分布的累积分布函数^［17-18］等拟合概率样本点，得到易损性曲线。

1.3　液化侧移计算案例

由于国内暂时缺乏较为完备的地震危险性分析数据，选择美国加州的一个微倾场地作为示例。场地位于紧邻美国第二大海岸风积沙丘Guadalupe-Nipomo Dunes的Oceano。在2003年12月22日发生的M_w6.5地震中，该场地发生了较大规模的液化侧移破坏^［19］，包含静力触探试验数据的场地典型地质剖面如图1a所示。图中，q_c1n为修正后静力触探锥尖阻力。

图1  案例场地

Fig.1  Illustrative example

a 地质剖面图 b 数值模型剖面图

采用多屈服面塑性砂本构模型P2PSand^［20］模拟场地的液化侧移。基于LEAP试验中动三轴及离心机试验标定的模型参数^［21］如表1所示，其中循环参数K_c根据El-Sekelly等^［21］提供的2组数据由相对密实度线性插值计算得到。

表1  数值模拟模型参数

Tab.1  Model parameters for numerical simulation

模型参数	取值
循环参数，Kc	线性插值
临界摩擦因数	33
剪胀因数	1
临界状态参数	9
弹性参数	772
退化因数	0.19
泊松比	0.1
塑性剪切比	0.4
渗透系数/（cm·s-1）	1.2×10-3

模型布局如图1b所示，表面倾角为5%。模型中心顶部（sand₀）的位移即为场地表面的液化侧移d_x。模型最大网格尺寸为1 m^［15，22］，小于地震波最小波长的1/10，可以较为准确地模拟地震波在模型介质中的传播^［23］。模型底部采用吸收边界，侧面边界采用绑定边界^［24-25］。模型底部施加了原始未旋转的两分量水平地震波，以考虑地震动方位角与场地倾向之间的不确定性。

增量动力分析时需要选择地震波数据集。由于数值模拟计算量较大，因此选择20组地震波，将20条地震波与20组土体参数样本随机匹配^［8］。从NGA数据库^［26］中选择30 m平均剪切波速v_s30在200~300 m·s^-1范围内的20条地震波，编号分别为31、126、158、160、162、183、430、437、439、461、507、515、719、720、766、992、996、1001、1088、1101，反应谱如图2所示。将上述地震波的二维水平分量输入模型，采用GMRotI50^［27］度量二维水平地震动的强度。选择地表峰值速度v作为目标缩放地震动强度参数^［28］。根据地表峰值速度v的年超越概率分布，选择1、2、5、10、20、30、40、50、60、80、100、120、140 cm·s^-1 13个等级作为v的缩放目标，相应的地表峰值速度v从0.007g至3.740g不等。

图2  本文采用的地震波反应谱

Fig.2  Response spectrum of selected seismic waves

1.4　传统增量动力分析法进行液化侧移易损性分析的不足

根据传统增量动力分析法计算得到的易损性曲线如图3所示，将观测失效概率样本点作为比较。失效概率由于仅由20个状态样本确定，因此具有较高的不确定性。失效概率的真实范围可用威尔逊置信区间^［29］表示，图3中灰色竖线代表了失效概率的90%威尔逊置信区间。

图3  方法1、2建立的液化侧移易损性曲线

Fig.3  Fragility curves for liquefaction-induced lateral displacement based on method 1，2

从图3可以看出，对于方法1建立的易损性曲线，不少观测概率点的90%置信区间与易损性曲线重合较小。此外，相比于方法2，方法1建立的易损性曲线在v较小时具有较小的破坏概率，在v较大时具有较大的破坏概率，曲线中段的斜率较大，说明式（3）中参数β较小，即易损性曲线考虑的地震动强度I不确定性较小。这可能与方法1在样条曲线插值过程中忽略插值的不确定性有关。由于方法2直接根据观测概率样本点拟合易损性曲线，因此易损性曲线与观测概率符合程度更高。另一方面，在观测点最小化拟合误差的过程中，观测点的微小变化可能会引起易损性曲线形状的突变。图4a展示了方法2侧移限值为0.15 m和0.20 m时的易损性曲线。可以明显地看到，2条易损性曲线发生交错，ln v=3时侧移超过0.20 m的概率甚至大于0.15 m的概率。事实上，当侧移超过0.20 m时则必然也超过0.15 m，上述现象违背了概率常识。值得注意的是，此时方法2仅使用增量动力分析法分析原始的13个地表峰值速度等级的失效概率样本点。为了获得更多的失效概率样本点来拟合易损性曲线，可对各条IDA曲线进行线性插值。通过插值获得了9个地表峰值速度等级的失效概率样本点。图4b比较了采用22个与13个样本点拟合的易损性曲线。可以看到，在地表峰值速度v较小时，采用22个样本点建立的易损性曲线预测失效概率相比13个样本点的易损性曲线显著提高，说明方法2受失效概率样本点的数量与位置影响较大。将图4b中22个样本点预测的0.20 m的易损性曲线与图4a中0.15 m的易损性曲线进行比较，ln v=3时侧移0.20 m的预测失效概率依然更高，说明增加样本点的数量也无法消除不同侧移易损性曲线交错的现象，无法保证符合概率逻辑的小位移超越概率大于等于大位移超越概率。

图4  基于方法2建立的液化侧移易损性曲线

Fig.4  Fragility curves of liquefaction-induced lateral displacement based on method 2

上述2种方法虽然建立了易损性曲线，但是无法量化参数的不确定性，也无法考虑易损性曲线的不确定性。因此，提出了一种能充分量化易损性曲线不确定性的易损性分析方法。

2 稳健的易损性分析

2.1　稳健的易损性曲线与超越概率曲线

数值分析或者模型试验中通常只会获得有限数量的样本，若样本数量较小，则基于有限样本建立的易损性曲线通常具有较大的不确定性。稳健的易损性分析由Jalayer等^［10］于2015年提出，考虑了有限样本产生的估计不确定性及地震波、土体参数等其他不确定性的影响。令D为有限的观测样本，稳健的易损性分析采用贝叶斯方法，计算易损性曲线参数χ的概率分布f（χ|D），对于易损性曲线P（d_x≥D_x|I，χ）在参数空间Ω_χ内按参数概率分布f（χ|D）积分，得到稳健的易损性曲线与标准差分别为^［10，30］：

$$P(d_x\ge D_x|I,\mathit{D})=\underset{{\Omega }_{\chi }}{\int }P(d_x\ge D_x|I,\chi )f(\chi |D)\mathrm{d}\chi =\underset{{\Omega }_{\mathit{\chi }}}{\mathrm{E}}(P(d_x\ge D_x|I,\mathit{D},\mathit{\chi }))$$

（5）

$$\begin{array}{l}{\sigma }^{\mathrm{2}}(P(d_x\ge D_x|I,\mathit{D}))=\underset{{\Omega }_{\chi }}{\int }(P(d_x\ge D_x|I,\chi )-P(d_x\ge D_x{|I,D))}^{\mathrm{2}}f(\chi |D)\mathrm{d}\chi =\\ \underset{{\Omega }_{\chi }}{\int }P(d_x\ge D_x{|I,\chi )}^{\mathrm{2}}f(\chi |D)\mathrm{d}\chi -{\underset{{\Omega }_{\mathit{\chi }}}{\mathrm{E}}}^{\mathrm{2}}(P(d_x\ge D_x|I,\mathit{D},\mathit{\chi }))=\end{array}$$

$$\underset{{\Omega }_{\mathit{\chi }}}{\mathrm{E}}(P^{\mathrm{2}}(d_x\ge D_x|I,\mathit{D},\mathit{\chi }))-P^{\mathrm{2}}(d_x\ge D_x|I,\mathit{D})$$

（6）

对式（5）、（6）加以拓展，可以建立稳健的超越概率曲线。对式（2）中d的年超越概率λ（d）按曲线参数的概率分布f（χ|D）积分，可得到稳健的超越概率曲线与标准差：

$$\lambda (d_x|\mathit{D})=\underset{{\Omega }_{\chi }}{\int }(\int P(d_x\ge D_x|I,\mathit{\chi })\mathrm{d}\lambda (I))f(\chi |D)\mathrm{d}\chi =$$

$$\underset{{\Omega }_{\mathit{\chi }}}{\mathrm{E}}(\int P(d_x\ge D_x|I,\mathit{\chi })\mathrm{d}\lambda (I))$$

（7）

$${\sigma }^{\mathrm{2}}(\lambda (d_x|\mathit{D}))=\underset{{\Omega }_{\chi }}{\int }(\int P(d_x\ge D_x|I,\mathit{\chi })d\lambda (I)-\lambda (d_x{|D))}^{\mathrm{2}}f(\chi |D)\mathrm{d}\chi =$$

$$\underset{{\Omega }_{\mathit{\chi }}}{\mathrm{E}}({(\int P(d_x\ge D_x|I,\mathit{\chi })\mathrm{d}\lambda (I))}^{\mathrm{2}})-{\lambda }^{\mathrm{2}}(d_x|\mathit{D})$$

（8）

2.2　直接回归法进行稳健的易损性分析以及不足

Jalayer等^［10］还提出采用直接回归法进行稳健的易损性分析。直接回归法是云图法中针对未经缩放地震波模拟数据的处理方法，将数据视为独立的数据点，忽略数据内部的相关性，对全部数据采用简单回归的方式建立d_x和I的关系。Straub等^［31］认为，云图法中具有相关性（如震级等参数相似）的观测值不宜采用传统的简单回归法计算，因为可能会低估参数的不确定性。在传统增量动力分析法中，每条IDA曲线内样本的地震波波形完全一致，具有极强的相关性。根据上述理论，采用直接回归法处理增量动力分析数据将极大地低估参数的不确定性。

直接回归法通常假定d_x服从对数正态分布，在对数空间采用线性函数拟合d_x与I的经验关系：

$$\mathrm{l}\mathrm{n}{d}_x=a\mathrm{l}\mathrm{n}I+b+\epsilon$$

（9）

式中：a、b为控制线性函数形状的参数；ε为拟合残差。

图5展示了线性函数的拟合结果与函数的90%预测区间（90%PI）。可以看到，线性函数能较好地拟合IDA曲线，但图中$\mathrm{l}\mathrm{n}{d}_x$与$\mathrm{l}\mathrm{n}v$的关系依然呈现一种S型的非线性形态。因此，还选择了可以刻画S形非线性关系的Logistic函数：

$$\mathrm{l}\mathrm{n}{d}_x=\frac{c}{\mathrm{1}+{\mathrm{e}}^{-(a\mathrm{l}\mathrm{n}I-b)}}+h+\epsilon$$

（10）

式中：h为S型曲线的渐进下限；c为渐进上限到下限的距离。根据Youd等^［32］对于液化侧移预测的建议，侧移小于0.01 m视为未发生液化侧移，当预测侧移超过6.00 m时认为可能会发生液化流滑。本案例中，只有少量数据样本小于0.01 m或者大于6.00 m，因此选择0.01 m作为渐进下限，选择6.00 m作为渐进上限，即h=ln 0.01=-4.60，c=ln 6.00-ln 0.01=6.40。从图5可以直观看出，S型非线性函数的拟合效果相比线性函数更好。采用S型非线性函数拟合残差的标准差为0.72，小于线性函数的0.80。

图5  对数空间的IDA曲线与线性函数、S型非线性函数的拟合结果

Fig.5  IDA curves in the logarithmic space and fitted with linear function and non-linear sigmoid function

根据直接回归法建立的d_x与I关系，式（5）中易损性曲线P（d_x≥D_x|I，χ）可以计算如下：

$$P(d_x\ge D_x|I,\mathit{\chi })=\Phi \left(\frac{\mathrm{l}\mathrm{n}{d}_x-\mathrm{l}\mathrm{n}{D}_x}{{\sigma }_{\epsilon }}\right)$$

（11）

式中：σ_ε为残差标准差；易损性曲线参数χ=（a，b，σ_ε）。

图6展示了采用线性函数式（9）和S型非线性函数式（10）建立的稳健的易损性曲线。可以看到，当侧移较小时（见图6a、b），线性函数的易损性曲线倾向于高估概率；当侧移较大时（见图6c、d），线性函数倾向于低估概率，易损性曲线的中段落在了观测概率90%置信区间之外。在图6a、b、d中S型非线性函数的易损性曲线与观测概率符合程度较高，在图6c中S型非线性函数也低估了损伤概率，但相对线性函数低估的程度更小。上述现象说明，直接回归法建立的稳健的易损性曲线存在与观测概率不符的情况，采用S型非线性函数相对于线性函数具有更好的拟合效果。

图6  基于直接回归法采用线性函数和S型非线性函数建立的稳健的易损性曲线

Fig.6  Robust fragility curves based on simple regression method with linear function and non-linear sigmoid function

3 基于贝叶斯分层模型的稳健的易损性分析

3.1　方法与步骤

2节讨论了采用传统增量动力分析法和直接回归法进行稳健的易损性分析的不足，本节中提出了一种基于贝叶斯分层模型的稳健的易损性分析方法。Vamvatsikos等^［2］建议，可以采用统计分析每条IDA曲线拟合参数的方式标定IDA曲线的分布。采用贝叶斯分层模型分析每条IDA曲线的拟合参数，量化IDA曲线的组间和组内变异性，标定IDA曲线的分布并建立稳健的易损性曲线。

图7所示为贝叶斯分层模型框架，该模型也是一种贝叶斯网络。边代表变量之间的条件依赖关系，箭头代表依赖关系的方向。例如，由节点χ_hyper指向节点χ₁的箭头代表χ₁对χ_hyper的条件依赖性，表示χ₁服从χ_hyper作为参数的概率分布。模型主要分为标定与预测2个部分，图8中左半边为标定部分，基于各条IDA曲线数据标定拟合曲线参数χ_i、σ_εi和超参数χ_hyper、σ_hyper的分布。各条地震波的χ_i与σ_εi在数值上不完全相同，但服从χ_hyper、σ_hyper作为参数的概率分布，可以考虑IDA曲线的组间变异性。各条地震波具有独立的σ_εi，可以考虑IDA曲线的组内变异性。图8中右半边为预测部分，根据超参数的分布预测任意地震波产生的侧移，再根据侧移分布建立稳健的易损性曲线。

图7  贝叶斯分层模型框架

Fig.7  Framework of hiearchical Bayesian model

图8  侧移的后验样本分布

Fig.8  Distribution of posterior samples of lateral displacement

建立稳健的易损性曲线可分为以下6个步骤：

步骤1   选择合适的IDA曲线拟合函数F（I|χ），其中χ既是拟合曲线参数又是基于拟合函数的易损性曲线的参数。

步骤2   假设不同地震波的IDA曲线参数χ_i以及拟合误差σ_εi分别服从χ_hyper、σ_hyper作为参数的概率分布。根据贝叶斯定理，曲线参数χ_i、拟合误差σ_εi、超参数χ_hyper、σ_hyper的后验概率密度计算式为

$$f({\mathit{\chi }}_i,{\mathit{\chi }}_{\mathrm{h}\mathrm{y}\mathrm{p}\mathrm{e}\mathrm{r}},{\sigma }_{\epsilon i},{\mathit{\sigma }}_{\mathrm{h}\mathrm{y}\mathrm{p}\mathrm{e}\mathrm{r}}|\mathit{D})\propto f({\mathit{\chi }}_i,{\mathit{\chi }}_{\mathrm{h}\mathrm{y}\mathrm{p}\mathrm{e}\mathrm{r}},{\sigma }_{\epsilon i},{\mathit{\sigma }}_{\mathrm{h}\mathrm{y}\mathrm{p}\mathrm{e}\mathrm{r}})l(\mathit{D}|{\mathit{\chi }}_i,{\mathit{\chi }}_{\mathrm{h}\mathrm{y}\mathrm{p}\mathrm{e}\mathrm{r}},{\sigma }_{\epsilon i},{\mathit{\sigma }}_{\mathrm{h}\mathrm{y}\mathrm{p}\mathrm{e}\mathrm{r}})$$

（12）

式中：f（χ_i，χ_hyper，σ_εi，σ_hyper）为先验概率密度函数； l（D|χ_i，χ_hyper，σ_εi，σ_hyper）为似然函数。f（χ_i，χ_hyper，σ_εi，σ_hyper）可以写成如下形式：

$$f({\mathit{\chi }}_i,{\mathit{\chi }}_{\mathrm{h}\mathrm{y}\mathrm{p}\mathrm{e}\mathrm{r}},{\sigma }_{\epsilon i},{\mathit{\sigma }}_{\mathrm{h}\mathrm{y}\mathrm{p}\mathrm{e}\mathrm{r}})=$$

$$f({\mathit{\chi }}_{\mathrm{h}\mathrm{y}\mathrm{p}\mathrm{e}\mathrm{r}})f({\mathit{\sigma }}_{\mathrm{h}\mathrm{y}\mathrm{p}\mathrm{e}\mathrm{r}}){\displaystyle \prod _{i=\mathrm{1}}^k}f({\mathit{\chi }}_i|{\mathit{\chi }}_{\mathrm{h}\mathrm{y}\mathrm{p}\mathrm{e}\mathrm{r}}){\displaystyle \prod _{i=\mathrm{1}}^k}f({\sigma }_{\epsilon i}|{\mathit{\sigma }}_{\mathrm{h}\mathrm{y}\mathrm{p}\mathrm{e}\mathrm{r}})$$

（13）

l（D|χ_i，χ_hyper，σ_εi，σ_hyper）等于各条IDA曲线在给定χ_i和σ_εi的情况下观察到已知侧移的概率之积，可以写成如下形式：

$$l(\mathit{D}|{\mathit{\chi }}_i,{\mathit{\chi }}_{\mathrm{h}\mathrm{y}\mathrm{p}\mathrm{e}\mathrm{r}},{\sigma }_{\epsilon i},{\mathit{\sigma }}_{\mathrm{h}\mathrm{y}\mathrm{p}\mathrm{e}\mathrm{r}})={\displaystyle \prod _{i=\mathrm{1}}^k}{\displaystyle \prod _{j=\mathrm{1}}^{n_i}}f(d_{ij}|{\mathit{\chi }}_i,{\sigma }_{\epsilon i})$$

（14）

式中：f（d_ij|χ_i，σ_εi）为在给定χ_i和σ_εi的情况下观察到已知侧移d_ij的概率；n_i是第i条IDA曲线上数据点的数量。将IDA曲线数据代入式（12），即可得到参数χ_hyper、σ_hyper的后验概率分布f（χ_hyper，σ_hyper|D）。

步骤3   场地内任意潜在地震波对应的IDA曲线服从χ_hyper、σ_hyper作为参数的概率分布。从上述概率分布中，抽取曲线参数样本χ${}^{\mathrm{*}}$和拟合残差标准差样本σ${}_{\epsilon }^{\mathrm{*}}$，根据残差标准差样本σ${}_{\epsilon }^{\mathrm{*}}$抽取残差样本ε${}^{\mathrm{*}}$。

步骤4   将χ${}^{\mathrm{*}}$和ε${}^{\mathrm{*}}$代入式（10），生成不同地震动强度I下的土体侧移样本d${}_x^{\mathrm{*}}$。曲线参数及拟合误差的不确定性将传递到侧移样本d${}_x^{\mathrm{*}}$中，这些样本将代表场地内所有潜在侧移的分布。

步骤5   对于每个地震动强度I，计算侧移样本d${}_x^{\mathrm{*}}$中超过损伤限值D_x的样本数量n，再利用式（4）计算出当前地震动强度I下的损伤概率。由于抽取的d${}_x^{\mathrm{*}}$数量不受限，因此可以选择较大的样本数。根据威尔逊置信区间的计算方法，当选取的样本数量N足够大时，威尔逊置信区间将收敛，从而近似消除了计算概率的不确定性。将不同地震动强度I下的损伤概率连接起来，即可得到给定χ_hyper、σ_hyper情况下的易损性曲线P（d_x≥D_x|I，χ_hyper，σ_hyper）。

步骤6   由步骤2获得的f（χ_hyper，σ_hyper|D）是概率分布。对于χ_hyper、σ_hyper的不同取值，重复步骤3~5，将易损性曲线P（d_x≥D_x|I，χ_hyper，σ_hyper）与参数的后验分布 f（χ_hyper，σ_hyper|D）代入式（5）、（6），即可得到稳健的易损性曲线与易损性曲线的标准差。

3.2　案例场地的计算流程与结果

（1）选择合适的拟合函数。根据2.2节直接回归法的比较结果，S型非线性函数相对线性函数具有更好的拟合效果。为了检验S型非线性函数对单条IDA曲线的拟合效果，采用最小二乘法拟合各条IDA曲线，拟合参数如表2所示。可以看到，曲线参数a、b存在较大差异。此外，从σ_ε可以看出，S型非线性函数整体上较好地拟合了IDA曲线，编号为3、4、8、9、10、12、13、15、17、19的IDA曲线拟合误差非常小，但编号2、11存在较大的拟合误差。这说明：一方面，IDA曲线形状具有较大差异，即曲线形状参数存在较大的组间变异性；另一方面，拟合误差存在较大差异，即存在组间变异性。采用如式（10）所示的S型非线性函数建立考虑曲线形状参数和拟合误差组间变异性的贝叶斯分层模型。

表2  不同IDA曲线的拟合参数

Tab.2  Fitted function coefficients for different IDA curves

编号	a	b	σε	编号	a	b	σε
1	1.14	1.65	0.40	11	1.19	3.69	0.98
2	1.00	3.05	1.11	12	1.03	2.24	0.26
3	0.95	2.43	0.15	13	1.21	3.20	0.27
4	0.90	1.97	0.28	14	1.63	4.35	0.37
5	1.57	2.87	0.35	15	1.72	4.39	0.24
6	1.28	3.54	0.53	16	1.93	4.36	0.43
7	0.92	2.44	0.29	17	1.32	3.19	0.17
8	1.39	2.88	0.22	18	0.71	1.75	0.62
9	1.36	2.41	0.23	19	0.96	1.93	0.21
10	1.00	1.91	0.25	20	2.57	4.64	0.68

（2）标定曲线参数的后验分布。对于S型非线性函数式（10），待标定的参数为控制曲线形状的参数a、b和拟合误差σ_εi。经Kolmogorov-Smirnov（KS）检验^［33］，表2中各曲线参数a、b服从正态分布假设的P值分别为0.691、0.674，可以认为a、b服从正态分布。假设参数a_i、b_i服从以［μ_a，μ_b］为均值、Σ为协方差的二元正态分布。Σ可以写成如下形式：

$$\mathit{\Sigma }=\left[\begin{array}{cc}{\sigma }_a^{\mathrm{2}}& \rho {\sigma }_a{\sigma }_b\\ \rho {\sigma }_a{\sigma }_b& {\sigma }_b^{\mathrm{2}}\end{array}\right]$$

（15）

式中：σ_a、σ_b分别为a_i、b_i的标准差；ρ为相关系数。

拟合误差σ_ε恒为正值，且经过KS检验，表2中ln σ_ε服从正态分布假设的P值为0.471。可以认为，拟合误差σ_εi服从对数正态分布，即ln σ_εi服从以ξ_ε为均值、β_ε为标准差的正态分布。于是，易损性曲线参数χ_hyper=（μ_a，μ_b，σ_a，σ_b，ρ），σ_hyper=（ξ_ε，β_ε），对于先验分布，采用弱先验假定，假设μ_a、μ_b、ξ_ε均服从均值为0、标准差为1 000的正态分布，标准差σ_a、σ_b、β_ε服从形状系数为1的半柯西分布，相关系数ρ服从Lewandowski-Kurowicka-Joe（LKJ）分布^［34］，代入式（12）、（13）即可得模型的先验概率密度函数。后验概率密度函数式（14）采用基于马尔科夫链的蒙特卡罗模拟求解。

（3）建立稳健的易损性曲线。对于一组超参数χ_hyper=（μ_a，μ_b，Σ），σ_hyper=（ξ_ε，β_ε），场地内任意潜在的IDA曲线的参数将服从以μ_a、μ_b为均值、Σ为协方差的二元正态分布，拟合误差σ_εi也服从以ξ_ε为均值、β_ε为标准差的对数正态分布。根据μ_a、μ_b、Σ抽取曲线参数样本a${}^{\mathrm{*}}$、b${}^{\mathrm{*}}$，根据ξ_ε、β_ε抽取拟合误差样本σ${}_{\epsilon }^{\mathrm{*}}$，再根据σ${}_{\epsilon }^{\mathrm{*}}$抽取残差样本ε${}^{\mathrm{*}}$。将样本代入式（10），得到侧移的后验样本d${}_x^{\mathrm{*}}$，如图8所示。对图8中侧移的后验样本，计算每个侧移样本中达到损伤限值D_x的样本数n，采用式（4）计算达到损伤状态的概率，将不同地震动强度I下的损伤概率连接起来即可得到易损性曲线。对参数μ_a、μ_b、Σ、ξ_ε、β_ε后验分布中的其他取值重复上述步骤，并与参数的后验分布共同代入式（5）、（6）中，得到稳健的易损性曲线和易损性曲线的标准差。

4 验证与比较

4.1　侧移的后验样本分布与稳健的易损性曲线

图8还展示了贝叶斯分层模型预测的侧移后验样本90%分布区间，即贝叶斯分层模型的90%预测区间。可以看到，贝叶斯分层模型90%预测区间很好地覆盖了IDA曲线，260个增量动力分析样本点中仅有9个落在90%预测区间外，说明贝叶斯分层模型较好地模拟了IDA曲线的分布。直接回归法的90%预测区间也展示在图8中。由于直接回归法考虑恒定的拟合误差，因此ln v较小时90%预测区间相比贝叶斯分层模型更宽、ln v较大时预测区间更窄，这与图8中观察到IDA曲线在ln v较小时离散程度更低、ln v较大时离散程度更高的现象相矛盾。上述现象说明，贝叶斯分层模型相比直接回归法对IDA曲线分布的预测更加准确。

图9a~d分别展示了贝叶斯分层模型与直接回归法建立的临界侧移为0.10、0.20、0.50、0.90 m时的稳健的易损性曲线。可以看到，相比于直接回归法，贝叶斯分层模型建立的曲线在地表峰值速度较小时具有更高的损伤概率、在地表峰值速度较大时具有更低的损伤概率，中段斜率也显著低于直接回归法，说明贝叶斯分层模型考虑的不确定性更大。从图9c、d可以看出，直接回归法的易损性曲线在ln v=1.6、2.2时倾向于低估侧移的超越概率。从实际工程应用的角度，工程师更关心损伤概率较低时预测的准确度。贝叶斯分层模型在ln v较大时出现了低估损伤概率的现象，损伤概率接近于1.0。贝叶斯分层模型在低损伤概率时预测准确、在高损伤概率时预测不准确，这对实际工程决策影响较小。直接回归法在损伤概率较小时低估了失效风险，严重降低了工程决策的可靠性。

图9  基于直接回归法与贝叶斯分层模型建立的稳健的易损性曲线

Fig.9  Robust fragility curves based on simple regression method and hierarchical Bayesian model

4.2　易损性曲线的标准差与置信区间

将易损性曲线与曲线参数的后验样本代入式（6），得到易损性曲线标准差、易损性曲线的90%置信区间（90%CI），如图10所示。当ln v=2时，侧移限值为0.50 m的观测概率点落在易损性曲线90%置信区间外，且距置信区间边缘较远，观测概率的90%置信区间也不与易损性曲线90%置信区间相交；当ln v较大时，侧移限值为0.50 m和0.90 m的观测概率出现了部分不为1.0的情况，90%置信区间无法包络这些观测点。易损性曲线的标准差整体较小，当ln v较大时0.50 m和0.90 m的观测概率发生了多次异动，但标准差为零，明显低估了真实的不确定性。上述现象说明，直接回归法的易损性曲线与观测概率相符程度较低，其90%置信区间偏小，预测标准差也偏小，印证了前文直接回归法可能会低估参数不确定性的论述。

图10  易损性曲线90 %置信区间与标准差

Fig.10  90 % confidence interval (CI) and standard deviation of fragility curves

图10d~f中，贝叶斯分层模型的易损性曲线90%置信区间大致覆盖了所有观测概率点。当损伤概率小于0.7时，易损性曲线90%置信区间全部被观测概率的置信区间覆盖，说明当损伤概率较小时贝叶斯分层模型的预测较为准确，并且准确反映了易损性曲线的不确定性。虽然贝叶斯分层模型在损伤概率接近于1.0时低估了观测概率，但易损性曲线的置信区间依然与所有观测概率的置信区间相交。总体而言，贝叶斯分层模型建立的易损性曲线较为准确地反映了真实损伤概率的不确定性。

4.3　液化侧移稳健的超越概率曲线

将4.2节得到的易损性曲线与地表峰值速度的年超越概率代入式（7），即可得到场地液化侧移稳健的年超越概率曲线，如图11所示。可以看到，贝叶斯分层模型与直接回归法的超越概率曲线整体上较为接近，但发生了多次交错，直接回归法的超越概率曲线变化相对较小。这可能是由于直接回归法采用样本独立和相等拟合误差等简化假设，而贝叶斯分层模型则可以模拟超越概率曲线可能具有的更加丰富的阶段特征。

图11  液化侧移稳健的年超越概率曲线

Fig.11  Robust exceedance curve for liquefaction-induced lateral displacement

图12展示了2种方法建立的年超越概率曲线的标准差和90%置信区间。与稳健的易损性曲线相同，直接回归法建立的超越概率曲线具有较小的置信区间和标准差。然而，超越概率曲线没有观测数据作为比较依据，考虑到贝叶斯分层模型能较为准确地量化易损性曲线的不确定性，传导到超越概率曲线中的不确定性也应当较为准确。贝叶斯分层模型的超越概率曲线90 %置信区间虽然更大，但是属于实际工程适用的范围，具有较高的工程实用性。

图12  液化侧移年超越概率曲线的90 %置信区间与标准差

Fig.12  90 % confidence interval and standard deviation of exceedance curves for liquefaction-induced lateral displacement

5 结 语

本文提出了一种适用于增量动力分析法的稳健的易损性分析方法。该方法可以较好地量化各个环节的不确定性，将不确定性传递到易损性曲线中，有助于更好地实现基于性能的抗震设计。该方法采用贝叶斯分层模型量化不同IDA曲线间的差异，预测侧移的潜在分布，建立液化侧移稳健的易损性曲线和超越概率曲线。贝叶斯分层模型可以准确地模拟IDA曲线的分布，建立的稳健易损性曲线在失效概率小于0.7时与观测失效概率符合程度较高。

作者贡献声明

葛一荀：数值仿真，数据分析，论文起草与校对。

张 洁：论文总体设计，方法指导，论文起草与校对。

黄宏伟：论文总体设计，方法指导，项目管理。