一种基于能量函数的室外点云场景中树木提取方法与流程

本发明属于人工智能计算方法技术领域，具体涉及一种基于能量函数的室外点云场景中树木提取方法。

背景技术：

随着虚拟现实、人工智能以及计算机视觉技术的发展，对各种三维室外场景的分析已经成为一个具有重要研究意义和应用价值的研究课题。而树木作为城市场景中的重要元素之一，已被广泛应用于城市规划建设、3d树木建模和城市树木检测等领域。近年来，研究人员已经提出了许多不同的方法从lidar数据中提取树木，主要就是一个分割和分类的过程，绝大部分算法都基于这种思想实现。目前，基于分割和分类的树木提取方法主要分为六大类：

(1)基于区域增长的树木提取：包括分割建筑结构和单棵树木。该方法提取的单棵树精度高，在相对简单的室外环境下表现出良好的性能。基于区域的方法分为两类：种子区域(或自下而上)方法和非种子区域(或自上而下)方法。基于种子的区域分割通过选择多个种子点来开始做分割，从这些种子点为起始点，通过添加种子邻域点的方式逐渐形成点云区域。这种方法处理时间长，且容易受噪声点的影响。后续提出了将坡度、法向量、曲率等特征作为生长准则来检测点云的新的改进方法，这类方法提高了分割的鲁棒性，对离群点和噪声点不敏感，但是如果初始种子点以及生长规则选取不当，也会导致过分割或欠分割等问题的出现。而非种子区域方法是一种基于自上而下的方法，首先把所有点归于同一个区域，然后对区域内的点进行细化处理，逐步分割成小的区域，直至划分完成，非种子区域方法的分割结果取决于细分的位置和方式，由于需要大量的先验知识比如对象模型，区域数量等，使得这一方法的使用也受到了很大的限制。

(2)基于属性特征的树木提取：该方法基于点云数据的属性进行分割，具有很好的鲁棒性。首先计算点云数据的属性，包括纹理特征，法向量，颜色信息等，然后对这些特征进行聚类，最终将不同特征的点云分割出来。这一方法依赖于点之间邻域的定义和点云数据的点密度，所以随着点云数据密度的增大，存在时间复杂度高、处理效率低的缺点，所以此类方法不适用于点云数据规模较大的处理。

(3)基于图优化的树木提取：基于图优化的方法是将点云数据当作顶点，利用点的空间邻域关系构造边，利用邻域点的相似性对连接边进行加权，构造成图，利用图割方法实现点云的分割。众所周知的图优化方法是基于图的分割算法(fh算法)，该方法简单且高效，用于在图中查找最小生成树。基于图优化的方法通常采用全局优化的策略求解，能够获得全局的优化结果，并且不受场景复杂程度和数学模型的影响。但是由于构造图的过程比较复杂，当点云数据规模较大时，优化处理通常无法实时运行。

(4)基于模型拟合的树木提取：主要用于提取具有特定形状的物体，包括平面，圆柱体和球体，需要检测点云对模型参数的拟合程度，对特定形状的对象数据中属于目标物的数据进行保留，再通过合并呈现出一棵完整的树。当使用基于模型拟合的方法时，每个提取的树簇都需要一个细化过程。这种方法具有纯粹的数学原理，快速且强大。其局限性在于对噪声敏感，所以在数据庞大、场景复杂的情况下，实验结果不是很理想。

(5)基于聚类的树木提取：将簇看成是数据空间中被低密度区域分割开的高密度区域。指导思想是：只要一个区域中的点的密度大过某个阀值，就把它加到与之相近的聚类中去。该方法被广泛应用于从mls数据中提取个体树，可以快速、成功地分离单棵树，但往往会存在聚类不足或过聚类的问题，导致树木提取不够完整。

(6)基于体素的树木提取：该方法在对点云数据进行预处理的时候先将点云转换为体素格式，以体素为单位进行邻域判断，通过对颜色，强度和点的空间位置的比较，实现点云数剧的分割，提取出候选树木；虽然也取得了不错的效果，但是其性能受点云分区不同引起的体素粒度和边界模糊的严重限制，导致损失分割精度。

