联合光学和SAR遥感影像的山区公路滑坡易发性评价方法

余绍淮, 徐乔, 余飞

(中交第二公路勘察设计研究院有限公司,武汉 430056)

0 引言

随着我国公路建设事业的发展,建设重点正逐步由东部向中西部地区、由平原微丘区向困难复杂的重丘山岭区转移,快速准确地对山区公路路线走廊区域进行滑坡易发性评价,可为后续道路的选线提供重要的参考资料; 但艰险山区地形、地质等基础资料严重匮乏,且气候条件复杂,地面数据获取困难,野外作业风险高。遥感技术的快速发展,为复杂山区的滑坡易发性评价提供了新的技术手段,当前遥感技术已广泛应用于区域性滑坡易发性评价中[1-4]。

滑坡灾害易发性评价是对滑坡灾害的时空分布和发生概率进行预测,可为滑坡灾害风险管理提供重要的决策依据。当前应用较多的易发性评价方法有2大类: 第一类是基于知识驱动的评价方法,如层次分析法[5]、专家打分法[6]等经验模型; 第二类是基于数据驱动的评价方法,包括逻辑回归[7]、信息量法[8]、确定系数法[9]等统计分析模型,及人工神经网络[10]、支持向量机[11]、随机森林[12]等机器学习模型。与经验模型和统计分析模型相比,基于数据驱动的机器学习方法能更好地分析评价因子和滑坡之间的非线性关系。随着对地观测技术的飞速发展,地表观测数据越来越丰富,基于数据驱动的机器学习模型已广泛应用到滑坡灾害易发性评价模型中[13-14]。

基于长时间序列的SAR影像获取的地表形变速率可直接反映滑坡运动状态,有助于有效识别不稳定区域,对滑坡监测具有重要意义[15-17]; 但仅利用地表形变速率进行滑坡易发性评价,而忽略诱发滑坡灾害的地形、地貌、地质等内在因素,会使评价结果易受地表形变误差影响。目前,基于机器学习的滑坡易发性评价模型中较少使用SAR影像,更多使用高程、坡度、坡向等地面静态数据,忽略了地表形变速率等地表动态数据,致使评价方法精度不高,实用性不强,难以推广到其他区域[18-19]。因此,在评价方法中增加利用地表形变速率,可使评价结果与区域滑坡实际分布情况更相符,进一步提高结果的准确性与可靠性。

为此,针对上述不足,本文在综合利用光学和SAR遥感影像基础上,提出一种结合地表形变速率的滑坡易发性评价方法,该方法先利用光学遥感影像和地形数据提取多种滑坡灾害静态因子; 然后利用随机森林(rondom forest, RF)模型计算滑坡易发性风险初始等级; 最后利用长时间序列SAR影像计算的地表形变速率对易发性风险初始等级进行修正,从而得到更加准确可靠的滑坡易发性评价分区图,为后续山区公路地质选线提供数据支撑。

1 基于多源遥感数据的滑坡易发性评价模型构建

1.1 滑坡易发性评价因子提取

滑坡灾害的发育是一个复杂的非线性过程,受多种因素综合影响,包括地形地貌、地质构造、地层岩性、岩体结构及地下水等内在地质因素,也有降雨、人工开挖、地震等外在诱发因素; 因此,充分提取滑坡孕灾环境信息对滑坡灾害易发性评价至关重要。为快速准确建立路线区域的滑坡易发性评价模型,利用高分遥感影像及DEM数据提取多种滑坡灾害因子; 其中,基于高分遥感影像的滑坡灾害因子提取主要是通过对区域地质、地形地貌及生态景观环境信息进行遥感解译识别,提取因子包括地层岩性、地质构造、地形地貌等; 基于DEM数据的因子提取主要是通过GIS分析获取,主要有坡度、坡向、坡长、水流流向、水流累积量等。

本文所用的滑坡易发性评价因子如表1所示。由于DEM数据和遥感影像的分辨率不同,致使基于二者所提取的滑坡灾害因子尺度不同,考虑到可免费获取的DEM数据的分辨率为30 m×30 m,本方法采用30 m×30 m的栅格单元作为滑坡易发性评价单元,并在此基础上提取滑坡灾害因子。

表1 滑坡易发性评价因子表

1.2 基于随机森林的滑坡易发性风险等级划分

随机森林(random forests,RF)是由Breiman L于2001年提出的一种以决策树为基础分类器的集成机器学习算法。该算法具有分类表现优异、人工干预少、运算非常快等优点,在遥感图像分类、特征选择等方面应用广泛[20-21]。

