机载激光雷达技术在滑坡调查中的应用
——以三峡库区张家湾滑坡为例

杜磊， 陈洁，2， 李敏敏， 郑雄伟， 李京， 高子弘

(1.中国自然资源航空物探遥感中心，北京 100083； 2.中国科学院遥感与数字地球研究所，北京 100101)

0 引言

滑坡作为一种全球范围内普遍存在而且频繁发生的地质灾害现象，对其规律性进行研究和预警一直是地学科研热点之一。目前滑坡调查的手段主要有野外实地勘查和遥感调查2种，前者作为传统的工作方法，可近距离现场研究滑坡形态，圈定滑坡要素，但其工作耗时久、效率低，抵达调查区域也受限较多； 卫星遥感技术支持下的滑坡调查能快速获取大尺度、多时相的数字遥感数据，现已形成一套较为完善的“数字滑坡”识别技术[1]。但该技术的不足之处在于受天气影响较大、且对被植被覆盖的滑坡识别能力有限，源数据生成的数字高程模型(digital elevation model，DEM)精度较低，难以刻画地表微小变形的精细特征。

机载激光雷达(light detection and ranging,LiDAR)作为一种主动式航空遥感技术，理论上可不受天气限制，其发射的激光点云能够穿透一定厚度的植被并获取地表信息。与传统航空摄影相比，机载LiDAR数据采集方式能在很大程度上减弱地形切割阴影的影响，其构建的DEM垂直精度可达10 cm，可以精确表达微地貌地形特征[2]，具备提升滑坡体识别精度的条件。近年来，美国地质调查局及欧洲太空局等机构均开展了基于机载LiDAR衍生DEM的滑坡识别研究，取得了一定的进展[3-4],包括： 基于高精度DEM计算滑坡坡度、坡向、地表粗糙度、半方差和分形等特征； 结合傅里叶变换、数据小波分析等多种算法开展地表地形有关参数精细特征研究[5]，并探索建立滑坡等地质灾害分布与上述有关因子的联系。以上研究充分显示了机载LiDAR技术在滑坡调查中的应用潜力。

本文基于机载LiDAR技术获取了研究区点云数据，通过对点云数据先后开展检校、滤波和分类等数据处理，构建研究区精细的DEM产品，再利用这些DEM产品对研究区的滑坡体进行识别以及重要参数提取，取得了较好的应用效果。

1 研究区概况和数据源

三峡库区一直以来是我国滑坡等地质灾害的频发区和重灾区。张家湾滑坡是库区较为典型的一处滑坡群。它由6个倾向长江的单体滑坡组成，位于童庄河右岸，隶属湖北省宜昌市秭归县郭家坝镇头道河村一组，与入江口相距约1.5 km。滑坡体总体前缘高程约为120 m，后缘高程约290 m，水位以上总面积约29×104m2，对应总体积约550×104m3。滑坡区地势整体东高西低，后侧山顶高程约500 m，滑坡前缘为三峡库区水面，中部及前缘地带坡面相对平缓，坡度为15°～25°，后缘及后侧山坡较陡峭，坡度为25°～35°。滑坡区内出露的地层主要为侏罗系砂页岩及三叠系中统巴东组(T2b)砂岩、泥灰岩等半坚硬岩类组成。

选择2009年由中国自然资源航空物探遥感中心获取的巴东—秭归地质灾害易发多发重点区1 000 km2LiDAR数据作为遥感数据源，其真实记录了三峡库区自移民迁建、175 m水位蓄水和2008年汶川震后的坝体微变化，具有开展库区库岸稳定性、地质灾害和生态地质环境变化调查等研究的潜力。

2 技术流程

针对现有滑坡调查技术在复杂地质条件下难以识别滑坡，以及传统基于DEM数据识别滑坡方法存在的问题，本文提出了一种从机载LiDAR构建的精细DEM中表征滑坡地表特征参数的技术方法，实现了滑坡边界的自动识别和相应参数提取。新方法步骤如下：

