CN102819023B - 基于LiDAR的复杂地质背景区滑坡识别的方法及系统 - Google Patents

Abstract

本发明公开了一种基于LiDAR的复杂地质背景区滑坡识别的方法及系统，其中方法包括以下步骤：S1、对LiDAR数据进行处理生成消除植被影响后裸地表的DEM数据；S2、提取传统的地貌特征参数以及计算纹理特征参数，生成特征参数文件；S3、确定最优特征参数组合；S4、获得满足预设精度条件的平衡系数；S5、计算平均用户精度、平均生产者精度和总体精度；S6、若满足精度要求，则使用边缘检测算子计算滑坡边界，实现滑坡识别。本发明的实现对数据需求少，可充分挖掘LiDAR-DEM数据在滑坡地形分析中的应用潜力；模型分类精度很高，能够实现复杂地质背景区滑坡边界的自动识别。

Description

基于LiDAR的复杂地质背景区滑坡识别的方法及系统

技术领域

本发明涉及地质环境遥感技术领域，尤其涉及一种基于LiDAR的复杂地质背景区滑坡识别的方法及系统。

背景技术

滑坡是当今世界除地震以外造成巨大经济损失的地质灾害。近些年来，在全球极端气候增多背景下，滑坡灾害事件频发，严重威胁人民群众的生命财产安全。滑坡预测预警是滑坡灾害防治的重要工作内容之一，而滑坡的精准识别又是滑坡预测预警的关键技术问题之一。因此，开展滑坡识别研究可以为滑坡预测预警及灾害防治工程提供技术支撑，具有重大的理论意义和工程价值。

目前，滑坡识别方法主要分为3类：（1）传统的野外实地调查方法。该方法能够近距离观察滑坡的基本地质和地貌条件，但是工作效率低，并且对于地形切割强烈、植被覆盖度高的复杂地质背景区，无法近距离到达进行实地调查；（2）光学遥感技术尤其是高空间分辨率遥感技术支持下的滑坡识别方法。该方法在复杂地质背景条件下，能够识别出已经发生滑动的滑坡几何属性，但是其成像特征决定了其无法穿透植被、难以获取山体阴影下的地表信息，所以难以识别出缓动型滑坡，并且数据采集受天气因素影响较大，因此限制了光学成像技术用于复杂地质背景条件下的滑坡识别；（3）干涉雷达测量技术。目前，干涉雷达卫星数据获取难度大，数据成本高，并且在复杂地质背景区受植被、大气等干扰因素太多，少有成功的案例。

LiDAR（Light Detection And Ranging，机载激光雷达）是近十年来飞速发展起来的遥感技术，通过多次回波能够穿透中等覆盖度以下的植被，获取大面积、高密度山体阴影下的地面点云数据，经过处理后能够生成消除植被影响的、裸地表1m级的DEM（Digital Elevation Model，数字高程模型），利用此DEM可提取出多种精细地形参数。因此，当面对复杂地质背景区滑坡识别研究时，机载LiDAR具有明显的技术优势。

但是，面对高分辨率LiDAR DEM，传统的DEM滑坡分析方法无法充分发挥LiDAR DEM的数据优势和技术优势。如何从LiDAR DEM中挖掘更多的、能够表征复杂地质背景区滑坡的地表特征参数，需要进一步研究。

发明内容

本发明要解决的技术问题在于针对现有技术中复杂地质背景区滑坡识别面临的困境，以及传统的基于DEM数据的滑坡识别方法存在的问题，提出一种数据需求少，充分挖掘了LiDAR DEM数据，实用性强，能够实现复杂地质背景区滑坡边界的自动识别的基于LiDAR数据的复杂地质背景区滑坡识别方法及系统。

本发明解决其技术问题所采用的技术方案是：

提供一种基于机载激光雷达数据的复杂地质背景区滑坡识别的方法，包括以下步骤：

S1、对LiDAR数据进行处理生成消除植被影响后裸地表的DEM数据；

S2、从所述DEM数据中提取传统的地貌特征参数以及根据所述地貌特征参数计算纹理特征参数，生成特征参数文件；

S3、根据选定的研究区的部分已知滑坡像元和非滑坡像元，利用智能分类算法对所述特征参数文件中的特征参数进行重要性计算，确定最优特征参数组合；

S4、选择所述研究区的部分已知滑坡像元、非滑坡像元的数据集作为训练集，并不断改变所述训练集中的元素，结合所述最优特征参数组合，利用智能分类算法预测已知滑坡像元、非滑坡像元的数据集并进行相关精度计算，获得满足预设精度条件的平衡系数，所述平衡系数为满足预设精度条件时训练集中的非滑坡像元数目与训练集的滑坡像元数目的比值；

