基于时序InSAR技术的中巴走廊 盖孜河谷段地表形变分析*

汪钰晴，唐伶俐，王新鸿†，周增光，李子扬，李传荣

(1 中国科学院空天信息创新研究院 中国科学院定量遥感信息技术重点实验室， 北京 100094； 2 中国科学院大学电子电气与通信工程学院， 北京 100049) (2020年6月12日收稿； 2020年7月6日收修改稿)

中巴经济走廊(简称“中巴走廊”)北起中国新疆喀什，南至巴基斯坦瓜达尔港，全长约3 000 km，连接着北部的“丝绸之路经济带”和南部的“21世纪海上丝绸之路”，是贯通南北丝路的一条重要贸易走廊。中巴走廊沿线滑坡、泥石流、崩塌、溜石坡等地质灾害频发，不仅给沿线人民生活造成极大不便，也给走廊的工程建设和运营带来了很大的挑战。然而，由于该地区地形地质条件复杂、自然环境恶劣、安全状况差、交通不便，不适合开展长期的野外调查活动，导致该区域相关研究工作进展缓慢。

20世纪70年代发展起来的卫星合成孔径雷达干涉测量技术(interferometric synthetic aperture radar，InSAR)是一种极具应用潜力的微波遥感技术，具有监测范围大、空间分辨率高、几乎不受云雨天气制约等优点，已在地震形变[1-2]、火山运动[3]、冰川漂移[4]、城市沉降[5-6]、山体滑坡[7-8]以及线性工程形变[9]等研究领域得到广泛应用。得益于星载SAR的发展，在相同研究区可积累大量的重复轨道SAR数据，这为开展长时间序列InSAR地表形变反演研究提供了可能。自本世纪初意大利米兰理工大学Ferretti提出永久散射体方法后，一系列基于点目标的时间序列分析方法先后被提出，包括永久散射体(permanent scatterers，PS)方法[10]，小基线集(small baseline subsets，SBAS)方法[11]，IPTA(interferometric point target analysis)方法[12]，StaMPS(stanford method for permanent scatterers)方法[13]等，这些方法通过分析时间序列SAR影像，提取出相位和幅度信息在长时间范围内保持稳定的离散点，即PS点或相干目标候选点(例如城市建筑、道路、裸露的岩石等地物)，利用这些点目标的相位特征，进行长期缓慢地表形变的反演。时序InSAR技术不仅可以克服传统差分合成孔径雷达干涉测量技术(differential InSAR，D-InSAR)易受时间、空间去相干影响的问题，还能在一定程度上估计大气影响和优化初始数字高程模型(digital elevation model，DEM)精度，达到抑制大气延迟相位和地形误差相位的效果。大量研究表明，时序InSAR技术能提取大范围、高精度的长时间缓慢地表形变信息，可广泛应用于城市地面沉降[14]、潜在滑坡形变监测[15-16]、地质灾害识别[17]等领域。

本文选择中巴走廊沿线地质灾害最为发育的盖孜河谷段为研究区，利用2015年4月5日至2018年12月9日共30景Sentinel-1A干涉宽幅模式下的雷达影像，分别采用PS-InSAR技术和SBAS-InSAR技术对研究区进行地表形变信息提取，比较两种时序InSAR技术所得形变序列结果，对地质灾害形变序列特征进行初步分析。

1 研究区概况与数据

1.1 研究区概况

中巴走廊穿过帕米尔高原腹地，地质构造上属于印度-巴基斯坦板块与欧亚板块强烈碰撞和挤压的区域，是地壳隆起迅速、地震活动强烈、沟谷侵蚀剧烈的典型地区之一[18]。因构造运动剧烈、地势陡峭、冰川广泛发育，中巴走廊沿线地质灾害频发，总体上具有灾种多、规模大、频度高、分布广的特点[19-20]。

盖孜河谷横切西昆仑山脉，北靠盖昆山(最高海拔6 670 m)，南依公格尔山(最高海拔7 649 m)，河谷大致呈北东东(NEE)走向，平均海拔2 500 m，最宽处不到2 km，最窄处不足50 m，是中巴走廊的咽喉地段。盖孜河谷位于欧亚大陆腹地，降水相对稀少而蒸发较为强烈，属于典型的大陆性气候。由于气候寒冷干旱、冰川广泛发育、风化作用强烈，该区域植被稀少、山体裸露。滑坡、崩塌、泥石流等地质灾害的频繁发生，给走廊建设和人民生活带来很大威胁[21]。

使用星载SAR数据进行时序InSAR处理，考虑到整景SAR影像覆盖范围较大，为提高数据处理效率，在进行时序InSAR处理前先对SAR影像做空间位置裁剪，以提取感兴趣的盖孜河谷段区域，裁剪范围为74°58′～75°30′ E，38°40′～38°53′ N(约45 km×25 km)，所得研究区地理位置如图1所示。

