2017年11月12日，伊朗伊拉克边境附近的萨尔波勒扎哈卜（Sarpol-e Zahab）地区发生Mw7.3级强震，震中位于34.886°N，45.941°E，震源深度为19 km^[1]，此次地震造成了严重的人员伤亡^[2]。由于地震发生在岩石层而不是松散的土壤厚层，因此地震波在传播中被放大得相对较小，故地震发生后，传到地表的震动也较少，没有导致发震断层在地表的破裂^[2]。为研究发震断层性质，众多学者对此次地震进行了研究。Kobayashi等^[3]和Wang等^[4]主要研究了地震与滑坡之间的关系，使用合成孔径雷达干涉测量（interferometry synthetic aperture radar，InSAR）方法获取了地震的形变场，但没有详细介绍地表破裂；He等^[5]、Ding等^[6]和Vajedian等^[7]基于合成孔径雷达（synthetic aperture radar，SAR）数据分析了地震造成的地表破裂特征并获取了断层面解。因此，使用基于高精度的大地测量方法获取的形变场可以对断层参数进行反演，能更好地分析地震构造环境。

目前，一般使用InSAR技术获取地震形变场，并基于同震形变场对断层参数进行反演分析。常用的反演方法有模拟退火法^[8]、最速下降法^[9]和随机搜索粒子群优化（particle swarm optimization，PSO）算法^[10]，前两种方法发展较早，但算法相对复杂，对参数要求较高。PSO算法可以基于InSAR、GPS结果进行约束，再选择最适合的反演路线对断层参数进行反演。本文使用ALOS-2和Sentinel-1卫星的升降轨数据分析地震同震形变场，基于PSO算法对震源机制进行反演，获得了萨尔波勒扎哈卜地震的详细震源参数，并分析了地震的断层几何特征；利用多源SAR数据的同震形变场获取了同震三维形变场，直观地展示了此次地震导致的地壳变形，为理解地震运动构造特征和断层性质等提供基础信息。

1 地震构造背景

伊朗西南部的扎格罗斯造山带是世界上发生地震最多的造山带之一，扎格罗斯造山带的形成与阿拉伯板块和欧亚大陆（包含伊朗）最北部的板块缓慢碰撞相关。扎格罗斯造山带是目前地球上最年轻的陆陆碰撞造山带，属于新特提斯碰撞造山带的一部分^[11-12]。扎格罗斯造山带位于伊朗西部，延伸到伊拉克北部，是长达1 500 km的褶皱冲断带，平均海拔约1 000~2 000 m。2017年11月12日的萨尔波勒扎哈卜地震位于扎格罗斯造山带西北部的逆冲断层（图 1）。由于阿拉伯板块以20~25 mm/a的速度向欧亚大陆北移^[13-14]（图 1中玫红色箭头），两个板块在该次地震附近的西北碰撞板块边界汇合，推动伊朗扎格罗斯山脉抬升，碰撞影响范围覆盖伊朗全境。板块碰撞的能量沿着或靠近这些边界造山带，以地震形式进行释放，故导致了本次地震的发生^[15-17]。如图 1所示，HZF（high Zagros fault）表示高扎格罗斯断层，MZF（main Zagros fault）表示主扎格罗斯断层，MFF（mountain front fault）表示山前断层。

虽然伊朗、伊拉克为地震常发区，但在过去20多年里，该区域未发生震级M≥4.5的地震，自1900年以来，只在1957年和1958年发生过两次震级为M6.7和M6.8的大地震，前两次大地震与本次地震震中相距约200 km，因此2017年11月12日的Mw7.3级地震为扎格罗斯山脉有史以来最大的一次地震^[18]。

2 InSAR同震形变场 2.1 数据处理

萨尔波勒扎哈卜地震发生后，在欧空局网站可以搜索到Sentinel-1卫星在地震前后获取的SAR影像，其中包括T72升轨影像和T79降轨影像。同时为了更精确地研究地震的震源机制，从日本宇宙航空研究开发机构申请了此次地震ALOS-2卫星的升降轨数据。图 1显示了本文使用SAR影像覆盖的区域及震中位置，可以看出，两种卫星的SAR数据都可以完整地覆盖震中区域。表 1为本文使用的SAR数据的具体参数。

获取SAR数据后，对ALOS-2和Sentinel-1数据进行差分干涉处理，获得地震同震形变场。使用GAMMA软件为InSAR处理平台，利用分辨率为90 m的数字地形高程模型先进行去地形处理，然后进行滤波和相位解缠，并对轨道误差进行校正，最后得到完整的差分干涉相位图，进行地理编码后可以获得同震形变场图。

