基于地震动参数的房屋震害矩阵曲线化分析方法*

陈 波，温增平，赵文哲

(中国地震局地球物理研究所，北京 100081)

0 引言

我国是世界上陆地地震灾害最严重的国家之一，特别是随着经济建设和城市化进程的不断加快，地震灾害潜在风险水平不断提高。易损性分析作为灾害风险分析的3个主要部分(地震危险性分析、承灾体的易损性分析、灾情损失评估)之一，其核心是建立建筑结构出现或超越各种破坏状态的概率与地震强度之间的关系。开展建筑结构地震易损性分析对于地震破坏和损失预测以及减轻地震灾害具有十分重要的意义。

在实际地震现场和震害预测的工作中，研究人员积累了大量基于烈度的建筑结构震害经验关系，即用烈度-震害程度的比例或概率构成的易损性矩阵，然而烈度的获得是基于人的感觉、器物反应、房屋为主的工程结构破坏、地表破坏等4类宏观现象来评定(胡聿贤，2006；中国地震烈度表，GB/T 17742—2008；明小娜等，2017)。这也意味着随着我国建筑抗震设防水平的提高，不同时期烈度本身的度量尺度也会发生变化，所以基于烈度的易损性矩阵具有较大的不确定性(Lagomarsino，Cattari，2014)。然而以往大多数震害预测中关于房屋建筑易损性的评估，都是基于烈度进行，即表示结构在不同烈度下出现各破坏状态的概率(韩淼等，2009)。为了准确预测建筑结构的抗震性能，更好地进行结构抗震设计，完成对现有建筑物的加固和维修，加快地震保险在全国范围内的推行，基于烈度的易损性矩阵很难满足实际需要，特别是对于土木结构、砖石结构等简易房屋，难以通过结构模拟分析获取其易损性曲线，如何完善现有的易损性矩阵，并将基于烈度的易损性矩阵转化为更合理的基于地震动参数的易损性矩阵显得尤为重要(胡少卿等，2007；刘如山等，2009；李静等，2012；马玉宏等，2015)。

本文采用最大似然估计的曲线方法对震害资料进行分析。选取PGA作为地震动参数，根据烈度与地震动参数之间的对应关系，将结构基于烈度的不同破坏状态的概率矩阵转化为基于PGA的不同破坏状态的概率矩阵，得到连续的基于地震动参数的易损性曲线。

1 易损性矩阵曲线化方法

地震易损性是指结构在不同地震作用下，出现各种破坏状态的条件概率，它从概率的角度定量表征结构的抗震能力，表达式为：

F(IM=x) =P(IM≥IMC/IM=x) =P(IMC≤x)

(1)

式中：IM为地震动强度指标；IMC是结构物达到破坏状态时对应的地震动参数；F(IM)为易损性函数。

建筑结构易损性主要通过地震易损性曲线和震害矩阵来表示。一般易损性曲线横坐标为地震的强度指标，如烈度、峰值加速度、峰值速度、地震动反应谱值等，纵坐标为达到不同破坏等级的超越概率。易损性曲线的获得可以通过经验法(Colombietal，2008)、解析法(Lagomarsino，Cattari，2014)及混合法(Aldemiretal，2013)。在易损性的研究进程中对易损性函数的选用，不同学者提出了不同的数学模型，如：逻辑回归函数(Singhal，Kiremidjian，1998)、威尔逊分布函数(Shinozukaetal，2007)、对数正态分布函数(Shinozukaetal，2000；Yamazakietal，2000)等，其中对数正态分布最为常用。

易损性曲线的拟合方法有很多种，如多项式拟合、最小二乘法拟合、最大似然估计法拟合等(Shinozukaetal，2000；Baker，2013；陈力波等，2012；吴子燕等，2014；林庆利等，2017)。本文采用最大似然估计法，通过MATLAB程序做具体计算和绘制易损性曲线，获取基于地震动参数的易损性曲线的具体流程，如图1所示。

图1 震害矩阵的曲线化计算方法

1.1 建筑物破坏等级分类

建筑物破坏等级通常可分为5类(GB/T 24335—2009)：DSk(k=0，1，…，4)分别表示基本完好、轻微破坏、中等破坏、严重破坏、毁坏5种破坏状态。LSs(s=1，2，3，4)为5种破坏状态之间的阈值，分别代表“轻微破坏”“中等破坏”“严重破坏”“毁坏”4种破坏状态的极限值，如图 2所示。

图2 易损性超越概率曲线不同破坏状态示意图

1.2 烈度与地震动参数对应关系

