CN103020637A - 一种基于K-plane算法的建筑物顶面点云数据分割方法 - Google Patents

Abstract

本发明公开了一种基于K-plane算法的建筑物顶面点云数据分割方法。利用格网化方法，得出建筑物顶面点云数据的高程图像；通过对建筑物顶面点云数据的高程图像中各个像素的邻域法向量进行分类分析，将高程图像分割为一系列与建筑物顶面内各个几何平面相对应的区域，将这些区域的最小二乘拟合平面作为K-plane算法的初始聚类中心平面；利用K-plane算法，通过计算建筑物顶面点云数据中各个激光脚点与其各自所被分配的类的聚类中心平面之间的距离之和的最小值，得出点云数据的分类结果，实现对建筑物顶面点云数据进行分割。

Description

一种基于K-plane算法的建筑物顶面点云数据分割方法

【技术领域】

本发明涉及遥感测绘领域，尤其是一种基于K-plane算法的建筑物顶面点云数据分割方法。

【背景技术】

建筑物作为三维地理信息系统中的重要组成部分，对其实现快速、高精度探测已经成为遥感测绘领域内的重大需求。虽然采用传统测绘技术及方法，可以得到大多数建筑物的相关数据，但是其获取数据的速度慢，过程繁琐，很难满足对大规模建筑物进行快速、高精度探测的需求。因此，目前在城市区域内多采用激光雷达对大范围建筑物进行快速探测。如何对所得到的建筑物点云数据进行恰当的分类处理，从中准确地提取出建筑物的相关信息，是采用激光雷达进行建筑物探测过程中的一个至关重要的步骤。

在蒋晶珏，姚春静和马洪超申请的发明专利(申请号：200910272643.0)所述方法中，利用区域生长的方法，对激光雷达点云数据进行基于分割段的自动分类。该发明所述方法仅适用于对点云数据进行地物分离处理，通过该方法，可将点云数据分类为：裸露地面、建筑物、植被、噪声和其它非地面对象五类。该方法不适用于对提取出的建筑物点云数据进行进一步细化分割处理。在李传荣，周梅，苏国中，黎荆梅，唐伶俐和夏冰申请的发明专利(申请号：200810224910.2)所述方法中，聚类运算的方法被引入点云数据分割领域。但由于该发明所述方法以点间距离作为约束对候选点进行聚类，这样的操作同样仅适用于对建筑物顶面点云数据进行整体提取，而不适用于对获得的建筑物顶面点云数据进行进一步的分割处理。

目前，基于聚类算法的建筑物顶面点云数据分割方法通常依靠采用K-means算法对点云数据内的激光脚点进行聚类运算来实现。由MacQueen提出的K-means算法是一种基于距离划分对特征空间进行聚类分析的经典方法，其算法简洁，具有很强的数据搜索能力，适合处理点云数据这种数据量大的分类对象。由于K-means算法采用点模型作为聚类中心，利用该算法对建筑物顶面点云数据进行分割时，需要计算点云数据中各个激光脚点的邻域法向量作为聚类运算的数据对象。目前，计算激光脚点的邻近点主要有以下三种方法：同半径圆覆盖方法、K-nearest neighborhood(K-nn)算法和Voronoi网格方法。其中，同半径圆覆盖方法对点云数据的要求较高，在被处理的点云数据中不能存在空洞，而且激光脚点的分布必须较为均匀。K-nn算法通过求解出与各个激光脚点距离最近的K个邻近点，构建该最近邻域的最小二乘拟合平面，实现法向量提取。在实际应用中，如何选取合适的K值是基于K-nn算法的激光脚点法向量计算方法中所存在的主要问题。当所选取的K值较小时，邻域内包含激光脚点的数目较少，所构成的拟合平面受到测量误差的影响较大，计算出的法向量误差较大；当选取的K值较大时，邻域内包含激光脚点的数目较多，所构成的邻域拟合平面的细节表征能力较弱，也会导致计算出的法向量存在较大误差。同时，与同半径圆覆盖方法类似，K-nn算法也对被处理点云数据的分布均匀性有一定要求，当点云数据中激光脚点的分布密度不均匀时，采用K-nn算法进行法向量计算的结果并不理想。Voronoi网格方法凭借Voronoi网格的构造唯一性和良好的空间方向性，在实际应用中通常可以得到比较好的激光脚点邻域法向量计算结果。可是，Voronoi网格方法也有其自身存在的问题，首先，对数目巨大的点云数据构造Voronoi网格的过程十分复杂；其次，在进行邻域计算时，由于Voronoi网格方法所具有的良好的空间方向性，在获得较高的邻近点计算精度的同时，所引起的负面问题是所得出的每个激光脚点的邻近点数目很多，从而导致计算各个激光脚点邻域的拟合平面和进一步求取其法向量的运算量巨大；最后，由于建筑物顶面自身空间结构的特点，建筑物顶面内的各个几何平面的法向量之间的差距很小，同时由于激光雷达测量误差和建筑物表面粗糙性的存在，即使采用Voronoi网格方法也不能保证计算所得的法向量数据一定具有良好的可分割性。

