李杨先，张慧春，杨旸

(南京林业大学机械电子工程学院，南京 210037)

植物表型是指能够反映植物结构和组成，或能反映植物生长发育过程和结果，由基因型与环境互作所产生的部分或全部可辨识植物物理、生理和生化特征及性状[1]。传统表型测量方法具有样本量小、效率低、误差大、适应性差等缺点，已成为制约植物表型发展的重要因素[2]。

面对植物表型测量技术发展滞后的现状，基于生物、传感器、机器视觉技术以及高性能计算机的表型分析平台能够为植物表型的精确化和高效化测量铺平道路。目前，针对农作物的表型监测技术已得到了较好的发展。与农作物相比，林木较高大、多枝叶密、生长周期长、根系发达和株型较大等特点决定了对其进行表型采集和分析的难度较大。但从林木的经济价值及生态价值考虑，有必要通过研究林木的表型，培育出有利于提高经济价值和提升人类生存环境的树种，从而完善表型信息采集系统在植物领域的应用。倪超等[3]设计了一套基于多目立体视觉的非接触式马尾松苗木形态学参数提取系统，利用该系统提取了苗木的根系体积参数，为苗木质量的评价提供了更加精确的形态指标。束义平等[4]采用车载二维激光扫描仪获取树木的点云数据，并基于点云数据检测分割树冠、识别树干，实现了树冠体积的在线测量；同时，李秋洁等[5]实现了基于移动二维激光扫描的单木三维绿量测定，真实反映了树冠内部体积及空隙。南玉龙等[6]搭建了一套植物冠层超声回波信号检测系统，实现了植物冠层密度超声量化。张慧春等[7]设计了一套植物表型测量系统，利用该系统可以获取植物的三维形态特征并建立植物时序生长可视模型。林木的株高、基径和分枝数等表型参数能够直接反映林木的生长发育情况，而叶片面积和分枝角对光合作用的强弱具有直接影响，但是目前鲜见有关林木这方面形态表型参数的研究。

簸箕柳(SalixsuchowensisCheng)作为模式树种杨树的姊妹种，已完成全基因组测序，且个体相对较小、幼龄期相对较短[8-9]，易于开展大规模田间试验，非常适合作为林木表型研究的对象。本研究以簸箕柳为研究对象，设计了一种基于图像处理技术的植物形态表型参数获取系统，能够对形态表型参数进行快速、精确、非破坏性测量，与传统测量方法相比优势突出，为簸箕柳等林木的生长发育状况研究提供了更精确的形态性状指标。

1 材料与方法

1.1 试验材料

选取簸箕柳为试验植株，剪取15段直径(4±0.2)mm、高(12±0.3)cm的簸箕柳枝条(亲本从江苏省新沂市采集)进行扦插式种植，基质选用南京林业大学林区有机土壤。过筛后装入规格相同的圆柱形黑色陶瓷花盆(口径10 cm、高11.5 cm)中。簸箕柳枝条扦插后，手动浇水至土壤饱和，并置于RXZ型人工气候箱(宁波江南制造厂)中培养。人工气候箱的环境参数设置为相对湿度80%、温度28 ℃、光照强度278 μmol/(m2·s)、昼夜时间比为16∶8。在此环境条件下，簸箕柳存活率高，生长速度较快，有利于试验的进行。

1.2 图像采集平台设计与构建

