北京时间2015年4月25日14点11分26秒，临近我国西藏自治区的尼泊尔境内发生了M_W7.9强烈地震.UGSG测定震源位置为28.147°N，84.708°E，震源深度15 km.震中位于尼泊尔博克拉市，距离首都加德满都西北77 km处，是自1950年M_W8.7察隅地震以后尼泊尔近70年来受到影响最大的地震(Mittal and Kumar，2015).此次尼泊尔地震发生在青藏高原南部与印度板块的边界处，位于全球陆地碰撞和垂直变形最剧烈的区域，受印度板块与欧亚板块陆陆俯冲碰撞和推挤(图 1).在此区域，历史上曾经发生过多次大地震.其中距离此次地震最近的是1833年8月26日尼泊尔加德满都北部M_S8.0地震和1934年1月15日加德满都M_W8.0地震(李培等，2015).截至2015年6月11日，地震至少造成8786人死亡，22303人受伤，中国西藏、印度、孟加拉国、不丹等地均出现人员伤亡.

图 1

Fig. 1

图 1 2015年4月25日尼泊尔MW7.9地震序列(USGS).主震由红色五星表示，黄色圆点代表M>3.0的余震 Fig. 1 Map of the aftershock distribution of the April 25,2015,Nepal earthquake from USGS.The main shock epicenter is indicated by red star.Yellow points represent the aftershocks of M>3.0

众所周知，主震后会伴随很多余震，这是因为主震会改变地壳的应力场从而引发余震甚至触发周边主震.Dieterich(1994)曾用均匀应力降的不规则体来模拟断层从而研究余震发生率的空间分布.Console和Catalli(2006)利用均匀应力降和非均匀应力降下的矩形位错模型模拟了依赖于时间的地震发生率.从物理的角度来看，在不考虑塑形形变的情况下，均匀应力降下的位错模型所表现出的无穷大的应力变化可能导致无穷的地震发生率；同时，均匀应力降下的位错模型同样也表明在地震滑移区内不应该发生余震.研究表明，余震大多发生在主震断层面大滑移区的外面或者边界(Benioff，1955;Das and Scholz，1981;Eaton et al.，1970;Henry and Das，2002).而Das和Henry(2003)发现一些小的地震通常会发生在断层面大滑移区内，且余震常常发生在应力增加的区域和滑动位移增速较快的区域.近期，Woessner等(2006)研究了余震与主震震源属性的空间相关性.研究表明，应力变化对于余震的空间分布相关性不大，这就意味着主震后静态库仑应力变化不是产生余震的主要因素.Helmstetter和Shaw(2006)研究了应力的各向异性同滑动速率与状态相依赖的摩擦定律下的地震活动性速率的关系，发现余震多分布在低滑动地区，且与简单的滑动位错模型相比空间分布较为分散.

在本研究中，我们基于张勇等(2015)以及Wang和Fialko(2015)给出的断层面滑动位移反演结果，将Ripperger和Mai(2004)提出的傅里叶变换法应用于2015年尼泊尔M_W7.9地震中，从而推出二维断层面上的静态应力分布，并与余震序列的空间分布相比较，从而获得更为合理的震源滑动模型机制.

在一个无限均匀各向同性弹性介质中，假设断层在直角坐标系下的x₂=0平面上，滑动位移D(x₁，x₃)是沿着x₁的方向.Ripperger和Mai(2004)给出了在波数域中估计断层面上的静态应力分布的方法.该方法能够将断层面上的滑动位移分布转换为平行于位移方向和垂直于位移方向的应力状态分布，如公式(1)(Andrews，1978，1980)：

(1)

其中，k是波数，D是静态滑动位移分布，Δσ₁和Δσ₃分别是平行于滑动位移方向和垂直于滑动位移方向的剪切应力变化.K₁和K₃分别为静态刚度矩阵(Ripperger and Mai，2004).我们将拉梅常数用杨氏模量E和泊松比ν代替，K₁和K₃可表示为(式中杨氏模量可抵消)：

(2)

(3)

其中，k₁和k₂分别为与滑动方向平行向和垂直向的波数.可见，当滑动位移被傅里叶变换到波数域后，只要给定泊松比，我们就可以计算出沿着滑移方向和垂直于滑移方向的静态应力变化.假设泊松比在地壳中是一致的，那么K₁和K₃仅与剪切模量有关.

