基于标尺图像识别的作物株高测量

孙 翔，裴晓芳，周 望，朱 平

(1.南京信息工程大学 电子与信息工程学院，江苏 南京 210044；2.南京信息工程大学滨江学院 电子与信息工程学院，江苏 无锡 214105；3.南京信息工程大学 江苏省大气环境与装备技术协同创新中心，江苏 南京 210044；4.航天新气象科技有限公司，江苏 无锡 214127)

作物株高是指作物冠层的高度，是衡量作物生长情况的重要指标，能直接反映出作物的生长速度。目前测量作物株高的方式主要有3类：(1)人工测量。该方法耗时费力且受测量人员的主观影响较大；(2)图像处理。通过对拍摄的作物的图像进行处理计算出株高；(3)通过激光雷达[1-3]、RGB-D相机[4-6]等获取作物和地面的点云数据，统计分析得出株高。

农情观测系统中需要采集大量图片，而图片可以反映大量的作物信息，图像处理的方法在该系统中有广泛的应用。图像处理在株高测量上分为双目视觉和单目视觉两种方法。双目立体视觉法[7-10]不需要参照物，但对环境光照非常敏感，匹配算法复杂，计算量大，商用推广难度较大，因此本文采用单目视觉的方法。单目视觉在测量长度时需要一个参照物，本文引入了一个红蓝相间的定制标尺，标尺的每个色块长度均为0.1 m。将标尺竖直放置于作物主茎秆后方，通过颜色识别出标尺未被遮挡的部分，进行滤波后通过连通域分析法计算出各个连通域的像素高度，最后按比例计算出作物的株高。该方案可应用于农作物自动观测或景观植物监测中，通过测量的株高信息反馈植物生产态势，指导人工干预实现增产增收的目标。本方案设备简单、成本低，只需要一个定制标尺和一架相机配合PC端进行处理即可，适合大规模农情观测。由于作物在生长过程中株高的变化较为缓慢，因此图像采集的频次可以根据需求适当降低，能有效减小系统功耗。此外，图像能反映出很多信息，采集的图片可以配合病虫害识别等算法进一步提升应用价值。

1 标尺图像采集

由于田间拍摄背景复杂，红色和蓝色相对易于区分，有利于图像分割，故本文选用红蓝相间的定制标尺。参考标尺精确到分米级，通过像素高度和实际高度的比例计算厘米级。

由于镜头的透视效应，拍摄的照片会出现近大远小的效果。为降低该影响，利用水平尺将标尺竖直放置且保持相机水平拍摄。算法的关键是识别作物冠层顶部遮挡标尺的位置，为保证图片中该位置高度与作物实际株高近似相等，标尺要放置在作物主茎秆后方，且标尺和作物主茎秆之间距离尽量小。保持相机中心与标尺中心对齐并能完整捕捉标尺的图像，调整相机位置至标尺在图像正中间位置。本文选用的标尺是1.2 m，因而量程为0～1.2 m，相机位置过远会导致每个像素点对应的实际高度变大，从而降低了精确度；相机位置过远会导致标尺色块丢失。因此，要使标尺占图像高度的2/3以上，但不能超出图像。

2 株高测量算法

2.1 算法流程

本文算法的核心是通过颜色阈值分割出未被作物遮挡的标尺并计算其在图像中的像素高度。株高测量算法流程图如图1所示。第1步要分割标尺色块和复杂的背景，将图像转换到HSV颜色空间进行阈值分割，保留符合标尺颜色特征的像素点并将背景置黑，这样就初步获得了标尺色块。第2步要对标尺色块进行优化去噪，滤波过程采用中值滤波和去除最小连通域的方法，随后进行一次闭运算防止色块出现断裂。第3步计算株高，采用连通域分析法计算未被遮挡的标尺色块的像素高度，通过标记连通域获取各个色块的位置和大小信息，筛选出未被遮挡的所有标尺色块进行计算。

图1 株高测量算法流程图

2.2 误差补偿

由于光线是直线传播的，所以在图像上反映出的株高会有一定的误差ΔH。如图2所示，标尺和作物主茎秆之间距离为a，固定相机位置后测量相机到标尺之间的距离为b，未被遮挡的色块总个数为n，通过标尺图像计算出作物遮挡位置到标尺中间位置的距离H，遮挡部位色块的像素高度为n。由于作物遮挡位置高于标尺中心点时计算出的株高偏高，此时H为负值，可得

