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.
-
引言
地磁场是地球最重要、最基本的物理场之一，包含丰富的地球物理信息（Němec et al，2008）。当前地磁场观测台站建设成本高、占地面积大，且由于城市扩张，高铁、地铁发展迅速以及基础建设需求较多，许多现有台站面临着严重且难以辨析的观测噪声干扰。台站仪器可能遭受雷击无法获得数据（姚休义等，2016）。与地面观测相比，井下温度波动较小、能滤除大部分地表干扰、观测精度可提高1—2个数量级（徐纪人，赵志新，2009），选址方便，背景噪声辨析度高（胡星星等，2011），此外还具有观测投入成本比建设地面观测台站要小得多以及几乎不占用地表土地面积等优点。因而，井下地磁观测系统逐渐成为国际上的先进观测方式（任烨等，2012）。
光泵磁力仪常被用作井下地磁观测的绝对测量。本文涉及的是一种分辨率为0.2 nT的跟踪式氦光泵磁力仪（胡睿帆，2017）。在光泵浦和磁共振的共同作用下（McGregor，1987; Chéron et al，1996; Gilles et al，2001），射频场的频率与待测磁场间构建了线性关系，从而可以通过测定射频场频率得到待测磁场的强度（张振宇等，2011a，b），在信号电路中通过捕捉基波和二次谐波的峰值进行频率的快速定位（张振宇等，2011c，d）。本文所用到的氦光泵磁力仪的测频公式（胡睿帆，2017）为
$$ f {\text{＝}} {B_0}\dfrac{{{\gamma _{\rm{s}}}}}{{2{\text{π}} }} {\text{＝}} 28.02{B_0}{\text{，}} $$ (1) 式中：γs为氦原子的磁旋比率，是物理常数；B0为待测磁场。考虑到地磁场范围介于3×104—7×104 nT，根据式（1），则射频场频率介于840.70 kHz—1.96 MHz，1 Hz的频率测量误差会导致0.036 nT的磁场测量误差，若磁力仪分辨率低于该测频误差一个量级，则认为频率测量精度较高。
频率测量有繁多的成熟技术。从硬件上划分，包含但不限于模拟电路测频、数字电路测频、数字信号处理结合现场可编程逻辑门阵列（digital signal process plus field programmable gate array，缩写为DSP+FPGA）测频、单片机（microcontroller，缩写为MCU）测频。井下光泵磁力仪的测频模块有以下三个必备的属性需求：一是体积小。井眼半径有限，信号处理电路不应过大，以光泵磁力仪探头大小为基准较为合适；二是对温度不敏感。井下温度范围宽，要求测频模块应在较大温度变化范围内始终能高精度测量；三是功耗小。井下环境窄小，井眼位置偏僻，从长期地磁观测或短期勘探勘察等方面考虑，为适应无法接入市电的地区，应减少整机功耗，提高电池续航能力，测频模块作为磁力仪的一部分，更应该尽可能减少功耗。传统的模拟或数字频率测量电路体积较大，且随着测量频率的提高，需要更多的外围电路以满足高精度测频要求，其中采用锁相环技术与直接数字频率合成（direct digital synthesis，缩写为DDS）技术的专业频率计均可以实现宽频带、高分辨率、高精度的频率测量，但体积和功耗较大；DSP＋FPGA测频相比于单片机测频在同体积条件下，对光泵磁力仪频率信号范围内（840.70 kHz—1.96 MHz）的测量精度较低，且FPGA电路功能较为单一，而单片机频率计不仅测量精度高且能够以更小体积的电路实现更强大更灵活的控制性能。鉴于此，基于单片机的频率测量是满足条件的较优解决方案。
本文使用MCU定时器外部时钟模式计数测频，体积小、功耗低。在840.70 kHz—1.96 MHz范围内，频率精度较高，误差稳定在−1 Hz，折算至地磁场，误差为−0.036 nT。此误差来源于MCU调用中断服务函数（interrupt service rountines，缩写为ISR）的时长大于待测信号1个周期，在上述频段范围内是稳定的−1 Hz误差，易通过软件进行补偿。
1. 单片机频率测量的实现方法
许军（1998）对直接测频法、测周法、变闸门测频法、细分测频法的频率误差及适用情况进行了分析，指出直接测频法适用于对频率较高信号的测量。马献果和焦阳（2004）给出了测周法和测频法的频率分界点的计算方法，即
$$ {f_{\rm{C}}} {\text{＝}} \sqrt { {\frac{{{f_{\rm{S}}}}}{{{T_{\rm{G}}}}}} } {\text{，}} $$ (2) 式中，fC为频率分界点，fS为标频信号频率，TG为闸门时间。对于给定的标频信号fS和闸门时间TG，总存在一个频率分界点fC，当待测频率小于频率分界点fC时，适合使用测周法，当待测频率大于频率分界点fC时，适合使用直接测频法。
光泵磁力仪的基波标频信号是1 MHz左右，我们取标频信号fS＝1 MHz，闸门时间TG＝1 s，从而频率分界点fC＝1 kHz，而待测频率范围为840.70 kHz—1.96 MHz，所以选择直接测频法。
如引言所说，考虑到项目中使用的光泵磁力仪需要小型化设计，且需要能够在井下无人操作且温度未知的环境中进行持续可靠的高精度测量，频率测量及数据处理设备不应较大，故排除大型频率计及传统的模拟电路和数字电路。同时考虑到DSP＋FPGA对高频信号测频误差较大，故选择使用单片机直接测频。
单片机利用定时器实现直接测频有多种方式：输入输出（input-output，缩写为IO）接口外部中断服务函数直接计数测频；定时器输入捕获模式测频；定时器外部时钟模式计数测频。
1.1 IO口外部中断服务函数计数测频
单片机定时器1进行1 s定时，待测信号在芯片IO口处的每一次跳变产生一次中断，定时器2对每一次产生的中断计数。单位时间内的中断次数就是待测信号的频率。使用这种方法测量的低频信号精度很高。但IO口外部中断服务函数计数测频占用CPU资源很大，并且中断函数调用与运算有一定耗时，对于频率较高的信号会产生漏计，从而产生较大误差。本文待测信号的频率范围是840.70 kHz—1.96 MHz，远高于这种方法的适用范围，不适合使用IO外部中断服务函数计数测频。
1.2 定时器输入捕获模式频率测量
图1给出了定时器的输入捕获模式。从图中可以看出，定时器用作计数器，待测信号ftest频率较低，定时器计数器频率fTIM较高，为已知的设置值。在待测信号的两次上升沿处记录定时器的计数器（time’s counter，缩写为CNT）寄存器的值，两次计数器的差值为ΔCNT，从而，待测频率为ftest＝fTIM/ΔCNT。输入捕获模式对中低频信号的频率测量比较准确，此时定时器时钟频率比待测信号大很多，ΔCNT的±1误差影响很小。当待测信号频率较高时，ΔCNT的±1误差的影响很大，导致最终测量误差增大，对于100 kHz及以上频率的待测信号，最终测量误差较大。
1.3 定时器外部时钟模式计数测频
基于Cortex-M3内核的STM32F103ZET6芯片的通用定时器可以使用外部信号作为时钟信号。因此，把待测信号作为定时器的时钟信号，使定时器工作在计数器模式，就可以进行频率测量。与IO外部中断服务函数计数相比，省去了大量的函数调用与计算时间，减少芯片的资源占用；与输入捕获模式相比，高频段±1误差影响较小。因而在840.70 kHz—1.96 MHz频率范围内，这种方法精度很高。
用一组定时器作为定时模块，定时1 s时间闸门。用另一组定时器作为计数模块，使用外部时钟模式对待测信号进行计数。单位时间的计数值就是频率。如图2所示，计数模块不断对信号源进行计数，同时计数模块在一定程度上受到定时模块的控制与调用。定时模块在每个时间闸门结束时给计数模块一个信号，这个信号首先读取并存储此时计数模块的计数值，然后使计数模块的定时器清零，从而开始下一次测量。读取、存储、清零这一系列工作在定时模块的中断服务函数中进行，并在这之后进行频率计算和打印输出。
本文的频率计使用高级精简指令计算机进行设计，满足了小型化的需要并为无人自动测量与数据传输作准备，且从后续测试结果可以看出，在840.70 kHz—1.96 MHz频率附近，频率计的精度较高。
单片机主频的稳定性对于处理器计算速度与精度非常重要，而影响主频晶振稳定性的主要原因是温度（孙敏等，2013）。本文设计的单片机频率计使用有源温补晶振（temperature compensated oscillator，缩写为TCXO）为处理器芯片提供更加稳定的主频。实际电路中对32 MHz的有源温补晶振进行四分频后供给芯片，从而满足主频8 MHz的需求。该晶振在常温25 ℃下的频差最大值是±1×10−7，在−40—95 ℃范围内的最大频差为±1×10−6，性能优秀，满足要求。
图3的光泵磁力仪流程图中给出了处于闭环反馈回路中的射频场发生器的工作流程，该发生器的工作频率可根据地磁场日变化不断调整，主要工作频段在840.70 kHz—1.96 MHz （对应地磁场范围为30 000—70 000 nT）。单片机需要测量的频率就是射频场发生器发出的频率。
2. 频率计的两种实现方法
2.1 二定时器法
二定时器法是指一个定时器完成计数功能，另一个定时器完成定时功能。
二定时器法框图（图4）显示，定时器TIM2处于定时器模式，用于产生1 s的时间闸门，定时器TIM3处于计数模式，用于对待测信号进行计数。两个定时器的主要配置参数列于表1。
表 1 二定时器法主要配置参数Table 1. Main configuration parameters of two timers method寄存器名称 TIM3 TIM2 TIMx_ARR 65536-1 60000-1 TIMx_PSC 0 1200-1 TIM3处于外部时钟模式，对待测信号不断计数，在TIM3中断服务函数中设置变量T3ova保存1 s内的溢出次数；TIM2进行1 200预分频，则TIM2计数60 000次时产生1 s闸门信号（闸门信号频率为60 kHz）；TIM3接收到1 s闸门信号，将此刻计数值T3c和溢出值T3ova发送给TIM2中断服务函数，TIM2中断服务函数计算待测频率并输出，频率计算公式为
$$ f {\text{＝}} 65\;536 {T_{{\rm{3ova}}}} {\text{＋}} {T_{{\rm{3c}}}}. $$ (3) 2.2 四定时器法
四定时器法是指对于计数模块与定时模块分别使用两个定时器完成对应功能。
四定时器法框图（图5）显示，高级定时器TIM8和通用定时器TIM2处于定时器模式，构成定时模块，用于产生1 s的时间闸门，通用定时器TIM3和TIM4处于计数器模式，构成计数模块。
定时模块设定为主从同步模式（external1模式），TIM8处于主模式，TIM8不断计数，其ARR寄存器重载时产生更新中断，使TIM2计数器自增1。计数模块也设定为主从同步模式（external1模式），即主定时器的更新中断作为从定时器的时钟信号。TIM3处于主模式，在外部时钟模式下对待测信号不断计数，ARR寄存器重载时产生更新中断，使TIM4计数器自增1。四个定时器的主要配置参数列于表2。
表 2 四定时器法主要配置参数Table 2. Main configuration parameters of four timers method寄存器名称 TIM3 TIM4 TIM8 TIM2 TIMx_ARR 65536-1 65536-1 50000-1 20-1 TIMx_PSC 0 0 72-1 0 可以看出，定时模块相当于给出了一个1 MHz的闸门时间信号（50 000×20 Hz），这个频率与待测信号相近，远高于二定时器法（60 kHz）。这减少了定时模块的中断更新次数，增加了1 s定时的精度。定时模块最终通过TIM2给计数模块一个闸门时间信号，接受到这个信号后，计数模块把计数值T4c和溢出值T3ova回传给TIM2。在TIM2中断服务函数中记录用于计算的计数值和溢出值，并重置TIM3和TIM4的计数器，准备再次测量计数。计算结果在TIM2的中断服务函数中打印输出。测量结果由式（4）确定，即
$$ f {\text{＝}} 65\;536 {T_{{\rm{4c}}}} {\text{＋}} {T_{{\rm{3ova}}}}. $$ (4) 3. 实验
3.1 系统误差实验
针对第二部分的两种实现方法，进行系统误差实验分析。使用同源测试来测试频率模块精度。待测信号由为芯片提供主频的TCXO分频提供。同源测试能够排除由于晶振自身不准确所带来的误差，从而确定来自系统本身或是程序本身的系统误差。
特别地，本文表格中的绝对误差是指测量频率与输入频率差值的绝对值，单位是Hz。相对误差为绝对误差除以输入值，是一个无量纲的小数。实验结果列于表3。
表 3 同源实验结果Table 3. Experiment result using the same signal source输入信号/kHz 二定时器法 四定时器法 测量值/Hz 绝对误差/Hz 相对误差 测量值/Hz 绝对误差/Hz 相对误差 62.5 62 500 0 0 62 500 0 0 250 250 000 0 0 25 0000 0 0 500 499 999 1.0 2.0×10−6 500 000 1.0 2.0×10−6 1000 999 997 3.0 3.0×10−6 999 999 1.0 1.0×10−6 2000 1 999 993 7.0 3.5×10−6 1 999 999 1.0 5.0×10−7 3.2 应用测试实验及不同频率计的实验结果比较
实际应用中，无法保证信号源与主频信号同步，多为异源测试。应用测试实验均使用四定时器法。例如使用异源有源温补晶振作为信号源，模拟实际测量。实验结果列于表4。或以高精度信号发生器作为信号源，为确保结果的准确性，并确定稳定的系统误差范围，使用频率精度为±1×10−6的UNI-T公司的UTG2062A高精度信号发生器作为异源信号源进行测试。实验结果列于表5。张前毅（2007）介绍了使用大规模可编程逻辑器件（complex programming logic device，缩写为CPLD）和高性能单片机制作的超高精度微波频率计，其部分实验结果列于表6；王永超（2012）介绍了一种基于CPLD的定闸门频率计，部分实验结果列于表7；谢尚豪等（2018）介绍了一种基于FPGA的中低频频率计，实验结果列于表8。
表 4 独立有源温补晶振作为信号源时的实验结果Table 4. Experiment results using independent active TCXO as signal source输入信号/kHz 实际测量值/Hz 绝对误差/Hz 相对误差 1 000 999 999 1 1.0×10−6 100 000 0 0 2 000 1 999 999 1 5.0×10−7 3 000 2 999 998 1 3.3×10−7 2 999 999 2 6.7×10−7 表 5 高精度信号发生器作为信号源时的实验结果Table 5. Experiment results using high-precision signal generator as signal source输入信号/kHz 四定时器法测量/Hz 绝对误差/Hz 相对误差 输入信号/kHz 四定时器法测量/Hz 绝对误差/Hz 相对误差 500 500 000 0 0 2 200 2 199 999 1.0 4.5×10−7 600 599 999 1.0 1.7×10−6 2 300 2 299 999 1.0 4.3×10−7 700 699 999 1.0 1.4×10−6 2 400 2 399 999 1.0 4.2×10−7 800 799 999 1.0 1.3.×10−6 2 500 2 499 999 1.0 4.0×10−7 900 899 999 1.0 1.1×10−6 2 600 2 599 999 1.0 3.8×10−7 1 000 999 999 1.0 1.0×10−6 2 700 2 699 999 1.0 3.7×10−7 1 500 1 499 999 1.0 6.7×10−7 2 800 2 799 999 1.0 3.6×10−7 2 000 1 999 999 1.0 5.0×10−7 2 900 2 899 999 1.0 3.4×10−7 2 100 2 099 999 1.0 4.7×10−7 3 000 2 999 999 1.0 3.3×10−7 表 6 微波频率计超高频段测量的实验结果Table 6. Experiment results of microwave frequency meter in ultra-high frequency measurement输入频率/Hz 实测频率/Hz 绝对误差/Hz 相对误差 5 000 000 000 5 000 000 373.12 373.12 7.5×10−8 5 100 000 000 5 100 000 326.17 326.17 6.4×10−8 5 270 000 000 5 269 999 878.63 121.37 2.3×10−8 5 520 000 000 5 520 000 061.08 61.08 1.1×10−8 5 603 000 000 5 603 000 047.03 47.03 8.4×10−9 表 7 CPLD频率计中频测量的实验结果Table 7. Experiment results of frequency memter of CPLD in intermediate frequency measurement输入频率/Hz 实测频率/Hz 绝对误差/Hz 相对误差 56 000.000 56 000.00 0 0 104 523.687 104 522.00 1.687 1.6×10−5 504 258.741 504 252.03 6.711 1.3×10−5 774 519.638 774 508.00 11.638 1.5×10−5 1 000 000.000 999 986.00 14.000 1.4×10−5 表 8 FPGA频率计中低频段测量的实验结果Table 8. Experimental results of FPGA frequency memter in low frequency range measurement输入频率/Hz 实测频率/Hz 绝对误差/Hz 相对误差 1 000 999 1.0 1.0×10−3 10 000 10 000 0 0 20 000 20 001 1.0 5.0×10−5 200 000 200 004 4.0 2.0×10−5 3.3 与其它频率计的比较
通过对比，可以看到，尽管微波频率计在超高频段（5 GHz）有很高的精度，但微波频率计不适用于光泵磁力仪测量磁共振频率。这是因为式（1）决定了光泵磁力仪信号频率范围是840.70 kHz—1.96 MHz。基于FPGA或CPLD的频率计在1.5 MHz—3 MHz的绝对误差都超过了14Hz，相对误差的数量级在10−5。另外，由于测频电路主体是相对传统的模拟电路，这类频率计的体积相对较大。在光泵磁力仪所需要的频率范围840.70 kHz—1.96 MHz内，本文设计的单片机频率模块，测量误差均在10−6之下，同时实现了较高的精度和较小的体积（长、宽均为5 cm）。从测量精度、电路尺寸、功耗等角度来看，单片机频率模块都有很大优势。
4. 误差分析及系统误差软件校正
4.1 中断服务函数带来的系统误差
系统误差实验结果表明，四定时器法能够有效地提高测频精度，但仍在500 kHz—3 MHz频段内有1 Hz的测量误差。考虑到同源实验中主频与待测信号同步，理论上不应存在误差。故认为这部分误差来自测量系统本身，属于仪器的系统误差。
图6给出了ISR产生系统误差的过程。由于定时模块结束时需调用ISR来完成读取计数模块计数值、计数模块清零、清除更新中断标志位、计算频率值等一系列指令时，中断服务函数本身会占用时间，从而导致高频待测信号在计数器清零过程中产生漏计。这也解释了在频率低于500 kHz时，同源实验没有测量误差，在500 kHz—3 MHz范围内，绝对误差是1 Hz的情况。可以预见，当待测信号频率更高时，误差会逐渐增大。当然，更高频率的待测信号不会出现在光泵磁力仪频率信号范围内。
4.2 直接测频法原理误差
直接测频法原理误差来自于闸门信号与待测信号不同步的±1量化误差ε1和闸门信号自身的标频误差ε2 （王永超，2012）。两者在一次实际测量中所带来的相对误差分别可通过式（5）和式（6）计算得到，即
$$ {\varepsilon _1} {\text{＝}} {\text{±}} \frac{1}{T_{\rm{G}}}, $$ (5) $$ {\varepsilon _2} {\text{＝}} {\text{±}} \frac{\Delta {f_{\rm{S}}}}{f_{\rm{S}}}, $$ (6) 式中：TG为闸门时间，本文中为1秒；fS为标频信号频率，本文采用72 MHz有源晶振分频后的1 MHz信号。
由式（1），地磁场范围对应的磁共振频率为840.70 kHz—1.96 MHz。结合式（5），|ε1|＝5.10×10−7—1.19×10−6。标频信号的误差主要来自有源温补晶振的精度，本文使用的32 MHz有源温补晶振在常温下（−10 ℃—60 ℃）精度可达±0.5×106，分频不影响晶振精度，故|ε2|＝5.0×10−7。综上所述，由直接测频法原理引入的相对误差的模值为|ε|＝5.10×10−7—1.19×10−6。
4.3 可行的系统误差软件校正
如4.1所述，在500 kHz—3 MHz频段内，系统本身有1 Hz的漏计的测量误差。考虑到地磁场在3×104—7×104 nT范围内波动，结合式（1），光泵磁力仪的磁共振频率为840.70 kHz—1.96 MHz。考虑到以上两点，实际应用中可对测量结果进行简单的软件校正以提高实际测量精度。
在中断服务函数中，若频率测量的范围在500 kHz—3 MHz内，则在中断服务函数中对测量值＋1后输出给函数返回值。软件矫正后以独立的有源温补晶振作为信号源的实验结果列于表9。
表 9 软件校正后以独立的有源温补晶振作为信号源Table 9. Independent active TCXO as signal source with software calibration输入信号/kHz 实际测量值/kHz 绝对误差/Hz 相对误差 1 000 1 000.000 0 0 2 000 2 000.000 0 0 3 000 2 999.999 0 0 3 000.000 1.0 3.3×10−7 可以看出，软件矫正后频率计在磁共振范围内，即840.70 kHz—1.96 MHz，实现了零误差。
5. 讨论与结论
井下氦光泵磁力仪的频率采集部分首先要满足高精度，其次是集成设备要足够小，功耗足够低，能够适用于井下工作环境。本文在比较了三种频率测量方法（传统模数电路、DSP＋FPGA、单片机）的优劣势后，选择单片机电路设计并实现了频率测量模块。比较了三种单片机的测频方法的优劣势后，选择使用四定时器法的外部时钟模式计数测频。
本文设计的频率计体积小，精度高，系统误差小且稳定。通过简单的软件校正可以实现高精度中高频段的频率测量。
中断服务函数本身及其附带的指令所占用的时间会造成高频信号测量结果的定量误差。如若从硬件考虑来消除这一误差，可以考虑使用直接存储器访问（direct memory access，缩写为DMA）获得计数模块定时器的计数值，每次计算频率用到的计数值直接通过DMA控制器实现读取与传输，不通过CPU调用资源，提高CPU利用率，减少系统延迟与误差。
-
![]()
图 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
计量
- 文章访问数: 251
- HTML全文浏览量: 119
- PDF下载量: 72
