地球物理学报  2016, Vol. 59 Issue (4): 1359-1370   PDF    
2015年尼泊尔Mw7.8地震震源机制InSAR反演及强地面运动模拟
李永生1, 申文豪1, 温扬茂2, 张景发1, 李振洪3, 姜文亮1,4, 罗毅1    
1. 中国地震局地壳应力研究所(地壳动力学重点实验室), 北京 100085;
2. 武汉大学测绘学院, 武汉 430079;
3. 英国纽卡斯尔大学土木工程与地球科学学院, 英国纽卡斯尔, NE1 7RU;
4. 中国地震局地质研究所, 北京 100029
摘要: 2015年4月25日,在尼泊尔中部发生了Mw7.8地震.本文利用ALOS-2和SENTINEL-1A宽幅数据获取了该地震大范围的同震形变场,并反演了该地震断层破裂的几何特征及运动机制,继而以此为约束资料反演地震强地面运动.InSAR结果显示本次地震造成了巨大的地表形变,LOS向最大抬升量达到1.3 m,最大下沉量达到0.7 m.震源机制反演得到的最优的滑动分布模型表明,断层的走向为291°,倾角为7.6°,倾滑主要分布在深度为12~18 km范围,主倾滑分布范围在长度上达到了140 km,该范围内的平均倾滑角为95°.本次地震最大倾滑量达到5.3 m,位于深度15 km处.累计释放地震矩达 6.5×1020N·m,约合矩震级Mw7.8.该地震发生在印度与欧亚板块俯冲逆冲界面之间,发震构造推断为主喜马拉雅逆冲断裂,属于典型的喜马拉雅型——低角度逆断层型强震.以该滑动分布模型参数为基础利用随机振动的有限断层模型进行尼泊尔地震的强地面运动模拟,结果显示最大地震烈度为Ⅸ度,烈度分布的范围及烈度等级与USGS模型结果对比具有很高的符合度.
关键词: 尼泊尔地震     主喜马拉雅冲断裂     宽幅InSAR     强地面运动    
Source parameters for the 2015 Nepal Earthquake revealed by InSAR observations and strong ground motion simulation
LI Yong-Sheng1, SHEN Wen-Hao1, WEN Yang-Mao2, ZHANG Jing-Fa1, LI Zhen-Hong3, JIANG Wen-Liang1,4, LUO Yi1    
1. Key Laboratory of Crustal Dynamics, Institute of Crustal Dynamics, China Earthquake Administration, Beijing 100085, China;
2. School of Geodesy and Geomatics, Wuhan University, Wuhan 430079, China;
3. COMET, School of Civil Engineering and Geosciences, Newcastle University, NE1 7RU, UK;
4. Institute of Geology, China Earthquake Administration, Beijing 100029, China
Abstract: On April 25 2015, an Mw7.8 earthquake occurred in Nepal, which is the largest since the 1934 Bihar Earthquake. 2015 Nepal earthquake occurred in subduction thrust interface due to collision of Indian Plate and Eurasian Plate, which is the main seismogenic fault structure named Main Himalaya Thrust fault (MHT).The earthquake released a huge amount of energy, causing huge losses of people and infrastructure.
Using a combination of ALOS-2 ScanSAR wide mode data and SENTINEL-1A Interferometric Wide Swath data, we constructed maps of what happened on and below Earth's surface during the Mw7.8 earthquake in Nepal. Wide scan SAR data is very suitable for obtaining a wide range of deformation fields, especially for Nepal earthquake. The earthquake size is often described as slip over a larger fault area. The primarily model implies that the dimension of the fault is about 120 km×80 km in size (length×width). The conventional Strip Map mode with maximum 100 km coverage is unable to obtain the complete deformation region of Nepal earthquake. In this paper, we use the ALOS-2 and SENTINEL-1A wide scan data to obtain the wide coverage of the coseismic deformation field, then a two-step inversion strategy is employed to determine fault geometry and slip distribution, at last the seismic intensity surface motion is inverted from the distributed slip model. InSAR results show that the maximum uplift of ~1.3 m and the maximum subsidence of ~0.7 m, suggests that earthquake caused a huge ground deformation. The best-fit slip inversion model suggests that the major seismogenic fault is a thrust fault with a strike of ~291°, a dip of ~7.6° and an average rake angle of ~95°. The slip is mainly distributed of 12~18 km in the depth and of 140 km in length. The maximum amount of slip is up to 5.3 m in the depth of 15 km. The accumulative seismic moment was up to 6.5×1020 N·m, equivalent to a magnitude of Mw7.8.
On the basic of the distributed slip model, we use the finite stochastic fault model to simulate the strong ground motion of Nepal earthquake. The maximum seismic intensity is of IX degree. The range and levels of intensity distribution in this paper are compared with USGS model results, and they have a high degree of consistency.
Key words: Nepal Earthquake     Main Himalaya Thrust fault     Wide Scan InSAR     Strong ground motion    
1 引言

