CN113205529B - 一种基于机载LiDAR点云分割建筑物顶面的方法 - Google Patents

Abstract

本发明提出一种基于机载LiDAR点云分割建筑物顶面的方法，很好的实现了房屋顶面的高精度提取，并且对于不同的数据类型具备很强的鲁棒性。本发明首先对原始建筑物点云进行欧式聚类，以分割出每栋建筑物的点云；然后分别对每一栋建筑进行基于L0梯度最小化的多尺度分割，得到初始平面分割结果。针对房屋顶面提取工作的特殊性我们对L0梯度最小化算法进行了针对性改进，包括排序和区域扩张的策略的应用以及附加的几何约束条件，很好的改进了L0梯度最小化算法在平面分割中的精度和鲁棒性。最后采用基于图割的后处理优化算法，将初始结果进行精化，同时解决初始结果中可能存在的过度分割以及边缘锯齿状的问题以得到最终的房屋顶面提取结果。

Description

一种基于机载LiDAR点云分割建筑物顶面的方法

技术领域

本发明涉及一种基于机载LiDAR点云利用L0梯度最小化和图割全局优化进行房屋顶面提取的方法，可用于LiDAR点云分类分割、三维建筑物模型重建等。

背景技术

建筑物的三维重建是计算机视觉、摄影测量与遥感等领域内的重要研究课题。随着激光探测及测距(Light Detection and Ranging)***的快速发展，迅速且准确的获得大规模城市三维点云成为可能，这为三维建筑物的重建工作带来了极大的便利。现有的基于机载LiDAR点云的建筑物三维建模的方法大致可以分为两类：模型驱动方法以及数据驱动方法。模型驱动方法是一种自上而下的方法并且认为建筑物模型都是存在于预设的模型库中的。而数据驱动方法则是自下而上进行建筑物建模，首先需要从点云中提取建筑物顶面，然后分析建筑物顶面的拓扑结构来最终构建完整的建筑物模型。由于模型驱动方法无法准确还原模型库中所不存在的房屋模型结构，因此从数据出发的数据驱动方法更为有效。而对于数据驱动方法，如何准确且高效地从点云中提取建筑物的顶面成为了三维建模工作的关键。针对建筑物顶面提取问题，目前国内外相关研究层出不从，但是由于不同数据点云密度不同，房屋复杂程度不同，现有的方法无法获得稳定且高精度房屋顶面提取结果，因此需要设计一种高精度且具备鲁棒性的房屋顶面提取算法来进一步推进建筑物模型重建工作。

发明内容

本发明针对现有技术的不足，提供一种基于L0梯度最小化算法以及图割优化算法的房屋顶面提取方法，该方法在面对不同类型数据时都可以得到高精度的房屋顶面提取结果。

本发明的技术方案为一种基于机载LiDAR点云分割建筑物顶面的方法，包括如下步骤：

步骤1，对原始建筑物点云进行欧式聚类，以分割出每栋建筑物的点云；

步骤2，分别对每一栋建筑进行基于区域扩张的L0梯度最小化算法的多尺度分割以得到初始平面分割结果；

步骤3，采用基于图割的后处理优化算法将初始平面分割结果进行精化，提取高精度建筑物顶面。

进一步的，步骤2中所述的基于区域扩张的L0梯度最小化算法是对L0梯度最小化算法的改进，其具体原理如下；

基于区域扩张的L0梯度最小化算法的能量函数定义为：

式中，M代表输入点云的长度，I_i和S_i分别代表点云中第i个点p_i输入和输出的法向量；N_i代表p_i的邻域，为k个最邻近的邻域点或者是三维空间中的以R为半径的球内的邻域点，i和j分别表示邻域点的索引，参数λ为L0范数的权重；

对公式(2)进行近似化求解，首先以一对相邻点p_i和p_j为例，它们对于整体能量函数F的贡献值f为：

为了使f最小，采取两种不同的策略：不融合p_i与p_j，此时S_i≠S_j或者融合p_i与p_j，此时S_i＝S_j，并且两个相邻点之间的L0范数被消除；

对于S_i≠S_j的情况，具体做法为：

对于S_i＝S_j的情况，具体做法为：