S5、利用满足所述平衡系数的训练集以及所述最优特征参数组合，训练智能分类模型，预测已知滑坡像元、非滑坡像元的数据集，并计算平均用户精度、平均生产者精度和总体精度；

S6、若所计算的平均用户精度、平均生产者精度和总体精度满足精度要求，则利用满足所述平衡系数的训练集以及所述最优特征参数组合，训练智能分类模型，并预测整个研究区的滑坡像元、非滑坡像元数据集，使用边缘检测算子计算滑坡边界，实现滑坡识别。

本发明所述的基于机载激光雷达数据的复杂地质背景区滑坡识别的方法中，所述智能分类算法为随机森林算法，所述智能分类模型为随机森林模型。

本发明所述的基于机载激光雷达数据的复杂地质背景区滑坡识别的方法中，所述传统的地貌特征参数包括高程、坡度、坡向以及像元检测窗口内的高程、坡度、坡向的平均值和标准差。

本发明所述的基于机载激光雷达数据的复杂地质背景区滑坡识别的方法中，所述像元检测窗口为3×3像元检测窗口。

本发明所述的基于机载激光雷达数据的复杂地质背景区滑坡识别的方法中，所述纹理特征参数包括相关性、对比度、二阶矩、熵、同质性在四个纹理方向的平均值，以及基于坡向方向的相关性、对比度、二阶矩、熵、同质性。

本发明所述的基于机载激光雷达数据的复杂地质背景区滑坡识别的方法中，所述四个纹理方向包括西-东方向、西南-东北方向，南-北方向和东南-西北方向。

本发明所述的基于机载激光雷达数据的复杂地质背景区滑坡识别的方法中，步骤S3具体为：

选取已知滑坡像元、非滑坡像元数据集中所有的滑坡像元和等数量的非滑坡像元作为训练集，并以该训练集为已知数据利用智能分类算法进行计算，得到所述特征参数文件中所有特征参数重要性的值，并对其从高到低进行排序，将前x个特征参数作为最优特征参数组合，其中x为袋外误差最小时对应的特征参数的个数。

本发明所述的基于机载激光雷达数据的复杂地质背景区滑坡识别的方法中，步骤S4中获得所述平衡系数具体过程为：

根据选择的像元确定已知滑坡像元、非滑坡像元数据集中非滑坡像元与滑坡像元的比值m；

选取已知滑坡像元、非滑坡像元的数据集中20%的滑坡像元和一定数量的非滑坡像元作为训练集，训练集中的非滑坡像元数目为训练集的滑坡像元数目的k倍，k=1+0.1*n，n为整数；

n从0开始取值，n每次增加1，训练集随着n值的变化而不断变化，用该训练集训练智能分类模型，预测整个已知数据集，并计算平均用户精度和平均生产者精度，直至k小于或等于m；

以k值为横轴，平均用户精度、平均生产者精度为纵轴作曲线，获得曲线中平均用户精度和平均生产者精度之差绝对值最小时的k值，将此时的k值作为平衡系数。

本发明解决其技术问题所采用的另一技术方案是：

提供一种基于机载激光雷达数据的复杂地质背景区滑坡识别系统，包括：

LiDAR数据处理模块，用于对LiDAR数据进行处理生成消除植被影响后裸地表的DEM数据；

特征参数文件生成模块，用于从所述DEM数据中提取传统的地貌特征参数以及根据所述地貌特征参数计算纹理特征参数，生成特征参数文件；

最优特征参数组合确定模块，用于根据选定的研究区的部分已知滑坡像元和非滑坡像元，利用智能分类算法对所述特征参数文件中的特征参数进行重要性计算，确定最优特征参数组合；

平衡系数计算模块，用于选择所述研究区的部分已知滑坡像元、非滑坡像元的数据集作为训练集，并不断改变所述训练集中的元素，结合所述最优特征参数组合，利用智能分类算法预测已知滑坡像元、非滑坡像元的数据集并进行相关精度计算，获得满足预设精度条件的平衡系数，所述平衡系数为满足预设精度条件时训练集中的非滑坡像元数目与训练集的滑坡像元数目的比值；

精度计算模块，用于根据所述研究区的部分已知滑坡像元和非滑坡像元以及所述最优特征参数组合，训练智能分类模型，预测已知滑坡像元、非滑坡像元的数据集，并计算平均用户精度、平均生产者精度和总体精度；

滑坡识别模块，用于在所计算的平均用户精度、平均生产者精度和总体精度满足精度要求时，利用满足所述平衡系数的训练集以及所述最优特征参数组合，训练智能分类模型，并预测整个研究区的滑坡像元、非滑坡像元数据集，使用边缘检测算子计算滑坡边界，实现滑坡识别。

