引言

2017年11月12日当地时间21时18分伊拉克哈莱卜杰（Halabjah）地区发生了M_W7.3强震，震中坐标为（34.911°N，45.959°E），震源深度为19 km （USGS，2017），伊朗受损程度尤为严重。据法新社报道，截至北京时间14日20时，地震已经造成600余人死亡，多达6 000人受伤，7万人流离失所。

2017年哈莱卜杰M_W7.3地震是继2011年土耳其凡城M_W7.1地震之后发生在伊朗高原扎格罗斯造山带的又一次强震（图1）。此处地质构造复杂，恰好位于伊拉克、伊朗交界处的扎格罗斯造山带南侧，属于阿拉伯—欧亚大陆碰撞带。由于阿拉伯板块以23 mm/a的扩张速率向欧亚板块运动，导致该板块东北缘快速隆升并形成了扎格罗斯造山带。该造山带平均海拔约为3 km，由伊朗西部自南东向北西延伸至伊拉克北部，是目前阿拉伯板块—欧亚板块碰撞区域地震活动性最强的区域之一（Agard et al，2006；Elliott et al，2013；Madanipour et al，2013）。由于阿拉伯板块和欧亚板块分别以走滑和逆冲的运动形式向北东斜向汇聚碰撞，使得扎格罗斯造山带迥异于典型的对称式结构，表现为不对称式结构的构造岩浆带（张洪瑞，侯增谦，2015）。两大板块相互挤压碰撞造成岩石内部的应力积累，当应力积累到临界状态就会使岩石发生错动，从而引发地震。此次哈莱卜杰地震的主震位于扎格罗斯主前缘断层（Tavani et al，2018）东侧区域，大部分余震集中分布在扎格罗斯主前缘断层附近区域，因此推测哈莱卜杰地震的发震断层为扎格罗斯主前缘断层。地震发生后，美国地质调查局（USGS，2017）和哈佛大学（GCMT，2017）均对这次地震的震源参数进行了反演，两个机构公布的震源机制解结果（表1）均反映出该地震的运动特征为低角度逆冲兼走滑。

表  1  2017年11月12日伊拉克哈莱卜杰M_W7.3地震的震源机制解

Table  1.  The focal mechanism solutions of the M_W7.3 Halabjah，Iraq，earthquake on November 12，2017

机构	MW	东经/°	北纬/°	走向/°	倾角/°	滑动角/°	深度/km
USGS （2017）	7.3	45.959	34.911	351	16	137	19.0
GCMT （2017）	7.4	45.840	34.830	351	11	140	17.9

下载: 导出CSV 

| 显示表格

图  1  哈莱卜杰地区地质构造背景图

Figure  1.  Tectonic setting map of Halabja region

下载: 全尺寸图片 幻灯片

InSAR技术作为一项新的对地观测技术，因高时空分辨率、大范围、高精度等优点，已然成为大面积区域形变监测的主要手段，并被广泛应用于地表形变提取和震源机制研究中（张国宏等，2011；许才军等，2012；Hetland et al，2012；单新建等，2017；赵强等，2017）。本文拟利用欧洲空间局（European Space Agency，缩写为ESA）升、降轨Sentinel-1A卫星数据快速获取震区的同震形变场，并基于Okada模型反演得到发震断层面的滑动分布，以对哈莱卜杰地震造成的形变场特征和运动性质进行深入分析。

1.   InSAR数据处理及形变场分析

1.1   数据选取与处理

哈莱卜杰地震发生之后，欧洲空间局陆续发布了覆盖地震区域的Sentinel-1A卫星宽幅干涉模式（interferometric wide swath，缩写为IW）观测影像。宽幅干涉模式采用渐进的条带扫描方式（terrain observation with progressive scans synthetic aperture radar，缩写为TOPSAR）来获取3个子条带，进而合成一景大的图像。相对于传统的条带扫描模式而言，TOPSAR工作模式不仅能获取覆盖面积广阔的地面影像，还能提高影像分辨率。鉴于Sentinel-1A卫星宽幅模式数据具有覆盖范围广、分辨率高、重访周期短、轨道精度高、时空基线间隔小等优点，在很大程度上保证了干涉图的相干性，因此选取覆盖震区的两对宽幅模式下的升、降轨影像进行干涉处理。具体的影像数据信息列于表2。

表  2  雷达影像对参数

Table  2.  Parameters of SAR image pair

