一种从机载激光雷达点云中快速提取平面片的方法与流程

本发明涉及摄影测量与遥感，尤其是一种从机载激光雷达点云中快速提取平面片的方法。

背景技术：

相较于传统的人工测绘或立体影像(像对)等，利用机载激光雷达(airbornelidar)所获取的高精度三维点云(pointcloud)，可以更高效、更精确、更智能化地重建数字城市三维模型(3dmodel)。建筑物(构筑物)是城市(城镇)中最常见、最主要的地物特征(以下统称为建筑物)，其最主要的几何形态是多面体(由多个平面形成的封闭体)；此外，对非平面表面，也可以通过多个小平面片(planarpatch)进行近似；也即，平面(平面片)是多面体建筑模型的基本单元之一，是城市三维激光点云数据中的基础几何形体。从三维点云中快速、精确、准确地提取出平面片，是数字城市三维模型重建的首要步骤。

平面区域分割过程是将属于同一平面片的lidar点聚类的过程；聚类则通常基于空间邻接性(proximity)及属性相似性(similarity)等特征。已有的从机载lidar数据中提取平面片的方法可分为两大类，即，空间域分割方法，或参数域分割方法。在空间域进行平面片聚类或分割，常用的方法是从种子点出发进行区域增长，但该方法严重依赖于种子点的选择，也即鲁棒性不足。为了改进依赖种子点的问题，也有改进利用ransac(randomsampleconsensus)，在物方空间中随机采样，选择种子点并估计平面参数，但ransac类的方法计算过程相对复杂。

在参数域进行平面聚类分割，无需选择种子点。参数域方法需要合理选择属性参数以及对这些参数进行估值。对机载lidar点云而言，其三维点云中的点仅是空间表面的某种采样，在恢复出由这些采样点所代表的空间表面的属性时，需要基于某种局部邻域对这些参数值进行估计。最常用的平面参数(属性)是法向量，lidar点的法向量通常在其某种(由不同邻接关系所定义)邻域内计算。计算出采样点的属性值后，hough变换常被用来检测三维空间中的平面，并估计相应的平面参数。参数域平面分割方法的主要局限是，无法区分共面但不邻接的平面片，另外，在数据量大的情况下，其计算效率低。

参数域及空间域结合的平面片提取方法，在一定程度上综合了上述两种方法的优点，但同时也面临种子点选择的性能，及提取效率较低的问题。

技术实现要素：

本发明要解决上述现有技术的缺点，即，依赖于种子点的选择(鲁棒性不足)、计算效率低、以及难以处理大范围、海量点云数据等问题(城市范围的点云数据通常为数十亿量级或更大)，提供一种鲁棒性好、效率高的的从机载激光雷达点云中快速提取平面片的方法，以此基础生成建筑物三维模型，可应用于智慧城市或数字城市的数据基础设施建设(三维建模、数据实时更新等)，服务于城市科学化、精细化管理与运营。

本发明解决其技术问题采用的技术方案：这种从机载激光雷达点云中快速提取平面片的方法，包括以下步骤：1)对机载激光雷达点云数据集，利用delaunay剖分三角形的法向量，作为平面片的属性参数；2)生成三角形单位法向量直方图，并检测直方图的峰值进而确定平面片种子点；3)选择任意一个与直方图局部峰值具有相似属性的三角形，开始平面片的区域增长过程；4)平面方程用聚类后三角形集合中的所有三角形的顶点进行最小二乘拟合，计算出平面片的边界。

作为优选，其中平面片的提取过程包括以下步骤：a)从直方图上选择任一未访问的峰值，以其属性构造相似属性子集s1；b)从s1中任选一点作种子点；c)将该种子点从s1移到平面片候选集c中；d)对c中的每个成员，s1搜索其空间邻接三角形，并将他们移到c中，直到s1中找不到邻接三角形为止；e)重复步骤b)，c)，d)，直到s1为空为止；f)标记该直方图峰值为已访问，重复上述过程，直到所有峰值处理完毕。

发明有益的效果是：本发明的创新主要体现如下：一、将平面法向量由三维属性空间直接映射到一维空间，并从中选择具有统计显著性的种子点，二、在空间域进行区域增长的过程中，利用属性相似性约束，将区域增长范围限制在具有相似属性的子集中。以上创新有效地减少了计算机内存占用，同时也显著地提升了平面片提取方法的时间性能，对大范围、海量机载lidar点云的三维重建具有重要的实际意义。

附图说明

图1是本发明处理流程框图；

图2是空间域平面片提取流程图；

