基于动态K阈值的苹果叶片点云聚类与生长参数提取

刘 刚 张伟洁 郭彩玲

(1.中国农业大学现代精细农业系统集成研究教育部重点实验室， 北京 100083；2.中国农业大学农业农村部农业信息获取技术重点实验室， 北京 100083； 3.唐山学院机电工程系， 唐山 063000)

0 引言

叶冠是果树光合作用的主要场所，其空间形态结构直接影响冠内的光照分布、光合生产力及枝叶器官的生理状态，进而影响果实品质和产量[1-2]。叶片作为叶冠的主要组成部分，获取其空间分布规律及生长参数对获取精确的叶冠结构具有重要的研究价值，也是进一步研究冠层光照分布计算的基础[3]，可为果树的自动化修形剪枝提供技术支持[4-6]。

叶片的识别方法多局限于二维空间中，虽然存在运算速度快等优点[7-8]，但难以进一步提取生长参数。地面三维激光扫描仪在农业领域的应用[9-12]使得快速获取果树冠层三维空间信息成为可能，进而可实现冠层叶片的自动分割及生长参数的提取，能在较大程度上减少人力，提高作业效率[13-15]。

国内外研究人员对点云分割方法进行了大量研究，主要分为6类[16]：基于区域边界的分割方法适用于边界明显、分布规则的数据[17-18]；基于区域增长的分割方法适用于有明显分布差异的点云[19]；基于图的分割方法和基于模型匹配的分割方法对规则点云有很好的聚类效果[20]；基于机器学习的分割方法对训练样本的要求较高[21]；基于聚类的分割方法采用特征聚类，根据不同的点云分布特征，选择适当的特征检测方法[22-23]。果树叶片空间分布相对密集，枝叶连续生长难以区分，因此利用特征聚类实现叶片分割较为合理。

叶片生长参数提取研究中，吴升等[24]通过计算叶片点云法向量与XOY坐标平面的夹角，实现点云在XOY坐标平面的正投影，进而将叶片几何参数求解降维到二维空间上；郑一力等[25]采用主成分分析与线性评判分析相结合的方法，实现了叶片的多特征降维计算。降维法计算叶片生长参数简单快速，但忽略了叶片的三维空间特性。针对以上问题，本文提出一种叶片点云聚类分割算法，融合局部凹凸性算法(Locally convex connected patches，LCCP)并改进K-means算法，动态获取K值，进而研究叶片的生长参数提取方法。

1 材料与方法

1.1 点云数据获取方法

随机选取中国农业大学苹果树采摘机器人试验基地(北京市昌平区南口镇辛力庄村)的1棵富士苹果树作为研究对象，树高3.1 m，树龄4年，采用美国Trimble 公司TX8型地面三维激光扫描仪(具体参数如表1所示)进行点云扫描。由于试验果树周围还生长有其他果树，遮挡较多，为了获得高质量的点云，数据获取试验采用5站扫描，如图1所示。其中的s1～s5为5个扫描点位，扫描仪的工作位置与被测果树树干的距离为3.5～5 m，c1～c6为用于配准定位的6个标靶球(标准直径100 mm)，分布于试验果树周围，保证在每站扫描能识别3个以上标靶球，以便多站点云的配准，扫描实地场景如图2a所示。本试验不考虑环境(如风速、温湿度、光照度等)的影响。

表1 Trimble TX8型扫描仪产品参数Tab.1 Product parameters of Trimble TX8

图1 扫描仪布置示意图Fig.1 Scanner layout diagram

图2 点云数据采集Fig.2 Point cloud acquisition

利用基于标靶球提取的算法[26]配准多站点云，点云数据获取过程中由于外界环境干扰和设备自身的误差，产生大量的噪点，采用前期研究成果[27]——基于距离均值计算的方法，剔除噪点。点云密度大且分布不均匀，采用均匀网格法完成抽稀处理，最终所得的苹果树冠层点云数据如图2b所示，图2a中白色小球为标靶球。

1.2 点云处理方法

本研究算法由C++语言编写实现，点云处理库为点云库(Point cloud library，PCL)。获取的点云数据以三维坐标(x，y，z)的形式存储，数据量庞大，约有115万点。为了简化果树叶片点云的搜索和存储，创建新的叶片点云PCD存储格式，每条记录存储一个叶片节点，包括的属性有叶片中心点Li(xi,yi,zi)、叶片法向量Ni、叶长Lleni和叶宽Lwidi。

本文研究方法流程主要分为4部分(图3)：数据获取、点云预处理、叶片点云聚类分割和叶片生长参数提取。数据获取试验需在晴朗无风条件下进行；点云预处理主要包括点云配准、去噪和抽稀；采用基于动态K阈值的算法，融合LCCP算法和改进的K-means算法，随机选取冠层的一个枝条，作为研究对象，进行叶片点云聚类分割；叶片生长参数提取方法采用PCA方法和边界提取方法，分别实现叶片生长角度和叶长叶宽参数值的获取。