2015年4月25日,尼泊尔发生了Mw7.8地震,这是尼泊尔1934年8.0级地震以来最大的一次地震.在地震主震发生后的一个小时内发生一系列Mw6~7级的余震,5月12日在主震的东北又发生一个更大的余震,震级达到Mw7.3,见图 1a.美国地质调查局(USGS)给出的震源机制解显示此次地震的震中位置为84.708°E,28.147°N,节面I:走向为295°,倾角为11°,滑动角为108°;节面II:走向为96°,倾角为79°,滑动角为87°.其他的研究者也利用不同的方法对尼泊尔地震进行研究,主要包括地震波反演(Grandin et al.,2015;Yagi and Okuwaki,2015;Avouac et al.,2015;Denolle et al.,2015;张贝等,2015王卫民等,2015张勇等,2015刘静等,2015),强地面运动反演(Grandin et al.,2015),高频GPS数据(Grandin et al.,2015;Galetzka et al.,2015),静态GPS数据(Grandin et al.,2015王卫民等,2015张勇等,2015苏小宁等,2015)以及InSAR数据(Grandin et al.,2015;Avouac et al.,2015;Galetzka et al.,2015)等,上述的研究成果表明本次地震是一次典型的低倾角逆冲推覆型强震,尼泊尔地区正处于该挤压的前缘地带,这是造成地震发生的主要原因.地震破裂是从尼泊尔中部廓尔喀地区(Gorkha)开始,深度为15 km,该破裂是和低倾角向北逆冲断裂的特征一致.尼泊尔Mw7.8地震属于浅源地震,断裂面范围约150 km×90 km,破裂传播到加德满都东约150 km,不过断裂并没有出露地表,这和地表并未发现清晰的主破裂带的特征一致.

图 1 (a) 2014年4月25日尼泊尔地震的构造背景场. 黑色三角表示加德满都位置,白色五角星表示震中位置,白色虚线框表示ALOS-2和SENTINEL-1A数据空间覆盖范围,红色圆点表示主震及余震位置,红线表示断层的位置, 蓝色实线表示剖面线P1-P2. (b) 图(a)中剖面线P1-P2的板块构造图(改自Ponraj et al., 2010), 其中红 色五星表示震中在剖面线上的投影位置,蓝色粗线表示地震发生的断层位置Fig. 1 (a) Tectonic background of the April 24, 2015 Nepal Earthquake overlaid on topographic relief. Black triangle indicate the location of Kathmandu, white star means the epicenter of the mainshock of April 24 quake. White dashed line means spatial coverage of the ALOS-2 and SENTINEL-1A images. Red circle-doted points means the location and mainshock and aftershock. Red lines imply the location of the seismic fault. (b) The plate structure profile across the central Nepal, which is exhibited in (a) (modified after Ponraj et al., 2010). Red star is the epicenter projection in profile for main shock and blue bold line denotes the location of the seismic fault