1)针对目标区域制定作业计划，使用定位定向系统(position orientation system,POS)辅助机载LiDAR系统对目标区域开展航空遥感数据采集。

2)对飞行获取的机载LiDAR点云数据先后进行检校、滤波和分类等数据处理，得到消除植被覆盖信息的纯地表精细DEM。

3)从DEM中分析、提取出地貌特征并建立纹理特征参数，形成特征参数文件。

4)利用计算机智能分类算法计算目标区域中已知滑坡像元以及非滑坡像元特征参数，结合上一步获取的特征参数文件，确定最优特征参数。

5)选取目标区域内的部分已知滑坡像元和非滑坡像元作为算法训练集，不断优化训练集各元素，综合确定最优特征参数组合。利用智能分类算法预测训练集中的像元属性同时进行精度评价，获得符合预设精度条件的平衡系数。

6)利用获取平衡系数的训练集以及最优特征参数组合，训练分类模型，预测已知滑坡像元和非滑坡像元的数据集，并分别计算数据集的平均用户精度、平均生产者精度和总体精度。

7)若所计算的精度满足设定阈值要求，则使用满足该平衡系数和最优特征参数组合的智能分类模型预测整个研究区，再采用边缘检测算子计算滑坡边界，最终实现滑坡识别。

8)对识别出的滑坡像元集，通过对应的精细DEM提取出滑坡坡度和地表粗糙度等滑坡参数，并完成定量分析。

3 数据处理与DEM制作

3.1 数据获取

2009年4—6月，使用运5-B型飞机作为航空遥感平台，搭载Leica Geosystem公司制造的ALS50-II机载LiDAR系统，获取了包含张家湾滑坡在内约1 000 km2研究区综合遥感数据资料，包括机载LiDAR点云数据、数字航空影像、机载GPS数据、机载IMU数据和地面GPS参考站定位数据等。已获取的机载LiDAR点云数据包含激光距离、角度和多次回波强度等重要信息。通过对获取的原始数据进行数据解算、系统误差检校等预处理，得到了含地物三维坐标信息的地表激光点云数据。

此次ALS50-II机载LiDAR系统飞行相对航线高度约为1 000 m，航线旁向重叠度为30%～39%，飞行地速为150～180 km/h，机载激光空中扫描角度为45°，其最大扫描频率设置为45 Hz，机载LiDAR平均点云密度为3.2个/m2。

3.2 点云解算

受系统误差和偶然误差影响，机载LiDAR系统获取的原始LiDAR点云坐标存在系统偏差，需先通过数据检校消除误差。机载LiDAR点云数据检校一般指检校侧滚角、俯仰角和航偏角3个安置角和测距误差。

3个安置角的检校方法如下： 首先,对检校场的数据开展人工初检校，分别获取较为准确的3个安置角检校参数； 然后,对初检校的机载LiDAR数据开展进一步精检校，获取检校参数的改正值; 最后，将上述2次检校参数相加得到安置角最终的检校参数。

完成安置角检校后，利用实测GPS点对测距误差进行检校，获取测距检校参数。利用上述安置角和测距检校参数，对研究区全部机载LiDAR点云数据进行检校，以消除研究区点云数据系统误差。经验证，研究区域各航带点云匹配效果较好，高程方向相对高差优于10 cm，平面位置相对差优于25 cm，结果稳定可靠。

对于数据检校解算依然无法较好消除的残余航带系统误差，本文采用目前已相对较为成熟的Burman航带平差方法进行改善，通过TerraSolid商业软件系列中的TerraMatch模块实现。

生成DEM前，先对机载LiDAR点云数据进行滤波分类处理，以区分地面点和非地面点，再利用地面点数据制作所需的高精度DEM。LiDAR点云分类处理采用了美国Bentley公司TerraSolid软件中TerraScan和TerraModeler等模块。其中，TerraScan模块地面点分类是基于Axelsson提出的不规则三角格网(triangulated irregular network，TIN)加密改进算法实现[6]。TerraScan模块提供了地面点分类的Ground工具命令，该方法利用迭代算法建立TIN表面模型来区分地面点和非地面点。