为获取植株二维图像序列，构建了一套植物形态表型图像采集平台(图1)，该平台主要由密闭采集环境和光源、旋转平台和竖直移动平台、图像采集模块三部分组成。整个暗箱环境框架选用铝型材构建，高角度打光的LED光源[10]和覆盖在框架上的黑色植绒布能够尽量减小拍照的噪音干扰。相机旋转平台和竖直移动平台、图像采集模块的三维示意图见图2a，旋转平台实现相机围绕簸箕柳进行圆周运动，从而获取簸箕柳的二维图像序列；竖直移动平台通过步进电机控制丝杠导轨改变相机在竖直方向上的高度(相机位置见图2b，根据植株的高度不同取H=45～70 cm、h=20～40 cm)，在每个相机位置调节相机角度(相机轴线与丝杠导轨轴线的夹角为15°～75°)，保证簸箕柳和标定板完全处于相机视野范围内，以获取2组不同高度下较为完整的图像信息；图像采集模块由高分辨率GO-5000C-PGE型RGB相机和LM8HC型镜头组成。簸箕柳形态表型参数获取平台需要满足能够自动从2个高度位置获取簸箕柳形态信息且具有较高图像质量的要求，因此，本平台使用LabVIEW 2016软件设计了一个具有相机启停、对焦显示窗口、图像自动存储、数量调整、拍摄时间间隔等功能的相机控制程序，图像采集控制界面程序如图3所示。本试验控制相机旋转一周的时间为195 s，每个相机高度拍摄30张簸箕柳二维图像，共60张。因此，控制相机每隔 6.5 s拍摄一张簸箕柳图像。

1. 棋盘格标定板；2. 暗箱框架；3. 竖直移动导轨；4. RGB相机；5. 暗箱采集背景；6. LED光源；7. 轨道板；8. 步进电机；9. 旋转圆台；10. 图像采集控制系统；11. 电机控制系统；12. 电机驱动系统。图1 簸箕柳形态表型图像采集平台Fig. 1 Morphological phenotype image acquisition platform of Salix suchowensis Cheng

1. 相机镜头；2. 相机采集卡；3. 角度调节装置；4. 丝杠导轨滑块；5. 丝杠导轨步进电机；6. 滑轮；7. 丝杠导轨；8. 定滑轮；9. 旋转台；10. 旋转台步进电机。图2 相机移动平台三维示意图Fig. 2 Three-dimensional diagram of camera moving platform

1. 图像路径存储；2. 图像命名；3. 相机选择；4. 拍摄间隔设置；5. 拍摄张数设置；6. 时间显示；7. 取消；8. 显示窗口。图3 使用LabVIEW设计的簸箕柳图像采集界面Fig. 3 Image acquisition interface of Salix suchowensis Cheng designed by LabVIEW

平台工作流程如下：首先基于图像采集平台进行标定板图像的采集，然后采用“张正友标定法”[11]基于MATLAB软件对相机进行标定，得到相机的内外参数及畸变参数。相机标定完成后，将簸箕柳摆放在平台内相机拍摄圆周圆心处，通过中控计算机的相机控制系统设置拍摄图像数量、拍摄间隔时间、存储路径等参数。启动图像采集平台，由图像采集平台搭载相机围绕簸箕柳旋转，相机控制系统控制相机拍摄，通过平台获取2组不同高度的图像；基于相机标定得到的标定参数对平台采集的簸箕柳二维图像序列进行畸变纠正。畸变纠正后，基于HSV(hue, saturation, value)模型进行阈值分割，目的是将植物和标定板从图像中提取出来。使用SFM(structure from motion)算法对采集平台获取的二维图像序列重建出三维点云，在此基础上使用PMVS(patch-based multi view system)算法重建出稠密三维点云。最后利用PCL(point cloud library)点云库、MATLAB软件等对生成的簸箕柳点云信息进行操作和计算，分别提取出簸箕柳株高、基径、叶面积、分枝数和分枝角参数。本平台可以稳定地获取簸箕柳二维图像，且图像质量清晰，符合平台的预期功能。平台获取的一幅细节图见图4。

图4 平台获取的图像Fig. 4 Image acquired by the platform

1.3 植株点云获取方法

1.3.1 图像预处理

图像预处理是为了将感兴趣的部分与背景等其他干扰像素分割，以此提高图像的匹配效率和后期的测量精度。由图4可知，所需保留的部分为簸箕柳和棋盘格标定板，花盆、泥土、背景幕布等信息为干扰因素。本研究使用HSV颜色模型并结合OpenCV库对图像进行分割[12]。棋盘格标定板的作用为求解点云坐标系中距离与世界坐标系中真实距离的转换。

