基于车载三维激光雷达的玉米点云数据滤波算法

张 漫 苗艳龙 仇瑞承 季宇寒 李 寒 李民赞

(1.中国农业大学现代精细农业系统集成研究教育部重点实验室， 北京 100083;2.中国农业大学农业农村部农业信息获取技术重点实验室， 北京 100083)

0 引言

作物表型是特定的基因在一定环境条件下的外在表现，是内在基因和外部环境共同作用的结果。作物表型与基因信息之间的关系是作物育种的重要研究课题，随着现代作物基因技术飞速发展，传统的作物表型测量技术已经严重制约了育种技术的发展[1-2]。因此，使用立体相机、深度相机和激光雷达等传感器测量作物表型参数研究得到高度重视。数字植物可以实现农业生产的数字化和信息化，在株型设计和种植规划等方面具有广泛的应用前景[3]。

作物表型测量和数字植物面临的首要问题是使用传感器采集作物三维数据信息，并对数据进行滤波处理。目前获取作物三维数据的方法主要有使用立体相机、深度相机和高精度三维激光扫描仪等。在室内环境下，立体视觉技术可以快速、高效、准确地获取作物的表型信息[4-7]，但在田间应用时会因光照和风等环境因素的影响而无法工作。另外，表型信息采集的过程虽然是无损的，但将样本采集到实验室却是破坏性的，并且在运输过程中不能保证作物表型参数前后的一致性。

Kinect深度相机是一款价格便宜的体感设备，可以获取作物的深度信息，生成三维点云数据。仇瑞承等[8]使用Kinect相机采集玉米的彩色和深度信息，并测量玉米的株高。何东健等[9]通过Kinect相机采集玉米和茄子的三维点云数据，并提出了一种基于密度分析和深度数据双边滤波的方法。夏春华等[10]使用Kinect相机在室内采集了番茄的三维点云数据，并采用改进密度分析算法和双边滤波算法去除离群点。以上研究多以室内为主，存在受外界光照和风的影响较大和分辨率低等问题。

高精度三维激光扫描仪受外界光照的影响小，可以提供作物的三维结构信息，具有精度高、扫描速度快等优点[11-16]。地面三维激光雷达价格昂贵，点云数据量大，并且覆盖样本量少，无法在运动中采集数据。为解决这一问题，机载激光雷达被用于采集作物的三维点云信息。苏伟等[17]对机载雷达采集的数据进行TIN算法滤波和随机采样一致性算法分割，提取了植被树木和电力塔等飞行障碍物。针对航空三维点云数据的滤波算法，还有一些学者取得了一定的成果[18-22]。机载激光雷达测量范围广，但是测量精度低，稳定性不足。

车载多层三维激光雷达可以避免光照的影响，快速获得精度较高的三维点云数据，且测量范围广，同时设备价格适中，是作物表型测量和数字植物的研究热点。EHLERT等[23-24]使用车载4层三维激光雷达测量了冬小麦的株高。获取高精度的三维点云数据是提高表型测量的前提，在获取的原始点云数据基础上进行滤波处理，滤除地面点、离群点和设备自身原因产生的噪声，获得高精度的三维点云数据。滤波处理的结果直接影响之后的特征提取、匹配和三维重建等处理。目前对车载三维激光雷达运动过程中采集的点云数据滤波的研究相对较少。激光雷达点云属于线状点云，点云类型不同，在数据滤波及离群点数据去除方法等方面也有所不同。

本文以大喇叭口期的玉米为研究对象，使用车载16层三维激光雷达采集点云数据，开展点云数据滤波算法研究。在对三维点云数据进行预处理的基础上，针对不同的离群点，提出一种基于统计分析去除三维点云数据离群点的两次滤波算法。对非冠层残缺的叶片点云和冠层叶片边缘的噪声进行分步去除，为后续作物表型参数测量和植物数字化提供高精度三维点云数据。

1 材料与方法

1.1 试验材料

选用Velodyne LiDAR公司的VLP-16型三维激光雷达作为信息采集设备，VLP-16型三维激光雷达每秒可以采集30万个点，采集频率5～20 Hz，水平视角范围360°，分辨率0.1°～0.4°；垂直视角范围30°，分辨率2°；测量范围0.5～100 m。使用VeloView软件采集玉米生长期的三维点云数据。在Windows 7操作系统下，以Visual Studio 2013为平台，安装点云库PCL1.8.0(Point Cloud Library)和Cmake3.8.0，使用C++语言编程实现。

