Sentinel-1数据在西南山区水库变形斜坡InSAR监测中的适用性评价:以溪洛渡水库为例

李凌婧姚 鑫周振凯王德富

1.中国地质科学院地质力学研究所,北京 100081;

2.中国地质大学(北京)水资源与环境学院,北京 100083;

3.自然资源部活动构造与地质安全重点实验室,北京 100081;

4.中国地质调查局新构造与地壳稳定性研究中心,北京 100081;

5.中国地质大学(北京)工程技术学院,北京 100083;

6.自然资源部第三地理信息制图院,四川 成都 610000

0 引言

自1963年意大利瓦伊昂特大滑坡 (the Vaiont Landslide)发生以来,岸坡地质灾害一直是峡谷区水库关注的主要问题之一 (Hendron and Patton, 1987;郑嘉豪等,2020;李滨等,2020),库区斜坡活动性大小主要受地质条件、地震活动、降水及库水位波动等影响(中村浩之,1990;唐晓松等,2013;张淼等,2014)。近年来,合成孔径雷达(Synthetic Aperture Radar,SAR)技术在斜坡地质灾害形变识别、监测等方面有广泛应用(李晓恩等,2021),作为一种主动式微波传感器,具有全天候、全天时成像的特点,尤其适用于水库区这样的多云多雨地区(廖明生等,2012)。随着国内三峡水库区、金沙江四个梯级电站的相继建成及运行,InSAR (Interferometric Synthetic Aperture Radar,合成孔径雷达干涉)技术在水库区斜坡变形识别、结合库水位变化的时序变形监测等方面都有成功的应用(王桂杰等,2011;廖明生等,2012;Shi et al.,2016)。2001年美国阿拉斯加大学地球物理研究所与中国地震局地震研究所使用欧洲空间局的 ERS(European Remote Sensing Satellite)数据,利用D-InSAR(Differential Interferometric Synthetic Aperture Radar,合成孔径雷达差分干涉)技术监测三峡地区蓄水引起的地表形变(游新兆等,2001),展示了InSAR技术在国内水库建设中的应用前景。在此之后,Gao and Zeng(2007)利用2003~2005年11景万州地区ASAR数据进行时序InSAR处理,认为数据量和大气效应是限制InSAR数据在该地区应用的主要原因。随着高分辨率、短重访周期的SAR卫星相继升空,库区InSAR滑坡监测有了进一步发展。廖明生等(2012)将TerraSAR-X数据应用于三峡库区滑坡和泥石流等地质灾害的监测和预警工作中,取得了一定成效。

现有SAR系统有星载SAR(Kankaku et al., 2014)、机载SAR(Scott et al., 2016)和地基SAR(Herrera et al., 2009)3种,但由于机载SAR系统获取数据不稳定和地基SAR设备的高花费,使得星载SAR的应用更为普及。截至2021年,根据地质灾害InSAR监测技术指南附录C(中国地质灾害防治工程行业协会, 2018)中所示,可用的星载SAR卫星数据可达10种。随着2014年4月Sentinel卫星的发射,其超高的辐射分辨率、稳定的轨道系统、较大的覆盖能力、较短的重返时间、可免费下载的数据,使InSAR技术在斜坡灾害识别、监测、预警的应用上又迈进了一步(欧阳伦曦等,2017;徐帅等,2020)。但是,由于雷达卫星的重访时间、观测角度、分辨率、波长等因素的不同,对于地质条件复杂的山区滑坡识别,在数量和重合度上均有差异(刘星洪等,2018),因此,分析不同SAR数据在山区滑坡识别监测方面的适用性,具有重要意义。适用性是相对的,当变形斜坡所处的地质环境条件出现较大的改变时,彼处的适用可能就变成了此处的不适用,所以,适用性的分析是一个理论结合实际的过程。

文章以金沙江上游溪洛渡水库区变形斜坡(主要为滑坡、部分危岩体)为例,基于Sentinel-1、PALSAR-2和TerraSAR-X卫星SAR数据,采用D-InSAR、SBAS-InSAR方法提取变形斜坡信息,结合研究区的地质条件和库水位波动数据,评价Sentinel-1卫星SAR数据在山区水库变形斜坡InSAR监测中的适用性,以理论结合实际结果分析Sentinel-1数据是否可以在一定条件下替代其他两种商业数据,为相关的行业应用提供参考。

