升降轨D-InSAR技术提取新疆伽师地震同震形变场

2021-05-07丛康林

矿山测量 2021年2期

张迪，丛康林，董超

(1.山东农业大学信息科学与工程学院，山东泰安 271018；2.南京信息工程大学遥感与测绘工程学院，江苏南京 210044)

合成孔径雷达差分干涉测量技术(D-InSAR，Differential Interferometric Synthetic Aperture Radar)是以InSAR技术为基础并进行延伸，差分干涉技术是指对同一地区的两幅不同时相的影像进行差分干涉处理，从形成的干涉图中去除地球曲率等因素的影响[1]，利用外部DEM信息消除地形因素的影响，据此探测到地表形变信息[2]。D-InSAR技术较传统的测量方法，具有较高的精准度、灵敏度和分辨率、全天时全天候、不受烟云雨雾的限制等优点，且不需要布设控制点，所以方便在危险地区成像，已经广泛应用于监测地震引起的地表形变[3]。北京时间2020年1月19日，在新疆维吾尔自治区伽师县发生了Ms6.4地震，震中位置为39.83°N、77.21°E，震源深度为16 km。本文通过对该次新疆伽师地区震前震后的升降轨影像分别处理，利用ENVI+SARscape软件，运用“二轨法”，结合精密轨道数据，获取雷达LOS方向的形变，再将升降轨数据结合成水平和垂直的形变量，最终反演出新疆伽师地区此次地震的同震形变场，对断裂结果进行位错分析，发现发震断层，探测隐藏的地震活动，定量地进行灾后评估，对于灾情分析与灾后重建具有重要意义[4]。

1 研究区概况

伽师县位于新疆维吾尔自治区西南部喀什地区，地形为盆地，地势西高东低，南高北低，属于温带大陆性气候，冬冷夏热，湿度低，昼夜温差较大。伽师一面环山，一面环湖，处于西昆仑山地震带和南天山地震带的交界地带，为地震多发区，五年来该地区发生的三级以上的地震多达180次，自二十世纪中叶起，震级超过六级以上的多达20余次，最强的一次地震震级高达8.25级。伽师地区的地理位置如图1所示。

图1 新疆伽师县地理位置图

2 数据

2.1 Sentinel-1数据

哨兵一号(Sentinel-1)是欧洲航天局哥白尼计划中的对地观测卫星，轨道的高度约为693 km，由位于同一轨道的A、B两颗星组成，相位差为180°，哨兵一号的IW模式包含三个sub-swaths，而每个sub-swaths又包含若干个burst，每个burst与相邻burst都在距离向与方位向上有一定的重叠度[5]。哨兵一号携带12 m长的C波段合成孔径雷达天线，可以向全球用户提供连续不间断的图像，其卫星影像数据可用于地震等自然灾害监测、地表形变监测、海洋环境监测、海上安全等方方面面。本文选用的哨兵一号SAR数据信息如表1所示，空间基线远远小于临界基线，适合于反演同震形变场[6]。

表1 SAR像对参数

2.2 精密轨道数据

精密轨道数据(POD，Precise Orbit Ephemerides)是D-InSAR技术中至关重要的信息，贯穿于整个处理过程，例如图像配准干涉图的生成以至于最后形变图的生成，可以修正轨道的误差，有效地去除系统误差，使差分干涉的结果更加精确。精密轨道数据的精度很高，但是只在GNSS下行21天之后才可以正常获取，定位精准度优于5 cm。

2.3 SRTM 90 m DEM数据

数字高程模型(DEM，Digital Elevation Model)，作为一种实体地面模型，是将地面划分为等间隔的格网，取每个网点的高程值作为地面高程值[7]，广泛地应用于等高线的生成、正射纠正等。SRTM(Shuttle Radar Topography Mission)是一种对地球表面进行遥感测量的航天飞机雷达地形测绘任务，SRTM数据是由美国太空总署(NASA，National Aeronautics and Space Administration)和国防部国家测绘局(NIMA，National Imagery and Mapping Agency)联合测量[8]，由美国发射的“奋进”号航天飞机搭载的SRTM系统完成，此次测图任务从2000年2月11日开始，至22日结束，仅仅利用航天飞机便可获得全球大部分的分辨率为30/90 m的外部DEM数据，绝对垂直精度为16 m[9]，相对水平精度为10 m，为D-InSAR技术提供了精确可靠的DEM数据[10]。DEM的数据分辨率主要有30 m、90 m和1 km，本次实验选用分辨率为90 m的DEM数据。

3 方法

