CN111709430B - 基于高斯过程回归的室外场景三维点云的地面提取方法 - Google Patents

Abstract

本发明涉及三维点云的地面提取方法，一种基于高斯过程回归的室外场景三维点云的地面提取方法，包括以下步骤：(1)获取室外场景的三维点云，(2)构建三维点云的邻域关系，(3)计算三维点云的协方差矩阵和法向量，(4)根据邻域形状对三维点云进行初步分类，(5)提取初步地面G_s，(6)分割初步地面，(7)二维高斯过程回归，(8)寻找局部地面邻域

(9)提取最终地面G_e。本发明利用逐层提取的思想和高斯过程回归的构建，从室外场景三维点云中准确完整地提取出了地面点云数据，有效地解决了由于室外场景复杂、地面零碎、起伏不定等因素而造成的地面提取不完整、不准确的问题。

Description

基于高斯过程回归的室外场景三维点云的地面提取方法

技术领域

本发明涉及一种三维点云的地面提取方法，更具体地说，涉及一种基于高斯过程回归的室外场景三维点云的地面提取方法。

背景技术

随着三维扫描测距技术的发展，三维点云在逆向工程、工业检测、自主导航等领域的应用越来越为广泛。三维点云处理技术作为实现上述应用的基础，发挥了至关重要的作用。在三维点云处理技术中，三维点云的特征提取是一个非常关键的技术，尤其是室外场景三维点云的地面特征提取，其对于室外场景的分割识别、移动机器人的路径规划等后续处理，都有着极其重要的作用。

在移动机器人自主导航领域，室外场景三维点云的地面提取是移动机器人进行路径规划的前提，完整的三维点云地面为移动机器人提供了准确的可行区域，提高了移动机器人的空间通过能力，保障了移动机器人在行进过程中的自身安全。在室外场景分析领域，由于室外场景极其复杂，会涉及到各种不同的物体，例如建筑、树木、车辆、人员等，为了便于进行场景分析，就需要对室外场景的三维点云进行有效的分割，而地面作为整幅场景的大背景，其准确完整的提取有助于将地面上的各类物体从空间中相互剥离，以便于后续的物体分割和场景分析。

目前，较为常见的三维点云地面提取方法为随机抽样一致性算法(RANSAC算法)，该方法将地面当作所处理场景中最大的平面而直接获取，这种方法对于比较平整大块的地面，具有较好的效果，但对于室外场景较为复杂、地面比较零碎、起伏不定的情况，这种方法就不能保证地面提取的完整性和准确性。

发明内容

为了克服现有技术中存在的不足，本发明目的是提供一种基于高斯过程回归的室外场景三维点云的地面提取方法。该方法针对一个室外场景，首先利用激光扫描测距仪获取室外场景的三维点云，其实质为三维空间中的一个点集，然后通过一定的点云提取方法，从室外场景三维点云中，准确完整地提取出地面点云数据。该方法解决了由于室外场景复杂、地面零碎、起伏不定等因素而造成的地面提取不完整、不准确的问题，以提高室外场景地面提取的准确性和完整性，具有较好的地面提取效果。

为了实现上述发明目的，解决已有技术中所存在的问题，本发明采取的技术方案是：一种基于高斯过程回归的室外场景三维点云的地面提取方法，包括以下步骤：

步骤1、获取室外场景的三维点云：利用激光扫描测距仪，获取室外场景的三维点云；

步骤2、构建三维点云的邻域关系：采用KD-Tree算法构建全体点云的结构树，根据点云的坐标将全体点云划分到不同的空间区域，在邻域构建时即可利用空间地址信息实现邻近点的搜索，以快速构建给定点p＝(x,y,z)的邻域N＝{p_i＝(x_i,y_i,z_i)|1≤i≤n_n}，其中：p_i为给定点邻点，i为给定点邻点的序号，n_n为给定点邻点的个数；

步骤3、计算三维点云的协方差矩阵和法向量：选取三维点云中的任意一点为给定点p＝(x,y,z)，利用其邻域N＝{p_i＝(x_i,y_i,z_i)|1≤i≤n_n}，构建协方差矩阵M，求解协方差矩阵M的特征值λ₁、λ₂、λ₃和特征向量v₁、v₂、v₃，以及给定点p的法向量n，具体包括以下子步骤：

