考虑参数不确定性的无筋砌体结构地震易损性分析

蒋亦庞，苏 亮，黄 鑫

(浙江大学空间结构研究中心，浙江，杭州 310058)

在历次地震中，砌体结构的大量破坏与倒塌造成了非常严重的人员伤亡和财产损失[1-2]，而无筋砌体结构由于未设置合理的抗震构造措施，其地震破坏则尤为严重。由于无筋砌体结构仍是我国广泛存在的建筑结构形式之一[3]，对其进行准确的地震易损性分析具有非常重要的工程实际意义。

结构地震易损性分析需要考虑不确定性的影响。以往的研究认为地震易损性分析的不确定性主要源于地震动作用，其影响大小则通常是通过对多条地震动输入的增量动力分析(incremental dynamic analysis, IDA)进行计算考察[4-5]。然而，随着近年来不确定性研究的不断深入，越来越多的研究成果表明[6-10]，当结构参数不确定性较大或结构响应非线性程度较高时，结构参数的不确定性将会带来无法忽视的影响。

目前对结构参数不确定性的研究主要集中于钢筋混凝土结构、钢结构以及桥梁结构等[6-10]，国内外对砌体结构参数不确定性影响的研究则尚处于起步阶段。Parisi等[11]利用静力推覆法研究了8个材料参数不确定性对一2层无筋砌体结构的抗震能力影响，分析结果显示结构变形能力产生了较大的离散性，表明材料参数不确定性对砌体结构抗震能力的影响较大，但文章没有考虑地震动不确定性的影响。Lu等[12]基于IDA分析，采用简化的多自由度层剪切模型，利用蒙特卡罗法和一次二阶矩(first-order second-moment, FOSM)方法对加筋砌体结构的区域震害预测进行了结构参数敏感性研究。研究表明，FOSM方法能够在较少的计算量下，得到与蒙特卡罗法相似精度的结果；结构参数的不确定性对单个结构的地震易损性较大，而在区域震害分析时则可以被忽略。当然由于文章仅采用了简化的多自由度层剪切模型，其精确程度值得进一步计算验证。

鉴于无筋砌体结构参数的高离散性及其地震响应的强非线性[11]，在其地震易损性分析中考虑结构参数的不确定性显得尤为必要。正是基于这样的工程背景，本文以4幢不同层数的无筋砌体结构为研究对象，在OpenSees中建立其等效框架有限元模型，采用IDA法和FOSM法考察了地震动以及包括荷载、材料和阻尼在内的9个参数不确定性对结构地震易损性的影响，并给出了相应的工程建议。

1 无筋砌体结构模型

1.1 无筋砌体结构参数

采用层数为3层~6层的4幢典型无筋砌体结构为研究对象[13]，其中3层结构模型示意图如图1(a)所示。除层数不同之外，其余3幢结构的参数均与3层结构相同。不考虑结构平面外破坏及扭转效应，在此仅取一榀平面无筋砌体墙作简化分析。单层单跨墙的尺寸为层高3 m、跨长3.6 m，平面外结构跨度3.9 m。单层单跨墙的门窗洞口尺寸为宽1.2 m、高1.2 m，位于墙体中心。墙厚370 mm，墙材料取MU20的烧结普通砖和M10的砂浆。楼板假定为钢筋混凝土材料，厚120 mm，计算时仅考虑其荷载作用。

1.2 结构建模方法

结构模型采用等效框架法建模，在OpenSees上进行有限元分析。等效框架法是目前被广泛采用的砌体结构建模方法，用于砌体结构的整体抗震分析，兼具高计算效率和高计算精度的优点[14-15]。等效框架法的基本思想是将墙肢和窗群梁等效成框架柱和框架梁单元，墙肢与窗群梁相交的部分由于在地震中一般不会破坏则等效成刚域，如图1(b)所示。

图1 3层无筋砌体结构及其等效框架模型Fig.1 3-storey unreinforced masonry structure and its equivalent frame model

等效框架模型在OpenSees中实现如下[14]：

1) 单元：墙肢和窗群梁选用基于柔度法的梁柱单元建模，刚域处采用足够大的弹性模量本构来模拟；

2) 等效高度：墙肢和窗群梁的等效高度选用Dolce[16]提出的准则，即外侧墙肢的等效高度为相邻洞口两顶点延伸出来的30°倾斜线所包围的高度，而内侧墙肢和所有窗群梁的等效高度与相邻洞口尺寸相同；

3) 截面：墙肢和窗群梁截面采用纤维模型，各个截面的单个纤维尺寸约为40 mm×40 mm，沿截面长、宽方向分别等分；