试验地点为北京市中国农业大学上庄实验站。试验时间为2017年7月13日09:00—11:00，玉米品种为农大84和京农科728，其中，农大84为舒展型品种，京农科728为紧凑型品种。行间距为90 cm，株间距为30 cm，生长期为大喇叭口期。VLP-16型三维激光雷达安装在车载测量平台顶端，如图1所示。垂直采集玉米三维点云数据，其与玉米冠层顶端的距离约为60 cm。车载测量平台以1 m/s左右的速度向前行进，VeloView软件同时采集三维点云数据，存储为.pcap格式，用于后续处理。本次试验共计采集了60株玉米。

1.2 三维点云数据处理

对采集的玉米三维点云数据首先进行预处理，主要是去除其他玉米植株和地面点云数据。预处理后的点云数据还存在非冠层残缺叶片的点云数据和冠层边缘的噪声。因VLP-16型激光雷达采集的数据在水平和垂直方向数据点的密度不同，常用的滤波算法难以适用于这种形式的数据。因此，根据点云及去除点的特征，提出一种基于统计分析去除离群点的两次滤波算法。

1.2.1三维点云预处理方法

使用VeloView软件在田间获取的三维点云数据为.pcap格式，且连续播放。点云库PCL1.8.0无法读取.pcap格式的数据，需要先将数据转换为.pcd格式数据，并且拆分为每帧都可以单独读取。在每帧三维点云数据中，非试验行玉米植株、测量平台和地面的点云数据占有比较大的比例，对后续处理无用，且影响处理时间。因此，采用简单快速的直通滤波法去除非试验行玉米植株、测量平台和地面的点云数据。在田间采集试验时，测量VLP-16型雷达的安装位置与其到地面的距离，通过这些数据确定玉米的三维空间范围，分别设置点云坐标X、Y、Z的阈值范围。进行预处理后，再进行后续滤波处理。

1.2.2精度评价标准

引用分类算法中的评价指标来优选滤波算法中k和α(k为临近点个数，α为标准差倍数)的值。对预处理之后的数据人工标定为离群点和保留点，即实际类别中的实际正例和实际负例。通过滤波算法判断的离群点和保留点为预测类别中的正例和负例。选用精确率(Precision)和召回率(Recall)来定量评价滤波的品质[25]，其计算公式为

P=TP/(TP+FP)

(1)

R=TP/(TP+FN)

(2)

式中P——精确率R——召回率

TP——正确划分为正例的个数

FP——错误划分为正例的个数

FN——错误划分为负例的个数

通过设置两次滤波的精确率和召回率阈值，缩小k和α的取值范围。在第1次滤波试验中，选取几个k和α组合，使用这些组合的保留点进行第2次滤波试验。最终选取3种组合：第1组为平均精确率与平均召回率和最大的k和α组合；第2组为平均召回率大于0.85，平均精确率最大的k和α组合；第3组为平均精确率大于0.85，平均召回率最大的k和α的组合。

1.2.3三维点云数据两次滤波算法

经过预处理的玉米三维点云数据仍然具有非冠层残缺叶片点云和叶片边缘离群点两类噪声。VLP-16型激光雷达采集的三维点云数据是由16根线状的雷达点云组成，在水平方向数据点非常密集，而在垂直方向却相当稀疏，只有16个数据点。不同于地面激光雷达和Kinect相机，这两个传感器采集的三维点云数据是阵列式的，水平和垂直方向数据点的密度相同。常用的基于三角网的滤波法、数学形态学的滤波算法和曲率统计的三维点云滤波算法，算法复杂，且很难适用线状稀疏三维点云数据。根据点云数据中存在的这两类噪声，提出一种基于统计分析的两次滤波算法，通过设置和修改统计分析算法的参数去除非冠层残缺叶片点云和叶片边缘的离群点。

统计分析算法[26-27]首先输入点云；设置点云中每个点的邻域大小，即邻域中临近点的个数k；计算每个点pi(xi,yi,zi) (i=1, 2,…,n)到它邻域中每一个临近点pj(xj,yj,zj) (j=1, 2,…,k)的平均距离di，计算式为

(3)

每个点的平均距离如果符合正态分布，全局邻域平均距离的均值μ和标准差σ计算公式为

(4)

(5)

式中n——整个点云个数

μ——邻域平均距离均值

σ——标准差

