基于SVF地形可视化方法的滑坡识别<br/>——以四川省丹巴县城典型滑坡为例

基于SVF地形可视化方法的滑坡识别
——以四川省丹巴县城典型滑坡为例

2021-12-22董秀军刘小莎佘金星

成都理工大学学报（自然科学版） 2021年6期

郭晨, 许强, 董秀军, 潘星, 刘小莎, 佘金星

(1.地质灾害防治与地质环境保护国家重点实验室(成都理工大学)，成都 610059； 2.成都理工大学地球科学学院，成都 610059； 3.四川测绘地理信息局测绘技术服务中心，成都 610081)

中国是地质灾害频发的国家之一，其中西南山区是中国地质灾害的高发地区。近年来发生的多起重大地质灾害事件愈发凸显了滑坡隐患识别的重要性[1-4]，而传统的地质灾害巡排查手段已经远不能满足当前防灾减灾的需求。随着遥感技术的不断发展，借助于高精度光学遥感技术、合成孔径雷达干涉测量技术(interferometric synthetic aperture radar, 简称InSAR)、激光雷达测量技术(light laser detection and ranging, 简称LiDAR)、无人机(unmanned aerial vehicle, 简称UAV)摄影测量技术等现代高精度对地观测技术的滑坡隐患识别方法得到了广泛应用[5-8]。其中机载激光雷达作为一种新型主动航空遥感技术，能够直接获取观测区三维激光点云，利用激光多次回波的特性结合相应的滤波算法可以有效去除地表植被、建筑物等，从而获得真实的地面高程信息，进而生成分米级甚至厘米级的高分辨率数字高程模型，为滑坡的识别与测绘带来了前所未有的机遇。其穿透植被直达地表的能力，更是为植被密集覆盖的山区地质灾害的识别提供了新的解决方案。从2000年开始，奥地利、美国、比利时、土耳其等国家先后利用机载激光雷达技术开展滑坡的调查识别等工作，取得了丰硕的成果[9-11]，推动了机载LiDAR技术在地质灾害防治中的应用。且意大利、日本以及中国香港、中国台湾地区已经完成了全域机载LiDAR飞行，为地质灾害防治提供了基础数据和科技支撑[12-13]。

利用机载激光雷达技术进行滑坡识别首先要对高分辨率DEM进行三维可视化处理，进而通过观察地貌特征进行滑坡的识别解译。山体阴影是最常用的DEM三维可视化方法，它通过模拟太阳光照射引起的明暗变化从而对地形进行渲染，使二维地形看起来有三维立体的效果，经验丰富的地质工作者可以根据滑坡展现的各种地貌特征来进行识别解译[14]。但是由于不同的太阳方位角和高度角会对山体阴影图的效果造成极大影响，例如太阳光垂直照射斜坡的坡面时，会造成斜坡上没有阴影；而坡面与太阳光平行时会造成斜坡上阴影过重，从而使得潜在滑坡的微小特征难以被发现；当光线方向与斜坡走向大致平行时，滑坡特征通常最明显。如图1所示的滑坡区山体阴影图，在方位角为315°时效果最佳；但仍有部分滑坡特征受光照角方向影响而不易被发现，尤其是在方位角为45°时由于光照过量，滑坡的细节特征几乎被掩盖了。由此可见不同角度的光照山体阴影图会在一定程度上影响解译人员识别滑坡的微小地貌特征，进而影响滑坡识别的准确性和可靠性。为了减轻单一光源对可视化影像的影响，国内外学者相继利用DEM的衍生物如晕染图、坡度图、粗糙度图等手段进行滑坡识别工作[15-16]，但是普适性欠佳，效果不甚理想。

图1 不同太阳方位角下的滑坡山体阴影图Fig.1 Landslide shadow map at different solar azimuth angles

