基于三维点云的大豆植株器官分割及表型分析

肖奕同， 刘帅， 侯晨连， 刘琦， 李富忠， 张吴平

（山西农业大学软件学院，山西 太谷 030801）

植物表型是指在基因和环境相互作用下产生的能够直观反映植物生长状态的特征及性状[1]。植物表型参数的分析与育种息息相关[2]，然而表型参数的测量结果受人员、环境和设备等因素影响较大，极大限制了大规模遗传育种筛选的效率[3]。因此，开发高效、准确的表型参数测量方法对提高植物表型测量效率、促进植物育种至关重要。基于二维图像的植物表型测量方法因设备要求低、易于使用等优势逐渐应用于农业工程领域[4-6]。但是，二维图像在表示成簇植株的冠层时有较大缺陷[7]，冠层结构中被遮挡的部位在图像中完全缺失。三维点云不仅具有色彩信息，同时弥补了二维图像深度信息不足的缺陷，在很大程度上避免了因叶片粘连、遮挡等问题而难以准确测量表型的问题。史蒲娟等[8]基于从不同角度拍摄的图像序列，对油菜植株进行三维重建，提取了株高、叶柄长及叶长等表型参数。李杨先等[9]采集不同视角的二维图像并结合运动恢复结构算法实现簸箕柳植株稀疏点云的重建，提取了株高、基径、叶片面积、分枝数和分枝角等表型参数。近年来，随着对植物表型分析需求的日益增多及三维点云表型研究的深入，对表型参数测量的精准度也有了更高的要求[10]，精准的表型分析需要依靠植株器官水平或具体部位的信息[11-13]。已有研究均是基于完整的植株点云数据进行表型分析，没有实现植株器官尺度的分割。朱超等[14]基于点云骨架和最优传输距离实现了玉米点云茎叶分割，分割算法平均总体准确率达到0.92。阳旭等[15]基于棉花植株的多时序点云数据，采用基于随机采样一致的区域生长算法完成冠层提取及分割。然而，相较于玉米及棉花植株，大豆植株叶片成簇且叶片间相互粘连现象严重，通过骨架提取[14]及柱体拟合[15]的方法均难以快速实现大豆植株茎叶分割。对植株点云进行高效、精确的器官分割既是表型分析的前提也是植物表型研究的重大难题[16]，为了解决多分枝作物叶片相互粘连带来的器官分割及高通量表型测量困难问题，本研究基于大豆植株的整体稠密三维点云，利用点云法向量、曲率差异，结合基于欧式聚类的区域生长算法实现大豆植株器官分割，并提取叶面积、叶宽、叶长、叶倾角、茎粗等表型参数，为多分枝作物茎叶分割、单叶分割以及表型参数的测量提供技术支持。

1 材料与方法

1.1 数据获取

本研究选取10 株分枝期盆栽大豆为研究对象。将植株置于光照均匀且无风的场景，利用RGB 相机（GoPro，HERO8 Black）以半球形状环绕大豆植株进行多视角拍摄（图1）。根据株型差异，采集30~50 张多视角图像组成图像序列。图像序列中相邻两幅图像间相互重叠区域达到70%以上。最后，分别测量各植株花盆口外径。

图1 相机位置参照Fig. 1 Camera position reference

1.2 点云重建与滤波

基于SFM（structure from motion）+MVS（multi view stereo）方法，利用Agisoft Photoscan 软件重构大豆植株稀疏点云（图2A）及稠密点云（图2B）。基于人工测量所得的真实花盆口外径，利用软件中的标尺工具实现植株点云尺度变换。

图2 重建的大豆植株三维点云Fig. 2 Reconstructed 3D point clouds of soybean plant