2.2 同震形变场特征

图 2为ALOS-2和Sentinel-1卫星SAR传感器获取的降轨和升轨同震形变场，可以看出，两种形变场结果的一致性较高。升降轨中，干涉图的椭圆形变区不同主要是由卫星方向不同所导致。对比同为升轨（A）的ALOS-2 T180和Sentinle-1 T72，可以发现，T180左盘干涉条纹中存在失相干区域，可能是由于ALOS-2 T180数据的时间基线较长导致，但两个升轨数据的形变干涉区域和形变量都相差不大。同为降轨（D）的T71和T79，其左右两盘形变干涉区域则具有较高的一致性。同时，ALOS-2 T180获取的形变场远场有部分误差存在，可能是由于差分干涉处理中，电离层延迟误差没有完全去除，但远场误差与形变场整体量级对比相对较小，对后续震源机制反演等分析处理的影响可以忽略。通过分析不同卫星升降轨的形变场可以发现，萨尔波勒扎哈卜地震的地面形变主要有两个形变区域，卫星视线向（line of sight，LOS）的最大位移量达88.8 cm，距离断层两侧位置越远，LOS向形变量逐渐减小，远场的形变量较小，这种变化特征与地震活动物理机制一致。

3 PSO模型反演 3.1 震源机制反演方法

本文以ALOS-2和Sentinel-1卫星数据获取的萨尔波勒扎哈卜地震的同震形变资料为反演的约束条件，对地震的断层参数和滑动机制进行详细的分析研究。为提高反演效率，需要对InSAR形变干涉图进行降采样处理，获取大小合适的InSAR图集，其降采样方法为数据分辨率约束的四叉树方法^[19]。降采样之后，ALOS-2升轨T180和降轨T71的点数分别变为241和198，Sentinel-1升轨T72和降轨T79的点数分别变为137和151。采用二步反演算法估计断层的破裂几何参数以及倾滑分布^[20]，具体流程为：首先假设断层为均一断层模型，计算出断层的空间几何参数；然后使用分布式断层模型计算断层面上的分布式滑动量。利用PSOKINV软件反演震源参数，该软件使用改进的群体协作的PSO算法^{[10, 21-23]}，该算法改进了适配函数，简化了非线性求解空间的维度，通过改进最优解获取方法，提高了断层参数反演中的收敛率。

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

断层的几何参数主要包括断层位置（经度和纬度）、顶部埋深、走向、倾角等；断层运动参数主要有滑移量和滑动角，断层几何参数和运动参数是理解地震同震形变的主要依据之一。本文首先使用均一断层模型获取断层的几何参数，以求得理论形变量与真实观测值的最优拟合。在均一断层模型中，使用单断层模型进行反演，基于Okada弹性位错模型^[24]使用PSO非线性优化算法自动对模拟结果进行比较，从而搜索出最优参数，求得可以使适配函数达到最小值的解。本文使用的适配函数定义为：

$ \sigma = \sqrt {\frac{{{{\left( {\mathit{\boldsymbol{W}}\left( {\mathit{\boldsymbol{D}} - \mathit{\boldsymbol{GS}}} \right)} \right)}^2}}}{N}} $

(1)

式中，G为系数矩阵，表示均匀断层倾滑量为1 m时引起的地表运动响应；S表示倾滑矢量；D为地表形变观测值；W为每个数据集的相对权重；N为形变观测值的个数。将获取的InSAR同震形变场作为约束条件，使用PSO方法在全参数域内进行搜索，找到可以使适配函数为最小值的解。初步断层几何参数反演的结果为：断层的最优拟合断层参数的经纬度分别为45.748°和34.716°，断层走向为347°，断层的倾角是15°，震级达到Mw7.36，详见表 2。

为进一步确定地震滑动的空间分布，将由均一断层模型确定的断层沿走向和倾向方向进行扩展，利用分布式断层模型分析地震滑动空间分布。该算法考虑了模型的空间分辨率，基于深部相关算法，浅层子断层使用较小的单元尺寸，深部子断层的单元尺寸随着深度线性增长。本文选取1 km×1 km的单元大小。滑动分布反演主要是线性求解问题，因此需要使用非负约束的最小二乘方法进行求解^[21]。使用二次差分拉普拉斯算子约束滑动粗糙度，可以防止子断层块之间存在过大抖动的滑动^[28]，基本的关系方程可以表示为：

