众所周知，大震后通常伴随了大量的余震活动.而对余震的发生频度和大小的描述则来自于Omori定律和Bath定律^{[1， 2]}.Omori定律和Bath定律都是由观测统计总结出的经验关系.前者描述了主震后余震个数随时间推移而逐步衰减的变化过程，最终达到区域内背景场地震活动性的程度，即从余震活动性转变到背景场地震活动性.目前，由于缺乏统一的判别准则，仍然无法准确地给出余震转至背景场的时间尺度.而Bath定律则指出了主震后余震震级一般比主震震级小1~1.2个震级.目前的研究成果表明，余震的诱发(触发机制)分为应力变化的静态触发和动态触发.前者可由主震同震位移造成主断层以及周边应力场变化.库仑应力转移模型通常用于对余震分布空间的预测.而动态触发机制则来源于地震波能量辐射传播所造成的主断层周边断裂的失稳.动态应力波触发可导致远距离的地震诱发，如Gomberg等^[3]认为断层前缘受破裂方向性的影响，可触发余震.

近期Parsons和Velasco^[4]通过对相似震级爆炸源的分析和模拟，认为近断层的余震触发过程主要来源于主震断层内部形成的应力变化(转移).如果这一结论正确，近断层，尤其是主震断层内部余震的产生则与同震位移所产生的应力变化或静态应力降相关.很明显，目前的研究大都集中在余震的触发机制，对于余震的大小则没有基于物理过程的明确解释.虽然Bath定律^[2]指出了最大余震震级一般比主震震级小1~1.2个震级，但同样缺乏相应的物理基础.

陈学忠等^[5]则通过对中国大陆近期主要地震的视应力(即地震波辐射能量同地震矩之比)与余震强度的统计关系，给出了一次主震后高视应力的地震一般对应了后续较大余震的产生，并尝试为余震危险性评估提供帮助.但是，由于作者忽略了造成视应力大小的震源力学过程或断层破裂过程中动摩擦机制的影响，所得结论存在很大的不确定性，甚至结论同观测有相违之处.2001年M_S8.1昆仑山口西地震和2008年M_S8.0 汶川地震发生于同一构造单元，即巴颜喀拉地块.昆仑山口西地震破裂长度约为500km, 汶川地震的破裂长度超过300 km.由USGS给出的震源参数解表明，昆仑山口西视应力约为5.83MPa, 而汶川地震视应力仅为0.55MPa.主震后两个地震造成的余震大小和时空分布有着明显的区别.中国地震台网中心数据显示，对昆仑山口西地震而言，在主震后两年内发生大于M_S5.0的余震仅为6个，余震的最大震级仅为M_S5.7，余震一般分布在主震断层的两侧；而对于汶川地区，两年内发生了46个大于M_S5.0的余震.至今为止，最后一次M_S5.0的余震发生在2010 年5 月25 日，所有余震的最大震级为M_S6.4，较大震级的余震一般都分布在主震断层带内.由此，我们推测，断层内部的破裂过程，或同震过程中的应变能的变化可能同余震强度相关联.

在这篇文章中，我们采用2008年汶川地震震源物理参数、井下应力测量数据和近断层参数，从简单的断层破裂模型出发，应用能量分配原理和库仑摩擦准则，估算和推断应力变化(静态应力降的大小)对断层内部余震强度的影响.同时，在平均意义下，给出了主震震源的辐射效率和地震效率的估计，并同以往的一些结果进行对比.最后，我们给出了对本工作中存在的一些关键问题的讨论.

地震时，断层周围的介质通过断层的错动释放出所贮存的能量，即应变能

其中τ₁是初始应力，τ₂是最终应力，Δu为最终平均滑动位移，$A = \int_\Sigma {dS} $为断层面积.这些应变能一部分为应变能变化量E_eff, 另一部分通过摩擦以热(Eheat = (-ΔUe)-E_eff)的形式释放^[6].在滑动弱化模型中，应变能变化一部分用于断层尖端的继续延展，称为破裂能E_f；一部分用于辐射地震波，即地震波辐射能E_S.

Savage和Wood^[7]提出了在不考虑破裂能的情况下，将震源动态摩擦过程的描述分为三种主要的模型:完全应力降^{[8， 9]}、应力上调^[7]和应力下调^{[8， 10~12]}.这三种模型分别对应于摩擦应力等于、大于、小于断层面上的终止剪切应力的情况.动摩擦应力上调现象表征了错动过头的物理现象，这一点可由早期Madariaga^[13]圆盘断层动态破裂模型得到描述.同位错模型相比，半径一定时，动态破裂模型给出的平均滑动位移是其1.52倍.静态应力降则为位错模型的5.5倍.动摩擦应力上调给出静态应力降Δσ_s 与动态应力降Δσ_d 的比值为

