基于机载LiDAR属性信息检测离群点云

李沛婷，赵庆展，田文忠，马永建，陈洪

(1石河子大学信息科学与技术学院/国家遥感中心新疆兵团分部/兵团空间信息工程技术研究中心,新疆 石河子 832003; 2 石河子大学机械电气工程学院,新疆 石河子 832003)

机载激光雷达(light detection and ranging,LiDAR)因其有非接触主动式和可以快速获取物体高精度点云等优点,已成为获取点云的一种新颖方式,其在智慧城市、资源调查和基础测绘等领域发挥着越来越重要的作用[1]。机载LiDAR集成了激光扫描仪、动态差分全球定位系统(global positioning system,GPS)、惯性导航系统(inertial navigation system,INS)数码相机等先进技术[2],其中定位定姿系统(positioning and orientation system,POS)控制确定飞行平台的空间位置和姿态,通过激光扫描仪获取光学中心与地物目标之间的距离信息,利用数码相机捕捉图像信息[3]。因为机载LiDAR采集的点云具有数据量大、冗余高、点密度不均匀等特点,所以研究点云数据处理是后期将其应用在不同领域的基础[4],其中点云处理包括预处理和后处理,预处理主要是点云滤波,目的是较准确分离地面点与非地面点,且利用地面点云生成数字高程模型(digital elevation model,DEM)。DEM是地理信息系统发展的基础数据,其在很多应用领域上都是重要的基本需求和地形产品[5],如生成荒漠区高精度DEM,是开展荒漠区植被、地形变换和生态系统分布等研究的重要基础。王智等[6]结合亚洲数字高程数据和新疆土地利用数据,实现对新疆山地-绿洲-荒漠中植被覆盖变化时空特征研究;李向婷等[7]结合机载LiDAR测量得到的DEM数据和MODIS影像数据,提取新疆荒漠稀疏植被覆盖度。

目前,国内外学者一方面直接使用点云的三维坐标属性,或者通过三维坐标获取点云法向量等间接特征实现滤波,另一方面也逐渐通过对点云回波次数和强度的深入研究来提高点云滤波精度。Pirotti F等[8]结合地面激光扫描仪的回波次数,实现点云滤波得到DEM;Reymann等[9]针对不同试验区的点云数据,得到结合点云回波次数和回波强度属性信息以实现点云分类,比仅使用点云的几何信息更加精确;Ghamisi P等[10]采用改进的支持向量机,结合点云高度变化、回波强度、激光扫描回波和拟合的残差实现点云分类;龚亮等[11]和曾静静等[12]利用点云强度、回波次数和高度三种属性信息,提取机载LiDAR点云的道路,结合属性信息可提高点云分类精度;刘凯斯[13]根据首末次回波高度值去除植被点,不仅提升点云的滤波效率,而且以回波次数、高度值和回波强度三个属性信息作为地物分类的约束条件构建决策树实现点云分类;樊敬敬[14]提出了直接利用首尾两次回波强度之差提取植被信息,结果表明该方法提取植被的效果较好,避免了使用两次回波高度差提取过程中的植被和建筑物边缘混合问题。

由上述研究可见,结合点云回波次数、回波强度和点云高度等属性信息,可以快速高效地实现点云滤波和点云地物分类等处理。在实际应用中,如果直接使用原始点云的回波属性进行滤波处理,将会大大降低运算速度和效率,因此，本文在进行点云滤波之前,首先采用箱形图对点云属性的离群情况进行分析,在去除大量离群点云的基础上,再通过反复建立三角网来实现点云滤波。

1 材料与方法

1.1 无人机平台

以瑞士AeroScout公司提供的Scout B1-100单旋翼油动无人直升机作为飞行平台,其空机质量50 kg,有效载荷18 kg,对应的长宽高为3.3 m×1.0 m×1.3 m,飞行时间可达90 min。

POS系统是由Oxford Technical Solutions(OXTS)公司提供的Survey+2,其采用GNSS / INS紧密耦合的方法实现飞行器空间位置和姿态的确定,刷新频率为100 Hz,定位精度为2 cm,横滚和俯仰精度为0.03°。

以Riegl公司提供的VUX-1全波段激光雷达为飞行器,激光扫描仪的扫描方式为线性,波长是近红外,扫描速度范围10～200 scan/s,最大扫描角度为330°,激光脉冲频率550 kHz,测量精度15 mm,不仅可以记录16 bit的回波强度,而且可以记录无限次的回波次数。

在数据获取过程中,因为激光扫描仪VUX-1和Survey+2是固定连接的,所以点云精度和POS系统关系紧密。

1.2 数据获取

1.2.1 研究区概况

研究区位于新疆第八师150团5连周边的荒漠植被区(85°58′35′E～85°59′4′E,44°57′29′N～44°58′0′N),主要以旱生和沙生灌木为主,如梭梭、骆驼蓬、麻黄、碱蓬草和驼绒藜等。于2017年7月31日,用无人机搭载传感器在荒漠植被区上空飞行采集数据,飞行速度5 m/s,飞行高度60 m,共3条航线,每条航线扫描面积86 400 m2。