(a)、利用步骤2所构建的三维点云的邻域关系，来快速构建给定点p＝(x,y,z)的邻域N＝{p_i＝(x_i,y_i,z_i)|1≤i≤n_n}；

(b)、构建给定点p的邻域N的协方差矩阵M，通过公式(1)进行描述，

式中：T为向量转置符号，其将列向量转置为行向量；

(c)、求解协方差矩阵M的特征值λ₁、λ₂、λ₃(λ₁<λ₂<λ₃)，以及对应的特征向量v₁、v₂、v₃；

(d)、将最小特征值λ₁对应的特征向量v₁单位化，即得到给定点p的法向量n；

(e)、对三维点云中的每一个点，重复步骤3子步骤(a)-(d)，继而求解每一个点的协方差矩阵的特征值、特征向量和法向量；

步骤4、根据邻域形状对三维点云进行初步分类：利用给定点p的协方差矩阵M的特征值λ₁、λ₂、λ₃之间的大小关系，来判断其邻域形状，并以此将三维点云分为三大类，即点性点集合C_p、线性点集合C_l和面性点集合C_s，具体包括以下子步骤：

(a)、若协方差矩阵M的特征值λ₁≈λ₂≈λ₃，即λ₃/λ₂≤8和λ₂/λ₁≤8，则给定点p与其邻点p_i呈散乱状分布，将给定点p分类为点性点；

(b)、若协方差矩阵M的特征值λ₁≈λ₂＜＜λ₃，即λ₃/λ₂>8和λ₂/λ₁≤8，则给定点p与其邻点p_i呈直线状分布，将给定点p分类为线性点；

(c)、若协方差矩阵M的特征值λ₁＜＜λ₂≈λ₃，即λ₃/λ₂≤8和λ₂/λ₁>8，则给定点p与其邻点p_i呈平面状分布，将给定点p分类为面性点；

(d)、对三维点云中的每一个点，重复步骤4子步骤(a)-(c)，将整个三维点云初步分为三大类，即点性点集合C_p、线性点集合C_l和面性点集合C_s；

步骤5、提取初步地面G_s：将所有面性点集合C_s的法向量映射到单位球上，以构建法向量球S，并利用Mean-Shift聚类算法，在法向量球S上对所有面性点集合C_s的法向量的顶点进行聚类，将面性点集合C_s分为若干个平面区域F_j，并从中提取出初步地面G_s，具体包括以下子步骤：

(a)、将面性点集合C_s内每一个面性点的法向量映射到单位球上，以构建法向量球S，法向量球S上的点即为面性点的法向量的顶点；

(b)、利用Mean-Shift聚类算法，在法向量球S上，对所有面性点的法向量的顶点进行聚类，进而可将面性点的法向量分为若干类，最终可将面性点分为若干个平面区域F_j，1≤j≤m，其中：j为平面区域的序号，m为平面区域的个数；

(c)、计算每一个平面区域F_j的平均高程

和平均法向量

若某个平面区域F_j的平均高程

和平均法向量

满足条件：

和

其中：

为平均法向量

与垂直方向的夹角，则认为该平面区域F_j为初步地面区域的一个组成部分，利用该方法对每一个平面区域F_j进行判断和筛选，则可获取整个的初步地面G_s；

步骤6、分割初步地面：使用K-Means聚类算法，根据初步地面G_s中的初步地面点p_t到各个局部地面中心点p_k之间的距离远近，将初步地面G_s细分为K个局部地面LG_k＝{p_ks＝(x_ks,y_ks,z_ks|1≤s≤n_k)}，1≤k≤K，其中，p_ks是局部地面点，k是局部地面序号，K是局部地面个数，s是局部地面LG_k中的点序号，n_k是局部地面LG_k包含的点数，具体包括以下子步骤：

(a)、确定局部地面的个数K，用最远点采样法选取K个点作为初始的局部地面中心点p_k，1≤k≤K；

(b)、计算初步地面点p_t(1≤t≤n_s)与各个局部地面中心点p_k之间的距离，把每个初步地面点p_t分配给距离它最近的局部地面中心点p_k，其中，t是初步地面点序号，n_s是初步地面点个数；

