0 引言

地震预测是世界公认的科学难题，因地表观测手段受限、地壳构造复杂且地震发生规律难重复，导致地震预测是小概率事件且是经验性的，但科学家从未停止对其探索（周硕愚等，2019）。地震预测有长中短临渐进式预测，其思路为：长期预测（未来10年）、中期预测（未来1—3年）、短期预测（未来3个月）和临震预测（未来10天）（江在森等，2020）。本文重点关注短临预测。短临预测核心是寻找突变异常，如前兆异常突跳、平静等。当地震前兆b值出现异常时，可在一定程度上预测地震的发生（曾宪伟等，2020；姜佳佳等，2021；Xie et al，2022）。

古登堡—里克特（G-R）震级—频度关系式为：lgN = a - bM（Gutenberg et al，1942），式中N为震级M≥M_c（最小完整性震级）的地震累计次数，a表示区域内地震活动水平，b表示区域内不同震级地震的比例关系。20世纪60年代，Mogi（1962）和Scholz（1968）开展岩石破裂实验，发现b值大小与介质应力水平、岩石性质有关。Mori等（1997）、Amitrano（2003）、Schorlemmer等（2005）通过岩石实验进行了验证。介质的非均匀程度、有效剪应力等参数会影响b值变化，b值越大，介质的不均匀性越显著（李纪汉，1987；Wiemer et al，1998）。b值与应力成反比关系，b值越低，表明应力越高（Scholz，1968；Wyss，1973；Urbancic et al，1992），低b值区具有更高应力累积，特别是岩石破裂前或地震发生前，b值下降趋势明显（刘子璇等，2020；Li et al，2024）。

在我国，李全林等（1978）、马鸿庆（1978）最早通过b值变化预测地震的发生。当某一区域出现b值下降时，反映介质内部应力较高，发生大破裂的可能性增加，发生强震概率高于其他区域。b值时空扫描结果可用来判定未来强震的活动水平，揭示和推断活动断裂在不同阶段应力积累的相对水平，判定活动断裂的强震危险性。强震前b值的时间变化是地震预测的重要前兆，但事实上，强震前b值时间变化过程复杂，并非简单的上升、下降，说明b值异常作为地震前兆具有多样性（孟昭彤等，2021）。本文基于川西地区地震目录，研究强震前后b值变化，并推测该区未来地震发生的可能性。

1 研究背景

川西地区（100°—103.1°E，29°—32°N）位于四川省西部、青藏高原东缘，地处巴颜喀拉块体和川滇块体，区内断裂构造发育，主要分布NE—NNE向、NW—NNW向和近NS向3组（徐锡伟等，2008）断裂带（图 1），其中鲜水河断裂带和龙门山断裂带是四川省构造活动活跃且具有潜在危险性的断裂带。区内地震活动强烈，据四川地震台地震目录，2010—2023年，川西地区共发生M_S≥5.0地震20次（M_S≥6.0地震4次），为该区中强地震前后b值变化研究提供了详实数据。

b值空间变化可用来判别强震危险区（易桂喜等, 2005, 2008）。以2013年4月20日雅安M_S 7.0地震、2014年11月22日康定M_S 6.3地震、2022年6月1日芦山M_S 6.1地震和2022年9月5日泸定M_S 6.8地震（图 1）为例，获取震前4—5年及震后8—9月地震目录，分析川西地区大震前后b值变化，并鉴于短临预测，选取4次强震发生前3个月的地震目录，根据b值空间变化特征圈定潜在地震危险区，分析强震是否发生在该危险区，以验证划定潜在地震危险区规律的有效性。

2 资料及方法 2.1 资料选取

据四川地震台2010年1月1日—2023年10月31日地震目录，川西地区共发生M≥0.1地震68 887次，其中M_S≥6.0地震4次，M_S 5.0—5.9地震16次，M_S 4.0—4.9地震100次，M_S 3.0—3.9地震765次，M_S 2.0—2.9地震5 950次，M_S 1.0—1.9地震32 289次，M_S 0.1—1.0地震29 763次（下文中震级标度统一使用M表示）。采用Uhrhammer（1986）提出的基于地震事件时间窗、距离窗的算法，去除前震和余震，获得所需主震目录，并基于完整性震级范围法（EMR）（Ogata et al，1993），确定最小完整性震级M_c。