开展数据迭代运算前需首先选定若干初始地面点，通过Max building size参数调节初始点的选择。实验表明： 当Max building size取值30 m时，不仅可以很好地去除植被和房屋并保留微地形特征，而且还可以减少手工编辑的工作量。

通过分类前后对比(图1)可以看出，经滤波分类后的DEM中大部分建筑物和植被信息被自动滤除。少量建筑物和植被残余信息则用人工编辑的方法予以剔除。

(a) 分类前(b) 分类后

图1分类前后DEM渲染对比

Fig.1ComparisonofDEMshadingmapsbeforeandaftercloud-pointclassification

3.3 DEM制作及评价

先利用机载LiDAR系统获取的海量LiDAR点云数据快速开展目标区域地面三维数字地形模型(digital terrain model，DTM)采集，再对原始点云数据进行系统误差纠正和点云数据分类，有效滤除了点云中的非地面点，进而获得目标区域高精度DEM。为了评价DEM产品的精度，在研究区内均匀布设30个地面检查点，通过DEM读取记录地面检查点的三维坐标后，再和实测同名点三维坐标进行对比计算，完成研究区DEM的精度质量分析(表1)。结果表明，生成的DEM高程中误差满足1: 1 000比例尺规范精度要求[7]，可应用于滑坡识别调查。

表1 DEM成果精度评价Tab.1 DEM precision evaluation (m)

4 滑坡调查应用与分析

4.1 解译标志建立

滑坡的解译标志主要为“箕”状色调带和纹理异常带。滑坡体一般位于地质较稳定的自然斜坡凸突的负地形中，滑坡后壁与滑体交接处往往形成洼地，中部则有多级垂直滑动方向的台坎。图2是张家湾滑坡的数字正射影像图(digital orthophoto map, DOM)，它是通过机载LiDAR同步获取的航空遥感影像制作而成。由图2可以看出，张家湾滑坡是一倾向长江的单斜顺向滑坡群，影像中滑体一般呈浅色调，通常具有一定的挤压、扰动或松脱等特征，岩(土)体破碎，结构松散，耕地、园地及林地被扰动或者破坏比较严重，滑体纹理与背景呈突变关系。

图2 张家湾滑坡区域数字正射影像图Fig.2 DOM of Zhangjiawan Village landslide area

本文的野外工作包括前期野外踏勘和解译成果野外验证2部分。前期野外踏勘主要是对研究内容进行地面初步感性认识，结合相应遥感影像特征，建立遥感解译标志； 解译成果野外验证则是对室内目视解译结果在实地开展查验，并对解译标志做进一步修正，从而达到理性认识，使遥感解译成果与实际情况更加吻合。

4.2 技术方法

通过机载LiDAR点云数据制作的张家湾滑坡区域DEM可精细表达该地区微地形地貌信息，本文通过对其进行定性及定量分析完成了滑坡调查有关应用研究。

定性分析是通过张家湾DEM制作不同视角山体阴影图完成。山体阴影图是通过模拟光线计算DEM每个栅格单元的照明值。通过对比分析不同角度山体阴影图，不仅可较好地掌握张家湾滑坡区域的立体形态，还可以有效地提取出某些地形遮蔽信息，如滑坡、断层崖等线性特征。制作张家湾滑坡区域山体阴影图需确定3个重要参数，分别是： 太阳方位角、太阳高度角和表面灰度值。考虑到研究区内干流及支流走向，选择以相同的太阳高度角(45°)，用180°,135°和90°这3个不同的方位角，制作山体阴影图，表面灰度值设定为0～255。