(c)、当初步地面点p_t均被分配完毕后，每一个局部地面中心点p_k会根据局部地面LG_k中的局部地面点p_ks被重新计算，通过公式(2)进行描述，

式中，p_k表示计算后新的局部地面中心点；

(d)、这个过程将不断重复，直到满足没有局部地面中心点p_k再发生变化为止；最后，将每一个局部地面中心点p_k和分配给它们的局部地面点p_ks看成一个簇，一共分为K簇，以此完成局部地面的划分，将初步地面G_s细分为K个局部地面LG_k(1≤k≤K)；

步骤7、二维高斯过程回归：构建局部地面LG_k的二维高斯过程回归模型，并以局部地面LG_k作为训练样本来训练二维高斯过程回归模型，用共轭梯度法优化方法求出该模型的超参数(l_k，

和

)的最优解，具体包括以下子步骤：

(a)、局部地面LG_k上定义离散函数z_ks＝f_d(t_ks),t_ks＝[x_ks,y_ks]^T，该函数在t_ks处的值f_d(t_ks)是高斯过程的随机变量；

(b)、通过均值函数m(t_ks)和协方差函数

指定高斯过程为，

于是，可以进一步将高斯过程GP写为

(c)、考虑到数据存在噪声，z_ks＝f_d(t_ks)+ε，在噪声

的假设条件下，得到观测值

的先验分布为，

其中：

K(T_k,T_k)＝(k_st)是输入样本T_k的一个n×n阶对称正定的协方差矩阵，矩阵元素k_st用来度量t_ks和t_kt之间的相关性，计算k_st的公式为平方指数协方差函数，

式中，l_k是方差尺度，

是信号方差，1≤s,t≤n_k,s≠t；

(d)、观测值z_k和预测值f_*的联合先验分布为，

K(T_k,t_*)＝K(t_*,T_k)^T是测试点t_*与训练输入T_k之间的n×1阶的协方差矩阵，K(t_*,t_*)是测试点t_*自身的协方差，I是n阶单位矩阵；由此可以计算出预测值f_*的后验分布，即要构建的局部地面LG_k的二维高斯过程归回数学模型为

f_*|T_k,z_k,t_*～N(m(f_*),cov(f_*)) (7)

式中，

(e)、用局部地面LG_k中的点作为训练样本，首先建立训练样本条件概率的负对数似然函数，并令其对超参数(l_k，

和

)求偏导；然后采用共轭梯度法优化方法对偏导数进行最小化，从而得到超参数的最优解；

步骤8、寻找局部地面邻域

对于某一个局部地面LG_k，逐个判断整个室外场景P＝{p_a＝(x_a,y_a,z_a)|1≤a≤n_a}中的场景点p_a，只要场景点p_a距离局部地面LG_k中任意一点p_ks(p_ks∈LG_k)小于距离阈值r_k，那么该点就归属于局部地面LG_k的一个邻域点，其中，a是场景点序号，n_a是整个室外场景的总点数，具体包括以下子步骤：

(a)、求出每一个局部地面点p_ks的地面点邻域

其中，n_ks是地面点邻点个数，

是LG_k中距离p_ks最近的点；

(b)、将步骤8子步骤(a)中得到的每一个地面点邻域

求并集，并将结果作为局部地面LG_k的邻域

记

其中，

是局部地面邻域点，s′是局部地面邻域点序号，

是局部地面邻域点个数；

步骤9、提取最终地面G_e：将每一个局部地面LG_k的邻域

代入到局部地面的二维高斯过程归回模型中，计算出每个

的预测均值

和方差

若预测均值和方差满足下列条件：

且

则

被认为是地面点，通过对每个局部地面LG_k使用上述方法检查其邻域

中的点，得到最终的地面G_e，其中

和

分别表示方差阈值和马氏距离阈值，

本发明有益效果是：一种基于高斯过程回归的室外场景三维点云的地面提取方法，包括以下步骤：(1)获取室外场景的三维点云，(2)、构建三维点云的邻域关系，(3)、计算三维点云的协方差矩阵和法向量，(4)、根据邻域形状对三维点云进行初步分类，(5)、提取初步地面G_s，(6)、分割初步地面，(7)、二维高斯过程回归，(8)、寻找局部地面邻域