点云中平均距离di大于μ+ασ取值范围内的点定义为离群点。

根据玉米三维点云试验结果知临近点个数k和标准差倍数α决定离群点的数量和分布的空间位置。在离群点数量上，首先由α决定离群点的基数，k决定变化趋势，如图2所示。当k不变时，离群点的基数数量随着α的增大而减少。当α不变时，且α取值范围在[0.1,1.6]内，离群点的数量基本随着k的增大而增多，只有在k值为1～2或256～512时，部分离群点数量出现减少。当α不变时，且α取值范围在[1.7, 1.9]内，离群点数量则呈现先减少后增加又减少的变化趋势。在离群点分布的空间位置上，由k决定离群点分布的轮廓形状，α决定轮廓形状的宽度。当α固定时，k值较小时离群点零散地分布在点云的各个位置，随着k值增大，离群点会聚集在边缘小块的点云；k值越大，被认为是离群点的点云块越大，即离群点分布空间大、数量多；当k值固定时，α越小，轮廓的宽度越大，如图3所示。

图2 离群点数量与k和α关系Fig.2 Relationship of number of outliers and k and α

图3 离群点空间分布与k和α关系Fig.3 Relationship of spatial distribution of outliers and k and α

根据统计分析算法两个参数的不同取值会滤除不同的离群点，契合玉米三维点云数据中两类噪声的特点，可以选取不同的参数进行两次滤波。第1次滤波去除的噪声点是非冠层残缺叶片点云。这类点云呈块状分布，每块点云的数量在几十个。点云在块状内部彼此之间的距离很近，但是块状点云之间的距离较远，尤其是距离最大的块状点云(玉米冠层点云)距离远。因此从较大的k值中选取适合点云块的值，并寻找合适的α值决定轮廓的宽度，避免去噪不完整或过度去噪。第2次滤波去除的噪声点是叶片边缘离群点。这类点云零散分布在玉米冠层叶片点云的附近，它们与其最近的几个点云之间的距离较远，然而整体与叶片点云块距离较近。因而从较小的k值中选取合适的值，并通过比较选取合适的α值，去除叶片边缘的离群点。

基于统计分析的两次滤波算法流程如图4所示。

图4 基于统计分析的两次滤波算法流程图Fig.4 Flow chart of two times filtering algorithm based on statistical outlier removal

2 结果与分析

2.1 三维点云预处理

通过测量得VLP-16型激光雷达安装位置距离边缘车架大于0.5 m，雷达距离地面为2.9 m。设定点云坐标X、Y、Z的阈值分别为X[-0.5, 0.5]、Y[0.5, 2.0]、Z[-0.5, 0.5]。对两个玉米品种的各25株三维点云数据分别进行直通滤波预处理。如图5所示，直通滤波前有非试验行玉米植株、测量平台和地面的点云数据，滤波后这些非研究对象点云被很好地去除。点云的数量级从12 000降至1 700，减少了点云数量，节约后续处理的时间。

图5 玉米点云直通滤波处理Fig.5 Maize point cloud pass filter processing

2.2 基于统计分析的第1次滤波

对标定后的20株京农科728玉米点云数据进行一次滤波参数选取试验。在参数选取试验中需要根据精确率和召回率阈值人工选取结果，确定k和α的取值区间。根据基于统计滤波算法的经验，k初始值设置为50，每次增加0.6倍，上限为326；α初始值设置为0.1，每次累加0.1，上限为2；每一个k值都有20个α值与之对应。通过设置P>0.8，同时R>0.5，进行一次滤波数据处理，选取k的取值区间为[50,160]，α的取值区间为[0.8,1.3]。为了进一步缩小k和α的取值区间，k初始值设置为50，每次增加10，上限为160；α初始值设置为0.8，每次累加0.1，上限为1.3。设置P>0.9，同时R>0.6，再进行一次滤波数据处理，选取k的取值区间为[90,110]，α的取值区间为[0.8,1.0]。求得取值区间内所有k和α组合的平均精确率和平均召回率，如表1所示。

表1 第1次滤波平均精确率与平均召回率Tab.1 The first filtered average precision and recall

