基于三维重建的秸秆地单株水稻生长形态检测研究

魏天翔,汪小旵*,施印炎,王延鹏,王凤杰,张先洁,卜俊怡,杨四军

(1.南京农业大学工学院/江苏省现代设施农业技术与装备工程实验室,江苏 南京 210031;2.江苏省农业科学院农业资源与环境研究所,江苏 南京 210014)

长江下游稻麦轮作区,普遍推行着基于秸秆全量还田的全程机械化生产方式,但小麦秸秆还田后1个月内会经历快速腐解期[1],引发微生物夺氮的同时产生还原性毒素,不利于下茬水稻的早期生长发育,体现为返青缓慢、分蘖减少、叶尖黄化甚至僵苗等负面现象[2],直接影响着最终的产量和品质[3]。诸多研究围绕促进水稻早期生长发育展开,包括调整还田方式、增施氮肥和优化水肥管理等[4];而评估这些措施的效果主要通过对比分析株高、叶长、分蘖数等水稻生长形态参数,高效的参数测量方法可以给相关研究提供强有力帮助。

水稻的生长形态参数测量经历过人工测量、二维图像分析[5]和专业数字化仪器测量[6]等阶段,这些方法或效率较低,或成本较高,或具有破坏性,难以推广使用。基于多视角RGB图像的摄影测量方法可以以较低成本实现三维重建,使用普通数码相机即可实现,将其用于植物形态参数的测量是当下研究热点[7]。根据三维重建原理的不同可分为运动恢复结构(structure from motion,SfM)[8]和空间雕刻(structure form space carving,SSC)[9]两大算法。运动恢复结构算法在提取并匹配多视角图像序列间的特征点后,通过对极约束计算目标物各点的三维空间坐标,是基于计算机视觉的重建技术中最先进的方法之一,但用于田间水稻仍存在噪声较大和重建时间较长等瓶颈问题。空间雕刻算法以目标物多视角图像及二值分割图像作为输入,将目标区域离散化后的体素投影至多视角图像,根据其在图像上是否属于植株部分而选择保留或剔除,最终保留体素构成目标的三维模型,因其具有高效稳定和易于实现等优点被研究者们用于水稻、小麦、高粱等作物的三维重建[10-12]。目前的研究多在温室内完成,可以通过挡板简化背景分割的难度,而田间水稻的生长环境中存在杂草、浮萍、倒影等干扰因素,使得基于颜色、纹理等特征的传统图像分割方法难以准确分割出田间水稻[13];而基于深度学习的方法可以在训练过程中实现自动特征提取,在田间图像的应用上能表现出优秀的抗干扰能力,显著提升分割准确度与效率[14-15]。

本研究旨在结合深度学习图像分割技术和空间雕刻技术探索出一种田间水稻单株尺度三维重建方法,以期为不同秸秆还田方式下水稻生长参数测量提供相应的技术与方法。

1 材料与方法

1.1 秸秆还田方案

试验地点位于江苏省农业科学院泗洪稻麦科技综合示范基地(118.27°E,33.36°N),属东亚季风区亚热带和暖温带的过渡区,土壤类型为砂姜黑土,作业全程机械化。前茬作物为小麦,2020年6月3日使用久保田688Q联合收割机高留茬收割,同时秸秆粉碎后全量还田。6月4日耕作整地,动力机械为久保田M1004R轮式拖拉机,旋耕作业悬挂1GQNGK-250型旋耕机(河北春耕机械制造有限公司),翻耕作业悬挂三铧犁(淮安淮犁机械有限公司)。6月7日蓄水泡田。供试水稻品种为江苏省农业科学院自主选育的‘优香梗’,2020年5月14日育苗,6月14日使用久保田SPV-8C插秧机插秧,每公顷栽插25.5万穴,每穴4～6苗,配施基肥(复合肥405 kg·hm-2,尿素81 kg·hm-2),田间管理常规。

