喜马拉雅是全球大陆地震最活跃的地区之一，仅过去的一个世纪就曾发生多起8级左右的特大地震(Bilham et al.，2001; Kumar et al.，2010)，进入新世纪来，曾出现过2005年克什米尔M_W7.6地震(Pathier et al.，2006)，而2015年4月25日尼泊尔中东部郎唐地区(Lamjung)的M_S8.1级地震则是自1950年察隅地震后，喜马拉雅地区震级又一次达到8级的地震.据全球台网观测，该次地震发生于UTC时间06:11:26.2，破裂从博卡拉下方开始，在此后约80 s时间，破裂(图 1)向东偏南方向单侧扩展，延伸到加德满都以东，并引发60余次(M_W>4)余震，其中两次为M_W6.7、6.8强震.尽管极震区主要位于喜马拉雅南缘的尼泊尔境内，但地震也波及到喜马拉雅北缘的西藏定日、聂拉木、吉隆三地，地震当日境内曾诱发出一次M_S5.9地震.

图 1

Fig. 1

图 1 尼泊尔MW7.9地震的滑动分布及GPS揭示的永久同震位移、动态形变波 GPS水平同震位移以蓝、绿、红边框箭头及95%的置信区间表示的不同尺度的形变量组成，垂直位移以点位处的不同颜色柱体表示，图左上方的插入图像为图中白色方框区域，标明了此处的GPS位移矢量；动态形变波形以归一化尺度绘制于点位四周，波形附近的同色数字表示最大振幅值，单位为m；黑色粗线为发震断层；灰色细线为绘图区域的主要断裂；黑色方框为本文断层模型边界在地表处投影；红色及粉色圆圈为余震分布；红色五角星为震源位置；灰色圆点表示城市或城镇位置；图右上角插入的地球模型上的品红色圆点为联合反演的远震地震台站分布. Fig. 1 Slip distribution and surface deformation including permanent offset and kinematic wave inferred from GPS measurements The horizontal displacements are depicted using differently scaled arrows with 95% confidence ellipses, as well as the vertical displacements are inferred by the rectangles fitted with different color. The upper left insert shows the GPS displacements in the white rectangle in detail. The locations of the XZAR, XZRK, XZZF, LHAZ, and XZZB station are surrounding the corresponding kinematic deformation waves, respectively, and the max amplitudes are annotated nearby with same colors. Black heavy line indicates the causative fault. Gray thin line shows the main faults in mapping area. Black solidline box outlines the profile of the fault model. Red and pink circles correspond to the aftershocks with magnitude greater than MW6.6 and smaller than MW6.0, respectively. Red star shows the epicenter. Gray dots with white circle infer the locations of cities and towns. Magenta dots on the earth model inserted on the right upper depict the distribution of teleseismic waves used in this paper.

喜马拉雅地区处于大陆板块边界，其高海拔、陡峭地形和高山地貌是印度-欧亚板块长时间大幅度汇聚缩短、地壳增厚和山体隆升的结果(Harrison et al.，1992; Avouac et al.，2003).该地区近来人口快速增长，但严峻的现实是一旦大震降临，灾害效应也急剧扩大(Bilham et al.，2001)，伤亡往往格外惨烈，2005年克什米尔地震即是显证.本次地震也不例外，截至5月5日已导致7300余人死亡，大量房屋毁损，震区数百万民众不同程度受灾.因此对该地区大震活动性的研究需求十分迫切，而此次地震对研究板缘地带的大震成因机制提供了一次难得机遇.我们利用“中国大陆构造环境监测网络”(简称“陆态网络”)西藏境内GPS基准站和尼泊尔境内加州理工学院布设的GPS连续观测台网(Ader et al.，2012)，结合全球地震台网资料反演本次地震的破裂过程，特别是通过近场强地面运动GPS观测展示大陆板缘地带一次逆冲型特大地震的震源时空分布.