1 地质背景

1.1 地形地貌

溪洛渡水库位于四川省雷波县和云南省永善县接壤的金沙江峡谷段,正常蓄水位600 m,死水位540 m,库尾距坝里程约196 km。库区地处云贵高原和四川盆地两大地貌单元过渡地段,属强侵蚀高山、中山地貌类型,河谷狭窄,岸坡坡度30°~60°,最大可达80°。地势西高东低,山脉走向与构造线大体一致,以近南北及北东向为主,主要活动断裂有峨边-金阳断裂、金阳断裂及五莲峰断裂(图1a、1b)。

Sentinel-1搭载C波段SAR传感器,具有4种工作模式,其中,干涉宽幅模式 (Interferometric Wide Swath,IW)数据可用于差分干涉,数据幅宽250 km,在溪洛渡库区需南、北2景拼接处理,其他星载SAR数据均需要2~3景才可全面覆盖。库区整体走向大致与卫星升轨方位向垂直,与降轨方位向平行,单一轨道无法满足两岸的识别。

1.2 灾害发育情况

在地层上,除石炭系、泥盆系、古近系及新近系外,库区地层从元古界到第四系均有出露,其中碎屑岩及碳酸盐岩地层出露较广,坝区、近坝段出露的主要地层为二叠系的阳新灰岩及峨眉山玄武岩,河谷、支沟两侧的山间谷地内零星分布第四系的冲洪积、崩坡积等(邓宏艳和王成华,2011)。蓄水前灾害主要发育在这些碎屑岩及碳酸盐岩中,大型古滑坡堆积体变形则受软岩控制(高改萍和杨建宏,2003)。2013年5月水库蓄水后,陈廷照和戴福初(2018)利用光学遥感数据圈定了库区主河道及支流范围内共计213处滑坡点,其中特大型及大型滑坡占比74%;在涉水滑坡中,按照滑坡类型,基岩顺层滑坡和堆积层顺基覆界面滑坡占比最大。库区大型及以上灾害较为发育的特点,使得大多数星载SAR数据分辨率均满足一定的探测条件。

此次研究所用SAR数据范围如图1c所示(仅标识原始数据拼接裁剪后的实际处理范围)。

图1 研究区地理位置、地质背景及SAR数据范围Fig.1 The location and the geological background of the study area, and the SAR data ranges(a) The location of the study area; (b) The elevation and active faults; (c) The SAR data types, and InSAR processing ranges after images being spliced and clipped

2 数据与方法

2.1 基础数据

针对覆盖研究区的存档SAR数据,选取了3种不同波长的SAR数据源,分别是Sentinel-1(C波段)、ALOS PALSAR-2(L波段)及TerraSAR-X(X波段),详细参数见表1。此次InSAR处理,Sentinel-1多视数为9×1,多视后的分辨率为距离向×方位向=10.48 m×13.99 m;PALSAR-2多视数为2×2,分辨率为2.86 m×4.43 m;TerraSAR-X多视数为3×3,分辨率为2.73 m×5.89 m。

表1 星载SAR传感器基本参数及特征表Table 1 Basic parameters and characteristics of Spaceborne SAR sensors

利用InSAR方法计算的形变结果是SAR影像坐标系,需要投影到地理坐标系统才利于分析,一般需利用外部数字高程模型 (Digital elevation model,DEM)数据配合投影转换(刘国祥等,2019)。文章利用的DEM数据分辨率为6m,为SPOT-7卫星立体像对制作,其分辨率优于Sentinel-1,与PALSAR-2和TerraSAR-X基本持平,在重采样过程中可基本还原位置。

Sentinel系列是欧洲全球环境与安全监测系统项目——“哥白尼计划”的成员,其于2014年4月3日发射首颗卫星Sentinel-1A,每12天覆盖全球一次,2016年4月25日Sentinel-1B发射成功,双星重访周期缩短至6天(欧阳伦曦等,2017)。PALSAR-2重访周期为14天,TerraSAR-X为11天。除Sentinel-1为免费数据外,其他两种SAR数据均为商业数据,需编程采购,因卫星执行任务、经费等原因,在数据连续性上较差,因此,从数据选取方面,Sentinel-1性价比更高。

