CN106570468A - 一种重建LiDAR原始点云建筑物轮廓线的方法 - Google Patents

Abstract

本发明提供了一种重建LiDAR原始点云建筑物轮廓线的方法，包括在LiDAR离散点云的基础上，按照Delaunay三角网的原则构建三角网格；将Delaunay三角网格在三维空间存在长边的三角形中高程值比较大的两个点提取出来，并将其以点对形式存储；然后将点对作为“数据链条节”，通过连接“数据链条节”构成环形数据链条，完成每个突起物边缘点的聚类；最后统计每个突起物点集内点的数量，筛除非建筑物边缘点集合；按照连接顺序将每个点逐次加入拟合数据中重新拟合，以新加入的点到前一步拟合直线的距离作为判断标准，最终将相邻拟合边的交点连接起来得到建筑物的轮廓线；本发明方法速度快、效率和准确度高，解决了LiDAR数据的离散性和非特征性以及轮廓线重建的规则性问题。

Description

一种重建LiDAR原始点云建筑物轮廓线的方法

技术领域

本发明涉及一种提取建筑物外边界轮廓线的方法，特别涉及重建机载LiDAR获得的离散点云数据中的建筑物屋顶轮廓的方法。

背景技术

建筑物轮廓信息在军事和民用许多相关领域有着越来越广泛的应用，特别是城市数字模型建立和规划管理、交通线路的设计、电力线路的勘测、地形的检测与分析、森林监管等方面，以及现代作战中精确目标打击和毁伤效果评估领域有着重要作用。

LiDAR是利用全球定位系统(GPS,Global Position System)和惯性测量装置(IMU,Inertial Measurement Unit)进行激光扫描获取目标三维信息的激光雷达探测系统。LiDAR系统探测得到的数据以离散点的形式存在，其中含有每个点的空间三维信息和激光回波强度信息。在LiDAR技术的运用过程中，建筑物轮廓信息的提取是最为关键及困难的一步。

徐景中等在“LiDAR点云中多层屋顶轮廓线提取方法研究”中是先从点云数据中滤波将屋顶点提取出来，而后进一步提取边缘点，这是许多建筑物重建数据驱动方法的基本流程，但是该类方法建筑物边缘点的提取精度大大依赖于从点云数据中提取建筑物的滤波精度；经典的Alpha shape算法能够很好的提取建筑物的边缘点，但是只能在单独的建筑物点集中运用、提取边缘点，无法实现LiDAR数据量大的场景中不同建筑物点的聚类。曾齐红等在“机载LiDAR点云数据的建筑物重建研究”中提出的利用Delaunay三角网公共边的方法也只能适用于单一的建筑物点集，不能实现不同建筑物点的聚类，实现数据的批量处理。程效军等在“融合航空影像和LiDAR点云的建筑物探测及轮廓提取”中分别提取机载点云和航空影像中的部分建筑物轮廓，将航空影像的光谱信息赋予机载点云，将光谱信息作为特征向量进行 聚类，利用高程信息从光谱信息相似的地面道路和建筑物中分离出建筑物，提取建筑物轮廓线，完成建筑物轮廓的探测。该类建筑物提取方法需要利用配准技术将航空影像信息与机载点云信息结合起来，虽然提取精度高，但计算量大，对信息源要求过高，需要探测同一场景的LiDAR数据及航空影像数据，没有充分利用机载点云数据信息。而另一类模型驱动的建筑物提取及重建方法需要建立类型种类健全的建筑物模型库，不仅工作量大，计算繁杂，并且难以涵盖实际情况中的所有建筑物模型。

LiDAR技术能够快速地获取目标空间三维信息，对获取的目标三维离散点信息进行轮廓线的拟合和区域的划分，其中存在的问题主要有三点：

1、三维信息的离散性问题。LiDAR技术获取的目标空间三维信息为离散点数据，点与点之间没有密切的联系，本应关联性最大的邻近点之间，在数据中无法确定相互位置关系。

2、目标数据的非特征性问题。对原始离散点数据很难利用现有的方法快速提取特征得到目标区域。

3、处理结果的规则性问题。从点云数据中提取建筑物边缘点后，对建筑物轮廓线进行规格化的关键就是确定其角点的位置，采用何种拟合方式和拟合精度确定角点的位置也是必须要解决的重要问题。

发明内容

为了克服现有技术的不足，本发明提供一种重建LiDAR原始点云建筑物轮廓线的方法，包括：

步骤A：在LiDAR离散点云的基础上，按照Delaunay三角网的构建原则建立三角网格；

步骤B：将在步骤A构建的三角网格中存在长边的三角形中高程值大的两个点提取出来，并将其以点对形式存储；然后将该点对作为数据链条节，根据数据链条节的连接关系形成数据链条节集合，完成每个突起物边缘点的聚类，形成突起物边缘点集合，最后统计每个突起物边缘点集合内点的数量，筛除非建筑物边缘点集合；

步骤C：构建每个建筑物边缘点集的环形数据链条，按照环形数据链条顺序将点加入拟合数据集合，进行直线拟合，当新加入的点到当前拟合直线的距离小于设定的阈值，则将该点加入当前拟合数据集合，重新拟合当前直线；当新加入的点到当前拟合直线的距离大于设定的阈值，则将该点视作下一个边的起始点，以该点与环形数据链条顺序中的下一点拟合新的直线；求出相邻两条拟合直线的交点并顺次连接，最终得到建筑物的轮廓线。

