CN116883754A - 一种面向地面LiDAR点云的建筑物信息提取方法 - Google Patents

Abstract

本发明公开了一种面向地面LiDAR点云的建筑物信息提取方法，包括：S1、利用体素滤波对原始点云进行降采样；S2、进行点云裁剪并综合利用统计、半径滤波进行去噪；S3、使用布料模拟滤波初步分离地面点与非地面点；S4、引入法线微分算子对初步获取的地面点过滤，得到最终地面点与非地面点；S5、针对非地面点构建Kd‑tree，并按高程值均匀分层，对每层点云进行密度聚类；S6、根据投影密度分隔聚类结果中建筑物与粘连的非建筑物；S7、对分隔开的点云进行区域增长；S8、通过点云总个数、高差、分层投影面积比和平面拟合残差四类特征对建筑物点云簇进行识别分类。本发明有效提高了传统建筑物提取方法的精度，为从地面LiDAR点云中提取建筑物提供了新的处理思路。

Description

一种面向地面LiDAR点云的建筑物信息提取方法

技术领域

本发明涉及激光扫描数据处理与信息提取领域，具体涉及一种面向地面LiDAR点云的建筑物信息提取方法。

背景技术

城市建设项目竣工验收阶段“多测合一”作业机制正在全国范围内推广应用，该机制使得竣工测绘时限缩短，因此对统一测绘的作业效率及成果质量提出更高要求。地面激光扫描技术(Light Detection and Ranging,LiDAR)以其高精度、高效率、高分辨率、全要素采集等特点，能够满足“多测合一”技术和时限要求，近年来在城市测量领域中得到广泛关注。但目前内业点云数据的处理、信息提取仍主要依靠技术人员与数据的交互，手工利用相应的软件进行，重复操作多，内业处理工作量和复杂度大。建筑物作为城市竣工测绘主要关注的测绘要素之一，如何从海量的点云中准确提取出建筑物点云信息一直是相关研究的热点、重点。

目前从LiDAR点云中提取建筑物的研究已有较多的成果，但在实际的建筑物点云提取中，需要结合不同的点云数据特征、场景特性选择合适的方法，传统的建筑物提取方法大都是针对机载LiDAR点云设计的，未顾及到地面LiDAR与机载LiDAR点云间的差异，导致利用传统方法获取到的地面LiDAR建筑物点云底部与其他地物粘连及高层建筑物顶部不完整的现象。因此，亟需新的一种面向地面LiDAR点云的建筑物信息提取方法以解决上述问题。

发明内容

针对现有技术中的上述不足，本发明提供了一种面向地面LiDAR点云的建筑物信息提取方法。

为了达到上述发明目的，本发明采用的技术方案为：

一种面向地面LiDAR点云的建筑物信息提取方法，包括以下步骤：

S1、利用体素滤波对原始点云进行降采样；

S2、进行点云裁剪并综合利用统计、半径滤波进行去噪；

S3、使用布料模拟滤波初步分离地面点与非地面点；

S4、引入法线微分算子对初步获取的地面点过滤，得到最终地面点与非地面点；

S5、针对非地面点构建Kd-tree，并按高程值均匀分层，对每层点云进行密度聚类；

S6、根据投影密度分隔聚类结果中建筑物与粘连的非建筑物；

S7、对分隔开的点云进行区域增长；

S8、通过点云总个数、高差、分层投影面积比和平面拟合残差四类特征对建筑物点云簇进行识别分类。

进一步地，步骤S1中利用体素滤波对原始点云进行降采样。其具体方法为：

S11、将原始地面LiDAR点云数据按照一定空间大小三维体素栅格进行划分；

S12、计算三维体素栅格的重心，具体公式如下所示：

式中，(x_h,y_h,z_h)为三维体素栅格的重心坐标，(x_i,y_i,z_i)为三维体素栅格中每个点对应的坐标，n为该三维体素栅格中所包含的点个数；

S13、在每个体素栅格内，用体素栅格的重心代表该体素栅格中所有点，过滤掉除重心点以外的剩余点，得到降采样后的点云数据。

进一步地，步骤S2进行点云裁剪并综合利用统计、半径滤波进行去噪。其具体方法为：

S21、在X、Y、Z轴上设置可接受的裁剪范围，对范围外的点云进行过滤，实现对非目标区域点云的裁剪；

S22、假设点集中点的总个数为N，首先计算点集中每个点p到其k个邻近点的平均欧式距离d_i，据此分别计算出点集中的全局平均距离和标准偏差σ：

式中(x_p，y_p,z_p)为每个点p的三维坐标，(x_j,y_j,z_j)为点p邻近点的三维坐标；