$ \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{G}}_\delta }}\\ {{\alpha ^2}\mathit{\boldsymbol{L}}} \end{array}} \right]\mathit{\boldsymbol{S}} = \left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{D}}\\ 0 \end{array}} \right] $

(2)

式中，L为二阶微分算子，用于评价滑动粗糙度；${{\alpha ^2}}$为平滑因子，用来平滑解析解；${{\mathit{\boldsymbol{G}}_\delta }}$与式（1）中的定义相同，表示通过弹性位错模理论计算出来的倾角δ的理论位移场；S为每个子断层的倾滑矢量；D为形变观测值。均一断层模型可以确定最优拟合解，但在滑动分布反演过程中，最优倾角会有一些细微的变化。因此，滑动模型使用迭代方法以进一步优化确定最优倾角。

根据模型定义，残差ξ和模型的粗糙度ψ可以分别定义为^{[20, 29]}：

$ \xi = \sqrt {\frac{{{{\left( {\mathit{\boldsymbol{W}}\left( {\mathit{\boldsymbol{D}} - \mathit{\boldsymbol{GS}}} \right)} \right)}^2}}}{N}} $

(3)

$ \psi = \sum\limits_{i = 1}^n {\left| {{p_i}} \right|} /\left( {2n} \right) $

(4)

式（3）中，ξ为残差，其他变量与式（1）相同。式（4）中，p_i为p的分量，，p = LS；n为子断层的个数。模型进行归一化后，粗糙度曲线是单调递减的函数，残差曲线为单调递增的函数，但通过这两条曲线无法获取全局最小值，因此构建对数函数$f\left( {\delta , \alpha } \right) = {\rm{lg}}\left( {\psi + \xi } \right)$，通过以下4步获取最优倾角和平滑因子：

1) 使用有界变量最小二乘法估计给定倾角下α²的解^[30]，从而恢复出拟合残差（ξ）和模型粗糙度（ψ）。需要注意，α²的振幅与式（2）中的L有关，可以设定一系列间隔为0.5、[0.1, 10]范围内的α²。

2）给定倾角范围为[10°, 20°]，采用迭代的方法，以1°为间隔，计算该倾角范围内的拟合残差（ξ）和模型粗糙度（ψ）。

3）使用形如$\left( {\left\{ \cdot \right\} - {\rm{min}}\left\{ \cdot \right\}} \right)/\left( {{\rm{max}}\left\{ \cdot \right\} - {\rm{min}}\left\{ \cdot \right\}} \right)$的公式对ξ和ψ进行归一化。其中，{·}表示一系列的ξ和ψ，且这两个变量都是无量纲的。

4）任意给定(δ，α)，计算f(δ，α)，并绘制f曲线，从而确定最优倾角和平滑系数。

3.3 地震震源反演结果

图 3为最优滑动分布结果，从反演结果可以看出，断层走向为347°，倾角为15°，倾滑主要分布在深度为16~19 km范围内，断层倾滑以走滑兼逆冲为主，最大倾滑量为4.5 m。反演结果与全球矩心矩张量（global centroid moment tensor catalog，GCMT）^[25]及美国地质勘探局（United States Geological Survey，USGS）^[1]的结果一致，但GCMT和USGS只是快速给出震源机制的主要参数，无法确定断层面的大小，使用二步法进行震源机制反演时，是依据USGS给出的初步震源参数获取更详细的震源机制解和较为精确的断层参数，为后续的形变分析和发震机理分析等提供基础。详细震源参数见表 2。

本文分别基于最优拟合滑动模型给出ALOS-2和Sentinel-1卫星降轨和升轨数据的拟合结果。图 4为各组数据基于InSAR结果反演模拟得到的干涉形变场，图 5为相应的残差分布。可以看出，ALOS-2 T71-D模拟的干涉形变场存在部分误差，可能是反演模型导致的误差残留，由于残差相对于整体地震造成的同震形变场来说相对较小，所以其对模型结果的影响可以忽略。ALOS-2 T180-A模拟的干涉形变场中，主要的残差为远场的形变场误差，是由于原始SAR形变场中的远场误差所导致，因为远场误差相对较小，所以不会对模型结果造成很大影响。通过分析观测值和模拟值的残差值可以看出，震中区域反演的最大残差小于5 cm。总体而言，InSAR模拟干涉图能够对干涉形变场进行很好的解释，证明本文确定的滑动区间是稳定可靠的。

