引言

北京时间2021年5月22日02：04，青海果洛州玛多县发生M_S7.4地震（34.59°N，98.34°E），震源深度17 km。震中位于玛多黄河乡，距县城38 km，距西宁380 km，震中海拔高度4258 m。地震表现出无明显前震、余震较多的特征，主震后1小时内发生M_S3.0以上余震10次，最大余震为M_S4.3，一天内发生M_S3.0以上余震36次，最大余震为M_S5.1，为主震-余震型地震。中国地震台网中心对震源机制解进行反演得出本次地震为走滑型事件，地震破裂持续45 s，震后30 s内释放大量能量，最大滑移量4.5 m，地表破裂超过160 km（邓文泽等，2021）。本次地震为汶川8.0级地震以来，中国境内发生的震级最高、地表破裂最长的地震事件。本文将对此次地震震前b值进行研究，寻求该区域b值与应力的变化关系，为该区域后续地震预报提供科学依据和参考。

Gutenberg等（1944）发现震级M与地震频次N之间存在统计学关系logN =a−bM，a表示研究区域地震活动水平，b表示研究区域不同震级的相对分布。Kun等（2013）通过数值模拟研究验证了震级、频次的数学关系。研究人员在进行岩石实验时发现介质的应力水平不断增加时，b值会呈现逐渐减小的趋势（Mogi，1962；Scholz，1968；Rivière等，2018）。Mousavi等（2017）对整个板块进行研究，同样发现b值的变化与介质的应力变化存在一定关系。因此，研究学者使用ｂ值研究区域应力的大小（Nuannin等，2012；Scholz，2015）。易桂喜等（2006）对龙门山断裂带b值进行分析，得到断裂带区域的强震危险性；刘艳辉等（2015）对青藏高原东南缘b值进行时空扫描，发现在大地震发生前后b值会出现较大变化，地震发生前b值长时间处于低值；Tormann等（2015）研究得到全球平均b值约为1，中强震多发地区b值小于1，小地震多发地区b值大于1； Nanjo等（2018）利用b值与应力关系，得到海槽地震的初始破裂部位；谢卓娟等（2019）利用b值与应力关系研究首都圈地震活动性；曾宪伟等（2021）通过b值分析得到宁夏吴忠—灵武地区的强震危险性。前人已开展大量区域应力变化与b值之间的关系研究，在此基础上Nievas等（2020）提出可以利用b值变化研究区域应力变化，寻求b值变化与地震发生的关系。本文对玛多M_S7.4地震震前地震序列进行分析，研究地震前b值的变化特征，为研究区域地震预测和地震发生风险提供参考。

1.   区域地质背景

新生代以来，青藏高原不断受到欧亚板块和印度板块的挤压，使青藏高原呈现扩张趋势，表现出板块运动活跃、地震多发的特点，21世纪以来发生7级以上地震10次（Molnar等，2009）。巴颜喀拉块体是青藏高原的次一级构造单元，该块体被东昆仑、甘孜-玉树-鲜水河、龙门山3条大型断裂带包围，成为中国目前地震活动最为强烈的地区之一（张培震等，2003；Rui等，2016；闻学泽，2018）。其中，巴颜喀拉块体东缘最为活跃，十年内相继发生了2008年汶川8.0级地震、2013年芦山7.0级地震、2017年九寨沟7.0级地震（图1）。

图 1  青藏高原主要活动断裂及历史地震分布

Figure 1.  Simplified distribution of main active faults and historic earthquakes in the Qinghai-Tibet plateau

下载: 全尺寸图片 幻灯片

东昆仑断裂位于巴颜喀拉块体北边界，EW走向，长约2000 km，是一条巨型左旋走滑断裂带，倾角55°～85°。该断裂自西向东走滑速率不断降低，玛沁段为12.5 mm/a，至玛曲段降为5 mm/a，至塔藏段已不足3 mm/a（Kirby等，2007；Ren等，2013）。有记录以来，该断裂发生7级以上地震6次，最大震级为2001年M_S8.1地震，地表破裂达426 km（徐锡伟等，2002；邓起东等，2003；Xu等，2006）。

甘孜-玉树-鲜水河断裂带位于巴颜喀拉块体南边界，为左旋走滑断裂带。学者对该断裂带滑动速率有不同见解，闻学泽、石峰等认为甘孜-玉树段滑动速率为8～14 mm/a，但 Chevalier等认为甘孜-玉树段滑动速率仅有5～8 mm/a（闻学泽等，2003；石峰等，2013；Chevalier等，2018）。学者们普遍认为鲜水河段滑动速率为10～13 mm/a（Zhang，2013；Chen等，2016）。有记录以来，该断裂发生7级以上地震17次，最大震级为1997年M_S7.9地震。