S23、通过自定义设定标准差倍数λ(一般取3)求取一个合适的距离阈值，所有大于等于这个阈值的点可以被视为偏离模型的点即离群点，并将其从数据集中剔除，以下二式为距离值的定义和离群点的判定方法：

式中t为设定的距离阈值，outlier表示离群点；

S24、计算点集中每个点在半径r范围内的临近点数k_i，如果其少于最小邻近点数k_min(与点云密度相关，其范围一般在20-100)，那么就认为该点的邻近点相当少，将其从点集中去除。

进一步地，步骤S3中使用布料模拟滤波初步分离地面点与非地面点。其具体方法为：

S31、转换点云几何坐标，将经降采样和去噪后的点云倒置；

S32、初始化布料网格，通过格网分辨率确定网格质点数量，记录格网对应的激光脚点高程HIV；

S33、对于每个可以移动的格网质点，计算每个粒子仅受到重力作用下的位移，其位移具体计算公式如下所示：

其中△t为时间步长；G为重力常数；X(t+△t)为下一步粒子所在位置；X(t)为当前粒子所处位置；X(t-△t)为上一步粒子所处位置。比较粒子当前所处位置高程与其对应激光脚点高程HIV，如果节点高程小于或等于激光脚点高程HIV，则将该节点位置替换到对应的激光脚点位置，并将其标记为不可动点；

S34、计算每个粒子仅受到内力作用下的位移，若2个相邻粒子都是可移动点，则其中一个向上移动一个向下移动；若2个相邻点中一个为可移动点，另一个为不可移动点，则移动可移动点，移动距离与布料硬度相关，布料硬度越大，移动位移越小，具体位移量可以通过下式计算：

式中，为粒子位移量，当粒子可移动时，b＝1，不可移动时，b＝0，/>分别表示粒子与其近邻粒子的当前位置，/>为点标准化到垂直方向上的单位向量；

S35、重复步骤S33和S34，当各个质点的位移变化开始小于阈值或者迭代次数达到最大次数的时候，模拟过程终止；

S36、对点云位置进行判断，如果该点的激光脚点和布料粒子下落最终位置间的距离小于预设的高度阈值，则将其归类为地面点，否则归类为非地面点。

进一步地，步骤S4中引入法线微分算子对初步获取的地面点过滤，得到最终非地面点。其具体方法为：

S41、设定法线估计参数中的较小支撑半径r₁和较大支撑半径r₂

S42、S36所得的初步地面点集为P＝{P₁,P₂,...,P_n}，n为初步获取的地面点个数，对于其中某点P_i，通过r邻域半径区域搜索到邻域集合N(P_i)，邻域集合点个数为k，计算该邻域集合的重心

S43、计算每一个P_i点的协方差矩阵C，其具体计算公式如下所示：

式中，λ_j是协方差矩阵的第j个特征值，是第j个特征向量。对矩阵C进行求解，得到最小特征值对应的特征向量即为点P_i在邻域半径r下的法向量/>

S44、计算不同尺度下的法线微分算子具体计算公式如下：

S45、设定向量域的阈值，过滤计算所得的向量域，剔除初步地面点中错分的非地面点，得到最终的地面点和非地面点。

进一步地，步骤S5中针对非地面点构建Kd-tree，并按高程值均匀分层，对每层点云进行密度聚类。其具体方法为：

S51、对S45获取的非地面点构建Kd-tree，并按照高程值均匀分层，分层结果如下式：

Range＝{(Z_min,Z₁)(Z₁,Z₂),……,(Z_f-2,Z_f-1)(Z_f-1,Z_max)}

式中，f为分段数，Z_min为非地面点高程最小值，Z_max为高程最大值。根据各层点云数量将其分为密集分层区和稀疏分层区，密集分层区一般为靠近地面的1-2层。

S52、分别设定不同点云分层区域密度聚类的两个参数，邻域半径Eps和密度阈值MinPts，当为密集分层区时，Eps值尽量设小一些(0.1-1.5m)，当为稀疏分层区时，Eps值尽量设大一些(2-10m)。MinPts与点云密度相关，其值一般设定在20到200间；

S53、从点云中选取一个尚未被检测的点P_j，根据给定的阈值Eps和MinPts判断P_j是否为核心点，若不为核心点，则标记P_j已检测并重新选取未被检测的点；若P_j为核心点，则创建新簇Q，并将点P_j和Eps邻域内的点加入簇Q中，同时标记P_j为已检测；

S54、遍历Q中所有未被检测的点P′,若P′在给定阈值Eps和MinPts下仍为核心点，则将点P′的Eps邻域内中所有点加入簇Q，并标记P′为已检测；

S55、重复步骤S54，直至没有新的点加入簇Q；