图1 研究区地理位置Fig.1 Geographical location of the study area

1.2 SAR数据

选择Sentinel-1A卫星SAR数据作为数据源。Sentinel-1A卫星于2014年4月3号发射，采用太阳同步轨道(轨道高度693 km)，具有条带、干涉宽幅、超宽幅及波浪4种成像模式[22]，主要参数如表1所示。

以中巴走廊盖孜河谷段部分区域为研究区，共选取30景Sentinel-1A卫星SAR影像(干涉宽幅(IW)模式、升轨、VV极化)。影像获取时间如表2所示，时间范围为2015年4月5日至2018年12月9日。

表1 Sentinel-1A卫星主要参数Table 1 Main parameters of Sentinel-1A

表2 SAR数据时间分布表Table 2 Date distribution of the SAR data

此外，Sentinel-1A卫星发布了精密定轨星历参数(precise orbit ephemerides，POD)用于修正轨道信息，定位精度优于5 cm，将它作为辅助数据可以有效去除轨道误差对形变测量的影响。

1.3 DEM数据

采用空间分辨率为30 m的ASTER GDEM V2数据(下载网址https:∥www.gscloud.cn)作为InSAR处理的辅助数据。该数据集根据NASA的EOS-Terra卫星的全球观测结果制作而成，其数据覆盖南北纬83°之间的所有陆地区域。该高程数据集的垂直精度为20 m，水平精度为30 m，可在InSAR处理过程中较好地模拟地形相位。

2 地表形变提取方法

分别采用PS-InSAR技术和SBAS-InSAR技提取研究区地表形变信息，术总体流程如图2所示。

PS-InSAR技术提取地表形变的主要步骤包括：主影像选择、差分干涉处理、PS候选点选择、估算候选点的平均形变速率及DEM校正系数、大气相位的估计和去除、形变时间序列估计等。本文选取2017年3月13日的SAR影像作为公共主影像，其他29幅影像作为辅影像。图3(a)为干涉对的时间基线和空间基线连接图，从中可见主辅影像干涉像对空间基线分布于-100～105 m，时间基线范围为-708～636 d。选取公共主影像之后，将其余29幅辅影像配准重采样到主影像的像素空间上并进行主辅影像干涉。利用Sentinel-1A卫星精密定轨星历参数与辅助DEM数据模拟并去除参考椭球面相位及地形相位后，得到29幅差分干涉图。

采用振幅离差指数方法，将振幅离差指数阈值设置为0.4，对研究区内PS点进行识别。PS点的干涉相位包含线性形变、高程误差、非线性形变、大气延迟和噪声等相位分量。首先用线性模型从所有差分干涉图中估算线性形变速率和残余高程信息，再在时域和空域上做相应的滤波处理将非线性形变相位分离出来，将它与线性形变分量相加，得到每个PS点的总地表形变信息。最后，将视线向形变进行地理编码，获得地理坐标系下的地表时序形变信息。

图2 研究区地表形变提取流程图Fig.2 Flow chart of surface deformation extraction in the study area

图3 PS及SBAS处理中SAR数据对的时空基线连接图Fig.3 Temporal-spatial baseline connection diagram of SAR data pairs in PS and SBAS processing

SBAS-InSAR技术提取地表形变的主要步骤包括：连接图生成、差分干涉处理及相位解缠、轨道精炼与重去平、大气相位估计与去除、形变时间序列估计等。本文选取2015年8月27日的SAR影像作为公共主影像。经过主辅影像配准得到像空间坐标系一致的30幅SLC影像后，设置时空基线阈值。本文使用的SAR数据总时间跨度为1 344 d，为保证干涉质量，将空间最大临界基线百分比设置为45%，时间基线阈值设置为180 d，一共生成94个干涉对，最长的空间基线约为210 m。图3(b)为生成的时空基线连接图，这里每幅影像平均有6.27个连接，连接得比较均匀，且每个时相都与其他时相建立了连接，表明数据配对良好。

对满足时空基线阈值的所有干涉像对进行差分干涉处理，基于卫星精密定轨星历参数与辅助DEM数据，去除参考椭球面相位及地形相位，然后对差分干涉图进行滤波处理，改善干涉条纹的清晰度，之后通过相位解缠得到2次观测期间同一目标的真实相位差。最后经过轨道精炼与重去平、形变速率与DEM校正速率估计、大气效应相位估计与去除等处理，得到所有PS点的年平均形变速率及时间序列形变信息。本文采用ENVI软件的SARscape模块完成上述时序InSAR处理过程。

3 结果与分析

3.1 PS-InSAR与SBAS-InSAR形变结果比较分析

