【www.zhangdahai.com--其他范文】
陈丽娟 陈学忠 李艳娥 龚丽文
1)重庆市地震局，重庆 401147
2)中国地震局地球物理研究所，北京 100081
Gutenberg等(1944)提出了重要的地震频度-震级定律，即logN=a-bM，其中M为地震震级，N是震级≥M的地震累积频度，式中的b值是反映地震活动性的重要参数。Suyehiro(1966)首次发现了前震b值比余震低的现象，并认为低b值异常为地震前兆，其他研究人员也观察到同样的现象(Knopoffetal.，1992；
Caoetal.，2002)。根据岩石模拟实验的研究结果，Mogi(1962)认为介质的不均匀性是导致b值变化的主要因素，Scholz(1968)则认为b值及其变化是由岩石所处的应力状态决定的。中国一些学者的系列岩石实验研究表明b值与应力状态和介质性质均相关，并不由单一的因素决定(吴开统等，1986)。有研究认为b值的变化可用来反映断层上的应力状态(Goebeletal.，2012；
易桂喜等，2013)，b值的降低被解释为主震来临前应力的增加(Wyss，1973；
Mainetal.，1989；
Urbancicetal.，1992；
Hainzletal.，1999)，在一些大地震前b值随时间推移呈下降变化(Nuannin，2005；
Nakaya，2006；
Nanjoetal.，2012；
薛艳等，2012)。
2008年5月12日四川省汶川(31°N，103.4°E)发生了MS8.0 地震，地震由汶川映秀开始破裂，沿龙门山断裂带向NE扩展，形成总长约240km的地表破裂带(徐锡伟等，2010)。学者们针对汶川地震前的b值开展了大量研究，利用的方法有最大似然法(Zhaoetal.，2008；
王辉等，2012；
刘艳辉等，2015；
史海霞等，2018；
Chenetal.，2020)、最小二乘法(易桂喜等，2011)、改进的最小二乘法(刘雁冰等，2017)和改进的贝叶斯方法(Zhangetal.，2015)等，得到的结论也不尽相同。不同学者通过空间扫描得到的汶川地震前的低b值区有较大差异，如汶川地震起始破裂点附近的西南区域(Zhaoetal.，2008)、龙门山断裂带中北段的绵竹—茂县与江油—平武段(易桂喜等，2006，2011)、映秀附近的彭灌杂岩区(刘雁冰等，2017)等，而王辉等(2012)给出的大面积空间扫描结果显示汶川地震前龙门山区域为高b值背景。大部分学者的时序分析结果认为汶川地震前b值有明显的下降变化，但下降的形态不一。解朝娣等(2010)的结果显示2006—2008年b值稳定波动，在2008年初急剧下降至最低值。史海霞等(2018)的结果认为b值从2002年初开始出现趋势性降低，地震前半年快速显著下降。刘艳辉等(2015)计算的b值在汶川地震前呈现正常—低值—高值—回落—发震—正常的变化。一些结果显示b值在震前逐渐降低，在汶川地震前达到最低值(Zhangetal.，2015；
刘雁冰等，2017)。然而也有结果显示汶川地震前b值的下降并不明显(Zhaoetal.，2008)，甚至在临震前有所上升(Chenetal.，2020)。这些结果主要是基于构造带、震源区、主破裂区或龙门山地区选取研究区域进行b值的时序分析。
研究区域的选取在b值时序分析工作中十分关键，以上研究主要是基于汶川地震的发震位置和地震破裂区域选取震中周边的区域分析b值的时间变化，选取的区域中可能包含与地震孕育过程无关的地震样本，这些无关的样本会影响统计特征，导致真正的孕震信息被忽略。此外，空间扫描出的低b值区是研究时段内区域的绝对b值，而忽略了这个低b值特征可能是由区域本身介质的强度或构造差异导致的，与应力的动态变化过程无关。大地震可能发生在低b值区，但大地震前的绝对b值也不一定明显偏低。例如，2004年苏门答腊9.0级地震发生前震中区的b值在1.0以上(Nuannin，2005；
Nanjoetal.，2012)；
2011年日本9.0级地震前b值为0.96(Nanjoetal.，2012)；
汶川地震前b值为1.08(史海霞等，2018)。这些大地震发生前绝对b值并没有出现明显的低值异常，但在震前均出现明显的下降变化过程。b值下降说明应力水平上升，反映了孕震的动态演化过程。本文中，我们关注的是b值的下降过程和下降幅度的空间分布。基于b值与应力的负相关关系，只有b值出现下降变化的区域才是应力水平上升的“真实”异常区，那么就要考虑研究时段内b值的下降幅度，这个量越大，说明应力上升的幅度越大。汶川地震前，较大范围内的应力是如何积聚的？是否只有龙门山断裂带的b值出现下降变化过程？如何在地震发生前在较大范围内准确地判断出未来强震的孕育区？以上都是本文想要探讨的问题。着眼于b值的下降变化过程，本研究通过格点搜索法，在较大空间范围内找出b值出现下降变化的区域，分析不同幅度的b值相对变化量对应的区域与地震震中的关系，进而分析出现b值下降变化区域内的b值随时间的变化特征，通过b值的时空演化直观反映汶川地震前应力的积累演变过程。
本文以汶川地震的起始破裂点为中心选取研究区域(26°～36°N，98°～108°E)，包含整个四川地区及部分邻省地区，远大于汶川地震的震源区，如图1 所示。图中，蓝色实心圆为2000年1月—2008年4月期间发生的ML≥1.0地震，绿色实线为汶川地震的主破裂带。这些地震的G-R关系如图2b 所示，目录完整性震级可达ML1.5，考虑不同时段Mc的差异，利用最大曲率法计算了区域内Mc随时间的变化，如图3 所示。最终选取研究区内ML≥1.7的地震目录作为b值格点搜索的研究资料。图2a 为震级随时间的变化。
图1 研究区及2000年1月—2008年4月ML≥1.0地震分布图Fig.1 The study area and earthquakes(ML≥1.0)from Jan.2000 to April 2008.蓝色实心圆为地震，红色震源球为汶川 MS8.0 地震，震源机制解来自GCMT(Global Centroid-Moment-Tensor)目录，绿色实线为汶川地震的主破裂带(修改自史海霞等，2018)，白色实线为断层
图2 研究区域内M-T图(a)和G-R关系图(b)Fig.2 The chronological orders of earthquakes(a)and magnitude-frequency relation(b)in the study area.
图3 研究区域内Mc随时间的变化Fig.3 Variation of Mc with time in the study area.
本文采用最大似然法计算b值(Aki，1965)，计算公式为
(1)
b值95%置信度的标准差为
(2)
将研究区域(26°～36°N，98°～108°E)划分为0.1°×0.1°的格点，以每个格点为中心点，选取1°×1°的空间子区，计算每个子区内的b值随时间的变化曲线，计算时使用等事件数的样本窗口，对窗口以相等的事件数为步长进行滑动。挑选出震前b值出现明显下降变化的子区，并将该子区的中心点作为一个异常点，多个异常点集中的区域可看作b值异常区。然后选取这个异常区内的地震计算b值随时间的变化，从而得到整个b值异常区的范围及其b值随时间的变化特征。
b值的计算与样本量的大小有重要关系。李世杰等(2018)的研究表明，最大似然法需要的最低样本量为200。为保证b值的可靠性，本文中对每个子区以200个地震样本作为1个时间窗计算出1个b值，以10个地震为步长向后逐步滑动样本窗，计算出每个子区的b值随时间的变化曲线。选取时间曲线大于10个数据点的子区，共计得到2410个子区的b值时间曲线。图4 中，灰色实心圆即为计算了b值时间曲线的格点，图中蓝色空心圆为ML≥1.7的地震，可以看出计算了时间曲线的格点分布与地震的分布有关，地震分布较稀疏的区域由于样本量不足，未计算b值的时间曲线。
图4 计算了b值时间曲线的格点分布与所选取地震的震中分布Fig.4 The distribution of selected earthquakes and the grids calculated b-value time curves.蓝色空心圆为ML≥1.7地震，灰色实心圆为计算了时间曲线的格点，红色震源球为汶川地震，绿色实线为汶川地震的主破裂带
我们通过程序在每个子区的b值时间变化曲线上取定 2002年1月—2008年4月间的n个b值数据，对于相邻2个b值，计算后一个b值与前一个b值之差，得到n-1个b值差，再计算出最后一个b值与第1个b值的相对变化量Δb，代表b值的下降幅度。设定挑选b值呈下降变化形态的条件为：1)n-1个b值差值中有55%为负值，即多数b值呈下降变化；
2)b值的下降幅度Δb分别≥5%、10%、15%、20%、25%和30%。
3.1 b值下降幅度的分布特征
当b值的下降幅度取不同值时，其对应的格点分布有何特征？图5a—f展示了b值的下降幅度Δb分别≥5%、10%、15%、20%、25%和30%时b值异常格点的空间位置，即图中的彩色圆点，分别有633、444、284、177、99、39个格点。图5a 显示，当Δb≥5% 时，异常格点分布在龙门山断裂带、安宁河断裂带、马尔康附近以及云南北部，范围较广，最大的异常区位于龙门山断裂带。随着Δb的增大，异常范围逐渐减小，当Δb≥20% 时，异常范围集中收敛至龙门山断裂带及马尔康地区。图5f 显示，当Δb≥30% 时，共有39个异常格点，但其中的29个异常格点分布在龙门山断裂带的都江堰—北川段，该段位于汶川地震主要破裂区的南段，汶川地震的起始破裂点位于异常区的东南缘。综上，汶川地震前b值下降幅度最大的地区位于汶川地震主破裂区的南段。
图5 不同b值下降幅度Δb对应的异常格点分布Fig.5 The grid points distribution with decreasing b value obtained by grid point search method.绿色实线为汶川地震的主破裂带，彩色实心圆为所在子区内b值时间曲线出现下降变化形态的格点，图a—f为b值下降幅度Δb分别≥5%、10%、15%、20%、25%和30%时的格点分布
图6 b值异常区的选取(a)及异常区内b值的时间曲线(b)Fig.6 The abnormal area(a)and the time curve of b-value(a).彩色实心圆为异常格点，蓝色线段为选取的异常区，蓝色震源球标示汶川地震的起始破裂点，绿色线段为汶川地震的主破裂带
3.2 异常区内的b值随时间的变化特征
选取比b值下降幅度Δb≥25% 的格点分布区稍大一些的区域作为汶川地震前的b值异常区，即图6a 中蓝色线段所围区域。根据已有研究可知，汶川地震震源区的最小完整性震级可达ML1.5(史海霞等，2018)，本文选取图6a 中蓝色线段所围区域内的ML≥1.5地震，以500个地震样本作为1个时间窗计算出1个b值，以10个地震为步长向后滑动样本窗，得到异常区内的b值随时间的变化曲线，如图6b 所示，图中灰色部分为b值95%置信度的标准差。可以看出在汶川MS8.0 地震前，异常区内的b值时间曲线呈3个阶段的变化过程：第1阶段b值持续下降，自2002年5月—2005年末从1.5下降至1.1，在约3.5a内的下降幅度约为27%；
第2阶段b值较为平稳，2006—2007年末b值保持在1.1左右，维持了约2a的平稳期；
第3阶段b值快速下降，2008年初b值在临震前从1.15下降至1.05，在短短4个月时间内达到约9%的下降幅度，之后发生了汶川MS8.0 地震。b值在震前短时间的快速显著下降过程与前人得到的结果(解朝娣等，2010；
史海霞等，2018)相似。
4.1 b值下降异常区与震源区的关系
不同的b值下降幅度对应的异常格点分布表明，在汶川地震前研究区内不止龙门山断裂带的b值出现了下降变化，但龙门山断裂带的b值下降幅度最大。地震最终发生在b值下降幅度最大的区域，说明较大地震前应力的加载并不局限于震源区，而是在较大范围内呈现上升变化，龙门山断裂带在汶川地震前承受了最大的应力加载，导致其最终破裂。汶川地震的破裂过程反演结果显示，最大的滑动量集中在汶川—映秀一带下方(张勇等，2008)。震源破裂过程的研究表明，汶川地震由可分辨的5次7.3级以上地震组成，这5次子事件的空间展布恰好与本文得到的b值异常范围吻合(赵翠萍等，2009)。马瑾等(2013)认为汶川地震的孕震区主要分布在映秀-北川断裂和后山断裂之间，孕震区的范围小于发震区。以上研究得到的汶川地震的最大滑动量、主要破裂范围、孕震区等与本文得到的b值下降幅度最大的区域较为相似，说明b值下降幅度较大的异常区能够反映强震的孕育范围。这为强震发震地点的判断提供了一种思路，即b值下降幅度最大的区域未来发生强震的可能性较大。
图7 2006年1月—2008年4月的b值空间扫描结果Fig.7 Spatial scanning distribution of b-value from Jan.2006 to April 2008.
4.2 b值下降异常区与b值空间扫描结果的对比
b值下降幅度最大的区域是否与空间扫描出的低b值区一致？通过图6b 可知，b值下降异常区内的b值于2006年下降至较低水平，我们将研究区域同样划分为0.1°×0.1°的格点，以每个格点为中心点，选取0.5°×0.5°的空间子区域，对2006年1月—2008年4月的地震目录进行空间扫描，计算每个子区内的b值，若子区内地震<50个，则不计算b值。图7a 为b值的空间扫描结果，可看出在2006年1月—2008年4月期间，研究区内甘川交界、道孚—石棉段、鲁甸—盐津等区域均表现为低b值。图7b 为去掉b≥1的结果，可看出都江堰—北川段只有少部分区域b<1，并没有表现出明显的低b值异常，在北川附近有少许显著的低b值异常，对比图5 可发现都江堰—北川段是b值下降幅度最大的区域。王辉等(2012)的结果也表明汶川地震震源区的b值高于整个川滇地区的背景，汶川地震前在整个川滇地区有很多低b值区，但在龙门山断裂带看不出明显的低b值异常，其统计结果也显示川滇地区大部分强震与低b值区域的关系并不显著，说明低b值并非大地震发生的必要条件，但强震前b值出现下降变化，则反映应力的动态演化过程，这种动态的应力积累是有利于发震的必要条件。在强震发生前如何在较大范围内找出未来发震的地点是地震研究者广为关注的问题，而常规的b值空间扫描结果很难找到这个未来孕震区。本文通过搜索b值出现下降变化的格点，找到下降幅度最大的区域，并发现其与未来的震中破裂区重合。因此，利用格点搜索法找出b值下降幅度最大的区域，可为强震的震中判断提供有益参考，格点搜索法的意义就在于此。
4.3 b值时序变化的物理意义
b值在趋势下降的前3.5a中下降了约27%，下降速率为每年7.7%。之后经2a的平稳变化状态，又于临震前4个月在平稳变化的背景上继续下降了约9%，下降速率为每年27%，临震前的下降速率较快。这种变化形态的物理意义是什么？岩石力学及声发射试验结果表明，岩石内部大尺度裂纹的比例增加时b值减小，微裂纹的破坏状态比较稳定时b值几乎不变(曾正文等，1995；
葛振龙，2019)，在岩石临近破坏时b值会大幅度下降(王春来等，2018；
赵建军等，2019)。低围压下声发射b值出现骤降特征为岩石破坏的前兆，而高围压下岩石破坏前b值并无明显的前兆特征(张黎明等，2015)。结合以上研究结果分析，汶川地震前b值随时间呈现“长时间持续下降—平稳变化—临震前再次下降”的变化形态或许反映了地下岩石裂纹发展的3个阶段：第1阶段内，在构造应力的作用下不断形成新的破裂，b值持续下降；
第2阶段内，构造应力继续增加，岩石内部的裂纹呈现渐进式稳定扩展，大、小尺度破裂事件的比例基本不变，b值基本不变；
第3阶段内，b值快速下降进入临震阶段。汶川地震前b值在临震前的快速下降或许能为短临预报提供参考。
本文采用格点搜索法将研究区域划分为0.1°×0.1°的格点，以每个格点为中心点，选取1°×1°的空间子区域，计算每个子区的b值随时间的变化曲线，并基于给定的挑选条件，选出b值呈现下降变化形态的子区，将其认定为汶川MS8.0 地震前的b值异常区，对异常区内的b值时序变化进行分析，得到以下结论：
汶川地震前在初始破裂点NE向形成了大面积的b值下降异常区。在空间分布上，当设定b值下降幅度Δb≥30% 时，异常区集中在都江堰—北川段，与地震孕震区吻合。异常区内的b值随时间呈现出“趋势下降—平稳变化—临震前快速下降”的变化形态，可能反映了地壳岩石裂纹发展的阶段性特征，与岩石声发射试验结果具有一定的相似性，临震前b值的短期快速下降可作为发震前兆。
格点搜索法的优点在于能在强震发生前，在大范围内搜索出b值具有下降变化形态的格点，根据选取的下降幅度搜索出应力积聚最强的区域，从而判定未来强震的孕震区。在实际的预测研究工作中，将b值和其他学科相结合，可做为判断孕震区的有益参考。
致谢中国地震台网中心为本研究提供了地震目录；
文中部分图件利用Generic Mapping Tools(GMT)绘制；
审稿专家为本文提出了诸多建设性的修改意见。在此一并表示感谢！