S56、重复步骤S53至S55，直至点云中的所有点都被检测。

进一步地，步骤S6中根据投影密度分隔聚类结果中建筑物与粘连的非建筑物。其具体方法为：

S61、将点云投影至XOY平面内并构建规则格网，设定格网尺寸Grid，将试验区划分为N×M个矩形格网，并对其编号，为避免在格网边界计算，扩大格网边界，格网数量计算公式具体如下：

其中x_max、x_min、y_max、y_min分别表示点云坐标在X、Y轴方向上的最值，INT为取整；

S62、统计每个格网内点数，即点云投影密度值，通过点云投影密度值分布情况判断地物粘连现象，设置合理点数阈值(一般在30-200间)，若格网内点数大于此阈值则保留此格网；

S63、将满足要求的格网全部取出来，即可将与建筑物粘连着的非建筑物点云分隔开。

进一步地，步骤S7中对分隔开的点云进行区域增长。其具体方法为：

S71、计算当前区域内点云数据的曲率大小，根据点的曲率对点进行排序，然后将曲率最小点设置为初始种子点并加入种子序列W，设置空聚类簇U；

S72、对当前种子点的邻近点进行搜索，计算其邻近点的法线，若当前种子点的法线与其邻近点的法线之间的夹角θ小于角度阈值，则将这些邻近点添至聚类簇U中；

S73、逐个检测这些邻域点的曲率值，如果曲率值小于曲率阈值，则将该点加入到种子序列W中，再把当前点删除，通过迭代，一直到种子序列清空为止；

S74、重复S72和S73，直到遍历完所有的数据点，区域增长结束。

进一步地，步骤S8中通过点云总个数、高差、分层投影面积比和平面拟合残差四类特征对建筑物点云簇进行识别分类。其具体方法为：

S81、统计密度聚类及区域增长获得的每个点云簇中的点云总个数。通过此参数可以快速的去除行人、垃圾、离散点等小体量点云簇；

S82、计算剩余点云簇中最高点的高程与最低点高程之差，根据高差与设置的阈值相比较，区分出低矮地物(通常认为此类地物高差小于1.5m)；

S83、计算剩余点云簇的水平投影面积，进而得到不同地物的分层投影面积比，根据分层投影面积比的差异(分层投影面积比远大于1)将树木、路灯等地物点云簇分离；

S84、计算剩余点云簇的平面拟合残差，将有多个面的对象和仅有单个面的对象区分开(如围墙、广告牌一般属于单一平面，其平面拟合残差通常较小，建筑物则属于多平面，其平面拟合残差通常较大)，最终得到建筑物点云。

进一步地，步骤S83具体为：

S831、在平面投影点集中选取Y坐标最小的一点G₀当作基点，如果存在多个点的Y坐标都为最小值，则选取X坐标最小的一点；

S832、对剩余点按G₀点的逆时针方向进行极角排序，当极角相同时，距离G₀比较近的点排在前面；

S833、根据“凸多边形的各顶点必在该多边形任意一条边同一侧”的原则对投影点集相继三个顶点连线进行一次有序扫描，删去投影边界序列中的非凸包顶点，最终剩余的点即为凸包的顶点；

S834、当点云数据被分为F层时，自下而上，第f层的投影面积记为S_f(f＝1,2,...,F)，设第f层的凸包多边形有L个顶点，则投影面积及分层投影面积比△S_F计算公式如下：

式中，x_i、y_i表示第i个顶点的坐标，x_i+1、y_i+1表示第i+1个顶点的坐标，max S_f为点云数据分层后最大的投影面积，min S_f为最小的投影面积。

进一步地，步骤S84具体为：

S841、将点云簇中所有点作为最小二乘平面拟合得到的平面方程ax+by+cz+d＝0的四个参数a、b、c、d；

S842、计算任一点D_v到此平面的欧式距离，进而得到平面拟合残差δ，具体计算公式如下：

式中，n为点云簇中点云总数，(x_v,y_v,z_v)为点D_v的三维空间坐标。

本发明具有以下有益效果：

本发明针对地面LiDAR点云数据，引入法线微分算子对布料模拟滤波提取到的地面点进行二次滤除，准确地将地面点和非地面点进行了分离。针对获取的非地面点，采用分层密度聚类的方法有效保证了提取结果高层建筑物顶部的完整性。针对密度聚类提取到的建筑物点云底部与其他地物粘连的现象，使用点云投影密度特征将粘连着的非目标点云分隔开，再利用区域增长算法实现建筑物点云的精提取，显著地改善了建筑物点云底部与其他地物粘连的现象。最后通过点云总个数、高差、分层投影面积比和平面拟合残差四类特征对建筑物点云簇进行识别分类，得到了完整、准确的建筑物点云，为从地面LiDAR点云中提取建筑物提供了相应的处理思路。

