基于图像处理和聚类算法的待考种大豆主茎节数统计

王跃亭 王敏娟 孙 石 杨 斯 郑立华

(1.中国农业大学现代精细农业系统集成研究教育部重点实验室， 北京 100083；2.中国农业科学院作物科学研究所, 北京 100081)

0 引言

为了熟悉新培育出来的大豆品种，研究人员需要对大豆植株性状进行详细的观察和记录(即考种)[1-2]。传统的人工考种存在操作误差大、效率低下、成本过高等问题。为了优化作物考种流程，研究者进行了大量研究：机器视觉、深度学习等技术应用于作物籽粒性状检测，提升了检测速率与准确性[3-5]；三维视觉技术的引入为获取作物茎秆信息，以及作物产量、抗倒伏等性状分析提供了支持[6-8]；高光谱技术与图像分析技术的结合为作物叶片性状的采集与监测提供了新的方式[9-11]。

豆荚多生长于主茎上的茎节点处，主茎节数关系到大豆的产量，故主茎节数是大豆考种过程的重要指标之一。文献检索发现，目前利用机器学习和图像处理方法进行大豆主茎节点和主茎节数统计的研究鲜有报道，但其他作物有关茎节点和茎节的研究均基于图像处理展开。文献[12-13]通过对甘蔗图像进行图像分割、局部像素点统计等对甘蔗茎节点进行识别；MOSHASHAI等[14]根据茎节处粗大的形态特点，通过检测甘蔗边缘信息，将局部最大值确定为茎节。上述方法虽然在甘蔗茎节点识别中获得不错的效果，但因大豆植株存在豆荚、分枝等因素的干扰，局部像素点较多处和边缘粗大处等不一定为大豆茎节点。葡萄、荔枝等果实采摘点与大豆豆荚生长习性类似，果实串采摘点的确定对大豆主茎节点的识别有一定借鉴意义。熊俊涛等[15-16]在识别、分割葡萄、荔枝等作物果实串的基础上，确定果梗感兴趣区域，并辅以角点检测等算法确定采摘点，为农业采摘机器人的实际应用提供了技术支持。上述方法应用于大豆茎节点识别时，需要解决大豆因色彩相近、遮挡严重等原因而造成的豆荚识别困难，以及豆荚与主茎连接处较短、不易识别等问题。

综合分析相关研究可知，对于大豆植株存在的豆荚、分枝干扰问题，可利用骨架提取等技术进行图像细化，配合角点检测可进一步实现对主茎节点的确定。针对因拍摄角度造成大豆茎节点无法有效识别的问题，可利用聚类算法[17]对多角度下的图像信息进行统计和分析，以实现大豆植株茎节点的识别、统计和分析。其中，HDBSCAN (Hierarchical density-based spatial clustering of applications with noise)[18-21]结合了密度聚类、层次聚类算法的优点，无需设定初始聚类簇大小、领域阈值等参数，能对噪声点进行有效识别和剔除，适用于未知节数情形下的大豆主茎节数统计研究。

本文在上述研究基础上，提出一种基于图像处理和聚类算法的待考种大豆主茎节数统计方法。以不同拍摄角度下的大豆图像作为输入，引入图像处理与聚类算法对大豆主茎节数进行逐级筛选、提纯，最终实现大豆主茎节数的统计。

1 材料与方法

1.1 图像采集和处理

实验对象为已脱叶、处于考种阶段的高产大豆品种“中黄30”，由中国农业科学院作物科学研究所提供。采集图像时，将种植于花盆内的待测植株放置于转盘上，背景采用黑色幕布，使用三脚架将相机架设到适宜高度，数据采集场景设计如图1a所示。拍摄时，转动转盘实现不同角度的大豆植株图像获取，图像采集实验场景如图1b所示。

图1 数据采集实验场景Fig.1 Scene of data acquisition

