地球物理学报  2016, Vol. 59 Issue (6): 2080-2093   PDF    
引用本文
谭凯, 赵斌, 张彩红 等. 2016. GPS和InSAR同震形变约束的尼泊尔MW7.9和MW7.3地震破裂滑动分布. 地球物理学报, 59(6): 2080-2093, doi: 10.6038/cjg20160614.  
TAN K, ZHAO B, ZHANG C H, et al. 2016. Rupture models of the Nepal MW7.9 earthquake and MW7.3 aftrershock constrained by GPS and InSAR coseismic deformations. Chinese J. Geophys. (in Chinese), 59(6): 2080-2093, doi: 10.6038/cjg20160614.  
GPS和InSAR同震形变约束的尼泊尔MW7.9和MW7.3地震破裂滑动分布
谭凯1 , 赵斌1 , 张彩红1 , 杜瑞林1 , 王琪2 , 黄勇1 , 张锐3 , 乔学军1     
1. 中国地震局地震研究所, 地震大地测量重点实验室, 武汉 430071;
2. 中国地质大学(武汉)地球物理与空间信息学院, 行星科学研究所, 武汉 430074;
3. 地壳运动监测工程研究中心, 北京 100036
收稿日期 2016-01-18, 2016-06-14 收修定稿
基金项目: 地震行业科研专项(201308009),中国地震局地震研究所所长基金(IS201506220),国家自然科学基金(40974012,41274027)联合资助.
作者简介: 谭凯,博士,研究员,主要从事大地测量与地球动力学研究.E-mail:whgpstan@163.com
通讯作者: 王琪,博士,教授,主要从事大地测量与地球动力学方面的研究.E-mail:wangqi@cug.edu.cn
摘要: 2015年4月25日尼泊尔爆发MW7.9地震,继而引发5月12日MW7.3级余震,GPS、InSAR监测到震源区及周边大范围同震形变.本文以国内外的GPS和InSAR同震形变为约束,考虑喜马拉雅断裂带岩石圈垂向分层和横向差异的影响,反演主喜马拉雅逆冲断裂在这次主震和余震中破裂面形状和滑动分布.结果显示,主震从USGS确定的震中位置向东偏南延伸100 km以上,破裂地面迹线与主前缘逆冲断裂迹线基本一致.破裂面倾角约7°~11°,大部分破裂集中在深度8~20 km,同余震分布深度一致.主震最大滑动量约6.0~6.6 m,位于14 km深处.余震破裂集中在震中附近30 km范围内,填补了主震东部破裂空区,最大滑动约3.6~4.6 m,位于13 km深.深度20 km以下基本没有破裂.地壳介质不均匀性对破裂滑动分布的影响较大,介质不均匀模型的观测值不符值比各向同性弹性半空间模型降低10%以上.本文地震破裂模型特征与地震反射剖面、以及根据震间期大地测量数据反演的喜马拉雅深部蠕滑剖面极其相似.跨喜马拉雅断裂剖面的震间形变量与地震破裂滑移量直接相关.以此推算,尼泊尔中部大震原地复发周期在300年以上.
关键词: 尼泊尔地震      GPS      InSAR      同震形变      破裂滑动分布     
Rupture models of the Nepal MW7.9 earthquake and MW7.3 aftrershock constrained by GPS and InSAR coseismic deformations
TAN Kai1, ZHAO Bin1, ZHANG Cai-Hong1, DU Rui-Lin1, WANG Qi2, HUANG Yong1, ZHANG Rui3, QIAO Xue-Jun1     
1. Key Laboratory of Earthquake Geodesy, Institute of Seismology, Chinese Earthquake Administration, Wuhan 430071, China;
2. Institute of Geophysics & Geomatics, China University of Geosciences, Wuhan 430074, China;
3. National Earthquake Infrastructure Service, Beijing 100036, China
Abstract: During the Nepal MW7.9 earthquake on 25 April 2015 and the MW7.3 aftershock on 12 May, large surface deformations were recorded by seismographs and geodetic instruments. The inversion of teleseismic waveforms enabled a quick gain of knowledge about spatio-temporal processes of rupture including the onset, rise time, speed of rupture, as well as the location, geometry and seismic moment, providing relatively imprecise solutions of distribution of slip. InSAR and GPS data collected in the near field can be used to constrain fault geometry and slip distribution at a finer resolution. However, previous geodetic models usually provided distributions of main shock slip in either an elastic half space or a spherically symmetric Earth. These solutions varied significantly in slip pattern. It is assumed that radial stratification and lateral heterogeneity of Earth medium might introduce considerable biases onto the geodetic inversions of slip distribution. This paper reanalyzes all available measurements of ground deformation associated with the main shock and its largest aftershock based on a radially layered and laterally heterogeneous Earth model. We used GPS-derived co-seismic displacements at a total of 33 sites in Nepal and south Tibet, together with two swaths of the ALOS-2 InSAR interferograms for constraining fault geometry and slip distribution. Firstly, using a genetic algorithm, we solved for rupture parameters (centroid location, depth, length, width, strike, dip, rake and slip magnitude) for a uniform slip model and estimated their statistic confidence intervals. Then, we divided the best-fitting fault plane further into numerous sub-fault patches and estimated slip values of these patches adopting a nonnegative least squares algorithm, which revealed in detail the heterogeneous features of slip on this optimal rupture plane. In the inversions of slip distribution, we minimized misfits to surface displacements while maintaining smoothness of slip across neighboring patches based on a trade-off curve of data misfits vs. slip roughness. The modeled surface displacements associated with the main shock are all directed to the epicenter, consistent with a simple thrusting mechanism. The co-seismic horizontal offset was 1.89 m at the KKN4 GPS station in Nepal where the maximum observed displacement was documented. In south Tibet, 100~400 km away from the epicenter, the modeled co-seismic horizontal displacements varied from several millimeters to approximately half a meter, consistent with the largest horizontal offset of 54 cm recorded by the J041 station. The modeled offsets of the largest aftershock match the observed ones up to 2 cm at GPS sites as far as approximately 200 km away from the epicenter and InSAR line-of-sight (LOS) range offsets of up to 1 m are fitted within their formal errors. Under the same constraint conditions, the slip amount inverted in the elastic half space (the homogeneous model) is usually larger than that for a layered crust (the inhomogeneous model), although they appear to be little changes in the distribution of slip. Furthermore, the homogeneous model yields an overall misfit to the data better that the inhomogeneous model does. The former misfits are estimated only 67% for GPS, and 84% for InSAR of the latter ones. Compared with the homogeneous model, the inhomogeneous model reduces, in particular, the postfit residuals by more than 10% of the GPS displacements on the hanging wall north of Kathmandu. These GPS sites including those in Tibet are sensitive to the effect of laterally inhomogeneous media on slip behavior. The modeled fault geometry shows a reverse fault striking 285° from the north, dipping 7°~11° to the north, and extending east-southeasterly over 100 kilometers from the USGS epicenter. The geometry of the Main Himalayan Thrust inferred from the coseismic deformation is similar to that as illustrated by a seismic reflection profile obtained in the INDEPTH project and is in agreement with what is adopted for a creeping segment of the Main Himalayan Thrust beneath the Himalaya in fitting GPS velocities. The modeling of coseismic and interseismic deformation shows that the India plate plunges into the Tibetan Plateau at a sub-horizontal dip angle. Our model of slip distribution of the main shock reveals a main patch with local peak slip about 6.0~6.6 m at 14 km depth. The main shock geodetic moment is estimated to be 7.65×1020 N·m, corresponding to the moment magnitude MW7.9, in accordance with the GCMT solution. Most of slip induced by the largest aftershock is localized within a circular region of 30 kilometers in radius around its epicenter, filling in the gap (area of no any slip) left by the main shock. This aftershock is characterized by a main patch with local peak slip of about 3.6~4.6 m at 13 km depth and its geodetic moment is estimated to be 1.12×1020 N·m, corresponding to MW7.3. In the 2015 event, most of slip is confined to the brittle crust at depths of 8 to 20 km, and slip below 20 km depth is barely observed inconsistent with the depth distribution of aftershocks. We analyzed the stress loading triggered by the mainshock according to our slip model. The calculation of Coulomb stress distribution indicates that the majority of the aftershocks is restricted to the areas that have experienced an enhanced stress loading of 0.2~0.8 MPa due to the main shock. Our model suggests that the 2015 earthquake in Nepal has released stresses accumulated in interseismic period due to the underthrusting of the India plate, and hypothesizes that the recurrence of such a large earthquake is more than 300 years in central Nepal where the slip rate of the seismogenic thrust fault is about 20 mm·a-1..
Key words: Nepal earthquake      GPS      InSAR      Coseismic deformation      Rupture slip distribution     
1 引言

