唐山及其邻区地质构造复杂，张家口-蓬莱地震构造带与华北平原地震构造带在这里交汇，多条断裂遍布周围区域（图1）。1976年唐山地区先后发生唐山7.8级地震、滦县7.1级地震和宁河6.9级地震，大地震虽已经过去45年，但在这期间唐山及邻近地区仍为小震频发的态势，近期又先后发生了2020年7月12日古冶5.1级地震、2021年4月16日滦州4.3级地震等中强地震。现今唐山及邻近地区发生的地震是否属于唐山地震余震？这些中小地震受到了唐山历史大地震怎样的影响？仍是各专家学者研究的热点。

图  1  唐山及邻区断裂分布图

前人从多个方面研究了唐山地震与余震的对应关系。其中，万永革等分别用简单破裂模型及基于分段破裂模型结合粘弹性力学模型，研究了唐山主震对滦县地震、宁河地震及3次大震对唐山地区后续地震的库仑破裂应力作用情况，得到了唐山地震对滦县地震及宁河地震有触发作用、3次大震对后续地震有触发作用的结论^[1-2]；刘桂萍等运用3段主震破裂断层模型，研究了唐山主震对3个余震区的弹性触发作用，发现主震有助于后续地震的发生^[3]；蒋长胜等利用Etas模型对2010年以来发生的3次4级以上地震进行分析研究，得到了唐山现今地震为背景地震的可能性较大的结论^[4]；王辉等对唐山地区2012年发生的2次M4.0以上地震进行震源机制讨论，结合对30余年唐山地区小震活动和地壳水平应变率计算，认为唐山强震区的余震数十年后仍在持续，震后粘性松弛可能是华北地区长时间余震持续的主要原因^[5]；朱琳等基于Burgers流变模型，模拟了1976年唐山强震群引起的震后形变场以及同震和震后库仑应力变化，分析认为强震群应力调整过程已基本稳定，古冶5.1级地震是1976年唐山大地震的余震的可能性不大^[6]。

前人通过不同方法理论对唐山地区现今地震是否为历史大震余震进行分析得到的结论不尽相同，静态库仑破裂应力的计算分析是用于判断不同的地震及其余震活动分布的一个重要可靠方法。本文将在万永革等基于粘弹性力学模型计算唐山、滦县及宁河地震在唐山地区的库仑破裂应力基础上^[2]，采用万永革等用形变资料计算1976年唐山地震序列的破裂分布后得出的唐山地区非均匀破裂模型^[7]，对1976年唐山3次大地震发生后的库仑破裂应力的演变过程进行计算，并将库仑破裂应力分布同地震时空特征分布对比分析。同时，收集唐山及邻近地区2002—2020年部分中强地震震源机制解及近期发生的古冶5.1级、滦州4.3级地震的震源机制解，通过计算其受3次大地震产生的库仑破裂应力情况，进一步探讨历史大震对现今地震的触发情况。

根据库仑破裂准则，当断层面剪切应力达到摩擦强度时，断层面将发生剪切破坏。库仑破裂应力变化$\Delta {\sigma _f}$定义为

式中：$\Delta {\tau _s}$为断层面上沿滑动方向的剪切应力；c为介质摩擦系数；$\Delta {\sigma _n}$和$\Delta p$分别为断层面上的正应力和地下的孔隙流体压力，张开为正^[8]。按照万永革等^[2]对唐山地区粘弹性力学模型地震序列应力触发研究中的做法，摩擦系数c取值为0.4。如果$\Delta {\sigma _f} {\text{＞}} 0$，则该应力有利于后续地震发生；否则，不利于后续地震的发生。

当大地震发生后，短期的库仑破裂应力体现在弹性形变上，而随着时间的积累粘弹性特征会越发明显。现有的认知地壳的介质为孔隙流体介质同弹性体介质的混合体。比较符合真实的介质模型为麦克斯威尔体^[9]，即一个弹性体和一个粘性体的串联介质。在短时间内呈现的是弹性体的特征，在长时间尺度呈现的是粘性介质特征。

麦克斯威尔体的总应变可表示为

引入剪切模量和粘性系数的表达式为：

式中：$ \mu_{m} $、$ \eta_{m} $分别为弹性介质的剪切模量和粘性介质的粘性系数，该方程为麦克斯威尔体的本构方程^[10]。

