基于无人机多光谱影像的番茄冠层SPAD预测研究

田军仓 杨振峰 冯克鹏 丁新军

(1.宁夏大学土木与水利工程学院， 银川 750021； 2.宁夏节水灌溉与水资源调控工程技术研究中心， 银川 750021)

0 引言

番茄是我国北方主要栽培蔬菜之一，其市场需求大、经济效益高[1-2]。叶绿素是植被营养胁迫、叶片氮素含量、植物病虫害、农作物估产的重要反映因子[3-4]，精确预测番茄冠层不同垂直位置叶绿素含量可为番茄生产提供指导。近年来，已有学者应用无人机对田块尺度的精准农业进行了相关研究[5]。无人机[6]遥感监测具有灵活性高、成本低、操作简单、作业周期短、时效性强等特点[7]，在农业领域，其应用从传统的作物种植面积、区域总产、品种单产向作物健康、作物品质、土壤湿度、农业污染控制等更深层方向发展[8-9]，无人机搭载各类光谱传感器技术成为精准农业的重要技术支撑。研究表明，光谱反射及透射特性与叶片及冠层叶绿素含量具有明显的相关性[10-13]，光谱特性成为预测农作物叶绿素含量的重要因素。

植被净光合速率、叶绿素荧光参数及叶面积指数因叶片位置不同而呈现规律性变化[14]。叶片特性受各种条件的影响，在冠层垂直方向上表现出差异性(如光合作用)[15]。随着生长期的变化，营养物质在冠层垂直方向上转移、聚集，以维持代谢平衡[16]。TAWANDA等[17]对多种植物冠层不同垂直位置(茎下部、中部、上部)的叶片性状及光谱特性进行研究，结果表明，冠层垂直位置对单株性状的叶片光谱特性有显著影响。谢传奇等[18]利用高光谱成像仪对灰霉病胁迫下番茄叶片叶绿素进行检测。GARA等[19]通过采集不同时期的树冠叶片反射特性，研究了垂直树冠剖面(冠层上部、下部)生物化学和形态在PROSPECT模型反演叶绿素含量时影响模型性能的因素。丁永军等[20]利用光谱辐射仪采集玻璃温室环境下番茄叶片光谱信息，提取敏感波段，建立叶绿素预测模型。SCHLEMMER等[21]基于遥感数据利用近红外绿边和红边光谱信息，分析了玉米叶片和冠层氮素与叶绿素之间的关系。陈鹏等[22]通过融合光谱与纹理信息筛选特征变量，对生育期马铃薯冠层叶绿素含量进行预测建模。

在农作物叶绿素遥感反演方法的研究中，基于物理模型的反演方法通过建立冠层光谱反射率与农作物结构参数及生化参数之间的辐射传输方程[23-25]来预测叶绿素含量，但因输入参数组合的不同，得到的冠层反射率存在差异，此外，建模复杂程度高，其应用受到限制[26]。以敏感波段及植被指数与叶绿素含量建立回归关系的经验统计模型[27-28]具有操作简单、可行性强的特点，因此，被广泛应用于农作物叶绿素定量遥感监测中。

冠层不同位置的叶片特性存在差异，冠层整体叶绿素预测研究忽略了冠层不同垂直位置的叶绿素含量不同，使预测结果不能精确反映冠层垂直方向上叶绿素的变化情况，而在实际生产中，冠层不同位置的叶绿素精准预测是指导田间管理的重要依据。因此，本文基于无人机多光谱技术[29]，获取番茄种植基地遥感影像，建立番茄冠层上部(TU)、中部(TM)、下部(TL)位置及整个冠层(TC)叶绿素含量预测模型。结合以往叶绿素建模方法研究[30-32]，综合考虑建模可行性及预测效果，选取偏最小二乘法、支持向量机及BP神经网络3种方法进行预测建模。

1 材料与方法

1.1 研究区概况