秸秆还田方案如图1所示,每个田块面积约400 m2,田块之间使用PVC覆膜隔断区分。基于全程机械化作业的模式,设有“泡田+旋耕(SR)”“旋耕+泡田(RS)”和“翻耕+旋耕+泡田(PRS)”3种秸秆还田方式,每种方式均设置秸秆不还田的空白对照组(CK)。“泡田+旋耕(SR)”的具体操作为前茬小麦收割后,先蓄水泡田40 h以上后进行旋耕。“旋耕+泡田(RS)”即前茬小麦收割后先旋耕再蓄水泡田。“翻耕+旋耕+泡田(PRS)”的具体操作为前茬小麦收割后先犁翻,再旋耕,最后蓄水泡田。

图1 田块分布Fig.1 Distribution of the field

1.2 水稻生长形态参数测量方法

1.2.1 图像采集为了便于图像分割和水稻单株尺度参数计算,在插秧完成后立刻在田间设置单株的观察组,参考生物统计学5点采样法在田间选取位置,以每穴1苗的方式栽插水稻,每个田块设置15株。此外用同样的方式另设置部分单株水稻作参数验证组,用于和人工测量数据对比作算法精度验证。

图像采集分为2类:一类用于三维重建和参数计算;另一类用于制作深度学习数据集。图像采集设备为iPhone 8数码手机原生相机,成像分辨率为3 024×4 032。采集方式如图2所示,首先在每一个单株观察点选取合适位置放置标定板,然后依次拍摄植株顶部视角和侧方视角图片,侧方视角根据植株生长阶段和叶片遮挡情况环绕植株拍摄4～6张照片。用于空间雕刻三维重建的照片首先要保证信息量足够,对输入图像的视角没有严格要求[16],拍摄时可根据实际情况微调拍摄视角与高度,以确保清晰完整地拍摄到植株和标定板。用于深度学习的数据集囊括水稻的不同生长时期,拍摄场景包括不同的光照和天气条件,在注重类别均衡的情况下,共采集500张照片用作深度学习训练数据集。

图2 图像采集示意图Fig.2 Schematic diagram of image acquisition

1.2.2 基于深度学习的田间图像分割在深度学习领域,自全卷积网络(fully convolutional network,FCN)[17]被提出以后,为了实现像素级的端到端(end-to-end)图片语义分割,已有诸多亮点技术被提出。空洞卷积(atrous convolutions)[18]技术可以实现扩大感受野的同时不增加计算量。编码器-解码器的结构(Encoder-Decoder)旨在融合不同维度的特征,实现空间信息和全局信息同时利用,同时在尽可能减少信息损失的情况下完成同尺寸输入输出,代表性网络有SegNet[19]、U-Net[20]等经典模型。在PSPNet[21]、DeepLab系列[22]等亮点模型中出现空间金字塔模块(spatial pyramid pooling,SPP),通过引入多尺度信息,融合不同感受野的特征图,可以更好地挖掘全局上下文信息。而DeepLabv3+[22]将空洞卷积空间金字塔模块(atrous spatial pyramid pooling,ASPP)应用到编码器-解码器架构上,实现了前述亮点技术的融合,优化了物体边缘细节分割效果,目前在多个数据集上保持着最先进的分割效果。

DeepLabv3+网络结构如图3所示,网络的编码器部分由主干网络和ASPP模块组成,可以接受任意尺寸的图片输入。网络的解码器部分采用了简单但高效的方式,将编码器输出的特征经过上采样后和解码器中同分辨率的特征进行Concat融合,经过3×3卷积细化特征后上采样至原图尺寸,经过Softmax分类层得到像素级的分割图。

图3 DeepLabv3+网络结构Fig.3 Architecture of DeepLabv3+

网络采用逐像素的交叉熵作为损失函数,图片上每个像素的Softmax分类器的输出为:

(1)

式中:x、y为像素坐标;K为总类别数;α(x,y)表示Softmax输出的该位置像素对应第k个类别的值;pk(x,y)表示该像素属于第k类的概率。整个网络的损失可表示为:

(2)

式中:wl表示类别l的损失权重;pl(x,y)为该点像素属于类别l的真实概率。

