CN111209828B - 从机载激光雷达点云中提取建筑物屋顶点的方法和系统 - Google Patents

Abstract

本发明公开了一种从机载激光雷达点云中提取建筑物屋顶点的方法和系统，其中提取建筑物屋顶点的方法通过三维欧几里得聚类将测区的机载激光雷达点云分割为多个聚类，对其中面积较大的聚类，计算DSM的拉普拉斯算子，如果拉普拉斯算子小于平坦度阈值的像素占比较高，将聚类内的点分为建筑物屋顶点；对面积较小的聚类，采用改进的RANSAC算法提取建筑物屋顶上的点。该方法能够克服提取建筑物平面屋顶点时计算复杂度高、执行效率低的问题，还可以避免将非建筑物的树木点分类为建筑物屋顶点，提高了建筑物屋顶点提取的效率和准确度。

Description

从机载激光雷达点云中提取建筑物屋顶点的方法和系统

技术领域

本发明属于机载激光雷达点云数据处理技术领域，具体涉及一种从机载激光雷达点云中提取落在建筑物屋顶上的点的方法和系统。

背景技术

机载LiDAR(Light Detection And Ranging，激光雷达)是当前测绘领域效率最高、发展最快的大面积测绘手段之一。通过使用激光器发射和接收高能激光脉冲来测距，GNSS(Global Navigation Satellite System，全球导航卫星系统)接收机给出激光器实时位置，INS(Inertial Navigation System，惯性导航系统)给出激光器实时三维姿态，可利用向量公式和坐标转换计算出散射面的三维坐标。借助机械扫描装置，机载LiDAR可以实现大面积面状测量。目前，机载LiDAR的发射频率普遍高达数百KHZ以上乃至2000KHZ。这意味着，在不考虑多回波和无回波的情况下，一秒钟激光雷达可采集数十万乃至数百万点。因此，机载LiDAR可以在短时间内获得海量的点，称为点云。

经过初步处理的机载LiDAR点云中，除点坐标外，还提供的点属性包括回波强度、第几次回波、数据采集时间等。但是，缺乏语义信息，不能给出激光脉冲打到的散射面的物理性质，不知道测量是地面、建筑物还是植被、鸟等。而在构建DEM(Digital ElevationModel，数字高程模型)时必须采用位于地面上的点，在构建三维建筑物模型时必须采用位于建筑物上的点，在进行森林制图时必须采用位于树木上的点。所以，对点云进行分类，确定点位于哪种基本地物类型上，是机载LiDAR数据处理中最重要的基础工作之一。

在机载LiDAR点云，尤其是覆盖居民区的点云中，往往含有落在建筑物上的点。提取出落在建筑物上的点，可以用于建筑物的三维建模、DSM(Digital Surface Model，数字表面模型)生成、城市规划等。建筑物点提取是点云数据处理的热门问题。

目前研究人员已经提出了很多从机载LiDAR点云中提取建筑物点的方法，既有基于点云的，也有基于图像的。基于图像的方法(参见徐文学，杨必胜，魏征，等.多标记点过程的LiDAR点云数据建筑物和树冠提取[J].测绘学报，2013，42(1)：51-58.)将点云内插为图像，而后通常借助一些图像分析方法来提取含有建筑物的像素，进而实现建筑物点的提取和边界的追踪。基于点云的方法通常是利用RANSAC(随机采样一致性，参见专利CN105139379B)算法或Hough变换(参见袁晨鑫，基于Hough变换的建筑物点云数据特征提取及三维建模，2019，东华理工大学硕士论文)或基于种子点的区域生长算法(参加专利CN104036544B和CN 107944383A)从离散点云中提取位于局部平面上的点作为屋顶点。

但是，基于图像和基于点云的方法都有明显的问题。基于图像的方法往往会损失原始点云中的信息，尤其是点云的垂直分布所反映的地物垂直结构信息，对于提取树木遮挡的建筑物的效果并不好。另外人工修剪过的平顶植被容易被误分为屋顶点。基于区域生长的算法需要给出法向量夹角阈值和曲率阈值，而这两个阈值的取值受到数据质量和屋顶平整度的影响，难以简单确定。更重要的是，该方法提取到的屋顶点并不全，经常存在大量漏检的屋顶点。Hough变换提取平面需要三个参数，而三参数的Hough变换非常复杂，需要在三维Hough空间进行分块，用作累加器的分区，以进行投票。空间分块的粒度直接影响后续平面检测的效果，而且在检测结果中经常存在错误。在点云量大时计算复杂度很高，需要大量存储空间。基于RANSAC方法的基本思路是从点云中随机抽取3个点构建平面模型，测试其它点是否可以很好地用该平面进行拟合。上述过程重复执行很多次。假设点云共有N个点，则随机取3个点共有

