1 引言

大豆作物在生长过程中所体现出的形态特征对选育大豆优良品种具有重要作用^[1]。株高的动态变化不同程度地影响大豆产量,研究该表型特征不仅有助于基因-环境相互作用的定量分析,而且对提升大豆品质至关重要^[2]。

传统株高测量方法是利用直尺、手持式激光测距仪等设备手动测量冠层最高点到植株根部的垂直长度,这种方法不仅费时费力,而且结果易受测量者主观因素的影响。为了弥补传统测量技术的不足,信息化的测量技术被广泛地应用到农作物表型特征的研究中。

目前植物表型参数测量研究大多局限于实验室环境或是基于虚拟植株三维重建,无法反映田间植株真实的生长情况^[3-4]。基于田间的高通量表型参数测量手段是制约植物基因组研究和作物改良的主要因素^[5]。目前农作物冠层信息获取技术主要包括手持式传感器、基于双目相机的立体视觉技术、手持式激光扫描技术和三维激光扫描及雷达技术等^[6-7]。手持式传感器采用的是接触式测量方法,会导致农作物发生形变。基于双目相机的立体视觉技术^[8],可对作物冠层图像进行非接触式无损采集,且不破坏作物生长形态,但采集的图像质量受光照影响较大,进而降低了表型参数计算的准确性。手持式激光扫描技术(如FastSCAN型激光扫描仪和Artec EVA型激光扫描仪)虽然能够高精度地重构作物冠层的三维结构形态,但其采集速度较慢^[9]。近年来,在地物表型参数分析的基础上,研究人员提出了空间上高通量的测量方法,如网络相机远程表型测量^[10],这种方法虽然实现了农作物表型特征的快速采集,但其测量结果误差波动范围较大。三维激光扫描及雷达技术^[11-12](如FARO Focus3D 120型地物激光扫描)能快速精确地获取植物冠层点云信息,但获取的点云信息中存在大量的冗余背景信息,需要进行点云筛选等操作才能精确构建冠层三维结构形态,且昂贵的价格进一步限制了该设备的广泛应用^[13-15]。上述技术手段已经在大田作物(如玉米和水稻)和树木等的表型参数计算中取得了较多的研究成果,但针对大豆的表型研究相对较少。

本文提出了一种价格相对低廉且精确度相对较高的农作物冠层信息获取方法。Kinect传感器具有较精确和快速获取空间信息的优势,已经被广泛用于农产品检测、农作物实时监控和植物表型参数获取等研究中^[16-18]。在此构建了基于Kinect 2.0传感器的大豆冠层图像采集平台,以快速高效地获取自然环境下大豆冠层彩色及距离点云数据,并以此开展大豆株高的计算方法研究。该研究成果为田间农作物高通量表型技术的研究提供了一种有力的获取手段,为大豆优良品种的选育提供了精确的数据支持^[19-20]。

2 材料与方法

2.1 实验材料

大豆种植及其冠层信息的采集工作是在黑龙江八一农垦大学农学院卓越农业人才创新创业训练园开展的。在自然光强为0~1200 μmol·m^-2·s^-1、温度为20~34 ℃的室外条件下,选取抗线9号、抗线13号和富豆6号这3个品种的大豆作为研究对象,其中抗线9号和抗线13号大豆的抗线虫能力较强,富豆6号富含丰富的油脂。将种子经过精选后进行消毒、催芽,播种于花盆中。具体操作为:首先将中等土块铺平于直径为0.3 m、高为0.18 m的聚氯乙烯(PVC)材料的盆底,装入筛选好的非盐碱细土,直至盆重为 4.5 kg,均匀平铺复合肥后继续在盆中放入细土至7.5 kg。在2018年5月21日种植第一批上述三个大豆品种,分别在9 d后和18 d后种植第二批和第三批大豆。每个品种的大豆各15盆,三批大豆共计45盆。

2.2 图像同步采集系统构建