图3 叶片点云分割及生长参数提取流程图Fig.3 Flow chart of apple tree point cloud segmentation and growth parameters extraction

1.2.1叶片点云的聚类方法

叶片点云聚类分割，采用基于动态K阈值的算法，融合LCCP算法，并采用改进的K-means算法，实现K阈值的动态获取。该算法主要包含以下过程：超体聚类、LCCP聚类[28]、K阈值获取和K-means分割。

超体聚类是一种过分割方法，过分割涉及的参数有体素距离Rv、晶核距离Rs。八叉树初始化点云，过分割生成独立的超体素块pi=(xi,ni,Ei)，其中pi为超体素索引、xi为体素中心、ni为体素法向量、Ei为邻接体素边集合，计算邻接体素相似度完成超体聚类

(1)

式中D——体素合并的概率

wc、ws、wn——颜色、空间、法向量的信息权重

Dc、Ds、Dn——颜色、空间、法向量容差

以超体聚类所得体素集合为研究对象，建立体素索引向量，采用CC(Extended convexity criterion)和SC(Sanity criterion)方法判定邻接体素凹凸性关系完成LCCP聚类。

(2)

(3)

β=∠(n1,n2)

(4)

式中 CCb(pi,pj)——邻接体素块凹凸性

n1、n2——邻接体素块法向量

βT——法向量角度偏差阈值

i、j——体素索引号

当邻接体素面不连接时，CC的判定准则不适用，此时需要使用SC准则对邻接体素关系进行判定，SC判定依据为

(5)

其中θ(pi,pj)=min(∠(d,s),∠(d,-s))=

min(∠(d,s),180°-∠(d,s))

(6)

(7)

s=n1×n2

(8)

式中 SCb(pi,pj)——邻接体素凹凸性

conv(pi,pj)=CCb(pi,pj)∧SCb(pi,pj)

(9)

式中 conv(pi,pj)——pi、pj的凹凸性

动态K阈值的获取主要分为两个步骤：一方面计算K=4时，K-means聚类后各点集三维质心间欧氏距离；另一方面计算叶片中心点到拟合直线的距离均值，作为距离阈值。动态获取K阈值的具体方法如下：

(1)根据LCCP所得类簇，设定K=4，进行K-means聚类分割，并计算4个点群的中心点Ckm(xm,ym,zm)。

(2)计算4个中心点之间的6组欧氏距离(其中m、n为类簇中心点编号)

(10)

(3)手动截取枝条叶片点云，计算所有叶片点云中心点Cs(xs,ys,zs)，根据所得点云中心点拟合直线L，计算叶片中心点到直线L的距离均值，作为距离阈值Dt，此阈值为统计结果，不需要重复计算。

(4)将步骤(2)所得的中心点欧氏距离，与距离阈值Dt进行比较，大于、小于距离阈值的个数分别为a、b，根据a与b的大小关系，确定此点群需要分割的点群个数K，其中K值的选定准则为