2.2 分析流程

以往分析主要是根据SAR卫星本身的参数(王志勇和张金芝,2013;Hu et al.,2014)或技术本身的误差 (大气、植被等;Xu et al.,2011;廖明生和王腾,2014;占文俊等,2015)进行的,此次研究基于影响地质灾害InSAR识别的相关因素进行分析,分析流程图见图2。采用从整体到局部的评价方法,区域识别上主要是时间(相干性、重访周期、数据量情况)和空间(理论及实际的有效识别结果),单体识别上主要为灾害单体特征,如范围、形变和准确度。

图2 分析流程图Fig.2 Analysis flow chart

3 区域性水库变形斜坡识别适用性评价

3.1 时序识别应用

与其他的在轨观测雷达卫星相比,Sentinel卫星最大的特点是具有优良的覆盖性能和重访性能,适合与库水位变化相关性极强的水库滑坡的监测工作。为了开展评价分析,在研究区编程购买了PALSAR-2及 TerraSAR-X 数据,但是,由于PALSAR-2卫星经常执行突发任务,使得编程拍摄一度中断,造成数据不连续,而TerraSAR-X数据由于波长较短,在植被较茂盛的夏季干涉效果很差,性价比低。因此,在实际应用过程中,如不考虑编程SAR数据,在只用存档数据的前提下,Sentinel是最优选择。

基于Sentinel丰富且连续的存档SAR数据(2014年10月至2021年2月),使得区域性的变形斜坡时序识别成为可能,此次在研究区成功解译的变形斜坡约200处,类型有滑坡、危岩体和塌岸3种。

溪洛渡水库水位在每个水文年内的6月份至10月末为水位上升阶段,于次年6月再降至540 m左右(陈廷照和戴福初,2018),仅针对大的蓄—排水周期,不考虑水位的小范围波动,结合Sentinel-1存档数据情况,将SAR数据分为6个完整的水位上升(下降)周期,每个阶段的SAR数据数量在图3和图4的括号中标示。2020年10月至开展此文工作时未满足一个水位下降期。

图4 水位上升(下降) 周期内降轨Sentinel-1数据存档情况(括号内数字为存档SAR数据数量)Fig.4 Archive of the Sentinel-1 data from the descending orbit in periods of rising (falling) water level (The numbers in brackets are the amount of archived SAR data)

由于季节、植被等因素的影响,相干性无法保证每一个像对D-InSAR结果都具有可识别性。在溪洛渡水库区,每年6月至11月初的水位上升期,周边植被发育较好,数据相干性较差。因此,为了最大化利用SAR数据,并获取水位上升(下降)周期的变形斜坡分布情况,将同期内D-InSAR数据做并集,提取水位上升(下降)周期内变形斜坡形变最大范围,以揭示整个周期内灾害发育情况。随着双星的运行,2017年后数据量逐渐增多,即2016年“水位下降3”之前的数据较少(图3),会引入较大的误差,因此只针对该期以后的数据进行分析。

图3 水位上升(下降) 周期内升轨Sentinel-1数据存档情况(括号内数字为存档SAR数据数量)Fig.3 Archiving of the Sentinel-1 data from the ascending orbit in periods of rising (falling) water level (The numbers in brackets are the amount of archived SAR data)

由于人工圈定范围的不确定性,按水位周期解释时采用 “点”的方式,分为 “变形”或 “非变形”点。此次研究规定,边坡表面发生变化的面积小于解译斜坡最大边界的10%时,这种变化不足以反映边坡的整体变形,边坡被视为不变形。由于解释方法的局限性,不能排除人工解译中的漏识。研究区时序解译滑坡点位解译标准示意图如图5所示。图5a所示InSAR结果为边坡的最大变形范围如图5a所示的InSAR变形区,红色线为最大变形边界线。从图5b—5h可以看出,后续边坡的变形量均没有图5a所示大,变形减小。在图5f和5h中,红色形变像元占比为3%和4%,不满足此次设定的10%,因此为该滑坡的未变形时期。

图5 研究区变形斜坡不同时段干涉相位图及斜坡变形情况Fig.5 Interferometric phase of the deformed slope at different periods in the study area and the deformation status