考虑到麦克斯韦粘弹性介质模型的力学特征，平衡方程的傅立叶变换解可以用剪切模量和体变模量给出的解来表示，即

式中：$ \mu (z) $，$ \kappa (z) $，$ \eta (z) $分别为剪切模量、体变模量和黏度系数^[11]。基于此关系对式（1）进行计算，可得出麦克斯韦粘弹性介质模型中地震产生的库仑破裂应力变化。

震源模型：万永革等用形变资料计算1976年唐山地震序列的破裂分布，在计算过程中将唐山地区地震断层分为4个大段，共计108个小断层，各个子断层的长度和宽度大约为5 km×5 km，最终得到了由108个子断层组成的非均匀的破裂分布结果，将该结果作为震源模型来进行计算^[7]。

地壳模型：选取万永革等^[2]基于粘弹性力学模型计算唐山地震序列库仑破裂应力触发的地壳结构模型（表1）。其地壳模型综合分析了刘昌铨等^[12]、曾融生等^[13]、徐锡伟等^[14]地震探测研究结果以及于湘伟等^[15]的反演计算结果，结合张学民等得到的体波速度结构模型^[16]得到唐山地区的综合地壳模型；地壳密度模型采用郑天愉等^[17]运用近震记录确定唐山地震震源过程时所给出的唐山地区地层密度；粘度模型分成了三大部分，其中上地壳以上部分（表1中1～4层）选取了一个较大的粘滞系数1.0×10³⁰ Pa·s，下地壳部分（表1中5～7层）和上地幔及以下（表1第8层）采用孙荀英等^[18] 给出的华北板块下方深部物质的粘性结构，粘滞系数为7.1×10¹⁸ Pa·s和2.1×10¹⁹ Pa·s。

接收断层参数：在计算3次大地震对唐山地区的库仑破裂应力分布过程中需将总库仑破裂应力投影到可能的最优破裂面及滑动方向上。前人对唐山地区现今应力场进行了诸多研究，本文选取了王晓山等^[19]将京津冀地区分成了1°×1°的网格后基于标量断层类型值分类方法，采用MSATSI软件进行反演得到的结果，同时将本文研究区域分成4个区域，选取了其中对应的4个计算结果来进行库仑破裂应力计算，具体范围及接收断层参数结果见表2。

地震数据：自1976年来唐山地区发生的地震数量较大，地震重定位可以更准确地反映出地震发震位置，而双差定位方法是一种目前应用较为广泛的相对定位方法，本文选取了双差重定位后的地震分布结果同库仑破裂应力分布进行对比分析。但鉴于1979年前的地震观测资料收集较困难，故将研究区域内1979年以来双差重定位后的地震数据同库仑破裂应力分布进行对比分析。

计算时间段划分：张素欣等通过逐月对比分析了唐山及邻区地震分布图像、分布频次及强度等特征，通过对唐山老震区的余震时空分布特征的研究，得到了1976—2015年中5个时空分布特征较为一致的时间段^[20]。本文选取了其中自1979年开始的后4个时间段，并增加了2016—2021年4月这一时间段，共计5个时段来进行库仑破裂应力分布情况计算。

通过计算，1976年唐山M_S7.8地震发生后的地震深度平均值为10.64 km，故选取了地下11 km的深度进行库仑破裂应力分布计算。基于库仑破裂应力计算的理论方法及所选取的模型参数，计算得到3次大震在指定点的总应力变化，然后在选取的唐山地区可能的最优破裂面上进行投影，从而得到大地震在唐山地区库仑破裂应力演变过程。

研究区域内地震双差定位结果来自2个部分：1979—2008年定位结果由王晓山提供，该时间段共计5 906个地震重定位结果；2009—2021年重定位结果通过双差重定位程序进行重定位后得到。Waldhauser等在提出双差重定位方法后，进一步开发了双差定位程序^[21]。本文选取了分布在河北、北京、天津等地的168个地震台站经纬度信息，以及经河北遥测台网地震数据库下载到的5 723个地震震相到时资料，通过双差定位程序进行双差重定位后得到4 865个地震定位结果；定位均方根残差为0.21 s，数值较小；在水平和垂直向的平均相对误差分别为0.82 km、0.92 km、1.15 km，重定位前后震中分布及误差分布见图2～3。得到的重定位震中分布结果显示，震中分布沿断裂分布更加收敛。