M_c是指地震目录中该震级以上地震是完整且无遗漏的，是进行地震目录统计的基本参数（谢卓娟等，2019）。Woessner等（2004）、韩立波等（2012）分析发现，完整性震级范围法（EMR）在稳定性和可靠程度上优于MBS（稳定b值法）（Cao et al，2002）、GFT（拟合优度法）（Woessner et al，2004）等方法。在EMR方法中，对于M_c的不确定度δM_c使用自举（Bootstrap）方法的蒙特卡罗近似来估计。文中在计算b值时，将Bootstrap选为200。

川西地区4次6级以上地震频度—震级分布见图 2，图中矩形表示累积频度—震级分布，三角形表示非累积频度—震级分布。通过计算其分布曲线的一阶导数，将一阶导数最大值，即曲率首次发生突变点作为M_c的值，得到雅安M 7.0地震序列M_c = 1.6（选取2010—2013年地震目录）、康定M 6.3地震序列M_c = 1.5（2010年1月1日—2015年6月1日地震资目录）、泸定M 6.8地震序列M_c = 1.3（2018年1月1日—2023年4月1日地震目录）、芦山M 6.1地震序列M_c = 1.3（2018年1月1日—2022年12月31日地震目录）。

2.2 b值计算方法

b值计算方法有多种，如归一化和标准化方法（段华琛等，1995）、稳健回归方法（杨马陵等，1999）、最大似然法（Maximum likelihood method，MLE）（Aki，1965；Utsu，1965）和最小二乘法（Least square method，LS）（Zöller et al，2002），其中最小二乘法和最大似然法较为常用。

最小二乘法以离散点到拟合直线的距离为权重，每个震级档之间的地震数据权重是不同的，受样本扰动大，当参与统计的地震数目较少时，受个别大地震随机性的影响大，其更重视含有丰富信息量较大地震的作用，低震级的地震权重占比小，换而言之，出现几个较大地震就会使得b值波动较大，其计算公式（Okal et al，1995）如下

$ b=\frac{n \sum M_i \lg N_i-\sum M_i \lg N_i}{\left(\sum M_i\right)^2-n \sum M_i^2} $

(1)

式中，将所有地震分为n档，M_i（i = 1, 2, …, n）包含M_i -ΔM/2≤M≤M_i + ΔM/2的所有震级档，在i档中M_i≤M的地震数目记为N_i。

最大似然法将每个地震事件的贡献看作是相同的，对所有地震的震级用同样权重求平均，也就是将大、小地震的作用视为相同，其计算公式（Ustu，1966）如下

$ b=\frac{\operatorname{lge}}{\bar{M}-(\Delta M / 2)} $

(2)

式中：M_c为最小完整性震级；ΔM为震级分档，文中震级间隔取0.1；$\bar{M}=\left(\sum\nolimits_{i=1}^N M_i\right) / N$，为震级平均值；e为自然常数。

采用上述2种方法计算b值，无论样本量大小，基于最大似然法计算的b值平均值均高于最小二乘法的计算结果，且随着样本量的增加，b值平均值均趋向于理论值，表明最小二乘法适用于计算区域整体b值，也就是适用于样本较为充裕的地震目录，但最大似然法在样本量较小时即能得到准确的估计值（谢卓娟，2020）。文中选取的地震目录样本量较小，且低震级的地震权重占比大，因此采用最大似然法计算研究区b值。

3 地震活动演化

应变能的积累和释放过程可用来描述研究区内地震活动的演化过程。应变释放曲线图即$\Sigma \sqrt{E}-T$关系曲线，可由震级能量公式lgE = 4.8 + 1.5M（赖智华等，2023）求得（震级已知）。基于所选地震目录，计算得到川西地区2010—2023年应变释放曲线，见图 3(a)，可知在4次强震孕育到发生前的时间段内，应变能增长平缓，处于能量积累阶段，表明在孕震区附近地震活动性往往表现出相对平静的特征，即处于相对平静期，强震发生后，应变能得到快速释放，处于地震活动活跃期，之后相对平静，进入未来强震孕育阶段。地震活动的活跃期和平静期通常称为地震活动周期（卢显等，2023）。

