一种多尺度自适应机载LiDAR点云建筑物单体化分割方法与流程

本发明涉及计算机三维重建技术领域，尤其涉及一种多尺度自适应机载lidar点云建筑物单体化分割方法，主要应用于建筑物的建模、城市规划、市政管理、数字城市建设等多个领域。

背景技术：

机载lidar可以快速获取地物三维表面信息，是建筑物三维建模的重要数据来源。我国所建设数字城市系统中的三维模型是缺乏语义信息的模型或者是由人工手动单体分割建立的模型，并不能满足真正数字城市建设与应用的需求，数字城市的建设亟需高精度的单体化三维模型作为基础。建筑物点云的单体分割可以并行模型重建、高效管理与分析各个单体建筑物，生成可查询、可分析的建筑物三维模型。然而，由于点云数据的立体性、层次性和海量性等并且场景中有大量老城区的密集建筑、物现代市区的裙楼以及由于建筑材料等导致建筑物点云密度不均匀等，其结构复杂和空间分布不规则，如何快速、高效、自动的实现建筑物的单体分割仍然是一个具有挑战的课题，也是一个亟需解决的难题。

目前对建筑物单体分割主要用八邻域聚类算法、欧式距离聚类算法、密度聚类算法、区域增长(或移动窗口法)等，但对于建筑物单体分割的研究大部分都是在建筑物建模过程中对建筑物点云在二维空间内聚类，对建筑物结构复杂和裙楼建筑物和点云密度不均匀建筑物缺少深入研究，无法实现单体分割。与现有建筑物点云单体化相比，本发明的技术优势在于抗噪声强、鲁棒性好，通过多尺度与自适应性实现了对密集城区的建筑物、裙楼结构的建筑物和点云密度不均匀建筑物点云的单体分割。

技术实现要素：

本发明要解决的技术问题在于规避现有方法单一尺度非自适应的问题，提供一种多尺度自适应机载lidar点云建筑物单体化分割方法，将现有的聚类方法与建筑物单体化评估方法相结合，充分考虑建筑物结构复杂和空间分布不规则的特点，以解决老城区的密集建筑物、现代市区的裙楼以及由于建筑材料等导致建筑物点云密度不均匀的数据建筑物单体分割问题，对多种数据具有自适应性，将过分割的单体进行合并，实现了多尺度自适应的建筑物单体分割。

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

本发明提供一种多尺度自适应机载lidar点云建筑物单体化分割方法，该方法包括以下步骤：

步骤1、获取建筑物的点云数据，并根据点云数据计算建筑物的点云间距；

步骤2、通过二维多尺度密度聚类算法对建筑物点云数据进行单体分割，将其分割为单体结构和非单体结构；

步骤3、去除非单体结构建筑物立面点云，将其在z轴方向上进行缩放变换，对变换后的点云运用三维多尺度密度聚类算法进行单体分割，将其分割为屋顶单体结构和非屋顶单体结构；

步骤4、对非屋顶单体结构的建筑物点云在未缩放变换的尺度上再次运用三维多尺度密度聚类算法，将其分割为屋顶单体结构；

步骤5、将屋顶单体结构投影到xoy平面内，进行区域增长，获得三维单体结构和三维单体结构之外的其它点云；

步骤6、将三维单体结构之外的其它点云运用密度聚类算法进行聚类，对屋顶细节结构进行识别并合并到对应的屋顶单体结构中，对于最小外接矩形长和宽大于最小建筑物阈值的非屋顶细节结构簇识别为小建筑物单体结构，否则，该簇为建筑物附属物点云；

步骤7、将建筑物附属物点云根据最近邻域内所属单体建筑类别，合并到对应的单体结构中，完成建筑物单体分割。

进一步地，本发明的步骤2中所述的通过二维多尺度密度聚类算法对建筑物点云进行单体分割的具体方法为：

步骤2.1、运用二维多尺度密度聚类算法对建筑物点云进行聚类分割；

步骤2.2、根据最小外接矩形旋转法求外接矩形的面积，利用面积的单调性，使点云向面积减小的方向旋转，求面积的最小值来改进最小外接矩形算法；