使用labelme标注工具对数据集进行人工标注得到掩膜标签,在本试验中标签类别为水稻和背景 2类,按照4∶1的比例随机划分为训练集和验证集。为验证DeepLabv3+网络在田间水稻分割方面的优越性,与PSPNet、SegNet两代典型语义分割网络以及植物分割常用的传统图像分割算法ExG+Otus[11]、Grabcut[23]进行比较。采用交并比(intersection over union,IoU)、像素精度(pixel accuracy,PA)、F值3种指标来评估分割效果。交并比是语义分割领域最常用的评价指标,为真实值与预测值之间的交集和并集之比,可以判断目标的捕获程度和分割的精确程度。像素精度为标记正确的像素占总像素的比例。F值为精准率(precision,P)和召回率(recall,R)的调和平均数,在样本不均衡时可以更客观反映模型的准确率。虽然上述指标定义不同,但均可以通过表1所示的混淆矩阵(confusion matrix)计算得到。其中:TP代表真实标签是植株且被识别为植株的像素个数;FP代表真实标签是背景但被识别为植株的像素个数;FN代表真实标签是植株但被识别为背景的像素个数;TN代表真实标签是背景且被识别为背景的像素个数。各指标的计算公式见式(3)—(7)。

表1 分割结果混淆矩阵Table 1 Confusion matrix of segmentation

(3)

(4)

(5)

(6)

(7)

1.2.3 基于空间雕刻的三维重建空间雕刻技术是一种基于多视角图像的物体外型恢复方法,该方法将待重建的目标空间按照需求分辨率划分为多个大小相同的立方体,即计算机图形学中的体素化。根据每个立方体体素在多视角图像上的投影坐标,来判断其是否需要删除或保留。其具有精度较高、易于实现的特点,具体流程如图4。

图4 三维重建流程图Fig.4 Flow chart of three-dimensional reconstruction

1)相机标定:相机标定是三维重建的重要步骤,通过相机标定确定的投影矩阵是重建时映射实际三维空间点与二维图像像素点的必要前提,该映射关系如式(8)。

(8)

式中:Xw、Yw、Zw为三维空间中某点在世界坐标系下的坐标值;u、v为该点在图像平面上的坐标值;Zc为该点在相机坐标系下的Z轴坐标,标定后可计算得到;fx与fy、cx与cy分别为相机焦距的像素单位表达和相机主点在图像平面上的坐标,这4个参数由相机本身的制造工艺决定,共同组成相机的内参数矩阵;R和t分别为旋转矩阵和平移矩阵,也称为外参数矩阵,用于描述相机在世界坐标系下的位姿。由内外参数矩阵相乘得到投影矩阵P。为了保证精度,使用外置标定物的方法,采用12×9规格的棋盘格标定板,每个小方格边长5 mm,采用经典标定法[24]计算每个视角对应的参数。

2)计算包围盒与体素化:包围盒作用在于确定需要进行体素重建的空间区域,为了提升算法的计算效率,首先采用每个体素4 mm3的空间分辨率进行体素雕刻,以此计算出目标空间包围盒的大致范围,再以每个体素1 mm3的空间分辨率执行算法。

3)空间雕刻与模型着色:当包围盒体素化完成后,取每个体素的中心位置作为其世界坐标系下的三维坐标,例如原点处的体素坐标为(0,0,0)。然后遍历每个视角,计算该视角下每个体素的投影坐标,若在该视角下投影坐标属于植株部分,则保留体素,反之删除。该过程形似实体的木块或陶土的雕刻,因此得名空间雕刻算法。体素雕刻完成后,将保留体素根据投影矩阵计算其在原RGB图上的投影坐标,以多视角图像的RGB均值作为该体素的颜色值,以此获得效果更逼真、信息量更丰富的植株三维模型。

1.2.4 基于三维模型的株型参数提取提取株高、叶长、分蘖数等物理形态参数是三维模型的优势,将这些参数与人工测量值进行对比可以较客观地评估三维模型精度,同时验证其在田间水稻上应用的可行性。本试验中的参数测量示意图见图5。