图2 误差补偿原理

H=0.1n-0.7+h0

(1)

通过三角形相似即可计算出ΔH。

(2)

2.3 阈值分割及滤波

由于标尺图像在完整图像中只占中间的一小部分，所以在相机固定后，本文裁剪出图像宽度的1/3，取中间部分进行图像处理。这样做既能减少不必要的运算，也能剔除部分无关因素的影响，在一定程度上提高了算法的工作效率和鲁棒性。在MATLAB上可使用imcrop函数完成图像裁剪。

相机采集的图片为RGB格式，RGB颜色空间最大的优点是直观，但R、G、B3个分量高度相关难以分离。本文将图片转换到HSV(Hue，Saturation，Value)颜色空间中进行颜色分割[11-13]。Hue分量(H)表示色相，取值范围为0°～360°。Saturation分量(S)表示饱和度，即颜色接近光谱色的程度，可以理解为混入白色的量，当Saturation分量为零时，只有灰度。Value分量(V)表示明度，即颜色明亮的程度，可以理解为混入黑色的量。在HSV颜色空间中，Hue分量相对独立，在颜色分割中易于实现，在HSV颜色空间中进行颜色分割可在一定程度上降低光照强度对分割效果的影响。Hue分量0°附近为红色，240°附近为蓝色。

RGB图像转HSV可表示为

(3)

在MATLAB上使用函数rgb2hsv将RGB图像转换到HSV颜色空间。

完成颜色空间转换后进行图像分割，确定分割阈值，统计出HSV图像中H分量的直方图。将用来标定的样本图片转换到HSV颜色空间后所统计的H分量的直方图如图3所示。在MATLAB中，H、S、V3个分量都是归一化的，取值范围是(0，1)。对色相环(0°，360°)进行归一化，0°附近为红色，240°附近为蓝色，归一化后H分量在0的附近和1的附近为红色，在2/3附近为蓝色。

图3 H分量直方图

对多张图片的直方图进行对比，反复进行分割实验，观察分割出的标尺色块的完整度及非标尺色块的数量，最终确定红色分量阈值为(0，0.04)∪(0.96，1)，蓝色分量阈值为(0.5，0.75)，此时分割效果最好。遍历像素点，保留目标阈值范围内的图像，将超出阈值范围的像素点置为黑色，使得复杂的背景转化为纯净的黑色背景。如图4所示为标尺红色色块的分割效果，图中可见存在一定的噪声。

图4 标尺色块分割效果

在MATLAB上采用for循环遍历各个像素点两次，分别保留红色和蓝色像素点，将其余像素点全部置黑，保存为两个图片分开处理。

接下来需要对图像进行二值化，本文选用Ostu算法进行二值化操作。Ostu算法，即最大类间方差法，是利用阈值将图像进行二值化的高效算法。对于待处理图像I(x,y)，二值化分割阈值记为T，目标和背景像素点数占图像比例分别记为ω0和ω1，平均灰度记为τ0和τ1，类间方差记为D，图像平均灰度记为τ，则有

(4)

阈值T为灰度值，取值范围为(0，255)。遍历阈值T，取类间方差D最大时T的值为分割阈值。

在MATLAB中用函数rgb2gray将分割出的图像转换成灰度图像，函数graythresh即是利用Ostu算法求灰度图像的分割阈值T，使用函数im2bw以T为分割阈值对图像进行二值化操作，得到标尺色块为白色的二值图像。

由于阈值分割后的图像中存在许多椒盐噪声，本文选用中值滤波进行处理。中值滤波的基本原理是把数字图像中一个点的像素值用该点邻域中各点像素值的中值代替，是基于排序统计理论的非线性滤波。中值滤波消除椒盐噪声效果显著，同时它能够保护信号的边缘，使之不被模糊。本文中体现标尺色块的长和宽的边缘信息是计算作物株高的关键，中值滤波正符合本文需求。