步骤2.3、利用改进的最小外接矩形算法，根据最小外接矩形长和宽的阈值判断该簇是噪点还是待评估的结构，同时根据该簇点云数与该簇最小外接矩形面积的比值判断该簇是噪点还是待评估的结构；

步骤2.4、对待评估的结构利用提出的建筑物单体分割评估方法进行评估，判断该簇点是单体结构还是非单体结构；

步骤2.5、对非单体结构，通过减小搜索半径阈值再次进行聚类，直至搜索半径减小到最小阈值，完成二维多尺度密度聚类算法的单体分割。

进一步地，本发明的步骤2.2中所述的改进的最小外接矩形算法的具体方法为：

步骤2.2.1、对该簇点云旋转α＝0度，求该簇点云的最大值与最小值，通过最大值与最小值计算该簇点的面积为area1；

步骤2.2.2、如果α＜90，对该点进行旋转α＝α+θ，求旋转后该簇点云的最大值与最小值，通过最大值与最小值计算该簇点的面积为area2，比较area1与area2的大小，如果area1大于area2，则将area2的值赋值给area1，则对点云簇继续旋转α＝α+θ，直至area1小于area2，返回α＝α-θ角度和旋转该角度对应的顶点坐标；

步骤2.2.3、如果α＜90，对该点进行旋转α＝α+θ，求旋转后该簇点云的最大值与最小值，通过最大值与最小值计算该簇点的面积为area2，比较area1与area2的大小，如果area1小于area2，将area2的值赋值给area1，则对该点进行旋转α＝α-θ，如果α＞-90，求旋转后该簇点云的最大值与最小值，通过最大值与最小值计算该簇点的面积为area2，如果area1大于area2，则对点云簇继续旋转α＝α-θ，直至area1小于area2，返回α＝α+θ角度和旋转该角度对应的顶点坐标；

步骤2.2.4、根据最小外接矩形顶点坐标，求最小外接矩形的长、宽和面积。

进一步地，本发明的步骤2.4中所述的对待评估的结构利用提出的建筑物单体分割评估方法进行评估，判断该簇点是单体结构还是非单体结构，其具体方法为：

步骤2.4.1、对该簇点云求最大值最小值，根据点云的最值，计算格网的行数和列数，建立二维格网，将点云填充到对应的格网内；

步骤2.4.2、对格网内点数不为零的格网标记为true，否则标记为fault，统计有点格网的个数number2dt；

步骤2.4.3、利用步骤2.4.1的最值，计算格网的行数nx、列数ny和层数nz，建立三维格网；

步骤2.4.4、对三维格网从零行零列的最上层开始往下查找格网内点数不为零的格网，将其标记为true，对应行列号下面所有层均标记为true，一直到底层，并对每个访问过得格网标记已访问；

步骤2.4.5、按照顺序，从最上层开始往下查找格网内点数不为零的格网，按照步骤2.4.4对所有的三维格网进行查找，统计三维格网标记为true的个数number3d；

步骤2.4.6、根据建筑物是柱体结构，利用二维格网内有点云的格网数与三维格网的层数，可以根据公式计算理论上标记三维格网标记为true的格网个数number3dt；其公式为：

number3dt＝number2dt*nz

式中，number3dt表示三维格网内理论上标记为true的格网个数；number2dt表示二维格网标记为true的格网个数；nz表示三维格网的层数；

步骤2.4.7、根据三维格网实际标记为true的格网个数与理论上标记为true的格网个数比值，判断该簇点云为单体结构还是非单体结构：

ratio＝number3d/number3dt

式中，ratio是三维格网内实际标记为true的格网个数与理论上标记为true的格网个数比值，比值越接近1越有可能是单体结构；number3d表示三维格网内实际标记为true的格网个数。

进一步地，本发明的步骤3中所述的去除非单体结构建筑物立面点云，将其在z轴方向上进行缩放变换，对变换后的点云运用三维多尺度密度聚类算法进行单体分割，将其分割为屋顶单体结构和非屋顶单体结构，其具体方法为：