构建基于Kinect 2.0的大豆冠层三维数据同步采集系统,动态采集大豆冠层彩色图像和深度信息的多源数据。该系统由Kinect 2.0相机、笔记本电脑、线圈和可上下左右调节的铁架组成。通过云台把Kinect 2.0相机安装在可调拍摄架的固定位置处,通过USB 3.0接口连接Kinect 2.0相机和笔记本电脑。Kinect 2.0由彩色摄像头(RGB Camera)、深度(红外)摄像头和红外发射器组成^[21-22]。其中彩色摄像头用于拍摄彩色图像,其分辨率为1920 pixel×1080 pixel。深度(红外)摄像头用于获取红外光谱并创建可视范围内的物体深度图像,其分辨率为512 pixel×424 pixel,帧速为30 frame/s。深度数据获取范围为0.5~4.5 m,水平角度为70°,垂直角度为60°。为保证大豆正常生长,数据采集工作均在室外自然光照下进行。图1为大豆植株采集方式示意图。

图 1. 大豆植株采集方式示意图

Fig. 1. Schematic of acquisition method for soybean plants

下载图片 查看所有图片

2.3 采集方法

为减少光照对大豆冠层图像的影响,并保证所采集的实验图像的清晰度,选择垂直方式获取大豆冠层图像信息,并基于深度信息开展大豆株高计算方法研究。数据采集周期为7~10 d。为保证深度信息的有效性,将Kinect 2.0相机置于可调铁架中央,通过调节铁架高度寻找最佳拍摄高度。大豆株高人工测量时间在设备采集的前一天或后一天进行,测量工具为直尺,单位为cm,将所获取的测量值与计算值进行对比。

3 大豆三维点云快速重建

为构建具有颜色信息的单盆及群体大豆冠层三维结构形态,利用搭建的同步采集系统,以顶视法采集富豆9号、抗线9号和13号大豆冠层的彩色图像及其深度图像,获取的图像分别包括红、绿、蓝(R、G、B)通道的颜色分量值及x、y和z轴的距离值。由于彩色相机和深度相机的视场和分辨率不同,颜色信息和深度图像的距离值不能直接一一对应,因此,需要对二者进行标定。标定方法为:求解图像坐标系和空间坐标系旋转矩阵,利用所求得的旋转矩阵将彩色及深度图像转换到同一坐标系下,实现彩色图像与深度图像的配准^[23]。

以单株抗线9号大豆为例,具有颜色信息的大豆冠层三维重建流程如图2所示。图2(a)为图像同步采集系统;图2(b)和图2(c)分别为在自然光环境下利用图像同步采集系统垂直拍摄得到的大豆冠层深度图像和彩色图像;图2(d)为构建的原始三维点云;图2(e)为彩色图像灰度直方图,用于设定颜色分割阈值;图2(f)为原始点云距离图;图2(g)为构建的具有颜色信息的原始点云;图2(h)为分割后大豆冠层HSV (Hue, Saturation, Value) 距离颜色图;图2(i)为分割后具有颜色信息的大豆冠层距离图。

图 2. 具有颜色信息的大豆冠层三维(3D)重建流程图。(a)图像同步采集系统;(b) 深度图像;(c) 彩色图像;(d)原始三维点云;(e) 彩色图像灰度直方图;(f) 原始点云距离图;(g) 具有颜色信息的原始点云距离图;(h) 分割后大豆冠层距离图;(i) 分割后具有颜色信息的大豆冠层距离图

Fig. 2. Flow chart of three dimensional (3D) reconstruction of soybean canopy with color information. (a) Synchronous acquisition system; (b) depth image; (c) color image; (d) raw 3D point cloud; (e) gray histogram of color image; (f) distance image of raw data; (g) distance image of raw data with color information; (h) distance image of canopy after segmentation; (i) distance image of canopy with color information after segmentation

下载图片 查看所有图片

3.1 Kinect相机标定