地震烈度是地震引起的地面震动及其影响的强弱程度。鉴于早期监测条件的限制，详细的强震记录很难获得，但是通过烈度评定，地震工作者积累了丰富的宏观烈度资料。胡聿贤和张敏政(1984)提出了缺乏强震观测资料地区地震动参数的估算方法；田启文等(1986)进行了根据烈度资料估算我国地震动参数衰减规律的一系列研究；特别是《中国地震烈度表》(GB/T 17742—2008)作为我国地震烈度的评定标准，给出了烈度与峰值加速度(PGA)和峰值速度(PGV)之间的对应关系。本文采用《中国地震烈度表》(GB/T 17742—2008)给出的烈度与地震动参数之间的对应值，从表1中可以看出，每一个烈度对应的PGA和PGV。

表1 地震烈度与峰值加速度PGA和峰值速度PGV的对应关系

1.3 基于最大似然估计法确定易损性双参数

在易损性研究的进程中，对数正态分布数学模型作为易损性的概率模型已经被验证是一种成熟、理想的模型，本文也采用对数累计正态分布函数作为易损性函数。易损性函数定义为：

(2)

式中：P(C|IM=x)是地震动强度达到IM=x时，结构物到达某一极限状态的概率；Φ是标准累积正态分布函数；θ是易损性函数的中位数，β为lnIM的标准差。

对于超越破坏状态的结构，由式(2)定义的易损性函数，可以得到在IM=xi处的似然函数的密度函数：

(3)

相应的，对于未超越破坏状态的结构，在IM=xi处的似然函数的密度函数为：

(4)

根据不同的破坏状态有不同的中位数θ和对数标准差β，定义似然函数为：

(5)

式中：N为统计结构总数；xi表示第i栋结构对应的峰值加速度；ti为伯努利事件，如果达到破坏状态，ti=1，否则，ti=0。

双参数θ和β计算公式为：

(6)

2 震害矩阵曲线化实例

2.1 震害资料的收集和统计

地震现场作为大的原型实验场，工程结构的真实震害经验是地震工程学发展的重要支柱，人类如何抗御地震灾害的认识很大程度也来源于此(温增平等，2009)。四川省地质构造复杂，地震活动频繁，房屋结构类型多样，历史上发生过多次7级以上特大地震，特别是2008年汶川8.0级地震、2013年芦山7.0级地震造成了严重的经济损失和人员伤亡。孙柏涛等(2014)采用不同分级和分区域的方法研究了南北地震带各类房屋建筑的分布特征和抗震能力，以四川省为例，将其分为：I区，西部地区；II区，中北部地区；III区，中南部地区。结合6个地级市204个调查点的23 964栋房屋建筑调查结果和汶川地震等历史震害经验对其房屋建筑进行调查，给出了分区域各类结构房屋的基于烈度的易损性矩阵。

2.2 土木结构房屋易损性曲线

本文首先以孙柏涛等(2014)统计给出的III区土木结构易损性矩阵为例，计算结构超越概率易损性曲线。图 3为分别以散点图和直方图来直观地表达III区土木结构房屋在不同烈度作用下的破坏比结果。

通过基于烈度的土木结构房屋易损性矩阵，从破坏状态DSi的破坏概率累积计算出超越每个极限破坏状态LSi的累计概率，得到相应土木结构房屋的超越概率矩阵，详见表2。

通过最大似然估计法，结合表 1，利用MATLAB软件即可求得土木结构房屋超越概率，计算得到其对数正态分布的易损性函数双参数θ和β，如表3所示。

表2 III区土木结构房屋超越概率矩阵(%)

图3 III区土木结构房屋震害矩阵散点图(a)及直方图(b)

易损性曲线双参数轻微破坏LS1中等破坏LS2严重破坏LS3毁 坏LS4θ0.0730.1370.2630.482β0.8390.7060.7790.711

采用易损性曲线双参数θ和β即可绘制出III区基于地震动参数PGA土木结构房屋的超越概率易损性曲线，如图4所示。超越概率的离散点是通过实际震害资料调查计算得到的，与最大似然估计法拟合的易损性超越概率校核曲线一致，可以得到理论分析的最大似然估计曲线是合理的、准确的。由此通过拟合的易损性曲线可以得到不同PGA(单位：g，1 g=9.8 m/s2)对应的不同破坏状态的超越概率(表4)。

图4 土木结构房屋易损性超越概率曲线

PGA/g轻微破坏LS1中等破坏LS2严重破坏LS3毁 坏LS40.01 0.89 0.01 0.00 0.000.10 64.66 32.78 10.72 1.350.20 88.55 70.40 36.28 10.810.30 95.41 86.66 56.74 25.260.40 97.88 93.55 70.51 39.670.50 98.91 96.67 79.55 52.060.60 99.40 98.18 85.54 62.100.70 99.65 98.96 89.58 70.010.80 99.79 99.38 92.36 76.190.90 99.86 99.62 94.30 81.011.00 99.91 99.76 95.69 84.76