技术实现要素：

本发明的目的是提供一种基于能量函数的室外点云场景中树木提取方法，解决了现有提取方法中分离准确率低、数据庞大运行慢问题。

本发明所采用的技术方案是，一种基于能量函数的室外点云场景中树木提取方法，包括以下步骤：

步骤1，选取lidar点云数据，剔除lidar点云数据中的平面点，得到候选树木点；

步骤2，将步骤1中候选树木点利用欧几里得聚类算法进行聚类处理，将所有聚类得到的簇保存在簇集c中；

步骤3，对步骤2得到的簇集c进行相似性矩阵计算；添加唯一性约束和条件性约束，随后计算簇集c的能量函数，根据能量函数对簇集c中的簇进行合并；

步骤4，根据步骤3中合并之后的簇集，对其进行筛选优化，对每个簇的点数进行判断，对点数不大于1000的簇进行删除，提取结束。

本发明的特征还在于，

步骤1具体为：

步骤1.1，在lidar点云数据中随机选取不少于三个且不共线的点，拟合呈一个平面，通过法向量计算平面的四个系数；

步骤1.2，计算lidar点云数据中剩余每个点至步骤1.1平面的距离，根据阈值判断距离平面较近的点，记为模型的局内点，记录模型的局内点数量；

步骤1.3，重复迭代，当满足迭代退出条件，则退出循环，得到整个迭代过程中最合理的解，即为平面点集；

步骤1.4，将步骤1.3中的平面点集剔除，得到候选树木点。

步骤1.3中迭代退出条件具体为：

步骤1.3.1，设w为随机抽到局内点的概率，n′表示计算模型参数需要选取的数据个数，得到随机抽取n′个点都是局内点的概率w^n′；计算至少存在一个局外点的概率，记为1-w^n′；计算迭代k′次都至少存在一个局外点的概率(1-w^n′)^k′；则迭代k′次至少可以获得一组正确模型参数的概率记为1-(1-w^n′)^k′；

步骤1.3.2，计算置信度p＝1-(1-w^n′)^k′；

步骤1.3.3，对步骤1.3.2中置信度p取对数，反推得到k′＝log(1-p)/log(1-w^n′)，k′值即为迭代终止条件，当实际迭代次数不小于k′时，退出循环。

步骤2具体为：

步骤2.1，利用kd-tree对步骤1中的候选树木点建立拓扑组织结构；

步骤2.2，选择步骤2.1中拓扑组织结构中的任意某点p，通过kd-tree近邻搜索算法找到k个距离p点最近的点，距离小于设定阈值的聚类到集合q中；

步骤2.3，选取集合q中除过p点的其余点点，通过kd-tree搜索k近邻，针对小于距离阈值的点继续将其添加到集合q中，重复直至q中元素的数目不在增加，划分完成一个簇；

步骤2.4，继续选取在集合q以外的点，重复步骤2.2和步骤2.3，直至拓扑组织结构中所有的点都被划分到某一个簇里面，将所有聚类得到的簇保存簇集c中。

步骤3具体为：

步骤3.1，设邻近矩阵的第i行第j个元素，簇ci的质心记作pi，簇cj的质心记作pj，计算欧几里得距离dis(pi,pj,ci,cj)、方向差值dir(pi,pj,ci,cj)、相似值w(ci,cj)如公式(1)-(3)：

dis(pi,pj,ci,cj)＝||pi-pj||(1)

w(ci,cj)＝λdis(pi,pj,ci,cj)+βdir(pi,pj,ci,cj)(3)

公式(1)-(3)中，参数v表示pi或pj的主方向，通过最小二乘线性拟合方法对pi或pj来局部拟合平面，然后求出最大特征值所对应的方向向量即为主方向；其中参数λ+β＝1，λ处于0到1之间且λ的取值取决于簇ci的密度，其中，参数p*属于簇ci中的任意一点，k″代表了簇ci中总点数；

