2019年巴基斯坦新米尔普尔MW6.0地震的同震形变场与断层滑动分布反演

贾 蕊 张国宏 解朝娣 单新建 张迎峰 李成龙 黄自成

1)中国地震局地质研究所，地震动力学重点实验室，北京 100029 2)云南大学地球科学学院，地球物理系，昆明 650500

0 引言

2019年9月24日11时01分(Universal Time Coordinated，UTC)，巴基斯坦北部距新米尔普尔(New Mirpur)7km处发生MW6.0地震。美国地质调查局(United States Geological Survey，USGS)利用远场地震波数据进行测定，给出的地震震中位于(33.078°N，73.794°E)，震源深度为11.5km(表1)。

表1 不同机构给出的2019年12月24日巴基斯坦新米尔普尔MW6.0地震震源机制解Table1 Source parameters of the MW6.0 earthquake in New Mirpur，Pakistan，24 December，2019 from different agencies

本次巴基斯坦MW6.0地震很可能是由于印度板块和欧亚板块碰撞系统内的逆冲断裂造成的，巴基斯坦北部、印度和阿富汗邻近地区的地震活动性和活断层都是印度板块与欧亚板块碰撞系统的直接结果。印度板块持续以45～52mm/a的速率向N移动，从而在该碰撞区引起地壳缩短和变形，这使得包括巴基斯坦西部和北部在内的喜马拉雅地区成为世界上地震活动性最活跃的地区之一(邓起东等，2014)。喜马拉雅造山带作为印度板块与欧亚板块碰撞形成的产物，主要由4条断裂带组成，从北向南依次为藏南拆离系(Southern Tibetan detachment system，STDS)、主中央逆冲带(main central thrust，MCT)、主边界逆冲带(main boundary thrust，MBT)以及主前缘逆冲带(main frontal thrust，MFT)(杨晓平等，2016)。本次地震的震中靠近喜马拉雅主前缘逆冲次级断裂(图1)。MFT是地球上最大的活动碰撞造山带中最长的活动收缩构造带，大体呈EW走向。自全新世以来，MFT作为板块边界断裂带，吸收了两大板块间至少 1/3 的碰撞缩短量。如此大规模的高速率会聚活动使得MFT断裂带上的地震活动频度高、强度大，历史上曾发生过多次7级以上大地震，如1897年印度西隆8.7级地震、2005年巴基斯坦克什米尔7.8级地震以及2015年尼泊尔8.1级地震等，这些地震分别发生在MFT的不同段落上(吴中海等，2016)。陆地上的活动断层造成了多次大型破坏性地震，产生了地表断裂，促进了构造地貌的发育。同震地表破裂的详细观测不仅为了解地震本身和估计地震重现周期提供了基本信息，而且为解释其他地区构造地貌特征提供了帮助。然而，逆断层地表断裂比走滑断层或正断层更少见。2005年克什米尔地震是喜马拉雅地区历史上第1次地表断裂事件，也是该地区已知的最大地震。在2019年巴基斯坦6.0地震后，至今还没有学者使用InSAR技术观测研究此次地震。研究此次地震为评估活动逆断层的地震危险性以及了解巴基斯坦北部的地震活动性提供了难得的机会。

图1 2019年12月24日巴基斯坦新米尔普尔MW6.0地震周边构造背景图Fig. 1 Tectonic background of the MW6.0 earthquake in New Mirpur，Pakistan，24 December，2019.a 区域构造背景图； 红色实线代表活动断层，红色沙滩球代表本研究所得的震源机制解，黑色沙滩球代表GCMT确定的本次地震的震源机制解，黄色沙滩球代表USGS确定的本次地震的震源机制解，黄色圆圈代表1919年至今4级及以上地震震中，黑色箭头代表板块的运动方向，黑色小方框表示SAR数据的覆盖范围。b 红色五角星指示本次地震震中(USGS)，白色 方框内为轨道T100和T107的Sentinel-1A影像范围

如图1b所示，该地区在历史上发生过多次中强地震。在全球范围内，平均每年将发生约10次7级以上破坏性地震及超过100万次破坏性有限甚至没有地表破裂但仍能被地震仪器记录到的小地震(震级范围为2.5～4.5级)，而在这两者之间，更有成百上千次中强地震(震级范围为5.5～6.5级)发生，破坏性显著。中强地震往往容易被忽视，只有极少数由于震源深度很浅而造成人员伤亡及结构破坏，得以被人们关注。构造应力如何沿着断层累积并转化为大地震引发能量的突然释放，这仍然是一个悬而未决的问题。对大量中强地震的震源破裂以及破裂过程开展研究，可以更详细地了解构造载荷如何从静止期发展至活跃期，对分析地震活动性以及大地震的突然爆发之间的关系也有一定参考价值。

