国道213汶川—松潘段滑坡隐患遥感识别与易发性评价

黄艳琴， 李为乐， 许 洲， 李鹏飞， 铁永波

(1.成都理工大学地质灾害防治与地质环境保护国家重点实验室，四川 成都 610059；2.四川峨眉山四零三建设工程有限责任公司，四川 乐山 614200； 3.中国电建集团贵阳勘测设计研究院有限公司，贵州 贵阳 550081；4.中国地质调查局成都地质调查中心，四川 成都 610081)

0 引言

目前，我国公路交通建设正处于快速发展阶段，山区公路建设的投入也在不断加大[1]。山区地质环境复杂、山高坡陡、沟谷纵横，在地震或强降雨作用下滑坡灾害频发，严重威胁着山区公路的建设和运营安全。近年来，随着山区大量高速公路的建设和使用，每年都会出现公路滑坡灾害发生的情况，造成公路损毁、交通瘫痪甚至人员伤亡。例如： 2014年7月17日，国道213线K774+600 m处突发山体滑坡，塌方量达3 000 m3，造成国道213线约100 m的道路被掩埋阻断[2]； 2018年7月26日，国道213线K852 +600 m处发生山体滑坡，造成交通临时中断[3]。

山区公路滑坡隐患一般具有高位、高隐蔽性、突发性强等特点[4]，传统地面调查手段无法进行大范围提前识别。利用光学遥感技术可大范围识别出形态特征和变形迹象显著的滑坡[5]。而采用合成孔径雷达干涉测量(interferometric synthetic aperture radar，InSAR)技术可获取滑坡地表形变，对目前正在变形的活动滑坡进行识别。其中，差分干涉图叠加(Stacking-InSAR)技术相比于其他InSAR技术，能够有效抑制大气效应并减少DEM误差的影响，快速获取地表形变速率，被广泛应用于滑坡隐患的早期识别[5]。光学遥感和InSAR技术的结合可提高滑坡隐患识别的准确性。

滑坡易发性是对滑坡空间发生概率的定量评价[6]。随着遥感技术的发展，国内外学者依托遥感技术对滑坡易发性评价开展了大量研究工作，评价方法也由定性评价转为定量评价。常见滑坡易发性定量评价模型有Logistic回归模型、随机森林模型、支持向量机模型、信息量模型、神经网络模型等[7]。

国道213汶川—松潘段位于强震山区，构造活动强烈，滑坡灾害频发，造成公路被堵断或损毁，影响公路交通的正常运营。本文以国道213汶川—松潘段为研究区，结合光学遥感技术和Stacking-InSAR技术识别研究区滑坡隐患，并采用Logistic回归模型进行易发性评价，为提前预防研究区山体滑坡，保障国道213线安全运行，并为我国公路规划建设、防灾减灾提供理论参考和技术支撑。

1 研究区概况

国道213汶川—松潘段位于四川强震山区，自四川省汶川县映秀镇开始，途径茂县，向北延伸至松潘县川主寺镇，总长约230 km(图1)。研究区整体地势北西高、南东低，海拔在800～4 700 m之间，相对高差500～1 200 m，属于中-高山地貌。岷江是研究区内的主要水系，由南向北贯穿全境，受内部水动力和外部地质作用的影响，岷江河床不断向下深切，形成了如今的高山峡谷地貌。受高海拔西北风气流和印度洋西南季风的影响，研究区呈现青藏高原季风气候特点，冬季干冷，夏季湿润凉爽，日照量大，气温日差大、年差小。研究区内地质构造活动强烈，主要发育的断层有龙门山断裂带、岷江断裂带、虎牙断裂带、松坪沟断层、汶川—茂汶断层、较弧场构造等[8]，2008年汶川Ms 8.0级地震和2017年九寨沟Ms 7.0级地震均造成研究区公路不同程度的损毁。区内出露地层有新元古界震旦系，古生界寒武系、奥陶系、志留系、泥盆系、石炭系、二叠系，中生界三叠系以及新生界第四系，岩体质量较差、结构破碎。

图1 研究区位置和影像覆盖范围

2 滑坡隐患综合遥感识别

