地震同震破裂是发生在地球岩石圈内的天然力学实验，在断层未破裂区域以及下地壳、上地幔深度产生应力扰动，并激发可以被大地测量技术捕获的震后形变。空间大地测量技术GPS、InSAR可以观测到6级以上地震的震后形变^[1]。这些空间密集、时间连续的大地测量资料为研究震后形变机理、岩石圈流变结构提供了重要约束^[2-3]。震后形变的机理主要包括同震应力驱动的震后余滑、中下地壳上地幔黏弹性松弛和地壳浅部的孔隙弹性回跳^[4-5]，不同的震后形变同时发生，难以完全区分。

2015-04-25喜马拉雅弧中东段发生M_W7.9地震，地震自西向东沿着脆、韧转换边界的狭窄区域破裂，但未破裂至地表^[6-8]。同震应力触发了大量余震，包括当天发生的中国藏南定日M_S5.9地震，以及18 d后发生的M_W7.3强余震。地震发生后，不同研究团队采用GPS或InSAR资料研究了该地震的早期震后余滑，结果显示，余滑主要发生在主喜马拉雅逆冲断裂同震区域的下倾区域^[9-10]，并且分布广泛，30 km深度仍大量存在。由于缺乏中国藏南震后形变资料的约束，且没有考虑震后黏弹性松弛作用的影响，这些结果高估了震后余滑的释放率。

本文首先介绍GPS观测的尼泊尔地震震后1 a的三维形变场，然后分析是否扣除黏弹性松弛效应对震后余滑反演结果的影响，并探讨断层浅部未破裂区域的地震危险性。

1 震后GPS观测资料处理 1.1 数据处理

采用GAMIT/GLOBK软件对尼泊尔境内和中国藏南的连续GPS站进行统一处理。采用IGS发布的精密星历，估计包括测站坐标、卫星轨道、大气延迟等参数，得到无基准的自由网解。将区域单日松弛解与全球单日松弛解(SOPAC)绑定，并通过坐标变化将解算结果归算至ITRF2008参考框架，由此得到原始坐标时间序列：

$ D = {V_t} + {d_{{\rm{period}}}} + {d_{{\rm{offset}}}} + {d_{{\rm{postseismic}}}} $

(1)

式中，V_t为测站在ITRF2008框架下的长期运动，d_period为季节性周年半周年变化，d_offset为地震同震形变或其他因素引起的阶跃形变，d_postseismic为震后形变。为获得震后形变，需要对上述非震后形变信号予以改正。

由于部分GPS站震前没有观测，无法直接估计其长期运动速率，因而本文采用内插估计方法计算每个GPS站的线性运动速率。首先对研究区域内震前的连续站、区域站速度场进行融合，得到ITRF2008框架下的震前背景速度场；再采用应变分析方法估计每个测站的长期运动速率^[11]；最后从原始坐标序列中扣除线性速率的影响。

尼泊尔以及西藏南部是东亚地区陆地水负荷形变最为显著的区域之一，垂直方向的周年季节性变化幅度(波峰至波谷)达15~20 mm^[12-13]。本文采用2002~2015年JPL公布的GRACE二级月产品计算每个GPS观测站的地表垂直形变时间序列，拟合得到周年、半周年的振幅和相位，并对原始坐标时间序列予以改正。此外，还对原始GPS时间序列中同震阶跃(2015-05-12 M_W7.3余震)予以改正^[14]。

1.2 震后形变场

为将不同观测时段的GPS时序归算到震后1 a的累积形变，对改正后的时间序列采用对数函数进行拟合：

$ D = F\ln \left( {1 + \Delta t/\tau } \right) $

(2)

式中，D是观测位移(单位mm), F为震后形变幅度(单位mm)，Δt为震后逝去时间(单位d)，τ为震后松弛时间(单位d)。本文取τ=30 d拟合每个GPS站的观测分量，拟合的平均水平误差为3.1 mm，平均垂直误差为6.8 mm。与同震形变相比，震后1 a的GPS观测结果显示，震后形变主要集中在尼泊尔北部区域，靠近主前缘断裂带(MFT)的区域震后形变并不明显(见图 1，图中橙色五角星为USGS测定的尼泊尔地震主震和最强余震的震中；误差椭圆表示95%的置信区间；黑色实线为断层线^[15]，MFT为主前缘断裂)。

2 震后余滑模型

震后余滑是由于地震同震破裂对断层破裂面产生了应力扰动，从而引起同震破裂面扩展区域(向上、向下或双侧)的蠕滑位错。本文以GPS观测震后1 a累积的形变为约束，采用弹性位错理论反演震后余滑的空间分布。设置反演的目标函数为：

$ F\left( {s, \beta } \right) = {\left\| {\mathit{\boldsymbol{W}}\left( {Gs-d} \right)} \right\|^2} + {\beta ^2}{\left\| {{\nabla ^2}s} \right\|^2} $

(3)

式中，G为格林函数，取决于断层的几何产状与GPS测站的分布；s为滑动量；d为观测值；W为观测值协方差矩阵；∇为有限差分拉普拉斯算子；β为光滑因子，通过断层滑动模型粗糙度和数据拟合度曲线折中选取。