步骤3.2，设步骤2中簇集c中包含n个簇，二维矩阵m表示每个簇之间的关系且为对称矩阵，取值为0或1，其中合并的原则为对二维矩阵m中值为1的簇进行合并。

步骤3.2具体为：

步骤3.2.1，基于欧几里得距离初始化二维矩阵m中的元素，若距离不大于0.2m，则将对应的二维矩阵m中的数据初始化为1，反之，将二维矩阵m中的数据初始化为0；

步骤3.2.2，添加唯一性约束如公式(4)，唯一性约束表示通过能量函数公式更新之后的二维矩阵m中每一行仅有一个元素被赋值1；使用mij来表示二维矩阵m中(i，j)的值，若能量函数更新之后mij＝1，则mji＝1，表示对应的簇ci和cj将进行合并；

步骤3.2.3，添加条件性约束如公式(5)，条件性约束表示根据两个簇的相似性矩阵，以最大的相似性值来找到最适合合并的簇；

其中参数sij,is是一个矩阵，其中sij,is中“ij”表示求取当前第i行第j列的值，“is”表示簇ci与其他簇的值，sij,is就表示簇ci和cj的权重，则为：

其中若j＝s，则γ＝0，反之γ＝1；

步骤3.2.4，结合步骤3.2.2的唯一性约束和步骤3.2.3的条件性约束，可得出能量函数，如公式(7)，

经对唯一性约束公式化简后，则更新前的能量函数pe'为公式(12)，

步骤3.2.5，更新之后能量函数peo和更新之前的能量函数pe'关系如公式(13)所示：

ed＝peo-pe'(13)

若由于j＝s，则s'pq,pq＝0，可得：

将公式(15)带入公式(14)中，可得：

由于ed≤0，提出更新的规则如公式(17)所示：

步骤3.2.6，通过步骤3.2.5中公式(17)来对二维对称矩阵m进行更新，根据更新的结果进行合并；

步骤3.2.7，对经步骤3.2.6处理后的二维矩阵中值为1的簇进行合并，重复执行以上步骤直至没有可以合并的簇。

本发明的有益效果是：本发明一种基于能量函数的室外点云场景中树木提取方法，通过有效地进行簇的合并，解决聚类不足的问题；将聚类完成之后不足成为一棵树的簇当做噪声点剔除，完成非树物体的去除提高了单棵树木提取的准确性，计算过程快，有很好的实用价值。

附图说明

图1为本发明一种基于能量函数的室外点云场景中树木提取方法的流程图；

图2为本发明树木提取方法中将原始数据映射二维平面的示意图；

图3为本发明树木提取方法中欧几里得聚类结果的示意图；

图4为本发明树木提取方法中能量函数方法聚类一次结果的示意图；

图5为本发明树木提取方法中能量函数方法最终聚类结果的示意图；

图6为本发明树木提取方法中优化处理结果的示意图。

具体实施方式

下面通过结合附图和具体实施方式对本发明进行详细说明。

本发明一种基于能量函数的室外点云场景中树木提取方法，如图1所示，包括以下步骤：

步骤1，如图2所示，选取lidar点云数据投射到二维平面上，剔除lidar点云数据中的平面点，得到候选树木点。初始的lidar点云数据分为地面点和非地面点，具体为：

步骤1.1，在lidar点云数据中随机选取不少于三个且不共线的点，拟合呈一个平面，通过法向量计算平面的四个系数；选取组成平面的点的个数要依据模型确定，一般来说选择适应于模型的最小数据个数，因为要拟合一个平面，所以最少需要选取三个点且不共线；

步骤1.2，计算lidar点云数据中剩余每个点至步骤1.1平面的距离，根据阈值判断距离平面较近的点，记为模型的局内点，记录模型的局内点数量；

步骤1.3，重复迭代，当满足迭代退出条件，则退出循环，得到整个迭代过程中最合理的解，即为能够最大程度拟合模型的平面点集；其中迭代退出条件具体为：