2.1 光学卫星遥感识别

利用2022年2月获取的高分2号卫星影像和Google Earth在线三维卫星影像数据，采用人工目视解译方法。光学遥感识别的滑坡隐患类型主要分为2类： 已发生整体失稳的老(古)滑坡隐患和正在变形的潜在滑坡隐患[9]。不同滑坡隐患类型的识别标志不同： 对于已经发生过整体失稳的古(老)滑坡，主要基于形态特征进行识别，整体平面形态呈圈椅状或舌形，后缘可见滑坡壁，前缘可见滑坡舌挤压河流，滑坡体上有群居居民或耕地[10](图2)； 对于正在变形的潜在滑坡隐患，主要识别变形迹象特征明显的区域，识别标志是坡体裂缝的形成过程以及前缘是否出现局部垮塌现象，针对这一类型滑坡隐患的识别，相对快速有效的方法是充分收集多时相光学遥感影像，通过对比不同时期的遥感影像特征解译滑坡隐患，从而提高遥感识别滑坡隐患的准确性(图3)。

图2 典型老(古)滑坡隐患卫星影像(左)与Google Earth卫星影像(右)

(a) 2011年 (b) 2016年 (c) 2019年 (d) 2021年

2.2 时序InSAR技术活动滑坡探测

InSAR技术是一种监测地表缓慢形变的遥感技术[11]，具有全天候、全天时、精度高、周期性重复观测等优点。本文选取2017年1月至2021年12月覆盖研究区域的288景Sentinel-1A降轨数据，采用Stacking-InSAR技术对研究区活动滑坡进行探测。

2.2.1 Stacking-InSAR技术原理

Stacking-InSAR技术由Price等[12]提出，其基本理论为： 假设地形相位误差、大气延迟相位误差等具有随机性，联合多期解缠后的差分干涉相位图建立形变速率和干涉相位间的线性函数模型，从而监测时间段的线性形变速率。利用该技术可以削弱InSAR技术中存在的随机轨道误差、大气延迟相位误差和地形误差等的影响，提高形变计算精度。Stacking-InSAR技术的计算方法可表示为[13]

。

(1)

式中：V为雷达卫星视线方向平均形变速率,mm/a；λ为波长,m；φi为第i个解缠后的差分干涉相位；M为干涉图数量；ti为第i个差分干涉相位的时间间隔,d。

相应的误差传播公式为

。

(2)

式中：ΔVdisp为雷达卫星视线方向形变速率误差值,mm/a；E为假定的单幅干涉图相位误差；M为干涉图数量；λ为波长,m；ti为第i个差分干涉相位的时间间隔,d。

2.2.2 Stacking-InSAR技术处理流程

通过Stacking-InSAR技术获取研究区地表形变速率图，可以大面积探测活动滑坡隐患。运用该技术首先对SAR影像数据进行处理，之后选取主从影像并配准生成干涉图。结合外部DEM与卫星成像参数，将DEM和SAR影像精确配准，并将DEM数据编码到SAR影像坐标系，生成差分干涉图。对差分干涉图进行空间滤波和相位解缠处理，解缠相位加权叠加后生成地表形变速率图。通过地表形变速率图筛选形变区域，从而确定滑坡隐患位置并识别滑坡隐患，数据处理流程见图4。

2.3 野外调查验证

结合现场调查与无人机航摄技术进行滑坡野外调查验证，主要工作包括判断室内解译的滑坡隐患是否为滑坡、滑坡边界的圈定是否正确并判定滑坡类型； 查明坡体裂缝的发育情况，如位置、方向、深度、宽度等(图5)。研究区共验证滑坡点185处，占滑坡总数的59.7%。其中有163处属于滑坡，排除滑坡22处，野外验证正确率为88.1%。

(a) 滑坡隐患边界范围 (b) 坡体变形及裂缝发育调查

2.4 滑坡隐患识别结果