(9)、提取最终地面G_e。与已有技术相比，本发明采用了逐层提取的思想：首先，通过协方差矩阵的特征值分析，来提取室外场景三维点云中的整体平面区域(面性点集合C_s)；然后，通过构建面性点C_s的法向量球S及其上的法向量聚类，将整体平面区域(面性点集合C_s)分为若干个平面区域F_j；接着，通过法向量信息和高程信息的结合，从若干个平面区域F_j中提取初步地面区域G_s；然后，用K-Means聚类算法将初步地面分割成若干个更紧凑的局部地面LG_k，并找出每个局部地面的邻域点；最终，通过局部地面区域的高斯过程回归，来获取完整地面区域G_e。逐层提取的思想能够使地面三维点云的提取更加完整和准确，尤其是当室外场景比较复杂时，将初步地面分割成若干个紧凑的局部地面，然后分别进行高斯过程回归，最后通过“或逻辑”获得完整地面，这对提取零碎、离散、起伏地面具有很大的帮助。因此，本发明所提方法有效地解决了由于室外场景复杂、地面零碎、起伏不定等因素而造成的地面提取不完整、不准确的问题，具有较好的地面提取效果。

附图说明

图1是本发明方法步骤流程图。

图2是室外场景三维点云显示图。

图3是本发明面性点提取结果图。

图4是本发明法向量球构建结果图。

图5是本发明法向量顶点聚类结果图。

图6是本发明初始地面图。

图7是本发明初始地面分割效果图。

图8是局部地面及其邻域图。

图9是本发明最终地面提取结果图。

具体实施方式

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

如图1所示，一种基于高斯过程回归的室外场景三维点云的地面提取方法，包括以下步骤：

步骤1、获取室外场景的三维点云，利用激光扫描测距仪，获取室外场景的三维点云。如图2所示，整幅室外场景由大约10万个点组成，其中包括地面、树木、草丛、建筑、车辆、人员等。

步骤2、构建三维点云的邻域关系：采用KD-Tree算法构建全体点云的结构树，根据点云的坐标将全体点云划分到不同的空间区域，在邻域构建时即可利用空间地址信息实现邻近点的搜索，以快速构建给定点p＝(x,y,z)的邻域N＝{p_i＝(x_i,y_i,z_i)|1≤i≤n_n}，其中：p_i为给定点邻点，i为给定点邻点的序号，n_n为给定点邻点的个数；

步骤3、计算三维点云的协方差矩阵和法向量：选取三维点云中的任意一点为给定点p＝(x,y,z)，利用其邻域N＝{p_i＝(x_i,y_i,z_i)|1≤i≤n_n}，构建协方差矩阵M，求解协方差矩阵M的特征值λ₁、λ₂、λ₃和特征向量v₁、v₂、v₃，以及给定点p的法向量n，具体包括以下子步骤：

(a)、利用步骤2所构建的三维点云的邻域关系，来快速构建给定点p＝(x,y,z)的邻域N＝{p_i＝(x_i,y_i,z_i)|1≤i≤n_n}；

(b)、构建给定点p的邻域N的协方差矩阵M，通过公式(1)进行描述，

式中，T为向量转置符号，其将列向量转置为行向量；

(c)、求解协方差矩阵M的特征值λ₁、λ₂、λ₃(λ₁<λ₂<λ₃)，以及对应的特征向量v₁、v₂、v₃；

(d)、将最小特征值λ₁对应的特征向量v₁单位化，即得到给定点p的法向量n；

(e)、对三维点云中的每一个点，重复步骤3子步骤(a)-(d)，继而求解每一个点的协方差矩阵的特征值、特征向量和法向量；

步骤4、根据邻域形状对三维点云进行初步分类：利用给定点p的协方差矩阵M的特征值λ₁、λ₂、λ₃之间的大小关系，来判断其邻域形状，并以此将三维点云分为三大类，即点性点集合C_p、线性点集合C_l和面性点集合C_s，具体包括以下子步骤：