步骤3.1、对非单体结构点云创建kd-tree，对非单体结构表面法向量进行估计，通过点云的法向量删除建筑物立面点云；

步骤3.2、将删除立面的非单体结构点云将z值进行缩放变换，x、y值不变；

步骤3.3、对缩放变换的点云，运用密度聚类算法进行聚类；

步骤3.4、利用改进的最小外接矩形算法，根据最小外接矩形的长和宽的阈值判断该簇是噪点还是待评估的结构，同时根据该簇点云数与该簇最小外接矩形面积的比值判断该簇是噪点还是待评估的结构；

步骤3.5、对待评估的结构利用提出的建筑物单体分割评估方法进行评估，判断该簇点是屋顶单体结构还是非屋顶单体结构；

步骤3.6、通过减小搜索半径阈值再次进行聚类，直至搜索半径减小到最小阈值，完成三维多尺度密度聚类算法的单体分割。

进一步地，本发明的所述的将屋顶单体结构投影到xoy平面内，进行区域增长，获得三维单体结构和三维单体结构之外的其它点云，其具体方法为：

步骤5.1、求步骤2中的非单体结构点云的最大值与最小值，根据格网的大小，计算格网的行数、列数，将点云填充到二维格网内；

步骤5.2、建立步骤5.1同样的格网，将某个单体结构填充到建立的格网内，根据格网内点数为非零的格网的id找到步骤5.1中对应的格网id，将该格网内的点添加到区域增长的类中，完成该单体的区域增长；

步骤5.3、对步骤3、步骤4所得的屋顶单体结构，重复步骤5.2，实现屋顶单体结构的区域增长，获得三维单体结构。

进一步地，本发明的所述的将三维单体结构之外的其它点云运用密度聚类算法进行聚类，对屋顶细节结构进行识别并合并到对应的屋顶单体结构中，对于最小外接矩形长和宽大于最小建筑物阈值的非屋顶细节结构簇识别为小建筑物单体结构，否则，该簇为建筑物附属物点云；其具体方法为：

步骤6.1、对三维单体结构之外的其它点云运用密度聚类算法聚类，获得各个簇结构；

步骤6.2、对每一个单体结构求最值，对步骤6.1聚类的簇求最值，通过下面公式判断该结构是否为屋顶细节结构；

maxsx＞maxcx

maxsy＞maxcy

maxsz-maxcz＞buffer_zmin

minsx＞mincx

minsy＞mincy

maxsz-maxcz＜buffer_zmax

式中，单体结构求最大值(maxsx，maxsy，maxsz)，最小值(minsx，minsy，minsz)，步骤6.1聚类簇的最大值(maxcx，maxcy，maxcz)，最小值(mincx，mincy，mincz)，buffer_zmin是屋顶单体细节结构的最大值与屋顶单体结构最大值差值的最小值；buffer_zmax是屋顶单体细节结构的最大值与屋顶单体结构最大值差值的最大值；

步骤6.3、满足步骤6.2的单体结构为屋顶细节结构，将该屋顶细节结构点云合并到对应的单体结构中，实现单体结构的合并；

步骤6.4、对于不满足步骤6.2的并且最小外接矩形的长和宽大于最小单体建筑物阈值的为小建筑物单体结构，否则，该簇为建筑物附属物点云。

本发明产生的有益效果是：本发明的多尺度自适应机载lidar点云建筑物单体化分割方法运用多尺度自适应的算法，相比于现有的技术有以下明显优势：

(1)本发明具有抗噪声强、鲁棒性好，通过多尺度与自适应性实现了对密集城区的建筑物、裙楼结构的建筑物和点云密度不均匀建筑物点云进行单体分割，取得了较好的效果。

(2)本发明公开的多尺度方法，可以解决建筑物点云密度不均匀的问题，尤其由于建筑物屋顶材料反光，部分点云密度较小，现在的算法对密度不均匀的建筑物点云，设置参数较大时存在多个建筑物无法分割成单体，参数较小时密度较小的区域分割为多个单体；