本发明所述的基于机载激光雷达数据的复杂地质背景区滑坡识别系统中，所述传统的地貌特征参数包括高程、坡度、坡向以及像元检测窗口内的高程、坡度、坡向的平均值和标准差；

所述纹理特征参数包括相关性、对比度、二阶矩、熵和同质性在四个纹理方向的平均值，以及基于坡向方向的相关性、对比度、二阶矩、熵和同质性。

本发明产生的有益效果是：本发明利用机载激光雷达技术能够获取米级、高密度地表点云数据、消除植被影响的DEM数据的特征，提取纹理特征参数，生成特征参数文件，再通过计算获取最优特征参数组合，利用最优特征参数组合和智能分类算法预测已知滑坡像元、非滑坡像元的数据集，并计算平均用户精度、平均生产者精度和总体精度；若所计算的平均用户精度、平均生产者精度和总体精度满足精度要求，则利用满足平衡系数的训练集以及最优特征参数组合，训练智能分类模型，并预测整个研究区的滑坡像元、非滑坡像元数据集，使用边缘检测算子计算滑坡边界，实现滑坡识别。本发明的实现对数据需求少，从DEM数据中提取传统的地貌特征参数以及从地貌特征参数中提取纹理特征参数，生成特征参数文件，可充分挖掘LiDAR-DEM数据在滑坡地形分析中的应用潜力；采用智能分类算法使得模型分类精度很高，实用性强，能够实现复杂地质背景区滑坡边界的自动识别。

附图说明

下面将结合附图及实施例对本发明作进一步说明，附图中：

图1为本发明实施例基于LiDAR的复杂地质背景区滑坡识别的方法的流程图；

图2为本发明较佳实施例基于LiDAR的复杂地质背景区滑坡识别的方法的具体流程图；

图3为本发明较佳实施例坡向与纹理方向的关系图；

图4为本发明实施例的灰度共生矩阵示意图；

图5为本发明实施例的特征参数选择图；

图6为本发明实施例不断改变所述训练集中的元素，利用随机森林算法预测已知滑坡像元、非滑坡像元的数据集时所计算的平均用户精度和平均生产者精度和k值的关系图；

图7为本发明实施例长江某地区的滑坡识别结果；

图8为本发明实施例基于机载激光雷达数据的复杂地质背景区滑坡识别系统结构示意图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。

本发明利用机载激光雷达技术能够获取米级、高密度地表点云数据、消除植被影响的DEM数据的特征，提出新的纹理特征参数，能够实现复杂地质背景区滑坡边界的自动识别，实用性强，可提升“数字滑坡”的定量化水平。

本发明实施例中利用长江某地区的基于LiDAR数据获取的DEM数据进行滑坡自动识别，如图1所示，本发明实施例基于机载激光雷达数据的复杂地质背景区滑坡识别的方法，包括以下步骤：

S1、对LiDAR（Light Detection And Ranging，机载激光雷达）数据进行处理生成消除植被影响后裸地表的DEM（Digital Elevation Model，数字高程模型）数据，可以简称为LiDAR-DEM数据；

S2、从DEM数据中提取传统的地貌特征参数以及根据地貌特征参数计算纹理特征参数，生成特征参数文件；

S3、根据选定的研究区的部分已知滑坡像元和非滑坡像元，利用智能分类算法从特征参数文件中确定最优特征参数组合；

S4、选择研究区的部分已知滑坡像元和非滑坡像元作为训练集，并不断改变训练集中的元素，获得满足预设精度条件的平衡系数，该平衡系数为满足预设条件时训练集中的非滑坡像元数目与训练集的滑坡像元数目的比值；

S5、根据研究区的部分已知滑坡像元和非滑坡像元以及最优特征参数组合，训练智能分类模型，预测已知滑坡像元、非滑坡像元的数据集，并计算平均用户精度、平均生产者精度和总体精度；

S6、若所计算的平均用户精度、平均生产者精度和总体精度满足精度要求，则利用满足平衡系数的训练集以及最优特征参数组合，训练智能分类模型，并预测整个研究区的滑坡像元、非滑坡像元数据集，使用边缘检测算子计算滑坡边界，实现滑坡识别。

在本发明实施例中，传统的地貌特征参数包括高程、坡度、坡向以及像元检测窗口内的高程、坡度、坡向的平均值和标准差。本发明实施例中采用的像元检测窗口为3×3像元检测窗口。在其他实施例中也可以选择如5×5、7×7像元检测窗口等。

在本发明实施例中，智能分类算法采用森林随机算法，相应地智能分类模型为随机森林模型。可以理解的是，也可以采用其他的智能分类算法，如神经元网络法、支持向量机法等。