利用光学卫星影像共识别滑坡隐患308处，基于InSAR技术共探测出活动滑坡隐患27处。活动滑坡主要集中分布于汶川—茂县段，其中25处为光学卫星遥感与InSAR技术共同识别，另外2处(B4、B16)为InSAR技术单独识别(图6、图7)。野外验证排除滑坡隐患22处，最终研究区共识别滑坡隐患288处(图8)。其中，汶川—茂县段滑坡隐患发育较多，分布较密集，该区域内地质构造活动强烈，地层岩性差异大，岩体结构较破碎，受公路或河流切割明显，临空条件较好。在地震或恶劣天气条件诱发作用下，发生滑坡的可能性增大。

图6 基于Stacking-InSAR技术的活动滑坡隐患形变探测

图7 Stacking-InSAR技术探测滑坡隐患放大效果图

图8 研究区滑坡隐患分布

3 滑坡易发性评价

3.1 评价因子的选取

通常易发性评价因子体系主要由内在因素和外在因素两部分构成[14]。内在因素一般侧重于斜坡本身的岩土体性质或地形地貌条件，其中地形地貌包括斜坡坡度、坡向和坡高等，岩土体性质包括地层岩性、斜坡及其周围环境的地质构造等。外在因素则被定义为除斜坡本身成灾环境条件之外诱发滑坡发生的因素，即外部环境条件对滑坡的影响，一种是自然环境的诱发因素，如地震或降雨； 另一种是人类工程活动的诱发因素，如交通建设、城镇建设、水利水电工程建设等。

根据上述对易发性评价因子体系的分析以及对研究区滑坡特点、分布规律、影响因素等的研究，选取高程、坡度、坡向、地表曲率、工程地质岩组、归一化植被指数(normalized difference vegetation index，NDVI)、距道路距离、距河流距离以及距断层距离共计9个评价因子，对研究区进行了滑坡易发性评价。

3.2 评价因子数据来源

评价因子数据来源主要分为4类： DEM数据来源于30 m GDEMV2数字高程数据(地理空间数据云)，用于提取高程、坡度、坡向、地表曲率因子； 地质数据来源于1∶20万地质图，用于获取岩性与断层分布； 水文数据来源于1∶20万地质图，主要提取河流与道路矢量； NDVI数据来源于2020年12月2日的全球植被覆盖Landsat 8卫星影像数据。

3.3 评价单元的确定

研究区采用的评价单元为栅格单元，栅格单元是滑坡易发性评价中应用最广泛的评价单元，把研究区按照一定尺寸划分为规则的网格，具有计算机处理方便、可操作性强等优点。栅格单元划分的大小由多种因素决定，常用专家经验公式求取[15]。同时，需综合考虑研究区地貌特征以及归一化植被指数栅格数据的空间分辨率，因此研究区采用的栅格单元大小为30 m×30 m，共划分出2 759 358个栅格单元。公式为

Gs=7.49+6×10-4×S-2.0×
10-6+2.9×10-15×S2。

(3)

式中：Gs为栅格单元建议大小；S为高程数据的比例尺倒数，高程数据比例尺为1∶50 000。

3.4 评价因子分级

分析各评价因子对滑坡发育的影响，统计各评价因子类别中滑坡发育的数量(图9)，从而对各评价因子进行分级。

(1)高程。高程反映区域的宏观地貌特征，表现为岩土体特征和地下水性质的差异，从而间接影响滑坡的发育条件。研究区处于中-高山区，海拔在828～4 694 m之间，将高程按300 m等间距分为10个区间，统计分析各区间内滑坡的分布。结果表明，滑坡在[1 300,1 600) m区间分布比较集中，该区间内滑坡较发育，易发性贡献率大，而在高程<1 300 m和≥3 400 m区间滑坡发育相对最弱，易发性贡献率小。根据滑坡在研究区高程的发育程度，将研究区高程分为： <1 300 m、[1 300,1 600) m、[1 600,2 500) m、[2 500,2 800) m、[2 800,3 400) m以及≥3 400 m共6个区间。