利用区域或远场地震数据可以反演出中强地震的震源机制，用于分析目标地震的断层特征。震源机制解提供2个节面，然而由于中强地震往往不会产生地表破裂，故很难从2个节面中区分出发震断层。尽管地震台站的仪器时间采样精度非常高，但其空间覆盖率往往不足以区分出真实的震源机制解。综上所述，针对中强地震震源断层的进一步研究具有实际意义，需要利用卫星及大地测量等传统地震学以外的更独立的研究方法。作为应用最广泛的卫星大地测量技术之一，干涉SAR(InSAR)卫星影像是获取与地震有关的同震变形的常用手段。InSAR虽然时间采样率低，视卫星重访时间长短而定，但具有非常高的空间采样率。InSAR的同震变形精度可达2～3cm，能够得到清晰的中强地震引起的地表变形结果。InSAR形变数据具有上升和下降等多种观测模式，可协助识别地表破裂及震源断层面。近年来，使用InSAR技术反演盲断层地震的发震断层参数和滑动分布已成为研究盲断层地震的重要手段。此外，使用InSAR重建的同震形变场能够获得高精度的同震观测结果，该方法已成功应用于2015年MW6.4皮山地震(Fengetal.，2016)、2016年门源地震(Lietal.，2016)以及2017年MW6.5九寨沟地震(单新建等，2017)等地震的研究中。

从表1 可看出，2个机构利用远场地震波数据给出的巴基斯坦MW6.0地震的震源机制解差异较大，并且无法得知此次地震真实的发震断层。因此，使用InSAR技术重建巴基斯坦新米尔普尔MW6.0地震的同震形变场，并反演这次地震的发震断层参数与滑动分布、深入了解其破裂特征与发震构造，对了解该区域内的逆冲断裂带活动性与发震构造机制有重要意义。此外，在新米尔普尔北部的吉拉姆河(Jhelum)上建有一座曼格拉大坝(Mangla)，该水库位于上新世—更新世西瓦利克(Siwalik)系砂岩和黏土层交替形成的低矮山丘上。西瓦利克堆积物是由河流沉积的淡水沉积物，在中新世时期开始的一次大的抬升使其恢复了活力。西瓦利克系下的默里(Murree)地层为同一地槽中新世时期的咸淡水和淡水沉积。较老的结晶岩形成了海拔较高的内喜马拉雅山脉，而较年轻的默里和西瓦利克系形成了外喜马拉雅山脉或山麓丘陵。该水库距离巴基斯坦新米尔普尔MW6.0地震的震中不到10km，因此该水库对该地震是否有触发作用是值得注意的一点。在本研究中，我们将结合水库周围的地质构造以及地震活动性对该问题加以讨论。

本研究中，我们利用欧洲航空局(European Space Agency，ESA)获取的升降轨Sentinel-1A SAR数据重建了此次地震的InSAR同震形变场，并在升、降轨形变场的联合约束下反演获得了发震断层参数和同震滑动分布，从而对这次地震造成的形变场特征、断层运动性质和发震构造机制进行分析与探讨，为进一步研究巴基斯坦北部的地震构造活动性提供一定依据。

1 InSAR同震形变观测

我们获取了2019年9月24日巴基斯坦MW6.0地震的Sentinel-1A卫星升、降轨SAR影像(TOPS模式)。震前Sentinel-1A升轨SAR影像的观测日期为2019年9月16日，震后升轨影像的观测日期为2019年9月28日； 震前Sentinel-1A降轨SAR影像的观测日期为2019年9月5日，震后降轨影像的观测日期为2019年10月11日。SAR影像干涉对的垂直基线分别为-44m(升轨)和-9m(降轨)(表2)。

表2 Sentinel-1A卫星影像参数信息Table2 Details of Sentinel-1A interferograms used in this study

基于D-InSAR技术，我们使用了GAMMA软件平台进行干涉处理。采用美国宇航局发布的SRTM 3sec数字高程模型(SRTM DEM 90m)消除地形相位误差。处理时对干涉图进行了10×2(距离向×方位向)多视处理以降低干涉图的噪声。采用最小费用流算法进行相位解缠(Werneretal.，2000； Gallandetal.，2016)，先对高质量的区域进行解缠，得到可靠的解缠相位值，然后以这些可靠相位值构建相位模型，并以此为基础对其他低相干区域进行解缠(Goldsteinetal.，1998； Farretal.，2007)。最后，将解缠后的干涉影像地理编码至WGS84坐标系下，得到升、降轨同震InSAR干涉条纹和LOS向形变场，如图2 所示。

