基于机载激光雷达的地质灾害识别关键技术及应用研究

彭艺伟，董 琦，田 冲,陈 刚，孟小军,曾 斌*

(1.中国地质大学(武汉)环境学院，湖北 武汉 430078;2.浙江久核地质生态环境规划设计有限公司，浙江 湖州 313002;3.江苏省地质调查研究院，江苏 南京 210018)

机载激光雷达是遥感手段中的激光成像测量系统，以飞机(本次研究采用无人机)或车辆为载具对地面沿着航线进行垂直扫描，记录距离数据和其他信息，生成高精度的激光点云数据以及强度信息和光谱信息。国内学者利用激光雷达测量(LiDAR)技术开展了多个领域的应用研究，如：陈鼎等利用LiDAR技术对铁路边坡形变进行了监测，可为后期铁路边坡的维护和加固提供较为准确的资料，预防边坡坍塌事件；
陈洁基于机载LiDAR技术分析和总结了断裂构造的点云特征，利用断裂的线性延伸、高程异常和坡度异常的特点，运用合适的算法自动提取了断裂构造信息；
何柯璐等采用LiDAR技术对采煤沉陷盆地进行地面扫描和建模，提出了基于格网高程排序的类地面点云提取算法，构建了精度更高的矿区开采沉陷模型；
杜磊等基于LiDAR技术获取了三峡库区张家湾地区的点云数据，并通过对点云数据进行加工处理，成功地识别出滑坡体，并提取了滑坡体重要的参数。在国外很多学者基于机载LiDAR技术开展了针对滑坡、泥石流的地质灾害调查，如：Mezaal等基于LiDAR技术，提出了基于相关的特征选择(CFS)算法用来识别滑坡并判别滑坡类型；
Pourghasemi等通过人工神经网络算法，基于机载激光点云数据进行了滑坡易发性分析与预警；
Booth等利用LiDAR重复扫描获取了一处滑坡的变形区域，并分析了滑坡的发生机理。虽然LiDAR技术已经在铁路边坡监测、地质构造调查以及滑坡体识别等方面有了一定的研究，但是对于植被覆盖区的地质灾害或地质灾害隐患识别的研究仍有待进一步深入。

激光雷达点云具有部分穿透植被的多次回波功能，可以快速且精准地构建研究区的真实数字高程模型(DEM)，从而进行数字地形分析并识别出研究区内植被以下的地质灾害或地质灾害隐患。传统的激光雷达多以直升机载或者车载，其优点是单次扫描任务面积大，但由于激光雷达距离地面较远且行进速度大，得到的激光点云较为稀疏，通常只有2～3 pts/m，滤波处理后得到剥离植被后的地面点云则密度更小，难以达到对小型地质灾害进行识别的目的。

针对江苏省山区植被覆盖茂盛、地质灾害隐蔽性强的特点，本文采用无人机搭载激光雷达，在离地80 m左右的位置对研究区进行扫描，得出的地面点云密度可达400～500 pts/m，能够精准地识别出江苏省植被覆盖区小型地质灾害的一系列参数，可为后期防灾减灾工作提供重要依据。

1.1 机载LiDAR技术的优势与特点

一直以来，传统的光学遥感技术，如光学卫星、无人机倾斜摄影三维建模等，在地质灾害调查中发挥了较大的作用，但是这些技术方法不能穿透植被，无法快速获取植被以下的真实地面信息。而机载LiDAR是通过仪器主动发射的脉冲信号，可部分穿透植被到达地表，通过点云去噪、滤波，剔除植被，可快速构建高精度地形地貌，免地面控制点数据相对精度可达厘米级。

机载激光雷达测量(LiDAR)系统主要由激光雷达扫描仪、GNSS/INS定位测姿系统(POS)、存储控制单元三大部分组成。其中,激光扫描仪包括扫描仪、测距单元和控制单元，主要用于测量扫描仪和地物之间的距离；
GNSS与INS主要用于测定激光雷达在空中的位置及五轴飞行姿态；
存储控制单元用于快速存储激光雷达获取的各项原始数据。三者的有机结合就能根据载具的位置反算出激光点云的精确三维坐标。