使用HSV模型分割首先要将图像转换为HSV模型再进行分割操作，HSV颜色模型中绿色的阈值区间为[35，43，46]～[77，255，255]，白色的阈值区间为[0，0，221]～[188，30，255]。由于光照等因素的影响，簸箕柳和标定板不完全是标准色的阈值区间。试验结果表明，簸箕柳选取阈值区间[30，120，130]～[60，190，255]，标定板选取阈值区间[0，0，100]～[90，70，255]得到的分割效果较好。基于HSV颜色模型的分割效果见图5。由图5可知，基于HSV的分割效果干扰信息少，簸箕柳信息保留完整。

图5 阈值分割图像Fig. 5 Threshold segmentation of the image

1.3.2 基于图像序列的植株三维点云获取与处理

图6 簸箕柳稀疏点云重建效果图Fig. 6 Rebuilding effect map of sparse point cloud of Salix suchowensis Cheng

植株三维点云使用基于多视觉立体运动中恢复结构(structure from motion with multi-view stereo，SFM-MVS)的开源软件VisualSFM生成。VisualSFM封装了基于SIFT图像特征匹配、相机参数确定及稀疏点云重建等算法[13]，通过导入2组不同视角的植株二维图像实现植株稀疏点云的重建。其中，每组图像数量分别为30张，共60张。2019年3月19日采集的编号为1的簸箕柳稀疏点云重建效果图见图6，其中，虚线局部放大图表示每一张二维图像的拍摄位置和角度，即拍摄时相机所在位置。

利用SFM重建得到的是稀疏三维点云，为获得更好的重建效果，需对目标重建出稠密三维点云。使用基于片面的三维多视角立体视觉算法(patch-based multi view system，PMVS)生成的植株稠密点云效果见图7。

从图7中可以看出，获得的簸箕柳三维稠密点云包含植株相应的纹理和颜色信息，表明此方法在簸箕柳形态构建上具有较好的效果，能够较真实反映植株形态。

图7 生成的稠密点云效果图Fig. 7 Generated dense point cloud renderings

将生成的三维稠密点云进行滤波预处理以去除噪声点和离群点，获取较为平滑的稠密点云；之后利用标定板进行比例缩放，标定板中每个方格尺寸为50 mm×50 mm，总体尺寸为150 mm×150 mm。在世界坐标系中，原点O坐标为(0,0,0)，标定板平面中心白方格右上角q1点坐标为(0,50,0)，中心白方格左下角q2点坐标为(50,0,0)，中心白方格右下角q3点坐标为(50,50,0)。在生成的点云坐标系中，原点Q在点云坐标系中的坐标为(x,y,z)，选取点p1(x1,y1,z1)、p2(x2,y2,z2)和p3(x3,y3,z3)。通过式(1)计算出坐标位置的缩放关系K，利用K即可进行点云坐标系中距离与世界坐标系中真实距离的转换。

(1)

2 簸箕柳表型三维形态参数提取

为研究簸箕柳的生长发育和光合作用等情况，需要选择合理的形态参数用于分析，本研究选取了株高、基径、叶面积、分枝数和分枝角等表型形态参数进行测量。利用手工测量方式选取的基径、叶面积、分枝角测量点和测量部位如图8所示，其中，去除顶部第一个展开的叶片，从上至下选取3个叶片进行叶面积测量。

图8 形态参数手工测量部位Fig. 8 Manual measurement location of morphological parameters

2.1 株高与基径参数提取

植物株高的生长速度能够反映植物长势情况，杨柳科属植物的株高一般是指植株与土壤交接点到植株顶部叶片自然伸展至最高处的垂直高度。传统株高测量方式是使用直尺量取生长点到植株顶部竖直方向的距离，该方法属于接触式测量，容易对植物造成损伤[14]。本研究中簸箕柳株高参数的提取方法如下：首先利用主成分分析法(principal component analysis，PCA)大致确定簸箕柳植株的主径方向(图9中箭头方向)，然后基于此方向做簸箕柳植株的最小包围盒，包围盒高度即为簸箕柳株高。PCA的实现方法如下：