图4(a)和4(b)分别为利用PS-InSAR及SBAS-InSAR技术获得的研究区PS点分布及其年平均形变速率图，结合光学影像(这里的背景底图)可以看出，两种方法提取PS点的分布范围较为一致，主要分布在海拔较低的中巴公路沿线、西部盆地及坡度较缓山区，在冰雪覆盖区基本无PS点，受Sentinel-1A卫星雷达波入射角和入射方向影响，由于阴影效应，雷达波束照射不到的北东(NE)向陡峭坡面上也没有形成PS点。

PS-InSAR处理时全区共提取了305 853个PS点，所测形变为传感器观测方向的形变，即雷达视线向(line of sight，LOS)形变量(地面向卫星运动方向为正，远离卫星运动方向为负)，PS点的LOS形变速率主要分布于-98～98 mm/a区间范围内。PS点的形变速率直方图大致为正态分布，如图5(a)所示，均值约为-0.32 mm/a，标准差约为5.25 mm/a。从图4(a)可以看出，研究区西部盆地及中巴公路沿线大部分PS点年平均形变速率都较小，表示这些区域在监测时间段内(2015年4月—2018年12月)未发生明显形变；而在中巴公路沿线部分路段及公路两侧冰川前缘，分布有年平均形变量较大的PS点(图4(a)中矩形区域)，表示这些区域在监测时间段内可能发生了由于地质灾害而引发的形变。此外，提取出的PS点中还存在一些孤立的年平均形变速率绝对值异常大的点簇(图4(a)中椭圆形区域)，由于研究区地处高寒干旱山区，高程落差大且气候变化明显，这些LOS形变异常高值PS点可能是在PS-InSAR数据处理过程中由于大气相位或DEM残差相位分离不当造成的。

基于相同的30幅雷达影像，在同一区域采用SBAS-InSAR技术共提取了670 757个PS点，测得PS点的LOS年平均形变量主要位于-46～38 mm/a区间范围内。PS点的形变速率直方图接近正态分布，如图5(b)所示，均值约为0.15 mm/a，标准差约为2.86 mm/a。从图4(b)可以看出，研究区西部盆地及中巴公路沿线大部分PS点在监测时间段内保持稳定，仅在公路沿线部分斜坡及两侧冰川前缘分布有形变量较大的PS点，这与PS-InSAR所得结果较为一致，两种时序InSAR技术可以相互印证。

图4 PS-InSAR及SBAS-InSAR技术提取的盖孜河谷LOS形变速率图Fig.4 LOS deformation rate map of the Gaizi valley by PS-InSAR and SBAS-InSAR techniques

图5 两种时序InSAR技术所提取PS点的形变速率分布图Fig.5 Distribution of LOS velocity obtained by the two time series InSAR techniques

值得注意的是，PS-InSAR技术是通过计算每个像素点的振幅离差指数来选择PS候选点，进而通过设定多时间相干系数阈值得到所有PS点的，这里多时间相干系数表示有多少形变趋势符合线性模型，该阈值设置越大，所得PS点越少(本文将此阈值设定为0.8)。与PS-InSAR技术利用强度信息确定PS点不同，SBAS-InSAR技术是通过设定相关系数阈值来获取高相干点的，本文处理过程中将相关系数阈值设置为0.35，低于该阈值的像素以空值输出。由于PS点的识别方法不同，这两种时序InSAR技术所获得的PS点在数量及分布上有较大差距。

相较而言，SBAS-InSAR技术更充分地利用了小基线数据子集内的良好相干性，提取出的时序形变结果在空间上更为连续，减少了出现形变异常结果的可能性；此外，在测量形变类型上，PS-InSAR主要适用于线性形变，而SBAS-InSAR对线性及非线性形变均适用。因此，下面将基于SBAS-InSAR提取结果对研究区部分潜在地质灾害进行形变特征分析。

3.2 研究区地质灾害形变特征分析

结合PS点正态分布图、PS点误差等综合分析，可将形变速率划分为7个区间，其中极小值-46～-20 mm/a与极大值20～38 mm/a用来识别LOS形变异常高值点，中间区间-5～5 mm/a可以看作背景值，表示该区域地表较为稳定，其余4个区间：-20～-10、-10～-5、5～10、10～20 mm/a可以用来标记产生一定程度形变的区域。

从图6形变速率划分结果可发现，研究区西部盆地与中巴公路沿线大部分区域年平均形变速率均不超过±5 mm/a，表示研究区大部分区域在监测时间段内保持稳定，未发生明显地表形变；而在中巴公路沿线两侧斜坡上，识别出较多年平均形变量大于±5 mm/a的PS点，由于公路沿线两侧山体植被稀少、风化作用强烈，易产生崩塌、溜石坡等斜坡灾害；在布伦口南北两侧冰川前缘区域，分布有大量年平均形变速率较大PS点，部分PS点形变速率超过±20 mm/a，这是由冰川运动或冰川融水引发泥石流所造成。