个不重复的组合。随着N的增加，组合的数目迅速增加。如果N达到上万，

会是个天文数字。而随机采样的次数是有限的，这就导致找到稳健模型的概率明显降低。即使有限度地增加随机采样次数，也不会明显解决问题，反而会造成计算时间显著延长。从这个角度讲，RANSAC方法并不适合点数较多的情况。此外，从理论上讲，即使对于随机分布的点，RANSAC都可能得到一些模型。例如，森林内的点云，也可以用RANSAC探测出一些平面出来，而且平面上的点可能足够多。在这种情况下，传统RANSAC难以确保平面来自建筑物而不是非建筑物。

为了克服RANSAC的这些缺陷，迫切需要对RANSAC方法进行改进，使得其更适合对海量机载LiDAR点云进行平面提取，区分森林和建筑物屋顶，提高提取效率和准确率。

发明内容

发明目的：本发明旨在解决在机载LiDAR点云中使用RANSAC方法提取建筑物平面屋顶点时计算复杂度高、执行效率低的问题，避免将非建筑物的点分类为建筑物屋顶点，以提高建筑物屋顶点提取的效率和准确度。

技术方案：本发明一方面公开了一种从机载激光雷达点云中提取建筑物屋顶点的方法，包括：

(1)获取测区的机载激光雷达点云；对可能包含建筑物屋顶的点类进行三维欧几里得聚类；计算每个聚类内点云的平均离地高度，对平均离地高度大于高度阈值的每一个聚类按步骤2-4进行处理，提取建筑物屋顶点；

(2)生成当前聚类C_p对应的DSM，获取C_p在DSM中有效像素的数目N_valid，计算C_p在XY平面的投影面积A_C；

(3)如果C_p的投影面积A_C大于等于面积阈值，计算DSM中拉普拉斯算子L_p小于平坦度阈值L_T的像素占DSM有效像素总数N_valid的比值P_Laplacian，如果所述比值大于平顶比例阈值P_{Laplacian_T}，则当前聚类C_p中的点全部为建筑物屋顶上的点；

(4)如果C_p的投影面积A_C小于面积阈值，按如下步骤判断：

(4.1)平移当前聚类C_p得到C_tmp，C_tmp的重心在坐标原点；初始化当前采样次数N_S＝0；初始化C_tmp中所有点的状态为未抽取，清空有效模型序列和最优模型序列；

(4.2)如果当前采样次数N_S达到最大采样次数N_{S_T}，跳转至步骤4.5；

从C_tmp中随机抽取状态为未抽取的点p₀，将p₀的状态修改为已抽取，寻找C_tmp内距离p₀最近的三个邻域点p₁、p₂、p₃；如果p₁、p₂、p₃距p₀的距离都足够小，且四点能够用平面拟合，计算拟合平面P_fit；否则当前采样次数N_S加一，重新执行步骤4.2；

(4.3)计算C_tmp中每个点到平面P_fit的距离；将所述距离小于平坦度阈值L_T的点组成内点集合C_inlier；如果集合C_inlier中点数小于聚类最小点数阈值，当前采样次数N_S加一，重新执行步骤4.2；否则对集合C_inlier再次进行三维欧几里得聚类；

(4.4)如果步骤4.3只得到一个聚类，将平面P_fit加入有效模型序列；

当前采样次数N_S加一，重新执行步骤4.2；

(4.5)对有效模型序列进行遍历，寻找内点数目最多的有效模型作为最优模型，加入最优模型序列；从C_tmp中移除最优模型所对应的内点集合中的点，如果剩余点数大于聚类最小点数阈值，则初始化当前采样次数N_s＝0；初始化C_tmp中所有点的状态为未抽取，清空有效模型序列和最优模型序列，跳转至步骤4.2；

(4.6)计算C_p内的每个点与最优模型序列中平面的最小距离d_min，将d_min小于平坦度阈值L_T的点组成集合C_{inlier_all}；如果C_{inlier_all}的点数占C_p中点数的比例大于阈值P_{inlier_all_T}，则将C_{inlier_all}中的点分为建筑物屋顶点。

另一方面，本发明还公开了实现上述提取建筑物屋顶点的方法的系统，包括：

机载激光雷达点云聚类模块，用于获取测区的机载激光雷达点云；采用三维欧几里得聚类对所获取的点云进行聚类；计算每个聚类内点云的平均离地高度，获取平均离地高度大于高度阈值的聚类；DSM生成模块，用于对机载激光雷达点云聚类模块获取的每一个聚类生成对应的DSM，并获取所述聚类在DSM中的有效像素数目，计算所述聚类在XY平面的投影面积A_C；并根据A_C将聚类分为大面积聚类和小面积聚类；