1)求取中心点。假设输入点集为P={pi|i=1,2,…,n}，点云中点的数量为n，则点云中心点(pm)为：

(2)

2)求特征协方差矩阵。通过式(2)求得的pm求取协方差3×3矩阵Cp：

(3)

3)求协方差矩阵的特征值和特征向量。由于Cp是对称实矩阵，可得到3个非负特征值λ0、λ1、λ2，从而计算出相应的特征向量e0、e1和e2：

Cpei=λiei，i∈{0,1,2}

(4)

其中，e0代表点云中最密集方向(植株主径方向)。

图9 植株高度最小包围盒Fig. 9 Minimum bounding box of plant height

通过人工获取的株高参数记录了15株簸箕柳快速生长期内株高的变化，采集时间为2019年2月20日至2019年3月16日，每隔6 d记录一组数据，记录5次共75组数据。由簸箕柳形态表型参数获取平台提取的簸箕柳株高参数，分别选取3组数据的平均值作为一次测量结果，测量日期和间隔都与手工测量相同。

平台提取基径的方法为：首先截取簸箕柳的生长枝条，通过确定簸箕柳生长点处基径边缘坐标(x1,y1,z1)与水平方向上另一边缘处坐标(x2,y2,z2)，计算两点之间的欧式距离。

通过手工测量的方法记录15株簸箕柳快速生长期内基径的变化。由于扦插后第一周萌芽的枝条较短，基径较小，手工测量和平台都不易提取基径参数，因此，选择簸箕柳扦插后的第二周开始测量，采集时间为2019年3月1日—25日，每隔6 d记录一组数据，记录5次共75组数据。

2.2 叶片面积提取

叶片是植物进行光合作用合成有机物的重要器官，叶片面积大小直接决定了光合作用的强弱，对农作物产量具有重要影响[15]。传统叶片面积测量方法有方格法、描形称质量法和仪器测定法等[16]，本试验中手工测量数据使用YMJ-D型叶面积仪(浙江托普云农科技股份有限公司)测量获得。由于SFM-MVS算法生成点云的规则特征和植株的形态差异，获取的点云数据通常分布不均且含有一定误差[17]。本研究采用最小二乘法对叶片的原始点云数据进行处理，以获取较好的曲面拟合效果。最小二乘法是一种数学优化算法，即利用最小化误差的平方和寻找数据的最佳函数匹配[18]。叶片拟合前后曲面重建的对比效果见图10。

图10 叶片拟合前后效果对比Fig. 10 Comparison of the effects before and after leaf fitting

计算叶片面积须将叶片点云剖分为三维三角形网格，所有三角剖分后三角形网格的面积之和即为叶片面积。由于分割出的各个叶片坐标系存在差异性，因此需要调整叶片点云，使之处于统一测量坐标系。定义的叶片标准测量坐标系以叶片的点云重心为坐标原点，X、Y、Z轴方向分别与叶片的长度、宽度、厚度方向一致，X、Y、Z轴的正方向符合右手法则。统一测量坐标系的目标是使叶片点云在空间中具有相同的朝向，利用主成分分析法调整叶片点云坐标系至本研究定义的标准测量坐标系，使叶片点云具有相同的朝向。

叶片点云经过坐标转换和最小二乘法拟合后，利用Delaunay三角剖分算法(一种特殊的三角剖分算法，符合空圆特征和最大化最小角特征2个重要准则)生成三角形网格，三角剖分后任一三角形的外接圆范围内不会有其他点存在[19]。Delaunay三角剖分算法主要通过以下2个步骤完成：首先，将叶片点云投影到XY平面上，由于对叶片点云进行了坐标转换，因此点云投影到XY平面不会造成点的重合，对投影过的二维点进行Delaunay三角剖分；然后，根据投影前x、y对应的z值将三角剖分后的三角形网格返回到三维坐标系中。上述2个步骤主要通过MATLAB软件中的Delaunay和trisurf函数(用于创建三角剖分曲面图)完成。

