一种融合多尺度特征的激光扫描数据物理平面自动化提取方法

摘要

本发明涉及一种融合多尺度特征的激光扫描数据物理平面自动化提取方法。本发明解决的现有问题主要包括：①发展了基于信息熵的最佳邻域自适应确定方法，克服了点密度变化、噪声、数据缺失等因素对三维点局部几何特征计算的影响；②综合利用基于点的特征（维数特征，法向量等）和基于区域的特征（粗糙度、紧凑度、尺度、长宽比等）进行物理平面区域提取，提高了物理平面提取的准确性；③扩展了现有平面分割方法的适用范围，本发明适用于机载、车载、固定站等多种类型激光扫描数据的物理平面提取。发明的整体技术流程图如下图。

说明书

技术领域

本发明属于计算机视觉和激光扫描数据处理的交叉领域，尤其涉及激光扫描 数据物理平面特征提取与自动化模型重建方法。

背景技术

激光扫描系统可以自动化的获取地物表面高精度的3维坐标信息，已成为 一种快速的空间数据获取手段，广泛运用于基础测绘、数字城市、交通运输、 森林资源调查、环境保护、文物保护与重建等领域。同时，激光扫描数据具有 数据量大，点云密度分布不均，场景目标多样(建筑物、道路、树木、车辆、 交通标志牌、交通信号灯等)，细节结构丰富等特点，难以直接对整个场景的激 光扫描数据进行三维重建。为了实现基于海量的激光扫描数据进行建筑物几何 模型三维重建，必须首先对激光扫描数据进行分割，识别出属于建筑物的面片 区域，然后基于提取的面片进行建筑物3维几何模型的重建。因此，建筑物平 面面片地自动、快速、准确提取和识别是建筑物三维重建的基础和前提，也是 本发明关注的重点。

目前从大范围城市复杂环境的激光扫描数据中快速提取物理平面的方法主 要包括：基于区域生长的分割方法，基于聚类的分割方法、基于辅助图像的分 割方法、基于随机采样一致性的分割方法等四类。BeltonandLichti(2001) 采用局部表面属性(点的邻近性、点到拟合平面的距离、平滑表面法向量夹角 等)作为相似性测度，然后将空间上邻近并且表面属性相似的点分割到同一区 域中。该分割方法对噪、点密度变化、数据遮挡缺失等比较敏感。Filin(2002) 首先计算激光扫描数据中所有的法向量和点到拟合平面的距离等特征，然后把 每个点从坐标空间转换到特征空间，最后利用模糊C均值方法在特征空间聚类 从而得到代提取的平面。该分割方法稳健性好，对邻域范围的大小不敏感，但 该方法需要人为指定面片个数的初始信息，因此自动化程度会受到影响。Van Hansenetal.(2006)利用了随机采样一致性方法来检测平面特征，该方法首 先将空间点云划分成3D空间体素，然后在每个体素内分别利用随机采样一致性 方法检测存在的平面，最后通过一定的合并规则将相邻体素中相似的平面特征 进行合并。杨必胜,魏征(2010)首先将激光扫描数据根据距离或者光谱信息 投影到特定的平面上形成深度图像或者彩色图像，然后利用图像处理中成熟的 平面检测算法进行平面提取，再将分割后的图像逆投影回激光扫描数据中，最 终实现激光扫描数据的分割。这种方法把3维的激光扫描数据转换成2维的图 像数据进行处理，加快了数据处理的效率，但是在转换的过程中会造成精度的 损失，同时还会受到图像投影误差、图像格网分辨率等因素的影响。

总体而言，从大范围城市复杂环境的激光扫描数据中快速、准确地提取平 面仍然存在：1)局部几何特征计算对点密度变化、噪声等影响比较敏感，导致 特征计算精度较低；2)平面提取的正确率和准确性相对较低，难以满足自动化 提取的需要；3)只适用于固定站、车载、机载扫描数据中某一单一类型数据 的平面提取，缺乏广泛的适用性，在实际生产中自动化程度低等问题。

发明内容

本发明针对现有技术的不足，提出一种融合多尺度特征的激光扫描数据物 理平面自动化提取方法。本发明解决的技术问题主要包括：1)发展了基于信息 熵的最佳邻域自适应确定方法，克服了点密度变化、噪声、数据缺失等因素对 三维点局部几何特征计算的影响；②综合利用基于点的特征(维数特征，法向 量等)和基于区域的特征(粗糙度、紧凑度、尺度、长宽比等)进行物理平面 区域提取，提高了物理平面提取的准确性；；③扩展了现有平面分割方法的适 用范围，本发明适用于机载、车载、固定站等多种类型激光扫描数据的物理平 面提取。