综上,根据利用Sentinel-1数据解译了不同蓄—排水周期内变形斜坡的发育变化情况(图6),可以得到:①水位下降期变形斜坡数量较水位上升期更多;②2017~2020年4年间变形斜坡的整体数量呈递减趋势;③2020年10月至2021年2月不完整水位下降期内数据解译情况显示,变形斜坡数量小幅度回升。

图6 变形岸坡数量时序解译情况示意图Fig.6 The annual interpretations of deformed slopes

3.2 空间识别可行性

由于雷达成像是根据微波回波返回时间长短记录的,因此为侧视系统,而星载SAR一般为右侧视。现有SAR卫星为近极地轨道,绕地球南北飞行,由于地球自转,卫星在同一地表位置的上空会有两种方式经过,即升轨(由南向北飞行)和降轨(由北向南飞行)。由于雷达为侧视观测,当地面对着微波方向的坡度过大时,回波会反射到同一成像单元,相应的SAR图像会出现叠掩现象;当地面背着微波方向的坡度过大时,较高地区会对背波的区域形成遮挡,雷达接收不到地表反射,SAR图像会出现阴影现象(任云,2013;张芳,2018)。阴影叠掩范围的影响因素主要为卫星入射角与地形之间的关系(张同同等,2019)。

溪洛渡水库区为高山峡谷区,叠掩、阴影现象更为严重。选取3种SAR数据升轨、降轨的交集区域,计算Sentinel-1升降轨、PALSAR-2升降轨及TerraSAR-X降轨数据在该地区的有效观测情况,计算结果见图7。雷达入射角为:Sentinel-1升轨35.26°,Sentinel-1降轨37.75°,PALSAR-2升轨39.66°,PALSAR-2降轨38.74°,TerraSAR-X降轨26.58°(叠掩最严重)。它们的有效观测比例如图7a—7e左下黑字所示,分别为:70.3%、68.9%、70.4%、67.6%和52.5%。在不考虑分辨率的情况下,仅对比有效观测面积,通过对雷达常用观测模式下的SAR数据比较可以看出,在库区Sentinel-1数据(图7a、7b)与其他两种SAR数据观测能力相比持平(图7c、7d)或更优秀(图7e)。

图7 研究区不同SAR卫星在SAR数据交集区域的有效观测比例Fig.7 Proportion of effective observations by different SAR satellites in the intersection area of SAR data(a) The effective observation area of the Sentinel-1 ascending image; (b) The effective observation area of the Sentinel-1 descending image; (c) The effective observation area of the PALSAR-2 ascending image; (d) The effective observation area of the PALSAR-2 descending image; (e) The effective observation area of the TerraSAR-X descending image; (f) The hill shade

对于实际解译情况,选取3种降轨数据,范围为它们的交集,分别进行变形斜坡解译,解译情况如图8所示,Sentinel-1解译58个点,PALSAR-2解译48个点,TerraSAR-X解译37个点。由于数据获取时间、分辨率、波长、入射角等因素影响,解译的位置不能完全一一对应,所识别的岸坡变形情况也不同,因此,不能片面的认为Sentinel-1更准确。但分辨率较低的Sentinel-1之所以能解译更多的点,主要原因在于其12天重访高频观测的特点,较好的连续性可提高干涉像对的相干性。但是,对于面积较小的变形斜坡,则需要高分辨率的SAR数据来补充。因此,如需全面获得库区的岸坡形变位置,需要多种SAR数据联合解译。

图8 三种SAR卫星降轨数据在同一地区的变形斜坡解译情况示意图Fig.8 Schematic diagram of deformed slope interpretations of three satellite SAR data on descending orbits in the same area

4 单体滑坡识别适用性评价

4.1 范围

此次处理的Sentinel-1数据分辨率为距离向×方位向=10.48 m×13.99 m,在变形范围识别上,变形斜坡大小需在地表形变空间尺度超过一定数量的分辨率单元。但由于地质条件复杂,变形干扰因素较多,解译过程中还需要针对变形地质体特征进行人机交互与现场校核。结合现场踏勘结果,研究区3种SAR数据解译的最小变形斜坡投影面积约为1300 m2,其中,Sentinel-1数据解译的最小变形斜坡投影面积约2400 m2,是约35 m(长)×77 m(宽)大小的圈椅状斜坡单元,共有16个变形像元聚集。

