多视角InSAR数据解算2017两伊地震三维同震形变场

韩 鸣，张永志，程 冬，尹 鹏

(长安大学地质工程与测绘学院，陕西 西安 710054)

传统的大地测量技术包括GPS、水准测量等采集的数据都是离散分布的，并且受时空条件的影响，难以准确和及时获取大面积地形测量数据[1]。差分合成孔径雷达干涉测量(DInSAR)技术是一种快速发展的大地测量技术，能够全天候获取地面高精度、大面积、低成本的地表变形信息，在地质灾害监测和评估方面已经取得一系列成果。

在各种地质灾害中，地震和火山运动及地表沉降和山体滑坡等均为需要利用测绘技术研究的问题[2]。DInSAR技术通过处理两幅或两幅以上的雷达影像，提取地表形变相位差值，计算得到厘米级甚至毫米级的地表沿视线方向(line-of-sight,LOS)的形变信息。利用覆盖发震区域的震前和震后两幅雷达影像，通过DInSAR结合数字高程模型数据DEM处理就能获取视线向的同震形变场。在伊拉克和伊朗边境发生的Mw7.3级(USGS发布)地震，美国地质调查局(USGS)研究认为，该地震的震源机制主要是逆冲兼具轻微右旋走滑地震，为了更深入地研究该地震，本文利用InSAR技术及升降轨3个视角的雷达影像数据，重建了此次地震的三维同震形变场，为更深入研究和认识该地震发震机制提供必要的基础。

1 地质构造背景

2017年11月12日在伊拉克的苏莱曼尼亚省和伊朗的克尔曼沙阿省之间发生了Mw7.3级地震，造成了至少530人死亡，超过7200人受伤。震中位置为如图1所示的五角星处(34.911°N，45.959°E，USGS)[3]。研究区属于伊朗高原，伊朗高原由许多地形起伏较大的陆地块体边界组成，阿拉伯板块以相对于欧亚板块20 mm/a的速度向北移动，使阿拉伯板块和欧亚板块持续碰撞，形成了扎格罗斯造山带，如图1所示。图1中的曲线是阿拉伯板块与欧亚板块的边界线[4]。碰撞带边界的斜向应力聚集引起了此次地震，发震区域处于阿拉伯-欧亚板块碰撞带的中心，紧邻扎格罗斯山，USGS公布的震中位于扎格罗斯山前断层(ZMFF)附近，这是一条逆冲断层,图1中的小圆圈为余震分布，大部分余震也都发生在扎格罗斯山前断层(ZMFF)附近[5]。

图1 研究区震中与地质构造背景

2 三维同震形变场解算方法

视线向形变场由两幅SAR影像干涉处理得到，同一地区由于卫星的飞行方向和雷达照射方向的不同，处理得到的视线向形变场也是不同的，甚至是相反的[6]。卫星雷达传感器成像模式分为升轨和降轨，如图2所示。三维形变分别为北向形变，东向形变和垂直向形变，任何变形都可以通过三维形变来合成。但是3个方向的形变对视线向形变的贡献各不相同[7]，Sentinel1-A卫星的视线方向为垂直于卫星飞行方向的右下方，因此雷达传感器对垂直向运动最敏感；而卫星的飞行方向是近南北向的，如果地震断层引起的地表运动是近NS向的，则将对InSAR观测结果非常不利，甚至观测不到任何形变信息，比如巴姆地震和美国的Hector Mine地震[8-9]。根据雷达成像的几何关系，影像干涉得到的视线向形变用dLOS表示，约定地表接近卫星的dLOS值为正，远离卫星的dLOS值为负，三维形变(垂向dU,北向dN,东向dE)与视线向形变dLOS之间的几何关系用公式表达为

dLOS=dUcosθ-dNsinθcosα-dEsinθsinα

(1)

图2 卫星轨道方向

式中,dLOS为InSAR视线向观测值；dU为地表垂向形变；dN为地表南北向形变；dE为地表东西向形变；θ为雷达入射角；α为方位视线方向，即距离向与北向的夹角，该式同时适用于升降轨的投影计算。由式(1)可知，要获取真实的地表三维形变,至少需要3个不同视线向的观测量。发生地震的区域有3对不同成像视角的雷达影像时，即可利用这3对影像数据得到的视线向形变场进行直接解算地表三维形变[10]，根据式(1)有

(2)

令

(3)

3 视线向同震形变场提取与分析

3.1 数据选取与处理流程

Sentinel1-A是由欧空局2014年4月发射的对地观测卫星，在近地太阳同步轨道上运行，重访周期为12 d，选取Sentinel1-A卫星的干涉宽幅(interferometric wide swath,IW)模式，对地表进行递进扫描得到3个字条带，该成像模式采用中等分辨率(5 m×20 m)，对地面数据提取的范围达到240 km[11]。本文选取了升降轨不同的3对影像，共6景数据，覆盖范围如图1中黑色虚线框所示，基本信息见表1。选取的DEM数据为SRTM航天飞机雷达地形测量任务数据，地面分辨率为90 m。