附图说明

图1所示为本发明提供的面向地面LiDAR点云建筑物提取方法流程图。

图2所示为本发明实施例提供的经降采样及裁剪去噪后的结果图。

图3(a)所示为本发明实施例提供的经地面滤波获得的地面点结果图。

图3(b)所示为本发明实施例提供的经地面滤波获得的非地面点结果图。

图4(a)所示为LiDAR 360软件提取建筑物点云结果图。

图4(b)所示为密度聚类算法提取建筑物点云结果图。

图4(c)所示为区域增长算法提取建筑物点云结果图。

图4(d)所示为本发明实施例提供的建筑物点云提取结果图。

具体实施方式

下面对本发明的具体实施方式进行描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。

本发明实施例提供了一种面向地面LiDAR点云建筑物提取方法，如图1所示，包括以下步骤S1～S8：

S1、利用体素滤波对原始点云进行降采样。

本发明实施例中，利用体素滤波对原始点云进行降采样，具体方法如下：

S11、将原始地面LiDAR点云数据按照一定空间大小三维体素栅格进行划分；

S12、计算三维体素栅格的重心，具体公式如下所示：

式中，(x_h,y_h,z_h)为三维体素栅格的重心坐标，(x_i,y_i,z_i)为三维体素栅格中每个点对应的坐标，n为该三维体素栅格中所包含的点个数；

S13、在每个体素栅格内，用体素栅格的重心代表该体素栅格中所有点，过滤掉除重心点以外的剩余点，得到降采样后的点云数据。

S2、进行点云裁剪并综合利用统计、半径滤波进行去噪。

步骤S2包括以下分步骤S21～S24：

S21、在X、Y、Z轴上设置可接受的裁剪范围，对范围外的点云进行过滤，实现对非目标区域点云的裁剪；

S22、假设点集中点的总个数为N，首先计算点集中每个点p到其k个邻近点的平均欧式距离d_i，据此分别计算出点集中的全局平均距离和标准偏差σ：

式中(x_p,y_p,z_p)为每个点p的三维坐标，(x_j,y_j,z_j)为点p邻近点的三维坐标；

S23、通过自定义设定标准差倍数λ(一般取3)求取一个合适的距离阈值，所有大于等于这个阈值的点可以被视为偏离模型的点即离群点，并将其从数据集中剔除，以下二式为距离值的定义和离群点的判定方法：

式中t为设定的距离阈值，outlier表示离群点；

S24、计算点集中每个点在半径r范围内的临近点数k_i，如果其少于最小邻近点数k_min(与点云密度相关，其范围一般在20-100)，那么就认为该点的邻近点相当少，将其从点集中去除。经降采样及裁剪去噪后的结果如图2所示。

S3、使用布料模拟滤波初步分离地面点与非地面点。

步骤S3包括以下分步骤S31-S36:

S31、转换点云几何坐标，将经降采样和去噪后的点云倒置；

S32、初始化布料网格，通过格网分辨率确定网格质点数量，记录格网对应的激光脚点高程HIV；

S33、对于每个可以移动的格网质点，计算每个粒子仅受到重力作用下的位移，其位移具体计算公式如下所示：

其中△t为时间步长；G为重力常数；X(t+△t)为下一步粒子所在位置；X(t)为当前粒子所处位置；X(t-△t)为上一步粒子所处位置。比较粒子当前所处位置高程与其对应激光脚点高程HIV，如果节点高程小于或等于激光脚点高程HIV，则将该节点位置替换到对应的激光脚点位置，并将其标记为不可动点；

S34、计算每个粒子仅受到内力作用下的位移，若2个相邻粒子都是可移动点，则其中一个向上移动一个向下移动；若2个相邻点中一个为可移动点，另一个为不可移动点，则移动可移动点，移动距离与布料硬度相关，布料硬度越大，移动位移越小，具体位移量可以通过下式计算：

式中，为粒子位移量，当粒子可移动时，b＝1，不可移动时，b＝0，/>分别表示粒子与其近邻粒子的当前位置，/>为点标准化到垂直方向上的单位向量；

S35、重复步骤S33和S34，当各个质点的位移变化开始小于阈值或者迭代次数达到最大次数的时候，模拟过程终止；

S36、对点云位置进行判断，如果该点的激光脚点和布料粒子下落最终位置间的距离小于预设的高度阈值，则将其归类为地面点，否则归类为非地面点。

S4、引入法线微分算子对初步获取的地面点过滤，得到最终非地面点。

步骤S4包括以下分步骤S41～S45：