在获取大豆植株点云数据时，由于设备精度、环境因素等影响，点云数据中不可避免地出现一些离植株点云主体较远的离群点和噪声点[17]。除了这些随机误差产生的噪声点之外，受作物本身生长限制，如叶片交叉、叶片共面、叶片重叠等因素的影响，点云数据中往往存在一些贴近植株边缘的冗余点。为了更好地进行器官分割、表型分析等后续处理，基于PCL（point cloud library）[18]对植株点云进行滤波处理，包括以下几个步骤：①直通滤波过滤掉地面及花盆；②颜色滤波剔除与植株本身颜色差异较大的点，例如土壤及叶片边缘的杂色冗余点；③统计滤波剔除离群点；④降采样降低点云密度，减少计算时间。经过滤波处理后的大豆植株点云如图3所示。

图3 滤波处理后的大豆植株点云Fig. 3 Point clouds of soybean plant after filtering

1.3 器官分割

1.3.1茎叶分割 大豆植株叶片相对扁平，同一点在合理尺度范围的法线方向变化不大，但其茎秆近似于圆柱形，不同尺度法向量会产生较大差异，若法线差大于设定阈值，则认为该点属于茎秆，反之该点属于叶片部分。因此，可以利用不同尺度下法向量特征的差异性提取出大豆植株茎秆部分。利用法线微分差异（difference of normals，DoN）[19]实现茎叶分割，具体步骤为：①计算输入点云数据的平均距离（Rm），作为较大支持半径（Rl）以及较小支持半径（Rs）的取值依据；②对于输入点云数据中的每个点（P），分别利用较大的支持半径（Rl）与较小支持半径（Rs）计算其法向量，根据式（1）进行差分运算得到该点差分法向量；③过滤向量域，分割出大豆植株茎秆点云；④基于欧式距离聚类，合并阈值内的冠层点云。部分叶片上的点受成像精度的影响以及阈值的选取可能被错误地归为茎秆（图4C）而被移除。同时，一些孤立的茎部，特别是靠近叶片的部分，可能被错误地归类在叶片点云（图4B）中。为了提高后期表型参数测量的精确度，利用欧式聚类移除多余的茎部点云（图4D）并且回填漏掉的叶片点云（图4E），得到较为完整的冠层结构（图4F）。

图4 大豆植株茎叶分割Fig. 4 Stem and leaf segmentation of soybean plant

式中，P表示点云数据表示大半径及大半径下的法向量表示小半径及小半径下的法向量，∆n̂表示差分法向量。

1.3.2单叶分割 由于叶片点云密度不均匀且存在遮挡、粘连等情况，去除茎秆的大豆植株冠层叶片点云仍需进一步分割才能实现高效准确的表型参数分析。单叶分割的难点在于如何识别并断开叶片粘连部位。对于每个点，根据式（2）计算其协方差矩阵(C)的特征向量和特征值，根据式（3）进一步运算估计得到曲率（δ）。

式中，Pi表示点云数据，k表示Pi邻近点的数目表示邻近点的三维质心，T表示点云矩阵的转置，λj是协方差矩阵的第j个特征值是第j个特征向量；λ0、λ1、λ2为每个点云数据协方差矩阵的3个特征值。

叶片粘连部位的点分布相对混乱且具有较大的曲率值（图5A），剔除曲率值大于某个阈值的点，粘连叶片能顺利地相互分离。去除叶片粘连点可能会导致叶片边缘结构缺损，从而影响分割之后得到的单片叶片（图5B）。因此，本研究利用欧式聚类回填部分粘连点，综合考虑曲率、法向量夹角以及点间距离，采用改进的区域生长算法实现冠层叶片分割，得到较完整的单片叶片（图5C）。

图5 大豆植株冠层叶片分割Fig. 5 Leaves segmentation of soybean canopy

1.4 表型参数提取

1.4.1叶面积 基于拉普拉斯平滑（laplace smoothing）的德洛内三角剖分（delaunay triangulation）法（简称LDT 法）计算大豆叶片叶面积。拉普拉斯平滑是一种网格平滑算法，其原理是在一定的松弛度（r）下将每个点用其邻域点的中心来代替（式4），通过不断迭代，得到较为光滑的网格。随着平滑过程中迭代次数及松弛系数增大，平滑程度越大，网格会不断向中心收缩。为了避免过度收缩造成的误差，本研究在r=0.2、iteration=20的情况下对叶片进行平滑处理。平滑后利用三角剖分将散乱的点云数据剖分为一系列三角网格（图6E），通过计算所有三角网格面积之和（式5）来代替叶面积。