天空视域因子(sky view factor，SVF)是对地表形态开阔度的定量描述，反映了视线范围内天空被遮挡的程度[17]，是研究城市热辐射、热岛效应的重要参数[18-20]。在考古学中，研究人员通过机载LiDAR获取高分辨率DEM，结合天空视域因子进行考古地形可视化，从而通过解译地面甚至地下遗迹的影像特征来探测与发现考古遗存[21-22]。如M.A.Canuto等[23]利用2000多平方千米的机载激光雷达数据结合SVF等多种可视化技术揭示了危地马拉北部密林深处的古玛雅文明。然而目前将SVF可视化方法应用于滑坡识别的研究还相对较少。

本文基于SVF三维地形可视化技术的滑坡识别方法，以获取的丹巴县机载LiDAR数据为基础，通过识别出的3处典型滑坡为例，比较常规山体阴影图与SVF可视化影像在表征滑坡等灾害特征上的差异，探讨最大搜索半径与搜索方向对SVF影像的影响。

1 数据与方法

1.1 数据

丹巴县位于四川省甘孜州东北部，属于青藏高原东南缘横断山区邛崃山脉，境内地质环境脆弱，崩塌、滑坡、泥石流等灾害频发，素有“地质灾害博物馆”之称。本次机载LiDAR数据由四川测绘地理信息局于2019年2月13日至2月19日利用小松鼠高原加强版直升机挂载“天眼”SE-J1500激光雷达设备获取，飞机飞行速度约120 km/h，激光发射器最大视场角为50°，激光脉冲频率为50～550 kHz，获取的调查区平均点密度大于10 m-2。在飞行获取机载LiDAR点云数据的同时，和激光扫描仪同轴挂载了8 000万像素的光学相机，获得的航空影像分辨率优于0.2 m。数据处理采用TerraSolid软件，流程主要包括点云去噪、航带平差、点云分类、DEM数据构网与栅格数据输出。因数据量比较大，通过编写宏命令程序来实现点云自动分类，进而结合人工精细化分类将获取的激光点云数据分为地面点和植被点。对处理后的机载激光点云数据利用ArcGIS 10.3软件进行处理，生成了0.5 m分辨率的数字地面模型(DSM)以及去除植被和建筑物的数字高程模型(DEM)，经检验DEM水平和高程精度均优于0.5 m。

1.2 Sky View Factor计算原理

天空视域因子利用漫反射的方法解决了单一光源照射下的阴影问题，其含义是某一点处对天空的可见度，是一个无量纲值[24]。利用DEM计算SVF要基于以下假设：①整个半球的光线是均匀的；②半球内没有其他光源；③忽略半径10 km内地球曲率的变化。

照射到DEM上某个点上的光量通常与该点对天空的可见度有关。例如陡峭的山脊或山峰比山谷更明亮，因为它们可以从周围的天空获得更多的照明。计算SVF的最简单方法是测量它的立体角(图2-A)，如果在 DEM 上选取一个观测点，那么立体角定义为: 以观测点为球心，构造一个单位球面，任意物体投影到该单位球面上的投影面积，即为该物体相对于该观测点的立体角。立体角是单位球面上的一块面积，它表示在某一点上可以包含的最大范围(以空间为单位观察者上方半球的投影面积)

图2 天空视域因子Fig.2 Calculation principle of sky view factor(根据K.Zakšek等[24] 修改)

(1)

式中φ、λ分别表示半球内的纬度和经度。由于天空可见度受到地物的限制，假设在水平面以上，水平方向上的方位角具有相同的仰角，则立体角为

(2)

式中γ表示相对于水平地表上的仰角。由于实际观测点在各个方向的地表高度不同，因此仰角也不相同，所以可以通过在水平面上选的一定数量的方位角，逐一计算其水平面垂直仰角γi来计算立体角(图3)。

(3)

式中:n为搜索方向的数量;γi为不同方向的高度角。则天空视域因子(SVF)的计算公式为

(4)

计算出的SVF值在0到1之间，主要与2个参数有关，即水平搜索方向n和最大搜索半径R。SVF=1代表观测点上空整个半球都是可见的，地表面以上通视良好，图像上灰度值表现为亮值；而SVF=0意味着观测点上方几乎没有可通视的天空，图像上显示为暗值。

