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

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

背景技术：

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

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

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

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

但是，基于图像和基于点云的方法都有明显的问题。基于图像的方法往往会损失原始点云中的信息，尤其是点云的垂直分布所反映的地物垂直结构信息，对于提取树木遮挡的建筑物的效果并不好。另外人工修剪过的平顶植被容易被误分为屋顶点。基于区域生长的算法需要给出法向量夹角阈值和曲率阈值，而这两个阈值的取值受到数据质量和屋顶平整度的影响，难以简单确定。更重要的是，该方法提取到的屋顶点并不全，经常存在大量漏检的屋顶点。hough变换提取平面需要三个参数，而三参数的hough变换非常复杂，需要在三维hough空间进行分块，用作累加器的分区，以进行投票。空间分块的粒度直接影响后续平面检测的效果，而且在检测结果中经常存在错误。在点云量大时计算复杂度很高，需要大量存储空间。基于ransac方法的基本思路是从点云中随机抽取3个点构建平面模型，测试其它点是否可以很好地用该平面进行拟合。上述过程重复执行很多次。假设点云共有n个点，则随机取3个点共有个不重复的组合。随着n的增加，组合的数目迅速增加。如果n达到上万，会是个天文数字。而随机采样的次数是有限的，这就导致找到稳健模型的概率明显降低。即使有限度地增加随机采样次数，也不会明显解决问题，反而会造成计算时间显著延长。从这个角度讲，ransac方法并不适合点数较多的情况。此外，从理论上讲，即使对于随机分布的点，ransac都可能得到一些模型。例如，森林内的点云，也可以用ransac探测出一些平面出来，而且平面上的点可能足够多。在这种情况下，传统ransac难以确保平面来自建筑物而不是非建筑物。

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

技术实现要素：

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

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

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

(2)生成当前聚类cp对应的dsm，获取cp在dsm中有效像素的数目nvalid，计算cp在xy平面的投影面积ac；

(3)如果cp的投影面积ac大于等于面积阈值，计算dsm中拉普拉斯算子lp小于平坦度阈值lt的像素占dsm有效像素总数nvalid的比值plaplacian，如果所述比值大于平顶比例阈值plaplacian_t，则当前聚类cp中的点全部为建筑物屋顶上的点；

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

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

(4.2)如果当前采样次数ns达到最大采样次数ns_t，跳转至步骤4.5；

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

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

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

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

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

(4.6)计算cp内的每个点与最优模型序列中平面的最小距离dmin，将dmin小于平坦度阈值lt的点组成集合cinlier_all；如果cinlier_all的点数占cp中点数的比例大于阈值pinlier_all_t，则将cinlier_all中的点分为建筑物屋顶点。

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

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

大面积聚类分类模块，用于判断大面积聚类中的点是否为建筑物屋顶上的点：判断方法为：对大面积聚类计算dsm，统计拉普拉斯算子lp小于平坦度阈值lt的像素占有效像素总数nvalid的比例plaplacian，如果所述比值大于平顶比例阈值plaplacian_t，则将当前聚类cp中的点全部分为建筑物屋顶点；

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

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

(4.2)如果当前采样次数ns达到最大采样次数ns_t，跳转至步骤4.5；

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

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

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

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

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

(4.6)计算cp内的每个点与最优模型序列中平面的最小距离dmin，将dmin小于平坦度阈值lt的点组成集合cinlier_all；如果cinlier_all的点数占cp中点数的比例大于阈值pinlier_all_t，则cinlier_all中的点为建筑物屋顶上的点。

有益效果：与现有技术相比，本发明公开的从机载激光雷达点云中提取建筑物屋顶点的方法具有以下优点：1、本发明采用三维欧几里得聚类，可以在大部分情况下有效地将建筑物屋顶与临近的树木分到不同的聚类，在提取到平面模型内点之后对平面模型内点进行欧几里得聚类可以将绝大部分非屋顶的树木点排除，提高了屋顶点提取的准确率；；2、本发明采用大小屋顶分而治之的策略，避免了对大屋顶也采用ransac进行平面检测带来的高计算量和潜在的漏检可能性，同时维持了高准确度；3、对小面积聚类，采用改进的ransac算法；设某聚类中点数为n，则潜在的平面模型数目从降到n。如果n小于指定运行次数，则随机取点实际上只需要运行n次；4、本发明理论简单易行，计算效率高。对于300万点，一般可以在十秒内完成。可以多次运行，以更全面地提取屋顶点。

附图说明

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

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

图3是构建dsm的示意图；