(a)、若协方差矩阵M的特征值λ₁≈λ₂≈λ₃，即λ₃/λ₂≤8和λ₂/λ₁≤8，则给定点p与其邻点p_i呈散乱状分布，将给定点p分类为点性点；

(b)、若协方差矩阵M的特征值λ₁≈λ₂＜＜λ₃，即λ₃/λ₂>8和λ₂/λ₁≤8，则给定点p与其邻点p_i呈直线状分布，将给定点p分类为线性点；

(c)、若协方差矩阵M的特征值λ₁＜＜λ₂≈λ₃，即λ₃/λ₂≤8和λ₂/λ₁>8，则给定点p与其邻点p_i呈平面状分布，将给定点p分类为面性点；

(d)、对三维点云中的每一个点，重复步骤4子步骤(a)-(c)，将整个三维点云初步分为三大类，即点性点集合C_p、线性点集合C_l和面性点集合C_s，室外场景中面性点集合C_s的提取结果如图3所示。

步骤5、提取初步地面G_s：将所有面性点集合C_s的法向量映射到单位球上，以构建法向量球S，并利用Mean-Shift聚类算法，在法向量球S上对所有面性点集合C_s的法向量的顶点进行聚类，将面性点集合C_s分为若干个平面区域F_j，并从中提取出初步地面G_s，具体包括以下子步骤：

(a)、将面性点集合C_s内每一个面性点的法向量映射到单位球上，以构建法向量球S，如图4所示，法向量球S上的点即为面性点的法向量的顶点；

(b)、利用Mean-Shift聚类算法，在法向量球S上，对所有面性点的法向量的顶点进行聚类，如图5所示，进而可将面性点的法向量分为若干类，最终可将面性点分为若干个平面区域F_j，1≤j≤m，其中：j为平面区域的序号，m为平面区域的个数；

(c)、计算每一个平面区域F_j的平均高程

和平均法向量

若某个平面区域F_j的平均高程

和平均法向量

满足条件：

和

其中：

为平均法向量

与垂直方向的夹角，则认为该平面区域F_j为初步地面区域的一个组成部分，利用该方法对每一个平面区域F_j进行判断和筛选，则可获取整个的初步地面G_s，如图6所示。

步骤6、分割初步地面：使用K-Means聚类算法，根据初步地面G_s中的初步地面点p_t到各个局部地面中心点p_k之间的距离远近，将初步地面G_s细分为K个局部地面LG_k＝{p_ks＝(x_ks,y_ks,z_ks|1≤s≤n_k)}，1≤k≤K，其中，p_ks是局部地面点，k是局部地面序号，K是局部地面个数，s是局部地面LG_k中的点序号，n_k是局部地面LG_k包含的点数，具体包括以下子步骤：

(a)、确定局部地面的个数K，用最远点采样法选取K个点作为初始的局部地面中心点p_k，1≤k≤K；

(b)、计算初步地面点p_t(1≤t≤n_s)与各个局部地面中心点p_k之间的距离，把每个初步地面点p_t分配给距离它最近的局部地面中心点p_k，其中，t是初步地面点序号，n_s是初步地面点个数；

(c)、当初步地面点p_t均被分配完毕后，每一个局部地面中心点p_k会根据局部地面LG_k中的局部地面点p_ks被重新计算，通过公式(2)进行描述，

式中，p_k表示计算后新的局部地面中心点；

(d)、这个过程将不断重复，直到满足没有局部地面中心点p_k再发生变化为止；最后，将每一个局部地面中心点p_k和分配给它们的局部地面点p_ks看成一个簇，一共分为K簇，以此完成局部地面的划分，将初步地面G_s细分为K个局部地面LG_k(1≤k≤K)，如图7所示。

步骤7、二维高斯过程回归：构建局部地面LG_k的二维高斯过程回归模型，并以局部地面LG_k作为训练样本来训练二维高斯过程回归模型，用共轭梯度法优化方法求出该模型的超参数(l_k，

和

)的最优解，具体包括以下子步骤：

(a)、局部地面LG_k上定义离散函数z_ks＝f_d(t_ks),t_ks＝[x_ks,y_ks]^T，该函数在t_ks处的值f_d(t_ks)是高斯过程的随机变量；

