1996年2月3日19时14分,云南丽江发生M_S7.0地震,美国国家地震信息中心测定为M_W6.5,震中位于丽江县金沙江边的大具乡。这次地震是滇西北地区自1976年宁蒗6.7级地震以来20年间发生的最强烈的一次地震。地震有感范围很广:北起四川甘孜,南至云南思茅,东至昆明。受灾面积达18 720 km²,丽江地区、迪庆藏族自治州、大理白族自治州和怒江傈僳自治州境内都遭到不同程度的震害,丽江县城及附近地区约20%的房屋倒塌。受灾乡镇51个,受灾人口达107.5万,重灾民有30多万。人员伤亡人数为17 221人,其中309人丧生,3 925人重伤。电力、交通、通讯以及水利等设施也遭到了严重破坏,全县直接经济损失42.665 7亿元人民币。

此次地震为主震—余震型,截至1996年10月1日,共发生1.0级以上余震5 493次,其中6级1次,5.0～5.9级地震6次,4.0～4.9 级地震39次,3.0～3.9级地震210次。余震活动主体沿丽江盆地西部呈近南北向展布,带长约60 km,带宽20 km,余震有由北向南迁移的趋势(皇甫岗,1997)。主震距丽江县城40 km,最大余震南迁约30 km,位于县城北北东方向10 km处。

大地震的震源破裂过程蕴含了地震发生时能量释放的强弱,有助于解释地壳构造学领域的研究结果,而且能使人们更直观地理解断裂摩擦性质和地震发生的过程。若能在大震发生后快速得出震源破裂模型,就能迅速掌握到震区受灾程度,从而清晰地指引震后应急救援工作有序展开。