图  2  重定位前后震中分布图

图  3  水平及垂直方向重定位平均相对误差分布

将研究区域内对应时间段双差重定位后的地震分布同库仑破裂应力分布对比分析得到图4。

由图4可以看出，在区域I、II中地震分布同库仑破裂应力分布相关性更佳，而区域III、IV中地震分布同库仑破裂应力正区域整体分布在东北和西南方向有一相关性，但在细节部分对应度不够，尤其是大震破裂迹线周围的地震分布十分密集但库仑破裂正应力分布却无法与之完全吻合，出现此现象的原因可能是这2个区域中的断裂分布较为复杂。由于断层错动的应力释放作用对地震的发生起到了主导作用，故剔除了3次大震的破裂迹线附近区域发生的地震后，分别对5个时段对比图中正负应力区域地震分布个数进行统计（表3）。

从3次大震产生的库仑应力分布情况演变过程中可以看出，4个小研究区域拼接得到的整个研究区域中库仑破裂应力情况较为连续。总研究区域的库仑破裂应力较清晰地展示出：分布大致呈正负应力交替出现的蝴蝶形特征，在震后10年内粘性因素造成的库仑破裂应力变化较小，在发震10年后粘性因素造成的应力变化趋于明显，尤其是库仑应力为正的区域在逐渐扩大，这一过程与麦克斯威尔体的粘弹性介质松弛过程一致。从地震分布演变过程可以发现，在大地震发生后的10年内，地震主要集中在3次大地震的破裂迹线周围。而此时间段内库仑破裂应力的变化也较为缓慢，且地震的发生可能与断层的应力释放更相关，同张素欣等对唐山地区的时空演化特征中同时间段为正常的衰减过程结论一致^[20]。从1986年开始正应力区域开始向外扩张，其中西北和东南的正应力花瓣变化较为明显，1996年后西北和东南的正应力花瓣中正应力的变化速率更加明显。根据前人对“库仑破裂应力对地震触发作用”的研究，触发阈值为0.01 MPa^[22-24]，在库仑破裂应力变化的过程中处于受触发区域地震发生的个数也明显多于影区数目。在1986—2015年间，随着正应力区域的扩大，触发区域地震发生的个数表现为显著增多特征，这说明历史大震对唐山地区后续地震的发生具有一定的触发作用。值得注意的是，在北西触发区域的地震要多于东南区域，这与万永革等唐山大地震对唐山地区余震触发研究中给出的“大地震有利于唐山地区余震地震发生且在东南和西北花瓣区域余震数目更多”^[2]的结论一致。2016—2021年4月地震分布图结合表3的统计结果可以看出，在研究区域北西和东南的正应力分布区域中仍有较多地震发生，但是发生在影区内的地震比例较1986—2015年有了明显的上升，该现象的出现可能预示着历史大震对唐山地区地震发生的影响在逐渐减弱。

基于较为平均的唐山地区断层接收参数计算库仑破裂应力结果可能不够精细，而使用前人计算出的单个地震的震源机制参数来计算其受库仑破裂应力情况则更加精确。万永革提供了其对研究区域2002—2020年部分地震震源机制解计算结果。通过对所有震源机制解2个节面受库仑破裂应力情况进行计算后，首先筛选出了2个节面受库仑破裂应力作用正负相同的结果；然后对受库仑破裂应力正负不同的两节面，参考万永革等^[25]利用小震分布计算出的对应区域可能的断层面参数（表4），保留与其所得断层面参数较为接近的节面计算结果。经筛选后的结果见表5，受正负应力作用地震分布见图5。

图  5  受正负应力作用的地震分布图

由图5可以看出，收集到的震源机制解结果全部分布在研究区域的III和IV区域，主要集中在唐山断裂、蓟运河断裂周围。分布在III区域的地震，受库仑破裂正应力个数较多；分布在IV区域的地震多发生在唐山断裂的北段，收集到的中强地震受历史大震库仑破裂应力作用情况见表5。

