唐山MS7.8地震前b值异常特征

陈学忠，李艳娥，陈丽娟

1 中国地震局地球物理研究所，北京 100081 2 重庆市地震局，重庆 401147

0 引言

1976年7月28日在河北唐山发生的MS7.8巨大地震给国家和人民造成了巨大灾难.地震预报，特别是短临地震预报，可以显著减轻地震造成的人员伤亡和财产损失.但是，由于科技水平的局限性、对地震产生的成因以及震前观测到的一些异常现象的认识，还有待深入.至今，地震预报仍然难以实现.对已发生的强震进行解剖，深入研究强震发生前的异常现象，将有助于深化对地震孕育发生过程的认识，同时对推进地震预报深入研究也有裨益.唐山MS7.8地震发生之后，不少研究者对唐山地震前出现的地震学方面的异常进行了分析研究(朱传镇等，1977,1981；李全林等，1978；林怀存等，1987；钱兆霞，1988；薛艳和梅世蓉，1999；隗永刚等，2018；李赫等，2020)，研究内容主要涉及区域地震活动图像、震群活动、b值、微震震源参数以及地震活动与地球自转的关系等.关于唐山地震前b值的研究，李全林等(1978)研究了唐山震区及其邻区b值随时间的变化，发现震前b值存在下降过程，同时还对京津唐渤张区域进行了b值空间扫描分析，发现低b值区由西向东迁移的过程.目前，在b值分析时，研究者们主要在两个方面开展工作：一个是特定时段的空间分布分析，另一个是选定空间范围内的时间变化分析.b值的空间分布可以反映各个地区b值的差异性，仅仅根据b值的空间分布难以辨别出某一地区b值随时间的变化情况.在b值的空间分布图上，常常出现b值高低区域交错分布，无规律可循，难以将其与地震孕育相联系.同样，在选定某空间范围进行b值随时间变化分析时，选取区域时没有合理的标准，具有任意性，难以保证所选取的区域完全包含了b值异常区，因而不能反映全部b值异常区的b值随时间的变化情况.假设地震发生前在震源附近的一个空间区域内存在b值异常变化，如果选取的空间范围位于这个空间区域以外，与地震孕育过程无关，这时，b值随时间的变化与地震的孕育之间难以建立有效的联系；如果选取的空间范围只是包含了这个空间区域的一部分，此时b值随时间的变化也不能全面反映地震的孕育过程.一些震例研究表明，地震发生前b值存在数月-数年的下降变化过程(Imoto,1991; Jaumé et al.,1999; Enescu and Ito,2001; Nuannin et al.,2005; Chan et al.,2012; Nanjo et al.,2012; 史海霞等,2018).我们着眼于b值下降变化过程，通过格点搜索把b值下降变化的区域搜索出来，然后再分析这个区域内b值随时间的变化特征，就可以全面地揭示b值异常的空间分布和时间变化特征.本文将根据上述思路对唐山地震前的b值异常特征进行研究.

1 资料和研究区域

图1a给出了1976年唐山MS7.8地震及其余震的震中位置和研究区域，取图中所示的经纬度范围(37.5°N—42°N，113.5°E—121°E)为格点搜索范围，即研究区.这个区域显然远远大于唐山MS7.8地震余震分布区域.对研究区内于1970年1月—1976年6月期间发生的地震，根据G-R关系得到的地震目录完整性震级可达ML1.5.因此，在本文分析中，选取研究区内1970年1月—1976年6月ML≥1.5的地震目录作为研究资料，共2214次地震.地震目录来自于中国台网中心汇编的全国小地震目录库.图1b给出了这些地震的空间分布，其震级随时间的变化见图2.

图1 (a)唐山MS7.8地震震中与余震分布,灰色圆圈为余震分布；(b)选择的研究区域范围与1970年1月—1976年6月期间发生的地震(ML≥1.5)震中分布,灰色圆圈为地震震中.红色六角星为1976 唐山MS7.8地震.震源机制解来自于GCMT目录

图2 ML≥1.5地震震级随时间变化

2 唐山地震前b值变化特征

2.1 分析方法

采用极大似然法计算b值(Aki，1965):

(1)

b值的95%置信度的标准差为

(2)

先将研究区域划分成0.1°×0.1°的格点，以每个格点为中心点，选取2°×2°的空间子区域，计算每个子区内地震的b值随时间的变化曲线.然后，挑选出震前b值呈明显下降变化形态的子区，并以该子区的中心点作为一个异常点，可以把多个异常点集中的区域看成是b值异常区.然后，选取这个异常区内的地震计算b值随时间的变化，这样就可以得到整个b值异常区的范围及其b值随时间变化的特征.

计算每个子区b值随时间的变化时，使用的是等事件数的样本窗口，将该窗口以相等事件数的增量进行滑动.因b值是一个统计量，一般会受样本量的影响，选取等事件数量的窗口，目的是为了确保分析结果不受样本大小差异的影响.

2.2 结果