得到震源破裂模型的方法有很多:(1)通过近场GPS和InSar观测,这一方法可提供较高的空间分辨率解,但往往需要较长时间,且大多数情况下缺乏此类资料;(2)用远震波形开展有限元破裂过程反演,该方法能得到破裂过程中比较详细的细节,但受很多基于破裂过程中的动力学参数的限制;(3)近年来发展起来的反投影远震P波法(XU et al,2009; 徐彦等,2011; Ishii et al,2005; D'Amico et al,2010),该方法可以对在破裂过程中的动力学参数进行最小限度约束的基础之上得到震源破裂的时空变化,并为有限元法提供动力学参数的约束范围。与传统的有限元反演方法相比,反投影远震P波法更快更直接,因为该方法只需要很少的信息:一维速度模型和震中信息。其中全球范围内一维速度模型是现有的,如PREM(Dziewonski,Andeeson,1981),AK135和IASP91(Kennett,Engdahl,1991,1995); 目前地震台网的监测能力已经能够在震后10 min内给出震中信息。本文运用反投影远震P波方法研究1996年2月3日的丽江地震,得到该地震的震源破裂过程。

丽江盆地为狭长形断陷盆地,南北长40 km,东西宽3～11 km,可分为南北两部分,其中北部的断陷深度在中海一带可达1 200 m,南部的第四纪断陷深度在连汪—丽江城区之间最深处可达700 m(韩竹军等,2004)。由于印度板块与欧亚板块的碰撞作用,青藏高原东部地区地块整体强烈隆起,在此背景上一些断裂出现差异性抬升,形成有明显反差的高山峡谷与断陷盆地地貌。丽江地震震区地质构造十分复杂,南北﹑北西﹑北东走向的断裂纵横交错(图1)。这些断裂控制着规模不等的第四纪断陷盆地的发育,位于本次地震震中附近的大具乡是北西和北东两条断裂的交汇点; 位于极震区的丽江盆地东西侧均受到南北向断裂的控制,其东北和南缘均有断裂通过。该区及其附近地区历史上也曾发生过强烈地震。1966年9月28日曾在大具乡西北约20 km处发生M_S6.4地震。

1996年丽江7.0级地震地表破裂带、余震基本上都是沿着南北走向的玉龙雪山东麓断裂分布。该断裂是丽江第四纪断陷盆地与玉龙—哈巴雪山之间的边界断裂,与地表地质特征相对应,丽江地震震源破裂也表现出较大的倾滑分量,这与滇西地区以走滑为主的震源破裂方式存在明显差异(秦嘉政等,1997)。张建国等(1997)实际调查,表明此次地震地表破裂带的规模及其位错幅度都明显偏小,地震地表破裂带在黑水断层谷东西两侧的正断层附近都有所反映(韩竹军等,2004),地震断层的平面展布及其力学性质与玉龙雪山东麓断裂十分吻合。

笔者首先通过IRIS http://www.iris.edu.网站选取并下载了丽江7.0级地震的75条宽频带垂直向记录,其震中距范围在30°～95°,因为在此范围内的地震波主要是在介质相对均匀的下地幔中传播,避免了地震波在上地幔和核幔边界传播时因介质的非均匀性导致波的复杂性,从而使由传播路径所造成的波的复杂性达到最小化。然后使用SAC软件逐条查

图1 研究区域及1976年以来M≥4.0地震震源机制解
(F1:怒江断裂; F2:澜沧江断裂; F3:维西—巍山断裂;
F4:金沙江断裂; F5:龙蟠—乔后断裂; F6:洪门口断裂;
F7:大具断裂; F8:玉龙雪山东麓断裂; F9:鹤庆—洱源
断裂; F10:小金河断裂; F11:博科—木里断裂;
F12:金河—箐河断裂)
Fig.1 Map of the study area with focal mechanisms of M>4 earthquake since 1976 (F1:Nujiang River Fault; F2:Lancangjiang Fault; F3:Weixi-Weishan Fault; F4:Jinshajiang River Fault; F5:Longpan-Qiaohou Fault; F6:Hongmenkou Fault; F7:DajuFault; F8:Eastern Piedmont Fault of Yulong Snow Mountain; F9:Heqing-Eryuan Fault; F10:XiaojinRiver Fault; F11:Boke-Muli Fault; F12:Jinhe-Jinghe Fault)

台站的分布情况对反投影结果的影响可通过台站响应函数(ARF)(Yan et al,2009)来体现,ARF能反映出由于台站分布不均而导致慢度域上地震波能量的扩散和泄露。若台站分布是均匀的,在发震时刻最大能量应在震中位置呈现二维Delta函数状分布,也就是最大能量应位于震中位置,并且向四周逐渐衰减。ARF的表达式为

ARF(θ,φ,h,f)=|1/N∑^N_j=1ω_je^{2πifΔt(θ,φ,h)_j}|².(1)

式中,Δt(θ,φ,h)_j=t(θ,φ,h)_j-t(hypocenter); θ是经度; φ是纬度; h是深度; f是频率; N是台站数; t(θ,φ,h)_j是从震中附近一点(θ,φ,h)到第j个台站的传播时间; ω_j是最大值为1的正数,代表各台站的权重。台站响应函数是一个正实数,理论上最大值位于震中。

由于目前全球地震台网分布不均匀,在丽江地震震中距30°～95°范围内大部分的地震台站集中在欧洲、北美和澳大利亚。笔者选取远震P波较敏感的频率(1 Hz)作为中心频率,计算由42个台站所组成的全球地震台网的ARF,结果显示出由于台站分布的间断和非均匀造成了能量的泄露。计算得出的ARF最大能量位于震中位置,但在图示区域范围内还有很多局域最大值位于震中位置之外。这些局域最大值会很大程度地影响到反投影远震P波的结果,这种现象称为旁瓣效应。可通过对台站进行重新选取来减小旁瓣效应的影响,通过对台站赋予不同的权重来对台站进行重新选取。为了找到最近似于Delta函数状ARF的台站组合,本文采用了类似于重采样的方式来改善ARF。经过选择,最后采用的赋值方法是把地球表面平均地分成700 km×700 km的块,从每个块中随机选取1个台站,把该台站的权重值赋为1,而块内其余台站的权重值均赋为0。经过重新选取,笔者得到27个台站,由27个台站所得到的ARF有了明显改善,局域最大值明显减少,旁瓣效应明显减弱(图2)。

用反投影远震P波研究震源破裂过程是基于波形相关性的方法。在本文中我们运用MCCC法(Multichannel cross correlation)(VanDecar,Crosson,1990)计算了所选的27个台站记录的P波波形相似度。根据MCCC方法,一个台站记录到的P波波形相似度值为该台站与其余台站形成的台站对的波形相似度的平均值。图3a为27个台站的P波相似度值,由图可见,有26个台站的P波相似度在0.6以上。我们把这27个台站的未经滤波的原始P波记录画在一起(图3b),同样可以直观地看到这27个台的P波记录包含有这相同的信息,即27条宽频带垂直向记录承载了相同的震源信号。

图2 选取的42个(a)和27个(b)全球地震台网的台站记录所得到的台站响应函数
Fig.2 ARF of waveform recorded by the selected 42(a)and 27(b)stations in Global Seismic Network

图3 所选用的27个台站的P波相似度值(a)和未滤波的P波记录(b)
Fig.3 Similarity values of P wave(a)and unfiltered P records(b)of the selected 27 stations

除了台站的选取,反投影方法还有一个重要技术点在于准确地得到从震中附近一点到远震距离上的地震台站的P波传播时间。径向地球速度模型因为三维地球结构的变化会造成传播时间的差异,为了减少这一差异的影响,对所有台站记录到的P波的前10 s记录进行校准,允许每一台站的P波前后移动,从而使所有台站的P波在同一时刻到达。假设P波移动时间变化很小,并且到时校正值被用于震源区域的所有网格点。反投影研究中重要的参数是振幅,在测量P波到时校正值的同时,我们得到了P波振幅的归一化参数,归一化参数可有效地去除台站场地、波的几何扩散、仪器放大值不同以及地震波辐射方向的影响。

反投影方法是在某一特定时间,通过对与某一可能的震源位置所对应的波形进行叠加来抵消噪音和传播路径中次生波的影响,从而突出从震源传出的信号,然后把叠加所得到的能量投影到与之相对应的震源位置。在对可能震源区域所有可能位置都进行了能量反投影后,得到该时间的震源图像。之后把这一处理过程运用到从震前到震后的一个连续时间段上,从而得到该地震全时间段的震源破裂过程。在进行反投影时,考虑一个均匀的四维空间(即经度、纬度、深度、时间)围绕在震源和发震时刻周围。对每一网格点都计算它与每一台站相对应的理论到时,然后对用带宽0.5～1.5 Hz的Buterworth滤波器滤波后的波形运用4阶方根叠加算法(XU et al,2009; 徐彦等,2011)进行叠加:

B'(t)=1/M∑^M_j=1|b_j(t)|^1/n·sign{b_j(t)},(2)

B(t)=|B'(t)|^Nsign{B'(t)}.(3)