结合图5和表5结果分析，在III区域蓟运河断裂附近，受库仑破裂应力触发的地震较多，而在IV区域特别是唐山断裂北段受正应力作用个数要少于位于影区的地震个数。其中在2002—2009年内正应力区域地震略多于应力影区地震，在2010—2020年间应力影区的地震远大于正应力区域地震。造成该现象的原因可能是IV区域中断裂分布较为复杂，且唐山7.8级地震震中附近因多次大震余震发生导致该区域地下介质较为破碎，不易积累大的应力，同时破碎会发生蠕滑消耗能量，所以断层活动的作用力向NE方向传递，进一步引起了该区域新的活动^[20]，故该区域受到历史大震库仑破裂应力的调制作用较小。

2020年7月12日发生了古冶5.1级地震，2021年4月16日发生了滦州4.3级地震。本文选取万永革等综合各家结果得到的震源机制中心解作为接收断层参数，即：古冶地震（走向：238.96°，倾向：76.76°，滑动角：−172.58°）、滦州地震（走向：295°，倾向：75°，滑动角：−22°），进行其所受库仑破裂应力计算，得到结果为古冶地震受历史大震库仑破裂的应力作用，约为−0.05 MPa（图6），滦州地震受历史大震库仑破裂应力作用约为−2.0 MPa（图7），计算结果显示这2次地震受到了抑制作用。

图  6  古冶地震受大震库仑破裂应力图

图  7  滦州地震受大震库仑破裂应力图

破裂模型的选取是计算库仑破裂应力中的一个较为关键的环节。本文在库仑破裂应力计算中选取的破裂模型将唐山地区的破裂结构分为了5 km×5 km的共计108个子破裂分布，较前人的破裂模型精细了许多。万永革等2008年对唐山大震库仑破裂应力的计算中选取的破裂模型较前人的计算引入了10 km范围内外的修正^[2]，已将结果精确性进行了提升。而本文中选取的模型子破裂结构的范围宽度为5 km，远小于其对断层附近的修正范围，故本文中的计算结果可能更接近真实情况。

基于较为精细的破裂模型计算出的库仑破裂应力分布结果同重定位后的地震分布对比，理论上可以更加清晰地展示出地震分布同库仑破裂应力分布的对应关系。但就所得结果可以看出，地震分布同库仑破裂正应力分布的大方向上有一定的对应关系，而在应力影区的地震数目却仍然较大，尤其是地震在历史大震的破裂迹线周围更加集中。在该地震集中区域内，张素欣等^[20，26]研究发现在1996年之前余震分布较均匀，1996年以来余震在唐山-古冶断裂的东北端点附近较集中，但是库仑破裂应力的变化与这些特征对应关系不够明显。这可能是唐山断裂附近区域断裂分布多且复杂，断裂间相互的吸纳能力和阻隔作用及地壳介质较破碎发生蠕滑对能量积累消耗产生影响^[20]造成的。

通过进一步对2002—2020年部分中强地震受库仑破裂应力作用计算结果分析，得到了历史大震对唐山断裂及蓟运河断裂附近区域的触发情况。由于部分震源机制解结果因2个节面受力的正负情况不一致或与计算所得断层面结果不相近而进行剔除，因此用于计算的震源机制解数量较少，所得结果有误差存在。就该区域库仑破裂应力作用情况，后续需综合更精确的断层接收参数、断裂间的相互作用、地介质能量变化等因素进行计算。

1）地震分布同库仑破裂应力分布对比分析显示，在正负应力花瓣分布区域中正应力区域地震分布的地震较为密集且数量较多，地震分布与正应力花瓣形分布方向有一定的相关性，在北西和东南方向触发作用较为明显。从时间上分析，历史大震发生10年后对唐山及邻近地区地震发生的触发作用更加凸显。近年来发生在库仑破裂应力影区的地震占比显著增大，可能预示着历史大震对唐山地区地震发生的影响在逐渐减弱。

2）通过对2002—2020年部分地震震源机制结果进行受库仑破裂应力作用计算，认为历史大震可能对唐山断裂北段地区地震发生的触发作用较小，对蓟运河断裂附近区域地震发生的触发作用较大。

3）通过对2020年7月12日古冶5.2级地震、2021年4月16日滦州4.3级地震受历史大震库仑破裂应力作用计算分析，认为这2次地震受到了历史大震的抑制作用。

致谢　审稿专家提出了宝贵意见，本文使用了Waldhauser F.和Ellsworth L提供的双差定位程序，Wessel P.和W. H. F. Smith提供的GMT4绘图程序，在此一并致谢。