本文选取1950—2022年中国大陆地区发生的356次5级以上历史震例(表1),数据主要来源于Li等(2021)编制的1950—2018年中国大陆综合破坏性地震目录(MCCDE-CAT)、中国地震台网中心以及“中国大陆地震灾害损失评述”系列文章(郑通彦,郑毅,2012,2014,2015; 郑通彦等,2015; 陈通,郑通彦,2016; 文鑫涛,郑通彦等,2018; 文鑫清等,2021)。其中人口数据来自MCCDE-CAT和美国橡树岭国家实验室的人口栅格数据https://landscan.ornl.gov.。

表1 1950—2022年中国大陆地区部分5级以上历史震例基本信息
Tab.1 Details of some historical earthquake cases in Chinese mainland during 1950-2022

本文对所选数据进行筛查,最终选取中国大陆356次历史震例的基本信息,包括极震区烈度、发震时间、人口密度、死亡人数以及受灾人口(表1)。一般极震区烈度越大,造成灾害规模及建筑破坏等级越高,相应的人员死亡率越高。发震时间对人员在室率有直接影响(周中红等,2021; 谷国梁等,2021),在夜间,人员熟睡时意识和应急反应能力相对较低,无法采取紧急避险措施,更易造成人员死亡。人口密度能够反映人口分布状况,相同破坏程度的地震发生在人口密度较高的地区往往会造成更多人员死亡。极震区烈度、发震时间、人口密度这3个因素对地震人员死亡有直接的影响,因此,本文将其作为评估结果的指标。

由于烈度为Ⅰ～Ⅴ度时造成人员死亡的可能性较小,根据《中国地震烈度表》(GB/T 17742—2008)中宏观调查评定地震烈度,本文只考虑极震区烈度为Ⅵ度、Ⅶ度、Ⅷ度、Ⅸ度及以上4个级别。根据中国大陆地区作息习惯将发震时间分为白天06:00—19:59和夜间20:00—05:59。最后在极震区烈度及发震时间的分级基础上,考虑人口密度对地震人员死亡的影响,将人口密度分为≤100人/km²和>100人/km²。由于以下3种情况案例较少,本文暂不考虑:①极震区烈度为Ⅷ度、发震时间为晚上且人口密度>100人/km²; ②极震区烈度为Ⅸ度及以上、发震时间为白天且人口密度>100人/km²; ③极震区烈度为Ⅸ度及以上、发震时间为晚上且人口密度>100人/km²。按照上述条件在考虑极震区烈度和发震时间的基础上,考虑人口密度的分级,进而构建基于极震区烈度、发震时间及人口密度的三级场景分析框架,见表2。

表2 场景分析框架
Tab.2 Framework of the scenario analysis

将历史震例按照上述场景进行划分后,发现多数场景下的震例为小样本,因此采用Bootstrap抽样方法对各场景进行扩充。Bootstrap方法由Efron(1979)提出,是一种再抽样统计方法,主要用于解决小样本的评估问题,其基本思想是:在原始数据的范围内进行有效放回再抽样,样本含量仍为n,每个观察单位每次被抽到的概率相等,均为1/n,所得样本称为Bootstrap样本,如此重复X次,就可以得到X个Bootstrap样本,针对该X个样本逐一计算所关心的统计量,得到该统计量的N个观测值。通过对这些观测值进行直方图分析,可得到Bootstrap统计量的经验分布,基于该经验分布可进行统计推断及分析。

以场景S1为例,将场景S1中的死亡率通过Bootstrap抽样扩充为1 000个数据集,每个数据集中包含100个数据,再将死亡率均值x作为随机变量,拟合1 000个数据集死亡率x的频率直方图(图1),经过Bootstrap抽样扩充后各个场景死亡率x的均值、方差见表3。

图1 场景S1经过Bootstrap抽样扩充后各数据集均值x的频率直方图
Fig.1 Scenario 1 expands the frequency histogram of the mean x of each data set by Bootstrap sampling

为验证Bootstrap抽样扩充后死亡率数据的可用性和可靠性,将Bootstrap抽样扩充后死亡率的均值、方差(表3)与原始死亡率数据(表4)进行对比。结果表明,Bootstrap抽样扩充后的死亡率均值与原始死亡率均值特征基本一致,表明Bootstrap抽样方法能够保持样本统计量的可靠性和精确性,不会引入明显的偏差。同时死亡率方差相较于原始数据普遍偏小,该特征说明Bootstrap抽样方法有效地减少了数据的不确定性,提高了数据的稳定性和可靠性。因此Bootstrap抽样方法扩充后的虚拟数据为后续的分析和预测工作提供可靠的基础。