进一步地，

步骤A实现从原始LiDAR点数据中直接提取突起物的边缘点；

步骤A中输入为原始点云数据、阈值h₁和阈值h₂；步骤A中输出为突起物边缘点的总集合，步骤A包括：

步骤A1：数据预处理，去除奇异点；

步骤A2：采用Delaunay三角网格，构建点数据之间关联关系；

步骤A3：计算Delaunay三角网中三角形每个边的边长和平均边长；

步骤A4：提取突起物边缘点；

步骤B实现从突起物点集合中提取建筑物点并将同一建筑物的点聚类；

步骤B输入为突起物边缘点的总集合、阈值h₃和阈值h₄；步骤B输出为每个建筑物边缘点的集合，步骤B包括：

步骤B1：在边缘点集A中建立Deluanay三角网格；

步骤B2：建筑物边缘点聚类；

步骤B3：滤除突起物边缘点集合中的非建筑物边缘点集合；

步骤C实现将离散的边缘点拟合规格化，得到连续的建筑物轮廓线；步骤C输入为每一个建筑物边缘点的集合和残差阈值r，步骤C输出为每一个建筑物的规则轮廓线，步骤C包括：

步骤C1：形成每个建筑物边缘点集的闭合环路；

步骤C2：按照环形数据链条的顺序进行直线拟合；

步骤C3：计算残差；

步骤C4：确定拐点；

步骤C5：生成规则轮廓线。

进一步地，

步骤A1包括：

步骤A1-1：设定搜索距离；设三维空间x-y-z坐标系每个轴上的搜索距离为△x、△y和△z；

步骤A1-2：建立每个点的搜索长方体：以每一个点为中心，向X、Y和Z轴的正向和负向分别扩展△x、△y和△z的搜索距离，形成每个点的搜索长方体；

步骤A1-3：设定扩展长方体内点数量的阈值h₁，h₁＝1；

步骤A1-4：删除奇异点：统计落在每个点扩展长方体中点的数量a_i；如果a_i大于阈值h₁，则保留，如果a_i小于阈值h₁，则将该点从数据集中删除；最终得到新的点集O1；

选取作为△x、△y和△z的值，为平均相邻离散点间距；

步骤A2包括：

步骤A2-1：建立Delaunay三角形：在点集O1中任意寻找一点，然后查找距离此点最近的点，连接起来作为初始线；寻找一个离初始线最近的点，作为第三点；生成一个Delaunay三角形；

步骤A2-2：建立Delaunay三角网：以Delaunay三角形新生成的两条边，作为新的基线；重复步骤A2-1、步骤A2-2，直到所有的基线处理完毕，得到Delaunay三角网格；

步骤A4包括：

步骤A4-1：设定边长阈值h₂：设定Delaunay三角网格的平均边长为阈值h₂；

步骤A4-2：提取符合特征的三角形的顶点：将h₂与每个三角形的边长c进行比较，如果同一三角形内有两条边长超过阈值h₂，则将该三角形三个顶点(设为i，j，k)的坐标提取出来；

步骤A4-3：计算最大与最小高程点的高程平均值：取出三个点中高程最大点i和高程最小点k，并求这两个点高程的平均值

步骤A4-4：边缘点判定：比较第三个点的高程h_j与高程平均值的大小，如果则将i，j提取出来，并以点对的形式存储在新的集合E中；否则，判断下一个三角形，直到Delaunay三角网格内的三角形都判断完毕；

步骤A4-5：去除重复数据：将提取出来重复的点合并，得到突 起物边缘点的集合A；

步骤B2包括：

步骤B2-1：提取数据链条节：求三角网中每个三角形的边长 Δx′_i，Δy_i′，Δz′_i分别为构成三角形边长的两个端点的x轴、y轴、z轴方向上的差值；设定阈值h₃，h₃＝h₂；将边长长度大于h₃的边的两个端点提取出来，并以点对的形式存储，放在集合E中；

步骤B2-2：构建数据链条节集合，实现不同突起物边缘点的聚类：将起始点对放入到搜索集合L1中；以集合E的第一个点对为起始，放入到集合RB_i中，其中下标i代表第i个突起物边缘点集合；搜索集合E中内包含集合L1点中任意一点的点对，将搜索的结果以点的形式放入到临时集合L2中；删除集合E中搜索得到的点对；将L1里面的点更新为集合L2中集合RB_i中不存在的点；将更新后L1中的点放入集合RB_i中；重复上述步骤直到点集RB_i不再发生变化，即第i个突起物的点集不再更新；重复上述步骤直到将集合E变为空集，即所有点都分配到各自的突起物边缘点集合中；

步骤B3：滤除突起物边缘点集合中的非建筑物边缘点集合：统计每一个突起物点集RB_i中点的数量sum_i；运用Otsu法计算确定阈值h₄；如果点集RB_i中点的数量sum_i大于阈值h₄，则保留点集；如果点集RB_i中点的数量sum_i小于阈值h₄，则删除整个点集；

步骤C1包括：

步骤C1-1：分别建立每个点集RB_i的Delaunay三角网格；