龙门山断裂位于巴颜喀拉块体东边界，为逆冲右旋走滑断裂带，由多个近平行断裂组成。21世纪以来，该断裂活动加剧，2008年汶川M_S8.0地震、2013年芦山M_S7.0地震都发生在该断裂上。

巴颜喀拉块体除四周存在3条大型断裂带外，其内部也发育着一系列活动断裂。玛多-甘德断裂带位于东昆仑断裂带以南70 km，为全新世活动的左旋走滑断裂，整体呈NW-SE走向，长约650 km，水平滑动速率小于6～8 mm/a，活动速率明显小于位于其南北两边的甘孜-玉树-鲜水河断裂带和东昆仑断裂带（熊仁伟等，2010）。

2.   数据及方法

2.1   地震目录数据

本次研究地震目录由中国地震台网提供，目录选取时间为2009年1月1日至2021年4月30日。青海玛多M_S7.4地震研究范围选取震中所在的青藏高原东北部巴颜喀拉块体内部，南北范围选取东昆仑断裂带、甘孜-玉树-鲜水河断裂带之间区域（32°N～37°N），东西范围覆盖玛多-甘德断裂带且包含本次地震所有余震（96°E～100°E）。

由于巴颜喀拉块体内部发育着众多活动断裂，导致区域内5级以上地震较为分散。其中玛多M_S7.4地震及其余震沿玛多-甘德断裂带分布，呈NW-SE走向，空间分布长约150 km（詹艳等，2021；张喆等，2021；潘家伟等，2021）（图2）。

图 2  研究区域中强震（M_S≥5.0）及玛多地震余震（M_L≥3.0）分布图

Figure 2.  Simplified distribution of major earthquake（M_S≥5.0）and Maduo earthquake aftershocks（M_L≥3.0）in the study area

下载: 全尺寸图片 幻灯片

在b值计算时，M_S震级适用于大范围、长时间、高震级（一般指一级板块以上范围、100年以上时间、6级以上震级）范围内b值计算（陈培善等，2003）。本次研究所选区域较小，时间尺度为12年，地震目录99.8%震级小于5级，因此本次研究采取M_L震级进行计算。

地震活动性研究中，余震是否应该剔除尚未有统一答案。众多学者认为中强地震引发的余震会造成研究区域地震频次变化，因此在b值变化研究时，应对地震目录内的余震进行剔除（Keilis-Borok等，1980；Reasenberg，1985；陈凌等，1998；韩晓明等，2016）。也有学者认为余震去除会影响地震目录完整性，在研究时应保持地震目录的原始性，选取最小完整性震级以上地震进行计算即可（史海霞等，2018）。为研究余震序列对b值计算结果的影响，本文采用完整地震目录和去余震目录 2种地震目录研究b值的时空变化，为准确对比2种地震目录得到的b值差异，除地震目录外其它参数均相同，其中余震剔除时采用G-C方法，即对5级以上地震余震时间窗915天、余震空间窗半径 $ {R}_{0} $ 内的余震进行剔除（陈凌等，1998）：

2.2   地震b值计算方法

目前b值计算方法主要有2种，最小二乘法和最大似然法。最小二乘法是直接通过震级-频次的关系确定b值的大小，该方法简单易行，但其稳定性受地震频次影响，当地震较少时结果容易出现偏差（Weichert，1980；Sandri等，2007）。最大似然法作为概率方法具有更好的稳定性，可降低少数地震对区域整体b值结果的影响。由于研究区域台站间距较大，导致区域内最小完备震级较大，地震频次较少，因此，采用最大似然法计算本区域b值，该方法最早由Aki（1965）提出：

式中，M为震级；M_c为最小完备震级； $ \Delta M $ 为震级分档，本次地震目录震级分档为0.1。b标准误差可由以下公式进行估算（Aki，1965）：

式中，N为样本量，即地震目录中大于等于 $ {M}_{\mathrm{c}} $ 的所有地震的个数。可以看出 $ N $ 与 $ \delta b $ 成反比，因此提高样本数量可降低b值误差。该方法先将研究区域网格化，然后计算每一个网格的b值。计算b值时需要震级为连续随机变量，而实际震级计算时，会对震级进行归档，使震级只保留一位小数。因此震级分档 $ \Delta M $ 越小，b值的计算结果越准确，本次地震目录 $ \Delta M $ 为0.1。

2.3   最小完备震级

最小完备震级通过最大曲率法确定，该方法认为震级-频次曲线的一阶导数最大值对应的震级即最小完备震级（Wiemer等，2000），研究区域最小完备震级为M_L2.1（图3）。最小完备震级反映该地区的地震监测水平，青海地区虽地震频发，但人口稀少、台站密度较小，因此最小完备震级较大。为避免因地震监测能力差异导致地震目录不完整，影响b值的时空变化研究，本次研究将选取M_L2.1以上地震对玛多7.4级地震震源区进行b值研究。

