LiDAR点云生成DEM的水面置平方法研究与实现

张 弓同，李四海，焦红波，曲 辉

(国家海洋信息中心，天津300171)

一、引 言

LiDAR(激光雷达测量)是美国航天局(NASA)1970年开发的三维激光扫描测量技术，它通过位置、距离、角度等观测直接获取对象表面点三维坐标，直接获取真实地表的高精度三维信息［1］，已广泛应用于道路工程、海岸线海岛礁、特殊地形测绘、数字城市等领域，被认为是测量领域的又一次技术革命。

由于激光光束能够接收多次回波，通过一定的数据后处理算法，可以有效地剔除植被的影响，获得高精度数字高程模型(DEM)［2］，与数字正射影像图(DOM)叠加，可以快速获得研究区域地面的三维信息。利用LiDAR产生的DEM，在三维城市建模、铁路勘察设计、地形景观的显现等领域有很好的应用。在灾害调查与环境监测中，可以迅速获得灾害现场的情况，对制定灾害的预防及补救方法、灾后评估等有着重要的意义［3］。

但是，由于LiDAR光束无法穿透水体，导致水体内部的激光点数量很少，因此在拟合数字高程模型等产品时，水面会出现空洞或崎岖不平的现象，影响了数字高程模型等产品的使用效果。本文就水面置平问题展开研究。

目前，关于水面置平有不同的方法和技术。要解决这个问题，方法之一就是水体信息的提取。在国外，Brzank等分析了水体在LiDAR数据中所呈现出的特征，通过监督分类区分水体和湿泥地，提取水体和结构线［4］。在国内，王宗跃等采用双层格网模式，基于LiDAR点云数据提取较窄的水体，并以朝向水体的边界点作为拟合轮廓线的关键点提取更精确的轮廓线，取得较好的效果［5］;张永军等提出将LiDAR数据对水体的敏感性与航空影像的高分辨率特征相结合的水体自动提取方法，并用数学形态学的方法进行边缘信息区域生长运算，最终获得水体区域［6］。

本文以LiDAR激光点云数据为基础，使用相邻点旋转角最小原理，对水体轮廓信息进行有效识别之后，再对水面进行合理的填充和置平处理，有效提高了数字高程模型等产品的应用效果。在数据处理中，笔者使用C#语言编程序实现。

二、水面置平方法的建立

1．LiDAR数据的格式

在对LiDAR数据进行地面点提取时，首先要分析LiDAR数据的格式［7］，实现读写。

2．算法过程

首先提取LiDAR数据末次滤波的地面点，对LiDAR数据空白处考虑为水面点，提取基本的(x，y，z)坐标，生成ArcGIS文件，按照平面边界点的旋转角最小原理，提取水面和岛屿边界，按照填补水面坐标点的方法转成DEM，达到水面置平的目的。本文的岛屿指的是水面里的点集，可能是真实的岛屿或水里的地物。具体流程如图1所示。

(1)提取数据整个边界

对整个数据的边界范围进行提取，生成ArcGIS文件，对边界进行编辑修改，得到自己感兴趣部分的范围。

(2)数据存储方法

首先对激光点云分块处理，本文采用等格网分割方法［8］对海量点云数据进行分块处理。

图1 算法流程

(3)水面和岛屿边界提取方法

由于水面为空数据点，因此适当调整网格宽度在水面部分就会产生一个空数据网格。将所有为空数据网格组成一个集合，从集合中的第一个空网格开始，依照网格向左查找，直到网格内部不为空为止，如图2所示。

在找到这个网格k后，取网格k内最右边的点a0，从这一点开始，以角度为0°、旋转半径为r，顺时针旋转，遇到第一个点即为下一点a1，其旋转角θ为最小值，如图3所示。这是得到第一个点的方法。

图3 第二个点的选择方法

以后点的选取如图4所示，以点am为圆心，按照am到am－1方向开始，依据选好的半径r，沿着顺时针的方向，遇到第一个点即为下一点am+1，其旋转角θ为最小值。依次类推，得到一系列的点，每个水面构成一个点集，得到一个点集合，称这个点集合就是一个水面或岛屿。

图4 第m+1个点的选择方法

(4)具体分析

对于当前得到的am点可能有以下几种情况:

1)噪声点:当找到一系列点后，当前点am在半径为r的圆范围内查找点时发现，除了am－1外无法找到其他的点，或都是已经找到的点，则称am点为噪声点，将am点删除或另行处理。使用图4所示方法可得到新的am点，依次找到水面的所有点。