机载LiDAR系统可获取目标对象的丰富属性信息,其数据形式称为点云[15]。本文研究以WGS84为地理坐标系统,以UTM_Zone_44 N为投影坐标系获取点云数据,以LAS格式存储数据[16]。在获取LiDAR数据过程中,由于外界环境的影响导致系统容易产生噪声点云[17],这些噪声点多为孤立点,表现为明显高于地物或者低于地表的点云[18],所以本文研究在采用激光扫描仪配套软件调整点云数据的扫描姿态、方位和粗差等处理后,再手动去除大量噪声点云。

1.2.2 点云可视化

截取研究区部分点云作为本文原始数据,其包含研究区主要地物类型,每个点的属性主要包括点云三维坐标、强度、回波次数、扫描角度、GPS时间等信息[19]。

为方便后期描述,本文将点云三维坐标、回波强度和回波次数分别命名为[X,Y,Z,I,Re]。截取区域的覆盖面积约为80 m×77 m,点云密度为128.224 pts/m2,点云三维坐标范围为{X=[20.042 5,100.697 5],Y=[80.507 5,157.875 0],Z=[275.290 0,284.039 2]},其中X和Y分别表示东和北方向的位置,Z表示点云高度,单位为m。回波强度范围为{I=[5 767,65 535]},扫描视场角度范围为{Angle=[-22,33]}。

图1是以点云高度和回波强度属性信息可视化原始激光雷达点云数据。

图1 试验区点云高度和回波强度图Fig.1 Elevation and echo intensity maps of original point clouds in test area

不同激光扫描仪系统,可以记录不同的回波次数,包括单次回波和多次回波。在多次回波中,接受第1个和最后一个回波信号的分别被称为首次和末次回波,除了首次和末次回波外的回波信号为中间次数的回波[20]。

图2a是以全部回波次数可视化试验区点云,其中深棕色表示第1次回波点云,深蓝表示第2次回波点云,绿色表示第3次回波点云。分析图2a可知:

第1次回波对应的点云个数为777 915,占点云总数98.49%,点云高度范围为{Z=[275.29,284.039 2]}。

第2次回波获取的点云个数为11 581,占点云总数1.49%,高度范围为{Z=[275.35,280.17]}。

第3次回波获取点云个数为339,占点云总数0.04%,点云高度范围为{Z=[275.49,278.44]}。

第4次和第5次对应的点云个数分别为20和4,高度范围分别为{Z=[275.81,278.21]}和{Z=[275.64,278.11]}。

上述分析表明,第1次与第2次回波对应点云的高度范围相差很小,但第1次与第3、4、5次回波的点云高度范围相差较大,且目视判别第3、4、5次回波对应的点云,发现其均为噪声点云。

对第2次回波点云的强度进行可视化,结果如图2b所示,图中不同颜色代表不同强度值。从图2b可知,第2次回波主要包括地面点云、电线的阴影点云和微量低矮植被点云。因为本文最后通过对比点云滤波结果分析检测离群点云结果,即只保留地面点云,所以本文后期仅使用箱形图方法分析第1次回波对应的点云,再合并检测结果和第2次回波点云,作为点云滤波处理的输入数据集。

图2 试验区点云回波次数可视化Fig.2 Showed by point clouds echo times of test area

1.3 检测离群点云

1.3.1 技术路线

利用Matlab R2017软件编程提取点云三维坐标、回波强度和回波次数,在分析每次回波对应点云的分布情况后,保留第1次和第2次回波对应的点云。首先，采用箱形图方法对第1次回波点云的高度属性进行离群分析,得到高度离群点云数据集D1；再在分离点云数据集D1的基础上,对剩余点云进行强度属性离群分析,得到强度离群点云数据集D2、分离高度及强度离群点云后对应的点云数据集D3；然后,合并数据集D3和第2次回波对应的点云,得到点云数据集D4,并结合软件ENVI 5.3对数据集D4进行滤波；最后,与直接对原始点云滤波结果对比,其技术流程如图3所示。

图3 结合点云属性和箱形图检测离群点云技术流程图Fig.3 Flowchart of detecting of point clouds′ outliers based on combining point clouds′ attribute information and box-plot method

1.3.2 箱形图原理

检测离群点云的主要目的是识别不同地物类别,常用方法主要是基于统计、邻近度、密度和聚类来实现检测[21]，但是,在对单变量数据集进行离群探测分析时,箱形图因其简单、直观的优点,成为最常选择的方法之一[22]。箱形图是描述最小值、第一四分位数Q1、中位数、第三四分位数Q3、最大值和异常值的一种统计方法,其中心位置为中位数,箱子长度表示四分位数的间距IQR,离群值是根据异常值截距线和极端值截距线得到,其中异常值截距线和极端值截距线的定义如式(1)所示。

IQR=Q3-Q1，

B1=Q3+1.5×IQR,B2=Q1-1.5×IQR，

BE1=Q3+3×IQR,BE2=Q1-3×IQR。

(1)

式(1)中，B1和B2为异常值截距线对应的值,BE1和BE2为极端值截距线对应的值。

本文离群点云是由属性大于B1和小于B2,或者大于BE1和小于BE2的点云组成。

