CN102520401B

CN102520401B - 一种基于LiDAR数据的建筑物区域提取方法

Info

Publication number: CN102520401B
Application number: CN201110432421.8A
Authority: CN
Inventors: 程亮; 李满春; 刘永学; 陈焱明; 赵威; 张雯; 童礼华; 李�真; 蔡文婷; 杨康; 冯逍
Original assignee: Nanjing University
Current assignee: Nanjing University
Priority date: 2011-12-21
Filing date: 2011-12-21
Publication date: 2013-05-08
Anticipated expiration: 2031-12-21
Also published as: CN102520401A

Abstract

本发明公开了一种基于LiDAR数据的建筑物区域提取方法，属于LiDAR数据提取建筑物的方法领域。其步骤包括：LiDAR数据前期处理；对原始LiDAR数据重采样；反向迭代数学形态学滤波；分离建筑物与密集树木。本发明直接对3D点云数据进行处理，而不是将点云转化为深度图像，避免了转化过程中信息的损失和转化过程中增加的计算量；同时通过反向渐进迭代使用不同窗口进行数学形态学滤波操作，基本消除了地形起伏对数学形态学滤波中建筑物提取的影响，能取得较高的提取精度；本发明能够从大范围的LiDAR数据中快速、准确地提取出属于建筑物区域的LiDAR点，能够为城市的三维建模提供可靠的数据支持。

Description

一种基于LiDAR数据的建筑物区域提取方法

技术领域

本发明涉及一种从LiDAR数据提取建筑物的方法。尤其是一种采用反向迭代数学形态学方法从LiDAR点数据中提取建筑物点的方法。

背景技术

机载激光雷达(Airborne Light Detection And Ranging，简记LiDAR；又称：机载激光扫描，Airborne Laser Scanner，简记：ALS)是近十几年来快速发展的一种新型的测量技术，被广泛用于地面数据探测和模型的恢复、重建等应用中，显示出了巨大的前景，并逐渐成为三维城市数据模型获取的一种重要方法。

随着LiDAR数据质量的不断提高、越来越易获得、价格不断下降，从LiDAR数据中自动检测、重构三维地物的需求也将日益增长，而建筑物的检测和重构是城市3D建模中最受关注的。实时、准确地更新3D建筑物模型数据在城市规划、交通、游览等领域中十分重要。而如何从LiDAR点云中获得可靠的建筑物点是进行建筑物三维重建的关键步骤，建筑物点的分割效率、精度直接决定着三维重建的可行性和准确性。面对越来越广泛的LiDAR数据使用范围和使用需求，深入研究如何有效地从原始LiDAR数据中准确、有效地获取建筑物点具有重要的理论价值和实际价值。

从LiDAR数据中检测建筑物，即将打在建筑物上的点与打在其他地形表面的点(如地面、树木、道路)区分开来，其本质是一个分割过程。

在LiDAR研究领域，分割是最活跃的一个方向，大量的研究和应用与分割相关，因此也产生了许多的分割方法。根据所使用的数学方法不同，可以简单地分为五类：

1.边缘检测法：主要有(Heath，Sarkar et al.1996)，(Jiang and Bunke 1994)，(Sappa and Devy 2001)。由于在计算机视觉领域已有大量针对图像分割的边缘检测算法，所以LiDAR数据被转换为深度图像以适合图像边缘检测算法。然而将3D点云数据转换为2.5D深度图像的操作不可避免地带来的信息的丢失。对于机载LiDAR数据，由于重叠表面的存在(多层屋顶、桥梁、屋顶树枝等)，使用边缘检测方法进行分割将产生较多的误分割。对于地面LiDAR数据，由于点云是在不同位置扫描得到的，从3D转为深度图像时将导致大量的信息丢失。

2.表面生长算法：主要有(Gorte 2002)，(Lee and Schenk 2002)，(Rottensteinerand Briese 2003)，(Pu and Vosselman 2006)，(Rabbani，van Den Heuvel et al.2006)。对于LiDAR点云的表面生长相当于图像处理中的区域生长。首先，一些平面或非平面的区域被标记为种子，这些种子可以用比较成熟的最小平面调整、霍夫变换法得到。然后使这些种子区域不断向外扩展，扩展时一般是简单地考虑周围点的接近程度，坡度，曲率等。虽然在LiDAR点云分割中，表面生长算法使用广泛，但种子点的选择是难点，因为难以判断一组种子点比另一组可靠，而不同的种子点选择又会导致不同的分割结果。