其中,B(t)是最终的叠加值; B'(t)是叠加中间值,b'(t)是第j条记录的振幅; M是波形总数。当N=1时,为线性叠加。由于线性叠加虽运算速度较快,但却不能很好地提高信号能量并抑制噪音的干扰,本文选用N=4时的4阶方根叠加算法是而运用高阶方根叠加则能很好地达到这个目的。笔者采用20 s的窗长,1 s的滑动窗来对总长100 s的数据进行研究,叠加值的能量被反投影到与之相对应的网格点上,最大能量值所在的时间和空间区域为破裂区域。

笔者通过ARF和MCCC对所下载的75条宽频带垂直向记录进行筛选后,最终运用P波波形相

图4 归一化的时间能量变化图
Fig.4 Normalized variation graph of beam power as a function of time

图4为丽江地震能量随时间变化曲线。在图中可以看到能量值在0即发震时刻之前就有逐渐上升的现象,这是由于所采用的滑动时间窗窗长为20 s,可以使时间窗中心点在发震时刻之前的时间

窗就包含有后面的能量,导致能量的明显上升发生在0时刻之前。能量最大值随时间的变化曲线上显示出,在震后30 s能量值稳定于0.05附近波动。虽然这一能量波动范围高于震前的能量水平,但这一现象来自于大震后尾波能量的干扰,通常大地震后背景噪音至少需要30 min才能恢复到震前的水平(Kaema,Ringdal,1999),因此,本文得到丽江地震震源破裂时间总长度约为30 s。丽江地震震源破裂包含了两个主要能量释放点(图5),

图5 (a)发震时刻,(b)震后5 s,(c)震后15 s和(d)震后28 s的震源破裂过程截图
Fig.5 Snapshots of source rupture process at the origin time(a),5s(b),15s(c),28s(d)after the origin time