S41、设定法线估计参数中的较小支撑半径r₁和较大支撑半径r₂

S42、S36所得的初步地面点集为P＝{P₁,P₂,...,P_n}，n为初步获取的地面点个数，对于其中某点P_i，通过r邻域半径区域搜索到邻域集合N(P_i)，邻域集合点个数为k，计算该邻域集合的重心

S43、计算每一个P_i点的协方差矩阵C，其具体计算公式如下所示：

式中，λ_j是协方差矩阵的第j个特征值，是第j个特征向量。对矩阵C进行求解，得到最小特征值对应的特征向量即为点P_i在邻域半径r下的法向量/>

S44、计算不同尺度下的法线微分算子具体计算公式如下：

S45、设定向量域的阈值，过滤计算所得的向量域，剔除初步地面点中错分的非地面点，得到最终地面点与非地面点，结果如图3(a)和图3(b)所示。

S5、针对非地面点构建Kd-tree，并按高程值均匀分层，对每层点云进行密度聚类。

步骤S5包括以下分步骤S51～S56：

S51、对S45获取的非地面点构建Kd-tree，并按照高程值均匀分层，分层结果如下式：

Range＝{(Z_min,Z₁)(Z₁,Z₂),……,(Z_f-2,Z_f-1)(Z_f-1,Z_max)}

式中，f为分段数，Z_min为非地面点高程最小值，Z_max为高程最大值。根据各层点云数量将其分为密集分层区和稀疏分层区，密集分层区一般为靠近地面的1-2层。

S52、分别设定不同点云分层区域密度聚类的两个参数，邻域半径Eps和密度阈值MinPts，当为密集分层区时，Eps值尽量设小一些(0.1-1.5m)，当为稀疏分层区时，Eps值尽量设大一些(2-10m)。MinPts与点云密度相关，其值一般设定在20到200间；

S53、从点云中选取一个尚未被检测的点P_j，根据给定的阈值Eps和MinPts判断P_j是否为核心点，若不为核心点，则标记P_j已检测并重新选取未被检测的点；若P_j为核心点，则创建新簇Q，并将点P_j和Eps邻域内的点加入簇Q中，同时标记P_j为已检测；

S54、遍历Q中所有未被检测的点P′,若P′在给定阈值Eps和MinPts下仍为核心点，则将点P′的Eps邻域内中所有点加入簇Q，并标记P′为已检测；

S55、重复步骤S54，直至没有新的点加入簇Q；

S56、重复步骤S53至S55，直至点云中的所有点都被检测。

S6、根据投影密度分隔聚类结果中建筑物与粘连的非建筑物。

步骤S6包括以下分步骤S61～S63：

S61、将点云投影至XOY平面内并构建规则格网，设定格网尺寸Grid，将试验区划分为N×M个矩形格网，并对其编号，为避免在格网边界计算，扩大格网边界，格网数量计算公式具体如下：

其中x_max、x_min、y_max、y_min分别表示点云坐标在X、Y轴方向上的最值，INT为取整；

S62、统计每个格网内点数，即点云投影密度值，通过点云投影密度值分布情况判断地物粘连现象，设置合理点数阈值(一般在30-200间)，若格网内点数大于此阈值则保留此格网；

S63、将满足要求的格网全部取出来，即可将与建筑物粘连着的非建筑物点云分隔开。

S7、对分隔开的点云进行区域增长。

步骤S7包括以下分步骤S71～S74：

S71、计算当前区域内点云数据的曲率大小，根据点的曲率对点进行排序，然后将曲率最小点设置为初始种子点并加入种子序列W，设置空聚类簇U；

S72、对当前种子点的邻近点进行搜索，计算其邻近点的法线，若当前种子点的法线与其邻近点的法线之间的夹角θ小于角度阈值，则将这些邻近点添至聚类簇U中；

S73、逐个检测这些邻域点的曲率值，如果曲率值小于曲率阈值，则将该点加入到种子序列W中，再把当前点删除，通过迭代，一直到种子序列清空为止；

S74、重复S72和S73，直到遍历完所有的数据点，区域增长结束。

S8、通过点云总个数、高差、分层投影面积比和平面拟合残差四类特征对建筑物点云簇进行识别分类。

步骤S8包括以下分步骤S81～S84：

S81、统计密度聚类及区域增长获得的每个点云簇中的点云总个数。通过此参数可以快速的去除行人、垃圾、离散点等小体量点云簇；

S82、计算剩余点云簇中最高点的高程与最低点高程之差，根据高差与设置的阈值相比较，区分出低矮地物(通常认为此类地物高差小于1.5m)；