簸箕柳叶面积参数的测量时间选择为簸箕柳扦插后的第2周，此时簸箕柳叶片生长较好，也达到了叶面积测量仪所能测量的最小值。采集时间为2019年3月1日—25日，每隔6天记录一组数据，记录5次共75组数据，其中，每组包含3片叶子的测量值。同时，使用叶面积测量仪记录了15株簸箕柳快速生长期内叶面积的变化，每次测量去除顶部第1个展开的叶片，从上至下依次选取3个叶片。

2.3 分枝数与分枝角提取

簸箕柳分枝数是单株簸箕柳从主干上生长出的枝条数。簸箕柳的分枝方式为单轴分枝，主干单一且较为明显，簸箕柳分枝数等于自身叶片数量，每一个分枝上仅有1个叶片，因此提取簸箕柳分枝数即对簸箕柳叶片数进行提取。将背景、花盆与主干部分进行分割，只留下叶片部分，并利用分割后的叶片进行三维重建，对于主干分割不理想的情况则后期在点云中进行人为分割处理。在三维空间中，每一个叶片点云之间无交集，可将每个叶片点云看做一个点云团，且每个点云之间的距离小于叶片点云团之间的距离。因此，本研究采用基于欧式距离聚类分割的思想对点云团进行计数，并设定阈值，将噪点信息排除。

由于簸箕柳前期分枝不明显，因此，选择扦插后第2周开始对15株簸箕柳进行分枝数的测量，采集时间为2019年3月1日—19日，每隔6天记录一组数据，记录4次共60组数据。

分枝角对树形结构的形成具有决定性作用，分枝角的大小对簸箕柳枝条间的生长和光照效果有一定影响[20]。传统的分枝角测量方式为使用量角器测量主干与分枝之间的角度θ。本研究提取分枝角的方法为提取两个向量之间的夹角，在点云中选取主干上两点设置为向量a，分枝上选取两点设置为向量b。分枝上选取的第1个点位于分枝与主干的交接处，第2个点与第1个点的连线应尽可能沿分枝初始伸长的切线方向，在保证上述选点的条件下，两点间的距离应尽可能远。利用式(5)求得分枝角θ：

(5)

簸箕柳顶部分枝角基本呈锐角，底部分枝角基本呈钝角，因此，分别提取顶部和底部两部分的分枝角进行测量。顶部分枝角的选择为去除顶部3个分枝后的第1个分枝，底部分枝角的选择为去除底部3个分枝后的第1个分枝。选择扦插后第3周开始对15株簸箕柳进行分枝数的测量，采集时间为2019年3月7日—19日，每隔6天记录一组数据，记录3次共45组数据。

2.4 试验数据对比

簸箕柳形态表型参数提取值与人工测量值的比较如图11所示。簸箕柳形态表型参数提取效果通过使用形态表型参数提取的结果误差来评价，估计误差由决定系数(R2)、误差均方根(RMSE)和平均绝对百分比误差(MAPE)衡量。本系统提取的簸箕柳各表型参数与人工测量值的吻合程度较好，能够较真实地反应簸箕柳的形态属性。株高、基径、叶片面积、分枝数以及分枝角的人工测量值与本系统测量值的RMSE分别为 8.91 mm、0.54 mm、63.65 mm2、1.36和 5.81°；MAPE分别为7.69%，13.94%，9.99%，5.32%和9.30%。