随机森林的基本原理是先利用Bagging重抽样方法从原始样本中抽取多个样本集; 然后对样本集构建多个分类回归树(classification and regression tree,CART),其中每个CART树的内部节点分裂时,对特征集进行一次随机抽样,并在抽取的特征集中进行最优选择; 最后,利用所有CART树的计算结果,通过均值或投票得出最终计算结果。随机森林模型的构建过程见图1。

图1 随机森林模型构建示意图

由图1可知,随机森林模型中引入了样本随机抽样及特征随机抽样,有效降低模型对样本噪声和异常值的敏感度,提高评价模型的准确率和稳定性。为此,本文基于滑坡灾害因子集,利用滑坡样本与非滑坡样本,训练基于随机森林的滑坡易发性评价模型,完成滑坡易发性风险初始等级的划分,主要包括以下步骤:

1)从N个原始滑坡训练样本中以有放回的方式重复取样N次,得到一个训练样本集; 重复上述过程k次,得到k个训练样本集;

2)针对每一个训练样本集,通过随机选取滑坡评价因子集作为分裂特征集,并选择最优的分裂方式进行内部节点的迭代分裂,且分裂过程中不做减枝处理,迭代完成后得到k棵CART树;

3)将生成的k棵CART树组合成随机森林,对输入的滑坡评价因子集进行分类预测,统计每棵CART树的分类结果,计算出滑坡易发性概率;

4)依据得到的滑坡易发性概率,利用自然断点法进行分级处理; 将研究区域划分低易发区、中易发区、高易发区和极高易发区4个等级,完成滑坡易发性风险初始等级的划分。

1.3 结合形变特征的滑坡易发性风险等级修正

上述滑坡易发性风险初始等级的划分是基于地面相关静态参数如高程、坡度、坡向等计算得到,而地表形变速率反映了地表动态变形,可直接反映滑坡的运动状态; 因此利用地表形变速率对滑坡易发性风险初始等级进行修正,可提高结果的准确性。为此,本文利用长时间序列SAR影像得到地表形变速率对易发性风险初始等级进行修正,从而得到更加准确有效的滑坡易发性评价分区图。

PS-InSAR[22]与SABS-InSAR[23]是2种应用最为广泛的时序InSAR技术,考虑到SABS-InSAR技术采用更加自由的时、空基线阈值来构成干涉对组合,削弱了空间失相关和大气延迟的影响,可得到长时间缓慢地表的形变规律,更适用于自然地表形变速率的监测。因此,本文采用SABS-InSAR技术获取该区域的形变速率,具体流程如图2所示。

图2 SABS-InSAR处理流程图

考虑到SABS-InSAR技术只能探测沿雷达视线方向的形变信息,而滑坡多沿斜坡面进行滑动,本文利用SAR成像信息及DEM数据,将上一步获取的形变速率转化为沿坡度方向的形变速率,以反映斜坡面的地表形变信息。斜坡面的形变速率(Vslope)的计算表达式为:

Vslope=VLOS/C

,

(1)

式中:VLOS为雷达视线方向的形变速率;C为坡面位移和雷达视线位移的比例系数,公式为:

,

(2)

hLOS=cosθ

,

(3)

,

(4)

,

(5)

式中:φ为地形坡度;α为坡向;θ为SAR卫星成像的入射角;ω为卫星轨道方向和正北方向的夹角;δ为地理正西方向与SAR卫星成像视线向的夹角。

SBAS-InSAR技术获取的坡向形变速率反映了区域坡面的动态变形状态,将坡面形变速率加入到滑坡易发性评价中,可提高模型的敏感度与准确性。为此,按照评价单元大小,将所提取的坡面形变速率图重采样至30 m×30 m,并依据评价单元所对应的形变速率大小,将评价单元划分为低形变区([0,15) mm/a)、中等形变区([15,30) mm/a)、高形变区([30,45) mm/a)及极高形变区(≥45 mm/a)4类; 然后,将滑坡易发性风险等级与形变等级关联,建立风险等级更新矩阵,实现对滑坡易发性风险初始等级的修正,得到最终的滑坡易发性评价分区图。滑坡易发性风险等级更新矩阵如表2所示。

表2 滑坡易发性风险等级更新表

表2中的数值表示滑坡易发性风险初始等级修正的级别,其中0表示不修正,+1,+2,+3分别表示将风险等级提高1级、2级、3级。例如,某个评价单元易发性风险初始等级为低易发区,形变等级为高形变区,根据表2,该评价单元的风险等级需提高2级,则该评价单元的最终风险等级为高易发区。

2 工程实践与分析

2.1 研究区概况

