引言

G-R关系式描述了区域地震活动频次与地震震级大小的关系，是地震学中最重要的统计关系式之一（Gutenberg，Richter，1954）：

式中，N表示震级M以上地震的频度，a表示区域内的地震活动水平，b则反映区域内不同震级地震的比例关系. 不同研究人员通过岩石实验等结果讨论了b值的物理意义，认为b值是表征应力状态和介质特征的指标（Scholz，1968；Mori，Abercrombie，1997；Amitrano，2003；Schorlemmer et al，2005），与环境应力大小成反比，低b值区具有更高的应力积累。此外，b值大小还与该区域地下介质的脆性、弹性、塑性、碎裂程度等岩石特性有关，而且也常常与地震断层的分形维数和自组织临界性（self-organization criticality，简写为SOC）模型中的标度常数联系起来（Bak，Tang，1989；Ito，Matsuzaki，1990）.

鉴于b值的物理意义和实用意义，地震学家对其进行了深入研究并将其应用于地震预测研究和地震危险性判断。大量的试验和实际震例表明，在岩石破裂或大地震前，b值均有明显的变化（李全林等，1978；马鸿庆，1978）. 其机理可解释为：强震往往发生在活动断裂带上的高应力凹凸体内（Aki，1984）或闭锁段内（Wiemer，Wyss，1997；Wyss，Wiemer，2000；王熠熙等，2015），地震破裂一旦成核，将会在高应力环境下扩展，从而使凹凸体内的b值降低（Wyss，2004），这一现象均已在实验室条件下的岩石破裂实验（Scholz，1968）、与流体抽取有关的地震活动（Wyss，1973）以及地下矿山岩石的破裂（Urbancic et al，1992）中观测到. 在所有的影响参数中，b值对地壳中的有效剪应力最为敏感，并与应力成反比关系（Scholz，1968；Wyss，1973；Urbancic et al，1992），因此，低b值可以作为高有效剪应力的表征指标（Wyss，Wiemer，2000）. Wiemer和Wyss （1997）提出通过小尺度（5—15 km）b值空间扫描，由异常低b值（b＜0.7）勾画高应力凹凸体、由异常高b值（b＞1.1）（Wyss，2004）勾画蠕滑断裂段的方法，并将其用于活动断裂带不同段落的强震危险性分析（Wyss，Stephen，1998；Wyss，Wiemer，2000；易桂喜等，2010）. Wiemer和Schorlemmer （2007）提出基于凹凸体的概率模型（asperity-based likelihood model，简写为ALM），通过赤池信息准则（Akaike information criteria，简写为AIC）系统地验证了b值的空间变化对未来地震活动性的预测意义. Kamer和Hiemer （2015）使用基于惩罚性概率和泰森多边形图的方法讨论了加州地区的b值，其结果显示，在加州检测能力较强的地区，b值介于（0.94±0.04）与（1.15±0.06）之间，使用传统的固定半径或者邻近法则所获得的b值空间变化可能会被高估，且会受到主观参数选择的影响，而基于概率准则可以构建和改进预测模型. 总之，b值的物理意义及其对地震活动的预测意义已被证实，计算方法的选择和计算参数的选取对b值计算结果的影响也有不同程度的探讨。为了探讨这些结论在不同研究区域、不同观测条件下的适用性，本文拟将b值空间扫描方法应用于鄂尔多斯地块区域。

鄂尔多斯地块北缘的新构造运动强烈，新生代沉积厚度可观，属典型的断陷盆地，具有发生强烈地震的构造条件，历史上曾发生过公元前7年和公元849年两次河套大地震，近几十年，中强地震频频发生，成为华北地区第四个地震活动期破坏性地震的主要活动区域之一（孙加林，1985）。作为独立的稳定构造单元，鄂尔多斯地块具有相对稳定、统一的应力场作用模式.

中国地震台网地震目录包含了1965年以来相对完整的地震事件资料，为b值研究提供了可靠的数据基础，特别是1970年至今多次中强以上地震，为b值研究提供了震例基础. 本文拟通过时空扫描方法计算鄂尔多斯地块不同时间和空间的b值分布，并通过震例回溯、R值评分法和Molchan图表法来检验b值的预测效能，以期对鄂尔多斯地块北缘b值的时空变化特征进行分析，探讨b值对于中强以上地震的预测效能.

