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

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 典型老(古)滑坡隐患卫星影像(左)与Google Earth卫星影像(右)Fig.2 Satellite images of typical old landslide hazard (left) and Google Earth satellite image (right)

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 正在变形的滑坡隐患卫星影像Fig.3 Satellite images of landslide hazard under deformation

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

2.2.1 Stacking-InSAR技术原理

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

$V=-\frac{\lambda }{4\pi }\frac{\sum\limits_{i=1}^{M}{{{\varphi }_{i}}}}{\sum\limits_{i=1}^{M}{{{t}_{i}}}}$。 (1)

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

相应的误差传播公式为

$\Delta {{V}_{\text{disp}}}=-\frac{\lambda \cdot \sqrt{M}\cdot E}{4\pi \cdot \sum\limits_{i=1}^{M}{{{t}_{i}}}}$。 (2)

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

2.2.2 Stacking-InSAR技术处理流程

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

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 Stacking-InSAR处理流程示意图Fig.4 Stacking-InSAR processing flow chart

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

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 滑坡隐患野外验证Fig.5 Field verification of landslide hazard

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

图6

Fig.6

	Figure Option ViewDownloadNew Window
	图6 基于Stacking-InSAR技术的活动滑坡隐患形变探测Fig.6 Landslide deformation detected by Stacking-InSAR technology

图7

Fig.7

	Figure Option ViewDownloadNew Window
	图7 Stacking-InSAR技术探测滑坡隐患放大效果图Fig.7 Landslide amplification effect map by Stacking-InSAR technology

图8

Fig.8

	Figure Option ViewDownloadNew Window
	图8 研究区滑坡隐患分布Fig.8 Distribution of landslide hazard in the study area

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

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

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