卫星	震前日期	震后日期	时间间隔/d	垂直基/m	轨道	模式
Sentinel-1A	2017−11−07	2017−11−19	12	36	6	升轨
Sentinel-1A	2017−11−11	2017−11−23	12	−52	72	降轨

下载: 导出CSV 

| 显示表格

获取卫星影像数据以后，以震前影像为主影像，震后影像为辅影像，以瑞士的GAMMA数据处理软件为平台，采用二轨差分法对影像数据进行干涉处理得到哈莱卜杰震区的视线向同震形变场。二轨差分法的基本思想是：先利用一对雷达影像数据进行干涉处理，得到包含形变相位和地形相位的干涉图，再将已知的数字高程模型（digital elevation model，缩写为DEM）数据模拟成雷达坐标系下的地形相位，然后从干涉图中去除地形相位，得到只包含形变相位的差分干涉图，最后进行相位解缠和地理编码，得到哈莱卜杰地震的同震形变场。

由于获取的Sentinel-1A卫星数据为单视复数影像（single look complex，缩写为SLC），可将其直接导入GAMMA软件进行干涉处理。在确定主从影像后，将多视比设为10:2，得到去平后的初始干涉图；然后将地理坐标系下已有的90 m分辨率的DEM （Farr et al，2007）数据模拟成雷达坐标系下的地形相位图，从初始干涉图中去除地形相位得到差分干涉图；为了滤除干涉图中的噪声以提高干涉图的质量，采用Goldstein滤波方法（Goldstein，Werner，1998）进行两次滤波，滤波窗口分别设为64和16；解缠时相干性阈值设为0.5，用最小费用流法（Pepe，Lanari，2006）避开相干性低的区域，生成掩膜图后提取高相干性的点，建立狄洛尼三角网对滤波后的干涉图进行相位解缠以得到真实的相位信息；最后根据相位与形变量之间的关系将相位图转为形变图，并进行地理编码，得到地理坐标系下的同震形变场，如图2和图3所示。

图  2  2017年哈莱卜杰M_W7.3地震降轨视线向同震形变场

（a） 降轨干涉条纹图；（b） 降轨干涉形变图

Figure  2.  The LOS descending coseismic deformation fields of the 2017 M_W7.3 Halabjah earthquake

（a） The descending interferogram；（b） The descending displacement map

下载: 全尺寸图片 幻灯片

图  3  2017年哈莱卜杰M_W7.3地震升轨视线向同震形变场

（a） 升轨干涉条纹图；（b） 升轨干涉位移图

Figure  3.  The LOS ascending coseismic deformation fields of the 2017 M_W7.3 Halabjah earthquake

（a） The ascending interferogram；（b） The ascending displacement map

下载: 全尺寸图片 幻灯片

1.2   同震形变场结果分析

从图2和图3中可以看出升、降轨形变场的覆盖范围大约为60 km×70 km，具有多个不对称的椭圆形变区，基本沿扎格罗斯主前缘断层东侧展开分布。降轨形变场有两个椭圆形变区域，东北盘呈现3个半条纹变化，最大形变值为−37 cm，表现为沉降趋势；西南盘呈现5个半条纹变化，最大形变值为55 cm，表现为隆升趋势。而升轨形变场表现为4个形变区域，西南盘出现了9个同心圆条纹，远比东北盘3个形变区域的条纹要密集，说明西南盘的形变梯度大，最大形变值为88 cm，呈抬升状态；东北盘总体存在3块条纹区域，条纹最多的区域位于东北盘最上方，存在1个半同心圆条纹，其最大的地表形变值为−16 cm，总体呈下降状态，与降轨呈现的两个形变区域的差异是由不同的入射角和方位角造成的。综合以上分析认为升降轨形变结果均呈现出不对称性，且地震造成的地表隆升形变值远大于沉降值，结合断层运动的特征，可以判断哈莱卜杰地震是由逆断层引起的。

为了进一步分析哈莱卜杰地震同震形变场的空间分布特征，分别从升、降轨同震形变场提取一条穿过形变中心的剖面线XX ′，如图4所示。结果表明升降轨的形变剖面线XX ′均有明显的抬升中心，最大形变值分别约为52 cm和80 cm，远离形变最大值的两侧区域，升降轨LOS向形变值均逐渐减小。不同成像模式会导致形变结果在空间分布和量级上的不同，但升降轨形变剖面线变化趋势基本一致，均以隆升为主，体现了逆冲型地震造成的地表形变特征。

图  4  2017年哈莱卜杰M_W7.3地震升降轨视线向同震形变剖面

