基于地形的LiDAR数据缺失区域填充方法研究

谷延超，范东明

(西南交通大学 地球科学与环境工程学院，四川 成都 611756)

基于地形的LiDAR数据缺失区域填充方法研究

谷延超，范东明

(西南交通大学 地球科学与环境工程学院，四川 成都 611756)

针对机载激光点云数据栅格化过程中出现的空白区域进行研究，在对邻近填充和最低点填充分析的基础上，提出基于地形的填充方法。对缺失数据边界进行一维形态学滤波得到其边界各点地形高，利用边界点构建不规则三角网进行地形内插填充。实验结果表明，地形填充法可适应不同原因造成的区域数据缺失，能够有效保证高程的连续性并可提高滤波精度。

机载LiDAR；栅格化；数据填充；不规则离散点；形态学滤波

机载激光雷达(Airborne Light Detection and Ranging)作为一种主动的新型测量技术，能够快速获取大范围高精度的空间地理信息，得到了广泛的应用[1-3]。LiDAR的原始三维激光点云作为数字表面模型，不仅包含地面点，还包含树木、建筑、车辆等地物信息[4]。滤波是从数字表面模型中提取地面点，进而生成数字地面模型(DEM)，是LiDAR数据处理的重要步骤，同样也是后续众多应用的前提。

数字形态学作为图像处理的重要部分，已经广泛应用于点云滤波，并取得了较好的效果[1-3 , 5-7]。数字形态学的滤波算法可以是基于原始LiDAR 点云数据，也可以是基于栅格数据进行，实验表明，后者无论在时间还是精度上都优于前者。由于点云数据是不规则分布的散点，因此在滤波前须对数据进行栅格化[2,6]。栅格化数据是形态学滤波算法的基本数据，将影响着滤波精度以及数据分层渲染的效果等，尤其是大面积数据缺失填充方法对滤波精度的影响，但至今未见此方面详细的对比分析。本文详细介绍了数据栅格化的基本流程，在对现有大范围数据缺失的填充方法分析的基础上提出了一种新的大面积数据缺失填充的方法——地形填充法，并对各种填充方法进行详细的评价。

1 点云栅格化

LiDAR获取的数据为一系列离散的空间三维坐标及反射强度，点云的位置和间隔在空间上都表现为离散性和不规则性。基于规则格网的形态学滤波可直接利用现有的形态学滤波方法进行大数据的处理[1]，故需要对原始的不规则离散数据进行栅格化。

规则栅格的格网大小应根据点云数据的空间间隔确定，一般情况取平均点云间隔以达到点云数据利用最大化。栅格化过程中，如果一个栅格内落入一个或几个激光脚点，把高程最小值记录在单元格内，若栅格内无激光脚点，则需对单元格进行填充[2,5,6]。在栅格化过程中，需要生成栅格高程影像和栅格掩膜影像。高程影像记录落入栅格内激光脚点的高程最低值；而掩膜影像记录栅格有无激光脚点落入，存在激光脚点的掩膜图像栅格值为1，否则为0。

2 空白栅格填充

空白栅格大致可分为两类：①由数据稀疏造成的点状或者小面积的数据缺失(Ⅰ类数据缺失)；②由高吸收性地物或者航带无重叠度造成的大面积数据缺失(Ⅱ类数据缺失)。

当激光脚点的点间距大于影像的空间分辨率时，会产生Ⅰ类数据缺失[2]，可利用限制搜索半径的邻近激光点对高程影像和掩膜影像进行填充，搜索半径可根据设备的扫描方式以及栅格大小来设定，通常可选取1～2倍的栅格大小。

产生Ⅱ类数据缺失的原因主要可归为两类：①相邻的扫描带无重叠区域；②强吸收性的地物对电磁波的吸收(以水体为主)，填充时都应加以考虑。针对Ⅱ类数据缺失通常有两种策略：邻近填充和最低点填充[5]。

2.1 邻近填充