表1 Sentinel1-A(来源ESA)卫星数据基本信息

常规二通差分DInSAR处理流程为影像配准、影像干涉、去平地相位、干涉图滤波、相位解缠、地理编码[12]。Sentinel1-A数据覆盖范围较宽，本文裁剪了比震区范围稍大的部分影像进行处理，在处理过程中使用地面分辨率为90 m的SRTM数据去除干涉图中的地形相位分量，研究区处于植被覆盖率低的中东地区，忽略由地表植被引起的干涉失相干影响[13]。为了尽可能消除干涉图的噪音，先对干涉影像进行距离向4视数，方位向1视数的多视处理，得到距离向分辨率为16.621 8 m，方位向分辨率为13.968 m，利用Goldstein滤波方法进行滤波,其滤波器是可变的，提高了干涉条纹的清晰度，减少了由空间基线或时间基线引起的失相干的噪声[14-15]。采用最小费用流法(MCF)进行相位解缠，设定相干性阈值，对相干性低于0.2的区域进行掩膜处理。

3.2 视线向同震形变场提取分析

利用最终的解缠结果，转换为视线向的形变信息，获取升降轨3种视角的干涉图和视线向形变信息。如图3—图5所示，分别为差分干涉条纹图和视线向形变图，图中F1—F2虚线为USGS推测出的发震断层顶部，3幅影像的重叠区没能完全覆盖形变区，但仍涵盖大部分形变区域，黑色虚线框为选取的三维形变解算区。升轨-A72的同震形变场视线向隆升中心集中在西南部，最大隆升量为84.7 cm，视线向沉降值最大为-26 cm；降轨-D79的同震形变场呈两个扇形沿NW向对称分布，西南部为视线向隆升区，最大隆升量为58.55 cm，东北部为视线向沉降区,最大沉降量为-41.88 cm；降轨-D6的同震形变场的形变分布范围和D79一致，隆升量与沉降量都比降轨-D79的结果稍小。

图3 升轨-A72视线向形变

图4 降轨-D79视线向形变

图5 降轨-D6视线向形变

图3—图5中设置线条AB为视线向形变值剖面线，剖面提取信息如图6所示,升轨-A72的视线向隆升量和沉降量与两个降轨相比差别较大，是由于成像视角不同造成的，但3种视角下的视线向上升量与沉降量总体趋势基本一致。3对影像的震后成像时间差值最大不超过5 d，可以忽略震后余滑的影响，认为是对同一同震形变场的3种不同的观测成果。

图6 沿AB剖线的升降轨视线向形变

4 三维形变结果分析

利用3种视角的视线向观测值，通过式(3)重建了研究区三维同震形变场，计算过程中选用表1中所列的平均入射角和平均视线方位角，计算结果如图7所示。图7中主要的形变区分为两个部分，西南部的主要形变区域称为C区，东北部的主要形变区域称为D区。图7(a)显示，C区向上隆升，最大上升量达到了84 cm，D区沉降，最大沉降量为-34.5 cm；图7(b)显示，C区和D区都朝西运动，最大形变量为-50 cm;图7(c)显示，形变区整体向南运动，最大形变量为-91 cm，运动主要集中在形变区南部。结合三维形变场特征，得出整个形变区呈向上抬升且在水平方向上向西南方向运动的特点。通过上文所述，本文可以通过几条线索来确认发震主断层：第一，视线向同震形变场和三维同震形变场的形变边界约束位置紧邻扎格罗斯山前断层(ZMFF)；第二，主震和余震都发生在扎格罗斯山前断层(ZMFF)周围；第三，主震具有低倾角及纯逆冲运动的特征，与扎格罗斯山前断层的运动特征一致。结合地质构造和形变分析，推测此次地震的发震断层为扎格罗斯山前断层(ZMFF)。

图7 三维同震形变场

5 结 语

本文利用Sentinel1-A卫星的3对升降轨SAR数据，通过二轨差分DInSAR技术处理得到2017年两伊地震区域的3种视角视线向同震形变场，将3幅干涉图重叠区域裁剪为三维解算区，利用直接解算法得到研究区的三维形变场。通过试验得到以下结论：①形变区域的视线向同震形变场中，升轨的最大隆升和最大沉降分别为84.7和-26 cm，降轨的最大隆升和最大沉降分别为58.55和-41.88 cm，且3种视线向的形变场总体形变趋势基本一致；②解算得到的三维形变场影响范围大致为75 km×80 km，分为C区和D区两个主要形变区，最大的垂直向、东西向、南北向的形变值分别为84、50、-91 cm，重建的三维形变场的显示研究区呈现整体抬升并朝西南向运动的特点；③综合地质构造背景和视线向同震形变场以及三维形变场的形变边界和趋势等因素，本文推测发震断层为扎格罗斯山前断层。