本发明所述方法采用K-plane算法对建筑物顶面点云数据进行分类运算。与K-means算法类似，K-plane算法也是一种基于距离划分的聚类算法，其思想是采用平面模型代替K-means算法中的点模型作为聚类中心，并通过最小化数据对象与其各自所属类的聚类中心平面的距离之和的方法进行分类。因此，在对点云数据进行分割时，K-plane算法直接采用点云数据内各个激光脚点的三维坐标值作为聚类运算的数据对象。由于不需要计算各个激光脚点的邻域法向量，采用K-plane算法对建筑物顶面点云数据进行分割可以减少分类运算的计算时间，同时，可以避免由于激光脚点邻域法向量计算错误而导致错误分类出现的可能性。一直以来，阻碍K-plane算法在建筑物顶面点云数据分割领域应用的难题是K-plane算法的初始聚类中心平面很难被合理地设定。若所选取的初始聚类中心平面不合理，在聚类运算过程中，会出现由于迭代次数过多而引起计算时间过长和聚类结果落入局部最优解等情况。本发明所述方法通过对建筑物顶面点云数据的高程图像中各个像素的邻域法向量进行分析，将高程图像分割为一系列与屋顶的几何平面相对应的区域，利用这些区域的最小二乘拟合平面作为K-plane算法的初始聚类中心平面，进而减少聚类运算的迭代次数，缩短聚类运算的时间，同时降低聚类结果落入局部最优解的可能性。在实际应用中，与基于K-means算法的建筑物顶面点云数据的分割方法相比，本发明所述方法具有速度快、精度高等优点。

【发明内容】

本发明是一种基于K-plane算法的建筑物顶面点云数据分割方法。利用格网化方法生成建筑物顶面点云数据的高程图像。利用高程图像中各个像素点及其最小邻域内的邻域像素点的三维坐标值计算各个像素点邻域的最小二乘拟合平面。利用所得到的最小二乘拟合平面的法向量对高程图像中各个像素点邻域的法向量进行表示。利用阈值滤波和尺度变换的方法对高程图像中各个像素点的邻域法向量的三个坐标轴参数分别进行滤波处理并将其映射至图像灰度值的定义域内，对得到的映射结果进行取整处理，得出分别对应于法向量的三个坐标轴参数的灰度值矩阵，组合这三个灰度值矩阵生成高程图像像素点邻域法向量彩色图像。利用阈值分类方法对三个灰度值矩阵分别进行分类处理。利用得到的分类结果对高程图像进行分割，计算高程图像内各个分割区域的最小二乘拟合平面并将得到的拟合平面设置为K-plane算法中的初始聚类中心平面。利用K-plane算法对建筑物顶面点云数据进行聚类运算，将得到的聚类结果作为点云数据的分割结果进行输出。

本发明的技术方案是这样实现的：

本发明的一种基于K-plane算法的建筑物顶面点云数据分割方法，步骤如下：

步骤一、对建筑物顶面点云数据进行格网化处理，生成高程图像。

在笛卡尔三维xyz坐标系内，将建筑物顶面点云数据沿z轴方向映射到xoy平面上，选取建筑物顶面点云数据的扫描线平均间距和扫描线内扫描点平均间距两者之中的较大值作为格网尺寸，构建正方形格网对点云数据进行分割。当格网内仅包含一个激光脚点时，采用该激光脚点的高程值作为格网的高程值；当格网内包含一个以上激光脚点时，采用这些激光脚点的高程平均值作为格网的高程值；当格网内不包含激光脚点时，采用该格网的最小邻域内所有格网的高程平均值作为其高程值。将所得到的包含高程值的格网作为图像的像素点，得出建筑物顶面点云数据的高程图像。