K=4(a=6)a+12(b

(11)

式中 ⎣·」——向下取整

计算每个类簇的K值，再次进行K-means聚类，得到单叶片和枝干点云。

1.2.2叶片的生长参数提取方法

将叶片点云三维空间信息作为输入参数，利用PCA法构建协方差矩阵，计算叶片点云的法向量，即计算由叶片点云三维坐标组成的协方差矩阵的特征值和特征向量。构建的协方差矩阵为

(12)

Avj=λjvj(j=0,1,2)

(13)

式中n0——叶片点云总数

Pi——叶片点坐标

λj——第j个特征值

vj——λj对应的特征向量

计算所得的最大特征值max_eigenvalue对应的特征向量max_eigenvector()作为叶片法平面方向向量，可根据计算所得结果，进一步计算叶方位角、叶倾角等参数。

图4为叶片生长参数求解示意图，根据前期研究成果，提取叶片边界点云，计算边界点中点与叶片中心点的最大距离，作为叶片的叶长。如图4b所示，N为计算所得叶片生长方向向量，O为三维质心，l1和l2为叶片长轴端点，与l1、O、l2所在平面近似垂直的最大欧氏距离为叶宽，w1和w2为叶片宽轴端点。

图4 叶片点云边界提取及参数求解示意图Fig.4 Sketch of leaf points boundary extraction and parameter solution

图5 传统方法聚类效果图Fig.5 Sketches of segmentation by traditional methods

2 结果与分析

2.1 叶片点云聚类效果与精度

叶片的聚类分割方法主要考虑了叶片点云的分布特征。叶片与枝干连接且分布密集，导致叶片与叶片相互区分困难、叶片与枝干相互分割困难，若采用传统的聚类算法则会有较多不足，图5采用3种传统的点云聚类方法，对同一单枝点云数据进行处理，图中可见，采用DBSCAN算法无法将连续的枝、叶点云分割开；另外在点云获取过程中由于遮挡等原因，会导致局部点云缺失，若仅利用K-means算法或Region-growing算法，会出现大量的单叶片被分割成两个到多个聚类的情况，导致所得单叶片点云不完整。

如上分析，需要找到一种既能成功分割单叶片点云，也能保证其完整性的方法。如图6所示，采用超体聚类和LCCP算法可以顾及点云局部凹凸性，将整体枝条点云分割为多簇，即将整体长枝条分割为较小的具有多个完整叶片的点云类簇，如图6b中的黄色立方体框为点云类簇的包围盒，此步骤能较好地保证叶片的完整性，避免一个叶片被分割到不同的类簇中。

图6 超体聚类与LCCP聚类效果图Fig.6 Sketches of super-voxel and LCCP clustering

图7 动态K阈值分割结果Fig.7 Sketches of dynamic-K-threshold segmentation

超体聚类中涉及的参数设置主要有Rv=10.0,Rs=0.3，由于点云无颜色信息，故ws=1.0,wn=10.0,Ds=1.0,Dn=1.0,wc=0,Dc=0。由图6可看出，叶片分割最大的难点是将叶片点云与连接的枝干点云进行分割。LCCP聚类结果能够保证叶片点云的完整性，但不能直接将叶片、枝干分割开来，而K-means算法将距离作为相似性度量的标准，叶片、枝干点云间具有明显的区域分布性，根据不同的LCCP聚类结果，设定不同的K值完成K-means聚类，即将叶片点云进行有效分割。因此，设定K-means算法中的K=4，计算所得4组聚类中两两的中心点距离，与距离阈值进行比较，根据动态K阈值判定准则可以动态获取K值，再次利用K-means方法聚类，得到枝叶分离的点云聚类结果。主要包含叶-叶(图7b)、枝-叶(图7c)、叶-枝-叶(图7d)以及叶-枝-叶-叶(图7e)4种情况。图8a为枝条叶片最终分割结果，图8b为最终苹果树冠层叶片分割结果。

图7为不同K阈值所得分割结果，图7a为动态K=1的分割结果，只包含1个叶片的点云数据；图7b、7c为K=2的分割结果，图7d为K=3的分割结果，图7e为K=4的分割结果，图7所包含的结果只代表不同的动态K值的情况，不包含所有的叶片分割结果。识别叶片和枝干点云，需要借助点云可视化的方法。为了验证本文提出的叶片分割方法的有效性，将本文所得冠层叶片分割结果与多种方法所得结果进行比较，如表2所示。

分析表2数据可得，DBSCAN方法所得聚类叶片数量较少，是由于单个聚类点集中往往会包含1个以上叶片。相反地，K-means方法和Region-growing方法所得叶片相对较多，是由于点云中若存在少量空洞或局部缺失则会导致数据不连续，算法执行时会将此区域聚类为独立的点集。本文提出的叶片聚类方法，叶片识别的正确率达到90%以上，另外，动态K阈值的获取方法也减少了人工干预，较大提高了运算效率。

2.2 生长参数提取效果与精度

叶片点云作为输入参数，利用PCA法构建协方差矩阵并计算叶片点云法向量，作为叶片生长方向向量；利用边界提取方法获取叶片边界点云，计算边界点与叶片中心点的欧氏距离，从而获取叶长、叶宽的参数值。为进一步验证本文提出的生长参数提取方法的可行性，随机选取10组叶片，每组中包含20个叶片，定义平均叶长、叶宽作为评价指标，分别记录利用实测方法、投影法、本文提出的方法所得结果，如表3所示。可以看出，利用本文提出的生长参数提取方法，叶片生长角度误差小于5°的叶片数量占叶片总数量的90%以上，叶长、叶宽参数值的平均准确率均达到91%以上，相较于降维法所得结果，本文提出的方法能够进一步提高提取精度。

图8 叶片点云分割效果图Fig.8 Sketches of leaf points segmentation

聚类方法叶片分割数量/个正确率/%本文方法110290.77K-means139485.18DBSCAN100782.95Region-growing143881.55实测值1214

3 结论

(1)根据苹果树叶片点云分布特征，提出了基于动态K阈值的叶片聚类方法，实现了冠层枝叶点云的分离，进而完成了叶片生长参数的提取。基于动态K阈值的算法，克服了传统算法中人工干预过多的缺点，可以有效提高生长参数的提取精度。

表3 叶片生长参数提取的精度对比Tab.3 Comparison of accuracy of leaf growth parameters extraction

(2)将本文方法与传统聚类方法的试验效果进行对比，结果显示，本文方法能够正确识别出的叶片数占冠层总叶片数的90%以上；本文提出的叶片生长参数提取方法所得结果，生长角度误差小于5°的在90%以上，叶长、叶宽平均准确率均在91%以上，与投影法提取准确率比较表明，生长参数的提取精度得到提高。

(3)本文提取的叶片生长参数种类较少，未实现枝叶点云的自动化分类识别，在后续的研究中，需要借助可视化的方法对此进行完善和实现。