1.   构造背景和数据

本文研究范围为鄂尔多斯地块北缘（39.5°N—41.5°N，105°E—115°E）. 鄂尔多斯地块处于华北克拉通西部，于始新世开始形成，至今仍在活动，其内部变形较小，周缘边界断裂带变形强烈（徐锡伟等，1994）. 鄂尔多斯地块北缘的主要构造有阴山隆起区、鄂尔多斯隆起区以及位于其间的河套断陷盆地. 自华北克拉通形成以来，该区域经过多期不同性质的构造活动，内部形成构造特征不同的次级构造单元，新构造运动强烈，中强地震活动较为频繁（曹刚，2001）. 其中河套断陷带为EW向剪切拉张带，由呼包、白彦花和临河等3个次级盆地右阶斜列组成，这几个盆地均是北深南浅、不对称，分别受北侧大青山、乌拉山和色尔腾山山前正断裂所控。河套断陷带南缘受鄂尔多斯北缘断裂控制，是河套不对称断陷盆地带在形成过程中产生的次生断裂，并在地表形成断层崖，显示了左旋正走滑断裂的活动特征（邓起东等，1999）. 历史上该区域曾发生公元前7年河套大地震、公元849年包头西M7左右地震，1970年以来发生了1976 年和林格尔M6.3、1979年五原M6.0以及1996年包头西M6.4等中强地震。鄂尔多斯地块北缘地区则是近年来鄂尔多斯周缘中强地震最为活跃的地段，1970年至今，共发生M5以上地震17次，其中独立地震事件14次. 本文选用1970年至2017年9月的中国地震台网统一地震目录，原始目录中包含地震事件19 664条，图1给出了研究区域的构造背景和地震分布.

图  1  鄂尔多斯地块北缘的构造背景和地震分布

Figure  1.  Tectonic settings and earthquake distribution in the northern margin of Ordos block

下载: 全尺寸图片 幻灯片

1970年以来，伴随鄂尔多斯地块北缘M6地震的活跃，该区域的小震活动也经历着复杂的起伏变化，地震活动的群集特征会对b值计算的背景性、区域性和真实性产生影响，因此，需要从这些群集地震事件中将“独立”地震事件剥离出来进行研究（韩晓明等，2016）。进行b值计算之前，运用伪随机时空窗方法（Reasenberg，1985）对研究区地震目录中的余震予以删除，图2a给出了余震删除前、后的地震频次时序统计，图2b给出了相应的震级频次统计，删除余震后地震事件共计14 498条.

图  2  余震删除前后的地震频次时序统计（a）和震级频次统计（b）

Figure  2.  Seismic frequency time series statistics （a） and magnitude frequency statistics （b） before and after aftershock deletion

下载: 全尺寸图片 幻灯片

受到台网监测能力的影响，小震级段地震的缺失会造成拟合得到的震级-频次曲线的斜率较真实的低，因此在b值计算中利用G-R关系时，确定最小完整性震级很重要. 本文使用完整性震级范围方法（entire-magnitude-range，简写为EMR）（Woessner，Wiemer，2005）对鄂尔多斯地块北缘的地震目录进行完整性评估，计算最小完整性震级M_c，再使用 Bootstrap方法的蒙特卡罗近似来估计M_c的不确定度δM_c. 图3给出了鄂尔多斯地块北缘地震完整性震级的时空分布，可见M_c在空间和时间上均表现出不均匀性：空间上鄂尔多斯地块东北缘的监测能力较高（图3a），M_c可达到1.5左右，而地块西北缘的监测能力较弱；从时间上看，鄂尔多斯地块北缘的监测能力随着时间有较大提升（图3c），1970—1975年M_c处于M_L2.3—2.5之间，1995年之后稳步提高，2016年后达到M_L1.0左右.

图  3  鄂尔多斯地块北缘的最小完整性震级计算结果

（a） M_c空间分布；（b） δM_c空间分布；（c） M_c时序分布

Figure  3.  Minimum completeness magnitude of the northern margin of Ordos block

（a） M_c spatial distribution；（b） δM_c spatial distribution；（c） M_c time-series distribution