由于SVF是一个代表天空可视度的指标，往往用来研究城市建筑物接受太阳辐射率或阳光散射等问题，一般来说当某一点处的天空大部分可见时，地表温度变化更快，因此被广泛用于城市热岛效应、地表热平衡等研究中[25-26]。我们利用数字地形分析软件Saga7.0实现了基于SVF的三维地形可视化。与其他可视化方法比较，SVF在保持总体地形地貌特征的同时，忽略单一光源影响，能较好地显示微小地貌特征，从而利于解译人员从地貌上观察滑坡变形破坏后在地表留下的滑坡壁、滑坡台坎、鼓丘、洼地、局部隆起等微地貌特征，以此来识别滑坡。

2 结果分析

图3-A所示为利用机载激光雷达数据识别出的丹巴县甲居乡2处滑坡，由于该滑坡区域植被覆盖较密，从光学影像上很难识别这两处滑坡。图3-B所示为在225°/45°的太阳方位角及高度角下生成的山体阴影图(hillshade)，滑坡整体特征被阴影掩盖，尤其是滑坡后缘陡坎和两侧边界的辨识度较低(图3-B中黄色箭头位置)。图3-C为SVF地形可视化影像，滑坡圈椅状特征极为明显，尤其是山体阴影图中较难辨识的滑坡后缘陡坎及滑源区两侧边界。另外，滑坡堆积区整体呈扇状位于滑坡前缘，与滑源区分界较为明显(图3-F)。对这两处识别出的滑坡进行野外复核证实了识别结果的准确性。图3-E中左侧滑坡滑源区边界清晰可见，受公路开挖影响存在局部垮塌现象，现场调查发现该滑坡现今仍处于变形中，其中后缘景区公路受滑坡变形影响出现了一条长约20 m、宽约0.1 m、最大高差约0.15 m的弧状裂缝(图3-D)。

图3 丹巴县聂拉村滑坡Fig.3 Pictures showing the Niela Village landslide in Danba County(A)正射影像; (B)山体阴影图; (C)SVF可视化影像; (D)滑坡现场照片; (E)滑源区UAV影像; (F)滑坡特征图

图4为丹巴县五里牌新区滑坡，正射影像显示斜坡的坡脚存在大量建筑物，且均在滑坡影响范围之内(图4-A)。在最佳光照角的山体阴影图上虽然也能较好地识别出滑坡整体特征，但是部分区域还是偏亮，且未能较好地显示滑坡体上堆积的块碎石以及局部细节特征(图4-C)。相比之下，基于SVF生成的地形可视化影像，能非常清晰地再现滑坡的地貌特征(图4-B)，尤其是对于局部高差变化明显的区域，如滑坡后缘右侧边界处，受滑坡滑动影响拉裂的破碎岩体(图4-D)，在SVF影像上清晰可辨，经无人机航拍验证该处确实存在较为明显的岩块(图4-E箭头所示)。滑坡前缘左侧受局部开挖影响存在一处局部垮塌现象。相比于山体阴影图(图4-C箭头位置)，SVF影像能够非常清晰地展示该处崩塌的边界特征(图4-F)。图4-G为现场验证拍摄照片，现场调查发现在滑坡前缘修筑的挡墙已出现局部变形(图4-F)。

图4 丹巴县五里牌滑坡Fig.4 Pictures showing the Wulipai landslide in Danba County(A)正射影像; (B)SVF可视化影像; (C)山体阴影图;(D)滑坡后缘破碎岩体SVF影像;(E)滑坡后缘破碎岩体UAV影像; (F)滑坡特征图; (G)滑坡前缘垮塌现场照片; (H)前缘挡墙开裂照片