(b)、通过均值函数m(t_ks)和协方差函数

指定高斯过程为，

于是，可以进一步将高斯过程GP写为

(c)、考虑到数据存在噪声，z_ks＝f_d(t_ks)+ε，在噪声

的假设条件下，得到观测值

的先验分布为，

其中：

K(T_k,T_k)＝(k_st)是输入样本T_k的一个n×n阶对称正定的协方差矩阵，矩阵元素k_st用来度量t_ks和t_kt之间的相关性，计算k_st的公式为平方指数协方差函数，

式中，l_k是方差尺度，

是信号方差，1≤s,t≤n_k,s≠t；

(d)、观测值z_k和预测值f_*的联合先验分布为，

K(T_k,t_*)＝K(t_*,T_k)^T是测试点t_*与训练输入T_k之间的n×1阶的协方差矩阵，K(t_*,t_*)是测试点t_*自身的协方差，I是n阶单位矩阵；由此可以计算出预测值f_*的后验分布，即要构建的局部地面LG_k的二维高斯过程归回数学模型为f_*

f_*|T_k,z_k,t_*～N(m(f_*),cov(f_*)) (7)

式中，

(e)、用局部地面LG_k中的点作为训练样本，首先建立训练样本条件概率的负对数似然函数，并令其对超参数(l_k，

和

)求偏导；然后采用共轭梯度法优化方法对偏导数进行最小化，从而得到超参数的最优解；

步骤8、寻找局部地面邻域

对于某一个局部地面LG_k，逐个判断整个室外场景P＝{p_a＝(x_a,y_a,z_a)|1≤a≤n_a}中的场景点p_a，只要场景点p_a距离局部地面LG_k中任意一点p_ks(p_ks∈LG_k)小于距离阈值r_k，那么该点就归属于局部地面LG_k的一个邻域点，其中，a是场景点序号，n_a是整个室外场景的总点数，具体包括以下子步骤：

(a)、求出每一个局部地面点p_ks的地面点邻域

其中，n_ks是地面点邻点个数，

是LG_k中距离p_ks最近的点；

(b)、将步骤8子步骤(a)中得到的每一个地面点邻域

求并集，并将结果作为局部地面LG_k的邻域

记

其中，

是局部地面邻域点，s′是局部地面邻域点序号，

是局部地面邻域点个数，局部地面邻域如图8所示。

步骤9、提取最终地面G_e：将每一个局部地面LG_k的邻域

代入到局部地面的二维高斯过程归回模型中，计算出每个

的预测均值

和方差

若预测均值和方差满足下列条件：

且

则

被认为是地面点，通过对每个局部地面LG_k使用上述方法检查其邻域

中的点，得到最终的地面G_e如图9所示。其中

和

分别表示方差阈值和马氏距离阈值，

本发明优点在于：在没有任何函数假设的情况下，例如线性或二次，高斯过程回归通过考虑函数值与给定观测值的相关性，可以很好地逼近未知函数值，即使是在突变的情况下。这一特点有利于灵活处理非常不规则的地面情况。本发明利用逐层提取的思想和高斯过程回归的构建，从室外场景三维点云中准确完整地提取出了地面点云数据，有效地解决了由于室外场景复杂、地面零碎、起伏不定等因素而造成的地面提取不完整、不准确的问题，具有较好的地面提取效果。

Claims (1)

1.一种基于高斯过程回归的室外场景三维点云的地面提取方法，其特征在于包括以下步骤：

步骤1、获取室外场景的三维点云，利用激光扫描测距仪，获取室外场景的三维点云；

步骤2、构建三维点云的邻域关系，采用KD-Tree算法构建全体点云的结构树，根据点云的坐标将全体点云划分到不同的空间区域，在邻域构建时即可利用空间地址信息实现邻近点的搜索，以快速构建给定点p＝(x,y,z)的邻域N＝{p_i＝(x_i,y_i,z_i)|1≤i≤n_n}，其中：p_i为给定点邻点，i为给定点邻点的序号，n_n为给定点邻点的个数；