下载: 全尺寸图片 幻灯片

考虑到数据资料的有限性，如果按照最小完整性震级M_c值来截取地震样本，将大大减少每组参与b 值计算的样本数，并最终影响计算结果的精度和可信度. 根据目录完整性的评估结果，并借鉴韩晓明等（2016）关于该区域截止震级的选取方式，本文选取截止震级为M_L2.0，将截止震级以上的6 975条地震作为b值计算的基础数据. 图4给出了目录中全部地震与截止震级以上地震的地震频次的时序对比图（时间窗长为半年），可见：目录中全部地震事件的半年频次在时间上呈不均匀分布，1970年至1995年间地震频次变化较平稳，1995年后地震频次大幅增加，2008年后也有明显的升高，这两个时段对应于中国地震台网“九五”和“十五”的改造，因此这种地震频次的突变不是真实地震活动变化所引起，而是台网监测能力提高，使得记录地震次数增加所致，最小完整性震级M_c的时序图（图3c）也反映了这种改变. 由图4还可以看出，截止震级以上地震的频次变化平稳，符合实际地震发生率的变化，这说明对于M_L2.0以上地震的记录能力比较平稳，选取M_L2.0作为截止震级下限是合理的.

图  4  参与计算的地震目录（M_L≥2.0）与全部地震目录的频次对比

Figure  4.  Comparison of frequencies of earthquakes used in calculation （M_L≥2.0） with those of all earthquakes

下载: 全尺寸图片 幻灯片

2.   b值时空扫描

2.1   b值计算方法及参数选取

破坏性大地震是在相当大的时空范围内，经历长时间的弹性能量积累后而发生，并不是突然发生，因此对b值变化进行动态时空扫描，可以作为观察应力集中、转移以及监测破坏性地震孕育过程的一种手段（李全林等，1978）. 常用的b值估计方法包括最大似然法和最小二乘法，其中最大似然法计算简便，且不受个别较大地震的影响，对所有地震的震级用同样的权重求平均，相当于给众多小地震信息加权，可在一定程度上排除因个别较大地震的发生所引起的b值突变，因此，本文选用最大似然法计算b值。

若N_i服从泊松分布，考虑到对数概率分布和最大真值原理，宇津德治提出最大似然估计b值（Utsu，1978）的方法：

式中，M₀为起算震级，$\bar M$为M₀以上地震的平均震级，lge＝0.434 3，当n为地震总次数时，95%置信水平的b值标准偏差$\sigma$为

另一方面，b值搜索半径r的设置应该与地震数据的集中分布程度相当，这样才可能使扫描得到的b值更加可靠，但是中小地震分布的集中程度具有一定的不均匀性，一般情况下r取值较小时，结果较可靠，而当r值较大时，所求b值涉及的范围较大，但对b值参数间细微变化的反映不敏感（张琳琳，唐兰兰，2015）. 参考已有研究（易桂喜等，2011；李正芳，周本刚，2014），本文在时空扫描中，按0.2°×0.2°划分格点，每个格点的b值使用搜索半径r≤0.5°范围内且预测相应时间段之前10年内的地震事件计算，最小地震次数设为30.

2.2   结果和分析

经过b值时空扫描，得到不同时段的b值空间分布图。图5a给出了b值的时空扫描结果，b值范围为0.3—2.4；图5b为95%置信水平的标准差δb的时空分布，可见绝大多数区域95%置信水平的标准差达到0.2以内. 根据我国经验，将b值低于0.7作为异常（梅世蓉，1993）。从空间分布上看，1970年至今低b值区域主要分布在鄂尔多斯地块西北缘的阿拉善地区和东北缘的晋冀蒙地区，表明两地区为应力积累区域. 从时间上看：1970年至今低b值区域的空间分布有局部变化；1996年包头M_S6.4地震后，包头地区的低b值区域逐渐消失，反映了包头地震后该区域的应力水平降低；鄂尔多斯地块西北缘从1975年至今均表现为低b值，而地块东北缘的b值近年来有升高的趋势，现阶段低b值特征不明显.

图  5  鄂尔多斯地块北缘b值 （a）和δb （b）的时空扫描结果