图6 大豆表型参数测量Fig. 6 Measurement of soybean phenotypic parameters

式中，x′表示点的新坐标位置，x表示点的当前坐标位置，j表示当前点的邻域点个数，xi表示第i个邻域点的坐标位置，S表示叶片点云面积，n表示三角网格个数，Si表示单个三角网格面积。

1.4.2叶宽 分割后的大豆叶片仍处于植株坐标系下（图6A），此时通过每个坐标轴最大点和最小点建立起来的包围盒称为轴对齐包围盒（axisaligned bounding box，AABB），其空间冗余大，对叶片的表型特征描述误差就大。为了更加贴近叶片真实的形状，基于主成分分析构建叶片有向包围盒（oriented bounding box，OBB）（图6B），提取包围盒的宽代表叶片的宽度（图6C）。具体步骤为：①利用主成分分析获得叶片点云的质心及3 个主方向，将点云转换至叶片坐标系下；②基于获得的主方向和质心，将输入的叶片点云转换至原点，建立变换到原点的叶片点云包围盒；③计算y方向上最大点和最小点间的距离。

1.4.3叶长 本研究采用提取点集以拟合叶片中脉的方法，能够提取叶片点云上近似于中脉的最短曲线（图6F）。具体步骤为：①在叶片坐标系下，找到叶片点云空间中距离最远的2 个点并分别定义为起点和终点；②将起点作为基点，基于最邻近（K-NearestNeighbor，KNN）算法搜索其邻近点，找出邻近点中距基点与终点距离之和最小的点作为下个基点；③重复搜索，直到基点与终点重合，此时找出的所有基点构成的曲线长度就是叶片长度。

1.4.4叶倾角 叶倾角是叶片腹面的法线与天顶轴的夹角，实际上它也是叶面与XOY平面的夹角。在得到叶片有向包围盒后，进一步将其逆变换至原始植株坐标系下，此时包围盒Z轴与天顶轴之间的夹角即为叶倾角（图6C）。

1.4.5茎粗 本研究统一将花盆口上方5 cm 处的平面上茎秆点云欧式距离最远的2 个点定义为最大点与最小点，将2 点水平方向上的差值作为大豆植株的茎粗。利用直通滤波截取目标茎秆点云，茎粗测量如图6D所示。

1.4.6表型参数真实值测量 针对叶片表型参数，随机选取10株大豆植株中的30片叶片用来测试所提出的表型参数测量方法。叶面积、叶宽、叶长和叶倾角均为30个样本，茎粗为10个样本。对涉及到的表型参数真实值作以下约定。

叶面积：将利用标定法[20]对获取的叶片二维图像进行像素转换得到的值作为叶面积真实值。

叶宽：将用卷尺测量平铺叶片的最大宽度距离值作为叶宽真实值。

叶长：将用卷尺测量平铺叶片叶枕到叶尖的距离值作为叶长真实值。

叶倾角：将叶片腹面与水平线间夹角值作为叶倾角真实值。

茎粗：将用游标卡尺测量花盆口平面上方5 cm处宽度距离值作为茎粗真实值。

1.5 评价指标

茎叶分割的主要目的在于分离大豆植株茎秆与冠层叶片，以分割后冠层叶片点云数量评价其分割效果。冠层叶片点云数量用Nc表示，分割前点云总数用N表示，基于以上统计根据式（6）计算分割后冠层叶片点云分割率（Rc）。

单叶分割的主要目的在于分离大豆植株冠层中粘连的叶片，以分割后叶片聚类点云总数评价其分割效果。叶片聚类点云总数用Nl表示，分割前点云总数即冠层叶片点云数Nc，基于以上统计根据式（7）计算分割后单叶点云分割率（Rl）。