研究区位于宁夏回族自治区贺兰县常信乡五渠村(北纬38°38′11″,东经106°16′49″)，属中温带干旱气候带，干旱少雨，光照条件充足，蒸发强烈。研究区分为Ⅰ、Ⅱ区，番茄均采用沟垄相间的种植方式，实验品种为2个水平：Ⅰ区为英冠218；Ⅱ区为英冠128。番茄种植行距为80 cm，株距为40 cm，株高在1.6～2.5 m之间，种植时间为2019年5月，土壤类型均为壤土，灌水量、施肥类型及施肥量等田间管理方式均无特定实验处理，Ⅰ区种植行为东西走向，Ⅱ区种植行为南北走向(图1)。

图1 研究概况Fig.1 Study area location

1.2 无人机影像数据获取

无人机影像数据获取时间分别为2019年7月12日(坐果期)、2019年8月6日(结果初期)和2019年9月15日(结果晚期)，选择天气晴朗无风、太阳光照强烈的12:00—14:00采集图像。采用eBee无人机智能系统(SenseFly公司,瑞士)，搭载Sequoia五通道多光谱传感器，参数如表1所示。飞行前用灰度板辐射定标，定标反射率为0.18，垂直地面飞行，设置飞行高度为120 m,航向重叠率80%，旁向重叠率70%，质量控制点数量为6个，获取正射影像图。

表1 Sequoia多光谱传感器参数Tab.1 Sequoia multispectral camera parameters

获取的单幅影像数据通过Pix4D软件的农业多光谱模板进行拼接；将带有POS数据的Sequoia多光谱图像剔除后加载到软件中，采用UTM投影方式，输出坐标系为CGCS2000；导入预先布设的控制点坐标信息进行手动刺点；通过点云加密及纹理特征匹配，最终生成完整影像。

1.3 番茄叶绿素测定

研究表明，叶绿素含量与SPAD(Soil and plant analysis development)值具有良好的相关性，SPAD测量是叶绿素无损测量的有效方法[33]。在无人机影像获取的同时，采用SPAD-502Plus型手持式叶绿素测量仪进行同步测量,测量样本点的株高H，最小值为160 cm，最大值为246 cm，平均值为198.1 cm，标准差为0.052 cm。

依番茄植株特性，取上部(2H/3～H)、中部(H/3～2H/3)、下部(0～H/3)各5片新旧程度存在差异、均匀分布的叶片，使用SPAD-502Plus型手持式叶绿素测量仪，每个叶片测量3组SPAD值，每个部位共获得5组数据，求平均值作为测量数据，对每个样本点的TU、TM、TL测量数据求取平均值作为TC的SPAD测量值。用RTK工具记录样本点位置，样本点均匀布置于研究区内。共设70个实测样本点，每个样本点实测番茄TU、TM、TL、TC 4组数据，获取开花坐果期、结果初期、结果晚期3个生育期共840组数据。整个实验期内番茄SPAD值箱线图如图2所示。

图2 番茄SPAD值箱线图Fig.2 SPAD value box map of tomato

1.4 多光谱辐射定标

eBee无人机在飞行准备时，搭载Sequoia多光谱传感器进行室外辐射定标。辐射定标时，Sequoia多光谱传感器与灰度板垂直距离为1.5～1.8 m，并确定无阴影遮挡灰度板，所获取的单通道反射率图像如图3a所示。在Pix4D软件合成多光谱影像时，取距灰度板外边缘1/3图像部分的数字量化值(DN),番茄多光谱DN值转换为反射率的公式为

(1)

式中R1——目标地物反射率

DN1——目标地物数字量化值

DN2——灰度板数字量化值

R2——灰度板反射率，取0.18

图3 单通道反射率图像Fig.3 Single channel reflectivity images

1.5 反射率提取