步骤C1-2：选取数据链条连接起点：以点集RB_i中一个点P_i开始，并将其存到新的点集RR_i中；

步骤C1-3：寻找相同的链条节点：搜索三角网格中存在点P_i的三角形，将三角形的顶点放入集合G中；将集合G中仅出现一次的点放进集合H中；将H中的两个点插到点集RR_i中位置的第一个和最后一个；

将新插入的点作为新的搜索点，按照步骤C1-3的方法进行搜索，得到新的两个点；

步骤C1-4：链条连接：将点集RR_i中第一个点进行搜索得到的点插入到点集RR_i中第一个位置，将点集RR_i中最后一个点进行搜索得到的点插入到点集RR_i中最后一个位置；

步骤C1-5：重复步骤C1-3，步骤C1-4直到点集RR_i中第一个点与最后一个点相同；

步骤C3中定义点p_j+2到前面拟合直线的距离为残差d_j，点p_j+2坐标为(p_x，p_y，p_z)，得到的拟合直线方程为Lax+Lby+Lcz+Ld＝0，其中La、Lb、Lc、Ld是直线方程的参数，x,y,z是直线方程的变量，残差d_j计算公式为：

步骤C4中设定逐点拟合的残差阈值r，判断r与残差d_j的大小，如果r＜d_j，则将点p_j+2加入当前拟合数据集合，返回步骤C2重新拟合当前直线，如果r＞d_j，则将点p_j+2作为下一条边的起始，以点p_j+2与环形数据链条顺序中的下一点返回步骤C2拟合新的直线；

步骤C5中按照拟合得到直线的顺序，依次求得直线的交点X₁，X₂，X₃，……X_i；将X₁，X₂，X₃，……X_i依次连接，得到最终的轮廓线。

进一步地，

C1中闭合环路的形成方法为选取第一个点为起始点，不断搜索三角网格内由该点构成的边长，该边长只存在于一个三角形中；将搜索到的端点再定为搜索点，继续搜索下一个点，直到将整个点集的点连接成一个闭合的环路。

进一步地，

轮廓线规则化后加入建筑物边缘点高程信息。

本发明方法重建建筑物轮廓速度快、效率和准确度高，解决了LiDAR数据的离散性和非特征性以及轮廓线重建的规则性问题。

下面结合附图和实施例对本发明作详细说明。

附图说明

图1是本发明实施例中原始LiDAR数据图；

图2是本发明方法流程图；

图3是本发明实施例中提取突起物边缘点效果图；

图4是本发明实施例中建筑物数据链条聚类效果图；

图5是本发明实施例最终轮廓线效果图；

图6是本发明实施例中提取突起物边缘点运行时间图；

图7是本发明实施例中建筑物数据链条聚类运行时间图。

具体实施方式

以下将结合具体实施例对本发明提供的技术方案进行详细说明，应理解下述具体实施方式仅用于说明本发明而不用于限制本发明的范围。

本发明一个实施例中的原始LiDAR数据如图1所示。

本发明在现有技术和需要解决的关键技术问题基础上，提出一种提取LiDAR原始点云建筑物轮廓线的新方法，如图2所示，包括：

步骤A：边缘点提取。在LiDAR离散点云的基础上，按照Delaunay三角网的原则构建三角网格，建立数据二维空间的关联关系；

步骤B：建筑物提取和聚类。将在步骤A构建的三角网格中存在长边的三角形中高程值大的两个点提取出来，并将其以点对形式存储；然后将点对作为数据链条节，根据数据链条节的连接关系形成数据链条节集合，完成每个突起物边缘点的聚类；最后统计每个突起物点集内点的数量，筛除非建筑物边缘点集合；

步骤C：轮廓线规格化。构建每个建筑物边缘点集的环形数据链条，按照环形数据链条顺序将点加入拟合数据集合，进行直线拟合，当新加入的点到当前拟合直线的距离小于设定的阈值，则将该点加入当前拟合数据集合，重新拟合当前直线；当新加入的点到当前拟合直线的距离大于设定的阈值，则将该点视作下一个边的起始点，以该点与环形数据链条顺序中的下一点拟合新的直线；求出相邻两条拟合直线的交点并顺次连接，最终得到建筑物的轮廓线。

其中，步骤A实现从原始LiDAR点数据中直接提取突起物的边缘点。步骤A中输入为原始点云数据、阈值h₁和阈值h₂；步骤A中输出为突起物边缘点的总集合。具体实施过程如下：

步骤A1：数据预处理，去除奇异点：本发明针对LiDAR扫描时 产生的噪声或者回波，原始点云数据中形成脱离建筑物或地面的孤立点，称为奇异点，在后面边缘点提取的步骤中会造成误判，在边缘点提取之前利用三维空间中点的分布特征将这些奇异点剔除。具体实施过程为：

步骤A1-1：设定搜索距离：设三维空间x-y-z坐标系每个轴上的搜索距离为△x、△y和△z；

步骤A1-2：建立每个点的搜索长方体：以每一个点为中心，向X、Y和Z轴的正向和负向分别扩展△x、△y和△z的搜索距离，形成每个点的搜索长方体；

步骤A1-3：设定扩展长方体内点数量的阈值h₁，h₁＝1；

步骤A1-4：删除奇异点：统计落在每个点扩展长方体中点的数量a_i；如果a_i大于阈值h₁，则保留，如果a_i小于阈值h₁，则将该点从数据集中删除，得到新的点集O1；