步骤1.3.1，设w表示随机抽到局内点的概率，n′表示计算模型参数需要选取的数据个数，计算随机抽取的n′个点都是局内点的概率w^n′；计算至少存在一个局外点的概率，记为1-w^n′；计算迭代k′次都至少存在一个局外点的概率(1-w^n′)^k′，当k′足够大时，概率(1-w^n′)^k′会很小，意味着模型可信度越高；则迭代k′次至少可以获得一组正确模型参数的概率记为1-(1-w^n′)^k′；

步骤1.3.2，计算置信度p＝1-(1-w^n′)^k′，表示ransac算法在运行后提供有用结果的期望概率；

步骤1.3.3，对置信度p取对数，反推k′＝log(1-p)/log(1-w^n′)，k′值即为迭代终止条件，当实际迭代次数不小于k′时，退出循环。

步骤1.4，将步骤1.3中的平面点集剔除，得到候选树木点。

步骤2，假设步骤1中任意一个候选树木点代表一个簇，利用欧几里得聚类算法对候选树木进行聚类处理，将所有聚类得到的簇保存在簇集c中。

如果每个候选树木点代表一个簇，直接对候选树木通过能量函数来进行聚类求解，这样在存储每个数据点与其余的数据点之间的相似性数据时，容易导致内存不足的情况。因此利用欧几里得聚类算法对候选树木进行聚类处理，把候选树木点分为多个簇，将初始的簇数量减少。

步骤2.1，利用kd-tree对步骤1中的候选树木点建立拓扑组织结构；具体为：

首先，划分维度，对所有维度数值计算方差，选择方差最大的维度d0；随后，选择切割标准的数据项，找出维度d0上所有数据项的中位数m，按m划分数据集，一分为二，记两个数据子集为dl和dr，建立树节点，存储本次划分；再次，对两个数据子集为dl和dr重复上述划分，并且将新生成的树节点设置为上一次划分的左右孩子；最后，充分划分至不能再划分为止，将对应的数据保存至最后的节点中，最后的节点也就是叶子节点。

步骤2.2，随意选择步骤2.1中拓扑组织结构中的任意某点p，通过kd-tree近邻搜索算法找到k个距离p点最近的点，距离小于设定阈值的便聚类到集合q中；

步骤2.3，在集合q中选取p点以外的点，通过kd-tree搜索k近邻，对于小于距离阈值的点再将其添加到集合q中，重复这一过程直至q中元素的数目不在增加，一个簇便划分完成；

步骤2.4，继续选取在集合q以外的点，重复步骤2.2和步骤2.3，直至拓扑组织结构中所有的点都被划分到某一个簇里面，将所有聚类得到的簇保存在簇集c中。

如图3所示，将图2中的点通过欧几里得聚类分为aa\bb\cc\dd\ee\ff\hh\ii\jj\gg\kk共计11部分。

步骤3，对步骤2得到的簇集c进行相似性矩阵计算；随后计算簇集c的能量函数，根据能量函数对簇集c中的簇进行合并。

通过对室外场景点云数据进行欧几里得距离聚类得到簇集c，目的其实就是将比较相似的数据样本合并到一起，然后通过计算相似性矩阵的值，找到两个最相似的簇并组合形成新的簇。但是上述方式来寻找最佳组合的簇的时候，只能体现出局部最佳的组合情况。因此，在对相似性求解的时候添加了唯一性约束和条件性约束来组成能量函数，通过能量函数来分析两个簇是否需要合并。基于能量函数最小化的原则，当一棵树被聚类成为多个簇时，能够有效地进行簇的合并。

步骤3.1，设邻近矩阵的第i行第j个元素，簇ci的质心记作pi，簇cj的质心记作pj，计算欧几里得距离dis(pi,pj,ci,cj)、方向差值dir(pi,pj,ci,cj)、相似值w(ci,cj)如公式(1)-(3)：

dis(pi,pj,ci,cj)＝||pi-pj||(1)

w(ci,cj)＝λdis(pi,pj,ci,cj)+βdir(pi,pj,ci,cj)(3)