2)水面点:按照步骤(3)方法得到水面范围，将水面范围内(不包括边界)所有的点提取出来构成一个集合list1，然后从list1中找到坐标点x方向坐标最大的点b0，从b0开始，依次考虑是噪声还是岛屿，一直做到list1中的点全部解决完为止。然后，对查找到的噪声点和岛屿点分别处理。

3)岛屿点:使用步骤(3)算法，很容易得到岛屿的外围。在得到岛屿点的集合并保存后，将岛屿内部及边界所有的点从网格内部删除，这样就不影响岛屿所在的水面的第一个边界点的选取，直到将所有的水面都画出来为止。

4)边界点:使用步骤(3)算法，如果遇到整个数据边界，就要边界加点依次画出水面的外围。

这样就生成了水面和岛屿的ArcGIS文件，如图5所示。

三、水面置平的技术实现

图5中的ArcGIS文件是参考Esri发布的Shapefile Technical Description格式，只保存水面和岛屿点的x、y、z坐标(如图6所示)。

图5 含有噪声点、岛屿点、水面点和边界点生成的面文件

图6 生成的水面岛屿文件

画出水面和岛屿之后，实现水面内部的加点，要求所加点的z值与水面边界点及其内部点的最小z值一致。由于LiDAR数据在水面部分会有少量的点返回，因此这样得到的z值即为水面的高度。添加的点生成新的ArcGIS点文件。水面部分在输出到DEM时水面部分z值相同，而且为水面所围范围的高度，保证了水面部分是平的，也就是达到水面置平的目的。

将原始LiDAR数据生成的ArcGIS文件(如图7所示)与水面加点后生成的新的ArcGIS文件(如图8所示)相比，很容易看出水面密集的加点效果。

图7 原始的LiDAR数据

图8 水面均匀加点后的数据

四、LiDAR数据水面置平试验

原始LiDAR数据构建TIN时水面崎岖不平的状况很明显，如图9所示。经过水面置平技术处理后构建TIN时水面部分已经平坦，而且加氧机也显现出来了，说明水面置平效果很明显，如图10所示。

图9 原始LiDAR数据生成的TIN

图10 LiDAR数据处理后生成的TIN

原始LiDAR数据转成DEM时水面有空洞现象，如图11所示;水面置平处理后LiDAR数据转成DEM时水面平坦，水面内部的加氧机很明显，如图12所示。可以从区域的遥感影像中看出实际的地形情况，可以看出水面置平效果很显著，如图13所示。

图11原始LiDAR数据生成的DEM(有空洞)

图12 LiDAR数据处理后生成的DEM

图13 研究区域的遥感影像图

五、结束语

本文使用的LiDAR数据为2006年12月21日海南省老高原地区沿海数据，中心经纬度为108．97°E、18．395°N，数据大小约为2．7 MB，范围约为0．54 km×0．53 km，在此范围内的水面展开研究，建立了一种LiDAR数据提取水面范围的方法，通过将LiDAR数据资料应用到上述方法中，所得到的置平结果比没有处理过的结果有明显改进，表明了该水面置平方法的有效性。

［1］ 周淑芳．基于机载雷达LIDAR与航空像片的单木树高提取研究［D］．哈尔滨:东北林业大学，2007．

［2］ 张良，马洪超，邬建伟．联合机载LIDAR数据和潮汐数据自动提取潮位线［J］．遥感学报，2012，16(2):411-416．

［3］ 赖旭东．机载激光雷达基础原理与应用［M］．北京:电子工业出版社，2010．

［4］ BRZANK A，HEIPKE C，GOEPFERT J，et a1．Aspects of Generating Precise Digital Terrain Models in the Wadden Sea from LiDAR-water Classification and Structure Line Extraction［J］．ISPRS Journal of Photogrammetry and Remote Sensing，2008(5):510-528．

［5］ 王宗跃，马洪超，徐宏根，等．基于LIDAR点云数据的水体轮廓线提取方法研究［J］．武汉大学学报:信息科学版，2010，35(4):432-435．

［6］ 张永军，吴磊，林立文，等．基于LIDAR数据和航空影像的水体自动提取［J］．武汉大学学报:信息科学版，2010，35(8):936-940．

［7］ 刘春，姚银银，吴杭彬．机载激光扫描(LIDAR)标准数据格式(LAS)的分析与数据提取［J］．遥感信息，2009(4):38-42．

［8］ 刘昌军，丁留谦，孙东亚，等．海量激光点云的分块处理及植被自动过滤技术研究［C］∥第一届全国激光雷达对地观测高级技术研讨会论文集．北京:［s．n．］，2010．