定量方法则是指通过目标区域的DEM数据提取精细微地形地貌参数，分析滑坡坡度、地表粗糙度、半方差和分维等滑坡要素信息。其中，地表粗糙度是指在某确定范围内地表表面积与其在水平方向上的投影面积之比，是反映地形起伏和被侵蚀程度的一个重要指标，可有效反映地表形态的宏观变化。半方差是指滑坡体相对于周围地形表现在DEM数据值上的突变，它反映地表性质的不同距离观测值之间的变化。当DEM设置为半方差变量时，半方差函数则能够衡量地形的空间自相关性。相同尺度下，半方差值越大，对应空间地形变化越大，滑坡体与非滑坡体识别则越容易。分形是指部分与整体以某种形式相似的形； 分维是指分形的维数，它不仅可以衡量地形粗糙度，还可以体现地形的复杂度和空间自相关性。相同条件下，目标区域分维数值变化越大，其对应地形越复杂，地表则越粗糙[3]。

4.3 滑坡识别与定量分析

通过张家湾地区高精度DEM制作不同视角下的山体阴影图，结合张家湾滑坡的坡度图和地表粗糙度图完成滑坡识别； 通过计算张家湾滑坡体半方差函数和分维完成对滑坡的定量分析。图3(a)是张家湾滑坡野外实地勘查照片，红线是其中2处滑坡的边界线； 图3(b)是张家湾滑坡DEM彩色晕渲图，通过对其进行数值分析可了解滑坡地貌形态特征。

(a) 野外实地勘查照片(b) 滑坡DEM晕渲图

图3张家湾滑坡野外实地勘查照片及其滑坡DEM晕渲图

Fig.3LandslidefieldphotoanditsDEMshadingmapofZhangjiawanVillage

图4是张家湾滑坡2009年山体阴影系列图，图4(a)，(b)，(c)分别是太阳方位角设置为90°，135°和180°时制作出的山体阴影图(太阳高度角均设置为45°，表面灰度值设定为0～255)。对比分析3幅山体阴影图可知： ①与DEM数据栅格图像相比，3幅山体阴影图均可较好地反映滑坡地形的立体形态，对冲沟等线性地物特征的表现效果尤为突出； ②对滑坡体局部地形地貌特征的表现，3幅山体阴影图则差异明显，比如： 滑坡冲沟等线性地物在太阳方位角为90°和180°的山体阴影图上明显比135°的山体阴影图要清晰可见； ③太阳方位角为90°和180°的山体阴影图反映出的局部地貌特征相似但又不完全相同，但2幅影像信息呈互补关系。综上分析，通过高精度DEM数据制作的不同太阳方位角山体阴影系列图，不仅能有效地识别出张家湾6个单体滑坡的滑动范围，还可较为准确地确定其滑坡后缘、侧缘和滑舌等滑坡要素信息。

(a) 方位角为90°(b) 方位角为135°(c) 方位角为180°

图4张家湾滑坡山体阴影系列图

Fig.4SlidemountainshadowseriesmapsofZhangjiawanVillage

图5是张家湾滑坡的坡度图和地表粗糙度图。由张家湾滑坡坡度图分析可知： ①滑坡侧缘或滑坡两侧坡度跳跃变化明显，且跳跃变化幅度大体一致，可在坡度图圈定出滑坡侧缘边界(图5(a)中红线所示)； ②滑坡壁与滑舌的坡度值较为集中，两者色彩较均匀，且在形态上表现为面状分布，亦可较好地将6个不同的滑坡体区分开； ③在滑坡体与滑舌、滑坡壁与滑坡体的分界处，滑坡侧缘及其附近地形地貌处，坡度变化明显并出现跳跃式变大或缩小，且跳跃变化幅度大致相同。综上，通过滑坡坡度图也可以有效圈定出张家湾滑坡的3大要素。此外，结合图4和图5可知，在张家湾滑坡群上方，地形轮廓清晰且呈封闭状，“双沟同源”现象比较明显，表明此处具有诱发滑坡地质灾害的地形特征。结合山体阴影系列图，可圈出滑坡隐患区域(图5中蓝线所示区域)。

(a) 滑坡坡度图(b) 地表粗糙度图

图5张家湾滑坡坡度图和地表粗糙度图

Fig.5LandslideslopemapandlandslidesurfaceroughnessmapofZhangjiawanVillage