2015年4月25日，尼泊尔发生MS8.1级地震(http://www.cea.gov.cn/publish/dizhenj/468/553/101803/index.htm )．5月12 日，在主震东偏南方向发生MW7.3余震．重烈度区从震中向东延伸，不少房屋倒塌，公路、通信等基础设施损坏．受灾范围包括尼泊尔、印度北部、巴基斯坦、孟加拉和中国藏南等地区．美国地质调查局(USGS)给出主震震源机制解的地震位置在84.73°E，28.23°N，深度10 km，走向295°，倾角10°( http://earthquake.usgs.gov/earthquakes/search/)．Global CMT目录开始给出MW7.9地震位置为85.37°E，27.77°N，深12 km，走向293°，倾角7°；在2015年底修改为走向287°，倾角6°(http://www.globalcmt.org/CMTsearch.html )．美国喷气推进实验室公布了主震和余震在尼泊尔境内的GPS同震形变(http://aria-share.jpl.nasa.gov/events/20150425-Nepal_EQ/GPS/ )，加州大学圣迭戈分校也发表了主震和余震的InSAR同震形变结果(Lindsey et al.，2015)．

张勇等(2015)基于USGS的震源机制解，主要利用地震波数据反演MW7.9地震破裂过程，破裂以逆冲为主，最大滑动量在3.2～5.2 m间．地震波数据一般位于远场或者中远场，以此约束获得的地震破裂位置、几何形状和滑动分布的空间分辨率较低，破裂子断层长度一般在十几到二十公里．单新建等(2015)基于InSAR和尼泊尔境内GPS数据，假设弹性半空间的地球模型，获得MW7.9主震更为精细的破裂滑动分布.这些模型一般仅给出MW7.9主震滑动分布结果，缺少最大余震的破裂分布信息，因此此次地震完整的滑动分布尚有待约束．另外，喜马拉雅断裂带岩石圈垂向分层和横向物质差异显著，可能会对研究结果产生显著影响(Xu and Xu，2015)．

我们基于“中国大陆构造环境监测网络”(简称陆态网络)计算获得了尼泊尔震中周围500km范围内GPS连续观测站和流动观测站主震和余震的静态同震形变．本文以国内外密集的GPS和InSAR同震形变数据为约束，考虑喜马拉雅断裂带岩石圈垂向分层和横向物质差异的影响(Wang et al.，2003)，反演子断层长度约3 km的精细破裂模型，对尼泊尔MW7.9主震和MW7.3余震断层几何形状、运动学特征进行更为准确的约束，以了解尼泊尔震区最后静态滑动的总体分布，为青藏高原动力学研究和危险性评估提供参考．

2 构造地震背景与尼泊尔地震形变资料

尼泊尔位于喜马拉雅山中段南麓，构造上系欧亚板块与印度板块的交界区，分布有喜马拉雅主前缘逆冲断裂、主边界逆冲断裂、主中央逆冲断裂和西藏南部拆离断裂系．GPS观测表明，印度—欧亚板块以36~37 mm·a-1的速率汇聚(Bilham et al.，1997； Ader et al.，2012；Wang et al.，2014)，印度次大陆沿北北东方向下插挤入青藏高原(Zhao et al.，2015c；Nelson et al.，1996；Hauck et al.，1998)，尼泊尔中部跨喜马拉雅山的地壳缩短汇聚速率约为20 mm·a-1(Bilham et al.，1997；Ader et al.，2012)，吸收了约一半的板间汇聚速率．

强烈的构造运动孕育了频繁的地震活动．地中海－喜马拉雅地震带释放地震能量占全球地震释放总能量的24%，过去的一个多世纪曾发生多起8级左右大地震(Bilham et al.，2001；Kumar et al.，2010)．在尼泊尔，西部在1505年曾发生MW8级以上地震，东部在1833年发生MW7.8级地震，1934年又发生MW8.2级地震(Bilham et al.，2001)．据全球地震目录，1976年以来尼泊尔中部地区长期缺失强震．2015年MW7.9级地震爆发在1833年和1934年大震的西侧，填补了前两次大震西侧100多公里的地震空区．

我们收集了“陆态网络”全国260个GPS基准站和西藏17个GPS区域站数据，其中区域站进行了多次观测(李强等，2012)，4月25日地震后“陆态网络”立即进行应急复测.UNAVCO公布了尼泊尔境内16个GPS连续站数据．我们将这些GPS数据综合处理(赵斌等，2015a)，给出主震和最大余震同震形变场(图 1)，并估算其方差－协方差矩阵．其中连续观测GPS站同震形变中误差一般在2 mm以内，流动观测站点误差在4 mm以内．

图 1主图部分显示MW7.9主震引起的33个GPS站点同震水平位移，符合逆冲断层破裂的形变特征，位移矢量从南北两侧指向震中．GPS站最大位移1.89 m发生在尼泊尔境内距主震GCMT震中仅9.5 km的KKN4站．尼泊尔境内距离震中约50 km的两个测站观测到1.3～1.5 m的同震位移．在藏南地区，GPS最大位移54 cm处在毗邻尼泊尔的聂拉木县J041区域站．距离震中约140 km的吉隆县J040和J339两个区域站位移约13～15 cm．离开震中约243 km的仲巴和昂仁位移约2 cm．距离震中400 km以外的其他测站同震形变一般小于5 mm．图 1内插图则显示尼泊尔和藏南地区GPS站(约25点)MW7.3级余震同震位移，在震中北东方向至我国珠峰测站之间产生2 cm以上的变形，影响范围在200 km左右.

图 1 尼泊尔地震构造背景与2015年MW7.9和MW7.3地震同震形变场 沙滩球：2015年尼泊尔地震序列的震源机制解； 灰色五角星：历史大震； 灰色点：历史地震； 箭头：GPS同震位移矢量； 彩色点：InSAR降采样点的LOS形变值； 红色矩形框：MW7.3余震InSAR覆盖区域； 左下角彩色点：余震InSAR降采样点的LOS形变值； 红细线：喜马拉雅主前缘逆冲断裂(Styron et al.，2011； Ader et al.，2012)； 黑细线：断裂； 左下角插图是MW7.3余震同震形变场． Fig. 1 The tectonic setting and coseismic displacements of the 2015 Nepal MW7.9 and MW7.3 earthquakes Beach balls: the focal mechanism solutions; gray star: historical great earthquakes; gray solid dots: historical seismicity; arrows: GPS displacements;colored dots: InSAR range offset; red rectangle: InSAR coverage for the MW7.3 aftershock; red curve: the Main Frontal Fault;black curve:fault. The left lower insets show geodetic data for the MW7.3 aftershock.

美国加州大学圣迭戈分校公布了4景日本ALOS-2卫星L波段InSAR主震同震形变结果，InSAR视线向形变值最大1 m左右．2015年2月22日—2015年5月3日升轨宽幅的InSAR资料(图 1)覆盖范围最广，本文用该幅InSAR资料去趋势项后的结果进行四叉树法降采样，获得4437个采样点．根据雷达入射角和卫星轨道方位角计算获得各观测点的视线向单位矢量大小范围是：[0.429792～0.749677，-0.12467～-0.084656，0.650093～0.898951]．将GPS三维形变投影到视线向上并与InSAR观测值比对，在尼泊尔境内两者差值在5 cm以内，与InSAR观测的标称误差一致(Qiao et al.，2015)；在藏南地区，特别是InSAR图幅的东北角，两者差值一般在7～10 cm，InSAR观测结果可能受到更多误差因素影响，或者受到后处理中去趋势项的消极影响．因此，将InSAR图幅东北角的一小块数据与其他区域的数据分开，在后面的反演中用不同的改正参数来估算其误差影响．

2015年5月3日—5月17日宽幅InSAR资料给出了余震同震形变．考虑到余震影响范围较小，离震中100 km外的InSAR形变可能误差所占比例较大，舍弃不用．降采样中，尽量选取近场样点，为削弱远场粗差的影响，远场以等间隔降采样．如此获取总共2405点InSAR观测值(图 1内插图)，各样点视线向单位矢量范围是：[0.495218～0.683239，-0.119019～-0.095208，0.72043～0.863536].InSAR与GPS计算的视线向LOS形变差值大部分在1～5 cm内，少部分在6～9 cm内．为避免观测数据不完全独立对反演结果带来误差，我们用Sun等(2008)的方法构建主震和余震InSAR观测值离散采样点的方差-协方差矩阵．

3 地震破裂模型反演

震源矩形断层位错引起的地表变形可用弹性位错模型计算，主要与断层几何特征7参数(长、宽、深、走向、倾向、水平坐标)和滑动参数(走滑量、倾滑量、张性分量)有关，以形变观测值为约束可以反演断层破裂参数．如果假设破裂由少数几个面积比较大的断裂组成，则反演的目的是求得破裂的大体几何形状和平均意义的滑动量．如果将破裂细分为更多的子断层，则反演的目的是为了获得精细的破裂滑动分布．两种反演方法都应该满足观测数据拟合度和滑动分布粗糙度最小，即

(1)

式中d是形变观测值，包括GPS位移(本文采用精度较高的水平位移)及InSAR形变观测值．W是观测值的权矩阵，与观测值方差-协方差矩阵D的关系为D=W－1W．G是形变格林函数，本文主要基于半无限空间地壳分层的弹性位错模型(Wang et al，2003)计算格林函数．s是待求的未知数矢量，一般指子断层滑动矢量．但是在InSAR数据参与的反演中，为消除轨道误差的影响，每个数据块应该引入3个未知改正参数：2个与距离成正比的双线性平面改正参数，1个相位常数差改正参数(Wang et al.，2011)．在MW7.9主震反演中，我们根据InSAR与GPS LOS差值的分区特征，将InSAR数据分为2块，所以引入6个InSAR未知改正参数，与子断裂滑动参数同时解算．MW7.3余震InSAR重采样区域较小，反演中只引入3个未知改正参数．β是子断裂滑动量平滑因子，单位为m－1，平衡了不符值与粗糙度的权重．根据观测值不符值和粗糙度的折衷(trade-off)曲线来选取平滑因子，平滑因子选取的原则是不能不恰当地增加观测值不符值．L是拉普拉斯二阶差分算子，沿走向和倾向的第(i，j)个子断层滑动量s的差分公式为

(2)

其中Δx、Δy分别是沿走向和倾向的相邻子断层距离(可以设为1)．对于破裂边缘的子断层，采用拉普拉斯一阶差分算子，例如对于最上一层的第(i，j)子断层，沿倾向(2式右边第二项)只能做一阶差分： ．

我们考虑了地壳垂向分层对形变的影响，地壳分层采用根据地震波数据、冰层及沉积岩深度数据、人工地震及钻深数据构建的crust 2.0模型( http://igppweb.ucsd.edu/~gabi/rem.htmlhttp://igppweb.ucsd.edu/~gabi/rem.html)．由于同震形变测点主要分布于28°N两侧，南北两侧地壳物质差异显著，地壳厚度差异在15 km以上，因此以28°N为分界考虑南北两侧地壳横向差异的影响，28°N以南测点将使用南部地壳分层模型计算格林函数，28°N以北使用北部地壳分层模型(表 1)．

表 1 地壳结构模型 Table 1 Crustal structure models

地表位移是断层几何形状参数和滑动参数的非线性函数，可通过非线性优化方法求解．如果假定断层几何形状参数已知，则地表位移是滑动参数的线性函数，可用线性反演方法求解滑动参数，滑动分布反演结果主要取决于破裂面先验几何特征．尼泊尔地震破裂没有完全出露地表，断裂位置、走向、长度、倾角、宽度等具有较大的不确定性．本文先用非线性反演的遗传算法，确定断层几何形状参数，即将断裂用数个断裂段表示，用遗传算法确定破裂面的几何形状，获得合理的、最优的几何形状参数如断裂倾角、走向等；然后，在最佳几何形状模型的基础上，将破裂面离散化为更多的子断层，用线性反演方法，如非负最小二乘法确定最终的破裂滑动分布模型．

3.1 主震破裂面几何形状参数

本节用单个矩形断裂表征MW7.9主震破裂面的几何形状，用并行遗传算法求解其7个几何形状参数，对于每组候选模型，以主震同震形变为观测值约束，用非负最小二乘法求解单断裂均匀滑动量和InSAR改正参数．目标函数采用(1)式的观测数据拟合度和滑动分布粗糙度之和．

由于地质调查没有发现破裂出露地表，所以需要在较大的范围内确定MW7.9主震破裂几何参数．通过多次遗传算法试验，根据不符值和粗糙度的折衷曲线选取β=4，反演获得MW7.9主震单断层破裂最优模型(表 2)的破裂长度约80 km，宽度约25 km，以逆冲为主，兼有少量走滑．

表 2 MW7.9主震单断裂模型参数的搜索上下界、最优模型、置信区间 Table 2 Searching ranges，optimum models and confidence intervals of rupture parameters of the MW7.9 main shock

单断裂模型的倾角在7°～11°变化，最优值倾向于取较小的角度，这与GCMT震源机制解小倾角7°对应．如果以单断裂模型的位置和倾角扩展断裂至地表，则地表迹线位于主前缘逆冲断裂以南约30 km．但是余震目录显示，余震向南终止于主前缘逆冲断裂(参见图 8)，而没有出现在主前缘逆冲断裂的南部．因而我们判断，破裂面倾角可能较低，接近地表倾角逐渐变高，最后终止于主前缘逆冲断裂．这与地质剖面(Hauck et al.，1998；Nelson et al.，1996；Zhao et al.，1993)特征符合．为了使最优破裂模型的剖面同时与深部和浅部吻合，我们取平均倾角9°作为最优断层模型的破裂倾角．

单断裂模型走向一般在280°～290°间，与Styron等(2011)地质调查给出的喜马拉雅主前缘逆冲断裂尼泊尔中东段走向基本一致，并与InSAR形变升高区和降低区的分界线大致平行．

模拟结果显示，GPS观测值残差均值为6.0 mm，InSAR观测值残差均值为9 cm．我们认为单断裂模型反映了破裂基本特征，与地质资料和震源机制解比较吻合，可以做为后续研究的基础．由于走向活动范围较大，在后续滑动分布反演中，还需要进一步精化验证．

3.2 主震破裂滑动分布

由于走向变动范围稍大，本节使走向在283°～295°间依次变化1°，并假设前面反演所得的地表迹线、倾角等其他破裂几何形状参数为已知，构建13个单断裂破裂模型(即13个单断裂候选模型)．考虑破裂可能超出单断裂模型范围，新建模型的破裂面向两侧及深处延长并超过余震分布范围(长220 km，深30 km)．将13个单断裂候选模型都离散化为3 km×3 km的4672个矩形子断层，构建了13个离散化的断裂滑动分布候选模型，以获得更精细的滑动分布．如此则地表位移是子断层滑动参数的线性函数，遵循(1)式观测数据拟合度和滑移分布粗糙度最小的原则，采用非负最小二乘法解算所有子断层的走滑量、倾滑量和InSAR改正参数．在相同的平滑因子下(例如β为30或35)，反演获得13个破裂滑动分布候选模型的观测值不符值与走向的关系曲线(如图 2)．由于走向285°在上节的单断裂模型反演结果中出现概率较大、在本节的滑动分布反演中观测值不符值也较小，所以选取285°做为最终破裂模型的走向．

图 2 走向变化引起的观测值不符值曲线 Fig. 2 Observation misfits WRSS curve vs. strikes

采用走向285°、倾角9°的破裂几何形状模型，进一步确定主震最终的滑动分布．使平滑因子在5～150间变化，获得不符值和粗糙度的折衷曲线(图 3)．兼顾观测值不符值最小和破裂光滑原则，平滑因子可以在20～60之间选取，则最大滑动量置信区间为6.0～6.6 m．平滑因子对其他破裂特征量的影响较小，USGS震源位置的滑动量都在1.0 m以上，GCMT震源位置的滑动量在2.5 m以上．GPS后验中误差在4～8 mm，InSAR后验中误差为5～10 cm，与先验中误差一致．矩震级一般为MW7.86～7.87．因此，主震破裂模型反演的稳健性较好．

图 3 MW7.9主震观测值不符值和滑动粗糙度的折衷曲线图 曲线右侧的数值是平滑因子，纵坐标轴上的虚线标出了最小不符值min WRSS位置、以及1.5倍、2倍和3倍最小不符值位置. Fig. 3 Trade-off curve of observation misfits WRSS vs. slip roughness of the MW7.9 mainshock The smoothing factors are labelled along the curve，and the WRSS minimum and 1.5，2 and 3 times values are indicated with dashed lines.

选取折衷曲线拐点的平滑因子β=35，在此光滑约束条件下反演获得主震最终滑动分布如图 4a所示，以逆冲为主，破裂集中在由西向东的狭长的范围内，最大滑动峰值高达6.37 m，位于13.83 km深处．USGS起始破裂点的滑动量1.48 m，向东扩展，在GCMT震源处的滑动量为2.74 m，加德满都深部滑动1.48 m．最大滑动量距离加德满都29 km．依据该模型计算的标量地震矩Mo=7.65×1020N·m，矩震级MW=7.86，与GCMT震级一致．

图 4 MW7.9主震和MW7.3余震在弹性分层和横向差异地壳中的同震滑动分布模型 (a)MW7.9主震滑动分布；(b)MW7.3余震滑动分布；(c)总滑动分布.黑五角星：USGS震中； 红五角星：GCMT震中； 滑动量等值线：1～6 m； 小红方框：加德满都； 红曲线：喜马拉雅主前缘逆冲断裂(Styron et al.，2011). Fig. 4 Coseismic slip distributions in a multi-layered and laterally heterogeneous brittle crust associated with the MW7.9 main shock and MW7.3 aftershock respectively (a)The MW7.9 main shock;(b)The MW7.3 aftershock;(c)Combined slip. Black star: USGS epicenter; red star: GCMT epicenter; slip contour line:1~6 m; little red square: Kathmandu; red curve: the Main Frontal Fault(Styron et al.，2011).

平滑因子为35的最佳模型获得较好的拟合效果(图 5)．形变最大的GPS点KKN4(形变1.88 m，残差0.8 cm)、NAST(形变1.33 m，残差2.3 cm)、CHLM(形变1.42 m，残差3.8 cm)的模拟残差都在观测值的3%以内．中尼边境两个站J040(形变0.16 m)、J041(形变0.54 m)的拟合残差也在3 cm左右，其他测站残差一般在2 cm以内．拟合的InSAR LOS值平均残差为5 cm，测区中部的残差一般在10 cm以内，震中附近部分点残差接近20 cm，在南北两侧的远场部分点残差接近20 cm．GPS和InSAR后验中误差与先验误差一致．

图 5 GPS和InSAR观测值模拟残差 黑箭头：MW7.9 GPS形变拟合残差； 红箭头：MW7.3 GPS形变拟合残差；颜色点：InSAR LOS形变值拟合残差；
红色矩形框：MW7.3地震InSAR观测区；左下角小图：MW7.3观测值残差. Fig. 5 Postfit residuals of GPS displacements and InSAR range offsets Arrows: GPS residuals; Coloured dots: InSAR residuals; Rectangle: InSAR coverage for the MW7.3 aftershock;
The main figure shows the mainshock and the left lower inset displays the MW7.3 aftershock.
3.3 余震破裂滑动分布

由于MW7.3余震位于主震破裂传播方向上，两者的震源机制很接近，并且Styron等(2011)给出的断裂地表迹线也很平直，为了简便，假设MW7.3余震沿着主震破裂面扩展，两者具有相同产状．使平滑因子在5～150间变化，以余震同震形变为约束，遵循(1)式观测数据拟合度和滑移分布粗糙度最小的原则，采用非负最小二乘法解算所有子断层的走滑量、倾滑量和InSAR改正参数．获得观测值不符值与粗糙度的折衷曲线(图 6)，平滑因子可以在22～45之间选取，则最大滑动量置信区间为3.6～4.6 m．其拐点位于平滑因子35附近，取β=35，获得对应的最佳滑动分布(见图 4b)．余震滑动以逆冲为主，集中在震中30 km范围内，从较深的USGS震源向较浅的GCMT震源扩展，最大滑动峰值3.97 m，位于13.4 km深．依据该模型计算的标量地震矩Mo=1.12×1020N·m，矩震级MW=7.3，与GCMT震级一致．余震形变拟合效果见图 5左下角小图．远场GPS点符合很好，只有近场两点残差较大，不到1 cm．拟合的InSAR观测值残差都在10 cm以内，平均残差为3.7 cm．GPS和InSAR后验中误差与设置的先验误差一致．

图 6 MW7.3余震观测值不符值和滑动粗糙度的折衷曲线 曲线右侧的数值是平滑因子，纵坐标轴上的虚线标出了最小不符值min WRSS位置、以及1.5倍、2倍和3倍最小不符值位置. Fig. 6 Trade-off curve of observation misfit WRSS vs. slip roughness of the MW7.3 aftershock The smoothing factors are labelled along the curve and the WRSS minimum and 1.5，2 and 3 times values are indicated by dashed lines.

余震破裂填补了主震东部靠近中国边境的一块破裂空区，两者的破裂区和空区具有很好的互补性．最大总滑动量为6.6 m，与主震最大滑动量基本一致．

3.4 地壳垂向分层、横向差异与各向同性弹性半空间对反演结果的不同影响

我们还用各向同性弹性半空间模型(Okada，1985)反演了MW7.9主震的滑动分布(图 7)．介质不均匀模型(同时考虑地壳垂直分层和横向差异)比弹性半空间模型更接近震区实际情况，在相同的约束条件下，观测值拟合度提高．介质不均匀模型的GPS不符值只有弹性半空间模型的67%，InSAR不符值是84%．弹性半空间模型设置的剪切模量一般与上地壳剪切模量相当(见表 1)，比浅部沉积层剪切模量大，比深部中下地壳剪切模量小，因此弹性半空间模型反演获得的浅部滑动比不均匀介质模型的稍大，而对滑动分布形态影响很小，两者总体分布形态很相似．考虑地壳垂直分层和横向差异影响后，GPS观测值拟合度改善很大，只有边境附近的J041拟合度提高较少．可能是J041位于InSAR监测区，附近的InSAR点较多且误差较大，使得J041对反演结果的贡献降低．如果要继续提高GPS拟合程度，可能需要用GPS对InSAR观测数据做进一步的校正，或者在反演中根据InSAR误差分区特征分几个不同区域，加入InSAR改正参数一起反演．如果只考虑地壳垂直分层而忽略横向差异影响，则28°N南侧地壳模型的拟合效果比28°N北侧模型的好，其不符值只有北侧模型的75%～90%，可能与震中和大部分测点位于28°N南侧有关．考虑介质横向差异后，28°N北侧观测值拟合程度有所提高．不过需要注意的是，本文以28°N区分南北两侧地壳模型是一种近似方法，可能会给计算结果带来一定的近似误差．

图 7 MW7.9主震在弹性半空间的滑动分布模型 Fig. 7 Coseismic slip distribution model in an elastic half-space due to the MW7.9 mainshock Black star:USGS epicenter; red star:GCMT epicenter; black contours:1~6 m slip; red square: Kathmandu ; red curve: the Main Frontal Fault(Styron et al.，2011).
4 讨论

尼泊尔地震发生后，USGS利用远场P波、SH波和长周期面波反演了地震破裂过程，最大滑移量为3.0～4.0 m(表 3).张勇等(2015)以地震波资料约束反演主震，获得其最大滑移量为3.2 m．随着近场GPS和InSAR形变资料的加入，破裂滑移量得到较好约束．苏小宁等(2015)以GPS为约束，反演获得6.8 m最大滑移量．JPL利用远震波形、GPS和InSAR数据联合，其最大滑移量超过6.0 m(http://aria.jpl.nasa.gov/node/43 )．这些模型反演一般以震源机制解(位置、走向、倾角)为先验模型，由于地震破裂没有出露地表，所以先验模型的位置、走向和倾角具有较大的不确定性．

表 3 尼泊尔地震破裂参数 Table 3 The Nepal earthquake focal parameters

单新建等(2015)认为主边界断裂是此次地震地表迹线，而图 8、图 9展示，余震分布往南终止于Styron等(2011)给出的主前缘断裂.另外，喜马拉雅断裂带岩石圈垂向分层和横向物质差异显著，也对研究结果产生一定影响．目前公布的模型一般较粗，子断层长10～20 km、深3～4 km，对近场地震应变应力变化和地震危险性评估的准确程度有一定的影响．

图 8 尼泊尔地震破裂剖面、库仑应力变化、INDEPTH剖面 (a)垂直于断裂方向的震间GPS水平速度剖面．红色小圈：垂直于断裂方向的GPS点水平速度；蓝色曲线:震间速度模型曲线(赵斌等，2015b)．(b)垂直于断裂方向的破裂剖面.MFT：主前缘逆冲断裂；MBT：主边界逆冲断裂；MCT:主中央逆冲断裂；STD：藏南拆离层；从MFT往北延伸的黑线：本文断裂剖面；其下方和左侧的曲线：沿剖面的主震和余震总的最大滑动量；下方绿色的五角星：USGS震中；上方黄色的五角星：GCMT震中；下方品红色五角星：USGS定位的MW7.3余震；库仑应力是MW7.9主震引起的库仑应力，计算深度:10 km，接收断层走向285°，倾角9°，滑动角90°，摩擦系数0.8. Fig. 8 Slip distribution，Coulomb stress and aftershocks distribution of the Nepal earthquake across the profile (a)Interseismic GPS horizontal velocity profile across the Himalaya.Red circle: GPS horizontal velocity across the Himalaya; blue line: interseismic dislocation model.(b)Rupture along a profile across the Himalaya. MFT: Main Frontal Thrust; MBT: Main Boundary Thrust; MCT: Main Central.Thrust; STD:Southern Tibet Detachment; thick black line extending from MFT: the rupture along a profile inferred from this paper; the bottom curve and the left curve: the maximum of total slip; green star: USGS epicenter; yellow star:GCMT epicentre; magenta star:USGS epicentre of the MW7.3 aftershock; the Coulomb stress were caused by the MW7.9 main shock at depth of 10 km，the receive fault is assumed striking 285°，diping 9° with a rake of 90°．The friction coefficient 0.8 is used.
图 9 尼泊尔MW7.9主震库仑应力与余震分布(库仑应力接收断层走向285°) 库仑应力是MW7.9主震库仑应力，计算深度10 km; 接收断层走向285°，倾角9°，滑动角90°，摩擦系数0.8； 灰色圆点：所有3级以上的余震； 虚线曲线：主震滑动量为1～6 m的等值线； Fig. 9 The MW7.9 mainshock Coulomb stress and the aftershocks distribution The Coulomb stress were caused by MW7.9 main shock at depth 10 km，the receive fault striking 285°，diping 9°，rake 90°，friction coefficient 0.8; gray dots: >M3 aftershocks; dashed curve: 1~6 m slip contour lines of the main shock; black plus:USGS main shock epicentre，GCMT main shock epicentre; USGS aftershock epicentre.

此外，这些破裂滑动模型一般反映MW7.9主震滑动，缺少余震滑移的反演结果(Avouac et al.，2015；Galetzka et al.，2015；Lindsey et al.，2015；张勇等，2015；单新建等，2015；刘刚等，2015；苏小宁等，2015)．

我们充分发挥近场InSAR和GPS资料高精度、广覆盖作用，首先直接约束破裂位置、走向、倾角，然后将断裂划分为3×3 km尺度的子断层，反演获得更精细的破裂滑动模型．近场形变资料约束反演获得的断层走向约为280°～290°，与InSAR观测值升高区与沉降区的分界线平行，与地质资料给出的尼泊尔中部喜马拉雅逆冲断裂(Styron et al.，2011)走向基本一致.破裂倾角约7°～11°，小倾角模型比大倾角模型的拟合度更好．

本文主震最大滑动量大约位于USGS初始破裂点和GCMT矩震中的连线上，靠近GCMT矩震中．余震最大滑动量介于USGS初始破裂点与GCMT矩震中位置之间，离GCMT矩震中较近．GCMT是利用全球数字地震台网(GDSN)体波资料，包括多重反射波和P、S波的转换波，使用的是地幔波的资料，其矩心矩张量解可以理解为在更平均化尺度下的震源体的某种“核心”(高原等，1997)．地表静态形变反演的最大滑动区域可能与GCMT震源的“平均核心”含义更接近．

破裂模型剖面与INDEPTH剖面(图 8)在几何特征上具有一致性．INDEPTH揭示了喜马拉雅山和青藏高原南部的地壳精细结构和深部构造(Hauck et al.，1998；Nelson et al.，1996；Zhao et al.，1993)，沿主喜马拉雅逆冲断裂MHT发育5～ 6 km厚度的低速层，推断是印度大陆上地壳叠加增厚的主滑移面，与本文的地震破裂面相对应．主前缘逆冲断裂MFT、主边界逆冲断裂MBT以及主中央逆冲断裂MCT向深部延伸成为铲式断裂，汇聚于主喜马拉雅逆冲断裂MHT．在地震反射剖面上，MHT表现为一条强地震反射带，它从喜马拉雅山脊下26 km深度(相对于印度大陆20 km深，与本文破裂模型的深度一致)以低缓倾角9°～ 10°向北延伸到雅鲁藏布缝合带下面42km深度，总长超过150 km．

沿剖面的地震破裂分段特征与震间形变积累分段特征具有较好对应关系．在图 8，沿倾向分别给出了随水平向和深度变化的最大滑动量曲线．最大滑动量在水平距离和深度上都对应于震中位置，2 m以上的滑动量一般对应于水平距离60～130 km处(深度8～20 km)．水平距离0～60 km(深度0～8 km)，水平距离130～190 km(深度20～30 km)的滑动量都小于2 m．震间速度剖面(Bilham et al.，1997)显示跨喜马拉雅山速率变化20 mm·a-1，水平距离0～60 km的速率变化2～4 mm·a-1，130～240 km的速率变化只有3～5 mm·a-1，大部分震间能量积累在中间段60～130 km，速率变化为9～12 mm·a-1，这一段也是地震破裂滑动量大于2 m的区间，对应深度8～20 km．震间期大地测量结果反演获得的喜马拉雅逆冲断裂深部滑动倾角约10°，浅部闭锁深度大于20 km(Ader et al.，2012；赵斌等，2015b)，与本文反演的破裂模型相符．

喜马拉雅山地震主要受控于印度—欧亚板块碰撞边界断裂闭锁段应变积累.地震释放的矩张量与破裂时的滑移量及面积有关(Feldl and Bilham，2006)．根据震间地壳缩短率、地震破裂最大滑移量、以及地震破裂长度，可以大致推算原地大震复发周期．如果考虑尼泊尔在经度84.5°～86°的至少100 km区段，以破裂最大滑动量6.6 m，按20 mm·a-1速率积累估算，则需要330 a原地复发一个MW7.9大震．

破裂滑动可引起区域库仑应力变化，进而影响地震活动(盛书中等，2015)，因而破裂滑动分布结果可与余震活动分布相互验证．为此，根据尼泊尔MW7.9主震最优破裂模型计算的静态库仑应力(Toda et al.，2005)分布，将静态库仑应力、破裂滑动分布、余震分布分别投影到垂直剖面(图 8)和水平面(图 9)上．在图 8为了与震间形变比对，展示了主震和余震总的滑动量，而图 9只展示主震滑动量．库仑应力变化具有三维空间的变化，鉴于篇幅关系，本文只展示出在深度10 km库仑应力水平变化，以及在断裂中线上的剖面变化．从图 8可知，大部分余震分布在8～20 km 深度内，与本文模型总滑动大于2 m 的区域一致．图 8和图 9都显示，库仑应力变化与余震分布有很好的一致性．比如在图 9，余震主要分布在高滑动区外部主震库仑应力升高的地方，并且MW7.3余震诱发的余震与其库仑应力相关联．因而，可以认为本文滑动分布模型与余震活动相一致，余震分布基本都落在最大滑动量等值线外部库仑应力增加的区域，与2008年汶川地震和2013年芦山地震中破裂与余震活动特征类似(Wang et al.，2011；谭凯等，2015)．

5 结论

本文顾及喜马拉雅地区岩石圈垂向分层和横向介质差异，联合反演InSAR和GPS同震位移，约束尼泊尔MW7.9主震和MW7.3余震破裂产状和精细的破裂滑动分布．结果显示，地震破裂迹线与喜马拉雅主前缘逆冲断裂迹线基本一致，破裂面倾角约7°～11°．主震最大滑动量约6.0～6.6 m，余震最大滑动量约3.6～4.6 m，主震和余震滑动区和破裂空区具有互补关系．大部分破裂集中在深度8～20 km，同余震分布深度一致．深度20 km以下基本没有破裂，与地震反射INDEPTH剖面、以及根据震间期大地测量数据反演的喜马拉雅深部蠕滑剖面相对应．与各向同性弹性半空间破裂模型相比，介质不均匀模型浅部滑动较小，观测值不符值降低10%以上，特别是在北侧高海拔地区．喜马拉雅北侧高精度GPS资料，包括西藏境内GPS资料，对横向介质差异影响敏感，藏南地区GPS加密观测对喜马拉雅震源过程研究具有重要价值．横跨喜马拉雅断裂震间形变量与地震破裂滑移量直接相关．以此推算，尼泊尔中部大震原地复发周期在300 a以上．

致谢

“中国大陆构造环境监测网络”提供的GPS观测资料，UNAVCO提供了尼泊尔地区GPS观测资料，IGS数据中心提供的全球IGS跟踪站资料，加州大学圣迭戈分校提供了InSAR结果，匿名专家给了很有价值的建议，在此一并致谢！

参考文献
Ader T, Avouac J P, Zeng J L, et al. 2012. Convergence rate across the Nepal Himalaya and interseismic coupling on the Main Himalayan Thrust: Implications for seismic hazard. J. Geophys. Res. , 117(B4): B04403.
Avouac J P, Meng L S, Wei S J, et al. 2015. Lower edge of locked Main Himalayan Thrust unzipped by the 2015 Gorkha earthquake. Nature Geosci. , 8(9): 708–711.
Bilham R, Larson K, Freymueller J. 1997. GPS measurements of present-day convergence across the Nepal Himalaya. Nature , 386(6620): 61–64.
Bilham R, Gaur V K, Molnar P. 2001. Himalayan seismic hazard. Science , 293(5534): 1442–1444.
Feldl N, Bilham R. 2006. Great Himalayan earthquakes and the Tibetan plateau. Nature , 444(7116): 165–170.
Galetzka J, Melgar D, Genrich J F, et al. 2015. Slip pulse and resonance of the Kathmandu basin during the 2015 MW7.8 Gorkha earthquake, Nepal. Science , 349(6252): 1091–1095.
Gao Y, Zhou H L, Zheng S H, et al. 1997. Preliminary discussion on implication of determination on source depth of earthquake. Earthquake Research in China (in Chinese) , 13(4): 321–329.
Hauck M L, Nelson K D, Brown L D, et al. 1998. Crustal structure of the Himalayan orogen at~90° east longitude from Project INDEPTH deep reflection profiles. Tectonics , 17(4): 481–500.
Kumar S, Wesnousky S G, Jayangondaperumal R, et al. 2010. Paleoseismological evidence of surface faulting along the northeastern Himalayan front, India: Timing, size, and spatial extent of great earthquakes. J. Geophys. Res. , 115(B12): B12422.
Li Q, You X Z, Yang S M, et al. 2012. A precise velocity field of tectonic deformation in China as inferred from intensive GPS observations. Sci. China Earth Sci. , 55(5): 695–698.
Liu G, Wang Q, Qiao X J, et al. 2015. The 25 April 2015 Nepal MS8.1 earthquake slip distribution from joint of teleseismic, static and high-rate GPS data. Chinese J. Geophys. (in Chinese) , 58(11): 4287–4297.
Lindsey E O, Natsuaki R, Xu X H, et al. 2015. Line-of-sight displacement from ALOS-2 Interferometry: MW7.8 Gorkha earthquake and MW7.3 aftershock. Geophys. Res. Lett. , 42(16): 6655–6661.
Nelson K D, Zhao W J, Brown L D, et al. 1996. Partially molten middle crust beneath southern Tibet: Synthesis of project INDEPTH results. Science , 274(5293): 1684–1688.
Okada Y. 1985. Surface deformation due to shear and tensile faults in a half-space. Bull. Seismol. Soc. Amer. , 75(4): 1135–1154.
Qiao X J, Wang Q, Yang S M, et al. 2015. The 2008 Nura MW6.7 earthquake: A shallow rupture on the Main Pamir Thrust revealed by GPS and InSAR. Geodesy and Geodynamics , 6(2): 91–100.
Shan X J, Zhang G H, Wang C S, et al. 2015. Joint inversion for the spatial fault slip distribution of the 2015 Nepal MW7.9 earthquake based on InSAR and GPS observations. Chinese J. Geophys. (in Chinese) , 58(11): 4266–4276.
Sheng S Z, Wang Y G, Jiang C S, et al. 2015. Preliminary study on the static stress triggering effects on China mainland with the 2015 Nepal MS8.1 earthquake. Chinese J. Geophys. (in Chinese) , 58(5): 1834–1842.
Styron R, Taylor M H, Murphy M A. 2011. Oblique convergence, arc-parallel extension, and the role of strike-slip faulting in the High Himalaya. Geosphere , 7(2): 582–596.
Su X N, Wang Z, Meng G J, et al. 2015. Pre-seismic strain accumulation and co-seismic deformation of the 2015 Nepal MS8.1 earthquake observed by GPS. Chin. Sci. Bull. (in Chinese) , 60(22): 2115–2123.
Sun J B, Shen Z K, Xu X W, et al. 2008. Synthetic normal faulting of the 9 January 2008 Nima (Tibet) earthquake from conventional and along-track SAR interferometry. Geophys. Res. Lett. , 35(22): L22308.
Tan K, Wang Q, Ding K H, et al. 2015. Rupture models of the 2013 Lushan earthquake constrained by near field displacements and its tectonic implications. Chinese J. Geophys. (in Chinese) , 58(9): 3169–3182.
Toda S, Stein R S, Richards-Dinger K, et al. 2005. Forecasting the evolution of seismicity in southern California: Animations built on earthquake stress transfer. J. Geophys. Res. , 110(B5): B05S16.
Wang Q, Qiao X J, Lan Q G, et al. 2011. Rupture of deep faults in the 2008 Wenchuan earthquake and uplift of the Longmen Shan. Nature Geosci. , 4(9): 634–640.
Wang R J, Martín F L, Roth F. 2003. Computation of deformation induced by earthquakes in a multi-layered elastic crust-FORTRAN programs EDGRN/EDCMP. Computer & Geosciences , 29(2): 195–207.
Wang W, Wang D J, Zhao B, et al. 2014. Horizontal crustal deformation in Chinese Mainland analyzed by CMONOC GPS data from 2009-2013. Geodesy and Geodynamics , 5(3): 41–45.
Xu B, Xu C J. 2015. Numerical simulation of influences of the earth medium's lateral heterogeneity on co- and post-seismic deformation. Geodesy and Geodynamics , 6(1): 46–54.
Zhang Y, Xu L S, Chen Y T. 2015. Rupture process of the 2015 Nepal MW7.9 earthquake: Fast inversion and preliminary joint inversion. Chinese J. Geophys. (in Chinese) , 58(5): 1804–1811.
Zhao B, Du R L, Zhang R, et al. 2015a. Co-seismic displacements associated with the 2015 Nepal MW7.9 earthquake and MW7.3 aftershock constrained by Global Positioning System Measurements. Chin. Sci. Bull. (in Chinese) , 60(28-29): 2758–2764.
Zhao B, Tan K, Zhang C H, et al. 2015b. Interseismic coupling and coseismic displacements of the 2015 MW7.9 Nepal earthquake. Journal of Geodesy and Geodynamics (in Chinese) , 35(5): 734–737.
Zhao B, Huang Y, Zhang C H, et al. 2015c. Crustal deformation on the Chinese mainland during 1998-2014 based on GPS data. Geodesy and Geodynamics , 6(1): 7–15.
Zhao W J, Nelson K D, Che J, et al. 1993. Deep seismic reflectionevidence for continental underthrusting beneath southern Tibet. Nature , 366(6455): 557–559.
高原, 周蕙兰, 郑斯华, 等. 1997. 测定震源深度的意义的初步讨论. 中国地震 , 13(4): 321–329.
李强, 游新兆, 杨少敏, 等. 2012. 中国大陆构造变形高精度大密度GPS监测-现今速度场. 中国科学:地球科学 , 42(5): 629–632.
刘刚, 王琪, 乔学军, 等. 2015. 用连续GPS与远震体波联合反演2015年尼泊尔中部MS8.1地震破裂过程. 地球物理学报 , 58(11): 4287–4297.
单新建, 张国宏, 汪驰升, 等. 2015. 基于InSAR和GPS观测数据的尼泊尔地震发震断层特征参数联合反演研究. 地球物理学报 , 58(11): 4266–4276.
盛书中, 万永革, 蒋长胜, 等. 2015. 2015年尼泊尔MS8.1强震对中国大陆静态应力触发影响的初探. 地球物理学报 , 58(5): 1834–1842.
苏小宁, 王振, 孟国杰, 等. 2015. GPS观测的2015年尼泊尔MS8.1级地震震前应变积累及同震变形特征. 科学通报 , 60(22): 2115–2123.
谭凯, 王琪, 丁开华, 等. 2015. 近场位移数据约束的2013年芦山地震破裂模型及其构造意义. 地球物理学报 , 58(9): 3169–3182.
张勇, 许力生, 陈运泰. 2015. 2015年尼泊尔MW7.9地震破裂过程: 快速反演与初步联合反演. 地球物理学报 , 58(5): 1804–1811.
赵斌, 杜瑞林, 张锐, 等. 2015a. GPS测定的尼泊尔MW7.9和MW7.3级地震同震形变场. 科学通报 , 60(28-29): 2758–2764.
赵斌, 谭凯, 张彩红, 等. 2015b. 2015年尼泊尔MW7.9地震震前断层闭锁与同震位移特征分析. 大地测量与地球动力学 , 35(5): 734–737.