表3 场景S1～S13经Bootstrap抽样扩充后死亡率的均值及方差
Tab.3 Mean and variance of death rate of Scenario 1-13 expanded by Bootstrap sampling

表4 场景S1～S13原始死亡率数据均值及方差
Tab.4 Mean and variance of the raw death rate data of Scenario 1-13

场景S1的死亡率均值服从 [FK(W][BG(][BHDWG*2,WK*2W]x[BG)W][FK)]～N(3.427 5×10^-4, 4.201 4×10^-8),将该正态分布作为先验分布,可得超参数μ₀=3.427 5×10^-4,σ²₀=4.201 8×10^-8。为了不让先验均值和方差的影响过大,取k₀=v₀=1。

由于该数据集均值总体μ,σ²均未知,故采用正态—逆伽马分布N—IG(μ_n,k_n,v_n,σ²_n)对其进行后验修正。根据联合后验密度公式可以求出后验参数值μ_n=3.571 7×10^-7,σ²_n=4.208 3×10^-8,k_n=v_n=1 001,即后验正态—逆伽马分布为N—IG(3.571 7×10^-7, 1 001, 1 001, 4.208 3×10^-8),则该正态—逆伽马分布的联合后验密度为:

S1场景联合后分布如图2所示。由图可见,后验概率密度函数的曲面为单峰对称。根据正态近似法将其后验分布近似为正态分布,再对此正态分布进行截断分析。

图2 S1场景联合后验密度分布
Fig.2 The joint posterior density distribution of Scenario 1

由于死亡率≥0,所以截断正态分布的范围为[0,+∞),经式(10)计算取得α=-0.001 7, β=+∞,将其带入上述所服从的截断正态分布,可得φ(1.638 1×10^-4,1.530 8×10^-8,0,+∞,x ),该截断正态分布如图3所示。

当该截断正态分布的逆累积分布函数为95%时,对应的死亡率为4.023 1×10^-4,因此得出烈度为Ⅵ度、发震时间在白天且人口密度≤100时,95%的死亡率区间为[0,4.023 1×10^-4)。

将收集的历史震例按照基础场景分级后,利用各场景下的死亡率数据结合死亡率评估模型,按照S1计算过程计算各场景下的死亡率均值以及95%的死亡率区间,见表5。由表可知,各场景下的死亡率均值及95%区间估计基本符合在极震区烈度相等并且人口密度相同的条件下,夜间的死亡率大于白天的死亡率,但在极震区烈度为Ⅵ度时,人口密度>100人/km²的死亡率均值及95%区间的上限值小于人口密度≤100人/km²的死亡率均值及95%区间的上限值。

图3 场景S1截断正态分布
Fig.3 Truncated normal distribution of Scenario 1

表5 各场景下的死亡率均值及95%死亡率区间
Tab.5 The average death rate and its 95% interval for each scenario

模型的准确性和可靠性,需要通过实际发生的震例进行验证,随机选取中国大陆地区13次地震的人员死亡数据,将其代入本文模型得到预测人员死亡数量,并与实际死亡人数进行对比。利用表5估算各场景下的死亡率与95%死亡人数区间,结果见表6。

表6 本文模型估算的中国大陆13次地震人员死亡率和95%死亡人数区间
Tab.6 The averge death rate and its interval for 13 earthquake in Chinese Mainland evaluated by the model

由表6可见,模型总体评价效果较好,大多数震例的实际死亡人数与模型评价结果在同一数量级且被包含在区间内,评价结果较为准确。由于涉及众多不确定性因素,个别震例评估结果与实际情况存在一定出入。如模型估算的2011年云南盈江5.9级地震95%死亡人数上限值与实际死亡人数相差不大,相对来说较为合理; 而2008年四川攀枝花凉山交界6.1级地震的评估结果与实际死亡人数悬殊较大,此次地震人员伤亡较多的主要原因是大量房屋倒塌和受损(张雪飞等,2008),进而导致评估结果与实际死亡人数相差较大。