本发明采用理论分析法确定步骤A1-1中搜索距离△x，△y和△z以及A1-3中点数量的阈值h₁：

Δx、Δy、Δz分别是每个点在x轴、y轴、z轴上的搜索距离，以点(x_i,y_i,z_i)为例，建立(x_i±Δx，y_i±Δy,z_i±Δz)的三维长方体。Δx、Δy、Δz的设定要大于平均相邻离散点间距使生成的长方体能够在x、y、z方向分别包含对应轴上正负方向的一个相邻点。本发明步骤A1-1中选取作为Δx、Δy和Δz的值，即平均相邻离散点间距取决于LiDAR系统扫描特性，本实施例中的值设定为1m。

奇异点的出现概率比较低，一般是以独立点的形式存在的，因此将阈值h₁设定为数量“1”。

步骤A2：构建点数据之间关联关系：本发明针对LiDAR点云数据的离散性特征，为了在数据库存储中，解决点与点之间无关联性的问题，该问题导致无法快速查询某一点在空间中的邻近点，并且由于LiDAR点云数据量极大，需要根据每个点及其邻域点的空间位置关联对数据进行处理。本发明采用Delaunay三角网格来建立一种有效的数据结构，以提高数据的搜索效率及可靠性。Delaunay三角网格的数据结构建立比较快，并且其表示方法没有数据冗余，更便于找到每个点的邻近点，并能够求得每个点的旋转度以及邻近点三角面片的法 向量等特性。具体实施过程为：

步骤A2-1：建立Delaunay三角形：在点集O1中任意寻找一点，然后查找距离此点最近的点，连接起来作为初始线；寻找一个离初始线最近的点，作为第三点；生成一个Delaunay三角形；

步骤A2-2：建立Delaunay三角网：以Delaunay三角形的两条新边(新生成的两条边)，作为新的基线；重复步骤A2-1、步骤A2-2，直到所有的基线处理完毕，得到Delaunay三角网格；

步骤A3：计算Delaunay三角网中三角形每个边的边长和平均边长：Delaunay三角网中三角形每个边的边长公式：Δx_i，Δy_i，Δz_i分别为构成三角形边长的两个端点的x轴、y轴、z轴方向上的差值绝对值；统计Delaunay三角网中三角形边的数量n；计算Delaunay三角网中三角形的平均边长

基于二维空间建立Delaunay三角网格并求Delaunay三角网格在三维空间中的边长，使得点距特征的应用比较明确，可以在原始点云基础上直接对突起物的边缘点进行提取。方法的目的性强，处理步骤直接，具有较快的反应速度，即良好的响应性。

步骤A4：提取突起物边缘点：本发明根据在突起物边缘处，相邻点的高程差Δz较大，没有突起物的时候，相邻点的高程差Δz接近于0的特征，比较三角形中的边长，如果同一个三角形里面有两个边长超过阈值h₂，并且三角形三个顶点里面有两个点的高程比较高，一个顶点的高程比较低，则将三个顶点里面高程比较高的两个点提取出来，以点对的形式存储，两个高程较高的点即为突起物的边缘点，并且是邻接点，将其提取，具体实施过程为：

步骤A4-1：设定边长阈值h₂：设定Delaunay三角网格的平均边长为阈值h₂；

步骤A4-2：提取符合特征的三角形的顶点：将h₂与每个三角形的边长c进行比较，如果同一三角形内有两条边长超过阈值h₂，则将三角形的三个顶点(设为i，j，k)坐标提取出来；

步骤A4-3：计算最大与最小高程点的高程平均值：取出三个点中高程最大点i和高程最小点k，并求这两个点高程的平均值

步骤A4-4：边缘点判定：比较第三个点的高程h_j与高程平均值的大小，如果则将i，j提取出来，并以点对的形式存储在新的集合E中；否则，判断下一个三角形，直到Delaunay三角网格内的三角形都判断完毕；

步骤A4-5：去除重复数据：将提取出来重复的点合并，得到突起物边缘点的集合A。

本发明采用理论分析法确定步骤A4-1中边长的阈值h₂：

在突起物边缘处，相邻点的高程差Δz较大，没有突起物的时候，相邻点的高程差Δz接近于0。当突起物出现的时候，边长比较大，即长边；没有突起物出现时候的边长比较小，即短边。设短边的数量为a，平均边长为长边的数量为b，平均边长为则平均边长其中有：由于a＞＞b，平均值略大于平均点距又低于有突起物出现时长边的平均边长 因此，通过平均边长做阈值，能够明确地检测出突起物的边缘点。

如图3所示，建筑物边缘点提取结果在用黑色的点表示。

步骤B实现从突起物点集合中提取建筑物点并将同一建筑物的点聚类。步骤B输入为突起物边缘点的总集合、阈值h₃和阈值h₄；步骤B输出为每个建筑物边缘点的集合。具体实施过程如下：

步骤B1：在边缘点集A中建立Deluanay三角网格：实施方法与步骤2相同；