步骤二、计算高程图像内各个像素点的邻域法向量。

将高程图像中各个像素点的高程值作为其z轴坐标，利用各个像素点及其最小邻域内各像素点的x轴、y轴和z轴坐标计算出各个像素点邻域的最小二乘拟合平面，将所得到的最小二乘拟合平面的法向量作为各个像素点邻域的法向量。

步骤三、对高程图像内像素点的邻域法向量的各坐标轴的参数分别进行图像化表示。

对高程图像内像素点的邻域法向量的三个坐标轴的参数分别进行阈值滤波处理，利用尺度变换的方法将滤波后的三个坐标轴参数分别映射至图像灰度值的定义域[0，255]内。通过对所得到的参数进行取整处理，得到分别对应法向量x轴、y轴和z轴的三个灰度值矩阵。将这三个矩阵作为彩色图像的三原色矩阵，对其进行组合，得出建筑物顶面点云数据高程图像中像素点邻域法向量的彩色图像。

步骤四、对图像化表示后的法向量进行阈值分类，根据分类结果对建筑物顶面的高程图像进行分割，将分割得到的各个区域的最小二乘拟合平面设置为K-plane算法中的初始聚类中心平面。

利用阈值分类方法对图像化表示后的法向量三个坐标轴的参数分别进行分类，利用各个坐标轴参数的分类结果对高程图像进行分割，利用分割得到的各个区域内像素点的三维坐标计算各个区域的最小二乘拟合平面，将得到的拟合平面设置为K-plane算法中的初始聚类中心平面。

步骤五、利用K-plane算法对建筑物顶面点云数据进行聚类运算，得出点云数据分割结果。

利用K-plane算法对建筑物顶面点云数据进行聚类运算，采用平面作为聚类中心模型，比较激光脚点与各个聚类中心平面之间的欧式距离，将该激光脚点分配给距离其最近的聚类中心平面所代表的类。将点云数据中所有的激光脚点都按照这样的方法进行分配，当点云数据中所有激光脚点都被分配后，将各个类的聚类中心平面更新为各个类内激光脚点的最小二乘拟合平面。再次比较各个激光脚点与各个更新后的聚类中心平面之间的距离，重新对点云数据中的各个激光脚点进行分配，待分配结束后，再次更新聚类中心平面，按照这样的方法进行迭代运算。当所得的分类结果不再随迭代运算发生变化或者迭次次数已经满足要求，将聚类运算的分类结果作为点云数据的分割结果输出。

本发明的有益效果：本发明所述方法采用平面作为聚类中心模型，通过比较激光脚点的三维坐标与各个聚类中心平面之间的欧式距离进行分类运算，因此，在采用本发明所述方法进行点云数据分割过程中，避免了采用基于K-means算法的点云数据分割方法中对点云数据内各个激光脚点的邻域法向量进行计算的步骤，这样即缩短了点云数据分割的整体运算时间，又可以避免由于激光脚点邻域法向量计算的误差而导致点云数据被错误分类的可能性。本发明所述方法中的K-plane算法的初始聚类中心平面是通过对点云数据的高程图像中各个像素的邻域法向量进行分析得到的，相对于传统的随机初始化方法，采用本发明所述的初始化方法可以有效的减少点云数据分割过程中聚类运算的迭代次数，缩短分类运算的时间，并且，可以降低得到的分类结果陷入聚类运算的局部最优解的可能性。

【说明书附图】

图1是本发明所述的一种基于K-plane算法的建筑物顶面点云数据分割方法的示意图；

图2是利用本发明所述方法对一真实建筑物顶面点云数据进行分割的具体实施例的示意图。

【具体实施方式】

本发明的一种基于K-plane算法的建筑物顶面点云数据分割方法，步骤如下：

步骤一、对建筑物顶面点云数据进行格网化处理，生成高程图像。

在笛卡尔三维xyz坐标系内，将建筑物顶面点云数据沿z轴方向映射到xoy平面上，选取建筑物顶面点云数据中的扫描线平均间距和扫描线内扫描点平均间距两者之中的较大值作为格网尺寸，构建正方形格网对点云数据进行分割。当格网内仅包含一个激光脚点时，采用该激光脚点的高程值作为格网的高程值；当格网内包含一个以上激光脚点时，采用这些激光脚点的高程平均值作为格网的高程值；当格网内不包含激光脚点时，采用该格网的最小邻域内所有格网的高程平均值作为其高程值。将所得到的包含高程值的格网作为图像的像素点，得出建筑物顶面点云数据的高程图像。