空间物体上任意一点(p)都与相机拍摄所得的图像存在对应关系,该p点在二维图像(O₁-u_Lv_L和O₂-u_Rv_R)中的位置经过针孔模型可还原成三维。Kinect相机标定就是求解此对应关系,将三维物体的空间坐标经过矩阵的旋转及平移运算后,映射到相机的二维成像平面(O_L-x_Ly_L和O_R-x_Ry_R)上。由于实际情况中的参照物不同,一般需要在统一的世界坐标系(O_W-X_WY_WZ_W)下表示不同相机坐标系(O_cL-X_LY_LZ_L和O_cR-X_RY_RZ_R)下的空间信息,该Kinect相机视觉模型^[24]如图3所示。

图 3. Kinect相机视觉模型

Fig. 3. Visual model of Kinect camera

下载图片 查看所有图片

为精确标定彩色相机和深度相机,首先,分别获取Kinect 相机视场下的12组不同位置的棋盘格彩色图像和红外图像,由于红外图像与深度图像具有统一的视角和分辨率,因此可用红外图像代替深度图像作为彩色图像标定的对象。利用 MATLAB标定工具箱,分别对两个相机进行标定^[25]。标定得到的棋盘格与相机的相对位置关系如图4所示。

图 4. 棋盘格与相机的相对位置关系。(a) 彩色相机拍摄的棋盘格角点检测;(b) 深度相机拍摄的棋盘格角点检测

Fig. 4. Relative position relationship between checkerboard and camera. (a) Corner detection of checkerboard by color camera; (b) corner detection of checkerboard by depth camera

下载图片 查看所有图片

世界坐标系中的坐标点经过旋转平移后,可转化到相机坐标系O_C-X_CY_CZ_C下,其数学表达公式如下所示:

XCYCZC=RXWYWZW+T=r11r12r13r21r22r23r31r32r33XWYWZW+t1t2t3,(1)

式中:R为3×3的旋转矩阵;r_ij(i=1,2,3;j=1,2,3)为R的矩阵元;T为 3×1的平移矩阵,是相机的外部参数;t₁、t₂、t₃为T的矩阵元^[26]。将(1)式转化为齐次坐标,即

XCYCZC1=RT0T31XWYWZW1,(2)

式中:0^T=(0,0,0)。 通过上述求解运算,可得

R=0.8923-0.09240.44180.43360.4474-0.7821-0.12530.88950.4393,T=-230.95365110.57616-132.76179。(3)

3.2 具有颜色信息的三维点云重建

利用上述标定关系,将彩色相机坐标系转换到深度相机坐标系下,即可得到深度图像中距离值与彩色图像中颜色值的对应关系。利用MATLAB的scatter3函数实现大豆冠层三维点云的重建,重建后的三维模型仅包括形态结构信息^[27]。为了重建具有完整颜色信息的大豆三维模型,对采集的点云颜色数据进行了均值融合处理。最后,根据相机间的标定关系,将点云数据的RGB颜色信息赋值到对应的深度图像中,最终重建了具有颜色信息的大豆冠层三维结构形态。

3.3 具有颜色信息的大豆冠层点云提取方法

为提高后续大豆株高计算的速度及准确性,需要对第3.2节重建后的三维点云进行分割,即将大豆冠层点云从复杂的背景点云中提取出来。在此采用最大类间方差阈值分割法确定分割阈值T^[28]。设(x,y,z)为距离图像中的任意一点,f(x,y,z)为相应的灰度级。若灰度值大于阈值,标记为大豆植株(用1表示),否则标记为背景(用0表示)。标记值h(x,y)可表示为

h(x,y)=1,f(x,y,z)>T0,f(x,y,z)≤T。(4)

上述T的选取对于大豆冠层的有效提取至关重要,其基本思路是以图像直方图中的某一灰度为阈值,将图像分成两组并计算其方差。当两组之间的方差最大时,将此灰度值确定为所需的阈值。图像概率分布的直方图可表示为

pi=niN,i=0,1,2,…,L-1,∑i=0L-1pi=1,(5)

式中:p_i为概率分布;N为图像的像素;i为该图像灰度级,i∈[0,L-1];n为对应灰度级i的像素。