研究区采用的评价单元为栅格单元, 栅格单元是滑坡易发性评价中应用最广泛的评价单元, 把研究区按照一定尺寸划分为规则的网格, 具有计算机处理方便、可操作性强等优点。栅格单元划分的大小由多种因素决定, 常用专家经验公式求取^[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。

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

图9

Fig.9

	Figure Option ViewDownloadNew Window
	图9 评价因子与滑坡数量关系Fig.9 Relationship between evaluation factors and landslide quantity

(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个区间。

(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.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为滑坡发生概率; X₁, X₂, …, X_n为各评价因子; B₀为常量, 取值-10.756; B₁, B₂, …, B_n为各评价因子回归系数。

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

Wij=Sijz/Sij, (7)

${{M}_{ij}}={{W}_{ij}}/\sum\limits_{j=1}^{m}{{{W}_{ij}}}$。 (8)

式中: S_ij和 Sijz分别为评价因子下的分级区间面积和滑坡面积; W_ij为 Sijz和S_ij的比值; M_ij为评价因子的归一化值。各评价因子的归一化结果见表1。

表1

Tab.1

表1(Tab.1)

表1 评价因子分级和因子归一化值 Tab.1 Grading and normalization of evaluation factor 评价 因子 因子分级 Sijz/km2 Sij/km2 Wij Mij 高程 < 1 300 m 1.464 61.506 0.024 0.211 [1 300.1 600) m 7.381 143.005 0.052 0.457 [1 600, 2 500) m 14.538 710.411 0.020 0.181 [2 500, 2 800) m 3.313 286.770 0.012 0.102 [2 800, 3 400) m 3.846 755.989 0.005 0.045 ≥ 3 400 m 0.086 525.738 0.000 0.001 坡度 < 20° 5.385 401.807 0.013 0.220 [20° , 30° ) 9.120 682.284 0.013 0.220 [30° , 40° ) 9.671 837.629 0.012 0.190 [40° , 50° ) 5.158 443.041 0.012 0.191 ≥ 55° 1.293 118.661 0.011 0.179 坡向 [0° , 22.5° ) 4.698 418.685 0.011 0.156 [22.5° , 112.5° ) 4.427 313.880 0.014 0.196 [112.5.5° , 202.5° ) 3.750 626.467 0.006 0.083 [202.5° , 247.5° ) 4.339 320.897 0.014 0.188 [247.5° , 337.5° ) 12.098 655.020 0.018 0.257 ≥ 337° 1.284 143.555 0.009 0.124 地表 曲率 < -0.5 7.675 693.828 0.011 0.299 [-0.5, 0.5) 14.173 1 071.537 0.013 0.357 ≥ 0.5 8.779 718.057 0.012 0.330 评价 因子 因子分级 Sijz/km2 Sij/km2 Wij Mij 工程 地质 岩组 松散岩体 0.471 29.797 0.016 0.304 较软弱岩体 6.080 421.116 0.014 0.278 较坚硬岩体 19.101 1 464.630 0.013 0.251 坚硬岩体 4.968 569.552 0.009 0.168 NDVI < 0.17 6.593 802.214 0.008 0.125 [0.17, 0.36) 7.145 587.581 0.012 0.181 [0.36, 0.58) 7.097 407.890 0.017 0.264 [0.58, 0.82) 6.287 348.445 0.018 0.273 [0.82, 1] 3.506 325.593 0.011 0.163 距道 路距 离 < 100 m 3.181 60.353 0.053 0.205 [100, 200) m 4.112 60.273 0.068 0.265 [200, 300) m 4.189 62.168 0.067 0.262 [300, 500) m 7.102 111.732 0.064 0.247 ≥ 500 m 12.044 2 188.875 0.006 0.021 距河 流距 离 < 100 m 3.399 104.659 0.032 0.233 [100, 300) m 8.058 196.316 0.041 0.294 [300, 500) m 7.838 194.169 0.040 0.289 [500, 800) m 6.030 265.799 0.023 0.162 ≥ 800 m 5.302 1 722.458 0.003 0.022 距断 层距 离 < 1 km 8.554 486.823 0.018 0.293 [1, 2) km 4.309 342.610 0.013 0.210 [2, 4) km 3.584 374.387 0.010 0.160 [4, 10) km 3.302 455.369 0.007 0.121 ≥ 10 km 10.878 824.212 0.013 0.220

表1 评价因子分级和因子归一化值 Tab.1 Grading and normalization of evaluation factor

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

Tab.2

表2(Tab.2)

表2 Logistic回归分析结果 Tab.2 Logistic regression analysis results 评价因子 回归系数 标准误差 Wals值 显著性 高程 5.423 0.764 10.013 0.007 坡度 8.153 1.029 37.521 0.000 坡向 5.624 1.328 15.097 0.003 地表曲率 3.334 1.202 5.305 0.016 工程地质岩组 6.603 1.354 13.482 0.001 归一化植被指数 1.621 1.805 2.135 0.032 距道路距离 5.025 1.530 9.522 0.009 距河流距离 2.329 0.407 4.229 0.021 距断层距离 4.069 1.018 6.128 0.013 常数项 -10.756 0.786 48.763 0.000

表2 Logistic回归分析结果 Tab.2 Logistic regression analysis results

将各因子回归系数代入式(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为距断层距离因子归一化值。

建立逻辑回归模型后, 利用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

Fig.10

	Figure Option ViewDownloadNew Window
	图10 研究区滑坡易发性等级分区Fig.10 Vulnerability level partition of landsides in the study area

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

图11

Fig.11

	Figure Option ViewDownloadNew Window
	图11 Logistic回归模型的ROC曲线Fig.11 ROC curve of logistic regression model

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

3.7.1 显著性检验

Logistic回归模型显著性检验指的是因变量与自变量之间相关关系的检验。采用Wald检验法^[23], 若各评价因子的显著性均小于0.05, 表示自变量通过95%的显著水平检验, 即逻辑回归模型的检验效果好^{[20, 21, 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)。