步骤B2：建筑物边缘点聚类：本发明依据提取出的边缘点构建的Delaunay三角网，在Delaunay三角网中，根据同一突起物边缘点之间构成的三角形的边，其长度较短的特征，将同一突起物的边缘点区分开。把每一个点对看作一个“链条节”，将具有相同节点的链条节连接起来形成一个链条，称为数据链条，点对中的每一个点就是数据链条中的节点。每一个完整的数据链条就是一个突起物边缘点的类。

Delaunay三角网构成的点的连接具有很强的关联性，基于此的 数据链条连接方法充分利用了这一特征，得到的边缘点环形链条可靠性强，聚类效果好。具体实施过程为：

步骤B2-1：提取数据链条节：求每个三角网的边长 Δx′_i，Δy′_i，Δz′_i分别为构成三角形边长的两个端点的x轴、y轴、z轴方向上的差值。设定阈值h₃，h₃＝h₂；将大于h₃两个端点提取出来，并以点对的形式存储，放在集合E中；

步骤B2-2：构建数据链条节集合，实现不同突起物边缘点的聚类：将起始点对放入到搜索集合L1中；以集合E的第一个点对为起始，放入到集合RB_i中(其中下标i代表第i个突起物的点集)；搜索集合E中内包含集合L1点中任意一点的点对，将搜索的结果以点的形式放入到临时集合L2中；删除集合E中搜索得到的点对；将L1里面的点更新为集合L2中集合RB_i中不存在的点；将更新后L1中的点放入集合RB_i中；重复上述步骤直到点集RB_i不再发生变化；重复上述步骤直到将集合E变为空集；

步骤B3：滤除突起物边缘点集合中的非建筑物边缘点集合：统计每个数据链条节点个数，即统计每一个突起物集合RB_i中点的数量sum_i；运用Otsu法计算确定阈值h₄；如果点集RB_i中点的数量sum_i大于阈值h₄，则保留点集；如果点集RB_i中点的数量sum_i小于阈值h₄，则删除整个点集。

本发明采用理论分析法确定步骤B2-1中边长的阈值h₃：

在提取出的边缘点点集基础上建立的Delaunay三角网中，同一突起物上的相邻边缘点之间构成的边长与平均点距近似，且与非相邻边缘点之间构成的边长长度差别很大，因此，本发明设置h₃＝h₂，将小于阈值h₃的边长保留，其端点以点对形式存储，大于阈值h₃的边长删除。

通过Delaunay三角网格中点的属性特征得到了突起物的边缘点和边缘点的连接关系，突起物中包含建筑物点与非建筑物点两种。本发明根据先验知识：建筑物的面积一般要大于其它突起物的面积，如树木，汽车等，将面积反映到边长上，则面积比较大的周长较大，建筑物轮廓线的周长要大于其它突起物的周长。在点距相当的前提下，建筑物具有更多的边缘点，即建筑物边缘点的数据链条更长，节点更多，利用这个特征，通过边缘点的连接关系，能够得到每个突起物边 缘点连接情况和边缘点数量，从而很好地过滤掉突起物中的植被点跟小型突起物，以此将建筑物边缘点集合与非建筑物边缘点集合区分开。

本发明根据上述先验知识及Otsu算法确定步骤B2-3中建筑物点集数量的阈值h₄：

聚类得到的边缘点点集是建筑物边缘点集和非建筑物边缘点集之和，本发明这一步目的是为了剔除掉非建筑物边缘点集，分类的依据为边缘点集内点的数量，将边缘点集分为两类。运用Otsu方法能够得到保证最大类间距离的单个阈值，很好地将其分为两类，因此在阈值选取的过程中，本发明采用Otsu方法得到阈值h₄。

滤除非建筑物边缘点集合之后并进行聚类，其结果如图4所示，用不同形状的点表示不同建筑物上的边缘点。

步骤C基于最小二乘法，采用了逐点拟合的策略，按照数据链条的顺序，逐个点求残差判断拐点的位置，使得准确率更高。

步骤C实现将离散的边缘点拟合规格化得到连续的建筑物轮廓线。步骤C输入为每一个建筑物边缘点的集合和残差阈值r，步骤C输出为每一个建筑物的规则轮廓线。具体步骤如下：

步骤C1：形成每个建筑物边缘点集的闭合环路，即生成数据链条连接顺序，具体实施过程如下：

步骤C1-1：分别建立每个边缘点集RB_i的Delaunay三角网格：实施方法与步骤2相同；

步骤C1-2：选取数据链条连接起点：以点集RB_i中一个点P_i开始，并将其存到新的点集RR_i中；

步骤C1-3：寻找相同的链条节点：搜索三角网格中存在点P_i的三角形，将三角形的顶点放入集合G中；将集合G中仅出现一次的点放进集合H中；将H中的两个点插到点集RR_i中位置的第一个和最后一个；

将新插入的点作为新的搜索点，按照步骤C1-3的方法进行搜索，得到新的两个点；

步骤C1-4：链条连接：将RR_i中第一个点进行搜索得到的点插入到RR_i中第一个位置，将RR_i中最后一个点进行搜索得到的点插入到RR_i中最后一个位置；

重复步骤C1-3，步骤C1-4直到RR_i中第一个点与最后一个点相同；

步骤C2：按照环形数据链条的顺序进行拟合：在RR_i中选取一个点p_j，将点p_j和下一个点p_j+1利用最小二乘法拟合；