假设T将图像中的像素按灰度值分为大豆冠层C₀∈[0,T]和背景区域C₁∈[T+1,L-1]两类,这两类具有不同灰度分布概率的整幅图像的均值可表示为

uT=∑i=0L-1ipi,(6)

式中:u_T为整幅图像的均值。

大豆冠层像素均值为u₀,背景区域像素均值为u₁,其数学表达式分别为

u0=∑i=0Tipiw0,u1=∑i=T+1L-1ipiw1,(7)

式中:w₀为大豆冠层像素均值出现的概率;w₁为背景区域像素均值出现的概率。w₀和w₁的数学表达式分别为

w0=Pr(C0)=∑i=1Tpi,(8)w1=Pr(C1)=∑i=T+1L-1pi=1-w0,(9)

式中:P_r(·)为所求区域(大豆冠层或背景区域)的概率分布。

最终求得两类间的方差u_L为

uL=w0u0+w1u1σ2=w0w1(u0-u1)2。(10)

当T的取值使得两类间的方差最大时,即

σB2(T)=max0≤T≤L-1{uT},(11)

此时T值为最佳阈值。

利用求解的T分割三维重建后具有颜色信息的三维点云,提取得到的大豆冠层如图5所示。图5(a)、(b)分别为分割前大豆植株的原始三维点云和具有颜色信息的原始点云,5(c)、(d)分别为分割后大豆冠层三维点云距离正视图和俯视图,5(e)、(f)分别为融合RGB颜色信息后大豆冠层的正视图和俯视图。

图 5. 具有颜色信息的大豆冠层三维重建。(a)原始三维点云;(b)具有颜色信息的原始三维点云;(c)分割后大豆冠层三维点云距离正视图;(d)分割后大豆冠层三维点云距离侧视图;(e)具有颜色信息的大豆冠层三维点云正视图;(f)具有颜色信息的大豆冠层三维点云俯视图

Fig. 5. 3D reconstruction of soybean canopy with color information. (a) Raw 3D point cloud; (b) raw 3D point cloud with color information; (c) 3D point cloud distance image of soybean canopy after segmentation under front view; (d) 3D point cloud distance image of soybean canopy after segmentation under side view; (e) 3D reconstruction of soybean canopy with color information under front view; (f) 3D reconstruction of soybean canopy with color information under top view

下载图片 查看所有图片

4 基于深度信息的株高计算方法研究

4.1 单盆大豆株高计算

假设地面相对平坦,即地面三维点云在z轴方向上波动较小。单盆大豆冠层三维图像采集方法如图6所示。为了提取大豆冠层三维点云图,首先将扫描得到的所有区域进行栅格划分,并将三维点云坐标投影到x-y、x-z或y-z的栅格平面上,其两两垂直的三条边分别与笛卡儿坐标系的三个坐标轴平行。根据z轴上的距离值设定阈值,去除三维点云中的背景部分后,即可提取出大豆植株冠层三维点云。由于Kinect相机获取的坐标信息是无序的,需通过改变x、y和z方向的阈值范围得到大豆冠层区域的三维点云集,进而计算每一单盆大豆的株高h,其计算公式分别如下:

xm∈X0|X1≤xm≤X2ym∈Y0|Y1≤ym≤Y2zm∈X0|Z1≤zm≤Z2,(12)

式中:(x_m,y_m,z_m)为在阈值范围内对某点(X₀,Y₀,Z₀)进行投影得到的点云坐标信息; [X₁,X₂],[Y₁,Y₂]和[Z₁,Z₂]分别为三个坐标方向上的大豆冠层阈值范围。则每一单盆大豆的株高h为

h=zmax-zmin,(13)

式中:z_max为z轴方向最大值;z_min为z轴方向最小值。

图 6. 单盆大豆冠层三维图像采集方法

Fig. 6. 3D image acquisition of individual soybean canopy

下载图片 查看所有图片

对单盆株高的计算总共进行了90次,求解平均值并将其作为最终株高值。

4.2 群体大豆株高计算