将两者情况相结合得到：

然后加入几何约束来进一步提升算法的稳定性，首先引入一个平面距离的概念，以一对相邻点p_i和p_j为例，它们之间的平面距离Pd_i,j定义如下：

Pd_i,j＝abs[p_ij·I_i] (6)

其中P_ij代表从p_i到p_j的向量，I_i代表p_i的法向量，abs[]为取绝对值；因此Pd_i,j代表着两个点在法向量I_i方向上的投影距离，当两个点之间的平面距离大于规定阈值T_pd，就不再惩罚它们之间的非零梯度；

因此公式6重写为：

接下来，把上述公式从相邻点扩展至相邻点集，以G_i代表一个点集，Y_i代表其法向量，点集中点数为w_i，点集G_i与G_j之间的连接点个数为c_i,j，结合上述几何约束条件相邻点集之间的能量公式为：

上述公式解为：

公式(9)即为基于区域扩张的L0梯度最小化算法的扩张条件。

进一步的，步骤2的基于区域扩张的L0梯度最小化算法的流程如下，

1)计算待分割点云中每个点的局部平滑特征，按照局部平滑特征值对点的操作顺序进行排序；

2)初始化点集并设定一个预设参数β，对于每一个点都看作一个只有一个点的点集；

3)引入辅助参数λ，并且每次迭代中λ的值增加k倍；辅助参数的初始化由下述公式完成：

其中median{}代表取中值；

4)在每一次迭代过程中，按顺序分别处理所有点集，对于点集G_i，根据扩张条件来确定是否将点集G_i的所有邻域点集G_j加入到G_i，如果G_j加入G_i，相对应的w_i将会被更新，G_j的邻域将会被加入N_i，完成一次扩张之后，G_j相关信息将会被删除，待处理点集个数减少1个；

5)循环处理所有点集，直至辅助参数λ大于预设参数β，迭代终止，剩余的点集则对应着分割出的平面，并且每个点集中的点都有着一致的平面法向量。

进一步的，所述局部平滑特征为该点与其邻域中其他邻域点法向量角度差异值的最大值。

进一步的，步骤2中多尺度分割的流程为，

(1)选取一个起始的分割尺度参数β，以欧式聚类得到的点云聚类为输入，其中分割尺度参数β与预设参数相同；

(2)对于每一个聚类首先进行平面验证，如果平面指数大于阈值，则直接输出为合法平面，否则即在当前分割尺度下通过基于区域扩张的L0梯度最小化算法进行平面分割以获得新的聚类；

平面指数的定义为：对一个聚类，先利用主成分分析法PCA来拟合平面，然后计算聚类中每个点到拟合平面的距离d_i，平面指数则为所有点到拟合平面的平均距离；

(3)当所有聚类完成(2)的处理，如果产生新的聚类，则以新生成聚类为输入进行下一次迭代，并且分割尺度参数减小；如果没有产生新的聚类即所有输入聚类都为合法平面，算法终止。

进一步的，步骤3中基于图割的后处理流程是以多尺度分割结果为初始平面通过图割算法来优化目标函数以得到每个点到平面的最优对应关系，所优化的目标函数由三部分组成，分别为：数据项，平滑项以及标签项，定义如下：

其中L为平面，P为平面L上的所有点的点集，L_p即代表点p所属的平面，公式中p,q∈N代表相邻的两点；

其中数据项是衡量点到平面的距离，计算公式为：

其中

为L_p的平面方程；

平滑项用来惩罚属于不同平面的相邻点的个数，δ()为指示函数，当L_p≠L_q，这个函数值为1，ω是权重；

最后一项是标签项，这一项用于惩罚最终所有点对应不同平面的个数，κ为权重。

进一步的，公式12中的能量函数最小化问题通过α-expansion求解。