结合研究区光学影像及前人相关研究成果和经验，根据地表LOS形变速率提取出研究区内多处不稳定斜坡及冰川运动，见图6。其中形变速率较大的不稳定斜坡主要分布在中巴公路两侧的坡面上(图6中椭圆区域)，存在大量不稳定PS点，个别不稳定斜坡与公路距离较近，可能发生损毁道路的安全隐患。公格尔山是研究区内海拔最高的山峰，冰川发育较为显著，由于冰川运动，在冰川前缘提取到较多形变速率相对高值PS点(图6中矩形区域)，此外，在公路北侧木吉县方向也有比较明显的冰川运动。

图6 盖孜河谷区域SBAS-InSAR处理所得PS点分布及形变速率Fig.6 Distribution and velocities of PS points obtained by SBAS-InSAR processing in the Gaizi valley

3.2.1 不稳定斜坡形变特征

研究区构造运动强烈、河流切割侵蚀、地震发生频繁，且气候高寒、冰川广泛发育，公路沿线两侧产生了大量陡坡，易发生滑坡灾害，对公路存在潜在危害。由地震、降雨、风化而引发的崩塌也是研究区主要地质灾害之一，规模一般较小，但对公路本身及过往人员车辆有较大威胁。同时威胁公路安全的还有由于岩石风化而形成的溜石坡。由于时序InSAR技术获取的地表信息有限，较难准确区分出滑坡、崩塌及溜石坡，这里将这些地质灾害统一称作不稳定斜坡。图7为中巴公路北侧西部一处不稳定斜坡的形变特征图，可以看到坡面上分布有大量负向形变量较大的PS点，累积最大负向形变达到-30 mm，且越靠近坡底形变量越小，在2016年6月之前未发生明显形变，之后大致呈线性趋势产生负向形变；坡底分布有少量正向形变PS点，累积正向形变达到20 mm，可能由碎石堆积引起，与不稳定斜坡的一般性特点相符。

图7 不稳定斜坡形变特征图Fig.7 Deformation characteristics diagram of unstable slope

3.2.2 冰川运动形变特征

由于盖孜河谷独特的地质地貌和气候条件，冰川运动及其引发的冰川融水泥石流是盖孜河谷区域的典型地质灾害之一，产生的形变相对较大，表现为年均形变量高值区域，地域上主要分布于各大沟谷中、国道、乡道公路两侧，时间上多集中于春夏融水强度较大的5—8月。研究区地势陡峭、冰雪融水丰富且具有丰富的松散固体源，冰川泥石流是盖孜河谷危害最为严重的地质灾害，但由于泥石流运动速度快，形变大，会引起严重失相干，在灾害点往往无法提取到有效的形变信息，但冰川泥石流会造成坡底碎石堆积，这在平均形变速率图上表现为正异常值集中区域。冰川前缘运动一般比较缓慢，可提取到较多有效PS点，故可根据冰川前缘冰碛物的形变特征合理推测冰川运动情况。由图6底图光学影像可看出，冰川在盖孜河谷南侧公格尔山及北侧盖昆山均广泛发育，图8为盖孜河谷北部一处冰川运动形变特征图，图8(a)为PS点分布及形变速率，图8(b)为矩形框和椭圆框内部PS点形变时间序列变化曲线，体现出冰渍物运动引发的地表形变，其中正负向累积形变分别达到+60 mm和-80 mm。由于该地区的土质多为季节性冻土，土质极不稳定，冰川运动可能会加剧不稳定斜坡的发生，甚至引发冰碛物滑坡。

图8 冰川运动形变特征图Fig.8 Deformation characteristics diagram of glacial movement

4 结语

本文以中巴经济走廊盖孜河谷段为研究区，利用PS-InSAR及SBAS-InSAR技术分别进行地表形变提取，比较两种时序InSAR技术所得结果发现，这两种方式获取的形变信息在形变速率的空间变化趋势上具有较好的一致性，其中SBAS-InSAR方法能更充分地利用小基线数据子集内的良好相干性，所获取的时序形变速率分布在空间上更为连续。根据SBAS-InSAR技术提取的结果可知研究区多处发生明显地表形变，主要分布于中巴公路两侧斜坡及布伦口以北的冰川前缘区域。结合光学影像识别出研究区内多处不稳定斜坡与冰川运动，并对其中典型实例进行了时序形变特征分析。不过，由于研究区历史上的实地勘测数据十分稀缺，所获取的形变结果尚难以进行完全的验证。本文所做的前期探索性研究初步展现了时序InSAR技术在高寒山区地质灾害监测中的应用潜力，有助于今后卫星数据与实地勘测同步开展条件下的进一步深入研究，并最终服务于盖孜河谷区域地质灾害易发性评价、基础工程设施建设与防护等工程应用。