(2)坡度。斜坡的应力分布受坡度的影响，坡度大小不同，则斜坡的应力与稳定性不同。研究区坡度在0°～87°之间，为了更好地分析滑坡在不同坡度范围内的发育数量和分布情况，将研究区的坡度因子以等间距5°划分为11个坡度区间。分析结果显示，滑坡在[35°,45°)区间内分布比较集中，该坡度区间滑坡较发育，易发性贡献率大，而坡度在<20°和≥50°区间内滑坡分布较稀疏，滑坡不发育，易发性贡献率小。根据滑坡在研究区坡度的分布情况，将研究区坡度分为<20°、[20°,30°)、[30°,40°)、[40°,50°)以及≥50°共5个区间。

(a) 高程-滑坡数量 (b) 坡度-滑坡数量 (c) 坡向-滑坡数量

(d) 地表曲率-滑坡数量 (e) 工程地质岩组-滑坡数量 (f) NDVI-滑坡数量

(g) 距道路距离-滑坡数量 (h) 距河流距离-滑坡数量 (i) 距断层距离-滑坡数量

(3)坡向。坡向是日照时数和太阳照射强弱的客观反映[16]。研究区坡向在0°～360°之间，将其按方位划分为8个方位向，并统计滑坡在各方位向的数量和分布情况。分析统计结果可知，研究区东朝向[67.5°,112.5°)、西南向[202.5°,247.5°)以及南朝向[247.5°,292.5°)滑坡分布较为密集，易发性贡献率较大，而坡向[0,22.5°)与[337°,360°)区间的滑坡分布少，易发性贡献率相对最低。根据滑坡在各朝向的分布情况，将研究区坡向划分为[0°,67.5°)、[67.5°,112.5°)、[112.5°,202.5°)、[202.5°,247.5°)、[247.5°,337.5°)以及≥337.5°共6个区间。

(4)地表曲率。地表曲率的大小指示了斜坡的形态特征。在不同斜坡形态条件下，滑坡的发生具有一定的差异性。根据滑坡发育形态特征进行分类： 地表曲率<-0.5为凹型坡，地表曲率[-0.5,0.5)的为直线型坡，地表曲率≥0.5为凸型坡。

(5)工程地质岩组。斜坡稳定性受自身物质组成的影响，物质组成不同，岩土体的物理参数不同，则岩土体的物理性质不同。研究区岩组的分类标准参考《GB 50218—2014工程岩体分级标准》[17]，共分为4个等级，分别为松散岩体、较软弱岩体、较坚硬岩体和坚硬岩体。统计分析各类岩组中的滑坡发育数量，可知松散岩体中滑坡数量最少，相对最不发育； 较坚硬岩体和坚硬岩体中滑坡发育最多，可能与岩体遭受风化剥蚀有关，风化越严重，岩体破碎程度越明显，斜坡越不稳定，发生滑坡灾害的可能性越大。

(6)归一化植被指数。植被对滑坡发生的影响主要表现为植被的根茎深入地底，可以起到根固的作用，还可以缓冲坡面水流流速，从而极大程度地减少山体水流对坡面的破坏，减少滑坡灾害的发生。NDVI能反映某一区域内植被的覆盖程度，利用Landsat8遥感影像求取NDVI，计算公式为[18]

NDVI=(IR-R)/(IR+R)。

(4)

式中：IR为近红外波段；R为红外波段。NDVI的取值区间为[-1，1]，其值越大表示植被覆盖程度越高，生长状况良好。结合研究植被分布规律，本研究采用自然间断法，将NDVI划分为<0.17、[0.17,0.36)、[0.36,0.58)、[0.58,0.82)以及[0.82,1]共5个区间。

(7)距道路距离。在道路建设过程中，坡脚开挖、坡顶荷载等人类工程活动会不同程度地破坏岩土体完整性和天然状态下的自然结构。按照100 m间距划分出11个区间，得出距道路距离与滑坡数量的关系。在距道路100 m范围内滑坡分布最多，道路对滑坡发生的影响最大。根据滑坡在不同距道路距离的分布情况，将研究区划分为<100 m、[100,200) m、[200,300) m、[300,500) m以及≥500 m共5个区间。