图2 2019年12月24日巴基斯坦新米尔普尔MW6.0地震升、降轨同震形变场Fig. 2 The coseismic deformation field of the MW6.0 earthquake in New Mirpur，Pakistan，24 December 2019 observed by Sentinel-1A.红色沙滩球为USGS给出的震源机制解，红色五角星为此次地震的震中，天蓝色为水库水体的范围

根据升、降轨同震形变场的长轴方向可确定发震断层的走向为NW-SE，其结构呈非对称性，分布于喜马拉雅前缘断裂的次级逆冲断裂带北缘。断层SW盘的干涉条纹明显密集于NE盘(图2a，c)，且SW盘的形变值为正、NE盘为负(图2b，d)。同一地区的升、降轨形变量相同，图2 显示出该发震断层的运动性质以逆冲运动为主，观测获得的卫星视线向(LOS)最大隆升形变量约为0.1m(升轨)和0.05m(降轨)，最大沉降形变量约为0.03m(升轨)和0.05m(降轨)。由于形变中心北部有一座水库，震中附近失相干严重，导致形变场条纹不连续，对获取对此次地震的震源性质产生了一定影响。

2 断层滑动模型的反演

同震观测数据过多会增加计算量且降低数据的信噪比，使反演结果难以收敛。对于形变梯度较大的区域，采用均匀采样法能有效降低其对整体结果的影响。因此，在反演之前，为避免数据量过大导致计算效率降低，我们采用均匀采样方法对升降轨InSAR形变场进行降采样处理。最终，分别得到4i195(升轨)和3i311(降轨)个形变数据点用以反演发震断层面上的滑动分布。

巴基斯坦新米尔普尔MW6.0地震属于未破裂到地表的盲断层事件，发震区域的地质与地震活动性研究较少，难以从现有资料中精准地确定发震断层的位置。分析InSAR形变场可知(图2)，形变区域大致沿NW-SW向延伸，再结合震中区域的地质构造特征猜想本次地震的发震断层有2种可能，即倾向SW的断裂带或倾向NE的断裂带(MFT的次级断裂)。因此，我们建立了2种初始断层模型，并对初始断层模型的几何参数进行多次测试，搜索断层参数的最优解以反演断层面上的精细滑动分布。然后，根据观测和模拟的拟合情况以及不同发震断层面的滑动情况来推测这次地震的可能发震断层。

2.1 模型1(倾向SW)断层滑动模型及滑动分布反演

我们将升、降轨形变场的观测数据作为约束条件，根据模型1获取断层面上的精细滑动分布。首先，固定模型1的断层迹线位置，再测试不同的断层倾角对数据的拟合情况，得出当倾角为90°时拟合效果相对较好，拟合度达89%。升、降轨数据的均方根误差约为0.01m，拟合残差≤2cm。图3 为模型1的观测-模拟-残差图，由图3c、f可看出拟合度较好。模型1反演得到的最佳断层参数详见表3。得到的最佳倾角为90°，但这种情况在自然界中实际并不存在，逆冲型地震的倾角一般为30°～50°，现已知发生过的逆冲型地震的最大倾角为65°，该模型不符合逆冲型地震的客观规律，因此排除了该模型。

图3 模型1(倾向SW)的升、降轨观测值、模拟值和残差Fig. 3 Observation，simulation and residual from ascending(a—c)and descending(d—f)InSAR data and inversion by model 1.黑色虚线为模型1的断层迹线，红色五角星为此次地震震中(USGS)

表3 均匀滑动反演发震断层的几何和运动学参数(模型1)Table3 The source parameters for the Pakistan earthquake(model 1)

2.2 模型2(倾向NE)断层模型及滑动分布反演

将InSAR同震形变场作为约束条件，根据发震断层的几何参数，采用模拟退火法反演发震断层的同震滑动分布。模拟退火法是一种非线性直接反演的优化算法，能够避免结果为局部极大值而非全局最大值的情况(齐浪等，2011)。同样以升、降轨形变场的观测数据作为约束，获得模型2断层面上的精细滑动分布。采用LOG函数确定最佳倾角和平滑因子(冯万鹏等，2010)。首先，根据USGS给出的震源机制解固定NE倾的断层倾角，并不断平移断层的位置，将拟合残差最小的断层位置作为最佳发震断层的位置，然后再固定该发震断层，在不断调整断层倾角的同时改变平滑因子，并计算每次反演结果的均方根误差(root-mean-square error，RMSE)(式(1))。设定的倾角范围为2°～60°，以2°为步长，设定平滑因子的范围为0～0.3，以0.05为步长，所得结果如图4 所示。