D-InSAR技术的方法主要有三种：“二轨法”、“三轨法”、“四轨法”。已知一个地区的DEM数据，将跨越地表形变的两幅不同时相的影像进行干涉处理，结合DEM数据去除地形因素，再结合精密轨道数据修正卫星轨道参数信息，即可获得地表的形变量，这种差分干涉处理方法称为“二轨法”[11]，处理流程如图2所示。

图2 “二轨法”处理流程

4 数据处理分析

4.1 数据处理

SARscape软件是由sarmap公司研发的专业雷达图像处理软件，可以用于处理Sentinel-1数据。关键步骤如下：

(1)生成干涉图

干涉图的生成与差分的精度密切相关。要生成干涉图，首先要将两景影像基于轨道信息配准，然后用SLC数据所对应的像元复数共轭相乘即可得到干涉图[12]。此时得到的干涉图是包含平地效应的，再结合SRTM 90 m DEM数据，基于精密轨道数据进行去平处理，得到去除平地效应的干涉图。

(2)滤波与生成相干系数

去除平地效应后的干涉图包含各种渠道所产生的噪声，影响解缠的精确度。Goldstein滤波方法可以提高干涉条纹的清晰度，有效地减少了失相干噪声。由于研究区位于我国西部，植被覆盖率较低，所以此次实验选取的两对像对的相干性都很高，接近于1。

(3)相位解缠

相位解缠的主要目的是为了将去平和滤波后产生的相位去除。主要的相位解缠方法有区域增长法(Region Growing)，最小费流法(Minimum Cost Flow)，分解等级法(Unwrapping Decomposition Level)[13]。本次实验采取最小费流法，其优势在于所有的像元都可以被考虑，尤其是在低相干范围较大的区域。

(4)相位转形变以及地理编码

把经过绝对校准和解缠的相位，结合成合成相位，转化为形变数据以及地理编码到制图坐标系统，默认得到的LOS方向的形变如图3所示。将形变量进行统计，得到如图4所示的地表变化分布曲线，分别沿着抬升区、沉降区平行于断层方向(a1、a2、b1、b2)，整个形变区垂直于断层的方向(c1、c2)绘制剖面线，得到如图5～图7所示的剖面图。

图3 LOS方向地表形变图

图4 地表形变统计分布曲线

图5 抬升区剖面图

图6 沉降区剖面图

图7 形变区南北向剖面图

4.2 结果分析

对同一个地区，采用相同的方法进行D-InSAR处理，升降轨影像得到的形变场在范围和量级上基本一致，由于航向和视线方向的不同，以及对地形的敏感程度不同，导致视线向的形变存在微小的差异。如图3～图7所示，此次形变主要集中在震中的东北方向的两瓣区域，南瓣表现为抬升，北瓣表现为沉降。就抬升区来说，升降轨的LOS向形变量位置和范围大小大致相同，最大抬升量约为6 cm。而对于沉降区，升轨的处理结果要比降轨的处理结果范围大得多，并且在数值上也有1.4 cm的差异。由图4地表形变分布曲线可知，绝大多数的像素形变量分布在±1 cm之间，少数的像素点发生较大的形变。如图5～图6所示，抬升区的剖面曲线相对比较光滑连续，而沉降区的剖面图存在很多“跳点”，分析可知，南部形变区呈现均匀抬升的趋势，而北部沉降区呈现不均匀沉降的趋势，造成地表破碎现象。

将升降轨LOS向结合成垂直和水平方向的形变场，如图8、图9所示。并沿着垂直于断层方向由南向北(d1、d2)绘制剖面线，得到如图10、图11所示的剖面曲线图。从垂直和水平方向的同震形变场可以得到如下结论：(1)垂直形变场的剖面图(图10)表明，沿着垂直于断层方向呈正负分布，上盘表现为下降，下盘表现为抬升，可以判定属于正断层，且下盘上升的趋势要远远大于上盘下降的趋势。(2)水平形变场的剖面图(图11)表明下盘向西运动，上盘向东运动，且有一个陡坡，侧移量先减少后缓慢增大。(3)由垂直形变场(图8)可知，整个垂直形变场由两个近似相同的椭圆呈南北对称分布，东西方向长度约为50 km，南北方向约为35 km，下盘表现为抬升，抬升范围为0～7.67 cm，上盘表现为沉降，沉降范围为0～-3.79 cm，颜色过渡较为平缓。(4)由水平形变场(图9)可知，上盘向东运动，下盘向西运动，且拉张的范围大致相同，约为3 cm。(5)在抬升区与沉降区之间形成了一个断层，方向为“西北-东南”，长度大约为28 km。地震波及的两处形变区颜色过渡比较平滑，结合图10跨越断层的剖面线可以发现，垂直位移分布较连续，没有出现明显的条纹断裂或者失相干现象，这表明此次地震并没有造成明显的地面破裂现象，即本次地震破裂并未到达地表。