首先取一个3×3的模板，将该模板遍历整个图像，对模板中所有元素进行排序，取中间值替换模板的中心位置。本文中对二值图像进行中值滤波，由于二值图像中只有0和1，所以排序操作可以简化为统计1的个数i。在3×3模板中，如果i≥5则取1放入模板的中心位置，反之取0。如图5所示为举例对比灰度图像和二值图像的中值滤波。

(a) (b)

在MATLAB中可使用函数medfilt2(I，[3，3])对图像I进行中值滤波，其中[3，3]指滤波模板的大小为3×3。

本文采用标记连通域的方式获取标尺色块位置信息，为防止标尺色块因光线影响出现断裂从而影响连通域标记，本文对图像进行一次闭运算。闭运算是对二值图像先膨胀再腐蚀，作用是填充物体内的小空洞，连接邻近的物体，平滑其边界的同时不改变面积。因为本文中标尺色块的边界是平滑的，所以平滑边界的影响较小。

在MATLAB中使用函数bwmorph(BW，operation)对二值图像BW进行形态学处理。operation是一个字符串，用于指定进行的形态学处理类型，闭运算operation 为'close'。本文采用默认结构元对图像进行操作，实际应用中可能出现未将断裂的色块连接起来的情况，此时可根据需要相应地调整结构元的大小，改用imclose函数进行闭运算。

2.4 连通域分析

连通区域分析是指将图像中的各个连通区域找出并标记[14-15]。彼此连通的像素点形成的一个区域即为连通域，标记连通域是二值图像分析的基础，通过对白色目标像素的标记，给每一个单独的连通域打上不同的标签，进而可以获得每个连通域的轮廓、面积、外接矩形等参数。图像中每个像素点周围有8个相邻的像素点。常见的邻接关系有4邻接和8邻接。4邻接为目标像素点上、下、左、右4个相邻的像素点；8邻接为目标像素点周围邻接的所有像素点。如果像素点A与B邻接，则称A与B连通，连通域连通关系类似于等式。如果A与B连通，B与C连通，则A与C连通。

在二值图像中，目标标尺色块面积较大，经过闭运算后一些小的噪点会连成一些较小的连通域，进行去除小连通域操作，剔除小连通域。

在MATLAB中使用函数bwareaopen(BW，P，conn)删除二值图像BW中面积小于P的对象。conn指的是连通域邻接方式，默认情况下使用8邻域。连通区域标记函数bwlabel中使用的算法是：一次遍历图像，并记下每一行(或列)中连续的团(run)和标记的等价对，然后通过等价对对原来的图像进行重新标记。返回一个和原图像大小相同的矩阵，包含了标记原图像中每个连通区域的类别标签，这些标签的值为1，2，…，num(连通区域的个数)。

首先逐行扫描图像，把每一行连续的白色像素点组成的序列称为一个团，并记录下起始、结束位置和行号，合并相邻两行有重叠区域的团并提取等价对。举例说明，如图6所示为1个4行16列的二值图像。

图6 4行16列的二值图像

第1行有两个团，即[3，8]和[12，14]，标记为1和2。第2行有两个团，即[5，9]和[14，16]，与上一行的团有重叠区间，所以沿用上一行的团编号。第3行有一个团[9，14]，与上一行的两个团都有重叠区域，所以给他标记为较小的团编号1，然后将(1，2)写入等价对。第4行有一个团[1，5]，与上一行的团没有重叠区域，单独标记为3。

将等价对装换为等价序列，例如某二值图像有12个团编号，等价对有(1，2)，(1，3)，(2，7)，(7，5)，(4，6)，(6，9)，(10，1)，(11，9)，(11，12)。将1～12都看做节点。第1条等价序列从1开始，有(1，2)，(1，3)和(10，1)，则1、2、3、10为等价序列，再从2开始，有(2，7)，7又通向5，可知1、2、3、5、7、10为等价序列。同理可得第2条等价序列：4、6、9、11、12。查找完毕后，8未出现在等价对中，说明8区域没有等价区域，所以8单独形成1条等价序列。

