CN110726998B - 一种激光雷达扫描测定矿区采煤塌陷盆地的方法 - Google Patents
一种激光雷达扫描测定矿区采煤塌陷盆地的方法 Download PDFInfo
- Publication number
- CN110726998B CN110726998B CN201911018042.7A CN201911018042A CN110726998B CN 110726998 B CN110726998 B CN 110726998B CN 201911018042 A CN201911018042 A CN 201911018042A CN 110726998 B CN110726998 B CN 110726998B
- Authority
- CN
- China
- Prior art keywords
- area
- cloud data
- group
- point
- point cloud
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S17/00—Systems using the reflection or reradiation of electromagnetic waves other than radio waves, e.g. lidar systems
- G01S17/88—Lidar systems specially adapted for specific applications
- G01S17/89—Lidar systems specially adapted for specific applications for mapping or imaging
Abstract
本发明公开了一种激光雷达扫描测定矿区采煤塌陷盆地的方法，该方法包括以下步骤：一、点云数据的获取；二、一次激光雷达点云数据的滤波；三、二次激光雷达点云数据的滤波；四、选取参照区及参照区域点云数据的导出；步骤五、参照区较低点集的获取；步骤六、参照区栅格DEM高程误差的计算；步骤七、全区域点云数据处理及矿区采煤塌陷盆地的测定。本发明方法步骤简单，设计合理，减少在人为干预的误差，处理过程简单，避免了非地面点噪声和DEM建模与叠加算法而产生的误差，确保矿区采煤塌陷盆地测定的准确性，实用性强。
Description
技术领域
本发明属于矿区采煤塌陷盆地测定技术领域，尤其是涉及一种激光雷达扫描测定矿区采煤塌陷盆地的方法。
背景技术
矿区采煤塌陷盆地的测定，传统的方法是在矿区采煤塌陷盆地布设一定数量的监测点，通过测量这些点的坐标变化来拟合和近似描述矿区采煤地表塌陷的主要特征，不能真实地获取地表塌陷盆地的三维精细特征，且测量工作的效率低下；而采用无人机摄光学影测量方法，虽然可以通过空三建模来得到地表塌陷盆地的大致形态，但在高程和平面上的精度均不能满足沉陷监测的基本要求。无人机激光扫描技术能够快速并准确地获取矿区塌陷盆地的三维点云数据，该技术具有速度快、精度高且非接触测量的特点，所获取的三维坐标数据完整且真实可靠。基于无人机激光扫描技术测量矿区塌陷盆地的基本原理是，利用地表塌陷前、后无人机激光扫描所获取的地面点云数据，经过去噪处理并生成DEM，再对两次生成的DEM进行叠加从而得到矿区塌陷盆地。但上述过程难以保证点云数据处理的精度，可能因包含大量非地面点或因DEM算法不合理，而导致所获得的地表塌陷盆地偏离其真实状态。
发明内容
本发明所要解决的技术问题在于针对上述现有技术中的不足，提供一种激光雷达扫描测定矿区采煤塌陷盆地的方法，其方法步骤简单，设计合理，减少在人为干预的误差，处理过程简单，避免了非地面点噪声和DEM建模与叠加算法而产生的误差，确保矿区采煤塌陷盆地测定的准确性，实用性强。
为解决上述技术问题，本发明采用的技术方案是：一种激光雷达扫描测定矿区采煤塌陷盆地的方法，其特征在于：
步骤一、点云数据的获取：
步骤101、采用无人机搭载激光雷达对待测定矿区进行一次激光雷达扫描，获取一次激光雷达点云数据；其中，一次激光雷达点云数据包括第1个待测定矿区扫描点坐标数据、第2个待测定矿区扫描点坐标数据，...，第i个待测定矿区扫描点坐标数据，...，第n′个待测定矿区扫描点坐标数据；其中，i和n′均为正整数，n′≥100，且1≤i≤n′，第i个待测定矿区扫描点坐标数据记作Pi(xi,yi,zi)，xi表示第i个待测定矿区扫描点在WGS-84坐标系下的X坐标，yi表示第i个待测定矿区扫描点在WGS-84坐标系下的Y坐标，zi表示第i个待测定矿区扫描点在WGS-84坐标系下的Z坐标；
步骤102、采用无人机搭载激光雷达对矿区采煤塌陷盆地进行二次激光雷达扫描，获取二次激光雷达点云数据；
步骤二、一次激光雷达点云数据的滤波：
采用计算机利用TerraSolid软件，对一次激光雷达点云数据进行滤波，得到一次扫描滤波后点云数据，并将一次扫描滤波后点云数据保存为ground1.las文件；
步骤三、二次激光雷达点云数据的滤波：
采用计算机利用TerraSolid软件，对二次激光雷达点云数据进行滤波，得到二次扫描滤波后点云数据，并将二次扫描滤波后点云数据保存为ground2.las文件；
步骤四、选取参照区及参照区域点云数据的导出：
步骤401、采用计算机将ground1.las文件和ground2.las文件导入到TerraSolid软件，得到点云覆盖区域；其中，点云覆盖区域为矩形区域；
步骤402、采用计算机利用TerraSolid软件中Place Block工具，并将Place Block工具中Method设置为Rotated，在点云覆盖区域外画矩形ABCD，并将AB边的中点标记为a，BC边的中点标记为b，CD边的中点标记为c，DA边的中点标记为d；
步骤403、采用计算机利用TerraSolid软件中Draw Vertical Section工具设置剖面的宽度为1m，分别沿AC、BD、ac和bd方向画剖面线，得到第一剖面图、第二剖面图、第三剖面图和第四剖面图；
步骤404、采用计算机利用TerraSolid软件中Place Block工具，将Place Block工具中Method设置为Orthogonal，将第一剖面图中稳定区域框选，得到第一稳定区域和第二稳定区域；
采用计算机利用TerraSolid软件中Place Block工具将第二剖面图中稳定区域框选，得到第三稳定区域和第四稳定区域；
采用计算机利用TerraSolid软件中Place Block工具将第三剖面图中稳定区域框选，得到第五稳定区域和第六稳定区域；
采用计算机利用TerraSolid软件中Place Block工具将第四剖面图中稳定区域框选，得到第七稳定区域和第八稳定区域；其中，第一稳定区域、第二稳定区域、第三稳定区域、第四稳定区域、第五稳定区域、第六稳定区域、第七稳定区域和第八稳定区域中一次扫描滤波后点云数据和二次扫描滤波后点云数据的高程相同；
步骤405、采用计算机利用TerraSolid软件中Place Point or Stream Curve工具将第一稳定区域、第二稳定区域、第三稳定区域、第四稳定区域、第五稳定区域、第六稳定区域、第七稳定区域和第八稳定区域靠近点云覆盖区域中心的内边界点用曲线连接；其中，曲线内区域为非稳定区域，曲线外与外界矩形ABCD之间区域为稳定区域；
步骤406、采用计算机利用TerraSolid软件中Place Tile工具从稳定区域中框选四个参照区；其中，每个参照区为32m×32m的正方形区域，且四个参照区分别为1参照区、2参照区、3参照区和4参照区，1参照区位于A顶点、d中点、a中点和曲线之间围成的区域，2参照区位于B顶点、a中点、b中点和曲线之间围成的区域，3参照区位于C顶点、b中点、c中点和曲线之间围成的区域，4参照区位于D顶点、c中点、d中点和曲线之间围成的区域；
步骤407、采用计算机利用TerraSolid软件将一次扫描滤波后点云数据和二次扫描滤波后点云数据关闭，将一次扫描滤波后点云数据ground1.las文件加载，利用TerraSolid软件中Delete point工具中的Outside fence工具，分别得到1参照区中点云数据、2参照区中点云数据、3参照区中点云数据和4参照区中点云数据，并将1参照区中点云数据、2参照区中点云数据、3参照区中点云数据和4参照区中点云数据分别记作一次第一组点云数据、一次第二组点云数据、一次第三组点云数据和一次第四组点云数据；
步骤408、采用计算机利用TerraSolid软件将一次扫描滤波后点云数据关闭，将二次扫描滤波后点云数据ground2.las文件加载，利用TerraSolid软件中Delete point工具中的Outside fence工具，分别得到1参照区中点云数据、2参照区中点云数据、3参照区中点云数据和4参照区中点云数据，并将1参照区中点云数据、2参照区中点云数据、3参照区中点云数据和4参照区中点云数据分别记作二次第一组点云数据、二次第二组点云数据、二次第三组点云数据和二次第四组点云数据；
步骤五、参照区较低点集的获取：
步骤501、采用计算机并将一次第一组点云数据按照X坐标从小到大，Y坐标从小到大顺序进行排序；
步骤502、采用计算机将一次第一组点云数据中X坐标最小值记作
一次第一组点云数据中X坐标最大值记作
一次第一组点云数据中Y坐标最小值记作
一次第一组点云数据中Y坐标最大值记作
采用计算机对一次第一组点云数据进行遍历，将一次第一组点云数据划分为多个遍历间隔为α的点集组；其中，任一个遍历间隔为α的点集组中点的X坐标在
区间，且第n+1个遍历间隔为α的点集组中点的Y坐标在
区间，n和m均为自然数，且n的取值范围为0≤n≤N1-1，m的取值范围为0≤m≤N2-1，且
α表示一次第一组点云数据的遍历间隔，ceil(·)为向上取整函数；
步骤503、采用计算机将多个遍历间隔为α的点集组分别记作第1个遍历间隔为α的点集组、第2个遍历间隔为α的点集组，...，第j个遍历间隔为α的点集组，...，第M个遍历间隔为α的点集组，j和M为正整数，且1≤j≤M，M表示遍历间隔为α的点集组的总数；
步骤504、采用计算机从第j个遍历间隔为α的点集组中提取较低点集的过程如下：
步骤5041、采用计算机将第j个遍历间隔为α的点集组中点集按照Z坐标从小到大进行排序；
步骤5042、采用计算机从第j个遍历间隔为α的点集组中筛选出Z坐标较小点集；其中，Z坐标较小点集中点的数量占第j个遍历间隔为α的点集组中点的总数的比例为b％；
步骤5043、采用计算机将从第j个遍历间隔为α的点集组中选出的Z坐标较小点集记作第j个遍历间隔为α选取比例为b％的较低点集；
步骤5044、多次重复步骤5041至步骤5043，直至完成第M个遍历间隔为α的点集组的较低点集的筛选，得到一次第一组遍历间隔为α选取比例为b％的较低点集；
步骤505、按照步骤501至步骤504所述的方法得到一次第二组遍历间隔为α选取比例为b％的较低点集、一次第三组遍历间隔为α选取比例为b％的较低点集和一次第四组遍历间隔为α选取比例为b％的较低点集；
步骤506、按照步骤501至步骤505所述的方法，得到二次第一组遍历间隔为α选取比例为b％的较低点集、二次第二组遍历间隔为α选取比例为b％的较低点集、二次第三组遍历间隔为α选取比例为b％的较低点集和二次第四组遍历间隔为α选取比例为b％的较低点集；
步骤六、参照区栅格DEM高程误差的计算：
步骤601、采用计算机利用ArcGIS软件中的添加XY数据工具，将一次第一组遍历间隔为α选取比例为b％的较低点集坐标导入；
步骤602、采用计算机利用ArcGIS软件中的自然邻域法插值工具，得到一次第一组遍历间隔为α选取比例为b％的点集数据栅格图；其中，输出像元大小为1m，一次第一组遍历间隔为α选取比例为b％的点集数据栅格图的尺寸为32m×32m；
步骤603、重复步骤601和步骤602，对二次第一组遍历间隔为α的点集组b％较低点集处理，得到二次第一组遍历间隔为α选取比例为b％的点集数据栅格图；
步骤604、采用计算机利用ArcGIS软件中的栅格计算工具，将二次第一组遍历间隔为α选取比例为b％的点集数据栅格图与一次第一组遍历间隔为α选取比例为b％的点集数据栅格图相减，得到遍历间隔为α选取比例为b％时的第一组差值栅格图；
步骤605、采用计算机利用ArcGIS软件中的栅格转点工具，将遍历间隔为α选取比例为b％时的第一组差值栅格图转换成遍历间隔为α选取比例为b％时的第一组差值点要素，并将遍历间隔为α选取比例为b％时的第一组差值点要素属性表导出为txt文件，得到遍历间隔为α选取比例为b％时的第一组差值点要素txt文件；其中，遍历间隔为α选取比例为b％时的第一组差值点要素txt文件包括差值栅格像元序号和差值栅格像元对应的高程差值；
步骤606、按照步骤601至步骤605所述的方法，对一次第二组遍历间隔为α选取比例为b％的点集数据栅格图和二次第二组遍历间隔为α选取比例为b％的点集数据栅格图处理，得到遍历间隔为α选取比例为b％时的第二组差值点要素txt文件；
对一次第三组遍历间隔为α选取比例为b％的点集数据栅格图和二次第三组遍历间隔为α选取比例为b％的点集数据栅格图处理，得到遍历间隔为α选取比例为b％时的第三组差值点要素txt文件；
对一次第四组遍历间隔为α选取比例为b％的点集数据栅格图和二次第四组遍历间隔为α选取比例为b％的点集数据栅格图处理，得到遍历间隔为α选取比例为b％时的第四组差值点要素txt文件；
步骤607、采用计算机根据公式
得到遍历间隔为α且选取比例为b％时的高程误差值；其中，
表示遍历间隔为α选取比例为b％时的第一组差值点要素txt文件中第f个差值栅格像元对应的高程差值，
表示遍历间隔为α选取比例为b％时的第二组差值点要素txt文件中第f个差值栅格像元对应的高程差值，
表示遍历间隔为α选取比例为b％时的第三组差值点要素txt文件中第f个差值栅格像元对应的高程差值，
表示遍历间隔为α选取比例为b％时的第四组差值点要素txt文件中第f个差值栅格像元对应的高程差值；其中，f为正整数，且1≤f≤1024；
步骤608、采用计算机按照步骤601至步骤607所述的方法，当α取1m、2m、4m、8m、16m或者32m，当b％取5％、10％或者15％时，分别得到遍历间隔为1m且选取比例为5％时的高程误差值e1,5、遍历间隔为1m且选取比例为10％时的高程误差值e1,10、遍历间隔为1m且选取比例为15％时的高程误差值e1,15；遍历间隔为2m且选取比例为5％时的高程误差值e2,5、遍历间隔为2m且选取比例为10％时的高程误差值e2,10、遍历间隔为2m且选取比例为15％时的高程误差值e2,15；遍历间隔为4m且选取比例为5％时的高程误差值e4,5、遍历间隔为4m且选取比例为10％时的高程误差值e4,10、遍历间隔为4m且选取比例为15％时的高程误差值e4,15；遍历间隔为8m且选取比例为5％时的高程误差值e8,5、遍历间隔为8m且选取比例为10％时的高程误差值e8,10、遍历间隔为8m且选取比例为15％时的高程误差值e8,15；遍历间隔为16m且选取比例为5％时的高程误差值e16,5、遍历间隔为16m且选取比例为10％时的高程误差值e16,10、遍历间隔为16m且选取比例为15％时的高程误差值e16,15；遍历间隔为32m且选取比例为5％时的高程误差值e32,5、遍历间隔为32m且选取比例为10％时的高程误差值e32,10、遍历间隔为32m且选取比例为15％时的高程误差值e32,15；
步骤609、采用计算机将步骤608中多个高程误差值按照从小到大顺序进行排序，得到最小高程误差值，并将最小高程误差值所对应的遍历间隔记作遍历间隔αmin，最小高程误差值所对应的选取比例记作选取比例bmin％；
步骤七、全区域点云数据处理及矿区采煤塌陷盆地的测定：
步骤701、采用计算机将ground1.las文件中点云数据按照X坐标从小到大，Y坐标从小到大顺序进行排序；将ground2.las文件中点云数据按照X坐标从小到大，Y坐标从小到大顺序进行排序；
步骤702、当选择遍历间隔为αmin且选取比例为bmin％，按照步骤502至步骤5044，对一次扫描滤波后点云数据处理，得到一次全区域点云数据的较低点集；
当选择遍历间隔为αmin且选取比例为bmin％，按照步骤502至步骤5044，对二次扫描滤波后点云数据进行处理，得到二次全区域点云数据的较低点集；
步骤703、按照步骤601至步骤604所述的方法，得到测定矿区全区域差值栅格图，则测定矿区全区域差值栅格图为矿区采煤塌陷盆地，实现矿区采煤塌陷盆地的测定。
上述的一种激光雷达扫描测定矿区采煤塌陷盆地的方法，其特征在于：步骤102中采用无人机搭载激光雷达对矿区采煤塌陷盆地进行激光扫描时，点云数据的间距小于0.15m，点云数据密度高于60个/平方米；且采用无人机搭载激光雷达对矿区采煤塌陷盆地进行激光扫描时，扫描区域为矩形；
采用无人机搭载激光雷达对待测定矿区进行一次激光雷达扫描和采用无人机搭载激光雷达对矿区采煤塌陷盆地进行二次激光雷达扫描之间的时间间隔为20days～60days。
上述的一种激光雷达扫描测定矿区采煤塌陷盆地的方法，其特征在于：步骤二中采用计算机利用TerraSolid软件，对一次激光雷达点云数据进行滤波，得到一次扫描滤波后点云数据，并数据保存为ground1.las文件，具体过程如下：
步骤201、将一次激光雷达点云数据导入TerraSolid软件中，利用TerraSolid软件中的Classify by class工具，将一次激光雷达点云数据全部放入Default类中；
步骤202、采用计算机利用TerraSolid软件中的Classify low points工具，将噪声点从Default类中分离出至low points类中；其中，Classify low points工具中MAXcount阈值范围为1～6，Classify low points工具中Search设置为Single poits，Classify low points工具中高度差阈值范围为0.3m～1m，距离搜索范围阈值为2m～8m；
步骤203、采用计算机利用TerraSolid软件中的Classify isolated points工具，将孤立点从Default类中分离出至low points类中；其中，Classify isolated points工具中距离搜索范围阈值为2m～8m；
步骤204、采用计算机利用TerraSolid软件中的Classify ground工具，将地面点从Default类中分离出至Ground类中；其中，Classify ground工具中Max buliding size阈值为10m～50m，Iteration angle即迭代角度阈值为4°～12°，Iteration distance即迭代距离阈值为0.5m～1.5m，Terrain angle即地形角为88°；
步骤205、采用计算机利用TerraSolid软件中的Classify height from ground工具，将地面以下点从Default类中分离出至Ground类中；其中，Classify height fromground工具中Max height阈值为零，Min height阈值为-999m；
步骤206、采用计算机将Ground类中点记作一次扫描滤波后点云数据，并数据保存为ground1.las文件。
上述的一种激光雷达扫描测定矿区采煤塌陷盆地的方法，其特征在于：步骤七中得到矿区采煤塌陷盆地，之后，具体执行如下：
采用计算机利用ArcGIS软件中的识别工具，查看矿区采煤塌陷盆地中各个栅格像元的对应的高程差值；
采用计算机利用ArcGIS软件中的栅格转点工具，将测定矿区全区域差值栅格图转换为测定矿区全区域点要素，并将测定矿区全区域点要素的属性数据表导出为txt文件，得到测定矿区全区域点要素txt文件；其中，测定矿区全区域点要素txt文件包括采煤塌陷盆地栅格像元序号和差值栅格像元对应的高程差值。
本发明与现有技术相比具有以下优点：
1、本发明方法步骤简单、设计合理且实现方便，精度高。
2、本发明首先利用分类算法去噪保留地面点而无需人为干预进行精细分类，一方面减小了在数据处理方面的工作量，一方面也避免了人为干预而产生的误差。
3、本发明首先利用TerraSolid软件对激光雷达点云数据进行滤波，之后，进行一次全区域点云数据的较低点集和二次全区域点云数据的较低点集的提取，实现对激光雷达点云数据进行二次滤波，减少在人为干预的误差，处理过程简单，避免了非地面点噪声而产生的误差。
4、本发明采用计算机通过点云覆盖区域中AC、BD、ac和bd方向画剖面线，得到第一剖面图、第二剖面图、第三剖面图和第四剖面图，从第一剖面图中框选得到第一稳定区域和第二稳定区域，第二剖面图中框选得到第三稳定区域和第四稳定区域，第三剖面图中框选得到第五稳定区域和第六稳定区域，第四剖面图中框选得到第七稳定区域和第八稳定区域，便于对点云覆盖区域进行非稳定区域和稳定区域的划分，从而便于四个参照区的选择，既保证了选取的参照区内地表较为稳定且选取的区域极具代表性，又不需要经过复杂的步骤找出沉陷边界，方法步骤简单易行，为后续全局点云数据的处理提供基础。
5、本发明采用计算机首先按照遍历间隔为α进行点集组的划分，得到多个点集组；然后对每个点集组按照选取比例为b％进行较低点集筛选，得到一次遍历间隔为α选取比例为b％的较低点集和二次遍历间隔为α选取比例为b％的较低点集，便于获取不同遍历间隔不同选取比例下的高程误差，便于得到应用于全局的最优较低点集选择所需的遍历间隔和选取比例，提高了全局较低点集获取的准确性，避免操作人员进行手动筛选，从而确保获取矿区采煤塌陷盆地的准确性。
6、本发明计算机将一次第一组遍历间隔为α选取比例为b％的较低点集和二次第一组遍历间隔为α的点集组b％较低点集、一次第二组遍历间隔为α选取比例为b％的较低点集和二次第二组遍历间隔为α的点集组b％较低点集、一次第三组遍历间隔为α选取比例为b％的较低点集和二次第三组遍历间隔为α的点集组b％较低点集以及一次第四组遍历间隔为α选取比例为b％的较低点集和二次第四组遍历间隔为α的点集组b％较低点集进行处理，得到遍历间隔为α选取比例为b％时的第一组差值栅格图、第二组差值栅格图、第三组差值栅格图和第四组差值栅格图，同时分别得到遍历间隔为α选取比例为b％时第一组差值点要素txt文件、第二组差值点要素txt文件、第三组差值点要素txt文件以及第四组差值点要素txt文件，便于获取各个差值栅格像元对应的高程差值，处理过程简单。
7、本发明根据获取的遍历间隔为α选取比例为b％时的各个差值栅格像元对应的高程差值，得到遍历间隔为α且选取比例为b％时的高程误差值，通过将多个高程误差值进行排序，得到最小高程误差值所对应的遍历间隔αmin，最小高程误差值所对应的选取比例，并将最小高程误差值所对应的遍历间隔αmin，最小高程误差值所对应的选取比例应用于一次全区域点云数据和二次全区域点云数据，从而获取一次全区域点云数据的较低点集和二次全区域点云数据的较低点集，提高了矿区采煤塌陷盆地测定的准确度，能够有效地避免了非地面点噪声和DEM建模与叠加算法而产生的误差。
综上所述，本发明方法步骤简单，设计合理，减少在人为干预的误差，处理过程简单，避免了非地面点噪声和DEM建模与叠加算法而产生的误差，确保矿区采煤塌陷盆地测定的准确性，实用性强。
下面通过附图和实施例，对本发明的技术方案做进一步的详细描述。
附图说明
图1为本发明的方法流程框图。
图2为本发明1参照区、2参照区、3参照区和4参照区的结构示意图。
图3为本发明第一剖面图的结构示意图。
图4为本发明矿区采煤塌陷盆地测定的结构示意图。
具体实施方式
如图1至图4所示的一种激光雷达扫描测定矿区采煤塌陷盆地的方法，包括以下步骤：
步骤一、点云数据的获取：
步骤101、采用无人机搭载激光雷达对待测定矿区进行一次激光雷达扫描，获取一次激光雷达点云数据；其中，一次激光雷达点云数据包括第1个待测定矿区扫描点坐标数据、第2个待测定矿区扫描点坐标数据，...，第i个待测定矿区扫描点坐标数据，...，第n′个待测定矿区扫描点坐标数据；其中，i和n′均为正整数，n′≥100，且1≤i≤n′，第i个待测定矿区扫描点坐标数据记作Pi(xi,yi,zi)，xi表示第i个待测定矿区扫描点在WGS-84坐标系下的X坐标，yi表示第i个待测定矿区扫描点在WGS-84坐标系下的Y坐标，zi表示第i个待测定矿区扫描点在WGS-84坐标系下的Z坐标；
步骤102、采用无人机搭载激光雷达对矿区采煤塌陷盆地进行二次激光雷达扫描，获取二次激光雷达点云数据；
步骤二、一次激光雷达点云数据的滤波：
采用计算机利用TerraSolid软件，对一次激光雷达点云数据进行滤波，得到一次扫描滤波后点云数据，并将一次扫描滤波后点云数据保存为ground1.las文件；
步骤三、二次激光雷达点云数据的滤波：
采用计算机利用TerraSolid软件，对二次激光雷达点云数据进行滤波，得到二次扫描滤波后点云数据，并将二次扫描滤波后点云数据保存为ground2.las文件；
步骤四、选取参照区及参照区域点云数据的导出：
步骤401、采用计算机将ground1.las文件和ground2.las文件导入到TerraSolid软件，得到点云覆盖区域；其中，点云覆盖区域为矩形区域；
步骤402、采用计算机利用TerraSolid软件中Place Block工具，并将Place Block工具中Method设置为Rotated，在点云覆盖区域外画矩形ABCD，并将AB边的中点标记为a，BC边的中点标记为b，CD边的中点标记为c，DA边的中点标记为d；
步骤403、采用计算机利用TerraSolid软件中Draw Vertical Section工具设置剖面的宽度为1m，分别沿AC、BD、ac和bd方向画剖面线，得到第一剖面图、第二剖面图、第三剖面图和第四剖面图；
步骤404、采用计算机利用TerraSolid软件中Place Block工具，将Place Block工具中Method设置为Orthogonal，将第一剖面图中稳定区域框选，得到第一稳定区域和第二稳定区域；
采用计算机利用TerraSolid软件中Place Block工具将第二剖面图中稳定区域框选，得到第三稳定区域和第四稳定区域；
采用计算机利用TerraSolid软件中Place Block工具将第三剖面图中稳定区域框选，得到第五稳定区域和第六稳定区域；
采用计算机利用TerraSolid软件中Place Block工具将第四剖面图中稳定区域框选，得到第七稳定区域和第八稳定区域；其中，第一稳定区域、第二稳定区域、第三稳定区域、第四稳定区域、第五稳定区域、第六稳定区域、第七稳定区域和第八稳定区域中一次扫描滤波后点云数据和二次扫描滤波后点云数据的高程相同；
步骤405、采用计算机利用TerraSolid软件中Place Point or Stream Curve工具将第一稳定区域、第二稳定区域、第三稳定区域、第四稳定区域、第五稳定区域、第六稳定区域、第七稳定区域和第八稳定区域靠近点云覆盖区域中心的内边界点用曲线连接；其中，曲线内区域为非稳定区域，曲线外与外界矩形ABCD之间区域为稳定区域；
步骤406、采用计算机利用TerraSolid软件中Place Tile工具从稳定区域中框选四个参照区；其中，每个参照区为32m×32m的正方形区域，且四个参照区分别为1参照区、2参照区、3参照区和4参照区，1参照区位于A顶点、d中点、a中点和曲线之间围成的区域，2参照区位于B顶点、a中点、b中点和曲线之间围成的区域，3参照区位于C顶点、b中点、c中点和曲线之间围成的区域，4参照区位于D顶点、c中点、d中点和曲线之间围成的区域；
步骤407、采用计算机利用TerraSolid软件将一次扫描滤波后点云数据和二次扫描滤波后点云数据关闭，将一次扫描滤波后点云数据ground1.las文件加载，利用TerraSolid软件中Delete point工具中的Outside fence工具，分别得到1参照区中点云数据、2参照区中点云数据、3参照区中点云数据和4参照区中点云数据，并将1参照区中点云数据、2参照区中点云数据、3参照区中点云数据和4参照区中点云数据分别记作一次第一组点云数据、一次第二组点云数据、一次第三组点云数据和一次第四组点云数据；
步骤408、采用计算机利用TerraSolid软件将一次扫描滤波后点云数据关闭，将二次扫描滤波后点云数据ground2.las文件加载，利用TerraSolid软件中Delete point工具中的Outside fence工具，分别得到1参照区中点云数据、2参照区中点云数据、3参照区中点云数据和4参照区中点云数据，并将1参照区中点云数据、2参照区中点云数据、3参照区中点云数据和4参照区中点云数据分别记作二次第一组点云数据、二次第二组点云数据、二次第三组点云数据和二次第四组点云数据；
步骤五、参照区较低点集的获取：
步骤501、采用计算机并将一次第一组点云数据按照X坐标从小到大，Y坐标从小到大顺序进行排序；
步骤502、采用计算机将一次第一组点云数据中X坐标最小值记作
一次第一组点云数据中X坐标最大值记作
一次第一组点云数据中Y坐标最小值记作
一次第一组点云数据中Y坐标最大值记作
采用计算机对一次第一组点云数据进行遍历，将一次第一组点云数据划分为多个遍历间隔为α的点集组；其中，任一个遍历间隔为α的点集组中点的X坐标在
区间，且第n+1个遍历间隔为α的点集组中点的Y坐标在
区间，n和m均为自然数，且n的取值范围为0≤n≤N1-1，m的取值范围为0≤m≤N2-1，且
α表示一次第一组点云数据的遍历间隔，ceil(·)为向上取整函数；
步骤503、采用计算机将多个遍历间隔为α的点集组分别记作第1个遍历间隔为α的点集组、第2个遍历间隔为α的点集组，...，第j个遍历间隔为α的点集组，...，第M个遍历间隔为α的点集组，j和M为正整数，且1≤j≤M，M表示遍历间隔为α的点集组的总数；
步骤504、采用计算机从第j个遍历间隔为α的点集组中提取较低点集的过程如下：
步骤5041、采用计算机将第j个遍历间隔为α的点集组中点集按照Z坐标从小到大进行排序；
步骤5042、采用计算机从第j个遍历间隔为α的点集组中筛选出Z坐标较小点集；其中，Z坐标较小点集中点的数量占第j个遍历间隔为α的点集组中点的总数的比例为b％；
步骤5043、采用计算机将从第j个遍历间隔为α的点集组中选出的Z坐标较小点集记作第j个遍历间隔为α选取比例为b％的较低点集；
步骤5044、多次重复步骤5041至步骤5043，直至完成第M个遍历间隔为α的点集组的较低点集的筛选，得到一次第一组遍历间隔为α选取比例为b％的较低点集；
步骤505、按照步骤501至步骤504所述的方法得到一次第二组遍历间隔为α选取比例为b％的较低点集、一次第三组遍历间隔为α选取比例为b％的较低点集和一次第四组遍历间隔为α选取比例为b％的较低点集；
步骤506、按照步骤501至步骤505所述的方法，得到二次第一组遍历间隔为α选取比例为b％的较低点集、二次第二组遍历间隔为α选取比例为b％的较低点集、二次第三组遍历间隔为α选取比例为b％的较低点集和二次第四组遍历间隔为α选取比例为b％的较低点集；
步骤六、参照区栅格DEM高程误差的计算：
步骤601、采用计算机利用ArcGIS软件中的添加XY数据工具，将一次第一组遍历间隔为α选取比例为b％的较低点集坐标导入；
步骤602、采用计算机利用ArcGIS软件中的自然邻域法插值工具，得到一次第一组遍历间隔为α选取比例为b％的点集数据栅格图；其中，输出像元大小为1m，一次第一组遍历间隔为α选取比例为b％的点集数据栅格图的尺寸为32m×32m；
步骤603、重复步骤601和步骤602，对二次第一组遍历间隔为α的点集组b％较低点集处理，得到二次第一组遍历间隔为α选取比例为b％的点集数据栅格图；
步骤604、采用计算机利用ArcGIS软件中的栅格计算工具，将二次第一组遍历间隔为α选取比例为b％的点集数据栅格图与一次第一组遍历间隔为α选取比例为b％的点集数据栅格图相减，得到遍历间隔为α选取比例为b％时的第一组差值栅格图；
步骤605、采用计算机利用ArcGIS软件中的栅格转点工具，将遍历间隔为α选取比例为b％时的第一组差值栅格图转换成遍历间隔为α选取比例为b％时的第一组差值点要素，并将遍历间隔为α选取比例为b％时的第一组差值点要素属性表导出为txt文件，得到遍历间隔为α选取比例为b％时的第一组差值点要素txt文件；其中，遍历间隔为α选取比例为b％时的第一组差值点要素txt文件包括差值栅格像元序号和差值栅格像元对应的高程差值；
步骤606、按照步骤601至步骤605所述的方法，对一次第二组遍历间隔为α选取比例为b％的点集数据栅格图和二次第二组遍历间隔为α选取比例为b％的点集数据栅格图处理，得到遍历间隔为α选取比例为b％时的第二组差值点要素txt文件；
对一次第三组遍历间隔为α选取比例为b％的点集数据栅格图和二次第三组遍历间隔为α选取比例为b％的点集数据栅格图处理，得到遍历间隔为α选取比例为b％时的第三组差值点要素txt文件；
对一次第四组遍历间隔为α选取比例为b％的点集数据栅格图和二次第四组遍历间隔为α选取比例为b％的点集数据栅格图处理，得到遍历间隔为α选取比例为b％时的第四组差值点要素txt文件；
步骤607、采用计算机根据公式
得到遍历间隔为α且选取比例为b％时的高程误差值；其中，
表示遍历间隔为α选取比例为b％时的第一组差值点要素txt文件中第f个差值栅格像元对应的高程差值，
表示遍历间隔为α选取比例为b％时的第二组差值点要素txt文件中第f个差值栅格像元对应的高程差值，
表示遍历间隔为α选取比例为b％时的第三组差值点要素txt文件中第f个差值栅格像元对应的高程差值，
表示遍历间隔为α选取比例为b％时的第四组差值点要素txt文件中第f个差值栅格像元对应的高程差值；
步骤608、采用计算机按照步骤601至步骤607所述的方法，当α取1m、2m、4m、8m、16m或者32m，当b％取5％、10％或者15％时，分别得到遍历间隔为1m且选取比例为5％时的高程误差值e1,5、遍历间隔为1m且选取比例为10％时的高程误差值e1,10、遍历间隔为1m且选取比例为15％时的高程误差值e1,15；遍历间隔为2m且选取比例为5％时的高程误差值e2,5、遍历间隔为2m且选取比例为10％时的高程误差值e2,10、遍历间隔为2m且选取比例为15％时的高程误差值e2,15；遍历间隔为4m且选取比例为5％时的高程误差值e4,5、遍历间隔为4m且选取比例为10％时的高程误差值e4,10、遍历间隔为4m且选取比例为15％时的高程误差值e4,15；遍历间隔为8m且选取比例为5％时的高程误差值e8,5、遍历间隔为8m且选取比例为10％时的高程误差值e8,10、遍历间隔为8m且选取比例为15％时的高程误差值e8,15；遍历间隔为16m且选取比例为5％时的高程误差值e16,5、遍历间隔为16m且选取比例为10％时的高程误差值e16,10、遍历间隔为16m且选取比例为15％时的高程误差值e16,15；遍历间隔为32m且选取比例为5％时的高程误差值e32,5、遍历间隔为32m且选取比例为10％时的高程误差值e32,10、遍历间隔为32m且选取比例为15％时的高程误差值e32,15；
步骤609、采用计算机将步骤608中多个高程误差值按照从小到大顺序进行排序，得到最小高程误差值，并将最小高程误差值所对应的遍历间隔记作遍历间隔αmin，最小高程误差值所对应的选取比例记作选取比例bmin％；
步骤七、全区域点云数据处理及矿区采煤塌陷盆地的测定：
步骤701、采用计算机将ground1.las文件中点云数据按照X坐标从小到大，Y坐标从小到大顺序进行排序；将ground2.las文件中点云数据按照X坐标从小到大，Y坐标从小到大顺序进行排序；
步骤702、当选择遍历间隔为αmin且选取比例为bmin％，按照步骤502至步骤5044，对一次扫描滤波后点云数据处理，得到一次全区域点云数据的较低点集；
当选择遍历间隔为αmin且选取比例为bmin％，按照步骤502至步骤5044，对二次扫描滤波后点云数据进行处理，得到二次全区域点云数据的较低点集；
步骤703、按照步骤601至步骤604所述的方法，得到测定矿区全区域差值栅格图，则测定矿区全区域差值栅格图为矿区采煤塌陷盆地，实现矿区采煤塌陷盆地的测定。
本实施例中，步骤102中采用无人机搭载激光雷达对矿区采煤塌陷盆地进行激光扫描时，点云数据的间距小于0.15m，点云数据密度高于60个/平方米；且采用无人机搭载激光雷达对矿区采煤塌陷盆地进行激光扫描时，扫描区域为矩形；
采用无人机搭载激光雷达对待测定矿区进行一次激光雷达扫描和采用无人机搭载激光雷达对矿区采煤塌陷盆地进行二次激光雷达扫描之间的时间间隔为20days～60days。
本实施例中，一次激光雷达扫描和二次激光雷达扫描之间的时间间隔为60days。
本实施例中，需要说明的是，
和相等。
本实施例中，步骤二中采用计算机利用TerraSolid软件，对一次激光雷达点云数据进行滤波，得到一次扫描滤波后点云数据，并数据保存为ground1.las文件，具体过程如下：
步骤201、将一次激光雷达点云数据导入TerraSolid软件中，利用TerraSolid软件中的Classify by class工具，将一次激光雷达点云数据全部放入Default类中；
步骤202、采用计算机利用TerraSolid软件中的Classify low points工具，将噪声点从Default类中分离出至low points类中；其中，Classify low points工具中MAXcount阈值范围为1～6，Classify low points工具中Search设置为Single poits，Classify low points工具中高度差阈值范围为0.3m～1m，距离搜索范围阈值为2m～8m；
步骤203、采用计算机利用TerraSolid软件中的Classify isolated points工具，将孤立点从Default类中分离出至low points类中；其中，Classify isolated points工具中距离搜索范围阈值为2m～8m；
步骤204、采用计算机利用TerraSolid软件中的Classify ground工具，将地面点从Default类中分离出至Ground类中；其中，Classify ground工具中Max buliding size阈值为10m～50m，Iteration angle即迭代角度阈值为4°～12°，Iteration distance即迭代距离阈值为0.5m～1.5m，Terrain angle即地形角为88°；
步骤205、采用计算机利用TerraSolid软件中的Classify height from ground工具，将地面以下点从Default类中分离出至Ground类中；其中，Classify height fromground工具中Max height阈值为零，Min height阈值为-999m；
步骤206、采用计算机将Ground类中点记作一次扫描滤波后点云数据，并数据保存为ground1.las文件。
本实施例中，步骤七中得到矿区采煤塌陷盆地，之后，具体执行如下：
采用计算机利用ArcGIS软件中的识别工具，查看矿区采煤塌陷盆地中各个栅格像元的对应的高程差值；
采用计算机利用ArcGIS软件中的栅格转点工具，将测定矿区全区域差值栅格图转换为测定矿区全区域点要素，并将测定矿区全区域点要素的属性数据表导出为txt文件，得到测定矿区全区域点要素txt文件；其中，测定矿区全区域点要素txt文件包括采煤塌陷盆地栅格像元序号和差值栅格像元对应的高程差值。
本实施例中，无人机为SZT-R250无人机机载移动测量系统，操作简单且工作效率高。
本实施例中，计算机首先按照遍历间隔为α进行点集组的划分，得到多个点集组；然后对每个点集组按照选取比例为b％进行较低点集筛选，得到一次遍历间隔为α选取比例为b％的较低点集和二次遍历间隔为α选取比例为b％的较低点集，便于获取不同遍历间隔不同选取比例下的高程误差，便于得到应用于全局的最优较低点集选择所需的遍历间隔和选取比例，提高了全局较低点集获取的准确性，避免操作人员进行手动筛选，从而确保获取矿区采煤塌陷盆地的准确性。
本实施例中，根据获取的遍历间隔为α选取比例为b％时的各个差值栅格像元对应的高程差值，得到遍历间隔为α且选取比例为b％时的高程误差值，通过将多个高程误差值进行排序，得到最小高程误差值所对应的遍历间隔αmin，最小高程误差值所对应的选取比例，并将最小高程误差值所对应的遍历间隔αmin，最小高程误差值所对应的选取比例应用于一次全区域点云数据和二次全区域点云数据，从而获取一次全区域点云数据的较低点集和二次全区域点云数据的较低点集，提高了矿区采煤塌陷盆地测定的准确度，能够有效地避免了非地面点噪声和DEM建模与叠加算法而产生的误差。
综上所述，本发明方法步骤简单，设计合理，减少在人为干预的误差，处理过程简单，避免了非地面点噪声和DEM建模与叠加算法而产生的误差，确保矿区采煤塌陷盆地测定的准确性，实用性强。
以上所述，仅是本发明的较佳实施例，并非对本发明作任何限制，凡是根据本发明技术实质对以上实施例所作的任何简单修改、变更以及等效结构变化，均仍属于本发明技术方案的保护范围内。
Claims (4)
1.一种激光雷达扫描测定矿区采煤塌陷盆地的方法，其特征在于，该方法包括以下步骤：
步骤一、点云数据的获取：
步骤101、采用无人机搭载激光雷达对待测定矿区进行一次激光雷达扫描，获取一次激光雷达点云数据；其中，一次激光雷达点云数据包括第1个待测定矿区扫描点坐标数据、第2个待测定矿区扫描点坐标数据，...，第i个待测定矿区扫描点坐标数据，...，第n′个待测定矿区扫描点坐标数据；其中，i和n′均为正整数，n′≥100，且1≤i≤n′，第i个待测定矿区扫描点坐标数据记作Pi(xi,yi,zi)，xi表示第i个待测定矿区扫描点在WGS-84坐标系下的X坐标，yi表示第i个待测定矿区扫描点在WGS-84坐标系下的Y坐标，zi表示第i个待测定矿区扫描点在WGS-84坐标系下的Z坐标；
步骤102、采用无人机搭载激光雷达对矿区采煤塌陷盆地进行二次激光雷达扫描，获取二次激光雷达点云数据；
步骤二、一次激光雷达点云数据的滤波：
采用计算机利用TerraSolid软件，对一次激光雷达点云数据进行滤波，得到一次扫描滤波后点云数据，并将一次扫描滤波后点云数据保存为ground1.las文件；
步骤三、二次激光雷达点云数据的滤波：
采用计算机利用TerraSolid软件，对二次激光雷达点云数据进行滤波，得到二次扫描滤波后点云数据，并将二次扫描滤波后点云数据保存为ground2.las文件；
步骤四、选取参照区及参照区域点云数据的导出：
步骤401、采用计算机将ground1.las文件和ground2.las文件导入到TerraSolid软件，得到点云覆盖区域；
步骤402、采用计算机利用TerraSolid软件中Place Block工具，并将Place Block工具中Method设置为Rotated，在点云覆盖区域外画矩形ABCD，并将AB边的中点标记为a，BC边的中点标记为b，CD边的中点标记为c，DA边的中点标记为d；
步骤403、采用计算机利用TerraSolid软件中Draw VerticalSection工具设置剖面的宽度为1m，分别沿AC、BD、ac和bd方向画剖面线，得到第一剖面图、第二剖面图、第三剖面图和第四剖面图；
步骤404、采用计算机利用TerraSolid软件中Place Block工具，将Place Block工具中Method设置为Orthogonal，将第一剖面图中稳定区域框选，得到第一稳定区域和第二稳定区域；
采用计算机利用TerraSolid软件中Place Block工具将第二剖面图中稳定区域框选，得到第三稳定区域和第四稳定区域；
采用计算机利用TerraSolid软件中Place Block工具将第三剖面图中稳定区域框选，得到第五稳定区域和第六稳定区域；
采用计算机利用TerraSolid软件中Place Block工具将第四剖面图中稳定区域框选，得到第七稳定区域和第八稳定区域；其中，第一稳定区域、第二稳定区域、第三稳定区域、第四稳定区域、第五稳定区域、第六稳定区域、第七稳定区域和第八稳定区域中一次扫描滤波后点云数据和二次扫描滤波后点云数据的高程相同；
步骤405、采用计算机利用TerraSolid软件中Place Point or Stream Curve工具将第一稳定区域、第二稳定区域、第三稳定区域、第四稳定区域、第五稳定区域、第六稳定区域、第七稳定区域和第八稳定区域靠近点云覆盖区域中心的内边界点用曲线连接；其中，曲线内区域为非稳定区域，曲线外与外界矩形ABCD之间区域为稳定区域；
步骤406、采用计算机利用TerraSolid软件中Place Tile工具从稳定区域中框选四个参照区；其中，每个参照区为32m×32m的正方形区域，且四个参照区分别为1参照区、2参照区、3参照区和4参照区，1参照区位于A顶点、d中点、a中点和曲线之间围成的区域，2参照区位于B顶点、a中点、b中点和曲线之间围成的区域，3参照区位于C顶点、b中点、c中点和曲线之间围成的区域，4参照区位于D顶点、c中点、d中点和曲线之间围成的区域；
步骤407、采用计算机利用TerraSolid软件将一次扫描滤波后点云数据和二次扫描滤波后点云数据关闭，将一次扫描滤波后点云数据ground1.las文件加载，利用TerraSolid软件中Delete point工具中的Outside fence工具，分别得到1参照区中点云数据、2参照区中点云数据、3参照区中点云数据和4参照区中点云数据，并将1参照区中点云数据、2参照区中点云数据、3参照区中点云数据和4参照区中点云数据分别记作一次第一组点云数据、一次第二组点云数据、一次第三组点云数据和一次第四组点云数据；
步骤408、采用计算机利用TerraSolid软件将一次扫描滤波后点云数据关闭，将二次扫描滤波后点云数据ground2.las文件加载，利用TerraSolid软件中Delete point工具中的Outside fence工具，分别得到1参照区中点云数据、2参照区中点云数据、3参照区中点云数据和4参照区中点云数据，并将1参照区中点云数据、2参照区中点云数据、3参照区中点云数据和4参照区中点云数据分别记作二次第一组点云数据、二次第二组点云数据、二次第三组点云数据和二次第四组点云数据；
步骤五、参照区较低点集的获取：
步骤501、采用计算机并将一次第一组点云数据按照X坐标从小到大，Y坐标从小到大顺序进行排序；
步骤502、采用计算机将一次第一组点云数据中X坐标最小值记作
一次第一组点云数据中X坐标最大值记作
一次第一组点云数据中Y坐标最小值记作
一次第一组点云数据中Y坐标最大值记作
采用计算机对一次第一组点云数据进行遍历，将一次第一组点云数据划分为多个遍历间隔为α的点集组；其中，任一个遍历间隔为α的点集组中点的X坐标在
区间，且第n+1个遍历间隔为α的点集组中点的Y坐标在
区间，n和m均为自然数，且n的取值范围为0≤n≤N1-1，m的取值范围为0≤m≤N2-1，且
α表示一次第一组点云数据的遍历间隔，ceil(·)为向上取整函数；
步骤503、采用计算机将多个遍历间隔为α的点集组分别记作第1个遍历间隔为α的点集组，第2个遍历间隔为α的点集组，...，第j个遍历间隔为α的点集组，...，第M个遍历间隔为α的点集组，j和M为正整数，且1≤j≤M；其中，M表示遍历间隔为α的点集组的总数；
步骤504、采用计算机从第j个遍历间隔为α的点集组中提取较低点集的过程如下：
步骤5041、采用计算机将第j个遍历间隔为α的点集组中点集按照Z坐标从小到大进行排序；
步骤5042、采用计算机从第j个遍历间隔为α的点集组中筛选出Z坐标较小点集；其中，Z坐标较小点集中点的数量占第j个遍历间隔为α的点集组中点的总数的比例为b％；
步骤5043、采用计算机将从第j个遍历间隔为α的点集组中选出的Z坐标较小点集记作第j个遍历间隔为α选取比例为b％的较低点集；
步骤5044、多次重复步骤5041至步骤5043，直至完成第M个遍历间隔为α的点集组的较低点集的筛选，得到一次第一组遍历间隔为α选取比例为b％的较低点集；
步骤505、按照步骤501至步骤504所述的方法得到一次第二组遍历间隔为α选取比例为b％的较低点集、一次第三组遍历间隔为α选取比例为b％的较低点集和一次第四组遍历间隔为α选取比例为b％的较低点集；
步骤506、按照步骤501至步骤505所述的方法，得到二次第一组遍历间隔为α选取比例为b％的较低点集、二次第二组遍历间隔为α选取比例为b％的较低点集、二次第三组遍历间隔为α选取比例为b％的较低点集和二次第四组遍历间隔为α选取比例为b％的较低点集；
步骤六、参照区栅格DEM高程误差的计算：
步骤601、采用计算机利用ArcGIS软件中的添加XY数据工具，将一次第一组遍历间隔为α选取比例为b％的较低点集坐标导入；
步骤602、采用计算机利用ArcGIS软件中的自然邻域法插值工具，得到一次第一组遍历间隔为α选取比例为b％的点集数据栅格图；其中，输出像元大小为1m，一次第一组遍历间隔为α选取比例为b％的点集数据栅格图的尺寸为32m×32m；
步骤603、重复步骤601和步骤602，对二次第一组遍历间隔为α的点集组b％较低点集处理，得到二次第一组遍历间隔为α选取比例为b％的点集数据栅格图；
步骤604、采用计算机利用ArcGIS软件中的栅格计算工具，将二次第一组遍历间隔为α选取比例为b％的点集数据栅格图与一次第一组遍历间隔为α选取比例为b％的点集数据栅格图相减，得到遍历间隔为α选取比例为b％时的第一组差值栅格图；
步骤605、采用计算机利用ArcGIS软件中的栅格转点工具，将遍历间隔为α选取比例为b％时的第一组差值栅格图转换成遍历间隔为α选取比例为b％时的第一组差值点要素，并将遍历间隔为α选取比例为b％时的第一组差值点要素属性表导出为txt文件，得到遍历间隔为α选取比例为b％时的第一组差值点要素txt文件；其中，遍历间隔为α选取比例为b％时的第一组差值点要素txt文件包括差值栅格像元序号和差值栅格像元对应的高程差值；
步骤606、按照步骤601至步骤605所述的方法，对一次第二组遍历间隔为α选取比例为b％的点集数据栅格图和二次第二组遍历间隔为α选取比例为b％的点集数据栅格图处理，得到遍历间隔为α选取比例为b％时的第二组差值点要素txt文件；
对一次第三组遍历间隔为α选取比例为b％的点集数据栅格图和二次第三组遍历间隔为α选取比例为b％的点集数据栅格图处理，得到遍历间隔为α选取比例为b％时的第三组差值点要素txt文件；
对一次第四组遍历间隔为α选取比例为b％的点集数据栅格图和二次第四组遍历间隔为α选取比例为b％的点集数据栅格图处理，得到遍历间隔为α选取比例为b％时的第四组差值点要素txt文件；
步骤607、采用计算机根据公式
得到遍历间隔为α且选取比例为b％时的高程误差值；其中，
表示遍历间隔为α选取比例为b％时的第一组差值点要素txt文件中第f个差值栅格像元对应的高程差值，
表示遍历间隔为α选取比例为b％时的第二组差值点要素txt文件中第f个差值栅格像元对应的高程差值，
表示遍历间隔为α选取比例为b％时的第三组差值点要素txt文件中第f个差值栅格像元对应的高程差值，
表示遍历间隔为α选取比例为b％时的第四组差值点要素txt文件中第f个差值栅格像元对应的高程差值；其中，f为正整数，且1≤f≤1024；
步骤608、采用计算机按照步骤601至步骤607所述的方法，当α取1m、2m、4m、8m、16m或者32m，当b％取5％、10％或者15％时，分别得到遍历间隔为1m且选取比例为5％时的高程误差值e1,5、遍历间隔为1m且选取比例为10％时的高程误差值e1,10、遍历间隔为1m且选取比例为15％时的高程误差值e1,15；遍历间隔为2m且选取比例为5％时的高程误差值e2,5、遍历间隔为2m且选取比例为10％时的高程误差值e2,10、遍历间隔为2m且选取比例为15％时的高程误差值e2,15；遍历间隔为4m且选取比例为5％时的高程误差值e4,5、遍历间隔为4m且选取比例为10％时的高程误差值e4,10、遍历间隔为4m且选取比例为15％时的高程误差值e4,15；遍历间隔为8m且选取比例为5％时的高程误差值e8,5、遍历间隔为8m且选取比例为10％时的高程误差值e8,10、遍历间隔为8m且选取比例为15％时的高程误差值e8,15；遍历间隔为16m且选取比例为5％时的高程误差值e16,5、遍历间隔为16m且选取比例为10％时的高程误差值e16,10、遍历间隔为16m且选取比例为15％时的高程误差值e16,15；遍历间隔为32m且选取比例为5％时的高程误差值e32,5、遍历间隔为32m且选取比例为10％时的高程误差值e32,10、遍历间隔为32m且选取比例为15％时的高程误差值e32,15；
步骤609、采用计算机将步骤608中多个高程误差值按照从小到大顺序进行排序，得到最小高程误差值，并将最小高程误差值所对应的遍历间隔记作遍历间隔αmin，最小高程误差值所对应的选取比例记作选取比例bmin％；
步骤七、全区域点云数据处理及矿区采煤塌陷盆地的测定：
步骤701、采用计算机将ground1.las文件中点云数据按照X坐标从小到大，Y坐标从小到大顺序进行排序；将ground2.las文件中点云数据按照X坐标从小到大，Y坐标从小到大顺序进行排序；
步骤702、当选择遍历间隔为αmin且选取比例为bmin％，按照步骤502至步骤5044，对一次扫描滤波后点云数据处理，得到一次全区域点云数据的较低点集；
当选择遍历间隔为αmin且选取比例为bmin％，按照步骤502至步骤5044，对二次扫描滤波后点云数据进行处理，得到二次全区域点云数据的较低点集；
步骤703、按照步骤601至步骤604所述的方法，得到测定矿区全区域差值栅格图，则测定矿区全区域差值栅格图为矿区采煤塌陷盆地，实现矿区采煤塌陷盆地的测定。
2.按照权利要求1所述的一种激光雷达扫描测定矿区采煤塌陷盆地的方法，其特征在于：步骤102中采用无人机搭载激光雷达对矿区采煤塌陷盆地进行激光扫描时，点云数据的间距小于0.15m，点云数据密度高于60个/平方米；且采用无人机搭载激光雷达对矿区采煤塌陷盆地进行激光扫描时，扫描区域为矩形；
采用无人机搭载激光雷达对待测定矿区进行一次激光雷达扫描和采用无人机搭载激光雷达对矿区采煤塌陷盆地进行二次激光雷达扫描之间的时间间隔为20days～60days。
3.按照权利要求1所述的一种激光雷达扫描测定矿区采煤塌陷盆地的方法，其特征在于：步骤二中采用计算机利用TerraSolid软件，对一次激光雷达点云数据进行滤波，得到一次扫描滤波后点云数据，并数据保存为ground1.las文件，具体过程如下：
步骤201、将一次激光雷达点云数据导入TerraSolid软件中，利用TerraSolid软件中的Classify by class工具，将一次激光雷达点云数据全部放入Default类中；
步骤202、采用计算机利用TerraSolid软件中的Classifylow points工具，将噪声点从Default类中分离出至low points类中；其中，Classify low points工具中MAX count阈值范围为1～6，Classify low points工具中Search设置为Single poits，Classify lowpoints工具中高度差阈值范围为0.3m～1m，距离搜索范围阈值为2m～8m；
步骤203、采用计算机利用TerraSolid软件中的Classifyisolated points工具，将孤立点从Default类中分离出至low points类中；其中，Classify isolated points工具中距离搜索范围阈值为2m～8m；
步骤204、采用计算机利用TerraSolid软件中的Classify ground工具，将地面点从Default类中分离出至Ground类中；其中，Classify ground工具中Max buliding size阈值为10m～50m，Iteration angle即迭代角度阈值为4°～12°，Iteration distance即迭代距离阈值为0.5m～1.5m，Terrain angle即地形角为88°；
步骤205、采用计算机利用TerraSolid软件中的Classify height from ground工具，将地面以下点从Default类中分离出至Ground类中；其中，Classify height from ground工具中Max height阈值为零，Min height阈值为-999m；
步骤206、采用计算机将Ground类中点记作一次扫描滤波后点云数据，并数据保存为ground1.las文件。
4.按照权利要求1所述的一种激光雷达扫描测定矿区采煤塌陷盆地的方法，其特征在于：步骤七中得到矿区采煤塌陷盆地，之后，具体执行如下：
采用计算机利用ArcGIS软件中的识别工具，查看矿区采煤塌陷盆地中各个栅格像元的对应的高程差值；
采用计算机利用ArcGIS软件中的栅格转点工具，将测定矿区全区域差值栅格图转换为测定矿区全区域点要素，并将测定矿区全区域点要素的属性数据表导出为txt文件，得到测定矿区全区域点要素txt文件；其中，测定矿区全区域点要素txt文件包括采煤塌陷盆地栅格像元序号和差值栅格像元对应的高程差值。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201911018042.7A CN110726998B (zh) | 2019-10-24 | 2019-10-24 | 一种激光雷达扫描测定矿区采煤塌陷盆地的方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201911018042.7A CN110726998B (zh) | 2019-10-24 | 2019-10-24 | 一种激光雷达扫描测定矿区采煤塌陷盆地的方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN110726998A CN110726998A (zh) | 2020-01-24 |
| CN110726998B true CN110726998B (zh) | 2020-08-07 |
Family
ID=69223091
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201911018042.7A Active CN110726998B (zh) | 2019-10-24 | 2019-10-24 | 一种激光雷达扫描测定矿区采煤塌陷盆地的方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN110726998B (zh) |
Families Citing this family (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111583302B (zh) * | 2020-05-06 | 2021-07-16 | 北京大学 | 一种基于三维激光点云的割煤顶板线提取方法 |
| CN111895911B (zh) * | 2020-08-01 | 2022-02-22 | 上海市地矿工程勘察(集团)有限公司 | 一种应用于浅部砂层地面塌陷隐患监测方法 |
| CN113086876B (zh) * | 2021-04-06 | 2022-06-28 | 上海建工五建集团有限公司 | 一种多层卷筒钢丝绳传动方式的卷筒钢丝绳故障检测方法 |
| CN115713531B (zh) * | 2023-01-05 | 2023-05-09 | 山东环宇地理信息工程有限公司 | 一种基于InSAR的地球表面图像数据处理系统 |
Citations (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101914890A (zh) * | 2010-08-31 | 2010-12-15 | 中交第二公路勘察设计研究院有限公司 | 一种基于机载激光测量的公路改扩建勘测方法 |
| CN103679655A (zh) * | 2013-12-02 | 2014-03-26 | 河海大学 | 一种基于坡度与区域生长的LiDAR点云滤波方法 |
| CN106952242A (zh) * | 2016-01-06 | 2017-07-14 | 北京林业大学 | 一种基于体素的渐进不规则三角网点云滤波方法 |
| CN106990401A (zh) * | 2017-05-24 | 2017-07-28 | 武汉大学 | 基于全波形机载激光雷达数据二类高程误差修正方法 |
| CN107092020A (zh) * | 2017-04-19 | 2017-08-25 | 北京大学 | 融合无人机LiDAR和高分影像的路面平整度监测方法 |
| CN107101683A (zh) * | 2017-06-15 | 2017-08-29 | 西安科技大学 | 一种基于激光雷达与速度信息的煤流量监测系统 |
| CN109143257A (zh) * | 2018-07-11 | 2019-01-04 | 中国地质调查局西安地质调查中心 | 无人机机载雷达矿山开采土地变化监测系统及方法 |
| CN109407113A (zh) * | 2018-11-19 | 2019-03-01 | 中南林业科技大学 | 一种基于机载激光雷达的林窗时空变化监测与量化方法 |
| CN109410265A (zh) * | 2019-01-22 | 2019-03-01 | 江苏省测绘工程院 | 一种基于往期dem辅助的tin滤波改进算法 |
Family Cites Families (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US8458188B2 (en) * | 2010-02-17 | 2013-06-04 | Lockheed Martin Corporation | Voxel approach to terrain repositories for modeling and simulation |
| CN102662179A (zh) * | 2012-05-18 | 2012-09-12 | 四川省科学城久利科技实业有限责任公司 | 基于机载激光雷达的三维优化选线方法 |
| CN104502919A (zh) * | 2015-01-13 | 2015-04-08 | 南京大学 | 利用机载激光雷达点云提取城市植被三维覆盖的方法 |
| US10127685B2 (en) * | 2015-12-16 | 2018-11-13 | Objectvideo Labs, Llc | Profile matching of buildings and urban structures |
-
2019
- 2019-10-24 CN CN201911018042.7A patent/CN110726998B/zh active Active
Patent Citations (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101914890A (zh) * | 2010-08-31 | 2010-12-15 | 中交第二公路勘察设计研究院有限公司 | 一种基于机载激光测量的公路改扩建勘测方法 |
| CN103679655A (zh) * | 2013-12-02 | 2014-03-26 | 河海大学 | 一种基于坡度与区域生长的LiDAR点云滤波方法 |
| CN106952242A (zh) * | 2016-01-06 | 2017-07-14 | 北京林业大学 | 一种基于体素的渐进不规则三角网点云滤波方法 |
| CN107092020A (zh) * | 2017-04-19 | 2017-08-25 | 北京大学 | 融合无人机LiDAR和高分影像的路面平整度监测方法 |
| CN106990401A (zh) * | 2017-05-24 | 2017-07-28 | 武汉大学 | 基于全波形机载激光雷达数据二类高程误差修正方法 |
| CN107101683A (zh) * | 2017-06-15 | 2017-08-29 | 西安科技大学 | 一种基于激光雷达与速度信息的煤流量监测系统 |
| CN109143257A (zh) * | 2018-07-11 | 2019-01-04 | 中国地质调查局西安地质调查中心 | 无人机机载雷达矿山开采土地变化监测系统及方法 |
| CN109407113A (zh) * | 2018-11-19 | 2019-03-01 | 中南林业科技大学 | 一种基于机载激光雷达的林窗时空变化监测与量化方法 |
| CN109410265A (zh) * | 2019-01-22 | 2019-03-01 | 江苏省测绘工程院 | 一种基于往期dem辅助的tin滤波改进算法 |
Non-Patent Citations (3)
| Title |
|---|
| Remote sensing monitoring of surface vegetation and soil moisture changes and the disturbance effect of coal mining subsidence in the Western mining area of China;TANG Fuquan and GU Jin;《2018 Fifth International Workshop on Earth Observation and Remote Sensing Applications》;20181231;正文第1-5页 * |
| 基于沉陷对称特征的近水平煤层开采InSAR三维位移反演模型;汤伏全 等;《煤炭学报》;20190131;第210-220页 * |
| 机载激光雷达森林冠层高度模型凹坑去除方法;段祝庚 等;《农业工程学报》;20141130;第209-217页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN110726998A (zh) | 2020-01-24 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN110726998B (zh) | 一种激光雷达扫描测定矿区采煤塌陷盆地的方法 | |
| CN107451982B (zh) | 一种基于无人机影像的高郁闭度林分树冠面积获取方法 | |
| CN103020342B (zh) | 一种从地面LiDAR数据中提取建筑物轮廓和角点的方法 | |
| CN106918311A (zh) | 基于车载激光点云数据的单株树树冠投影面积自动计算方法 | |
| CN113537141B (zh) | 一种堤坝管涌和滑坡病害快速检测方法及系统 | |
| CN104751479A (zh) | 基于tin数据的建筑物提取方法和装置 | |
| CN109636904A (zh) | 一种基于uav航测地形数据的噪声处理技术 | |
| Salah et al. | Evaluation of the self‐organizing map classifier for building detection from lidar data and multispectral aerial images | |
| CN114119902A (zh) | 一种基于无人机倾斜三维模型的建筑物提取方法 | |
| CN115564926A (zh) | 基于影像建筑物结构学习的三维面片模型构建方法 | |
| CN109961512B (zh) | 地形机载点云提取方法及装置 | |
| CN114862715A (zh) | 一种融合地形特征语义信息的tin渐进加密去噪方法 | |
| CN115131231A (zh) | 辅以多特征聚类的复杂地形区点云层次滤波方法 | |
| KR101078238B1 (ko) | 항공 레이저 측량의 점군 집합을 활용한 하천 지역 내 3차원 인공 제방선 추출 방법 | |
| CN116721228B (zh) | 一种基于低密度点云的建筑物高程提取方法及系统 | |
| CN111895907B (zh) | 一种电塔点云提取方法、系统及设备 | |
| Kim et al. | Adaptive morphological filtering for DEM generation | |
| CN115661398A (zh) | 一种用于实景三维模型的建筑物提取方法、装置及设备 | |
| Sajadian et al. | A data driven method for building reconstruction from LiDAR point clouds | |
| Mu et al. | Canopy lidar point cloud data k-means clustering watershed segmentation method | |
| Kim et al. | Building footprints extraction of dense residential areas from Lidar data | |
| CN113570621A (zh) | 一种基于高精度点云与影像的树木信息提取方法及装置 | |
| CN113496551B (zh) | 一种基于地质露头三维模型的地形剖面线绘制方法 | |
| CN117523124B (zh) | 一种基于无人机激光雷达的大流域地形变化切片方法 | |
| CN115031689B (zh) | 一种基于激光点云数据的电力输电塔倾斜状态识别方法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| CB03 | Change of inventor or designer information | ||
| CB03 | Change of inventor or designer information |
Inventor after: Tang Fuquan Inventor after: Yang Qian Inventor after: Liu Shiwei Inventor after: He Kelu Inventor before: Tang Fuquan Inventor before: Yang Qian Inventor before: He Kelu Inventor before: Liu Shiwei |
|
| GR01 | Patent grant | ||
| GR01 | Patent grant |