震级—时间分布图（M—T图）是描述地震活动时间进程的最简单直观的方法。根据M—T图显示的地震疏密程度和大小来划分地震活动的相对平静期和活跃期（蒋溥等，1993）。如图 3(b)所示，研究区4次强震（图中序号①、②、③、④）发生前一段时间内，均处于地震活动相对平静期，与图 3(a)中的平静期相对应，并在2022年泸定地震（图中序号④）发生后，进入地震平静期，即未来地震的孕育阶段。结合上文所得结论，推测该区域地震活动水平目前仍处于活跃期。

4 研究区b值时空特征 4.1 b值时间分布特征

为全面掌握川西地区b值变化，选取2010年1月1日至2023年10月31日地震目录（包含4次强震震例），采用200次的bootstrap重采样，以500个样本为一个时间窗，0.1为震级间隔，计算得到研究区b值分布在0.553—1.134，平均值为0.879，结果见图 4。

结合图 3(b)可知，研究区地震以小震为主、中强震较少，且中强震多在b值小于平均值的时段发生，小震多在b值处于平均值以上时段发生。4次强震发生前，b值均突然下降，后逐步升高。而且，2022年先后发生芦山、泸定2次6级以上地震，间隔时间仅3个月左右，受芦山地震影响，泸定地震b值并未表现出明显的时空变化规律。b值时间扫描结果表明，b值的上升、下降变化与中强地震发生有一定对应关系，可作为地震预测时间参考。

4.2 b值空间分布特征

通过b值空间扫描，可判别地震发生的潜在危险区，同时对研究区内b值处于异常状态区域做出解释（王辉等，2012；曾宪伟等，2020；孟昭彤等，2021）。在进行b值空间扫描时，应选择适当的研究范围，避免b值异常被遗漏。2013年雅安M 7.0地震震源深度13 km，震中（30.3°N，103°E），取（102°—103°E，30°—31°N）范围作为研究区域；2014年康定M 6.3地震震源深度16 km，震中（30.26°N，101.69°E），取（101°—102°E，30°—31°N）范围作为研究区域；2022年芦山M 6.1地震震源深度17 km，震中（30.37°N，102.94°E），取（102°—103°E，30°—31°N）范围作为研究区域；2022年泸定M 6.8地震震源深度16 km，震中（29.59°N，102.08°E），取（102°—103°E，29°—30°N）范围作为研究区域。

选取4次强震发生前4—5年的地震目录，采用0.01°×0.01°的空间间隔进行网格化，以每个网格节点的中心点为圆心，选取半径为50 km的圆形单元内M≥M_c的地震（M雅安≥1.6，M康定≥1.5，M泸定≥1.3，M芦山≥1.3），采用最大似然法（MLE）进行b值空间扫描，为保证样本足够多，少于50个地震的网格不参与计算，以空白区显示，得到b值空间分布图，见图 5，图中不同颜色代表不同b值，表示应力水平不同。

从b值空间扫描结果可知，4次强震震源区及邻近范围处于较低b值区域，且b＜1，与Tormann等（2015）的研究结论一致，表明该区域为高应力积累区，地震危险性较高。此外，低b值区与地震危险区的判定，在一定程度上依赖于对区域活动断裂的认识程度。图 5(a)与图 5(c)选择的区域范围相同，潜在地震危险区即震中所在的绿色区域，该区域有3条活动断裂穿过。活动断裂带上具有高应力积累的凹凸体或闭锁断裂段，是强震易发地（Aki，1984；Wiemer et al，1998；Wyss et al，2000），可用来判定潜在地震危险区。其中：凹凸体可以理解为地震发生前断层面上应力集中、强度较大的区域，在同震过程中将是破裂起点或者破坏最严重的地方（许才军等，2022），闭锁断裂段是应变能易积累或应力集中地段。由图 5(c)可见，芦山地震发生在低b值区域，即发生在潜在地震危险区，验证了本研究对潜在地震危险区划分的有效性。中长期异常有助于地震发生规律的经验总结，但这些规律往往不够直观和明确，此时更应关注短临异常特征。