机载激光雷达测量系统对地定位示意图见图1。若空间有一向量，其模为

S

，方向为(

φ

，

ω

，

κ

)，如能测出该向量起点

O

的坐标(

X

，

Y

，

Z

)，则另一端点

P

(

X

，

Y

，

Z

)可唯一确定。对机载LiDAR系统来说，起点

O

为激光雷达信号发射参考点，其坐标(

X

，

Y

，

Z

)可通过机载GPS的记录数据计算得到；
向量的模是由激光测距系统测定的机载激光测距仪参考中心到地面激光脚点间的距离，姿态参数(

φ

，

ω

，

κ

)由惯性测量装置(IMU)获得。此外，还必考虑各种系统误差，如相位中心的误差、安置角误差、各集成系统中心间的距离和偏心角等。

图1 机载激光雷达测量(LiDAR)系统对地定位示意图[2]Fig.1 Working principle diagram of LiDAR system

1.2 机载LiDAR数据处理技术流程

激光雷达数据测量作业的生产环节，主要包括航线设计、激光雷达数据采集、数据预处理、点云解算以及DEM、等高线和数字地形分析图件制作等，如图2所示。

图2 激光雷达(LiDAR)识别地质灾害技术流程图Fig.2 Technical flow chart of LiDAR identification of geological hazards

(1) 航线设计。本文采用在专业地图软件平台上勾画航线，并在飞行器遥控器上根据等高线设计飞行高度、飞行速度以及转弯半径的方法进行航线设计。设计飞行高度为相对高度80 m，航线起始端进行8字飞行，航带旁向重复率大于10%。

(2) 研究区激光雷达数据采集。本文采用大疆经纬M300 RTK无人机搭载GS-MID 40激光雷达扫描系统获取研究区点云数据，同时辅以千寻星矩SR3专业级GNSS接收机(地面基站)获取当地的坐标和高程数据。

(3) 数据预处理。本文使用GS-MID 40激光雷达扫描系统适用的软件对原始点云数据进行预处理，得到通用格式点云数据。

(4) 点云解算。本文使用Mircrostation+TerraSolid软件中的TerraScan模块对点云数据进行处理。首先进行航带校正、数据检校消除误差；
其次对校正后的点云进行滤波、分类，得到滤除植被后的地面点云。

(5) DEM制作及后处理。本文在ArcGIS平台中使用相应的工具生成滤除植被后的DEM，并生成山体阴影、地表坡度、地表粗糙度和地形起伏度等数字地形分析图件，最后对地质灾害体进行识别和分析。

2.1 研究区概况

本文以江苏省无锡市宜兴竹海风景区公路旁一处道路边坡为研究区，其地理坐标为东经119°42′35.39″、北纬31°09′54.30″，区内堆积体较厚，下伏基岩为泥盆纪观山组石英砂岩，道路位于半山腰处，南侧已发生过滑坡，目前已经得到治理。但由于植被覆盖密度大，传统的光学卫星、倾斜摄影等遥感手段都无法识别植被覆盖层以下的地质灾害。本次采用激光雷达探测发现研究区新滑坡体，其位于道路北侧植被茂密区，见图3。

图3 无人机航拍研究区图(由于植被茂密，普通光学 遥感无法发现滑坡体)Fig.3 Aerial photo of the study area taken by UAV (Due to dense vegetation,ordinary optical remote sensing cannot find landslide)

2.2 研究方法

首先利用无人机机载LiDAR低空测取研究区高密度激光点云数据；
然后通过Microstation+Terrasolid软件，对点云数据进行滤波处理，得到剥离植被后的地面点云；
再基于地面点云，在ArcGIS中生成研究区的DEM，并由DEM二次生成等高线图、多视角山体阴影图、地形坡度图、地表粗糙度图和地形起伏度图等数字地形分析图件；
最后利用点云剖面和数字地形分析图，完成植被覆盖区地质灾害体的识别和分析。

