Analysis of Seismic Statistical Characteristics Based on POT Model in Kunlun Mountain Area
-
摘要: 极值统计是研究较少发生但一旦发生即产生极大影响的随机事件的有效方法。本文以地震活动频繁的昆仑山地区作为研究区域，建立了基于广义帕累托分布的超阈值（POT）模型，并讨论了该地区若干地震活动性参数，包括强震震级分布、潜在震级上限、强震平均复发间隔、一定周期内的强震发震概率、一定时期内的重现水平和超定值重现震级。经统计分析得到：该地区震级阈值选定为MS5.5，超阈值期望震级为MS6.81，潜在震级上限高达MS9.08，MS8.0的平均复发间隔仅为66.8年，未来3年该地区发生MS5.5～MS6.5的概率在80%以上，百年重现水平即可达到历史最大震级MS8.1。
-
关键词:
- 广义帕累托分布 /
-
超阈值(
POT)模型 / - 潜在震级上限 /
- 重现水平
Abstract: Extreme value statistics is an effective method to study random events that rarely occur but can cause great impact once they occur. This article takes the Kunlun Mountains area with frequent seismic activities as the research area. We establish a peaks over threshold (POT) model based on the generalized Pareto distribution. Then we discuss several seismic activity parameters, including: strong earthquakes magnitude distribution, upper limit of potential magnitude, average recurrence period, probability of strong earthquakes in a certain period in the future, the recurrence level and expected recurrence magnitude in a certain period. According to statistical analysis, the magnitude threshold of the region is selected as MS5.5. The expected magnitude over the threshold is MS6.81 and the upper limit of potential magnitude is MS9.08. The average recurrence period of MS8.0 is only 66.8 years. The probabilities of MS5.5 ~ MS6.5 are all above 80%. The 100-year return period can reach the historical maximum magnitude MS8.1. -
引言
地震是一种破坏力极强的自然灾害，对地震危险性的认识尤为重要。地震危险性分析是工程地震工作的重要组成部分，对于地震区划、抗震设防、应急救灾等工作的实施具有重要指导意义。地震危险性分析方法主要分为确定性方法和概率性方法（李昌珑，2016）。概率性方法由Cornell于1968年提出，之后多位学者对泊松模型假设进行改进，提出了多种地震危险性模型。近年来，概率模型已成为国内外普遍应用的地震活动性模型（Cornell，1968；胡聿贤，1999；徐伟进等，2012），其采用某种地震动参数在一段时间内的超越概率作为评价地震危险性的指标。传统的概率模型假设地震的发生服从泊松分布，各震级档的地震发生率服从G-R关系（Gutenberg等，1944）。然而，基于大量实际地震资料的研究表明，用G-R关系模型拟合地震震级分布时，在震级高端或低端会出现明显的“掉头”现象（胡聿贤，1999）。由于实际地震资料对G-R关系的偏离，尤其是强震级段的偏离，影响了预测未来强震危险性可信度。因此，对强震危险性模型的改进和拓展一直是地震学者研究的重点。
极值统计模型主要研究较少发生，但一旦发生即产生极大影响的随机事件（Coles，2001）。近年来，极值统计模型被不断应用于地震危险性分析中。陈培善等（1973）在进一步分析地震过程和已有工作的基础上，结合地震震级应有上限的事实，对震级极值分布函数进行了修正，并将新的分布函数应用于我国部分地震带危险性评估中。根据Balkema等（1974）、Pickands（1975）的研究，对于充分大的阈值，随机变量的超出量分布可近似为广义帕累托分布（GPD）。相比于区组最大估计法，该方法能够充分利用大震级区段的信息，适用于历史地震记录时间长但低震级地震记录缺失的地区，便于构建强震活动模型。Pisarenko等（2003）利用广义帕累托分布分析了哈佛地震目录中18个地震区的浅源地震地震矩分布。Huyse等（2010）利用太平洋地震工程研究中心的地面峰值加速度数据和残差数据比较了对数正态分布与广义帕累托分布的拟合优度。钱小仕等（2012）基于广义极值分布（GEV）给出了地震危险性分析的相关公式与方法，研究了我国台湾地区大震发展规律，并在此基础上预测了未来几年发震的危险性。钱小仕等 （2013） 基于广义帕累托分布给出了若干地震活动性参数的估计公式，并讨论了云南地区地震活动性的相关参数。田建伟等（2017）利用广义帕累托分布估计了马尼拉海沟俯冲带潜在地震海啸源区的地震危险性。任梦依（2018）以1930—2016年龙门山地区MS≥4.5历史地震目录为基础，构建了该地区广义帕累托分布的超出量分布模型，估计了该地区震级上限。任晴晴等（2021）利用广义极值分布对巴颜喀拉地块及其周边地区的最大震级进行了极值统计分析，计算了相关的地震危险性参数，并基于Monte Carlo模拟方法验证了运用广义极值分布对研究区域进行极值统计分析的稳定性。
本文首先给出超阈值（Peaks Over Threshold，简称POT）模型的理论描述及基于POT模型的地震危险性参数度量公式；然后以地震活动频繁的昆仑山地区作为研究区域，对该地区1923—2019年历史地震目录数据进行统计分析，根据平均剩余寿命及置信区间等合理确定阈值，构建昆仑山地区超阈值模型，并与G-R关系计算的震级上限进行比较；最后讨论该地区若干地震活动性参数，包括强震震级分布、潜在震级上限、强震平均复发间隔、一定时期内的重现水平和期望重现震级及未来一定时期内的强震发震概率等。
1. 模型与方法
1.1 广义极值模型
设
${X_1},{X_2}, \cdots ,{X_n}$ 为地震震级随机变量序列，其相互独立且服从同一分布$F(x)$ ，记其最大震级为${M_n} = \max ({X_1},{X_2}, \cdots ,{X_n})$ 。若存在$ {a_n} > 0,{b_n} \in R $ 和非退化分布函数$H(x)$ ，使得$\mathop {\lim }\limits_{n \to \infty } P\left(\dfrac{{{M_n} - {b_n}}}{{{a_n}}} \leqslant x\right) = H(x)$ ，则称$H(x)$ 为极值分布，$F(x)$ 属于极值分布$H(x)$ 的最大值吸引场，记作$X \in MDA(H)$ 或$F \in MDA(H)$ 。Fisher等(1928)获得了极值分布的3种形式：I型${\text{Gumbel}}$ 分布$(\xi = 0)$ ，Ⅱ型${{{\rm{Fr}}{\dot {\rm{e}}}{\rm{chet}}}}$ 分布$(\xi > 0)$ ，Ⅲ型${\text{Weibull}}$ 分布$(\xi < 0)$ ，可统一为如下的广义极值分布：$$ {H_\xi }(x;\mu ,\sigma ,\xi ) = \exp \left\{ - {\left(1 + \xi \dfrac{{x - \mu }}{\sigma }\right)^{ - \tfrac{1}{\xi }}}\right\} {\text{，}}1 + \xi \dfrac{{x - \mu }}{\sigma } \geqslant 0{\text{，}}\sigma > 0{\text{，}}\xi \in {{\rm{R}}} {\text{，}}\mu \in {{\rm{R}}} $$ (1) 式中，
$\mu $ 为位置参数；$\sigma $ 为尺度参数；$\xi $ 为形状参数，是广义极值分布最重要的参数，也称为广义极值分布的极值指数（尾指数）。1.2 POT模型
若记
$ {x^ + } = \sup \{ x:0 < F(x) < 1\} $ ，为分布函数的右端点，则称$ {F}_{u}(x)=P\{X-u\leqslant x|X > u\}=\dfrac{F(x+u)-F(u)}{1-F(u)}(x\geqslant 0) $ 为地震震级x的超过阈值超出量的分布函数，简称超出量分布。Balkema等(1974)与Pickands (1975)指出，当震级分布F(x)属于极值分布H(x)的最大值吸引场时，有：$$ \mathop {\lim }\limits_{u \to {x^ + }} \mathop {\sup }\limits_{0 \leqslant x \leqslant {x^ + } - u} |{F_u}(x) - G(x,\tilde \sigma ,\xi )| = 0 $$ (2) 式中，
$G（x,\tilde \sigma ,\xi ） = 1 - {\left(1 + \xi \dfrac{x}{{\tilde \sigma }}\right)^{ - 1/\xi }},x \geqslant 0,1 + \xi \dfrac{x}{{\tilde \sigma }} \geqslant 0,\tilde \sigma = \sigma + \xi (u - \mu )$ 为两参数广义帕累托分布，即若最大震级Mn近似服从广义极值分布${H_\xi }（x,\mu ,\sigma ,\xi ）$ ，则对于充分大的阈值u，震级超出量X-u近似服从广义帕累托分布$ G（x,\tilde \sigma ,\xi ） $ ，且二者具有相同的形状参数$ \xi $ 。$ \xi $ 的取值反映了分布尾部的收敛性质，$ \xi $ 取值越大，尾部越厚，尾分布收敛速度越慢。对于分布$ G（x,\tilde \sigma ,\xi ） $ ，当$ \xi \geqslant 0 $ 时支撑集为$ [0, + \infty ） $ ，当$ \xi < 0 $ 时支撑集为有限集$ [0, - \sigma /\xi ） $ 。基于广义帕累托分布对超过某充分大的阈值所有观测数据进行极值统计建模，渐近地刻画分布的尾部特征模型，将其称为POT模型。1.3 模型估计
记
$\bar F(x) = 1 - F(x)$ ，表示震级分布$F(x)$ 的尾部，也称为生存函数；$ {N_u} = \displaystyle\sum\limits_{i = 1}^n {{I_{\left\{ {{X_i} > u} \right\}}}} $ 表示序列${X_1},{X_2}, \cdots ,{X_n}$ 中超过阈值u的次数；$ {\Delta _n}(u) = \left\{ {i:{X_i} > u} \right\} $ 表示超阈值u的观测值下标记集；${\bar F_n}(u) = \dfrac{1}{n}\displaystyle\sum\limits_{i = 1}^n {{I_{\left\{ {{X_i} > u} \right\}}}} = \dfrac{{{N_u}}}{n}$ 表示尾部$\bar F(u)$ 的经验分布函数；$ e(u) = E(X - u|X > u) $ 表示震级X超过阈值u时的平均超出量函数，基于广义帕累托分布$ G(x,\tilde \sigma ,\xi ) $ 的平均超出量函数为$ e(u) = \dfrac{{\tilde \sigma }}{{1 - \xi }} = \dfrac{{\sigma + \xi (u - \mu )}}{{1 - \xi }} $ ，即e(u)为u的线性函数；$ {e_n}(u) = \dfrac{1}{{{N_u}}}\displaystyle\sum\limits_{i \in {\Delta _n}(u)} {({X_i} - u)} ,u > 0 $ 表示样本平均超出量函数。首先给出POT模型中广义帕累托分布参数
$ \xi $ 和$ \tilde \sigma $ 的估计值。假设${x_1},{x_2}, \cdots ,{x_n}$ 为地震震级观测数据，u为阈值，$ N = {N_u} $ 为观测数据的超出量，则超出量$ {Y_1},{Y_2}, \cdots ,{Y_N} $ 近似服从参数为$ \tilde \sigma ,\xi $ 的广义帕累托分布，其对数似然函数为：$$ l(\tilde \sigma ,\xi ) = - n\log \tilde \sigma - \left(\dfrac{1}{\xi } + 1\right)\sum\limits_{i = 1}^N {\log \left(1 + \dfrac{\xi }{{\tilde \sigma }}{y_i}\right)} $$ (3) $ \tilde \sigma {\text{，}}\xi $ 的极大似然估计无解析解，通常借助于数值解法求得其极大似然估计$ \hat{\tilde{\sigma }}和\hat{\xi } $ 。类似于广义极值分布，当$ \xi > - 0.5 $ 时，极大似然估计的统计性质较好，根据史道济（2006）的研究，当$ n \to \infty $ 时，由观测信息矩阵的逆可得到参数$ \tilde \sigma ,\xi $ 的渐进协方差矩阵的估计为：$$ {{\boldsymbol{Cov}}_{\tilde \sigma ,\,\xi }} = \left[ {\begin{array}{*{20}{c}} {D(\tilde \sigma )}&{Cov(\tilde \sigma ,\xi )} \\ {Cov(\tilde \sigma ,\xi )}&{D(\xi )} \end{array}} \right] $$ (4) 在一定的正则条件下，当n充分大时，根据渐进正态性，
$ \hat{\tilde{\sigma }}，\hat{\xi } $ 的$ （\text{1}－\alpha ） $ 置信区间分别为：$$ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left( {\hat {\tilde \sigma} \pm {Z_{1 - \alpha /2}}\sqrt {{k_{11}}} } \right){\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left( {0 < \alpha < 1} \right) $$ (5) $$ {\kern 1pt} \left( {\hat \xi \pm {Z_{1 - \alpha /2}}\sqrt {{k_{22}}} } \right){\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left( {0 < \alpha < 1} \right) $$ (6) 式中，
${k_{11}},{k_{22}}$ 分别为估计的协方差矩阵$ {\boldsymbol{Co}}{{\boldsymbol{v}}_{\tilde \sigma ,\,\xi }} $ 主对角线上的元素。由Delta方法可得到下分位数
$ {x_p} $ 的估计$ {\hat x_p} $ 方差为：$$ Var（{\hat x_p}） \approx \nabla x_p^{\rm{T}}V\nabla {x_p} $$ (7) 式中，
$ \nabla {x_p} $ 为$ \left(\dfrac{\partial {x}_{p}}{\partial \overline{F}（u）}，\dfrac{\partial {x}_{p}}{\partial \tilde{\sigma }}，\dfrac{\partial {x}_{p}}{\partial \xi }\right) $ 在点$ （\hat {\bar F}（u）,\hat{ \tilde \sigma} ,\hat \xi ） $ 处的值，其中矩阵V为：$$ {\boldsymbol{V}} = \left[ {\begin{array}{*{20}{c}} {\hat {\bar F}(u)[1 - \hat {\bar F}(u)]/n{\kern 1pt} {\kern 1pt} }&0 \\ 0&{Co{v_{\tilde \sigma ,\, \xi }}} \end{array}} \right] $$ (8) $ {\hat x_p} $ 的$ （\text{1}－\alpha ） $ 置信区间近似为：$$ \left({\hat{x}}_{p}\pm {Z}_{\text{1-}\alpha \text{/2}}\sqrt{Var（{\hat{x}}_{p}）}\right)\text{，}0 < \alpha < 1 $$ (9) $\bar F(u + x) = \bar F(u){\bar F_u}(x)$ ，则分布函数$ F(x) $ 的尾部估计为：$$ \hat {\bar F}(u + x) = \hat {\bar F}(u){\hat {\bar F}_u}(x) = {\bar F_n}(u)\bar G(x,\hat {\tilde \sigma} ,\hat \xi ) = \dfrac{{{N_u}}}{n}{\left(1 + \hat \xi \dfrac{x}{{\hat {\tilde \sigma} }}\right)^{ - 1/\hat \xi }}{\text{，}}x \geqslant 0 $$ (10) 对应的超阈值震级分布函数的估计为：
$$ F(y) = 1 - \dfrac{{{N_u}}}{n}{\left(1 + \hat \xi \dfrac{{y - u}}{{\hat {\tilde \sigma} }}\right)^{ - 1/\hat \xi }}\text{，}y \geqslant u $$ (11) 1.4 地震危险性参数
（1）潜在震级上限
地震数据为日观测数据，一年按365天计算，若超阈值分布为广义帕累托分布
$ {G_u}(x,\tilde \sigma ,\xi ) $ ，则当形状参数$ \xi < 0 $ 时，随机变量X具有有限的支撑集，右端点即为潜在震级上限，其估计值为：$$ {\hat x^ + } = u - {{\hat {\tilde \sigma} } \mathord{\left/ {\vphantom {{\hat {\tilde \sigma} } {\hat \xi }}} \right. } {\hat \xi }} $$ (12) （2）平均复发间隔
以年为单位，震级x的平均复发间隔为：
$$ T(x) = \dfrac{1}{{365[1 - {G_u}(x;\tilde \sigma ,\xi )]}} $$ (13) （3）重现水平
T年重现期的重现水平相当于
$ p = 1 - \dfrac{1}{{365 T}} $ 下分位数$ {x_p} $ ，其估计值为：$$ {\hat x_p} = u + \dfrac{{\hat {\tilde \sigma} }}{{\hat \xi }}\left\{ {{{\left[ {\dfrac{n}{{{N_u}}}(1 - p)} \right]}^{ - \hat \xi }} - 1} \right\} $$ (14) （4）超定值期望震级
称
$ E（X|X > \tau ） $ 为震级超过$ \tau （u < \tau < {x}^{+}） $ 时的期望震级，该参数旨在确定当发生超过某一震级的实际震级，其值记为$ \bar \tau $ ，则有：$$ \bar \tau = E(X|X > \tau ) = \dfrac{1}{{\bar F(\tau )}}\int_\tau ^{{x^ + }} {x{\rm{d}}F(x)} = \tau + e(\tau ) $$ (15) 其估计值为：
$$ \hat {\bar \tau } = \tau + \dfrac{{\hat {\tilde \sigma }}}{{1 - \hat \xi }} + \dfrac{{\hat \xi (\tau - u)}}{{1 - \hat \xi }} $$ (16) （5）定周期内发震概率
未来T年内，发生超过x级地震的概率为：
$$ P(\min \left\{ {m:{X_m} \geqslant m} \right\} \leqslant 365T) = 1 - {[{G_u}(x,\tilde \sigma ,\xi )]^{365T}} $$ (17) 2. 昆仑山地区POT模型构建
2.1 数据资料选取
昆仑山及其周边地区地震活动频繁，中强震频发，为该地区地震危险性评估提供了丰富的数据资源。本文选取昆仑山地区（35oN～38oN，75oE～102oE）作为研究区域。地震震级采用面波震级MS，余震删除时采用C-S余震时空窗方法（陈凌等，1998）。经余震删除后，自1900年开始的MS2.0以上地震频数直方图如图1（a）所示。结合黄玮琼等（1994）对中国大陆大部分地区历史地震资料的研究结果及M-t图（图1（b））可知，昆仑山区域 MS≥4¾地震资料约在1923年后基本完整。本文截取1923—2019年MS≥2.0地震目录，共12 000余条地震记录作为基础研究数据。其中MS≥5.0地震共257次，MS≥6.0地震共41次，MS≥7.0地震共11次，MS≥8.0地震共1次。
2.2 阈值选取与模型估计
图2所示为昆仑山地区震级的平均剩余寿命及置信区间。由图2可知，当阈值u＞5.5时，平均剩余寿命近似线性，所以选取5.5为阈值。为进一步检验阈值u=5.5是否合适，选择一系列阈值，对不同阈值利用极大似然估计得到一系列参数值，包括修正的尺度参数
${\tilde \sigma ^*}{\text{ = }}\tilde \sigma - \xi u$ 和形状参数$\xi $ ，然后确定${\tilde \sigma ^*}$ 和$\xi $ 关于阈值u的图形及相应的置信区间，选择使估计值可保持为常数的最小值u作为阈值。如果参数估计值在所选阈值附近是稳定的，说明阈值合适，阈值应尽量大，以保证模型的正确性。在3.0～6.2选取33个值作为阈值，用广义帕累托分布模型进行估计，结果如图3所示。由图3可知，对于较大的阈值（>5.5），参数估计值的变化与平均剩余寿命的变化类似，但相对于抽样误差，参数估计值的扰动较小。因此，选取阈值u=5.5是合适的。在选取的昆仑山地区地震目录中，超阈值u=5.5的震级基本信息如表1所示。其中平均震级为MS6.20，最大震级为2011年11月14日发生在昆仑山口西的MS8.1巨震（徐昊等，2018）。表 1 超阈值震级基本信息Table 1. Basic information of over threshold magnitudes最小值 四分之一分位数 中位数 平均值 四分之三分位数 最大值 极差 标准差 5.59 5.80 6.00 6.20 6.43 8.10 2.51 0.58 根据上述地震目录数据，超阈值广义帕累托分布模型参数的极大似然估计为：
$$ （\hat{\tilde{\sigma }},\hat{\xi }）=（0.8625\text{，}-0.2410） $$ (18) 参数的协方差矩阵估计为：
$$ {{\boldsymbol{Cov}}_{\tilde \sigma ,\xi }} = \left[ {\begin{array}{*{20}{c}} {0.01843}&{ - 0.01228} \\ { - 0.01228}&{0.01169} \end{array}} \right] $$ (19) 昆仑山地区超出量广义帕累托分布拟合诊断结果如图4所示。由图4可知，样本数据符合广义帕累托分布；样本数据均在广义帕累托分布分位数估计的置信区间内，且重现水平为凸曲线；广义帕累托分布密度函数曲线与频率直方图的走势一致。综上所述，利用广义帕累托分布模型对昆仑山地区进行统计分析是合理的。
图 4 昆仑山地区震级的阈值超出量模型诊断结果Figure 4. Model fitting diagnosis chart of Kunlun mountains area计算得到尺度参数
$ \hat {\tilde \sigma} $ 的95%置信区间为$ （0.8625 \pm 0.2661） $ ，形状参数$ \xi $ 的95%置信区间为（$ - 0.2410 \pm 0.2119 $ ）。形状参数$ \xi $ 的95%置信区间均为负值，显示分布支撑集上端点为有限值，即存在有限的潜在震级上限，与实际情况吻合。在地震目录中，超阈值u=5.5的数据共72个，所以尾分布
$ \bar F（u） $ 的极大似然估计为：$ \hat{ \bar F}（u）{\text{ = 0}}{\text{.005940}} $ ，则$ \left(\hat{\overline{F}}（u），\tilde{\sigma }，\xi \right) $ 的协方差矩阵为：$$ {\boldsymbol{V}} = \left[ {\begin{array}{*{20}{c}} {4.8716 \times {{10}^{{{ - 7}}}}}&0&0 \\ 0&{0.01843}&{ - 0.01228} \\ 0&{ - 0.01228}&{0.01169} \end{array}} \right] $$ (20) 震级分布为：
$$ F(y) = 1 - 0.005964{[1 - 0.2794(y - u)]^{4.1499}},y \geqslant 5.5 $$ (21) Gutenberg等（1956）提出的震级-频度经验公式也称为G-R关系式，表达式为：
$$ \lg N(m) = a - bm,m \geqslant {M_{\min }} $$ (22) 式中，
$ N（m） $ 为震级≥$ m $ 的地震次数。G-R关系式已成为地震学中基本的统计公式之一。当
$ N（m） = 1 $ 时，对应的震级可视为研究区内的理论最大震级$ {M_{{\rm{theo}}}} $ ，即$ {M_{{\rm{theo}}}} = a/b $ 。图5所示为昆仑山1923—2019年震级-频度拟合关系。由图5可知，MS≥2.5的地震较好地符合G-R关系式，因此本文研究昆仑山地区地震危险性时震级下限取为MS2.5。根据地震目录数据，相关计算结果如表2所示。表 2 G-R关系计算结果Table 2. Calculation results of G-R relationship起始时间/年 $ {M_{\min }} $ a b ${M_{{\rm{theo}}} }$ 1923 2.5 12.956 5 1.536 9 8.43 Sun等（1995）给出了基于G-R关系的泊松（指数）分布模型，该模型被广泛应用于地震危险性评估中（苏有锦等，2011）。将本文POT模型分布与基于G-R关系的泊松（指数）分布进行比较，如图6所示。由图6可知，在中强震处，本文模型更接近于实际地震数据，明显优于泊松分布。这主要是因为本文的广义帕累托分布理论并不是对所有震级数据建模，而主要是针对超阈值的数据建模，根据相关理论，阈值往往较大，使本文模型的适用性主要体现在中强震级段。此外，指数分布无震级上限，这与实际情况不符，而本文广义帕累托分布模型拟合结果中形状参数
$ \xi $ 的95%置信区间均为负值，说明震级存在潜在上限，与实际情况吻合，这也是本文模型的优势。
图 6 本文模型与泊松（指数）分布比较Figure 6. Comparison between the model in this paper and the poisson （exponential） distribution3. 地震危险性分析
潜在震级上限的估计为：
$ {\hat x^ + } = u - {{\hat{ \tilde \sigma }} \mathord{\left/ {\vphantom {{\hat {\tilde \sigma} } {\hat \xi }}} \right. } {\hat \xi }}{\text{ = }}9.08 $ 。超阈值u=5.5的期望震级为MS6.81。表3所示为强震复发间隔的估计，其中MS6.0地震的平均复发间隔为0.9年，这与其他学者关于西部地区强震发震周期短的结论一致（张国民等，2005；任晴晴等，2013）。由表3可知，MS8.0及以上巨震的复发间隔为66.8年，考虑到该地区2011年已发生过MS8.1地震，所以未来50年左右该地区发生8.0级以上巨震的概率较小。表 3 不同震级复发间隔及发震概率Table 3. Recurrence cycle and occurrence probability of different magnitudes震级MS/级 5.5 6.0 6.5 7.0 7.5 8.0 8.5 平均复发间隔/年 0.5 0.9 1.8 4.4 13.8 66.8 882.9 1年内发震概率 0.886 4 0.687 5 0.427 1 0.203 6 0.070 1 0.014 9 0.001 1 3年内发震概率 0.998 5 0.969 5 0.812 0 0.494 9 0.196 0 0.043 9 0.003 4 5年内发震概率 1.000 0 0.997 0 0.938 3 0.679 7 0.304 8 0.072 2 0.005 6 10年内发震概率 1.000 0 1.000 0 0.996 2 0.897 4 0.516 7 0.139 1 0.011 3 由表3可知，7级以上大震发生的可能性均超过50%。未来5年MS5.5及以上地震的发震概率均超过93%，说明该地区短期内中强震和强震仍较活跃；7级以上的强震发震概率达68%，应引起关注。未来10年7级以上的强震发震概率超过90%。
表4所示为昆仑山地区不同重现期的重现水平估计结果，即给定年限内可能出现的最大震级，也是高分位数
$ {x_p} $ 的点估计。由表4可知，随着周期的增加，重现水平不断提高，因时间增加带来的不确定性增大，导致重现水平的95%置信区间长度增大，即精度降低，这与实际情况吻合；未来1年的重现水平为MS6.11，在95%置信水平下，估计的取值区间为（5.89，6.33）；未来5年的重现水平高达MS7.06，且超过MS7.06的平均震级为MS7.46；百年重现水平可达历史最大震级MS8.1，较基于G-R关系计算得到的理论最大震级MS8.43更接近于真实数据；百年重现水平的期望震级高达MS8.30。表 4 重现水平Table 4. Recurrence level项目 周期/年 1 2 5 10 20 50 100 重现水平（MS） 6.11 6.57 7.06 7.37 7.64 7.92 8.10 95%置信区间 （5.89，6.33） （6.29，6.85） （6.77，7.35） （7.03，7.71） （7.23，8.05） （7.37，8.47） （7.43，8.77） 超定值期望震级（MS） 6.69 7.05 7.46 7.71 7.92 8.15 8.30 4. 结论与讨论
一般情况下，地震目录中强震的记录较充分，POT模型可充分利用地震目录中强震级区段的数据资料，在小震资料不充足的情况下也可构建潜源区的地震活动模型，其统计结果可作为地震危险性分析和结构抗震等研究工作的参考。
本文首先给出POT模型中的相关理论和衡量地震危险性参数的公式；然后对昆仑山地区的地震目录进行POT建模和实证分析，在实证分析过程中，选定研究的起始时间为1923年，并结合平均剩余寿命及修正的尺度参数和形状参数随阈值变化情况，选定昆仑山地区地震震级阈值为MS5.5；最后对超阈值震级进行极值统计分析。
根据拟合结果可知，昆仑山地区的震级超出量服从广义帕累托分布。经模型计算，该地区超阈值期望震级为MS6.81，潜在震级上限为MS9.08。未来3年该地区发生MS5.5～MS6.5中强震和强震的概率较大，均超过80%，说明短期内该地区中强震和强震依然活跃。中强震及以上地震的复发间隔相对较短，MS8.0以上巨震的复发间隔仅为66.8年，因该地区2011年已发生了MS8.1巨震，故未来50年左右该地区发生MS8.0以上巨震的概率较小。未来5年的重现水平高达MS7.06，且超过MS7.06的平均震级为MS7.46。百年重现水平可达历史最大震级MS8.10，较基于G-R关系计算得到的理论最大震级MS8.43更接近于真实数据。
本文基于POT模型对昆仑山地区的地震危险性分析结果的可靠性和合理性需进一步研究，模型稳定性及前震、余震在POT模型中的影响也需进一步研究。同时可考虑将地质构造结合到地震危险性极值分析框架中，多学科的交叉应用会得到更合理的结果。
致谢 中国西部环境与生态科学数据中心和国家地震科学数据共享中心为本文提供地震目录数据，审稿专家对本文进行了专业、细致的评阅，并提出了宝贵的修改意见，在此一并表示感谢。
-
图 4 昆仑山地区震级的阈值超出量模型诊断结果
Figure 4. Model fitting diagnosis chart of Kunlun mountains area
图 6 本文模型与泊松（指数）分布比较
Figure 6. Comparison between the model in this paper and the poisson （exponential） distribution
表 1 超阈值震级基本信息
Table 1. Basic information of over threshold magnitudes
最小值 四分之一分位数 中位数 平均值 四分之三分位数 最大值 极差 标准差 5.59 5.80 6.00 6.20 6.43 8.10 2.51 0.58 
下载: 导出CSV
表 2 G-R关系计算结果
Table 2. Calculation results of G-R relationship
起始时间/年 $ {M_{\min }} $ a b ${M_{{\rm{theo}}} }$ 1923 2.5 12.956 5 1.536 9 8.43
下载: 导出CSV
表 3 不同震级复发间隔及发震概率
Table 3. Recurrence cycle and occurrence probability of different magnitudes
震级MS/级 5.5 6.0 6.5 7.0 7.5 8.0 8.5 平均复发间隔/年 0.5 0.9 1.8 4.4 13.8 66.8 882.9 1年内发震概率 0.886 4 0.687 5 0.427 1 0.203 6 0.070 1 0.014 9 0.001 1 3年内发震概率 0.998 5 0.969 5 0.812 0 0.494 9 0.196 0 0.043 9 0.003 4 5年内发震概率 1.000 0 0.997 0 0.938 3 0.679 7 0.304 8 0.072 2 0.005 6 10年内发震概率 1.000 0 1.000 0 0.996 2 0.897 4 0.516 7 0.139 1 0.011 3
下载: 导出CSV
表 4 重现水平
Table 4. Recurrence level
项目 周期/年 1 2 5 10 20 50 100 重现水平（MS） 6.11 6.57 7.06 7.37 7.64 7.92 8.10 95%置信区间 （5.89，6.33） （6.29，6.85） （6.77，7.35） （7.03，7.71） （7.23，8.05） （7.37，8.47） （7.43，8.77） 超定值期望震级（MS） 6.69 7.05 7.46 7.71 7.92 8.15 8.30
下载: 导出CSV
-
陈凌, 刘杰, 陈颙等, 1998. 地震活动性分析中余震的删除. 地球物理学报, 41（S1）: 244—252Chen L. , Liu J. , Chen Y. , et al. , 1998. Aftershock deletion in seismicity analysis. Acta Geophysica Sinica, 41(S1): 244—252. (in Chinese) 陈培善, 林邦慧, 1973. 极值理论在中长期地震预报中的应用. 地球物理学报, 16（1）: 6—24Chen P. S. , Lin B. H. , 1973. An application of statistical theory of extreme values to moderate and long interval earthquake prediction. Acta Geophysica Sinica, 16(1): 6—24. (in Chinese) 胡聿贤, 1999. 地震安全性评价技术教程. 北京: 地震出版社. 黄玮琼, 李文香, 曹学锋, 1994. 中国大陆地震资料完整性研究之二——分区地震资料基本完整的起始年分布图象. 地震学报, 16（4）: 423—432. 李昌珑, 2016. 时间相依的地震危险性区划研究及应用. 北京: 中国地震局地球物理研究所.Li C. L., 2016. Study and application of time-dependent seismic hazard zonation. Beijing: Institute of Geophysics, China Earthquake Administration. (in Chinese) 钱小仕, 王福昌, 曹桂荣等, 2012. 广义极值分布在地震危险性分析中的应用. 地震研究, 35（1）: 73—78 doi: 10.3969/j.issn.1000-0666.2012.01.013Qian X. S. , Wang F. C. , Cao G. R. , et al. , 2012. Application of the generalized extreme value distribution in seismic hazard analysis. Journal of Seismological Research, 35(1): 73—78. (in Chinese) doi: 10.3969/j.issn.1000-0666.2012.01.013 钱小仕, 王福昌, 盛书中, 2013. 基于广义帕累托分布的地震震级分布尾部特征分析. 地震学报, 35（3）: 341—350 doi: 10.3969/j.issn.0253-3782.2013.03.006Qian X. S. , Wang F. C. , Sheng S. Z. , 2013. Characterization of tail distribution of earthquake magnitudes via generalized Pareto distribution. Acta Seismologica Sinica, 35(3): 341—350. (in Chinese) doi: 10.3969/j.issn.0253-3782.2013.03.006 任梦依, 2018. 龙门山地区的地震活动性广义帕累托模型构建. 地震研究, 41（2）: 226—232 doi: 10.3969/j.issn.1000-0666.2018.02.010Ren M. Y. , 2018. The establishment of generalized Pareto distribution model of seismicity in Longmenshan region. Journal of Seismological Research, 41(2): 226—232. (in Chinese) doi: 10.3969/j.issn.1000-0666.2018.02.010 任晴晴, 钱小仕, 赵玲玲等, 2013. 中国大陆活动地块边界带最大震级分布特征研究. 地震, 33（3）: 67—76 doi: 10.3969/j.issn.1000-3274.2013.03.008Ren Q. Q. , Qian X. S. , Zhao L. L. , et al. , 2013. Characteristics of maximum magnitude distributions for active block boundaries in China’s mainland. Earthquake, 33(3): 67—76. (in Chinese) doi: 10.3969/j.issn.1000-3274.2013.03.008 任晴晴, 陆丽娜, 钱小仕等, 2021. 巴颜喀拉地块及其周边地震危险性分析. 地震, 41（3）: 144—156 doi: 10.12196/j.issn.1000-3274.2021.03.011Ren Q. Q. , Lu L. N. , Qian X. S. , et al. , 2021. Earthquake hazard analysis of the Bayankala block and its surroundings. Earthquake, 41(3): 144—156. (in Chinese) doi: 10.12196/j.issn.1000-3274.2021.03.011 史道济, 2006. 实用极值统计方法. 天津: 天津科学技术出版社. 苏有锦, 李忠华, 2011. 云南地区6级以上强震时间分布特征及其概率预测模型研究. 地震研究, 34（1）: 1—7. doi: 10.3969/j.issn.1000-0666.2011.01.001Sun Y. J. , Li Z. H. , 2011. Interval Distribution and Probability Model of the Strong Earthquakes with M≥6.0 in Yunnan. Journal of Seismological Research, 34(1): 1—7. (in Chinese) doi: 10.3969/j.issn.1000-0666.2011.01.001 田建伟, 刘哲, 任鲁川, 2017. 基于广义帕累托分布的马尼拉海沟俯冲带地震危险性估计. 地震, 37（1）: 158—165 doi: 10.3969/j.issn.1000-3274.2017.01.016Tian J. W. , Liu Z. , Ren L. C. , 2017. Seismic hazard estimation of the Manila trench subduction zone based on generalized Pareto distribution. Earthquake, 37(1): 158—165. (in Chinese) doi: 10.3969/j.issn.1000-3274.2017.01.016 徐昊, 孙玉军, 吴中海, 2018. 岩石圈结构对大地震震后形变的影响——以1976年唐山大地震和2001年昆仑山大地震为例. 地球物理学报, 61（8）: 3170—3184 doi: 10.6038/cjg2018L0637Xu H. , Sun Y. J. , Wu Z. H. , 2018. The effect of lithospheric structure on the seismic deformation—taking the 1976 Tangshan earthquake and 2001 Kunlunshan earthquake as an example. Chinese Journal of Geophysics, 61(8): 3170—3184. (in Chinese) doi: 10.6038/cjg2018L0637 徐伟进, 高孟潭. 2012. 根据截断的G-R模型计算东北地震区震级上限. 地球物理学报, 55（5）: 1710—1717Xu W. J., Gao M. T., 2012. Calculation of upper limit earthquake magnitude for Northeast seismic region of China based on truncated G-R model. Chinese Journal of Geophysics, 55(5): 1710—1717. (in Chinese) 张国民, 马宏生, 王辉等, 2005. 中国大陆活动地块边界带与强震活动. 地球物理学报, 48（3）: 602—610 doi: 10.3321/j.issn:0001-5733.2005.03.018Zhang G. M. , Ma H. S. , Wang H. , et al. , 2005. Boundaries between active-tectonic blocks and strong earthquakes in the China mainland. Chinese Journal of Geophysics, 48(3): 602—610. (in Chinese) doi: 10.3321/j.issn:0001-5733.2005.03.018 Balkema A. A. , de Haan L. , 1974. Residual life time at great age. The Annals of Probability, 2(5): 792—804. Coles S., 2001. An introduction to statistical modeling of extreme values. London: Springer. Cornell C. A. , 1968. Engineering seismic risk analysis. Bulletin of the Seismological Society of America, 58(5): 1583—1606. doi: 10.1785/BSSA0580051583 Fisher R. A. , Tippett L. H. C. , 1928. Limiting forms of the frequency distribution of the largest or smallest member of a sample. Mathematical Proceedings of the Cambridge Philosophical Society, 24(2): 180—190. doi: 10.1017/S0305004100015681 Gutenberg B. , Richter C. F. , 1944. Frequency of earthquakes in California. Bulletin of the Seismological Society of America, 34(4): 185—188. doi: 10.1785/BSSA0340040185 Gutenberg B. , Richter C. F. , 1956. Earthquake magnitude, intensity, energy, and acceleration: (second paper). Bulletin of the Seismological Society of America, 46(2): 105—145. doi: 10.1785/BSSA0460020105 Huyse L. , Chen R. , Stamatakos J. A. , 2010. Application of generalized Pareto distribution to constrain uncertainty in peak ground accelerations. Bulletin of the Seismological Society of America, 100(1): 87—101. doi: 10.1785/0120080265 Pickands J. , 1975. Statistical inference using extreme order statistics. The Annals of Statistics, 3(1): 119—131. Pisarenko V. F. , Sornette D. , 2003. Characterization of the frequency of extreme earthquake events by the generalized Pareto distribution. Pure and Applied Geophysics, 160(12): 2343—2364. doi: 10.1007/s00024-003-2397-x Sun J. C. , Pan T. C. , 1995. The probability of very large earthquakes in Sumatra. Bulletin of the Seismological Society of America, 85(4): 1226—1231. doi: 10.1785/BSSA0850041226 期刊类型引用(2)
1. 周仰新，李永强，李明吉，潘泓序，许言，陈波. 基于Folium地理信息可视化的地震数据统计分析. 河北工程大学学报(自然科学版). 2025(02): 89-93 . 
百度学术2. 许一涵，尤添革，温芫姚，宁静，肖扬岚，尤学敏. 基于POT模型的我国地震损失分布研究. 哈尔滨商业大学学报(自然科学版). 2023(04): 446-452 . 百度学术其他类型引用(1)
-
下载:
点击查看大图
计量
- 文章访问数: 228
- HTML全文浏览量: 41
- PDF下载量: 10
- 被引次数: 3