本发明的技术方案为一种融合多尺度特征的激光扫描数据物理平面自动化 提取方法，主要包含以下步骤：

一个激光扫描数据的步骤：利用机载、车载、固定站激光扫描仪获取原始 点云数据，该数据包括道路及道路两侧地物的高精度的3维坐标信息。

一个基于信息熵的最佳邻域确定和局部几何特征计算的步骤：针对激光扫 描获取原始数据确定局部几何特征计算的最大邻域半径r_max，最小邻域半径，并 计算该最佳邻域内每个点的几何特征，具体包括：

步骤1.1、确定局部几何特征计算的最大邻域半径r_max，最小邻域半径r_min， 半径的增量r_Δ，强度差阈值I_Δ，并初始化当前半径r_C＝r_min。

步骤1.2、逐步增大邻域半径r_C(r_C+＝r_Δ)，直到r_C≥r_max。对邻域r_C内与当 前点反射强度差小于I_Δ的数据进行主分量分析(PCA)得到点云数据分布的特征 值λ₁，λ₂，λ₃，(λ₁≥λ₂≥λ₃)，并定义维数特征：和熵函数E_f＝-a_1Dln(a_1D)-a_2Dln(a_2D)-a_3Dln(a_3D)。

步骤1.3、把函数E_f为最小值时的半径确定为最佳的邻域半径r_optimal，并计 算每个激光脚点的维数特征(a_1D，a_2D，a_3D)和法向量(N_x,N_y,N_z)。

一个基于支持向量机的逐点粗分类的步骤：根据计算的每个点的几何特征 利用支持向量机(SVM)对每个激光脚点分类为线状分布点、面状分布点和球状 分布点3类：具体方法是：

将每一个激光脚点的维数特征(a_1D，a_2D，a_3D)归一化到区间[-1,1]，手工 标记一部分点类别(线状分布点、面状分布点和球状分布点)并将这些已知类 别的点分为训练集和验证集；输入训练集进行训练(核函数采用高斯核函数)， 得到最优分类模型；用得到的模型对验证集和未知类别的点进行分类，分类为 线状分布点、面状分布点和球状分布点3类。

一个针对面状分布的点进行区域生长并提取初始平面的步骤：利用最小割 方法合并相邻的相似平面，具体包括：

步骤2.1、标记所有的面状分布的点为“未分割”，并从未分割的扫描点中 随机选取一个点作为生长的初始种子点；

步骤2.2、利用KD树搜索种子点的邻域点，如果邻域点法向量和种子点法 向量的夹角小于阈值T_a并且邻域点到种子点所在平面的距离小于阈值T_s，则认 为邻域点和种子点属于同一平面，并将该邻域点作为下次生长的种子点；

步骤2.3、重复步骤2.1至步骤2.3，直到所有的点都已分割完毕，得到初 始的分割区域。

步骤2.4、把初始的分割区域作为结点，构建无向加权图G＝(V,E,W)， V表示图中所有的节点，E表示节点间的边，W＝[w(i,j)]_n*n是权值矩阵，w(i,j)表 示节点i和j之间的相似性。假设将图G分为两个不相交的部分A与B，移去连 接A与B之间的边可以使该图一分为二，定义该图的一个割如下：

$cut(A,B)=\underset{i\in A,j\in B}{\Sigma }w(i,j)$

节点间的边反映了节点间的相似程度，权值越小，相似程度越低。每次选取图 的“最小割值”对图G进行二分，直至满足迭代终止条件。

一个基于语义知识的平面精细化提取的步骤，对提取的初始平面进行精化 处理，进一步剔除虚假的平面，具体方法是：计算每个初始片面的尺寸，包括 长、宽、高、面积、粗糙度其中l₁和l₃为对面片进行主分量分析(PCA) 得到的点云数据分布的特征值，紧凑度其中，area和perimeter为 区域的面积和周长，保留同时满足宽度大于阈值T_S、高度大于阈值T_H、面积 大于阈值T_β、粗糙度小于阈值T_r并且紧凑度大于阈值T_c的初始面片作为最终的 平面提取结果。