邻近填充利用相邻的地物性质一致进行填充，假设栅格高程影像为I，点云数据集为P，若数据缺失区域为L，则对L内的所有单元格逐一进行最近点搜索，步骤如下：

1)以空白单元格lij为中心，统计搜索半径r内的所有激光脚点集N。

2)若N为空，则扩大搜索半径r，重复1)直至脚点集N为非空为止。

3)计算lij与邻近激光脚点集N中所有脚点的距离，距离最小值所对应的高程即为待填充点的高程。

2.2 最低点填充

文献[5]中认为：水面往往接近其周边范围内的最低点，故提出了最低点填充。对消除Ⅰ类数据缺失后的掩膜影像进行空白区域探测，统计所有的空白区域，依次对Ⅱ类数据缺失利用边界最低点进行填充，步骤如下[5]：

1)对掩膜影像进行形态学腐蚀运算得到腐蚀后的掩膜影像，将原掩膜影像与腐蚀后的掩膜影像相减即可得到数据缺失的边界。

2)提取 Ⅱ 类数据缺失的边界中对应高程影像I的最小值，赋给高程影像中对应的 Ⅱ 类数据缺失区域。

3)依次对影像进行1)、2)步骤，直至所有的Ⅱ类数据缺失全部填充完成。

图1显示边界提取的过程。图1(a)为未经任何填充的原始掩膜影像；图1(b)为经过Ⅰ类数据缺失填充后，残留的大面积数据缺失即Ⅱ类数据缺失；图1(c)为数据缺失区域腐蚀后的掩膜影像；图1(d)为Ⅱ类数据缺失的边界。

图1 边界提取示意图

通常情况下，针对由航带无重叠造成的Ⅱ类数据缺失进行邻近填补为较优，针对水体吸收造成的Ⅱ类数据缺失进行最低点填补较优[5]，然而在现实数据处理过程中，无法通过空白区域判断其产生原因，故不能按照上述分类进行填充。

2.3 地形填充

邻近填充会改变地物本身的属性，例如尺寸等，而在形态学滤波中需要设定的最大滤波窗口与建筑物的尺寸有关，最大窗口必须大于滤波区域内最大建筑物的尺寸，这样才能保证在滤波过程中将建筑物过滤[1,5]。同时邻近填充可能在数据缺失区域中部形成明显的高程跳跃，对后续的滤波处理会造成一定的影响。

水体造成的缺失可分为由于湖泊和河流造成。由湖泊造成的数据缺失采用最低点填充方法往往与现实地形一致；但河流往往具有一定的地形起伏，尤其在大范围内的河流落差较大，此时利用最低点填充的数据并未与现实地形保持一致。

综上所述，最低点填充是为了更好地接近水体的高程即数字高程模型；而邻近点填充目的是更好地接近周围地物或地形，注意到如果邻近点为地形点，则是将填充点更接近数字高程模型，当其为地物时，也要在后续的滤波过程中进行剔除，以得到相应点的地面高程。故无论最低点填充还是邻近点填充均可理解为将空白区域填充到相应的地形上。

针对 Ⅱ 类数据缺失利用区域边界，通过一维形态学滤波得到其边界各点地形，然后利用边界地形基于TIN对空白区域进行填充。地形填充既能统一不同原因引起的 Ⅱ 类数据缺失，又能保证填充区域内部和外部地形的连续性。地形填充的具体流程如下：

1)对掩膜影像进行形态学腐蚀运算得到腐蚀后的掩膜影像，将原掩膜影像与腐蚀后的掩膜影像相减即可得到数据缺失的边界。

2)对提取的边界进行编码连接。由于提取的边界点无拓扑关系，而一维形态学滤波是根据周边地形确定各点高程，故将相邻的边界点进行编码连接构成一维向量，使得相邻的边界点在一维向量中同样紧邻。

