本发明涉及滑坡监测领域，室一种基于激光雷达技术的滑坡体变化检测方法。
背景技术：
伴随着经济社会发展的加速，经济活动的加剧，对自然资源的过度开发，人口数量的剧增，生态环境的恶化，使得我国自然灾害呈现日趋严重的发展趋势，抗灾救灾工作面临严峻挑战。
滑坡灾害是斜坡岩土体沿着惯通的剪切破坏面所发生的滑移现象，其机制是某一滑移面上剪应力超过了该面的抗剪强度所致。中国西南地区，特别是西南丘陵山区，最基本的地形地貌特征就是山体众多，山势陡峻，沟谷河流遍布于山体之中，与之相互切割，因而形成众多的具有足够滑动空间的斜坡体和切割面，广泛存在滑坡发生的基本条件，滑坡灾害相当频繁。另外，地震对滑坡的影响很大，首先是地震的强烈作用使斜坡土石的内部结构发生破坏和变化，原有的结构面张裂、松弛，加上地下水也有较大变化，特别是地下水位的突然升高或降低对斜坡稳定是很不利的。另外，一次强烈地震的发生往往伴随着许多余震，在地震力的反复振动冲击下，斜坡土石体就更容易发生变形，最后就会发展成滑坡。
此外，由于地质构造、降雨造成的滑坡也时有发生。滑坡灾害的特点是发生突然，破坏力巨大，特别是在滑坡救援过程中，摇摇欲坠的滚石、岩石边坡、土体边坡和危岩时刻对受灾群众和救援人员构成威胁。这些滑坡对建筑物、农田、道路及地下管线等产生危害或破坏，造成一系列的生态环境问题。
针对滑坡，国内外已经做了很多研究，传统检测滑坡的方法有：大地精密测量法、gps监测点布设，insar技术等，但具有一定优势的同时存在一定的缺点。例如，传统的大地测量方式和gps监测点布设为点式检测手段，难以获取空间的信息。而insar技术作为高精度监测地表形变的技术，在地表起伏较大时易产生相位失相关，从而影响计算精度。
作为空间信息技术之一的激光雷达技术，不受天气和周围自然环境的影响，而且采集数据精度和效率方面都明显优于其它方式，目前已经在滑坡等自然灾害动态检测，以及预报预测方面体现了非常重要的作用和优势。
技术实现要素：
本发明的目的是提供一种基于激光雷达的滑坡体变化检测方法，该方法在地震和滑坡灾害中，尤其对复杂地形下潜在的滑坡体进行变化检测，自动检测某个时间段滑坡体的位置和姿态发生的变化，一旦检测出潜在滑坡体发生的变化，并圈定变化范围和预测下一步滑坡体可能发生的变化，提醒现场工作人员注意，采取安全措施或者离开危险区域。
为达到上述目的，本发明的技术解决方案是：
一种基于激光雷达的滑坡体变化检测方法，对潜在的滑坡体进行的变化检测；其特征在于，包括具体步骤如下：
a)选取滑坡体研究区不同时间的两期或者两期以上激光雷达点云数据；
b)首先对步骤a)中选取研究区的激光雷达点云数据进行预处理，根据地形的复杂程度对激光雷达点云数据进行重采样，提高后期计算的高效性和精确性；
c)对步骤b)重采样后的点云数据进行坐标转换，激光雷达获得的点云数据坐标系是wgs84下的世界大地坐标系[xyz]，将重采样后的滑坡体点云数据转换成当地的地方坐标系[xyz]，以便统一；
d)将电影数据转换到地方坐标系以后，对其进行迭代线性最小二乘内插滤波，滤波后分类出地面点和非地面点，地面点即为剔除植被、建筑物后，真实地形表面点，滤波后的地面点可以更加真实地模拟滑坡体模型；
e)步骤d)滤波后得到的地面点，对地面点云数据构建滑坡体研究区的dem；
f)对滑坡体研究区的两期dem数据计算分别得到其每期dem的坡度模型、坡向模型，从整个潜在的滑坡区域按照坡度和坡向的大小分析危险滑坡区域；
g)对步骤e)处理后得到的滑坡体研究区dem，如果该滑坡体是自然地形形成的，利用dem得到研究区的等高线，如果该滑坡体是人工干预后形成的(比如露天矿的边坡、水坝的边坡)，可以根据地形的特点得到滑坡体的坡顶坡底的特征线、台阶的特征线和等高线；
h)对步骤e)、f)和g)中的滑坡体研究区的dem、坡度模型、坡向模型、等高线模型，坡顶坡底的特征线模型进行做差计算，得到某段时间内滑坡体的高度差变化模型、坡度差变化模型、坡向差变化模型、等高线差值变化模型、坡顶坡底的特征线差值变化模型、台阶的特征线差值变化模型；
i)对步骤h)中的高度差变化模型估算滑坡体的滑坡岩石土方量，为后期的滑坡清理运输提供正确的岩石土方量数据；
j)根据步骤h)中得到的滑坡体的各种变化模型，对滑坡体的变化进行变化检测。
所述的滑坡体变化检测方法，其特征在于，所述潜在滑坡体，为岩石边坡、土体边坡、危岩和滚石。
所述的基于激光雷达的滑坡体变化检测方法，其特征在于，所述a)步激光雷达点云数据可以是机载的雷达数据和地面的雷达数据的任意一种，数量上至少为两期扫描数据，本专利申请书采用了两期地面激光扫描的雷达数据。
所述的滑坡体变化检测方法，其特征在于，所述b)步，激光雷达扫描的数据量是巨大的，预处理重采样是根据实际地形的复杂程度决定滑坡体点云数据点密度，点密度过高则会影响处理的速度，点密度过低则不能很好的描述地表模型，因此需要选择合适的点密度，用曲率分析的方法进行点云数据的重采样，既要数据量不宜过大，还要最大程度地保留地形的特征点。
所述的滑坡体变化检测方法，其特征在于，所述c)步中滑坡体研究区的点云模型从wgs84坐标转换为地方坐标系，不管是星载、机载和地面激光雷达，其坐标系都为wgs84坐标系，转换为地方坐标系可以为其他系统融合到统一坐标系统中。
所述的滑坡体变化检测方法，其特征在于，所述d)和e)步中，滑坡体的地表会覆盖植被或者其他地物，通过迭代线性最小二乘内插滤波后得到滑坡体最贴近地表的模型，地物点的高程值比其周围地面点的高程值高，内插后激光脚点高程与拟合面的拟合残差不服从正态分布，地物点的拟合残差都为正值且偏差较大，地面点的拟合残差较小且可能为负值.该方法为迭代计算，首先应用较粗的格网内插生成初始拟合表面，实际上该表面是介于真实地面(dtm)和地物表面(dsm)之间的一个面，再进行内插处理，以此迭代下去得到最终的dem。
所述的滑坡体变化检测方法，其特征在于，所述e)步骤中生成的dem为高分辨率的dem，目前来说，通过激光雷达获取的dem的精度相对与其他方式而言是比较高的，f)和g)步骤中提取出的坡度模型、坡向模型、等高线模型、坡顶坡底的特征线模型、台阶的特征线模型都是高精度的滑坡体模型。
高精度滑坡体模型特征包括：
a)滑坡体坡度模型可以从滑坡体研究区整体观察到哪些是潜在滑坡区域，按照坡度的大小对滑坡体研究区域进行分类：坡度小于30度属于比较稳定的区域，坡度在30度和50度之间属于易发生滑坡的区域，坡度大于50度属于极易发生滑坡的区域；
b)根据滑坡体坡向模型可以观察到滑坡体在发生滑坡以后的滑坡走向，以及滑坡岩土容易堆积的区域；
c)滑坡体的等高线模型、坡顶坡底特征线模型、台阶特征线模型的生成主要是针对在不同机理下形成的滑坡体，有自然形成的，有人工填挖形成，采取其中一种或者几种滑坡体分析模型。
滑坡体变化检测方法特征包括：
a)滑坡体坡度差值变化模型是由两期滑坡体坡度模型做差得到，这段时间内滑坡体研究区的坡度变化，根据滑坡体坡度差值变化模型可以分析出该段时间内滑坡体的变化，哪块区域发生了滑坡；
b)滑坡体坡向差值变化模型是由两期滑坡体坡向模型做差得到，这段时间内滑坡体研究区的坡向变化，根据滑坡体坡向差值变化模型可以分析出该段时间内滑坡体的变化走向，以及下一步的滑坡堆积区；
c)滑坡体的高度差变化模型是由两期dem做差得到，从整体模型的h角度来分析，哪块区域高度增加，哪块区域高度降低；
d)滑坡体的等高线差值变化模型、坡顶坡底特征线的差值变化模型、台阶特征线的差值变化模型分别是滑坡体模型两期的等高线模型、坡顶坡底特征线模型、台阶特征线模型进行比对，从细微线性角度观察滑坡体的变化，在必要时可以对模型进行任意角度和位置的剖切，分段观察滑坡体的变化。
所述的滑坡体变化检测方法，其特征在于，i)步骤中的岩石土方量的计算牵扯到后期的滑坡处理，及时转运滑石，对整个工期的进度进行合理的科学安排。
所述的滑坡体变化检测方法，其特征在于，h)、g)、h)和i)步骤中的多种滑坡体参数模型综合得到了滑坡体变化检测的方法，提出了对于滑坡体的变化进行了从整体到局部的变化检测方法。
附图说明
图1为本发明的一种基于激光雷达的滑坡体变化检测方法流程图；
图2为本发明的滑坡体灾害研究区的两期卫星影像叠加全景图(蓝色传送带为滑动后的位置，灰色传送带为滑动前的位置，红色线框内为滑坡区域)；
图3为滑坡灾害救援现场；
图4a为预处理重采样后2013年第一期滑坡体点云数据；
图4b为预处理重采样后2014年第二期局部滑坡体点云数据；
图5为滑坡体局部滤波前后效果图；
图6为滑坡体滤波后得到的dem；
图7a为dem数据提取的等高线图；
图7b为提取的坡顶坡底的特征线图；
图7c为提取的台阶的特征线图；
图8a为滑坡体2014年激光点云的坡度模型图(坡度大小0-87.6度)；
图8b为滑坡体2014年激光点云的坡向模型图(坡向大小0-180度)；
图9a为带有坐标格网的两期数据叠放在一起的位置关系图(浅灰色模型为2013年激光点云模型，深灰色模型为2014年激光点云模型)；
图9b为高度差变化模型图(上图为整体情况，下图为南邦情况)；
图9c为坡度差变化模型图(坡度变化大小-65.9-87.3度)；
图9d为等高线变化模型图；
图9e为台阶特征线变化模型图。
具体实施方式
本发明是一种基于激光雷达的滑坡体变化检测方法，本研究方法的实施主要由不少于两期的同一区域不同时相的激光扫描数据和计算机系统构成。在整体上而言，该变化检测方法的功耗很低、计算精度高、计算效率高，有以下几大特点：
特点1：本发明的检测方法需要配套两个核心模块，即激光雷达点云数据处理模块和变化检测模块。激光雷达点云数据处理模块可以实现：
a)因激光雷达点云数据通常获得包括研究区和非研究区的大量冗余的数据，需要对数据进行预处理和重采样，剔除研究区以外的冗余点，而且一般激光雷达点云数据的点密度较高，根据地形的实际情况对数据按照曲率分析进行重采样，既最大程度地保留了地形特征点，考虑到了数据的精度要求，又压缩了数据量，提高后期变化检测的效率；
b)研究区点云数据地方坐标系的转换，激光雷达器自身带有在wgs84坐标系的gps定位模块，所有采集的数据均为wgs84坐标，为实现和当地坐标系的统一，采用布尔沙七参数转换模型对点云数据进行空间坐标转换；
c)点云滤波是为了实现地面点和非地面点的分离，通过迭代线性最小二乘内插滤波，该滤波方法是基于线性迭代法的改进，相对于其他滤波方法，该滤波方法不仅精度明显优于其他方式，而且滤波的效率非常高，点云的数据量是以数以百万计的，运算量是相当大，滤波后可以得到最真实地面dem。
变化检测模块由以下功能构成：
a)通过处理以后的dem得到滑坡体研究区的坡度模型、坡向模型、等高线模型、坡顶坡底特征线模型、台阶特征线模型，得到分析滑坡体一些基本模型；
b)由两期或者两期以上的a)中的坡度模型、坡向模型、等高线模型、坡顶坡底特征线模型、台阶特征线模型得到坡度差变化模型、坡向差变化模型、等高差变化模型、坡顶坡底特征线变化模型、台阶特征线变化模型；
c)滑坡体研究区的高度差变化模型，通过两期或者多期dem差值计算得到，根据高度差变化模型可以得到整体滑坡区域或者部分滑坡区域的岩土方量，为后期的滑坡处理提供了科学准确的参考数据。
特点2：本发明采用了“基于离散曲率分析的方法”对原始激光雷达点云数据的重采样和坡顶坡底特征线、台阶特征线的提起。激光雷达点云数据都是离散的、没有规则的点，在点曲率变化较小的地方，重采样时保留较少一些点，而在曲率变化大，地形变化大的地方，重采样时将会保留较多的点，可以保证科技高精度模拟地形表面模型，既压缩了数据量的大小，又在最大程度上保留了地形的特征，数据在压缩后将会提高滑坡体变化检测的效率。在这些离散点曲率发生突变的地方一定是地形特征发生突变的地方，用线连接这些曲率发生突变的特征点就生成了坡顶坡底特征线、台阶特征线。
特点3：本发明采用了“迭代线性最小二乘内插滤波的方法”进行滤波，该滤波是在线性迭代的基础上进行的改进，该方法不仅考虑的多次迭代拟合迭代面，首先应用较粗的格网内插生成初始拟合表面，实际上该表面是介于真实地面(dtm)和地物表面(dsm)之间的一个面，再进行内插处理，以此迭代下去得到最终的dem，而且处理效率非常高，不管再大的数据量，基本上滤波的处理时间都是以分/秒为时间单位。
特点4：本发明的核心技术是“基于滑坡体变化参数集”的滑坡体变化检测方法，改变以往的单种或者几种参数，或者单纯依靠个别一定量标定点的移动进行滑坡体的变化检测。“基于滑坡体变化参数集”的滑坡体变化检测方法从初步得到坡度模型参数、坡向模型参数、等高线模型参数、坡顶坡底特征线模型参数、台阶特征线模型参数进行滑坡体的检测分析，并进一步得到坡度差变化模型参数、坡向差变化模型参数、等高线差变化模型参数、坡顶坡底特征线变化模型参数、台阶特征线变化模型参数，高度差变化模型参数构成滑坡体变化检测参数集，对滑坡体的变化进行动态检测，分析滑坡体的变化。
特点5：滑坡体变化检测的方法最滑坡体发生滑坡以后，岩石在位置和姿态上会发生一定变化，在某一区域范围进行堆积，通过高度差变化模型估计岩石土方堆积的量，为下一步的及时处理提供科学准确的数据。
实施例：
以2013年-2014年国内某大型露天矿滑坡为例，利用本研究方法对该滑坡体(包括岩石边坡、土体边坡、危岩和滚石)进行变化检测。
a)激光雷达点云数据的预处理：针对地面激光扫描获得的两期滑坡体点云数据进行预处理，预处理主要是处理扫描数据中的一些飞点和杂点，并按照曲率分析的方法进行点云的重采样，地面激光扫描在获得滑坡体点云数据的同时，也会扫描到滑坡区以外的其他地物点，这些滑坡体以外的地物点会干扰滑坡体的分析，重采样的目的是保证滑坡体点云数据量的大小足够小，这样可以保证计算的速度，但同样必须要保证数据的精度要求，按照曲率分析的重采样方法，以最大保留地形特征的方式对点云数据进行重采样，预处理后的两期数据如图4a，4b；
b)点云数据的坐标转换：滑坡体研究区的点云模型从wgs84坐标转换为地方坐标系，地面激光雷达自身有定位系统，接收gps信号对自身进行定位，其扫描的数据坐标系都为gws84坐标系下的坐标，根据七参数坐标转换进行转换，转换为地方坐标系可以为其他系统融合到统一坐标系统中，并对转换后的精度进行评估，转换的转换精度误差在±0.001m内；
c)点云数据的滤波：通过迭代线性最小二乘内插滤波后得到滑坡体最贴近地表的模型(图5滑坡体局部滤波前后效果图)，地物点的高程值比其周围地面点的高程值高，内插后激光脚点高程与拟合面的拟合残差不服从正态分布，地物点的拟合残差都为正值且偏差较大，地面点的拟合残差较小且可能为负值.该方法为迭代计算，首先应用较粗的格网内插生成初始拟合表面，实际上该表面是介于真实地面(dtm)和地物表面(dsm)之间的一个面，再进行内插处理，以此迭代下去得到最终的dem(图6)；
d)滑坡体点云数据的等高线、坡顶坡底的特征线、台阶的特征线模型：根据滑坡体dem提取的等高线图7a，采用离散曲率分析的方法提取坡顶坡底的特征线、台阶的特征线，滑坡体的每个点都存在曲率，通常点与点之间的曲率变化相差都不会太大，当曲率发生大的变化地方一定是滑坡体边坡存在特征线的地方，按照曲率突变拟合的思路提取坡顶坡底的特征线、台阶的特征线图7b、7c；
e)滑坡体的坡度、坡向模型：对滑坡体dem矢量栅格化，栅格化后的像素大小为0.1m(可以根据实际情况自己设定像素的大小)，然后计算滑坡体的坡度、坡向数值得到滑坡体的坡度模型(图8a)、坡向模型(图8b)；
f)滑坡体的高度差变化模型、坡度变化模型、坡向变化模型、坡顶坡底的特征线变化模型、台阶特征线变化模型：滑坡体滑坡变化后的相对应的滑坡体高差位置自然会有相对应的变化(两期数据叠放在一起的位置关系图9a)，两期数据进行差值计算得到高度差变化模型(图9b)，坡度、坡向差变化模型(图9c)由两期坡度、坡向模型进行差值计算，坡顶坡底的特征线变化模型、台阶特征线变化模型由两期的得到的特征进行拓扑关系分析得到相对应位置的特征线变化特征(图9d)；
g)滑坡岩石土方量计算：根据滑坡体滑坡变化后，按照高差的变化计算局部或者整体的滑石土方量；
h)通过坡度、坡向、等高线、坡顶坡底特征线、台阶特征线，以及高差变化参数对滑坡体的变化进行检测，改变了以往用单种或者几种参数对滑坡体的变化进行检测。