S83、计算剩余点云簇的水平投影面积，进而得到不同地物的分层投影面积比，根据分层投影面积比的差异(分层投影面积比远大于1)将树木、路灯等地物点云簇分离；

S84、计算剩余点云簇的平面拟合残差，将有多个面的对象和仅有单个面的对象区分开(如围墙、广告牌一般属于单一平面，其平面拟合残差通常较小，建筑物则属于多平面，其平面拟合残差通常较大)，最终得到建筑物点云，结果如图4(d)所示。

进一步地，步骤S83包括以下分步骤S831～S834：

S831、在平面投影点集中选取Y坐标最小的一点G₀当作基点，如果存在多个点的Y坐标都为最小值，则选取X坐标最小的一点；

S832、对剩余点按G₀点的逆时针方向进行极角排序，当极角相同时，距离G₀比较近的点排在前面；

S833、根据“凸多边形的各顶点必在该多边形任意一条边同一侧”的原则对投影点集相继三个顶点连线进行一次有序扫描，删去投影边界序列中的非凸包顶点，最终剩余的点即为凸包的顶点；

S834、当点云数据被分为F层时，自下而上，第f层的投影面积记为S_f(f＝1,2,...,F)，设第f层的凸包多边形有L个顶点，则投影面积及分层投影面积比△S_F计算公式如下：

/>

式中，x_i、y_i表示第i个顶点的坐标，x_i+1、y_i+1表示第i+1个顶点的坐标，max S_f为点云数据分层后最大的投影面积，min S_f为最小的投影面积。

进一步地，步骤S84包括以下分步骤S841～S842：

S841、将点云簇中所有点作为最小二乘平面拟合得到的平面方程ax+by+cz+d＝0的四个参数a、b、c、d；

S842、计算任一点D_v到此平面欧式距离，进而得到平面拟合残差δ，具体计算公式如下：

式中，n为点云簇中点云总数，(x_v,y_v,z_v)为点D_v的三维空间坐标。

为进一步说明本发明的有效性，将本发明提取到的建筑物点云与LiDAR 360软件、密度聚类算法和区域增长算法提取到的建筑物点云进行对比，结果如图4所示。由图4(a)利用LiDAR 360软件提取建筑物点云效果很差，存在较多错分现象，如将树木、花坛、围墙、台阶等错分为建筑物，同时，针对高层建筑顶部的破碎立面，提取得到的建筑物点云不完整现象严重；由图4(b)和(c)可知密度聚类算法和区域增长算法提取建筑物点云效果较差，但效果比软件自动提取的要好些，两类算法提取的建筑物点云在底部均出现了不同程度的粘连现象，且两类算法提取的建筑物点云顶部均出现了一定的残缺；由图4(d)可知，本发明提取建筑物点云效果较好，提取到的建筑物点云底部均未出现粘连现象，且针对高层建筑物，提取到的建筑物点云顶部也较为完整。

本发明是参照根据本发明实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。

这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。

这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。

本发明中应用了具体实施例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处，综上所述，本说明书内容不应理解为对本发明的限制。

本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通技术人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体变形和组合，这些变形和组合仍然在本发明的保护范围内。

Claims (10)

1.一种面向地面LiDAR点云的建筑物信息提取方法，其特征在于，包括如下步骤：

S1、利用体素滤波对原始点云进行降采样；

S2、进行点云裁剪并综合利用统计、半径滤波进行去噪；

S3、使用布料模拟滤波初步分离地面点与非地面点；

S4、引入法线微分算子对初步获取的地面点过滤，得到最终地面点与非地面点；

S5、针对非地面点构建Kd-tree，并按高程值均匀分层，对每层点云进行密度聚类；

S6、根据投影密度分隔聚类结果中建筑物与粘连的非建筑物；

S7、对分隔开的点云进行区域增长；

S8、通过点云总个数、高差、分层投影面积比和平面拟合残差四类特征对建筑物点云簇进行识别分类。

2.根据权利要求1所述的面向地面LiDAR点云的建筑物信息提取方法，其特征在于，所述步骤S1中使用体素滤波对原始点云进行降采样，其具体方法为：

S11、将原始地面LiDAR点云数据按照一定空间大小三维体素栅格进行划分；

S12、计算三维体素栅格的重心，具体公式如下所示：

式中，(x_h,y_h,z_h)为三维体素栅格的重心坐标，(x_i,y_i,z_i)为三维体素栅格中每个点对应的坐标，n为该三维体素栅格中所包含的点个数；

S13、在每个体素栅格内，用体素栅格的重心代表该体素栅格中所有点，过滤掉除重心点以外的剩余点，得到降采样后的点云数据。

3.根据权利要求1所述的面向地面LiDAR点云的建筑物信息提取方法，其特征在于，所述步骤S2中进行点云裁剪并综合利用统计、半径滤波进行去噪，其具体方法为：

