在中国,山地和丘陵地区占据国土面积的65%,而滑坡是山地丘陵地区最常见的自然地质灾害。滑坡主要是由活动构造和历史地震影响的内部动力地质作用和降雨等外部自然因素两方面控制产生的^[1],滑坡一旦形成,就会对居住在山区的人们有着致命的影响,因此提前预警滑坡的形成对防止危险事故的发生至关重要,而斜坡的运动特征是滑坡最重要的指标之一,是滑坡发生的先兆。常规滑坡灾害点监测(如水准和全球定位系统(global positioning system,GPS)测量)要耗费大量的人力物力,十几个到上百个观测点的布置需要耗费大量的时间,地理条件偏远和自然环境恶劣的山区等地往往无法开展专门的地面观测。合成孔径雷达差分干涉测量技术(differential interferometric synthetic aperture Radar,DInSAR)的出现有效改善了这一难题。在滑坡灾害的相关研究中,范景辉^[2]利用合成孔径雷达(synthetic aperture Radar,SAR)数据、反射体及GPS测量对范家坪滑坡变形进行监测,综合分析滑坡影响因素; Liu等^[3]和王振林等^[4]用ALOS-2数据对中国西南部的滑坡变形进行监测,取得了不错的成果; 王振林等^[4]和Intrieri等^[5]利用Sentinel-1数据对滑坡变形进行时序性分析监测。

四川省茂县地区位于西南山区,地质灾害频发,滑坡灾害造成巨大人员伤亡。本文综合过往研究经验,获取了2015—2017年长波长L波段的ALOS-2和短周期观测的Sentinel-1两种不同分辨率的遥感卫星数据,并结合资源三号(ZY-3)卫星制作的高分辨率5 m数字高程模型(digital elevation model,DEM)数据和航天飞机雷达地形测绘(shuttle Radar topography mission,SRTM)30 m分辨率DEM数据,使用两轨法和Stacking方法探测茂县中部地区的滑坡形变,分析茂县滑坡的形变特征,为地区的防灾减灾提供相关参考。

四川省茂县位于四川省阿坝藏族羌族自治州,是中国西南地区典型的山区县(图1)。县域内跨越岷江和涪江上游高山河谷地带,发育着茂汶断层和龙门山断层,断裂顺着岷江河谷发育,存在着明显的破碎带和脆弱的地质条件。当地的年降雨量集中于5—10月,分布不均,最大日降雨量达到75.2 mm,瞬时降雨量大,很容易导致地质灾害的发生^[6]。地区内以千枚岩、板岩的变质岩为主,而在岷江流域的河谷冲刷沿岸坡上堆积了大量的松散土块,基础设施的防治条件较差,极易造成土块崩塌。沿岸半山腰一旦崩塌,地形较陡,落差较大,很容易发展成巨型、岩质或基岩滑坡^[1]。地区处于南北方向地震带即松潘和龙门山地震带的中部地区,同时九顶山华夏系构造带从南向北斜贯该地区,受构造因素运动影响较剧烈。地区内地震现象常有发生,1933年叠溪7.5级地震和2008年的汶川8.0级大地震更是减弱了茂县地区内陡坡的稳定性,在暴雨等外部诱发因素下,区内的滑坡运动极有可能加剧,需要特别防范。2017年6月,茂县叠溪镇新磨村突发山体高位垮塌,造成河道堵塞,100余人被掩埋。当地山高坡陡,位于流域上游而建设的公路、水电、水利设施等人工建设项目,更是有着巨大的隐患,而对山区游客、当地居民和通行车辆存在着潜在的人身安全威胁^[7,8,9]。在自然和人为的综合作用下,该区域为重要地质灾害隐患地区,区内的白布村滑坡群便是典型案例。研究区地形和数据覆盖情况如图1所示,图中红色框表示具体研究区域,水蓝色框表示Sentinel-1数据挑选的burst覆盖范围,黄色和深蓝色框分别代表ALOS-2 SM2级别和SM3级别数据覆盖范围。

ALOS-2卫星是由日本航空航天局(Japan Aerospace Exploratin Agency, JAXA)于2014年5月发射的高级陆地观测卫星。其上搭载PALSAR-2传感器,工作波段为L波段(1.2 GHz段),最高分辨率为: 3 m(距离向)×1 m(方位向),能在复杂大气条件下全天候工作,最大观测范围为490 km×350 km。相对于短波波段的SAR,ALOS-2能更好地去除植被在斜坡上的去相关效应,更加准确地对地表信息进行观察,及时监测到一些自然灾害,例如滑坡。本文使用的ALOS-2数据参数如表1所示。因为本次实验采集到的ALOS-2数据量较少,无法使用多时相InSAR技术,便采用常规DInSAR技术中的两轨法探测地表形变。