公式(1)-(3)中，参数v表示pi或pj的主方向，通过最小二乘线性拟合方法对pi或pj来局部拟合平面，然后求出最大特征值所对应的方向向量即为主方向；其中参数λ+β＝1，λ处于0到1之间且λ的取值取决于簇ci的密度，其中，参数p*属于簇ci中的任意一点，k″代表了簇ci中总点数。

步骤3.2，设步骤2中簇集c中包含n个簇，二维矩阵m表示每个簇之间的关系且为对称矩阵，取值为0或1，其中合并的原则为对二维矩阵m中值为1的簇进行合并。具体为：

步骤3.2.1，基于欧几里得距离初始化二维矩阵m中的元素，若距离不大于0.2m，则将对应的二维矩阵m中的数据初始化为1，反之，将二维矩阵m中的数据初始化为0；

步骤3.2.2，添加唯一性约束，如公式(4)，唯一性约束表示通过能量函数公式更新之后的二维矩阵m中每一行仅有一个元素被赋值1，具体为每一次迭代聚类的过程中只会出现一对一的关系；使用mij来表示二维矩阵m中(i，j)的值，若能量函数更新之后mij＝1，则mji＝1，表示对应的簇ci和cj将进行合并；

步骤3.2.3，添加条件性约束，如公式(5)，条件性约束表示根据两个簇的相似性矩阵，以最大的相似性值来找到最适合合并的簇；

其中参数sij,is是一个矩阵，由于要求取当前第i行与簇ci最相似的另一个簇，故添加当前簇ci与其他簇对于计算簇ci和cj时的影响的值；其中sij,is中“ij”表示求取当前第i行第j列的值，“is”表示簇ci与其他簇的值，sij,is就表示簇ci和cj的权重，则为：

其中若j＝s，则γ＝0，反之γ＝1；

步骤3.2.4，结合步骤3.2.2的唯一性约束和步骤3.2.3的条件性约束，可得出能量函数，如公式(7)，

为了能够更加方便直观的求出能量函数的值，分别对唯一性约束和条件性约束的公式进行化简。

对唯一性约束公式进行化简：

其中的推导过程如公式(9)所示：

通过化简之后，能量函数pe如公式(10)所示：

则更新前的能量函数pe'为公式(11)，

将公式(10)带入公式(11)中，得到

步骤3.2.5，为了最小化能量函数，需要从更新之后能量函数peo和更新之前的能量函数pe'来进行分析，如公式(13)所示：

ed＝peo-pe'(13)

基于更新过程能量变化的分析，在公式中将从mpq变成m'pq，如公式(14)：

由于j＝s，则s'pq,pq＝0，可得：

通过公式(14)可以得出

因为始终要保持ed≤0，在这个基础上，提出更新的规则如公式(17)所示：

步骤3.2.6，通过步骤3.2.5中公式(17)来对二维对称矩阵m进行更新，根据更新的结果进行合并；若公式(17)计算值大于0，则将簇ci和簇cj在二维数组中所对应的值记为0；若公式(17)计算值等于0，则将簇ci和簇cj在二维数组中所对应的值记为1；若公式(17)计算值小于0，则簇ci和簇cj在二维数组中对应的值不变；

步骤3.2.7，对经步骤3.2.6处理后的二维矩阵中值为1的簇进行合并，重复执行以上步骤直至没有可以合并的簇。

如图4所示，将图3中的通过能量函数方法聚类一次，得到a-b(aa和bb合并)，c-d(cc和dd合并)，e-f(ff和ee合并)，hh，i-j(ii和jj合并)，g-k(gg和kk合并)共计六部分；经过重复的合并，最终的合并结果如图5所示，分为a-b，c-d-e-f(c-d和e-f合并)，hh，i-j-g-k(i-j和g-k合并)共计四部分。

步骤4，根据步骤3中合并之后的簇集，对其进行筛选优化，对每个簇的点数进行判断，对点数不大于1000的簇进行删除，提取结束。如图6所示，剔除掉图5中不大于1000的簇，剩余c-d-e-f和i-j-g-k为最终的提取结果。