1976年发生的唐山7.8级大地震使百万人口的唐山市顷刻间变成废墟，死伤近40万余人.之后，地震科技工作者在过去的30多年中，对它进行了全方位的、大量的研究^[1~5]，包括前兆、地质结构、孕震过程等，以便探询地震发生机制.由于地震科学仍处于艰难的探索阶段，人类对地震孕育过程、前兆及前兆场的形成与演化图像至今知之甚少，这些研究结果尚不能全面地解释唐山地震的孕育、发展、发生和各类前兆响应，如唐山地震前一段时间震中周围及其附近地区小震乃至中强地震增强对唐山地震有何作用以及如何作用？唐山地震初始破裂后，为何要分别向北东、西南方向破裂，而不向一个方向破裂？

地震的孕育与发生造成应力在空间上重新分布或转移^[6]，应力重新分布或转移可以分为准静态与动态两种，静应力（static stress）对地震的影响国内外研究比较早，取得了大量成果^{[5, 7~10]}；但是，由于静应力变化量具有近程性^[7]、时滞性^[8]、无方向性^[10]等，限制了对地震机理的揭示和对地震现象的解释.国外一些学者于20世纪90年代初开始从动应力（dynamic stress）角度研究地震触发问题^{[7, 11, 12]}.而动应力为何触发地震，即机制方面人们对它还很不清楚，研究结果仍是推测：比如液化过程中饱和断层接触面孔隙压力的增加^[13]；周期性疲劳至损^[14]；调整地热系统扩散过程^[15].国内动应力相关研究工作刚开始不久，未曾全面、系统地展开，现在有很多学者关注此方面研究情况并开始作了一些相关研究工作^[16~18]；但还没有运用动应力作用及触发机制对一个地震进行完整研究的报道，本文运用三维非线性动态有限元方法（Non-linear Dynamic FEM）仿真模拟研究了相似唐山地震孕震体的菱形模型的动应力作用，结合有关地学理论、相关研究成果和地震资料，对唐山地震进行动态分析和研究，从动态角度探讨唐山地震的发震机制和过程.

本文在数值模拟计算中采用拉格朗日坐标构架.动态有限元原理如下：对于某一时刻t具有总体积V，包含有m种介质材料，内作用面S_in，位移边界S_u与外力边界S_T的相互作用系统，应满足虚功原理^[19]并可得到如下关系式：

(1)

式中σ_ij是应力；V_m是第m种介质体积；S_Tm是第m种介质边界面积，dV是体积微分元，dS是微面积元；ρ_i对应于第i种介质密度；χ_i（i=1，2，3）是位置坐标；是加速度；δχ_i和δχ_i，j分别是虚位移及其增量；F_i和T_i是体积力和面积力.

将总体积离散化为n个单元，为了保持四面体单元的抗畸变能力，同时又使计算精度得到保证，在4节点四面体单元的基础上，建立8节点六面体超级单元.每个单元的节点数为8，位置坐标χ_i（i=1，2，3）可用形函数N_j（j=1，2，…，8）与节点位置χ_i^j的插值表示为

(2)

其中，ξ，η，ζ是单元的局部坐标，t是时间变量，形函数N_j定义为

(3)

式中，ξ_j，η_j，ζ_j是单元局部坐标系下节点位置，取值限为（±1，±1，±1）；于是离散后的虚功原理可用矩阵表达为

(4)

式中，N是单元形函数N_i的集合矩阵；σ是应力矢量，F与T分别是给定单元的体积力与面积力载荷矢量；a是节点加速度矢量，与加速度的关系为；B是变形协调矩阵，定义为B=D^TN，D为矩阵算子，定义为

(5)

上面（4）式即为进行动态计算的基本公式.但是，由于动态系统中一般具有较强的非线性（几何、材料和内作用边界的耦合非线性），要求动态计算中时间积分增量适当小，从而能正确清楚地反映动态响应过程.因此，本计算中采用了显示积分方法，根据计算稳定性的要求，在计算过程中自动确定适当的时间增量步长.为了讨论简便起见，取以下形式的动力学平衡方程：

(6)

式中M、K、u、F、t分别为质量矩阵、刚度矩阵、位移、载荷、时间.两阶中心差分的格式为

(7)

(8)

由以上三式，可以得到每一时刻递推关系：

(9)