在上述一种融合多尺度特征的激光扫描数据物理平面自动化提取方法，在 基于信息熵的最佳邻域确定和局部几何特征计算的步骤中，确定局部几何特征 计算的最大邻域半径r_max，最小邻域半径r_min，半径的增量r_Δ，强度差阈值I_Δ，并 初始化当前半径r_C＝r_min。然后逐步增大邻域半径r_C，其中r_C+＝r_Δ，直到r_C≥r_max。 通过对不同邻域r_C内的数据进行主分量分析，发展了基于信息熵的最佳邻域自 适应确定方法，克服了点密度变化、噪声、数据缺失等因素对局部几何特征计 算的影响，从而提高了后续平面提取的完整性和正确率。

在上述一种融合多尺度特征的激光扫描数据物理平面自动化提取方法，在 基于支持向量机的逐点粗分类的步骤和针对面状分布的点进行区域生长并提取 初始平面的步骤中，利用支持向量机对每个激光脚点分类为线状分布点、面状 分布点和球状分布点3类，然后对分类为面状分布的点进行区域生长并提取初 始平面，并利用最小割方法合并相邻的相似平面。该方法提高了物理平面提取 的精度，同时降低了误提取的可能性。

本发明发展了基于信息熵的最佳邻域自适应确定方法，克服了点密度变化、 噪声、数据缺失等因素对局部几何特征计算的影响，提高了初始平面提取的精 度；同时，综合利用基于点的特征(维数特征，法向量等)和基于区域的特征 (粗糙度、紧凑度、尺度、长宽比等)进行平面区域提取，提高了平面提取的 准确性。

附图说明

图1是本发明实施例的基于知识的激光扫描数据平面自动化提取方法的整 体流程图。

图2是本发明实施例的基于信息熵的最佳邻域自适应确定方法的示意图。

图3是本发明实施例的基于SVM逐点粗分类的示意图。

图4是本发明实施例的利用“最小割”方法合并相邻权值构造的示意图。

具体实施方式

本发明主要基于计算机视觉和空间几何推理理论，提出一种融合多尺度特 征的激光扫描数据物理平面自动化提取方法。本发明综合利用基于点的特征(维 数特征，法向量等)和基于区域的特征(粗糙度、紧凑度、尺度、长宽比等) 进行平面区域提取，提高了物理平面提取的准确性，扩展了现有平面分割方法 的适用范围。

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实 施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅 用以解释本发明，并不用于限定本发明。本发明提供的方法能够用计算机软件 技术实现流程，整体技术流程图参见图1，包括以下步骤:

步骤1，为确定激光脚点的最佳邻域，首先需要设定最大邻域半径r_max＝0.5 米，最小邻域半径r_min＝0.1米，半径的增量r_Δ＝0.05米，强度差阈值I_Δ＝5，并 初始化当前半径r_C＝r_min，参见图2。实施例具体的实施过程说明如下：

逐步增大邻域半径r_C(r_C+＝r_Δ)，直到r_C≥r_max。利用邻域r_C内与当前点反射 强度差小于I_Δ的数据构建协方差矩阵M_3×3：

其中，k为邻域点个数，为每个邻域点坐标，为邻域 点坐标的平均值。把协方差矩阵特征值分解可以得到：

其中，λ₁,λ₂,λ₃(λ₁≥λ₂≥λ₃)为协方差矩阵M_3×3的特征值，为对应的特征值的 特征向量。并计算邻域点的维数特征：

$\begin{array}{ccc}{a}_{1D}=\frac{\sqrt{{\lambda }_1}-\sqrt{{\lambda }_2}}{\sqrt{{\lambda }_1}},& a_{2D}=\frac{\sqrt{{\lambda }_2}-\sqrt{{\lambda }_3}}{\sqrt{{\lambda }_1}}& a_{3D}=\frac{\sqrt{{\lambda }_3}}{\sqrt{{\lambda }_1}}\end{array}---\left(3\right)$

和熵函数：

E_f＝-a_1Dln(a_1D)-a_2Dln(a_2D)-a_3Dln(a_3D)(4)

选择使得熵函数E_f为最小值时的半径作为最佳的邻域半径r_optimal＝0.3米， 并把该邻域下计算得到的维数特征(a_1D，a_2D，a_3D)和法向量()用于步骤2 中。

步骤2，根据步骤1所计算的每个点的维数特征(a_1D，a_2D，a_3D)，利用支 持向量机(SVM)对每个激光脚点分类，把它们分类为线状分布点、面状分布点 和球状分布点3类，参见图3。实施例具体的实施过程说明如下：

