引言

地震危险性估计中，通常用地震活动性模型描述潜源区地震的时空分布、强度分布、频度分布等地震活动特征（蒋溥，戴丽思，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   地震活动性广义帕累托模型构建

设X₁， X₂， ···， X_n是一个独立且同分布的随机变量序列，具有边际分布函数，令M_n＝max{X₁，···，X_n}，若存在{a_n＞0，b_n∈R}和非退化分布函数G（x），使

则称G为极值分布。Fisher和Tippett （1928）证明了极值分布的三大类型定理，指出G必属于耿贝尔（Gumbel）分布、弗雷歇（Fréchet）分布、威布尔（Weibull）分布三种类型之一。广义极值理论将这三种分布类型统一为一个分布类型，由位置参数μ，尺度参数σ，形状参数ξ三个参数来表示（Coles，2001；史道济，2006），即

对于足够大的u，当X＞u，y＝（X－u）的条件分布函数可以近似服从广义帕累托分布，即

定义域为$\{y：y ＞ 0， （ 1＋\xi y/\tilde{\sigma } ） ＞ 0\}$，其中$\tilde{\sigma }＝\sigma ＋\xi （ u－\mu ） $。

在构建地震活动性广义帕累托模型时，首先需确定震级阈值u。阈值u的选取通过两种方法：① 依据震级平均超出量分布函数图，震级平均超出量与震级阈值应满足线性关系；② 随着阈值不断增大，形状参数$ \xi $和修饰的尺度参数$ \tilde{\sigma } $估计值应基本保持稳定（Coles，2001）。

N年重现水平为预期每N年至少会发生一次震级大于x_N的地震事件，如果每年地震发生次数为n_y，当ξ≠0时，N年重现水平为

式中，u为震级阈值，ζ_u为超过震级阈值的地震样本数与地震样本总数的比值。

考虑到潜在震源区的地震事件的震级必为有限值，因此只有当形状参数ξ＜0，广义帕累托分布具有有限上界，才符合震级存在上限的条件，震级上限估计值可表示为

利用最大似然估计方法可得到三个参数${\hat{\zeta }}_{u}，{\hat{\tilde{\sigma }}}，\hat{\xi }$的估计值。根据Delta定理可知，三个参数的最大似然估计值近似符合正态分布，因此根据方差协方差矩阵可以获得置信水平为1－α的震级上限的置信区间为

式中$\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_i＝V ［ E（y|x_i） ］ ，V_i，j＝V ［ E（y|x_i，x_j） ］ －V_i－V_j，V_i，j，k~V_{1，2，···，n}，以此类推。

定义x_i的一阶敏感性指标或主效应指标为

式中，S_xi表示单个参数独立作用的敏感性，依据S_xi的大小可以对输入参数独立作用的敏感性进行排序。

定义x_i，x_j的二阶敏感性指标为

式中，S_xi，xj表示输入参数x_i与x_j之间的交互作用对模型输出结果不确定性的影响，更高阶的敏感性指标的含义以此类推。

定义x_i的全效应指标为

式中，${S}_{{x}_{i}}^{{\rm{T}}}$表示输入参数x_i对模型输出结果的总体影响，包括单个参数的敏感性和该参数与其它参数间交互作用的敏感性。因此可以通过主效应指标与全效应指标的差异，来判断输入参数与其它输入参数是否有交互作用。

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）

式中，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

下载: 导出CSV 

| 显示表格

图  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

下载: 全尺寸图片 幻灯片

2.2   地震危险性估计

本文分析了1885—2020年青藏高原东北缘的地震目录，利用广义帕累托模型估算了青藏高原东北缘震级重现水平和震级上限。图3为青藏高原东北缘震级数据的平均超出量分布图，置信区间为95%。在平均超出量函数具有线性特征的区段选取对应的阈值。当u大于5.5时，平均超出量函数近似呈线性。另外，还要综合考虑样本数据的实际情况，通过观察选取不同阈值时，形状参数和修饰的尺度参数估计值的变化情况（图4）来进行取舍。应选取形状参数和修饰的尺度参数较为稳定呈近似水平直线的区段所对应的阈值。由图4容易看出，当阈值大于5.5时，两个参数的估计值随阈值的增大相对比较稳定。另外需要考虑的是，阈值取得过大，超阈值的样本量较少，会导致估计量的方差过大，反之阈值取得太小，使得超出量分布与广义帕累托分布出现较大偏差，会导致估计量的标准差过大。综合考虑超出量分布函数的线性要求和分布参数估计的稳定性要求，选取阈值u＝5.5。

