焦志平,, 江利明,, 牛富俊 ,, 郭 瑞,, 周志伟
(1.中国科学院精密测量科学与技术创新研究院大地测量与地球动力学国家重点实验室,湖北武汉 430077; 2.中国科学院大学,北京 100049; 3.中国科学院西北生态环境资源研究院冻土工程国家重点实验室,甘肃兰州 730000)
冻融作用是冻土区滑坡发育的重要外动力因素之一[1],冻融滑坡形态与常规降雨型滑坡相似,但影响因素和失稳机制更为复杂。位于川藏交界处的西藏自治区东部,季节性和多年冻土分布广泛、地质环境复杂、地形和气象条件恶劣,是我国滑坡地质灾害最严重的地区之一。近年来,随着青藏高原暖湿化加剧,该地区冻融滑坡灾害呈现多发频发的趋势,对人民生命财产安全造成严重威胁,制约了当地经济社会发展[2]。因此,加强冻土区滑坡灾害隐患的监测识别研究,对藏东地区冻融地质灾害早期识别及风险评估具有重要意义。
目前,滑坡灾害监测主要采用GNSS位移监测、水准测量、裂缝计、无人机与卫星遥感等实地测量与遥感解译相结合的方式[3]。但是,藏东地区滑坡灾害多位于高海拔、高寒、高陡等艰险山区,人员难以抵达,实地测量通常要耗费大量人力财力物力,监测范围及密度有限,且难以获取滑坡历史变化信息[4]。光学遥感手段以其大范围、高效率的优势在地质灾害调查识别中发挥了重要作用[5],但数据采取受云雾遮蔽影响、人工解译不确定性大,在滑坡形变监测中存在错检漏检、难以进行定量分析等问题[6]。
作为一种快速发展的空间大地测量技术,合成孔径雷达干涉测量(InSAR)具有全天候、高分辨率、高精度、大范围等特点[7],近年来在滑坡[8]、地震[9]、冻土[10]等自然地表形变监测中应用广泛。特别是随着时序InSAR 技术的发展,国内外学者在四川、云南、甘肃、三峡等地开展了大量滑坡监测研究,取得一系列重要成果[11-15]。但是,对于高山峡谷复杂地形条件下的滑坡监测而言,雷达侧视成像更易受叠掩、透视收缩、阴影等几何畸变的影响,单一轨道方向(升轨或降轨)InSAR在地形起伏地区导致数据空洞[16],而且对滑坡形变解译提出了挑战[17]。欧洲航空局(ESA)分别于2014年和2016年发射的Sentinel-1A/B SAR 卫星,具有短重访周期(6/12/24 天)、高分辨、大幅宽、多角度数据获取优势[18],能够获取更丰富的冻融变形过程信息,在高山峡谷区冻融滑坡隐患早期识别和风险评估方面具有很好的应用前景。
本文以藏东地区317 国道矮拉山段为例,利用小基线集(Small Baseline Subset,SBAS)时序In-SAR 技术,分别对2017 年3 月—2019 年7 月期间Sentinel-1A55 景升轨、69 景降轨SAR 数据进行分析,获取该区域高时空分辨率的地表形变速率和滑坡形变时间序列,并结合气象资料和现场调查分析了典型滑坡体的形变特征及影响因素。
研究区位于西藏自治区昌都市江达县岗托镇矮拉山地区(98.43°~98.58° E,31.6°~31.68° N),是经由317国道进入西藏地区的重要通行路段。平均海拔4 710 m,地势由东南向西北方向倾斜[19],属于青藏高原中高山剥蚀地貌,区域地震活动构造背景较强,崩塌、滑坡等自然灾害严重。地表分布2条河流,主要由大气降水和上游冰雪融化水补给[20]。该区域年均温差较大,太阳辐射强烈,气温范围为-15~28 ℃,平均气温4.5 ℃,每年7—8 月降雨量最多,日最大降雨量达50 mm以上,11月到次年3月为冻结期,一般冻土深度为0.85 m,海拔高于5 200 m时为多年冻土区[20]。
本文InSAR 处理选用欧洲空间局(ESA)于2014 年发射的C 波段Sentinel-1A 卫星的干涉宽幅模式SAR 数据,极化方式为VV,每景影像幅宽250 km,由3 个具有一定重叠度的子条带拼接而成,用于反演滑坡整体形变速率及形变序列。数据时段为2017 年3 月—2019 年7 月,包括升轨数据55 景(path99),降轨数据69景(path106)。采用ESA 发布的精密轨道数据提高Sentinel-1A 影像的轨道数据精度;30 m 分辨率的SRTM DEM 数据去除地形相位和地理编码。
选用与SAR 数据时间段相对应的气象资料用于分析形变趋势,该气象资料来源于距离矮拉山约20 公里的德格气象站。通过查询2017—2019 年温度及降雨历时数据,该地区一年中气温及每月累积降雨量呈现周期性变化,最低气温出现在1 月份,为-12 ℃,最高气温出现在7 月,达到29 ℃,年温差最大为41 ℃;6—8 月为雨季,2018 年7 月降雨量达到225.3 mm,为近三年最高;每年11月至来年4月,降雨量普遍低于20 mm,平均每月降雪天数小于5天。
小基线集(SBAS)InSAR 是由Berardino 等[21]提出的一种时序InSAR 分析方法。该方法优点在于采用小基线(时空基线较短)干涉对组合,提高了SAR 数据使用率和监测点密度,并利用奇异值分解(SVD)方法高效获取监测点的形变速率及形变时间序列[22]。SBAS-InSAR 方法详细介绍见参考文献[23],本文简要介绍其主要数据处理步骤:
表1 Sentinel-1A影像数据信息Table 1 Image parameters of Sentinel-1A
(1)计算SAR 影像的时空基线,选择合适时空基线阈值并生成干涉对;
(2)对选取的干涉对,进行差分干涉及相位解缠处理;
(3)采用SVD 方法对干涉图组成的相位方程进行形变参数解算。
本文数据处理流程如图2所示,在数据配准时,选取日期为20180524的影像作为主影像,减少因时间基线较长造成相关性降低导致配准失败。在基线连接时,结合已有研究[24],时间基线阈值为60天,空间基线阈值为±200 m,时空基线分布情况如图3所示(蓝色连接点代表观测影像,线段表示相应的干涉对)。
图2 SBAS方法技术流程Fig.2 Flow chart of SBAS
图3 干涉对时空基线分布:升轨(path99)(a);降轨(path106)(b)Fig.3 Temporal-spatial baselines of interferometric pairs:ascending(path99)(a);descending(path106)(b)
利用SBAS-InSAR 方法处理了55 景升轨及69景降轨Sentinel-1 SAR 数据,获取了研究区卫星视线向(Line of Sight,LOS)向2017—2019 年均形变速率和形变时间序列。参考目前基于滑坡形变量速度的分类[25],假设沿LOS 向形变在-10~10 mm·a-1之间形变区为稳定区域。本文年均形变速率结果表明(图4),大部分区域处于稳定状态,存在明显形变信号的区域分布于山谷两侧,最大形变速率达到-78 mm·a-1。
升、降轨InSAR 形变结果的对比分析表明(表2):(1)地形平坦及地表稳定区域的形变信号较为一致;(2)在部分地表形变信号明显的区域,监测点密度、范围以及形变量大小存在差异,甚至LOS 向的形变方向相反(如格普滑坡,图4);(3)升降、轨In-SAR 结果在监测范围方面具有较好的互补性,但增加了解译的难度。
表2 典型区域升、降轨InSAR形变统计结果对比Table 2 Comparison of deformation statistics of typical landslides along different radar LOS directions
根据3.1分析,下面将结合降轨形变结果、卫星影像(来源于天地图,图5)以及现场考察情况(图6),重点分析朱乔、格普、隧道出入口四处典型滑坡形变时空特征,并选取四处典型滑坡共8个点位(位置见图4)定量分析滑坡的形变趋势(图7)。
图4 基于SBAS-InSAR方法获取的雷达LOS向年均形变速率:升轨(a);降轨(b)Fig.4 Mean ground deformation along radar LOS direction derived by SBAS-InSAR method:ascending track(a);descending track(b)
(1)朱乔滑坡范围较小,长度为300 m,宽度为200 m,与山脚高差约560 m。InSAR 结果[图4(b)]表明,2017—2019 年期间,滑坡后缘形变速率高于前缘,LOS 向年均形变速率为23~33 mm·a-1,发生两次变化,持续时间较长,累计形变值约为100 mm[图7(a)]。变形坡体表面植被稀疏,两侧发育多条冲沟[图5(a)],直通山脚埃曲河。
(2)隧道入口处滑坡位于矮拉山隧道入口西侧,冲沟尽头,长度为240 m,宽度为120 m,坡向近似南北方向[图5(b)],滑坡中部LOS 形变速率达到24 mm·a-1,累计形变达到90 mm。该区域国道旁多处设有水泥防护墩,实地调查发现部分已受到冲击破坏。
图5 滑坡卫星影像:朱乔(a);隧道入口(b)Fig.5 Satellite image of landslides:ZhuQiao(a);tunnel entrance(b)
(3)隧道出口处滑坡为小型碎石滑坡,图6(a)表明,地表风化现象严重,植被以草甸为主,最大LOS 向形变速率约20 mm·a-1,具有突发性,但变形相对较小[图7(c)]。
(4)格普滑坡为牵引式堆积体滑坡,两处形变坡体呈“八”字分布,每一处长度为600 m,宽度为350 m,在InSAR 结果,最大LOS 向形变速率为74 mm·a-1。如图7(d)所示,该滑坡变化期和稳定期存在明显分界线,累计形变量最大超过150 mm。现场发现该区域滑坡后缘存在明显张拉裂缝,西侧伴有发育冲沟,在谷底滑坡前缘存在松散堆积体[图6(b)]。
星载SAR 侧视成像产生的几何畸变不仅导致在高山峡谷部分地区无法提取有效的结果,而且也给升、降轨InSAR LOS 向形变结果的解译带来挑战。如图8 所示,形变表现为靠近卫星还是远离卫星,取决于坡面的朝向及坡度与入射角的大小,因此同一目标的位移在升、降轨结果中LOS 方向形变中可能相同,也可能相反。另外,当使用单一轨道数据处理时,在阴影区会出现无法测量的情况,因此在山区有必要联合升、降轨InSAR结果进行分析。
图8 卫星入射角(LOS)观测方向与沿坡向形变关系[26]Fig.8 Relationships between the Line of Sight(LOS)and the downslope displacements for different slope orientations[26]
四处典型滑坡中,朱乔、隧道入口处滑坡的具有相同的变化方向,而隧道出口处、格普滑坡变化方向相反(图9)。从定量角度分析,以格普滑坡为例,坡向为西,坡度范围处于10°到30°之间。Sentinel-1A 升轨数据飞行方向由南向北,入射角为33.8°。由于坡面朝向卫星,且坡度角小于卫星入射角,因此沿坡面向下发生的形变在LOS 向表现为接近卫星。当使用降轨数据时,卫星由北向南飞行,坡面背向卫星,坡角小于卫星入射角余角,此时为最佳观测情况,沿坡面向下的形变在LOS 向表现为远离卫星[26]。
图9 典型滑坡升、降轨InSAR形变特征对比Fig.9 Comparison of deformation characters of typical landslides along different radar LOS directions
通常情况下,在青藏高原内陆多年冻土区[27]及城区[28]等地形平坦区域,同一系列卫星升、降轨In-SAR 结果具有较好一致性。但在藏东高山峡谷地区,由于形变大多位于具有一定坡度且坡向不同的斜坡处,在SAR 卫星侧视成像及几何畸变影响下,难以通过分析升、降轨形变结果的相关性对两者的一致性做出评判。取决于卫星成像参数与坡体方向之间的几何关系,不同卫星、不同轨道数据,适用坡体也不同,例如Sentinel-1 卫星,升轨数据更适用于监测东坡向形变,而降轨数据适用于监测西坡向形变。
对于稳定区及地形平坦的区域,仍然可以通过对升、降轨结果进行交叉验证。而对于失衡斜坡,通过联合升、降轨数据进行反演形变,可以减少因几何畸变造成的观测盲区,有利于判断滑坡实际变化方向及形变序列的提取和分析,在提高滑坡监测识别准确度和覆盖度方面具有一定的应用价值。
大量研究表明,斜坡失衡主要由降雨、地震、人为活动等外部因素诱发,而在冻土区,季节性冻融作用也是导致斜坡失衡的重要因素,且海拔越高,冻融作用将更加强烈[29]。
研究区地表形变普遍存在季节性变形信号特征,为了分析冻融作用对滑坡形变过程的影响,我们重点分析了格普、朱乔、隧道入口三处滑坡形变过程与温度、降水变化的关系。由图10 可知,随着季节交替,冻融作用反复发生,三处滑坡形变过程呈现“阶梯下滑”形状,存在明显的平稳期和失稳期,与多年冻土地表形变“周期性”变化特征不同[30]。
图1 研究区位置Fig.1 The location of study area
平稳期主要出现在每年10—11 月至来年3—4月,这段时间该地区处于冻结期,温度低于摄氏零度,降雨量普遍低于20 mm,低温导致地下排水受阻加之气候干燥,在无地震及人类活动影响下,地表保持稳定状态(图10)。
图10 三个典型滑坡形变过程与气温、降雨的关系(灰色线为零度上下分界线)Fig.10 Correlation between landslide deformation and meteorological factors(grey lines are the boundary of the zero degree)
图10表明失稳期主要出现在春夏季(4—9月),呈现两个时段的持续变形,均始发于每年的3 月下旬,最低温度开始从零下转为零上,此时降雨量低于30 mm,且降雪较少,因此排除降水和融雪作用,冻土融化是诱发滑坡变形的主要因素,其诱发机理为“季节性冻结滞水促滑效应”[31];冻结作用使斜坡内地下水聚集滞留,导致坡体内含水层扩大,增大了斜坡区静、动水压力,影响斜坡土体强度和稳定性[32];当气温持续上升进入解冻期,土体软化,地下水自由流通产生侵蚀作用,从而导致斜坡失衡破坏[33];随着降雨增加,一方面给地面带来额外热量,使冻土进一步融化,另一方面外部水分补给渗透,起到润滑作用,导致坡体形变速度加快,当累计降雨达到最大时,具有最大形变速度。
藏东高山峡谷区气候具有雨热同期的特点,在降雨和冻融共同作用下,相比青藏高原内陆冻土区,该地区滑坡变形失稳机制更为复杂[34],冻融作用对滑坡变形的影响还需进一步研究。另外,同一处滑坡每年的变化趋势并不相同,这也是此类滑坡预测分析的难点。
本文以青藏高原东部317 国道矮拉山地区为例,利用SBAS-InSAR 技术和55 景升轨、69 景降轨Sentinel-1A SAR数据,提取了2017—2019年该该地区升降轨InSAR 地表形变信息,据此获取了该地区冻融滑坡体隐患的空间分布情况,并分析讨论了滑坡历史形变演化特征及影响因素。主要结论如下:
(1)研究区整体处于稳定状态(-10~10 mm·a-1),但存在部分明显形变区,主要形变区集中于山谷河流两侧,其LOS 向形变速率极值达到-78 mm·a-1。
(2)对比分析升降轨InSAR 地表形变范围、速率及形变序列,发现多处滑坡隐患,分别位于格普、矮拉山隧道出入口、朱乔,通过光学影像判读和实地调查验证了InSAR结果的可靠性。
(3)冻融滑坡形变过程呈现平稳期和失稳期交替出现的季节性变化特征,主要受冻融与降雨的综合影响,但具体作用机制还需进一步研究。
本文验证了升降轨时序InSAR 在高山峡谷冻土区滑坡大范围监测识别方面的有效性,可为此类地区滑坡灾害隐患的早期识别与监测防治提供重要参考。但受植被和积雪的影响,InSAR 失相干导致部分地区的形变信号难以准确提取,随着今后国内外L 波段SAR 卫星(如我国“陆探一号”、美国NISAR)的发射,将为植被和积雪地区地表形变监测提供更加可靠、丰富的研究资料。未来有必要加强滑坡三维形变反演及冻融和降雨影响的研究,更好地揭示冻融滑坡形变模式和形成机理。