Figure  5.  Spatio-temporal scanning results of b value （a） and δb （b） in the northern margin of Ordos block

下载: 全尺寸图片 幻灯片

3.   b值时空扫描的预测效能检验

3.1   震例回溯

1970—2017年期间鄂尔多斯地块北缘共发生中强以上地震17次，其中3次为前震或余震，14次为独立地震事件. 对14次独立地震事件，使用主震前10年间的地震目录进行b值时空扫描，图6给出了14次地震前的b值空间分布，红色区域为低b值区域. 以0.7为阈值，将低于该阈值的区域及其邻近区域作为预测区域. 表1给出了b值的预测效果，可见10次地震发生在低b值区域的内部及边缘，4次发生在非低b值区域，表明主震前的低b值区域对中强地震的位置有较好的指示意义，即中强地震易发生在低b值区域的内部或边缘.

图  6  鄂尔多斯地块北缘14次震例的b值回溯检验

Figure  6.  Retrospective test for b values of the fourteen earthquake examples in the northern margin of Ordos block

下载: 全尺寸图片 幻灯片

表  1  鄂尔多斯地块北缘14次震例的回溯情况

Table  1.  Retrospective test for the fourteen earthquake examples in the northern margin of Ordos block

编号	年-月-日	北纬/°	东经/°	MS	地点	低b值
1	1976-04-06	40.20	112.10	6.2	和林格尔	×
2	1976-09-23	39.90	106.40	6.2	巴音木仁	√
3	1979-08-25	41.20	108.10	6.0	五原	√
4	1981-08-13	40.58	113.41	5.6	丰镇	√
5	1989-10-19	39.92	113.91	5.9	山西大同	√
6	1991-01-13	40.40	105.60	5.2	阿拉善左旗	×
7	1991-03-26	40.00	113.80	5.8	山西大同	√
8	1991-09-14	40.20	105.50	5.1	阿拉善左旗	×
9	1996-05-03	40.80	109.60	6.4	包头	√
10	1997-10-21	41.20	107.10	5.0	乌拉特后旗	√
11	1998-01-10	41.10	114.30	6.2	河北张北	×
12	1999-03-11	41.20	114.60	5.6	河北张北	√
13	1999-11-01	39.80	113.90	5.6	山西浑源	√
14	2015-04-15	39.80	106.30	5.8	阿拉善左旗	√

下载: 导出CSV 

| 显示表格

3.2   Molchan图表法检验

Molchan图表法（Molchan，1997）由于能够客观、科学地进行地震预测评估，目前已被广泛地应用于确定性和概率性预测的统计检验和效能评估中（蒋长胜等，2011，孙丽娜等，2012，Chen et al，2016；薛艳等，2017），而且在国际地震可预测性研究计划中已列为常规使用. 通过Molchan图表法可以对地震前兆的预报效能进行评估，并得到概率增益曲线，其间定义时空占有率τ和漏报率ν，当异常阈值由最高变化至最低时，时空占有率τ从0变到1，漏报率ν从1变到0，最佳的预测效果对应于最大的预测成功率和最小的时空占有率. Molchan图表法还可通过显著性水平α和概率增益值G有效地评价地震危险性概率预测模型的预测效能，也可通过漏报率$\upsilon $（τ）曲线向左包围的面积来判断不同预测模型整体预测效果的优劣（Zechar，Jordan，2010）.

使用某个指标预测时空范围的目标地震，设目标地震的次数为N，其中有≥h次事件“击中”的分布概率可写为二项式形式，则相应的显著性水平α值为

式中，N为地震次数，α为显著性水平，τ为时空占有率，h为击中次数。

以空间步长0.2°，时间步长2年将鄂尔多斯地块北缘区域（空间范围：39.5°N—41.5°N，105°E—115°E；时间范围：1970—2017）的时空预测范围划分为1万个时空格点，以空间上低于阈值的区域为预测区域，计算时间窗长为10年，预测标的为2年内该格点区域发生的M≥5.0地震. 由于选取的计算时间步长和预测时间均为2年，可以通过格点数目计算不同阈值下的漏报率、报准率和时空占有率，其中时空占有率τ为超过阈值的格点数与全部格点数之比，对应的漏报率$\upsilon $为未发生在低于阈值格点内的地震数与全部地震数之比. 图7给出了b值时空扫描的Molchan曲线，可见当阈值从0.3变化至2.4时，时空占有率从0增至0.8，漏报率从1降至0.1. 从图7还可看出，b值时空扫描的预测效能优于泊松分布，表明b值具有一定的预测效能，当阈值为0.8时，概率增益达到2，显著性水平达到25%.