图  3  震级平均超出量分布函数图

Figure  3.  Diagram of mean excess function of magnitude

下载: 全尺寸图片 幻灯片

图  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 ］

下载: 导出CSV 

| 显示表格

对广义帕累托分布拟合的情况进行诊断，通常根据实际数据与分布模型的一致性来判断模型的准确性。诊断结果由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_T为T年内累积频度为1的地震震级；A＝lg（v₀T），v₀为M≥4.0地震年平均发生率；T为所考虑时间段。从式中可以看到，$\left({1－{10}^{-A}}/{1－{10}^{-A＋b{M}_{T}}}\right)$大于1，则M＞M_T，也就是说，地震带的震级上限应大于该带历史最大地震震级。当然，由于构造规模和岩石强度都有极限，实际上M8.6地震已为数极少（胡聿贤，1999）。

任雪梅等（2012b）基于修正的G-R关系模型，对祁连山—六盘山地震带1970年以来M_L3.0以上地震的震级-频度关系进行拟合，得到最大震级为9.66。钱小仕等（2013b）基于广义帕累托分布的超阈值分布模型，对中国大陆活动地块边界带强震震级分布特征开展研究，以1921年为起始年，震级阈值为5.3，估计了海源—祁连带震级上限为9.95。

对于广义帕累托模型，由于当ξ逐渐趋近于0时，模型估计的震级上限值可能会变得不稳定，Pisarenko等（2008，2014）提出可以利用分位数来进一步量化尾部特征。钱小仕等（2013a，b）也提到在工程地震危险性分析中，潜在最大震级的确定往往与特定建筑结构的抗震设计年限相关，他们利用0.999 97高分位数作为最大震级估计，其值往往低于广义帕累托模型公式直接计算的震级上限估计值。本文旨在提出基于全域敏感性分析的地震危险性估计不确定性分析流程和方法，震级上限估计值未采用分位数的形式。

基于现有的地震活动性模型，进行地震危险性估计，都不可避免地存在不确定性。

2.3   不确定性分析

构建地震活动性广义帕累托模型，需要人为选取历史地震记录起始时间和震级阈值。这些输入参量选取的不确定性可能导致地震危险性估计结果的不确定性。本文采用E-FAST方法开展定量的全域敏感性分析。

地震活动性广义帕累托模型的输入参数为地震目录起始时间t_s和震级阈值u，输出结果为不同时期的震级重现水平、震级上限及其对应的置信区间。青藏高原东北缘地震目录起始时间范围的选取是基于对该地区地震目录完整性的研究（黄玮琼等，1994）；震级阈值的选择范围通过上述1.1节中阈值的两种选择方法确定。设地震目录起始时间范围为1880—1890年，震级阈值的取值范围介于5.5—5.9之间，并假设它们在其不确定性范围内服从均匀分布，利用E-FAST方法进行蒙特卡洛抽样，共采集输入参数的抽样样本194组。

对青藏高原东北缘地震活动性广义帕累托模型进行参数敏感性分析。首先，根据输入参数的样本组合，计算相应的震级上限，以及对应一定重现期的震级重现水平；其次，将输入参数和相应的地震危险性分析结果文件导入敏感性分析软件中，利用E-FAST方法，计算广义帕累托模型地震危险性估计结果对各输入参数的敏感性指标，即t_s和u的主效应指标S_i和全效应指标$S^{\rm{T}}_i $ （表3）。为了更直观地显示各输入参数的敏感性指标，绘制如图6—8所示的直方图。

表  3  地震危险性估计结果对地震目录起始时间t_s和震级阈值u的主效应指标S_i和全效应指标$S^{\rm{T}}_i$

Table  3.  Total and the first-order effects of the seismic hazard estimation on earthquake catalogue initial time t_s 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 

| 显示表格

图  6  地震危险性估计结果对地震目录起始时间t_s和震级阈值u的主效应指标S_i直方图

