北京时间2015年4月25日14时11分，在尼泊尔博克拉发生M_W7.9地震，震中位于28.1°N，84.7°E，震源深度15 km(由USGS(United States Geological Survey，美国地质调查局)给出).此次地震位于人口相对密集、建筑物抗震性能较差的地区，因此造成了尼泊尔境内的大量人员伤亡.中国西藏、印度、孟加拉国、不丹等地也均出现人员伤亡.此次地震发生在印度板块和欧亚板块的陆陆碰撞带上并且震级较大，所以也引发了雪崩、滑坡等次生灾害，造成重大人员伤亡和财产损失.本次地震发生在印度板块向北俯冲到欧亚板块的汇聚边界上，在其西部远处1505年发生过一次更大的地震，而紧邻本次地震的东部，在1934年发生过一次M8.0地震，该地震附近在1255年还发生过一次更大的地震.考虑到本次地震对加德满都乃至尼泊尔全境造成的严重损失，系统研究本次地震的破裂过程并准确给出能量释放极大的位置，对于揭示当地板块构造活动的规律和抗震救灾都有着非常重要的意义.

地震发生后，中国地震局地球物理研究所给出了有限断层反演的解(http://www.cea-igp.ac.cn/ tpxw/272110.shtml)，其结果表明主要能量释放区域在震中与加德满都之间的位置.不同于有限断层反演，在当地断层的几何形态未知的情况下，反投影方法是近十年来发展起来的获取大地震破裂过程的稳定有效的方法(Zhang and Ge，2010; Ishii et al.，2005; Koper et al.，2012).IRIS(Incorporated Research Institutions for Seismology，美国地震学研究联合会)网站给出了不同台阵的反投影结果(http://ds.iris.edu/spud/backprojection/9924467)，但由于单一台阵的方位角局限，这些结果表现出了较大的差异.为了提高破裂的分辨率，多重信号分类算法(Schmidt，1986; Goldstein and Archuleta，1991)也被引入了地球物理学领域来得到震源破裂的信息(Meng et al.，2011，2012).作为反投影方法和反演方法的结合，压缩传感算法(Donoho，2006; Candès et al.，2006; Malioutov et al.，2005)基于反演计算而不是直接的波形聚束，可以复现出稀疏分布的源而不依赖于任何的先验信息.这种方法已经被广泛使用在了核磁共振成像(Lustig et al.，2007)和计算摄影学(Duarte et al.，2008)等领域.在一个给定的时间，对于网格化的空间中，源的分布应该较为离散，数量较少，这一点符合稀疏性的前提条件，并满足震源的物理假设.因此，基于某一个台阵的大量数据，可以将压缩传感方法应用在地震学领域来对大地震的时空演化过程进行成像(Yao et al.，2011).但是，在传统的压缩传感方法中，为了数据之间更好的相关性，一般只会使用由一个台阵记录到的P波数据来进行反演.

考虑到一个台阵的方位角的局限性可能会带来一些误差，本文将传统的压缩方法进行改进，使用多个不同方位角的台阵的数据进行反演，并将其称为多台阵压缩传感方法.本文利用该方法，使用阿拉斯加、欧洲和澳大利亚三个台网共179个台站的远场P波垂直分量的数据来反演本次M_W7.9尼泊尔地震的破裂过程.

地震位错点源产生的位移场可以表示为(Aki and Richards，2002)

对于有限断层，考虑震源的有限破裂尺度效应，则(1)式可修改为

写成离散的形式，则有

其中下标k为子破裂元的编号，，τ_k，S_k分别为子破裂元k的位置坐标、破裂发生时间和破裂面积.由于考虑的是远场的P波垂直分量，并使用归一化的波形而非振幅的绝对值来计算每个子破裂元在频率域释放能量的相对强度，则意味着格林函数 G _ip，q中有关辐射花样和几何衰减因子的项可以忽略，从而格林函数 G _ip，q在频率域可简化为相移项e^iωt，其中t为从子破裂元到地震台站的走时.子破裂元的矩张量m_pqk=μS_ku_k.则(3)式可重新改写为矩阵形式

其中，x(ω)=(x₁，…，x_M)^T是具有复数值的子破裂元的向量，代表着M个待求解的子破裂元的相对滑动量值. b(ω)是频率域的在理论到时后的波形数据. A(ω)是由格林函数组成的矩阵. 在这里，由于每一条数据都已经归一化，因此只考虑相移对数据的影响，所以它的元素被设为A_nm=e^{iω(t_nm－t_n0)}.这里，t_nm是从震源格点 θ_m 到台站n的理论走时，t_n⁰是从震源到台站n的理论走时.震源破裂反演的最终目的是求解出 x(ω).

因为使用的是多个台阵的数据，因此对于基于单台阵的传统压缩传感方法的解的形式(Yao et al.，2011)就转化为

这里，j代表的是台阵的数量，i是某个台阵的序号.λ是一个正值的拉格朗日乘子.w_i是每一个台阵的权重因子，它的大小取决于这个台阵的孔径和数据质量.因此，这个反演问题的解是数据残差和解的稀疏性之间的折中.由于远震P波对深度几乎没有分辨，所以我们只在震源的深度平面内进行反演(Xu et al.，2009).

为了使用多台阵压缩传感方法反演此次地震的破裂过程，我们使用了来自三个台阵(欧洲 79台站，阿拉斯加 81台站，澳大利亚19台站)，共计179个台站的远场P波垂直分量的数据(图 1，图 2).首先，所有波形都进行0.2~5 Hz的带通滤波，然后分别根据每一个台站的数据最大值进行归一化.之所以选择这个频段进行滤波是因为过于低频的信息空间的分辨率太低，很难获得精细的细节信息；高频能量由于噪声过大可能导致结果的不稳定.然后，为了减少由于台站所在的接收区域的浅层速度结构的三维非均匀性对于走时的影响，我们先利用PREM模型(Dziewonski and Anderson，1981)计算出了由震源(28.147°N，84.708°E，15 km)到各个台站的理论到时，并通过计算最早的7 s数据的互相关函数，将同一台阵内的数据对齐到零时刻.随后，由于到达不同台阵的射线路径不同，也会影响到时，因此需要将每个台阵近似看作一个观测点，对齐所有来自不同台阵的数据.为了达到此目的，我们分别将每个台阵已经对齐好的数据分别叠加起来获得三条叠加的数据，从而提高信噪比，达到P波初动可以清晰识别的目的.最后，将叠加好的数据通过人工标定的P波初动时间作为参考时间来对齐并进行相应的极性翻转(图 2).

图 1

Fig. 1

图 1 尼泊尔地震的震中位置以及所使用的台阵数据的分布 黑色五角星代表震中的位置，蓝色圆形代表台站.Fig. 1 Map of the Nepal earthquake location and distribution of stations from three arrays Black star is the epicenter. Blue circles are the stations.

接下来以震源为中心，将震源所在区域划分为101×101个空间间隔为0.05°的格点，每个格点作为可能的子震源用于反演.在本文中，由于地震波形的数据资料质量非常好(信噪比约为29.4)，并且台阵的孔径相差不大，设权重因子w₁=w₂=w₃=1，并采用(Yao et al.，2011)，其中r为信噪比的倒数，N为台站总数，所以设λ=0.45.

为了获得此次地震震源区域不同位置能量释放的时间分布，我们使用了一个长度为10 s的滑动窗，滑动的步长为2 s.第一个时间窗被设定为只包含P波初动的波形，这是为了更好地在结果中识别出破裂发生的初期的性质.在反演计算的过程中，每 一个时间窗的数据都转换到频率域，并计算在0.25~1 Hz 频段范围内的解.之所以选择这个频段是由于数据的频谱决定的(图 2d).如果频段选择的过低，则无法看到破裂的细节；如果选择的过高，则会由于高频的能量较少信噪比较低而使得解不够精确.在反演过程结束后，对于每一个选定的时间窗，我们得 到了不同频率下能量释放的空间分布，并将这些分布叠加起来，得到在0.25~1 Hz频带范围内的能量分布.最终，将每一个时间窗的解进行二维的空间高斯平滑，得到了震源区域能量的释放位置随时间的变化过程，即此次地震的时空演化过程.

图 2

Fig. 2

图 2 本文计算中所使用的数据的信息 (a)欧洲，(b)阿拉斯加，(c)澳大利亚三个台阵上方的彩色图案代表地震波形数据及各台阵以P波初值对齐后叠加的地震图，箭头标示P波初动；(d)图为台网中代表性台站的数据的频谱的范例.Fig. 2 The information about used data (a)Europe，(b)Alaska，(c)Australia. In the figure of each array，the top colored image is the aligned waveforms. The lower image is the stacked data，while the arrow indicates the first pulse of P wave.(d)shows the frequency spectrum of three records from three respectively as an example.

由于多台阵压缩传感方法是基于反演的算法，所以数据中的噪声和一些反射波的震相(例如pP)可能使得反演的结果中有一些假象，因此就非常有必要在通过反演直接求出的所有解中挑选出真实的解.为了达到判定解是否真实可信的目的，设定了一个判据.首先，定义：

这里，对于第q个时间窗，P_mq是按由大到小排列的第m个区域极大值，M_q是相对残差，定义为：

本文取m=1，也就是每个时间窗所求得的解只取极大值作为待选.通过公式(6)，可以计算出α随时间窗变化的图像(图 3).然后，再定义一个鉴别因子(DF).这意味着，所有α值小于DF的时间窗的解将会被认为是假象而舍去.在本文研究中，DF=0.03，这与所使用的数据的信噪比(SNR=30)的倒数是几乎一致的.但是由于在震源开始的破裂阶段能量释放较小，因此使用此判据会丢失掉解中初期的信息，所以此判据对于刚开始破裂的前7 s并不适用.也就是说，虽然前7 s的解的α小于DF，但是依然保留作为最后的解.图 4和图 5分别给出的是筛选前和筛选后的结果.可以看到，经过筛选，一些能量相对较小，残差相对较大的解被删去，这有助于我们更精确地估计震源破裂持续时间和断层破裂的尺度，获取破裂特征信息.

图 3

Fig. 3

图 3 蓝线代表α随时间变化的曲线，即不同时间窗中的解的α的变化.红线代表本次选定的鉴别因子的 值，这意味着除前7 s外红线下方的解是被舍弃的.Fig. 3 α versus rupture time. The red line is the value of the chosen DF which means that the results below the red line are removed except the results of first 7 s.

图 4

Fig. 4

图 4 对结果进行筛选前的地震破裂过程 黑色五角星代表本次地震的震中位置，彩色圆圈代表能量释放位置，颜色代表时间，圆圈的大小代表释放能量的相对大小.Fig. 4 The raw rupture process of 2015 Nepal earthquake before filtering Black star is the epicenter of the 2015 earthquake. Colored circles denote the location of power radiation. Color means rupture time. The radius of circle is proportional to the relative value of power radiation.

图 5

Fig. 5

图 5 筛选后的2015尼泊尔地震的最终破裂过程 黑色五角星代表本次地震的震中位置，彩色圆圈代表能量释放位置，颜色代表时间，圆圈的大小代表释放能量的相对大小.灰色圆圈代表震后5天内的余震分布.红色五角星代表 1833年和1934年的震中位置. 虚框的椭圆为2015和1934两次地震的估计断层破裂范围.Fig. 5 The final rupture process of 2015 Nepal earthquake after filtering Black star is the epicenter of the 2015 earthquake. Colored circles denote the location of power radiation. Color means rupture time. The radius of circle is proportional to the relative value of power radiation. Grey circle denotes aftershocks within 5 days after the main shock. Red stars are the location of 1833 and 1934 Nepal earthquake. Dashed ellipses are the estimated rupture zones of 2015 and 1934 earthquake.

经过筛选，最终得到了这次地震的破裂过程(图 5).结果显示，本次地震从震中开始，向南东东方向破裂约105 km.其中，破裂分为两个阶段.第一个阶段是地震发生后的前15 s，破裂所释放的能量基本围绕在震中附近10 km以内，并没有向远处传播.第二个阶段从第16 s开始，破裂从震中东部约40 km处开始向南东东方向传播，总的破裂持续时间为59 s.值得注意的是，破裂到断层最东边的时间为51 s，此时破裂尚未结束.在破裂的最后阶段，破裂能量集中的位置又回到了此次地震能量释放最强的地方(28°N，85.4°E).这意味着这个位置从29 s开始破裂，直到59 s才结束，持续时间长达30 s，持续释放了较长时间的能量.并且可以判断此处应该是受灾较为严重，释放能量较大的区域.

从上述解中提取出破裂前沿的位置，并沿断层走向投影，可以进一步得到在破裂离开震中附近到达震中东部40 km处之后，破裂向南东东方向以约 1.9 km·s^-1的速度传播(图 6)，到达距离震中105 km 处终止.

图 6

Fig. 6

图 6 破裂前沿与震源之间的距离随破裂前沿传播的时间的变化 黑线给出的是线性拟合出的第二阶段中的破裂速度，为1.9 km·s-1.Fig. 6 The distance from the rupture front to the epicenter versus rupture time The black line is obtained by linear regression to fit time-distance relations of rupture process. Its slope is the mean rupture velocity with a value of 1.9 km·s-1.

对比震级极其相近且也是发生在陆地板块间的2008年汶川地震，可以发现两次具有巨大破坏力的地震有着明显的几个相同点：(1)破裂都为单侧破裂:汶川地震的破裂方向是从震源开始向东北方向沿龙门山断裂带发展(张勇等，2008；Zhang and Ge，2010)，尼泊尔地震是向东南东方向沿印度板块与欧亚大陆板块的俯冲带走向的反方向上发展，均没有向另一方向破裂的能量释放.(2)两次地震都存在明显的阶段性: 汶川地震在初始阶段的0~20 s，看不出明显的破裂方向(Zhang et al.，2014)，尼泊尔地震前15 s破裂也都围绕在震源附近.两次地震都是在初始的小部分能量在震源附近区域释放后的一段时间内才表现出明显的破裂方向性和主要能量释放特征.(3)能量释放最大的位置均不在震中: 汶 川地震能量释放的最大位置在震源东北方向约50 km 处(Zhang and Ge，2010)，尼泊尔地震在距震中70 km处.然而，这两次地震也有很明显的不同点：(1)破裂速度和尺度：汶川地震的破裂速度约为3 km·s^-1，破裂尺度约为300 km(张勇等，2008；Xu et al.，2009；Zhang and Ge，2010)；而尼泊尔地震的第二阶段有明显方向性的破裂速度为1.9 km·s^-1，长度为105 km.经过对比可以发现，尼泊尔地震无论是破裂尺度还是破裂速度均要小于汶川地震.这表明尼泊尔地震的能量释放更加集中.(2)地质资料和地震形成的地表破裂轨迹表明，汶川地震的破裂并不是发生在单一的断层上，而是发生在龙门山构造带的映秀—北川断裂和灌县—江油断裂，中间的能量释放在空间上有一个跳跃(王卫民等，2008).从大的地质背景来看，汶川地震属于欧亚大陆板块内部的地震.而尼泊尔地震并没有发现跨断层触发的明显证据，属于板间地震.

另一方面，可以看到在破裂的第一阶段和第二阶段空间位置的中间存在一片没有明显能量分布的空区.同时这一空区的位置也和本次地震后五天内余震的分布的空区相一致(图 5).这表明该部分地区无论是震时还是在震后在高频部分均没有较大幅度的能量释放，两者分布一致，证明了利用多台阵压缩传感算法反演得到的能量释放位置是可信的.而且，根据我们的结果，可以看到能量释放最集中的区域就在尼泊尔首都加德满都附近，而加德满都也是本次受灾最为严重的地区.由于造成房屋毁坏倒塌，破坏力较强的地震波为相对高频的波动，相比利用有限断层反演的低频的结果，高频能量辐射位置会与破坏较为严重的位置更为一致.因此，利用多台压缩传感方法来反演大地震的破裂过程可以快速得到能量释放最为集中、破坏程度较为严重的地区的地理位置，这给震后第一时间的快速救援提供了很大的帮助，同时体现出了在相对较高频率进行破裂过程反演在这一方面具有一定的优越性.同时，我们将本次地震在0.25~1 Hz频段释放的能量累加之后，可得到高频总能量释放的空间分布(图 7).多台阵压缩传感方法给出的结果与利用INSAR数据给出的破裂范围基本相同(http://e-decider.org/content/nepal-earthquake-april-2015)，这证明我们有关破裂尺度的估计是可信的.以上对比证明，通过多台阵压缩传感方法得到的高频能量释放区域的解与其他类型数据得到的结果具有一致性.

图 7

Fig. 7

图 7 本次地震在0.25~1 Hz能量释放的空间分布 红色五角星代表震中位置，虚框矩形为INSAR给出的断层破裂范围.Fig. 7 The spatial distribution of power radiation in frequency b and of 0.25~1 Hz The red star denotes the epicenter of 2015 Nepal earthquake. Dashed rectangle is the range of rupture observed by INSAR.

本次研究所用到的多台阵压缩传感方法使用了 来自三个台阵的数据，目的是减轻由于台阵方位角和孔径的局限性而导致的对于破裂细节的分辨率下降的问题.我们通过计算1 Hz时的台阵响应函数(ARF)(Rost and Thomas，2002；Xu et al.，2009)来判断这样的台站分布是否可以有效提升分辨率.由图 8可以发现，在选取同样的置信区间的时候(ARF＞0.7)，使用多台阵分布可以使解的范围更加精确，等值包络线范围最小，而选取其他单一台阵的时候的精度较低，包络范围较大.同时，我们发现，若只使用单一台阵，例如阿拉斯加或欧洲台阵的数据，椭圆的倾斜方向即长轴方向与破裂方向较为一致，因此可能对于本次地震的破裂细节的反演精度较低，在破裂方向上分辨率较差.这表明台阵中台站的分布以及台阵的方位角都可能影响到最终解的精确度.值得注意的是，在使用多台阵数据计算时，虽然解的精度有较为明显的提升，但是台阵响应函数的结果显示其分布并非是一个圆形，而在南北方向可能出现一些误差，但是考虑到此次地震的破裂方 向为南东东，因此对于整体破裂信息的获取影响不大.

图 8

Fig. 8

图 8 各个台阵在1 Hz的台阵响应函数(ARF)大于0.7的包络线 紫色五角星为震中位置.箭头代表指向台阵的方向.Fig. 8 The envelopes of array response functions(ARF)larger than 0.7 of three arrays，and their sum，at 1 Hz The purple star denotes the epicenter. The black arrow represents the direction toward array.

本次地震发生的位置是在从南向北俯冲的印度板块和在北边上覆的欧亚板块的俯冲带上，这一区所处的断裂带被称为主喜马拉雅逆冲断裂(MHT).通过GPS数据计算可以得到在MHT上，从地表到沿倾角向下的约100 km是锁定的区域，近20年的观测也表明这一地区几乎没有慢滑移事件，也就是说累积的能量都需要通过地震来释放(Ader et al.，2012).作为两个板块的边界，历史上这个地方也曾发生过一些大地震.在本次地震的发生区域的方圆250 km内曾经在过去的一个世纪发生过4次6级以上的大地震.最大的一次是震级大于M8.0，发生 在1934年(Hough and Roger，2008；Sapkota et al.，2013； Bollinger et al.，2014)，断层破裂范围至少从 85°52′E至87°17′E，最小破裂长度为150 km(图 5). 从我们的结果可以看出，2015地震破裂以及余震活动的分部刚好延伸至1934地震的破裂带的西缘，两个破裂带紧邻，所以2015地震可以被看作是1934地震所处断层的另一区域的能量释放.在更早的19世纪，1833年也曾有一次M_W7.6的地震(Ambraseys and Douglas，2004)发生在与2015年地震相同的发震断层(图 5).为了研究这两次地震之间的关系，我们发现由GPS数据得到的此区域的板块汇聚速率为17.8±0.5 mm·a^-1，但是由于MHT浅部的锁定，其滑移速率极慢(<5 mm·a^-1)，所以这一区域的滑移缺损累积速率至少为13 mm·a^-1(Bilham and Ambraseys，2005)；MHT在尼泊尔地区的震间的地震矩缺损累积速率为6.6±0.4×10¹⁹Nm·a^-1(Ader et al.，2012).那么根据USGS的有限断层反演结果(http://earthquake.usgs.gov/earthquakes/eventpage/us20002926 scientific_finitefault)得到本次地震的最大滑移量约为3.3 m，平均滑移量约为2 m，标量地震矩约为0.8×10²¹ N·m.可以推断出：本次地震的平均滑移量与自1833年积累至今的180年的缺损所积攒的滑移缺损总量几乎相等；但本次地震的标量地震矩仅为地震矩缺损累积总量的约十分之一，可以认为地震矩缺损已经释放了一小部分.因此从该断层破裂区域的小尺度范围看，我们认为本次地震可以被看作是1833年地震后，局部能量积累的再次释放.但是由于尼泊尔境内的MHT积累的地震矩缺损绝大部分能量还没有释放，因此整个尼泊尔地区在未来一段时间仍不排除有继续发生更大地震的可能.

本文利用台阵压缩传感方法给出了2015尼泊 尔M_W7.9地震的破裂过程和破裂细节. 结果表明此次破裂是一个典型的单侧破裂，破裂方向为南东东，破裂尺度约为105 km，破裂传播速度为1.9 km·s^-1，破裂累积能量最大值位于尼泊尔首都加德满都附近.都是造成加德满都震害严重的原因.

通过与震级相近的2008汶川地震的细节破裂特征进行比较，我们发现两次地震都是单侧破裂、存在明显的分段特征且最大能量释放位置并非震中；但是，与汶川地震相比尼泊尔地震的破裂尺度更短，说明其能量释放更加集中.

根据历史上曾经发生在尼泊尔境内的主喜马拉雅逆冲断裂(MHT)上的地震以及板块的相对运动速率，推测本次2015尼泊尔地震可能是继1833年M_W7.6地震后，在该断层上又一次小尺度的180年间积累的能量释放的过程，不排除整个尼泊尔地区发生更大地震的可能性.

致谢 本文地震波数据来自IRIS 数据中心，部分图形的绘制使用的是GMT绘图软件. 两位匿名审稿专家提出了非常中肯的修改意见，作者对他们致以诚挚的感谢！