其中，ρ₀ =1.52，C(v_R/β)为Kostrov函数，其中，v_R是破裂速度，β 是剪切波速度，在v_R/β = 0.9 时，C(v_R/β)等于0.82，Δσs/Δσd 约为1.27，即

在考虑动摩擦应力上调的情况下，断层附近的介质通过断层错动释放出所贮存的应变能，除了转化成摩擦热能、破裂能和地震波辐射能^[14]之外，还有一部分应变能变化分配给滑动停止后断层面上应力松弛所做的功，简称松弛功E_r^{[13， 15]}(如图 1).松弛功不仅体现了主震后的断层上应力转移过程，还表征了周围应力场的变化.Favreau和Archuleta^[15]对美国M_W6.5ImperialValley地震动力学模拟的结果表明，该地震能量分配比例E_r∶E_S∶E_f 约为1∶1.3∶3.5，说明松弛功和破裂能在整个应变能变化中占相当大的份额.事实上，松弛应力过程可能对应了黏弹性体能量逐步释放的过程，其释放率可能伴随了断层内部或周边的余震触发，并进一步影响到余震的持续时间.这些问题涉及较为复杂的断层动力学过程，超过了这篇文章的讨论范围，所以暂不予考虑.在这里，我们则主要以现有数据，计算出松弛功在所有能量中的份额，由此给出后续余震发生的震级上限.

图 1

Fig. 1

图 1 滑动弱化下的能量分配.τ1、τ2是初始剪应力和最终剪应力，τf为滑动摩擦应力，τr为裂纹尖端扩展时的阻碍力，τp为屈服应力或最大静摩擦力，Dc 为滑移弱化距离.地震波辐射能由ES =Eeff-Ef-Er 计算得到.(a)由弹性地壳形变产生的弹性功Eeff；(b)破裂尖端的破裂能Ef 和破裂阻止段的松弛功Er Fig. 1 Energy partitions based on the slip-weakening law.τ1 and τ2 are the initial and final shear stresses applied on thefault, respectively, τf and τr are the dynamical frictional and resisting stresses, respectively.τp is the yield stress whichindicates the strength offault, Dc is the slip-weakening distance.Radiated seismic energy is given by Es =Eeff -Ef -Er.(a)Ee, the elastic work done due to elastic crust deformation； (b) Ef, the fracture work spent at the rupture frontand Er, the relaxation work spent at the rupture arrest

主震过后，通常伴随余震的发生.余震的触发分为两类，一类是近断层的静态触发，是由于主震造成了近断层应力变化，一部分能量以松弛功的方式贮存在周围介质中，激活了近断层产生余震；另一类是主断层外的动态触发，是由主震产生的地震波带来的高频能量激活了断层.对于摩擦应力上调的汶川地震来说，余震非常频繁，而且大多分布在主断层附近，很可能是由于静态应力变化触发了余震^[4].

根据中国国家地震台网测定，2008 年5 月12日在我国四川省汶川县发生了M_S8.0地震，震中位于映秀镇(31.0°N，103.4°E)，震源深度约10km.而在5月6日，中国地震局地壳应力研究所等单位恰好在震区北端青川断裂带内部和断裂下盘的四川盆地里完成了四个测孔的水压致裂绝对应力测量，测孔深度400余米.震前，应力从四川盆地向龙门山断裂带逐渐增大，龙门山断裂带内观测点在408~418m深处附近的最大水平主压应力σ₁为21.8~21.9MPa, 最小水平主应力σ₃约为10.8 MPa.震后，发震断裂附近的观测点390~420m 深度段的主压应力σ₁为15.7~16.0 MPa^[16]，测孔位置如图 2 所示.作用于断层面上的最大主应力σ₁为水平方向，而最小主应力σ₃为垂直方向，汶川地震主断层断裂倾角θ 近地表较大，约为60°~70°，甚至直立，而随深度向下逐渐变缓^[17].由于汶川地震主断层是一次以逆冲为主、兼具小量右旋走滑分量的断层，所以取主断层上平均倾角θ 为33°^{[18， 19]}.因此作用于倾角33°的逆断层的剪切应力τ和正应力σ_N 可由公式τ=1/2(σ₁-σ₃)siN₂θ和σN =1/2(σ₁+σ₃)-1/2(σ₁-σ₃)cos2θ获得.其中σ₁和σ₃分别为最大水平应力和垂直应力，此时σ₃一般等于ρgh，其中ρ 为地壳介质密度，g为重力加速度，h为垂向深度.由库仑破裂准则τ ≥τ_p =μ_sσ_N，其中μ_s 为静摩擦系数，由此可得震前、震后剪应力分别约为τ₁=5.0 MPa和τ₂=2.3 MPa.而正应力则为σ_N=14.06 MPa.库仑破裂准则的临界条件τ =τ_p =μ_sσ_N 时，有