步骤C3：计算残差。定义点p_j+2到前面拟合直线的距离为残差d_j，点p_j+2坐标为(p_x，p_y，p_z)，得到的拟合直线方程为Lax+Lby+Lcz+Ld＝0，其中La、Lb、Lc、Ld是直线方程的参数，x,y,z是直线方程的变量，残差d_j计算公式为：

步骤C4：确定拐点：设定逐点拟合的残差阈值r，判断r与残差d_j的大小，如果r＜d_j，则将点p_j+2加入当前拟合数据集合，返回步骤C2重新拟合当前直线，如果r＞d_j，则将点p_j+2看作是下一条边的起始，以点p_j+2与环形数据链条顺序中的下一点返回步骤C2拟合新的直线；

步骤C5：规则轮廓线的生成：按照拟合得到直线的顺序，依次求得直线的交点X₁，X₂，X₃，……X_i；将X₁，X₂，X₃，……X_i依次连接，得到最终的轮廓线。

本发明上述步骤C1中闭合环路的形成，即数据链条的连接顺序设定：选取第一个点为起始点，不断搜索三角网格内由该点构成的边长，并且该边长只存在于一个三角形中；将搜索到的端点再定为搜索点，继续搜索下一个点，直到将整个点集的点连接成一个闭合的环路。本发明采取了一个“后补”起始点的策略：先随机选取一个点开始进行拟合，当出现第一个拐点时将其定为轮廓线拟合的起始点。因为轮廓线是一个闭合环路，最终的拟合会回到起始点。这样能够避免出现由于起始点选择不同导致最后拟合结果不同的问题，提高了拟合的准确性。

本发明采用理论分析法设定步骤C4-1中残差阈值r：

逐点拟合中残差阈值r为局部点允许波动的范围，如果r设置得比较大，轮廓线的拟合精确度将会降低；如果r设置得比较小，点数据局部的波动将会影响轮廓线规格化的效果。按照规则的矩形建筑物 边缘点分布情况来看，局部波动比较大的位置是直角边拐点处，此时点的波动距离约等于因此，波动的值r的设定要小于

本发明步骤C3中采用最小二乘法作为边缘点轮廓线拟合方法：按照数据链条的连接顺序将数据点依次加入拟合点集合，进行拟合，以下一点到拟合直线的距离作为每条边分段的标准，当下一点与当前拟合直线的距离超过设定的阈值r时，则将其视为另一条边的开始，上一个点则看作是上一条边的结束；当距离没有超过设定的阈值时，则仍看作是这条边上的点，并加入下次的直线拟合。对数据链条的连接顺序集合中每一个点按连接顺序依照上述准则进行判断，得到每条边上点的聚类。通过求两条拟合直线的交点得到建筑物轮廓线中的拐点，按照顺序依次连接起来，得到建筑物规格化的轮廓线。

轮廓线规则化后加入建筑物边缘点高程信息，其结果如图5所示，用黑色的线表示该场景中的建筑物轮廓线。

本发明能够利用特征快速地将突起物的边缘点从原始点云中提取出来，当点云数量变大的时候，本发明依然能够保持较快的速度。如图6所示，在本发明的一个实施例中，原始数据采用不同数量数据点集，如1000，3000，5000，10000等，对本发明方法的效果进行对比分析。当处理1000个点时，运行时间为1.0966秒，当处理3000个点时，运行时间为2.2362秒，当处理5000个点时，运行时间为3.4148秒，当处理10000个点时，运行时间为7.0448秒，当处理50000个点时，运行时间为70.3564秒，当处理该场景全部点52119个点时，运行时间为74.8244秒；由此可见运行时间与处理的点数量基本成正比。当突起物出现的时候，对边长起主要作用的为高程变化，因此当LiDAR的扫描点密度发生变化的时候，本发明的提取效果不会随之发生变化。由于边长阈值的选取是利用所有三角形边长求平均长度得到的，因此，只要LiDAR的扫描距离是相对稳定就能得提取到边缘点所需的合适阈值。

本发明边缘点聚类速度快。在本发明的另一个实施例中采用不同数量数据点集测试边缘点聚类速度，如图4所示，滤除非建筑物边缘点集合之后并进行聚类，用不同形状的点表示不同建筑物上的边缘点。如图7所示，对本发明方法的效果进行对比分析。当处理1000个点时，运行时间为1.3369秒，当处理2000个点时，运行时间为2.6999 秒，当处理3000个点时，运行时间为4.5934秒，当处理4000个点时，运行时间为6.5391秒，当处理5000个点时，运行时间为8.7179秒，当处理该场景全部边缘点对5277个点对时，运行时间为8.8169秒；由此可见运行时间与处理的点数量基本成正比，这是由于：1、数据链条连接所利用的聚类特征明确，即边缘点邻接点对的关联性极强，仅需要点对的搜索，无需进行多余的计算；2、对数据链条的连接进行了并行搜索的方式，搜索策略是对数据链条的两个端点同时进行搜索，使得每次搜索可得到两个数据链条节点，大幅提升了搜索效率。