得到等价序列后遍历图像重新标记，将等价序列标记为同一个团编号，并记录下该团在图像中最左、最右、最上和最下端的像素点坐标，这样就可以得到关于目标标尺连通域的坐标数据。将这些坐标数据作差即可轻松求得各个目标连通域的像素宽度和高度。将这些数据存为一个5列的矩阵，每一行为一个连通域的数据，第1列为连通域编号，第2列为连通域的列最小值，第3列为连通域的行最小值，第4列为连通域的像素宽度，第5列为连通域的像素高度。

获得连通域数据后需要筛选出未被作物遮挡的标尺色块，因为其中还混杂着少量非目标色块。方法如下：

(1)起始列位置滤波。由于标尺色块理论上是在一条竖直直线上的，反映在图像中则是标尺色块的起始列相近，在允许误差下设定阈值可以滤除位置偏差较大的非目标色块，该阈值X0可以通过设备架设后的标定过程获得。取矩阵的第2列与阈值X0对比，删除距离第X0列较远的元素所在行；

(2)宽度滤波。在不移动设备的前提下标尺色块的宽度相对固定，通过标定结果设定一个宽度阈值ΔX滤除宽度过小或过大的非目标色块。原理同上，取矩阵的第4列与宽度阈值ΔX对比，删除宽度过小或过大的元素所在行。

完成对图像中的非目标色块的滤除后，要在各标尺色块中筛选出作物遮挡部位的色块。由于该色块有一部分被遮挡，色块的像素高度必然比完整的色块小。假设红色和蓝色色块中像素高度最小的色块为作物遮挡部位的色块，并分析该假设可行性：(1)若该色块不是作物遮挡部位的色块，说明遮挡位置刚好在色块分界处，所有色块都是完整的，而后续处理会将该色块误认为是作物遮挡部位的色块，从而出现误差，但该误差根本来源在于前端采集和阈值分割部分，在允许范围内；(2)若该色块是作物遮挡部位的色块，假设成立。

综上，该假设对最终结果影响较小。根据目标色块像素高度对连通域数据进行降序排序，即以矩阵第5列元素为关键元素对连通域数据进行降序排序，在MATLAB中函数sortrows(A，column)表示对矩阵A按照第column列进行排序，其中参数column为正值时为升序排序，column为负值时为降序排序。统计色块的总个数n。再次遍历图像保留作物遮挡的色块并将其他位置置黑，即将该色块最小外包矩形以外的像素点置黑。由于目标色块已经提取出来，为减少计算量可再次对图像进行裁剪，只保留有用的部分。处理结果如图7所示。

图7 作物遮挡部位色块的提取

(5)

识别出对应0.1 m的标尺色块可取的像素高度最大值为样本Temp。有作物遮挡的色块高度hs为

(6)

标尺未被遮挡部分的总长度Δh为

Δh=0.1(n-1)+hs

(7)

作物株高为High为

High=120-Δh+ΔH″

(8)

结合式(1)、式(2)、式(5)～式(8)即可计算出株高High。

2.5 参数标定

前文中提到很多参数需要在设备架设后进行标定。操作时，使标尺竖直放置，保持摄像机水平并与标尺中心位置对齐，记录下标尺和作物主茎秆之间距离a以及相机到标尺之间的距离b。此时拍摄4～5张图片，运行算法，通过HSV色彩空间3个分量的直方图确定分割阈值。通过色块数据矩阵第2列数据可以检查标尺在图像中是否竖直放置，否则调整标尺或相机角度。在没有明显倾斜的情况下，记录此时色块数据矩阵第2列数据列最小值的平均值X0以及色块数据矩阵第4列数据色块像素宽度的平均值ΔS，然后分别计算出色块数据矩阵第4列和第5列数据的极差并记录下较大的极差ΔS，取ΔS为识别允许的最大误差。

2.6 基于MATLAB的GUI设计

图形用户界面(Graphical User Interface，GUI)是MATLAB的可视化操作功能，但在实际使用中MATLAB仿真不够便捷。利用前端相机采集图像，发送到PC端，利用GUI打开保存的图像并计算作物株高。将数据可转化为exe格式的可执行文件，即可脱离MATLAB环境运行，方便用户使用[16-19]。

3 实验及误差分析