(8)距河流距离。河流掏蚀作用会使岩体产生裂缝，越靠近河流的岩土体含水率越大，抗剪强度降低，斜坡稳定性随之降低。以100 m为间隔将距河流距离划分为11个区间，统计各区间的滑坡发育情况。分析表明，滑坡主要集中分布在距河流距离300 m的范围内，而距离河流大于800 m的区域滑坡总体发育较少。通过分析滑坡发育数量及分布与河流距离的关系，将研究区距河流距离划分为<100 m、[100,300) m、[300,500) m、[500,800) m以及≥800 m共5个区间。

(9)距断层距离。区域构造活动的强烈程度会影响节理裂隙发育和岩土体结构，可直接或间接影响滑坡的形成和破坏。以1 km为缓冲间隔将研究区的距断层距离划分为11个区间，统计分析各缓冲区域内滑坡的发育数量。结果表明，在[0，1) km区间内滑坡相对最发育，而在[2，8) km区间滑坡发育数量呈现逐步降低趋势，当距断层距离≥10 km时，滑坡发育数量呈增长模式，由于研究区区域的形状较长，断层分布跨度大，在此区间内断层相对集中，滑坡发育数量增多。统计不同断层缓冲区内滑坡的数量，将研究区断层缓冲区划分为<1 km、[1，2) km、[2，4) km、[4，10) km以及≥10 km共5个区间。

3.5 评价模型

3.5.1 Logistic回归模型

Logistic回归模型是一种预测性建模技术，其基本原理为利用回归分析方法研究因变量(目标)与自变量(预测器)之间的关系，由一个或多个滑坡影响因素(自变量)来判定一个滑坡(自变量)是否发生(0或1)的概率[19-20]。Logistic回归函数表达式为

Y=Logit(P)=ln[P/(1-P)]=
B0+B1X1+B2X2+…BnXn，

(5)

P=1/(1-exp(-Y))。

(6)

式中：Y为二元值(0或1)，表示滑坡发生与否，发生Y=1，不发生Y=0；P为滑坡发生概率；X1,X2,…,Xn为各评价因子；B0为常量，取值-10.756；B1,B2,…,Bn为各评价因子回归系数。

本文选取Ayalew[21]提出的评价因子量纲归一化方法获取归一化值。计算公式为

，

(7)

。

(8)

表1 评价因子分级和因子归一化值

(续表)

3.5.2 Logistic回归模型的构建

建立逻辑回归模型的步骤如下：

(1)将研究区矢量滑坡分布图转为栅格图层，发生滑坡的区域赋值为1，未发生滑坡的区域赋值为0。

(2)根据所求的归一化值，分别将9个评价因子图栅格化。

(3)在研究区域范围内，随机选取230个滑坡点(占总滑坡点的80%)，将剩余的58个滑坡点(占总滑坡点的20%)作为本次评价的检验点。同时，借助ArcGIS软件中的随机点生成工具，在研究区范围内,排除滑坡点以相同的点间距500 m生成随机点，创建与滑坡点数量相同的230个非滑坡点，共460个点作为本次评价的总样本点。

(4)将步骤(1)和步骤(2)中的值赋在460个样本点上，得到样本点属性表，将其导入SPSS软件中，利用回归分析工具进行总样本点的逻辑回归分析，得出回归分析结果(表2)。

表2 Logistic回归分析结果

将各因子回归系数代入式(5)和式(6)，建立Logistic回归模型为

P=1/[exp(-10.756+5.423a+8.153b+

5.624c+3.334d+6.603e+1.621f+

5.025g+2.329h+4.069i)] 。

(7)

式中：P为滑坡发生概率，取值范围为0～1；a为高程因子归一化值；b为坡度因子归一化值；c为坡向因子归一化值；d为地表曲率因子归一化值；e为工程地质岩组因子归一化值；f为NDVI因子归一化值；g为距道路距离因子归一化值；h为距河流距离因子归一化值；i为距断层距离因子归一化值。

3.6 评价结果与分析

建立逻辑回归模型后，利用ArcGIS的栅格计算器工具，将各评价因子回归系数的栅格图层进行叠加运算，得到研究区滑坡灾害易发性概率，概率区间为[0.008，0.963]，采用自然断点法[22]，按照易发性程度将易发性概率进行分级，可分为极低易发区[0.008，0.176)、低易发区[0.176，0.323)、中等易发区[0.323，0.481)、高易发区[0.481，0.688)、极高易发区[0.688，0.963]共5级，采用重分类生成研究区滑坡易发性等级分区，结果如图10所示。