本发明的一个实施例中，纹理特征参数包括相关性、对比度、二阶矩、熵和同质性在四个纹理方向的平均值，以及基于坡向方向的相关性、对比度、二阶矩、熵和同质性。其中，在本发明的一个实施例中，四个纹理方向包括西-东（W-E）方向、西南-东北（SW-NE）方向、南-北（S-N）方向和东南-西北（SE-NW）方向。

如图2所示，本发明较佳实施例中，步骤S2中从LiDAR-DEM数据中提取的地貌特征参数具体包括：

1）提取高程、坡度和坡向信息。

把高程、坡度、坡向图像每个栅格的值分别输入到三个矩阵中，即高程矩阵、坡度矩阵、坡向矩阵。图像中栅格的行列号与矩阵的行列号是一一对应的，即图像中最上面的一行为第一行，从上到下行数依次增大，图像的最左边的一列为第一列，从左到右列数依次增大。

2）提取3×3像元窗口内高程、坡度、坡向的平均值和标准差。

首先取3×3像元检测窗口，并新建六个矩阵，分别为高程平均值矩阵、高程标准差矩阵、坡度平均值矩阵、坡度标准差矩阵、坡向平均值矩阵、坡向标准差矩阵，且每个矩阵中的每个元素均预设为0；然后从高程图像的左上角开始，获取3×3像元检测窗口内的9个像元的高程值，计算9个像元值的平均值（Mea）和标准差（Stdv），计算公式如下：

$\mathrm{Mean}=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{P}_i---\left(1\right)$

$\mathrm{Stdv}=\sqrt{\frac{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{(P_i-\mathrm{Mean})}^2}{N-1}}---\left(2\right)$

其中P_i为检测窗口内像元的高程值，N为像元总数。

将计算出来的平均值和标准差赋值给检测窗口的中心像元，即分别赋给高度平均值矩阵和高度标准差矩阵的第二行第二列的元素；接着按从左到右、从上到下的顺序移动3×3像元检测窗口（每次移动一个像元），依次计算平均值和标准差，并赋值给高程平均值矩阵和高程标准差矩阵中相应的元素；最后按照上述方法分别计算坡度平均值矩阵、坡度标准差矩阵和坡向平均值矩阵、坡向标准差矩阵。

步骤S2中根据地貌特征参数计算纹理特征参数具体包括：

1）计算高程、坡度、坡向在四个纹理方向的平均相关性、对比度、二阶矩、熵、同质性等纹理特征参数，具体包括以下步骤：

首先取3×3像元检测窗口，并新建十五个矩阵，分别为平均的高程、坡度、坡向相关性、对比度、二阶矩、熵、同质性，且每个矩阵中的每个元素均预设为0。

然后从高程图像的左上角开始，获取3×3像元检测窗口内的像元值，并建立一个3×3的矩阵，通过该矩阵计算四个纹理方向（0°(W-E)、45°(SW-NE)、90°(S-N)和135°(SE-NW)）的灰度共生矩阵（GLCM），然后计算得到GLCM的对称、归一化矩阵。

GLCM及其对称、归一化矩阵计算方法如下：以图4为例，矩阵I有8个元素1-8，则其GLCM大小为8×8，GLCM中每个元素代表矩阵I中元素i、j（i，j=1，2，3，…，8）邻接的频数，从例图中可以看出：GLCM中第一行第一列的元素1表示矩阵I中两个元素1水平邻接个次数为1；GLCM中第一行第二列的元素2表示矩阵I中元素1和元素2水平邻接个次数为2。然后使GLCM加上其转置矩阵，即得到其对称阵。接着把GLCM对称阵中的每个元素均除以该矩阵所有元素之和，及得到GLCM的对称、归一化矩阵。

其中纹理方向如图3所示，图3中，实线箭头为坡向方向，虚线箭头为纹理方向，W-E为西-东方向，SW-NE为西南-东北方向，S-N为南-北方向，SE-NW为东南-西北方向。纹理方向实质上是矩阵I中元素的邻接方向。图4中椭圆内的元素对1、1和1、2为水平方向（0°）邻接，其纹理方向为0°(W-E)；椭圆内的元素对7、8为45°方向邻接，其纹理方向为45°(SW-NE)；椭圆内的元素对4、2为90°方向邻接，其纹理方向为90°(S-N)；椭圆内的元素对7、3为135°邻接方向，其纹理方向为135°(SE-NW)）。

由于纹理窗口为3×3像元，本方法中元素之间的邻接距离均设为1，即图4中显示的元素对之间在纹理方向上距离为1个像元。