Sentinel-1卫星由欧空局(European Space Agency, ESA)运营,搭载了C波段(5.6 cm)传感器,获取数据时间上具有稳定的连续性,其IW成像模式,幅宽250 km,地面空间分辨率5 m(距离向)×20 m(方位向),采用新型循序扫描地形观测(terrain observation with progression scans,TOPS)成像技术,数据量丰富,下载方便。研究采集了45景Sentinel-1升轨数据,日期为2015年11月26日—2017年12月21日,采用Stacking技术探测地表形变。Sentinel-1数据参数如表2所示。在本次研究中,DEM采用了两种不同分辨率的数据。一个为基于我国的ZY-3的两景立体像对通过ENVI和PCI软件生成的5 m DEM数据。ZY-3卫星配备三线阵测绘相机和多光谱相机,全色波段分辨率2.1 m,幅宽52 km^[10],具有广泛的覆盖范围,受大气和地形的影响较小,适合进行偏远地区的测图工作。在本次研究中用ZY-3卫星制作的5 m分辨率的DEM辅助ALOS-2数据实施图像配准和地理校正^[11],在InSAR差分干涉的过程中,高精度的高程数据能够提高配准精度、削弱高程残余相位,5 m的数字表面模型能更好地满足研究的需要。另一个的DEM地形数据是由美国国家地理空间情报局(National Geospatial-Intelligence Agency,NGA)和美国国家航空航天局(National Aeronautics and Space Administration,NASA)免费提供的30 m分辨率的SRTM-DEM数据,被用于辅助更大范围内Sentinel-1 数据的图像配准和地理校正。

DInSAR两轨法是对同一地区的两景SAR影像进行干涉差分。它需要一个外部DEM来模拟和消除地形的影响,外部DEM被转为地形相位后,地形因素被消除,从而得到形变相位。两轨法首先需要对外部DEM 和主影像进行精确匹配,把DEM 的高程模拟成地形干涉条纹形式(即模拟地形相位),然后就可以从干涉图中移除DEM 模拟的地形相位,经过滤波及相位解缠后,从而得到形变信息。

Stacking技术对生成的多景差分干涉图进行解缠后估算出其线性相位速率^[12],可以在较少的时序数据基础上获取年平均形变速率,一定程度上减少大气误差、提高形变精度,从而完成滑坡的探测。其原理是基于最小二乘法对N组观测值线性回归,估算公式为:

式中: Δt _j为第j组干涉图的时间基线; φ _j为第j组解缠的差分干涉相位; ph_rate为线性相位速率;N为生成的差分干涉图的数量。

InSAR数据处理的具体步骤如图2所示。第一步,ALOS-2数据配准是将其与高分辨率的DEM进行地理编码后生成初始查找表,计算轨道偏移量生成多项式精化配准查找表,多次重采样后进行亚像元级配准; Sentinel-1数据配准是将采集的多个SAR数据选出覆盖研究区范围的条带,按照研究区范围进行拼接、裁剪,随后选择一张图像作为主图像,利用DEM地形数据和精密的轨道参数对SAR数据进行初步配准,再建立主辅影像配准的查找表进行精确配准,最后通过强度互相关最大化的方法估算残余偏移量,精确配准后再基于频谱多样性进行重叠区配准。第二步,分别将配准后的两种数据进行影像多视处理,ALOS-2的多视比为4∶2,Sentinel-1的多视比为5∶1。第三步,将Sentinel-1数据通过限定时间基线和空间基线的原则生成67对干涉对,进行像对差分干涉处理(空间基线限定150 m,时间基线限定24 d); 选择相同级别的ALOS-2数据两两组合生成干涉对。第四步,将ALOS-2数据干涉对进行自适应滤波,以及去除大气延迟相位,选择相干图中相关系数较高、位于稳定区域的像元作为参考点以进行相位解缠,便得到了研究区的形变相位。第五步,针对Sentinel-1数据将多个差分条纹图叠加求平均后进行滤波,并利用最小费用流算法进行相位解缠,然后对多景解缠的差分干涉图根据式(1)计算像元的线性形变速率。相比于DInSAR技术,Stacking方法可以有效地抑制大气效应与DEM误差的影响。

ALOS-2数据和Sentinel-1数据的干涉对生成的差分干涉图如图3所示。从图3可以看出,ALOS-2的L波段比Sentinel-1的C波段,干涉相位受地表植被影响更小,克服时间去相干的能力更佳,干涉相位信息更加丰富。