选取k和α组合进行第2次滤波，分别是平均精确率与平均召回率和最大的(90,0.8)，平均精确率最大的(90,1.0)，平均召回率最大的(110,0.8)；每行每列中平均精确率与平均召回率和最大的组合：(90,0.8)、(100,0.8)、(110,0.9)、(90,0.8)、(100,0.9)、(100,1.0)。去掉重复组合，共计7个k和α组合：(90,0.8)、(90,1.0)、(100,0.8)、(100,0.9)、(100,1.0)、(110,0.8)、(110,0.9)；保存这些组合的保留点数据依次为数据组1～7，每组数据都有20株玉米。记录这些组合下的离群点个数、离群点中人工标定为离群点的个数和召回率。其中k和α组合为(90,0.8)的滤波效果如图6所示。白色点云为第1次滤波保留点云，红色点云为去除的非冠层残缺的叶片点云，由非冠层叶片和残缺叶片组成，呈块状分布，与冠层距离较远，每块点云数量在几十个，主要位于非冠层部位。经过该次滤波，平均点云数量从1 700降至1 400，去除离群点300个。

图6 k和α组合为(90,0.8)的滤波效果Fig.6 Filtered effect of combination of k and α as (90,0.8)

2.3 基于统计分析的第2次滤波

对第1次滤波后的7组数据进行第2次滤波试验。以k和α选取(90,0.8)为例说明试验过程，其他6组进行相同的试验。

读入k和α选取(90,0.8)滤波后保存的数据，设置k的初始值为2，每次增加0.5倍，上限为50；设置α的初始值为0.1，每次累加0.1，上限为2。通过第1次滤波记录的数据，计算整个算法的精确率和召回率，并设置P>0.85，同时R>0.85；进行1次滤波处理，选取k的取值区间为[3,13]，α的取值区间为[0.9,1.5]。为了进一步缩小k和α的取值区间，k初始值设置为3，每次增加1，上限为13；α初始值设置为0.9，每次累加0.1，上限为1.5。设置P>0.87，同时R>0.87，再进行1次滤波数据处理，选取k的取值区间为[4,9]，α的取值区间为[1,1.4]。计算取值区间内所有k和α组合的平均精确率和平均召回率，找到平均精确率与平均召回率和的最大值，并记录该k和α组合和该最大值。然后扩大k和α的取值范围，查找平均召回率大于0.85，平均精确率最大的k和α组合；以及平均精确率大于0.85，平均召回率最大的k和α的组合。

经过数据处理得到7组数据的平均精确率与平均召回率和的最大值及对应的k和α的组合；平均召回率大于0.85，平均精确率最大值及对应的k和α组合；平均精确率大于0.85，平均召回率最大值及对应的k和α的组合，如表2所示。

表2 第2次滤波k和α组合及其平均精确率与平均召回率Tab.2 The second filtered k and α combination and average precision and recall

确定平均精确率与平均召回率和最大值，其两次滤波k和α组合为(110,0.9)、(6,1.2)；确定平均召回率大于0.85，平均精确率最大值，其两次滤波k和α组合为(100,1.0)、(6,1.2)；确定平均精确率大于0.85，平均召回率最大值，其两次滤波k和α组合为(110,0.8)、(5,0.9)。两次滤波组合为(110,0.9)、(6,1.2)的第2次滤波效果如图7所示。白色点云为第2次滤波保留点云，红色点云为去除的冠层叶片边缘离群点，呈线状分布，在叶片边缘部位，数量在几个。经过该次滤波，平均点云数量从1 400降至1 300，去除离群点100个。

图7 组合为(110,0.9)、(6,1.2)的第2次滤波效果图Fig.7 The second filtering effect with combination of (110,0.9) and (6,1.2)

2.4 试验验证

验证数据集包括3部分。一部分是参与滤波参数选取的6株玉米的其他帧点云数据，为验证集A；其他两部分是没有参与滤波参数选取的数据，5株同一生长期京农科728玉米点云数据，为验证集B；5株同一生长期农大84玉米点云数据，为验证集C。验证集数据的两次滤波处理是自动完成的。

使用选取的两次滤波k和α的组合(110,0.9)、(6,1.2)，处理验证集数据分别得到其精确率和召回率，如表3所示。从表中数据可以计算出验证集A的平均精确率为0.937、平均召回率为0.930；验证集B的平均精确率为0.927、平均召回率为0.906；验证集C的平均精确率为0.904、平均召回率为0.947。可以看出，参与滤波参数选取的验证集A的平均精确率和平均召回率之和略高于没有参与滤波参数选取的验证集B和C；不过验证集B和C的平均精确率和平均召回率都大于0.9，满足滤波要求。且验证集C的平均精确率与平均召回率的和大于验证集B，证明该两次滤波参数组合在同一玉米品种内可以通用，使用紧凑型玉米品种选取的该参数与舒展型玉米品种可以通用。