图4是判断某个随机选取的点p0是否和周围3个邻域点p1、p2、p3可以近似地用平面拟合的示意图；

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

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

具体实施方式

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

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

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

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

搜索半径r的范围为(rmin,rmax)，r最小值rmin与激光脉冲脚点密度负相关，本实施例中，d是激光脉冲脚点的密度，可以由航飞设计方案获得。如果r小于rmin，可能会导致聚类点数过少，一个屋顶被分为多个聚类。r最大值rmax可以取相邻建筑物之间的最小距离，如1.5m，2.0m等。如果r超过rmax，可能会将多个屋顶、屋顶与树木等聚到一类，不能实现准确分割。

nmin、nmax与建筑物屋顶的最小、最大面积正相关。设建筑物屋顶最小面积为amin，最大面积为amax，则nmin、nmax的计算式为：

nmin＝amind-q

nmax＝amaxd+q

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

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

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

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

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

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

n为聚类内点的总数。

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

本实施例中，δhmean_t取值为2.5m。

步骤2、生成当前聚类cp对应的dsm，获取cp在dsm中有效像素的数目nvalid，计算cp在xy平面的面积ac；具体步骤为：

(2.1)计算cp中点的最小x坐标xmin、最大x坐标xmax、最小y坐标ymin、最大y坐标ymax；将以(xmin,ymin)、(xmax,ymax)为对角顶点的矩形区域划分为nrow行、ncol列、格网大小为r的栅格，其中：

nrow＝int(ymax-ymin)/r+1

ncol＝int(xmax-xmin)/r+1

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

(2.2)将cp中的点根据其x、y坐标分配到栅格中的像素，得到当前聚类cp对应的dsm；

第n个点映射到像素(i,j)中，其中(xn,yn,hn)为第n个点的坐标；

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

(2.4)对dsm逐像素进行遍历，统计有效像素的数目nvalid；cp在xy平面的投影面积ac为：ac＝nvalidr²。

如果ac大于等于指定的面积阈值ac_t，则认为认为cp为大面积聚类；否则cp为小面积聚类。

步骤3、对于大面积聚类，计算dsm中拉普拉斯算子lp小于平坦度阈值lt的像素占有效像素总数nvalid的比值plaplacian，如果所述比值大于平顶比例阈值plaplacian_t，则当前聚类cp中的点全部为建筑物屋顶上的点；

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

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

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

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

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

计算聚类cp的重心op(xop,yop,hop)：

ctmp中第n个点的坐标为：(xn-xop,yn-yop,hn-hop)，则ctmp的重心在坐标原点；

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

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

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

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

p1、p2、p3和p0构成三个三角形，即δp0p1p2、δp0p2p3、δp0p3p1。计算这些三角形的法向量，δp0p1p2的法向量n012为：n012＝(p1-p0)×(p2-p0)；对n012进行模标准化。δp0p2p3的法向量n023、δp0p3p1的法向量n031采用类似的方法进行计算和标准化。

计算n012和n023的夹角θ0：θ0＝cos^-1(n012·n023)；

n012和n031的夹角θ1、n023和n031的夹角θ2采用类似的方法进行计算。cos^-1为反余弦函数。

如果θ0、θ1、θ2均接近0°或180°，认为p0、p1、p2和p3可以用某平面pfit进行拟合；否则认为四个点不能用平面进行拟合；

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

a、b、a分别为平面法向量在三个坐标轴上的分量，d＝-(ax0+by0+ch0)，(x0,y0,h0)为p0的坐标。

(4.3)计算ctmp中每个点到平面pfit的距离，其中第n个点到pfit的距离dfitn为：

(xtmpn,ytmpn,htmpn)为ctmp中第n个点的坐标；

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

搜索半径为r，聚类最小点数阈值为1，聚类最大点数阈值为集合cinlier中的总点数ninlier；

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

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

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

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

(4.6)计算cp内的每个点与最优模型序列中平面的最小距离dmin，将dmin小于平坦度阈值lt的点组成集合cinlier_all；如果cinlier_all的点数占cp中点数的比例大于阈值pinlier_all_t，则认为cinlier_all为建筑物屋顶，将其中的点分为建筑物屋顶点。本实施例中，pinlier_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方法快速从小面积聚类点云中提取平面以及对模型内点使用聚类数目来判断是否为有效的模型。具体实现该技术方案的方法和途径很多，以上所述仅是本发明的优选实施方式。应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰。这些改进和润饰也应视为本发明的保护范围。本实施例中未明确的各组成部分均可用现有技术加以实现。