4.2 形变

水库运行期间,岸坡破坏表现为缓动型和剧动型两类(陈洪凯和周晓涵,2016)。在变形斜坡InSAR监测应用上,由于技术本身的特点,主要是针对缓动型(李晓恩等,2021)。这里的缓动型变形斜坡,指的是滑坡失稳前,其加速趋势可以被大于2期的SAR数据捕捉到 (Dong et al., 2018;Intrieri et al., 2018)。理论上,当差分干涉图中一个分辨率单元内发生沿雷达视线向超过二分之一波长的形变时,D-InSAR就无法监测到这种形变(Massonnet and Feigl, 1998)。进一步,数据处理时的多视数和相干性也对监测效果影响显著(蒋弥等,2009)。因此,当变形量过大时,应选取雷达波长更长的SAR数据,或以缩短重访周期、增加空间分辨率等方式改进。根据搜集的GNSS(Global Navigation Satellite System,全球导航卫星系统)资料,溪洛渡水库区部分滑坡形变最大可达1 m/月,且不会溃决入库,如此大的变形量,使得传统InSAR技术无法准确测量其实际变形量,只能监测其相对变形趋势。以上问题,也是InSAR技术在滑坡监测与预警上一直无法大面积应用的原因。

4.3 典型滑坡分析

在对整个库区的变形斜坡识别基础上,选取升轨、降轨数据均可识别的翌子村滑坡(滑向西—西南)进行单体适用性评价分析。全景照片见图9a,红色虚线为InSAR识别的实际变形区。

4.3.1 基本情况

翌子村滑坡处于金沙江右岸,该滑坡位于巨型古滑坡堆积体前缘,为一个次级堆积层滑坡,水库蓄水前已有变形 (Li et al., 2019a)。坡面形态较复杂,陡坎发育,总体上坡体后缘稍陡,中部相对平缓。位于该区段的道路开裂、房屋变形(图9b)。滑坡体中已形成连贯的滑动面,因此该滑坡常年变形。

图9 库区翌子村滑坡现场照片Fig.9 Field photos of the Yizicun landslide(a) Panoramic view of the landslide (towards NW); (b) Photo of a abandoned house

4.3.2 多源SAR数据D-InSAR结果

用Sentinel-1升轨、降轨及PALSAR-2升轨、降轨相似时间段的D-InSAR处理结果,和2017年相似时间段 TerraSAR-X降轨 D-InSAR结果(TerraSAR-X升轨在同期未编程购买数据)进行对比分析。为使图像获得最佳细节展示,采用0~π周期进行出图,图例已根据不同数据的波长进行实际变形大小的转换,见图10。Sentinel-1数据升降轨D-InSAR结果见图10a、10b,PALSAR-2数据升降轨D-InSAR结果见图10c、10d,TerraSAR-X数据降轨D-InSAR结果见图10e。光学遥感影像显示的滑坡边界,与InSAR识别的实际变形区有较大差别(图10f)。在相似时间段,分辨率高的SAR数据(PALSAR-2约4 m,TerraSAR-X约6 m;表1)D-InSAR结果显示出了形变的详细细节,如图10c、10e识别出了翌子村滑体中部的三段滑坡台阶 (Li et al., 2019a),结合现场调查,与InSAR识别情况相符。

4.3.3 多源SAR数据PS-InSAR结果

由于翌子村滑坡没有GNSS数据,通过对比Sentinel-1降轨、PALSAR-2降轨的PS-InSAR时序曲线,来分析Sentinel在该滑坡的监测适用性。方法是通过相同位置一定区域内(15 m直径)变形点差值,分析得到两种卫星识别的变形规律是否一致。由于变形较大的前部位置Sentinel-1时序结果数据缺失,因此选取对比区域位置在后缘(位置见图10f)。变形曲线见图11,图中十字点为PALSAR-2数据的PS-InSAR时序曲线,实心圆点为Sentinel-1时序曲线,他们的趋势线由2周期移动平均方法合成。PALSAR-2趋势线上标注的A、B、C三个线段,数据连续性较差(即无SAR数据)。