(3)本发明公开的自适应方法，可以通过二维多尺度密度聚类实现密度不均匀，各个单体在二维空间内非连通的可以实现较好的分割；去除立面点云，点云缩放变换后的三维多尺度密度聚类算法可以实现现代都市裙楼结构的建筑物单体分割；对于建筑物高度一致距离相近的再次用密度聚类算法，实现老城区密集区域建筑物单体分割；

(4)本发明公开的建筑物单体分割评估方法实现了算法的多尺度与自适应性；

(5)本发明公开的改进的最小外接矩形算法，根据面积的单调性，在最小外接矩形面积递减再递增的拐点，找到最小外接矩形面积的最小值，可以较快的通过旋转法求得该簇点云的最小外接矩形，与现有的方法相比，按照一度旋转一次，需要旋转90次，结合面积的单调性，最少旋转3次，最多旋转47次，想到与提高1.91～30倍；本发明对密集城区的建筑物、裙楼结构的建筑物和点云密度不均匀建筑物点云的单体分割可以具有较好的自适应性。

附图说明

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

图1为本发明的流程图；

图2为二维格网的标记示意图，(a)点云投影到xoy平面二维格网内，(b)搜索含有点云的格网并标记网进行标记，(c)对整列二维格网进行标记，(d)整个二维格网进行标记；

图3为三维格网的标记示意图,(a)点云填充到三维格网内,(b)从最上层开始搜索含有点云的格网,(c)将有点格网该行列号层下面所有的格网均标记为true,(d)根据(b)、(c)的方法标记所有格网；

图4为建筑物单体分割示例1，德国法伊英根(vaihingen)地区机载lidar点云数据，不同颜色标记为不同建筑物单体：(a)分类后的建筑物点云数据，(b)八邻域聚类算法单体分割结果，(c)密度聚类算法单体分割结果，(d)本文方法单体分割结果。

图5为建筑物单体分割示例2，宁波地区机载lidar点云数据，不同颜色标记为不同建筑物单体：(a)分类后的建筑物点云数据，(b)八邻域聚类算法单体分割结果，(c)密度聚类算法单体分割结果，(d)本文方法单体分割结果。

具体实施方式

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

为了解决密老城区的密集建筑物、现代都市的裙楼以及由于建筑材料等导致建筑物点云密度不均匀建筑物点云单体化问题，以下结合附图和实施例对本发明技术方案进行具体说明。

为了测试该技术方案的正确性，选取isprs上的点云数据和在宁波采集的机载lidar点云数据进行实验。该数据含有较为密集的建筑物点云数据、现代都市的裙楼以及由于建筑材料等导致建筑物点云密度不均匀建筑物点云如图4(a)、图5(a)所示，对本发明的过程进行详细讲解。一种多尺度自适应机载lidar点云建筑物单体化分割方法，其流程图如图1所示，包括以下步骤：

步骤1、获取建筑物的点云数据，并根据点云数据计算建筑物的点云间距；

步骤1.1、读取建筑物点云，遍历点云，求点云x、y的最值；

中部1.2、对于给定的格网边长l，统计落入格网内点云的数量，记录格网点数大于零并且最多格网数量对应的点数；

步骤1.3、根据最多格网数量对应的点数，计算每平米点云数量，根据如下公式计算建筑物点云的点间距；

式中，iinterval为点云的点间距，nnumber为每平米点云的个数。

步骤2、通过二维多尺度密度聚类算法对建筑物点云数据进行单体分割，将其分割为单体结构和非单体结构；

步骤2.1、运用二维多尺度密度聚类算法对建筑物点云进行聚类分割；

步骤2.2、根据最小外接矩形旋转法求外接矩形的面积，利用面积的单调性，使点云向面积减小的方向旋转，求面积的最小值来改进最小外接矩形算法；

步骤2.2.1、对该簇点云旋转α＝0度，求该簇点云的最大值与最小值，通过最大值与最小值计算该簇点的面积为area1；

步骤2.2.2、如果α＜90，对该点进行旋转α＝α+θ，求旋转后该簇点云的最大值与最小值，通过最大值与最小值计算该簇点的面积为area2，比较area1与area2的大小，如果area1大于area2，则将area2的值赋值给area1，则对点云簇继续旋转α＝α+θ，直至area1小于area2，返回α＝α-θ角度和旋转该角度对应的顶点坐标；

