基于POT模型的我国地震损失分布研究.pdf

第39 卷第4期2023年8 月Journal of Harbin University of Commerce(Natural Sciences Edition)哈尔滨商业大学学报(自然科学版)Vol.39 No.4Aug.2023基于POT模型的我国地震损失分布研究许一涵12,尤添革12,温芜姚12,宁静12,肖扬岚12,尤学敏3（1福建农林大学计算机与信息学院，福州350 0 0 2;2 福建省统计信息研究中心，福州350 0 0 2；3.国网信通亿力科技有限责任公司，福州350 0 0 2）摘要：确定地震灾害损失分布有利于防灾减灾及灾后赔偿工作的开展.运用超越阈值模型（POT）、马尔科夫链蒙特卡洛估计（MCMC)等研究方法，探索我国地震损失分布特征，并进行参数估计.POT模型能够有效运用极端数据.研究结果显示，我国地震损失具有“尖峰、厚尾”的特征，适合使用广义Pareto分布（GPD）拟合，且通过分布拟合图、残差QQ图等检验.经Hill图、均值超额函数图比较不同值下GPD分布的参数估计，得到阅值=3.4较合适并建立模型，使用极大似然法进行参数估计且通过检验.考虑到小样本情形下极大似然法可能失效的问题，使用蒙特卡洛法估计参数，并比较不同样本量下极大似然法和蒙特卡洛法的参数估计结果，得到小样本情形下蒙特卡洛法估计效果更理想.关键词：地震灾害;GPD分布；厚尾特征；MCMC法；极大似然估计中图分类号：0 2 12.8文献标识码：A文章编号：16 7 2-0 9 46(2 0 2 3)0 4-0 446-0 7Study of earthquake losses distribution in China based on POT modelXU Yihan-2,YOU Tiange*2,WEN Yuanyao-2,NING Jing2,XIAO Yanglan-2,YOU Xuemin(1.College of Computer and Information Sciences,Fujian Agriculture and Forestry University,Fuzhou 350002,China;2.Statistical Information Research Center of Fujian Province,Fuzhou 350002,China;3.State Grid Communication Yili Technology Co.Ltd.,Fuzhou 350002,China)Abstract:Determining the distribution of earthquake damage losses is conducive to disasterprevention and mitigation and post-disaster compensation.In this paper,the BeyondThreshold Model(POT)and the Markov Chain Monte Carlo Estimation(MCMC)were usedto explore the distribution characteristics of earthquake loss in China and estimate theparameters.POT models can effectively use extreme data.The results showed that Chinasseismic loss had the characteristics of“spike and thick tail,which was suitable forgeneralized Pareto distribution(GPD)ftting,and passed the distribution fiting plot,residual QQ plot and other tests.The parameter estimation of GPD distribution under收稿日期：2 0 2 2-10-2 1.基金项目：福建省社会科学规划项目（FJ2018B063）作者简介：许一涵（19 9 9-），女，硕士研究生，研究方向：资源与环境统计.E-mail:第4期different thresholds was compared by Hill plot and mean excess function plot,and thethreshold =3.4 was more suitable and the model was established,and the parameterestimation was carried out by the maximum likelihood method and passed the test.Considering the problem that the maximum likelihood method may fail in the small samplecase,the MCMC method was used to estimate the parameters,and the parameter estimationresults of the maximum likelihood method and MCMC method were compared under differentdata amounts,and the MCMC method was obtained to be more satisfactory in the smallsample case.Key words:earthquake disaster;GPD distribution;thick tails characteristics;MCMCmethod;maximum likelihood estimation许一涵，等：基于POT模型的我国地震损失分布研究447我国是一个极端灾害多发的国家，巨大自然灾害不但会影响民众的生命安全，还会对我国政府造成极大的经济负担，影响国家经济发展.2 0 0 3年非典，2 0 0 8 年四川8.0 级大地震，2 0 10 年西南大旱灾,2 0 2 1年郑州暴雨等等巨灾事件，都对国家经济造成巨大的损失.巨灾的特点是极少发生，但只要发生都将会对各方面或者某一方面造成巨大损失，因此，有必要对巨灾进行相关分析.巨灾事件通常都有“尖峰厚尾”的特性，此时极端数据比中间数据更具有分析价值，研究这部分数据对于研究巨灾风险更具实际意义，极值模型正是为此提供了一个合适的模型.极值模型在水文、降水、地震、工程、环境、金融及保险等方面都有重要应用.近年来，极值理论在巨灾风险建模方面得到了广泛的应用.在地震灾害方面，主要运用极值理论分析其损失数据、震级数据,计算地震最大震级的分布，复发规律并预测,研究其震级上限及重现水平2-5,分析地震灾害损失及其尾部特征6-7,并计算了VaR、ES 等值度量地震灾害风险8-9,以此得到地震保险纯保费和人均保费10.基于VaR值，建立了风险分散层级,还利用了CVaR、PM L这两种风险测度指标改进VaR的测度12,近年来，许多学者还采用混合模型分析我国地震损失13.以上对地震的研究大多基于其极端值，也有学者认为中间数据对于GPD模型的建立也十分重要建立了四种混合模型,并提出用贝叶斯法估计参数和阈值14本文将基于极值理论,综合运用MCMC、极大似然估计等研究方法，分析我国地震损失分布和巨灾保险政策，为研究我国地震损失提供理论基础1模型与方法1.1模型极值理论是统计建模中的一个重要理论，常常用于研究和分析样本数据中罕见情况.极值理论是由Fisher和Tippet15（19 2 8）提出的，随后Gnedenkol10（19 43)对此深入分析,之后 Gumbel 17(1958)将其标准化.极值模型主要由两种方法组成：一种是分块极值方法(Block Maxima Method,BMM）；另一种是超越阈值方法（PeakOverthreshold,POT）.PO T 模型由Smith18（19 8 6)提出，是一种更为有效、应用更加广泛的模型，主要分析数据尾部的渐近分布，也就是广义Pareto 分布(Generalized Pareto Distribution,GPD).相比BMM模型,POT模型更充分利用数据,计算也更准确,是BMM模型的改进.对某固定的大值小于*，成为阈值，若X,大于，则称它为超阈值，此时X，-表示超出量，得到F(x)=P,(X-x|X)=F(+)-F()1-F()一,0称F（x)为超出量分布19.对应的密度函数为fa(x)=x+u)F(u),x0Frm(x)=P,(Xx/X)=E(ca)=F(u),1-F(),xu(3)称为超阈值分布函数,对应的密度函数为(1)(2)448称为X的平均超出量函数.考虑它的极限分布,如果X的分布函数为C(x;,0,5)=1-(1+*_4)*,x,1+*=0则称X服从广义Pareto分布.其中：ER代表位置参数,0代表尺度参数,R代表形状参数.1.2研究方法1.2.1厚尾性检验通过绘制QQ图检验数据的尾部特征，先将数据升序排序，然后利用数据和标准正态分布的各个分位点绘制散点图.分位点计算通常使用：$-(i-)/(-)如果符合正态分布，图像中的大部分点应该围绕着一条直线波动；否则,图像将在一端或者两端有摆动.当经验分位数增速较理论分位数快时，图像将向上摆动，为厚尾分布；相反，图像将向下摆动，为短尾分布.1.2.2阈值的选取POT模型主要考察超阈值的样本数据，阈值的取值大小会影响研究结果.如果阈值太大，可用数据量太少，会导致方差增大；如果阈值太小，样本数据太多,将产生有偏的参数估计，不服从GPD分布.因此，合理选择阈值，既要保证阈值充分大，又要保证数据量足够.本文使用两种图示法确定阈值:1)Hill图法，绘制次序统计量与Hill 估计值的图像,Hill估计式为：(8)1(X(k)Hill图中尾部稳定区域开始时所对应的横坐标X就为所要选取的阈值；2）均值超额函数图法，取超出量为Y,=X,-,均值超额函数是阈值的线性函数,函数为：(9)在图中X时，如果大部分点围绕在一条直线附近波动，就可以选取这个值为阈值.哈尔滨商业大学学报（自然科学版）1.2.3参数估计fru()=)e()=E(X-lx)X(1)第39 卷(4)1）极大似然估计(5)设样本Xf(X,),在参数空间取值,若=T(X)是一个统计量,满足条件：f(X,0)=supf(X,0)则称=T(X)是的极大似然估计.(6)2)Gibbs 抽样设=（1,）是p维参数向量,（D)是观察到数据集D后的后验分布,则基本 Gibbs抽样方法如下：第0 步，任意选取一个初始点(）=（1.o，02.0,p,0),并置i=0;第1步，按下列方法生成(+)=（0 1.+1,2.+1,p.:+1):生成 0 1,+1 T(0 1 02,0p,i,D),(7)生成0 2,i+1 (02 01,i+1,Q3,i,0p,i,D),生成 0 p,+1 (0,01,+1,2,1+,p-1+1,D);第2 步,置i=i+1,并返回到第1步.在这个算法过程中，0 的每一个分量按照自然顺序生成，每一个循环需要生成p个随机变量.2实实证分析2.1数据来源与预处理由于地处板块的相互作用、自然环境的破坏和人类活动的增加，我国经常受到地震的侵扰.地震发生后,民众伤亡惨重、经济损失巨大，政府部门和保险机构都需要投人巨大的资金，例如2 0 0 8 年汶川8.0 级大地震，造成约七万人死亡，损失金额八千多亿.随着我国现代化进程的不断加快，交通、建筑等城市设施不断地优化，人口越来越密集，地震的发生只会使经济损失和人员伤亡越来越大以我国省份空间布局视角，截止2 0 19 年，我国32个省市自治区中有2 7 个发生过地震，其中四川、云南、新疆、青海省、西藏这五个省市自治区发生频率较高,是我国发生地震的主要省份.考虑到19 9 0 年以前数据误差较大，决定以中国大陆19 9 0 2 0 19 年地震数据作为样本数据.数据中震级均为4级以上；记相同地区短时间发生的全部地震和余震的总损失为该次地震经济损失，记(10)第4期最大震级为该次地震震级.本文中所有数据均是国家标准统计数据，来源为各年中国地震年鉴等，数据处理软件为R、O p e n b u g s.由于损失数据都是根据当年经济水平记录的，为了尽量减少经济发展对研究结果的影响，使数据更具说服力，本文按照2 0 19 年的GDP对数据进行调整,即每次的地震损失数据乘以2 0 19 年的GDP，再除以每年的GDP.2.2描述性分析通过表1描述性统计可以看到，数据的最小值与最大值相差巨大；7 5%分位数远小于均值，标准差较大；偏度系数背离正态分布0 特征值，右偏严重；峰度系数明显背离正态分布3特征值.由QQTable 1 Descriptive statistics of earthquake disaster damage data in mainland China(billion yuan)样本容量343图1地震灾害损失数据QQ图Figure 1QQ chart of earthquake disaster damage data90060030001990图2 地震灾害损失数据散点图Figure 2 Scatterplot of earthquake disaster damage data2.3模型估计与分析以上结果可以得出正态分布并不能很好的拟许一涵，等：基于POT模型的我国地震损失分布研究表1中国大陆地震灾害损失数据描述性统计（亿元）最小值最大值0.021 102.35可8004000-3-2-1 0 1 2 3理论分位数金20002010年份449.图1可见，图像尾部呈现下凸的形状，下尾偏离明显，初步认为数据有“尖峰、厚尾”特征.利用散点图2 能够更直观看到图中有一极强影响点发生于2013年四川省芦山县的7 级特大地震，本次特大地震直接经济损失110 2.35亿元，死亡人数196人，受伤人数130 19 人.此次地震的特点是震级较大；地震后余震发生十余次且震级较高，对震区造成持续破坏，交通、通信等设施都有较大程度的毁坏；相比同震级地震，虽然此次地震造成经济损失较大，但是人员伤亡情况较轻，主要由于2008年汶川特大地震灾后重建、加固工作较好，大部分建筑的主体结构并没有倒塌，人群密集的医院、学校等大型建筑毁坏、倒塌程度较轻.均值标准差20.8377.72合地震损失数据，考虑使用GPD分布来拟合.首先对数据进行Box-Cox变换得到变换参数入近似8为0,即对数据进行对数变换.通过图3中的Hill图和均值超额函数图确定阈值，比较不同阈值下1GPD分布参数（,）的极大似然估计及相应的95%置信区间，发现阈值=3.4较合理，此时该值附近呈现平稳状态，建立POT模型得到参数结果如表2 所示，此时超阈值样本数为48 个，数据量满足模型的需要，也就是说这48 个数据属于高损数据.并且通过GPD分布拟合图和残差QQ图检验GPD分布拟合情况，由图4可以看出拟合效果较好.使用极大似然估计法对广义帕累托分布的参数进行估计，得到结果表2、诊断图5,在P-P图和分位数图中只有几个强影响点偏离直线；重现水平图中几乎所有的点都落在了置信区间内；密度曲线的估计和直方图基本一致，说明极大似然法可以估计该模型参数.2020表2 POT模型参数估计表Table 2 POT model parameter estimation table超阅值阅值样本3.448偏度9.65(,)(1.174 6,-0.220 5)峰度116.30数超闯值比例0.139 9(0.222 3,0.126 5)75%分位数11.12标准误差 450.哈尔滨商业大学学报（自然科学版）第39 卷阅值4.52002.36001.24000.135 04420H155085125170215次序图3Hill图和均值超额函数图Figure 3 Hill chart and mean-excess function chart20-4-20246阅值u一0.80.8-x)n0.40.80.400Return Level Plot108641e-01Return period(years)图5累计概率分布图、分位数图、重现水平图和密度函数图Figure 5 Cumulative probability distribution plots,quantile plots,recurrence level plots and density function plots极大似然估计法适用于样本数据较多的情况，当样本数据较少时其参数估计结果不再稳定，因此调用Openbugs软件，使用MCMC估计法中的Gibbs抽样对参数进行估计.利用48 个超阈值数据，参数估计值=1.1746,=-0.2205,反复送代2 0 0 0 0 次，为避免起始点对收敛结果的影响，剔除前40 0 0 次数值进行估计，并通过核密度图、送0.404567810 x(on log scale)图4GPD分布拟合图和残差QQ图Figure 4 GPD distribution fit plots and residual QQ plotsProbability PlotQuantile Plot6.55.03.50.40.8Empirical1e+011e+034567810Ordered Data3.54.8ModelDensity Plot0.80.4034567X代轨迹图、自相关系数图判断迭代效果,如图6 所示，当反复迭代到一定程度时，参数、的核密度图近似正态；迭代轨迹在一个稳定的区间波动；自相关函数值趋向于0,因此可以认为迭代已经收敛，参数取值为1.16 3，取值为-0.146 5.比较极大似然法和MCMC法估计不同阈值时的参数，结果如表3所示，当超阈值样本较多时，两5.56.5第4期者估计结果相近；当超阈值样本较少时，极大似然法变得很不稳定，此时两者的参数估计值相差sigma sample:48 000(eueis)d2.01.000许一涵，等：基于POT模型的我国地震损失分布研究xi sample:48 0002.001.02.0sigma:451:0.8077，因此，小样本情况下选择MCMC法估计参数更加理想。3.0-1.0-0.50Xi0.51.01.50.5sigma1.00-1.00Figure 6 Nuclear density diagram,iterative trajectory diagram and auto-correlation graph表3极大似然法与MCMC法参数估计结果比较Table 3 Comparison of parameter estimation resultsbetween the maximum likelihood method and超闯值阅值极大似然估计(o,)样本数2.9593.4483.9334.4184.99最后，计算地震损失的高分位数点估计（估计式为x,=+%(),如表4所示，说明适合我国的地震保险制度，当地震损失金额小于8 4.4亿元时，保险公司承担大部分的赔偿金额，剩余部分由再保险承担或者国家扶持；当损失金额介于8 4.4 339.2 亿元时，再保险承担大部分的赔偿金额；当损失金额大于339.2 亿元，甚至大于9 13.9 亿元时，国家承担大部分的赔偿金额，0.5-0.51985019.900iteration501ag图6 核密度图、选代轨迹图和自相关系数图MCMC methodMCMC 估计(o,)(1.634 0,-0.344 0)(1.608 0,0.289 3)(1.174 6,-0.220 5)(1.163 0,0.146 5)(0.861 8,-0.089 1)(0.838 2,0.409 0)(0.842 2,-0.109 8)(0.783 2,0.188 6)(1.264 4,0.511 7)(1.001 0,0.296 0)199501985019.900iterationxi1.00-1.00因为此程度的地震通常会引起公共基础设施大面积损毁、建筑物倒塌等等，考虑这样的巨额损失靠单纯的保险和再保险很难承担，国家承担主要的理赔金额.表4我国地震灾害损失金额的高分位数点估计值（亿元）Table 4 High-quartile point estimates of the amount ofearthquake disaster damage in China(billion yuan)0.95分位数0.99分位数0.995分位数0.999分位数84.40339.203结 语本文基于GPD分布与 Gibbs 抽样法对我国地震损失数据展开实证研究，研究结果表明：1)我国地震损失数据具有明显的“尖峰、厚尾”特性，经分布拟合图、残差QQ图等验证了GPD分布能够很好拟合此特性，并通过Hill图、均值超额函数图等确定最优阈值=3.4;2)基于该阈值建立POT模型,使用极大似然法估计参数,并通过P-P图、分位数图、重现水平图等检验；3在参数估计中，考19.950501ag447.20913.90452虑到极大似然法在小样本情形下可能失效的问题，再次采用蒙特卡洛法进行参数估计，并通过调整阈值、控制样本量来比较不同数据量下两者的参数估计结果，发现在小样本情形下蒙特卡洛法比极大似然法估计效果更好；4）通过高分位数点估计值，根据我国实际国情，简单说明适合我国的巨灾保险制度.本文运用POT模型充分利用了样本数据，保留了巨灾风险的厚尾特征，能更好地把握巨灾数据，对于研究巨灾损失具有重要意义.参考文献：1 IRMELAZ,GABRIELE A,EMMANUEL V.BayesianEstimation of the Maximum Magnitude m max Based onthe Extreme Value Distribution for Probabilistic SeismicHazard Analyses J.Pure and Applied Geophysics,2020,177(12).2任晴晴,赵宜宾，钱小仕，等.基于POT模型的昆仑山地区地震统计特征分析J震灾防御技术,2 0 2 2,17(3):529 538.3张锟,任鲁川,田建伟,等基于广义极值理论的潜在地震海啸源震级上限及强震重现水平的估计一以琉球海沟俯冲带为例J中国地震，2 0 16，32(4):702-709.4田建伟,刘哲，任鲁川基于广义帕累托分布的马尼拉海沟俯冲带地震危险性估计J地震，2 0 17，37(1):158-165.5 JAN B,ANDRZEJ K,TOM R,et al.Estimating themaximum possible earthquake magnitude using extremevalue methodology:the Groningen case J.NaturalHazards,2019,98(3).6耿贵珍，王慧彦.基于POT-CPD模型的地震巨灾损失分布研究J自然灾害学报，2 0 16，2 5（3）：153 158.7 AANNE D.Estimation of Tail Distribution of the AnnualMaximum Earthquake Magnitude Using Extreme ValueTheoryJ.Pure and Applied Geophysics,2019,176哈尔滨商业大学学报（自然科学版）8巢文,邹辉文基于POT模型的巨灾损失VaR和CVaR估计J北京化工大学学报：社会科学版，2020(2):18 22.9孟生旺,李政宵。地震死亡人数预测与巨灾保险基金测算J.统计研究，2 0 18，35（10）：8 9-10 2.10 郭静，张连增基于Mixed Erlang-Pareto 组合分布的巨灾风险评估一一以中国地震灾害为例J统计与信息论坛，2 0 2 1，36（3）：119-12 8.11 郝军章，崔玉杰.基于POT模型的巨灾风险度量与保险模式研究一一以地震风险为例J数理统计与管理，2 0 16,35(1)：132-141.12 任婧，张节松基于POT模型的巨灾损失分布拟合及风险度量J.科技与管理，2 0 15，17（1）：75-80.13李云仙，董志伟，钱振伟基于混合模型对地震巨灾风险的分析J.数理统计与管理，2 0 17，36（4）：571-579.14YUNXIAN L,NIANSHENG T,XUEJUN J.Bayesianapproaches for analyzing earthquake catastrophic risk J.Insurance Mathematics and Economics,2016,68.15 RA F,L H C T.Limiting forms of the frequencydistribution of the largest or smallest member of asample J.Mathematical Proceedings of theCambridge Philosophical Society,1928,24(2).16(GNEDENKO B.Sur La Distribution Limite Du TermeMaximum D Une Serie Aleatoire J.Annals ofMathematics,1943,44(3).17GUMBEL E J.Statistics of ExtremesM.S.1.:Columbia University Press.1958.18 SMITH R L.Extreme value theory based on the rlargest annual events J.Journal of Hydrology.1986,86(1-2).19 5史道济实用极值统计方法M.天津：天津科学技术出版社,2 0 0 6.第39 卷(2).