汶川巨震发生在龙门山断裂带上,龙门山断裂带作为川—青地块与扬子地块的碰撞边界,位于青藏高原东缘,呈NE—SW走向,断层性质为逆冲兼右旋走滑(王二七等,2001)。它主要由3 条近乎平行的主干断裂组成(图1),3 条主断裂之下为朝NW方向缓倾的、统一的逆冲滑脱面,深度大约15～25 km(王椿镛等,2003)。龙门山断裂带是四川地区的主要活动断裂带,1657年汶川曾发生过6.5 级地震(唐荣昌等,1993)。我们采用美国地质调查局(USGS)地震目录和矩心矩张量解以及《2008年5月12日汶川特大地震震源特性分析报告》 陈运泰,许力生,张勇,等.2008.2008年5月12日汶川特大地震震源特性分析报告.中的相关数据,在图1中给出了汶川主震和主震后15天内发生的7次M_W≥5.5强余震的震源机制及震中位置分布图。

图1 汶川主震及其强余震震源机制及震中位置分布图

根据《2008年5月12日汶川特大地震震源特性分析报告》^①,汶川大地震是一次以逆冲为主、兼少量右旋走滑分量的地震。断层向WN方向倾斜,走向为229°。汶川主震能量辐射源跨度达600 km,但主要能量释放发生在前80 s,破裂开始于汶川县震中位置,然后以3.1 km/s的平均速度向EN方向传播了300 km,最后,破裂反向传播,在震中西南有微弱的能量辐射。破裂大致由3个主要的子事件组成(表1),第1个子事件位于0～15 s,这个时间内释放了整个过程16%的地震矩; 第2个子事件开始于第17 s,结束于第46 s,为本次地震过程中最大也是最主要的一次事件,释放了63%的地震矩; 最后一次事件在第48 s和第60 s之间,同样释放了16%的地震矩。第60 s以后,总体上整个破裂已经基本完成,仍有一些地方发生断断续续的破裂,这些破裂只占所有地震矩的5%。因此,此次地震的主要能量释放时间大约为60 s。本研究中用于计算的断层长度约为200 km,即采用了主要能量释放阶段的断层长度。

表1 汶川主震破裂过程3个子事件的参数

根据《2008年5月12日汶川特大地震震源特性分析报告》^①给出的相关参数和USGS提供的矩心矩张量解数据,同时参考王卫民、姚振兴研究小组的主震模型参数 王卫民,姚振兴.2008.汶川8.0级地震震源破裂过程. http://www.csi.ac.cn.,我们列出了汶川主震和7次强余震的震源位置和震源断层面参数(表2)。我们取节面走向与宏观破裂方向即余震震中展布方向较接近的节面为余震的发震断层面。

表2 汶川巨震及其强余震震源断层面参数

注:*表示陈运泰等的参数^①; **表示王卫民等的参数^②。

我们以二维情况说明DWN法的原理(Bouchon,2003)。无限均匀介质中一个线源的稳态辐射,可以用柱状波表示,也等价于均匀与非均匀平面波的连续叠加。设x表示水平坐标轴,z表示垂直坐标轴,则位移或应力可以表示为

F(x,z; ω)=e^iωt∫^∞∫_-∞f(k,z)e^-ikxdk.(1)

其中ω是频率,k是水平波数。当介质有限或垂直方向不均匀时,积分存在奇点,因此对水平波数的积分变得复杂而困难。DWN法首先用复合源问题代替上述单一源问题,则(1)式变为

G(x,z; ω)=∫^∞∫_-∞f(k,z)e^-ikx∑^∞_m=-∞e^ikmLdk.(2)

G(x,z; ω)即为格林函数,其中L是周期源的间隔。于是(2)式可进一步化为

G(x,z; ω)=(2π)/L∑^∞_n=-∞f(k_n,z)e^-ik_nx.(3)

式中k_n=(2π)/(L)n。

如果上式中的序列收敛,则可以表示为有限求和形式:

G(x,z; ω)=(2π)/L∑^N_n=-Nf(k_n,z)e^-ik_nx.(4)

DWN法的下一步是将已在频率域中解决的复合源问题恢复到单一源问题。如果能计算函数G的连续傅立叶变换,问题将变得较直接,但实际上仅能对某些频率计算出函数G,并使用离散傅里叶变换得到时间域的解。因此,一方面,要处理一个含无限时间响应的单一源问题,另一方面,要使用离散傅立叶变换得到有限持续时间为T=2π/Δω的信号。这些可通过在复频率域中作傅立叶变换实现:

g(x,z; t)=∫^{∞+iω _l}∫_-∞+iωlG(x,z; ω)e^iωtdω.(5)

式中ω_l表示频率的虚部,且要满足

e^ω_lT<<1.(6)