本发明提出一种基于机载LiDAR点云利用L0梯度最小化和图割全局优化进行房屋顶面提取的算法，很好的实现了房屋顶面的高精度提取，并且对于不同的数据类型具备很强的鲁棒性。本发明首先对原始建筑物点云进行欧式聚类，以分割出每栋建筑物的点云；然后分别对每一栋建筑进行基于L0梯度最小化的多尺度分割，得到初始平面分割结果。针对房屋顶面提取工作的特殊性我们对L0梯度最小化算法进行了针对性改进，包括排序和区域扩张的策略的应用以及附加的几何约束条件，很好的改进了L0梯度最小化算法在平面分割中的精度和鲁棒性。最后采用基于图割的后处理优化算法，将初始结果进行精化，同时解决初始结果中可能存在的过度分割以及边缘锯齿状的问题以得到最终的房屋顶面提取结果。

附图说明

图1是本发明的整体算法流程图。

图2是本发明的多尺度分割算法流程图。

具体实施方式

下面通过实施例，并结合附图，对本发明的技术方案作进一步的说明。

本发明采用的技术方案是：首先对原始建筑物点云进行欧式聚类，以分割出每栋建筑物的点云；然后分别对每一栋建筑进行基于L0梯度最小化的多尺度分割以得到初始平面分割结果；最后采用基于图割的后处理优化算法将初始结果进行精化，同时解决初始结果中可能存在的过度分割以及边缘锯齿状的问题，得到最终的建筑物顶面提取结果。

上述的欧式聚类方法为点云实例分割中常用方法，通过欧式聚类可以得到每一栋建筑物的点云，这么做提高了算法的时间效率并且提升算法的稳定性。

上述L0梯度最小化算法是我们的多尺度分割的基础方法，这种方法最初被应用于信号处理中并且用来控制相邻信号之间非零梯度个数。本发明根据发明目的将L0梯度最小化算法延伸扩至平面分割算法中，提出了一种基于区域扩张的L0梯度最小化算法。

L0梯度最小化算法的能量函数定义为：

I表示输入信号，S则表示输出信号。

表示输出信号的梯度。||·||₀代表L0范数。参数λ控制L0范数的权重从而决定最终输出信号的平滑程度。

公式1扩展到点云中去并且离散化之后可以表述为：

式中，M代表输入点云的长度，I_i和S_i分别代表点云中第i个点p_i输入和输出的法向量。N_i代表p_i的邻域，可以为k个最邻近的邻域点或者是三维空间中的以R为半径的球内的邻域点，i和j分别表示邻域点的索引，参数λ为L0范数的权重。

由于直接求取公式2的最小值是一个NP困难问题，因此需要对这个问题进行近似化求解。

首先以一对相邻点p_i和p_j为例，它们对于整体能量函数F的贡献值f为：

为了使f最小，可以采取两种不同的策略：不融合p_i与p_j，此时S_i≠S_j或者融合p_i与p_j，此时S_i＝S_j,并且两个相邻点之间的L0范数被消除。

对于S_i≠S_j的情况，具体做法为：

对于S_i＝S_j的情况，具体做法为：

将两者情况相结合，我们可得到：

同时考虑到三维平面分割问题的特殊性，我们加入几何约束来进一步提升算法的稳定性。这里引入一个平面距离的概念，同样以一对相邻点p_i和p_j为例，它们之间的平面距离Pd_i,j可以定义如下：

Pd_i,j＝abs[p_ij·I_i] (6)

其中P_ij代表从p_i到p_j的向量，I_i代表p_i的法向量,abs[]为取绝对值。因此Pd_i,j就代表着两个点在法向量I_i方向上的投影距离，当这个投影距离过大时我们认为他们不属于同一平面，尽管他们可能拥有相似的法向量。因此当两个点之间的平面距离大于规定阈值T_pd，我们就不再惩罚它们之间的非零梯度。

因此公式6可以重写为：

接下来，把上述公式从相邻点扩展至相邻点集。以G_i代表一个点集，Y_i代表其法向量，点集中点数为w_i,点集G_i与G_j之间的连接点个数为c_i,j；结合上述几何约束条件相邻点集之间的能量公式为：

对于相邻点集，平面距离的计算也是通过公式7，其中P_ij代表着点集G_i与G_j的重心点连接所形成的向量。

上述公式解为：

公式(9)即为我们基于区域扩张的L0梯度最小化算法的扩张条件。下面我们介绍L0梯度最小化算法的整体流程。