2 结果与分析

2.1 箱形图分析高度属性

利用回波次数去除噪声点云后,采用箱形图对第1次回波对应的点云高度进行离群统计分析,得到如图4所示的点云高度离群统计结果,其中纵坐标表示点云高度值。

图4 第1次回波点云对应的高度箱形图及高度离群点云可视化Fig.4 Box-plot and outlier point clouds of elevation values corresponding to the first echo point cloud data

图4a为第1次回波次数对应的高度箱形图,经过分析得出:点云高度Z为278.976 7是整个研究区高度临界点,且将高度范围在{Z=[278.976 7,284.039 2]}之间的点云归为高度离群点,对高度离群点云进行三维可视化得到灰度图4b,令该数据集为D1,其对应的回波强度范围为{I=[5 767,65 535]}。将高度范围在{Z=[275.290 0,278.976 7]}之间的点云归为分离高度离群后的点云,回波强度范围为{I=[7 143,58 369]},此时第一四分位数对应的点云高度为276.628 2,即第1次回波中包含25%的点云高度的范围为275.290 0～276.628 2;第三四分位数对应的点云高度为277.567 9,第1次回波包含75%的点云高度的范围为277.567 9～278.976 7。

2.2 箱形图分析强度属性

通过上述对点云高度属性数据进行箱形图分析,可以去除大量电线、车辆和微量高植被地物点云数据集D1，之后在分离高度离群点云的基础上,对剩余点云进行强度数据离群分析,得到分离高度离群点云后对应的强度离群点云数据集D2,再去除点云数据集D1和D2得到剩余点云数据为D3。

不同数据集可视化结果如图5所示。

图5 分离第1次回波高度离群后对应的点云强度分析结果可视化Fig.5 Visualization of analyzing point clouds echo intensity after separating elevation values outlier point clouds from first echo

图5a是分离高度离群点后点云强度箱形图,纵坐标表示点云回波强度值。由图5a可知,此时强度临界值是31 325,故将回波强度范围{I=[7 143,31 325]}归为分离高度离群后对应的强度离群点云数据集D2,点云高度范围为{Z=[275.490 5,278.972 0]};将回波强度范围{I=[31 325,58 369]}归为分离高度和强度离群后对应的点云数据集D3,其点云个数为766 380,点云高度范围为{Z=[275.290 0,278.976 8]}。另外，第一四分位数和第三四分位数分别对应的点云强度值是41 571和48 408,即分离高度离群后对应25%的点云强度的范围为7 143～41 571,对应75%的点云强度的范围为48 408～58 369。

图5b为强度离群点云数据集D2三维可视化结果,从图5b可知，点云数据集D2主要包括大量植被点云、微量道路点云和微量车辆点云。

图5c为分离高度和强度离群后的点云数据集D3的可视化结果,主要包含了研究区大量的地面点云。为了进一步提高点云滤波精度,得到研究区全部地面点云,将点云数据集D3和第2次回波点云合并得到数据集D4。

2.3 点云滤波

在软件ENVI 5.3 ENVI LiDAR模块的支持下,通过反复建立三角网分别对点云数据集D4和原始点云进行滤波,图6是滤波后地面点云可视化结果。

图6a是对数据集D4滤波得到的283 852个地面点云,点云高度范围为{Z=[275.291 5,278.312 7]},回波强度范围为{I=[8 322,58 369]},对应滤波消耗的时间为14.1 s。

图6b是对原始点云滤波获取的281 301个地面点云,点云高度范围为{Z=[275.294 5,278.363 0]},回波强度范围为{I=[7 470,57 976]},对应滤波消耗的时间为46.3 s。

对比两次滤波的结果(图6a、b)可知,基于点云回波次数、高度和强度属性信息,采用箱形图方法对原始激光雷达点云进行滤波前处理,有效保留了地面点云2551个,滤波消耗时间降低了32.2 s，而且提高了点云滤波速度。

图6 点云数据集D4和试验区原始点云滤波结果可视化Fig.6 Filtering of point cloud datasets D4 and original point cloud data in test area respectively

为了分析对比两次地面点云是否存在差异,本文将上述地面点云导入SPSS软件中,对点云高度进行统计分析,得到点云数据集D4与原始点云分别滤波后对应地面点云高度中位数、众数、标准差、方差、偏度相差和偏度标准误差都为276.900、276.512、0.610、0.372、0.003和0.005,即对2种点云数据集进行滤波得到的地面点云无明显差异,因此,本文采用简单易实现的箱形图方法对点云属性进行滤波前处理是可行且有效的。

3 结论

本文在点云滤波前,结合点云属性信息,利用箱形图对原始激光雷达点云进行滤波前预处理,通过对比处理前后的滤波结果,得到以下结论:

(1)结合点云回波次数,可以有效剔除多次回波点云363个。

(2)结合箱形图先后分析点云高度和点云强度的离群情况,且利用反复建立的三角网实现滤波,可获取地面点云766 380个,滤波消耗时间减少了32.2 s，表明通过本文方法可提高点云滤波速度。

(3)本文采用简单易实现的箱形图方法对点云属性进行滤波前处理是可行且有效的。