图5为丹巴县中路乡一处滑坡，该滑坡所处位置较高，最高位置的海拔高度>3.4 km，与坡脚河床的高差达1.5 km。该滑坡区位置植被发育，因而隐蔽性极强，光学影像上很难发现滑坡痕迹(图5-A)。在山体阴影图像上，整个斜坡区域受光照影响颜色较亮，未能清晰地反映实际的地形特征，图左下角的冲沟及右上角的山脊在山体阴影图上均较难辨识(图5-C黑色箭头所示)。在SVF影像上不仅山脊、冲沟这类特征较为明显的地貌可以清晰显示(图5-B)，滑坡后缘拉裂缝也可以被准确识别(图5-E)，从影像上可以看出受融雪、降雨等流水冲刷影响，滑坡区发育多条侵蚀细沟，且深浅不一，在SVF影像上颜色稍深的侵蚀沟表明其侵蚀深度相对较深。另外，滑坡体上多处次级滑动痕迹也非常明显(图5-E)，滑坡后缘崩塌(图5-D)和前缘一处相对较大的次级滑坡(图5-F)都较为清晰。上述特征在山体阴影图上的可辨识度均较低。

图5 丹巴县中路乡后山滑坡Fig.5 Pictures showing the Zhonglu landslide in Danba County(A)正射影像； (B)SVF可视化影像； (C)山体阴影图； (D)滑坡后缘局部崩塌SVF影像； (E)滑坡特征图； (G)滑坡次级滑动SVF影像

3 讨论

影响SVF地形可视化结果的主要参数是搜索方向的数目(n)、最大搜索半径(R)和DEM的空间分辨率。在空间分辨率一定的情况下，由于确定SVF的值需要大量的计算，在处理大型数据集时，搜索方向的数目和最大搜索半径的最优参数选择尤为重要。

3.1 搜索方向

利用空间分辨率为0.5 m、搜索半径为25 m(50像素)的DEM，对中路乡滑坡进行了可视化分析，分别对8、 16、 32和64个方向的SVF进行了目视定性和统计定量评估。可视化影像的视觉比较显示，在所有不同搜索方向的影像中都可清晰显示出滑坡后缘裂缝、细小冲沟等微小的地形特征；但是当选择8或64个方向时，在冲沟或裂缝的边缘上存在一些微小的差异。另外，图6-C中的散点图显示了在8和16个方向以及32和64个方向计算的SVF值，从图中可以看出，随着搜索方向数目的增加，SVF计算结果的差异不断减小，在32个和64个方向上计算的SVF值趋于一致，表明在更多方向(超过32)上计算SVF意义并不大，类似的，图6-B中的剖面图也显示了相同的结果。

图6 不同搜索方向下SVF可视化效果Fig.6 Visualization factor of SVF in different search directions (A)搜索半径为50像素，搜索方向分别为8、 16、 32、 64的SVF影像；(B)8、 16、 32、 64个搜索方向下的A-A’剖面图; (C)8、 16个搜索方向与32和64个搜索方向的SVF值散点图

因此，在利用DEM影像进行SVF可视化时，建议至少选择8个方向；但是为了提高计算效率，不建议大于32个搜索方向。在利用SVF进行辐射研究方面，J.Dozier等[27]研究者也得出了计算SVF的值时32个搜索方向就足够了的相似结论。

3.2 最大搜索半径

图7比较了当搜索方向为16时，搜索半径分别为5 m(10像素)、20 m(40像素)和50 m(100像素)时的SVF影像(与图6相同的区域)。从图7-A、B、C可以看出，较大的地形起伏特征在更大的搜索半径下才愈发明显；而当搜索半径增大到一定程度后，可视化影像之间的差异就逐渐减小。如图左下角的黄色箭头所示的冲沟，在100像素的搜索半径下比10像素的搜索半径下更明显；但是当搜索半径增大到100像素后，其与40像素的结果差异明显较小。

图7 不同搜索半径下SVF可视化效果图Fig.7 Visualization factor of SVF in different search radius (A)-(C)搜索方向为16时，搜索半径分别为10、 50、 100像素的SVF影像; (D)-(F)分别为搜索半径下的局部放大图; (G)10、 50、 100像素下的B-B’剖面图

而对于微小的地形特征，随着搜索半径增大，视觉上可见的特征愈发明显(图7-D、E、F中的箭头所示)；但类似的，搜索半径增大到40和100像素时，目视差异较小。图7-G中的剖面图也进一步印证了上述结果，且表明SVF主要在地形突变的地方发生急剧变化，如陡崖、裂缝、冲沟等。因此，对于类似数据，建议选用10～50像素的搜索半径。