图 3  研究区域2009年1月1日—2021年4月30日震级-频次图

Figure 3.  The frequency-magnitude distribution of earthquakes in the study area during January, 2009—April, 2021

下载: 全尺寸图片 幻灯片

3.   地震b值时空扫描特征

3.1   b值时间扫描特征

为研究时间段区域内b值随时间的变化规律，选用M_L2.1以上地震目录进行计算，为避免地震数量较少对结果准确性造成影响，以4个月为时间段进行b值时间扫描，分别计算余震剔除前后b值随时间的变化（图4）。由图4可知，2种地震目录得到的b值整体变化趋势相同，在5级地震发生前两者吻合的较好，但在5级以上地震发生后的一段时间内存在明显差异。图4中2条水平直线分别表示2种地震目录b值均值，其中完整目录b值均值为0.86、去余震目录b值均值为0.85，均值相差0.01。

图 4  b值时间曲线

Figure 4.  The temporal variation of b value

下载: 全尺寸图片 幻灯片

从b值时间扫描特征来看，是否剔除余震对b值时间变化影响不明显，仅在中强震发生后一段时间内存在差异。这种现象是由于研究区域中强震多表现为主震-余震型地震序列类型，余震主要集中在主震发生后，因此是否剔除余震对中强震发生后一段时间较为明显，对其他时间段影响较小。同时，虽然中强震余震持续时间较长，但随时间的推移余震不断减小，而该区域最小完备震级较大，不管是否剔除余震，震级小于M_L2.1的地震都被删除，后续较小余震对2种地震目录未产生影响，中强震发生后一段时间后2种地震目录的b值变化曲线又互相重合。因此，是否剔除余震对中强震余震窗内b值时间变化存在一定影响。

玛多M_S7.4地震前1年b值持续低于均值，并处于下降趋势，b值有上升趋势时地震发生，扩大到区域内其它5级以上地震，也发现此规律。同时，2种地震目录的b值曲线均存在3段较长的低b值时间段，且起止时间差异不大，分别为2009年5月—2011年12月、2014年5月—2016年4月、2019年9月至今，基本为中强震分割点，地震发生后b值明显上升，短时间内又下降到较低位置，并一直处于较低位置直至下次地震发生。这是由于地震发生前，该区域处于应力积累阶段，b值处于较低位置，地震发生后应力得以释放，区域内b值恢复到正常水平。

3.2   地震b值空间扫描特征

b值时间扫描结果表明，b值变化与中强震发生存在一定关系，可以作为地震预测时间上的参考。但研究区域较大，时间扫描特征仅能反映研究区域时间尺度上发生地震的可能，无法准确给出地震发生的空间区域。要实现通过b值分析给出研究区内中强震发生的空间预测范围，就要对b值进行空间扫描，分析研究区域内b值处于异常状态的地区。由b值时间扫描特征知，是否剔除余震对中强震余震窗内b值存在一定影响。为验证余震剔除是否对b值空间扫描特征同样具有影响，在b值空间扫描特征研究中同样计算M_L 2.1及以上震级去除余震地震目录及完整地震目录，分别计算2009年1月—2021年4月、2009年1月—2012年12月、2013年1月—2016年12月、2017年1月—2021年4月的b值空间特征，研究玛多地震前该区域b值空间变化特征及b值空间分布与中强地震的关系。玛多M_S7.4地震前去除余震目录及完整目录b值空间扫描结果分别如图5、图6所示。

图 5  去余震目录b值空间分布图像

Figure 5.  Spatial distribution of b value with deleting aftershock

下载: 全尺寸图片 幻灯片

图 6  完整目录b值空间分布图像

Figure 6.  Spatial distribution of b value without deleting aftershock

下载: 全尺寸图片 幻灯片

根据我国西部地震多发的特点，本文将研究区域经纬度扫描步长设为0.1°，并以此进行网格化（易桂喜等，2007）。由于研究区域地震较少，网格内地震扫描半径设为70 km，为保证结果的可信度，仅选取地震数量超过30个的网格进行计算，以每个网格点为单位计算b值，网格内的b值误差最大为0.18b，小于地震b值的变化范围，可信度较高。再根据每个网格点的b值做出研究区域b值空间分布，若网格点地震数量达不到要求则进行空白处理。

图5（a）、图6（a）分别为2009年1月1日—2021年4月30日去除余震目录及完整目录b值空间扫描结果，可以看出两者在玛多震中东南部低b值区域存在差别，图5（a）中该区域b值约为0.7，图6（a）中该区域b值约为0.8。区域内整体b值范围为0.5～2.0，玛多M_S7.4震中处于低b值区域内部，低b值区域以玛多M_S7.4震中所在的玛多-甘德断裂为中心向两侧扩展，该区域b值整体低于0.6，玛多M_S7.4震中所在区域为b值最低区域，b值约为0.5，表明该断裂震前构造应力处于较高状态。由图5（a）、图6（a）可知，b值的空间分布与中强地震的发生有较好的对应，为研究玛多M_S7.4地震前震中区域b值变化过程，本文将地震目录时间进行分割。