3)边界点一维形态学滤波。根据实验区域的先验知识，确定最大建筑物的尺寸。以最大建筑物的尺寸为滤波结构元的大小，依次得到各边界点的地形。图2所示为某一数据缺失区域提取的边界地形，一维形态学滤波较好地提取了边界各点的地形。

图2 某数据缺失区域边界地形提取示意图

4)边界点构TIN并进行内插填充。利用边界点的像素坐标构建不规则三角形，对空白区域的所有栅格进行判断其所属的三角形，然后利用不规则三角形顶点的地形高进行内插求得各空白栅格的地形。

5)依次对影像进行1)、2)、3)、4)步骤，直至所有的Ⅱ类数据缺失全部填充完成。

3 评价分析

为了充分对比和分析3种方法的填充效果以及其对滤波结果的影响，本文采用目视分析和定量分析相结合，定量分析采用ISPRS提供的精确分类的点云样本以及其评价体系[8]。

3.1 目视分析

ISPRS提供的FSite5中包含由水体吸收和航带无重叠度而造成的 Ⅱ 类数据缺失。图3中格网大小为2 m，图3(a)显示 Ⅱ 类数据缺失区域，图3(b)、图3(c)、图3(d)分别为最低点填充、邻近填充和地形填充的填充影像。区域1包含航带无重叠度而造成的 Ⅱ 类数据缺失，呈条带状长约1600 m，宽约70 m，且条带两侧高差较为明显，最低点填充与边界地形高程相差很大;邻近填充在数据缺失区域中线处存在较为明显的高程跳跃;而地形填充的效果较好，其高程的连续性以及与周边地形的吻合性均较为理想。区域2包含水体吸收而造成的 Ⅱ 类数据缺失，长约50 m，宽约20 m，最低点填充和地形填充能较好地对空白区域进行填充，填充部分与周边地形吻合较好;而邻近填充则将河流区域分块填充，填充视觉效果较差。

图3 数据缺失区域及其填充影像

综合上述，利用最低点填充对于水体效果最佳，但对由航带无重叠造成的数据缺失填充效果较差；邻近填充会在数据缺失区域中央产生高程不连续，同时还改变了地物的尺寸；地形填充法可以兼顾两者的特点，填充效果从整体而言优于前两种方法，其高程的连续性及与周边地形的吻合性都达到了最佳。

3.2 定量分析

在ISPRS提供的15个样本中有5个样本涉及到大面积数据缺失，分别为Sample21、Sample41、Sample51、Sample52、Sample61，其中Sample51、Sample52两个样本中包含由水体吸收和航带无重叠度而造成的数据缺失，其他样本均只包含航带无重叠度而造成的数据缺失。保证滤波参数一致，利用形态学滤波对3种填充方法获得的填充影像进行滤波，探讨不同填充方法对滤波结果的影响。

滤波质量通过三类误差进行定量评价，其中第1类误差(Type Ⅰ)是地面点被误分为地物点的百分比；第2类误差(TypeⅡ)是地物点被误分为地面点的百分比；第3类误差为总误差，即被错分的点占整个数据的百分比[8-9]。

表1为样本的三类误差统计，最低点填充TypeⅡ较小但TypeⅠ较大，主要是当边界存在较大高差时，其将边界周边内几乎所有的地面点和非地面点全部标记非地面点。Sample52 3种填充方法存在较大的差异，图4中(a)、(b)、(c)分别为其最低点填充、邻近填充以及地形填充的误差分布图。Sample52下半部由水体吸收造成的条带状的Ⅱ类数据缺失，其上部存在由航带无重叠造成的Ⅱ类数据缺失。从误差分布图可得出：3种填充方法在河流沿

岸无明显区别，在因航带无重叠度而造成的Ⅱ类数据缺失区域有明显的差异，主要是由于最低点填充和邻近填充产生的高程不连续造成的。其他样本地形填充与邻近填充相比，其TypeⅠ误差、TypeⅡ误差和总误差都保持一致。从所有样本汇总来看，地形填充效果最优，邻近填充次之，最低点填充较差。综上分析可得出：栅格数据的填充方法会对滤波精度产生影响，若空白区域两侧存在较大高差时，其影响更为明显。