文献[34-35]表明，绿光、红光、近红外及红边波段均与叶绿素值具有明显的相关性。在RTK测量样点位置之前，设置RTK测量仪坐标系为CGCS2000，标高为2 m，布设RTK矩形测点，在矩形测点内实测番茄SPAD值。将RTK记录样点坐标输出并在影像图中显示，如图4所示。根据显示位置，在ENVI 5.3支持下，通过RTK工具记录样点位置，绘制感兴趣区域(Region of interest, RoI)，RoI尺寸为20像素×20像素，利用New Region of Interest工具的计算统计功能提取各波段RoI反射率数据。

1.6 植被指数选取

植被指数(Vegetation index，VI)从数学角度出发能更多地反映植被特征，从而提高叶绿素反演精度[23]。结合以往研究成果，本研究选取的植被指数如表2所示。

1.7 多光谱植被指数变化敏感度分析

光谱敏感性分析(Sensitivity analysis，SA)是植被理化参量进行定量反演的基础[42]。为进一步研究所选取多光谱植被指数对番茄SPAD值变化的敏感度，假设多光谱植被指数为模型指标参数，实测番茄SPAD值为不确定因素，光谱参数敏感度指数计算公式为

图4 RoI区域的建立Fig.4 Creating of RoI regions

表2 多光谱植被指数Tab.2 Vegetation index

(2)

式中S——光谱参数敏感度指数

ΔSP——多光谱植被指数变化率

ΔSPAD——番茄实测SPAD值变化率

根据式(2)分别计算各植被指数在TU、TM、TL、TC位置上的敏感度指数S。S越接近1，说明光谱参量对番茄SPAD值的线性敏感程度越好；S越远离1，说明光谱参量对番茄SPAD值的非线性敏感程度越好。

1.8 反演建模

偏最小二乘法(Partial least squares，PLS)是WORD于1983年提出的主要用于研究因变量对于自变量回归的多元统计方法[43]。

支持向量机(Support vector machine，SVR)是一种运用统计学原理的机器算法[44]。本文采用径向基(RBF)函数，通过交叉验证最终确定惩罚因子C为1.66，参数δ为0.5。

BP神经网络基于误差逆向传播算法，可实现正向传播和误差反向传播过程[45]。文中采用输入、隐含、输出3层结构，最终确定输入层-隐含层-输出层结构为9-16-1，学习速率为0.01，迭代次数5 000次。

1.9 模型验证

为避免人为因素对选取模型和验证样本的影响，本文在Matlab平台支持下，应用Divide函数分别将TU、TM、TL、TC 4种冠层位置实测样本按照建模集∶验证集为3∶1的比例进行多次重复分割，直至分割后的建模集和验证集样本满足时空分布均匀。最终确定建模集为158个样本数据，验证集为51个样本数据，固定建模集和验证集样本数据，进行分析建模，采用决定系数R2(Coefficient of determination)、均方根误差(Root mean square error，RMSE)、平均相对误差(Mean relative error，MRE)评价模型反演精度。

图5 不同位置处番茄光谱特征参数与SPAD值相关系数矩阵Fig.5 Matrix diagram of correlation coefficient between spectral characteristic parameters and SPAD value

2 结果与分析

2.1 SPAD值与多光谱植被指数相关性分析

对表2所示的多光谱植被指数与番茄坐果期、结果初期、结果晚期在TU、TM、TL、TC 4种冠层位置上的SPAD值分别进行相关性分析,结果如图5所示。