地震发生后，国内外研究机构相继给出了基于远震地震波资料反演的震源模型(张勇，2015； 王卫民，2015; USGS，2015；Wei Shengji，2015；Yagi and Okuwaki，2015)，这些模型具有大致相同的破裂特征，但模型之间的差异也十分明显.张勇等(2015)在此基础上，利用境外两个GPS点的同震形变作为约束，初步联合反演了此次地震的破裂过程，对破裂模型进行了一定程度的修订.以往的实践表明，近场的强地面观测可在一定程度弥补由于远震纪录对破裂细节不敏感的缺陷，恰巧本次地震我国境内分布了一定数量的高频GPS连续观测站，本文利用连续GPS揭示的近场强地面运动、永久同震形变与全球均匀分布的远震P波，联合反演了此次地震的破裂过程，以优化现有的破裂模型细节特征.

本文所用14个GPS台站的地理分布如图 1所示，其中境内的5个GPS站距离震中相对较远，较近的3个台站(XZZB、XZAR、XZZF)分别位于震中的西北、东北和正东方位，震中距大约170~300 km.除位于拉萨的GPS站外，境内的4个站全部为基岩观测墩(李强等，2012)，采用Trimble接收机和扼流圈天线，其已稳定运行3年有余.这4个站以0.02、1及30 s采样三种方式记录数据，其中1 s高频采样数据每隔1 h上传到北京的“陆态网络”数据中心.境外9个连续观测GPS站均距震中较近，最近的CHLM、KKN4、NAST位于极震区内，震中距分别为60 km、70 km和80 km.境外GPS站公布了15 s采样数据，上传到美国UNAVCO的数据库(http://www.unavco.org)供自由下载.本文所用远震数据由国际地震学联合会(IRIS)地震数据中心提供，我们选取了全球地震台网中震中距30°~90°范围、信噪比较高且方位角均匀分布的43个台站的垂向P波纪录(图 1).数据的时间窗口为110 s，各窗口起始时刻设定为P波达到前10 s.

图 1中12个GPS站30 s采样数据被用于计算由于地震导致的永久性位移.数据处理采用GIPSY软件，以各站24 h观测数据为一个单元，首先计算各站震前、震后各三天在ITRF框架下的坐标位置，坐标解算选用了IGS发布的最终精密星历和卫星钟差，采用PPP静态处理策略和各种标准模型(Zumberge et al.，1997)，以前后三日坐标均值的差分估计同震位移(图 1)，其水平误差为4 mm，垂直误差约为10 mm.境内GPS站的位移大致在5~30 mm范围，全部指向地震震中.境外各站中，位于极震区的三个站南向移动1.3~1.9 m，西向位移0.2~0.5 m，隆升达0.6~1.3 m.

我们采用PANDA软件(Shi et al.，2010; Fang et al.，2013)，按PPP瞬时精密定位(Larson et al.，2003)处理1 s采样数据，获取各站的地震动态位移.数据处理选用IGS发布的快速星历和钟差.瞬时坐标的解算精度水平向优于2 cm，垂直向优于3 cm. 1 s采样GPS站动态效应得到充分展示，其中XZZF、XZAR、XZRK水平向动态幅度达12 cm左右(图 1)，是静态同震位移的6~10倍.甚至600 km以外的LHAS站的波形信号也十分清晰.不过距离较近的XZZB站的波形信号微弱，多路径效应相对突出，动态波形的信噪比不高，为此特别利用震前三日的数据对该站地震动态位移时间序列进行恒星日滤波(Choi et al.，2004).

本文采用有限断层位错模型反演破裂时空过程.破裂被限定在走向295°、倾角11°、长210 km、宽为160 km的单一矩形平面内，该断层面与USGS公布的震源机制解相一致，也与喜马拉雅主滑脱断层的产状相似(Schult-Pelkum et al.，2005).模型断层面被进一步划分为21×8=168个子断层，每个子断层的尺度为20 km×10 km.其中走向方向上第11个、倾向方向上第4个子断层设为初始破裂点，该点对应于USGS公布的震源位置(28.147°N，84.708°E，depth=15 km).

反演基于多时窗叠加的线性算法(Hertzell and Heaton，1983)，将子断层的震源时间函数表示为5个上升时间为2 s的三角函数的叠加，通过非负最小二乘法线性解算每个三角函数的幅度，最佳拟合GPS站同震三维静态位移、动态波形及远震P波垂直分量，进而得到子断层的地震矩和滑动矢量.远场体波的格林函数计算基于反射法(Kikuchi and Kanamori，1991)和一维分层的全球速度模型PREM，对波形数据和对应的格林函数进行0.002~1 Hz带通滤波；动态GPS波形的格林函数则基于青藏高原一维速度结构模型(Wang et al.，2007)，使用频率波数法计算得到(Zhu and Rivera，2002)，而GPS静态位移的格林函数计算基于同样的岩石区速度结构，采用分层弹性半空间的位错模型计算(Wang et al. 2003).

约束破裂模型的数据由于其不同的空间分布和数据类型，对震源破裂过程具有不同分辨力(Yue and Lay，2013).本文设计了三组棋盘测试来检验用于反演的三种数据单独模型及联合模型的分辨力：(1)走向倾向均分布两个凹凸体于断层面边界上(2×2，图 2左列)；(2)走向倾向分别分布两个与三个凹凸体于断层面边界与内部(2×3，图 2中间列)；(3)走向倾向分别分布三个与两个凹凸体于断层面内部(3×2，图 2右列).在棋盘测试中，正反演所采用的断层模型、站点分布及格林函数的计算均与此次地震的实际情况相同，滑动矢量设计为逆冲.测试结果(图 2)表明：(1)联合模型在走向上和倾向上均具备最佳分辨力，很好地恢复了凹凸体的分布范围及滑动量值；(2)高频GPS与远震体波均具备较好的分辨力，且走向优于倾向；(3)同震永久位移的分辨力相对较弱，且断层面的东南部分辨力优于西北部的，主要原因是此次地震以南北向逆冲为主，使得分布在断层面东西向的static GPS约束能力较弱，而具备有效约束的南北向的站点分布于断层面东南部且数量较少(仅为CHLM，KKN4与NAST).

图 2

Fig. 2

图 2 棋盘测试 左、中、右列分为为2×2、2×3、3×2的测试；第一行为设计模型(Input model)，第二行为远震体波单独约束模型(Teleseismic)，第三行为同震永久形变单独约束模型(static GPS),第四行为高频GPS揭示的强地面运动单独约束模型(hrGPS),第五行为联合约束模型(Joint). Fig. 2 Examples of checkerboard tests The rupture models of left, middle and right array are 2×2 test, 2×3 test and 3×2 test, respectively. The first line are input models, and the second line to the last line are the inversions using teleseismic P waves, static GPS, hrGPS separately and jointly.

联合反演的结果如图 1、图 3a所示，地震主体破裂发生在断层面的东南段，呈单侧破裂的特征，走向上破裂展布约为100 km，扩展至加德满都以东30 km处，倾向上约130 km，向北延伸至我国西藏吉隆县区域，深度范围为5~25 km.破裂面滑动量区间为0.5~4.3 m，平均破裂幅度达2.4 m，最大滑动发生在约11 km深部，距离震源位置约35 km，位于加德满都西北面约20 km处，量值达4.3 m.在5 km以上的浅层部位，滑动沿倾向向上迅速衰减，断层浅表处未发生明显破裂.在深度为20~30 km 范围，滑动量随深度迅速衰减，我国吉隆县境内的破裂可达0.5~1.0 m，吉隆镇下方破裂达1 m以上.联合模型的地震矩释放率函数曲线(图 3b)及滑动历史(图 4)表明此次地震持续时间为80 s，主要子事件发生在15~40 s区间，破裂极值发生约35 s前后，为此次地震的最大滑动量，位于震中东南侧35 km、深度11 km处；次要子事件发生于约40~55 s区间，破裂极值发生于45 s前后，位于震中东南侧60 km、深度18 km处；55~80 s区间则为破裂向东南及深部传播，并伴随部分破裂抵达浅部但未上及地表，部分破裂传播至我国吉隆县下部；此次破裂的平均速度约为2.5 km·s^-1.根据破裂模型，本次地震释放地震矩约为9.4×10²⁰ N·m，对应的矩震级为M_W7.9，大于USGS发布((7～8)×10²⁰ N·m)的计算结果.

图 3

Fig. 3

图 3 断层面滑动分布与地震矩释放率函数 上图为联合反演的结果，下图为远震体波单一约束的结果.图(a,c)为断层面的滑动分布，红色五角星为震源位置；图(b,d)的红色曲线为地震矩释放率函数. Fig. 3 Slip distributions and moment rate functions The upper part depicts results by joint inversion, and the lower part shows results by teleseismic-only inversion. The contours of (a) and (c) denote the slip model, the color bar implicates the slip amplitudes, red star shows the epicenter. The red curves of (b) and (d) show the moment rate functions.

图 4

Fig. 4

图 4 联合模型破裂的时空过程 Fig. 4 Space-time history of fault slip inferred by joint inversion

与单一用远震P波的资料约束的破裂模型相比(图 3c，d)，联合反演模型在走向上延展出更大破裂范围，地震矩心位置也相应向东南方向移动16 km，此外联合反演显示5~6 km深度的浅层也有1 m左右的破裂.

GPS同震位移观测值和联合模型的模拟值及残差如图 5所示，拟合残差水平向均值为0.71 cm，垂直向均值为2 cm，其中CHLM、KKN4和NAST水平向拟合残差在3 mm内，由于这三个GPS站正好位于主破裂的上方，表明浅层的破裂状况得到了较好约束.基于联合约束破裂模型计算的GPS波形与观测波形的拟合状况见图 6，两者具有较好一致性，其中垂向与东西向的拟合度略优于南北向的.合成的远场P波与观测波形拟合较好(见图 7)，其平均互相关系数约为0.85，该结果表明联合反演的破裂模型较好地反映远场P波的起伏变化.

图 5

Fig. 5

图 5 GPS同震位移观测值与模拟值的拟合及残差图 (a,c) 水平向观测值与模拟值的对比与残差图；(b,c) 垂直向观测值与模拟值的拟合与残差图； 图a, b中蓝色系的位移矢量为观测值，红色系为模拟值，矢量值与图 1中的相对应. Fig. 5 GPS data fitting and residuals between observed and synthetic displacements (a,c) The data fitting and residuals in horizontal element, respectively. (b, d) The data fitting and residuals in vertical element, respectively. The blue arrows of figure (a) and (b) are observed value, as well as red arrows are synthetic value. The vectors of displacements in figure (a) and (b) correspond with the ones in Fig.1.

图 6

Fig. 6

图 6 高频GPS模拟值与观测值的对比 红色为观测值，黑色为模拟值.曲线左端上下方的数字分别为震中距(单位(°))与方位角(单位(°))，右端上方数字为最大幅度值(单位为m) Fig. 6 Comparison of observed and synthetic waves of HRGPS Red lines are the observed waves and black lines are the synthetic waves. In each sub-graph, the numbers on the left top and on the left bottom are epicentral distance (unit:(°)) and azimuth (unit:(°)), respectively; the number on the right top shows the max amplitude of observed wave (unit: meter).

图 7

Fig. 7

图 7 远震P波的模拟值与实测值的拟合图 红色为观测值，黑色为模拟值.曲线下方数值为两种曲线的互相关系数，单位为无量纲. Fig. 7 Comparison of observed and synthetic waves of teleseismic wave Red lines are the observed waves and blue lines are the synthetic waves. In each sub-graph, the number on the left bottom shows cross correlation coefficient between observed wave and synthetic wave.

图 8

Fig. 8

图 8 尼泊尔MW7.9地震构造背景及历史地震分布 Fig. 8 Tectonic setting and history earthquakes distribution of Himalayas earthquake zone

联合反演破裂模型显示本次地震发生在喜马拉雅深部低角度滑脱断层上.GPS观测表明印度板块以大约18~21 mm·a^-1速率沿该断层插入欧亚板块(Bilham et al.，1997; Avouac，2003)，占印度-欧亚大陆汇聚速率的一半.在尼泊尔东部，地震深反射资料显示该滑脱断层向北延伸到藏南的雅鲁藏布一带(Zhao et al.，1993)，向南在恒河平原最北缘斯瓦利克山前出露地表(Schulte-Pelkum et al.，2005)，滑脱断层以上显示为楔状薄皮构造，以挤压褶皱变形为主(Lavé and Avouac，2000)，而断层面以下的印度板块基本不变形，平缓插入青藏高原之下，控制了区域内变形的样式和空间展布.

喜马拉雅大震活动则与构造楔的弹性应变状态有关.根据GPS和微震观测，在两次大震间，位于喜马拉雅主峰以北的滑脱断层则自由蠕滑，以南的滑脱断层及其向下延伸的高角度断坡完全闭锁，宽度可达90~100 km(Feldl and Bilham，2006).积累的弹性应变大多通过闭锁的滑脱断层局部和整体破裂向喜马拉雅最南端的斯瓦利克山前传递，如此前1905年印度康拉M_W7.8(Wallace et al.，2005)、1934年尼泊尔和印度交界的比哈尔M_W8.1(Sapkota et al.，2013)、1950年中国西藏察隅M_W8.3级地震(Chen and Molnar，1977)以及2015年郎唐地震都是如此.此前普遍认为，此类地震破裂应从喜马拉雅主峰下方闭锁与蠕滑的交接部位开始，低角度向上扩展，部分甚至可能上达地表，如1934年比哈尔地震；但大多数如1905年康拉、1950年察隅以及本次地震，可能在应变传递的途中湮灭，形成所谓盲地震.不过联合反演结果显示此类地震破裂可能深入青藏高原下方自由蠕滑的滑脱断层，与2008年汶川地震的反演结果十分相似(Wang et al.，2011).

本次地震极震区与1934年比哈尔地震震区相接或重叠(Kumar et al. 2010)，几乎是1833年M_W7.7和1866年M_W7.2 两次地震的原地再现(Szelig et al.，2010).此前普遍认为1833—1866年序列地震不可能发生在滑脱断层部位，而是位于其上的高角度逆断层，如此1934年比哈尔地震才可能在时隔101 a在震区下方的滑脱断层爆发(Rajendran and Rajendran，2011).联合反演模型显示本次地震平均2.4 m左右的滑动幅度，基本上等于从1833年以来其闭锁导致的滑移亏失，这说明1833—1866年序列地震可能与本次地震破裂方式类似，并由此我们推测地震在此次地震破裂区域复发的周期可能为150~200 a.鉴于本次地震破裂面南端至喜马拉雅前缘推覆断层间滑脱断层仍保持闭锁，此后数十年间再次爆发灾害性地震的可能性依然存在.

联合反演GPS静、动态位移和远震体波一致显示，2015年4月20日尼泊尔M_S8.1地震是一次发生在喜马拉雅深部滑脱断层上的特大地震，破裂以逆冲为主，喜马拉雅整体向南仰冲，导致喜马拉雅南缘极震区内GPS测站南向水平位移2 m以上，垂直隆升1 m以上.地震破裂没有上达地表，但向下延伸至藏南吉隆县，位于震中东北400 km的日喀则县高频采样GPS站纪录到持续时间60 s、幅度为10~15 cm的动态位移，由于地震向东偏南方向单侧破裂，位于震中西北450 km的仲巴县GPS站动态位移仅在2~4 cm.本次地震破裂抵消了自1833年以来积累的弹性应变，但极震区以南的深部滑脱断层仍保持闭锁，未来仍有导致灾害性大震的可能性.

致谢 本文使用的地震波数据来源于IRIS数据中心，尼泊尔境内GPS数据来源于UNAVCO，GPS数据解算利用了GIPSY软件，部分图件绘制使用了GMT绘图软件.感谢两位匿名审稿专家提出的宝贵修改意见和编辑部老师的帮助！