对于每个子区，以100个地震计算b值，5个地震滑动，计算b值随时间的变化曲线.共计算了2392个子区的b值随时间的变化.有的子区b值个数较少，其时间变化曲线难以反映出可靠的b值变化信息.我们选取具有15个及以上b值数据的子区进行分析，共1747个子区.图3a中给出了上述2392个计算了b值随时间变化的子区中心点位置，其中灰色圆点“”为b值数据小于15个的子区，红色圆点“”为具有15个及以上b值数据的子区.从1747个子区的b值随时间的变化曲线中找出震前b值呈下降变化形态的子区，靠人工是很费事的.我们在b值时间变化曲线上取定1973年1月—1976年6月间的m个b值数据，对于相邻两个b值，计算出后1个b值与前1个b值之差，一共得到m-1个b值差值.再计算出m个b值中最后1个b值与第1个b值之间的相对变化量.挑选b值呈下降变化形态的子区的条件设定为：(1)m-1个b值差值中有55%为负值，即多数b值是下降的；(2)相对下降变化量≥10%.根据这两个条件，可用程序自动挑选出满足设定条件的子区.图3b给出了按设定条件挑选出来的子区的中心点的空间位置，以绿色圆点表示，共391个点.这些绿色圆点中有376个点集中分布在京—津—唐一带地区.我们从376个点中随机抽取了10个点，图4中给出了它们所在的子区的b值随时间的变化，还给出了各子区的中心点的经纬度，竖直线“│”表示唐山MS7.8地震的发生时间，曲线末端的数字为10个子区的序号，这些序号所在的空间位置也用红色数字标注在图3b中.从图4中可以看出，这10个子区的b值在震前确实呈下降变化形态，说明前述设定的选择条件是可行的.选取图3b中黑色粗实线围成的区域，这个区域比376个绿色圆点集中区稍大一点，可以基本包含整个b值下降异常区，以100个地震计算b值，5个地震滑动，可以得到b值随时间的变化曲线(图5).可以发现，在唐山MS7.8地震发生前，这个区域内b值大约在1973年6月从1.1开始下降，到1974年10月达到最低值0.72，相对下降幅度约为35%.从1974年11月到1976年6月，b值维持在低位，没有明显变化.

图3 (a)b值时间变化的子区中心点位置分布.灰色圆“”为b值数据小于15个的子区，红色圆“”为具有15个及以上b值数据的子区；(b)b值时间变化呈下降变化形态的子区中心点位置分布(绿色圆点).六角星为唐山地震震中

图4 随机选取的10个格点为中心点，2°×2°的空间子区域内地震的b值随时间的变化

图5 b值(黑色粗实线)随时间的变化

综上所述，唐山地震前在京—津—唐一带出现大面积b值下降异常，从1973年中到1974年10月，b值呈趋势下降变化形态.之后，直到临近地震发生时，b值稳定在低值，没有明显变化.在空间上异常面积较大，在时间上表现为先下降，然后维持在低值平稳变化，这就是唐山地震前b值的变化特征.

3 结论与讨论

着眼于b值随时间变化过程中的下降变化形态，采用格点搜索法，将研究区域划分成0.1°×0.1°的格点，以每个格点为中心点，选取2°×2°的空间子区域，计算每个子区b值随时间的变化曲线.根据各子区b值随时间变化曲线，得到了唐山MS7.8地震前b值下降异常分布区域.b值下降异常集中分布在京—津—唐一带地区，位于(116°E—119°E，38.8°N—41°N)的经纬度范围内，所占面积大约为35000 km2.

在b值下降异常区内，唐山MS7.8地震发生前b值大约从1973年6月到1974年10月经历了17个月的下降过程，从1.1降至0.72，相对下降幅度约为35%.从1974年11月—1976年6月的20个月中，b值一直维持在低值，没有明显变化.一般把余震分布区作为地震的破裂区，唐山地震的破裂区约占b值异常区内的26%.震前破裂区内发生的ML≥1.5总共111个，显然难以给出破裂区内可靠的b值随时间变化曲线来.

为了考察计算b值的地震个数的影响，分别以150、120和80个地震为一组计算b值，以5个地震滑动得到b值随时间的变化.以同样的挑选b值呈下降变化形态的子区的设定条件进行挑选，得到的b值呈下降变化形态的子区的中心点分布如图6.当计算b值的地震个数增加时，满足挑选条件的子区数会减少，b值呈下降变化形态的子区的中心点也会有所减少；当计算b值的地震个数减少时，满足挑选条件的子区数会增多，b值呈下降变化形态的子区的中心点略有增加.总体上来讲，b值下降异常空间分布概貌没有本质的变化.

图6 b值时间变化呈下降变化形态的子区中心点位置分布(绿色圆点)

b值变化除了与应力有关外，还与介质的均匀度和热梯度有关.应力增加，b值降低(Scholz，1968；Wyss，1973)；介质越均匀，b值越低(Mogi, 1962)；热梯度降低，b值也降低(Warren and Latham, 1970).一般认为强震发生前存在区域应力增强过程，所以认为震前几个月或几年的b值下降变化与应力有关，这更容易理解.2011年日本本州MW9.1地震前，震区b值自2005年起下降，一直持续至临近地震发生前(Nanjo et al.，2012).同期，视应力呈上升趋势变化(李艳娥和陈学忠，2011).即日本本州MW9.1地震前，视应力与b值呈反向变化，即负相关.视应力上升，同时b值下降.视应力与b值之间的负相关关系是构造应力增强的结果.因此，尽管缺乏唐山地震前视应力变化结果，但仍然可以把b值的下降过程理解为构造应力的增强.

选取图3b中黑色粗实线围成的区域内的地震，以40天为时间窗口统计地震的应变释放，然后将统计的应变释放进行累加，可以得到累积应变释放曲线，如图7a.根据累积应变释放数据，以25个数据为窗口，1个数据滑动，通过直线拟合可以得到累积应变释放曲线的斜率K随时间的变化，如图7b.对于应变减速释放，斜率K逐渐减小；对于应变加速释放，斜率K逐渐增大.图7b中显示，1975年4月以前，K逐渐减小，应变释放是减速的；从1975年4月直到临近地震发生时，K逐渐增大，应变释放是加速的.当K逐渐增大时，b值维持在低值，没有明显变化.因此，低b值条件下的应变释放加速可能更具有预测意义.

图7 (a)累积应变释放曲线；(b)累积应变释放曲线斜率K随时间的变化.“↓”为唐山地震