图10 库区翌子村滑坡多源SAR数据D-InSAR结果对比Fig.10 Comparison of D-InSAR results based on the multi-source SAR data of the Yizicun landslide in the reservoir(a) The interference figure of the Sentinel-1 ascending image; (b) The interference figure of the Sentinel-1 descending image; (c) The interference figure of the PALSAR-2 ascending image; (d) The interference figure of the PALSAR-2 descending image; (e) The interference figure of the TerraSAR-X descending image; (f) The optical image and the deformation boundary

(1)监测面积。翌子村滑坡地面农耕植被较多,Sentinel-1时序监测结果在坡体大变形部位(中部、前部)结果丢失。PALSAR-2由于其波长较长,更适合植被密集区的大面积的时序变形的监测,在翌子村滑坡大范围监测上的适用性更高。

(2)趋势识别。如图 11所示,由于PALSAR-2在A、B、C三个区域,数据连续性较差(无数据),因此该时间段内的变形趋势无法有效识别。与此同时,由于长时间的数据缺失,会造成最近一对SAR数据间的变形超过其InSAR最大量程,解缠时丢失相位周期。因此,从图11两种SAR数据的变形趋势上,每当PALSAR-2长时间缺失数据时,后期两种变形结果的差值就越大。而Sentinel-1由于连续性较好,变形趋势较为连续。因此Sentinel-1更适合连续小变形的趋势识别。

图11 库区翌子村滑坡多源SAR数据PS-InSAR时序变形曲线结果对比Fig.11 Comparison of PS-InSAR time-series deformation results based on the multi-source SAR data of the Yizicun landslide in the reservoir

(3)形变测量。该时序分析的点位区域靠近滑坡后缘,水位对其的影响减弱,降雨对其的影响增强。对于变形较大的区域,Sentinel-1数据也会有相位丢失的现象,造成实际的监测值较小。同时,SAR数据的分辨率会影响实际形变的测量,多视处理后的一个像元内的形变会被平均,一些高值点会被弱化。像素偏移追踪测量的方法(Pixel Offset Tracking,POT)对于分米至米级的大变形的监测较为合适 (Strozzi et al., 2002; Singleton et al., 2014),但监测精度会有所降低(Shi et al., 2015; Li et al., 2019b; Jia et al., 2020)。

5 结论

造成InSAR测量结果与实际变形之间的差异的主要原因,主要有InSAR技术的视线向测量特点、波长对可测最大变形量的限制、分辨率控制的像元区域内形变被平均、夏季植被造成的失相干、SAR数据时序不连续造成的失相干等。文章以金沙江上游溪洛渡水库区为研究区,以上述影响因素为切入点,基于Sentinel-1、PALSAR-2、TerraSAR-X数据,采用D-InSAR、PS-InSAR方法,分析评价Sentinel-1数据在西南山区水库变形斜坡InSAR监测中的适用性,结论如下。

(1)Sentinel-1数据在研究区可解译的变形斜坡约200处,类型有滑坡、危岩体和塌岸;经现场核查,Sentinel-1数据解译的最小变形斜坡体投影面积约为2400 m2,为一个约35 m(长) ×77 m(宽)大小,16个变形像元聚集的斜坡单元。

(2)高山峡谷区叠掩、阴影现象严重,通过对常用雷达观测模式下的SAR数据的比较,在SAR数据交集区域,有效观测面积为Sentinel-1升轨70.3%,Sentinel-1降轨68.9%,PALSAR-2升轨70.4%,PALSAR-2降轨67.6%,TerraSAR-X降轨52.5%,在不考虑分辨率的情况下,在溪洛渡水库区Sentinel-1数据与其他2种SAR数据相比观测能力持平或更优秀。

(3)溪洛渡水库区在6月至11月初的水位上升期,周边植被发育较好,数据相干性也较差,2017年后Sentinel-1A(1B)双星拍摄获取SAR数据量增加,高频观测使相干性提高,利用2017年后该卫星数据可有效识别水库蓄—排水周期内的区域性变形斜坡发育变化情况。

(4)当长时间缺失SAR数据时,会造成最近一对SAR数据间的变形超过其InSAR最大量程,解缠时丢失相位周期,Sentinel-1由于连续性较好,变形趋势较为连续,因此Sentinel-1更适合连续小变形的趋势识别。

致谢:感谢评审专家在文章内容、结构及结论组织上提出的宝贵意见。感谢编辑的审阅,使本文在语言表述方面有了较大的改进。