2011, Vol. 54
2010年4 月14 日当地时间07 时49 分，青海玉树藏族自治区发生了7.1级地震.截止2010年7月13日，此次7.1级地震导致玉树地区2220人遇难，失踪70人(中国地震信息网).据中国台网中心测定，震后3个月时段内玉树地震序列共发生余震2847次，其中M6.0~6.9级1次，M5.0~5.9级2次，M4.0~4.9 级9 次，M3.0~3.9 级24 次(中国地震台网中心).这次玉树地震为中国大陆自2008年汶川地震后的又一次造成大量人员伤亡的破坏性地震，也是玉树地区近百年来的最大地震.
一个大地震发生后，地震监测部门的一项重要职责就是尽快地得到震源破裂模型.这一破裂模型所包含的区域最大震动信息能够指导地震应急救援工作，这是由于在大多数情况下，震中并不是破坏最重的区域，例如汶川地震，其破坏最重的地区和震中相距几十公里.除去震后快速应急救援工作对震源破裂过程研究的需求外，震源破裂过程还有助于解释地壳构造学领域的观测和研究结果，并能帮助加深对断裂摩擦性质和破裂过程的理解.
大地震的破裂模型通常可通过三种方法来得到：一是GPS或InSAR 观测；二是有限源反演远震波形；三是反投影远震P 波.然而在震后数分钟内就能开始的对大地震震源破裂的快速研究迄今为止仍来自于地震观测，这是由于当体波和面波传播到达分布于全球的地震台站，且被实时地传回数据中心时，大量的地震波数据就已经能被用来对震源进行研究了.
有限源反演能得到破裂过程中比较详细的细节，然而有限源模型能够给出较好波形模拟的基础在于对破裂过程中的动力学参数的限制[1].近年来发展起来的反投影研究破裂过程的方法[2~6]能够在最小限度内对破裂过程中的动力学参数进行限制的基础上得到震源破裂的时空变化.反投影方法可以为有限源方法提供动力学参数的约束范围.这是一种和传统的有限源反演相互补充的研究方法，但却相对更直接且更快(能在震后约30min的时间内得到地震的破裂模型).这是由于反投影远震P 波只需要很少的信息：一维速度模型和震中信息.全球范围的一维速度模型已知，如PREM[7]，AK135[8]和IASP91[9]，且我们台网的监测能力能够在震后10min内给出较为精确的震中位置，这都为运用反投影远震P波快速得到震源破裂过程创造了条件.本文将基于反投影远震P 波研究震源破裂模型的方法运用于4 月14 日玉树地震，研究玉树地震破裂过程.
2 研究区域青海玉树地区位于青藏高原东部(图 1).始于约5000万年前的印度板块和欧亚板块间的碰撞造就了平均海拔5000 m 的世界屋脊———青藏高原[10].印度板块持续北向运动的挤压力转换为了欧亚板块的内部变形[1~14]和大量的大型走滑断裂，如：阿尔泰断裂、昆仑断裂、卡拉昆仑断裂、海原断裂、鲜水河断裂.玉树地区平均海拔4500m，是两大水系(金沙江、澜沧江)的分水岭.本次玉树地震的发震断层为北西—南东走向的玉树断裂，该断裂位于鲜水河断裂系统的北西端，是这一重要的走滑断裂系统的组成部分.沿断裂发展的溪水错断、断层陡崖、拉张盆地、雁列式山脊都显示玉树断裂是一年轻的活动断裂[15].
|
图 1 研究区域及1976年以来的M>4级地震震源机制解 F1:喜马拉雅山前断裂；F2:石阶断裂；F3:红河断裂；F4:龙门山 断裂;F5 :玉树断裂；F6 :科多断裂；F7 :阿尼玛卿断裂；F8 :阿尔 金断裂.红色震源机制解为2010年4月14日玉树地震.黑色箭 头为印度板块移动方向. Fig. 1 Map of the study area with focal mechanisms of M>4 since 1976 Eight major faults are shown ： F1 ： Himalayan front fault； F2 ： Shijie fault； F3 ： Red river fault； F4 ： Longmen shan fault； F5 ： Yushu fault； F6 ： Keduo fault； F7 ： A’ nyemaqen fault； F8 ： Altun fault. Red beach ball represents April 14, 2010 Yushu earthquake. The black arrow is themoving direction of India plate. |
本次玉树地震为玉树地区自19世纪以来最大地震.1896年的7级地震为与玉树地震在同一断裂上的次大地震.中国地震局地球物理研究所刘超等(http://www.csi.ac.cn/manage/html/4028861611c5c2ba0111c5c558b00001/\- content/10\- 04/17/1271488288058.html)[2011-01]给出的矩张量解显示破裂发生在一走向东南(119°)，近乎垂直(倾角83°)的走滑断层面上.据地震现场地质灾害调查结果，玉树地震最大位移约1.8m，地表破裂带约为31km，走向310°.
3 数据和方法我们从IRIS(www.iris.edu/data)[2011-01]下载了381条宽频带垂直向记录.为避免波在上地幔和核幔边界传播时介质的非均匀性所造成的波的复杂性，我们所选取的381 条记录的震中距范围在30°~95°之间.这一震中距范围的地震波主要是在介质相对均匀的下地幔中传播，该条件最小化了传播路径所造成的复杂性.从IRIS下载了波形后，我们对所有波形逐一进行检查，去除波形记录不好的台站.经过这一波形筛查后，共206条宽频带垂直向记录被用于我们的反投影研究中.
台站响应函数(ARF)[6]给出了台站分布对反投影结果的影响.台站响应函数能够给出由于台站分布的非均匀性，在慢度域上地震波能量所显示出的扩散和泄露.如果台站分布均匀，在发震时刻最大能量应在震中位置呈现二维Delta函数状分布.即最大能量位于震中位置且能量向四周逐渐衰减.下式为台站响应函数(ARF)的表达式，
|
(1) |
其中，θ是经度，Φ 是纬度，h是深度，f是频率，N是台站数，t(θ，Φ，h)j 是从震中附近一点(θ，Φ，h)到第j个台站的传播时间.wj是一最大为1 的正数，代表各台站的权重.台站响应函数是一正实数，最大值位于震中.
由于目前我们的全球台网分布不均匀，在玉树地震震中距30°~95°的范围内大部分地震台站集中在欧洲、北美和澳大利亚，见图 2a.图 2 为运用式(1)得到的中心频率为1.0Hz的台站响应函数.我们选取中心频率为1.0 Hz是由于远震P 波对这一频率的信息较为敏感.这样的一个台站分布的台站响应函数直观地显示出由于台站分布的间断和非均匀所造成的能量分布的泄露.在图 2a中最大能量的确位于震中位置，但在图示区域范围内还有很多局域最大值位于震中位置之外.这些局域最大值会很大程度地影响到反投影远震P 波的结果.我们也可称这些局域最大值为旁瓣效应.
|
图 2 (a，c)381个全球宽频带垂直向记录所得到的台站响应函数.图a中的红色五角星为玉树地震震中，黑色点为台站 位置.台站响应函数中的绿色五角星为震中位置.（b，d)29个全球子台网的台站响应函数.图b中红色五角星为震中，黑 色点为台站位置，空心圈为6个权重为0的台站.台站响应函数中绿色五角星为震中位置 Fig. 2 (a，c) ARF of 381 global stations. Red star is theepicenter ofYushuearthquake，black dots are the stations.Green star m the ARF s the epicenter of Yushu earthquake. (b,d) ARF of 29 global subarray. Red star is the epicenter， black dots are the station，cicles are the stations with zero weight. Green star in the ARF s the epicenter |
然而台站响应函数可以通过对台站赋予不同的权重，即对台站进行选取，来进行改进.赋予台站不同的权重目的是为了找到能得到最近似于Delta函数状ARF的台站组合.我们采用类似于重采样的方式来改善台站响应函数.我们把地球表面平均分成700km×700km 的块，并从每个块中随机地选取1个台站，该台站的权重值被赋为1，块内的其他台站的权重均为零.经过这样的处理，我们得到76个台站，这76个台站所得到的台站响应函数得到了明显改善，局域最大值明显减少，旁瓣效应明显减弱.
反投影远震P 波研究震源破裂过程是基于波形的相关性上的方法.远震P 波的出射角度较小，例如一震源深度为10km 的地震在震中距为95°距离处的P 波的出射角度为13.8°，而震中距为1°的近场P波的出射角度为95°.远震P 波的能量近乎垂直出发，对于走滑型的地震，小角度出射的远震P波是从靠近两节面交叉的位置出发的，而两节面交叉的区域也是P波能量较弱的区域，因此在远震记录上信号相对较弱.而这也使得各台站记录到的P波的波形相似度较低.图 3a为运用MCCC(Multichannel crosscorrelation)[16]对76 个台站记录计算的P 波波形相似度和76条未经过滤波处理的P波记录.在运用MCCC 计算波形相似度时，我们首先截取P 波前5s到P波后10s共15s长的未滤波的波形.然后对76个台站所构成的2850个台站对的记录进行波形相似度计算.一台站记录到的P 波的波形相似度值为该台站与所有其他台站形成的台站对的波形相似度的平均值.这76个台站的P波波形相似度大于0.6的仅占2%.且这76个台站的P 波波形图也显示出了波形的低相似度，图 3把这76个台站的P波画在一张图上，并不能看到任何清晰的震相.这一结果正是由于玉树地震震源的走滑机制所造成.因而，我们只选取这76个台站中MCCC 计算的P 波波形相似度大于0.3 的35 个台站.由于35 个台站中有近一半位于澳大利亚，为避免由于某一方向台站过于密集所带来的台站响应函数的歧形，我们对这35个台站采用1000km×1000km 的地球表面块再次进行了重采样，图 2b为经过以上台站选取步骤后所用的29个台站和其台站响应函数.这29 个台站所组成的全球子台网显示最大能量位于震中位置，且旁瓣较少.除能得到一较好的近似Delta函数分布的台站响应函数外，这29个台站的P波波形相似度大于0.6的台站有20个，占总数的69%，见图 3b.
|
图 3 (A)76条未滤波P波记录（a)和MCCC波形相似度（c);B)29条未滤波P波记录（b)和MCCC波形相似度（d) Fig. 3 (A) 76unfiltedP records (a) and MCCCvalues (c) ; (B) 29unfiltedP records (b) and MCCCvalues (d) |
除了台站的选取，反投影的一个重要技术点在于准确地得到从震中附近的一点到远震距离上的地震台站的P波的传播时间.径向地球速度模型的精度大体上来说并不足以满足反投影的需要.这是由于三维的地球结构的变化会造成传播时间的异常.为了弥补三维速度结构对到时的影响，我们对所有台站记录到的P波的头10s的记录进行校准，允许每一台站的P波前后移动，从而使得所有台站的P波在同一时刻到达.我们假设这一P 波移动时间变化很小，这一P波到时校正值从而被用于震源区域的所有网格点.
在校正了P波到时后，另一需要加以关注的量是振幅.在测量P 波到时校正值时，我们还同时得到P波振幅的归一化参数.这一归一化参数可以有效地去除台站场地影响，波的几何扩散影响，仪器放大值不同的影响和地震波辐射方向的影响.
我们的反投影远震P 波研究过程的方法和Ishii等[3]运用于2004年苏门答腊9.3级地震的方法相似.反投影方法是在某一指定时间，通过对与某一可能的震源位置所对应的波形进行叠加来抵消噪音和传播路径中次生波的影响，从而突出从震源传出的信号.然后，把叠加所得到的能量投影到与之相对应的震源位置.在对可能震源区域所有可能位置都进行了能量反投影后，能够得到该时间的震源图像.之后，把这一处理过程运用到从震前到震后的一个连续时间段上，从而得到该地震全时间段的震源破裂过程.
在进行反投影时，我们考虑一均匀的四维空间(经度、纬度、深度、时间)围绕在震源和发震时刻周围.对每一网格点，我们都计算其与每一台站相对应的理论到时，然后对用带宽0.5~1.5 Hz 的Buterworth滤波器滤波后的波形运用4 阶方根叠加算法进行叠加.叠加公式如下：
|
(2) |
|
(3) |
其中，B(t)是最终的叠加值，B′(t)是叠加中间值，b(t)是第j条记录的振幅，M是波形总数.当N等于1时，为线性叠加.我们选用4 阶方根叠加算法(N=4)进行叠加是由于线性叠加虽然运算速度较快，但却不能很好地提高信号的能量并抑制噪音的干扰，而这两个目的在运用高阶方根叠加时能够很好的达到.在叠加后，叠加值的能量被反投影到与之相对应的网格点上.最大能量值所在的时间和空间区域为破裂的区域.叠加时采用2s窗长，1s滑动的滑动窗进行叠加.
4 结果与讨论通过对玉树地震运用反投影P 波研究其震源破裂特征，体现了在反投影研究中台站选取对结果的影响.由于全球台站分别不均，运用全球台网所得台站响应函数(ARF)的分布并不近于二维Delta函数，我们需对全球台站进行重新选取，选取的标准是全球子台网的台站响应函数接近于二维Delta函数分布.因而在对震源破裂尺度在几十至上百公里的地震进行反投影研究时，台站的选取是影响反投影结果的一个最重要的因素.
在对玉树地震震源破裂过程进行研究的过程中，我们通过台网响应函数和MCCC 对从IRIS 下载的381条宽频带垂直向记录进行筛选后，我们运用29个台站构成全球子台网记录对玉树地震的震源破裂过程进行研究.这29个台站有较高的P波波形相关性，且这29个台站组成的全球子台网的台网响应函数在震中位置能量为最大，旁瓣较少，能量分布接近于二维Delta函数分布.这一台站选取保证了反投影结果的可靠性.
我们运用4阶方根叠加法对总长100s，发震时刻前20s到震后80s的时间段的波形进行反投影，得出能量最大值随时间的变化曲线(图 4)，图 4 中的时间轴上，0代表的是发震时刻.能量值在0时刻前就有逐渐爬升的现象是由于我们采用的是滑动时间窗进行波形叠加.滑动时间窗会使时间窗中心点在发震时刻前的窗包含有后面的能量，从而导致能量的明显抬升发生在0 时刻前.能量最大值随时间的变化曲线显示在震后20s 后能量值稳定地在0.15附近波动，虽然此时能量没有恢复到震前的低能量水平，但这是一整体的能量值的抬高，受到大震尾波能量的干扰.Kærna[17]对震后全球台网监控能量的研究指出，一次大地震发生时全球的噪音水平被明显地提高了，这一背景噪音的提高先是表现在包含有两个主要能量释放，一个能量释放点是在震中位置附近，在震后约6s的时间；另一个是位于震中位置西南40km 处，发生在震后约12s的时间，此次能量释放也是整个破裂过程中能量最大的，见图 5.张勇等[18]对玉树地震震源破裂过程的反演结果同样给出约为20s长的破裂时间、两次主要的能量释放、且第二次能量释放为最大能量释放的研究结果.在和张勇等[18]的反演结果进行对比时，两种研究震源破裂过程的地震学方法得到的玉树地震震源破裂的特征基本一致.但也存在差异，那就是本研究得到的两次主要能量释放的时间为震后6s和12s，而反演结果的两次主要能量释放时间为震后2.7s和7.7s.造成这一能量释放时间差异的原因有可能是两种方法所研究的频率范围不同.Xu 等[6]对2008年汶川地震的研究表明，破裂过程的一些细节和所研究的频率范围有直接的联系.研究的频率不同时，所得到的震源破裂过程的主要特征不会受到影响，但在某些细节上会有所不同.两种方法所研究的频率范围的差异，有可能造成最大能量释放时间的差异.但并不改变破裂过程的主要特征：两个主要能量释放过程.
|
图 4 归一化的能量时间变化 Fig. 4 Normalized beam power as a function of time |
|
图 5 震源破裂过程截图.插人曲线为总破裂能量时间函数.在每一截图上红色线为相应时间，黑色星为玉树地震震中 Fig. 5 Snapshots of the rupture. In each case the red line mdicates the corresponding time，black star is the Yushu earthquake |
反投影远震P 波得到的玉树地震能量释放积累图见图 6.玉树地震的能量呈现北东-南西向分布，主要沿玉树断裂发展，破裂长度约为60km.两个主要的能量释放点，一个在震中附近，一个位于震中西南，靠近结古镇的位置.且靠近结古镇附近的能量释放是最大的，这也是结古镇的损失是本次地震中最重的地区的原因.这一结果和有限源反演[18]及InSAR 结果(http://web2.ges.gla.ac.uk/~zhli/Qinghai-Yushu-EQ.htm)一致.图 6所显示研究区域的玉树断裂分成了南北两条，反投影研究所得到的能量释放较为靠近玉树南断裂而不是北面的玉树断裂；张勇等[18]得到的结果更靠近北面的玉树断裂，造成这一差异的原因可能有两个：一是在研究时所用震中经纬度不同，我们用的是USGS提供的结果：北纬33.156°，东经96.529°;而张勇等[18]的反演研究所用的震中经纬度为：北纬33.1°，东经96.7°.反演研究所用的经纬度位置更靠近北面的玉树断裂；另一个原因可能是频率，如前所述，研究时所用频率范围不同，所得到的破裂过程的一些细节会有所不同.
由反投影结果得到的玉树震源破裂尺度为60km，震源破裂时间为20s，我们可以得到本次玉树地震震源破裂的传播速度，约为3km/s.
Faenza[19]指出由于台站分布的影响，由软件自动生成的震动图并不能充分地反映地面运动，因而反投影方法得到的破裂模型能成为震动图的一个补充，为震后的应急救援提供有用信息.我们的反投影远震P波和有限源反演结果[18]所给出的最大能量释放区域为震中西南，靠近结古镇，和震后实际震害调查的结果相一致，这一事实证明了地震学方法所得到的震源破裂模型是震动图的有效补充.
反投影P波技术多被运用于震级大于7.5的大震，其有着上百公里的破裂尺度.然而D′Amico等[2]运用本方法研究了2009年4月6日的意大利Mw6.3级地震.这一研究表明反投影P 波技术的应用范围除了有着上百公里破裂尺度的大震外，有着重要破裂特征的、破裂尺度在几十公里范围的中强地震同样可以运用反投影技术快速地得到震源的破裂过程.对玉树地震的研究也同样体现出反投影方法可以较好地得到破裂尺度在几十公里范围的中强地震的破裂过程.虽然本方法不能像反演方法那样得到详细的破裂图像，但由于反投影方法所需动力学参数限制少，且运算快，是一种能为震后快速应急救援工作服务的地震学方法.
5 结 论通过反投影远震P 波，我们得到2010 年4 月14日的玉树地震震源破裂时间为20s，破裂长度为60km，破裂主要沿北西—南东向的玉树断裂发展，破裂的传播速度为3km/s.玉树地震的能量释放主要是在两个时间和空间点上，一个能量释放的时间点是震后6s，能量释放的空间位置为震中附近；另一个能量释放的时间点是震后12s，释放的空间位置为震中东南方向约40多公里处，靠近此次玉树地震破坏最大的结古镇，这次能量释放是整个玉树地震震源破裂过程中最大能量释放的时空点.
反投影远震P 波是一对破裂过程中的动力学参数限制最少的研究震源破裂过程的方法.该方法和有限源反演相互补充，是一能在震后较短时间内得到震源破裂过程的地震学方法.
致谢我们感谢IRISDMC 提供的波形数据.我们用GMT(Wessel and Smith，1991)[20]作图.感谢两位审稿人对本文提出的非常有益的修改意见.
| [1] | Lay T, Ammon C H, Hutko A R, et al. Effects of kinematic constraints on teleseismic finite-source rupture inversions: Great Peruvian earthquakes of 23 June 2001 and 15 August 2007. Bull Seismol Soc Am , 2010, 100. DOI:10.1785/0120090274 |
| [2] | D'Amico S, Koper K D, Herrmann R B, et al. Imaging the rupture of the Mw6.3 April 6, 2009 L'Aquila, Italy earthquake using back-projection of teleseismic P-waves. Geophys Res Lett , 2010, 37: L03301. DOI:10.1029/2009GL042156 |
| [3] | Ishii M, Shearer P, Houston H, et al. Extent, duration and speed of the 2004 Sumatra-Andamann earthquake imaged by the Hi-Net array. Nature , 2005, 435: 933-936. |
| [4] | Krüger F, Ohrnberger M. racking the rupture of the Mw=9.3 Sumatra earthquake over 1150 km at teleseismic distance. Nature , 2005, 435: 937-939. |
| [5] | Walker K T, Ishii M, Shearer P M. Rupture details of the 28 March 2005 Sumatra Mw8.6 earthquake imaged with teleseismic P waves. Geophys Res Lett , 2005, 32: L24303. DOI:10.1029/2005GL024395 |
| [6] | Xu Y, Koper K D, Sufri O, et al. Rupture imaging of the Mw7.9 May 2008 Wenchuan earthquake from back projection of teleseismic P waves. Geochem Geophys Geosyst , 2009, 10: Q04006. DOI:10.1029/2008GC002335 |
| [7] | Dziewonski A M, Anderson D L. Preliminary reference Earth model (PREM). Phys Earth Planet Int , 1981, 25: 297-356. DOI:10.1016/0031-9201(81)90046-7 |
| [8] | Kennett B L N, Engdahl E R, Buland R. Constraints on seismic velocities in the earth from travel times. Geophys J Int , 1995, 122: 108-124. DOI:10.1111/gji.1995.122.issue-1 |
| [9] | Kennett B L N, Engdahl E R. Travel times for global earthquake location and phase identification. Geophys J Int , 1991, 122: 429-465. |
| [10] | Rowley D B. Age of initiation of collision between India and Asia: A review of stratigraphic data. Earth Planet Sci Lett , 1996, 145. DOI:10.1016/S0012-821X(96)-201-4 |
| [11] | England P, Housemann G. Finite strain calculations of continental deformation II: Comparison with the India-Asia collision zone. J Geophys Res , 1986, 91: 3664-3676. DOI:10.1029/JB091iB03p03664 |
| [12] | Royden L H, Burchfiel B C, King R W, et al. Surface deformation and lower crustal flow in eastern Tibet. Science , 1996, 276: 788-790. |
| [13] | Tapponnier P, Peltzer G, Le Dain, et al. Propagation extrusion tectonics in Asia: New insights from simple experiments with plasticine. Geology , 1982, 10: 611-616. DOI:10.1130/0091-7613(1982)10<611:PETIAN>2.0.CO;2 |
| [14] | Yin A, Harrison T M. Geologic evolution of the Himalayan-Tibetan Orogen. Annual Review of Earth and Planetary Sciences , 2000, 28. DOI:10.1146/annuerv.earth.28.1.211 |
| [15] | Wang S, Fan C, Wang G, et al. Late Cenozoic deformation along the northwestern continuation of the Xianshuihe fault systym, Eastern Tibetan Plateau. GSA Bulletin , 2008, 120. DOI:10.1130/B25833.1 |
| [16] | VanDecar J C, Crosson R. Determination of teleseismic relative arrival times using multi-channel cross-correlation and least squares. Bull Seismol Soc Am , 1990, 80: 150-159. |
| [17] | Krna T, Ringdal F. Seismic threshold monitoring for continuous assessment of global detection capability. Bull Seismo Soc Am , 1999, 89: 946-959. |
| [18] | 张勇, 许力生, 陈运泰. 2010年4月14日青海玉树地震破裂过程快速反演. 地震学报 , 2010, 32(3): 361–365. Zhang Y, Xu L S, Chen Y T. Fast inversion of rupture process for 14 April 2010 Yushu, Qinghai earthquake. Acta Seismologica Sinica (in Chinese) , 2010, 32(3): 361-365. |
| [19] | Faenza L, Lauciani V, Michelini A. ShakeMaps of the L'Aquila main shock. Eos Trans AGU , 2009, 89(53): Fall Meet Suppl, Abstract U23B-0036. |
| [20] | Wessel P, Smith W. Free software helps map and display data. Eos Trans AGU , 1991, 72: 441. |