3.扫描线法：主要有(Jiang and Bunke 1994)，(Sithole and Vosselman 2003)，(Khalifa，Moussa et al.2003)。扫描线法采用一种分裂-融合的策略。深度影像被分裂成沿指定方向的扫描线，比如每行视为一条扫描线。对于平的表面，任何一条三维平面上的扫描线都是一条直线。每一条扫描线都将被独立地分割成线段直到点到其对应直线的垂直距离小于某个固定的阈值。然后按照线段之间的关系融合成平面。扫描线法主要是为了提取平面而设计，然而在无结构的点云中实际上不存在扫描线。所以要将扫描线法引入到对点云的处理中需要决定预设方向的问题，并在点云的切片中进行扫描线构建，这就使得分割结果好坏取决于方向。

4.聚类方法：主要有(Roggero 2001)，(Filin 2002)，(Biosca and Lerma 2008)，(Chehata，David et al.2008)。在这些算法中，每个LiDAR点被赋予一个包含几个几何特征和光谱特征的特征矢量。然后LiDAR点就可以在特征空间通过K均值、最大似然、模糊聚类等聚类方法进行分割了。与其他方法不同，由于聚类是在特征空间进行的，所以可以按同样的方式处理点云、网格、TIN。聚类算法的效果依赖于特征矢量的和聚类方法的选取。

5.图形分区方法：主要有(Kim and Muller 1999)，(Fuchs 2001)。因为在同一分割中的点比不同分割中的点之间接触更加紧密，所以边界应当位于点之间间隔最大，即点最稀疏的地方。这个简单的想法被用来构造邻接图和邻域图。邻接图是在点云上的属性图G(V，E)。在一个小邻域里的每个或每组LiDAR点是集合V的一个节点，而集合E则由连接每对点的边构成，每条边有一个衡量点对间相似性的属性。通过图形分区算法寻找到同时满足最小化分割间点的相似性和最大化分割内部点的相似性的最优分割。分割可以通过递归分区或直接多路分区来实现。

现有算法较多是将LiDAR数据转化成深度影像，从而使用原有的图像分割算法进行处理，而非直接对3D点云数据进行处理。而数学形态学方法则是直接基于3D点数据进行处理的，不会因为数据转换而丢失信息。数学形态学是建立在严格的数学理论基础上的，从某种意义上讲，它是以几何学为基础的。它的主要思想是：利用被称之为结构元素(structuring element)的几何模板(例如：线、三角形、正方形、圆等)探测图像中物体的形状，具体通过数学形态学操作(如“开”或“闭”操作)来逼近地形表面。与其他方法相比，数学形态学方法概念简单，易于实现，并能取得较好的效果。但原有的数学形态学方法在结构元素的选择(形状和尺寸等)对处理后的结果有很大的影响，存在消除了地形起伏对区分地面点和建筑物的影响的缺陷，导致其精确性的不理想，需要根据具体的应用(如地面、植被、建筑物等)选择合适的结构元素，才能获得较为较为准确的提取结果。

发明内容

1.发明要解决的技术问题

本发明针对现有分割方法较多是将LiDAR数据转化成深度影像，从而使用原有的图像分割算法进行处理，而非直接对3D点云数据进行处理，以及对原有的数学形态学方法的缺陷，本发明提供一种基于LiDAR数据的建筑物区域提取方法，该方法主要是在数学形态学方法的基础上进行一定的改进，形成了一种新的反向迭代数学形态学滤波方法，此方法可以消除地形起伏对建筑物提取的影响，取得较好的提取效果。

2.技术方案

一种基于LiDAR数据的建筑物区域提取方法，包括以下步骤：

步骤1：LiDAR数据前期处理。主要是将原始LiDAR数据，由初始的数据格式(.xyz)转换为便于处理的矢量数据格式(.shp)，在此首先从xyz文件中读出每个点的三维坐标，存放到一个点集数组中，然后使用GDAL类库，将每个点的属性按顺序输出到一个shp文件中。