(6)式保证了在时间窗T内的衰减。因此,如果选择L足够大,以至在时间窗T内没有从最邻近的源到达接收点(x,z)的干扰,时间域的单一源的解f(x,z; t)便可以从频率域的复合源的解G(x,z; ω)恢复得到:

f(x,z; t)=e^-ω_lt∫^∞∫_-∞G(x,z; ω)e^iω_Rtdω_R.(7)

任何类型的弹性震源都可以用一组点源的组合来代替。地震学中最感兴趣的是剪切位错产生的辐射场,而一个双力偶震源就相当于一个弹性剪切位错源。设有一剪切断层,用(s_x,s_y,s_z)表示沿断层滑动方向的单位矢量,用(n_x,n_y,n_z)表示断层面法线方向的单位矢量,则相应的矩张量分量为

m_ij=-μslip(ω)ΔS(s_in_j+s_jn_i).(8)

其中ΔS是发生滑动的断层平面的面积元,slip(ω)是滑动矢量。将地震断层面离散成小面元上位错点源组成的二维阵列,则每个小单元面上的位错产生的位移场等价于在小单元面处双力偶点源所产生的位移场。每个点源以e^-iωt_r的相位延迟辐射,其中t_r表示破裂从小点源传播到断层某个确定位置所需要的时间,所有小点源的贡献在频率—波数域内求和。

分别计算汶川主震3次破裂在强余震断层面上产生的完全库仑破裂应力变化,然后将其叠加,可得到强余震断层面上总的完全库仑破裂应力变化。计算步骤如下:

首先采用离散波数法合成位移理论地震图计算程序,计算汶川主震在其余震接收点(余震震源)处及其邻近3个点处产生的地震波位移u_i(x,t); 然后应用差分法原理由位移计算出应变变化分量Δε_ij(x,t),根据虎克定律将应变转换成应力,可得到主震在余震接收断层处产生的应力变化分量Δσ_ij(x,t)。

设接收断层面的法向单位矢量为n(n₁,n₂,n₃),滑动方向单位矢量为s(s₁,s₂,s₃),由柯西公式可得接收点处断层面上的动态应力变化矢量:

ΔT_i(x,t)=Δσ_ij(x,t)n_j.(9)

将应力矢量分别投影到接收断层面的法线方向和滑动方向上,可得到正应力变化Δσ(x,t)和切应力变化Δτ(x,t):

Δσ(x,t)=Δσ_ij(x,t)n_jn_i,(10)

Δτ(x,t)=Δσ_ij(x,t)n_js_i.(11)

完全库仑破裂应力变化(Kilb等,2002)为

ΔCFS(x,t)=Δτ(x,t)+μ'Δσ(x,t).(12)

其中,μ'为接收断层面介质的视摩擦系数,一般取0.2～0.8,本研究取μ'=0.5。断层面上的正应力变化Δσ(x,t)拉张时为正、压时为负。

完全库仑破裂应力包括了动、静态库仑破裂应力(Kilb等,2002)。在以下表述中,动态库仑破裂应力变化采用符号ΔCFS(t)表示,静态库仑破裂应力变化采用符号ΔCFS表示。

图2给出了计算得到的汶川巨震在后续7个强余震断层面上产生的完全库仑破裂应力变化ΔCFS(t)的时程曲线。

本文取余震震中位置误差为0.1°、震源深度误差为5 km,计算了定位误差对库仑破裂应力变化产生的影响。在距离主震震源100 km以内,峰值动态库仑破裂应力最大绝对误差≤1.36 MPa(绝对值),静态库仑破裂应力最大绝对误差≤0.19 MPa; 在距离主震震源100 km以外,峰值动态库仑破裂应力绝对误差≤0.93 MPa,静态库仑破裂应力绝对误差≤0.06 MPa。

考虑不同破裂模型对计算结果的影响,需要综合考虑破裂模型各种参数,包括断层走向、倾角、滑动角、长度、宽度等的误差(变化)对库仑破裂应力变化计算结果所造成的误差。模型参数误差绝

图2 汶川主震在其各强余震断层面上产生的完全库仑破裂应力变化的时程曲线(a)1号强余震;(b)2号强余震;(c)3号强余震; (d)4号强余震;(e)5号强余震; (f)6号强余震;(g)7号强余震

由图2可以看到,汶川主震在其强余震断层面上产生的完全库仑破裂应力变化ΔCFS(t)的时程曲线由起伏振荡部分和尾部平稳部分组成,它们分别代表了动态库仑破裂应力变化和静态库仑破裂应力变化(Kilb等,2002)。图2显示,汶川主震在各后续强余震断层面上产生的动态库仑破裂应力变化绝大多数达到了10 MPa量级,静态库仑破裂应力变化绝大多数达到了兆帕量级。而Kilb等(2000)和Harris(1998)认为,有意义的动态库仑破裂应力触发和静态库仑破裂应力触发阈值分别为0.5 MPa和0.01 MPa。