尼泊尔M_W7.9地震发生后，很多学者运用不同方法和手段计算了尼泊尔地震破裂分布(刘刚等，2015;单新建等，2015;苏小宁等，2015;张勇等，2015;Wang and Fialko，2015)：刘刚等(2015)用连续GPS与远震体波联合反演的方法得到了尼泊尔地震的破裂过程，并推测了此次地震破裂区域的复发周期；单新建等(2015)利用日本ALOS-2和欧空局Sentinel-1A卫星获得的尼泊尔地震同震形变场，结合GPS同震位移数据，联合反演了断层滑动分布特征和空间展布；苏小宁等(2015)依据GPS观测到的同震位移场，利用非负最小二乘方法反演得到了震源断层面上的滑动分布；张勇等(2015)除了利用同震位移明显的两个台站的GPS位移数据外，还利用远震地震波数据、全球地震台网(GSN，Global Seismographic Network)震中距90°以内的宽频带地震数据，对破裂模型进行了反演，并根据反演结果与历史大地震的破裂分布情况，简单定性讨论了此次地震对加德满都附近区域地震活动性和危险性的影响；Wang和Fialko(2015)利用了包含干涉合成孔径雷达的表面形变测量数据，结合全球定位系统数据计算得到了2015尼泊尔M_W7.9地震的滑动位移反演结果.综合考虑上述反演过程，我们选取了利用台站的GPS位移数据、远震地震波数据、全球地震台网宽频带地震数据联合反演的张勇等(2015)的反演结果(图 2a)以及包含干涉合成孔径雷达的表面形变测量数据的Wang和Fialko(2015)的结果(图 2b)对尼泊尔地震进行讨论.

图 2

Fig. 2

图 2 (a)张勇等(2015)与(b)Wang和Fialko(2015)反演滑动位移结果 黑色五角星为主震位置，蓝色五角星为5月12日MW7.2地震位置. Fig. 2 Results of joint inversion by Zhang et al.(2015) and Wang and Fialko(2015) Black stars represent the main shock epicenters;blue stars represent the 5.12 MW7.2 earthquake epicenters.

张勇等(2015)反演得到的最大静态滑动量为5.2 m，断层面走向301°，倾角11°.Wang和Fialko(2015)反演得到的最大静态滑动量为5.8 m，断层面走向285°，倾角7°.

主震发生后短时间内发生的余震可以被认为是断层面上应力调整所导致，为了清晰地分辨出每个时间段内余震发生的空间位置，我们将余震分成12 h、12~24 h、24~48 h和48 h到5月12日M_W7.2地震之前四个阶段.此次尼泊尔地震序列发生在中国与尼泊尔交界地区，在两国均有余震分布，由于两国地震台站分布不均，且各自的地震目录均不完善，所以我们选用USGS公布的地震目录作为研究对象.根据USGS公布的地震目录，4月25日到5月12日M_W7.2地震前，共有M>3.0余震156次，其中12 h内发生72次，24 h内发生87次，48 h内发生110次余震.应用Ripperger和Mai(2004)的模型，假设泊松比为0.25，剪切模量为3×10⁴ MPa，通过波数域中的傅里叶变换，我们分别计算了上述两个模型沿断层走向和倾向两分量的同震应力变化(图 3a，3b和图 4a，4b)和应力变化梯度(图 3c，3d和图 4c，4d).该应力变化梯度是对a和b图中两方向的应力变化求梯度得到的.同时，由于此次地震断层面为浅角度逆断层，且从UGSG给出的地震目录来看，绝大多数余震都发生在10 km深度以内，所以我们近似将余震投影到地表.在反演结果中，由深及浅存在滑动角向南到向西的变化(图 2)，而主断层又是以逆冲为主的浅角度断层，因此我们在计算断层面上的应力分布时，仅考虑倾向方向的滑动位移带来的影响.

图 3

Fig. 3

图 3 (a)和(b)分别为利用张勇等(2015)模型得到的2015尼泊尔地震沿走向和倾向两个方向断层面上的静态应力变化以及余震在地表的投影，(c)和(d)分别为与(a)和(b)相对应的静态应力变化梯度以及发生在图(a)和(b)应力下降区域的余震在地表的投影.黑色、绿色、蓝色和白色圆圈分别为12 h、12~24 h、24~48 h以及48 h到5.12 MW7.2地震之前发生的余震；黑色五角星为主震位置，蓝色五角星为5月12日MW7.2地震，应力变化和应力梯度单位均为MPa Fig. 3 (a)and(b)are along-strike and down-dip components of static stress changes for the Nepal earthquake based on the kinematic slip distribution inversion form Zhang et al.(2015),while(c)and(d)are the shear-stress gradient respectively.Black,green,blue and white circles represent the aftershocks in 12,24,48 hours and from 48 hours to 5.12 MW7.2 earthquake of M>3.0,black stars represent the main shock epicenters,blue stars represent the 5.12 MW7.2 earthquake epicenters.The units of stress change and stress gradient are both MPa