图 1为唐山7.8级地震发震构造图，由图可知：唐山地震孕震体是被四个断裂，即宁河-昌黎断裂、丰台-野鸡坨断裂、滦县-乐亭断裂和蓟运河断裂所包围的北东东向的近似菱形块体^[1]，初始破裂后向东北传播70 km、向西南传播45 km，破裂深度34 km^[3]；同时实际地震初始破裂处不可能是点或线.为此我们建立如图 2所示的模型：长120 km、宽70 km、厚35 km的菱形模型，中间有8 km的孔，模拟地震初始产生的冲击性加卸载力的作用.考虑到对称性，只取一半定义对称面.考虑三维砖单元，由于载荷和形状的轴对称性，属于面应变问题，因此Z方向的尺寸对计算结果无影响.

图 1

Fig. 1

图 1 唐山7.8级地震发震构造图[1] 1.与发震构造有关的断裂; 2.发震断裂; 3.推测7.8级地震发生时的破裂发展方向; 4.地裂缝剪切带; 5.平移断层; .地震时地面上升区; 7.地震时地面下降区; 8.地震震中; 9.7.8级地震烈度线; 10.区域应力场中主压应力方向; ①丰台-野鸡坨断裂; ②蓟运河断裂; ③宁河-昌黎断裂; ④滦县-乐亭断裂. Fig. 1 Seismotectonic map of Tangshan M7.8 earthquake[1] 1.Fault related with Tangshan event; 2.Seismic fault; 3.Inferred rupture direction of the main shock; 4.Ground fissure; 5.Translational fault; 6.Rise area during the event; 7.Subsidence area during the event 8.Epicenter; 9.Seismic intensity line of M7.8 earthquake; 10.Direction of regional principal stress; ①Fengtai-Yejituo fault; ②JiYunhe fault; ③Ninghe-changli fault; ④Luanxian-Leting Fault.

图 2

Fig. 2

图 2 菱形模型 Fig. 2 Rhombic model

模型采用介质均匀弹塑性材料，密度2700 kg/m³，单元定义为三维弹塑性可破坏单元，弹塑性用双线性强化模型表征如图 3a，图中σ、ε、σ_a和σ_b分别为应力、应变、屈服应力和断裂应力.体变模量为1.5× 10¹¹ Pa，强化模量为1.0×10⁸ Pa，屈服应力为-2.0× 10⁹ Pa（压应力），过了屈服点，材料进入非线性状态，断裂条件为拉应力达到0.8×10⁹ Pa，此应力小于屈服应力，因此在拉应力状态下尚未进入非线性状态单元就已经发生断裂了，即压应力破坏条件定得很高，此处不会发生压引力破坏，这与地壳介质长期受高围压作用而不易破裂的情况一致.引入两个假定：（1）单元断裂后，裂纹不会延伸到邻近单元；（2）单元断裂后，不能承受也不传播任何应力.因此，在算法上使断裂的单元消失，而该单元的质量也同时被忽略.假定研究对象4个侧面为自由边界，冲击波在自由面上的反射会自然地发生“半波损失”，也就是说压缩波反射后成为拉伸波；与Z轴垂直的两个表面加以约束：限制Z方向的位移，即dz=0.由于载荷是轴对称的，问题本质上是平面应变问题，但仍采用三维有限元定义对称面、建模计算，即dy=0；从而整个模型在三维空间中不会发生刚性位移、不会有刚性转动.

图 3

Fig. 3

图 3 双线性强化模型(a)和载荷曲线(b) Fig. 3 Bilinear strengthen model (a) and loading curve (b)

在现实中，水平构造应力和垂直重力只引起一个静态的不变的应力场.因此，不加载水平构造应力和重力，定义密度，各点的惯性力都会根据有关节点的加速度自动被纳入计算.在中央圆孔上加载一个冲击压力，载荷因子为4.5×10⁹ Pa，载荷因子与载荷曲线相乘，得到加载的载荷-时间曲线.图 3b为载荷-时间曲线图，该曲线是归一化的，载荷在0时刻从0突变到1，然后持续0.07 s，所以这是一个冲击，相当于地震发生时在这个柱形孔中产生的冲击力.设定输出时间间隔，计算结果处理后，可以得到任意物理场任意视角的动画画面.

计算仿真结果如图 4所示.图 4a是压缩性冲击力产生的压缩波从加载面向外传播的某一瞬间，一部分刚刚触及边界.由于边界为自由面，压缩波达到边界后，发生半波损失，压缩波变成拉伸波.图 4b是压缩波反射的拉伸波从边界向内传播，在某些地方相遇，叠加或干涉使拉伸力增大，在短时间内达到破裂阀值，结果在某些地方出现了断裂，最终出现大量断裂的情况.此后，由于已有的断裂形成了新的边界，波的传播变得特别复杂.结果发现整个动态断裂过程、位置变得十分复杂而不规则，说明反射的拉伸波发生叠加、或干涉的位置与模型边界形状有关，不同形状的边造成反射的拉伸波发生叠加、或干涉的位置是不同的，即动态断裂与模型边界形状有关.