将σ₁和σ₃的比值定义为R，有

因此要使断层滑动R的值最小就要使dR/dθ =0，得到R 的最小值${R^\prime } = {\left[ {\sqrt {1 + \mu _s^2} - {\mu _s}} \right]^{ - 2}}$，此时σ₁与断层走向所夹的角θ=1/2tan^-1(1/μ^s).对于汶川地震，取θ=33°，摩擦系数μ^s=0.5，最小值${R^\prime } = {\left[ {\sqrt {1 + \mu _s^2} - {\mu _s}} \right]^{ - 2}}$= 2.6.断层面上任一点破裂前缘的屈服应力约为τ_p ≈μ_sσ_N =7.0MPa.假设主应力方向与断层走向近似垂直(如图 3)即α=90°，计算得到静态应力降Δτ_s 为2.7 MPa.

如果断层破裂到地表，对于正断层或逆断层，利用Starr^[20]的泊松体模型，静态应力降Δτ_s =8μD/(3πW)，其中W为断层宽度，约为22km^[18]，D为断层上的滑动位移，约为2.23m；采用王卫民等^[21]对汶川地震反演得到的平均位错量，取剪切模量μ=3.0×10⁴ MPa, 则可得静态应力降Δτ_s 约为2.6MPa.

根据一维瞬态破裂模型^[8]给出近断层的初始质点速度和加速度与动态应力降的关系可以得到:

(1)

(2)

其中，最大质点速度和加速度利用兰景岩等^[22]给出的绵竹清平强震记录，如图 2.绵竹清平台站位于绵竹市清平乡，距断层约5km.该台址地层较为复杂，主要由松散土以第四系全新统杂填土、第四系残坡积含碎石粉质粘土等构成，覆盖层厚度9.1m, 属中硬场地土^[22].测得该场地南北向和东西向峰值速度分别为52cm/s和76cm/s, 剪切模量β 取3.0×10⁵cm/s, C为覆盖层放大系数，取C=1.4.考虑到自由表面加倍效应，利用公式(1)求得动态应力降Δτ_d约为1.9MPa和2.7MPa, 平均值约为2.3MPa.而场地南北向和东西向峰值加速度ü 分别为802cm/s₂和824cm/s₂，介质密度ρ为2.7g/cm³，由于加速度波形主频大约在5 Hz左右，所以Δt取0.2s.当破裂速度v$\gg $β时，公式(2)近似为Δτ_{d_a} =ρβΔtü/(2C)，为动态应力降Δτ_d 的上限.考虑到自由表面加倍效应，得到动态应力降上限分别约为2.3 MPa和2.4MPa.综上所述，断层内部平均动态应力降Δτ_d 约为2.3 MPa.由此可见，不同的方法所得到的应力降相差不大，清平台站附近断层动态破裂符合应力上调动态应力降小于静态应力降的基本特点.

图 2

Fig. 2

图 2 地应力测孔位置和清平台站观测记录.黑色箭头所指为测孔位置 Fig. 2 The location of borehole strain and the observations in Qingping station.The solid arrow indicates the location of borehole strain

图 3

Fig. 3

图 3 汶川地震断层空间展布示意图.σ1在与xoy面平行的面内，σ1为最大主压应力，σ3为最小主压应力.θ为断层的倾角，α 为σ1与y轴的夹角 Fig. 3 Spatial distribution of earthquake fault.σ1is the maximum principal stress in a parallel plane of xoy, σ3 is theminimum principal stress.θ is the dip angle of thefault, α is the angle between σ1 and y

利用汶川地震作为研究实例，根据Brune^[8]，有效动态应力降Δτ_d =τ₁-τ_f, 其中τ₁是400m 深处初始剪应力，τ_f则为滑动摩擦应力.由此可得τ_f=2.7 MPa.表 1给出了汶川地震400m 深处震前、震后的应力τ₁和τ₂、屈服应力τ_p、摩擦应力τ_f以及静态、动态应力降Δτ_s 和Δτ_d 的值.