步骤二、计算高程图像内各个像素点的邻域法向量。

将高程图像中各个像素点的高程值作为其z轴坐标，计算各个像素点最小邻域(3*3邻域)内的邻域像素点。将高程图像中的像素点(i，j)记为P_(i，j)，其邻域切平面Plane_(i，j)的数学表达式记为：A(_i，j)x+B(_i，j)y+C(_i，j)z+D(_i，j)＝0，其中，x轴和y轴分别表示高程图像的横坐标轴和纵坐标轴，z轴表示高程值的坐标轴，A_(i，j)，B_(i，j)，C_(i，j)和D_(i，j)分别表示切平面Plane_(i，j)的相应参数。该像素点(i，j)及其在3*3邻域内的全部邻近点到拟合切平面Plane_(i，j)的距离的平方和表示为：

采用最小二乘拟合的方法进行处理，计算出当sumdist_(i，j)取最小值时所对应的拟合切平面Plane_(i，j)的相关参数A^* _(i，j)，B^* _(i，j)，C^* _(i，j)和D^* _(i，j)，从而得出像素点P_(i，j)的最佳拟合切平面。利用其最佳拟合切平面的方向信息对该像素点邻域的标准化法向量进行表示：

其中，

为标准正交化运算的公分母因子，NVx_(i，j)，NVy_(i，j)和NVz_(i，j)分别为标准化法向量NV_(i，j)在x轴，y轴和z轴上的参数。当像素点位于高程图像的边缘和角点位置处不满足3*3邻域要求时，分别利用其最邻近的5个像素点或3个像素点采用上述方法计算最佳拟合切平面，进而得出相应的法向量信息。

步骤三、对高程图像内像素点的邻域法向量的各坐标轴的参数进行图像化表示。