青海省沿黄公路共和至大河家段公路地处青藏高原东北缘,青海省东南部,起点位于青海省共和县,终点位于甘肃省积石山县,是青海省第一条沿黄河公路; 建设高质量的沿黄河公路,对于开发黄河上游水能资源,促进黄河谷地开发利用,完善公路网建设,适应地区经济快速发展的交通需求,具有十分重要的现实意义和社会意义。

工程地处青海省拉脊山断裂南侧,受拉脊山断裂带的控制作用,并随着青藏高原的急剧抬升和黄河不断下切,导致路线经过的龙羊峡至拉西瓦、李家峡库区、隆务峡至公伯峡等黄河峡谷段两岸地势陡峭,深切河谷发育,地形地貌复杂多样,泥石流、滑坡、崩塌等不良地质现象发育。为此,本文选取选择隆务峡至公伯峡段进行工程实践,以验证本文提出的联合光学和SAR遥感影像的山区公路滑坡易发性评价方法的准确性和有效性。研究区的范围与地形如图3所示。

图3 研究区范围与地形示意图

2.2 实验数据及处理

为准确评价研究区的滑坡易发性风险等级,选择QuickBird卫星影像、Sentinel-1A卫星影像及ASTER GDEM数据作为评价模型的基础数据。其中,QuickBird卫星影像用于提取滑坡易发性评价因子与历史滑坡区的解译,其分辨率为全色0.61 m,多光谱2.44 m。考虑到研究区内的地形起伏较大,先对QuickBird影像进行几何纠正与空间配准,以消除影像上存在的几何畸变; 然后对QuickBird卫星的多光谱与全色波段进行融合处理,以提高QuickBird影像的空间分辨率的同时保留其多光谱特性。在此基础上,利用处理后QuickBird影像提取出植被指数、地表覆盖等多种评价因子; 然后对QuickBird影像进行滑坡遥感解译,以得到研究区的历史滑坡区,并在此基础上选择训练样本。图4为研究区的历史滑坡遥感解译图。

图4 研究区历史滑坡遥感解译图

Sentinel-1A卫星影像用于提取研究区内的坡面形变速率,作为后续滑坡易发性风险初始等级调整的依据。实验选择覆盖研究区范围的37景2018年1月—2020年6月的Sentinel-1A升轨影像,影像的成像方式为干涉宽幅模式,选用VV极化方式的影像进行干涉处理,并利用精密轨道文件进行轨道误差校正,所用Sentinel-1A的影像信息如表3所示。

表3 Sentinel-1A 影像信息表

本文利用ENVI 5.3中的SARScape模块,对上述Sentinel-1A影像进行裁剪、配准、干涉、去平、滤波、解缠及地理编码等列干涉处理; 其中,所选择的超级主影像日期为2018-05-14,所设置的时间基线阈值与最大空间基线阈值分别为120 d与5%,所利用的干涉相对数目为165,平均绝对时间基线与绝对空间基线分别为70 d与53.88 m; 在此基础上,利用SBAS-InSAR方法得到研究区沿雷达视线方向的形变速率,结果如图5所示。

图5 研究区2018—2020年沿雷达视线方向形变速率分布图

利用上述结果,结合Sentinel-1A卫星成像信息及ASTER GDEM数据,利用前文提及的转换方法,将沿雷达视线方向的形变速率转换为沿坡面方向的速率。由于斜坡面和雷达视线方向的夹角接近90°时,会导致斜坡面形变速率值趋近无穷大,结合区域的实际情况与其他学者的研究经验[24-25],本文将比例系数C的阈值固定为±0.3,即:

。

(6)

另外,考虑到斜坡体的形变应沿着斜坡面向下运动,只有形变速率为负值的形变点才能反应斜坡体的实际形变状态,故剔除形变速率为正值的形变点,从而得到研究区沿斜坡面的形变速率分布图,结果如图6所示。

图6 研究区2018—2020年沿斜坡面方向形变速率分布图

2.3 结果分析与讨论

利用基于高分辨率QuickBird影像解译的历史滑坡区域,从中选取257个样本点作为训练样本,并从影像其他区域随机选择非滑坡点,二者比例保持1∶1,提取出各样本点的评价因子数据作为随机森林模型的训练样本。本文利用随机森林算法得到的研究区滑坡易发性初始评价分区图如图7所示。

图7 研究区滑坡易发性初始评价分区图