（a） 升轨形变剖面；（b） 降轨形变剖面

Figure  4.  The LOS deformation profiles of coseismic for the 2017 M_W7.3 Halabjah earthquake

（a） The ascending deformation profiles；（b） The descending deformation profiles

下载: 全尺寸图片 幻灯片

2.   断层滑动反演

2.1   形变场重采样

图2和图3获取的视线向同震形变场的相干性比较好，影像范围也完整地覆盖了形变区域，保证了反演的可靠性。考虑到InSAR数据分辨率高的特征，如果将形变场所有的观测数据都作为反演数据，需要计算超大数据量的同时也会影响反演结果的精度，因此需要对视线向形变场进行降采样。现在常用的降采样方法有均匀降采样和四叉树降采样，鉴于该区域的地表形变梯度大，更适合采用四叉树方法进行重采样。本文选用四叉树降采样方法（Lohman，Simons，2013），该方法在最大程度上保留原有形变场总体特征的同时还可以保证形变梯度大的区域采样点数量多于形变梯度小的区域。

2.2   均匀滑动反演

为了得到此次地震的震源机制解，以升、降轨LOS向同震形变场作为约束，首先基于Okada弹性半无限平面空间位错理论模型（Okada，1985）忽略地球介质分层的影响，把发震断层当作一个光滑的矩形平面，然后采用非线性方法利用均匀滑动模型来获取哈莱卜杰地震的断层几何参数，具体包括断层走向、倾角、经纬度、长度、宽度、深度、滑动量等。本文采用单纯型算法进行均匀反演，其目的是使模拟形变值与观测值达到最佳拟合，从而得到最佳的断层几何参数。在反演过程中需要给出断层长度、宽度以及深度等约束条件（上、下限）来确保这些参数具有物理意义，其余参数不进行约束。

为了确定断层的几何参数，本文采用psokinv程序（冯万鹏，李振洪，2010；Feng et al，2013，2014）进行均匀反演。该程序基于多峰值粒子群算法（multi-peak particle swarm optimization，缩写为MPSO）来搜索参数的最优解。多峰值粒子群优化算法的特点是利用个体在解空间中的随机速度来改变个体，该算法收敛速度快、稳健性明显。反演得到的断层走向为350.9°，倾角约为20°，滑动量为3.38 m。根据矩震级公式，得到该地震的地震矩为7.91×10¹⁹ N·m，矩震级为7.2。

根据均匀反演得到的发震断层几何参数进行正演，得到InSAR同震拟合形变场和残差图（图5）。从图5中可以看出，对单一断层进行均匀滑动反演得到的结果与观测值拟合得较好，升、降轨的拟合形变残差值大部分均在0左右，仅形变中心处少量残差达到了±10 cm，这是因为形变中心处形变梯度大，断裂情况复杂。

图  5  均匀滑动模型拟合的形变图和残差图

（a） 降轨形变；（b） 降轨拟合形变；（c） 降轨残差；（d） 升轨形变；（e） 升轨拟合形变；（f） 升轨残差

Figure  5.  The fitted and residual displacement diagram based on the uniform slip model

（a） The displacement diagram for descending；（b） The fitted displacement diagram for descending；（c） The residual displacement diagram for descending；（d） The displacement diagram for ascending；（e） The fitted displacement diagram for ascending；（f） The residual displacement diagram for ascending

下载: 全尺寸图片 幻灯片

2.3   分布式滑动反演

基于均匀反演的断层几何参数，将断层面划分为均匀、连续排列的矩形子单元进行分布式滑动反演，即可获得断层面上的精细滑动分布。

本研究利用的反演算法为基于最速下降法（steepst descent method，缩写为SDM）（Wang et al，2013），在断层的几何参数确定的情况下，断层面上的运动参数与形变观测数据就转化为一般线性问题，其求解公式为

式中，d为视线向InSAR观测数据，G为将位错模型与形变观测数据联系起来的格林矩阵；m为描述地下断层走向和倾向的滑动量，ε为形变数据观测误差。

本文根据均匀反演给出的断层位置、长度、宽度和走向等参数，将断层走向长度延至90 km，倾向长度延至80 km，然后将断层按5 km×5 km离散成18×16个子断层。考虑到均匀反演模型中得到的倾角并不是最佳倾角，这里将倾角范围设置为10°—25°，每隔1°进行多次试验，最终得出当倾角为17°时，残差最小（图6）。在求解过程中为避免滑动分布解出现震荡现象添加拉普拉斯平滑约束条件（Feng et al，2014），即