(1)

其中，zobs为观测值，zmod为模拟值(冯万鹏等，2010)。

根据反演结果来看，模拟退火法在基于InSAR数据计算发震断层参数方面的效果较好，所得结果基本满足了反演需要。

图4 倾角与RMS的关系图Fig. 4 Diagram of inclination and RMS.a 倾角与平滑因子的曲线关系； b 倾角与平滑因子的函数关系

图5 模型2(倾向NE)的升、降轨观测值、模拟值和残差Fig. 5 Observation，simulation and residual from ascending(a—c)and descending(d—f)InSAR data and inversion by model 2.黑色虚线为模型2的断层迹线，红色五角星为此次地震震中(USGS)

图6 模型2(倾向NE)断层面上的滑动分布Fig. 6 Slip distribution of the northeast-dipping fault model 2.

表4 均匀滑动反演发震断层的几何和运动学参数(模型2)Table4 The variable source parameters for the Pakistan earthquake(model 2)

由图4 可看出，当倾角约为10°～20°平滑因子约为0.05时拟合效果最好，拟合度达89.37%。图5 为升、降轨拟合效果图，无论是升轨还是降轨，拟合得到的形变场都有明显的一升一降盘，拟合误差≤4cm，符合一般逆冲型地震的特点。根据所确定的断层最佳参数，画出模型2的断层滑动分布图，如图6 所示，滑动整体呈椭圆形集中分布。滑动分布沿走向集中在10～14km范围内，沿发震断层面倾向主要集中在10～20km内，滑动分布表现为逆冲性质。该模型反演得到的最佳断层参数详见表4。分析认为，该断层模型与此次地震的发震断层特征更相符。但反演得到的矩震级与USGS给出的震级差距较大，这可能是由于此次地震的震级及形变区域的形变量太小，未能被InSAR测量捕捉。大气延迟、轨道、DEM等误差也对反演结果产生了影响。

3 讨论

基于此次新米尔普尔MW6.0地震的InSAR同震形变场数据，反演计算得到的震中位置为(33.09°N，73.75°E)，与USGS确定的结果相差4.48km，与GCMT所得结果相差68.65km。经分析认为，USGS和GCMT给出的发震地点和震源机制解是基于远场地震波资料反演得出的，但由于台站资料具有不确定性，导致其给出的震中位置差异较大。USGS等机构在确定震源深度时主要利用了地震波的走时，但这类方法对地震台的间距要求很高，只有当最小震中距与震源深度达到一定比例时，所得的震源深度才能足够准确。对于6级以上地震，远震体波能够较准确地测定出震源深度，但对于5级以下的中小地震而言，远场台站数据信噪比较低，不利于准确测定震源深度。

总体而言，传统的地震学方法在确定震源参数时对台站的分布密度及观测数据的质量等要求很高。对于中小地震，地震学方法费时费力，且结果误差较大，不适合进行震源参数研究。InSAR是一种新兴的大地观测技术，具有分辨率高、形变敏感度高且覆盖范围广等特点，以InSAR数据为基础可利用弹性位错模型获取震源参数，所反演的断层参数和滑动分布结果更具可靠性。中小地震的形变量较小，而InSAR测量技术的精度很高，能识别出mm级的形变。InSAR空间分辨率较高，能对发震位置起到很好的约束作用。此外，近年来发展后的时序InSAR技术能很好地处理震前、震后时间序列的形变场，在很大程度上提升了中小地震地面形变的观测精度，为研究中小地震断层面三维几何形态和力学属性提供更加可靠的地面位移数据，为研究发震区域的活断层的运动性质和地震周期提供更加清晰的认识。

图7 同震滑动矢量与地形(a)、重力(b)的分布关系图Fig. 7 Diagram of co-seismic sliding vector(a)and gravity(b)distribution.黄色沙滩球是本次地震震源机制解，蓝色和黑色箭头均为滑动矢量

新米尔普尔MW6.0地震也为研究青藏高原南缘的大规模构造提供了机会。震中区域处于欧亚板块与印度板块碰撞区，印度板块持续向N推挤，使得该碰撞区发生变形。由图7 可知，欧亚板块与印度板块之间的相对运动方位角虽然变化不大，但逆冲滑动方位角确实发生了变化。因此，我们推测除了驱动板块边界运动的力之外，重力也控制着逆冲运动的方向(Copley，2012)。当山地与低地之间存在重力势能差则意味着山地会变形以降低势能，从而导致横向扩展至刚性下边界之上，这种响应重力驱动力的横向运动被称为重力流(Copleyetal.，2015)，在重力驱动力的作用下地壳向垂直于地形等高线的下坡方向移动，且导致地壳沿青藏高原南缘的弧形山脉径向逆冲。因此，此次地震的滑动矢量结果表明，重力驱动力对于控制青藏高原南缘的弧形边缘而言十分重要(Copleyetal.，2007)。

