2021, Vol. 64
2. 中山大学测绘科学与技术学院, 广东 珠海 519082
2. School of Geospatial Engineering and Science, Sun Yat-Sen University, Zhuhai Guangdong 519082, China
2017年4月5日6时9分(协调世界时，UTC)，伊朗东北部桑塞菲德(Sang sefid)地区发生了MW6.1地震，震中位于35.776°N 60.436°E，震源深度13 km(https://earthquak-e.usgs.gov/earthquakes).据当地媒体报道，此次地震造成两人死亡，数十人受伤，靠近震中西北部的是伊朗第二大人口城市马什哈德(Mashhad).此次地震附近没有发现有明显位移的已知断层，分析构造背景认为此次地震属于一次板块间俯冲带地块的挤压造成的盲冲事件(Aflaki et al., 2019)，俯冲带的挤压对Binalud断层带的东南部和Kopeh Dagh断层带的东部(Hollingsworth et al., 2008)造成了一定的冲击；从大范围的构造机制来看，地震的发震原因可以归结为：阿拉伯板块相对于欧亚板块以31 mm·a-1北移(DeMets et al., 2010)，与南移的欧亚板块挤压碰撞造成板块间分布式缩短与剪切作用，同时阿拉伯板块的北移导致伊朗东部的Lut地块随之北向推进，进而与桑塞菲德地区坚硬岩石层挤压碰撞而造成桑塞菲德内部剪切应力作用发震(Jackson，1992).
伊朗桑塞菲德地区被一系列WNW-ESE走向逆冲走滑断层包围(图 1)，且由于阿拉伯板块的北移俯冲运动，造成了该区域存在强烈的地壳运动、构造变形与地震事件.其中桑塞菲德北部是伊朗东北部的主要构造变形区，主要活动断裂以一条WNW-ESE走向线型山脉Kopeh Dagh为主，这条由阿拉伯板块北向运动孕育的长约600 km的逆冲系统隔绝了伊朗与土库曼斯坦地区；桑塞菲德西部为著名的Binalud断层带和Mashhad断层带与伊朗东北部的一系列WNW-ESE走向逆冲走滑断层共同吸收了绝大部分的阿拉伯板块汇聚速率，这条NW走向范围约130 km的Binalud山脉是伊朗地区构造变形最严重的断层带之一(Shabanian et al., 2012)，它也是稳定的伊朗中部板块与活跃的东北部的边界，形成了伊朗东北部活跃的构造带(Shabanian et al., 2009, 2010)；南部是以Lut地块为主的N-S走向的剪切带(Foroutan et al., 2014)，它在阿拉伯板块的推动下逐渐北移与桑塞菲德地块挤压扭转造成Lut地块东西侧逆向运动，形成了Lut地块东西两侧N-S走向的右旋走滑断层以及北侧Khaf和Jangal逆冲系统，这也被认为是伊朗-阿富汗东部边界地震带的主要构造(Walker and Jackson, 2004).然而由于桑塞菲德地区的沙漠性气候且地处偏远，又未发生过大型事件导致了该地区构造资料的不足，也阻碍了学者们对该地区盲断层的识别.
|
图 1 伊朗桑塞菲德地震构造背景 (a) 伊朗地区相对于欧亚大陆参考框架二维GPS速度场; (b) 本文研究区域. 其中震源球为不同机构和学者提供的震源机制解，深黄色点为IRSC(Iranian Seismological Center)记录的桑塞菲德主震后两个月内的余震(http://irsc.ut.ac.ir)，红色线条来自USGS记录的伊朗已知断层(https://catalog.data.gov/dataset/major-faults-in-iran-flt2cg)，KDF为Kopeh Dagh断层带，MF为Mashhad断层带，BF为Binalud断层带，KF为Khaf断层带，JF为Jangal断层带，深蓝色为伊朗地区2006年至2015年相对于欧亚大陆参考框架二维GPS速度场(Khorrami et al., 2019). Fig. 1 Tectonic background of Sang sefid earthquake in IRAN (a) GPS velocities relative to a Eurasian reference frame in Iran; (b) The Area Of Interest of this study. The focal mechanisms indicate different results from different institutions and scholars, The dark yellow dots represent the aftershocks two months after the mainshock recorded by IRSC (Iranian Seismological Center; http://irsc.ut.ac.ir), The red lines represent faults from the USGS World Energy Project (https://catalog.data.gov/dataset/major-faults-in-iran-flt2cg), Abbreviations for faults are KDF: Kopeh Dagh fault, MF: Mashhad fault, BF: Binalud fault, KF: Khaf fault, JF: Jangal fault, The dark blue vectors are the GPS velocity field relative to a Eurasian reference frame for the period 2006-2015 in Iran. |
伊朗桑塞菲德地震发生后，一些国内外学者和专业机构对此次地震的构造背景和发震断层的倾向进行了研究(表 1)对于发震断层的倾向目前仍存在争议.Ghayournajarkar等(2020)使用哨兵一号和ALOS-2的InSAR数据获取了精确的同震形变场，并利用PSO(Particle Swarm Optimization)非线性优化反演方法推断该地震的盲断层为向东北倾斜的右旋逆冲断层.Xu等(2018)使用哨兵一号数据用于同震反演，通过分析恢复的二维形变场并研究了主震库伦应力的分布，推断该地震的发震断层为向西南倾斜的逆冲断层.Su等(2019)基于均匀断层模型采用拉普拉斯平滑约束的方式反演推断断层滑动分布主要为右旋剪切位移，进而判断该地震的发震断层为东北倾向的逆冲兼具右旋走滑断层.Aflaki等(2019)综合分析余震分布与震区现场地质特征考察结果，推断该地震发震断层为西北-东南走向且向东北倾斜的逆冲断层.综上，不同学者对于桑塞菲德地震的断层构造机制的认识充满了争议，其中的主要原因在于基本都是使用传统的二步法反演断层滑动分布特征，这也导致获取的滑动分布特征的单一，无法分析不同反演方法之间的各向异性结果，以得到更全面的印证.因此，本文采用相同的断层模型和地表观测数据作为约束，运用不同的反演方法更辩证地分析滑动分布的合理性，进而为该盲断层的构造机制提供更为全面的认识.
|
表 1 本文使用的Sentinel-1A数据参数 Table 1 Parameters of Sentinel-1A data used in this study |
本文首先根据震源机制解(USGS，GCMT，IRSC)，结合获取的InSAR升降轨同震形变场的空间变化形态及余震分布，假设了发震断层东北倾和南倾的两种倾向可能，并基于弹性半空间构建了单一平面断层模型，使用二维InSAR形变数据为约束，采用基于约束最小二乘的动力学反演方法(Wang et al., 2008)以及贝叶斯反演方法(Amey et al., 2018)结合解算得到潜在发震断层非线性断层参数和断层面上的滑动分布特征.最后，本文具体讨论了反演得到的断层参数和滑动分布特征之间的各向异性，并结合余震空间分布特征以及地质构造背景进行综合分析，对桑塞菲德盲冲断层的构造机制进行了解释.
1 InSAR数据及盲断层倾向确定 1.1 InSAR数据本研究采用两对Sentinel-1A宽幅数据和两对ALOS-2影像分别获取了桑塞菲德地震的同震形变场.其中，选取了哨兵一号T13升轨和T93降轨，ALOS-2 T171升轨和T64降轨数据.图 1为伊朗东北部构造背景以及桑塞菲德主震位置及余震分布，表 1为本文使用SAR影像详细参数.
首先利用精密轨道以及SRTM 1 sec数字高程模型采用几何配准的方式进行哨兵数据参考影像的配准，然后在进行差分干涉处理之前对配准好的堆栈进行增强谱分集(Enhanced Spectral Diversity，ESD)处理，以校正子带中每个burst重叠部分的相位不连续(Ma et al., 2019)，从而使得burst重叠区域的相位差小于3°；在进行精配准之后，利用分辨率为30 m的数字地形高程模型(SRTM 1 sec)进行去地形处理，然后采用传统的Goldstein滤波器(Goldstein and Werner, 1998)对干涉图进行滤波处理，以抑制相位噪声，提高干涉图信噪比，使干涉条纹更光滑，从而减小相位解缠的难度.最后使用SNAPHU方法(Chen and Zebker, 2000)对干涉相位进行二维解缠，ALOS-2数据采用的L波段受电离层与轨道误差影响较大，在解缠后采用二阶曲面去除电离层与轨道误差代表的长波信号，并将干涉相位地理编码为WGS84地理坐标，最终得到清晰的同震干涉条纹图和同震形变图(图 2-3).
|
图 2 基于InSAR数据获取的同震形变干涉图 (a)，(b) 哨兵一号的升(左)降(右)轨干涉图；(c)，(d) ALOS-2的升(左)降(右)轨干涉图.其中黑色和红色分别代表东北倾断层和南倾断层地表迹线及震源机制解. Fig. 2 The wrapped interferogram of the co-seismic deformation based on InSAR data (a), (b) The wrapped interferogram of Sentinel-1A ascending and descending track; (c), (d) The wrapped interferogram of ALOS-2 ascending and descending track. The black and red represent the fault plane and focal mechanisms of NE-dipping and South-dipping fault respectively. |
|
图 3 地震同震形变场LOS向解缠干涉位移 (a)，(b) 哨兵一号的升(左)降(右)轨; (c)，(d) ALOS-2的升(左)降(右)轨. 其中黑色和红色分别代表东北倾断层和南倾断层地表迹线及震源机制解. Fig. 3 Unwrapped Interferometric displacement in LOS direction of earthquake event (a), (b) The unwrapped interferometric displacement of Sentinel-1A ascending and descending track; (c), (d) The unwrapped interferometric displacement of ALOS-2 ascending and descending track. The black and red represent the fault plane and focal mechanisms of NE-dipping and South-dipping fault respectively. |
图 2为进行D-InSAR处理后的Sentinel-1A卫星和ALOS-2 SAR传感器获取的升轨和降轨同震形变干涉图，它清晰地反映出了桑塞菲德地震引起的地表形变，整个升降轨干涉条纹图空间形态特征明显，大致呈西北-东南走向椭球牛眼形分布，同震干涉图分布均匀，覆盖范围达到30 km×20 km.干涉条纹沿包裹中心向外扩散，较为稀疏，形变梯度较小，但也足够表明地震引起了明显的地表位移；其次，升降轨干涉图都为单叶干涉条纹，说明桑塞菲德地震引起的地表形变以沿视线向单一隆升或远离，表明此次地震发震断层以逆冲为主的断层性质.
1.2 盲断层模型确定当地震构造背景已知或者发震断层确定时，利用大地测量观测数据(InSAR，GNSS)进行单断层或多断层参数反演的研究已经趋于成熟(Kuang et al., 2019；Feng et al., 2020；屈春燕等，2017)，但当遇到中型地震(6≤MW≤6.5)且地表破裂情况未知，发震构造背景模糊时，确定发震断层参数尤其是倾向一直是一个难点，这也阻碍了盲断裂地带的地震危险性分析(Biggs et al., 2006；Materna et al., 2019).针对如何在多数地震情况未知时确定盲断层的倾向和位置这个问题，本文提出了相应的解决方案.众所周知，InSAR数据具有很高的空间分辨率，这也为确定断层的位置提供了十分宝贵的信息，通过InSAR观测数据和重新定位的余震信息可以假设出潜在的盲断层位置，然后由InSAR观测数据反演约束出最佳拟合断层面，进而获得最佳断层模型.
表 2为国内外学者和科研机构反演获得的桑塞菲德地震的震源机制解.其中，USGS，GCMT和IRSC认为该地震的发震断层具有两种潜在的可能断层，主要为东北倾和南倾.不同学者对该发震断层的倾向也有不同的认识，Xu等(2018)认为该盲断层向西南倾的可能性更大，断层为西北-东南走向，方位角和倾角分别为120°和40°；具有争议的是其他大多数学者都认为该盲断层为东北倾，且反演结果的滑动角差异较小，均反映了该断层逆冲为主且具有右旋分量的特性；但不同学者给出的方位角和倾角均有一定的差异，方位角最大差9.2°，倾角最大差9.9°，这也造成了该发震断层具体倾向的模糊.综合分析不同学者给出的结果，可以认为桑塞菲德地震的发震断层主要有东北倾和南倾两种潜在可能，关于本文的结果将做进一步分析和解释.
|
|
表 2 不同学者和机构反演的断层参数 Table 2 Fault parameters retrieved by different scholars and institutions |
通常，发震断层的地表迹线投影位于干涉条纹的边缘附近.因此，根据桑塞菲德地震获取的同震形变干涉图，本文假设了WNW-ESE走向和W-E走向两种潜在断层(图 2).同时，根据USGS和IRSC给出的震源机制解，设置了断层初始走向、倾角和滑动角搜索范围，并假设地球是一个均匀的弹性半空间，泊松比为0.25，基于Okada85弹性位错模型(Okada，1985)构建了东北倾和南倾两种单一平面断层模型，断层模型深度为15 km，面积为30 km×20 km，以覆盖整个滑动区域.
2 二维形变场分解与滑动分布反演 2.1 二维形变场分解通过目视判别桑塞菲德地震同震形变场(图 3)，可以观察到升降轨的同震形场在空间变化形态具有一定的相似性.从分布范围来看，升降轨同震形变场都呈西北-东南向分布，这与该区域北部一系列剪切断层(Binalud断层和Kopeh Dagh断层)平行，且形变场的分布区域也大致相似；从分布特征来看，升降轨同震形变场都呈椭球牛眼形状，且都展现了沿卫星视线方向(line of sight，LOS)方向隆升位移的特征；其中，升轨和降轨形变图各自展现出沿视线向靠近卫星的最大位移量达10.2 cm和13.4 cm；对于视线向位移，降轨的东向位移会对视线向位移产生正向作用，升轨的东向位移会对视线向位移产生负作用(Wang et al., 2017)，而降轨的视线向隆升要稍大于升轨，因此可知桑塞菲德地震引起的地表形变以东向为主.从地震活动物理机制来看，该地震升降轨同震形变场沿视线方向一致性隆起的特征，也更加印证了此次地震发震断层的逆冲特性.
因为升降轨视线向位移不具备极性特征，仅使用视线向位移无法确定桑塞菲德地震的发震断层，本文采用结合应力应变模型的方差分量估计方法(Liu et al., 2018)联合哨兵和ALOS-2数据恢复二维形变.根据卫星的成像模式可知，雷达视线向位移由东西向、南北向和垂直向位移构成，而由于哨兵一号和ALOS-2的近极轨道飞行(与经线差10°左右)，导致雷达观测对南北向位移极不敏感且计算矩阵奇异，因此为保证东西向和垂直向的分解精度，将南北向位移做忽略不计处理，此时将雷达视线分解如式(1)：
|
(1) |
式中DLOS为视线向位移，UEU为分解矩阵，λ为雷达入射角，θ为卫星升降轨卫星方位角，vE和vU为东西和垂直位移.通过以上的关系式，并结合应力应变模型(袁霜等，2020；He et al., 2019)即可解算出东西和垂直位移.由于桑塞菲德地震区大多以荒漠岩石层为主，相干性较高，因此在解算过程中，为提高解算精度，提取相干性高于0.4的像素点进行解算，得到二维形变场(图 4).
|
图 4 二维形变场 (a) 东西方向形变; (b) 垂直方向形变, 其中黑色和红色分别代表东北倾断层和南倾断层地表迹线及震源机制解. Fig. 4 2-dimensional displacement field (a) East and West displacement; (b) Vertical displacement. The black and red represent the fault plane and focal mechanisms of NE-dipping and South-dipping fault respectively. |
分析二维形变场可以发现，桑塞菲德地震引起了明显的东向位移和隆升位移，其中隆升位移主要位于震源中心附近(图 4b)，最大值达到13.4 cm，并且在隆升区域西南部和东北部展现出轻微的下沉；东向位移主要位于形变场的东北部(图 4a)，最大值达到6.8 cm，在形变场的西南部也显示了部分西向位移，最大值达到4.1 cm；其次，分析本文假设的两种潜在断层，水平形变场主导的东向位移即为两种潜在断层上盘的东向运动，可解释为东北倾盲断层的右旋运动或者南倾盲断层的左旋运动；而水平形变场西南部的西向运动量与范围相对于东北部的东向运动较小，分析局部构造认为这可能是由于震区南部Lut地块向北推挤造成的块体挤压缩短而导致的东西形变场的极性特征(图 1)，这与假设东北倾盲断层为发震断层时的右旋运动并不矛盾.两种潜在断层都可以很好的解释二维形变场，因此，仅从二维形变场仍然无法确定盲断层的倾向，还需进行反演分析.
2.2 滑动分布反演 2.2.1 约束最小二乘动力学反演本文采用解算升降轨一维形变获取的二维形变场数据为约束，并基于弹性半空间构建的东北倾和南倾两种矩形位错模型，沿长度和宽度将断层划分成30×20个矩形块，矩形位错大小设置为1 km×1 km，使用SDM(Steepest Descent Method)反演获得了与地表测量数据具有最佳吻合度的非线性断层参数和发震断层面上的最优滑动分布特征.InSAR的高空间分辨率数据虽然有利于描述断层细节，但也加大了反演难度.因此，在反演之前，为提高计算效率，减小因数据冗杂而产生的收敛问题，对获取的二维形变场数据进行四叉树降采样(Jónsson et al., 2002)处理，采用GBIS软件(Bagnardi and Hooper, 2018)提供的降采样包，设置形变梯度阈值为0.005 m，并根据相干性加权以提高反演精度，充分利用质量较好的数据点，最终降采样得到2645个数据点参与反演.
经过近10000次迭代达到收敛，得到两种潜在断层的非线性断层参数(表 3)和滑动分布特征(图 5a-b).分析滑动分布特性可以发现，东北倾和南倾断层破裂都未到达地表，这与野外地质考察结果一致.东北倾断层滑动较为集中，沿走向6~22 km处，主要盲冲作用在7 km深度(图 5a)，最大滑动量达到1.09 m；且反演得到的平均滑动角为128.80°(表 3)，显示该断层是以逆冲为主的兼具右旋走滑的逆断层，这与二维形变场分析的结果保持一致.南倾断层盲冲作用主要分布在3~11 km深度，沿走向8~20 km处，最大滑动量出现在7 km深处达到0.69 m(图 5b)；平均滑动角为77.50°(表 3)，表明该断层是以逆冲运动为主的逆断层，这也与二维形变场分析的结果一致.另外，反演得到震源参数与各机构(USGS，GCMT，IRSC)基本一致，其中各几何参数的误差主要是因为各机构给出的矩震级是利用主震发生时地震波反演而得，它描述的是主震发生时的瞬时强度，而桑塞菲德地震发生后的很短时间内在主震附近发生了至少三次MW5.0以上的余震(Aflaki et al., 2019)，这会导致主震的发震断层进一步破裂，而InSAR数据又是在震前和震后获取的静态数据，其中ALOS-2的升轨数据的时间基线达到154天，这必然也会涵盖了由于余震导致的震后位移以及震前的部分非构造信号；其次，ALOS-2数据在研究区右上角的不足也会导致InSAR空间分辨率的降低而造成一定反演结果的差异.
|
|
表 3 SDM反演获取的断层参数 Table 3 Fault parameters inverted by SDM |
|
图 5 SDM反演和贝叶斯反演的滑动分布 (a) SDM反演的东北倾断层滑动分布; (b) SDM反演的南倾断层滑动分布; (c) 贝叶斯反演的东北倾断层滑动分布; (d) 贝叶斯反演的南倾断层滑动分布. 其中蓝色实心点为USGS-NEIC记录的主震后三个月内(2017年4月5日至2017年7月5日)MW≥3.5的余震数据，黑色五角星为主震震中. Fig. 5 Slip distribution inverted by SDM and a Bayesian inversion method (a) Slip distribution of NE-dipping fault inverted by SDM; (b) Slip distribution of South-dipping fault inverted by SDM; (c) Slip distribution of NE-dipping fault inverted by Bayesian inversion method; (d) Slip distribution of South-dipping fault inverted by Bayesian inversion method. The blue dots represent the aftershocks three months (5 April 2017 to 5 July 2017) after the mainshock recorded by USGS-NEIC, the black star represents the epicenter of mainshock. |
图 6和图 7为反演获取的二维形变场的模拟值和残差值，从残差图可以看出东北倾断层和南倾断层模型与二维同震形变场都有着较好的适配，两种潜在断层模型的拟合残差都较小，其中，东北倾和南倾断层模型反演得到的二维形变值和模拟值的平均RMS(Root Mean Square)残差分别为东西向0.42 cm和0.50 cm，垂直向0.31 cm和0.49 cm，观测数据与模型的相关系数分别为0.997和0.992.可见SDM反演的结果更支持东北倾断层为可能的发震断层，但由于两种断层模型拟合性都较好，本文将基于此断层模型进行贝叶斯反演滑动分布，获取最优滑动分布特性的后验概率密度，并做进一步验证分析.
|
图 6 反演获得的东北倾断层模型的二维形变观测值、模拟值和残差 Fig. 6 2D observed, modeled and residual displacement for NE-dipping fault |
|
图 7 反演获得的南倾断层模型的二维形变观测值、模拟值和残差 Fig. 7 2D observed, modeled and residual displacement for South-dipping fault |
利用贝叶斯方法反演大地测量数据集可以快速地获取源断层模型参数的后验概率密度，进而获取断层模型上具有最大后验概率的滑动分布特征(Amey，2018).根据贝叶斯定理，利用似然函数和先验概率构建后验概率函数如式(2)：
|
(2) |
式中，m为包含滑动量和滑动角的断层模型参数，d为大地测量数据集(InSAR数据)，p(m|d)为断层模型参数的后验概率，p(d|m)是描述大地测量数据与断层模型参数匹配度的似然值，p(m)为断层模型参数的先验概率，p(d)为不考虑模型参数的标准化常量.反演之前，通过设置断层模型参数的可能范围来计算模型参数的先验概率，其中东北倾断层的滑动角的搜索范围设置在90°~180°之间，南倾断层的滑动角的搜索范围设置在0°~90°之间，然后使用蒙特卡罗方法进行先验样本采样，并通过对使用模型参数预测的观测值与实际观测值的拟合度加权计算似然函数如式(3)：
|
(3) |
式中，N为InSAR数据点的总个数，d为InSAR观测值，Σd为InSAR数据的方差-协方差，它描述的是InSAR观测值之间的相关性，用来对InSAR数据加权，s为滑动量，G为基于弹性半空间内由倾角和滑动角计算的设计矩阵(Bagnardi and Hooper, 2018)，模拟形变值由G与s计算而得.反演通过迭代的方式进行，在每次迭代后设计矩阵G会随着重新计算的滑动角而更新，进而会不断更新似然函数以计算断层参数的后验概率，以获取与观测值拟合度最高的断层滑动参数.
采用SDM反演得到的具有最佳拟合的断层模型作为贝叶斯反演的先验断层模型，通过迭代反演获取了具有最大后验概率的断层参数(表 4)和滑动分布特征(图 5c-d).其中，东北倾和南倾断层的平均滑动角分别为129.07°和58.94°，最大滑动量分别为1.02 m和0.66 m(表 4).滑动分布结果表明东北倾断层逆冲兼具右旋走滑以及南倾断层逆冲兼具左旋走滑的运动特性，东北倾断层的滑动主要集中在5~9 km深度范围，最大滑动量出现在7 km深处(图 5c)，与SDM反演结果基本一致；而南倾断层的滑动集中在3~9 km深度范围，最大滑动量分布在6 km深处(图 5d)，与SDM反演结果略有差异.其次，对比分析USGS-NEIC提供的震后的数次余震与滑动分布(图 5)，余震主要分布在东北倾断层主震滑动分布范围四周(图 5a，c)，与主震的滑动分布形成良好的互补，而在南倾断层上却显示出与主震滑动范围更多的重合(图 5b，d)，这与地震学余震机理不符，因此可以认为USGS给出的余震结果与东北倾断层更符合.综上所述，由于采用相同的断层模型与大地测量观测数据，贝叶斯反演的结果与SDM反演具有一定的相似性，但仍然可以看到差异性，关于差异以及断层倾向本文将做进一步分析和讨论.
|
|
表 4 贝叶斯反演获取的断层参数 Table 4 Fault parameters inverted by a Bayesian inversion method |
为了更好地讨论桑塞菲德地震发震断层的倾向，本文将约束最小二乘动力学反演方法与贝叶斯反演方法获取的结果进行对比与讨论.两种反演方法得到的东北倾断层参数基本一致(表 3-4)，其中滑动角仅差0.27°，滑动量相差0.07 m；而南倾断层参数的差异较大，滑动角相差18.56°，滑动量相差0.03 m.其次，两种反演方法获得的东北倾断层的滑动分布特征也具有较高的一致性(图 5a，c)，滑动都集中在7 km深度上下，且滑动角的方向都显示出了发震断层逆冲兼具右旋走滑的运动性质.对于南倾断层，利用SDM方法反演的滑动分布范围更大(图 5b)，较图 5d破裂范围要大2 km，且展现的几乎是纯逆冲的运动性质，而贝叶斯反演滑动分布结果却显示出了更多的左旋分量(图 5d)；此外在11 km至13 km处显示出部分余滑，与SDM结果有一定差异.利用地震学方法得出的震源机制解是研究地震性质，剖析断层模型的主要手段，大地测量数据与断层模型的拟合度也是作为确定盲断层模型决策的重要信息，反演残差结果也显示南倾断层反演获得东西和垂直向位移的残差0.50 cm和0.49 cm都要大于东北倾断层的残差0.42 cm和0.31 cm，并且从残差图也可以清晰的看到南倾断层在主要形变区较东北倾断层显示出了更多的系统性残差(图 6-7)，因此认为东北倾断层为此次地震的发震断层.
利用升降轨联合解算的桑塞菲德地震二维形变场可以观察到更多地表位移信息，东西和垂直形变场都显示出WNW-ESE走向的形变区域(图 4a)，并且东西向形变场中主要东向位移和部分西向位移也表明了该发震断层走滑特征，这与两种反演方法获取的东北倾断层的滑动分布特征一致，而SDM反演的南倾断层却显示出几乎纯逆冲的运动特性(图 5b)，与二维形变场特征不符；其次，东北倾向断层的WNW-ESE走向与伊朗东北部断层带已知的逆冲活动断层(Binalud断层带，Kopeh Dagh断层带)走向基本一致(图 1).
通过构建桑塞菲德主震后进行重新定位的59次余震分布(Aflaki et al., 2019)与两种潜在断层相对位置模型(图 8)，可以看到余震在反演获得的东北倾断层模型震源中心附近分布较为密集，且余震的震源中心有沿着WNW-ESE走向发展的趋势(图 8)，这与东北倾断层的走向基本一致.从两种潜在断层与余震的相对深度位置来看(图 8)，余震的深度有随着WNW-ESE走向断层倾向而逐渐增加的趋势，且余震的深度分布与东北倾向断层模型更吻合.另外，发现重新定位的余震深度都系统性地大于本文反演的主震深度，这可能是由于该地区地壳介质不均匀，岩石新生成且岩体破碎导致反演算法假设的均匀介质的地壳速度要大于实际地壳速度，从而使反演的主震深度偏高.
|
图 8 桑塞菲德地震余震分布 黑色表示东北倾断层，红色表示南倾断层. Fig. 8 The aftershocks of Sang sefid earthquake Black and red planes represent NE-dipping and South-dipping respectively. |
综合以上原因，本文认为伊朗桑塞菲德地震的发震断层为WNW-ESE走向，向东北倾斜的逆冲兼具右旋走滑分量的盲断层.同时也能说明对于构造背景模糊又缺乏大地震事件记录，现场考察地质特征不明显的盲断层事件，仅利用高空间分辨率的InSAR观测数据来判断盲断层的倾向仍然非常具有挑战性，未来随着InSAR技术观测精度的进一步提高，或者更丰富的观测数据如GNSS数据参与进行联合约束，以及分布密度更高的余震数据辅助判断，盲断层倾向的判断的可靠性将进一步提升.
3.2 盲断层地震构造意义伊朗桑塞菲德地震发生在伊朗东北部，震区周围分布着已知的活跃走滑断层和逆冲断层(图 1a)，但震区附近缺乏活动断层以及重大地震事件的记录造成了该震区构造背景的模糊.众所周知，断层位置的确定是进行地震危险性分析的首要任务，因此桑塞菲德地震的发生以及该盲断层的识别对于补充该地区地震构造背景具有重要意义.
此次地震发生后短时间内发生了数次MW5.0以上的余震，也充分表明了该盲断层的活跃性，并且本文利用InSAR观测数据反演的滑动分布结果也显示该断层在深度7 km处出现了大于1 m的明显滑动，但都未使该断层破裂至地表，这可能是由于桑塞菲德地区存在覆盖层较厚的沉积物以及坚硬的岩石层(Shabanian et al., 2010)，这也可能导致本文反演获得的震源深度存在高估以及地震震级存在低估，进而也解释了图 8重新定位的余震深度分布(10~20 km)要比反演获得的震源深度(7 km)更深的结果.其次余震整体空间分布呈广泛分散形态，符合盲断层活动断裂带余震分布的特征(Stein and Ekström，1992)，因此该地区可能普遍存在盲断层或次生断层而形成活动盲断裂带，这样一个未被识别的断裂带也对该地区的地震风险造成了隐患，将来应对构造背景模糊的桑塞菲德地区做进一步关注与研究，以丰富该地区地震活动机制并掌握此活动断裂的更多细节.
3.3 桑塞菲德地震破裂模式伊朗桑塞菲德地震被认为是阿拉伯板块北移与欧亚大陆板块碰撞挤压而造成的逆冲事件(Su et al., 2019)，局部来看认为是伊朗东部Lut地块北移与桑塞菲德地区坚硬的岩石层碰撞而孕育的剪切斜冲事件，因此该断层破裂模式可以认为是阿拉伯板块与欧亚大陆板块间俯冲带逆冲断层破裂(Hyndman and Wang, 1993).板块间俯冲带逆冲事件通常有三个特征：同震时的断层上盘前缘隆升、后端下沉和下盘后端下沉，这是由于地震破裂期间逆断层上下盘相向运动，上盘前缘受到下盘的俯冲推覆而缩短隆起，后端受到拉张而下沉，下盘后缘受到上盘的仰冲挤压而下沉，这种破裂模式造成的地表形变与本文获取的二维形变场中垂直位移空间分布形态一致，主要表现为震区的主要隆升与西南部和东北部的轻微下沉.
震前的大地测量学以及地质学研究表明伊朗东北部的构造变形主要集中在块体边界处，包括Lut地块边界以及Kopeh Dagh构造带附近(Mousavi et al., 2013)，GPS速度场的结果也表明Lut地块较欧亚大陆北移的速率已达到5.7 mm·a-1，这一北向作用几乎都被伊朗东北部的一系列WNW-ESE的走滑逆冲系统所吸收(Walpersdorf et al., 2014)，这也是Kopeh Dagh与Binalud构造带西北走向褶皱地貌与逆冲构造的主要发育成因(Lyberis and Manby, 1999)，与此相似的是，桑塞菲德地区也以西北走向褶皱地貌为主且逐渐展现出背斜上升的趋势(Aflaki et al., 2019)，此外桑塞菲德地区地处Lut地块北部，受其北移俯冲的直接推挤碰撞作用效果，因此可以认为桑塞菲德地区也正逐渐形成与Kopeh Dagh与Binalud构造带相似的构造带，即WNW-ESE方向走滑逆冲盲断层为主的构造系统，这也进一步印证了本论文认为的桑塞菲德地震WNW-ESE盲断层为发震断层的判断.由于发震断层未破裂至地表，本文只能通过地表地貌特征以及构造形变来分析桑塞菲德地震的破裂模式，将来随着该发震断层的构造背景逐渐明晰，桑塞菲德断层的破裂模式也将得到进一步解释.
4 结论本文利用哨兵一号和ALOS-2 SAR影像数据获取了伊朗桑塞菲德MW6.1地震高空间分辨率的同震形变场，并结合四幅干涉图和应力应变模型的方差分量估计方法恢复了二维形变场.通过构建东北倾和南倾两种潜在单一断层模型，基于弹性半空间以二维形变数据为约束利用约束最小二乘与贝叶斯反演算法获取了两种潜在断层的滑动分布特征，并结合余震数据和已知活动断层迹线判断了发震断层的倾向，得到以下几点认识：
(1) 获取的哨兵和ALOS-2升降轨同震干涉条纹都呈WNW-ESE走向椭球牛眼形分布，由于地震未破裂至地表，形变场相干性较好，条纹清晰，且升降轨都为单叶条纹表明此次地震断层主要逆冲的特性.解缠后的同震形变场都以LOS向隆升为主.利用多源数据融合恢复的二维形变场主要显示出震源附近的主要东向位移和隆升位移，以及西南部的少量西向位移，形变主要发生在断层上盘，解释为东北倾断层的逆冲右旋走滑运动以及南倾断层的逆冲左旋走滑运动，且水平方向的西向位移与震区南部Lut地块的北移俯冲推挤构造机理一致，二维形变场的空间分布特征也与俯冲带逆冲断层破裂模式吻合.
(2) 使用约束最小二乘动力学反演和贝叶斯反演的东北倾断层的滑动分布具有较高的一致性，结果显示滑动主要分布在地下5~9 km处并未破裂至地表，在地下7 km处孕育了超过1 m的最大滑动量，平均滑动角约128.8°，矩张量达1.427×1018N·m，相当于震级MW6.07的地震，显示出主要逆冲兼具右旋走滑分量的运动性质；而两种方法反演的南倾断层的滑动分布结果在破裂范围和滑动角具有一定的差异性；东北倾断层WNW-ESE走向也与伊朗东北部已知的活动断层平行(图 1a)，通过将余震在断层深度上的投影也发现余震的分布与东北倾断层滑动分布具有更好的符合，并且余震深度的扩展方向也与东北倾断层的倾向保持一致，反演残差也显示东北倾断层具有更好的拟合结果.综合分析东北倾断层与已知活动断层构造关系，地震学反演获取的震源机制解的重要性，余震分布和滑动分布合理性结果认为东北倾断层为桑塞菲德地震的发震断裂.
(3) 通过分析余震的空间分布形态，震前伊朗东北部地块构造位移以及桑塞菲德地区褶皱背斜隆升的地貌认为该地区可能普遍存在盲断层而形成盲断裂带，桑塞菲德地震的发生也让该活动断裂带受到关注，未来还应继续观察该地区的长期震后松弛效应以获取更细节性的构造特征，为获取该地区置信度更高的地震危险性评估为服务.本文着眼于对盲断层的识别，发现仅通过地表观测数据来识别盲断层依旧充满挑战性，而采用不同反演方法进行比较分析，再通过余震分布和已知活动断裂资料加以辅助有助于对盲断层构造机制的研究，也为进一步对盲断层的识别研究提供了新思路.
致谢 本文所使用的Sentinel-1数据由ESA/Copernicus提供.GFZ德国地学中心汪荣江老师提供了SDM约束最小二乘反演软件.Dr.Amey提供的SlipBERI贝叶斯滑动分布反演工具包.Dr.Bagnardi提供的GBIS贝叶斯反演软件工具包.中山大学地球科学与工程学院周宇教授和日本东北大学理学院地球物理系Dr.Ghayournajarkar对本文提出的宝贵建议.中南大学地球科学与信息物理学院刘计洪博士提出的关于应力应变模型的建议.本文中的全部图形由GMT 6.1.0绘制.在此一并感谢.
Aflaki M, Mousavi Z, Ghods A, et al. 2019. The 2017 MW 6 Sefid Sang earthquake and its implication for the geodynamics of NE Iran. Geophysical Journal International, 218(2): 1227-1245. DOI:10.1093/gji/ggz172 |
Amey R M J. 2018. The fractal nature of fault slip and its incorporation into earthquake slip inversions[Ph. D. thesis]. Leeds: The University of Leeds.
|
Amey R M J, Hooper A, Walters R J. 2018. A Bayesian method for incorporating self-similarity into earthquake slip inversions. Journal of Geophysical Research: Solid Earth, 123(7): 6052-6071. DOI:10.1029/2017JB015316 |
Bagnardi M, Hooper A. 2018. Inversion of surface deformation data for rapid estimates of source parameters and uncertainties: a bayesian approach. Geochemistry, Geophysics, Geosystems, 19(7): 2194-2211. DOI:10.1029/2018gc007585 |
Biggs J, Bergman E, Emmerson B, et al. 2006. Fault identification for buried strike-slip earthquakes using InSAR: The 1994 and 2004 Al Hoceima, Morocco earthquakes. Geophysical Journal International, 166(3): 1347-1362. DOI:10.1111/j.1365-246X.2006.03071.x |
Chen C W, Zebker H A. 2000. Network approaches to two-dimensional phase unwrapping: intractability and two new algorithms. Journal of the Optical Society of America A, 17(3): 401-414. DOI:10.1364/JOSAA.17.000401 |
DeMets C, Gordon R G, Argus D F. 2010. Geologically current plate motions. Geophysical Journal International, 181(1): 1-80. DOI:10.1111/j.1365-246X.2009.04491.x |
Feng W P, Samsonov S, Qiu Q, et al. 2020. Orthogonal fault rupture and rapid Postseismic deformation following 2019 Ridgecrest, California, earthquake sequence revealed from geodetic observations. Geophysical Research Letters, 47(5): e2019GL086888. DOI:10.1029/2019gl086888 |
Foroutan M, Meyer B, Sébrier M, et al. 2014. Late Pleistocene-Holocene right slip rate and paleoseismology of the Nayband fault, western margin of the Lut block, Iran. Journal of Geophysical Research: Solid Earth, 119(4): 3517-3560. DOI:10.1002/2013JB010746 |
Ghayournajarkar N, Fukushima Y. 2020. Determination of the dipping direction of a blind reverse fault from InSAR: case study on the 2017 Sefid Sang earthquake, northeastern Iran. Earth Planets Space, 72: 64. DOI:10.1186/s40623-020-01190-6 |
Goldstein R M, Werner C L. 1998. Radar interferogram filtering for geophysical applications. Geophysical Research Letters, 25(21): 4035-4038. DOI:10.1029/1998gl900033 |
He P, Wen Y M, Xu C J, et al. 2019. Complete three-dimensional near-field surface displacements from imaging geodesy techniques applied to the 2016 Kumamoto earthquake. Remote Sensing of Environment, 232: 111321. DOI:10.1016/j.rse.2019.111321 |
Hollingsworth J, Jackson J, Walker R, et al. 2008. Extrusion tectonics and subduction in the eastern South Caspian region since 10 Ma. Geology, 36(10): 763-766. DOI:10.1130/G25008A.1 |
Hyndman R D, Wang K. 1993. Thermal constraints on the zone of major thrust earthquake failure: The Cascadia Subduction Zone. Journal of Geophysical Research: Solid Earth, 98(B2): 2039-2060. DOI:10.1029/92jb02279 |
Jónsson S, Zebker H, Segall P, et al. 2002. Fault slip distribution of the 1999 MW7.1 Hector mine, California, earthquake, estimated from satellite radar and GPS measurements. Bulletin of the Seismological Society of America, 92(4): 1377-1389. DOI:10.1785/0120000922 |
Jackson J. 1992. Partitioning of strike-slip and convergent motion between Eurasia and Arabia in Eastern Turkey and the Caucasus. Journal of Geophysical Research: Solid Earth, 97(B9): 12471-12479. DOI:10.1029/92JB00944 |
Khorrami F, Vernant P, Masson F, et al. 2019. An up-to-date crustal deformation map of Iran using integrated campaign-mode and permanent GPS velocities. Geophysical Journal International, 217(2): 832-843. DOI:10.1093/gji/ggz045 |
Kuang J M, Ge L L, Metternicht G I, et al. 2019. Coseismic deformation and source model of the 12 November 2017 MW7.3 Kermanshah earthquake (Iran-Iraq border) investigated through DInSAR measurements. International Journal of Remote Sensing, 40(2): 532-554. DOI:10.1080/01431161.2018.1514542 |
Liu J H, Hu J, Li Z W, et al. 2018. A method for measuring 3-D surface deformations with InSAR based on strain model and variance component estimation. IEEE Transactions on Geoscience and Remote Sensing, 56(1): 239-250. DOI:10.1109/TGRS.2017.2745576 |
Lyberis N, Manby G. 1999. Oblique to orthogonal convergence across the Turan block in the post-Miocene. AAPG Bulletin, 83(7): 1135-1160. DOI:10.1306/E4FD2E97-1732-11D7-8645000102C1865D |
Ma Z F, Jiang M, Zhao Y, et al. 2019. Minimum spanning Tree Co-registration approach for time-series sentinel-1 TOPS data. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 12(8): 3004-3013. DOI:10.1109/JSTARS.2019.2920717 |
Materna K, Wei S J, Wang X, et al. 2019. Source characteristics of the 2017 MW6.4 Moijabana, Botswana earthquake, a rare lower-crustal event within an ancient zone of weakness. Earth and Planetary Science Letters, 506: 348-359. DOI:10.1016/j.epsl.2018.11.007 |
Mousavi Z, Walpersdorf A, Walker R T, et al. 2013. Global positioning system constraints on the active tectonics of NE Iran and the South Caspian region. Earth and Planetary Science Letters, 377-378: 287-298. DOI:10.1016/j.epsl.2013.07.007 |
Okada Y. 1985. Surface deformation due to shear and tensile faults in a half-space. Bulletin of the Seismological Society of America, 75(4): 1135-1154. |
Qu C Y, Zuo R H, Shan X J, et al. 2017. Coseismic deformation field of the Nepal MS8.1 earthquake from Sentinel-1A/InSAR data and fault slip inversion. Chinese Journal of Geophysics (in Chinese), 60(1): 151-162. DOI:10.6038/cjg20170113 |
Shabanian E, Bellier O, Siame L, et al. 2009. New tectonic configuration in NE Iran: Active strike-slip faulting between the Kopeh Dagh and Binalud mountains. Tectonics, 28(5): TC5002. DOI:10.1029/2008TC002444 |
Shabanian E, Bellier O, Abbassi M R, et al. 2010. Plio-Quaternary stress states in NE Iran: kopeh dagh and allah dagh-binaludmountain ranges. Tectonophysics, 480(1-4): 280-304. DOI:10.1016/j.tecto.2009.10.022 |
Shabanian E, Bellier O, Siame L, et al. 2012. The Binalud Mountains: A key piece for the geodynamic puzzle of NE Iran. Tectonics, 31(6): TC6003. DOI:10.1029/2012TC003183 |
Stein R S, Ekström G. 1992. Seismicity and geometry of a 110-km-long blind thrust fault 2. Synthesis of the 1982-1985 California Earthquake Sequence. Journal of Geophysical Research: Solid Earth, 97(B4): 4865-4883. DOI:10.1029/91JB02847 |
Su Z, Yang Y H, Li Y S, et al. 2019. Coseismic displacement of the 5 April 2017 Mashhad earthquake (MW6.1) in NEIran through Sentinel-1A TOPS data: New implications for the strainpartitioning in the southern Binalud Mountains. Journal of Asian Earth Sciences, 169: 244-256. DOI:10.1016/j.jseaes.2018.08.010 |
Walker R, Jackson J. 2004. Active tectonics and late Cenozoic strain distribution in central and eastern Iran. Tectonics, 23(5): TC5010. DOI:10.1029/2003TC001529 |
Walpersdorf A, Manighetti I, Mousavi Z, et al. 2014. Present-day kinematics and fault slip rates in eastern Iran, derived from 11 years of GPS data. Journal of Geophysical Research: Solid Earth, 119(2): 1359-1383. DOI:10.1002/2013JB010620 |
Wang H, Liu-Zeng J, Ng A H M, et al. 2017. Sentinel-1 observations of the 2016 Menyuan earthquake: A buried reverse event linked to the left-lateral Haiyuan fault. International Journal of Applied Earth Observation and Geoinformation, 61: 14-21. DOI:10.1016/j.jag.2017.04.011 |
Wang R, Motagh M, Walter T R. 2008. Inversion of slip distribution from coseismic deformation data by a sensitivity-based iterative fitting (SBIF) method. //EGU General Assembly 2008, 10, EGU2008-A-07971.
|
Xu G Y, Xu C J, Wen Y M. 2018. Sentinel-1 observation of the 2017 Sangsefid earthquake, northeastern Iran: rupture of a blind reserve-slip fault near the Eastern Kopeh Dagh. Tectonophysics, 731-732: 131-138. DOI:10.1016/j.tecto.2018.03.009 |
Yuan S, He P, Wen Y M, et al. 2020. Integrated InSAR and strain tensor to estimate three-dimensional coseismic displacements associated with the 2016 MW7.0 Kumamoto Earthquake. ChineseJournal of Geophysics (in Chinese), 63(4): 1340-1356. DOI:10.6038/cjg2020N0308 |
屈春燕, 左荣虎, 单新建, 等. 2017. 尼泊尔MW7.8地震InSAR同震形变场及断层滑动分布. 地球物理学报, 60(1): 151-162. |
袁霜, 何平, 温扬茂, 等. 2020. 综合InSAR和应变张量估计2016年MW7.0熊本地震同震三维形变场. 地球物理学报, 63(4): 1340-1356. DOI:10.6038/cjg2020N0308 |