步骤2：对原始LiDAR数据重采样，即用最近邻法对原始无规则LiDAR数据进行重采样，得到规则间距的LiDAR点阵。

步骤3：反向迭代数学形态学滤波，即对重采样后的等间距点进行反向迭代数学形态学的开操作，以将地面点和稀疏树木点滤除，得到仅含建筑物区域和密集树木区域的点集。

步骤4：分离建筑物与密集树木，即根据树木点区域和建筑物点区域之间粗糙度的差异进行数学形态学的开操作，以将密集树木点和建筑物点分离，以达到精确提取出建筑物区域的目的。

上述步骤1中的原始LiDAR数据是空载LiDAR数据集，是将航空激光雷达扫描器安装在飞机平台上，通过平台的移动及仪器的摆扫，获得大面积的地表三维数据。一般可从测绘部门购买获得。

上述步骤2中的LiDAR数据重采样方法如下：按照1m的点间隔，采用最近邻法进行重采样，具体来说就是在LiDAR数据四至坐标范围内按照指定的采样间隔(1m)选择采样位置，对每个位置点，在一定邻域范围内(大于采样间隔)寻找最近邻的LiDAR点的高程值作为其高程值。

上述步骤3中的数学形态学开操作方法如下：对于任意一个LiDAR点p(x，y，z)和指定的邻域窗口ω，则在该点处对高程z的侵蚀为：

e_{p} = \min_{(x_{p}, y_{p}, eω)} (z_{p}) - - - (1)

其中(x_p，y_p，z_p)代表在p点邻域窗口ω内所有的点坐标值(窗口可为一维线形窗口，也可为二维矩形或其他形状窗口)，所有点的高程最小值记为e_p，即点p的侵蚀结果，侵蚀的结果是获得p点邻域内最小的高程值。

对应的，LiDAR数据的膨胀操作为：

d_{p} = \max_{(x_{p}, y_{p}, eω)} (z_{p}) - - - (2)

膨胀的结果是获得p点邻域内所有点的最大高程值为目标点p的高程值。先进行侵蚀操作，再进行膨胀操作称为开操作。而先进行膨胀操作，再进行侵蚀操作称为关闭操作。本发明主要采用数学形态学的开操作进行处理，滤波窗口使用的是二维窗口。

上述步骤3中的反向迭代数学形态学滤波方法如下：

首先选择一个最大窗口和一个最小窗口，从最大窗口开始，选择以一个较小的步长逐渐减小窗口大小实现采取逐步缩小窗口尺寸，对每个窗口都使用迭代地对经过上一步处理后的LiDAR数据进行数学形态学的开操作进行滤波，同时在这个过程中不断对相邻的两个窗口滤波结果进行作差，根据差值大小就可以用阈值来获取建筑物点了。在这个过程中，由于步长很小，相邻两个窗口滤波结果中地面点的下沉相差很小，影响可以忽略，而随着窗口的不断减小，建筑物会按照从大到小的顺序依次从地面“冒出”，迭代到最小窗口时，所有建筑物便都“冒出”了，当然还有大片的密集树木“冒出”，对这些所有“冒出”的地物做上标记最终即可得到所有的建筑物点和密集树木点。

上述步骤4中的建筑物与密集树木点分离方法如下：首先对滤波得到的LiDAR点以一定大小的邻域窗口来计算每个点的粗糙度(邻域高程的方差)，并根据一定的阈值，即根据粗糙度进行初步分类，将树木区域点和建筑物点分离，树木区域被大量移除，剩余的树木区域多为十分破碎、细长的区域；由于数学形态学滤波操作对细小地物的去除效果较好，进而采用数学形态学对初步分类结果进行开操作，对细碎的树木区域进行滤波，细小的破碎部分被消除，并在此基础上采用面积阈值进一步去除树木区域残留的细小碎片，稍大的区块滤波后变得圆滑，而大块的建筑物得到完好的保留；在前面的操作中，部分建筑物屋顶面片之间的边界处由于粗糙度较大被错误地滤除了，所以最后根据滤除点被大片建筑物区域的包围关系来判断是否为建筑物面片间的边界区域点，如果是边界点，根据其与建筑物屋顶面片的位置关系来恢复，将其重新标记为建筑物点。这样就完成了建筑物区域与密集树木区域的分离。