表3 两次滤波k和α为(110,0.9)、(6,1.2)组合时验证集 数据的精确率和召回率Tab.3 Precision and recall of verification set data when two times filter k and α were combination of (110, 0.9) and (6, 1.2)

使用选取的两次滤波k和α的组合(100,1.0)、(6,1.2)，处理验证集数据分别得到其精确率和召回率，如表4所示。从表中数据可以计算出验证集A的平均精确率为0.949、平均召回率为0.851；验证集B的平均精确率为0.934、平均召回率为0.872；验证集C的平均精确率为0.910、平均召回率为0.924。可以看出验证集A和B都可以确保滤波的平均召回率在0.85附近，并且平均精确率更高；不过验证集C的平均精确率和召回率都大于0.90，偏离该两次滤波k和α组合的要求。所以该组合在同一玉米品种的同一生长期可以通用，使用紧凑型玉米品种选取的该参数与舒展型玉米品种不可以通用。

表4 两次滤波k和α为(100,1.0)、(6,1.2)组合时 验证集数据的精确率和召回率Tab.4 Precision and recall of verification set data when two times filter k and α were combination of (100, 1.0) and (6, 1.2)

使用选取的两次滤波k和α的组合(110,0.8)、(5,0.9)，处理验证集数据分别得到其精确率和召回率，如表5所示。从表中数据可以计算出验证集A的平均精确率为0.863、平均召回率为0.960；验证集B的平均精确率为0.847、平均召回率为0.933；验证集C的平均精确率为0.822、平均召回率为0.974。可以看出验证集A和B的平均精确率都在0.85附近，并且平均召回率更高；但是验证集C平均精确率偏离0.85较大，没有满足该两次滤波k和α组合的要求。所以该组合在同一玉米品种的同一生长期可以通用，使用紧凑型玉米品种选取的该参数与舒展型玉米品种不可以通用。

表5 两次滤波k和α为(110,0.8)、(5,0.9)组合时 验证集数据的精确率和召回率Tab.5 Precision and recall of verification set data when two times filter k and α were combination of (110, 0.8) and (5, 0.9)

两次滤波k和α的3种组合在不同品种、同一生长期的通用性有不同表现。原因可能是京农科728是紧凑型植株，农大84是舒展型植株。采集的玉米三维点云数据在空间分布上有相同点，也有明显差别，相同点是京农科728和农大84点云数据都主要分布在玉米冠层；差别是京农科728冠层数据占整个点云的比例低于农大84。因为有点云数据都主要分布在冠层的相同点，所以在平均精确率与平均召回率和最大的组合可以通用；然而冠层点云占整个点云的比例不同，会导致京农科728寻找到的边际参数组合不适用于农大84。

3 讨论

本文仅针对一个生长期的玉米三维点云数据进行滤波研究，玉米点云在点个数和点密度变化很小，使用固定的k和α组合可以满足基本要求。在后续研究中，为了进一步提高滤波的精确率、召回率和鲁棒性，在玉米的同一生长期，应根据点个数和点密度的微小变化对k和α值动态微调；不同生长期的玉米点云数据的点个数和点密度变化较大，k和α的组合应根据不同的生长期进行动态变化。

经过滤波处理后的数据，成功去除了两类离群点，这对后续进行的特征提取、点云配准、三维重建和可视化有着重要作用。未经滤波处理的点云直接进行特征提取，会导致点云匹配失败，进而无法获得多帧点云匹配后的数据。多帧点云匹配的数据与单帧点云数据相比，具有三维结构信息更丰富、精度更高等优点，这对于提高玉米株高、叶长、叶宽等表型参数的测量精度具有重要作用。

4 结论

(1)提出一种基于车载三维激光雷达采集的玉米点云数据的两次滤波算法，并求出平均精确率与平均召回率和最大的k和α组合为(110,0.9)、(6,1.2)；平均召回率大于0.85，平均精确率最大的k和α组合为(100,1.0)、(6,1.2)；平均精确率大于0.85，平均召回率最大的k和α的组合为(110,0.8)、(5,0.9)。

(2)验证试验发现，平均精确率与平均召回率和最大的k和α组合可以在京农科728和农大84之间通用，其他两组k和α的组合在京农科728和农大84之间不通用。