Global sensitivity analysis of the generalized Pareto distribution model for seismicity in the northeast Tibetan
-
摘要: 基于广义帕累托分布构建地震活动性模型，因其输入参数取值难以避免不确定性，导致依据该模型所得的地震危险性估计结果具有不确定性。鉴于此，本文选取青藏高原东北缘为研究区，提出了基于全域敏感性分析的地震危险性估计的不确定性分析流程和方法。首先，利用地震活动性广义帕累托模型，进行研究区地震危险性估计；然后，选取地震记录的起始时间和震级阈值作为地震活动性模型的输入参数，采用具有全域敏感性分析功能的E-FAST方法，对上述两个参数的不确定性以及两参数之间的相互作用对地震危险性估计不确定性的影响进行定量分析。结果表明：地震危险性估计结果（不同重现期的震级重现水平、震级上限及相应的置信区间）对两个输入参数中的震级阈值更为敏感；不同重现期的地震危险性估计结果对震级阈值的敏感程度不同；对不同的重现期而言，在影响地震危险性估计结果的不确定性上，两个输入参数之间存在非线性效应，且非线性效应程度不同。本文提出的不确定性分析流程和方法，可以推广应用于基于其它类型地震活动性模型的地震危险性估计不确定性分析。Abstract: Because the selected values of input parameters of generalized Pareto distribution （GPD） model are difficult to avoid uncertainty, the input parameters uncertainty of this model may lead to uncertainty in the seismic hazard estimation. In this paper, we selected northeastern Tibetan Plateau as the case studied area, and proposed an uncertainty analysis process and method of seismic hazard estimation based on the global sensitivity analysis. First, we used the GPD seismicity model to obtain the results of seismic hazard estimation. And then, we selected starting time of earthquake catalog and magnitude threshold to be the input parameters of seismicity model. The E-FAST method with global sensitivity analysis function was used to quantitatively analyze the influence of the uncertainties of the two parameters and the interaction between the two parameters on the uncertainty of seismic hazard estimation. The results show that the seismic hazard estimation of the GPD model is more sensitive to the magnitude threshold. With different return periods, the sensitivity degree of seismic hazard estimation to magnitude threshold is different. For different return periods, there are nonlinear effects between the two input parameters on the uncertainty of seismic hazard estimation, and the degree of nonlinear effects is different. The uncertainty analysis process and method proposed in this paper can be applied to the uncertainty analysis of seismic hazard estimation based on other seismicity models.
-
引言
地震危险性估计中，通常用地震活动性模型描述潜源区地震的时空分布、强度分布、频度分布等地震活动特征（蒋溥，戴丽思，1993）。已有的分析方法可分为确定性方法和概率性方法两个类别，概率性方法是目前构建地震活动性模型中经常使用的分析方法，该方法主要基于研究区的历史地震活动性构建模型，其中最常用的是利用经典的古登堡-里克特的震级频度关系（G-R关系）中对数线性关系的外推。然而对于震级上限而言，传统的G-R关系模型因为线性外推震级上限无限大，与实际情况不符。任雪梅等（2012a）对我国分区域进行了G-R关系的修正和震级极限值的确定，结果表明我国大陆存在非线性的震级-频度关系，修正的G-R关系模型运用于我国中强以上地震震级-频度关系的拟合更为适合。
利用极值理论分析地震活动性也属于概率性方法中常用的一种。近年来，一些学者探讨了依据广义极值理论基于广义帕累托分布（generalized Pareto distribution，缩写为GPD）构建地震活动性模型的方法，用以研究震级分布的尾部特征。Pisarenko和Sornette （2003）分析了哈佛CMT目录中1977—2000年12个地震带的浅层地震，发现广义帕累托分布可以较好地拟合震级分布的尾部特征；钱小仕等（2013a，b）利用广义帕累托分布分析了云南地区历史地震资料以及中国大陆各活动地块边界带的强震震级分布特征，提出了一种利用强震数据推断最大震级分布的可能途径；田建伟等（2017）利用广义帕累托分布估计了马尼拉海沟俯冲带潜在地震海啸源区的地震危险性；任梦依（2018）利用地震活动性广义帕累托模型估计了龙门山地区的震级上限；Dutfoy （2019）将广义帕累托分布与泊松分布相结合，分析了法国南部地区地震震级分布的尾部特征并计算了年最大地震震级。
构建地震活动性广义帕累托模型，利用的是历史地震记录中大震级区段的记录数据，假设大震级区段的震级分布符合广义帕累托分布。这种方法适合历史地震记录时间长、中小地震记录缺失的地区。然而，构建地震活动性广义帕累托模型，需要人为选取历史地震记录起始时间和震级阈值，这些模型输入参量选取的不确定性可能导致地震危险性估计结果的不确定性。
对地震危险性估计的不确定性分析，一直为国内外地震学家高度重视，也开展了一些相关研究。McGuire和Shedlock （1981）及McGuire （1987）是较早对地震危险性分析中模型输入参数和模型输出的不确定性进行研究的学者之一，在对美国中部及东部地区进行概率地震危险性评价时，就知识不完备导致的认识不确定性，采用了逻辑树的方法进行表达和处理；对于地震危险性分析模型的输入参数的不确定性满足连续分布（如正态分布）的情况，Lawrence Livermore National Laboratory （1989）及Cramer等（1996）提出了一种蒙特卡洛抽样结合逻辑树的不确定性表示方法；胡聿贤和鹿林（1990）在关于华北地区地震危险性研究中，分析了潜源区划分的不确定性及地震活动性参数不确定性对地震危险性分析结果的影响；王健和高孟潭（1996）强调，除了要研究单个参数的敏感性外，还必须弄清参数间的交互作用以及多个参数同时变化时结果的范围及其分布，他们利用因子设计中的正交设计法，分析了地震年平均发生率、b值、地震带震级上限、潜在震源区震级上限与地震空间分布函数等几个参数之间的交互作用，以及各参数取不同权重时对烈度概率分布的影响。美国公布的国家地震危险性图，充分考虑了震级上限的不确定性，并将其列为2008年国家地震危险性图编制工作的重大改进（Petersen et al，2008）。我国第五代地震区划图在构建地震活动性模型时，也强调了需充分认识地震活动性的不确定性，并加以综合分析和处理（潘华，李金臣，2016）。对于近些年应用的广义帕累托分布模型，Dutfoy （2019）利用自举法（bootstrap）得到了地震活动性广义帕累托模型参数的分布特征以及年度最大震级分布特征，分析讨论了形状参数概率密度分布的双峰特征导致地震危险性估计结果的不确定性。
然而，已有的相关研究大都仅应用局域敏感性分析方法（local sensitivity analysis method），只能分析当单一输入参量变化而其余参量保持不变时，基于地震活动性模型所得到的地震危险性估计结果的不确定性，并未系统地分析输入参数同时变化以及输入参数之间存在相互作用时，地震危险性估计结果的不确定性，即缺少全域敏感性分析（global sensitivity analysis method）的研究。
因此，本文提出基于全域敏感性分析的地震危险性估计的不确定性分析的流程和方法。以青藏高原东北缘为示例研究区，应用广义帕累托模型，估计研究区的强震重现水平、震级上限及其相应的置信区间，并将全域敏感性分析方法引入地震危险性估计的不确定性分析，利用拓展傅里叶幅度敏感性检验法（extended fourier amplitude sensitivity test，缩写为E-FAST）开展全域敏感性的定量分析。
1. 方法
1.1 地震活动性广义帕累托模型构建
设X1， X2， ···， Xn是一个独立且同分布的随机变量序列，具有边际分布函数，令Mn＝max{X1，···，Xn}，若存在{an＞0，bn∈R}和非退化分布函数G（x），使
$$ {P}_{r}\left\{ （ {M}_{n}－{b}_{n} ） ∕{a}_{n}{\text{≤}} x\right\}\to G （ x ） \text{，} n\to \infty \text{，} $$ (1) 则称G为极值分布。Fisher和Tippett （1928）证明了极值分布的三大类型定理，指出G必属于耿贝尔（Gumbel）分布、弗雷歇（Fréchet）分布、威布尔（Weibull）分布三种类型之一。广义极值理论将这三种分布类型统一为一个分布类型，由位置参数μ，尺度参数σ，形状参数ξ三个参数来表示（Coles，2001；史道济，2006），即
$$ G （ x ） ＝{{\rm{exp}}}\left\{-{\left[1＋\xi \left(\frac{x－\mu }{\sigma }\right)\right]}^{-1/\xi }\right\} {\text{.}} $$ (2) 对于足够大的u，当X＞u，y＝（X－u）的条件分布函数可以近似服从广义帕累托分布，即
$$ H （ y ） ＝1－\left(1＋\frac{\xi y}{\tilde{\sigma }}\right)^{-1/\xi } \text{，} $$ (3) 定义域为
$\{y：y ＞ 0， （ 1＋\xi y/\tilde{\sigma } ） ＞ 0\}$ ，其中$\tilde{\sigma }＝\sigma ＋\xi （ u－\mu ） $ 。在构建地震活动性广义帕累托模型时，首先需确定震级阈值u。阈值u的选取通过两种方法：① 依据震级平均超出量分布函数图，震级平均超出量与震级阈值应满足线性关系；② 随着阈值不断增大，形状参数
$ \xi $ 和修饰的尺度参数$ \tilde{\sigma } $ 估计值应基本保持稳定（Coles，2001）。N年重现水平为预期每N年至少会发生一次震级大于xN的地震事件，如果每年地震发生次数为ny，当ξ≠0时，N年重现水平为
$$ {x}_{N}＝u＋\frac{\tilde{\sigma }}{\xi } ［ （ N{n}_{y}{\zeta }_{u} ） ^{\xi }－1 ］ \text{，} $$ (4) 式中，u为震级阈值，ζu为超过震级阈值的地震样本数与地震样本总数的比值。
考虑到潜在震源区的地震事件的震级必为有限值，因此只有当形状参数ξ＜0，广义帕累托分布具有有限上界，才符合震级存在上限的条件，震级上限估计值可表示为
$$ \hat{x}＝u－\frac{\hat{\tilde{\sigma }}}{\hat{\xi }} {\text{.}} $$ (5) 利用最大似然估计方法可得到三个参数
${\hat{\zeta }}_{u}，{\hat{\tilde{\sigma }}}，\hat{\xi }$ 的估计值。根据Delta定理可知，三个参数的最大似然估计值近似符合正态分布，因此根据方差协方差矩阵可以获得置信水平为1－α的震级上限的置信区间为$$ \hat{x}{\text{±}} {Z}_{1－\frac{\alpha }{2}}\sqrt{{\rm{Var}} （ \hat{x} ） } \;\text{，} $$ (6) 式中
$\text{：}{\rm{Var}} （ \hat{{\boldsymbol{x}}} ） \sim\nabla {\hat{{\boldsymbol{x}}}}^{{\rm{T}}}{\boldsymbol{V}}\nabla \hat{x}$ ，$ V $ 是$ （ {\hat{\zeta }}_{u}，{\hat{\tilde{\sigma }}}，\hat{\xi } ） $ 参数估计的方差协方差矩阵，$ \nabla {\hat{{\boldsymbol{x}}}}^{{\rm{T}}} $ 为$ \dfrac{\partial x}{\partial {\zeta }_{u}}，\dfrac{\partial x} {\partial \tilde{\sigma }}，\dfrac{\partial x}{\partial \xi } $ 在$ （ {\hat{\zeta }}_{u}，{\hat{\tilde{\sigma }}}，\hat{\xi } ） $ 处的值；$ {Z}_{1－\frac{\alpha }{2}} $ 为标准正态分布1－α/2的分位数（Coles，2001）。1.2 定量的全域敏感性分析方法—E-FAST法
E-FAST法是Saltelli等（1999）结合Sobol法（Sobol，1993）的优点，对FAST （Fourier amplitude sensitivity test）法（Cukier et al，1973，1978）进行改进后得到的一种基于方差的全域敏感性定量分析方法（洪明理等，2014；宋明丹等，2014）。该方法将模型的敏感性概括为单个输入参数的敏感性及参数间交互作用的敏感性。模型输出y的总方差V（y）可以分解为
$$ V （ y ） ＝\sum\limits _{i＝1}^{n}{V}_{i}＋\sum\limits _{i\ne j}{V}_{i， j}＋\sum\limits _{i{\text{≠}} j{\text{≠}} k}{V}_{i， j， k}＋\cdots ＋ {V}_{\mathrm{1， 2}， \cdots ， n} \text{，} $$ (7) 式中，Vi＝V ［ E（y|xi） ］ ，Vi，j＝V ［ E（y|xi，xj） ］ －Vi－Vj，Vi，j，k~V1，2，···，n，以此类推。
定义xi的一阶敏感性指标或主效应指标为
$$ {S}_{{ x}_{i}}＝\frac{{V}_{i}}{V （ y ） } \text{，} $$ (8) 式中，Sxi表示单个参数独立作用的敏感性，依据Sxi的大小可以对输入参数独立作用的敏感性进行排序。
定义xi，xj的二阶敏感性指标为
$$ {S}_{{x}_{i}， {x}_{j}}＝\frac{{V}_{i， j}}{V （ y ） } \text{，} $$ (9) 式中，Sxi，xj表示输入参数xi与xj之间的交互作用对模型输出结果不确定性的影响，更高阶的敏感性指标的含义以此类推。
定义xi的全效应指标为
$$ {S}_{{x}_{i}}^{{\rm{T}}}＝{S}_{{x}_{i}}＋\sum\limits_{j{\text{≠}} i}{S}_{{x}_{i}， {x}_{j}}＋\sum\limits_{j{\text{≠}} k{\text{≠}} i}{S}_{{x}_{i}， {x}_{j}， {x}_{k}}＋\cdots ＋{S}_{{x}_{1}， {x}_{2}， \cdots ， {x}_{n}} \text{，} $$ (10) 式中，
${S}_{{x}_{i}}^{{\rm{T}}}$ 表示输入参数xi对模型输出结果的总体影响，包括单个参数的敏感性和该参数与其它参数间交互作用的敏感性。因此可以通过主效应指标与全效应指标的差异，来判断输入参数与其它输入参数是否有交互作用。E-FAST方法的优点是：① 适用于各类模型，对于模型本身的性质没有特殊的限定，可用于分析多个输入参数同时变化时，模型输出结果的不确定性；② 可用于定量地分析各个输入参量的不确定性，以及各个输入参量不确定性的相互作用对模型输出结果不确定性的影响。
2. 案例研究区—青藏高原东北缘
2.1 地震地质背景
青藏高原东北缘位于青藏高原与华北克拉通的交会地带，西南一侧是不断向北东推挤并不断增厚的青藏高原，北东一侧是相对稳定的鄂尔多斯地块，北西一侧是相对稳定的阿拉善地块。在印度板块的推动作用下，青藏高原整体向北东一侧运动，东北缘地区就处在北东向扩展的最前沿。
青藏高原东北缘的构造多形成于印度板块与欧亚大陆碰撞引起的SW-NE向挤压，该挤压始于50 Ma前（Yin，Harrison，2000；Taylor，Yin，2009）。全球定位系统（GPS）的测量结果显示，青藏高原东北缘约以每年10 mm的速率沿NE-SW向或NNE-SSW向水平缩短（Gan et al，2007）。长期以来，青藏高原东北缘构造变形活动强烈、大震多发（Li et al，2017）。该地区分布着几组具有不同活动性质的断裂带，包括左旋走滑断裂带，右旋走滑断裂带，端部逆冲断裂带（徐化超，2018）。这些大型活动断裂对该区地壳稳定性起着重要的调节作用，而在彼此的共同作用下，该区新构造运动强烈，一直以来地震活动十分活跃。
历史上青藏高原东北缘发生过M≥8.0地震5次，分别为1654年天水M8.0地震、1739年平罗M8.0地震、1879年武都M8.0地震、1920年海原M8.5地震和1927年古浪M8.0地震（图1）。本文所用地震目录主要源自国家地震科学数据共享中心的中国历史地震目录（顾功叙，1983）、国家地震台网地震目录、中国台网正式地震目录以及国家地震局震害防御司编制的中国历史强震目录（1995），共收集了青藏高原东北缘（95°E—107°E，32°N—40°N）1880—2020年M≥4.0历史地震1 222次，其中M≥8.0地震3次，7.0≤M＜8.0地震13次，6.0≤M＜7.0地震54次，5.0≤M＜6.0地震228次。黄玮琼等（1994）对中国大陆大部分地区的历史地震资料作了较详细地分析与研究，显示青藏高原东北缘M≥5.0地震资料在1885年之后基本完整。
![]()
图 1 青藏高原东北缘主要活动断裂及M≥4.0地震分布图Figure 1. Map of main active faults and M≥4.0 earthquakes distribution in northeastern Tibetan Plateau为了使震级样本尽量满足构建广义帕累托模型时独立同分布的条件，在进行地震危险性估计之前需要尽量删除地震目录中的余震。目前国内外地震危险性分析研究中，考虑到计算的简便和工程应用上的可操作性，普遍采用Gardner和Knopoff （1974）提出的时空窗法来删除余震（徐伟进，2012）。空间窗口的确定方法为（Knopoff，Gardner，1972）
$$ \mathrm{l}\mathrm{g}R＝aM＋b \text{，} $$ (11) 式中，R为窗口半径，a，b为固定参数值，在我国a，b的取值分别为0.5和−1.78 （刘杰等，1996）。
最终确定的时间窗口如表1所示。删除余震后，共收集到青藏高原东北缘1880—2020年M≥4.0历史地震864次，M≥5.0地震的M-t图如图2所示。
表 1 不同震级地震的余震时间窗（Gardner，Knopoff，1974）Table 1. Aftershock time windows of different magnitudes（Gardner，Knopoff，1974）主震震级 持续天数/d 主震震级 持续天数/d M4.0 42 M6.5 790 M4.5 83 M7.0 915 M5.0 155 M7.5 960 M5.5 290 M8.0 985 M6.0 510 M8.5 985 ![]()
图 2 青藏高原东北缘1880—2020年M≥5.0地震的M-t图Figure 2. M-t map of the M≥5.0 earthquakes from 1880 to 2020 in northeastern Tibetan Plateau2.2 地震危险性估计
本文分析了1885—2020年青藏高原东北缘的地震目录，利用广义帕累托模型估算了青藏高原东北缘震级重现水平和震级上限。图3为青藏高原东北缘震级数据的平均超出量分布图，置信区间为95%。在平均超出量函数具有线性特征的区段选取对应的阈值。当u大于5.5时，平均超出量函数近似呈线性。另外，还要综合考虑样本数据的实际情况，通过观察选取不同阈值时，形状参数和修饰的尺度参数估计值的变化情况（图4）来进行取舍。应选取形状参数和修饰的尺度参数较为稳定呈近似水平直线的区段所对应的阈值。由图4容易看出，当阈值大于5.5时，两个参数的估计值随阈值的增大相对比较稳定。另外需要考虑的是，阈值取得过大，超阈值的样本量较少，会导致估计量的方差过大，反之阈值取得太小，使得超出量分布与广义帕累托分布出现较大偏差，会导致估计量的标准差过大。综合考虑超出量分布函数的线性要求和分布参数估计的稳定性要求，选取阈值u＝5.5。
![]()
图 4 修饰的尺度参数（a）和形状参数（b）估计值随震级阈值选取的变化Figure 4. Variation diagram of estimators for modified scale parameter （a） and shape parameter （b） according to the magnitude threshold地震活动性广义帕累托模型，可以充分利用历史地震记录中大震级区段的信息，为研究震级分布的尾部特征提供统计分析方法。本文模型设置震级阈值为5.5，因而，如果所采用的地震目录小于M5.5的地震记录缺失，将不会对模型结果产生影响，该模型的特点也正是适合历史地震记录时间长、中小地震记录缺失的地区。
确定震级阈值后，得到研究区震级的超出量样本。依据式（3）进行广义帕累托分布的拟合，利用极大似然估计方法，得到修饰的尺度参数和形状参数的估计值分别为1.12和−0.32。由于形状参数的估计值为负值，因此广义帕累托分布为威布尔分布，右端点为有限值，意味着研究区存在震级上限。将阈值、形状参数和修饰的尺度参数估计值分别代入式（4）和式（5）进行研究区地震危险性估计，结果列于表2。
表 2 青藏高原东北缘地震危险性估计Table 2. Seismic hazard estimation for northeastern Tibetan Plateau重现期/a M 置信度95%的置信区间 20 7.40 ［ 7.10，7.69 ］ 50 7.80 ［ 7.49，8.11 ］ 100 8.03 ［ 7.69，8.37 ］ 200 8.22 ［ 7.84，8.59 ］ 500 8.41 ［ 7.97，8.84 ］ ∞ 8.95 ［ 8.03，9.87 ］ 对广义帕累托分布拟合的情况进行诊断，通常根据实际数据与分布模型的一致性来判断模型的准确性。诊断结果由4个图（概率图、分位数图、重现水平图、密度曲线图）表示（图5）。概率图根据样本数据的累积概率与模型分布的累积概率之间的关系绘制，分位数图依据样本数据分布的分位数与模型分布的分位数之间的关系绘制。当样本数据符合理论假设的分布时，概率图（图5a）和分位数图（图5b）应近似为直线。图5c代表的是对应于不同重现期的重现水平，样本数据应落在给定分布分位数估计置信区间内。由于ξ＜0，为威布尔分布，符合函数存在上限值的条件，因此重现水平曲线为凸曲线，且逐渐趋于有限值。图5d中拟合的密度曲线与数据的直方图也较为一致。综上，从图5可以看出，各散点数据基本紧密围绕各参考线分布，表明拟合状态良好，上述诊断结果不能否定利用广义帕累托分布分析青藏高原东北缘强震震级分布特征的合理性和适用性。
![]()
图 5 广义帕累托分布拟合诊断图（a） 概率图；（b） 分位数图；（c） 重现水平图；（d） 概率密度图Figure 5. Diagnostic plots of the generalized Pareto distribution fitted to magnitude（a） Probability plot；（b） Quantile plot；（c） Return level plot；（d） Probability density plot研究区震级上限估计值为8.95，置信度95%的置信区间为 ［ 8.03，9.87 ］ ，而研究区记录到的历史最大地震震级为8.5。
震级上限的含义，是指“在地震带震级-频度关系式中，累积频度趋于零的震级极限值”（胡聿贤，1999）。地震带的震级上限与带内历史上最大地震震级、地震带的地震活动水平、b值以及该带最大地震的重现期有关，即
$$ {M}＝{M}_{T}＋\frac{1}{D}\lg\left(\frac{1－{10}^{-A}}{1－{10}^{-A＋b{M}_{T}}}\right) \text{，} $$ (12) 式中：MT为T年内累积频度为1的地震震级；A＝lg（v0T），v0为M≥4.0地震年平均发生率；T为所考虑时间段。从式中可以看到，
$\left({1－{10}^{-A}}/{1－{10}^{-A＋b{M}_{T}}}\right)$ 大于1，则M＞MT，也就是说，地震带的震级上限应大于该带历史最大地震震级。当然，由于构造规模和岩石强度都有极限，实际上M8.6地震已为数极少（胡聿贤，1999）。任雪梅等（2012b）基于修正的G-R关系模型，对祁连山—六盘山地震带1970年以来ML3.0以上地震的震级-频度关系进行拟合，得到最大震级为9.66。钱小仕等（2013b）基于广义帕累托分布的超阈值分布模型，对中国大陆活动地块边界带强震震级分布特征开展研究，以1921年为起始年，震级阈值为5.3，估计了海源—祁连带震级上限为9.95。
对于广义帕累托模型，由于当ξ逐渐趋近于0时，模型估计的震级上限值可能会变得不稳定，Pisarenko等（2008，2014）提出可以利用分位数来进一步量化尾部特征。钱小仕等（2013a，b）也提到在工程地震危险性分析中，潜在最大震级的确定往往与特定建筑结构的抗震设计年限相关，他们利用0.999 97高分位数作为最大震级估计，其值往往低于广义帕累托模型公式直接计算的震级上限估计值。本文旨在提出基于全域敏感性分析的地震危险性估计不确定性分析流程和方法，震级上限估计值未采用分位数的形式。
基于现有的地震活动性模型，进行地震危险性估计，都不可避免地存在不确定性。
2.3 不确定性分析
构建地震活动性广义帕累托模型，需要人为选取历史地震记录起始时间和震级阈值。这些输入参量选取的不确定性可能导致地震危险性估计结果的不确定性。本文采用E-FAST方法开展定量的全域敏感性分析。
地震活动性广义帕累托模型的输入参数为地震目录起始时间ts和震级阈值u，输出结果为不同时期的震级重现水平、震级上限及其对应的置信区间。青藏高原东北缘地震目录起始时间范围的选取是基于对该地区地震目录完整性的研究（黄玮琼等，1994）；震级阈值的选择范围通过上述1.1节中阈值的两种选择方法确定。设地震目录起始时间范围为1880—1890年，震级阈值的取值范围介于5.5—5.9之间，并假设它们在其不确定性范围内服从均匀分布，利用E-FAST方法进行蒙特卡洛抽样，共采集输入参数的抽样样本194组。
对青藏高原东北缘地震活动性广义帕累托模型进行参数敏感性分析。首先，根据输入参数的样本组合，计算相应的震级上限，以及对应一定重现期的震级重现水平；其次，将输入参数和相应的地震危险性分析结果文件导入敏感性分析软件中，利用E-FAST方法，计算广义帕累托模型地震危险性估计结果对各输入参数的敏感性指标，即ts和u的主效应指标Si和全效应指标
$S^{\rm{T}}_i $ （表3）。为了更直观地显示各输入参数的敏感性指标，绘制如图6—8所示的直方图。表 3 地震危险性估计结果对地震目录起始时间ts和震级阈值u的主效应指标Si和全效应指标$S^{\rm{T}}_i$ Table 3. Total and the first-order effects of the seismic hazard estimation on earthquake catalogue initial time ts and magnitude threshold u重现期/a 参数 M 置信度95%的置信区间下端点 置信度95%的置信区间上端点 Si $S^{\rm{T}}_i$ Si $S^{\rm{T}}_i$ Si $S^{\rm{T}}_i$ 20 u 0.742 4 0.937 8 0.779 1 0.968 4 0.681 3 0.880 5 ts 0.059 5 0.184 0 0.030 2 0.158 9 0.112 5 0.229 4 50 u 0.749 5 0.868 0 0.740 6 0.982 8 0.402 3 0.797 9 ts 0.117 8 0.178 2 0.015 3 0.182 4 0.138 9 0.502 4 100 u 0.374 1 0.678 8 0.709 5 0.987 7 0.570 5 0.964 7 ts 0.228 6 0.522 3 0.008 8 0.204 5 0.011 8 0.329 7 200 u 0.487 6 0.936 6 0.685 8 0.988 9 0.600 7 0.978 9 ts 0.024 9 0.405 1 0.005 6 0.223 1 0.004 3 0.300 2 500 u 0.568 6 0.975 0 0.661 6 0.987 6 0.609 4 0.982 2 ts 0.004 8 0.329 7 0.004 0 0.244 3 0.003 3 0.292 1 ∞ u 0.460 0 0.958 8 0.345 0 0.929 2 0.386 0 0.940 2 ts 0.015 2 0.473 8 0.029 5 0.623 5 0.024 2 0.573 0 ![]()
图 6 地震危险性估计结果对地震目录起始时间ts和震级阈值u的主效应指标Si直方图横轴数字代表的含义下同Figure 6. Histogram of the first-order effect Si of seismic hazard estimation on earthquake catalogue start time ts and magnitude threshold uThe meanings of the number represent are the same below![]()
图 8 地震危险性估计结果对地震目录起始时间ts和震级阈值u的主效应指标和交互效应指标直方图Figure 8. Histogram of the first-order effect and interactions on earthquake catalogue start time ts and magnitude threshold u图6为广义帕累托模型的地震危险性估计结果对两个输入参数的主效应指标直方图。从图上可以看出，广义帕累托模型估计结果明显对u更为敏感。对于ts，除了100年震级重现水平对应的主效应指标略高于0.2，其余的主效应指标均不到0.2。u的主效应指标远高于ts的主效应指标，因此在建立模型时，首先要考虑震级阈值u选取的合理性和准确性。
结合表3综合分析，对应不同重现期的u的主效应指标有所差别，表明震级重现水平及震级上限对u的敏感性程度也有所不同。除了100年震级重现水平对应的u的主效应指标为0.374 1，其余重现期为20，50，200和500年的震级重现水平以及震级上限对应的u的主效应指标均高于0.45，分别为0.742 4，0.749 5，0.487 6，0.568 6和0.460 0。
图7为广义帕累托模型的地震危险性估计结果对两个输入参数的全效应指标直方图。主效应与全效应的大小差异，即交互效应，体现了输入参数之间对输出结果不确定性的影响存在非线性效应与否。图8主要体现了广义帕累托模型的地震危险性估计结果对两个输入参数的交互效应指标。图中柱状长度对应的数值代表了输入参数的全效应指标大小，其中，蓝色部分代表主效应指标在全效应指标中的占比，橘色部分代表交互效应指标在全效应指标中的占比。
![]()
图 7 地震危险性估计结果对地震目录起始时间ts和震级阈值u的全效应指标$S^{\rm{T}}_i$ 直方图Figure 7. Histogram of the total effect$S^{\rm{T}}_i$ of seismic hazard estimation on earthquake catalogue start time ts and magnitude threshold u从图7和8可以看出，两个参数对地震危险性估计结果不确定性的影响均存在非线性效应，即两个参数之间存在相互作用。对应不同的重现期，在影响地震危险性估计结果不确定性上，两个参数之间的相互作用程度不同。对比不同重现期的震级重现水平，输入参数u的全效应指标中，主效应占主要部分，交互效应的比重随重现期的增大而逐渐增大；而另一输入参数ts的全效应指标中，交互效应占主要部分。对于震级上限而言，两个输入参数的全效应指标中，交互效应占主要部分，表明这两个参数主要是以交互作用的形式对震级上限的不确定性产生影响。
3. 讨论与结论
本文选取青藏高原东北缘为案例研究区，提出了基于全域敏感性分析的地震危险性估计不确定性分析流程和方法。
基于广义帕累托分布，构建了研究区地震活动性模型，进行了地震危险性估计（包括震级重现水平和震级上限的估计）。鉴于构建模型需要人为选取ts和u这两个输入参数，其取值的不确定性，可能导致地震危险性估计结果的不确定性，因此，选取具有全域敏感性分析功能的E-FAST方法，定量分析了上述两个输入参数的不确定性对地震危险性估计结果不确定性的影响。结果表明：地震危险性估计结果，对两个输入参数中的u更为敏感；对应不同重现期，地震危险性估计结果对u的敏感性程度有所不同；两个输入参数对地震危险性估计结果不确定性的影响，存在非线性效应；对不同的重现期而言，在影响地震危险性估计结果不确定性上，两个输入参数之间的相互作用程度不同。
目前，确定某一构造或构造区的最大震级，常用的方法有概率统计方法和确定性方法。前者有基于G-R关系统计分析的方法或利用震级-断层破裂尺度之间经验关系统计分析的方法；后者有基于研究区历史最大地震估计的方法和构造类比法等。新近的研究，开始倾向于以形变观测等资料作为约束，结合构造或区域的历史地震信息，估计该区域的最大震级。
本文选用基于广义帕累托分布构建的地震活动性模型，进行地震危险性估计，所用的方法属于概率统计方法。这一方法，利用历史地震记录中大震级区段的信息，无需先验地选定震级下限，较为适合历史地震记录时间长但低震级地震记录缺失的地区，便于构建研究区强震活动性模型，和其它概率统计方法一样，也有其局限性，即未能考虑地震发生的物理机制和物理过程，因为仅利用历史地震记录中大震级区段的信息，有时统计所用的历史地震记录样本偏少。
必须指出，不论通过哪一种方法构建地震活动性模型进行地震危险性估计，都难以避免不确定性，因此有必要开展相应的不确定性和敏感性分析研究，分析地震活动性模型不确定性的来源和特征，即分析输入参数的不确定性如何影响和制约地震危险性估计的不确定性。这些分析结果，可以指导地震活动性模型构建时输入参数的选取和调整，从而使地震活动性模型构建得更为完善。
本文提出的定量全域敏感性分析方法，对地震活动性模型的性质没有特殊的限定，便于在地震危险性估计不确定性分析中推广应用。
-
![]()
图 1 青藏高原东北缘主要活动断裂及M≥4.0地震分布图
Figure 1. Map of main active faults and M≥4.0 earthquakes distribution in northeastern Tibetan Plateau
![]()
图 2 青藏高原东北缘1880—2020年M≥5.0地震的M-t图
Figure 2. M-t map of the M≥5.0 earthquakes from 1880 to 2020 in northeastern Tibetan Plateau
![]()
图 4 修饰的尺度参数（a）和形状参数（b）估计值随震级阈值选取的变化
Figure 4. Variation diagram of estimators for modified scale parameter （a） and shape parameter （b） according to the magnitude threshold
![]()
图 5 广义帕累托分布拟合诊断图
（a） 概率图；（b） 分位数图；（c） 重现水平图；（d） 概率密度图
Figure 5. Diagnostic plots of the generalized Pareto distribution fitted to magnitude
（a） Probability plot；（b） Quantile plot；（c） Return level plot；（d） Probability density plot
![]()
图 6 地震危险性估计结果对地震目录起始时间ts和震级阈值u的主效应指标Si直方图
横轴数字代表的含义下同
Figure 6. Histogram of the first-order effect Si of seismic hazard estimation on earthquake catalogue start time ts and magnitude threshold u
The meanings of the number represent are the same below
![]()
图 8 地震危险性估计结果对地震目录起始时间ts和震级阈值u的主效应指标和交互效应指标直方图
Figure 8. Histogram of the first-order effect and interactions on earthquake catalogue start time ts and magnitude threshold u
![]()
图 7 地震危险性估计结果对地震目录起始时间ts和震级阈值u的全效应指标
$S^{\rm{T}}_i$ 直方图Figure 7. Histogram of the total effect
$S^{\rm{T}}_i$ of seismic hazard estimation on earthquake catalogue start time ts and magnitude threshold u表 1 不同震级地震的余震时间窗（Gardner，Knopoff，1974）
Table 1 Aftershock time windows of different magnitudes（Gardner，Knopoff，1974）
主震震级 持续天数/d 主震震级 持续天数/d M4.0 42 M6.5 790 M4.5 83 M7.0 915 M5.0 155 M7.5 960 M5.5 290 M8.0 985 M6.0 510 M8.5 985 
下载: 导出CSV
表 2 青藏高原东北缘地震危险性估计
Table 2 Seismic hazard estimation for northeastern Tibetan Plateau
重现期/a M 置信度95%的置信区间 20 7.40 ［ 7.10，7.69 ］ 50 7.80 ［ 7.49，8.11 ］ 100 8.03 ［ 7.69，8.37 ］ 200 8.22 ［ 7.84，8.59 ］ 500 8.41 ［ 7.97，8.84 ］ ∞ 8.95 ［ 8.03，9.87 ］
下载: 导出CSV
表 3 地震危险性估计结果对地震目录起始时间ts和震级阈值u的主效应指标Si和全效应指标
$S^{\rm{T}}_i$ Table 3 Total and the first-order effects of the seismic hazard estimation on earthquake catalogue initial time ts and magnitude threshold u
重现期/a 参数 M 置信度95%的置信区间下端点 置信度95%的置信区间上端点 Si $S^{\rm{T}}_i$ Si $S^{\rm{T}}_i$ Si $S^{\rm{T}}_i$ 20 u 0.742 4 0.937 8 0.779 1 0.968 4 0.681 3 0.880 5 ts 0.059 5 0.184 0 0.030 2 0.158 9 0.112 5 0.229 4 50 u 0.749 5 0.868 0 0.740 6 0.982 8 0.402 3 0.797 9 ts 0.117 8 0.178 2 0.015 3 0.182 4 0.138 9 0.502 4 100 u 0.374 1 0.678 8 0.709 5 0.987 7 0.570 5 0.964 7 ts 0.228 6 0.522 3 0.008 8 0.204 5 0.011 8 0.329 7 200 u 0.487 6 0.936 6 0.685 8 0.988 9 0.600 7 0.978 9 ts 0.024 9 0.405 1 0.005 6 0.223 1 0.004 3 0.300 2 500 u 0.568 6 0.975 0 0.661 6 0.987 6 0.609 4 0.982 2 ts 0.004 8 0.329 7 0.004 0 0.244 3 0.003 3 0.292 1 ∞ u 0.460 0 0.958 8 0.345 0 0.929 2 0.386 0 0.940 2 ts 0.015 2 0.473 8 0.029 5 0.623 5 0.024 2 0.573 0
下载: 导出CSV
-
顾功叙. 1983. 中国地震目录: 公元前1831—公元1969年[M]. 北京: 科学出版社: 773–791. Gu G X. 1983. Catalogue of Earthquakes in China: 1831BC-1969AD[M]. Beijing: Science Press: 773–791 （in Chinese）.
国家地震局震害防御司. 1995. 中国历史强震目录[M]. 北京: 地震出版社: 496–499. Department of Earthquake Damage Prevention, State Seismological Bureau. 1995. Catalogue of Historical Strong Earthquakes in China[M]. Beijing: Seismological Press: 496–499 （in Chinese）.
洪明理,任鲁川,霍振香. 2014. 基于E-FAST法分析海啸波高对潜在海啸源参数的敏感性[J]. 地震学报,36(2):252–260. Hong M L,Ren L C,Huo Z X. 2014. Sensitivity analysis on maximum tsunami wave heights to the potential tsunami source parameters based on extended FAST method[J]. Acta Seismologica Sinica,36(2):252–260 (in Chinese).
胡聿贤, 鹿林. 1990. 地震活动性估计的不确定性[G]//地震危险性分析中的综合概率法. 北京: 地震出版社: 176–177. Hu Y X, Lu L. 1990. Uncertainty in seismicity estimation[G]//Synthetic Probability Method in Seismic Hazard Analysis. Beijing: Seismological Press: 176–177 （in Chinese）.
胡聿贤. 1999. 地震安全性评价技术教程[M]. 北京: 地震出版社: 227–228. Hu Y X. 1999. Seismic Safety Evaluation Technology Tutorial[M]. Beijing: Seismological Press: 227–228 （in Chinese）.
黄玮琼,李文香,曹学锋. 1994. 中国大陆地震资料完整性研究之二:分区地震资料基本完整的起始年分布图象[J]. 地震学报,16(4):423–432. Huang W Q,Li W X,Cao X F. 1994. Study on the integrity of seismic data in Mainland China Ⅱ:The initial year distribution images of the complete data on each seismic zone[J]. Acta Seismologica Sinica,16(4):423–432 (in Chinese).
蒋溥, 戴丽思. 1993. 工程地震学概论[M]. 北京: 地震出版社: 110. Jiang P, Dai L S. 1993. Introduction to Engineering Seismology[M]. Beijing: Seismological Press: 110 （in Chinese）.
刘杰,陈棋福,陈顒. 1996. 华北地区地震目录完全性分析[J]. 地震,16(1):59–67. Liu J,Chen Q F,Chen Y. 1996. Completeness analysis of the seismic catalog in North China region[J]. Earthquake,16(1):59–67 (in Chinese).
潘华,李金臣. 2016. 新一代地震区划图的地震活动性模型[J]. 城市与减灾,(3):28–33. doi: 10.3969/j.issn.1671-0495.2016.03.008 Pan H,Li J C. 2016. A seismicity model for a new generation of seismic zoning maps[J]. City and Disaster Reduction,(3):28–33 (in Chinese).
钱小仕,王福昌,盛书中. 2013a. 基于广义帕累托分布的地震震级分布尾部特征分析[J]. 地震学报,35(3):341–350. Qian X S,Wang F C,Sheng S Z. 2013a. Characterization of tail distribution of earthquake magnitudes via generalized Pareto distribution[J]. Acta Seismologica Sinica,35(3):341–350 (in Chinese).
钱小仕,蔡晓光,任晴晴. 2013b. 中国大陆活动地块边界带强震震级分布特征研究[J]. 地震工程与工程振动,33(1):212–220. Qian X S,Cai X G,Ren Q Q. 2013b. Characteristics of the great earthquake magnitude distributions for active tectonic boundaries in Chinese mainland[J]. Earthquake Engineering and Engineering Vibration,33(1):212–220 (in Chinese).
任梦依. 2018. 龙门山地区的地震活动性广义帕累托模型构建[J]. 地震研究,41(2):226–232. doi: 10.3969/j.issn.1000-0666.2018.02.010 Ren M Y. 2018. The establishment of generalized Pareto distribution model of seismicity in Longmenshan region[J]. Journal of Seismological Research,41(2):226–232 (in Chinese).
任雪梅,高孟潭,俞言祥. 2012a. 基于MGR模型修正我国大陆中强以上地震的震级-频度关系和确定震级极限值[J]. 地震学报,34(3):331–338. Ren X M,Gao M T,Yu Y X. 2012a. Modification of magnitude-frequency relation and magnitude limit determination based on MGR model for moderate-strong earthquakes in Chinese mainland[J]. Acta Seismologica Sinica,34(3):331–338 (in Chinese).
任雪梅,高孟潭,张纳莉. 2012b. 基于MGR模型的我国大陆地区各地震带1970年以来震级·频度关系和震级上限[J]. 中国地震,28(3):320–327. Ren X M,Gao M T,Zhang N L. 2012b. Magnitude-frequency relation and magnitude limit of seismic zones based on the MGR model in Chinese mainland since 1970[J]. Earthquake Research in China,28(3):320–327 (in Chinese).
史道济. 2006. 实用极值统计方法[M]. 天津: 天津科学技术出版社: 28–32. Shi D J. 2006. Practical Extremum Statistical Methods[M]. Tianjin: Tianjin Science and Technology Press: 28–32 （in Chinese）.
宋明丹,冯浩,李正鹏,高建恩. 2014. 基于Morris和EFAST的CERES-Wheat模型敏感性分析[J]. 农业机械学报,45(10):124–131. doi: 10.6041/j.issn.1000-1298.2014.10.020 Song M D,Feng H,Li Z P,Gao J E. 2014. Global sensitivity analyses of DSSAT-CERES-Wheat model using Morris and EFAST methods[J]. Transactions of the Chinese Society for Agricultural Machinery,45(10):124–131 (in Chinese).
田建伟,刘哲,任鲁川. 2017. 基于广义帕累托分布的马尼拉海沟俯冲带地震危险性估计[J]. 地震,37(1):158–165. Tian J W,Liu Z,Ren L C. 2017. Seismic hazard estimation of the Manila trench subduction zone based on generalized Pareto distribution[J]. Earthquake,37(1):158–165 (in Chinese).
王健,高孟潭. 1996. 地震危险性分析中的参数敏感性研究[J]. 地震学报,18(4):489–493. Wang J,Gao M T. 1996. Study on parameter sensitivity in seismic hazard analysis[J]. Acta Seismologica Sinica,18(4):489–493 (in Chinese).
徐化超. 2018. 青藏高原东北缘地区主要活动断裂带的运动学研究[D]. 北京: 中国地震局地震预测研究所: 15–41. Xu H C. 2008. Kinematics Study of Main Active Fault Zones in the Northeastern Qinghai-Tibet Plateau[D]. Beijing: Institute of Earthquake Forecasting, China Earthquake Administration: 15–41 （in Chinese）.
徐伟进. 2012. 地震危险性分析中地震时空统计分布模型研究[D]. 北京: 中国地震局地球物理研究所: 4–5. Xu W J. 2012. Study on Space-Time Statistical Distribution Models of Earthquakes in Seismic Hazard Analysis[D]. Beijing: Institute of Earthquake Forecasting, China Earthquake Administration: 4–5 （in Chinese）.
Coles S. 2001. An Introduction to Statistical Modeling of Extreme Values[M]. London: Springer-Verlag: 74–91.
Cramer C H,Petersen M D,Reichle M S. 1996. A Monte Carlo approach in estimating uncertainty for a seismic hazard assessment of Los Angeles,Ventura,and Orange counties,California[J]. Bull Seismol Soc Am,86(6):1681–1691. doi: 10.1785/BSSA0860061681
Cukier R I,Fortuin C M,Shuler K E. 1973. Study of the sensitivity of coupled reaction systems to uncertainties in rate coefficients I:Theory[J]. J Chem Phys,59(8):3873–3878. doi: 10.1063/1.1680571
Cukier R I,Levine H B,Shuler K E. 1978. Nonlinear sensitivity analysis of multiparameter model systems[J]. J Comput Phys,26(1):1–42. doi: 10.1016/0021-9991(78)90097-9
Dutfoy A. 2019. Estimation of tail distribution of the annual maximum earthquake magnitude using extreme value theory[J]. Pure Appl Geophys,176(2):527–540. doi: 10.1007/s00024-018-2029-0
Fisher R A,Tippett L H C. 1928. Limiting forms of the frequency distribution of the largest or smallest member of a sample[J]. Math Proc Camb Philos Soc,24(2):180–190. doi: 10.1017/S0305004100015681
Gan W J,Zhang P Z,Shen Z K,Niu Z J,Wang M,Wan Y G,Zhou D M,Cheng J. 2007. Present-day crustal motion within the Tibetan Plateau inferred from GPS measurements[J]. J Geophys Res:Solid Earth,112(B8):B08416.
Gardner J K, Knopoff L. 1974. Is the sequence of earthquakes in southern California, with aftershocks removed, Poissonian?[J] Bull Seismol Soc Am, 64（5）: 1363-1367.
Knopoff L,Gardner J K. 1972. Higher seismic activity during local night on the raw worldwide earthquake catalogue[J]. Geophys J Int,28(3):311–313. doi: 10.1111/j.1365-246X.1972.tb06133.x
Lawrence Livermore National Laboratory （LLNL）. 1989. Seismic hazard characterization of 69 nuclear plant sites East of the Rocky Mountains: Methodology, input data and comparisons to previous results for ten test site[J]. NUREG/CR-5250-V1.
Li Y H,Wang X C,Zhang R Q,Wu Q J,Ding Z F. 2017. Crustal structure across the NE Tibetan Plateau and Ordos block from the joint inversion of receiver functions and Rayleigh-wave dispersions[J]. Tectonophysics,705:33–41. doi: 10.1016/j.tecto.2017.03.020
McGuire R K,Shedlock K M. 1981. Statistical uncertainties in seismic hazard evaluations in the United States[J]. Bull Seismol Soc Am,71(4):1287–1308.
McGuire R K. 1987. Seismic hazard uncertainty and its effects on design earthquake ground motions[C]//Proceeding of International Seminar on Seismic Zonation. Guangzhou: State Seismological Bureau: 351–359.
Petersen M D, Frankel A D, Harmsen S C, Mueller C S, Haller K M, Wheeler R L, Wesson R L, Zeng Y H, Boyd O S, Perkins D M, Luco N, Field E H, Wills C J, Rukstales K S. 2008. Documentation for the 2008 Update of the United States National Seismic Hazard Maps[R]. Reston, Virginia: U.S. Geological Survey: 1–60.
Pisarenko V F,Sornette D. 2003. Characterization of the frequency of extreme earthquake events by the generalized Pareto distribution[J]. Pure Appl Geophys,160(12):2343–2364. doi: 10.1007/s00024-003-2397-x
Pisarenko V F,Sornette A,Sornette D,Rodkin M V. 2014. Characterization of the tail of the distribution of earthquake magnitudes by combining the GEV and GPD descriptions of extreme value theory[J]. Pure Appl Geophys,171(8):1599–1624. doi: 10.1007/s00024-014-0882-z
Saltelli A,Tarantola S,Chan K P S. 1999. A quantitative model-independent method for global sensitivity analysis of model output[J]. Technometrics,41(1):39–56. doi: 10.1080/00401706.1999.10485594
Sobol I M. 1993. Sensitivity estimates for nonlinear mathematical models[J]. Math Model Comput Exp,1(4):407–414.
Taylor M,Yin A. 2009. Active structures of the Himalayan-Tibetan orogen and their relationships to earthquake distribution,contemporary strain field,and Cenozoic volcanism[J]. Geosphere,5(3):199–214. doi: 10.1130/GES00217.1
Yin A,Harrison T M. 2000. Geologic evolution of the Himalayan-Tibetan orogen[J]. Annu Rev Earth Planet Sci,28:211–280. doi: 10.1146/annurev.earth.28.1.211
-
期刊类型引用(0)
其他类型引用(1)
计量
- 文章访问数: 243
- HTML全文浏览量: 118
- PDF下载量: 72
- 被引次数: 1