图 4

Fig. 4

图 4 模拟结果(σ单位:Pa) Fig. 4 Simulation results (unit of σ:Pa)

根据岩石力学的Coulomb-Mohr破裂准则可以对动应力作用模拟结果进行解释.根据岩石力学的Coulomb-Mohr破裂准则，当破裂面的剪切应力τ到达该物体抗剪强度τ_f时，破裂面将发生破裂而产生地震.库仑破裂应力^[20]CFS=τ-τ_f，抗剪强度τ_f=τ₀ +μσ，其中τ₀为接触面固有剪切强度，μ是内摩擦系数，σ为围压.从上面的破裂准则可看出：拉应力增大，围压σ降低，致使抗剪强度τ_f也降低，这样在剪应力τ不变的情况下，等于增大了库仑破裂应力CFS，从而发生破裂，形成地震.因此，要增大拉应力，就必须增大压缩性冲击力，但在静态情况下，破裂以及断层错动必须要增加足够大的剪切应力，否则不能发生.因为它不仅要克服静态增大的压力所产生的摩擦力，而且还要克服以前围压产生的摩擦力；而在瞬态增大压缩性冲击力的情况下，破裂以及断层错动结果大不一样，而且是压缩性冲击力增大得越大，在不改变剪切力情况下，破裂以及断层错动越容易.

地应力是最直接反映孕震体状态的物理量，图 5为唐山地震地应力观测曲线.图 5a是唐山地震震中区赵各庄站地应力1976年7月27日一段时间速率曲线^[1]，发现在震前出现了数小时明显而又持续的张应力活动^[1]；图 5b是离震源区较远的大连站地应力曲线^[1]，发现在震时及以后一段时间，应力均出现低值和波动；图 5c是震源区附近一些站地应力变化曲线^[1]；图 5d是震源区及附近简易地应力变化曲线^[1]，也具有这种特征.文献[1]对地应力研究后认为震中区震时及震后受力状态均以张应力为主，震后变化几乎都是在地震当天完成的.而大量破裂出现时的模拟结果（图 4b）所表现出来的应力变化与图 4a压缩性冲击力产生的压缩波从加载面向外传播的某一瞬间完全不同，应力变化十分复杂且在断裂处应力有所降低，这与实际地震时应力观测结果比较一致.其原因是在地震发生前后岩石产生裂纹、破裂并发生串通，其过程产生压缩性冲击力，压缩性冲击力所产生的冲击波在遇断层侧面后反射形成的动应力叠加到围压上，使围压出现较大的波动和下降.产生的压缩性冲击力加卸载存在时间差，所形成的加卸载冲击波相位不同，而应力波又在破裂面和周围断层之间来回反射，结果形成的动应力和围压叠动后致使围压剧烈波动^[17].唐山地震的发震机制和过程存在一种可能与静态模式不同的动态机理.

图 5

Fig. 5

图 5 地应力观测曲线图 (a)1976年7月27日唐山赵各庄站地应力速率曲线; (b)大连站地应力曲线; (c)北京昌平站地应力曲线; (d)简易地应力曲线. Fig. 5 Observational data curves of ground stress (a) Rate of ground stress recorded by Zhaoge village station in Tangshan on July 27, 1976;(b) Ground stress recorded by Dalian station by means of inductance; (c) Ground stress recorded by Changping station in Beijing by means of inductance; (d) Ground stress recorded by smple methods by means of electric current.

下面运用上述模拟结果，结合有关地学理论、相关研究成果和地震资料，对唐山地震进行动态分析，探讨唐山地震的动态发震机制和过程.

唐山地震发生前，唐山附近中强地震活动增强，并向未来大震震中-唐山附近集中^[1].图 6为1950年以来燕山带东段4级以上地震累加频度图，可以看出唐山附近地区4级以上中强地震在地震前一段时间出现明显增加的现象.强震前中强地震增强表征了一个正反馈过程，即中强地震产生的瞬态性冲击压力形成的应力波反射后发生半波损失，成为动态的拉伸应力，通过减小内摩擦致使岩石产生裂纹、局部破裂并发生串通，表现为某种链式反应，直到触发强震.从动态应力传播、演化特征看，动态应力容易在断裂带交汇点处周围聚集并呈现高量值^[16]，1976年唐山7.8级地震位于华北板块的北边界阴山东西向构造带南缘-燕山褶断带和冀中块体东边界带-沧东断裂的交汇部位，孕震体是一个被四个断裂包围着的菱形块体^[1]，所以唐山地震孕震体所处位置适合动态应力传播能量的聚集.因此，唐山地震前唐山附近中强地震活动增强所产生的动态应力不仅直接触发地震，而且也为后来孕育唐山大地震积累了能量.