图3是实施例1中某多双坡结构屋顶平面片提取(建筑#1)示意图：(a)lidar点云及边界，(b)边界约束的delaunay三角剖分，(c)三角形法向量在高斯球上的分布，(d)线性化的法向量直方图，(e)三角形法向量在平面上的分布(法向量相似则颜色相同)，(f)空间域平面片聚类(平行平面片通过空间邻接关系分离)，(g)平面片提取结果，(h)参考影像；

图4是实施例2中平屋顶的平面片提取过程(建筑#2)示意图：(a)边界约束的delaunay三角剖分，(b)法向量一维直方图，(c)平面片提取结果(其中的小点为原始lidar点)，(d)参考影像。

具体实施方式

下面结合附图对本发明作进一步说明：

本发明用于从机载lidar数据中自动提取平面片的方法，是参数域与空间域结合的方法，即，在参数域选择种子点，再用这些种子点在空间域中进行区域增长，从而自动提取各平面片。主要步骤包括，(1)属性参数的选择以及计算，(2)参数域聚类选择种子点，(3)空间域区域增长分割出平面片，以及，(4)最小二乘平面参数估计。图1给出本方法的主要步骤及处理流程：

(1)属性参数选择

数学上，平面是没有厚度、无限延伸的、平的二维表面。在3d欧氏空间中，通常用平面上的任一点，以及垂直于平面的法向量(normalvector)来描述一个平面，平面的法向量则是我们通常意义上的平面的“朝向”。

显然，平面具有两个属性，即，位置(位置向量)以及朝向(法向量)。但对于平面片来说，都有确定的边界，或者说是具有确定的空间分布及(有限的)面积。事实上，共面的平面片都具有相同的平面属性，因此，还必须通过空间邻接关系来区分共面的平面片。

对机载lidar数据而言，点云是未结构化的点的集合，要计算各个点的属性就必须对点集进行结构化(确定邻域关系等)。对给定的无结构化的点集，delaunay三角剖分将点集的凸包用三角形格网化，即，对无结构化的机载lidar点云数据集，可利用delaunay三角剖分进行结构化，确定各点的邻接关系。点云三角剖分中的每个三角形都有确定的平面属性，且同一平面上的三角形具有相同的平面属性。相比于由邻域拟合曲面所确定的点的法向量，三角形的法向量更直观，更易计算。因此，本发明选择点云delaunay剖分三角形的(单位)法向量，作为平面片的属性用以在属性空间聚类选择种子点。

已知三角形在3d空间中三个顶点(非共线点)，p0(x0,y0,z0)，p1(x1,y1,z1)，p2(x2,y2,z2)，由这三点所确定的平面方程为：

n·(r-p)＝0

其中，n＝(p1–p0)×(p2–p0)，即，法向量为三个顶点所确定的三个向量中任意两个向量的叉积，p是p0，p1，或p2中的任意一点。所计算的法向量单位化后作为平面片的属性(聚类)参数。

(2)种子点选择

一方面，lidar数据中含有噪声，且各采样点有一定的测量误差，另一方面则是地物的表面非光滑(有一定的粗糙度)，就某一个三角形的法向量而言，用其作种子点不够稳定。为了获取更鲁棒的平面片种子点，其选择必须具有统计意义的显著性。直观地看，当三角形的单位法向量映射到高斯球(gaussiansphere)上时，具有相似属性的三角形(共面或平行)，其法向量指向相同(或相近)的方向，进而在球面上形成一个显著区域。从统计直方图上看，这些显著区域对应直方图的峰值。因此，本发明在对lidar点集进行三角剖分后，生成三角形单位法向量直方图，并检测直方图的峰值进而确定平面片种子点。高斯球上的显著区域，或直方图上的局部峰值，表示平面片最可能的朝向(法向量)。

为构造三角形法向量直方图，向量的每个分量必须量化为离散值(discretevalue)。量化步长(或直方图单元格大小，bincellsize)需要合理选择，以平衡(或满足)平面检测的精度以及可靠性要求。步长越小，参数量化精度越高，所对应提取的平面参数就越精确，但也更容易受噪声或外点影响。相应地，步长越小，量化累加器数组也越大；对于3d直方图来说，其占用的计算机内存也越大。如果注意到，在将三角形法向量映射到3d直方图时，其三维累加器数组中有很多空(null)值，因此，可以采用类似稀疏矩阵的存贮方式，将3维数组映射到一维空间，也即，用一个四元组来表示一个直方图单元格：

θ(i,j,k,value)→f(x,y,z)

其中，i,j,k对应数组下标，且i,j,k可以直接映射到一维空间，或一维数组(即，3d直方图线性化)。

(3)空间域平面片聚类分割

考虑到数据获取精度，量化步长以及计算误差等，本发明使用chebyshev距离定义“属性相似性”(attributesimilarity)，其阈值与量化步长相关。

定义(相似性)：

由定义可以看出，若两向量属性分量差值的最大值小于给定阈值，则相似。法向量直方图上的给定点p，由相似性定义所确定的邻域内的所有点都与p相似，即，若3d属性空间中的点到p的“距离”小于等于给定的阈值dt，则与p有相似属性。例如，若dt＝1，则p与其26个邻点(记为b)形成一个3×3×3的立方体，max(|px-bx|,|py–by|,|pz–bz|)≤dt＝1，即这26个邻点与p相似。

定义相似属性后，所有的(未处理)三角形集合s按照其法向量的属性(分量)划分为两个集合，即，候选集s1与非候选集s2。s1中的每个元素与直方图的峰值(种子点)属性相似，其他三角形则构成非候选集。基于空间邻接性的平面片区域增长过程，则限制在候选集s1中，而非整个集合s中。三角形的空间邻接性定义如下，即，两个邻接三角形之间至少有一个公共点，但不重叠。

定义(邻接性)：其中，“^”表示并且(and关系)。

构建法向量直方图后，选择任意一个与直方图局部峰值具有相似属性的三角形，开始平面片的区域增长过程(分割提取平面片)。在直方图中，每个局部峰值表示参数空间中相似的平面属性，即法向量相似，或空间域中近似平行或共面的平面片；因此，还需要在空间域分割这些相似属性集合，进而提取相应的平面片。平面片提取过程如下(参见图2)：

(a)从直方图上选择任一未访问的峰值，以其属性造相似属性子集s1；

(b)从s1中任选一点作种子点；

(c)将该种子点从s1移到平面片候选集c中(标号)；

(d)对c中的每个成员，在s1中搜索其空间邻接点(三角形)，并将他们移到c中，直到在s1中找不到邻接点为止(区域增长)；

(e)重复步骤(b)，(c)，(d)，直到s1为空为止(分割相同属性的多平面片)；

(f)标记该直方图峰值为已访问，重复上述过程，直到所有峰值处理完毕。

(4)平面方程计算

经过区域增长聚类后的三角形集合，属于同一平面片的三角形具有相同的标号，这些三角形合并后生成3d空间中的“多边形”。平面方程则用所有这些三角形的顶点进行最小二乘拟合，拟合方程的形式如下(机载lidar数据，平面不竖直，法线z方向的分量不为0)：

z＝ax+by+d

其最小二乘方程为，

上述方程的解是拟合平面(方程)的最小二乘解，拟合误差是z方向误差。所提取平面片的边界，则是聚类分割过程提取的“多边形”在拟合平面上的投影多边形的边界。

以下以某区域的机载lidar样例数据为实施例加以说明。三维激光点云数据采集于2010年，飞行平台平均航高500米，数据采集视场角(fov)为45°,点云的平均密度约6.3点/米²。

实施例1：

将本方法应用于一幢多双坡结构屋顶(建筑#1)的平面片提取过程如图3所示。从法向量分布图(图3(c)，两个显著区域)，或从法向量直方图(图3(d)，两个明显的峰值)可以看出，屋顶结构包括两组平行(或共面)的平面片(图3(e),各组颜色相同)。当用种子点在空间域进行聚类时，基于空间邻接关系(这些平面片相离)，这些平行平面片被相应地分离开(图3(f))。图3(g)中，红圈所示的面片，连接两个平行平面片，对应竖直平面(墙)，在后续的重建过程中需要加以利用。

实施例2：

将本方法应用于从另一平屋顶机载lidar数据集中提取平面片的过程(建筑#2)如图4所示。整个平屋顶被女儿墙分割成三个不同的平面片，另有少量的如烟囱等非平面点(图4(d))。由于由激光扫描系统获取的采样点，很难精确覆盖线特征，因此，仅从lidar数据上看，这几个平面片在空间上是相互连通的。此外，由于存在少量墙面点，这些点与屋顶平面点之间构成几个倾斜面，也即平面提取结果中所包含的几个“伪”平面片(图4(c))。这些“伪”平面片都具有特定的形状参数(面积小，面积/周长比小)，在重建过程中可以方便地加以剔除。

除上述实施例外，本发明还可以有其他实施方式。凡采用等同替换或等效变换形成的技术方案，均落在本发明要求的保护范围。