大面积聚类分类模块，用于判断大面积聚类中的点是否为建筑物屋顶上的点：判断方法为：对大面积聚类计算DSM，统计拉普拉斯算子L_p小于平坦度阈值L_T的像素占有效像素总数N_valid的比例P_Laplacian，如果所述比值大于平顶比例阈值P_{Laplacian_T}，则将当前聚类C_p中的点全部分为建筑物屋顶点；

小面积聚类分类模块，用于提取小面积聚类中为建筑物屋顶上的点，提取方法为：

(4.1)平移当前小面积聚类C_p得到C_tmp，C_tmp的重心在坐标原点；初始化当前采样次数N_S＝0；初始化C_tmp中所有点的状态为未抽取，清空有效模型序列和最优模型序列；

(4.2)如果当前采样次数N_s达到最大采样次数N_{s_T}，跳转至步骤4.5；

从C_tmp中随机抽取状态为未抽取的点p₀，将p₀的状态修改为已抽取，寻找C_tmp内距离p₀最近的三个邻域点p₁、p₂、p₃；如果p₁、p₂、p₃距p₀的距离都足够小，且四点能够用平面拟合，计算拟合平面P_fit；否则当前采样次数N_S加一，重新执行步骤4.2；

(4.3)计算C_tmp中每个点到平面P_fit的距离；将所述距离小于平坦度阈值L_T的点组成内点集合C_inlier；如果集合C_inlier中点数小于聚类最小点数阈值，当前采样次数N_S加一，重新执行步骤4.2；否则对集合C_inlier再次进行三维欧几里得聚类；

(4.4)如果步骤4.3只得到一个聚类，将平面P_fit加入有效模型序列；

当前采样次数N_S加一，重新执行步骤4.2；

(4.5)对有效模型序列进行遍历，寻找内点数目最多的有效模型作为最优模型，加入最优模型序列；从C_tmp中移除最优模型所对应的内点集合中的点，如果剩余点数大于聚类最小点数阈值，则初始化当前采样次数N_S＝0；初始化C_tmp中所有点的状态为未抽取，清空有效模型序列和最优模型序列，跳转至步骤4.2；

(4.6)计算C_p内的每个点与最优模型序列中平面的最小距离d_min，将d_min小于平坦度阈值L_T的点组成集合C_{inlier_all}；如果C_{inlier_all}的点数占C_p中点数的比例大于阈值P_{inlier_all_T}，则C_{inlier_all}中的点为建筑物屋顶上的点。

有益效果：与现有技术相比，本发明公开的从机载激光雷达点云中提取建筑物屋顶点的方法具有以下优点：1、本发明采用三维欧几里得聚类，可以在大部分情况下有效地将建筑物屋顶与临近的树木分到不同的聚类，在提取到平面模型内点之后对平面模型内点进行欧几里得聚类可以将绝大部分非屋顶的树木点排除，提高了屋顶点提取的准确率；；2、本发明采用大小屋顶分而治之的策略，避免了对大屋顶也采用RANSAC进行平面检测带来的高计算量和潜在的漏检可能性，同时维持了高准确度；3、对小面积聚类，采用改进的RANSAC算法；设某聚类中点数为N，则潜在的平面模型数目从

降到N。如果N小于指定运行次数，则随机取点实际上只需要运行N次；4、本发明理论简单易行，计算效率高。对于300万点，一般可以在十秒内完成。可以多次运行，以更全面地提取屋顶点。

附图说明

图1为本发明公开的从机载激光雷达点云中提取建筑物屋顶点方法的流程图；

图2是欧几里得聚类的示意图；

图3是构建DSM的示意图；

图4是判断某个随机选取的点p₀是否和周围3个邻域点p₁、p₂、p₃可以近似地用平面拟合的示意图；

图5为本发明公开的从机载激光雷达点云中提取建筑物屋顶点的系统的组成示意图；

图6为本发明的实施效果图。

具体实施方式

下面结合附图和具体实施方式，进一步阐明本发明。

如图1所示，本发明公开了一种从机载激光雷达点云中提取建筑物屋顶点的方法，包括：

步骤1、获取测区的机载激光雷达点云；采用三维欧几里得聚类对所获取的点云进行聚类；计算每个聚类内点云的平均离地高度，对平均离地高度大于高度阈值的每一个聚类按步骤2-4进行处理，提取建筑物屋顶点；