3.1 数据获取

本次研究利用机载激光雷达系统和相应的地面基站，获取了研究区的遥感数据，包括：激光雷达数据，即GNSS、INS、测距扫描角、强度时间；
基站数据，即坐标、绝对高程和姿态数据。由于研究区位于低山丘陵区，海拔在150～240 m之间，区域内大部分面积为竹林所覆盖，故综合考虑地貌、气象、飞行安全及应用需求等因素，飞行高度设定为240～320 m之间，飞行速度设定为4 m/s，旁向重复率高于10%，数据量为100 000 pts/s，实际航测面积为66 208 m，激光点云密度为487.3 pts/ m，远高于以往研究所采用的有人机机载LiDAR。

3.2 数据处理

3.2.1 数据预处理

获取的研究区LiDAR数据并不能直接地进行解算，需要先在与GS-MID 40激光雷达扫描系统配套的数据处理软件中将原始点云数据转换为可以在TerraSolid软件中适用的格式，即Las格式点云文件。首先在配套的数据处理软件中将“.dat”格式的流文件分拣为解算所需的文件；
然后新建工程文件，导入相关文件，输入基站信息，差分后显示差分定位结果，输出GNSS文件；
再使用输出的GNSS文件并添加INS文件，处理得到GNSS/INS组合结果POS文件(见图4)；
最后采用配套软件处理获得的激光数据和GNSS/INS组合结果POS文件进行点云计算，并输出Las格式点云文件(见图5)。

图4 GNSS/INS组合结果POS文件Fig.4 GNSS & INS fusion results

图5 激光数据预处理结果Fig.5 Laser data preprocessing results

3.2.2 数据后处理

3.2.2.1 点云数据处理

(1) 点云数据噪声去除。在TerraSolid软件的TerraScan模块中将预处理得到的Las格式点云数据进行后续的点云解算，具体点云解算及后续的DEM图件加工流程见图2。由于受仪器误差和偶然误差的影响，获取的点云数据坐标会产生偏移，因此在滤波分类前需进行航带校正和数据检校，以消除误差，提高数据精度。然后裁剪掉与所需滑坡信息无关的点云，即缩小范围，避免在解算过程中由周围临近点云引起的精度不高的问题。

(2) 点云滤波获取地面点。在精确点云范围后，开始进行滤波分类。首先对点云数据采用手动编辑的方式去除明显的噪点和孤点；
然后采用点云滤波处理剔除由于植被影响的多次回波现象形成的不规则形状和因激光脉冲的折射、多路径效应等引起的点云

x

、

y

、

z

值异常而产生的噪点；
最后通过软件自动分类算法，将点云数据分成若干类别。由于滤波后的目标产品是DEM，故只需要得到准确的地面点云数据，因此将点云数据自动分成地面、低点、重叠点等几类。TerraScan模块地面点的分类是基于Axelsson提出的不规则三角网滤波算法实现的，可以得到高精度的地面点。该算法以局部高程最低点作为种子点构建TIN网格，通过计算待定点到三角形各节点的坡度值和到三角形平面的距离是否处于设定好的迭代阈值(见图6)，来判断该点是否为地面点，若是地面点则利用该点重新构建TIN网格，并迭代加入地面点直至没有新的地面点被加入，结束滤波。

图6 迭代阈值Fig.6 Iterative threshold

此外，根据渐进三角网滤波分类算法自动分出的地面点存在一定的误差，需要进行人机交互编辑，进行手动分类，将不正确的分类通过人工干预的方法得到正确的地面点。经过滤波处理滤除植被后生成的地面点云的密度可以达到10～40 pts/m，满足地质灾害解译的需求。

3.2.2.2 DEM的生成