1)计算待分割点云中每个点的局部平滑特征，按照局部平滑特征值对点的操作顺序进行排序，所述局部平滑特征为该点与其邻域中其他邻域点法向量角度差异值的最大值；

2)初始化点集并设定一个预设参数β，对于每一个点都看作一个只有一个点的点集；

3)引入辅助参数λ，并且每次迭代中λ的值增加1.2倍；辅助参数的初始化由下述公式完成：

其中median{}代表取中值。

4)在每一次迭代过程中，按顺序分别处理所有点集，对于点集G_i，根据扩张条件(公式10)来确定是否将点集G_i的所有邻域点集G_j加入到G_i，如果G_j加入G_i，相对应的w_i将会被更新，G_j的邻域将会被加入N_i，完成一次扩张之后，G_j相关信息将会被删除，待处理点集个数减少1个；

5)循环处理所有点集，直至辅助参数λ大于预设参数β，迭代终止，剩余的点集则对应着分割出的平面，并且每个点集中的点都有着一致的平面法向量。

算法流程如下：

由于单次L0梯度最小化算法的分割效果受参数λ影响很大，因此为了让算法具备自适应性以及更高的鲁棒性，本发明提出一种多尺度分割流程，首先引入平面指数概念，平面指数c用来衡量一个聚类是否为合法平面；对一个聚类，先利用主成分分析法(PCA)来拟合平面，然后计算聚类中每个点到拟合平面的距离d_i，平面指数则为所有点到拟合平面的平均距离；只有当一个聚类的平面指数大于预设阈值T_c时我们才会认为它是一个合法平面。

1)先进行欧式聚类，来获得不同房屋的点云聚类。

2)选取一个较大的分割尺度参数β为起始β，以欧式聚类得到的点云聚类为输入。

3)对于每一个聚类首先进行平面验证，如果为平面指数大于阈值，则直接输出为合法平面，否则即在当前分割尺度下通过基于区域扩张的L0梯度最小化算法进行平面分割以获得新的聚类。

4)当所有聚类完成步骤3的处理，如果产生新的聚类，则算法以新生成聚类为输入进行下一次迭代，并且分割尺度参数减小，参数β缩小1.5倍。如果没有产生新的聚类即所有输入聚类都为合法平面，算法终止。

多尺度算法对于起始λ具备鲁棒性，并且可以尽可能的检测到所有的潜在平面。为了精化多尺度分割得到地分割结果，我们采用一个基于图割的后处理流程来解决多尺度分割中存在的过度分割以及边缘锯齿状的问题。

基于图割的后处理流程主要是以多尺度分割结果为初始平面通过图割算法来优化目标函数以得到每个点到平面的最优对应关系。所优化的目标函数由三部分组成，分别为：数据项，平滑项以及标签项。定义如下：

其中L为平面，P为平面L上的所有点的点集，L_p即代表点p所属的平面，公式中p,q∈N代表相邻的两点；

其中数据项主要衡量点到平面的距离，计算公式为：

其中

为L_p的平面方程。

平滑项用来惩罚属于不同平面的相邻点的个数，δ()为指示函数，当L_p≠L_q，这个函数值为1，ω是这一项代价的权重。

最后一项是标签项，这一项代价主要惩罚最终所有点对应不同平面的个数，如果权重κ设置的比较小，那么最终的平面个数会更多。需要注意的是，这里最终产生的平面是输入平面的子集。

公式12中的能量函数最小化问题可以通过α-expansion求解。

后处理流程结束之后，我们最终得到精化的平面分割结果。

为了验证本发明在建筑物顶面提取中的作用，我们在两套数据集上对本发明方法进行评价，其中一套数据集来自于中国杭州市，另一套数据集来自德国Vaihingen。Vaihingen数据集具有三个测区，下面的表格详细介绍了两套数据集的具体情况。

表1两套数据集的详细情况说明

主要评价指标为：完整度(completeness，C_m)，准确度(correctness,C_r)，分割质量(quality,Q_l)，过度分割指数(over-segmentation rate,R_o)，分割不足指数(under-segmentation rate,R_u)，边缘精确度(boundary precision,B_p)，边缘回召(boundaryrecall,B_r)，边缘F指数(boundary F-measure,F_m)。这八个指标是平面分割中常用的评价指标。