为了得到处于不同破坏状态的概率，用PDSk表示在DSk状态下的概率；PLSs表示极限破坏状态的概率。结构不同破坏状态下的损害概率，通过以下公式求得：

PDS0(xi)=1-PLS1(xi)

(7)

PDSk(xi)=PLSk(xi)-PLSk+1(xi)=

(8)

表5 基于PGA的土木结构房屋易损性矩阵(%)

式中：k=1,2,3。

PDS4(xi)=PLS4(xi)

(9)

根据土木结构易损性超越概率曲线结果，结合表4，通过式(7)～(9)可以得到不同破坏状态下的易损性矩阵，如表 5所示。

通过式(7)～(9)绘制基于PGA土木结构房屋的易损性曲线，如图5所示，可直观看到不同破坏状态的损害概率随PGA增大的变化趋势，得到给定PGA下结构超越不同破坏状态的概率。

图5 基于PGA土木结构房屋震害矩阵易损性拟合曲线

2.3 多种房屋建筑的易损性曲线

与III区土木结构房屋易损性曲线计算一致，通过最大似然估计方法，利用烈度与PGA之间的对应关系，本文将提供不同结构处于不同破坏状态下的易损性函数的双参数θ和β，并绘制出其易损性曲线，如表 6所示。基于易损性曲线双参数绘制四川省3个分区不同类型结构的易损性曲线，如图 6所示。

按照烈度与PGV之间的对应关系，同样可以给出了基于PGV的不同结构房屋在不同破坏状态下的θ和β，如表7所示。

3 结论与讨论

本文提出一种房屋震害矩阵双参数曲线化方法，该方法基于烈度与地震动参数之间的对应关系，采用最大似然估计方法，将既有房屋易损性矩阵或实际震害统计的破坏比结果转化为双参数易损性曲线。具体认识如下：

(1)该方法弥补了基于烈度易损性曲线的不足，特别是针对难以进行模拟分析的结构类型，如土木结构、砖石结构等简易房屋，提供了更为简便的途径来得到相应的易损性曲线。

(2)以四川省房屋震害矩阵为例，基于中国地震烈度表中烈度与地震动参数的关系，给出了各类房屋结构的地震易损性曲线特征参数，分析结果可为四川省基于地震动参数的房屋建筑震害评估提供基础数据。需要注意的是，在实际工程应用中，针对不同的区域特征，应合理选择烈度-地震动参数(PGA，PGV)的关系模型。

图6 四川省分区域各类型房屋结构震害矩阵拟合曲线

分区及结构类型轻微破坏LS1中等破坏LS2严重破坏LS3毁 坏LS4θβθβθβθβ土木结构0.0700.6960.1260.6750.2390.7320.4670.694西部地区砖木结构0.0830.7030.1710.7390.3460.6780.7840.720(I区)穿斗木结构0.1020.6690.2120.7260.5190.8481.7300.926砖石结构0.0920.6890.1530.5470.2380.6130.4910.845土木结构0.0730.8390.1370.7060.2630.7790.4820.711中北部地区(II区)砖木结构0.0900.7320.1850.7370.3650.6870.8080.707穿斗木结构0.1090.6660.2240.6960.5580.8541.7500.851土木结构0.0730.8390.1370.7060.2630.7790.4820.711中南部及东部地区(III区)砖木结构0.0910.7230.1900.7420.3710.6880.8250.694穿斗木结构0.1050.6730.2130.7010.5420.8681.6960.853高层建筑(高设防)0.2980.5990.7370.6241.3650.4761.1580.077高层建筑(一般设防)0.2770.5310.6340.5951.2760.5141.1340.075钢混结构(高设防)0.2670.7850.5400.5480.8410.5061.6290.558钢混建筑(一般设防)0.2280.7380.4810.5840.6960.4661.5270.607设防砖混(高设防)0.1760.8610.3570.7340.6020.6001.6240.740设防砖混(一般设防)0.1390.8450.2920.7090.5100.6081.3720.828未设防砖混0.1190.7700.2230.6230.3960.6601.0310.786

表7 基于PGV的不同结构不同破坏状态下的易损性曲线双参数θ(单位：m/s)和β

(3)本文提出的易损性矩阵曲线化方法具有较高的可靠性和适用性，在进一步研究工作中可应用该方法对已有的易损性矩阵或曲线进行校验和完善，以便更好地服务于地震灾害快速评估、风险评价和保险估算等。

考虑到目前直接获取的强震记录均为未修正的结果，烈度与地震动参数对应关系还需要进一步完善，本文通过基于房屋震害矩阵曲线化方法得到的易损性曲线与真实情况相比存在一定误差。随着烈度与地震动参数对应关系的完善，本文方法可得到更加广泛的应用。