图5 植株参数测量示意图Fig.5 Schematic diagram of plant shape parameters extractionA.植株三维模型3D model of the plant;B.点云分割示意图Schematic diagram of point cloud segmentation;C.骨架提取示意图Schematic diagram of skeleton extraction;D.叶片轴线拟合示意图Schematic diagram of leaf axis fitting;E.叶长计算示意图Schematic diagram of leaf length calculation.

1)株高:通常作物的株高可以通过计算三维模型Z轴坐标极差来获得,但水稻的生长环境决定了会有部分植株没于水面下方而无法被相机采集到,因此需要在拍摄照片之前测量标定板平面相对于土壤层的高度(h0),在模型生成后,Z轴坐标的最大值(Zmax)为植株在标定板平面上方的高度数据,将h0和Zmax相加得到植株的完整株高参数。

2)分蘖数:根据水稻分蘖发生在近地面处的特性,对植株底部2 cm的体素进行删除后将模型转化为点云,然后使用基于欧式距离的点云聚类算法实现植株分蘖的分割,图5-B中不同颜色部分即为聚类得到的不同簇,簇的数量即为分蘖数。

3)叶片数和叶长:参考Fang等[11]的方法并加以改进,利用骨架提取的方法实现叶片数和叶长的计算。立方体体素在空间中根据公共点、公共边、公共面的不同,存在着26邻接、18邻接和6邻接3种关系,采用拓扑细化法可以得到拓扑结构不变的单连通骨架[25]。具体实现如下,将属于植株模型的点定义为前景点,反之定义为背景点,依次遍历植株模型中的每个点,同时满足如下4个条件则删除该点:①该点6邻域中至少有1个背景点;②该点的26邻域中必须有2个或以上前景点;③该点26邻域的前景点集合中,任意2个点可以通过集合中其他前景点形成26连通路径;④该点6邻域中的背景点,可以通过其 18邻域中的背景点形成6连通路径。迭代上述步骤至保留体素点数量不变。

骨架提取完成后,根据其26邻域单连通的特性,可以实现分支点和叶尖点的识别,满足26邻域内超过2个前景点的即为分支点,满足26邻域内仅有1个前景点的为端点,叶尖部分和植株底部都会出现端点,故还需要对每个分支点沿Z轴坐标向下搜索,将得到的端点删去后才得到最终该点叶尖点数量。此外,尽管拓扑细化法在理论上可以得到单连通骨架,但因模型边缘存在少量噪声,导致提取出的骨架偶尔出现伪分支和毛刺[16],进而引起特征点的错误识别,可人工筛选剔除错误点,最终效果如图5-C所示。

叶尖点的数量即为当前植株的叶片数。叶尖点和其直接连接的分支点的骨架为叶片骨架,对其进行降采样后可得到中心轴点,如图5-D和5-E所示,利用B样条(B-spline curve)曲线拟合叶片中心轴线后可通过曲线积分计算各叶片的长度。

2 结果与分析

2.1 背景分割结果

深度学习算法均基于Tensorflow(v1.14.0)+Keras(v2.3.1)框架实现,硬件配置为处理器Intel Core i9-9900k,内存32 GB,显卡NVIDIA RTX 2070s,操作系统Ubuntu 20.04。训练超参数设置如下:优化器选择自适应算法Adam(adaptive moment estimation);初始学习率为0.001,当验证集误差连续2个epoch不下降时,触发系数为0.5的指数衰减;epoch设置为100,同时设置验证集误差连续6个epoch不下降时触发早停(early stop);Batchsize设置为4。

5种方法的分割性能对比见表2,可以看出基于DeepLabv3+的分割算法性能最优,对本试验任务的测试图像平均分割性能为IoU值0.801,PA值0.986,F值0.822,均为所有方法中最高。

表2 不同网络性能比较Table 2 Performance comparison of different networks