利用公式将每一个激光脚点的维数特征(a_1D，a_2D，a_3D)归一 化到区间[-1,1]，在激光扫描数据中手工标记10000个点的类别，其中线状分 布点3000个、面状分布点4000个和球状分布点3000个，并将这些已知类别的 点分为训练集(线状分布点2000个、面状分布点3000个和球状分布点2000个) 和验证集(线状分布点1000个、面状分布点1000个和球状分布点1000个)两 部分；输入训练集进行训练(核函数采用高斯核函数)，得到支持向量机的最 优分类模型；用得到的分类模型对验证集点进行分类，获得线状分布点、面状 分布点、球状分布点分类的精度分别为95％，96％和94％；用得到的分类模型对 对为分类的进行分类，将它们分类为线状分布点、面状分布点和球状分布点3 类。

步骤3，对步骤2中分类为面状分布的点进行区域生长并提取初始平面，并 利用“最小割”方法合并初始的平面。实施例具体的实施过程说明如下：

步骤3.1:初始化所有的面状分布的点为“未分割”，并初始化队列Q为空。

步骤3.2:从“未分割”的扫描点中随机选取一个点作为生长的初始种子点 P_s，并初始化种子点P_s的区域标记L_C＝1。

步骤3.3:利用KD树搜索种子点P_s的邻域点P_j(j＝1,2…N)，N为邻域点的总 个数.如果同时满足条件：①邻域点P_j和种子点P_S之间的欧式距离小于阈值T_s(实施例中取值为0.5米)②邻域点P_j和种子点P_S的法向量夹角小于阈值(实 施例中取值为5°)③邻域点P_j到种子点P_S所在平面的垂直距离小于阈值T_h(实 施例中取值为0.2米)，则邻域点P_j和种子点P_S合并为同一区域，并把邻域点P_j压入队列Q的队尾，并赋值邻域点P_j的区域标记

步骤3.4:如果队列Q不空，则从队列Q的队首弹出一个点，并把该点作为 新的种子点P_S，并返回步骤3.3；如果队列Q为空，一个平面生长结束，返回步 骤3.2。

步骤3.5:重复步骤3.2-步骤3.4，直到所有的面状分布的点都已经被分割， 得到初始的平面分割区域。

步骤3.6:把步骤3.5获得的初始分割区域作为结点，构建无向加权图G＝ (V,E,W)，V表示图中所有的节点，E表示节点间的边，W＝[w(i,j)]_n*n是权值 矩阵，w(i,j)表示节点i和j之间的相似性。其中 $\omega \left(i,j\right)=e^{-\alpha \times {\mathrm{Pangle}}_{ij}}\times e^{-\beta \times d_{ij}}\times e^{-\gamma \times {\mathrm{PL}}_{ij}}\left(\alpha >0,\beta >0,\gamma >0,\alpha +\beta +\gamma =1.0\right),$分别为区 域i和区域j的法向量夹角值、相离距离(两个区域的最近距离)、区域i的中 心点到区域j的距离(点到面的距离)，α、β、γ(α＞0,β＞0,γ＞0；α+β+γ＝1)分别 为每个影响因子的权值调节参数，参见图4。

步骤3.7:每次选取图的“最小割值”对图G进行二分，直至满足迭代终止 条件(实施例中终止条件设为待合并的两个区域的相似度大于0.5)。其中，图 的“最小割值”为即区域A与区域B之间所有的边的权 值之和。

步骤4，对步骤3中提取的平面进行精化处理，进一步剔除虚假的平面。实 施例具体的实施过程说明如下：

步骤4.1:计算每个区域的尺寸(长L、高H、面积Area，周长Perimeter), H＝maxz-minz，Area＝L*H，Perimeter＝2*(L+H)，其中maxx， maxy，maxz，minx，miny，minz分别为区域的最大和最小的x,y,z值。

步骤4.2:计算每个区域的粗糙度(λ₁和λ₃计算参照步骤1)和紧 凑度$Compactness=\frac{Area}{Perimeter}.$

步骤4.3:遍历步骤3中所有的区域，保留同时满足宽度大于阈值T_S(实 施例中T_S取值为0.5米)、高度大于阈值T_H(实施例中T_H取值为0.5米)、面积 大于阈值T_β(实施例中T_β取值为0.5平方米)、粗糙度小于阈值T_r(实施例中T_r取值为0.1)并且紧凑度大于阈值T_c(实施例中T_c取值为0.3)的初始面片作为 最终的平面提取结果。

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