图10 研究区滑坡易发性等级分区

由图11可知，国道213汶川—松潘段滑坡分布特征如下： ①极高易发区和高易发区主要集中分布在汶川县映秀镇至茂县石大关乡一带，该区域地势陡峻，斜坡坡度较大，岩体破碎，易发生滑坡灾害。其中，极高易发区面积约106.75 km2，占研究区总面积的4.31%，发育滑坡92处，滑坡点密度为0.862处/km2； 高易发区面积约505.25 km2，占研究区总面积的12.15%，发育滑坡54处，滑坡点密度为0.179处/km2； ②滑坡的中易发区面积约492.54 km2，占研究区总面积的19.83%，发育滑坡37处，滑坡点密度约为0.075处/km2； ③低易发区和极低易发区主要集中在茂县石大关乡以北至松潘县城一带。其中，低易发区面积约505.25 km2，占研究区总面积的20.34%，发育滑坡25处，滑坡点密度为0.049处/km2； 极低易发区面积约1 077.23 km2，占研究区总面积的43.38%，发育滑坡22处，滑坡点密度为0.021处/km2。低易发区和极低易发区范围内地势相对平缓，斜坡坡度较小，坡度大部分约20°，发生滑坡的可能性较小。

由分析结果可知： 随着滑坡易发性等级的升高，易发性面积占研究区总面积的比例减少，发育滑坡数量增加，因而滑坡点密度不断增大，说明易发性分区是合理的。

3.7 评价结果检验

3.7.1 显著性检验

Logistic回归模型显著性检验指的是因变量与自变量之间相关关系的检验。采用Wald检验法[23]，若各评价因子的显著性均小于0.05，表示自变量通过95%的显著水平检验，即逻辑回归模型的检验效果好[20-22]。从表2可知，各评价因子显著性均小于0.05，自变量通过了95%的显著水平检验，模型显著性较高。

3.7.2 精度检验

易发性评价结果精度可以采用受试者工作特征曲线(receiver operating characteristic curve,ROC)进行检验[24]。ROC曲线是依据多种不同种类的二分类方式，以敏感度为纵坐标，以1-特异度为横坐标绘制的曲线，表示敏感度和特异的相互关系。ROC曲线利用曲线下面积(the area under the ROC curve,AUC)来评定评价结果的精度，AUC的取值范围为[0.5，1]，其值越大，检验越具有价值，评价结果的精度越高。检验结果得出AUC值为0.837，表明评价结果的精度较高，说明采用逻辑回归模型对国道213汶川—松潘段滑坡易发性进行评价是较为客观准确的(图11)。

图11 Logistic回归模型的ROC曲线

5 结论

本文以国道213汶川—松潘段为研究区，综合运用光学卫星遥感与InSAR变形探测技术，结合统计分析与逻辑回归评价方法，开展研究区滑坡隐患综合遥感识别和易发性评价研究，得出以下结论：

(1)采用光学遥感和InSAR变形观测技术对研究区进行滑坡隐患识别，共识别出滑坡隐患288处，其中InSAR技术探测出活动滑坡27处，说明综合遥感识别技术在该研究区具有较好的识别效果。

(2)极高易发区和高易发区主要集中分布在汶川—茂县一带，总面积约612 km2； 极低易发区和低易发区主要集中在茂县以北至松潘一带，总面积约1 582 km2。同时，采用ROC曲线进行易发性结果精度分析，AUC值为0.837，表明评价结果的精度较高，说明采用Logistic回归模型对国道213汶川—松潘段滑坡隐患进行易发性评价是较为客观准确的。

(3)分析9个评价因子的逻辑回归系数，结果表明坡度、高程、坡向、距道路距离4个评价因子对国道213汶川—松潘段滑坡易发程度的贡献相对较大，其中坡度[35°,45°)、高程[1 300,1 600) m、坡向[67.5°,112.5°)、距道路距离<100 m是滑坡最容易发生的地区。