图 4

Fig. 4

图 4 (a)和(b)分别为利用Wang和Fialko(2015)模型得到的2015尼泊尔地震沿走向和倾向两个方向断层面上的静态应力变化以及余震在地表的投影，(c)和(d)分别为与(a)和(b)相对应的静态应力变化梯度以及发生在图(a)和(b)应力下降区域的余震在地表的投影.黑色、绿色、蓝色和白色圆圈分别为12 h、12~24 h、24~48 h以及48 h到5.12 MW7.2地震之前发生的余震；黑色五角星为主震位置，蓝色五角星为5月12日MW7.2地震，应力变化和应力梯度单位均为MPa Fig. 4 (a)and(b)are along-strike and down-dip components of static stress changes for the Nepal earthquake based on the kinematic slip distribution inversion form Wang and Fialko(2015),while(c)and(d)are the shear-stress gradient respectively.Black,green,blue and white circles represent the aftershocks in 12,24,48 hours and from 48 hours to 5.12 MW7.2 earthquake of M>3.0,black stars represent the main shock epicenters,blue stars represent the 5.12 MW7.2 earthquake epicenters.The units of stress change and stress gradient are both MPa

图 3a，3b和图 4a，4b中应力为负代表了应力的下降，且应力分布显示出了极大的不均匀性.可以看出应力下降的区域比图 2中滑动量为正的区域明显缩小了.由于2015年尼泊尔地震是以逆冲为主，因此倾向方向的静态应力变化远远大于走向方向的应力变化.表 1给出图 3a，3b和图 4a，4b中两个分量应力变化的最大值和最小值.

表 1(Table 1)

表 1 2015年尼泊尔地震断层面上两分量应力极值 Table 1 The maximum and minimum stress changes of 2015 Nepal earthquake 走向分量 倾向分量 最大值 /MPa 最小值 /MPa 最大值 /MPa 最小值 /MPa 张勇等(2015) 0.98 -2.03 4.68 -28.02 Wang和Fialko(2015) 0.46 -0.84 5.54 -10.15

表 1 2015年尼泊尔地震断层面上两分量应力极值 Table 1 The maximum and minimum stress changes of 2015 Nepal earthquake

从图 3和图 4中可以看出，余震分布较为集中，主震周围余震较少，反而在断层面另一侧余震分布较多.表 2和表 3给出了余震在四个不同时间段内的分布情况.在距离主震较近的时间段内(12 h内)，在应力变化主导方向倾向分量上，分别在图 3和图 4中，有72%和69%的余震分布在应力增加的地方，随着时间的推移该百分比稍有变化，而且5月12日M_W7.2地震也发生在应力增加的区域.这说明主震过后，断层面上经历一个短暂的应力调整的过程，应力增加的地方，通过余震进行能量的释放.虽然有少部分余震发生在应力下降的区域，但从图 3c，3d和图 4c，4d中可以看到这些地震多发生在应力变化梯度较大的地方.由于图 3c，3d和图 4c，4d相比，梯度相差近一个数量级，所以图 3和图 4中，d图更能说明问题.忽略滑动位移反演结果的不唯一性，应力变化梯度大的地区可能在主震发生时没有完全破裂，而是发生了震后滑移(afterslip)现象.Perfettini和Avouac(2007)指出，尽管震后滑移与库仑应力变化对区域介质的影响机制不同，但也可能导致余震的发生.

表 2(Table 2)

表 2 余震在四个不同时段的分布情况统计(利用张勇等(2015)的模型) Table 2 Statistics of the distribution of the aftershocks in four periods(Zhang et al.(2015)model) 时段 总个数 走向分量 倾向分量 正值 负值 正值 百分比 正值 负值 正值 百分比 12 h内 72 54 18 75% 52 20 72% 24 h内 87 62 25 71% 61 26 70% 48 h内 110 78 32 71% 77 33 70% 48 h后 156 110 46 71% 107 49 69%