其中完整度反映真值中的平面被正确提取的百分比，准确度反映提取平面中为正确提取平面的百分比，分割质量是一个综合指标，可以看作是主要指标。过度分割指数以及分割不足指数用来衡量提取结果的过度分割的程度以及分割不足的程度。边缘F指数使边缘精确度以及边缘回召的综合指数，反映提取平面的边缘精度。

实验结果见下表：

表2实验结果

实施例：

参见图1和图2，本发明通过多尺度L0梯度最小化算法以及图割优化算法从机载LiDAR点云中高效且高精度的提取建筑物顶面。

以ISPRS公开数据集Vaihingen数据集一号测区为例。

步骤一：欧式聚类

首先我们计算点云的分辨率，计算方式是统计所有点的与最临近点的平均距离。根据点云分辨率设置欧式聚类的最小距离阈值，经过实验证明三倍点云分辨率可以很好的得到欧式聚类分割效果。欧式聚类的具体流程如下：

根据最小距离阈值建立Kd-tree，通过Kd-tree可以确定每个点的邻域点，而邻域点距离中心点的距离均不大于上述距离阈值；然后随机从一个点开始，这一个点构成一个初始的聚类，根据Kd-tree将这个点的邻域点加入聚类中，然后对于每个新加入的点，再通过上述操作将他们的邻域点进一步加入聚类，重复上述操作直至没有新的点加入为止，此时就可以得到一个聚类，将形成聚类的点从原始点云中剔除，再从随机选点开始重复上述流程直至所有的点都被加入某个聚类中。

步骤二：多尺度分割

参见图2，以欧式聚类的聚类结果为输入，根据数据特点进行参数设置，经过实验整明多尺度分割流程对于起始λ不敏感，为了减少不必要的迭代，设置起始λ为0.1。然后确定平面验证时的阈值T_c，实验证明T_c选取为0.05可以得到比较好的结果。对于L0梯度最小化算法，设定平面距离阈值T_pd为0.1。确定所有的参数设置之后，我们按照图2所述流程进行多尺度分割以获得初始平面分割结果。

步骤三：基于图割后处理

首先确定图割优化的参数设定，分别是平滑项和标签项的权重。针对数据特点，我们设定平滑项的权重为0.05,标签项的权重为5。然后根据初始平面分割结果获得初始标签，每一个平面模型作为一个标签，一个平面标签包括重心以及平面法向量。然后利用α-expansion算法进行全局优化问题求解，得到所有的点对应的标签，根据标签可以确定每个面所对应的点集，最终得到房屋顶面的提取结果。

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。

Claims (6)

1.一种基于机载LiDAR点云分割建筑物顶面的方法，其特征在于，包括如下步骤：

步骤1，对原始建筑物点云进行欧式聚类，以分割出每栋建筑物的点云；

步骤2，分别对每一栋建筑物进行基于区域扩张的L0梯度最小化算法的多尺度分割以得到初始平面分割结果；

步骤2中所述的基于区域扩张的L0梯度最小化算法是对L0梯度最小化算法的改进，其具体原理如下；

基于区域扩张的L0梯度最小化算法的能量函数定义为：

式中，M代表输入点云的长度，I_i和S_i分别代表点云中第i个点p_i输入和输出的法向量；N_i代表p_i的邻域，为k个最邻近的邻域点或者是三维空间中的以R为半径的球内的邻域点，i和j表示邻域点的索引，参数λ为L0范数的权重；

对公式(2)进行近似化求解，首先以一对相邻点p_i和p_j为例，它们对于整体能量函数F的贡献值f为：

为了使f最小，采取两种不同的策略：不融合p_i与p_j，此时S_i≠S_j或者融合p_i与p_j，此时S_i＝S_j，并且两个相邻点之间的L0范数被消除；

对于S_i≠S_j的情况，具体做法为：

对于S_i＝S_j的情况，具体做法为：

将两者情况相结合得到：

然后加入几何约束来进一步提升算法的稳定性，首先引入一个平面距离的概念，以一对相邻点p_i和p_j为例，它们之间的平面距离Pd_i,j定义如下：