新米尔普尔MW6.0地震北部约10km的吉拉姆河上建有一座中型水库及曼格拉大坝。该水库区域的地质构造为基底岩石上覆盖着厚约3km的西瓦利克系易碎砂岩和黏土，这些岩石对施加的荷载几乎没有抵抗作用，并且会逐渐变形以适应荷载，这样的地质构造特征导致该水库没有出现明显的地震效应，在水库发生过的几次大地震都与其东北部的断层有关。另外，在此次地震的发生地新米尔普尔建有地震观测站，对蓄水前、后记录到的地震事件的频率进行分析，发现蓄水导致地震活动性出现明显增加，然而经进一步分析发现，地震的增多实则是因为观测仪器的灵敏度提高了，而非水库蓄水触发了地震(Brown，1974)。Brown(1974)认为未来地震活动性可能由于构造力的累积再次增加，而非由于水库水位波动引起的荷载变化。综合水库区域的地质构造分析以及古地震与历史地震的研究成果，我们认为此次地震并非由水库蓄水触发的。

4 结论

本文以InSAR升、降轨数据为约束，重建了2019年9月24日巴基斯坦MW6.0地震同震形变场，通过对不同断层模型的反演试验，给出了发震断层可靠的几何参数和滑动分布特征，并得到了以下结论：

基于2019年巴基斯坦地震的升、降轨SAR数据，得到了该地震的同震形变场，其形态呈现明显的一升一降两盘分布，显示此次地震为逆冲型，兼具少量右旋走滑，LOS向最大隆升量为0.1m，最大沉降量为0.05m。巴基斯坦MW6.0地震没有出现地表破裂位置的失相干条带，这可能是由于此次地震发生在低倾角的印度板块和欧亚板块俯冲带上，而板块俯冲带上逆断层的运动和破裂方式都与板内的逆断层不同。断层面的倾角较缓，断层破裂时上盘向S仰冲，北盘向N俯冲，由于挤压推覆构造造成了上盘的前端隆升，而其后端受到拉张成为沉降区(Savage，1983； Hyndmanetal.，1993； 屈春燕等，2017)。仅根据InSAR数据无法准确得到发震断层的倾向，因此我们建立了倾向SW和倾向NE 2种断层模型，并分别进行多次反演实验，再结合USGS等机构给出的震源机制解及巴基斯坦北部区域的地质构造背景信息，最终认为倾向NE的断层模型与此次地震的发震断层特征更相符。该模型表明，2019年巴基斯坦地震的断层模型为主前缘逆冲断裂(MFT)的次级断裂。断层走向为304°，倾角为15°，最大滑动量为0.64m，平均滑动角为124.85°，震源深度为2.77km，反演得到的矩震级为MW6.0。

此次地震产生的破裂并未出露地表，是一次盲断层型的低角度逆冲活动。根据反演结果，分析认为该地震是由喜马拉雅主前缘逆冲断裂的次级断裂引起，再结合区域地貌特征推测此次事件是印度板块向欧亚板块持续低角度俯冲所造成的。东部西太平洋-菲律宾板块向W的俯冲作用、西部印度板块向N的俯冲碰撞作用使得中国大陆地壳强烈变形。汪素云等(1996)发现造成这种地壳强烈变形的主要动力来源于印度板块，印度板块与欧亚板块的碰撞作用使中国大陆地壳呈现出 “西强东弱”的形变特点，大多数的地震活动集中在青藏高原及其邻区(吴中海等，2016)。位于青藏高原南部的喜马拉雅造山带是印度板块与欧亚板块碰撞的产物，其规模巨大，且正持续高速运动，由此导致高频率、大强度地震的不断发生。因此，了解喜马拉雅造山带上的地震活动特征有利于研究青藏高原的地壳形变特征，同时也有助于探索中国未来的大地震潜在的危险性以及活动趋势，助力中国的防震减灾事业。

致谢欧洲航空局(ESA)为本研究提供了Sentinel-1A雷达卫星影像； 德国GMZ的汪荣江老师提供了断层滑动分布反演SDM程序包； 本文部分图件使用GMT软件绘制。在此一并表示感谢!