横轴数字代表的含义下同

Figure  6.  Histogram of the first-order effect S_i of seismic hazard estimation on earthquake catalogue start time t_s and magnitude threshold u

The meanings of the number represent are the same below

下载: 全尺寸图片 幻灯片

图  8  地震危险性估计结果对地震目录起始时间t_s和震级阈值u的主效应指标和交互效应指标直方图

Figure  8.  Histogram of the first-order effect and interactions on earthquake catalogue start time t_s and magnitude threshold u

下载: 全尺寸图片 幻灯片

图6为广义帕累托模型的地震危险性估计结果对两个输入参数的主效应指标直方图。从图上可以看出，广义帕累托模型估计结果明显对u更为敏感。对于t_s，除了100年震级重现水平对应的主效应指标略高于0.2，其余的主效应指标均不到0.2。u的主效应指标远高于t_s的主效应指标，因此在建立模型时，首先要考虑震级阈值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  地震危险性估计结果对地震目录起始时间t_s和震级阈值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 t_s and magnitude threshold u

下载: 全尺寸图片 幻灯片

从图7和8可以看出，两个参数对地震危险性估计结果不确定性的影响均存在非线性效应，即两个参数之间存在相互作用。对应不同的重现期，在影响地震危险性估计结果不确定性上，两个参数之间的相互作用程度不同。对比不同重现期的震级重现水平，输入参数u的全效应指标中，主效应占主要部分，交互效应的比重随重现期的增大而逐渐增大；而另一输入参数t_s的全效应指标中，交互效应占主要部分。对于震级上限而言，两个输入参数的全效应指标中，交互效应占主要部分，表明这两个参数主要是以交互作用的形式对震级上限的不确定性产生影响。

3.   讨论与结论

本文选取青藏高原东北缘为案例研究区，提出了基于全域敏感性分析的地震危险性估计不确定性分析流程和方法。

基于广义帕累托分布，构建了研究区地震活动性模型，进行了地震危险性估计（包括震级重现水平和震级上限的估计）。鉴于构建模型需要人为选取t_s和u这两个输入参数，其取值的不确定性，可能导致地震危险性估计结果的不确定性，因此，选取具有全域敏感性分析功能的E-FAST方法，定量分析了上述两个输入参数的不确定性对地震危险性估计结果不确定性的影响。结果表明：地震危险性估计结果，对两个输入参数中的u更为敏感；对应不同重现期，地震危险性估计结果对u的敏感性程度有所不同；两个输入参数对地震危险性估计结果不确定性的影响，存在非线性效应；对不同的重现期而言，在影响地震危险性估计结果不确定性上，两个输入参数之间的相互作用程度不同。

目前，确定某一构造或构造区的最大震级，常用的方法有概率统计方法和确定性方法。前者有基于G-R关系统计分析的方法或利用震级-断层破裂尺度之间经验关系统计分析的方法；后者有基于研究区历史最大地震估计的方法和构造类比法等。新近的研究，开始倾向于以形变观测等资料作为约束，结合构造或区域的历史地震信息，估计该区域的最大震级。

本文选用基于广义帕累托分布构建的地震活动性模型，进行地震危险性估计，所用的方法属于概率统计方法。这一方法，利用历史地震记录中大震级区段的信息，无需先验地选定震级下限，较为适合历史地震记录时间长但低震级地震记录缺失的地区，便于构建研究区强震活动性模型，和其它概率统计方法一样，也有其局限性，即未能考虑地震发生的物理机制和物理过程，因为仅利用历史地震记录中大震级区段的信息，有时统计所用的历史地震记录样本偏少。

必须指出，不论通过哪一种方法构建地震活动性模型进行地震危险性估计，都难以避免不确定性，因此有必要开展相应的不确定性和敏感性分析研究，分析地震活动性模型不确定性的来源和特征，即分析输入参数的不确定性如何影响和制约地震危险性估计的不确定性。这些分析结果，可以指导地震活动性模型构建时输入参数的选取和调整，从而使地震活动性模型构建得更为完善。

本文提出的定量全域敏感性分析方法，对地震活动性模型的性质没有特殊的限定，便于在地震危险性估计不确定性分析中推广应用。