由图11a可知，与人工测量值相比，该系统提取的株高整体偏高，结合系统获取株高的方法进行分析可知，系统获取的株高值普遍高于人工测量值的原因在于所生成的簸箕柳三维点云数据受到噪点的影响，在植株顶端存在噪点干扰，导致生成的包围盒产生了一定误差；簸箕柳生长点处一般会萌发新叶片，生长点处会出现叶片遮挡情况，使生长点处点云重建效果不理想，影响包围盒的大小。由图11b可知，系统基径测量值普遍大于人工测量值，且误差随着基径的增大逐渐减小。结合试验分析，由平台获取的基径测量值平均绝对百分比误差较大，主要是因为簸箕柳基径尺寸较小，一般为2～7 mm。系统基径测量误差来源主要是基径选取点与人工选取点存在一定的偏差，此外，基径选取点的偏离也是导致误差产生的主要原因，而基径边缘噪点的干扰也会导致测量值普遍大于人工测量值。

图11 形态表型参数测量准确性评价Fig. 11 Accuracy evaluation of measurements of morphological phenotype parameters

由图11c可知，与人工测量值相比，该系统提取的叶片面积值整体偏低。簸箕柳的叶片表面有褶皱结构，利用SFM和CMVS算法生成的植株三维点云不能很好地体现这些细节，且通过最小二乘法进行拟合会平滑叶片点云，导致叶片面积的计算值偏小。由图11d可知，系统测量分枝数与人工测量分枝数大多重合或较为接近，对于分枝数较少的阶段有部分数据偏低，分枝数较多的阶段有部分数据偏高。经过对产生误差的数据分析可知，产生误差的主要原因在于图像的分割质量。虽然该方法对图像边缘要求较低，但图像内部不能存在缺失部分，如果存在缺失部分就会导致生成的点云之间存在较大间隙，叶片数会增大，从而导致分枝数增加。而在分枝数较少时，主要也是因为叶片存在重叠导致分割效果不理想，使叶片无法区分，叶片数减少，导致分枝数也减少。由图11e可知，系统测量得到的分枝角基本都比人工测量得到的分枝角大。结合实际试验对此结果进行分析可知，测量误差主要来源于主干向量与分枝向量的选取，对于点云重建效果不理想的分枝部位进行向量选取时会产生相对较大的误差；同时，对于分枝较细的部位也会产生较大误差，分枝较细时，线段的选取不够准确。

3 结 论

本研究提出了一种基于图像处理技术的植物形态表型参数获取方法，以模式树种杨树的姊妹种簸箕柳作为研究对象，设计了一套基于RGB相机传感器的形态表型参数获取平台，结合二维图像恢复物体几何外形的方法无损提取了簸箕柳的株高、基径、叶片面积、分枝数和分枝角表型参数。结果表明：

1)通过对相机高度的调节获取2组不同高度的图像序列，补充一部分被遮挡部分的信息，能够在一定程度上减少叶片等器官相互遮挡带来的影响，提高生成的植株点云数据的精度。

2)基于HSV阈值分割的方法能够较好地将目标植物与干扰背景进行分割，簸箕柳自身信息保留较为完整，有利于提高特征匹配速度和形态表型参数的提取。

3)基于二维图像恢复物体几何外形的方法可以较精确获得植株点云数据，提取的株高、基径、叶面积、分枝数和分枝角参数除基径外平均绝对百分比误差均在10%以内。基径由于数值较小，因此误差相对偏大，但能够满足植株实际测量和生长观测的需求， 同时，经过预处理后三维点云的生成速率比使用原图像提高20%左右。

本研究中基于图像处理技术的植物形态表型参数获取方法适用于林木表型信息测量，为林木生长监测、育种改良等提供了一种新的表型测量方法，今后将探讨该方法用于杨树等其他林木性状参数测量的可行性，为植物表型性状提取与测量提供新的技术手段。本研究设计的植物形态表型参数获取平台在自动化程度和参数提取效率上仍有提高的空间，在植株器官分割上未能实现自动分割。因此，提高平台的自动化程度和采集数据的多源性能够实现大规模的植物形态表型参数提取，为植物学提供有力的研究平台。