图  7  Molchan检验

（a） 不同概率增益G下的τ - υ曲线；（b） 不同显著性水平α下的τ - υ曲线；

Figure  7.  Molchan test

（a） τ - υ curve with different probabilistic gain G；（b） τ - υ curve with different significant level α

下载: 全尺寸图片 幻灯片

3.3   R值评分法检验

R值评分法由许绍燮（1989）提出和改进，主要考虑实际预测地震的有震报准率c和预报时间占有率d. R值定义为

式中后一项可推广为时空占有率，则

R值的统计含义是扣除了随机概率的预报成功率. 统计样本的数量会影响R值的置信水平，许绍燮（1989）编制了保证97.5%置信水平的最低R值（R₀）表，如果实际计算R值大于表中对应值，则认为该R值至少有97.5%置信水平. R值的大小可以反映指标的预测价值，对与异常指标有关的参数，例如判定异常的起始值、异常持续时间、异常范围、异常幅度等，以及诸如震级阈值、地震地点（与异常点的距离）、发震时间（与异常的间隔时间）等与所研究地震相关的参数的不断调整，多次计算R值，这样通过参数调整，可以寻找预报效能最佳的R值（郑兆苾等，2001）. 本文中，我们计算了不同阈值下的报准率和时空占有率，绘制了R值随阈值变化的曲线，如图8所示，同时给出了R₀值曲线作为参考. 由图8可看出，当阈值从0.3增至2.3时，R值在0.14—0.28范围内变化，当阈值为0.8时，有最高R值0.28，并且达到R₀.

图  8  R值随阈值的变化曲线

Figure  8.  R value versus threshold curve

下载: 全尺寸图片 幻灯片

4.   讨论与结论

选取鄂尔多斯地块北缘1970—2017年的地震目录，使用伪随机时空窗方法删除余震后，使用完整性震级范围方法对目录的完整性进行评估，其结果显示，最小完整性震级M_c在空间上和时间上均表现出不均匀性：空间上，鄂尔多斯地块东北缘的监测能力较高，M_c达到1.5左右，而该地块西北缘的监测能力较弱；时间上，M_c逐步降低. 权衡数据资料充足性和地震目录的完整性，选取M_L2.0以上的6 975条地震参与计算.

通过地震目录完整性评估确定起算震级和删除余震后，使用最大似然法对鄂尔多斯地块北缘的b值进行时空扫描，在扫描的范围内，绝大多数区域95%置信水平的标准差在0.2以内. b值在空间上分布不均匀，但在时间上有一定的稳定性，反映了中长期的应力环境及其变化，低b值区域主要分布在鄂尔多斯块体西北缘和东北缘. 鄂尔多斯块体西北缘从1975年至今均表现为低b值，低b值范围有一定变化；鄂尔多斯块体东北缘的b值有升高趋势，现阶段低b值特征不明显.

对14次独立的中强以上地震前的b值分布进行回溯性检验，其中10次地震发生在低b值区域的内部或边缘，4次发生在非低b值区域，表明主震前的低b值区域对中强地震的位置有较好的指示意义，即中强地震易发生在低b值区域的内部或边缘.

使用Molchan图表法和R值评分法评估不同阈值下的b值预测效能，其检验结果显示，b值时空扫描的预测效能优于泊松分布，表明b值具有一定的预测效能，当阈值为0.8时，概率增益达到2，显著性水平达到25%. 为了对比不同阈值的预测效果，绘制R值随阈值变化的曲线，结果显示当阈值从0.3增至2.3时，R值在0.14—0.28之间变化，当阈值为0.8时，有最高的R值0.28. 根据本文的分析，可将b值低于0.8作为鄂尔多斯地块北缘中强以上地震的预测指标.

审稿专家为本文的撰写提出了科学细致的修改建议，作者在此表示感谢。