相机分辨率为1 920像素×1 080像素，可获得彩色图像。为了避免单一拍摄角度下处于视觉盲区的茎节点丢失，实验时设定旋转基础步长为15°，转盘旋转一周可获得24幅不同角度的单株大豆植株图像。本文在Windows 10操作系统下，利用PyCharm集成开发平台，采用Python语言，实现大豆植株图像的植株分割、骨架提取、角点识别、主茎节点筛选、数据转换以及密度聚类统计等操作，进而完成大豆植株主茎的茎节统计。

1.2 大豆植株主茎节数统计算法

图2 大豆植株主茎节数统计算法流程图Fig.2 Flow chart of statistical algorithm for number of main stem nodes of soybean plants

为了有效统计大豆植株主茎节数，本文设计了如图2所示的基于图像处理和聚类算法的待考种大豆植株主茎节数统计算法。算法将采集到的大豆植株图像随机分为训练集和验证集，并进行对应的训练、验证操作。其中，训练的目的是获得主茎节数识别准确率最大值所对应的图像采集间隔以及聚类参数，据此表明对同一品种大豆，只需训练少量样本即可获得稳定、收敛的结果。验证集则可进一步验证上述算法的效果。

在训练阶段(图2a)，首先通过人工输入初始图像采集间隔和抽样步长，初始化算法参数。其次，根据上述参数抽取训练集样本，并经过图像分割、骨架提取、角点识别、主茎节点筛选等操作，获取待检测茎节点数据集。再次，对分布离散、距离较难计算的二维数据点进行数据降维，转化为更易聚类的一维空间数据。然后，利用HDBSCAN算法对不同采集间隔下的待检测大豆主茎节点进行聚类分析，完成当前大豆植株的主茎节数的识别、统计，并计算对应的大豆植株识别准确率。最后，依据抽样步长减小图像采集间隔，重复上述步骤直至训练集数据运算完毕，根据主茎节数识别准确率筛选出最优采集间隔与聚类参数。

在验证阶段(图2b)，根据训练阶段所获得的最优采集间隔与聚类参数，对验证集图像进行对应图像的处理、聚类操作，并对获得的大豆识别准确率进行统计，进一步验证算法的准确性、有效性。

1.2.1植株分割

在本文设计实验场景下采集的大豆植株图像包含花盆、转盘、幕布等其他背景信息以及干扰噪声。为了获得准确、有效的大豆植株图像，需要对图像进行分割以获得仅含有目标的图像。对此，本文利用HSV色彩模型制作掩模，并以此实现对大豆植株的图像分割，获取只含大豆植株的图像。其中，HSV色彩模型是一种面向视觉感知的颜色模型，相较于RGB色彩空间更加符合人的视觉特征。HSV色彩模型由3个要素组成，即色调(H)、饱和度(S)以及亮度(V)。

1.2.2骨架提取

植株茎节点是大豆主茎节数计数的标志点，通过识别、标识、统计植株茎节点即可确定主茎节数，进而实现大豆产量的评估。若直接进行茎节点识别，因存在豆荚遮挡、主茎分枝等干扰，将引起茎节点遗漏、重复识别、错误识别等问题(图3)。故本文利用骨架提取算法进行大豆植株图像细化，以保证主茎节点的准确、有效提取。

图3 大豆植株图像噪声分布示意图Fig.3 Noise distribution schematic of soybean plant

在图3中，编号1处的茎节点存在豆荚遮挡问题；编号2、3处因豆荚与主茎的交叉而引入2个噪声点；编号4、5处则因分枝与主茎相交而引入2个噪声点；编号6处的噪声点由豆荚末端引入；编号7处则为大豆主茎与成像边界相交而引入的噪声点。

本文首先利用大津法(Otsu)[22]对图像分割后的大豆植株图像进行二值化，然后选择由LEE[23]提出的骨架提取算法进行大豆植株的骨架提取。

1.2.3大豆主茎节点检测

