玉米果穗的表型性状是反映玉米果穗形态与结构所有特征的总和, 是玉米考种的常规测量项目。玉米果穗的表型性状可分为几何性状、数量形状与颜色性状, 其中数量性状主要包括穗行数、行粒数和总籽粒数等。长期以来, 果穗的数量性状主要通过人工统计来完成, 存在效率低、工作量大、主观误差大、精确性难以保证等问题。如何快速、准确测量果穗数量性状成为玉米自动化考种和产量预测的热门研究课题。

近年, 国内涉及玉米果穗数量性状的自动化考种技术主要有3个研究方向。第一, 采集截断果穗的截面图像, 计算各数量性状, 但对果穗产生损伤:韩仲志等^[1]通过采集果穗横断面图像, 依据边缘标记和离散曲率, 得到穗行数。第二, 采集旋转果穗的图像序列, 通过图像拼接技术^[2]获得整穗图像, 并直接测量整穗数量性状:王传宇等^[3]利用全景技术采集旋转果穗图像序列, 无缝拼接后计算行数与行粒数; 柳冠伊等^[4]基于线阵扫描技术获取旋转果穗表面图像, 采用OTSU算法分割籽粒, 通过像素累计识别穗行; 刘长青等^[5]分析旋转果穗局部区域的像素累计, 提取出穗行并得到行数与行粒数; 杜建军等^[6]利用步进电机旋转果穗, 通过畸变校正、图像拼接后得到完整穗粒分布图, 并提取果穗表型形状参数。第三, 采集单幅果穗图像, 通过建立数学模型估算整穗数量性状:梁秀英^[7]通过建立理想几何模型计算行数, 利用改进的阈值分割算法与分水岭算法获得行粒数; 周金辉等^[8]通过玉米颜色特征与生物特性, 建立估算模型得到行数与行粒数。在国外, 机器视觉在玉米果穗或种子上的应用主要集中在种子损伤与活性检测、黄曲霉污染检测以及种子品种分类与品质分级等方面^[9-11], 与果穗数量性状相关的文献较少。

本文应用机器视觉技术识别果穗, 提出了基于单幅图像的自动化玉米果穗数量性状测量方法, 该方法能较快速、精确地测量果穗行数、行粒数和总粒数, 为玉米自动化考种提供了基础。

1 设备与算法流程

图像采集设备如下:CMOS工业相机(维视数字图像技术有限公司提供, 型号MV-EM1400C, 1 400万像素, 图像分辨率3 288×4 608, 红、绿、蓝的比例为1.75：1：1.4), LED条形光源(型号AFT-WL21244-22W, 白光, 最大功率22 W), 光源控制器, PC机(x64双核处理器2.4 GHz, 4 GB内存, Windows 10操作系统), 图像采集设备如图 1-a所示。图像处理平台为MATLAB 2014a, 在处理开始前用邻近插值法缩小图像尺寸为1 096×1 536, 算法流程如图 1-b所示。

2 步骤与方法 2.1 果穗轮廓提取

YCrCb(亦称YUV)色彩空间将色彩同亮度分离, 适合处理光照条件复杂的图像, 且YCrCb与RGB之间的转换为线性, 算法简单^[12], 转换公式^[13]如下:

(1)

式中:R, G, B∈[0, 255]。本文将Cr改进后用于果穗轮廓提取, 公式如下:

(2)

对改进后的Cr′通道采用双峰法进行二值化^[14-15], 提取二值图像上最大连通域(代表果穗)的轮廓, 将原图像RGB各通道与轮廓作Hadamard乘积, 获得仅含果穗的彩色图像。

该轮廓提取算法对背景颜色、光照条件、背景阴影与反光等具有较强的适应性, 效果如图 2所示。

2.2 像素累计

现阶段识别穗粒的常规方法是通过籽粒的颜色特征进行图像分割, 但图像分割通常需要结合先验知识, 并没有一种普适最优的分割方法^[16], 而籽粒颜色受果穗品种、籽粒含水量、图像采集环境等影响, 因此难以保证一种方法准确分割大多数果穗。不同条件下的分割效果如图 3所示。

根据果穗籽粒纵横排列的特点, 可通过累计像素的方法定位穗行与行间隙^[4-5], 该方法避免了对穗粒进行图像分割, 所以具有较高的精确性与可靠性。

2.2.1 校正果穗倾斜角

为保证像素累计的准确性, 需要校正果穗倾斜角并使其竖直于图像中。本文采用Radon变换法求得果穗倾斜角^[17]。因为轮廓边缘有不规则突起, 在Radon变换前将轮廓进行过度的腐蚀处理, 使轮廓近似呈椭圆形。得到倾斜角后, 以质心为基准将图像旋转以完成校正。