Pd_i,j＝abs[p_ij·I_i] (6)

其中p_ij代表从p_i到p_j的向量，I_i代表p_i的法向量，abs[]为取绝对值；因此Pd_i,j代表着两个点在法向量I_i方向上的投影距离，当两个点之间的平面距离大于规定阈值T_pd，就不再惩罚它们之间的非零梯度；

因此公式(6)重写为：

接下来，把上述公式从相邻点扩展至相邻点集，以G_i代表一个点集，Y_i代表其法向量，点集中点数为w_i，点集G_i与G_j之间的连接点个数为c_i,j，结合上述几何约束条件相邻点集之间的能量公式为：

上述公式解为：

公式(9)即为基于区域扩张的L0梯度最小化算法的扩张条件。

2.如权利要求1所述的一种基于机载LiDAR点云分割建筑物顶面的方法，其特征在于：步骤2的基于区域扩张的L0梯度最小化算法的流程如下，

1)计算待分割点云中每个点的局部平滑特征，按照局部平滑特征值对点的操作顺序进行排序；

2)初始化点集并设定一个预设参数β，对于每一个点都看作一个只有一个点的点集；

3)引入辅助参数λ’，并且每次迭代中λ’的值增加k倍；辅助参数的初始化由下述公式完成：

其中median{}代表取中值；

4)在每一次迭代过程中，按顺序分别处理所有点集，对于点集G_i，根据扩张条件来确定是否将点集G_i的所有邻域点集G_j加入到G_i，如果G_j加入G_i，相对应的w_i将会被更新，G_j的邻域将会被加入N_i，完成一次扩张之后，G_j相关信息将会被删除，待处理点集个数减少1个；

5)循环处理所有点集，直至辅助参数λ’大于预设参数β，迭代终止，剩余的点集则对应着分割出的平面，并且每个点集中的点都有着一致的平面法向量；

步骤3，采用基于图割的后处理优化算法将初始平面分割结果进行精化，提取高精度建筑物顶面。

3.如权利要求2所述的一种基于机载LiDAR点云分割建筑物顶面的方法，其特征在于：所述局部平滑特征为该点与其邻域中其他邻域点法向量角度差异值的最大值。

4.如权利要求2所述的一种基于机载LiDAR点云分割建筑物顶面的方法，其特征在于：步骤2中多尺度分割的流程为，

(1)选取一个起始的分割尺度参数β，以欧式聚类得到的点云聚类为输入，其中分割尺度参数β与预设参数相同；

(2)对于每一个聚类首先进行平面验证，如果平面指数大于阈值，则直接输出为合法平面，否则即在当前分割尺度下通过基于区域扩张的L0梯度最小化算法进行平面分割以获得新的聚类；

平面指数的定义为：对一个聚类，先利用主成分分析法PCA来拟合平面，然后计算聚类中每个点到拟合平面的距离d_i，平面指数则为所有点到拟合平面的平均距离；

(3)当所有聚类完成(2)的处理，如果产生新的聚类，则以新生成聚类为输入进行下一次迭代，并且分割尺度参数减小；如果没有产生新的聚类即所有输入聚类都为合法平面，算法终止。

5.如权利要求1所述的一种基于机载LiDAR点云分割建筑物顶面的方法，其特征在于：步骤3中基于图割的后处理流程是以多尺度分割结果为初始平面通过图割算法来优化目标函数以得到每个点到平面的最优对应关系，所优化的目标函数由三部分组成，分别为：数据项，平滑项以及标签项，定义如下：

其中L为平面，P为平面L上的所有点的点集，L_p即代表点p所属的平面，公式中p,q∈N代表相邻的两点；

其中数据项是衡量点到平面的距离，计算公式为：

其中

为L_p的平面方程；

平滑项用来惩罚属于不同平面的相邻点的个数，δ()为指示函数，当L_p≠L_q，这个函数值为1，ω是权重；

最后一项是标签项，这一项用于惩罚最终所有点对应不同平面的个数，κ为权重。

6.如权利要求1所述的一种基于机载LiDAR点云分割建筑物顶面的方法，其特征在于：公式(12)中的能量函数最小化问题通过α-expansion求解。

Images