图 6

Fig. 6

图 6 燕山带东段地震(M≥4.0)累加频度曲线图[1] Fig. 6 Accumulative frequency curve of earthquakes (M≥4.0) in the eastern segment of Yanshan fault[1]

1976年唐山7.8级地震位于华北板块的北边界阴山东西向构造带南缘-燕山褶断带和冀中块体东边界带-沧东断裂的交汇部位^[1].用人工地震测深和转换波探测工作的结果都指出，在唐山地区确有切断地壳下界面的深断裂，在震区构造显著的特点是由宁河-昌黎深断裂、丰台-野鸡坨大断裂、滦县-乐亭断裂和蓟运河断裂的切割形成一个北东东向的菱形块体^[1].唐山地震就发生在这个菱形块中，P波资料测定震源参数结果表明^[2]，唐山地震不是简单的破裂而是不对称的双侧破裂，初始破裂在菱形中央偏北、走向是北东30°，但总体走向是北东49°，向北东传播70 km、向西南传播45 km，平均位错136 cm，最大位错173 cm，在地表造成长8 km、宽30 m的地表裂缝带（图 1）.分析其破裂原因是初次破裂产生的冲击波传播出去以后，遇孕震体四周的菱形块体的四个侧面，即宁河-昌黎深断裂、丰台-野鸡坨大断裂、滦县-乐亭断裂和蓟运河断裂，然后反射回来成为拉伸波，在菱形中部相遇，造成相遇处围压降低、摩擦力减小，结果是在强大的剪切应力作用下，发生破裂并分别向北东、西南传播（图 1）.

1976年唐山M_s7.8级地震不是一个简单的破裂过程，而是由北东向30°的8 km长的初始破裂以及分别向两侧传播的不对称双侧破裂，即向北东传播70 km、向西南传播45 km所组成，所以此次主震由三次破裂组成（图 1）.为何会如此发生呢？

从图 1可以看出，唐山地震孕震体是被菱形断层所包围，当初次破裂发生后，产生的瞬间冲击压力波便迅速向外传播，当遇周围断层后，应力波被反射，然后相遇，从图 1、4中可以看出，它们反射后相遇处与后来发生的实际断裂位置基本一致.而造成的破裂又产生新的瞬间冲击压力波，然后又被反射，因此一个应力波在破裂面和周围断层间可以来回反射，直到衰竭，所以破裂可以多次发生.只有在反射的拉伸波所形成的动应力叠加围压后所引起的摩擦强度大于剪切应力，破裂就终止发生.所以，主震的多次破裂实际上是瞬间冲击压力多次遇孕震体断层反射后所形成的动应力致裂后错动而形成的.

唐山地震在初始破裂的基础上分别向北东破裂70 km、向西南破裂45 km^[3]，其速度分别为2.4 km/s、3.0 km/s^[4]，破裂时间分别为29.17 s、15 s；而一次瞬态动应力有效作用时间，即动应力使压力降低的时间为3.755 s、错动滑动时间为6.555 s^[17]，两者比较接近，其差异可能是地壳为多种介质构成的各向异性黏弹性体以及地壳存在各种不同类型和产状的断层所致，所以，从破裂时间方面说明唐山7.8级地震的发生与动态应力作用有关.

唐山地震是一个发生在城市里的特大地震，震时观察到了很多现象，震后地震部门进行了考察和总结，其中之一就是震源区地面剧烈的颠动、摇动等振动现象.据文献[4]记载，在唐山地震发震时沿发震构造裂缝带及其附近出现明显上下颠动现象，这种运动共三次，造成独立砖烟囱从下部折断倒落；在断裂裂缝带附近的地面位移现象也很明显，在裂缝带东侧50 m和西侧200 m，两位目击者发现大地先向南北来回反复运动三次；在裂缝带的西盘及其附近地区，人和物均与地面作相对运动而向西、向南位移2~3 m.据当地群众反映^[4]，包括上、下颠动和水平位移摇动的全过程，约在10 s以内.这些现象说明地震时震区发生了明显的振动，这种振动正是动应力作用的结果：第一可以从动态应力作用时间（有效作用时间3.755 s、错动持续时间6.555 s）^[17]与震时观察到的振动时间（10 s以内）基本一致得到印证；第二可以从动态应力作用使围压呈现剧烈波动，致使发生振动式错动^[17]得到解释.在实际地震发生时将造成破裂、错动时颠动、水平摇动等振动现象.这时由于：破裂产生的冲击力加卸载时间不同、应力波在破裂面和孕震体周围断层之间的来回反射，所形成的动应力和围压叠动后致使围压剧烈波动.动态应力作用造成围压波动而引起地震破裂时颠动、水平摇动等振动又将有利于孕震体裂纹串通、介质松动、接触面变平，进而减小摩擦系数、降低固有剪切强度，从而有利于地震的发生.