步骤2.2.3、如果α＜90，对该点进行旋转α＝α+θ，求旋转后该簇点云的最大值与最小值，通过最大值与最小值计算该簇点的面积为area2，比较area1与area2的大小，如果area1小于area2，将area2的值赋值给area1，则对该点进行旋转α＝α-θ，如果α＞-90，求旋转后该簇点云的最大值与最小值，通过最大值与最小值计算该簇点的面积为area2，如果area1大于area2，则对点云簇继续旋转α＝α-θ，直至area1小于area2，返回α＝α+θ角度和旋转该角度对应的顶点坐标；

步骤2.2.4、根据最小外接矩形顶点坐标，求最小外接矩形的长、宽和面积。

步骤2.3、利用改进的最小外接矩形算法，根据最小外接矩形长和宽的阈值判断该簇是噪点还是待评估的结构，同时根据该簇点云数与该簇最小外接矩形面积的比值判断该簇是噪点还是待评估的结构；

步骤2.4、对待评估的结构利用提出的建筑物单体分割评估方法进行评估，判断该簇点是单体结构还是非单体结构；

步骤2.4.1、对该簇点云求最大值最小值，根据点云的最值，计算格网的行数和列数，建立二维格网，将点云填充到对应的格网内；

步骤2.4.2、对格网内点数不为零的格网标记为true，否则标记为fault，统计有点格网的个数number2dt；

步骤2.4.3、利用步骤2.4.1的最值，计算格网的行数(nx)、列数(ny)和层数(nz)，建立三维格网；

步骤2.4.4、对三维格网，从零行零列的最上层开始往下查找格网内点数不为零的格网，将其标记为true，对应行列号下面所有层均标记为true，一直到底层，并对每个访问过得格网标记已访问；

步骤2.4.5、按照顺序，从最上层开始往下查找格网内点数不为零的格网，按照步骤2.4.4对所有的三维格网进行查找，统计三维格网标记为true的个数number3d；

步骤2.4.6、根据建筑物是柱体结构，利用二维格网内有点云的格网数与三维格网的层数，可以根据公式计算理论上标记三维格网标记为true的格网个数number3dt；其公式为：

number3dt＝number2dt*nz

式中，number3dt表示三维格网内理论上标记为true的格网个数；number2dt表示二维格网标记为true的格网个数；nz表示三维格网的层数；

步骤2.4.7、根据三维格网实际标记为true的格网个数与理论上标记为true的格网个数比值，判断该簇点云为单体结构还是非单体结构：

ratio＝number3d/number3dt

式中，ratio是三维格网内实际标记为true的格网个数与理论上标记为true的格网个数比值，比值越接近1越有可能是单体结构；number3d表示三维格网内实际标记为true的格网个数。

步骤2.5、对非单体结构，通过减小搜索半径阈值再次进行聚类，直至搜索半径减小到最小阈值，完成二维多尺度密度聚类算法的单体分割。

步骤3、去除非单体结构建筑物立面点云，将其在z轴方向上进行缩放变换，对变换后的点云运用三维多尺度密度聚类算法进行单体分割，将其分割为屋顶单体结构和非屋顶单体结构；

步骤3.1、对非单体结构点云创建kd-tree，对非单体结构表面法向量进行估计，通过点云的法向量删除建筑物立面点云；

步骤3.2、将删除立面的非单体结构点云将z值进行缩放变换，x、y值不变；

步骤3.3、对缩放变换的点云，运用密度聚类算法进行聚类；

步骤3.4、利用改进的最小外接矩形算法，根据最小外接矩形的长和宽的阈值判断该簇是噪点还是待评估的结构，同时根据该簇点云数与该簇最小外接矩形面积的比值判断该簇是噪点还是待评估的结构；

步骤3.5、对待评估的结构利用提出的建筑物单体分割评估方法进行评估，判断该簇点是屋顶单体结构还是非屋顶单体结构；