本发明还能够准确重建出建筑物规则轮廓线，这是由于：1、采用逐点拟合的策略，能够准确地找到轮廓线的拐点，即下一条边的起始点，实现了轮廓线每条边上点的聚类；2、采用最小二乘法拟合，在每条边上点的聚类的基础上，最小二乘法能够使得拟合残差和最小，精度高；3、轮廓线角点定位准确，通过求相邻拟合直线的交点，可以得到轮廓线的角点，形成一个直线构成的闭合轮廓线，比较符合建筑物的真实情况。如图5所示，轮廓线规则化后加入建筑物边缘点高程信息，黑色的线表示该场景中的建筑物轮廓线。

综上所述，本发明的有益效果有如下几点：

(1)本发明算法重建建筑物轮廓速度快，效率高，准确度高。

(2)基于Delaunay三角网格在二维空间中建立数据之间的联系，由于本发明利用的特征为点的空间特征，相比于其它数据结构，Delaunay三角网格更便于对数据的搜索，效率更高，解决了LiDAR点云数据的离散性问题；

(3)在Delaunay三角网格基础上，本发明利用了三维空间中Delaunay三角网的边长特征和数据链条的连接关系，即利用相邻点高程的突变性和建立的数据链条节点关系实现对建筑物边缘点的提取，精度高，解决了目标数据的非特征性问题；

(4)由于采用了逐点拟合的策略，在寻找拐点的过程中能够比较准确的检测出拐点的位置，实现对轮廓线边的拟合，确定出轮廓线角点的位置，解决了轮廓线重建的规则性问题；

(5)相比于其它轮廓线重建的方法，本发明无需进行前期对突起物点提取的滤波处理，能够直接在LiDAR原始点云中将建筑物的边 缘点自动提取出来。本发明采用的聚类方法不用提前设定聚类的种类数量，在种类数量不明确的前提下，本发明依靠于Delaunay三角网格得到的连通点对关系，采用了数据链条连接的方法对各个点进行聚类，找到突起物边缘点的数据链条并统计其节点的数量。为下一步非建筑物点的剔除做准备。实现每个建筑物边缘点的聚类，通过拟合最终得到建筑物完整的轮廓线。

Claims (5)

1.一种重建LiDAR原始点云建筑物轮廓线的方法，包括：

步骤A：在LiDAR离散点云的基础上，按照Delaunay三角网的构建原则建立三角网格；

步骤B：将在步骤A构建的三角网格中存在长边的三角形中高程值大的两个点提取出来，并将其以点对形式存储；然后将该点对作为数据链条节，根据数据链条节的连接关系形成数据链条节集合，完成每个突起物边缘点的聚类，形成突起物边缘点集合，最后统计每个突起物边缘点集合内点的数量，筛除非建筑物边缘点集合；

步骤C：构建每个建筑物边缘点集的环形数据链条，按照环形数据链条顺序将点加入拟合数据集合，进行直线拟合，当新加入的点到当前拟合直线的距离小于设定的阈值，则将该点加入当前拟合数据集合，重新拟合当前直线；当新加入的点到当前拟合直线的距离大于设定的阈值，则将该点视作下一个边的起始点，以该点与环形数据链条顺序中的下一点拟合新的直线；求出相邻两条拟合直线的交点并顺次连接，最终得到建筑物的轮廓线。

2.如权利要求1所述的一种重建LiDAR原始点云建筑物轮廓线的方法，其特征在于：

步骤A实现从原始LiDAR点数据中直接提取突起物的边缘点；

步骤A中输入为原始点云数据、阈值h₁和阈值h₂；步骤A中输出为突起物边缘点的总集合，步骤A包括：

步骤A1：数据预处理，去除奇异点；

步骤A2：采用Delaunay三角网格，构建点数据之间关联关系；

步骤A3：计算Delaunay三角网中三角形每个边的边长和平均边长；

步骤A4：提取突起物边缘点；

步骤B实现从突起物点集合中提取建筑物点并将同一建筑物的点聚类；

步骤B输入为突起物边缘点的总集合、阈值h₃和阈值h₄；步骤B 输出为每个建筑物边缘点的集合，步骤B包括：

步骤B1：在边缘点集A中建立Deluanay三角网格；

步骤B2：建筑物边缘点聚类；

步骤B3：滤除突起物边缘点集合中的非建筑物边缘点集合；

步骤C实现将离散的边缘点拟合规格化，得到连续的建筑物轮廓线；步骤C输入为每一个建筑物边缘点的集合和残差阈值r，步骤C输出为每一个建筑物的规则轮廓线，步骤C包括：

步骤C1：形成每个建筑物边缘点集的闭合环路；

步骤C2：按照环形数据链条的顺序进行直线拟合；

步骤C3：计算残差；

步骤C4：确定拐点；

步骤C5：生成规则轮廓线。

3.如权利要求2所述的一种重建LiDAR原始点云建筑物轮廓线的方法，其特征在于：

步骤A1包括：

步骤A1-1：设定搜索距离；设三维空间x-y-z坐标系每个轴上的搜索距离为Δx、Δy和Δz；

步骤A1-2：建立每个点的搜索长方体：以每一个点为中心，向x、y和z轴的正向和负向分别扩展Δx、Δy和Δz的搜索距离，形成每个点的搜索长方体；

步骤A1-3：设定扩展长方体内点数量的阈值h₁，h₁＝1；

步骤A1-4：删除奇异点：统计落在每个点扩展长方体中点的数量a_i；如果a_i大于阈值h₁，则保留，如果a_i小于阈值h₁，则将该点从数据集中删除；最终得到新的点集O1；