根据GLCM对称、归一化矩阵按照相关性（Cor）、对比度（Con）、二阶矩(Asm)、熵(Ent)和同质性(Hom)的公式分别计算四个方向的高程相关性、高程对比度、高程二阶矩、高程熵、高程同质性，并对它们取平均值，分别得到平均的高程相关性、平均的高程对比度、平均的高程二阶矩、平均的高程熵、平均的高程同质性。然后把它们赋值给平均的高程相关性矩阵、平均的高程对比度矩阵、平均的高程二阶矩矩阵、平均的高程熵矩阵、平均的高程同质性矩阵中相应的元素。

接着按从左到右、从上到下的顺序移动3×3像元检测窗口（每次移动一个像元），依次计算平均的高程相关性、平均的高程对比度、平均的高程二阶矩、平均的高程熵、平均的高程同质性，并赋值给相应矩阵的相应元素；最后按照上述方法分别计算得到平均的坡度相关性矩阵、平均的坡度对比度矩阵、平均的坡度二阶矩矩阵、平均的坡度熵矩阵、平均的坡度同质性矩阵和平均的坡向相关性矩阵、平均的坡向对比度矩阵、平均的坡向二阶矩矩阵、平均的坡向熵矩阵、平均的坡向同质性矩阵。

$\mathrm{Cor}=\underset{i,j=1}{\overset{N}{\mathrm{\Σ}}}{P}_{\mathrm{ij}}\left[\frac{\left({i-\mathrm{\μ}}_i\right)(j-{\mathrm{\μ}}_j)}{\sqrt{\left({{\mathrm{\σ}}_i}^2\right)\left({{\mathrm{\σ}}_j}^2\right)}}\right]---\left(3\right)$

$\mathrm{Con}=\underset{i,j=1}{\overset{N}{\mathrm{\Σ}}}{P}_{\mathrm{ij}}{(i-j)}^2---\left(4\right)$

$\mathrm{Asm}=\underset{i,j=1}{\overset{N}{\mathrm{\Σ}}}{P_{\mathrm{ij}}}^2---\left(5\right)$

$\mathrm{Ent}=\underset{i,j=1}{\overset{N}{\mathrm{\Σ}}}{P}_{\mathrm{ij}}\left({-\ln P}_{\mathrm{ij}}\right)---\left(6\right)$

$\mathrm{Hom}=\underset{i,j=1}{\overset{N}{\mathrm{\Σ}}}\frac{P_{\mathrm{ij}}}{1+{(i-j)}^2}---\left(7\right)$

其中，P_ij为GLCM中的像元值，i、j为矩阵元素的行列号，N为矩阵的行列总数，μ_i，μ_j，σ_i ²，σ_j ²的计算公式如下：

${\mathrm{\μ}}_i=\underset{i,j=1}{\overset{N}{\mathrm{\Σ}}}i\left(P_{\mathrm{ij}}\right)---\left(8\right)$

${\mathrm{\μ}}_j=\underset{i,j=1}{\overset{N}{\mathrm{\Σ}}}j\left(P_{\mathrm{ij}}\right)---\left(9\right)$

${{\mathrm{\σ}}_i}^2=\underset{i,j=1}{\overset{N}{\mathrm{\Σ}}}{P}_{\mathrm{ij}}{(i-{\mathrm{\μ}}_i)}^2---\left(10\right)$

${{\mathrm{\σ}}_j}^2=\underset{i,j=1}{\overset{N}{\mathrm{\Σ}}}{P}_{\mathrm{ij}}{(j-{\mathrm{\μ}}_j)}^2---\left(11\right)$

2）计算坡向方向的高程、坡度和坡向相关性、对比度、二阶矩、熵、同质性。

首先，根据坡向值确定纹理方向，其中正北方向的坡向值为0，坡向值按顺时针方向逐渐增大。坡向值和纹理方向的关系见附图3。当坡向值在区间[022.5]、(157.5202.5]和(337.5360]内时，把纹理方向定为90°(S-N)；当坡向值在区间(22.567.5]和(202.5247.5]内时，把纹理方向定为45°(SW-NE)；当坡向值在区间(67.5112.5]和(247.5292.5]内时，把纹理方向定为0°(W-E)；当坡向值在区间(112.5 157.5]和(292.5 337.5]内时，把纹理方向定为135°(SE-NW)。然后，按照上述1)中的方法计算由坡向方向确定的纹理方向的高程、坡度和坡向相关性、对比度、二阶矩、熵、同质性，最后可以得到：基于坡向方向的高程相关性矩阵、基于坡向方向的高程对比度矩阵、基于坡向方向的高程二阶矩矩阵、基于坡向方向的高程熵矩阵、基于坡向方向的高程同质性矩阵，基于坡向方向的坡度相关性矩阵、基于坡向方向的坡度对比度矩阵、基于坡向方向的坡度二阶矩矩阵、基于坡向方向的坡度熵矩阵、基于坡向方向的坡度同质性矩阵和基于坡向方向的坡向相关性矩阵、基于坡向方向的坡向对比度矩阵、基于坡向方向的坡向二阶矩矩阵、基于坡向方向的坡向熵矩阵、基于坡向方向的坡向同质性矩阵。

