基于运动恢复结构的玉米植株三维重建与性状提取
2020-06-29梁秀英周风燃许锡晨杨万能
梁秀英 周风燃 陈 欢 梁 博 许锡晨 杨万能
(1.华中农业大学工学院， 武汉 430070； 2.华中农业大学作物遗传改良国家重点实验室， 武汉 430070；3.华中农业大学植物科学技术学院， 武汉 430070)
0 引言
玉米是世界三大粮食作物之一[1-4]。在玉米生长过程中，株高、茎粗、叶面积等形态学参数决定植株对养分的竞争和获取，进而影响玉米的产量[5-7]。叶投影面积、叶夹角等参数直接影响玉米植株对光能的吸收[8]，叶片数直接反映了整株玉米的地上部生物量[9]，单株最小包围盒体积反映了玉米植株的空间紧凑度，单片叶最小包围盒体积、叶周长、叶投影宽度、叶投影长度等参数直接影响叶片的空间分布[10]，而玉米叶片的空间分布状态变化能够较大程度地影响玉米植株的光能利用率，进而影响玉米的生物量和产量[11]。因此，及时准确获取玉米植株的株高、茎粗、叶面积、叶片数、叶夹角等形态学参数是研究玉米植株的必要条件。传统的玉米植株形态学参数测量主要依靠人工或二维图像测量。人工测量速度慢、成本高、主观误差较大[12-13]，并且对植株有损伤，所以不适合玉米植株的连续测量[14]；基于二维图像的测量往往由于植株间的遮挡导致精度不高，甚至由于维度的限制不适合用于测量叶面积等形态学参数[15]。
近几年，为获得更高的测量精度，基于三维重建的方法被用于目标性状提取[16-20]。运动恢复结构(SfM)是一种基于多视图几何基本原理的三维重建技术，具备自校正、受环境约束较少、室内室外均能使用等优势，因此，在三维重建中使用非常普遍[21-23]。冯林等[24]基于无人机采集的倾斜影像与部署的地面控制点，采用SfM算法构建了高精度侵蚀沟表面模型，对其建模精度与数字高程模型、正射影像等进行分析，并与传统正射航图建模结果进行比较，证明该方法在侵蚀沟的高精度建模与检测方面具有显著优势。陈钰等[25]用简单便携的非测量智能手机获取砾岩露头的图像序列，基于SfM算法进行三维重建、生成点云数据，结果表明，该方法可以快速准确地对地质体进行重建，为野外地质调查及施工等工作提供一种新的思路。在植物表型测量方面，基于SfM算法的植株三维重建具有重建精度较高，能够动态、无损、大批量重建研究对象等优点，为植株性状提取和表型参数测量提供了一种新的方法[26-28]。张慧春等[29]构建了一套机器视觉系统，用于采集拟南芥植株的二维图像，利用SfM算法生成三维点云数据，测量拟南芥的叶片宽度、长度、主茎长度、叶片面积和叶片夹角等表型参数，与传统人工接触式测量值相比，该系统测量的拟南芥叶片宽度、长度、主茎长度、叶片面积、叶片间夹角的平均相对误差分别为9.830%、10.100%、1.070%、4.090%和4.370%，效率高、速度快，可满足拟南芥形态表型分析的需要。HUI等[30]基于RGB相机获取茄子、青椒和黄瓜3种植物不同生长时期的图像序列，利用SfM算法重建三维点云模型，基于得到的点云数据，提取3种植物的叶片数、叶长、叶宽和叶面积等表型参数，并与利用激光扫描仪得到的三维点云重建结果进行对比，结果表明，基于图像序列的三维模型与对应的三维激光扫描模型差别较小，其三维重建精度较高。HE等[31]搭建了一套多视角的立体成像系统，用于采集草莓果实360°图像数据，并基于SfM算法重建三维模型，利用自定义的软件进行三维点云分割，提取草莓果实的高度、长度、宽度、体积、花萼大小、颜色和瘦果数7个性状参数，与人工实测数据对比表明，提取的7个性状参数相关性良好、准确性较高。
本文基于SfM算法重建玉米植株三维点云模型,提出一种自动分割、提取玉米单株、茎秆、叶片的方法。使用前期研制的履带式拍照小车在户外自动获取多视角玉米植株二维图像，采用SfM算法进行三维重建，生成点云数据，提出玉米三维点云自动分割和性状提取方法，以实现户外玉米株高、单株最小包围盒体积、茎粗、叶片数、叶周长、叶面积、单片叶最小包围盒体积、叶投影面积、叶投影宽度、叶投影长度、叶夹角等参数的动态、准确、无损测量。
1 材料与方法
1.1 实验材料与数据获取
本实验共选择3个玉米品种(晋单73、大丰30、郑单958)，每个品种有30株(分3个重复种植，每个重复单独种植一行，一行每个品种各10株)，共90株玉米植株样本，玉米行间距为90 cm,每行内两株玉米的初始距离为30 cm。为了数据后期的处理方便，实验样本采用盆栽种植，避免植株之间的遮挡。用于种植玉米的土壤先在太阳下晒干，然后将晒干的土壤初步碾碎，过6 mm网筛去除土壤中的石头和杂草，保证土壤的均一性，最后将尿素、磷酸二氢钾、氯化钾肥料按照4∶3∶1的配比搅拌均匀。在玉米生长过程中提供正常的水肥管理，待玉米播种60 d后(玉米抽穗期)开始进行玉米株高、茎粗、叶面积等参数的动态无损测量及人工对比验证。
由于天气等外在原因，最终用于实验的玉米植株样本共80株。人工测量玉米植株株高、茎粗、叶面积和地上部生物量的具体步骤为：① 3名工作人员分别用塑料长尺以盆沿为零刻度线基准测量同一株玉米植株高度，取3名工作人员读数的平均值作为该玉米植株的最终人工实测株高。② 3名工作人员用游标卡尺(0～200 mm开式，精度0.01 mm)测量同一株玉米植株的相同部位(距根部10 cm处)并读取数值，取3名人员读数的平均值作为该玉米植株的人工实测茎粗。③用叶面积仪(YMJ-C/YMJC型)测量玉米叶片面积。④用电子天平(LT2002T型)称量玉米植株的地上部生物量。
为全自动地在户外获取玉米植株二维图像，本实验使用前期研制的能按规定路线自动行走的履带式拍照小车采集不同视角下的玉米植株图像。小车主要由计算机系统、图像采集系统、可编程逻辑控制器(Programmable logic controller,PLC)系统、磁导航系统、供电系统、履带式底盘结构等组成，如图1所示。为达到比较好的重建效果，要求相邻两幅图像的重叠率在75%以上[13]，这要求相机的移动速度(即履带式拍照小车的速度)不能太快，通过测试选取小车速度为0.1 m/s。相邻两幅图像的拍摄间隔为5～6 cm。图像采集系统由两个可见光相机组成，分别用来采集侧视图像和俯视图像，俯视相机镜头与水平方向夹角为90°，与地面的垂直距离为150 cm,侧视相机镜头与水平方向夹角为30°，与地面的垂直距离为110 cm，相机型号均为Canon EOS 77D,镜头为EF-S18-200mm f/3.5-5.6 IS,分辨率为6 000像素×4 000像素，图像格式为jpg；计算机内存为3.5 GB，主频3.5 GHz。相机获取图像经过USB线传送到计算机，最终获取的俯视图和侧视图各700幅。
图1 履带式小车户外工作现场图Fig.1 Work scene of autonomous crawler in outdoor1.磁条 2.背景布 3.俯视相机 4.太阳能板 5.侧视相机 6.小车框架 7.控制箱 8.磁导航传感器
实验所用软件为Microsoft Visual Studio 软件、Robot Operating System下的跨平台开源C++点云库(Point cloud library,PCL)和三维重建开源软件Visual SFM,性状提取系统流程图如图2所示。履带式小车在户外以0.1 m/s的速度连续工作4 h获取玉米植株二维图像；将获取的图像导入到Visual SFM软件进行三维重建生成点云数据；将点云数据进行下采样、去噪、坐标校正等预处理；根据预处理后的点云数据，使用直通滤波、条件欧氏聚类、直径拟合等算法实现自动去除背景点云、分割单株、提取茎秆和叶片等功能；基于分割出来的单株、茎秆和叶片，利用距离最值遍历、三角面片化、条件欧氏聚类和随机采样一致性等算法自动测量株高、单株最小包围盒体积、茎粗、单片叶最小包围盒体积、叶投影宽度、叶投影长度、叶周长、叶面积、叶投影面积、叶片数和叶夹角等参数；与人工测量的株高、茎粗、叶面积作对比，来评估SfM重建玉米植株的精度和鲁棒性。
图2 性状提取系统流程图Fig.2 Flowchart of traits extraction system
1.2 玉米植株三维模型重建
玉米植株三维模型重建是在Visual SFM软件中完成，重建的基本原理是利用大量重叠度较高的图像通过SfM算法重建物体三维模型[32]。SfM算法的核心主要包括特征点提取、立体匹配、姿态估计和参数优化等[33-35]。具体重建过程如下：
(1)特征点提取与匹配：通过尺度不变局部特征描述子(Scale-invariant feature transform,SIFT)提取图像特征，用K维空间二叉树(kd-tree)模型计算两幅图像特征点之间的欧氏距离来进行特征点的立体匹配。
(2)稀疏点云重建：用集束调整算法(Bundle adjustment,BA)减少因图像增多而累计的误差，并估计每幅图像的相机参数，生成稀疏点云。
(3)稠密点云重建：用多视角立体集群算法(Cluster multi-view stereo,CMVS)对图像集进行聚簇，以减少重建过程的数据量，提高运算速度和重建精度，然后用多视角拼接算法(Patch-based multi-view stereo,PMVS)通过匹配、膨胀、过滤，在局部光度一致性和全局可见性约束下完成稠密点云的重建[36]。
将履带式拍照小车获取的700幅图像导入到Visual SFM软件进行三维重建并生成三维点云数据。图3为玉米植株三维重建过程。
图3 玉米植株三维重建过程Fig.3 Three-dimensional reconstruction process of maize plants
1.3 玉米点云数据预处理
由于玉米三维重建后点云数据很大，夹杂着许多噪声背景点云，且玉米三维点云模型与实际植株在标准三维空间方向、尺度上不一致，因此需要进行下采样、降噪、坐标校正等预处理。
1.3.1点云数据下采样
由于用Visual SFM软件重建出来的原始玉米三维点云数据很稠密，为了节约算法运算时间，采用下采样方法来精简点云数据。下采样是在不改变点云外轮廓的前提下大幅度减少点云的数量，从而提高算法运算效率。下采样的方法采用三维体素栅格法[37-38]，通过将输入的点云数据创建三维体素栅格，然后在每个体素内，用体素中所有点的重心近似代替该体素内的所有其他点，这样所有体素内的点最终都由该体素的重心所代替，从而达到精简点云数据、减少运算时间、提高程序运行速度的目的。如图4b所示，点的数量降至点云原图(图4a)的30%，但单株的外轮廓几乎没有变化，不影响后面性状参数的提取。
图4 玉米点云下采样和去噪效果Fig.4 Effects of downsampling and denoising
1.3.2点云数据降噪
由于受RGB相机精度、实验操作者经验、外界环境等因素的影响，点云数据中不可避免地会出现一些噪声点和局外点。这些噪声点和局外点对性状提取、特征匹配和曲面重建都有不良影响[39]。本研究采用统计方法去除噪声点和局外点[40-41]。对于下采样后点云数据中的每个点，计算其与k近邻的平均距离(k近邻的值设置为50)，通过假设得到的结果服从高斯分布，具有一个平均值(μ)和标准差(σ),如果所有这些点的平均距离大于阈值Tthreshold，将被视为异常值,从数据集中删除，计算公式如下
Tthreshold=σC
(1)
式中C——用户所使用的常量，本研究取1
本研究也根据三维点云的颜色信息来去除噪声点，计算所有点的超绿分量值(ExG)[42]，将ExG小于0.36的点视为噪声点，将其从三维点云数据中去除，计算公式如下
(2)
式中R、G、B——红、绿、蓝颜色分量
ExG——超绿分量值
玉米三维点云去噪效果如图4c所示，其中点的数量降至下采样后单株点云数量的95%，与图4b相比，减少的点不多但都为噪声点(图4b中的红色椭圆框)。
1.3.3点云数据坐标校正
为了精确地提取玉米植株性状，需进行玉米三维点云的坐标校正。为了获取玉米三维点云模型和真实玉米植株的比例，需进行坐标比例校正，以花盆为基准计算出坐标比例，计算公式如下
(3)
式中hreal——花盆真实高度
hreconstructed——花盆重建模型的高度
r——坐标比例
由于SfM重建的三维玉米点云模型的主方向与PCL点云库中可视化下的三维坐标轴方向存在一定的偏差，所以需要使用旋转与平移矩阵进行坐标转换校正，将点云模型的质心点转换到坐标系原点，转换的主要步骤如下：①计算预处理后的三维点云模型质心点的三维坐标。②基于主成分分析法(Principal component analysis,PCA)[43]将点云模型进行质心化，然后计算点云模型的协方差，求解协方差矩阵的特征值和特征向量，特征向量即为玉米三维点云的主方向。③根据三维点云的质心坐标和主方向创建一个相对于空间坐标系原点的平移和旋转矩阵。④通过创建的旋转和平移矩阵，将原始点云转换为以原点为中心的新点云，空间坐标轴为新点云主方向的位置。
点云坐标转换公式如下
TA=MTTO
(4)
式中MT——平移旋转矩阵
TO——原始点云数据
TA——转换之后的点云数据
在可视化PCL点云数据库下，新位置的点云模型主方向与玉米实际株高的生长方向不一致，为了便于后续点云数据的处理，需要再次使用平移旋转矩阵使玉米植株生长方向与坐标轴z的正方向保持一致。点云坐标转换校正流程如图5所示。图中，红轴和蓝轴分别表示x轴和z轴,黄色十字交点为坐标原点。
图5 点云坐标转换校正流程图Fig.5 Flowchart of point cloud coordinate conversion correction
1.4 玉米植株三维点云分割
玉米植株三维点云分割主要分为分割玉米植株与花盆及以下点云、单株分割、茎秆提取以及叶片提取4部分。处理过程如图6所示，具体点云处理步骤如下：①预处理后的原始玉米点云如图6a所示。通过直通滤波[44]，以花盆真实高度为基准将目标点云分割为玉米植株和花盆及以下点云两部分，分割效果如图6b所示。图6中不同颜色代表不同类，绿色代表玉米植株点云，蓝色则代表花盆及以下点云。②将分割出来的玉米植株采用条件欧氏聚类算法[45-46]分割出单株玉米植株，分割效果如图6c所示。图6c中不同颜色代表不同的玉米植株点云。③提取出单株玉米点云，如图6d所示，基于随机采样一致性算法(Random sample consensus,RANSAC)[47]进行圆柱拟合分割出玉米茎秆，分割效果如图6e所示。同时提取出非茎秆点云，如图6f所示，将提取的非茎秆点云再次采用统计法去除噪声点和局外点，去除后的效果如图6g所示，被去除的点如图6f中的红色椭圆框所示，把去除噪声点和局外点的非茎秆点云再次采用条件欧氏聚类算法分割出单片叶片，分割效果如图6h所示，图6h中不同的颜色代表不同的类。
图6 玉米三维点云分割流程Fig.6 Flowchart of three-dimensional point cloud segmentation
图7 玉米性状参数提取流程Fig.7 Traits parameters extraction process of maizes
1.5 玉米植株性状提取
基于分割得到的玉米点云，计算玉米株高、单株最小包围盒体积、茎粗、叶片数、叶周长、叶面积、单片叶最小包围盒体积、叶投影面积、叶投影宽度、叶投影长度、叶夹角等参数，性状参数提取流程如图7所示，具体参数计算方法如下：
(1)玉米株高
本研究中玉米单株的最高点与盆沿平面在z轴方向的差值默认为玉米的株高。提取出单株玉米植株点云(图7a)，遍历所有点，由于校正后单株玉米的生长方向与z轴的正方向一致，且玉米的底部平面(盆沿平面)为oxy平面，所以要找到单株玉米植株z坐标的最大值和盆沿平面的z坐标值(当前为0)，其差值的绝对值即为单株玉米植株的高度(图7b)。
(2)单株玉米植株最小包围盒体积
提取出单株玉米植株，校正前的单株玉米植株如图7c所示，黄色线条所组成的长方体为包围盒，校正后的单株玉米植株已转换至其主方向上(图7d所示，黄色线条所组成的长方体即为最小包围盒)。要找到校正后单株玉米点云最大的x、y、z坐标值和最小的x、y、z坐标值，得到8个顶点，这8个顶点连接起来所形成的长方体体积即为单株最小包围盒体积。
(3)玉米茎粗
本研究中盆沿平面往上10 cm的平面内茎秆点云坐标最大值与最小值的差值默认为整株玉米的茎粗。该茎粗为本研究的定义，人工测量茎粗时也是按照这个标准进行测量。提取出茎秆点云(图7m)，首先用算法自动截取z=10 cm处平面内的点云数据(图7n)，遍历该平面内的所有点，找到该平面内坐标的最大值和最小值，它们的差值绝对值即为茎粗(图7o)。
(4)叶片数
提取去除噪声点和局外点的非茎秆点云(图7e)，进行条件欧氏聚类，分割出单片叶片(图7f，不同的颜色代表不同的类)，其中聚类出的不同类的数量即为叶片数。
(5)叶周长
分割出三维玉米叶片(图7g)，采用贪婪投影三角算法[48]进行三角面片化，面片化后的三维叶片模型由若干个空间小三角面片组成，其中每个小三角面片的3条边至少被使用一次，而最外围的边即叶片的边界仅被某个小三角面片使用一次，不会被其它三角面片所共用，所以找出所有仅被使用一次的边，计算该边的长度并求和即为叶周长(图7h，红色边界线的长度为叶周长)。
(6)玉米叶面积
分割出玉米叶片，首先采用贪婪投影三角算法进行三角面片化，面片化后的叶片模型由若干个空间三角面片组成，然后通过海伦公式计算单个空间三角面片的面积，最后通过面积求和公式计算单个叶片的面积，计算公式如下
(5)
(6)
式中pi——面片化三角形周长的一半
ai、bi、ci——面片化三角形各边边长
n——总面片数i——面片索引序号
Si——单个空间三角形面片的面积
S3D——单片叶片总面积
(7)玉米单片叶最小包围盒体积
分割出玉米单片叶片后，将叶片点云转换至主方向位置，转换前的玉米单片叶片如图7i所示，黄色线条所组成的长方体为包围盒，转换后的玉米单片叶片如图7j所示，黄色线条所组成的长方体即为最小包围盒。找到转换后单片叶点云最大的x、y、z坐标值和最小的x、y、z坐标值，得到8个顶点，这8个顶点连接起来所形成的长方体体积即为单片叶最小包围盒体积。
(8)玉米叶投影面积
将分割出的玉米叶片向oxy平面投影(图7p)，生成对应的投影叶片点云(图7q)，首先采用贪婪投影三角算法对投影叶片点云进行三角面片化，面片化后的投影叶片模型由若干个平面三角面片组成，然后通过海伦公式计算单个平面三角面片的面积，最后通过面积求和公式计算单片叶片的叶投影面积，计算公式如下
(7)
(8)
式中pj——面片化三角形周长的一半
aj、bj、cj——面片化三角形各边长
m——总面片数j——面片索引序号
Sj——单个平面三角形面片投影面积
S2D——单片叶片总投影面积
(9)叶投影宽度
将分割出的玉米叶片向oxy平面投影，生成对应的投影叶片点云，把生成的投影叶片点云转换至主方向位置(图7r)，计算宽度方向坐标最大值和最小值，其差值的绝对值默认为本研究中的叶投影宽度(图7r)。
(10)叶投影长度
将分割出的玉米叶片向oxy平面投影，生成对应的投影叶片点云，把生成的投影叶片点云转换至主方向位置，计算长度方向坐标最大值和最小值，其差值的绝对值默认为本研究中的叶投影长度(图7r)。
(11)叶夹角
分割出单片玉米叶片后，由于整片玉米叶片弯曲幅度比较大且每片玉米叶片弯曲的幅度不一致，将其用随机采样一致性算法进行平面拟合，求取所得到的平面与茎秆的一个方向向量的夹角，该方法所求取的叶夹角可信度比较低，因此本研究统一截取叶片点云的一部分进行叶夹角求取。截取的规则为：以靠近茎秆的一端为起始端，在向叶片延伸的方向上截取某个阈值内的叶片点云(图7k所示，绿色为完整的单片叶片点云，红色为截取的点云，蓝色为茎秆所在位置)，将截取的叶片点云进行平面拟合，拟合出的平面与茎秆的一个方向向量夹角默认为叶夹角(如图7l所示，红色为截取的点云，黄色平面为截取的点云所拟合出来的平面，蓝色为茎秆所在位置)。
2 结果与分析
2.1 计算评估
将自动测量的玉米植株株高、茎粗、叶面积与人工测量值进行对比来评价本文方法的精度。测量结果精度用平均绝对百分比误差(Mean absolute percentage error,MAPE)、均方根误差(Root mean square error,RMSE)和决定系数R2(Determination coefficient)来评估。
2.2 性状参数测量结果与分析
本实验对抽穗期的80株玉米样本用SfM进行三维重建，自动分割玉米的单株、茎秆和叶片点云，并进行测量分析。将自动测量的株高、茎粗、叶面积和叶投影面积与人工测量值进行对比，结果如图8所示。
图8 玉米植株性状的人工测量值与算法测量值的对比图Fig.8 Comparison plots of algorithm and manual measurements of maize plant traits
由图8a可看出，R2=0.970,MAPE为3.163%,RMSE为3.557 cm,算法测量株高的准确率为96.837%。由图8b可知，R2=0.842，MAPE为4.760%,RMSE为1.540 mm,算法测量的茎粗准确性为95.240%。由图8c可看出，R2=0.901, MAPE为19.102%, RMSE为48.163 cm2,算法测量的叶面积准确率为80.898%。由图8d可看出，R2=0.728,MAPE 为69.321%, RMSE 为95.493 cm2。相比图8c，投影叶面积测量值与叶面积测量值相差较大，说明叶面积测量方法比单个方向(投影至oxy平面)投影叶面积测量方法更接近实际玉米叶面积。图8结果表明，本文方法具有较高的准确性，算法测量值与人工测量值具有较好的一致性。
单株最小包围盒体积、叶片数、单片叶最小包围盒体积、叶周长、叶投影宽度、叶投影长度和叶夹角7个性状参数的测量结果如图9所示。在单株性状参数上，晋单73品种玉米的单株最小包围盒体积和叶片数总体上高于大丰30和郑单958品种，而在单片叶性状上，3个品种的5个性状间(单片叶最小包围盒体积、叶周长、叶投影宽度、叶投影长度和叶夹角)没有非常显著的差异。
2.3 单株性状参数显著性分析
将80株玉米植株分为低地上部生物量和高地上部生物量两类，通过SPSS软件将单株玉米植株的株高、最小包围盒体积、茎粗和叶片数4个性状做单因素方差分析，显著性差异P值分别为0.000 3、0.000 4、0.317 0和0.241 5。当P≤0.01时，说明因素对试验指标的影响非常显著；当0.01
3 结束语
针对传统的玉米植株性状测量方法存在主观性强、劳动强度大、有损伤等问题，以80株抽穗期玉米植株为研究对象，研究了基于运动恢复结构的玉米植株三维重建方法和性状提取方法。结果表明，与人工测量值相比，株高、茎粗、叶面积的平均绝对百分比误差(MAPE)分别为3.163%、4.760%、19.102%，均方根误差(RMSE)分别为3.557 cm、1.540 mm、48.163 cm2,决定系数(R2)分别为0.970、0.842、0.901。同时，实验验证了株高和单株最小包围盒体积可以显著区分低地上部生物量玉米植株和高地上部生物量玉米植株。
图9 玉米叶片性状参数测量结果Fig.9 Measurement results of trait parameters in maize leaves
图10 单株性状参数显著性差异分析Fig.10 Analysis of significant differences in trait parameters of single plant