表1 样本的三类误差统计 %

综上所述，不同填充方法将会对滤波的精度产生影响，填充方法主要在数据空白区域周边产生影响；3种填充方法对小范围水体吸收区域都较为有效；对航带无重叠造成的数据缺失，邻近填充和地形填充精度相当，但在空白区域两侧存在较大高差时地形填充效果优于邻近填充。

图4 Sample52不同填充方法误差分布

4 结束语

本文介绍了点云数据栅格化的基本原理，详细讨论了栅格化过程中存在的空白栅格及其处理方法，在对不同类型数据缺失的填充方法分析的基础上提出了基于地形填充法。利用ISRPS提供的测试数据，通过目视分析和定量分析对3种填充方法进行了测试和分析，得出：①最低点填充和邻近填充都会产生高程不连续现象，而地形填充效果从整体而言优于前两种方法，其高程的连续性及与周边地形的吻合性均较理想；②填充方法将影响滤波的精度，填充方法主要在数据空白区域周边产生影响；③地形填充能够较好地解决由水体吸收和航带无重叠等不同原因造成的数据缺失，可进一步提高滤波精度。

[1]PINGEL T, CLARKE K, MCBRIDE W. An improved simple morphological filter for the terrain classification of airborne LIDAR data [J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2013, 77: 21-30.

[2]沈晶，刘纪平，林祥国．用形态学重建方法进行机载LiDAR 数据滤波[J]．武汉大学学报:信息科学版，2011，36 (2)：167-176.

[3]李峰，崔希民，袁德宝，等．改进坡度的LiDAR 点云形态学滤波算法[J]．大地测量与地球动力学，2012，32 (5)：128-132.

[4]周晓明，马秋禾，许哓亮，等．LIDAR点云滤波算法分析—以ISPRS测试实验为参考 [J]．测绘工程，2011，20 (5)：36-39.

[5]CHEN Q, GONG P, BALDOCCHI D, et al. Filtering airborne laser scanning data with morphological methods[J]. Photogrammetric Engineering and Remote Sensing, 2007, 73(2): 175-184.

[6]罗伊萍，姜挺，王鑫,等．基于数学形态学的LiDAR 数据滤波新方法[J]．测绘通报，2011(3)：15-19.

[7]隋立春，张熠斌，柳艳，等．基于改进的数学形态学算法的LiDAR点云数据滤波[J]．测绘学报，2010，39(4)：390-396.

[8]SITHOLE G, VOSSELMAN G. Experimental comparison of filter algorithms for bare-earth extraction from airborne laser scanning point clouds[J]． ISPRS Journal of Photogrammetry and Remote Sensing，2004，59(1-2):85-101.

[9]程晓光，黄先锋，张帆.机载LiDAR数据的城区树木点提取方法[J].测绘科学，2014，39(3)：52-56.

[责任编辑：刘文霞]

Research on filling aggregated missing data of LiDAR with topographic method

GU Yan-chao, FAN Dong-ming

(School of Geoscience and Environment Engineering, Southwest Jiaotong University, Chengdu 611756, China)

Through studying on aggregated missing data of LiDAR in the procedures of rasterizing, a method which fills aggregated missing data with terrain of boundary is presented based on analyzing lowest method and nearest method. Terrain of boundary is obtained by one-dimensional morphological filter, and triangulate irregular network (TIN) is constructed. Interpolation is carried out based on Terrain of boundary and TIN. Experiments show that the proposed method adapts to aggregated missing data caused by different reasons to effectively ensure the continuity of the elevation, and improves the filtering accuracy.

airborne LiDAR; rasterizing; data filling; irregular distribution points; morphological filter

2013-08-18,2014-07-28补充更新

谷延超(1989-)，男，硕士研究生.

P237

：A

：1006-7949(2014)10-0023-04