为实现高通量大豆株高计算,开展了群体的大豆株高计算方法研究。群体数据采集实验场景如图7所示。利用搭建的采集平台获取了大豆生长生殖过程中出枝期、开花期、结痂期和鼓粒期4个时期的彩色图像和深度图像,利用单盆大豆的处理方法,重建了具有颜色信息的群体大豆冠层三维结构形态,并有效地提取了大豆冠层区域,用于开展基于群体的大豆株高计算方法研究。大豆生长茂密时期植株叶片会出现遮挡现象,因此需要对相机安装高度和群体大豆植株摆放的盆间距做出动态调整。

图 7. 群体大豆株高计算方法

Fig. 7. Calculation method of plant height for grouped soybean plants

下载图片 查看所有图片

群体大豆株高的计算方法如图7所示。其中H为Kinect相机安装位置距离地面的高度;h_p为盆高(h_p=0.18 m);d为两植株中心点投影到地面的距离(图8);S为Kinect相机距离群体大豆中任意一盆大豆植株中心点顶端间的距离;h_a为除去大豆冠层的地面部分(包括地面和盆)相机距离大豆冠层的距离;θ为相机扫描点与地面的水平夹角。含下标1的参数(S₁、θ₁、h_a1和h₁)表示最左侧一盆大豆的计算参数,含下标2的参数(S₂、θ₂、h_a2和h₂)表示左侧第二盆大豆的计算参数。根据测量Kinect相机距离地面的高度确定点C(x₁,y₁,z₁),通过设定群体中每盆大豆冠层区域的范围来确定每盆大豆冠层顶端最高点三维坐标A(x₂,y₂,z₂),计算得到点C与点A之间的距离为

S(C,A)=(x2-x1)2+(y2-y1)2+(z2-z1)2。(14)

图 8. 盆间距示意图

Fig. 8. Schematic of spacing between pots

下载图片 查看所有图片

根据三角模型距离计算原理计算得到相机距离大豆冠层的距离为

ha=dtanθ=dsin(arccosθ)cosθ,(15)

其中:

d=r+0.03,(16)cosθ=Sd。(17)

从而得到群体图像中每盆大豆植株的株高为

h=H-ha-hp。(18)

5 大豆株高计算方法有效性验证

5.1 大豆冠层三维重建效果

在标定彩色图像和红外图像的基础上,遍历深度图像的所有三维点云,根据求解得到的旋转矩阵和平移矢量,将彩色图像坐标系转换到深度图像坐标系,即将颜色信息赋值于对应的点云,构建具有颜色信息的大豆冠层三维结构形态。

为清晰观察到大豆植株三维重建效果,在此以2018年7月19日获取的15盆大豆(三个品种)的冠层彩色和点云数据为例,构建其三维结构形态,三维重建效果如图9所示。图9(a)为逐一获取15盆大豆冠层数据后,对其结构形态进行重建的单盆重建图;图9(b)为大豆冠层结构形态的群体重建图。从重建效果上看,由于单盆拍摄时不存在冠层相互遮挡,因此冠层重建更加完整。从重建时间上看,重建单盆大豆三维结构形态需要0.587449 s,而群体重建15盆大豆需要1.312698 s。因此,所建立的采集平台及三维重建方法适用于重建大豆冠层三维结构形态的重建。

图 9. 2018年7月19日获得的大豆冠层三维重建效果。(a)单盆重建;(b)群体重建

Fig. 9. 3D reconstruction effect of soybean canopy on July 19, 2018. (a) 3D reconstruction of individual soybean plant; (b) 3D reconstruction of grouped soybean plants

下载图片 查看所有图片

重建采集的所有单盆大豆冠层图像,三维重建的大豆冠层能够反映生殖生长期的大豆冠层生长状态的动态变化,且重建后的冠层点云包含颜色信息,通过提供连续的数据记录能较为真实地还原植株形态和生长状态^[29]。以大豆9-5号为例,其2018年6月19日(0619)—2018年9月5日(0905)生殖生长期的冠层重建效果如图10所示。

图 10. 大豆9-5号的2018年6月19日(0619)—2018年9月5日(0905)生殖生长期的冠层3D重建效果