3.有益效果

本发明公开了一种基于LiDAR数据的建筑物区域提取方法数学形态学，与原有提取方法相比，有以下优势：

第一，本方法直接对3D点云数据进行处理，而不是将点云转化为深度图像，避免了转化过程中信息的损失和转化过程中增加的计算量；

第二，本方法通过反向渐进迭代使用不同窗口进行数学形态学滤波操作，基本消除了地形起伏对数学形态学滤波中建筑物提取的影响，因而能取得较高的提取精度；

综上所述，本发明直接基于LiDAR点云数据进行处理，同时采用反向迭代数学形态学方法滤除地面点和稀疏树木点，再根据高程分布的不同来区分树木区域和建筑物区域的点，能够取得较好的提取效果。应用本发明能够从大范围的LiDAR数据中快速、准确地提取出属于建筑物区域的LiDAR点，能够为城市的三维建模提供可靠的数据支持；而实时、准确地获取建筑物信息在城市规划、交通、旅游、土地资源监控等领域也有十分重要的应用价值。

附图说明

图1反向迭代数学形态学滤波方法技术流程图；

图2按高程着色的原始LiDAR数据图；

图3反向迭代数学形态学滤波结果图；

图4LiDAR点的粗糙度值图；

图5按照粗糙度初步分类树木点和建筑物点结果图；

图6对LiDAR点的粗糙度进行数学形态学滤波的结果图；

图7区域标号方法的示意图；

图8移除较小的树木区块结果图；

图9对被移除建筑物点的恢复结果图；

图10再次滤除剩余植被区域结果图；

图11建筑物点的分割结果图；

图12样区10栋建筑的分割结果图；

图13对样区影像手工数字化提取的10栋建筑物轮廓图；

具体实施方式

实施例1

本发明所提出的方法采用Visual Studio 2008开发平台C#语言编程实现，要处理的LiDAR数据为点状数据，在此选择采用开源地理数据格式转换类库OGR完成对点云数据的读写与其他处理操作。

参照技术方案的四个技术步骤和图1反向迭代数学形态学滤波方法技术流程图所示来完成建筑物点的提取，具体处理步骤如下：

步骤1：LiDAR数据前期处理

本发明采用大范围LiDAR数据进行反向迭代分割实验，实验区大约2000m×2000m，原始LiDAR数据如图2所示(按高程着色)，平均点间距约为1m，高程精度约为15cm，平面精度约为30cm，实验区内包含了不同屋顶结构、不同走向、不同复杂程度的多个建筑物，可以充分验证算法的鲁棒性。

原始LiDAR数据是.xyz格式的，为了对点进行标记，添加属性等处理操作，首先要将其转换为shapefile格式即转换为便于处理的矢量数据格式(.shp)，xyz格式是逐行记录每个点三个坐标信息的文件，我们逐行读取每组坐标，并用每行中的三个坐标构建一个三维点，最终使用OGR类库将这些三维点以shapefile文件形式输出。

步骤2：对原始LiDAR数据进行重采样

由于原始LiDAR数据是无序分布的点云，为了便于进行数学形态学操作，首先对原始数据进行重采样。因为原始数据平均点间距为1m，在此选择1m的间隔，用最近邻法进行重采样。重采样得到的规则间距LiDAR点存放到一个点矩阵中。

步骤3：反向迭代数学形态学滤波操作

进行反向迭代滤波操作首先需要选择一个最大滤波窗口大小和一个最小滤波窗口大小以及一个窗口渐变的步长。最大窗口是要用来滤除所有地物的，要求是比最大地物的尺寸更大一些，所以尽量选择大的尺寸即可；最小窗口是能为了够滤除稀疏的树木而保留建筑物，所以选择单棵树木的尺寸上限基本就可以了；窗口减小的步长对滤波结果影响不大，步长小则迭代次数多，计算时间长，步长大则迭代次数少，计算时间短，当然步长不宜超过20m。本文所采用的迭代滤波方法对窗口尺寸的选择要求十分宽松，窗口选择的差异对迭代滤波结果影响不大。在本实施例中，选择最大窗口为106m，最小窗口为6m，窗口减小步长10m，共迭代滤波11次。