图  6  不同倾角与残差的权衡曲线

Figure  6.  Trade-off curve between the dip and misfit

下载: 全尺寸图片 幻灯片

式中，d为视线向InSAR观测数据，s为地下断层面的滑动量，α为用于权衡观测数据拟合度和反演结果粗糙度的平滑因子，H为拉普拉斯算子。

为了进一步保证计算结果的合理性，需要引入合理的平滑因子α。通过模型反演结果的粗糙度和数据拟合度之间的折中曲线（图7），最终确定α为0.06。

图  7  模型粗糙度与残差的权衡曲线

Figure  7.  Trade-off curve between the model roughness and misfit

下载: 全尺寸图片 幻灯片

经SDM方法反演计算，最终获取的断层面滑动分布结果如图8所示。反演结果表明，断层滑动主要集中在沿倾向20—42 km和沿走向25—50 km处，最大滑动量为3.34 m，位于地下约21 km处，平均滑动角为142.1°，平均滑动量为0.78 m，断层整体呈以逆冲为主兼少量右旋走滑的运动特征。图8a右下角可以看到一个明显的滑动区域，这可能是由于本研究采用单一断层模型所致。反演得到的地震矩为1.68×10²⁰ N·m，矩震级M_W为7.4，与表1两个机构给出的矩震级非常接近。图9是由反演的滑动结果模拟得到的升、降轨形变残差图，从图中可知，InSAR观测形变值与模拟值在空间分布和量级上基本一致，拟合残差大部分集中在±5 cm左右，说明本文采用的断层模型可靠，降轨形变中心附近的残差大部分达到了±15 cm，升轨形变中心附近的残差远小于降轨，仅在边缘处少量达到了±15 cm，可能是由升轨影像的飞行方向与断层走向几乎平行，而降轨影像的飞行方向与断层走向交叉所造成，另外也与大气误差、断层模型的简单化有关。

图  8  断层滑动分布显示

（a） 滑动分布模型；（b） 地理坐标系下的滑动分布

Figure  8.  Fault slip distribution from InSAR inversion

（a） Slip distribution model；（b） Slip distribution in geographic coordinates

下载: 全尺寸图片 幻灯片

图  9  滑动分布模型拟合的形变图和残差图

（a） 降轨形变；（b） 降轨拟合形变；（c） 降轨残差；（d） 升轨形变；（e） 升轨拟合形变；（f） 升轨残差

Figure  9.  The fitted and residual displacement diagram based on the slip distribution model

（a） The displacement diagram for descending；（b） The fitted displacement diagram for descending；（c） The residual displacement diagram for descending；（d） The displacement diagram for ascending；（e） The fitted displacement diagram for ascending；（f） The residual displacement diagram for ascending

下载: 全尺寸图片 幻灯片

3.   讨论与结论

本文利用二轨差分法（D-InSAR）对Sentinel-1A卫星获取的两对雷达影像进行干涉处理，分析了2017年伊拉克哈莱卜杰M_W7.3地震的同震形变场特征，同时基于弹性半平面空间矩形位错模型反演了此次地震的断层滑动分布，得到以下结论：

1） 哈莱卜杰震区同震形变场的分布与余震的分布一致，均集中在扎格罗斯主前缘断层的东侧，说明哈莱卜杰地震与该断层有关，形变场影响范围约为60 km×70 km，升轨视线向形变最大值为88 cm，降轨视线向形变最大值为55 cm，整体表现为隆升，符合逆冲为主的地震运动特征。

2） 基于均匀滑动和分布滑动模型模拟得到的形变场与真实形变场拟合良好，残差基本控制在±5 cm以内，均匀滑动模型得到的发震断层线与扎格罗斯主前缘断层接近平行，说明此次地震的发生与扎格罗斯主前缘断层的构造运动密切相关。

3） 联合升、降轨InSAR形变数据来反演哈莱卜杰地震的滑动分布。反演结果显示此次地震的地震矩为1.68×10²⁰ N·m，矩震级为M_W7.4，发震断层的走向为350.9°、倾角为17°，且以逆冲为主，兼具少许右旋走滑，最大滑动量为3.34 m，位于地下21 km处，滑动量主要集中在沿断层倾向20—42 km和走向25—50 km处。

欧洲空间局（ESA）提供了Sentinel-1A影像数据，德国地学中心的汪荣江教授提供了SDM反演程序，作者在此一并表示感谢。