机载LiDAR数据加工DEM产品技术研究

2018-10-11张辉

福建地质 2018年3期

张辉

(福建省地质测绘院，福州，350011)

机载激光雷达(又称LiDAR)测绘技术是集激光扫描、全球定位系统(GPS)和惯性导航系统(INS)三大技术于一体，继GPS全球定位系统后一次重大技术革命[1]。目前已渗透国土资源、林业、测绘、遥感领域中广泛应用。

机载LiDAR技术能够快速、准确地获取地物三维坐标数据，生成高精度、高密度的数字高程模型，将复杂的地形(地貌)精准并形象地展现出来，具有传统航空摄影测量和地面常规测量技术无法取代的优势，被广泛应用于民用工程测量[2]。笔者以福建连江可门港经济开发区的机载LiDAR数据加工实验为例，研究机载LiDAR数据滤波处理工艺流程及关键技术。

1 LiDAR数据获取

机载LiDAR通过航空平台搭载的激光扫描设备进行空对地扫描探测，实现激光对地表或地物沿航线方向纵向扫描及通过转镜的横向扫描；同步GPS定位系统提供飞机精确方位数据和INS惯性导航系统提供的飞机姿态数据(航向、侧滚、俯仰和加速度)，从而获得一个大范围带状区域内的地物点云数据集。

1.1 原始POS数据获取

数据采集前，首先制定详细的飞行计划，完成测区资料的搜集；根据测区的情况，确定扫描的区域及飞机航线分区与敷设；根据项目数据比例尺精度要求，确定使用设备型号、飞行的高度及速度、地面与机载GPS接收机性能匹配、地面GPS基站位置情况等。测量前，需在2个以上GPS基准站上架设高精度GPS信号接收机，应与机载POS设备同步记录。测量时，开启基准站接收机，POS设备记录数分钟的静态GPS数据，进行测量初始化和IMU姿态初始化；然后进入测区按预定航线进行扫描。

1.2 POS数据解算

将原始的POS数据，通过POSPAC软件平台的Extract模块提取出GPS数据、IMU数据和辅助传感器数据。将GPS采集数据与地面GPS基准站同步的数据进行差分拟合计算，得到曝光瞬间摄影机投影中心精确的GPS定位坐标。然后将GPS定位坐标数据、IMU数据和辅助传感器数据加载到POSProc模块进行解算，得到6个方位元素。最后设置大气改正、角度改正、扫描仪改正、定位定姿偏移、高程偏差、回波强度值改正等，以及与外业测量时仪器状态和客观环境相关的测量区域相符合的椭圆参数、中央子午线投影等参数。最后解算生成标准LAS格式点文件[3]。标准LAS格式点云文件记录包含航带编号、XY坐标、高程值、回波次数、回点序列号、回波强度等6项信息。POS数据解算流程如(图1)所示。

图1 POS数据解算流程Fig.1 POS data processing diagram

1.3 数据精度与误差

在机载LiDAR扫描所生成的点云数据中，影响点云坐标测量精度主要因素是激光测距仪误差、差分GPS定位误差、IMU姿态测量误差和扫描角误差[4]。这些误差的存在导致POS系统对地定位结果产生误差，从而影响成果的精度。必须制定严密的飞行计划，严格仪器检校、合理敷设航线与地面基站点，达到最佳飞行方案，以此保障获取数据的精度。

1.4 机载LiDAR信息处理流程

对点云数据进行加工生产时，通过点云分类滤波算法研究，实现有效点云数据的提取。机载LiDAR数据信息处理流程可分为原始LiDAR数据获取与POS处理工序、LiDAR数据滤波分类处理工序及高精度测绘产品生产工艺。

2 点云数据滤波分类研究实验

2.1 实验区地貌特征

实验区位于福建连江可门港经济开发区，北纬 26°18′00″～26°27′00″，东经119°36′00″～119°43′15″，区内主要为大片滩涂、平地和丘陵，海拔0～400 m。东北面为罗源湾，分布大片浅滩涂地，西南面为高山地，202 省道、高铁路线南北贯穿测区，交通状况一般。

2.2 LiDAR数据质量评述

选择的实验区面积约131.96 km2，采用机载LiDAR系统航空遥感数据，获取时间2008年3～4月。相机幅面4 k×4 k彩色CCD数码相机，相机焦距55 mm；航摄1∶12 000；航速240 km/h，相对航高为1 200 m；相机像元地面分辨率为0.20 m，航向重叠度≥60%；旁向重叠度≥30%，像片倾斜角均小于2°，飞行旋偏角≤6°；点云点距≤2.0m，扫描角40°，脉冲频率22 kHz、衰减系数0.7。航飞获取实验区点云集数据总量为5 171万点，点云密度0.37～0.48点/m2符合 1∶2 000密度1点/m2的要求，数据精度符合要求。