表 1(Table 1)

表 1 汶川震前、震后的应力，屈服应力，摩擦应力以及静态、动态应力降(MPa) Table 1 The initial and final stress, yield stress, frictional stress, static and dynamic stress drop of Wenchuan earthquake (MPa)

表 1 汶川震前、震后的应力，屈服应力，摩擦应力以及静态、动态应力降(MPa) Table 1 The initial and final stress, yield stress, frictional stress, static and dynamic stress drop of Wenchuan earthquake (MPa)

事实上，断层深部τ₁、τ₂和τ_f 都会比400 m 深处的值大的多，变化可能更为复杂.为使问题简化，我们假定Δτ_s、Δτ_d、τ_p-τ₁在不同深度均为常数.这样现有的结果(表 1)则可用于计算远场视应力.由E_r/E_S = (Δτ_s-Δτ_d)/τ_a (如图 1)可得松弛功E_r, 由E_f =E_eff-E_S-E_r 可得破裂能E_f 和平均滑动弱化距离D_c.根据地震波辐射效率η_S和地震效率η 定义，η_S=E_S/E_eff 和ητ=σ_a, 其中τ=1/2(τ₁+τ₂)，可以得到辐射效率η_S和地震效率η.具体参数如表 2所示.

表 2(Table 2)

表 2 汶川地震能量分配 Table 2 Source parameters of the 2008 Ms8.0 Wenchuan earthquake

表 2 汶川地震能量分配 Table 2 Source parameters of the 2008 Ms8.0 Wenchuan earthquake

从表 2中看出，滑动弱化距离D_c 约为0.47m, 占平均滑动位移2.23 m 的21%.地震辐射效率ηS的值约为36.6% ，地震效率η近似等于13.7% .但是考虑到断层的平均震前、震后剪应力比400m深处震前、震后剪应力大，尽管理论上地震波辐射能E_S 和应变能变化E_eff随着深度的增加没有变化，但是在断层深部摩擦产生的热能比400m 处区域产生的热量多.那么按照地震效率的公式η = E_S/(E_eff +E_heat)，我们在400m 处求得的地震效率应为实际地震效率的上限.在给定上述两种不同地震矩M0 的情况下，采用E_f= (τ f-τ_f)DA，其中D为平均滑动位移，τf 为平均摩擦应力，A为断层面面积，那么有E_f/E_S = (τ f-τ_f)/σa, 由此可得断层从破裂到滑动全过程的平均摩擦应力τf 约为3.2 MPa.McGarr^{[24， 25]}给出了地震效率的公式η = (α-1)(1-2γ)/[1+α(1-2γ)]，其中，α =τ1/τf 表征了静态摩擦系数与动态摩擦系数之比μs/μd；γ= (τ f-τ2)/(τ1-τ2)为动摩擦应力上调系数，有α ≈1.6，γ ≈0.3；由此可得地震效率为η=14.6%，这与我们得到的13.7%的地震效率上限结果基本一致.对比1999年M_W7.6中国台湾集集地震，Wang^[26]利用不同的方法得到了该地震的地震效率，最高值达26.2%.因此，我们得到的13.7% 的地震效率上限也是可以接受的.同时，也说明汶川地震是一个中低效率的地震.根据Zller等^[27]给出的动态摩擦应力上调率:B= (τ1-τ2)/(τ₁-τ_f)，求得B=1.17，接近于Zller等^[27]在各向异性介质的蠕变区域对余震研究过程中使用的B=1.25的上调率.也就是说汶川地震约有17%的应力上调.那么，总的应变能变化E_eff有一部分通过地震波的形式辐射到地表，一部分用于破裂尖端的向前扩展，还有接近于地震波辐射能E_S 的能量作为松弛功E_r 贮存在断层周围的介质中.Stein 和Liu^[28]将美国加州LomaPrieta(板块边界上)地震、内华达州DixieValley和FairviewPeak地震(板块边界带)同美国中部New Madrid地震带(内陆)的余震进行对比，利用滑移速率和状态变化的摩擦定律(rateandstatevariablefrictionlaws)^[29]发现:在板块边界的断层上，板块移动速率较大(多大于10mm/yr)，余震持续时间为几十年；在板块边界带内，板块移动速率较小(一年几个毫米)，但具有较长的余震序列；发生在大陆板块的余震序列持续时间最长，其板块移动速率小于1 mm/yr, 例如美国1812年的New Madrid地震，余震序列能够持续几百年或更长时间.由于汶川地震发生在青藏高原北部，是沿龙门山断裂向四川盆地推覆逆冲造成的，GPS数据给出大地形变率小于3mm/yr^[30]，这就造成了断层的加载非常缓慢.如果Stein 和Liu^[28]对于余震序列持续时间的结论是合理可信的，那么发生在大陆板块内部的低板块移动速率区(小于3mm/yr)的汶川地震余震序列也能够延续几十年甚至上百年.而储存在断层周围的松弛功就通过几十年甚至上百年的余震或缓慢蠕变最终得到释放.假使E_r 取表 1中的平均值1.13×10¹⁶J，余震的能量释放遵从同主震一样的分配规律，余震的地震效率上限约为13.7%，那么全部余震所释放出的地震波辐射能为1.13×10¹⁶×13.7% ≈1.55×10¹⁵J.我们得到地震能量分配比例E_r∶E_S∶E_f 约为1∶1.2∶1.1.根据Richter^[31]给出的震级-能量关系式，lgE_S =1.5M_S+4.8，其中M_S 为面波震级，可知面波震级相差一级，地震波辐射能相差10^1.5倍(31.6倍).当汶川主震的全部松弛功都用于发生一个余震时，可以产生一个约为M_S6.9的余震.由于主震后不可能只有一个余震，况且我们使用的13.7%是地震效率上限，所以可以预计在汶川地震主断层上，不会发生大于M_S6.9 的余震.而根据目前统计，汶川地震最大余震震级仅为6.4级，还有16个大于等于M_S5.5的地震(如图 4).除了一个5.7级的地震离断层较远外，其余15个M_S5.5以上的地震均分布在断层带内，由此推断这些余震的产生属于静态触发，而且余震震级与我们的计算结果非常吻合.余震多发生在主震过后几天时间内，这体现了主震过后周围介质的应力调整，同时也是松弛功通过余震的形式耗散的过程.