如图2所示，欧几里得聚类的结果是使不同聚类之间的最小距离d大于等于搜索半径R。欧几里得聚类算法的参数包括：搜索半径R、聚类最小点数阈值N_min、聚类最大点数阈值N_max，这些参数值影响到聚类的结果。本发明中，上述参数的设置如下：

搜索半径R的范围为(R_min,R_max)，R最小值R_min与激光脉冲脚点密度负相关，本实施例中，

D是激光脉冲脚点的密度，可以由航飞设计方案获得。如果R小于R_min，可能会导致聚类点数过少，一个屋顶被分为多个聚类。R最大值R_max可以取相邻建筑物之间的最小距离，如1.5m，2.0m等。如果R超过R_max，可能会将多个屋顶、屋顶与树木等聚到一类，不能实现准确分割。

N_min、N_max与建筑物屋顶的最小、最大面积正相关。设建筑物屋顶最小面积为A_min，最大面积为A_max，则N_min、N_max的计算式为：

N_min＝A_minD-Q

N_max＝A_maxD+Q

其中Q为误差调节常数，本实施例中Q＝50。使用Q是为了容纳A_minD、A_maxD和屋顶真实的最小、最大点数之间的误差。

计算聚类内点云的平均离地高度，包括：

(1.1)获取聚类中每个点的地面高度；

一个点的地面高度可以通过对现有已知的地面点构建二维索引，进行最邻近查询获得，也可以通过查询其他来源的DEM数据获得。如果使用已知地面点计算，需要先对已知地面点利用其XY坐标构建二维索引T_grd。T_grd既可以是KD树索引，也可以是八叉树索引。设待计算的聚类中第n个点的坐标为(X_n,Y_n,H_n)，通过使用T_grd查询XY平面上距离(X_n,Y_n)最近的地面点G的序号，进而得到G的三维坐标，将G的高程作为聚类中第n个点对应的地面高度H_Tn；

(1.2)计算聚类中每个点的离地高度：△H_n＝H_n-H_Tn；

聚类内点云的平均离地高度：

N为聚类内点的总数。

鉴于建筑物屋顶一般至少离地有一定高度，比如2.5m以上，所以如果ΔH_mean小于给定的高度阈值ΔH_{mean_T}，则认为聚类不代表建筑物，如果ΔH_mean大于等于ΔH_{mean_T}，按步骤2-4进行处理，确定其中属于建筑物屋顶上的点；

本实施例中，ΔH_{mean_T}取值为2.5m。

步骤2、生成当前聚类C_p对应的DSM，获取C_p在DSM中有效像素的数目N_valid，计算C_p在XY平面的面积A_c；具体步骤为：

(2.1)计算C_p中点的最小X坐标X_min、最大X坐标X_max、最小Y坐标Y_min、最大Y坐标Y_max；将以(X_min,Y_min)、(X_max,Y_max)为对角顶点的矩形区域划分为N_row行、N_col列、格网大小为R的栅格，其中：

N_row＝int(Y_max-Y_min)/R+1

N_col＝int(X_max-X_min)/R+1

R为三维欧几里得聚类的搜索半径，int(·)为取整函数；

(2.2)将C_p中的点根据其X、Y坐标分配到栅格中的像素，得到当前聚类C_p对应的DSM；

第n个点映射到像素(i,j)中，

其中(X_n,Y_n,H_n)为第n个点的坐标；

(2.3)对栅格内的像素进行遍历；如果一个像素有点落入，则标识该像素为有效像素，取落入点的高度指标作为DSM的值，无点落入的像素的DSM值设为H_invalid；如图3所示，白色像素为无点落入的像素，灰色像素为有点落入的像素；所述高度指标为：落入点的高程最小值，或落入点的高程最大值，或落入点的高程平均值；

(2.4)对DSM逐像素进行遍历，统计有效像素的数目N_valid；C_p在XY平面的投影面积A_C为：A_c＝N_validR²。

如果A_C大于等于指定的面积阈值A_{C_T}，则认为认为C_p为大面积聚类；否则C_p为小面积聚类。

步骤3、对于大面积聚类，计算DSM中拉普拉斯算子L_p小于平坦度阈值L_T的像素占有效像素总数N_valid的比值P_Laplacian，如果所述比值大于平顶比例阈值P_{Laplacian_T}，则当前聚类C_p中的点全部为建筑物屋顶上的点；

DSM在像素(i,j)处的拉普拉斯算子L_p(i,j)为：

平坦度阈值L_T取值与点云的相对误差正相关，本实施例中取0.10-0.20m。设拉普拉斯算子L_p值小于L_T的像素共N_plane个，则：