LiDAR源数据分区。源数据基本分块编号：block01-80～block27-80共计27块；成果按分区(片)处理，分区编号：A～H区共计8区块；LiDAR源数据量：A区775万点、B区865万点、C区749万点、D区504万点、E区666万点、F区431万点、G区749万点、H区431万点。

2.3 点云空间信息

实验区点云空间基本描述地形表面、建筑物、林地、植被、路网等构建物面的测量点精确坐标信息，水域(水库、池塘、滩涂)无激光雷达反射波，多为点云空白区。依相对高程高度(相对于地形表面)可划分为地面点(0 m)、低植被(0.15～1.0 m)、中植被(1.0～2.5 m)、高植被(2.5～30.0 m)、建筑物(2.5～50.0 m)、高空中物(30.0～80.0 m)等不同点云高度层。上述针对实验区点云数据集的属性、密度及空间分布规律分析，结合实验区地形(地貌)特点，可作为制定“滤波分类方案”的重要科学依据。实验区点云密度大，点云数据总量超过5 000万，海量数据严重约束“滤波分类”计算能力。要解决海量点云处理技术瓶颈，只有采用合理分区设计方案，分区分片进行数据滤波处理。

2.4 实验区 LiDAR数据处理方案

基于TerraSolid平台对LiDAR点云数据的滤波分类技术方法，研究LiDAR数据加工DEM产品工艺技术。点云数据中包含真实的地面点集、植被点集(树木、草丛)、建筑物点集(房屋、桥梁、电力线等)及由各种噪声造成的异常点集。若要从庞大多要素、复合点云数据集中提取有效地面(或地表)高程点，并最终加工生产DEM测绘产品，则必须研究如何从这些点云中分离出地面点与非地面点，再从非地面点中分解出植被点、建筑物点、空中悬挂点和其他非地面点。LiDAR点云数据处理设计流程如(图2)所示。

图2 LiDAR 点云数据处理设计流程图Fig.2 Flow chart of LiDAR point cloud data processing design

点云滤波处理基本步骤：将点云噪声点滤除，使点云数据达到合理空间分布状态。现实中地表、地物表面激光雷达反射形成点云聚合体很复杂，无法依托一种算法将所有的地表与植被、地物(各类构造物)的点集提取出来[5]。因而，LiDAR数据处理滤波分类方法根据测区地形(地貌)条件，以及地表与植被、构建物之间依存空间关系，制定一系列滤波方案，逐级分解点云数据集，采用多组算法联合方可实现最终分类目的。根据不同的地物地貌选择不同的滤波、分类算法，试验不同的参数实现地物类别的有效提取。具体的滤波分类算法如(图3)所示。

图3 LiDAR点云数据分类流程图Fig.3 Flow chartt of LiDAR point cloud data classification

严格按照“机载激光雷达数据处理技术规范”[6]作相关点云分类生产设计，首先考虑机载 LiDAR 数据应用用途与目的。基于实验区机载 LiDAR 数据加工生产 1∶1 000～1∶5 000数字高程模型DEM产品重点是地面点、近地表点的点云信息提取，而低、中、高植被与建筑物点云顺带提取(表1)。

2.5 噪声点滤除

原始LiDAR数据包含噪声造成的异常点，偏差较大。如果直接进行分类处理，会因为某一个异常点的值造成整批数据分类失败。噪声异常点应先行滤除，使点云集处于较合理的空间位置状态下，以利后续数据批量自动处理[7]。

噪声点数据滤波算法相对较容易实现。根据实验区地形(地貌)海拔变化规律，实验区海拔变化0～400 m，因而0～400 m范围以外的点数据可判别为异常点。有些异常点在整个区域内并未超出最大、最小高程，则此时可用相对高差阈值法来处理。点云数据可理解为带坐标的空间点数据集，根据一定的方向(如X轴)成连续、离散的排列；设计一个高差的阈值，比对该点云与周边的点云相对高差来判断其是否为异常点。

表1 LIDAR点数据类型

2.6 地面点提取

原始LiDAR数据集首先进行地面点与非地面点2大类滤波提取，再从非地面点云集中将某种地物类别的点集分离出来，即为分类算法。分类算法即利用基于反射强度、回波次数、地物形状等的算法或算法组合对点云数据进行批量自动分类[8]。裸露地表处只有一次回波对应的发射点，此次回波对应的反射点即为地面点如实验区丘陵、山地区域植被覆盖区域可能对应多次回波；正常的地面点是最后一次回波对应的发射点，如丘陵、山地区域。

渐近加密地面点提取算法主要原理是选取一定数据窗口内较低点作为初始地面点，基于格网窗低洼点建立初始TIN三角网模型，反复构建地表TIN三角网模型方式，迭代将邻近点加入到地面点，直到满足迭代终止条件为止[9]。算法在开始时选择一些低点，认为它们是位于地表处，通过“Max building size”参数来控制初始点的选择。如果建筑物的最大边长是60 m，应用程序认为每隔60 m至少存在一个位于地表处的点，也就意味着该点就位于地表处。初始模型三角网构TIN大多数低于地面，只有最高点接触到地表。算法通过“Iteration angle”“Iteration distace”筛选控制参数，反复加入新的地面点开始向上扩建模型，每个加入的点使模型更加贴近地表面。

