2019, Vol. 39
2017-11-12(UTC)伊拉克苏莱曼尼亚省东南部哈莱卜杰地区发生MW7.3地震(34.9° N，45.75° E)，震源深度15.5 km(USGS)，地震位于扎格罗斯褶皱冲断带。秦四清等[1]分析该地区大地震孕育过程，并判断其未来地震情况。杨百存等[2]和左荣虎等[3]应用InSAR技术分析该地震，认为其为巴格达地震区第3锁固段向峰值强度点演化过程中发生的一次显著地震，并预测该地区未来有再次发生地震的趋势。
本文收集了覆盖该区域的Sentinel-1A卫星多个轨道SAR数据，分别选择不同轨道震前与震后的SAR影像组成不同轨道的差分干涉对。在所有产生的干涉对中选取干涉效果最优的差分干涉对对本次地震进行成像处理，提取地震形变量，反演出地震在该区域的高质量、连续性较好的D-InSAR同震形变场。在此基础上，应用降采样之后的同震形变场作为地震参数反演约束条件，基于弹性半空间模型(Okada)，采用非线性参数反演方法获取地震发震断层的几何参数，并将得到的结果与USGS及CMT公布的结果进行比较，研究此次地震同震形变场和发震断层参数，对认识该区域地震断层和地质构造移动有一定现实意义，也为地震监测和预防提供了一定的理论支持。
1 地质构造背景扎格罗斯褶皱冲断带(Zagros fold-thrust belt，ZFTB)是扎格罗斯造山带的前陆褶皱冲断带，后者是阿拉伯板块与欧亚板块的碰撞造山带，是新特提斯碰撞造山带的一部分，也是地球上最年轻的陆-陆碰撞造山带，至今构造作用活跃[4]。扎格罗斯褶皱冲断带既是扎格罗斯造山带的组成部分，同时也是波斯湾周缘前陆盆地的楔顶带，其从西北部的土耳其延伸到东南部的霍尔木兹海峡，长达2 000 km，在西北部宽约200 km，到东南部增加到400 km。
扎格罗斯褶皱冲断带形成始于晚白垩纪。由断裂分布指示可知，研究区主要受北东-南西向的挤压应力作用，阿拉伯板块向北俯冲到欧亚板块之下，褶皱断层构造从北东部缝合带向南西方向伸展。扎格罗斯褶皱冲断带在早古生代时以碳酸盐岩和碎屑岩沉积为主，二叠纪至第三纪以碳酸盐岩为主，阿尔卑斯运动使得第三纪和中生代的沉积发生褶皱，形成许多背斜构造[5]。研究区内发育一系列的断裂构造，包括扎格罗斯造山带(从南至北由扎格罗斯褶皱冲断带为板块缝合线，其走向和长度都与褶皱带一致)、萨南达杰-锡尔詹岩浆变质带(SSZ)、乌尔米耶-达克塔尔火山岩浆带(UDMA)和伊朗中部地块等4个构造单元[6]，具体位置如图 1所示。
|
MRF:主近断层; MZRF:主扎格罗斯逆断层; MZPPS-Z:丹帕拉米断层系统; DBF:达什太断层; DF:多罗伊断层; BTZ:比特里斯逆断带; NAFZ:北安纳托利亚断裂带; EAFZ:东安纳托利亚断裂带; NF:纳带断裂区; GF:高尔克断层; SSZ:西斯坦缝合带; KTF:北Tabriz断层 图 1 扎格罗斯造山带主要断裂分布 Fig. 1 Distribution of major faults in the Zagroth orogenic belt |
1957年至今，伊拉克和伊朗边界处共发生6次6级以上地震。其中1957-12-13、1958-08-26和1963-03-24克尔曼沙汗(Kermanshah)MW 6.5、MW6.7和MW6.0地震发生时间间隔很短[7-8]。其余3次为1967-01-11迪亚拉(Diyala)MW6.1地震、2014-08-18伊拉姆(Ilam)MW6.2地震和2017-11-13苏莱曼尼亚(Sulaymaniyah)MW7.3地震。6次地震均发生在扎格罗斯褶皱冲断带。该断层属于逆冲断层，在靠近震中区域，断面走向为西北-东南方向，上盘位于断层西北方向，下盘位于东北方向。此褶皱带由于其特殊的地质构造，使中东地区富含油气资源，也使伊拉克成为一个少震的国家。
2 InSAR数据及处理 2.1 Sentinel-1 SAR数据Sentinel-1由2颗卫星组成，载有C波段合成孔径雷达，不受天气等因素的限制，可提供连续图像[9]。
伊拉克地震前后，Sentinel-1卫星在此区域共有4个轨道的SAR数据，其中降轨数据2个，轨道号为6和79；升轨数据2个，轨道号为72和174，每个轨道分别拍摄三景数据。处理与比较所有数据，选择覆盖地震发生区域的降轨数据为6，时间分别为2017-10-26和2017-11-19；升轨轨道为174，时间分别为2017-10-25和2017-11-18。
对于升轨174，由于数据覆盖地震影响区域的面积较小，故选择同一时间内有部分重叠区域的数据，对数据镶嵌后，再选择合适的覆盖区域进行裁剪，选取的SAR影像覆盖区域如图 1所示。对SAR数据进行简单处理后，选择轨道6和174的SAR数据进行研究，经基线处理后形成干涉对的详细信息如表 1。
|
表 1 Sentinel-1A干涉像对参数 Tab. 1 Parameters of the Sentinel-1A interferogram |
D-InSAR技术的基本思想是基于2幅或多幅雷达复数影像的相位信息，通过干涉的方法提取地面目标三维空间信息和地表微小形变。处理过程中借助90 m分辨率SRTM DEM数据去除差分干涉对的地形效应，经过地理编码获得LOS向同震形变场，结果如图 2。
|
图 2 InSAR地表形变干涉图 Fig. 2 InSAR surface deformation interferogram |
图 2为2个轨道形变干涉纹图，其中降轨干涉纹图覆盖了整个地震形变区域，且干涉纹图以高扎格罗斯断层对称分布，干涉条纹密集程度和干涉条纹数量代表断层产生的形变值，下盘干涉条纹稀疏，但条纹呈现出圆环状；上盘干涉条纹密集，除中间有圆环，边缘区域有非闭合条纹。对于升轨，干涉纹图只覆盖了主要形变区域。升轨干涉纹图和降轨不同，由于卫星拍摄角度等多种因素，虽然干涉条纹在上盘和下盘均有呈现，但表现为非对称分布。上盘干涉条纹表现密集且闭合，下盘干涉条纹表现稀疏，同样闭合，这是逆冲断层地震的明显表现。另外，在下盘的东南方向也呈现出一个不明显干涉条纹，近似为闭合状。
为了更加直观地分析研究同震形变场，将升、降轨产生的同震形变场应用不同的分层显示方法得到差值色显示。由于2个干涉对产生的同震形变场经过卫星拍摄后的差异较大，为避免出现较大空值，建立独立分层显示图例，最大程度上对地震同震形变场进行还原，结果如图 3。
|
图 3 InSAR同震形变场结果 Fig. 3 InSAR coseismal deformation field |
由图 1可知，高扎格罗斯断层走向为NW-SE，而震中位于该断层附近。结合图 3升、降轨同震形变场，由震中位置和断层位置可以得到，本次地震使地面发生较为严重的变形，上盘发生较大抬升，下盘发生较大下沉，上盘和下盘2个区域以断层为对称轴对称分布。对图 3(b)升轨同震形变场进行精细化分析后发现，除了与降轨得到相似的形变外，在下盘的东南方向产生了一个微小抬升区；对比干涉纹图，在下盘的东南方向同样有一个不明显的干涉条纹。结合地震发震机制可知，该微小形变为下盘走滑产生的地面抬升。
为了更加直观地分析地震产生的地表形变，以断层为基准，连接最大抬升区和最大沉降区作A-B方向的剖面线，如图 3(a)、(b)，生成的剖面图如图 3(c)、(d)。降轨形变场在A-B方向的剖面线显示，地震产生了不同程度的形变值，在断层与剖面的相交点，地震发生的形变值近似为0，在上盘区域发生巨大抬升(最大抬升值为52 cm)，在下盘位置发生明显下降(最大下降值为40 cm)。升轨形变场虽然有相同的形变趋势，但是上盘和下盘的形变值相差很大，在断层的上盘区域抬升87 cm，而下盘仅下降28 cm。
升轨和降轨的结果有很大的差异，其主要原因与卫星轨道数和卫星的拍摄方式有关。InSAR系统观测到的是电磁波入射地球表面后反射(后向散射)的雷达脉冲的强度和相位信息，这些信息编码到雷达坐标系统下(即斜距坐标系)被记录下来，通过地理编码将处理后的InSAR干涉图从斜距坐标系转到地理坐标系，以获得LOS向的形变图。由于卫星升、降轨的视角差异，即降轨时卫星运动方向由北偏东至南偏西，升轨时卫星运动方向由南偏东至北偏西，导致相同的地表三维形变呈现不同的LOS向形变特征。除此之外，SAR卫星对距离向更加敏感，也是造成升、降轨结果差异的一个因素，断层西部上盘在LOS向远离降轨，也远离升轨，表现为较大的抬升；东部下盘在LOS向远离降轨，靠近升轨，表现为降轨的下沉值大于升轨的下沉值。
3 伊拉克地震参数反演 3.1 Okada模型算法Okada[10]提出通用点源及有限面元弹性位错模型，可用于走滑断层、倾滑断层、张裂断层、各向同性及各向异性介质，可以计算弹性半空间条件下地表走滑、倾滑、张裂三分量及梯度变化，从而提供一个通用的定量计算框架。Okada模型非线性反演主要采用Levemberg-Marquardt非线性最小二乘算法进行计算[11]，其公式为：
| $ \begin{array}{*{20}{c}} {{\rm{ CostFunction = }}}\\ {\sum\limits_{j = 1}^M {\left( {\frac{1}{N}\sum\limits_{i = 1}^{{N_j}} {{{\left( {\frac{{{s_i} - {{\tilde s}_i}}}{{{\sigma _i}}}} \right)}^2}} \frac{{{a_j}}}{{\sum\limits_{j = 1}^M {{a_j}} }}} \right)} } \end{array} $ | (1) |
式中，M为数据子集的个数，j为数据子集索引号，i为jth数据子集中单个观测数据i的索引号，αj为jth数据子集的权，Nj为jth数据子集中的测量数据个数，si、ŝi和σi为jth数据子集中的观测位移。预测位移和相应数据的中误差断层面上的错动矢量与笛卡尔坐标系的几何关系如图 4所示。
|
图 4 断层模型的坐标系统 Fig. 4 Coordinate system of fault model |
选取图 4所示的坐标系统，z≤0的部分为弹性空间。取x为断层的走向方向，并且定义U1、U2、U3单位位错分别为走向、倾滑和张裂方向[12]，计算得到面元引起的3个分量的大小，由此可以得到由点源位错引起的位移场和应变场。对于一个长为L、宽为W的矩形断层，作一些替换并进行积分
应用InSAR技术获取升、降轨的同震形变场，分析2个轨道同震形变场后，选择降轨同震形变场反演震源参数。反演震源参数分布主要分为2步：1)对InSAR同震形变场降采样；2)求解非线性反演断层的几何参数，包括经度、纬度、走向、倾向、滑动角、深度、断层的长度、宽度以及滑动量。
通过InSAR求解获得的同震形变场是连续的，其像元点数可以达到6×107。对于Okada模型而言，过多的数据不仅无法提供更多的细节信息，反而会增加计算机负担。因此，在进行参数反演前，对InSAR产生的同震形变场进行降采样处理[14]。
在InSAR获取的同震形变场中选取2个区域(2个区域均覆盖同震形变场)，一个区域为模型形变场，另一个为震源参数反演场。采用数据分辨率约束的四叉树方法[15]，以InSAR获得的同震形变场、视线方位角、视线入射角以及地理编码作为降采样参考系的约束条件，获得约9 040个代表同震形变场的数据点。
在反演断层精细滑动分布时，可通过几个不同参数对震源参数进行描述，包括经度、纬度、走向、倾向、滑动角、深度、断层的长度、宽度、滑动量、震级以及地震矩。获取如此多未知解时，可将反演过程转换为非线性参数求解问题。在非线性反演时，利用CMT官网推荐的地震编号及地震参数推荐区间反复迭代，最后求解地震参数最优值。经过对伊拉克地震几何参数的非线性反演，得到地震震源深度及断层走向、倾角、长度、宽度等参数，结果如图 5、6。
|
图 5 地震参数反演结果 Fig. 5 Seismic parameter inversion results |
|
图 6 Sentinel-1得到的最优断层分布 Fig. 6 Optimal fault distribution of Sentinel-1 |
由图 5、6最优断层滑动模型结果可知，经过反演获取的震源参数和发震机制同由InSAR得到的同震形变场获取的发震机制相符，震中位于扎格罗斯褶皱冲断带，发震机制为逆冲走滑断层，地震形成的破裂面主要位于震中东南方向，产生了较大破裂带，破裂带长度为44.5 km、宽度为16.1 km；破裂带集中在断层两侧，上盘发生上移，下盘有相应下降。此外，地震还伴有走滑，断层的滑动角为140.68°，倾角为8°。将反演得到的震源参数与CMT及USGS公布的震源参数进行比较，如表 2。
|
|
表 2 伊拉克地震最优断层震源参数 Tab. 2 Optimal fault geometric parameters with Iraq earthquake |
从表 2看出，本实验获得的参数与CMT公布的结果有较好的一致性，但与USGS公布的震源深度和倾角有一定的差异，主要原因是研究地震所用的地震波参考标准不同。综合分析可知，本次地震产生长度为44.5 km、宽度为16.1 km的破裂带，震源深度为15.5 km，滑动角为140.68°，断层倾角为8°，震级为MW7.3，累计释放地震矩为0.97×1020 Nm。
分析地震海滩球可以清晰直观地了解地震机制。由表 2中的海滩球可知，公布的地震参数和本实验得到的结果相似。伊拉克地震所在的断层为扎格罗斯主前缘断层，该断层的走向为NW-SE，由此可以断定，本次地震的发震机制为逆冲断层，且下盘同时发生走滑，滑动量为3.83 m，走滑方向由NW到SE，而通过升轨同震形变场的形变分析可以得到相同的结论。
4 结语本文基于Sentinel-1升、降轨SAR数据，利用D-InSAR技术得到LOS向同震形变场，并用四叉树方法对降轨同震形变场降采样作为反演数据，基于Okada模型非线性反演得到发震断层参数，最后结合扎格罗斯造山断裂带断层分布，分析本次地震震源参数。结果显示，伊拉克地震由扎格罗斯前缘断层引起，发震机制为逆冲走滑断层，产生滑动量为3.83 m，从降轨视角得到上盘最大抬升量约为50 cm，下盘最大沉降量约为40 cm。同时由基于Okada反演的震源参数可知，本次地震产生的破裂面长度为44.5 km、宽度为16.1 km，震源深度15.5 km，震级MW7.3，断层走向140.68°，倾角8°，地震矩0.97×1020 Nm。
| [1] |
秦四清, 杨百存, 薛雷, 等. 欧亚地震带主要地震区主震事件判识[J]. 地球物理学进展, 2016, 31(2): 559-573 (Qin Siqing, Yang Baicun, Xue Lei, et al. The Identification of Mainshock Events for Main Seismic Zones in the Eurasian Seismic Belt[J]. Progress in Geophysics, 2016, 31(2): 559-573)
(  0) |
| [2] |
杨百存, 秦四清, 薛雷, 等. 2017年伊拉克MW7.3地震的类型界定及其震后趋势分析[J]. 地球物理学报, 2018, 61(2): 616-624 (Yang Baicun, Qin Siqing, Xue Lei, et al. Identification of Seismic Type of 2017 Iraq MW7.3 Earthquake and Analysis of Its Post-Quake Trend[J]. Chinese J Geophys, 2018, 61(2): 616-624)
( 0) |
| [3] |
左荣虎, 屈春燕, 张国宏, 等. 基于Sentinel-1A数据获取美国纳帕MW6.1地震同震形变场及断层滑动反演[J]. 地震地质, 2016, 38(2): 278-289 (Zuo Ronghu, Qu Chunyan, Zhang Guohong, et al. Coseismic Displacement and Fault Slip of the MW6.1 Nape Earthquake America Revealed by Sentinel-1A InSAR data[J]. Seismology and Geology, 2016, 38(2): 278-289 DOI:10.3969/j.issn.0253-4967.2016.02.004)
( 0) |
| [4] |
McQuarrie N. Crustal Scale Geometry of the Zagros Fold-Thrust Belt, Iran[J]. Journal of Structural Geology, 2004, 26(3): 519-535 DOI:10.1016/j.jsg.2003.08.009
( 0) |
| [5] |
文磊, 张光亚, 李曰俊, 等. 扎格罗斯褶皱冲断带构造变形特征[J]. 地质科学, 2015, 50(2): 653-664 (Wen Lei, Zhang Guangya, Li Yuejun, et al. Structure-Deformation Features of the Zagros Fold and Thrust Belt[J]. Chinese Journal of Geology, 2015, 50(2): 653-664 DOI:10.3969/j.issn.0563-5020.2015.02.020)
( 0) |
| [6] |
张洪瑞, 侯增谦. 伊朗扎格罗斯造山带构造演化与成矿[J]. 地质学报, 2015, 89(9): 1560-1572 (Zhang Hongrui, Hou Zengqian. Tectonic Evolution and Metallogeny of Zagros Iran[J]. Acta Geologica Sinica, 2015, 89(9): 1560-1572 DOI:10.3969/j.issn.0001-5717.2015.09.002)
( 0) |
| [7] |
Ambraseys N, Moinfar A, Peronaci P. The Seismieity of Iran the Farsinaj (Kermanshah) Earthquake of 13 Decemher 1957[J]. Annals of Geophysics, 1973, 26(4): 679-692
( 0) |
| [8] |
Niazi M, Asudeh I, Ballard G, et al. The Depth of Seismicity in the Kermanshah Region of the Zagros Mountains (Iran)[J]. Earth and Planetary Science Letters, 1978, 40(2): 270-274 DOI:10.1016/0012-821X(78)90097-3
( 0) |
| [9] |
杨亚夫, 朱建军, 王永哲, 等. 利用Sentinel-1A数据、D-InSAR和沿轨干涉技术获取2016年高雄MS6.7地震三维形变场[J]. 大地测量与地球动力学, 2017, 37(4): 339-343 (Yang Yafu, Zhu Jianjun, Wang Yongzhe, et al. 3D Displacements Field of 2016 Kaohsiung MS6.7 Earthquake form D-InSAR and along-Track Interferometry with Ascending and Descending Sentinel-1A Timages[J]. Journal of Geodesy and Geodynamics, 2017, 37(4): 339-343)
( 0) |
| [10] |
Okada Y. Surface Deformation Due to Shear and Tensile Faults in a Half-Space[J]. Bulletin of the Seismological Society of America, 1985, 75(4): 1135-1154
( 0) |
| [11] |
张明.利用空间大地测量资料反演玉树地震带区域形变场特征[D].西安: 长安大学, 2014 (Zhang Ming. Inversing Deformation Characteristics of Yushu Seismic Zone with the Space Geodesy Information[D]. Xi'an: Chang'an University, 2014) http://cdmd.cnki.com.cn/Article/CDMD-10710-1014070896.htm
( 0) |
| [12] |
任小冲.基于InSAR的台湾集集地震形变场提取和震源参数确定[D].长沙: 中南大学, 2008 (Ren Xiaochong. Using InSAR to Extract the Deformation Field and Determine the Source Parameters of Chi-Chi Earthquake[D]. Changsha: Central South University, 2008) http://cdmd.cnki.com.cn/Article/CDMD-10533-2008167701.htm
( 0) |
| [13] |
凌勇, 曾祺明, 罗扬, 等. 巴姆地震变形场和应力场: Ⅰ.用差分干涉雷达和Okada方法求解[J]. 岩石学报, 2006, 22(9): 2367-2374 (Ling Yong, Zeng Qiming, Luo Yang, et al. Deformation and Stress Field of Bam Earthquake:Ⅰ.InSAR and Okada Calculation[J]. Acta Petrologica Sinica, 2006, 22(9): 2367-2374)
( 0) |
| [14] |
Lohman R B, Simons M. Some Thoughts on the Use of InSAR Data to Constrain Models of Surface Deformation: Noise Structure and Data Downsampling[J]. Geochemistry, Geophysics, Geosystems, 2005, 6(1)
( 0) |
| [15] |
Jónsson S, Zebker H, Segall P, et al. Fault Slip Distribution of the 1999 MW7.1 Hector Mine, California, Earthquake, Estimated from Satellite Radar and GPS Measurements[J]. Bulletin of the Seismological Society of America, 2002, 92(4): 1377-1389 DOI:10.1785/0120000922
( 0) |