如果P_Laplacian大于等于给定平顶比例阈值P_{Laplacian_T}，则认为C_p是一个由平顶构成的建筑物，将C_p内的点全部分为建筑物点。面积越大、结构越简单的屋顶，P_{Laplacian_T}可以取得越大，本实施例中取值为0.6。

步骤4、如果聚类投影区域的面积A_C小于面积阈值，采用改进的RANSAC算法判断C_p中点的类别，步骤如下：

(4.1)平移当前聚类C_p得到C_tmp，使C_tmp的重心在坐标原点；

计算聚类C_p的重心O_p(X_Op,Y_Op,H_Op)：

C_tmp中第n个点的坐标为：(X_n-X_Op,Y_n-Y_Op,H_n-H_Op)，则C_tmp的重心在坐标原点；

初始化当前采样次数N_S＝0；初始化C_tmp中所有点的状态为未抽取，清空有效模型序列和最优模型序列；

(4.2)如果当前采样次数N_S达到最大采样次数N_{S_T}，跳转至步骤4.5；最大采样次数N_{s_T}为预设的值，其取值最大不超过当前聚类中的点数N，本实施例中取值为100；

从C_tmp中随机抽取状态为未抽取的点p₀，将p₀的状态修改为已抽取，寻找C_tmp内距离p₀最近的三个邻域点p₁、p₂、p₃；如果p₁、p₂、p₃距p₀的距离都足够小，且四点能够用平面拟合，计算拟合平面P_fit；否则当前采样次数N_S加一，重新执行步骤4.2；

本实施例中，判断p₁、p₂、p₃和p₀的欧几里得距离是否均小于R，如果否，则认为距离不够小，邻域关系不够稳健，这四个点不能用于稳健地估计局部平面；如果是，判断四个点是否可以近似地用平面进行拟合，判断的方法既可以采用最小二乘拟合出平面再判断p₁、p₂、p₃到平面的最大距离是否够小，也可以基于平面夹角进行判断。鉴于基于平面夹角的方法简单易行效率高，下面结合图4介绍该方法。

p₁、p₂、p₃和p₀构成三个三角形，即Δp₀p₁p₂、Δp₀p₂p₃、Δp₀p₃p₁。计算这些三角形的法向量，Δp₀p₁p₂的法向量N₀₁₂为：N₀₁₂＝(p₁-p₀)×(p₂-p₀)；对N₀₁₂进行模标准化。Δp₀p₂p₃的法向量N₀₂₃、Δp₀p₃p₁的法向量N₀₃₁采用类似的方法进行计算和标准化。

计算N₀₁₂和N₀₂₃的夹角θ₀：θ₀＝cos^-1(N₀₁₂·N₀₂₃)；

N₀₁₂和N₀₃₁的夹角θ₁、N₀₂₃和N₀₃₁的夹角θ₂采用类似的方法进行计算。cos^-1为反余弦函数。

如果θ₀、θ₁、θ₂均接近0°或180°，认为p₀、p₁、p₂和p₃可以用某平面P_fit进行拟合；否则认为四个点不能用平面进行拟合；

计算P_fit的参数，设P_fit的方程为：ax+by+ch+d＝0；

a、b、a分别为平面法向量在三个坐标轴上的分量，d＝-(ax₀+by₀+ch₀)，(x₀,y₀,h₀)为p₀的坐标。

(4.3)计算C_tmp中每个点到平面P_fit的距离，其中第n个点到P_fit的距离d_fitn为：

(x_tmpn,y_tmpn,h_tmpn)为C_tmp中第n个点的坐标；

将所述距离小于平坦度阈值L_T的点组成内点集合C_inlier；如果集合C_inlier中点数小于聚类最小点数阈值N_min，当前采样次数N_s加一，重新执行步骤4.2；否则对集合C_inlier再次进行三维欧几里得聚类，此处欧几里得聚类的参数设置如下：

搜索半径为R，聚类最小点数阈值为1，聚类最大点数阈值为集合C_inlier中的总点数N_inlier；

(4.4)如果步骤4.3只得到一个聚类，则平面P_fit为有效模型，将其加入有效模型序列；

当前采样次数N_s加一，重新执行步骤4.2；

(4.5)对有效模型序列进行遍历，寻找内点数目最多的有效模型作为最优模型，加入最优模型序列；如果内点数目最多的模型有多个，则分别计算集合C_inlier中的内点到这几个模型的距离的标准差，取其中标准差最小的模型为最优模型，加入最优模型序列；

从C_tmp中移除最优模型所对应的内点集合中的点，如果剩余点数大于聚类最小点数阈值N_min，则初始化当前采样次数N_s＝0；初始化C_tmp中所有点的状态为未抽取，清空有效模型序列和最优模型序列，跳转至步骤4.2；