唐山地震孕震体四周被四条断裂包围，在其外围还存在大量规模不一的断裂，如：滦县-卢龙断裂、汉沽断裂等^[1].唐山地震发生过程中产生的压缩性冲击力所产生的应力波一部分被包围的四个断裂所反射，一部分穿过断裂继续向前传播并在传播过程中衰减.当向前传播的压缩性冲击力所产生的应力波再遇其他断裂时，又发生半波损失，成为拉伸波，产生拉伸力，从而使相遇处围压减小，增大了库仑破裂应力CFS，从而发生破裂，形成地震即余震.从唐山地震主震破裂看，产生的压缩性冲击力方向是北东、西南方向，所以产生的余震也主要是分布在北偏东的一个近似长方形范围内（图 7）；而动态应力在传播过程中衰减^[16]，所以远离主破裂区，产生的余震应该小，这与实际余震发生结果比较一致（图 7）.

图 7

Fig. 7

图 7 唐山7.8级地震和强余震分布图(1976-07-28~1976-12-23)[1] Fig. 7 Distribution map of Tangshan M7.8 earthquake and its strong aftershocks (1976-07-28~1976-12-23)[1]

当唐山地震产生的压缩性冲击力卸载时，相当于加了一个拉伸性冲击力，它产生的冲击波遇断层发生半波损失成为压缩波，所产生的压缩力使围压增大，结果唐山地震破裂和错动很快停止.因此，在唐山地震产生的压缩性冲击力的卸载作用和动态应力传播能量减小的共同影响下，唐山地震产生的余震越靠近主震，发生的越多，以后逐渐减少，这是一个负反馈过程，其结果与唐山地震实际发生的余震情况相吻合.

通过本文分析与研究发现，动态应力在唐山地震孕育、发展、发生及前兆产生中起了比较重要的作用：

（1） 唐山地震前中强地震增强表征了一个正反馈过程，即中强地震产生的瞬态性冲击压力形成的应力波反射后发生半波损失，成为动态的拉伸应力，通过减小内摩擦而引起岩石产生裂纹、局部破裂并发生串通，直到触发强震；而唐山地震孕震体所处位置适合动态应力传播能量的聚集.因此，唐山地震前中强地震活动增强所产生的动态应力不仅直接触发地震，而且也为后来孕育唐山大地震积累了能量.

（2） 地震产生的压缩性冲击力加卸载存在时间差，所形成的加卸载冲击波相位不同，结果形成的动应力和围压叠动后致使围压剧烈波动和下降，导致地应力在地震前后一段时间均出现低值和波动.

（3） 唐山地震初始破裂发生后，产生的动态应力使围压低值波动，导致摩擦强度降低，使初始破裂分别向北东、西南方向传播.

（4） 唐山地震产生的冲击力卸载，相当于产生拉伸性的应力波，在遇断层时发生半波损失，成为动态的压缩应力，使围压增大并导致断层破裂和错动停止，而动态应力在传播过程中衰减.因此，余震的发生是一个负反馈过程.

虽然对唐山地震进行了动态仿真模拟研究，探讨了唐山地震的发震机制和过程，但由于地壳为非均匀性结构，地震孕育、发展、发生是一个极其复杂的过程；而研究模型为简单的地震孕震环境且与实际构造存在差距，其研究是初步的、探索性的，需要在模型设计上更接近实际以便进行进一步的深入细致的研究.然而，动态应力作用及触发机制是一个崭新的、起步较晚的且研究不成熟的领域，特别是对于国内地学界来说，更是一个陌生的研究领域，对于极其复杂的地震孕育、发展、发生过程，运用各方面知识从不同角度进行不断的、艰难的探索，有助于对地震机理的深入认识和合理解释一些地震现象.

感谢本文匿名审稿者的有益意见和建议.