图 4

Fig. 4

图 4 汶川地震主震和MS4.5以上余震分布图(USGS公布数据).黄色线为断层出露地表位置，黄色圆圈为汶川地震主震，红色圆圈代表MS5.5以上余震分布，黑色圆圈代表小于MS5.5余震分布 Fig. 4 Map of seismicity for the Wenchuan earthquake.Yellow circle indicates epicenter of main shock.Red circles arethe epicenters of aftershocks above Ms5.5，and black circles are the aftershocks below Ms5.5

地震是一个复杂的动态应力变化过程，其能量的求取过程也有很多不确定因素.由于不同学者对于断层的形态和展布有不同的认识，造成同一个地震静态应力降的估计有所差别.由于远场地震波辐射能的估计不够精确，导致了对视应力的估计也存在误差，从而影响了对于辐射效率和地震效率的估计.此外，对青川断裂带的四个测孔的应力测量，并不能完全等同于断层面上的应力状况，尽管有些研究认为断层带是脆弱的，在0~20km 深度的平均应力降在约为10 MPa以下，但是还有一些观点认为断层的强度很大，应力降随着深度线性增加，可以达到100 MPa^[32]；尽管绵竹清平台站距离断层仅仅5km, 但是观测得到的强地面运动地表速度和加速度并不能完全代表断层面上质点速度和加速度状况.这些都增加了对于地震能量估计的不确定性.但总体上说，我们利用不同的模型得到了相似的应力降结果，也定量化地估计了汶川地震的能量分配.

利用龙门山断裂带内400m 深钻孔的数据，求得了静态应力降Δτ_s 为2.7 MPa, 与位错模型求得的静态应力降2.6 MPa相似，且略大于由绵竹清平站强地面运动速度和加速度求得的动态应力降的平均值2.3 MPa.大地测量方法可以被用来测量近断层应变变化，取最具代表性的10^-4，乘以地壳弹性模量约为30GPa, 得到静态应力降约为3 MPa^[32].这与我们求得的2.6 MPa相似.这充分表明，汶川地震主震服从应力上调机制.此外，还计算了在不同的地震矩M0 下的破裂能和松弛功，得到了地震辐射效率和地震效率上限.并且定量地计算了全部余震的总地震波辐射能.假设余震与主震的能量分配方式完全相同，汶川地震的余震释放的全部能量相当于一个M_S6.9的地震.由此可推测汶川地震主断层上不会发生M_S6.9以上的余震.

感谢各位审稿人对于本文提出的宝贵意见和建议，使本文在问题的描述上有了更大的改进，加深了我们对于科学问题的理解.