表 2 余震在四个不同时段的分布情况统计(利用张勇等(2015)的模型) Table 2 Statistics of the distribution of the aftershocks in four periods(Zhang et al.(2015)model)

表 3(Table 3)

表 3 余震在四个不同时段的分布情况统计(利用Wang和Fialko(2015)的模型) Table 3 Statistics of the distribution of the aftershocks in four periods(Wang and Fialko(2015) model) 时段 总个数 走向分量 倾向分量 正值 负值 正值 百分比 正值 负值 正值 百分比 12 h内 72 36 36 50% 50 22 69% 24 h内 87 43 44 49% 59 28 68% 48 h内 110 56 54 51% 75 35 68% 48 h后 156 79 77 51% 107 49 69%

表 3 余震在四个不同时段的分布情况统计(利用Wang和Fialko(2015)的模型) Table 3 Statistics of the distribution of the aftershocks in four periods(Wang and Fialko(2015) model)

滑动位移图案多数情况下可以被看作是一个非均匀应力变化的圆形破裂，最大的滑动位移与应力降(负应力变化)都发生在圆盘模型的中间部位(在张勇等(2015)模型中更为突出).我们将半径为a的圆盘断层模型中的静态滑动位移表示为

(4)

其中Δu_o是最大滑动量，H(x)是赫维赛德函数，m是正常数(对于均匀应力降模型，m=0.5)，r是圆心到断层上某点的距离.在波数域中，我们将其表示为(Singh，1977)

(5)

其中，J_n(x)是n阶贝塞尔函数，k_r是径向波数.对于一个很小的x，如果n给定，x趋于零时，J_n(x)≈ ，而对于一个很大的x，在n给定的情况下，.所以，对于非常小的k_r来说，Δu(k_r)是一个常数；而对于很大的k_r来说，Δu(k_r)~1/(k_ra)^m+1.5.对圆盘断层模型，我们遵循Keilis-Borok滑动模型，且令m=0.5，Δu就近似为(1-r²/a²)^1/2，因而在面积为A(A=πa²)的模型中得到的应力变化就是一个常数.如果要让圆盘模型边界上(r=a)的应力变化为有限值，m就必须大于1，这就意味着Δu(r)的一阶导数在r=a的地方是连续的(Singh，1977).

我们暂且利用张勇等(2015)的模型，计算出了滑动位移谱和静态应力变化谱，得到滑动位移谱近似按照k^-3(k^-3.0621)衰减，而与之相对应的应力变化谱是按照k^-2(k^-2.0453)衰减的(图 5a).图 5a中的蓝点代表了滑动位移和应力状态在波数域中不同波数下的谱值，红色点线代表了平均值，绿色点线是红色点线的包络线.从公式Δu(k_r)~1/(k_ra)^m+1.5中我们发现，滑动位移方程可近似成m=3/2，而应力状态在破裂面内成二次衰减.Singh(1977)指出，尽管Brune(1970，1971)和Dahlen(1974)的模型是用不同方法得到的，但是假定破裂速率等于剪切波的传播速率，对于S波来说，m=3/2时，Dahlen(1974)给出的圆盘模型的半径a与Brune(1970，1971)给出的非常相似.所以，m=3/2可能与真实情况更为接近.

图 5

Fig. 5

图 5 2015年尼泊尔地震滑动位移和应力状态谱，红色点线代表了平均值，绿色点线是红色点线的包络线；(b)基于圆盘模型得到的一维归一化的滑动位移(虚线)与静态应力(实线)曲线,其中re是有效半径，在该点应力改变正负号 Fig. 5 (a)Spectral representation of the slip displacement and stress changes on the fault of the Nepal earthquake modeled by Zhang et al.(2015).The red points represent the average value,and the green lines represent the envelopes of the red points.(b)The normalized slip(dashed line)and static stress(real line)corresponding to the circular slip distribution.Here,re is the equilibrium radius,which indicates where the stress drop changes sign