图6显示了6种田间图像场景在不同方法下的分割效果,图6-A—C中的植株在拍摄前进行过人工的杂物清理,目的是简化背景以检验各算法在水稻分蘖期不同阶段的分割效果;图6-D—F为实际水稻生长过程中常见的场景,在水稻插秧后至除草前,田间会常见如图6-D所示的杂草,而除草期间需要保水至一定高度,又会出现图6-E所示的浮萍和图6-F倒影等干扰图像分割的因素。图中的人工分割结果由labelme标注后输出,用作各方法分割性能指标的计算。可以看出:本文选取的DeepLabv3+方法在水稻的不同生长时期都可以达到较好的分割效果,且应对田间各种干扰因素也能表现出较好的稳定性。SegNet作为经典的语义分割网络能实现一定程度的水稻植株分割,但在边缘和细小连接处等细节部分处理较为粗糙,例如部分茎秆和叶片连接处出现断裂,整体效果逊色于拥有为细节和小物体分割而设计的ASPP结构的DeepLapv3+网络[26]。PSPNet得益于其金字塔池化模块的设计,对物体边缘的处理效果要优于SegNet,但其网络结构中处理特征图采用的是池化操作,而不是DeepLabv3+中采用的空洞卷积,虽融合了上下文信息但造成了一定程度的降采样[22],应用于水稻这类以细长叶片和茎秆为主的目标时,会导致背景容易过拟合,即过多地将属于前景的水稻错误识别为背景。

相较于深度学习的方法,传统的植物图像分割方法最大的缺点在于难以应对环境的干扰;相比于深度学习的自动特征提取,传统算法往往需要人为的设计特征,而田间环境复杂多变,难以设计出通用的方法[27]。如温室植物图像分割常用的ExG+Otus算法,其核心在于图像灰度化时使用超绿分量突出植物的绿色特征,导致执行分割任务时不仅会保留背景中的杂草、浮萍等绿色干扰物,还会将植物体上的非绿色部分当作背景舍弃。而基于图论的Grabcut算法虽然在黄瓜、烟叶等阔叶植物叶片的分割上有良好的效果,但应用于水稻这样的窄叶多分蘖植物就会暴露出其对复杂边缘处理效果差的缺点,且其需要交互式人工抠出待分割物体的大致轮廓,在水稻图像的分割上仍需要不小的工作量[28]。综上,采用DeepLapv3+网络可以很好实现复杂环境下田间水稻图像的分割任务,节约大量的人工成本,显著提升整体算法的自动化程度。

2.2 三维模型重建结果

图7展示了不同侧方视角数量下水稻的3个生长时期的重建结果。当视角数量小于4的时候,模型较为粗糙,重影噪声较多,同时存在相邻叶片黏连等现象。这是由于植株叶片间存在遮挡,少量的视角难以满足重建要求的信息量,造成部分背景体素未被删去;视角数量过多时则会造成模型失真,具体表现为在植株较细的叶片或茎秆处出现断裂,这是因为相机的标定难以做到完全零误差,视角数量过多会造成误差的累积,误删除属于植株的体素。故视角数量的确定要根据植株结构的复杂程度而定,从试验效果看,在水稻生长的苗期和分蘖中期,选用4个侧方视角的方案,在分蘖中期过后,选用6个侧方视角的方案效果最佳。

图7 不同视角数量的水稻植株三维重建结果(n为侧方视角数量)Fig.7 Reconstructed results of rice plant using different number of views(n is the number of side views)

2.3 模型精度验证

图8为4个水稻株型参数与人工测量的对比结果。在株高和叶长散点图中,各数据的散点分布主要集中在参考直线附近,表明人工测量值与模型输出值具有较强的相关性。此外,可以发现大部分样本的株高模型提取值低于人工测量值,造成该现象的可能原因是侧方视角的标定误差累积,使空间雕刻过程中误删了部分植株顶端的细小部分,导致株高输出值的误差,但从整体上看本方法重建的植株模型具有较高精度。

图8 模型精度验证Fig.8 Accuracy validation of models