大豆主茎节点是植株主茎上两节茎秆的连接部分，而相邻两节茎秆的直径、方向不同等使得其在该点存在属性差异(即角点)。为了获取大豆植株主茎上待检测的茎节点，本文选择Shi-Tomasi特征点检测算法[24]对大豆植株骨架进行角点检测。

经过角点检测后，骨架端点处仍存在因豆荚、分枝以及成像限制等产生的噪声点。为了剔除噪声点，本文首先设定端点识别阈值，选定卷积模板；然后，根据设定的卷积模板尺寸构建卷积模板，对骨架提取后的大豆植株图像各像素点进行卷积操作；最后，将各卷积结果与端点识别阈值进行比较，小于阈值的点即为所寻找的端点。通过对20组不同植株在任意角度下拍摄的图像进行主茎节点提取试验，观察发现当卷积模板尺寸设置为3×3、端点阈值设置为4像素时，端点识别效果较好，能够识别出所有非遮挡的骨架节点，平均准确率达90.39%。

1.2.4基于空间距离的数据转换

大豆植株生长过程中的弯曲、地面不平整等，使得不同拍摄角度的大豆主茎节点分布无序，从而使得后续的聚类算法缺乏合理的聚类依据。故需将空间分布无序的大豆茎节点进行数据转换，得出便于数据空间聚类的数据集。

本文选择MDS降维算法[25]将不同分布的待检测二维空间点映射至一维空间内，简化后续的聚类操作。MDS算法是一种基于空间距离的数据降维算法，其要求在低维空间内样本距离能够保持原始空间中的各样本之间的距离[26]。

算法假定m个样本在原始空间(其维度为d)的距离矩阵为D∈Rm×m，其第i行、第j列的元素dij为样本xi到xj的距离。样本在d′维空间的标识形式为Z∈Rd′×m,d′≪d，且任意两个样本zi、zj在d′维空间中的欧氏距离等于原始空间中的距离，即‖zi-zj‖=dij。

令B=ZTZ∈Rm×m，其中B为降维后样本的内积矩阵。对矩阵B做特征值分解，B=VΛVT，其中V为特征向量矩阵；Λ=diag(λ1,λ2,…,λd)为特征值构成的对角矩阵，λ1≥λ2≥…≥λd。假定其中有d*个非零特征值，它们构成对角矩阵Λ*=diag(λ1,λ2,…,λd*)，令V*表示相应的特征向量矩阵，则Z可表达为

(1)

1.2.5HDBSCAN聚类算法

HDBSCAN算法是一种密度聚类与层次聚类相结合的聚类算法。其只需确定每个簇中样本数的最小值(即最小簇大小)即可实现样本集的可视化聚类。其在未知大豆主茎节点的情况下，只需考虑待检测茎节点在总的抽样样本集中出现的频次即可，更加适合大豆主茎茎节点的聚类统计。算法具体步骤如下：

(1)根据数据点分布密度进行空间变换。为确定不同簇的聚类依据，算法引入了核心距离与互达距离的概念。核心距离ck(x)=d(x,Nk(x))指样本与第k个最邻近样本点的距离。a、b两点的互达距离dmreach-k(a,b)计算公式为

dmreach-k(a,b)=max{ck(a),ck(b),d(a,b)}

(2)

式中d(a,b)——a、b两点的直线距离

对于空间中存在的点A、B、C，当k=5时，A、B两点和A、C两点的互达距离分别为c5(x)、d(A,C)。互达距离示意图如图4所示。

图4 互达距离示意图Fig.4 Diagram of reaching distance

(2)依据互达距离，建立最小生成树。为了能够获取数据集内所有数据点的信息，HDBSCAN算法利用Prim算法，以数据点为顶点、任意两点间互达距离为权重，构建最小生成树。

(3)最小生成树转换为簇层次结构。簇状结构能够对不同层次的数据集进行划分，故需将最小生成树转化为簇层次结构。算法以互达距离为权重，对树的边按升序进行排序、遍历，为每条边创建一个新的合并簇。