图 5b表示了满足圆盘滑动位移Δu=(1-r²/a²)^3/2(Segall，2010)的一维归一化滑动位移(虚线)和应力状态(实线).在r=a的位置，滑动位移为零，而应力却正向最大，其绝对值与r=0的位置相等.有效半径r_e代表了应力变化反转符号的位置，r_e约等于0.7a.在r_e与a之间应力变化为正，而r_e点的滑动位移Δu(r_e)约为圆盘中心点滑动位移Δu(0)的40%.进一步，如果真实的滑动位移满足k^-3衰减(m=3/2)而非均匀应力降模型中m=1/2，那么圆盘模型中应力增加的部分应该占整个圆盘破裂面积的51%.这次2015年尼泊尔M_W7.9地震的破裂过程近似遵循滑动位移k^-3衰减的现象可以部分解释余震的空间分布规律.事实上，Yamashita(1978)强调了非均匀应力降模型在理解断层动态破裂过程中的重要作用.Yamashita(1978)指出，如果地震是由初始应力累积造成的，且其滑动速率小于破裂速度，那么震源的滑动位移方程可以表示为Δu(r)~Δu_o(1-r²/a²)^3/2.Δu(r)在破裂尖端是有限的，且与之相对应的非均匀应力降模型中的初始应力满足Δσ_S(1-r²/a²)^1/2，其中Δσ_S是静摩擦应力.因此，2015年尼泊尔地震的滑动位移和应力变化图案说明了有限的应力变化造成的滑动位移模型可能比均匀应力降模型或Keilis-Borok模型更加接近真实情况.

基于张勇等(2015)、Wang和Fialko(2015)给出的断层面滑动位移反演结果，结合Ripperger和Mai(2004)提出的傅里叶变换法，我们得到了2015年尼泊尔M_W7.9地震二维断层面上的静态应力分布，通过与余震序列的空间分布相比较，获得了更合理的震源滑动模型机制.然而，在计算过程中，存在以下几方面不确定性因素.

首先，应力变化是依赖滑动位移的反演结果的，不同的反演模型可能导致不同的应力变化图案，尽管不同的反演模型给出的滑动位移具有相似的图案，但是不同模型之间仍然存有差别.我们最终选取了张勇等(2015)、Wang和Fialko(2015)的结果，是因为他们不仅仅利用了远震地震波数据和全球地震台网(GSN，Global Seismographic Network)的宽频带地震数据，还运用了同震位移明显的两个台站的GPS位移数据，在最初的反演结果上有了更新和改进；第二，由于此次尼泊尔地震的余震目录较少，而且余震发生的位置分布于尼泊尔和中国国界线附近，对地震目录的统计可能不是很完备，同时也未对这些余震进行重定位，所以余震的分布可能会不够准确；第三，尽管发震断层是低角度逆冲兼走滑，但是在处理过程中，将断层和余震投影到地表可能会带来一定的误差.尽管有上述因素，但利用同震滑动位移计算断层面上的两分量应力变化仍然是一种有效的手段.

通过计算2015年尼泊尔M_W7.9地震的应力状态分布，加深了我们对地震破裂过程的理解.在主断层面上的静态应力变化与滑动位移分布图案有一些相似之处：滑动位移分布和应力变化分布都相对集中，且滑动位移较大的地区，应力下降也快.结合2015年尼泊尔M_W7.9地震的余震目录，我们统计了12 h内、24 h内、48 h内和48 h到5月12日M_W7.2地震之前四个阶段内M>3.0的余震分布情况.我们发现，大多数余震分布在远离主震的断层面的另一端，在主要滑动方向(倾向方向上)大约70%的余震分布在应力变化为正的区域，也就是主震同震位移造成的区域静态应力的增加导致了余震的发生.此外，在应力降较大的区域，也发生了少许的余震，但这些余震大多发生在应力变化梯度较大的地区，忽略滑动位移反演结果的不唯一性，应力变化梯度大的地区可能在主震破裂过程中没有完全破裂.为了得到一个更符合实际的滑动模型来解释余震的触发机制，我们利用张勇等(2015)的模型计算了波数域中滑动位移和应力状态的傅里叶谱.滑动位移和应力状态近似满足k^-3和k^-2衰减.在模型中，有效半径r_e约等于0.7a，这说明圆盘模型中应力增加的部分应该占整个圆盘破裂面积的51%.将本次尼泊尔M_W7.9地震的应力状态分布与滑动位移分布进行对比，可以看出，应力状态为负的地方比滑动位移为正的区域明显缩小了.这说明有很大一部分区域是处在有滑动位移但应力却增加的区域，即地震可以发生在有滑动位移的区域.上述现象证实了有效半径r_e的存在以及滑动位移满足k^-3衰减模型的结论.该结论说明，非均匀应力降模型比均匀应力降模型更加接近真实的震源破裂过程.

我们感谢张勇老师和Kang Wang老师提供的滑动位移联合反演数据以及对本项研究给予的大力支持.感谢评审老师的耐心评阅.