3）将上面得到的39个矩阵整合到一个文本文件中，即得到研究区39个特征参数的特征参数文件，以下步骤中的计算和分析均基于该特征参数文件。

本发明较佳实施例中，步骤S3中确定最优特征参数组合具体为：

选取已知滑坡像元、非滑坡像元数据集中所有的滑坡像元和等数量的非滑坡像元作为训练集，并以该训练集为已知数据利用随机森林算法进行计算，得到特征参数文件中所有特征参数重要性的值，并对其从高到低进行排序，将前x（x指图5中袋外误差最小时对应的特征参数个数，x为自然数）个特征参数作为最优特征参数组合。

本发明实施例中，如图5所示，根据计算得到袋外误差最小时对应的特征参数个数为10，则前十个重要的特征参数为：高程值，高程平均值，高程标准差，坡度值，平均的坡度对比度，坡度平均值，坡度标准差，坡向值，坡向平均值，坡向标准差，即最优特征参数组合。其中袋外误差为：在构建分量分类器时，未被选中的样本组成袋外（out-of-bag，OOB）数据集，用袋外数据进行测试得到的误差被称之为袋外误差（out-of-bag error，OOB Err）。

本发明的一个实施例中，步骤S4中获得平衡系数具体过程为：

根据选择的像元确定已知滑坡像元、非滑坡像元数据集中非滑坡像元与滑坡像元的比值m，m＞0；

选取已知滑坡像元、非滑坡像元的数据集中20%的滑坡像元和一定数量的非滑坡像元作为训练集，训练集中的非滑坡像元数目为训练集的滑坡像元数目的k倍，k=1+0.1*n，n为整数；

n从0开始取值，n每次增加1，训练集随着n值的变化而不断变化，用该训练集训练随机森林模型，预测整个已知数据集，并计算平均用户精度和平均生产者精度，直至k小于或等于m；

以k值为横轴，平均用户精度、平均生产者精度为纵轴作曲线，获得曲线中平均用户精度和平均生产者精度之差绝对值最小时的k值，将此时的k值作为平衡系数。其中生产者精度为已知数据集的分类结果经与地面实况核对的判对率；用户精度是正确分类的某类像元数目与分为该类的采样个数的比值。

本发明实施例中，如图6所示，已知数据集中非滑坡像元数量与滑坡像元数量之比为4.3，计算得到的平均用户精度和平均生产者精度之差绝对值最小时的k值为2.7。采用k值为2.7的训练集，以及最优特征参数组合，训练随机森林模型，然后预测已知滑坡像元、非滑坡像元数据集。计算得到平均用户精度为80.5%，平均生产者精度为80.4%，总体精度为87.9%。

由上述精度评价结果可知，滑坡像元和非滑坡像元分类精度很高。在该前提下，进行滑坡自动识别，可以保证自动识别结果的可靠性和正确性，具体步骤如下：

选取已知滑坡像元、非滑坡像元数据集中所有的滑坡像元和一定数量的非滑坡像元作为训练集，训练集中的非滑坡像元数目与训练集的滑坡像元数目的比值为平衡系数k值。采用该训练集，以及最优特征参数组合，训练随机森林模型，然后预测整个研究区的滑坡像元、非滑坡像元数据集。根据分类结果，使用Canny边缘检测算子可以运算得到滑坡边界，从而实现滑坡自动识别。

利用边缘检测算子自动提取滑坡边界信息，设置较大的阈值，从而剔除零散、破碎且区域很小的滑坡像元区域，得到了图7所示的结果。自动识别的边界与实地调查结果基本一致。

由上述本发明实施例可知本发明基于LiDAR的复杂地质背景区滑坡识别的方法：（1）数据需求少：仅仅利用LIDAR-DEM数据及研究区内很少部分的滑坡像元和非滑坡像元；（2）充分挖掘了LiDAR-DEM数据在滑坡地形分析中的应用潜力：在常用地形参数的基础上，结合复杂地质背景区滑坡识别的瓶颈和高精度LiDAR-DEM的数据特征，提出坡度、坡向纹理信息及基于坡向的纹理特征，并通过其进行滑坡自动识别；（3）模型分类精度很高，实用性强，能够实现复杂地质背景区滑坡边界的自动识别。

本发明实施例基于机载激光雷达数据的复杂地质背景区滑坡识别系统用于实现上述方法，其中智能算法采用随机森林算法，相应地智能分类模型为随机森林模型。如图8所示，包括：