(4)压缩簇层次结构。为了剔除数据集中的噪声点、压缩簇层次结构，需要人工设定最小簇的值。在每次簇层次结构分割时，比较新分割出的簇点数与最小簇大小：当新簇不少于最小簇时，新分割的簇保留，否则该簇被剔除且保留较大的父簇；最终获得一个拥有较少分枝的层级树结构。

(3)

式中cluster——聚类簇集合

将所有子簇均声明为选定簇，通过反向拓扑排序遍历比较子簇稳定度与父簇稳定度：当子簇稳定度之和大于父簇稳定度时，则将簇稳定度设置为子簇稳定度之和；否则，将父簇声明为选定簇，并取消选定其子簇；当到达根节点时，当前选定簇集合即为最终提取的稳定簇集合。

1.2.6识别准确率计算及寻优

为了衡量训练集大豆主茎节数的识别效果，需要设计、制定合理的评价标准。因单个植株主茎节数较小(“中黄30”大豆多为12～18节)，即使较小的误差(1～3节)也将影响每一样本的精度统计结果。故针对不同采集间隔下的识别准确率P，其计算公式为

(4)

式中NI——全部训练集识别的主茎总节数

NT——实际主茎总节数

在不同采集条件下对不同植株的识别准确率进行记录，将最高识别准确率的采集条件作为最优采集间隔，以供剩余植株的主茎节数识别、统计。

2 结果与分析

2.1 大豆植株图像分割

为了方便大豆植株的分割，在植株后方布置灰色幕布，因此HSV色彩模型选定灰色色彩空间(Hmin=0，Hmax=180；Smin=0，Smax=30；Vmin=46，Vmax=220)进行背景色彩描述；然后，通过取反操作制作掩模与拍摄图像进行与操作完成大豆植株的图像分割。图5为4株已脱叶的待考种“中黄30”大豆在转角0°时的原始图像及目标分割后的图像。

图5 大豆植株原始图像及分割后图像Fig.5 Images of original and image segmentation

2.2 大豆植株骨架提取与角点检测

2.2.1大豆植株骨架提取

在获取到只含有大豆植株的图像后，首先利用大津法对图像进行二值化处理，将大豆植株所在区域灰度设置为255(即白色)，背景区域灰度设置为0(即黑色)；然后，利用骨架提取算法进行二值图像的骨架提取，对图5进行处理的结果如图6所示。

图6 大豆植株二值化图像与骨架提取效果Fig.6 Images of binarization and skeleton extraction

通过骨架提取将大豆植株细化为宽度为1像素的二值化图像，减少同一豆荚、分枝因其空间不同而引入的噪声点数量。此外，在宽度为1像素的植株图像中更容易确定茎节点的位置。

2.2.2大豆植株角点检测

为保证获得所有可能茎节，针对考种阶段的大豆植株，通过前期试验设定最大识别角点数为100、识别等级为0.01、最小间隔为30像素进行大豆植株骨架的角点检测，图7为图6的处理结果，并将对应的图像坐标信息存入Excel表格内备用。

图7 不同视角图像的待检测茎节点识别结果Fig.7 Recognition results of nodes from different perspectives

图8 大豆植株骨架端点识别及筛选效果Fig.8 Recognition and filtering of skeleton endpoint and screening

由图7可知，通过角点检测获得的待检测茎节点中存在由豆荚骨架末端以及茎秆与图像边界相交所引入的噪声点，需要进行大豆植株主茎节点的筛选和提纯。

2.3 大豆植株主茎节点筛选

为了将骨架端点处的噪声点剔除，本文对骨架图像端点进行了识别、记录，并将其从待检测茎节点数据集中剔除，图7中各目标的处理结果如图8所示。

观察发现，通过图像的卷积操作能够保留植株主茎上的茎节点，并滤除植株、分枝、豆荚等部位的末端噪声点干扰。因已经滤除主茎末端噪声干扰，现存的茎节点个数即可应用于后续聚类分析。