根据上述分析,我们对汶川巨震在各强余震断层面上产生的完全库仑破裂应力变化的作用做如下讨论:

1号强余震(M_W5.7)相对其它强余震,其断层面上的动态和静态库仑破裂应力幅值最大(图2a),而且均为正值,远远超过动态和静态触发阈值,因此,1号强余震充分受到了汶川主震动态和静态应力的触发作用。

2号强余震(M_W6.1)断层面上正值部分的动态库仑破裂应力和静态库仑破裂应力均超过了相应的触发阈值(图2b),也受到了汶川主震动态和静态应力的触发作用。但其断层面上的动态和静态库仑破裂应力相对较小,这可能是因为2号余震的位置背于汶川主震破裂的主要传播方向。

3号强余震(M_W5.3)断层面上的静态库仑破裂应力幅值相对较小(图2c),约为0.7 MPa,刚刚超过静态应力触发阈值。动态库仑破裂应力低频振荡明显,正值部分的动态库仑破裂应力超过了相应的触发阈值。因此,3号强余震受到了汶川主震静态和动态库仑破裂应力的触发作用。

4号强余震(M_W5.6)断层面上的动态和静态库仑破裂应力数值较大(图2d),仅次于1号强余震,数值均为正值,远超过触发阈值,因此该强余震受到了汶川主震动态和静态应力的触发作用。另外,4号强余震断层面上的动态库仑破裂应力曲线显示出了汶川主震3次破裂产生的3次动态应力脉冲序列的特点。

5号强余震(M_W5.7)、6号强余震(M_W6.0)和7号强余震(M_W5.5)的震中位置距离汶川主震较远,其断层面上的动态和静态库仑破裂应力幅值均相对较小(图2e、f、g),但正值部分的动态库仑破裂应力均超过了动态应力触发阈值,静态库仑破裂应力也均超过了相应的触发阈值,因此这3个强余震均受到了汶川主震动态和静态库仑破裂应力的触发作用。

正的静态库仑破裂应力变化对断层有促进滑动的作用,负的静态库仑破裂应力变化对断层有阻碍或延迟滑动的作用。对于快速的动态库仑破裂应力振荡,从效果上讲,负的动态库仑应力变化对断层滑动没有大的促进作用,但也没有阻碍或延迟作用,或者说,断层对负的动态应力脉冲不敏感,正的动态库仑破裂应力变化比负的同样大小的应力变化对断层的力学作用要大得多。这是因为作用在断层上的瞬态的应力会引起断层的不可逆弱化,即使破裂发生愈合,但由于愈合过程太慢,以至于不能阻止或抵消如亚临界裂隙生长这样的弱化结果。岩石的滑动速率和状态依赖性本构特征控制着断层摩擦滑动不稳定性或成核过程,动态应力可以通过改变断层成核区的特征滑动距离来影响滑动速率及状态,快速的正的动态应力变化对断层的弱化作用比负的动态应力变化对断层的作用要大得多(Kilb等,2000; Gomberg等,1997)。

断层一旦受到这种力学作用其力学性质会发生改变,但并不代表破裂会立即或在一定的时间内发生,这就意味着在触发与被触发的地震之间存在不确定的时间延迟(Kilb等,2000),被触发地震的时间延迟不能用应力传播来解释,但可以用与断层摩擦滑动或成核过程相关的岩石速率—状态依赖性本构模型来解释。由于地震波可以通过改变断层成核区的特征滑动距离来影响滑动速率及状态,因此对伴随地震波产生的动态库仑破裂应力的延迟触发机理,能够用速率及状态依赖性摩擦滑动不稳定模式进行充分说明(Kilb等,2000; Gomberg等,1997; Parsons,2005)。Hill等(1993)和Anderson等(1992)认为,动态应力触发的最大延迟时间为15天,超过15天发生的后续地震受到的动态库仑破裂应力触发作用是无效的。

综上所述我们可得出以下几点结论:

(1)汶川主震后15天内发生的7次M_W≥5.5强余震全部受到了动态库仑破裂应力和静态库仑破裂应力的综合触发作用。

(2)与仅有单次破裂过程的一般强震不同,汶川巨震在其各强余震断层面上产生的完全库仑破裂应力变化的时程曲线,基本上都反映了主震破裂过程3个阶段的特征,动态库仑破裂应力的最大峰值多出现在第2阶段。巨震产生的库仑破裂应力,无论是动态的还是静态的,都远远超过以往研究中所总结出的触发阈值。

(3)在短时间内,巨震产生的动态库仑破裂应力对后续地震的触发效应可能大于静态库仑破裂应力,高强度、高频率的动态库仑破裂应力主要依靠改变(降低)后续地震断层面上的摩擦滑动强度来触发地震。