对比图6与图7可知,部分形变速率较大的区域被划分为低易发区,而长时间存在形变的区域往往发生滑坡的风险较高,因此仅考虑高程、坡度、坡向等地面静态因子对研究区的滑坡易发性进行评价,会导致评价结果的可靠性不高,与区域实际情况存在较大差异。为此,本文利用图6所示的斜坡面方向的形变速率对图7的滑坡易发性风险初始等级进行修正,得到最终的滑坡易发性评价分区图,结果如图8所示。经统计,图8所示的滑坡易发性评价分区图中,极高易发区、高易发区、中易发区及低易发区的面积占比分别为7.73%,13.67%,20.34%和58.26%,其中高易发区及极高易发区占比超过20%,表明研究区内大片区域受到滑坡的潜在威胁。通过前文解译的历史滑坡区可知,大部分滑坡区位于黄河峡谷及黄河两岸,与图7所示的高易发区、极高易发区的分布基本一致。此外,选取部分重点区域,将形变速率图、初始评价分区图与最终的评价分区图进行对比分析,以验证本文方法的有效性,对比图见图9。

图8 研究区滑坡易发性评价分区图

(a) 形变速率图 (b) 滑坡易发性初始评价分区图 (c) 结合形变特征的易发性评价分区图

由图9(a)可知,该区域存在大量形变速率较高的区域,其中A和B区域形变速率最高,具有较大的滑坡潜在风险,应具有极高的易发性风险等级; C,D,E区域也存在明显形变,应具有较高的易发性风险等级。但图9(b)所示的初始评价分区图将上述区域主要划分为低易发区,少量区域划分为中易发区与高易发区,与研究区形变速率较高的现状不符,而图9(c)所示的结合形变特征的评价分区图将上述区域主要划分为极高易发区与高易发区,与研究区的实际状况一致。对比分析结果表明,仅利用地表静态因子得到的滑坡易发性评价分区图精度有限,增加利用地表形变特征可显著提高评价结果的准确性与实用性。

频率比是指滑坡易发性等级内滑坡样本分级占比与易发性分级面积占比的比值,可用于定量评价滑坡易发性分级的准确性[26-27]。为进一步验证本文方法的准确性,利用研究区的历史滑坡遥感解译图,采用频率比检验的方式,对所得到的滑坡易发性评价分区图进行定量评价,评价结果如表4所示。

表4 滑坡易发性评价精度统计表

由表4可知,上述2种方式所得到的滑坡易发性评价分区图中,高易发区与极高易发区的分级面积占比均不超过研究区总面积的22%,但却分布有 37%以上的滑坡样本,远超过其他易发区等级,且频率比值从低易发区到极高易发区均显著增大,极高易发区频率比值也远大于其他易发性等级,表明上述2种方式能有效地评价研究区的滑坡易发性。

利用形变特征对滑坡易发性初始评级结果进行修正后,有26.73%的滑坡样本分布于极高易发区,22.28%的滑坡样本分布于高易发区,均高于初始分级图中的21.88%与15.18%,改正后分区图中的极高易发区的面积占比从6.99%提高到7.73%,频率比值从3.130 2提高到3.458 0,精度提升明显,表明增加利用地表形变特征可显著提高评价结果的准确性,尤其是具有较高风险等级的区域。修正后的分区图中低易发区的滑坡样本比例也显著降低,频率比值从0.749 1降低到0.555 6,但仍有32.37%的滑坡样本区域位于低风险区,这主要是由于使用的DEM数据精度与分辨率不高及滑坡样本点包含部分处于稳定状态的古滑坡所造成。总的来说,本方法得到的滑坡易发性评价分区图与研究区实际情况基本一致,验证了本文联合光学和SAR遥感影像对山区公路走廊进行滑坡易发性评价的有效性。

3 结论与展望

本文采用联合光学和SAR遥感影像进行山区公路滑坡易发性评价的方法,准确提取了多种静态与动态滑坡灾害因子; 基于随机森林算法综合利用地质、水文、地形地貌及地表覆盖等多种因子,实现了滑坡易发性风险初始等级的快速评估; 利用SABS-InSAR技术获取的地表形变因子,对滑坡易发性风险初始等级进行精准修正,并有效开展艰险山区公路走廊范围内的滑坡易发性评价工作。通过在青海省沿黄公路隆务峡至公伯峡段中的应用实践,表明本文方法可快速准确地得到山区公路走廊范围内的滑坡危险性评价分区图,适用于山区公路遥感地质勘察,可为后续公路地质选线提供数据参考。

但是,本文提出的滑坡易发性评价方法仍存在以下不足: ①滑坡样本需人工选择,影响工作效率与评价精度; ②采用的DEM数据精度与分辨率有限,降低了滑坡易发性评价分区图的准确性。因此,下一步工作应改进滑坡易发性评价算法,减少人工干预,采用高精度的DEM数据,提高滑坡易发性评价工作的效率与准确性。