InSAR形变结果通常采用水准仪或GPS观测数据进行验证。研究区位于陡峭的山区,地理条件恶劣,山高险峻,无法提供精准的地面测量数据来定量验证InSAR形变结果。本研究将采用光学影像特征和已有历史调查资料作为辅助验证数据,对ALOS-2与Sentinel-1两种数据的InSAR探测结果进行比较和交叉验证,以间接验证InSAR探测滑坡形变的可靠性。两种数据集探测的形变结果如图4所示。图4(a)显示了基于ALOS-2数据集获取的雷达视线向(LOS)年均地表形变信息,形变速率结果集中于-80~-10 mm/a。如图4(b)所示,基于Sentinel-1A数据集提取出的点目标分布较均匀密集,形变速率结果分布集中于-90~-10 mm/a。对比两种数据集探测的形变结果可以得到初步结论,研究区出现了多个滑坡变形中心,主要分布于研究区的中部,山谷的两侧。形变方向主要以两岸向河流滑动为主,形变速率空间分布不均一,从边缘到中心逐渐增加,呈现圆形或椭圆形。由于本次实验采用的是升轨SAR数据在河谷西侧形变信息比东侧更为丰富,主要原因是河谷东侧斜坡受压缩、叠掩等几何畸变影响导致可靠结果难以获取^[13]。综合分析两种数据集的形变速率结果,圈定了8个活动滑坡,分别为黄草坪滑坡a,白布村滑坡b和c,窝窝店滑坡d和e,苏家坪滑坡f,水草沟滑坡g,立角村滑坡h。两种数据的探测结果在其中7个活动滑坡处的空间分布能互相对应,分别是a,c,d,e,f,g和h区域。从图4(b)可以看出,a和b区域的形变量较大,形变速率最大绝对量达到了200 mm/a。综合以上可以得出,采用少量长波段ALOS-2影像和DInSAR两轨法,能够在具有一定植被覆盖度的山区探测到较为明显的滑坡地表形变; 而基于长时序Sentinel-1影像和Stacking时序分析方法,探测效果比采用ALOS-2影像的两轨法更佳。

对比两种数据得到的探测结果,可以看出圈定的活动滑坡空间分布能互相对应,侧面验证了本实验的结果准确性。不过因为不同传感器、入射角、以及时间、空间失相关等影响,两个数据集得到的形变结果在形变值有略微差异。ALOS-2数据结果未探测到a和b活动滑坡区的原因可能是因为数据集具有较长的空间基线,区域内年变化较大,数据相干性差,导致区域内形变点太少甚至没有点。c,d,e,g,h等区域Sentinel-1A的形变结果值分布较小,这可能是由于采集45景2015年11月26日—2017年12月21日之间的线性年平均形变速率小于ALOS-2数据采集时段的线性形变速率。本研究对圈定的a,b,c,d,e,g等6处进行了野外地质调查。调查发现,a,b,c,d,e等遥感识别的活动滑坡存在明显的地表形变,滑坡堆积体、山坡小范围坍塌、公路裂缝、陡坎等滑坡形变特征在多处出现,对于一些地势起伏较大的区域g,野外调查时只能远眺,发现河沟对岸山坡小范围滑塌,有大块板状碎石。比如,从图5(a)可以看出在圈定活动滑坡区域黄草坪a的右边缘存在拉裂缝、堆积物,从图5(b)可以看出在圈定的活动滑坡区域窝窝店滑坡e前缘,存在不断滑动的迹象,右侧拉裂缝明显,有潜在危险,整体判断为古滑坡复活。

滑坡运动是滑坡灾害发生的先兆,为防治灾害,探测早期的滑坡形变分布至关重要。而不同的滑坡发育阶段也有着不同的形变特征,形变速度也各有不同,所以就需要运用多种数据进行组合探测。

1)本文利用高分辨率的ALOS-2数据和5 m分辨率的DEM数据探测茂县白布村的滑坡形变,高分辨率的DEM能够辅助实现InSAR影像的高精度配准,并削弱高程残余误差,为应用DInSAR方法获得更为可靠的形变速率图提供了有利条件。而重访周期较短的Sentinel-1数据基于Stacking技术生成的形变速率图与ALOS-2数据表现了空间分布上的一致性。野外调查显示,在依据两种InSAR结果圈定的典型活动滑坡上,存在实地变形特征。

2)结果表明,数据集在时间基线不同的情况下,其探测能力和优势是不同的。在保持一定相干性的基础上,时间基线越长,越有利于识别出具有较小形变速率的活动滑坡; 时间基线越短,则更有利于识别具有较大形变速率的活动滑坡。长、短时间基线数据结合,能对滑坡的慢速运动和快速运动实现组合识别。

3)结果显示,采用少量长波段ALOS-2影像和DInSAR两轨法,能够在具有一定植被覆盖度的山区探测到滑坡地表形变; 当应用C波段的Sentinel-1影像时,通过积累连续的多景数据,并应用时序分析方法,探测效果比采用ALOS-2影像的两轨法效果更佳。

志谢: 衷心感谢JAXA EO-RA2项目(编号: P3073002)提供ALOS-2数据,ESA提供Sentinel-1数据。