选取作为Δx、Δy和Δz的值，为平均相邻离散点间距；

步骤A2包括：

步骤A2-1：建立Delaunay三角形：在点集O1中任意寻找一点，然后查找距离此点最近的点，连接起来作为初始线；寻找一个离初始线最近的点，作为第三点；生成一个Delaunay三角形；

步骤A2-2：建立Delaunay三角网：以Delaunay三角形新生成的两条边，作为新的基线；重复步骤A2-1、步骤A2-2，直到所有的基线处理完毕，得到Delaunay三角网格；

步骤A4包括：

步骤A4-1：设定边长阈值h₂：设定Delaunay三角网格的平均边长为阈值h₂；

步骤A4-2：提取符合特征的三角形的顶点：将h₂与每个三角形的边长c进行比较，如果同一三角形内有两条边长超过阈值h₂，则将该三角形三个顶点(设为i，j，k)的坐标提取出来；

步骤A4-3：计算最大与最小高程点的高程平均值：取出三个点中高程最大点i和高程最小点k，并求这两个点高程的平均值

步骤A4-4：边缘点判定：比较第三个点的高程h_j与高程平均值的大小，如果则将i，j提取出来，并以点对的形式存储在新的集合E中；否则，判断下一个三角形，直到Delaunay三角网格内的三角形都判断完毕；

步骤A4-5：去除重复数据：将提取出来重复的点合并，得到突起物边缘点的集合A；

步骤B2包括：

步骤B2-1：提取数据链条节：求三角网中每个三角形的边长Δx′_i，Δy′_i，Δz′_i分别为构成三角形边长的两个端点的x轴、y轴、z轴方向上的差值；设定阈值h₃，h₃＝h₂；将边长长度大于h₃的边的两个端点提取出来，并以点对的形式存储，放在集合E中；

步骤B2-2：构建数据链条节集合，实现不同突起物边缘点的聚类：将起始点对放入到搜索集合L1中；以集合E的第一个点对为起始，放入到集合RB_i中，其中下标i代表第i个突起物边缘点集合；搜索集合E中内包含集合L1点中任意一点的点对，将搜索的结果以点的形式放入到临时集合L2中；删除集合E中搜索得到的点对；将L1里面的点更新为集合L2中集合RB_i中不存在的点；将更新后L1中的点放入集合RB_i中；重复上述步骤直到点集RB_i不再发生变化，即第i个突起物的点集不再更新；重复上述步骤直到将集合E变为空集，即所有点都分配到各自的突起物边缘点集合中；

步骤B3：滤除突起物边缘点集合中的非建筑物边缘点集合：统计每一个突起物点集RB_i中点的数量sum_i；运用Otsu法计算确定阈值h₄；如果点集RB_i中点的数量sum_i大于阈值h₄，则保留点集；如果点集RB_i中点的数量sum_i小于阈值h₄，则删除整个点集；

步骤C1包括：

步骤C1-1：分别建立每个点集RB_i的Delaunay三角网格；

步骤C1-2：选取数据链条连接起点：以点集RB_i中一个点P_i开始，并将其存到新的点集RR_i中；

步骤C1-3：寻找相同的链条节点：搜索三角网格中存在点P_i的三角形，将三角形的顶点放入集合G中；将集合G中仅出现一次的点放进集合H中；将H中的两个点插到点集RR_i中位置的第一个和最后一个；

将新插入的点作为新的搜索点，按照步骤C1-3的方法进行搜索，得到新的两个点；

步骤C1-4：链条连接：将点集RR_i中第一个点进行搜索得到的点插入到点集RR_i中第一个位置，将点集RR_i中最后一个点进行搜索得到的点插入到点集RR_i中最后一个位置；

步骤C1-5：重复步骤C1-3，步骤C1-4直到点集RR_i中第一个点与最后一个点相同；

步骤C3中定义点p_j+2到前面拟合直线的距离为残差d_j，点p_j+2坐标为(p_x,p_y,p_z)，得到的拟合直线方程为Lax+Lby+Lcz+Ld＝0，其中La、Lb、Lc、Ld是直线方程的参数，x,y,z是直线方程的变量，残差d_j计算公式为：

步骤C4中设定逐点拟合的残差阈值r，判断r与残差d_j的大小，如果r＜d_j，则将点p_j+2加入当前拟合数据集合，返回步骤C2重新拟合当前直线，如果r＞d_j，则将点p_j+2作为下一条边的起始，以点p_j+2与环形数据链条顺序中的下一点返回步骤C2拟合新的直线；

步骤C5中按照拟合得到直线的顺序，依次求得直线的交点X₁，X₂，X₃，……X_i；将X₁，X₂，X₃，……X_i依次连接，得到最终的轮廓线。

4.如权利要求2所述的一种重建LiDAR原始点云建筑物轮廓线的方法，其特征在于：C1中闭合环路的形成方法为选取第一个点为起始点，不断搜索三角网格内由该点构成的边长，该边长只存在于一个三角形中；将搜索到的端点再定为搜索点，继续搜索下一个点，直到将整个点集的点连接成一个闭合的环路。

5.如权利要求1至4中任意权利要求所述的一种重建LiDAR原始点云建筑物轮廓线的方法，其特征在于：轮廓线规则化后加入建筑物边缘点高程信息。