对张家湾地区地表粗糙度图进行解译，结果显示： ①滑坡侧缘地表粗糙度与两侧地表粗糙度存在明显跳跃，地表粗糙度跳跃幅度也基本相同，据此从地表粗糙度图上可圈定出张家湾地区6个不同滑坡体滑坡侧缘的边界(图5(b)中红线所示)； ②对于张家湾地区滑坡图像灰阶分布，其滑坡体后部与滑坡壁的地表粗糙度值较集中，色彩也较匀滑，形态上表现为面状分布，据此亦可以将它们与滑坡体其他部分区别开； ③在滑坡壁与滑坡体、滑坡侧缘与其附近地貌处，地表粗糙度明显出现跳跃式变大或缩小，且跳跃变化幅度大致相同。综上，通过地表粗糙度图亦可有效圈定出张家湾滑坡的3大要素。

分别选择滑坡后缘、滑坡体和滑舌3部分作为半方差函数计算范围。张家湾地区滑坡后缘、滑坡体和滑舌的半方差曲线见图6。

图6 张家湾滑坡要素半方差曲线Fig.6 Landslide elements semivariogram curve of Zhangjiawan Village

由图6可知： ①滑坡后缘、滑坡体和滑舌半方差曲线斜率存在明显区别，这是滑坡要素识别的数学理论基础； ②滑坡后缘、滑坡体和滑舌的半方差函数曲线随变程增加均变大，表明张家湾地区各滑坡要素的地形地貌相关性较强； ③滑坡后缘半方差值最大，表明滑坡后缘地形空间变化最大； 滑坡体半方差值次之，表明滑坡体地形空间变化幅度介于滑坡后缘与滑舌之间； 滑舌半方差值最小，表明滑舌部分对应的空间地形变化也最小。虽然3种滑坡要素的半方差函数曲线均随着变程的增大而增大，但也存在明显区别： 在相同变程条件下，滑坡后缘值最大，滑坡体次之，而滑舌最小。

张家湾3种滑坡要素分维函数的计算范围与半方差相同，其分维曲线与同区域2006年分维曲线对比可知： 滑坡后缘的分维值变化最大，表明滑坡后缘对应地形空间的变化最大； 滑坡体的分维值变化次之，表明滑坡体对应地形空间的变化幅度介于滑坡后缘与滑舌之间； 滑舌的分维值变化最小，说明滑舌空间地形变化最小(如图7所示)。滑坡要素分维值的计算结果与半方差计算结果具有较好的一致性。

(a) 滑坡后缘分维变化(b) 滑坡体分维变化(c) 滑舌分维变化

图7张家湾滑坡要素分维曲线对比

Fig.7ComparisonoflandslidefractaldimensioncalculationresultsofZhangjiawanVillage

5 结论

机载激光雷达技术作为一种新的航空遥感调查与监测手段，在地质灾害调查与监测领域已获得较广泛关注和应用。与卫星和航空影像相比，机载激光雷达具有一定的植被穿透能力，多次回波可以用来精确测量和表达微地貌特征信息，在雾气较大或阴雨天气时也能快速获取真实地形信息，具备提升滑坡体识别精度的条件。具体表现在2方面： ①与边坡的整体变化相关的特征，如存在可以被机载激光雷达点云数据派生的DEM识别的斜坡凹陷和突变； ②结合机载激光雷达点云数据派生的精细DEM，结合地表粗糙度变化，可以有效识别某些特殊地表特征，如地表内部结构异常、地裂隙或滑痕等。

本文基于机载激光雷达系统获取的航空遥感数据构建了精细的DEM，通过对DEM产品进行深入挖掘，进行综合对比分析，参考有关地质资料，结合一些必要的野外验证，有效地获得了三峡库区张家湾滑坡群的滑坡信息，应用效果良好，表明采用机载激光雷达技术对地质灾害进行调查与监测是一种直观、快速而经济的方法，尤其对于灾害应急反应和灾后救援具有重要意义。