4) 单元的轴向与弯曲行为：截面纤维的单轴本构模型采用图2(a)所示的Kent-Park模型[17]，以模拟单元的轴向及弯曲行为。图中极限抗压强度fu和极限压应变εu取值分别为0.88倍抗压强度fm及1.6倍抗压强度处压应变εm[18]；

5) 单元剪切强度：单元剪切强度Vu取代表对角破坏的剪切强度Vu1和代表剪切滑移破坏的剪切强度Vu2的较小值。Vu1参考Turnšek等[19]提出的计算公式：

式中：ft为砌体材料抗拉强度；σ0为截面平均压应力(σ0=N/lt)，l和t分别为单元截面的宽度和厚度，系数b由单元的长宽比确定，如下式：

式中，h为单元长度。Vu2的计算公式如下[20]：

式中，fv0为砌体抗剪强度。

图2 模型所使用的两种本构模型Fig.2 Two constitutive laws used in model

6) 单元剪切行为：选用如图2(b)所示的Pinching4单轴本构模型[21]来描述单元的剪切行为，并以截面组合的形式加到纤维截面上，组合后单元响应为弯曲与剪切本构的共同作用结果。模型骨架曲线可由点1~4的坐标位置(γi，Vi)确定，骨架曲线正向与反向对称，其中γi与Vi分别代表第i点对应的位移与荷载，其计算详见式(4)；模型卸载与再加载曲线则由a、b、c、d、e与f的坐标位置确定，其中b点坐标(γm,Vm)代表历史加载的最大位移γm及对应的荷载Vm；a点坐标为(0.5γm, 0.25Vm)；c点横坐标由卸载刚度确定，其纵坐标为-0.05V3。d、e与f点分别与a、b、c点关于原点对称。其余控制滞回反应中逐渐增加的卸载刚度退化、再加载刚度退化以及强度退化的16个参数均使用Lowes等[22]的建议值。

式中：G为砌体材料剪切模量；AS为单元截面积。

1.3 结构不确定性参数

选取9个无筋砌体结构不确定性参数，分别为：砌体材料密度ρm、钢筋混凝土板材料重度γ、活荷载q、砌体抗压强度fm、砌体抗拉强度ft、砌体抗剪强度fv0、砌体弹性模量E、砌体剪切模量G、阻尼比ξ。由于目前针对砌体结构不确定性参数之间的相关性研究仍不充分，现有文献不足以获得各参数间的相关系数，因此本文作简化分析，假设各参数之间相互独立。同时，本文也不考虑模型参数在结构空间分布上的不确定性。各参数详细统计特征见表1。

表1 结构不确定性参数详细统计特征Table 1 Detailed statistical characteristics of structural uncertainty parameters

2 地震动参数及结构极限状态

2.1 地震动的选择

既有研究表明，在IDA分析中选取20条地震动记录足以考虑地震动的不确定性影响[4,9]。考虑到IDA分析比较耗时，本文从文献[24]建议的100条远场地震动记录中选取20条作为地震动输入考虑地震动的不确定性。所选地震动均来自美国太平洋地震研究中心PEER的强震数据库[26]，如表2所示。所选地震波的震级分布在5.65级~7.6级，峰值加速度(peak ground acceleration, PGA)分布在0.05g~0.81g。同时，来自同一场地震事件的地震动记录不超过两条。

2.2 地震动强度指标的选取

地震易损性分析中常用的地震动强度参数(intensity measure, IM)包括地面峰值加速度PGA和结构基本周期对应的加速度谱值Sa(T1, 5%)，T1为结构基本周期。由于结构层数和结构不确定性参数的改变都会使得结构的基本周期成为一个不确定值，如采用Sa(T1, 5%)作为IM参数会使分析变得更为复杂，故选取PGA作为IM参数。在进行IDA分析时，分别将20条地震动记录的PGA分别调幅为0.05g~0.8g(间隔0.05g)，对结构进行弹塑性时程分析。

表2 20条天然地震动记录Table 2 20 natural ground motion records

2.3 结构极限状态的定义

选取最大层间位移角为结构抗震性能指标，将无筋砌体结构的极限状态(limit state, LS)分为三种，即轻微破坏(LS1)、中等破坏(LS2)、严重破坏(LS3)，各极限状态的层间位移角限值如表3所示[27]。

表3 层间位移角限值Table 3 Allowance for story drift angle

3 考虑结构不确定性的地震易损性分析