2.4 基于空间距离的数据转换

完成单幅图像的主茎节点筛选后，需要将不同视角下的主茎节点进行汇总，并将数据转换为更适宜聚类的一维空间数据。因不同角度下采集的图像尺寸(1 920像素×1 080像素)相同，且每幅图像具有相同的像素坐标系，故不同视角下的主茎节点可放于同一文件内直接进行比较。

图9为从同一株大豆植株的3组图像中提取的主茎节点组合以及经过MDS数据转换的茎节点空间分布效果图。这3组图像分别为拍摄角度相隔180°的2幅图像(图9a)、相隔120°的3幅图像(图9b)、相隔90°的4幅图像(图9c)。其中，左图为茎节点图像组合，右图为降维转换后的数据(为了便于数据点描述，X坐标统一置零，Y坐标为转换后的数值)。

通过对数据转换后的图像进行分析可知，当采集间隔较大、输入图像较少时，转换后的各待检测茎节点分布较为疏散，但数据量较少；随着采集间隔减小、输入图像增多，数据点分布逐渐稠密，在多幅图像中出现的茎节点重叠使得图像中该点颜色加深。

图9 大豆植株主茎节点数据降维效果Fig.9 Effect diagrams of dimensionality reduction of node data

2.5 茎节点聚类

在HDBSCAN算法聚类过程中将最小簇设置为2(在2幅图像中均稳定出现的茎节点即认定为主茎节点)，部分聚类效果如图10所示，其中被圈出的是稳定的聚类簇，即被认定为大豆植株主茎上的茎节点。

由图10可知，当图像采集间隔较大时，输入的待检测茎节点较少，样本量不足以及部分样本点相距较远，导致识别的茎节点数少于真实值。随着图像采集间隔的不断减小，输入的待检测茎节点总数不断增加，更多的待检测茎节点被判定为茎节点，但也更容易受到噪声点的影响。

2.6 最优采集间隔选取与识别效果分析

为了确定最优的图像采集间隔，本文选择21株待考种的“中黄30”大豆植株作为训练集，以180°采集间隔为初始采集间隔，设定抽样步长为30°，分别

对采集间隔为180°、150°、120°、90°、60°、30°时获得的主茎节数识别准确率进行记录、统计和对比分析，结果如表1所示。

由表1可知，当图像采集间隔为90°时，其识别准确率为91.11%，识别效果最优，因此将其作为剩余样本的最优图像采集间隔。本文对7种不同主茎节数的42株验证集大豆植株(即每种主茎节数包括6株样本)进行主茎节数识别、统计，并计算其识别准确率，结果如表2所示。

由表2可知，在图像采集间隔为90°的条件下，不同主茎节数的大豆植株的平均识别准确率可达98.25%，单株考种大豆主茎节数识别准确率大于88%，能够满足大豆实际考种需求。

3 结论

(1)根据大豆植株的生长期及其骨架角点分布情况，通过实验设定卷积模板尺寸，并进行卷积计算，实现大豆植株主茎骨架上的待检测数据点提取。

图10 大豆植株主茎节点聚类效果Fig.10 Clustering of main stem nodes of soybean plants

表1 训练集不同采集间隔下的主茎节数识别结果Tab.1 Results of soybean main stem nodes from different collection intervals

表2 采集间隔为90°下的主茎节数识别结果Tab.2 Results of soybean main stem nodes at 90°

(2)基于空间距离的数据转换算法解决了因大豆植株自然生长、地势起伏等原因而引起数据点分布离散的问题，为聚类算法提供较好的样本数据集。

(3)在未知茎节点个数的情况下，本文方法能够统计出大豆植株主茎节数，并可剔除分枝中茎节点的干扰。实验发现，对于42株待考种“中黄30”大豆植株,利用90°采集间隔获得的图像的主茎节数平均识别准确率达98.25%。