由图5a可知，TU位置上，REG、NIR、NDVI、OSAVI、CIred-edge、TVI与SPAD值的相关系数在0.60～0.68之间(P<0.01)，而RVI、CVI、MCARI与SPAD值相关系数在0.51～0.55之间(P<0.01)。由图5b可知，TM位置上，OSAVI与SPAD值相关系数最高，为0.68，CIred-edge最低，为0.53(P<0.01)。由图5c可知，TL位置上，REG、NIR、NDVI、TVI与SPAD值相关系数在0.51～0.53之间，OSAVI最高，为0.58(P<0.01)，MCARI、CIred-edge与SPAD值相关系数较低，分别为0.40和0.41(P<0.01)；由图5d可知，TC位置上，除CIred-edge外，其他植被指数与SPAD值相关系数在0.62～0.72之间(P<0.01)，REG、OSAVI较高，分别为0.72和0.71。综合比较，所选植被指数在TC位置上整体相关性较好，其次为TU和TM位置上，在 TL位置上相关程度最差；OSAVI、REG、NIR、TVI与SPAD值相关系数在4种冠层位置下最稳定；CIred-edge在TU位置上与SPAD值相关性程度优于TM、TL和TC位置，CVI、MCARI、RVI在TC位置与SPAD相关性程度较好。

2.2 多光谱植被指数敏感度分析

对所选取建模的9个多光谱植被指数进行SPAD值敏感度分析，如图6所示。由图可知，S变化范围在0.44～1.69之间。OSAVI 在TU、TM、 TL、TC位置上的S变化幅度最小，说明OSAVI 对SPAD值变化线性敏感程度较好，NDVI次之；MCARI在TU、TM、 TL、TC位置上的S变化幅度最大，说明MCARI 对SPAD值变化非线性敏感程度较高，CIred-edge次之。TL位置上，各植被指数的S变化范围均高于TU、TM、TC位置。综合来看，各植被指数在各冠层位置上线性敏感度由优到差的顺序依次为OSAVI、NDVI、REG、NIR、TVI、RV、CVI、CIred-edge、MCARI。结合相关性分析，OSAVI、REG、NIR、TVI与SPAD值的相关系数在4种冠层位置最稳定。因此，TU、TM、TC位置上，本研究选取OSAVI、NDVI、REG、NIR 4种植被指数作为 PLS建模的自变量，TL位置上，各植被指数的S变化幅度较大，通过减少模型输入变量以提高线性回归程度，选取OSAVI和NDVI两种植被指数作为PLS建模的自变量。考虑到BP、SVR建模时为非线性映射，因此以全变量(9种植被指数)作为自变量进行建模。

图6 多光谱植被指数敏感度指数变化范围Fig.6 Sensitivity range of spectral parameters

2.3 番茄的SPAD值反演建模

基于植被指数的番茄冠层不同垂直位置的SPAD值预测建模及验证结果如表3所示。由表3可知，在TU位置上，SVR模型预测结果最优，R2、RMSE、MRE分别为0.68、2.84和0.39。在TM位置上，SVR模型的R2比PLS和BP模型分别高0.11和0.04；RMSE比PLS模型降低了0.26，比BP模型高0.18；MRE介于两者中间。TL位置上，SVR模型R2最高，BP模型RMSE最低。TC位置上，R2由大到小依次为SVR、PLS、BP模型，PLS模型RMSE最小，为2.22，BP模型MRE最小，为0.05。同一预测模型纵向比较，TU和TC位置上R2高于TM和TL位置。在TC、TU、TM、TL位置上PLS和BP模型R2依次减小。SVR模型在TU位置上的R2最高，在TM位置上MRE最小，为0.04。PLS模型，MRE值在TC位置上比TU、TM、TL位置分别低0.30、0.09和0.18，BP模型的MRE从小到大依次为TM、TU、TC、TL。

图7为不同预测模型的验证结果。由图可知，TU、TM、TL、TC位置上，SVR模型验证集R2均高于PLS、BP模型。在TU、TL位置上，BP模型RMSE值均高于PLS和SVR模型，在TM和TC位置上，PLS模型RMSE值最高。同一模型的验证集在不同位置比较结果为，TU位置上验证效果最优，TL位置上验证效果最差。不同模型之间比较，基于全变量的SVR模型在建模集和验证集的拟合结果最优，模型较稳定。

2.4 基于最优预测模型的番茄SPAD值可视化制图