采用线性回归分析方法评价表型参数（叶面积、叶宽、叶长、叶倾角和茎粗）测量结果，以本研究算法测量值与人工实测值之间的平均相对误差（mean relative error，MRE）、均方根误差（root mean square error，RMSE）和决定系数R2为指标，量化各参数的精度与误差。

2 结果与分析

2.1 器官分割结果

2.1.1茎叶分割结果 研究表明，当Rs>Rm且Rl≤2Rs的情况下才能体现不同尺度下法向量差异。由图7 可知，当尺度比例较小时（Rl=2Rs），冠层叶片与茎秆差异不明显，无法分离出冠层叶片；尺度比例较大时（Rl=5Rs），大量叶片点云错误归类为茎秆点云；当Rl=3Rs时，叶片部分与茎秆的差异最为明显。在该条件下，当阈值较小时（threshold=0.03），冠层叶片点云过分割；当阈值较大时（threshold=0.10），冠层叶片点云欠分割。为了能够最大程度地去除茎秆并且保留冠层叶片部分，试验均选取阈值为0.08的分割结果进行后续处理。

图7 不同参数下茎叶分割结果Fig. 7 Results of stem and leaf segmentation under different parameters

为了验证茎叶分割方法的准确性和有效性，将本研究提出的方法（DoN）与随机采样一致性法（RANSAC）[15]进行比较。表1 为不同方法下植株茎叶分割结果统计，可知本研究方法冠层叶片点云分割率整体高于随机采样一致性法，冠层叶片点云平均分割率分别为84.24%和78.34%。

表1 植株茎叶分割结果统计Table 1 Statistics of stem and leaf segmentation results of plants

2.1.2单叶分割结果 由图8 可知，当阈值设置过大（number of neighbours = 60，smoothness threshold=15/180π，curvature threshold=1）时，叶片不能完全分割并且有大量点归类于未分类的红色点集中（图8A）；随着阈值减小（number of neighbours=30， smoothness threshold=12/180π，curvature threshold=0.05），大量未分类的点能够重新归类，叶片分割效果最好（图8B），在该条件下基于欧式聚类回填部分点云进而得到完整的叶片点云（图8C）。表2为植株叶片分割结果统计，利用本研究算法实现叶片分割后，叶片单叶点云分割率均达到95%以上，单叶点云平均分割率为96.39%。

表2 植株叶片分割结果统计Table 2 Statistics of leaf segmentation results of plants

图8 不同参数下叶片分割结果Fig. 8 Results of leaf segmentation under different parameters

2.2 表型参数测量结果评价分析

为了验证叶面积计算方法的准确性和有效性，将本研究提出的方法（LDT）与投影法（projection）[21]和贪婪投影三角化法（GPT）[22]进行比较。将大豆叶片叶面积测量值与真实值进行拟合，图9 显示了3 种叶面积测量方法的回归比较，其中本研究方法的回归线最接近标准回归线y=x，回归方程为y=0.971 1x+0.271 7，决定系数R2=0.987 9, 算法测量值与真实值接近，表明本研究方法可以应用于大豆叶片叶面积的测量。通过表3 可知，在测量叶片面积时，投影法平均相对误差最大，决定系数R2和均方根误差分别为0.978 7和1.251 1 cm2；贪婪投影三角化法平均相对误差低于投影法但高于本研究方法，决定系数R2和均方根误差分别为0.979 1 和0.763 3 cm2；本研究提出的方法平均相对误差最小，决定系数R2和均方根误差分别为0.987 9 和0.541 7 cm2，表明本研究方法整体优于其他2种方法。

表3 不同方法测量叶面积结果评价Table 3 Result evaluation of leaf area by different methods

图9 大豆叶面积不同方法测量值与真实值比较Fig. 9 Comparison of soybean leaf area by different method