(4.6)计算C_p内的每个点与最优模型序列中平面的最小距离d_min，将d_min小于平坦度阈值L_T的点组成集合C_{inlier_all}；如果C_{inlier_all}的点数占C_p中点数的比例大于阈值P_{inlier_all_T}，则认为C_{inlier_all}为建筑物屋顶，将其中的点分为建筑物屋顶点。本实施例中，P_{inlier_all_T}的值为0.50。

跳转至步骤2对下一个聚类进行判断处理，直至步骤1中划分的所有聚类都处理完毕。

本实施例还公开了实现上述从机载激光雷达点云中提取建筑物屋顶点方法的系统，如图5所示，包括：

机载激光雷达点云聚类模块，用于根据步骤1获取测区的机载激光雷达点云；采用三维欧几里得聚类对所获取的点云进行聚类；计算每个聚类内点云的平均离地高度，获取平均离地高度大于高度阈值的聚类；

DSM生成模块，用于根据步骤2对机载激光雷达点云聚类模块获取的每一个聚类生成对应的DSM，并获取所述聚类在DSM中的有效像素数目，计算所述聚类在XY平面的投影面积；

大面积聚类分类模块，用于根据步骤3判断大面积聚类中的点是否为建筑物屋顶点；

小面积聚类分类模块，用于根据步骤4判断小面积聚类中的点是否为建筑物屋顶点。

图6为本发明的实施效果图，其中图6-(a)为某地的高分辨率光学影像，图6-(b)为输入的点云数据，其中点云已经进行了初步分类，图6-(c)为本发明实施分类后的点云数据。对比图6-(a)和6-(c)，可以看到本发明的方法对大建筑物和小建筑物、平顶和非平顶、是否有树木遮挡的建筑物屋顶点都有良好的提取效果。

本发明提供了一种基于分而治之策略从机载LiDAR点云中提取建筑物屋顶点的方法，尤其是利用改进的RANSAC方法快速从小面积聚类点云中提取平面以及对模型内点使用聚类数目来判断是否为有效的模型。具体实现该技术方案的方法和途径很多，以上所述仅是本发明的优选实施方式。应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰。这些改进和润饰也应视为本发明的保护范围。本实施例中未明确的各组成部分均可用现有技术加以实现。

Claims (10)

1.从机载激光雷达点云中提取建筑物屋顶点的方法，其特征在于，包括：

(1)获取测区的机载激光雷达点云；采用三维欧几里得聚类对所获取的点云进行聚类；计算每个聚类内点云的平均离地高度，对平均离地高度大于高度阈值的每一个聚类按步骤(2)-(4)进行处理，提取建筑物屋顶点；

(2)生成当前聚类C_p对应的DSM，获取C_p在DSM中有效像素的数目N_valid，计算C_p在XY平面的投影面积A_C；

(3)如果聚类C_p投影面积A_C大于等于面积阈值，计算DSM中拉普拉斯算子L_p小于平坦度阈值L_T的像素占有效像素总数N_valid的比值P_Laplacian，如果所述比值大于平顶比例阈值P_{Laplacian_T}，则将当前聚类C_p中的点全部分为建筑物屋顶点；

(4)如果聚类投影面积A_C小于面积阈值，按如下步骤进行判断：

(4.1)平移当前聚类C_p得到C_tmp，C_tmp的重心在坐标原点；初始化当前采样次数N_S＝0；初始化C_tmp中所有点的状态为未抽取，清空有效模型序列和最优模型序列；

(4.2)如果当前采样次数N_S达到最大采样次数N_{S_T}，跳转至步骤(4.5)；

从C_tmp中随机抽取状态为未抽取的点p₀，将p₀的状态修改为已抽取，寻找C_tmp内距离p₀最近的三个邻域点p₁、p₂、p₃；如果p₁、p₂、p₃距p₀的距离都足够小，且四点能够近似用平面拟合，计算拟合平面P_fit；否则当前采样次数N_S加一，重新执行步骤(4.2)；

(4.3)计算C_tmp中每个点到平面P_fit的距离；将所述距离小于平坦度阈值L_T的点组成内点集合C_inlier；如果集合C_inlier中点数小于聚类最小点数阈值，当前采样次数N_S加一，重新执行步骤(4.2)；否则对集合C_inlier再次进行三维欧几里得聚类；

(4.4)如果步骤(4.3)只得到一个聚类，将平面P_fit加入有效模型序列；

当前采样次数N_S加一，重新执行步骤(4.2)；