根据表2中各植被指数计算公式，得到每个像素点的植被指数特征。根据2.3节可知，基于全变量构建的SVR模型预测能力最优。因此，将每个像素点的植被指数特征作为模型参量代入构建的SVR预测模型中，从而得到每个像素点的SPAD值，将栅格分类显示，得到番茄SPAD值可视化制图，如图8所示，从左到右依次为开花坐果期、结果初期、结果晚期。由图可知，番茄冠层各部位叶片SPAD值在开花坐果期、结果初期、结果晚期总体呈现衰减趋势。总体来看，与实测番茄SPAD值箱线图(图2)规律一致。在开花坐果期，番茄植株生长茂盛，营养物质向上聚集，因此，上层叶片SPAD值是番茄纵向生长监测的关键；在结果期，果实生产量大，且集中在中层部位，因此，中层SPAD值监测是及时防治番茄结果期病虫害的关键。通过番茄冠层各部位SPAD值可视化填图，更清晰直观地说明番茄冠层不同部位叶片SPAD值的差异，对更加及时准确地掌握番茄植株营养状况分布情况以及给灌水、施肥和病虫害监测等具有重要意义。

表3 不同估算模型的SPAD值建模与验证结果Tab.3 Modeling and verification results of different estimation models for SPAD value

2.5 讨论

(1)本文研究结果显示，无人机多光谱植被指数OSAVI、REG、NIR、NDVI与番茄冠层不同垂直位置上的SPAD值均具有良好的相关性及线性敏感度，并且随着冠层位置降低，其与对应的SPAD值相关性程度及线性敏感度也依次减小，原因可能是番茄株高较高，行距、株距较小，随着冠层高度降低，叶片透光能力减弱，使得下层叶片对光谱信息的贡献较小。后续研究可以考虑用倾斜摄影等方法实现光谱信息在垂直位置上的分层，进而增强垂直地面获取的光谱信息对冠层不同垂直位置叶片生理特性预测的精确性。

(2)SVR模型将低维的非线性关系在高维空间中线性映射，因而具有较好的鲁棒性[46]。在番茄冠层不同位置上用SPAD值预测建模，SVR模型较BP神经网络和PLS模型的R2高，主要原因是SVR模型通过交叉验证得到的惩罚因子C和参数δ使得模型拟合情况最优；BP模型受到隐含层不确定性和样本数量(建模集158个样本，预测集51个样本)的影响，在一定程度上影响了BP模型学习及预测能力；对于多元回归，即使添加一个不显著的自变量，模型R2都会增大，产生过拟合现象[31]，本文通过敏感度分析筛选少数几个PLS模型建模变量，建模变量的数量成为PLS模型的R2较低的主要原因。

3 结论

(1)开花坐果期，番茄冠层上层叶片的SPAD值高于中层和下层叶片，结果初期和结果晚期，番茄中层叶片的SPAD值高于上层和下层叶片。

(2)所选取的光谱植被指数与番茄SPAD值的相关性在TU和TC位置上优于TM和TL位置，OSAVI、REG、NIR、TVI与SPAD值的相关性在不同冠层位置上均最高。

(3)光谱参数线性敏感度分析表明，各光谱植被指数在TU、TM、TL、TC冠层位置上线性敏感度由优到差依次为OSAVI、NDVI、REG、NIR、TVI、RVI、CVI、CIred-edge、MCARI。

(4)同一模型不同冠层位置上的SPAD值预测的R2由大到小依次为TU、TM、 TL，表明利用冠层光谱信息预测冠层上、中、下层SPAD值的精确性依次降低。

(4)采用筛选变量的PLS模型、全变量的SVR模型及BP模型对番茄SPAD值进行预测建模，结果表明：SVR模型结果最优，在TU、TM、TL、TC位置上，R2分别为0.68、0.64、0.55、0.65，RMSE分别为2.84、2.75、3.51、2.26。