为了使高程图像中像素的邻域法向量信息更易观测和分析，对其进行图像化表示。针对法向量的x轴、y轴和z轴三个坐标轴上的参数，分别构建三个与之相对应的色彩矩阵M_R，M_G和M_B。利用阈值滤波和尺度变换的方法将高程图像中各个像素点的标准化法向量信息分别映射到这三个色彩矩阵中。由于机载激光雷达探测是一个自上而下的扫描过程，所得的激光脚点都来自于地面或者地物的上表面，因此，在高程图像中，对应各个像素点的法向量z轴的参数NVz_(i，j)不可能为负数。由此可知，对应法向量x轴的参数NVx_(i，j)的分布区间为[-1，1]；对应y轴的参数NVy_(i，j)的分布区间为[-1，1]；而对应z轴的参数NVz_(i，j)的分布区间为[0，1]。为强化法向量各个坐标轴参数的色彩对比度，对其进行阈值滤波处理。以对应法向量x轴的参数NVx_(i，j)为例，计算出NVx_(i，j)负值部分和正值部分的平均值NVx^- _(i，j)和NVx⁺ _(i，j)分别作为其下限和上限滤波阈值，即： ${{\mathrm{NVx}}^{*}}_{(i,j)}=\left\{\begin{array}{cc}{\mathrm{NVx}}_{(i,j)}^{+},& \mathrm{if}{\mathrm{NVx}}_{(i,j)}\≥{\mathrm{NVx}}_{(i,j)}^{+}\\ {\mathrm{NVx}}_{(i,j)}^{-},& \mathrm{if}{\mathrm{NVx}}_{(i,j)}\≤{\mathrm{NVx}}_{(i,j)}^{-}\\ {\mathrm{NVx}}_{(i,j)},& \mathrm{else}\end{array}\right..$ 采用尺度变换的方法将阈值滤波后的法向量x轴参数NV^*x_(i，j)映射至色彩矩阵M_R中： $M_{R(i,j)}=255\·({{\mathrm{NVx}}^{*}}_{(i,j)}-{\mathrm{NVx}}_{(i,j)}^{-})/({\mathrm{NVx}}_{(i,j)}^{+}-{\mathrm{NVx}}_{(i,j)}^{-}).$ 通过这样的操作，色彩矩阵M_R中参数显示的对比度得到提高，并可以抑制建筑物顶面的表面突起物和激光雷达测量误差引起法向量参数出现跳变的情况。采用同样的方法对法向量y轴的参数NVy_(i，j)进行阈值滤波，并将其映射至对应的色彩矩阵M_G中。由于法向量z轴的参数NVz_(i，j)为非负形式分布，因此计算其正值部分的平均值作为下限滤波阈值NVz^- _(i，j)，不再设定上限滤波阈值，对其进行滤波处理： ${{\mathrm{NVz}}^{*}}_{(i,j)}=\left\{\begin{array}{cc}{\mathrm{NVz}}_{(i,j)}^{-},& \mathrm{if}{\mathrm{NVz}}_{(i,j)}\≤{\mathrm{NVz}}_{(i,j)}^{-}\\ {\mathrm{NVz}}_{(i,j)},& \mathrm{else}\end{array}\right..$ 将滤波后的法向量z轴的参数映射至M_B中： $M_{B(i,j)}=255\·(\mathrm{NV}{z^{*}}_{(i,j)}-{\mathrm{NVz}}_{(i,j)}^{-})/\left(\max \left({\mathrm{NVz}}^{*}\right){-\mathrm{NVz}}_{(i,j)}^{-}\right),$ 其中，max(NVz^*)表示滤波后的法向量z轴参数的最大值。通过上述运算，得到分别与法向量三个坐标轴参数相对应的色彩矩阵MR，M_G和M_B，将这三个色彩矩阵作为三原色矩阵进行组合，可得出高程图像像素点邻域法向量的彩色图像M_RGB：M_RGB＝(M_R，M_G，M_B)。

步骤四、对图像化表示后的法向量进行阈值分类，根据得到的分类结果对建筑物顶面的高程图像进行分割，将分割后各个区域的最小二乘拟合平面设置为K-plane算法中的初始聚类中心平面。

当高程图像的像素点P_(i，j)处于建筑物顶面的边缘位置，该像素点的邻域切平面在笛卡尔三维xyz空间坐标系内近似垂直于xoy平面，因此，滤波后的法向量z轴参数NVz^* _(i，j)应趋近于0。另一方面，当高程图像的像素点P_(i，j)处于建筑物顶面的内部区域时，该像素点的邻域切平面在三维xyz空间坐标系内应近似平行于xoy平面，因此，滤波后的法向量z轴参数NVz^* _(i，j)应该趋近于1。采用二值化方法对M_B矩阵中的建筑物顶面边缘区域进行检测： $M_{B(i,j)}^b=\left\{\begin{array}{cc}255,& \mathrm{if}{M}_{B(i,j)}\&GreaterEqual;255/2\\ 0,& \mathrm{if}{M}_{B(i,j)}<255/2\end{array}\right..$ 当M^b _B(i，j)等于0时，其所对应的高程图像像素点P_(i，j)为建筑物顶面边缘区域的像素。色彩矩阵M_R和M_G中参数的定义域为[0，255]，通过对定义域三等分的方法，得出用于分类处理的上限阈值TH_P＝2/3·255和下限阈值TH_N＝1/3·255。采用阈值分类方法将M_R中的参数分为三类： ${M^{*}}_{R(i,j)}=\left\{\begin{array}{cc}255,& \mathrm{if}{M}_{R(i,j)}>{\mathrm{TH}}_P\\ 0,& \mathrm{if}{M}_{R(i,j)}<{\mathrm{TH}}_N\\ M_{R(i,j)}& \mathrm{else}\end{array}\right..$ 当M_R(i，j)大于TH_P时，该参数被重新赋值为255，其所对应的高程图像中的P_(i，j)像素点被判断为是属于建筑物顶面中法向量在x轴上处于正方向的几何平面的。当M_R(i，j)小于TH_N时，该参数被重新赋值为0，其所对应的高程图像中的P_(i，j)像素点被判断为是属于建筑物顶面中法向量在x轴上处于负方向的几何平面的。除此以外，当M_R(i，j)处于[TH_N，TH_P]区间内，该参数值保持不变，其所对应的高程图像中的P_(i，j)像素点被判断为是属于建筑物顶面中法向量与x轴近似垂直的几何平面的。根据色彩矩阵M^* _R中参数的分类结果，该色彩矩阵所对应的图像被分割为一系列区域。采用相同的方法将M_G中的参数分为三类： ${M^{*}}_{G(i,j)}=\left\{\begin{array}{cc}255,& \mathrm{if}{M}_{G(i,j)}>{\mathrm{TH}}_P\\ 0,& \mathrm{if}{M}_{G(i,j)}<{\mathrm{TH}}_N\\ M_{G(i,j)}& \mathrm{else}\end{array}\right..$ 当M_G(i，j)大于TH_P时，该参数被重新赋值为255，其所对应的高程图像中的P_(i，j)像素点被判断为是属于建筑物顶面中法向量在y轴上处于正方向的几何平面的。当M_G(i，j)小于TH_N时，该参数被重新赋值为0，其所对应的高程图像中的P_(i，j)像素点被判断为是属于建筑物顶面中法向量在y轴上处于负方向的几何平面的。除此以外，当M_G(i，j)处于[TH_N，TH_P]区间内，该参数值保持不变，其所对应的高程图像中的P_(i，j)像素点被判断为是属于建筑物顶面中法向量与y轴近似垂直的几何平面的。根据色彩矩阵M^* _G中参数的分类结果以及已经获得的M^* _R图像的分割结果，对色彩矩阵M^* _G所对应的图像进行分割。将所得到的M^* _G分割结果与高程图像像素邻域法向量的彩色图像M_RGB中的色彩分布进行比较，若M^* _G的分割结果与M_RGB中的色彩分布相一致，利用M^* _G的分割结果对点云数据的高程图像进行分割。利用高程图像中分割所得各个区域内的像素的三维坐标值，计算出各个区域的最小二乘拟合平面，将得到的最小二乘拟合平面作为K-plane算法的初始聚类中心平面。

步骤五、利用K-plane算法对建筑物顶面点云数据进行聚类运算，得出点云分割结果。

在处理建筑物顶面点云数据分割问题时，K-plane算法中的数据对象为点云数据中N个激光脚点的三维坐标数据，设K-plane算法共设置M个聚类中心平面，其中，第i个激光脚点表示为v_i＝{x_i，y_i，z_i}，1≤i≤N，第j个聚类中心平面记为ccp_j，1≤j≤M。第j个聚类中心平面ccp_j的数学表达式为：w_j*(x，y，z)^T+r_j＝0，其中，w_j和r_j分别表示聚类中心平面的单位法向量和过原点垂线上的截距，T表示矩阵的转置。第i个激光脚点与其中的第j个聚类中心平面之间的距离为：d(v_i，ccp_j)＝|w_j*(x_i，y_i，z_i)T+r_j|。按照这样的方法计算出第i个激光脚点与其余M-1个聚类中心平面之间的距离。通过比较得到的M个距离数值，提取出与第i个激光脚点距离最近的聚类中心平面ccp_k：

将第i个激光脚点分配给聚类中心平面ccp_k所代表的类，记为C_k。按照相同的方法将其余N-1个激光脚点分别分配至M个类内。计算各个类内所包含激光脚点的最小二乘拟合平面，将其更新为各个类新的聚类中心平面。将N个激光脚点重新分配至更新后的各个类内，再次更新聚类中心平面，再次将N个激光脚点重新分配，以此类推，进行迭代运算。K-plane算法的目标函数表示为：

在迭代运算中，每次N个激光脚点分配完毕后都计算目标函数J_K-plane的数值，随着迭代运算的进行，当J_K-plane收敛为一稳定不变的最小值时，将所得到的最优聚类结果v_i∈C_k(1≤i≤N，1≤k≤M)作为点云数据的分割结果输出。

通过应用本发明所述方法对一真实建筑物顶面点云数据进行分割的具体实施例对本发明所述方法的实施方式做进一步说明：

应用本发明所述方法对一真实建筑物顶面点云数据进行分割。该点云数据中建筑物顶面包含三个几何平面如图2(a)所示，该点云数据所对应的扫描区域的遥感影像如图2(b)所示。对该点云数据进行滤波处理，将其中所包含的建筑物顶面点云数据提取出来，如图2(c)所示。利用格网化的方法得到该点云数据的高程图像，如图2(d)所示。对高程图像中各个像素点的最小邻域法向量进行计算。按照本发明所述方法对法向量x轴，y轴和z轴上的参数分别进行滤波处理，并分别映射至其各自所对应的色彩矩阵M_R，M_G和M_B中。得到的对应滤波后的法向量z轴坐标参数的色彩矩阵M_B如图2(e)所示，利用二值化处理方法从中得出建筑物顶面的边界信息。进而利用下限阈值TH_N和上限阈值TH_P对色彩矩阵M_R中的参数进行分类处理，利用得到的分类结果对M^* _R图像进行分割，如图2(f)所示。利用阈值TH_N和TH_P对色彩矩阵M_G中的参数进行分类处理，利用得到的分类结果和M^* _R图像的分割结果对M^* _G的图像进行分割，如图2(g)所示。通过组合M_R，M_G和M_B得到高程图像像素邻域法向量的色彩图像M_RGB，将M^* _G的分割结果与M_RGB中的色彩分布进行比对，如图2(h)所示。通过对图2(h)观察可知，M^* _G的分割结果与M_RGB中的色彩分布相符合。按照M^* _G的分割结果对高程图像进行分割，如图2(i)所示。计算高程图像中各个被分割区域的最小二乘拟合平面，将得到的拟合平面设置为K-plane算法的初始聚类中心平面。以点云数据中各个激光脚点的三维坐标作为数据对象，利用K-plane算法进行聚类运算。当聚类运算结束，将聚类结果输出作为点云数据的分割结果，如图2(j)所示。通过利用遥感图像进行比对可知，利用本发明所述方法，该建筑物顶面点云数据被合理地进行了分割。

以上所述，仅为本发明具体实施方法的基本方案，但本发明的保护范围并不局限于此，任何熟悉本技术领域的人员在本发明公开的技术范围内，可想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求的保护范围为准。所有落入权利要求的等同的含义和范围内的变化都将包括在权利要求的范围之内。

Claims (1)

1.一种基于K-plane算法的建筑物顶面点云数据分割方法，其特征在于包括以下五个步骤：

步骤一、对建筑物顶面点云数据进行格网化处理,生成高程图像；

在笛卡尔三维xyz坐标系内，将建筑物顶面点云数据沿z轴方向映射到xoy平面上，选取建筑物顶面点云数据的扫描线平均间距和扫描线内扫描点平均间距两者之中的较大值作为格网尺寸，构建正方形格网对点云数据进行分割，当格网内仅包含一个激光脚点时，采用该激光脚点的高程值作为格网的高程值，当格网内包含一个以上激光脚点时，采用这些激光脚点的高程平均值作为格网的高程值，当格网内不包含激光脚点时，采用该格网的最小邻域内所有格网的高程平均值作为其高程值，将所得到的包含高程值的格网作为图像的像素点，得出建筑物顶面点云数据的高程图像；

步骤二、计算高程图像内各个像素点的邻域法向量；

将高程图像中各个像素点的高程值作为其z轴坐标，利用各个像素点及其最小邻域内各像素点的x轴、y轴和z轴坐标计算出各个像素点邻域的最小二乘拟合平面，将所得到的最小二乘拟合平面的法向量作为各个像素点邻域的法向量；

步骤三、对高程图像内像素点的邻域法向量的各坐标轴的参数进行图像化表示；

对高程图像内像素点的邻域法向量的三个坐标轴的参数分别进行阈值滤波处理，利用尺度变换的方法将滤波后的三个坐标轴参数分别映射至图像灰度值的定义域[0,255]内，通过对所得到的参数进行取整处理，得到分别对应法向量x轴、y轴和z轴的三个灰度值矩阵，将这三个矩阵作为彩色图像的三原色矩阵，对其进行组合，得出建筑物顶面点云数据高程图像中像素点邻域法向量的彩色图像；

步骤四、对图像化表示后的法向量进行阈值分类，根据分类结果对建筑物顶面的高程图像进行分割，将分割得到的各个区域的最小二乘拟合平面设置为K-plane算法中的初始聚类中心平面；

利用阈值分类方法对图像化表示后的法向量三个坐标轴的参数分别进行分类，利用各个坐标轴参数的分类结果对高程图像进行分割，利用分割得到的各个区域内像素点的三维坐标计算各个区域的最小二乘拟合平面，将得到的拟合平面设置为K-plane算法中的初始聚类中心平面；

步骤五、利用K-plane算法对建筑物顶面点云数据进行聚类运算，得出点云分割结果；

利用K-plane算法对建筑物顶面点云数据进行聚类运算，将获得的聚类运算结果作为点云数据的分割结果进行输出。

Images