LiDAR数据处理模块10，用于对LiDAR数据进行处理生成消除植被影响后裸地表的DEM数据；

特征参数文件生成模块20，用于从DEM数据中提取传统的地貌特征参数以及从地貌特征参数中提取纹理特征参数，生成特征参数文件；

最优特征参数组合确定模块30，用于根据选定的研究区的部分已知滑坡像元和非滑坡像元，利用随机森林算法从特征参数文件中确定最优特征参数组合；较佳的最优特征参数组合确定过程在上文已有详细描述，在此不赘述。

平衡系数计算模块40，用于选择研究区的部分已知滑坡像元和非滑坡像元作为训练集，并不断改变训练集中的元素，获得满足预设精度条件的平衡系数，该平衡系数为满足预设条件时训练集中的非滑坡像元数目与训练集的滑坡像元数目的比值；较佳的平衡系数计算过程在上文已有详细描述，在此不赘述。

精度计算模块50，用于根据研究区的部分已知滑坡像元和非滑坡像元以及最优特征参数组合，训练随机森林模型，预测已知滑坡像元、非滑坡像元的数据集，并计算平均用户精度、平均生产者精度和总体精度；具体计算过程在上文已有详细描述，在此不赘述。

滑坡识别模块60，用于在所计算的平均用户精度、平均生产者精度和总体精度满足精度要求时，利用满足平衡系数的训练集以及最优特征参数组合，训练随机森林模型，并预测整个研究区的滑坡像元、非滑坡像元数据集，使用边缘检测算子计算滑坡边界，实现滑坡识别。

进一步地，本发明实施例中传统的地貌特征参数包括高程、坡度、坡向以及像元检测窗口内的高程、坡度、坡向的平均值和标准差；

纹理特征参数包括相关性、对比度、二阶矩、熵、同质性以及在四个纹理方向的平均值，基于坡向方向的相关性、对比度、二阶矩、熵、同质性；

应当理解的是，对本领域普通技术人员来说，可以根据上述说明加以改进或变换，而所有这些改进和变换都应属于本发明所附权利要求的保护范围。

Claims (9)

1.一种基于机载激光雷达数据的复杂地质背景区滑坡识别的方法，其特征在于，包括以下步骤：

S1、对LiDAR数据进行处理生成消除植被影响后裸地表的DEM数据；

S2、从所述DEM数据中提取传统的地貌特征参数以及根据所述地貌特征参数计算纹理特征参数，生成特征参数文件；

S3、根据选定的研究区的部分已知滑坡像元和非滑坡像元，利用智能分类算法对所述特征参数文件中的特征参数进行重要性计算，并按重要性从高到低对特征参数进行排序，将前x个特征参数作为最优特征参数组合，其中x为袋外误差最小时对应的特征参数的个数； 

S4、选择所述研究区的部分已知滑坡像元、非滑坡像元的数据集作为训练集，并不断改变所述训练集中的元素，结合所述最优特征参数组合，利用智能分类算法预测已知滑坡像元、非滑坡像元的数据集并计算平均用户精度和平均生产者精度，获得二者之差的绝对值最小时训练集中的非滑坡像元数目与训练集的滑坡像元数目的比值，该比值为满足预设精度条件的平衡系数；

S5、利用满足所述平衡系数的训练集以及所述最优特征参数组合，训练智能分类模型，预测已知滑坡像元、非滑坡像元的数据集，并计算平均用户精度、平均生产者精度和总体精度；

S6、若所计算的平均用户精度、平均生产者精度和总体精度满足精度要求，则利用满足所述平衡系数的训练集以及所述最优特征参数组合，训练智能分类模型，并预测整个研究区的滑坡像元、非滑坡像元数据集，使用边缘检测算子计算滑坡边界，实现滑坡识别；

其中，步骤S4中获得所述平衡系数具体过程为：

根据选择的像元确定已知滑坡像元、非滑坡像元数据集中非滑坡像元与滑坡像元的比值m，m>0；

选取已知滑坡像元、非滑坡像元的数据集中20%的滑坡像元和一定数量的非滑坡像元作为训练集，训练集中的非滑坡像元数目为训练集的滑坡像元数目的k倍，k=1+0.1*n，n为整数；

n从0开始取值，n每次增加1，训练集随着n值的变化而不断变化，用该训练集训练智能分类模型，预测整个已知数据集，并计算平均用户精度和平均生产者精度，直至k小于或等于m；

以k值为横轴，平均用户精度、平均生产者精度为纵轴作曲线，获得曲线中平均用户精度和平均生产者精度之差绝对值最小时的k值，将此时的k值作为平衡系数。