本文采用FOSM方法考虑结构的不确定性，首先对FOSM方法作简单介绍。假设随机变量y是随机向量的函数y=f(X)，FOSM方法将该函数在随机向量X的均值μX处近似展开为一阶泰勒级数形式，从而估计随机变量y的均值μy和方差分别如式(5)及式(6)所示：

式中：ρij为xi和xj的相关系数；σxi为xi的标准差。由下式计算：

上述公式详细的推导过程可见文献[12]及其他相关文献。此外，值得一提的是，上述方法实际上是平均值一次二阶矩法(mean value first-order second-moment, MVFOSM)，属于FOSM方法的一种。然而在很多文献中通常将该方法直接称为FOSM方法[12,28]，为简化起见，本文也直接将其称为FOSM方法。

在本文中，随机向量X由9个结构不确定性参数确定，如下式所示：

结构的地震易损性[24]定义为给定地震强度水平时结构达到或超过某种极限状态的条件概率：

式中：D≥C表示结构达到或超过某种极限状态，其中D为地震需求，C为抗震能力；FR(x)为称为地震易损性函数，一般符合对数正态分布。采用基于地震动强度(IM-based)的地震易损性函数来展开研究，针对某一特定极限状态LSi，其表达式为：

式中：Φ[·]是标准正态概率分布函数；Ci为极限状态LSi时结构抗震能力限值；mR和βR分别为结构地震易损性函数的中位值和对数标准差。

针对某一特定极限状态LSi的结构地震易损性曲线，可由下述步骤获得[29]：

1) 选择NRTR条地震动记录对结构进行IDA分析；

2) 统计当地震动强度IM=im时，NRTR条地震动记录中导致结构超过该极限状态的地震动数量Ni，采用基于频率的统计方法获得结构超过该极限状态的失效概率为：

3) 将获得的失效概率用式(10)拟合，可获得结构对应该极限状态的易损性曲线，如图3所示。

将随机变量y定义为结构达到某一极限状态的平均地震动强度mR的对数，即：

图3 地震易损性曲线的拟合Fig.3 Fitting of a fragility curve

由文献[12, 30]可知，当只考虑地震动的不确定性时，式(10)的对数标准差为βR=βRTR，表征地震动不确定性的影响；当同时考虑地震动和结构参数的不确定性时，由结构参数不确定性贡献的对数标准差βMDL可由FOSM方法通过式(6)获得，此时定义新易损性函数的总对数标准差βTOT可由式(13)获得[30]：

由上述原理可知，当进行考虑结构参数不确定性的地震易损性分析时，若结构数量为n个，结构不确定性参数为m个，则需要进行的IDA分析次数为：NIDA=n×(2×m+1)。

4 结果分析

以表2中的20条地震动记录作为输入，按照第3节所介绍的计算方法，对层数为3层~6层的4幢无筋砌体结构进行4×(2×9+1)=76次IDA分析和地震易损性分析。其中IDA分析在主频为2.8 GHz的Intel i9-7960X处理器配置电脑中进行，共耗时约180 h。

图4为4幢无筋砌体结构在不同极限状态下的地震易损性曲线，其中实线表示只考虑地震动不确定性时的计算结果，虚线为同时考虑地震动和结构参数不确定性时的计算结果。为考察地震动和结构参数不确定性对无筋砌体结构地震易损性的影响程度，表4进一步列出了βRTR、βMDL和βTOT三个分别考虑不同不确定性参数组合的易损性曲线对数标准差值。

由地震易损性曲线特性可知，曲线的倾斜程度表征了不确定性对地震易损性的影响程度。以图4(c)中的6号易损性曲线为例，5层无筋砌体结构在PGA处于0.35g~0.7g时，其达到LS3极限状态的概率为20%~80%，这表明不确定性导致了结构地震响应存在很大离散性，从而使得结构达到LS3极限状态所需要的地震强度也变得较为离散。因此，图4的结果表明不确定性对无筋砌体结构的地震易损性影响是不可忽视的。而这种影响在结构极限状态LS3时表现的更为明显，体现在图4中即LS3对应的易损性曲线倾斜度明显大于LS1和LS2。究其原因，是因为结构极限状态LS3所对应的地震动输入增强，使得结构的非线性程度增大，放大了不确定性对无筋砌体结构地震易损性的影响。与图4的分析结果相呼应，表4中每幢结构极限状态LS3所对应的βTOT值均明显大于极限状态LS1和LS2所对应的βTOT值。

图4 考虑和不考虑结构参数不确定性的易损性曲线Fig.4 Fragility curves with and without structural parameter uncertainty