数字高程模型(DEM)是表现地表形态最常用的方式之一。本次研究采用易于分析与计算的规则格网DEM来表达斜坡的地表信息。生成DEM的过程是利用ArcGIS转换工具中的Las数据集转栅格工具对Las数据进行内插,得出规则格网DEM。内插方法包括反距离加权法、克里金插值法和自然领域插值法3种，其中反距离加权法插值的精度最高，故本次采用反距离加权法内插生成高精度的DEM，作为等高线、山体阴影、地表坡度、地表粗糙度等辅助地质灾害解译后续图件加工的基础。

4.1 基于地面形态要素的滑坡识别原理

利用传统光学遥感工具进行滑坡识别，是基于滑坡在遥感图像中的梨形、舌形、簸箕形等典型坡面形态进行的。在遥感解译中，根据这些图形特征便可以在遥感影像中大致圈定滑坡范围。利用机载LiDAR测得的点云数据具有厘米级精度，同时拥有滤除地表植被的功能，能够在植被覆盖的区域分离出地面点云。基于机载LiDAR的滑坡识别，除了能够准确地识别传统光学遥感工具所能识别的典型地貌特征外，还增加了对滑坡微地貌的识别，能够清晰地识别其滑坡舌、滑坡壁、滑坡台阶和滑坡裂缝等地貌形态，达到精确判别滑坡类型和识别滑坡边界的效果。

通常情况下滑坡具有如图7所示的全部或部分空间形态要素，遥感解译中，则主要是基于这些形态要素的识别来确认滑坡。

图7 滑坡及其空间形态要素示意图Fig.7 Schematic diagram of landslide and its spatial form elements

4.2 基于山体阴影图和数字地形分析图的滑坡识别

利用遥感技术对地质灾害进行全面的调查，是地面调查的有效补充和辅助工具，而利用传统光学遥感卫星图像进行的遥感解译，其结果因受到植被覆盖的影响和图像精度的限制会存在较大的误差，在这种情况下，本文引入LiDAR技术对地质灾害进行识别。

4.2.1 基于山体阴影图的滑坡范围识别

采用LiDAR 数据生成的高精度DEM(见图8)和利用ArcGIS软件的3DAnalyst扩展模块生成研究区一系列不同太阳方位角的山体阴影(Hillshade)，作为地质灾害隐患识别的定性参数。

图8 研究区DEM图Fig.8 DEM map of the study area

本次研究的山体阴影图是利用方位角分别为270°、315°和360°，高度角为45°的假想光源进行相邻像元的照明值计算得到的，可较好地显示出研究区地形的立体形态，有利于提取研究区滑坡等地形信息。图9为太阳方位角为360°时研究区的山体阴影图，与实际航拍照片(见图3)相比，山体阴影图可以还原出真实地表。在山体阴影图(见图9)中，由于滑坡与周边环境的灰度值呈现明显的差异，并且其圈椅状边界也符合滑坡发育形态。

图9 研究区山体阴影图(太阳方位角为360°)Fig.9 Mountain shadow map of the study area (solar azimuth 360°)

4.2.2 数字地形分析

数字地形分析(Digital Terrain Analysis，简称DTA)是在DEM上进行地形属性计算和特征提取的数字信息处理技术，是对DEM应用范围的拓展和延伸。本次研究采用地表坡度图、地表粗糙度图和地形起伏度图作为分析对象。

地表坡度是指地表面上某一点的切面和水平面所成的夹角。地表坡度反映了地势的陡峭程度，在地形表面分析中起到重要的作用。本次研究采用2D笛卡尔坐标系对投影平面计算地表坡度。

地表粗糙度是指一定范围内地表表面积

S

与其在水平方向上的投影面积

S

之比。地表粗糙度反映了地表形态的宏观变化，是反映地表起伏变化和侵蚀程度的一个重要指标。其数学表达式为

式中:

R

为地表粗糙度(无量纲)；
angle表示地表面与水平面的坡度(°)。

地形起伏度是指某一确定面积内最高点与最低点的高程之差。地形起伏度是反映地形起伏的宏观地形因子，在区域研究中，能够定量地描述地貌形态，是划分地貌类型的一项指标。其数学表达公式为

RA

=

Z

max-

Z

min式中:

RA

为某区域的地形起伏度(m)；

Z

max为该区域内的最大高程值(m)；

Z

min为该区域内的最小高程值(m)；

i

为自然数。

研究区地表坡度图、地表粗糙度图和地形起伏度图，见图10至图12。

图10 研究区地表坡度图Fig.10 Surface slope map of the study area

由图10至图12中可以从不同角度识别出研究区滑坡体在平面上的空间形态及微地貌特征：①地表坡度值域在0～89°之间，在滑坡周界范围内地表坡度有明显的数值升高，表明该局部区块地表坡度因滑坡滑动后发生变化而与周边稳定斜坡体有所区别；
②地表粗糙度值域在1～58之间，在滑坡后缘部分地表范围粗糙度明显增加，表明这个局部区块因滑坡滑动过程中岩土体的位移及拉裂破碎而造成了比面积的增加；
③地形起伏度值域在0～6 m之间，而在滑坡周边则地形起伏度值域为3～6 m之间，表明滑坡滑动后地形起伏程度突变，滑坡后缘变化明显。

由于研究区域范围较小，数字地形分析图整体数值变化不大，仅靠目视难以发现其细微变化，但在滑坡区域有较为明显的变化，即滑坡后缘的地表坡度值、地表粗糙度值和地形起伏度值发生了明显的跳跃性变化，跳跃幅度基本相同，因此从数字地形分析图中能够直接圈定滑坡后缘边界。由图11可以看出，滑坡坡脚处地表粗糙度值也明显呈跳跃变化，由此可以判定滑坡前缘。同时，从图10至图12中图像颜色来看，滑坡后缘和滑坡前缘颜色与临近区域有明显的差异，可以圈定出滑坡的范围。

图11 研究区地表粗糙度图Fig.11 Surface roughness diagram of the study area

图12 研究区地形起伏度图Fig.12 Topographic relief diagram of the study area

综上数字地形分析图可以判断出该滑坡面积为930 m，宽为19 m、长为91 m，滑坡整体地表坡度为35°，滑坡体主滑方向为南偏西6°。

4.3 基于点云地形剖面图的滑坡细节信息分析

基于高密度点云的地形剖面图能直观地表现出地面的起伏、地势的变化和地表坡度的陡缓，还能形象地显示出一个地区的地形类型及其特征。由4.2.2节综合数字地形分析图件，可实现对滑坡的初步识别与判定(见图13),滑坡地形剖面图见图14至图17。

图13 滑坡解译结果及结构要素图Fig.13 Landslide interpretation results and structural elements

图14 滑坡地形剖面A—A′图(未剥离植被)Fig.14 A—A′ profile (vegetation unremoved)

图15 滑坡地形剖面A—A′图(剥离植被后)Fig.15 A—A′ profile(vegetation removed)

图16 滑坡地形剖面B—B′图(剥离植被后)Fig.16 B—B′ profile(vegetation removed)

图17 滑坡地形剖面C—C′图(剥离植被后)Fig.17 C—C′ profile(vegetation removed)

由于斜坡地表被植被覆盖，地面的微地貌特征信息无法识别(见图14)。而剥离植被后，可根据滑坡地形剖面进一步刻画滑坡体的详细特征。如：滑坡地形剖面

A

—

A

′(剥离植被后)上，滑坡后缘距公路的高差为8.7 m，滑坡后缘距滑坡前缘坡脚的水平距离为11.7 m，滑坡后缘可见明显的陡坎(高差约3.1 m，坡度约68°)，滑坡中部可见滑坡台阶(宽约3.5 m)，滑坡前缘地表坡度约30°；
滑坡地形剖面

B

—

B

′(剥离植被后)上，滑坡后缘距公路的高差为10.1 m，滑坡后缘距滑坡前缘坡脚的水平距离为16.5 m，滑坡后缘可见滑坡壁(高差约1.9 m)，滑坡中部可见滑坡台阶(宽约6.2 m)，滑坡前缘坡度约36°，该剖面为主滑剖面；
滑坡地形剖面