2.7 非地面点滤波分类

根据点的高度及点云分布的形状、密度、坡度等特征，非地面点云包括植被、交通路网及附设施、水体及附设施、构筑物等要素点云数据识别与分类提取。非地面点中滤波分类主要依托地面高程点三角构网模型为参考框架，设计低、中、高植被点的相对高差阀值，即可将植被点分离出来。由于算法的适用性会将建筑物点合并在植被点集中，对已分离的植被点集应用数学形态学算法处理，可将建筑物点二次分离出来。

2.8 交互式人工辅助分类

交互式人工辅助分类是自动分类方法的补充，通过交互剖面法、坡度法、影像信息等可视环境，对错误分类的点云进行重新分类或更正分类，以提高分类的准确率[10]。一些精度要求较高的项目，还要对LiDAR点云数据进行人工交互式的分类检查。在芬兰TerraSolid-SCAD软件平台下，可通过剖面法、坡度法、影像色彩法等方式直观地判别点云数据分类的效果，采用人工交互式的分类二次处理，将异常点或错分点分离(图4)。

图4 TerraSolid软件平台分离各类地物点前后对比图Fig.4 Comparison of before and after separating all kinds of object points on TerraSolid software platform

3 专题测绘平台数字测绘成果加工

3.1 基于 GeoTIN平台DEM格网计算处理

GeoTIN模块具海量数据处理能力。将LiDAR点云数据的地面点集提取为X、Y、Z点坐标的明码格式数据网格化计算，便可直接被GeoTIN平台导入[11]。在GeoTIN平台中读取离散的点集数据，加载地形地貌特征点(山峰、低谷)、特征线(山脊、山沟谷)为约束条件，经过TIN三角网构建，设定比例尺格网尺寸后，进行高精度网格化数值计算，即可生成DEM数据产品。

3.2 基于Surfer平台DEM产品生成

在Surfer平台下导入LiDAR地面点云数据,由离散点高程计算网格化点数字高程模型DEM。实验区选择1 m×1 m网格尺寸，通过Surfer平台由DEM数据再现三维地形地貌晕渲图。Surfer平台中提供TIN三角构网法、加权反距离法、克里格、最近邻近点、多项式回归方法等12种将离散点网格化的算法，不同网格化算法对DEM最终产品具有不同效应(图5)。

图5 TIN三角构网法生成的地形地貌DEM渲染效果图Fig.5 DEM rendering rendering of terrain and landform generated by TIN triangulation method

3.3 基于ArcGIS平台DEM产品生成

ArcGIS由ArcMap、ArcCatalog、ArcScene及ArcToolbox 4大功能程序组成。①DEM网格化算法：ArcToolbox平台中转换工具(Conversion Tools)提供离散点到栅格数据处理能力，ASCII to Raster菜单驱动可实现LAS格式点云数据网格计算，并生成DEM产品数据。TIN三角构网网格化算法具有计算速度快，数值精度高特点。②可视化立体晕渲图：采用自主研发LDS-System系统处理，导出ArcGIS平台数据，按X、Y、Z坐标值排列。在ArcGIS平台Add XY Date模块，选取导入的数据文件后，自动读取X、Y、Z三列的数据值，生成点图层。在构建三角网处理后，在ArcGIS Layers对TIN的三角构网进行渲染显示，可根据所构建的TIN数据生成DEM数据(图6)。

图6 ARCGIS平台DEM渲染图及DEM效果图Fig.6 Diagrams of DEM rendering and DEM effect on ARCGIS platform

3.4 DEM产品质量检验

以实验区DEM成果数据为例，对DEM产品进行检验。①XY平面坐标精度检验通过GeoTIN 3.16软件平台逐幅回放DEM等高线图结合1∶1 000 DLG、DOM叠加检验，通过观测道路、山峰等地貌特征部位比对，进行平面精度检测。②高程精度检验依据实验区控制点、大比例尺地面离散高程点资料，基于GeoTIN 3.16软件平台质量检查程序功能，按离散点法进行高程精度的检查，抽样统计中误差平地、山丘<0.64 m，山地<1.7 m，符合“CH/T 8023-2011机载激光雷达数据处理技术规范”要求(表2)。

表2实验区DEM产品高程精度检查统计

Table2StatisticaltableofhighprecisioninspectionofDEMproductsinexperimentationarea

图幅号地面检查点>3 rms样品Block14-805368(0.5%)Block15-8045821(1.1%)Block16-80621(1.6%)Block17-8030524(2.7%)Block19-807641(0.1%)Block20-8090313(3.3%)Block21-8070724(1.4%)

注：据1∶1 000地形图高程点资料。