尼泊尔地区的地震活动绝大多数是浅源地震,并且其现代地震活动带多出现在高喜马拉雅山山前,在尼泊尔东部(85°E—88°E)和西端(82°E以西)的地震尤其活跃,中部(82°E—85°E)不太活跃,而这次Mw7.8地震则发生在中部地区,表明这次地震的发生是长期应力聚集的结果(http://www.cags.ac.cn/YWJX/2015/0504-1.html).该地震发生在印度和欧亚两个巨型板块的边界,印度板块和欧亚板块以~40 mm·a-1的速度汇聚,过去一二十年的大地测量数据表明印度和藏南之间的相对缩短速率大致为17~21 mm·a-1(刘静等,2015).在漫长地质时间尺度上,长期的碰撞和高速率汇聚过程造成喜马拉雅山体在不断抬升.先前的地壳地震成像研究,结果表明主喜马拉雅逆冲断裂(MHT,Main Himalaya Thrust)是一个连续的断层系统,发育在中上地壳的滑脱面,由断坡—断坪组成.从Nepal地区的浅层断裂一直延伸到西藏南部的中等深度的断裂系统(Brown et al.,1996;Nábělek et al.,2009).在尼泊尔地区从南到北发育了3条主要的逆冲断裂,分别是喜马拉雅主前缘断裂(MFT,Main Frontal Thurst),喜马拉雅主边界断裂(MBT,Main boundary Thrust)和喜马拉雅主中央断裂(MCT,Main central Thrust),西藏南部发育了拆离断裂系(STD,South Tibetan Detachment Fault system),见图 1b.目前研究发现该地震发生在印度与欧亚板块俯冲逆冲界面之间,即发震构造是主喜马拉雅逆冲断裂(Avouac et al.,2015;赵文津,2015).已有研究者使用GPS测量手段显示通过尼泊尔的喜马拉雅山主逆冲断裂在深度为15~20 km处被闭锁,倾滑面宽度为100 km(Ader et al.,2012).利用InSAR方法得到加德满都以西150 km同样存在闭锁现象(Grandin et al.,2012).

目前众多合成孔径雷达在轨运行,其中2014年发射升空的两颗新的卫星系统SENTINEL-1A和ALOS-2备受关注,它们具备的宽幅成像能力在此次地震同震形变场获取过程中也得到了检验.利用差分干涉测量技术可以获取震区大范围、高精度的地表形变信息,本文利用干涉测量方法获得尼泊尔地震的InSAR同震形变场并确定断层几何特征,并对尼泊尔地震滑动机制进行反演.对于理解尼泊尔地震发生的运动学特征,认识大陆型板块碰撞俯冲的边界低角度逆冲推覆型地震动力作用、发震断层力学特性和大震孕育机理提供了条件.本文利用尼泊尔地震的滑动分布信息反演得到确定强地面运动信息,由于该地区严重缺乏足够的强地面运动记录,而模拟得到的强震数据可以补充该地区强震数据的不足,这对于未来地震动预测、工程场地地震安全性评价和地震区划工作将具有重要意义.

2 InSAR数据及同震形变场

地震发生后,日本宇航局(JAXA)新投入使用的ALOS-2卫星以及欧洲空间局去年发射的SENTINEL-1A卫星对尼泊尔地震区域进行大量成像,由于此次地震释放的能量巨大,地表的影响范围在东西向达到近200 km,传统的SAR条带模式无法一次成像能完整获取这次大范围的地表形变场.由于ALOS-2和SENTINEL-1A卫星的传感器都支持宽幅成像,非常适合本次尼泊尔地震成像.利用新开发的ALOS-2宽幅数据读取模块并嵌入到JPL/Caltech ROI_PAC软件中,经过后续的干涉处理得到尼泊尔地震的InSAR同震形变场.

ALOS-2卫星是于2014年5月24日发射升空,用于接替已经退役的ALOS-1卫星.ALOS-2卫星SAR传感器依旧采用波长较长的L波段.本文采用2015-04-05—2015-05-03干涉对,时间基线为28天,详细参数见表 1.ALOS-2具备宽幅成像的能力,在一次成像中共有5个子条带,所以比SENTINEL-1(IWS模式共3个子条带)具有更加广阔的覆盖范围,见图 1a.得益于其载荷的长波长,ALOS-2干涉图能保持非常高的相干性.从图 2a中可以看出ALOS-2宽幅模式能完整地获取得到4月25日尼泊尔地震的形变区域,充分显示出该卫星在未来地震行业应用的前景.该结果与Lindsey等发布的ALOS-2同震形变结果一致(Lindsey et al.,2015).

表 1 2015年尼泊尔地震干涉图详细信息 Table 1 Interferograms used for the 2015 Nepal Earthquake

图 2 (a) ALOS-2 SAR卫星地表形变干涉图; (b) SENTINEL-1A SAR卫星地表形变干涉图; (c) 表示图(a)、(b)中A-B剖面线在视线向上的位移比较,红线表示ALOS-2,蓝线表示SENTINEL-1AFig. 2 (a) Original ALOS-2 coseismic interferogram. (b) Original SENTINEL-1A coseismic interferogram. (c) Comparison of modelled LOS surface deformation along the profile A-B. Red line denotes ALOS-2 and blue line means SENTINEL-1A

2014年8月24日,美国加利福尼亚州的纳帕峡谷(Napa Valley)发生Mw6.1地震,SENTINEL-1A第一次获取的该地震的同震形变场,至此揭开该卫星在地震等地质灾害领域应用的序幕(Wei et al.,2015;李永生等,2015).在尼泊尔地震中,SENTINEL-1A卫星采用全新的干涉宽模式,特别适合于大空间范围(>100 km)地表形变监测.本文中为监测尼泊尔地震,采用干涉对为2015-04-17—2015-04-29,时间基线只有12天,干涉对参数见表 1.虽然SENTINEL-1A的载荷为C波段,而且地震震中区域植被茂密、地壳起伏较大,但是从干涉结果中可以看出在震中位置表现出较高的相干性.本文采用Gamma软件处理SENTINEL-1A IWS模式数据,并采用由两个相邻的Frame拼接得到最终的干涉图.虽然该干涉对由于SENTINEL-1A成像视角的原因没有完整获取到震中形变场,但也可以看出大部分形变场特征,见图 2b.

从两种卫星平台独立获取的同震形变场(见图 2a2b),发现地震InSAR形变场在空间上分布基本一致,在整个极震区呈现两个形变中心.在南形变中心表现出LOS向抬升现象,最大的抬升量达到了1.3 m.该形变中心位于加德满都北侧,从地震发生的范围来看,加德满都的建筑物因此受到极大破坏.北形变中心表现出地表沉降现象,在LOS向下沉量达到0.7 m.从尼泊尔地区的地质构造背景以及本次地震释放出的能量推断,南北两个形变中心位于断层的上盘,这也符合低倾角逆断层强震的变形特征.为了对形变干涉场精度进行评估,在图 2a2b中分别选择一条剖面线,两条剖面线横跨两个形变中心,其代表的空间位置是一致的.从图 2c中看出不同平台独立获取的形变量在空间分布上是一致的,从侧面也验证了本文的干涉形变场的精确性.

3 InSAR震源机制反演 3.1 震源机制反演方法

本文以ALOS-2以及SENTINEL-1A卫星数据获取得到的尼泊尔地震的同震形变资料作为反演约束,对本次地震的断层参数以及滑动机制进行详细研究和分析.在获取InSAR形变干涉图之后进行降采样处理,降采样采用数据分辨率约束的四叉树方法(Lohman and Simons,2005),降采样后ALOS-2和SENTINEL-1A数据点的个数分别为291和379个,这极大地降低反演的计算量.本文采用二步反演法对断层的破裂几何参数以及倾滑分布进行估计(Feng et al.,2013),首先假设断层是均一断层模型(Uniform fault model),计算出断层的空间几何参数.再采用分布式断层模型(Distributed fault model),计算出断层面上的分布式滑动量.本文利用PSOKINV软件进行反演(冯万鹏和李振洪,2010Li et al.,2011;Feng et al.,2013),其采用粒子群优化算法(Particle Swarm optimization,PSO)进行震源参数反演.PSO是一种基于迭代的优化算法,初始化为一群随机解并通过迭代方法寻找最优解.冯万鹏等分析了InSAR形变场资料在震源参数反演过程中的参数特点,采用PSO震源参数反演策略在较少参数约束下即可确保获得较高成功率(冯万鹏和李振洪,2010).

3.2 二步法最优滑动分布反演

使用均一断层模型(Uniform fault model)进行反演主要目的是确定断层的位置(经纬度)、顶部埋深、走向、倾角等参数.本文中所使用的适配函数定义为

其中系数矩阵G表示均匀断层上倾滑量达到1 m时引起的地表运动响应,即为格林函数.S表示倾滑矢量,W表示每个数据集的相对权重,D表示地表形变观测值,N表示形变观测值的个数(Feng et al.,2013).以ALOS-2和SENTINEL-1A同震形变场为约束资料,采用PSO方法在全参数域内寻找到使得适配函数达到最小的解.利用均一断层模型后得到断层最优拟合断层参数位于(85.43°E,27.93°N),深度为15 km,断层的走向为291°,倾向为8.9°,倾滑角为94°,震级达到Mw7.75,见表 2.

表 2 基于InSAR同震干涉图形变资料反演获取纳帕地震最优断层震源参数 Table 2 Optimal fault geometric parameters determined with InSAR coseismic interferogram

在均一断层模型的基础上进一步确定地震滑动的空间分布,采用分布式断层模型(Distributed fault model).断层面在走向和倾向上分别扩展到200 km×200 km.在反演过程中离散子断层的大小选择为6 km×6 km,并采用二次差分拉普拉斯算子对滑动粗糙度进行约束(Harris and Segall,1987).由于基于均一断层模型的假设对于分布式滑动模型并不是最优,所以在分布式滑动反演研究中,倾角与平滑系数(α2)需要进一步的优化.具体反演算法见公式(2):

其中S表示每个子断层的倾滑矢量, D表示地表形变观测值,L表示二阶微分算子,用于滑动粗糙度评价,α2为平滑系数,用于平滑解析解(Feng et al.,2013).

本研究中通过固定不同倾角值,调整不同的平滑系数对模型粗糙度和残差的变化趋势进行分析.以均方根误差(Root Mean Square Error)代表残差ξ,定义同公式(1).模型的粗糙度ψ的定义如下:

其中p=LSn表示子断层个数.构建log函数模型f(δ,α)=log(ψ+ξ)用于获得反演的最优粗糙度ψ和残差ξ.在归一化后,粗糙度曲线表现为单调递增函数,而残差曲线则表现为单调递减函数,见图 3a.从图中可以看出,α2≈2.1时可以获得最优平滑系数.在图 3b中可以看出对于给定的倾角范围〖5.5° 9.5°〗和平滑系数α2范围〖1.0 5.0〗,获得其全局最小点即为最优倾角以及平滑系数,见图 3b中的五角星,得到最优倾角为7.6°.

图 3 (a) 模型残差曲线和模型粗糙度折中曲线图,其中的虚线表示归一化后的模型粗糙度,点画线表示归一化后的模型残差曲线,实线表示归一化后的log(ξ+ψ); (b)表示log函数模型log(ξ+ψ)的轮廓线图,该轮廓线是由倾角Dip 和平滑系数α2所决定的,白色五角星表示全局最小解Fig. 3 (a) A trade-off curve line for residuals roughness. The dashed and dash-dotted black lines show the trends of model roughness and the residuals of modeled simulations after normalizing ([ξ,ψ]), respectively, while the solid grey line represents normalizing log(ψ+ξ). (b) Contour map of log(ψ+ξ) with variations of dips and hyperparameters (α2). White star indicates the point of global minimum

3.3 地震震源反演结果图 4a为最优的滑动分布模型,从反演结果中可以看出,断层的走向为291°,倾角为7.6°,倾滑主要分布在深度为12~18 km范围,断层倾滑以逆冲为主.主倾滑分布范围在长度上达到140 km,该范围内的平均倾滑角为95°,详细参数见表 2.此次地震,累计释放地震矩达6.5×1020N·m,约合矩震级Mw7.8,见图 4b.基于蒙特卡洛算法分析100组加入扰动的采样点得到断层面倾滑量的不确定性,见图 4c,不确定性最大值约为0.5 m,约为同震最大滑动量的10%.本文反演得到本次地震最大倾滑量达到5.3 m,位于深度15 km处.美国JPL(Jet Propulsion Laboratory)联合远震波形数据以及GPS和InSAR等近场形变数据反演得到的最大倾滑量超过了6.0 m(http://aria.jpl.nasa.gov/node/43),USGS利用远场P波、远场SH波和长周期面波反演得到的最大滑移量为3.0~4.0 m.本文中倾滑分布范围和方法与上述机构公布的结果基本一致,但本文的倾滑量要大于USGS的结果,主要原因是本文使用的反演约束数据主要是静态的近场形变场,而且InSAR同震形变场可能还包含了部分震后形变信号以及强余震影响(李永生等,2015苏小宁等,2015).

图 4 InSAR反演断层空间滑动模型
(a) 基于InSAR反演的滑动分布模型; (b) 表示对应模型的地震矩累计分布; (c) 滑动分布模型的不确定性.
Fig. 4 The slip models determined with InSAR datasets
(a) InSAR based slip model; (b) The accumulative seismic moment distribution along the depth for slip model; (c) Uncertainties of the slip model.

图 5给出由最优拟合滑动模型给出的拟合结果.图 5a5d表示ALOS-2和SENTINEL-1A InSAR观测值,图 5b5e表示分别基于ALOS-2和SENTINEL-1A形变场反演模拟得到的干涉形变场,图 5c5f分别表示图 5b5e两种情况下的残差.从图 5c5f的残差图中可以看出,在主形变场的边缘存在一些残余形变信号,推测原因可能是同震引起的山体滑坡造成的.总体而言,InSAR模拟干涉图能很好地对干涉形变场进行解释,证明本文确定的滑动区间稳定可靠.反演过程中我们通过对InSAR观测数据和同震模拟形变资料进行相关系数对比,当选取最优倾滑模型时,相关系数达到0.973.另外从观测值和模拟值的残差值可以看出(图 5c5f),震中区域的最大反演残差约为8 cm.从而证明了本文结果的可靠性.

图 5 InSAR资料约束下同震形变场(a,d),最优拟合模型预测位移场(b,e)以及残差分布(c,f), (a)、(b)及(c)表示ALOS-2卫星,(d)、(e)及(f)表示SENTINEL-1A卫星Fig. 5 Original (a,d), modelled (b,e) and residual interferograms (c,f) for InSAR models. (a),(b),(c) mean ALOS-2, (d),(e),(f) mean SENTINEL-1A
4 强地面运动模拟 4.1 模型参数确定

基于InSAR资料反演得到的震源机制解为约束进行强地面运动模拟研究,本文采用随机振动有限断层模型进行强地面运动模拟,在该模型中断层面被分成N个大小相等的矩形子断层,每一个子断层即为一个点源,通过设定震源模式和破裂传播速度,就可以得到子源破裂的时间顺序,然后根据子源与场地的几何关系,计算每个子源对场地的影响,即可叠加成场地的地震动时程.针对地震动随机合成中的不足之处,Motazedian等提出应用动力学拐角频率保证地震矩和断层辐射的总能量均是守恒的(Motazedian and Atkinson,2015),在高频部分引入保证远场高频辐射能守恒的标度因子,并引入脉冲子断层(Pulsing subfault)百分比的概念,断层上子断层的划分不受大小控制,从而避免子震的多次触发这一非物理现象,至此该方法理论框架已基本完善,并且已经广泛用于世界各地的地表运动研究中(Atkinson and Boore,2006;Ameri et al.,2011).

发震断层震源模型首先需要确定断层的全局震源参数,该类参数主要用来描述发震断层的宏观几何特征和平均破裂过程,包括断层破裂面积、长度、宽度、埋藏深度、走向、倾角、地震矩、应力降、Q模型、几何衰减模型等.应力降是随机有限断层模型最关键的参数之一,尤其对近场强地面运动有着至关重要的影响.Allmann等计算全球范围内1990—2007年间约2000个mb≥5.5中强地震的应力降值,结果显示应力降值在0.3~50 MPa之间,且平均应力降的值并未随地震矩变化(Allmann and Shearer,2009),根据其对喜马拉雅地区的研究结果,本文中我们取应力降值为2.0 MPa.Sharma和Wason(1994)通过对一系列喜马拉雅地区的中小震研究发现该地区地震应力降水平相对较低,这说明该地区地壳承受累积应变能的能力相对较低,这与喜马拉雅处在俯冲带区域有关系.断层的长度、走向、倾角、埋深、位置等几何参数分别来自USGS和InSAR数据反演,而Q模型参数的确定来自Mukhopadhyay等(2006)的研究,取Q=113f1.01.几何衰减模型我们采用R-1模型,R是震源距.就平均破裂速度而言,大量的反演结果表明,一般在剪切波速的0.6~0.9倍之间,对本研究的所有发震断层我们取平均破裂速度为剪切波速的0.8倍.表 3总结了模型计算参数.

表 3 随机有限断层模型模拟尼泊尔Mw7.8地震强地面运动参数表 Table 3 Simulating strong motion parameters for Nepal Mw7.8 earthquake by Finite Stochastic Model
4.2 模拟结果

Somerville等(1999)的研究指出:反演过程中的断层尺度通常被设定为至少大到足以容纳整个断层破裂,因此反演过程中通常高估了破裂区域的实际尺寸,我们使用Somerville等(1999)的标准判据减小矩形USGS和InSAR反演断层滑动模型的尺寸,修整了滑动模型的边缘(Somerville et al.,1999).地震发生后USGS基于远场地震波进行了震源破裂过程的快速反演,给出了地震破裂的滑动分布模型.我们分别基于USGS模型和InSAR反演模型作为滑动分布的输入值,利用随机振动有限断层模型进行尼泊尔地震的强地面运动模拟,见图 6.

图 6 基于两种滑动模型得到的尼泊尔Mw7.8地震模拟强地面运动
(a)、(b)、(c)分别为USGS滑动模型得到的PGA,PGV和烈度分布; (d)、(e)、(f)分别为 InSAR反演滑动模型得到的PGA,PGV和烈度分布.蓝色五角星为震中位置.
Fig. 6 Strong motion simulation for Nepal Mw7.8 earthquake based on two slip models
(a),(b),(c) indicate PGA, PGV and MMI distribution estimated from USGS slip model respectively; (d),(e),(f) indicate PGA, PGV and MMI distribution estimated from InSAR inversion slip model respectively.

本文针对尼泊尔地震影响范围广、受灾地区多的特点将研究区域划定在(25.5°N—30°N,83°E—87.5°E),区域网格大小为4.5°×4.5°,网格精度为0.05°,计算每个网格点的峰值加速度(PGA)、峰值速度(PGV)以及网格点到断层的最短距离.我们利用Boore和Atkinson建立衰减关系时对浅层平均横波速度Vs30的处理方法加入了Vs30对模型模拟结果的影响,使最终结果更加合理(Boore and Atkinson,2007).从图 6中看到,两种模型模拟结果在运动水平及分布范围方面有诸多相似点:两种模型模拟得到的PGA、PGV均沿断层走向方向分布;由于受到西南部平原区域低值浅层横波速度影响,模拟结果均体现出不对称性:在断层西南侧运动水平要高于断层东北侧水平.而相同PGA或者PGV水平情况下,USGS模型模拟结果对应的范围要大于InSAR反演模型结果的,例如PGA≥300 cm·s-2,PGV≥40 cm·s-1对应的范围在图 6a6b中要略大于图 6d6e中,两种模型模拟结果的差异性可能与滑动分布模型的不均匀性有关.

地震发生后USGS基于其“Shakemap”系统给出PGA,PGV和烈度分布示意图,5月1日中国地震局根据灾情、实地调查数据、观测数据、卫星遥感震害解译和全球有关机构数据,参照《中国地震烈度表》(GB/T 17742-2008),对此次地震烈度分布进行评估(中国地震局,2015).通过与这两个机构发布的结果对比,我们发现本文模拟结果在烈度分布的范围以及烈度等级方面与公布结果符合度很高:高烈度区(≥VIII)主要集中在震中东南侧区域及加德满都北部区域,见图(6c,6f).震中西部区域烈度明显高于其东部区域,说明此次地震破裂方式以单向破裂为主,破裂自震源开始沿东南向传播,这一点与汶川地震破裂方式非常类似(王卫民等,2008).

5 讨论 5.1 地震背景及能量释放

尼泊尔是世界上最易发生大型地震的地区之一,原因在于尼泊尔位于两个大陆板块之间,并且两板块在剧烈地碰撞挤压.该挤压汇聚在印度北部以及西藏南部之间存在一个显著的调节过程,在尼泊尔地区的形变速率为18~20 mm·a-1(Ader et al.,2012),常常认为在一个单一的断层结构(主喜马拉雅逆冲断层,Main Himalyan Thrust)由浅层断裂向中层断裂调整,见图 1b.该挤压蓄积的能量其中一部分加载到地壳断层中,并以弹性应力方式存储.在地震发生时,地壳存储的能量会瞬间释放,造成地表巨大的破坏.但根据GPS监测和地质学研究的结果,部分地震专家认为以这次地震的震源深度(15 km)判断此次地震并没有完全释放加德满都地区地壳运动积累的能量(http://bonerjea.org/2015/04/30/bigger-earthquake-might-be-coming-to-nepal/).先前研究结果表明至少需要10~15 m的地表形变量才能完全释放地壳积累的能量,但是本次地震地表形变只有约3 m(Oskin,2015),因此该地区未来依然需要重点关注.本次尼泊尔地震造成非常明显的地表形变,足以说明该次地震释放的能量巨大,但是发震断层破裂带并未出露地表,其原因可能是断层倾角较小(约7.6°),沿断层面向地表破裂过程无法持续近100 km,具体原因需要进一步研究.中国的天山两侧、祁连山山前都具备发生低角度逆冲推覆型地震的条件,对这类地震深入研究对于中国的防震减灾和地震预测非常重要(刘静等,2015).

5.2 强地面运动模拟

利用随机有限断层模型基于USGS和InSAR反演滑动模型对尼泊尔Mw7.8地震进行强地面运动评估,将本文的模拟结果与USGS公布的Shackmap结果和中国地震局调查烈度图(中国地震局,2015)进行对比,结果显示在极震区分布特征和衰减特征方面两者符合较好.由于该地区缺乏足够的强地面运动记录,无法利用实际记录到的地震动资料直接得到地震动衰减关系,而模拟得到的强震数据可以补充该地区强震数据的不足,利用这些强震数据可以完善该地区衰减关系曲线(Atkinson and Boore,1995),这对于未来地震动预测、工程场地地震安全性评价和地震区划工作将具有重要意义.

在喜马拉雅山脚前平原地区出现明显的烈度增强,PGA、PGV在该区域存在明显放大现象,这主要是因为在该平原区及山前盆地区域沉积层较厚,浅层平均横波速度Vs30较低,直接导致地震波的放大效应.尼泊尔当地房屋抗震能力较低,这种放大效应更是加剧山前盆地的加德满都的震害.模拟结果与调查结果的烈度衰减趋势上的一致说明本文使用的几何衰减模型以及Q值衰减模型是合适的.在模拟过程中我们忽略地形因素的影响,事实上在山麓地区地形对地震波的影响也很重要,此次地震我们也注意到有雪崩和泥石流现象的发生,这些现象是否跟地形坡度有关、如何科学有效增加地形因素的影响也是我们未来需要关注的.

目前利用大地测量数据进行反演的最大缺陷在于其时效性不足,往往需要几天甚至十几天才能获得形变数据,这使得反演结果的实时效率大打折扣.而利用远场地震波反演虽无时效性缺陷却存在因其覆盖点少导致的反演结果空间分辨率不足的问题,地震断层破裂面上的滑移分布直接控制着强地震波各种频率分量的激发,若反演结果空间分辨率不足则会影响强地面运动模拟结果的频率带宽准确性.随着InSAR技术发展,数据产出效率不断提高,结合形变数据与震波数据进行联合反演将是未来趋势,联合反演在数据采样率方面互为补充,确保信号的更全面覆盖.另外,在本文中我们只采用了单一的Q衰减模型,事实上,许多研究已经表明Q值是依赖于震中距变化的(Tsang et al.,2010),因此今后工作中我们需要建立随距离变化的Q值模型来改善模拟结果.另外,本研究在将地震动峰值参数转化为烈度时使用的公式是Atkinson和Sonley(2000)基于美国加州地震数据给出的经验公式,虽然喜马拉雅地区与美国加州地区都属于板块活跃的俯冲带区域,但由于各个地方地质构造和场地衰减特征不同,转换公式的系数也会有所差异,在实际应用中结果与实际情况也会存在一定差异.

6 结论

2015年4月25日发生的尼泊尔地震是印度板块与欧亚板块持续的南北向汇聚,印度板块向北俯冲在欧亚板块之下而产生的逆冲断裂作用的结果,地震造成地表向南剧烈运动,位于加德满都的前锋地表向南运动量达到3 m.断层破裂起始于西部并沿着断层向东延续150 km.基于ALOS-2和SENTINEL-1A两种卫星传感器获取的地震同震形变场进行地震破裂机制反演,采用二步反演策略,结果表明Nepal地震的发震断层的走向为291°,倾角为7.6°,断层倾滑以逆冲为主,最大倾滑量达到5.3 m,位于深度15 km处.主倾滑分布范围在长度上达到140 km,该范围内的平均倾滑角为95°,累计释放地震矩达6.5×1020 N·m,约合矩震级Mw7.8.在InSAR震源机制反演的基础上进行强地面运动模拟,InSAR模拟结果与USGS结果对比具有较高的吻合度.通过对这次地震的研究可以对发震断层的结构进行详细描述,并为进一步研究主逆冲断裂几何特征以及其主震滑动和震间断层的活动之间的空间关系奠定基础.

致谢 感谢两位匿名评审专家给出的宝贵修改意见,对本文的结构完整和逻辑严谨起到极大的推动作用.感谢JAXA提供的ALOS-2数据(1368),牛津大学为本文提供了SENTINEL-1A形变资料,本文部分图件使用GMT软件绘制,在此一并表示感谢!

参考文献
[1] Ader T, Avouac J P, Jing L Z, et al. 2012. Convergence rate across the Nepal Himalaya and interseismic coupling on the Main Himalayan Thrust: Implications for seismic hazard. Journal of Geophysical Research, 117: B04403, doi: 10.1029/2011JB009071.
[2] Allmann B P, Shearer P M. 2009. Global variations of stress drop for moderate to large earthquakes. Journal of Geophysical Research, 114: B01310, doi: 10.1029/2008JB005821.
[3] Ameri G, Emolo A, Pacor F, et al. 2011. Ground-motion simulations for the 1980 M6.9 Irpinia earthquake (Southern Italy) and scenario events. Bulletin of the Seismological Society of America, 101(3): 1136-1151.
[4] Atkinson G M, Boore D M. 1995. Ground-motion relations for eastern North America. Bull. Seismol. Soc. Am., 85(1): 17-30.
[5] Atkinson G M, Boore D M. 2006. Earthquake ground-motion prediction equations for eastern North America. Bull. Seismol. Soc. Am., 96(6): 2181-2205.
[6] Atkinson G M, Silva W. 2000. Stochastic modeling of California ground motions. Bulletin of the Seismological Society of America, 90(2): 255-274.
[7] Atkinson G M, Sonley E. 2000. Empirical relationships between modified Mercalli intensity and response spectra. Bulletin of the Seismological Society of America, 90(2): 537-544.
[8] Atkinson G M, Boore D M. 2006. Earthquake ground-motion prediction equations for eastern North America. Bulletin of the Seismological Society of America, 96(6): 2181-2205.
[9] 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 Geoscience, 8(9): 708-711, doi: 10.1038/ngeo2518.
[10] Boore D M, Atkinson G M. 2007. Boore-Atkinson NGA ground motion relations for the geometric mean horizontal component of peak and spectral ground motion parameters. PEER Report 2007/01. Berkeley, California: Pacific Earthquake Engineering Research Center.
[11] Brown L D, Zhao W J, Nelson K D, et al. 1996. Bright spots, structure, and magmatism in Southern Tibet from INDEPTH seismic reflection profiling. Science, 274(5293): 1688-1690, doi: 10.1126/science.274.5293.1688.
[12] China Earthquake Administration. 2015. "China Earthquake Administration released 8.1 of Nepal earthquake intensity map", http://www.cea.gov.cn/publish/dizhenj/464/515/201505012211 23823783663/index.html[2015-05-01].
[13] Denolle M A, Fan W Y, Shearer P M. 2015. Dynamics of the 2015 M7.8 Nepal earthquake. Geophys. Res. Lett., 42(18): 7467-7475, doi: 10.1002/2015GL065336.
[14] Feng W P, Li Z H. 2010. A novel hybrid PSO/simplex algorithm for determining earthquake source parameters using InSAR data. Progress in Geophysics (in Chinese), 25(4): 1189-1196, doi: 10.3969/j.issn.1004-2903.2010.04.007.
[15] Feng W P, Li Z H, Elliott J R, et al. 2013. The 2011 Mw6.8 Burma earthquake: fault constraints provided by multiple SAR techniques. Geophysical Journal International, 195(1): 650-660, doi: 10.1093/gji/ggt254.
[16] Galetzka J, Melgar D, Genrich J F, et al. 2015. Slip pulse and resonance of the Kathmandu basin during the 2015 Gorkha earthquake, Nepal. Science, 349(6252): 1091-1095.
[17] Grandin R, Doin M P, Bollinger L, et al. 2012. Long-term growth of the Himalaya inferred from interseismic InSAR measurement. Geology, 40(12): 1059-1062, doi: 10.1130/G33154.1.
[18] Grandin R, Vallée M, Satriano C, et al. 2015. Rupture process of the Mw=7.9 2015 Gorkha earthquake (Nepal): insights into Himalayan megathrust segmentation. Geophys. Res. Lett., 42(20): 8373-8382, doi: 10.1002/2015GL066044.
[19] Harris R A, Segall P. 1987. Detection of a locked zone at depth on the Parkfield, California, segment of the San Andreas fault. Journal of Geophysical Research: Solid Earth (1978—2012), 92(B8): 7945-7962.
[20] Li Y S, Feng W P, Zhang J F, et al. 2015. Coseismic slip of the 2014 MW6.1 Napa, California Earthquake revealed by Sentinel-1A InSAR. Chinese Journal of Geophysics (in Chinese), 58(7): 2339-2349, doi: 10.6038/cjg20150712.
[21] Li Z H, Elliott J R, Feng W P, et al. 2011. The 2010 MW6.8 Yushu (Qinghai, China) earthquake: constraints provided by InSAR and body wave seismology. Journal of Geophysical Research: Solid Earth, 116(B10): B10302.
[22] 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, doi: 10.1002/2015GL065385.
[23] Liu J, Ji C, Zhang J Y, et al. 2015. Tectonic setting and general features of coseismic rupture of the 25 April, 2015 Mw7.8 Gorkha, Nepal earthquake. Chin. Sci. Bull. (in Chinese), 60(27): 2640-2655.
[24] Lohman R B, Simons M. 2005. Some thoughts on the use of InSAR data to constrain models of surface deformation: Noise structure and data downsampling. Geochemistry, Geophysics, Geosystems, 6: Q01007.
[25] Motazedian D, Atkinson G M. 2005. Stochastic finite-fault modeling based on a dynamic corner frequency. Bulletin of the Seismological Society of America, 95(3): 995-1010.
[26] Mukhopadhyay S, Tyagi C, Rai S S. 2006. The attenuation mechanism of seismic waves in northwestern Himalayas. Geophysical Journal International, 167(1): 354-360.
[27] Nábělek J, Hetényi G, Vergne J, et al. 2009. Underplating in the Himalaya-Tibet collision zone revealed by the Hi-CLIMB experiment. Science, 325(5946): 1371-1374.
[28] Oskin B. 2015. Bigger Earthquake Coming on Nepal's Terrifying Faults. https://in.news.yahoo.com/bigger-earthquake-coming-nepals-terrifying-faults-104617764.html.
[29] Ponraj M, Miura S, Reddy C D, et al. 2010. Estimation of strain distribution using GPS measurements in the Kumaun region of Lesser Himalaya. Journal of Asian Earth Sciences, 39(6): 658-667, doi: 10.1016/j.jseaes.2010.04.037.
[30] Somerville P, Irikura K, Graves R, et al. 1999. Characterizing crustal earthquake slip models for the prediction of strong ground motion. Seismological Research Letters, 70(1): 59-80.
[31] 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, doi: 10.1360/N972015-00534.
[32] Sharma M L, Wason H R. 1994. Occurrence of low stress drop earthquakes in the Garhwal Himalaya region. Physics of the Earth and Planetary Interiors, 85(3): 265-272.
[33] Tsang H H, Sheikh M N, Lam N T K, et al. 2010. Regional differences in attenuation modeling for Eastern China. Journal of Asian Earth Sciences, 39(5): 441-459.
[34] Wang W M, Zhao L F, Li J, et al. 2008. Rupture process of the Ms8.0 Wenchuan earthquake of Sichuan, China. Chinese Journal of Geophysics (in Chinese), 51(5): 1403-1410, doi: 10.3321/j.issn:0001-5733.2008.05.013.
[35] Wang W M, Hao J L, He J K, et al. 2015. Rupture process of the Mw7.9 Nepal earthquake April 25, 2015. Science China: Earth Sciences, 58(10): 1895-1900, doi: 10.1007/s11430-015-5170-y.
[36] Wei S J, Barbot S, Graves R, et al. 2015. The 2014 Mw6.1 South Napa Earthquake: A unilateral rupture with shallow asperity and rapid Afterslip. Seismological Research Letters, 86(2A): 344-354.
[37] Yagi Y, Okuwaki R. 2015. Integrated seismic source model of the 2015 Gorkha, Nepal, earthquake. Geophys. Res. Lett., 42(15): 6229-6235, doi: 10.1002/2015GL064995.
[38] Zhang B, Cheng H H, Shi Y L. 2015. Calculation of the co-seismic effect of MS8.1 earthquake, April 25, 2015, Nepal. Chinese Journal of Geophysics (in Chinese), 58(5): 1794-1803, doi: 10.6038/cjg20150529.
[39] 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 Journal of Geophysics (in Chinese), 58(5): 1804-1811, doi: 10.6038/cjg20150530.
[40] Zhao W J. 2015. Geological background of the Nepal's Ms8.1 earthquake and its trend in the future. Chin. Sci. Bull. (in Chinese), 60(21): 1953-1957.
[41] 李永生, 冯万鹏, 张景发等. 2015. 2014年美国加州纳帕MW6.1地震断层参数的Sentinel-1A InSAR反演. 地球物理学报, 58(7): 2339-2349, doi: 10.6038/cjg20150712.
[42] 冯万鹏, 李振洪. 2010. InSAR资料约束下震源参数的PSO混合算法反演策略. 地球物理学进展, 25(4): 1189-1196, doi: 10.3969/j.issn.1004-2903.2010.04.007.
[43] 刘静, 纪晨, 张金玉等. 2015. 2015年4月25日尼泊尔Mw7.8级地震的孕震构造背景和特征. 科学通报, 60(27): 2640-2655.
[44] 苏小宁, 王振, 孟国杰等. 2015. GPS观测的2015年尼泊尔MS8.1级地震震前应变积累及同震变形特征. 科学通报, 60(22): 2115-2123.
[45] 王卫民, 赵连锋, 李娟等. 2008. 四川汶川8.0级地震震源过程. 地球物理学报, 51(5): 1403-1410, doi: 10.3321/j.issn:0001-5733.2008.05.013.
[46] 王卫民, 郝金来, 何建坤等. 2015. 2015年4月25日尼泊尔Mw7.9级地震震源过程. 中国科学: 地球科学, 45(9): 1421-1426.
[47] 张贝, 程惠红, 石耀霖. 2015. 2015年4月25日尼泊尔MS8.1大地震的同震效应. 地球物理学报, 58(5): 1794-1803, doi: 10.6038/cjg20150529.
[48] 张勇, 许力生, 陈运泰. 2015. 2015年尼泊尔Mw7.9地震破裂过程: 快速反演与初步联合反演. 地球物理学报, 58(5): 1804-1811, doi: 10.6038/cjg20150530.
[49] 赵文津. 2015. 尼泊尔大地震发生的构造背景及发展趋势. 科学通报, 60(21): 1953-1957.
[50] 中国地震局. 2015. 《中国地震局发布尼泊尔8.1级地震烈度图》http://www.cea.gov.cn/publish/dizhenj/464/515/201505012 21123823783663/index.html.[2015-05-01].