震后余滑模型的断层几何形状以同震破裂模型为基础^[9]，对断层的走向进行微调，使之与断层地表迹线一致。通过大量试算，最终选取断层的走向为289°，断层倾角为8°，断层顶部埋深取3 km。考虑到震后余滑通常发生在同震断层破裂面的延展区域，对断层的长度和宽度进行扩展，分别取260 km和200 km，并将其沿着断层走向和倾向均匀划分为20×20个子断层。由于尼泊尔地震机制为纯逆冲性质，格林函数计算时仅考虑倾滑分量。采用非负最小二乘算法对式(3)的震后余滑分布进行反演。

根据图 2曲线，最终取β=0.2以平衡反演的数据拟合度与断层滑动模型粗糙度，对应的观测值与模拟值的拟合误差为3.6 mm。

采用未经黏弹性松弛改正的GPS三维形变场反演的余滑模型(简称为“改正前的余滑模型”)见图 3(图中黑色等值线为同震破裂等值线，黑色实线为断层，MFT为主前缘断裂)。反演的最大滑动量为18 cm，平均余滑深度为20 km。震后余滑释放的地震矩为8.8×10¹⁹Nm，对应的矩震级为M_W7.3，占同震释放总能量的12%。图 3显示，震后余滑模型不但可以解释GPS观测的地表水平位移，而且还可以解释近场大幅度垂直位移。改正前的余滑模型显示，余滑主要分布在同震破裂的下倾延伸部分，分布范围广泛，甚至在模型断层底端亦有存在。这是由忽略黏弹性松弛效应引起的。

3 顾及黏弹性松弛的余滑模型

震后黏弹性松弛是震后形变中另一重要机理，是同震应力在下地壳上地幔扰动，使得深部低黏滞系数的物质流动引起的变形。黏弹性松弛效应与震后余滑相互耦合。在震后形变建模时，如果忽略黏弹性松弛效应的影响，通常会高估断层深部的震后余滑^{[5, 15-16]}。

尼泊尔-喜马拉雅地区的大地电磁、地震波层析成像、地热结构、有效岩石圈厚度等相关研究揭示了喜马拉雅两侧显著的结构差异^[17-20]。在此基础上，我们构建该地区的岩石圈结构(图 4)。根据青藏高原的有效弹性层厚度，选取青藏高原的上地壳厚度为35 km，地壳厚度取莫霍面平均深度70 km，之下为地幔。依据印度大陆板块的有效弹性厚度^[20]，取70 km作为弹性层与黏弹性层的分界面。青藏高原-印度大陆板块岩石圈的边界位置则参考大地电磁探测结果及地热结构^{[17, 19]}。青藏高原岩石圈黏滞系数参考青藏高原其他地区(龙门山地区、可可西里地区)的有效黏滞系数数值范围^{[16, 21]}。印度大陆板块上地幔黏滞系数参考同样具有克拉通属性的四川盆地的高黏滞系数^[16]。

采用图 4岩石圈流变结构以及谭凯等^[22]反演的同震破裂模型正演计算的黏弹性松弛效应可以解释藏南地区的GPS观测值，但不能解释近场的震后形变(图 5)，调整青藏高原下地壳黏滞系数亦不能同时拟合远场和近场观测的震后形变。因而我们认为，远场的震后形变是黏弹性松弛引起的，而近场则由震后余滑主导。首先采用侧向非均匀的岩石圈流变结构正演计算黏弹性松弛引起的理论地表形变，并从图 1所示的震后形变场中予以扣除，最后采用相同的方法反演震后余滑分布(图 6)。

图 6为扣除黏弹性松弛效应后的震后余滑分布(简称为“改正后的余滑模型”)。改正后的震后余滑主要集中在同震破裂的下倾部分，在浅部同震未破裂区域极其有限。反演的最大滑动量为19.7 cm。震后无震滑移释放的地震矩为5.7×10¹⁹ Nm，对应的矩震级为M_W7.2，占同震释放总能量的8%。改正后的余滑模型在空间分布上更加集中，平均余滑深度由改正前的20 km变为16 km。这种余滑分布特征与应力驱动的余滑分布特征更为接近^[10]。扣除黏弹性松弛后，GPS观测残差与余滑模型模拟结果的拟合误差减小为3.4 mm。

4 结语

1) 改正前的震后余滑主要集中在同震断层破裂面的下倾延伸段，分布范围比较大，平均余滑深度为20 km。震后余滑释放的地震矩为8.8×10¹⁹ Nm，对应的矩震级为M_W7.3，占同震释放总能量的12%。

2) 扣除黏弹性松弛效应后反演的震后余滑分布更加集中，平均余滑深度减小为16 km，震后1 a释放的地震矩为5.7×10¹⁹ Nm，等效于M_W7.2地震，占同震释放总能量的8%。顾及黏弹性松弛效应的震后余滑分布与应力驱动的余滑分布比较接近。

3) 同震未破裂的断层浅部受同震影响而库仑应力增加。震后该区域未观测到明显的震后余滑，说明断层浅部区域仍处于闭锁状态，地震危险性值得密切关注。