步骤3.6、通过减小搜索半径阈值再次进行聚类，直至搜索半径减小到最小阈值，完成三维多尺度密度聚类算法的单体分割。

步骤4、对非屋顶单体结构的建筑物点云在未缩放变换的尺度上再次运用三维多尺度密度聚类算法，将其分割为屋顶单体结构；

步骤5、将屋顶单体结构投影到xoy平面内，进行区域增长，获得三维单体结构和三维单体结构之外的其它点云；

步骤5.1、求步骤2中的非单体结构点云的最大值与最小值，根据格网的大小，计算格网的行数、列数，将点云填充到二维格网内；

步骤5.2、建立步骤5.1同样的格网，将某个单体结构填充到建立的格网内，根据格网内点数为非零的格网的id找到步骤5.1中对应的格网id，将该格网内的点添加到区域增长的类中，完成该单体的区域增长；

步骤5.3、对步骤3、步骤4所得的屋顶单体结构，重复步骤5.2，实现屋顶单体结构的区域增长，获得三维单体结构和三维单体结构之外的其它点云。

步骤6、将三维单体结构之外的其它点云运用密度聚类算法进行聚类，对屋顶细节结构进行识别并合并到对应的屋顶单体结构中，对于最小外接矩形长和宽大于最小建筑物阈值的非屋顶细节结构簇识别为小建筑物单体结构，否则，该簇为建筑物附属物点云；

步骤6.1、对三维单体结构之外的其它点云运用密度聚类算法聚类，获得各个簇结构；

步骤6.2、对每一个单体结构求最值，对步骤6.1聚类的簇求最值，通过下面公式判断该结构是否为屋顶细节结构；

maxsx＞maxcx

maxsy＞maxcy

maxsz-maxcz＞buffer_zmin

minsx＞mincx

minsy＞mincy

maxsz-maxcz＜buffer_zmax

式中，单体结构求最大值(maxsx，maxsy，maxsz)，最小值(minsx，minsy，minsz)，步骤6.1聚类簇的最大值(maxcx，maxcy，maxcz)，最小值(mincx，mincy，mincz)，buffer_zmin是屋顶单体细节结构的最大值与单体结构(即屋顶单体结构)最大值差值的最小值；buffer_zmax是屋顶单体细节结构的最大值与单体结构(即屋顶单体结构)最大值差值的最大值；

步骤6.3、满足步骤6.2的单体结构为屋顶细节结构，将该屋顶细节结构点云合并到对应的单体结构中，实现单体结构的合并；

步骤6.4、对于不满足步骤6.2的并且最小外接矩形的长和宽大于最小单体建筑物阈值的为小建筑物单体结构，否则，该簇为建筑物附属物点云。

步骤7、将建筑物附属物点云根据最近邻域内所属单体建筑类别，合并到对应的单体结构中，完成建筑物单体分割。

步骤7.1、将各个单体结构标记上label，合并到一个数组中，对该数据建立kd-tree；

步骤7.2、以建筑物附属物点为搜索中心，在单体结构中搜索最邻近的点，将最邻近点的单体结构label赋值给该点，并将该点添加到该簇中；

步骤7.3、重复步骤7.2，将步骤处理后剩余的点合并到对应的单体上，实现最终的建筑物单体分割。

图2是将点云填充到二维格网内，对格网内点数大于0的格网进行标记并统计个数。

图3是将点云填充到二维格网内，从三维格网的最顶层开始查找，对于格网点数大于0的格网，对应行列号下面所有层格网均标记为true。

以下以德国vaihingen区域和中国宁波的样例数据为实施例说明。本算法仅应用分类后的建筑点云数据(集)，其中德国vaihingen区域为样例1数据如图4所示。子图(a)为所提取的建筑区域根据高程进行显示的激光点云；(b)八邻域聚类算法单体分割结果；(c)密度聚类算法单体分割结果；(d)本文方法单体分割结果。中国宁波为样例2数据如图5所示。子图(a)为所提取的建筑区域根据高程进行显示的激光点云；(b)八邻域聚类算法单体分割结果，(c)密度聚类算法单体分割结果；(d)本文方法单体分割结果。

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