Fig. 10. 3D reconstruction of 9-5 soybean canopy in reproductive period from June 19, 2018 (0619) to September 5, 2018 (0905)

下载图片 查看所有图片

5.2 大豆株高计算值比较结果

以重建大豆冠层点云为研究目的,对大豆株高进行计算。分别对2018年6月6日—2018年7月19日的三个品种大豆植株冠层采用了单盆大豆株高测量方法,计算值与实测值的平均误差为0.14 cm(图11)。对2018年7月27日—2018年9月5日的大豆植株冠层采用了群体大豆株高测量方法,计算值与实测值的平均误差为0.54 cm。单盆大豆图片清晰,植株间无遮挡,冠层信息完整,但进行高吞吐量的株高计算的时间较慢,如15个单盆大豆株高所需的计算时间约为1.223895 s。而群体大豆株高测量方法可同时计算15盆大豆株高,所需的计算时间约为0.798435 s。

每种样本总数为60盆。在群体采集方式下,分别针对抗线9号、抗线13号和富豆6号三个品种进行株高计算,计算值与实测值的相关系数R²分别为0.9717,0.9730,0.9697。实验结果表明,在自然环境下基于深度信息的大豆株高测量方法能准确地进行高通量株高计算。

图 11. 三个品种大豆株高的测量值与计算值的比较。(a)抗线9号的对比结果;(b)抗线13号的对比结果;(c)富豆6号的对比结果

Fig. 11. Comparison of manual measurement values and sensor measurement values of plant height for three soybean varieties. (a) Comparison result of Kangxian-9; (b) comparison result of Kangxian-13; (c) comparison result of Fudou-6

下载图片 查看所有图片

在大豆生长末期,大豆株高的实测值与计算值均约为50 cm。基于时间序列的三个品种大豆株高的动态变化如图12所示。随着大豆生长时期的变化,观察到结痂期(2018年8月19日)叶片开始变黄脱落,大豆株高的变化呈下降趋势。结合具有颜色信息的重建效果可见,花期大豆植株长势增快、颜色青翠、叶片挺立,但分枝增多、植株间叶片遮挡严重^[30],且叶位不易确定等导致点云数据缺失,这是影响计算精确度的主要因素。

图 12. 基于时间序列的三个品种大豆株高动态变化图

Fig. 12. Dynamic change of mean value of plant height for three soybean varieties based on time sequence

下载图片 查看所有图片

由图12可知,计算值的动态变化相比实测值较平稳。但2018年8月19日的计算值与测量值误差较大。造成该测量误差的原因主要有:1)实验仪器测量与人工测量并未在同一天进行,两日自然天气相差较大;2)大风、突然降雨和作物自身水分变化情况等因素极易导致大豆冠层叶片萎蔫不够坚挺,使误差增大;3)地面局部不平坦和人为操控实验仪器导致测量系统并非是水平的^[31]。因此,在晴朗无风且地面水平环境下采集大豆冠层彩色及深度数据,能够获得更为精确的株高计算值。

6 结论

在自然环境下,利用所设计的大豆冠层图像同步采集系统,采集了抗线9号、抗线13号以及富豆6号三个寒地大豆品种的生殖时期冠层彩色图像和深度图像数据,标定了彩色相机和深度相机,重建了以大豆生长时间为序列的发芽期、出枝期、开花期、结痂期和鼓粒期5个阶段的具有颜色信息的大豆冠层三维结构形态。在此三维重建的基础上,计算得到的单盆和群体大豆株高平均误差分别为0.14 cm和0.54 cm。抗线9号、抗线13号和富豆6号的株高计算值与人工测量值的决定系数分别为0.9717,0.9730,0.9697。该研究成果为高通量计算农作物的其他表型参数提供了一种现实可行的数据获取手段。在未来工作中,将继续利用该采集平台及技术探索大豆植株其他重要的表型特征(如茎粗、冠幅和分枝角度等)的计算方法,为进一步实现大豆植株表型参数的高通量计算提供前期技术基础。