2.2.2 分区域累计像素

建立如图 4-a所示坐标系, 记P(x, y), x∈[x_ss, x_ee], y∈[y_ss, y_ee]表示果穗最小外接矩形图像上坐标为(x, y)的像素值, 将图像沿y轴等分为20个区域S₁, S₂, …, S₂₀, 对每个区域沿y方向累计像素, 计算公式^[5]如下:

(3)

累计像素得到的像素值分布函数记为f_y1(x), f_y2(x), …, f_y20(x)。以S₁₀为例, 其累计效果如图 4-b所示。

2.3 计数行粒数 2.3.1 定位穗行线

由图 4-b可知, 像素值分布函数的极大值点代表该区域的穗行, 将20个区域的极大值点(即穗行定位点)按一定顺序连接即可获得穗行线。本文仅定位果穗中间较清晰的若干穗行。

对f_y1(x), f_y2(x), …, f_y20(x)进行一维高斯滤波处理, 并提取极大值点, 得到集合D₁, D₂, …, D₂₀, 每个集合的元素代表该区域穗行的定位点。在相邻区域, 同一穗行的定位点在x方向上位置变化不大, 所以按距离最近原则对D₁, D₂, …, D₂₀内的元素分组, 使同一穗行的定位点在同一组内。对各组定位点进行线性插值, 即可获得所有的穗行线。定位效果如图 5所示。

2.3.2 计数方法

沿着穗行线在灰度图像上抽取像素值, 得到像素值分布函数m(x)。对m(x)进行一维高斯滤波处理, 并提取极大值点, 统计极大值点的数量即得行粒数。将中间4行的行粒数取平均值。从左往右, 图 5-a中第一条穗行的像素值分布函数如图 6所示。不考虑曲线上的毛刺, 图 6中有41个主要的波峰(极大值), 代表该穗行上有41颗籽粒。

2.4 估算穗行数 2.4.1 定位行间隙

由图 4-b可知, 像素值分布函数的极小值点代表该区域穗行间隙。本文对f_y8(x)、f_y10(x)、f_y12(x)进行分析。以f_y8(x)为例, 对其进行一维高斯滤波处理, 并提取极小值点, 得到集合E₈, 即区域S₈行间隙的位置。通过扫描轮廓左右边缘, 得到S₈横截面的半径R₈和圆心O₈, 取距离O₈不超过0.75R₈的穗行间隙。同理, 定位区域S₁₀和S₁₂的行间隙。定位效果如图 7所示。

2.4.2 估算模型

本文通过横截面几何模型估算整穗行数。如图 8所示, 假设果穗横截面为圆, 忽略投影线倾斜造成的半径误差(即拍摄高度远大于穗粗), O为截面圆心, R为截面半径, A和C是中间若干行的左右边缘(穗行间隙), a和c分别是它们的投影, P和Q分别是果穗左右边缘, OB⊥PQ, ∠AOB=α, ∠BOC=β。

计算中间若干行的圆心角∠AOC为:

(4)

若∠AOC包含n条穗行, 由比例式

(5)

可得整穗行数

(6)

2.4.3 误差分析与校正

进一步研究发现行间隙的凹陷使横截面并非圆形, 并产生误差, 如图 9-a所示。图 9-b中, 粗实线表示果穗真实轮廓, 细虚线表示轮廓外接圆(代表纸带), 带箭头实线表示射入镜头的光线, C表示间隙的实际位置, C′表示被定位的位置(光线与轮廓的切点), L表示直线OC′在水平方向上的投影长度, Δθ表示实际位置与定位位置之间的偏差角, W表示间隙宽度。可知该误差会使∠AOC偏小, 导致行数偏大。

分析后发现, 在间隙宽度W和截面半径R一定的前提下, 偏差角Δθ是水平距离L的函数, 记为Δθ=f(L)。可知Δθ=f(L)是单调递增函数, 且在果穗中心B点(L=0)取得最小值0, 在果穗的图上边缘Q′点(L≈R)取得最大值W/2R(此时C与C′之间夹1/2个间隙宽度)。于是函数可化为:

(7)

通过Auto CAD计算得到:

(8)

由图 4-b可知, 波形图中波谷的宽度反映了间隙宽度W, 但波谷上代表W的实际宽度MN较难确定(如图 9-c所示), 因此本文以波谷的中径为衡量标准间接度量间隙宽度W。波谷中径d_m的概念如图 9-c所示。对多根果穗的图 4-b曲线进行测量后得到W≈1.2d_m。再结合式(7)、(8)得到偏差角的计算公式:

(9)

根据式(9)在式(6)的分母中补偿误差, 得到整穗行数的最终计算式为:

(10)

式中:n为∠AOC包含的穗行数; d_m为距果穗中心最近波谷的中径。

3 试验验证与分析

本文选用‘登海605’和‘先玉335’2个品种共计165个果穗进行试验验证, 以证明本文所述方法的可行性和精确性。利用图 1-a所示设备采集果穗图像, 然后利用本文方法测量行数与平均行粒数, 行数的小数结果就近取偶(果穗实际行数必定为偶数), 行数与平均行粒数相乘得到总粒数。人工统计每个样本的行数与每行籽粒数, 将各行粒数取平均得到平均行粒数, 各行粒数求和得到总籽粒数, 以人工统计值作为真实值。使用式(11)计算各性状平均测量精度(P)^[6], 使用式(12)计算行粒数的平均绝对偏差(D)。

(11)

(12)

式中:m为果穗样本个数; x_i和x_i分别为果穗各性状的测量值与真实值。对于穗行数统计零误差率。

使用Matlab探查器统计测量速度, 得到平均测量速度为每穗1.936 s。在图像上对定位点进行标记, 人工检查后发现有9根果穗定位错误, 得到穗行定位准确率为94.55%。统计试验结果, 得到总体的穗行数、行粒数、总粒数的平均测量精度分别为98.84%、96.33%、95.67%, 行粒数的平均绝对偏差为每行1.24粒, 穗行数零误差率为90.91%。绘制统计直方图, 如图 10所示。

从直方图形状上看, 3种测量指标的分布大致呈偏态分布, 且分布趋势较偏向100%精度, 说明该测量方法具有得到正确结果的趋势, 方法具有可行性; 从样本分布区间上来看, 除少数结果精度为80%~90%, 大部分结果精度在90%以上, 说明该测量方法能较精确地得到测量结果, 方法具有精确性; 从整体统计数据上来看, 3种测量指标的整体统计结果均在95%以上, 而样本数量足够多(参考文献[1, 3-8]的试验样本均未超过100), 整体统计结果较为可信, 亦能说明本测量方法的精确性。

为了验证本文对行数计算误差的校正效果, 进行对照试验。操作方法为:试验组和对照组均按本文方法测量果穗行数, 除最终行数计算公式不同以外其余均相同。试验组采用式(10)计算穗行数(校正误差), 对照组采用式(6)计算穗行数(未校正误差), 测量的小数结果均不取偶, 计算测量结果与真实值的平均绝对偏差, 统计测量结果(表 1)并绘制偏差的分布直方图(图 11)。

从图 11可以看出, 统计结果的分布均大致呈正态分布, 但对照组的均值明显大于0, 而试验组的均值较为接近0, 说明该误差对行数测量结果有较明显的增大效果, 而本文的校正方法能有效抑制误差对结果的增大效果。从表 1可以看出, 对照组的大多数测量结果比真实值偏大, 而试验组偏大的样本数量约占总体的50%, 平均绝对偏差也比对照组显著减小, 亦说明本文校正方法能有效校正该误差。

在对测量误差较大的样本进行分析, 发现本测量方法的误差主要来自行粒数的测量。对于穗行严重倾斜或人眼无法判断穗行的果穗, 本文的方法识别效果较差, 当识别错误时行粒数的测量精度将无法保证。另外, 对于穗粒中部发生严重下凹的情形, 由于凹陷处形成阴影, 行粒数测量结果可能偏大。

4 讨论

综合以上试验结果, 与玉米自动化考种的现阶段研究相比较, 本方法在多个测量参数上具有优势。当前玉米自动化考种的两个主流研究方向分别为旋转果穗测量法和单幅图像测量法。目前旋转果穗测量法的主要测量参数可达到:行数精确度98.231%, 行粒数精确度94.351%, 总粒数精确度96.921%, 测量速度每穗12 s^[6]。与旋转果穗测量法的现阶段研究水平相比, 本方法的行数测量精度(98.84%)、行粒数测量精度(96.33%)较高, 测量速度(每穗1.936 s)较快, 且无需步进电机等设备。目前单幅图像测量法的主要测量参数可达到:行数零误差率在93%以上, 行粒数绝对误差在2粒左右, 测量速度达到每分钟30穗以上^[8]。与单幅图像测量法的现阶段研究水平相比, 本方法的行粒数平均绝对偏差(每行1.24粒)较小, 并且本文给出了总粒数的测量结果(平均精度95.67%)。另外, 本文提出一种出现在单幅图像中的使穗行目测位置偏向中部的成像误差, 并给出解决方案。