S21、在X、Y、Z轴上设置裁剪范围，对非目标区域点云的裁剪；

S22、假设点集中点的总个数为N，计算点集中每个点p到其k个邻近点的平均欧式距离d_i，计算方式为：

根据所得到的平均欧氏距离d_i分别计算出点集中的全局平均距离和标准偏差σ，具体计算方式为：

式中(x_p,y_p,z_p)为每个点p的三维坐标，(x_j,y_j,z_j)为点p邻近点的三维坐标；

S23、通过自定义设定标准差倍数λ求取距离阈值，所有大于等于该距离阈值的点被视为离群点，并将其从数据集中剔除，以下二式为距离值的定义和离群点的判定方法：

式中t为设定的距离阈值，outlier表示离群点；

S24、计算点集中每个点在半径r范围内的临近点数k_i，若少于最小邻近点数k_min，则将该点从点集中去除。

4.根据权利要求1所述的面向地面LiDAR点云的建筑物信息提取方法，其特征在于，所述步骤S3中使用布料模拟滤波初步分离地面点与非地面点，其具体方法为：

S31、转换点云几何坐标，将经降采样和去噪后的点云倒置；

S32、初始化布料网格，通过格网分辨率确定网格质点数量，记录格网对应的激光脚点高程HIV；

S33、对于每个可以移动的格网质点，计算每个粒子仅受到重力作用下的位移，其位移具体计算公式如下所示：

其中△t为时间步长；G为重力常数；X(t+△t)为下一步粒子所在位置；X(t)为当前粒子所处位置；X(t-△t)为上一步粒子所处位置，

比较粒子当前所处位置高程与其对应激光脚点高程HIV，若节点高程小于或等于激光脚点高程HIV，则将该节点位置替换到对应的激光脚点位置，并将其标记为不可动点；

S34、计算每个粒子仅受到内力作用下的位移，若2个相邻粒子都是可移动点，则其中一个向上移动，另一个向下移动；若2个相邻点中一个为可移动点，另一个为不可移动点，则移动可移动点，具体位移量可以通过下式计算：

式中，为粒子位移量，当粒子可移动时，b＝1，不可移动时，b＝0，/>分别表示粒子与其近邻粒子的当前位置，/>为点标准化到垂直方向上的单位向量；

S35、重复步骤S33和S34，当各个质点的位移变化开始小于阈值或者迭代次数达到最大次数的时候，模拟过程终止；

S36、对点云位置进行判断，如果该点的激光脚点和布料粒子下落最终位置间的距离小于预设的高度阈值，则将其归类为地面点，否则归类为非地面点。

5.根据权利要求1所述的面向地面LiDAR点云的建筑物信息提取方法，其特征在于，所述步骤S4中引入法线微分算子对初步获取的地面点过滤，得到最终非地面点，其具体方法为：

S41、设定法线估计参数中的较小支撑半径r₁和较大支撑半径r₂；

S42、根据地面点集P＝{P₁,P₂,...,P_n}，n为地面点个数，

对于其中任一点P_i，通过r邻域半径区域搜索到邻域集合N(P_i)，邻域集合点个数为k，

计算该邻域集合的重心计算方式为：

S43、对矩阵C进行求解，计算每一个P_i点的协方差矩阵C，得到最小特征值对应的特征向量即为点P_i在邻域半径r下的法向量计算方式为：

式中，λ_j是协方差矩阵的第j个特征值，是第j个特征向量；

S44、计算不同尺度下的法线微分算子具体计算公式如下：

S45、设定向量域的阈值，过滤计算所得的向量域，剔除初步地面点中错分的非地面点，得到最终的地面点和非地面点。

6.根据权利要求1所述的面向地面LiDAR点云的建筑物信息提取方法，其特征在于，所述步骤S5中针对非地面点构建Kd-tree，并按高程值均匀分层，对每层点云进行密度聚类，其具体方法为：

S51、对S45获取的非地面点构建Kd-tree，并按照高程值均匀分层，分层结果如下式：

Range＝{(Z_min，Z₁)(Z₁，Z₂)，……，(Z_f-2，Z_f-1)(Z_f-1，Z_max)}

式中，f为分段数，Z_min为非地面点高程最小值，Z_max为高程最大值，根据各层点云数量将其分为密集分层区和稀疏分层区，密集分层区一般为靠近地面的1-2层；

S52、分别设定不同点云分层区域密度聚类的两个参数，邻域半径Eps和密度阈值MinPts；