以上4次强震余震目录较少，在进行网格化时，将扫描半径缩小到30 km，不足30个地震的网格不参与b值计算，以空白区显示，其他参数设定与4次强震发生前4—5年地震目录参数相同，b值空间分布见图 6。结合强震前b值空间分布，可知在不同强震前后，震源区及邻近区域b值变化规律并不一致，其中：①雅安地震：发震前后b值无明显变化；②康定地震：震后b值上升；③芦山地震：震后b值上升；④泸定地震：震后b值下降。

震前4—5年地震数据仅能提供b值变化的大致范围，而震前3个月的数据则能反映更精确的b值变化。在发震前的短临阶段，因研究区域地震目录较少，在进行网格化处理时，将扫描半径缩小至10 km，不足10个地震的网格不参与计算，其他参数设定与前面设定相同，得到4个强震发生前3个月的b值空间分布，结果见图 7。

将4个强震孕育、发生不同阶段的b值空间分布进行对比，可知：与震前4—5年震源区b值空间分布（图 5）相比，临震阶段的b值（图 7）进一步下降，呈低b值变化，说明区域地震活动进入孕育阶段，震后震源区及邻近区域b值上升（图 6），呈较高b值状态。综上可知，b值的空间变化与时间变化结果吻合，均具有震前降低、震后升高的特征。需要注意的是，在进行b值研究时，要选择合适的研究区域及地震目录范围，以避免b值异常区不明显或被遗漏。

5 结论与讨论

川西地区地震活动强烈，为分析b值变化提供了良好的数据基础。基于四川省地震台提供的2010年1月1日至2023年10月31日的地震目录，分析该区发生的4次显著强震，即2013年4月20日雅安M_S 7.0地震、2014年11月22日康定M_S 6.3地震、2022年6月1日芦山M_S 6.1地震和2022年9月5日泸定M_S 6.8地震前后b值变化，得出以下结论：

（1）通过应变释放曲线图和时间分布图，可知4次强震从孕育到发生的前一段时间内，应变能增长平缓，处于应变能积累阶段，表明孕震区附近地震活动性往往表现出平静的特征，可将该阶段划分为相对平静期。地震的发生使得应变能快速释放，随后进入平静阶段，为下一次地震的孕育发生储备能量，以此循环往复。因此，地震平静现象可作为地震前典型的前兆异常（梅世蓉，1960）。

（2）b值时间分布显示，该区域b值分布在0.553—1.134，平均值为0.879，研究区内b值波动幅度较大，尤其是4次强震发生前后，地震活动性高，b值均呈突然下降，达低峰值后逐步升高到高峰值的异常特征。

（3）对4次强震发生前4—5年地震目录进行b值空间扫描，可知震源区及邻近区域b值呈较低值状态，表明该区域为高应力积累区，地震危险性较高；震前3个月的空间扫描结果显示，b值进一步下降，呈低b值状态，说明区域处于地震孕育阶段；震后b值空间扫描结果与震前4—5年数据相比，变化规律不明显，与震前3个月数据相比，b值呈上升特征，结合b值随时间的变化规律，可知震后b值呈上升趋势变化。

（4）要根据情况选择合适的地震目录参与b值计算，通过多次研究确定所需数据，并选择合适的研究方法，确保b值结果真实、可靠。而且，务必选择合适的研究区域及地震目录范围，避免b值异常区不明显或被遗漏。

（5）川西地区b值计算结果表明，该区域目前处于较高应力状态，具有发生中强震的危险，需密切关注区域地震活动。今后，将考虑把b值与其他地震活动性参数以及其他学科多种手段相结合开展综合研究，以便在短临预测、地震活动性和地震危险性分析中发挥更重要的作用。

基于MATLAB的Z-MAP软件包由瑞士地震服务中心提供，在此表示感谢。