通过表4 可知，叶宽、叶长及茎粗算法测量值与真实值的平均相对误差分别为2.68%、2.87%和3.99%，决定系数R2分别为0.961 3、0.962 6 和0.963 4，均方根误差分别为0.142 2、0.175 5 和0.047 5 cm。将叶宽、叶长及茎粗算法测量值与真实值进行拟合，拟合结果如图10A、B和C所示，回归方程分别为y=0.984 8x+0.086 3、y=0.897 0x+0.424 1 和y=0.830 2x+0.215 5。算法测量值与真实值接近，表明本研究方法可以应用于大豆叶片叶宽、叶长及茎粗的测量。

表4 叶宽、叶长及茎粗测量结果评价Table 4 Result evaluation of leaf width， leaf length and stem diameter

图10 大豆表型参数测量值与真实值比较Fig. 10 Comparison of soybean phenotypic parameters extracted with actual values

叶倾角算法测量值与真实值的平均相对误差为7.22%，决定系数R2和均方根误差分别为0.931 1 和3.279 6°。将叶倾角算法测量值与真实值进行拟合，拟合结果如图10D 所示，回归方程为y=0.986 0x+0.442 2。叶倾角算法测量值与真实值接近，但其平均相对误差相较于其余表型参数而言偏高，表明叶倾角测量算法在应用于大豆植株叶倾角测量方面仍有较大改进空间。

2.3 时效统计分析

本研究采用Intel Core i7-5930K 处理器以及NVIDIA Quadro P4000 GPU 进行试验，不同处理阶段时效统计结果见表5。每株大豆植株点云重建平均用时为24.17 min，滤波处理、茎叶分割、单叶分割以及表型参数提取平均用时分别为7.35、14.63、8.42 和46.72 s。点云重建受场景以及特征点数量的影响耗时较长，在匹配特征点时通过缩小目标区域能够进一步提高重建效率。

表5 时效分析结果Table 5 Results of processing time analysis

3 讨论

基于三维点云的植株表型参数测量为表型数据获取提供了有效的解决方案。相对于传统的常规测量方法，基于三维点云的测量技术具有非接触、速度快等诸多优点。相较于二维图像，三维点云数据包含的深度信息能够使植株与无关场景天然解耦，很大程度上避免了植株本身形变造成的误差。已有研究探讨了基于植株三维点云测量表型参数的实践[8-9,13-15]，并取得了良好的效果，但仍需要大量人为干预。本研究提出的器官分割及表型分析方法与已有研究[14-15]相比，在一定程度上解决了叶片相互粘连造成的茎叶分割及高通量表型测量困难等问题，提高了表型分析的效率和准确度。通过法向量差异实现茎叶分割，并结合点云曲率特征对冠层叶片实施了基于欧式聚类的区域生长分割，以上器官尺度的分割使植株表型参数（叶面积、叶宽、叶长、叶倾角及茎粗）测量更加自动、快速、可靠。结果显示，器官分割后冠层叶片点云平均分割率为84.24%，单叶点云分割率均高于95%，叶面积、叶宽、叶长、叶倾角和茎粗参数测量值与人工实测值的决定系数分别为0.987 9、0.961 3、0.962 6、0.931 1 和0.963 4，均方根误差分别为0.541 7 cm2、0.141 2 cm、 0.175 5 cm、3.279 6°和0.047 5 cm。每株大豆植株点云重建和滤波处理、茎叶分割、单叶分割以及表型参数提取平均用时分别为24.17 min 和7.35、14.63、8.42 和46.72 s。该结果表明，本研究方法能够快速地实现大豆植株器官分割，为多分枝作物高通量自动化表型参数测量提供高效的技术支持。

本研究方法可以类推应用于其他冠层密度大，叶片较为平坦，叶片呈近圆形、椭圆形的多分枝植株。在后续的研究中，一方面应考虑冠层密度小但叶片卷曲，叶片呈披针形的植株，提高器官分割及表型测量算法应用于不同形态植株上的精度。另一方面，针对大豆开花结荚期分枝和叶片较多、叶片遮挡严重的情况，未来应根据植株自相似性考虑从部分到整体的三维重建与表型分析。