步骤3、计算三维点云的协方差矩阵和法向量，选取三维点云中的任意一点为给定点p＝(x,y,z)，利用其邻域N＝{p_i＝(x_i,y_i,z_i)|1≤i≤n_n}，构建协方差矩阵M，求解协方差矩阵M的特征值λ₁、λ₂、λ₃和特征向量v₁、v₂、v₃，以及给定点p的法向量n，具体包括以下子步骤：

(a)、利用步骤2所构建的三维点云的邻域关系，来快速构建给定点p＝(x,y,z)的邻域N＝{p_i＝(x_i,y_i,z_i)|1≤i≤n_n}；

(b)、构建给定点p的邻域N的协方差矩阵M，通过公式(1)进行描述，

式中：T为向量转置符号，其将列向量转置为行向量；

(c)、求解协方差矩阵M的特征值λ₁、λ₂、λ₃(λ₁<λ₂<λ₃)，以及对应的特征向量v₁、v₂、v₃；

(d)、将最小特征值λ₁对应的特征向量v₁单位化，即得到给定点p的法向量n；

(e)、对三维点云中的每一个点，重复步骤3子步骤(a)-(d)，继而求解每一个点的协方差矩阵的特征值、特征向量和法向量；

步骤4、根据邻域形状对三维点云进行初步分类，利用给定点p的协方差矩阵M的特征值λ₁、λ₂、λ₃之间的大小关系，来判断其邻域形状，并以此将三维点云分为三大类，即点性点集合C_p、线性点集合C_l和面性点集合C_s，具体包括以下子步骤：

(a)、若协方差矩阵M的特征值λ₁≈λ₂≈λ₃，即λ₃/λ₂≤8和λ₂/λ₁≤8，则给定点p与其邻点p_i呈散乱状分布，将给定点p分类为点性点；

(b)、若协方差矩阵M的特征值λ₁≈λ₂＜＜λ₃，即λ₃/λ₂>8和λ₂/λ₁≤8，则给定点p与其邻点p_i呈直线状分布，将给定点p分类为线性点；

(c)、若协方差矩阵M的特征值λ₁＜＜λ₂≈λ₃，即λ₃/λ₂≤8和λ₂/λ₁>8，则给定点p与其邻点p_i呈平面状分布，将给定点p分类为面性点；

(d)、对三维点云中的每一个点，重复步骤4子步骤(a)-(c)，将整个三维点云初步分为三大类，即点性点集合C_p、线性点集合C_l和面性点集合C_s；

步骤5、提取初步地面G_s，将所有面性点集合C_s的法向量映射到单位球上，以构建法向量球S，并利用Mean-Shift聚类算法，在法向量球S上对所有面性点集合C_s的法向量的顶点进行聚类，将面性点集合C_s分为若干个平面区域F_j，并从中提取出初步地面G_s，具体包括以下子步骤：

(a)、将面性点集合C_s内每一个面性点的法向量映射到单位球上，以构建法向量球S，法向量球S上的点即为面性点的法向量的顶点；

(b)、利用Mean-Shift聚类算法，在法向量球S上，对所有面性点的法向量的顶点进行聚类，进而可将面性点的法向量分为若干类，最终可将面性点分为若干个平面区域F_j，1≤j≤m，其中：j为平面区域的序号，m为平面区域的个数；

(c)、计算每一个平面区域F_j的平均高程

和平均法向量

若某个平面区域F_j的平均高程

和平均法向量

满足条件：

和

其中：

为平均法向量

与垂直方向的夹角，则认为该平面区域F_j为初步地面区域的一个组成部分，利用该方法对每一个平面区域F_j进行判断和筛选，则可获取整个的初步地面G_s；

步骤6、分割初步地面，使用K-Means聚类算法，根据初步地面G_s中的初步地面点p_t到各个局部地面中心点p_k之间的距离远近，将初步地面G_s细分为K个局部地面LG_k＝{p_ks＝(x_ks,y_ks,z_ks|1≤s≤n_k)}，1≤k≤K，其中，p_ks是局部地面点，k是局部地面序号，K是局部地面个数，s是局部地面LG_k中的点序号，n_k是局部地面LG_k包含的点数，具体包括以下子步骤：

(a)、确定局部地面的个数K，用最远点采样法选取K个点作为初始的局部地面中心点p_k，1≤k≤K；

