2015年4月25日6时11分26.27 s(国际时)于 尼泊尔发生MS8.1(Mw7.8)地震,震中位于84.708°E,28.147°E,震源深度15 km. 大约半小时后,在主震附近发生了一次MS7.0(Mw6.6)的强余震(A1),一天后,即2015年4月26日7时9分11.01 s(国际时),在距主震~130 km处又发生了一次MS7.1(Mw6.7)的强余震(A2)(图 1).
这次MS8.1地震发生在印度板块与欧亚板块的边界带,位于青藏高原的南缘. 根据史料记载,这里发生过多次灾难性地震(Bilham et al.,2001),形成一系列低角度的逆冲断裂(Yin and Harrison,2000; Bilham et al.,2001). 距今最近的地震当属1950年的阿萨姆地震(Chen and Molnar,1977; Molnar,1990). 2015年4月25日发生的MS8.1地震及其强余震当属这类低倾角逆冲型事件,这已得到震源机制初步结果的证实(表 1,表 2,表 3).
地震发生不久,除震源机制外,国内外地震研究机构还很快发布了震源破裂过程的结果(张勇等,2015;USGS,2015; 王卫民等,2015). 这些结果之间有相同之处,均表明位错发生在起始破裂点的东侧,似乎表明为一单侧破裂,地震过程持续约100 s,可以分为两个阶段,即两次事件,第一次事件发生在前50 s,是主要的能量释放区间,第二次事件发生在后50 s,但释放的能量较小. 不过,更引人注意的是彼此之间不可忽视的差异.
震源破裂过程的反演本身是一项比较复杂的工作,初步结果之间存在差异在所难免. 况且,作为震后快速响应,使用的技术比较单一,使用的资料也比较有限,而且不同的人采用的技术和资料也不尽相同. 因此,对这次地震的震源过程进行多方面更深入的研究是非常必要的.
我们获取的地震记录中主要包含两种效应,一是地震的震源效应,二是路径的传播效应. 我们研究震源破裂过程,首要任务是去除路径的传播效应. 地震学的发展使我们能够通过理论计算得到路径效应(Kennett,1983; Wang,1999),即理论格林函数,但其精确和可靠程度很大程度上依赖于介质模型. 然而,介质模型和实际的介质之间的差别始终是存在的,因此理论格林函数是对路径效应的近似. 与理论格林函数相比,经验格林函数不依赖于介质模型,可以在相当程度上代表路径传播效应. 因此,利用经验格林函数可以更好地分离出震源效应(Hartzell et al.,1978; Mueller et al.,1985).
根据已有的震源机制解(表 1,表 2和表 3),余震A1和A2与主震都具有非常类似的震源机制. 余震A1的震源机制不但与主震相似,而且震源位置也非常接近,但是其发震时刻距主震太近,以致其信号在很多台站都无法与主震信号区分开来. 余震A2的震源机制与主震也非常接近,震源位置略远(二者相距约130 km),但仍在主震的震源区. 且考虑到我们所用的台站的震中距和信号的波长,余震A2仍不失为理想的经验格林函数事件.
已有的震源破裂过程多借助于理论格林函数获得(张勇等,2015;USGS,2015; 王卫民等,2015),所以,本研究以余震A2为经验格林函数事件,首先利用PLD(Projected L and weber Deconvolution)(Bertero et al.,1997;Piana and Bertero,1997;张勇等,2009)技术从遍布全球和全国的宽频带台站(图 1)记录的波形数据中提取主震的P波视震源时间函数和Rayleigh波视震源时间函数,分析主震的破裂方向性特征;然后借助于震源破裂过程的视震源时间函数反演技术(Chen and Xu,2000; Xu et al.,2002; 张勇,2008),利用提取的P波视震源时间函数和Rayleigh波视震源时间函数联合反演主震的时空破裂过程.
2 数据为了利用经验格林函数技术分析主震震源过程的复杂性,我们利用分布在全球(IRIS)和全国(郑秀芬等,2009)的宽频带台站记录的主震和余震A2的波形记录,并初步选择震中距10°~90°范围的垂直向波形记录. 为了使所用台站分布在空间上较为均匀,我们对选取的资料按照震中距和方位角分布进行了重新筛选,使震中距间隔和方位角间隔均为10°,最终用于本研究的台站分布如图 1b所示.
3 视震源时间函数 3.1 P波视震源时间函数首先我们将主震和余震A2的垂直向P波记录的采样率降至10sps,并采用0.01~0.05 Hz的带通滤波去除低频和高频噪声,然后分别截取主震和余震P波初动前10 s至后110 s信号,最后利用PLD技术提取视震源时间函数. 利用PLD技术提取各台站的视震源时间函数(即P波视震源时间函数)后,仍需进行一些后处理,例如,根据合成地震图(视震源时间函数与经验格林函数的褶积)和观测地震图的匹配情况去掉一部分匹配较差的台站的视震源时间函数,而且还要根据视震源时间函数的方位依赖性特征再去掉一些被认为奇异的视震源时间函数.造成这种情况的原因很可能是主震记录和/或经验格林函数记录受到了“瞬时”干扰. 图 2a展示了最后入选的所有P波视震源时间函数.
把从不同台站的P波视震源时间函数按照台站相对于震中的方位进行排列(图 2a),我们可以清楚地发现这些视震源时间函数的形状在随方位的变化而变化. 这种方位依赖特征表明,这次地震的震源具有相当的尺度且破裂具有明显的方向性(Lay and Wallace,1995; 许力生等,2014). 破裂的优势方向应该在100°左右,因为在此周围的震源时间函数的有效持续时间最短,且起始时间最早;相反,300°左右应该是破裂的背向,因为在此周围的震源时间函数的有效持续时间最长,且起始时间最晚. 另外,这次地震的主要能量应该释放在前50 s. 50 s后虽然有能量释放,但已相当少. 同时,从180°与270°之间的震源时间函数可以看出,即使前50 s,震源过程仍具有一定的复杂性,至少有两次大小相当的子事件.
利用经验格林函数技术提取的视震源时间函数只包含地震矩相对大小的信息,而没有地震矩绝对 大小的信息. 为了讨论问题的方便,我们采用USGS测定的标量矩 5.449×1020 N·m(http://earthquake. usgs.gov/earthquakes/eventpage/us20002926 scientific_ tensor:us_us_20002926_mwc〖2015-05-09〗)对各台站提取的视震源时间函数进行标定,让这些视震源 时间函数包含地震矩的绝对大小,并将其投影到极坐标中(图 3a).
视震源时间函数的极坐标更直观且更清楚地展示了视震源时间函数的方位依赖性. 从图 3a可以看出,地震矩释放的有效时间窗(高亮度带)大体呈椭圆形,在东南方向,高亮度带较窄,且离中心圆较近,这是破裂的优势方向;在西北方向,高亮度带相对较宽,且距中心圆较远,这是破裂的背向. 同时,我们还注意到,在东北方向和西北方向有一些高亮度点,是能量汇聚的方向,表明破裂有朝这些方向扩展的迹象.
为了定量估计地震破裂的优势方向,我们测量了所有P波视震源时间函数的峰值时间TP,并利用正弦函数对这些时间点进行拟合(图 3b). 根据拟合的结果,TP最小值位于142°,最大值位于334°. 这两个值分别对应于破裂的优势方向和优势方向相反的方向. 注意,142°的相反方向为322°,而不是334°,这是由于实际的破裂模型不是一个简单的线源模型的缘故,是实际模型的复杂性的表现.
3.2 Rayleigh波视震源时间函数与P波视震源时间函数相比,面波视震源时间函数具有更好的时间分辨能力(Lay and Wallace,1995; 许力生等,2014),因此,面波视震源时间函数的方位依赖性更强,更能够反映出破裂的优势方向. 为此,我们从Rayleigh波中也提取了视震源时间函数(即Rayleigh波视震源时间函数).
在提取Rayleigh波视震源时间函数时,我们采用了与处理P波资料类似的流程. 但考虑到面波频率较低,使用了0.005~0.0333 Hz的带通滤波,且采用2.5~4.2 km·s-1的群速度窗截取Rayleigh波资料.
图 2b展示了最终优选的Rayleigh波视震源时间函数. 正如所料,我们可以更清楚地看到这些视震源时间函数的形状随方位的变化. 再次表明,这次地震的震源具有相当的尺度且破裂具有明显的方向性.
不过,Rayleigh波视震源时间函数表明,破裂的优势方向在70°~80°之间,而不是在100°左右;破裂的背向在280°左右,而不是在300°左右. 这种差异是由于P波包含周期相对较短的破裂信息而Rayleigh波包含周期相对较长的破裂信息所致.
图 4a是Rayleigh波视震源时间函数在极坐标中的展示. 与P波视震源时间函数的极坐标展示相比,Rayleigh波视震源时间函数反映的方位依赖性更强烈、更清楚. 高亮度带形成更加明显但不很规则的椭圆形,在东北方向和东南方向,高亮度带较窄,且离中心圆较近,这是破裂的优势方向;在西北方向,高亮度带相对较宽,且距中心圆较远,这是破裂的背向.同时,在东北方向、西北方向和西南方向有一些高亮度点,表明这些方向也是能量汇聚的方向,这都是破裂过程复杂性的表现.
图 4b展示了所有Rayleigh波视震源时间函数的峰值时间TP测量值以及利用正弦函数对这些点的拟合曲线. 需要说明的是,在这里我们进行了两次拟合,第一次拟合考虑了所有的测量值,拟合结果如红色曲线所示,TP最小值位于71°,最大值位于285°. 第二次拟合没有考虑紫色的测量值,拟合的结果如蓝色曲线所示,TP最小值位于108°,最大值位于288°. 图 4b中紫色的测量值来自图 4a中东南方向那块特殊的区域,这块区域在图 3a中看不到,因此这部分能量还不能确认. 或者说,虽然是震源过程复杂性的表现,但不属于共性特征。 考虑到这种特殊情况,我们认为第二种拟合结果更可取,即破裂的优势方向在108°,而背向在288°.
4 震源破裂过程P波视震源时间函数包含相对高频的震源破裂信息,而Rayleigh波视震源时间函数包含相对低频的震源破裂信息. 为了获取主震比较完整的震源破裂过程图像,我们联合反演P波视震源时间函数和Rayleigh波视震源时间函数.
经过如前所述的优选后,有些台站既有P波视震源时间函数,又有Rayleigh波视震源时间函数,而有些则只有P波视震源时间函数或者Rayleigh波视震源时间函数,而且这些台站在空间上很不均匀. 太过不均匀的台站分布必然对反演结果造成影响,所以在反演之前,我们对台站进行了重新筛选,尽可能使得用于反演的视震源时间函数资料在空间上趋于均匀. 同时,为了在确保反演结果不失时间分辨能力的情况下减少计算量,我们把用于反演的视震源时间函数的采样率降至1 sps.
根据USGS发布的结果,主震的微观震中在28.147°N,84.708°E,震源深度为15 km;主震的震源机制如表 1所示. 我们以走向295°,倾角11°的节面为发震断层面,以断层面与地表的交线为断层上边界,自地表沿断层倾向180 km处为断层的下边界. 我们以USGS确定的震源位置为起始破裂点. 由起始破裂点沿断层面向西北80 km处作为断层的西北边界,由起始破裂点沿断层面向东南150 km处作为断层的东南边界. 将这个矩形区域分割成23×18个子断层,使子断层成为10×10 km的正方形.
为了稳定反演结果,我们不但引入了空间光滑约束,还引入了时间光滑约束(Yagi et al.,2004; 张勇,2008),同时,还引入了标量地震矩最小的约束以压制过低频的噪声(张勇,2008). 尽管引入了上述约束,但方程系统依旧是一个线性系统,因此,我们采用一直以来使用的共轭梯度法(Chen and Xu,2000; Xu et al.,2002; 张勇,2008)求解这个方程系统.
需要说明的是,为了保持线性的反演系统,我们不得不假设震源破裂的最大破裂速度,也不得不设定每个子断层的最大滑动时间. 经过多次尝试,3 km·s-1的最大破裂速度和40 s的子断层最大滑动时间不但能够满足方程系统而且能够最好地解释所有观测数据.
图 5和图 6展示了反演的结果. 破裂具有明显 的单侧破裂特征,主要的位错分布在起始破裂点左 下方,破裂区域大体呈三角形,水平方向长达~100 km,沿断层面向下延伸~80 km. 最大位错达~5.8 m,平均位错达~2.3 m. 图 6以快照的形式展示了破裂的传播过程. 可以看出,破裂过程是一个从起始点开始逐渐向左下方传播的过程. 换句话说,在水平方向上,从右向左传播;在垂直方向上,从上向下 传播.图 5b展示了与震源破裂过程对应的地震矩率随时间的变化过程,最大矩率为~2.5×1019 Nm·s-1. 需要说明的是,在图 5的左上角有一较小的破裂区,这一破裂区尺度小滑动弱,且出现在地震过程即将结束的时段(图 6),所以很可能是噪声所致.
图 7展示了观测视震源时间函数与合成视震源函数的对比,平均相关系数达0.88,这表明反演得到的同震位错模型能够很好地解释观测资料.
主震和余震的空间关系是我们关注的一个重要问题. 为此,我们计算了每个子断层释放的标量地震矩,并计算了相应的矩震级(Lay and Wallace,1995). 把每个子断层作为一次地震事件展示于图 8. 可以看出,主震事件和余震之间是空间上互补的 关系,而且余震大多发生在主震破裂的尾端. 同时,从图 8还可以更清楚地看出,破裂从西向东从浅至深的单侧破裂特征.
尼泊尔MS8.1地震和MS7.1余震震源机制的相似性允许我们借助于经验格林函数技术从一个新视角认识主震震源破裂过程的复杂性. 与理论格林函数相比,经验格林函数能更好地描述路径的传播效应. 面波视震源时间函数比体波视震源时间函数的方位依赖性更强. 体波视震源时间函数携带相对高频的破裂信息,而面波视震源时间函数携带相对低频的破裂信息. 正因为这些特点,我们选择余震作为经验格林函数事件对主震的震源过程进行反演分析.
为了尽可能准确认识尼泊尔MS8.1地震的震源复杂性,我们尽最大努力收集了国内外的宽频带地震记录. 经过主震和余震记录的配对筛选,视震源时间函数方位依赖性筛选,最终获得127条P波视震源时间函数和140条Rayleigh波视震源时间函数. 这些视震源时间函数最大程度地揭示了主震震源断层的有限性和破裂的方向性.
通过联合反演P波视震源函数和Rayleigh波视震源时间函数获得的震源破裂过程表明,尼泊尔MS8.1地震是一次单侧破裂事件. 破裂从起始点开始,沿断层面向东南方向扩展~100 km,与此同时,破裂沿断层面向深部扩展~80 km,破裂面呈三角形状,最大位错约5.8 m. 跟已有的结果相比(张勇等,2015;USGS,2015;王卫民等,2015),相同之处仅在破裂的方向性上,均为向东南方向扩展的单侧破裂;不同之处也十分明显,破裂面的形状各不相同,位错量也各不相同.
已有的结果表明(张勇等,2015;USGS,2015; 王卫民等,2015),尼泊尔MS8.1地震持续了约100 s,而且有两次事件,第一次在前50 s,是主要的事件,第二次在后50 s,是次要的事件. 我们的反演结果则清楚地显示,这次地震只有一次事件,发生在前50 s. 至于已有结果中出现的第二次事件,我们认为可能是青藏高原特殊的介质结构引起的特殊的路径效应的反映,也可能是观测资料不完善造成的结果. 为了确认尼泊尔MS8.1地震的破裂历史,将我们反演破裂过程时得到的震源时间函数、平均的P波视震源时间函数和Rayleigh波视震源时间函数以及二者平均后的视震源时间函数展示于图 9. 可以看出,在P波视震源时间函数和Rayleigh波视震源时间函数以及二者平均后的视震源时间函数上都可以看出50 s之后的第二次事件,但这是方位依赖性引起的虚假现象. 因此,我们可以确认,尼泊尔MS8.1地震的时间历史相对比较简单,这次地震是一次“一气呵成”的事件. 至于它的持续时间,这里的反演结果显示为40 s. 但如果考虑经验格林函数事件的持续时间(10 s左右),那么这次地震的持续时间应为50 s左右.
根据P波视震源时间函数的方位依赖性,破裂的优势方向在142°;根据Rayleigh波视震源时间函数的方位依赖性,破裂的优势方向在108°. 我们知道,体波视震源时间函数包含相对高频的破裂信息,而面波视震源时间函数包含相对低频的破裂信息,所以不难理解142°的破裂方向主要取决于浅部破裂;而108°的破裂方向主要取决于深部破裂. 这一特征与断层破裂向东南向深部扩展的特征一致. 如果考虑断层综合的破裂方向,我们不妨取二者的算术平均,即125°.
综上所述,尼泊尔MS8.1地震的震源几乎 是纯粹的单侧破裂,从破裂起始点开始,沿断层面向东南 方向扩展~100 km,同时沿断层面向深部扩展~80 km,形成破裂的优势方向125°. 地震的能量释放历史总体比较简单,属于一次非间断性扩展的事件,持续时间约50 s,形成最大位错 ~5.8 m.
致谢 本研究使用的波形资料来源于IRIS数据中心和中国地震局地球物理研究所“国家数字测震台网数据备份中心”.
[1] | Bertero M, Bindi D, Boccacci P, et al. 1997. Application of the projected Landweber method to the estimation of the source time function in seismology. Inverse Problems, 13: 465-486. |
[2] | Bilham R, Gaur V K, Molnar P. 2001. Himalayan Seismic Hazard. Science, 293: 1442-1444. |
[3] | Chen W P, Molnar P. 1977. Seismic moments of major earthquakes and the average rate of slip in Central Asia. J. Geophys. Res., 82(20), 2945-2969. |
[4] | Chen Y T, Xu L S. 2000. A time-domain inversion technique for the tempo-spatial distribution of slip on a finite fault plane with applications to recent large earthquakes in the Tibetan Plateau. Geophys. J. Int., 143: 407-416. |
[5] | Hartzell S H. 1978. Earthquake aftershocks as Green's functions. Geophys. Res. Lett., 5(1): 1-4. |
[6] | Kennett B L N. 1983. Seismic wave propagation in stratified media.Cambridge University Press. |
[7] | Lay T, Wallace T C. 1995.Modern global Seismology. Academic Press, INC. |
[8] | Mueller C S. 1985. Source pulse enhancement by deconvolution of an empirical Green's function. Geophys. Res. Lett., 12(1): 33-36. |
[9] | Molnar P. 1990. A review of the Seismicity and the Rates of active underthrusting and deformatation at the Himalaya. J. Himalayan Geol., 1, 131-154. |
[10] | Piana M, Bertero M. 1997. Projected Landweber method and preconditioning. Inverse Problems, 13: 441-463. |
[11] | USGS. 2015. http://earthquake.usgs.gov/earthquakes/eventpage/us20002926 scientific_fifinitefau[2015-05-09]. |
[12] | Wang R J. 1999. A simple orthonormalization method for stable and efficient computation of Green's functions. Bull. Seismol. Soc. Am., 89(3): 733-7410. |
[13] | Wang W M, Hao J L and Yao Z X. 2015. http://www.itpcas.ac.cn/xwzx/zhxw/201504/t20150426_4344080.html[2015-05-09]. |
[14] | Xu L S, Chen Y T, Teng T L, et al. 2002. Temporal-spatial rupture process of the 1999 Chi-Chi Earthquake from IRIS and GEOSCOPE long-period waveform data using aftershocks as Empirical Green's Functions. Bull. Seismol. Soc. Am., 92(8): 3210-3228. |
[15] | Xu L S, Zhang X, Yan C, et al. 2014. Analysis of the Love waves for the source complexity of the Ludian MS6.5 earthquake. Chinese J. Geophys.(in Chinese), 57(9): 3006-3017. |
[16] | Yagi Y, Mikumo T, Pacheco J, et al. 2004. Source rupture process of the Tecomán, Colima, Mexico Earthquake of 22 January 2003, determined by joint inversion of teleseismic body-wave and near-source data. Bull. ,Seismol. Soc. Am. 94(5): 1795-1807. |
[17] | Yin A, Harrison T M. 2000. Cenozoic evolution of the Himalayan-Tibetan orogeny. Annu. Rev. Earth, Planet Sci., 28: 211-280. |
[18] | Zhang Y. 2008. Study on the Inversion Methods of source rupture process [Ph. D. thesis]. Beijing: Peking University. |
[19] | Zhang Y, Xu L S and Chen Y T. 2009.PLD method for retrieving apparent source time function and its application to the 2005 Kashmir Mw7.6 earthquake.Chinese J. Geophys.(in Chinese), 52(3): 672-680. |
[20] | Zhang Y, Xu L S, Chen Y T. 2015. http://www.cea-igp.ac.cn/tpxw/272110.shtml[2015-05-09]. |
[21] | Zheng X F, Ouyang B, Zhang D N, et al. 2009. Technical system construction of Data Backup Centre for China Seismograph Network and the data support to researches on the Wenchuan earthquake. Chinese J. Geophys.(in Chinese), 52(5): 1412-1417. |
[22] | 王卫民, 郝金来, 姚振兴. 2015. http://www.itpcas.ac.cn/xwzx/zhxw/201504/t20150426_4344080.html[2015-05-09]. |
[23] | 许力生, 张旭, 严川等. 2014. 基于勒夫波的鲁甸MS6.5地震震源复杂性分析. 地球物理学报, 57(9): 3006-3017. |
[24] | 张勇. 2008. 震源破裂过程反演方法研究[博士论文]. 北京: 北京大学. |
[25] | 张勇, 许力生, 陈运泰. 2009. 提取视震源时间函数的PLD方法及其对2005年克什米尔Mw7.6地震的应用. 地球物理学报, 52(3): 672-680. |
[26] | 张勇, 许力生, 陈运泰. 2015. http://www.cea-igp.ac.cn/tpxw/272110.shtml[2015-05-09]. |
[27] | 郑秀芬, 欧阳飚, 张东宁等. 2009. "国家数字测震台网数据备份中心"技术系统建设及其对汶川大地震研究的数据支撑.地球物理学报, 52(5):1412-1417. |