其中反向迭代方法代码实现方式如下：

privatevoid iterate()

{

定义当前迭代窗口尺寸

for(当前迭代窗口尺寸从最大值按照步长减小到最小值)

{

用当前窗口对LiDAR数据进行侵蚀操作；

用当前窗口对LiDAR数据进行膨胀操作；

根据这一次开操作(侵蚀+膨胀)的结果与前一次开操作结果的高

程差来标记建筑物点；

}

其中侵蚀、膨胀操作均采用当前窗口大小的正方形窗口进行。侵蚀即对每一个点开一个窗口，用窗口区域内高程最小值去取代其原来的高程值；膨胀操作即对每一个点开一个窗口，用窗口区域内高程最大值去取代其原来的高程值。对于两次相邻窗口下开操作(侵蚀+膨胀)结果的高程差，用一个阈值来区分建筑物，由于在窗口减小的时候，当窗口小于建筑物尺寸时，建筑物点就不会被滤除而与前一窗口滤波结果有很大的高程差，所以此阈值的设置也十分宽松(建议2-5m)，在此选择3m作为阈值。

完成反向迭代数学形态学滤波操作后，地面点(建筑物所在地面上的平面点)和稀疏的树木点均被滤除，提取出了建筑物区域和密集树木区域。

步骤4：分离建筑物与密集树木

由前面的数学形态学迭代滤波得到如图3所示的建筑物和树木部分点集，空白部分表示被滤除的地面LiDAR点和稀疏树木点。树木区域虽然在高程上与建筑物接近，但其高程的分布与建筑物区域是有较大差异的，所以本发明采用高程表面的粗糙度来区分建筑物和树木。

(1)按粗糙度分类

对图3中的每个点计算粗糙度，以每个点为中心，开9m*9m的窗口进行粗糙度计算(高程的方差)得到如图4所示粗糙度分布图，图中白色表示粗糙度高的区域，黑色表示粗糙度低的区域，可见树木和建筑物屋顶面片边缘粗糙度较高，而且树木区域粗糙度变化频繁、无序。

以0.4为粗糙度阈值将所有点分为两类，得到如图5示所示的分类结果，可见，虽然树木部分有不少区域粗糙度很低，但都十分破碎，没有明显的轮廓，而建筑物每个屋顶面片都有大片的低粗糙度区域，轮廓清晰、明显。

(2)滤除破碎的树木区域

在按照粗糙度分类后，树木区域被大量移除，剩余的树木区域多为十分破碎、细长的区域，由于数学形态学滤波操作对细小地物的去除效果较好，所以本文采用数学形态学开操作对细碎的树木区域进行滤波，细小的破碎部分被消除，稍大的区块滤波后变得圆滑，而建筑物得到完好的保留，如图6所示(在此次滤波中使用的滤波窗口大小为5m*5m，滤波窗口较小是为了更好地保留建筑物的轮廓细节)。

由于滤波窗口较小，在数学形态学滤波的结果中(如图6)，仍有不少树木区域，但都是较小的碎片，其面积远小于建筑物区域的面积，故本发明进一步采用面积阈值来除去这些“树木区域碎片”。

由于所得结果均为点数据，要进行面积计算必须先按照同属性点的邻接关系进行分块，即将连为一体的同属性点标记为同一区块。

鉴于LiDAR数据是重采样成等间距的，本文采用一种效率很高的逐行扫描法来进行区域标号，该方法对所有LiDAR点进行两次逐行扫描即可完成区域标号：首先，从上往下逐行扫描各点，如果其上方、左侧、右侧相邻点已经被标号，则将当前点标记为与其一致的标号，否则将当前点标记上新的标号(如图7(a))。然而由于后面行连成一体的区块可能在前面行中并不相连，进行一次扫描后会存在同一区块被标记上不同标号的现象(如图7(b))。所以需要进行第二次扫描，通过从左往右逐列扫描各点，将相互邻接的标号不同的区块标记为同种区块。

在用逐行扫描法完成对连续的建筑物区块进行标号过程中，同时记录每个区块的大小。由于开操作后残留的植被区域较小，按照阈值300，将较小的区块移除，得到如图8所示结果。

(3)建筑物点的恢复

在前面的操作中，部分建筑物点被滤除，现根据滤除点被大片建筑物点包围和邻近情况对滤除点进行恢复，即根据其与建筑物屋顶面片的位置关系来恢复，具体来说是根据滤除点被大片建筑物点包围和邻近情况，遍历每个被标记为树木的点(即前面所滤除的点)，如果上下左右四个方向中有两个或以上方向有建筑物点且在距离阈值(设置为2m)内，则将此树木点标记为建筑物点。

遍历恢复完成后得到如图9所示结果，此时建筑物区域基本连接成较大区块，与误分的植被区域有较大面积差异。依据较大面积阈值600，将剩余植被区域重新滤除。得到最终的建筑物点分割结果，如10所示。

全区域建筑物分割结果如图11所示。图中共有建筑物173栋，最小面积139平方米，最大面积14211平方米，平均面积1475平方米。

从图12中可以看出，除了左下角的一些很小的建筑物之外，绝大多数建筑物都分割出来，遗漏较少，当然也有几个树木区域被误分为建筑物区域的。总体分割效果还是令人满意的。此区域LiDAR数据总共有280万个扫描点，属于大区域的数据，分割效果充分证明，本发明对大区域数据进行分割的鲁棒性，说明算法能够满足对大范围地面区域进行三维建模中的分割需求。

从全区域分割结果中选择一个较小区域单元来进行分割精度评估。图12为小区域内10栋完整建筑物的面状分割结果，右下角建筑物由于不完整，在此不进行精度评价。

精度参照数据是对此区域的航空影像中建筑物进行手工数字化的结果，由于影像分辨率极高(5cm)，所以数字化结果可靠性很强，可以用来作为真实数据如图13所示。

对所选区域内10栋建筑物的分割结果进行完整率和正确率评价得到结果如表1所示。

从表格中可以看出，10栋建筑物的分割正确率有9栋在90％以上，仅第10栋由于建筑物很小，而又紧邻树木区域所以正确率很低，而完整率都在85％以上，平均完整率和正确率都超过90％。

表1建筑物分割完整率和正确率

结合图12和图13可以看出，主要的误分割集中表现为两种情况，一是与建筑物连接在一起的树木，这样部分树木会被误分为建筑物；二是有些较细的建筑物由于某一方向的尺寸小于反向迭代滤波时最小的窗口尺寸，所以在滤波中没有提取出来。

Claims

1.一种基于LiDAR数据的建筑物区域提取方法，包括以下步骤：

步骤1：LiDAR数据前期处理，即将原始LiDAR数据转换为适合于处理的矢量数据格式；

步骤2：对原始LiDAR数据重采样，即用最近邻法对原始无规则LiDAR数据按照1m的数据点采样间距进行重采样，得到规则间距的LiDAR点阵；

步骤3：反向迭代数学形态学滤波，即对重采样后的等间距点进行反向迭代数学形态学滤波，以将地面点和稀疏树木点滤除，得到仅含建筑物区域和密集树木区域的点集，反向迭代数学形态学滤波的具体流程如下：

（1）选择最大滤波窗口、最小滤波窗口和窗口减小的步长；

（2）从最大窗口开始，到最小窗口结束，迭代地对重采样后的LiDAR数据进行数学形态学的开操作；

（3）迭代过程中，每次都将本次数学形态学的开操作结果与上一窗口下开操作结果进行作差，并根据这个高程差来初步区分建筑物区域的点和地面点；

步骤4：分离建筑物与密集树木，即根据树木点区域和建筑物点区域之间粗糙度的差异进行数学形态学的开操作，以将树木和建筑物点分离，这个分离过程的具体流程如下：

（1）对滤波得到的LiDAR点，计算每个点的粗糙度，并根据粗糙度进行初步分类即分成建筑物点和树木区域点；

（2）对初步分类结果进行数学形态学开操作，对细碎的树木区域进行初步滤除，并在此基础上采用面积阈值进一步去除树木区域残留的细小碎片；

（3）部分建筑物屋顶面片之间的边界处由于粗糙度较大，在滤波过程中被滤除，现根据其与建筑物屋顶面片的位置关系来恢复。