表4 不确定性对结构地震易损性的影响Table 4 Effect of uncertainty on the structural seismic fragility

由前文可知，βMDL与βRTR值的大小分别表征了结构参数及地震动不确定性对地震易损性的影响大小，其值越大表示影响越大。因此，计算了表4中βMDL与βRTR的比值并列于其中以对比结构参数以及地震动不确定性对无筋砌体结构地震易损性的影响大小，其最小值为6层LS3时的0.466，最大值为3层LS1时的1.191。同时，将表4中的βMDL与βRTR值图形化作成柱状图进行直观比较，如图5所示，左斜线柱状图表示βRTR值，右斜线柱状图表示βMDL值。可见，结构参数的不确定性相比于地震动的不确定性是不可忽略的，有时甚至会超过地震动不确定性的影响，如3层结构对应于极限状态LS1、LS2以及4层结构对应于极限状态LS1的易损性分析结果中，βMDL值均比与βRTR值大。从图4中也可以看出，考虑与不考虑结构参数不确定性的两条易损性曲线并不重合，考虑结构参数不确定性后易损性曲线倾斜度有明显增加。

图5 地震易损性曲线对数标准差对比Fig.5 Comparison of logarithmic standard deviation of seismic fragility curves

同时，从表4及图5可看出，βRTR值随着结构层数的增加总体上呈现了明显增长的趋势，例如极限状态LS2对应的βRTR值从3层的0.158增长到了6层的0.277。说明地震动不确定性的影响在大部分情况下会随着结构层数的增加而增大。而同时也能发现，βMDL值并没有随着结构层数增加而增加，反而有下降趋势。除了3层LS3的0.237外，其余的βMDL值均分布在0.138~0.191内，总体来看变化不大。这说明结构参数不确定性对不同层数无筋砌体结构地震易损性的影响基本相同。而这就导致，随着结构层数的增加，结构参数不确定性相比于地震动不确定性的影响呈现减小趋势，体现在图4中即随着层数的增加，考虑与不考虑结构参数不确定性的两条地震易损性曲线越来越趋近于重合。由此可见，结构参数不确定性在较低层无筋砌体结构地震易损性中的相对影响更大，在计算时尤其需要得到重视。

为了进一步研究无筋砌体结构地震易损性分析中9个结构不确定性参数的敏感性大小，以3层无筋砌体结构为例，在前文计算的基础上，再对每个参数分别取增加和减少2倍标准差进行地震易损性分析。此时，每个参数共有5种计算工况：以易损性曲线的地震强度中位值mR为评价指标，对每个变量的5种计算工况结果进行归一化处理：

式中：mR(·)为参数xi取相应值而其他参数取平均值时的地震强度中位值；(·)为对应的归一化结果。(·)偏离1越多，表示该参数敏感性越高。

结果如图6所示。从图6可见，除阻尼比外的8个参数变化时计算的易损性曲线地震强度中位值归一化结果在0.85~1.1之间变化，变化范围为0.25；而阻尼比变化时的归一化结果最小值与最大值分别达到0.4和1.4，变化范围为1，是前者的4倍。表明阻尼比的敏感性远高于其他8个参数，是结构地震易损性分析结果的控制性参数。在无筋砌体结构地震易损性分析当中，阻尼比的准确取值尤其需要被重视。同时，中位值mR随着阻尼比的增大而增大。这是由于随着阻尼比的增大，结构耗能增加，结构响应随之减小，从而导致结构易损性曲线的中位值不断增大。

图6 地震强度中位值随各结构参数的变化规律Fig.6 Variation of the median earthquake intensity with different structural parameters

5 结论

本文以4幢典型无筋砌体结构为例，采用IDA和FOSM方法同时考虑地震动及结构参数不确定性对其地震易损性的影响，得到以下结论：

(1) 在无筋砌体结构的地震易损性分析中需考虑地震动和结构参数的不确定性影响。对于同一结构，地震破坏程度越高这种影响程度越大；

(2) 结构参数的不确定性相比于地震动的不确定性是不可忽略的，有时甚至会超过地震动不确定性的影响，且结构的层数越小，结构参数不确定性相比于地震动不确定性的影响越大，因此有必要在无筋砌体结构地震易损性分析中同时考虑这两种不确定性；

(3) 以3层结构为例，对无筋砌体结构地震易损性分析中9个结构不确定性参数进行了敏感性分析。结果表明，阻尼比变化对易损性曲线地震强度中位值的影响可达到其他单个参数的4倍，其敏感性最高，同时随着阻尼比的增大，易损性曲线的地震强度中位值也随之增大。