(b)、计算初步地面点p_t(1≤t≤n_s)与各个局部地面中心点p_k之间的距离，把每个初步地面点p_t分配给距离它最近的局部地面中心点p_k，其中，t是初步地面点序号，n_s是初步地面点个数；

(c)、当初步地面点p_t均被分配完毕后，每一个局部地面中心点p_k会根据局部地面LG_k中的局部地面点p_ks被重新计算，通过公式(2)进行描述，

式中，p_k表示计算后新的局部地面中心点；

(d)、这个过程将不断重复，直到满足没有局部地面中心点p_k再发生变化为止；最后，将每一个局部地面中心点p_k和分配给它们的局部地面点p_ks看成一个簇，一共分为K簇，以此完成局部地面的划分，将初步地面G_s细分为K个局部地面LG_k(1≤k≤K)；

步骤7、二维高斯过程回归，构建局部地面LG_k的二维高斯过程回归模型，并以局部地面LG_k作为训练样本来训练二维高斯过程回归模型，用共轭梯度法优化方法求出该模型的超参数(l_k，

和

)的最优解，具体包括以下子步骤：

(a)、局部地面LG_k上定义离散函数z_ks＝f_d(t_ks),t_ks＝[x_ks,y_ks]^T，该函数在t_ks处的值f_d(t_ks)是高斯过程的随机变量；

(b)、通过均值函数m(t_ks)和协方差函数k(t_ks,t_kt),

指定高斯过程为，

于是，可以进一步将高斯过程GP写为

(c)、考虑到数据存在噪声，z_ks＝f_d(t_ks)+ε，在噪声

的假设条件下，得到观测值z_k＝[z_k1,z_k2,…,z_knk]的先验分布为，

其中：

K(T_k,T_k)＝(k_st)是输入样本T_k的一个n×n阶对称正定的协方差矩阵，矩阵元素k_st用来度量t_ks和t_kt之间的相关性，计算k_st的公式为平方指数协方差函数，

式中，l_k是方差尺度，

是信号方差，1≤s,t≤n_k,s≠t；

(d)、观测值z_k和预测值f_*的联合先验分布为，

K(T_k,t_*)＝K(t_*,T_k)^T是测试点t_*与训练输入T_k之间的n×1阶的协方差矩阵，K(t_*,t_*)是测试点t_*自身的协方差，I是n阶单位矩阵；由此可以计算出预测值f_*的后验分布，即要构建的局部地面LG_k的二维高斯过程归回数学模型为

f_*|T_k,z_k,t_*～N(m(f_*),cov(f_*)) (7)

式中，

(e)、用局部地面LG_k中的点作为训练样本，首先建立训练样本条件概率的负对数似然函数，并令其对超参数(l_k，

和

)求偏导；然后采用共轭梯度法优化方法对偏导数进行最小化，从而得到超参数的最优解；

步骤8、寻找局部地面邻域

对于某一个局部地面LG_k，逐个判断整个室外场景P＝{p_a＝(x_a,y_a,z_a)|1≤a≤n_a}中的场景点p_a，只要场景点p_a距离局部地面LG_k中任意一点p_ks(p_ks∈LG_k)小于距离阈值r_k，那么该点就归属于局部地面LG_k的一个邻域点，其中，a是场景点序号，n_a是整个室外场景的总点数，具体包括以下子步骤：

(a)、求出每一个局部地面点p_ks的地面点邻域

其中，n_ks是地面点邻点个数，

是LG_k中距离p_ks最近的点；

(b)、将步骤8子步骤(a)中得到的每一个地面点邻域

求并集，并将结果作为局部地面LG_k的邻域

记

其中，

是局部地面邻域点，s′是局部地面邻域点序号，

是局部地面邻域点个数；

步骤9、提取最终地面G_e，将每一个局部地面LG_k的邻域

代入到局部地面的二维高斯过程归回模型中，计算出每个

的预测均值

和方差

若预测均值和方差满足下列条件：

且

则

被认为是地面点，通过对每个局部地面LG_k使用上述方法检查其邻域

中的点，得到最终的地面G_e，其中

和

分别表示方差阈值和马氏距离阈值，

Images