Murat^[31]基于地震波进行反演，确定地震发震断层为倾斜的，与阿拉伯板块沿着扎格罗斯冲断带（Zagros thrust zone，ZTZ）斜向控制的源区地震构造一致，具有逆冲和右旋分量，地震主要影响ZTZ断裂结构的MFF和HZF，MFF是一个基底逆冲断层（北东向），HZF是一个具有小的右旋分量（北东向）的逆断层。通过同震滑动分布模型得到基于地震波的反演结果是破裂发生在20 km深兼具最大滑动量为9 m的滑动，主要破裂长度为50 km，宽度为15 km，其中有2/3向南（伊朗方向）传播，因此地震发生后，伊朗地区首先受到地震的破坏。该结果与本文用InSAR分析和反演的结果具有高度的一致性。

InSAR获取的是LOS向形变场，该形变场是地表东（E）、北（N）、垂直向（U）的形变分量在雷达脉冲入射方向投影的叠加结果。常规地震三维形变场的获取需要用不同卫星的升降轨数据或结合GPS等大地测量数据、地震破裂模型等^[32]对地表E、N、U方向的形变分量进行参数解算，但如果卫星参数比较接近或者没有GPS以及足够的SAR数据时，则无法进行精确的三维解算。

本文使用由InSAR数据反演的断层滑动分布重新计算地表位移场，同时得到三维形变场，该方法获得的三维形变场去除了失相干地区的影响，同时考虑了断层的滑动分布对形变结果的影响，使得结果更加合理、精确。三维形变分量的结果如图 6所示，最终获取的三维形变场结果如图 7所示。

图 6中，从东西向分量（图 6（a））的结果中可以看出，主要存在3个形变中心，西南和东北形变中心均为向东运动的形变特征，最大东向形变量为32.9 cm；东南形变中心表现为向西运动形变特征，最大西向形变量为16.5 cm；远场表现为接近于0的东西向运动。南北向分量（图 6（b））的结果表明，形变区域主要表现为向南水平运动，最大形变量约53.3 cm，远场主要表现为接近于0的运动。图 6（c）为垂直向形变场，即为图 7中的形变场底图。

图 7中，箭头表示东西向和南北向矢量的叠加，即水平形变场，定义北方向和东方向为正。垂直向的形变分量结果表明，地震导致左右两盘形变区域，左盘表现为上升形变特征，得到的垂直向最大隆升形变约为78.8 cm；右盘表现为沉降形变特征，垂直向的最大沉降量约31.4 cm；左盘隆升形变量远大于右盘沉降形变量。通过结合东西向和南北向形变场与垂直向形变场的结果，可以发现，发震断层面为自西北向东南延伸的倾斜面，发震断层为倾斜的逆冲断裂，兼具右旋特性。

He等^[5]使用SAR数据通过差分干涉处理、像素偏移跟踪、多孔径InSAR和脉冲重叠干涉（burst-overlap interferometry，BOI）方法分析了本次地震的三维形变场，并得到南北向的形变最为显著，形变量达76 cm，且InSAR结合BOI获取的三维形变场精度更高，本文获取的形变场结果与其得到的三维形变场吻合度较高，只有东西向形变场存在细微的差别，可能是获取三维形变场模型不一致所导致，但整体形变量和形变场的分布非常一致。

4 结语

2017年11月12日的伊朗萨尔波勒扎哈卜Mw7.3地震是板块碰撞产生倾斜逆冲断层作用的结果，地震造成地表在断裂带两侧产生了相对抬升和沉降现象，LOS向的最大位移量约为88.8 cm。基于ALOS-2和Sentinel-1两种卫星传感器获取的地震同震形变场进行震源机制反演，使用二步反演算法进行反演，反演结果表明，两伊地震的发震断层走向为347°，倾角为15°，断层倾滑以走滑兼逆冲为主，最大倾滑量达到4.5 m，其位于深度16~19 km处。累积释放地震矩约合矩震级为Mw7.36。InSAR震源机制反演结果与USGS结果有很好的吻合度。使用ALOS-2和Sentinel-1数据获取了同震形变场，结合反演结果获得了本次地震同震三维形变场。通过三维形变场结果可以看出，主要形变是垂直向的抬升，最大抬升量达78.8 cm，东西向和南北向的形变合成矢量后可以得出，发震断层为具有右旋性质的北北西向逆冲断层，右旋走滑相对较小，但也对形变场的形成有一定贡献，整体分析结果表明，本次地震的破裂面为东北倾斜的平面。通过对地震的震源机制进行反演分析，可以更好地了解和描述本次地震的发震断层，为进一步研究逆冲断裂的几何特征等奠定基础。