S53、从点云中选取一个尚未被检测的点P_j，根据给定的阈值Eps和MinPts判断P_j是否为核心点，若不为核心点，则标记P_j已检测并重新选取未被检测的点；若P_j为核心点，则创建新簇Q，并将点P_j和Eps邻域内的点加入簇Q中，同时标记P_j为已检测；

S54、遍历Q中所有未被检测的点P′，若P′在给定阈值Eps和MinPts下仍为核心点，则将点P′的Eps邻域内中所有点加入簇Q，并标记P′为已检测；

S55、重复步骤S54，直至没有新的点加入簇Q；

S56、重复步骤S53至S55，直至点云中的所有点都被检测。

7.根据权利要求1所述的面向地面LiDAR点云的建筑物信息提取方法，其特征在于，根据权利要求1所述的面向地面LiDAR点云的建筑物信息提取方法，其特征在于，所述步骤S6中根据投影密度分隔聚类结果中建筑物与粘连的非建筑物，其具体方法为：

S61、将点云投影至XOY平面内并构建规则格网，设定格网尺寸Grid，将试验区划分为N×M个矩形格网，并对其编号，为避免在格网边界计算，扩大格网边界，格网数量计算公式具体如下：

其中x_max、x_min、y_max、y_min分别表示点云坐标在X、Y轴方向上的最值，INT为取整；

S62、统计每个格网内点数，即点云投影密度值，通过点云投影密度值分布情况判断地物粘连现象，设置合理点数阈值，若格网内点数大于此阈值则保留此格网；

S63、将满足要求的格网全部取出来，将与建筑物粘连着的非建筑物点云分隔开。

8.根据权利要求1所述的面向地面LiDAR点云的建筑物信息提取方法，其特征在于，所述步骤S7中对分隔开的点云进行区域增长，其具体方法为：

S71、计算当前区域内点云数据的曲率大小，根据点的曲率对点进行排序，然后将曲率最小点设置为初始种子点并加入种子序列W，设置空聚类簇U；

S72、对当前种子点的邻近点进行搜索，计算其邻近点的法线，若当前种子点的法线与其邻近点的法线之间的夹角θ小于角度阈值，则将这些邻近点添至聚类簇U中；

S73、逐个检测这些邻域点的曲率值，如果曲率值小于曲率阈值，则将该点加入到种子序列W中，再把当前点删除，通过迭代，一直到种子序列清空为止；

S74、重复S72和S73，直到遍历完所有的数据点，区域增长结束。

9.根据权利要求1所述的面向地面LiDAR点云的建筑物信息提取方法，其特征在于，所述步骤S8中通过点云总个数、高差、分层投影面积比和平面拟合残差四类特征对建筑物点云簇进行识别分类，其具体方法为：

S81、统计密度聚类及区域增长获得的每个点云簇中的点云总个数，通过此参数可以快速的去除行人、垃圾、离散点等小体量点云簇；

S82、计算剩余点云簇中最高点的高程与最低点高程之差，根据高差与设置的阈值相比较，区分出低矮地物；

S83、计算剩余点云簇的水平投影面积，进而得到不同地物的分层投影面积比，根据分层投影面积比的差异将树木、路灯等地物点云簇分离；

S84、计算剩余点云簇的平面拟合残差，将有多个面的对象和仅有单个面的对象区分开，最终得到建筑物点云，具体为：

S841、将点云簇中所有点作为最小二乘平面拟合得到的平面方程ax+by+cz+d＝0的四个参数a、b、c、d；

S842、计算任一点D_v到此平面的欧式距离，进而得到平面拟合残差δ，具体计算公式如下：

式中，n为点云簇中点云总数，(x_v，y_v，z_v)为点D_v的三维空间坐标。

10.根据权利要求1所述的面向地面LiDAR点云的建筑物信息提取方法，其特征在于，所述步骤S83具体为：

S831、在平面投影点集中选取Y坐标最小的一点G₀当作基点，如果存在多个点的Y坐标都为最小值，则选取X坐标最小的一点；

S832、对剩余点按G₀点的逆时针方向进行极角排序，当极角相同时，距离G₀比较近的点排在前面；

S833、根据“凸多边形的各顶点必在该多边形任意一条边同一侧”的原则对投影点集相继三个顶点连线进行一次有序扫描，删去投影边界序列中的非凸包顶点，最终剩余的点即为凸包的顶点；

S834、当点云数据被分为F层时，自下而上，第f层的投影面积记为S_f(f＝1，2，...，F)，设第f层的凸包多边形有L个顶点，则投影面积及分层投影面积比ΔS_F计算公式如下：

式中，x_i、y_i表示第i个顶点的坐标，x_i+1、y_i+1表示第i+1个顶点的坐标，max S_f为点云数据分层后最大的投影面积，min S_f为最小的投影面积。