2.根据权利要求1所述的基于机载激光雷达数据的复杂地质背景区滑坡识别的方法，其特征在于，所述智能分类算法为随机森林算法，所述智能分类模型为随机森林模型。

3.根据权利要求2所述的基于机载激光雷达数据的复杂地质背景区滑坡识别的方法，其特征在于，所述传统的地貌特征参数包括高程、坡度、坡向以及像元检测窗口内的高程、坡度、坡向的平均值和标准差。

4.根据权利要求3所述的基于机载激光雷达数据的复杂地质背景区滑坡识别的方法，其特征在于，所述像元检测窗口为3×3像元检测窗口。

5.根据权利要求4所述的基于机载激光雷达数据的复杂地质背景区滑坡识别的方法，其特征在于，所述纹理特征参数包括相关性、对比度、二阶矩、熵和同质性在四个纹理方向的平均值，以及基于坡向方向的相关性、对比度、二阶矩、熵和同质性。

6.根据权利要求5所述的基于机载激光雷达数据的复杂地质背景区滑坡识别的方法，其特征在于，所述四个纹理方向包括西-东方向、西南-东北方向，南-北方向和东南-西北方向。

7.根据权利要求6所述的基于机载激光雷达数据的复杂地质背景区滑坡识别的方法，其特征在于，步骤S3具体为：

选取已知滑坡像元、非滑坡像元数据集中所有的滑坡像元和等数量的非滑坡像元作为训练集，并以该训练集为已知数据利用智能分类算法进行计算，得到所述特征参数文件中所有特征参数重要性的值，并对其从高到低进行排序，将前x个特征参数作为最优特征参数组合，其中x为袋外误差最小时对应的特征参数的个数。

8.一种基于机载激光雷达数据的复杂地质背景区滑坡识别系统，其特征在于,包括：

LiDAR数据处理模块，用于对LiDAR数据进行处理生成消除植被影响后裸地表的DEM数据；

特征参数文件生成模块，用于从所述DEM数据中提取传统的地貌特征参数以及根据所述地貌特征参数计算纹理特征参数，生成特征参数文件；

最优特征参数组合确定模块，用于根据选定的研究区的部分已知滑坡像元和非滑坡像元，利用智能分类算法对所述特征参数文件中的特征参数进行重要性计算，并按重要性从高到低对特征参数进行排序，将前x个特征参数作为最优特征参数组合，其中x为袋外误差最小时对应的特征参数的个数；平衡系数计算模块，用于选择所述研究区的部分已知滑坡像元、非滑坡像元的数据集作为训练集，并不断改变所述训练集中的元素，结合所述最优特征参数组合，利用智能分类算法预测已知滑坡像元、非滑坡像元的数据集并计算平均用户精度和平均生产者精度，获得二者之差的绝对值最小时训练集中的非滑坡像元数目与训练集的滑坡像元数目的比值，该比值为预设精度条件的平衡系数；

其中，获得所述平衡系数具体过程为：

根据选择的像元确定已知滑坡像元、非滑坡像元数据集中非滑坡像元与滑坡像元的比值m，m>0；

选取已知滑坡像元、非滑坡像元的数据集中20%的滑坡像元和一定数量的非滑坡像元作为训练集，训练集中的非滑坡像元数目为训练集的滑坡像元数目的k倍，k=1+0.1*n，n为整数；

n从0开始取值，n每次增加1，训练集随着n值的变化而不断变化，用该训练集训练智能分类模型，预测整个已知数据集，并计算平均用户精度和平均生产者精度，直至k小于或等于m；

以k值为横轴，平均用户精度、平均生产者精度为纵轴作曲线，获得曲线中平均用户精度和平均生产者精度之差绝对值最小时的k值，将此时的k值作为平衡系数；

精度计算模块，用于根据所述研究区的部分已知滑坡像元和非滑坡像元以及所述最优特征参数组合，训练智能分类模型，预测已知滑坡像元、非滑坡像元的数据集，并计算平均用户精度、平均生产者精度和总体精度；

滑坡识别模块，用于在所计算的平均用户精度、平均生产者精度和总体精度满足精度要求时，利用满足所述平衡系数的训练集以及所述最优特征参数组合，训练智能分类模型，并预测整个研究区的滑坡像元、非滑坡像元数据集，使用边缘检测算子计算滑坡边界，实现滑坡识别。

9.根据权利要求8所述的基于机载激光雷达数据的复杂地质背景区滑坡识别系统，其特征在于，所述传统的地貌特征参数包括高程、坡度、坡向以及像元检测窗口内的高程、坡度、坡向的平均值和标准差；

所述纹理特征参数包括相关性、对比度、二阶矩、熵和同质性在四个纹理方向的平均值，以及基于坡向方向的相关性、对比度、二阶矩、熵和同质性。