本文的测量对象为小叶黄杨和金叶女贞绿篱，均为某公司内景观植物，需要定期修剪整齐，本实验对其株高进行测量可以指导修枝的时机。选用的标尺的每段色块都是0.1 m，一共有12个色块，总长1.2 m。将标尺架设在灌木中，由于识别方法是通过颜色分割标尺色块，所以拍摄背景不可有大量红色和蓝色，天空蓝色较浅对实验影响不大，相机采集的图片大小为3 024×4 032。固定标尺和相机位置，首先拍摄3～4张图片完成参数标定过程，随后将标定的参数设置在MATLAB算法中。图像采集时间为2020年9月25日，天气晴，分别在9∶00、14∶00和17∶00拍摄。选择拍摄时间的原因是在9∶00和17∶00时光线较为柔和，在14∶00时光线较强，可以对比得出光线强弱对本方案的影响。由于灌木生长密集，相对容易测量，故可省略株高补偿过程。对比数据采用人工测量的方法，人工测量采用卷尺直接测量，可精确到0.001 m，估读1位。将标尺架设在同一地点，选择3个不同相机位置拍摄，分别拍摄3个不同机位的图片，实际拍摄距离为1.7～2.2 m。测得株高数据如表1所示。分别拍摄标尺架设在不同灌木丛中的图片，测得株高数据如表2所示。

表1 不同机位测得株高数据

表2 不同灌木丛测得株高数据

表1和表2测量结果表明在正常日光照射下测量结果不会出现太大偏差。表1为相同待测点不同拍摄位置测出的株高数据，测量结果表明按标定要求架设相机拍摄的图片经图像处理测得的株高数据是可信的，不同相机位置并不会带来显著的偏差。表2是标尺在不同待测灌木丛拍摄的图片测出的株高数据，由表1可知不同相机位置并不会带来太大的偏差。表2用标尺测不同灌木丛高度时，即可近似认为相机位置改变对结果无影响，此时人工测量结果表明该实验是可复制的，在不同待测点测得的数据误差均可接受。农业气象试点站的《作物气象自动观测站测试评估大纲》中对自动观测株高的标准如表3所示。

表3 作物冠层高度准确性评估评分标准

实验测得的株高数据与人工测量的数据相比偏小，误差在0.017 3 m以内，达到一级标准。由于标尺的存在，本文算法的误差在测量作物幼苗时期相对较大，但在测量较高的作物时表现较好。本文算法的误差主要来自于以下3个方面：一是相机镜头的畸变，虽然使用水平仪使得拍摄角度带来的透视效应尽可能的小，但是相机镜头本身带来的畸变无法调整；二是阈值分割时标尺色块边缘缺失，标尺色块的边缘信息是丈量色块长和宽的关键，在复杂的田间环境中对色块进行阈值分割时会出现边缘缺失，此时对测量结果影响较明显；三是在人工测量时，作物冠层参差不齐，读数时读取最高叶片顶部的高度，而在图像处理过程中是以作物遮挡色块的底部边界的平均高度作为株高，所以测得的数据与人工测量的数据相比偏小。

4 结束语

本文提出了一种图像处理配合标尺测量作物株高的方法。通过对田间架设的定制标尺的照片进行颜色识别提取出标尺色块，将图像二值化以后进行一系列滤波将非目标色块滤除，采用连通域分析法计算出图像中未被遮挡的标尺的像素高度，最后通过比例关系计算株高。经实验数据验证该方案测得的株高相对人工测量株高的误差小于0.017 3 m，可实现测量株高的要求。该方法采用MATLAB编写GUI进行图像处理，并转化为exe格式的可执行文件脱离MATLAB环境运行，方便用户使用。为减小测量误差，可以针对误差来源对算法进行改进，例如引入相机畸变矫正算法解决相机镜头自身带来的畸变问题。

本文的方案也有地方需要改进，在作物生长初期或作物相对稀疏时，标尺、作物与相机三者之间的相对位置对测量影响较大，目前解决办法是标尺要放置在作物主茎秆后方并尽量靠近作物主茎秆。虽然本方案对设备架设和标定过程要求较高，但本方案设备简单，只需要一套图像采集系统和一个定制标尺，节约了设备成本。此外，图像采集系统还可以配合其他农情观测系统使用，例如作物病虫害观测等。本方案后期投入少，标尺和相机架设完成后只需定期采集图片即可，节约了劳动力成本。