分蘖数的计算结果RMSE为0.82,与人工测量值进行线性拟合得到R2为0.89。进一步统计可得,误差为0的占总样本的58%,误差在±1以内的占总样本的88%,表明本方法在田间水稻分蘖数的自动获取上具有较高精度。叶片数的计算结果RMSE为1.39,与人工测量值进行线性拟合得到R2为0.95。其中,误差为0的占总样本44%,误差在±1以内的占总样本的74%,但存在少量(4%)误差在±4以上的样本,这可能是由于水稻生长中后期,存在部分凋零和浸泡于水中的叶片,导致模型难以精确分割和计数,但从整体上看仍具有较高精度。此外,在水稻田中人工测量参数时,难免因视角和环境干扰等因素产生主观误差,而本方法测量的数据更具客观真实性。

2.4 不同秸秆还田方式下水稻早期生长形态数据分析

使用上述方法对6块试验田中的单株观察组进行了为期4周的生长表型参数采集与测量,结果见图9。第1周3种秸秆还田方式下的株高、叶片数等参数增量均低于空白对照,部分田块甚至出现负增长,意味着水稻的缓苗和返青过程延长,原因在于泡水秸秆腐解的前10～20 d会出现微生物夺氮,同时会生成有机酸、硫化氢等还原性毒素,不利于水稻苗期的生长发育。此外,也可以看出配施基肥后,没有出现僵苗等严重后果。而在第2至4周的分蘖期中,3种秸秆还田方式下的水稻分蘖数均多于空白对照组,这有利于产量增加,可能原因是这一时段还田秸秆腐解逐步释放氮、磷、钾等营养物质,对土壤的肥力给予一定程度的补充,促进了水稻分蘖的早生快发。

图9 不同秸秆还田方式下水稻早期株型数据对比Fig.9 Comparison of early plant shape data of rice under different straw returning methodsSR:泡田+旋耕Soaking+rotary tillage;RS:旋耕+泡田Rotary tillage+soaking;PRS:翻耕+旋耕+泡田Ploughing+rotary tillage+soaking.CK:空白对照Black control.

特别值得注意的是,从生长参数的增量对比来看,水稻受秸秆还田的负面影响从小到大的处理依次为PRS、RS、SR。相比于RS和SR方式,PRS方式下的田块前期受到的抑制效果较小,具体体现在从第2周开始,秸秆还田组的株高、分蘖数、叶片数增量均要高于空白对照组,而其他方式要在第3周才出现该效果,其原因在于翻耕相对于旋耕有更深的耕深,可以将秸秆翻入更深的土层,在利于秸秆分解养分释放的同时,还可以有效增大土壤孔隙度,减小根系拓展阻力,促进水稻生长发育。而RS方式是在田块蓄水前干旋,相对于SR的水旋可以让秸秆得到更进一步的破碎,使秸秆腐解更快,从而使前2周的抑制效果要稍小。此外,叶长变化未看出明显规律,这可能和品种的基因性状有关,需要进一步挖掘分析,也体现出高通量表型技术给农艺性状信息挖掘带来的挑战性。

3 结论

本研究结合深度学习图像分割技术探索了一种田间水稻单株尺度三维重建方法,利用三维模型提取水稻早期生长形态参数,并应用于秸秆还田的效果评估,主要结论如下:

1)使用DeepLabv3+网络,通过深度学习的方法可以实现田间水稻图像的分割,其性能优于传统图像分割方法(ExG+Otus、Grabcut)和前代分割网络(SegNet、PSPNet),性能表现为IoU值0.801,PA值0.986,F值0.822。DeepLabv3+网络在水稻分蘖期的不同阶段均能实现较好的分割效果,且能克服田间的杂草、浮萍、倒影等干扰因素。

2)使用空间雕刻算法可以实现田间水稻的三维模型重建,根据植株三维模型计算得到株高、叶长、分蘖数和叶片数4个形态参数,与人工测量参数对比,决定系数(R2)分别为0.99、0.95、0.89、0.95,均方根误差(RMSE)分别为1.03、1.19、0.82、1.39,该方法结果可靠且操作简单,易于推广,可用于大田环境下的秸秆还田效果对比分析。

3)利用本文探索的方法对不同秸秆还田方式下的水稻进行了为期4周的参数测量,分析发现,不同处理方式对秸秆还田后的水稻生长有不同的影响,从生长参数的增量对比来看,PRS方式下田块前期受到的抑制效果较小。