(4.5)对有效模型序列进行遍历，寻找内点数目最多的有效模型作为最优模型，加入最优模型序列；从C_tmp中移除最优模型所对应的内点集合中的点，如果剩余点数大于聚类最小点数阈值，则初始化当前采样次数N_S＝0；初始化C_tmp中所有点的状态为未抽取，清空有效模型序列和最优模型序列，跳转至步骤(4.2)；

(4.6)计算C_p内的每个点与最优模型序列中平面的最小距离d_min，将d_min小于平坦度阈值L_T的点组成集合C_{inlier_all}；如果C_{inlier_all}的点数占C_p中点数的比例大于阈值P_{inlier_all_T}，则将C_{inlier_all}中的点分为建筑物屋顶点。

2.根据权利要求1所述的提取建筑物屋顶点的方法，其特征在于，步骤(1)中所述三维欧几里得聚类的参数设置为：

搜索半径R的范围为(R_min，R_max)，其中

D是激光脉冲脚点的密度，R_max为相邻建筑物之间的最小距离；

聚类最小点数阈值N_min＝A_minD-Q

聚类最大点数阈值N_max＝A_maxD+Q

其中Q为误差调节常数，A_min为建筑物屋顶最小面积，A_max为建筑物屋顶最大面积。

3.根据权利要求1所述的提取建筑物屋顶点的方法，其特征在于，步骤(1)中计算聚类内点云的平均离地高度，包括：

获取聚类中每个点的地面高度：第n个点的坐标为(X_n，Y_n，H_n)，对应的地面高度H_Tn为XY平面上距离(X_n，Y_n)最近的地面点的高程；

计算聚类中每个点的离地高度：ΔH_n＝H_n-HT_n；

聚类内点云的平均离地高度：

N为聚类内点的总数。

4.根据权利要求1所述的提取建筑物屋顶点的方法，其特征在于，所述步骤(2)具体包括：

(2.1)计算C_p中点的最小X坐标X_min、最大X坐标X_max、最小Y坐标Y_min、最大Y坐标Y_max；将以(X_min，Y_min)、(X_max，Y_max)为对角顶点的矩形区域划分为N_row行、N_col列、格网大小为R的栅格，其中：

N_row＝int(Y_max-Y_min)/R+1

N_col＝int(X_max-X_min)/R+1

R为三维欧几里得聚类的搜索半径，int(·)为取整函数；

(2.2)将C_p中的点分配到栅格中的像素，得到当前聚类C_p对应的DSM；

第n个点映射到像素(i，j)中，

其中(X_n，Y_n，H_n)为第n个点的坐标；

(2.3)对栅格内的像素进行遍历；如果一个像素有点落入，则取落入点的高度指标作为DSM的值，该像素为有效像素，无点落入的像素的DSM值设为H_invalid；所述高度指标为：落入点的高程最小值，或落入点的高程最大值，或落入点的高程平均值；

(2.4)对DSM逐像素进行遍历，统计有效像素的数目N_valid；C_p在XY平面的投影面积A_C为：A_c＝N_validR²。

5.根据权利要求4所述的提取建筑物屋顶点的方法，其特征在于，所述步骤(3)DSM在像素(i，j)处的拉普拉斯算子L_p(i，j)为：

6.从机载激光雷达点云中提取建筑物屋顶点的系统，其特征在于，包括：

机载激光雷达点云聚类模块，用于获取测区的机载激光雷达点云；采用三维欧几里得聚类对所获取的点云进行聚类；计算每个聚类内点云的平均离地高度，获取平均离地高度大于高度阈值的聚类；

DSM生成模块，用于对机载激光雷达点云聚类模块获取的每一个聚类生成对应的DSM，并获取所述聚类在DSM中的有效像素数目，计算所述聚类在XY平面的投影面积A_C；并根据A_C将聚类分为大面积聚类和小面积聚类；

大面积聚类分类模块，用于判断大面积聚类中的点是否为建筑物屋顶上的点：判断方法为：对大面积聚类的DSM，统计拉普拉斯算子L_p小于平坦度阈值L_T的像素占DSM有效像素总数N_valid的比值P_Laplacian，如果所述比值大于平顶比例阈值P_{Laplacian_T}，则将当前聚类C_p中的点全部分为建筑物屋顶点；

小面积聚类分类模块，用于提取小面积聚类中为建筑物屋顶上的点，提取方法为：

(4.1)平移当前小面积聚类C_p得到C_tmp，使C_tmp的重心在坐标原点；初始化当前采样次数N_S＝0；初始化C_tmp中所有点的状态为未抽取，清空有效模型序列和最优模型序列；

(4.2)如果当前采样次数N_S达到最大采样次数N_{S_T}，跳转至步骤(4.5)；