图5（b）、图6（b）分别为2009年1月1日—2012年12月31日去除余震目录及完整目录b值空间扫描结果，两者在称多西南部低b值区域存在差别，图5（b）中该区域b值约为0.6，图6（b）中该区域b值约为0.7。2010年玉树M_S7.4地震和2011年囊谦M_S5.2地震均发生于该区域，两者差异可能由是否剔除这2次地震余震引起。图5（b）、图6（b）中玛多M_S7.4地震震中位置空白，震中东部存在低b值区域，区域面积较小，b值在0.6附近。

图5（c）、图6（c）为2013年1月1日—2016年12月31日去除余震目录及完整目录b值空间扫描结果，两者没有差异。

图5（d）、图6（d）为2017年1月1日—2021年4月30日去除余震目录及完整目录b值空间扫描结果，两者在玛多震中东南部存在差别，差异区域与图5（a）、图6（a）一致。造成该差异的主要原因为2020年石渠M_S5.9地震发生于此区域，此次地震余震较多，去除余震会引起差异区域b值发生变化，因此该区域是否去除余震b值空间扫描结果呈现较为明显的差异性。

图5（b）、（c）、（d）显示随着时间推移，玛多震中附近低b值区域范围逐渐扩大，由正东不断向震中扩展。图5（b）中低b值区域位于东昆仑断裂带与玛多-甘德断裂带之间，图5（d）中低b值区域迁移至玛多-甘德断裂带，呈现南移现象。图5（d）中低b值区域贯通震中东南区域，该区域b值逐年降低，因此存在发生强震的危险性。

对比2种地震目录在不同时间段的b值空间扫描结果，整体差异不大，但在低b值区域存在差异。该区域去余震目录的b值相对较低，而完整地震目录b值较高，这将掩盖中强震区域的低b值特性。对比同一种地震目录在不同时间段的b值空间扫描结果，可发现玛多M_S7.4地震发生前低b值区域向震中的迁移现象，表明地震发生前震中附近应力不断集中。

4.   结论及讨论

本文根据中国地震台网提供的2009年1月1日—2021年4月30日地震目录，对玛多M_S7.4地震震源区震前地震活动性进行研究，分析区域内中强震发生前b值的变化特征。通过对完整地震目录和去余震地震目录分别进行b值时空扫描，发现是否剔除余震将影响中强震余震窗时间段内b值结果，空间扫描结果表明，在低b值区域完整地震目录b值高于剔除余震地震目录b值，从而掩盖了中强震区域低b值特性，而去余震地震目录可以较好凸显低b值区域，对中强震预测的指示考意义更强。因此，在进行b值时空变化研究时，尤其是通过空间扫描结果研究中强震的位置时，应剔除余震。对于文中结果整体差异较小的现象，主要是由于该地区最小完备震级较大，大量余震都处于最小完备震级以下，导致2种地震目录差异较小，因此所得结果差异也较小。虽然2种地震目录差异较小，但在较为关心的低b值区域仍存在一定差异，这将掩盖中强震区域的低b值特性。因此，后续研究可选则地震监测能力强、最小完备震级小的地区进行去余震地震目录与完整地震目录对比的研究，验证是否随着余震的增多，2种地震目录所得结果差异也随之增大。具体结论如下：

（1）玛多M_S7.4地震前1年b值开始低于均值，且不断下降，至b值有上升趋势时地震发生，扩大至区域内其它5级以上地震亦符合此规律，地震发生后b值明显上升，短时间内又下降至较低位置，并一直处于较低位置直至下次地震发生。这是由于地震发生前，该区域处于应力积累状态，b值处于较低位置，地震发生后应力被释放，区域内的b值即恢复到正常水平。

（2）从b值空间扫描结果看，玛多M_S7.4地震前震中位于低b值区域中心，该位置为b值最低处。

（3）不同时间段的b值空间扫描结果表明，玛多M_S7.4地震发生前低b值区域向震中逐渐迁移，表明地震发生前应力向震中附近不断集中，因此低b值区域的迁移现象可为地震预报提供参考。

（4）余震剔除可以有效避免中强地震余震序列对b值结果的干扰，b值空间扫描时完整地震目录掩盖了中强震区域的低b值特性，而去余震地震目录可以较好凸显低b值区域，对中强震预测的指示意义更强。因此，在b值时空扫描研究时，应剔除余震，力求更准确的表现ｂ值时空变化特征。