刘洪江 王 磊
(1.安徽理工大学空间信息与测绘工程学院,安徽 淮南 232001;2.安徽理工大学矿山采动灾害空天地协同监测与预警安徽普通高校重点实验室,安徽 淮南 232001;3.安徽理工大学矿区环境与灾害协同监测煤炭行业工程研究中心,安徽 淮南 232001)
滑坡是我国最常见的地质灾害之一,每年发生的滑坡数量高达数千起。为减少滑坡所造成的危害,国内外学者先后提出数十种滑坡监测技术与方法。早期由于科技设备的限制,对于滑坡的监测大多为定性监测,如地表及地物大幅度变形、地下水、动物、气味出现异常等。这些方法操作简单、成本低,但精度也低,且时效性无法保证。因此,定量的滑坡表面变形监测技术逐渐发展起来。自2008 年以来,徐进军[1]、刘圣伟[2]等众学者先后利用三维激光扫描技术、机载激光雷达技术、激光扫描点云[3]、GPS 观测技术[4]等对滑坡进行了监测,监测结果均表明上述技术的可行性。
相较于上述方法,InSAR(合成孔径雷达干涉)技术应用于监测滑坡的时间更早一些,1996 年Achache[5]等人首次使用InSAR 技术就法国南部的一个滑坡进行研究,验证了该技术在滑坡监测方面的可行性。随后InSAR 技术便广泛应用到早期的滑坡形变监测,得到一众学者的认可。朱建军[6]、李达[7]、Carlà[8]、陈兴芳[9]等众多国内外学者[10-11]利用单一的InSAR 技术手段对滑坡进行监测,取得了一些技术上的进展,表明该技术可以很好的实现对滑坡的监测。但是单一的InSAR 技术手段在滑坡监测中仍存在很多不足。例如滑坡表面往往覆盖着不同程度的植被,受植被影响,InSAR 技术会在一定程度上产生几何形变和时空失相干等问题,另外地形表面的起伏形态也会在一定程度上影响In-SAR技术的监测结果。
随着遥感技术的发展,滑坡监测的方法也由单一监测方法逐渐走向多元化手段,不同遥感技术及遥感数据的结合,正慢慢被更多学者所认可。2019年陆会燕等[12]使用高精度的光学数据及Sentinel-1A 数据联合的方法就金沙江沿岸的滑坡进行了一个宏观性的识别监测,监测结果表明多元数据的结合可以更好的实现对滑坡的监测;2020年刘筱怡[13]采用多源遥感数据和多种数据处理方法,就古滑坡的判别及发展规律进行了研究。2021 年吴绿川[14]基于InSAR 技术和光学遥感,大规模地对贵州省部分地区进行了地表形变监测和危险形变区识别,并依靠NDVI和滑坡发育要素对危险区进行调查。
江达县金沙江流域是滑坡易发区域,给当地人民和经济带来巨大的影响。鉴于此,本研究以白格滑坡为例,拟采用SBAS-InSAR 技术和植被覆盖度反演的综合遥感技术,开展滑坡灾害的时空演变(特别是滑坡前兆)规律研究,以期为该区域今后的地灾监测和防治提供有效科学判据。
白格滑坡位于四川甘孜白玉县与西藏江达县交界处的白格村,地处金沙江(北纬31°04'53.21″、东经98°42'14.1″),具体位置如图1 所示。金沙江河谷呈东西走向,平均海拔3 000~4 000 m,属于典型的高山峡谷,受高原的隆起及金沙江江水持续的下切侵蚀,让两岸形成了险峻、陡峭的“V”形的高山峡谷地貌。处在金沙江右岸的白格滑坡体整体呈楔形,滑面形态呈阶梯状,地势起伏较大。而且白格滑坡整体巨大,坡顶高程近乎3 720 m,滑坡前缘高程约2 880 m,高差将近840 m,体积约33 280 m3,主滑方向长约为1 600 m,平均宽度约500 m,主滑方向坡度约82°。就地质条件而言,白格滑坡各层区域岩性十分复杂,主要为石英片岩,沉积岩,强度有所不同,越靠近深层,岩性强度越强[15]。
图1 研究区位置示意
本研究所需数据有32 景Sentinel-1A 的SAR数据、3 景Sentinel-2 的光学数据和空间分辨率为30 m的DEM数据Sentinel-1A 属于主动微波遥感卫星,数据分辨率高且可以免费获取,因此被广泛应用到地表变形监测领域中。本研究涉及SAR影像参数见表1。
表1 Sentinel-1A SAR影像参数
Sentinel-2 属于高分辨率多光谱成像卫星,其数据是一种获取地表植被的重要光学数据源。Sentinel-2影像参数见表2。
表2 Sentinel-2影像参数
在进行SBAS-InSAR 试验中,仅依靠SAR 数据会产生地形相位影响。为消除地形相位影响,需要用到DEM 数据进行辅助,本研究所采用的DEM 数据为NASA(美国国家航天航空局)所主持的SRTM(航天飞机雷达地形测绘使命)项目中所获取的空间分辨率为30 m的高程数据。
相关研究表明[16-18],大型滑坡发生前伴随着许多变化,这些变化发生的同时也会影响周边的环境条件,如水和土壤等,从而影响滑坡体上及滑坡体周边的植被。而光学遥感技术凭借其高空间分辨率、高时间分辨率等优势,能够对坡体植被的生长状况进行监测。因此,可以利用光学遥感技术,监测坡体植被生长情况,将滑坡表面植被生长状况作为监测滑坡的一个指标,间接监测滑坡的演变。
FVC(植被覆盖度)是一种描述植被生长状况及变化的重要指数,通常用于评估生态,气候,土壤等环境状况,其值在0~1 之间。植被覆盖度的表达式为式(1)。
式中:Fc为植被覆盖度;NDVImax为纯绿色像元的NDVI值;NDVImin为纯裸土像元的NDVI值。
为获取白格滑坡表面的植FVC,首先要获取滑坡表面的NDVI(归一化植被指数)指数。NDVI是一种可以反映植被长势和健康状况的指数,范围在-1~1 之间,当NDVI值为正时表明该地区覆有植被,且越接近1 时表明该区域植被健康程度越好,植被生长越茂盛,其表达式为式(2)。
式中:ρNIR为植被在近红外波段的反射率;ρR为植被在红色波段的反射率。
由于NDVI指数在一定程度上受大气、土壤、地形和植被类型等条件的影响,不能直接表观滑坡表面的植被覆盖度,为获取坡体表面真实的植被生长状况需要克服上述影响,而像元二分模型[19]可以有效的解决这一问题,此模型可以在一定程度上削弱上述条件的影响。
像元二分模型是一种简单且实用的植被覆盖度估算模型,该模型将遥感影像上像元单纯的分为绿色植被像元与非绿色植被像元。植被覆盖度与像元的遥感信息有一定的关联,因此只要能够确定纯绿色植被像元所包含的信息和纯裸土像元所包含的信息这两个参数的大小,就可以估算像元乃至影像的植被覆盖度。为了获得两参数的大小,选择用NDVI指数确定,NDVI值的大小就可以看成是像元中绿色植被信息和裸土信息相互作用的结果,相当于是像元的遥感信息。首先利用式(2)获取遥感影像的NDVI值,并对NDVI值进行统计,得到整幅影像的NDVI值分布情况,然后根据前人的经验,设置5%~95%的像元值为置信区间,以累积概率为0.05 的NDVI值作为NDVImin,0.95 处的NDVI值作为NDVImax。最后根据式(1)可计算出白格滑坡区域及附近的区域的植被覆盖度。
大量研究表明[6-14]滑坡发生前就有明显的形变位移特征,也有学者将滑坡的整个形变过程分为4个阶段,蠕滑阶段、滑动阶段、剧滑阶段、稳定阶段,其中蠕滑阶段和滑动阶段属于滑坡发生前形变位移特征。在蠕滑阶段,滑坡内部会存在许多应力分布不均匀区域,部分坡体会因此发生变形,并发展为连通性的滑动。到了滑动阶段,滑坡内部许多裂隙渐渐联通,内部的薄弱部分形成一个移动面时,滑坡两侧会产生若干裂缝,形似羽毛,滑坡前部会不断出现一些鼓状裂缝和放射性裂缝,此时,滑坡体形成。在这两个阶段中,滑坡表面均会产生较大的形变位移,因此将滑坡表面形变情况作为滑坡前期监测预警的又一指标是可行且有效的。为获取白格滑坡的时序性的形变位移情况,本研究选取Sentinel-1A数据对滑坡进行了时序性的监测。
短基线集技术(SBAS-InSAR)是Berardino[20]等于2002 年提出的,是一种时间序列分析方法,能克服传统D-InSAR在时间、空间及大气的不足。相较于PS-InSAR又可以获取更为连续的空间形变序列。
本研究获取白格滑坡表面时序性形变位移的方法为SBAS-InSAR 技术。SBAS-InSAR 处理流程如图2所示。
图2 SBAS-InSAR处理流程
SBAS-InSAR 技术在处理SAR 影像的过程中会将所有影像生成连接并同时确定超级主影像,然后生成一系列具有相应关系的干涉像对。之后这些干涉像对会进行自动配对,配对后的像对需要进行去平、相干性计算、相位解缠等操作,生成干涉图与解缠图。在得到干涉图与解缠图之后要对每一个干涉像对的质量进行人工检验,对干涉质量和解缠效果差的像对进行手动剔除,此处共剔除31 对质量较差的像对。然后对剩余的干涉像对进行轨道精炼与重去平,将剩余的恒定相位和残余相位坡道去除,并进行第一次估算形变速率和残余地形,在经过第二次解缠之后便可得到形变速率。最后在得到形变速率的基础上进行大气滤波,估算和消除大气相位,可得到时间序列上的位移结果。
利用3 景不同年份同一季度的Sentinel-2 号光学数据对白格滑坡体发生前的表面植被覆盖度进行反演,并依据FVC 值对植被覆盖度进行分划,并以颜色区分。将实验划分结果与空间分辨率为10 m的真彩色遥感影像进行对比,最终确定表示滑坡表面的低植被覆盖区域的FVC 值在0~0.4 范围内,表示滑坡表面的中植被覆盖区域的FVC 值在0.4~0.6 范围内,表示滑坡表面高植被覆盖区域的FVC值在0.6~1范围内,(划分规则见表3)。植被覆盖度反演结果如图3所示。
表3 植被覆盖度分划规则
图3 滑坡区真彩影像及植被覆盖
通过对Sentinel-1A 数据进行SBAS-InSAR 技术处理,获取了白格滑坡体在灾前的形变信息,形变速率信息如图4 所示。由于本研究没有实地观测数据,参照了杨成业[21]利用SBAS-InSAR 技术对白格滑坡体的研究,在数据选取时间跨度接近的情况下(杨成业:2017/07/4—2018/12/2,本研究2017/07/04—2018/10/15),所得实验结果近似,因此认为数据处理结果可靠。
通过遥感目视解译的方法对白格滑坡表面的植被覆盖程度进行分析,发现早在2016 年,白格滑坡体表面的植被的生长状态就异于其周边地区,(见图3),从图3(a)的2016 年光学真彩影像可以看出滑坡表面有明显的四处近乎裸露的区域(由A、B、C、D 四个字母表示四处裸露的区域),而这四处区域随着时间的推移也逐渐扩张,扩张趋势明显。为了较为真实地反映滑坡体表面植被生长状态,图3(b)、图3(c)分别展示了从2016—2018 年滑坡发生前滑坡体整体的植被覆盖状况及滑坡体表面中低植被覆盖所占区域。从2016 年的3(c)图至2018 年的图3(c)中可以看出D 区域的低植被区域扩张最为明显,扩张趋势由滑坡边缘向滑坡中部靠近,与B区域的低植被区扩张趋势较为接近,A区域的低植被区扩张面积次之,其低植被区域由滑坡后缘向滑坡前缘的方向扩张,这与C 的低植被区的扩张趋势较为接近。即A 区域与C 区域的低植被区呈东西扩张趋势,B 区域与D 区域的低植被区呈南北扩张趋势。整体来看,白格滑坡体表面的植被生长状况在2016—2018 年,呈现着异常的生长状态,表现为低植被区域逐年扩张、高植被区域逐年缩减的趋势。为了更为精准地反映白格滑坡体表面的植被生长状态,本研究对白格滑坡体表面的植被覆盖区域进行了定量分析。采用像元统计的方法对不同覆盖程度的像元进行统计。统计结果见表4,植被覆盖度占比如图5 所示。统计结果显示2016—2018 年白格滑坡体表面的低植被覆盖区域由起初的10.5%增至滑坡发生前的26.3%;中植被覆盖区域由起初的13.1% 降至滑坡发生前的10.9%;高植被覆盖区域由起初的76.4%降至滑坡发生前的62.8%。该结果与遥感目视解译结果一致。
表4 像元统计表
图5 植被覆盖占比
为证明白格滑坡发生前其表面植被退化并非特例,同时证明利用光学遥感技术监测滑坡发生前的植被生长状况可以有效地确定滑坡的潜在发生区,本研究又利用Landsat-8 光学数据提取了四川省茂县叠溪镇新磨村滑坡发生前三年的滑坡体表面的植被生长状况。新磨村滑坡植被覆盖度分级规则与白格滑坡植被覆盖度划分规则一致,其结果如图6所示。
图6 新磨村滑坡前植被覆盖图及滑坡前后遥感影像
新磨村滑坡发生于2017年6月24日,从图6可以看出,新磨村滑坡后缘在2015 年也存在一处裸露的区域(低植被覆盖区),且该区域也随着滑期临近逐渐变大。滑坡前三年的植被覆盖图显示新磨村滑坡在发生前,滑坡后缘与滑坡前缘表面的植被覆盖度存在不同程度的退化,主要表现在中低植被区域逐渐扩张,高植被区域减少。新磨村滑坡前的滑坡体表面植被覆盖度监测结果表明利用光学遥感技术监测滑坡发生前的植被覆盖度,可以有效的确定部分滑坡的潜在发生区域。
对SAR数据处理的结果表明,白格滑坡在发生前就存在着一定程度上的形变。从图4 可以看出滑坡体前缘与后缘呈现明显的沉降趋势,年平均下降速率最高可达-85 mm/a,中部部分区域呈现抬升趋势,抬升原因可能是滑坡后缘产生碎石岩土在此堆积。为了统计滑坡表面部分点位的具体形变值,选取滑坡表面A、B 两点,并就A、B 两点进行了28期的形变值统计,并依据统计结果做出点线图,如图7 所示。由图7 可知白格滑坡发生的一年前,坡体就开始产生形变,且随着滑期临近,形变程度也愈发明显。
图7 A、B两点形变统计
由于滑坡表面的植被生长状和形变程度与该地区降雨量有着密切的关系,本研究同时获取白格滑坡发生前的日降水量,昌都市2017 年7 月4 日至2018年10月1日的日降雨情况统计如图8所示。
图8 滑坡发生前日降雨情况
由图8可知该地区的主要降水集中在7月到10月份之间,大量的降水会使该地区的土壤水分含量过高,含氧量降低,植被根部会因此缺氧,进而影响植被生长。推测滑坡表面植被退化可能与该地降水量也有一定关系。同时频繁的降水会让滑坡体上土体重量增加,当雨水渗透内部岩石还会导致斜坡体抗剪能力下降,进而加快坡体形变趋势。结合图7和图8可以看出,降雨频繁发生期间,滑坡表面形变速率逐渐加快,因此在降雨频繁发生期间要格外注意形变变化大的坡体,以防滑坡灾害发生。另外通过对比图3植被覆盖图与图4滑坡区形变速率图,发现白格滑坡体表面植被覆盖度低的区域与形变量高的区域较吻合。可以推测,在白格滑坡体表面上的植被对于斜坡的演化和稳定性具有一定的影响。随着滑坡体表面植被的退化,裸露区域变大,植被的护土能力会有所下降,滑坡表面蠕变速度也就慢慢加快。这也说明利用光学影像对坡体表面植被覆盖度进行时序性监测是一种间接且可行的监测方法,可用于寻找坡体植被生长异常的区域,以确定滑坡潜在发生区。但由于植被发育不是地质灾害形成的根本原因,因此针对植被异常的区域可以进一步进行变形监测,从而缩小滑坡潜在发生识别区域,才能更有效地对滑坡的发生进行预警。
本研究结合光学数据与SAR数据,以白格滑坡为例利用遥感植被覆盖度反演技术及SBASInSAR 技术对白格滑坡进行了滑坡表面植被覆盖度变化监测与形变监测,获得白格滑坡发生前的植被覆盖度变化及年平均形变速率,根据监测结果得出以下结论。
①白格滑坡发生前表面植被就表现出异于周边地区植被的现象,表现为白格滑坡表面植被呈逐年退化趋势,裸露地区逐年扩张。
②白格滑坡在发生前就已出现较大形变趋势,随着滑坡发生日期的临近加剧,且降雨越频繁,形变速率越大,滑坡发生的可能性更大。
③通过对试验结果的分析,证明了植被覆盖度反演技术与SBAS-InSAR 形变位移监测技术的结合对部分滑坡的早期识别与预警有着一定的可行性,在地质灾害实际监测中具有较好的发展前景。
由于光学遥感数据相较于SAR 数据更容易获取且相对更容易处理,因此利用光学遥感数据对监测区域进行宏观性监测,选取出边坡区域表现异常的范围区域,再利用SAR技术对边坡表面沉降信息进行监测将会更有效的实现对滑坡的监测,同时也能在一定程度上降低滑坡监测成本。