一个是在震中位置附近,震后约5 s的时间(图5b),该点也是整个破裂过程中能量释放最大的; 另一个在震后约15 s处(图5c)。从时间上看(图4),两次主要能量释放的时间差约为10 s; 从空间上看(图5),此次地震能量的释放有自北向南传播的特征,第二个主要能量释放点位于震中位置以南约20 km。许力生等(1997)用经验格林函数反卷积的办法估计了丽江主震的震源时间过程,得到两源的时间差约为10 s,马淑田等(1998)运用同时反演两个点源震源机制的原理,得到丽江主震由两次主破裂构成,两次破裂时间差约为12 s,相距约26 km。本文得到的两次破裂时间差为10s,空间相距20 km; 本文结果与许力生等(1997)的结果对比,反投影远震P波法得到的丽江地震震源破裂过程在时间上是一致的,两次破裂时间差都是10 s; 与马淑田等(1998)的反演结果进行对比,两种研究震源破裂过程的地震学方法得到的丽江地震震源破裂过程基本一致,但也存在一些细节上的差异,造成这些差异的原因可能是两种方法所研究的频率范围不同。XU等(2009)对2008年汶川地震的研究表明,破裂过程的一些细节与所研究的频率范围有着直接联系,研究选取的频率不同时,所得到的震源破裂过程主要特征虽不受影响,但在某些细节上会有所不同。

图6为运用反投影P波得到的丽江地震能量释放积累图,包含了所选取的100 s的时间长度的结果。由图中可以看到,在空间上,丽江地震的主要能量释放呈近乎北南向分布,沿玉龙雪山东麓断裂发展,破裂长度约为40 km,在中段10 km左右只有很少的能量释放,而南北两端存在着很强的能量释放点。这与张建国等(1997)考察得出的丽江地区地表破裂特征相符。值得注意的是,第二次能量释放过程中,除了震中以南这一能量释放点之外,在震中的西南方(主破裂点以西约15 km处)也几乎同时出现一个相对弱的能量释放点(最大释放能量只相当于最大能量释放的40%),由于这一能量释放区域正好位于龙蟠—乔后断裂,笔者认为,丽江地震在时间上有两个主要的能量释放点,而在空间上则有3个主要的能量释放,最大能量释放时空点位于震后5秒的震中位置,而震后15 s时,有两个主要的能量释放点在同时进行,位于玉龙雪山东麓的能量释放点所释放的能量约为龙蟠—乔后断裂上能量释放的2倍。

通过反投影远震P波,笔者得到发生在1996年2月3日的丽江地震震源破裂时间约为30 s,破裂长度约40 km,破裂主要沿玉龙雪山东麓断裂发展。此次地震的能量释放过程主要体现在两个时间点上,第一个能量释放点是在震后5 s,另一个能量释放点是在震后15 s。而在空间上,丽江地震有3个能量释放点,最大能量释放点为震中附近,次大能量释放点位于震中以南约20 km的玉龙雪山东麓断裂处,第3个能量释放点在时间上与次大能量释放同步,而位置位于震中西南的龙蟠—乔后断裂上。

反投影远震P波技术多被用于7.5级以上的大震,然而D'Amico等(2010)运用此方法研究了2009年4月6日的意大利M_W6.3地震,得到了较好的结果,表明反投影P波法除了能运用于有着上百公里破裂尺度的大震外,对有重要破裂特征,破裂尺度在数十公里范围的中强地震也同样适用。本文对于丽江地震的研究也同样体现出反投影P波法能够较好地得到破裂尺度在数十公里范围的

图6 丽江地震能量积累图
Fig.6 Cumulative beam power of Lijiang M7.0 earthquake

目前由软件自动定位生成的震动图不能充分反映出地面运动的信息,因而反投影P波法得到的震源破裂模型可以作为震动图一个补充。虽然此方法不能像反演方法那样得到很详细的破裂图像,但由于反投影远震P波法对破裂过程中的动力学参数限制很少,运算快,因此能在震后较短时间内得到震源破裂过程,从而能为一般速报只提供震中位置的监测手段作重要补充,继而为震后应急救援工作提供依据。