C

—

C

′(剥离植被后)上，滑坡后缘距公路的高差为7.9 m，滑坡后缘距滑坡前缘坡脚的水平距离为9.5 m，滑坡后缘可见滑坡壁(高差约2.4 m)，滑坡中部未见明显的滑坡台阶，滑坡总体地表坡度约为62°。

在Terrasolid中对滤波后的点云模型和地形剖面进行分析可知：公路边内侧斜坡上存在明显的呈圈椅状的滑坡，主滑方向正对道路；
斜坡体上存在滑坡壁、滑坡台阶等较为明显的滑坡要素，滑坡后壁高差最大为3.1 m，滑坡台阶最宽为6.2 m；
在滑坡后缘暂未发现明显的由滑坡引起的如张拉裂缝等变形迹象。但由于滑坡整体地表坡度、地形起伏度较大，故在人类工程活动以及降雨、风化等自然因素的影响下，该滑坡存在继续发生形变的可能性。

本文运用无人机机载LiDAR技术对江苏省宜兴市竹海风景区公路旁一处道路边坡66 208 m的研究区进行了高密度点云测取，最终测得原始点云密度为487.3 pts/ m，滤波处理剥离植被后的地面点云密度为13 pts/m，远高于传统有人机机载LiDAR。通过对点云数据的处理，获取了研究区内剥离植被后的DEM图、山体阴影图、地形起伏度图、地表坡度图和地表粗糙度图，并通过对地面点云多视角的观测和滑坡地形剖面图测量，结合山体阴影图和数字地形图分析，圈定了江苏省宜兴市竹海风景区公路旁植被茂盛区一处滑坡的范围，获取了滑坡的相关基本信息：滑坡体平面上呈圈椅状，面积为930 m，长为91 m、宽为19 m，整体地表坡度为35°，滑坡体主滑方向为南偏西6°。另根据滑坡地形剖面图可以看出，滑坡体的滑坡前、后缘高差为10.1 m，滑坡后壁高差为1.9～3.1m，滑坡台阶宽度为3.5～6.2 m。

在以往的地质灾害调查过程中，传统的光学遥感技术以及无人机低空航拍遥感不能穿透植被，不能获取植被覆盖区以下真实的地表(形变)信息，且数据的几何精度受地形地貌、传感器类型、天气因素等的影响较大。而机载LiDAR技术是利用激光仪器主动发射脉冲信号，可部分穿透植被到达地面，通过点云数据去噪、滤波，剔除植被点云层数据后，快速构建与重现高精度的地面形态，并提供地面变形迹象信息，且在无地面控制点情况下数据的相对精度可达厘米级。

本文利用无人机机载激光雷达可滤除植被的二次回波特点以及其点云高密度的特性，在低山丘陵植被覆盖区进行斜坡体变形迹象的地面信息读取，从而达到对植被茂密区地质灾害隐患早期识别的目的，极具实用性和创新性，能够为地质灾害的地面详查提供极为有效的技术支撑。

猜你喜欢 坡度激光雷达植被 基于双轴加速度的车辆坡度优化算法研究汽车电器(2022年6期)2022-07-02激光雷达实时提取甘蔗垄间导航线农业工程学报(2022年4期)2022-04-24法雷奥第二代SCALA?激光雷达汽车观察(2021年8期)2021-09-01追踪盗猎者课外生活·趣知识(2021年2期)2021-05-24融合激光雷达与超声波数据的障碍物检测方法科技研究·理论版(2021年20期)2021-04-20第一节 主要植被与自然环境 教学设计学校教育研究(2021年24期)2021-03-28Ouster发布首款全固态数字激光雷达计算机与网络(2020年19期)2020-12-04Aqueducts考试与评价·高二版(2020年4期)2020-09-10杭锦旗植被遥感分析农业与技术(2020年12期)2020-07-04放缓坡度 因势利导 激发潜能——第二学段自主习作教学的有效尝试福建基础教育研究(2020年3期)2020-05-28