从C_tmp中随机抽取状态为未抽取的点p₀，将p₀的状态修改为已抽取，寻找C_tmp内距离p₀最近的三个邻域点p₁、p₂、p₃；如果p₁、p₂、p₃距p₀的距离都足够小，且四点能够近似用平面拟合，计算拟合平面P_fit；否则当前采样次数N_S加一，重新执行步骤(4.2)；

(4.3)计算C_tmp中每个点到平面P_fit的距离；将所述距离小于平坦度阈值L_T的点组成内点集合C_inlier；如果集合C_inlier中点数小于聚类最小点数阈值，当前采样次数N_S加一，重新执行步骤(4.2)；否则对集合C_inlier再次进行三维欧几里得聚类；

(4.4)如果步骤(4.3)只得到一个聚类，将平面P_fit加入有效模型序列；

当前采样次数N_S加一，重新执行步骤(4.2)；

(4.5)对有效模型序列进行遍历，寻找内点数目最多的有效模型作为最优模型，加入最优模型序列；从C_tmp中移除最优模型所对应的内点集合中的点，如果剩余点数大于聚类最小点数阈值，则初始化当前采样次数N_S＝0；初始化C_tmp中所有点的状态为未抽取，清空有效模型序列和最优模型序列，跳转至步骤(4.2)；

(4.6)计算C_p内的每个点与最优模型序列中平面的最小距离d_min，将d_min小于平坦度阈值L_T的点组成集合C_{inlier_all}；如果C_{inlier_all}的点数占C_p中点数的比例大于阈值P_{inlier_all_T}，则C_{inlier_all}中的点为建筑物屋顶上的点。

7.根据权利要求6所述的从机载激光雷达点云中提取建筑物屋顶点的系统，其特征在于，所述机载激光雷达点云聚类模块和小面积聚类分类模块中的三维欧几里得聚类的参数设置为：

搜索半径R的范围为(R_min，R_max)，其中

D是激光脉冲脚点的密度，R_max为相邻建筑物之间的最小距离；

聚类最小点数阈值N_min＝A_minD-Q

聚类最大点数阈值N_max＝A_maxD+Q

其中Q为误差调节常数，A_min为建筑物屋顶最小面积，A_max为建筑物屋顶最大面积。

8.根据权利要求6所述的从机载激光雷达点云中提取建筑物屋顶点的系统，其特征在于，所述机载激光雷达点云聚类模块中计算聚类内点云的平均离地高度，包括：

获取聚类中每个点的地面高度：第n个点的坐标为(X_n，Y_n，H_n)，地面高度H_Tn为XY平面上距离(X_n，Y_n)最近的地面点的高程；

计算聚类中每个点的离地高度：ΔH_n＝H_n-H_Tn；

聚类内点云的平均离地高度：

N为聚类内点的总数。

9.根据权利要求6所述的从机载激光雷达点云中提取建筑物屋顶点的系统，其特征在于，所述DSM生成模块生成聚类C_p对应的DSM、获取C_p在DSM中的有效像素数目并计算聚类C_p在XY平面的投影面积，以及将聚类分为大面积聚类和小面积聚类的步骤包括：

(2.1)计算C_p中点的最小X坐标X_min、最大X坐标X_max、最小Y坐标Y_min、最大Y坐标Y_max；将以(X_min，Y_min)、(X_max，Y_max)为对角顶点的矩形区域划分为N_row行、N_col列、格网大小为R的栅格，其中：

N_row＝int(Y_max-Y_min)/R+1

N_col＝int(X_max-X_min)/R+1

R为三维欧几里得聚类的搜索半径，int(·)为取整函数；

(2.2)将C_p中的点分配到栅格中的像素，得到当前聚类C_p对应的DSM；

第n个点映射到像素(i，j)中，

其中(X_n，Y_n，H_n)为第n个点的坐标；

(2.3)对栅格内的像素进行遍历；如果一个像素有点落入，则标识该像素为有效像素，取落入点的高度指标作为DSM的值，无点落入的像素的DSM值设为H_invalid；所述高度指标为：落入点的高程最小值，或落入点的高程最大值，或落入点的高程平均值；

(2.4)对DSM逐像素进行遍历，统计有效像素的数目N_valid；C_p在XY平面的投影面积A_C为：A_c＝N_validR²；

如果聚类投影面积A_C大于等于面积阈值，则聚类C_p为大面积聚类；

如果聚类投影面积A_C小于面积阈值，则聚类C_p为小面积聚类。

10.根据权利要求6所述的从机载激光雷达点云中提取建筑物屋顶点的系统，其特征在于，所述大面积聚类分类模块中计算DSM在像素(i，j)处的拉普拉斯算子L_p(i，j)为：

Images