无人机多光谱影像的马铃薯地上生物量估算

刘 杨，孙 乾，黄 珏，冯海宽，王娇娇，杨贵军

1. 农业部农业遥感机理与定量遥感重点实验室，北京农业信息技术研究中心，北京 100097 2. 山东科技大学测绘科学与工程学院，山东 青岛 266590 3.南京农业大学国家信息农业工程技术中心，江苏 南京 210095 4.国家农业信息化工程技术研究中心，北京 100097

引 言

地上生物量(above-ground biomass，AGB)是作物活力和净初级生产力的基础，能够表征作物光合作用强度和评估作物营养状况以及长势情况[1]。故快速无损地获取AGB信息可监测作物长势变化，这对于指导农业生产管理和精准施肥至关重要。传统上，AGB测量是采用破坏性取样方法，需要人工收割作物、 称重和记录，这使得大面积、 长时间的监测作物变为困难[2]。

随着遥感平台的多样性和影像空间分辨率的提高，遥感技术被当今认为是监测和估算大面积作物AGB最有效的技术手段[3]。其中低空无人机遥感平台具有机动灵活，能够提供更高的时间、 空间和光谱分辨率的遥感数据等优势，成为了精准农业中定量观测使用最为频繁的技术工具[4]。当前，无人机平台搭载的传感器主要有3种(数码、 多光谱和高光谱)。数码相机虽然价格低廉，但是波段数目较少，难以解译作物冠层光谱中包含的复杂信息[5]。高光谱传感器在高精度表征光谱响应能力方面虽表现突出，但价格昂贵，数据后续处理较为复杂[6]。而多光谱传感器由于经济上更适宜，且包含监测农学参数信息重要的红边波段，因此在农业定量遥感中得到广泛关注[7]。国内外学者基于无人机多光谱遥感技术开展了作物理化参数监测的大量研究，如Han等[4]使用6种植被指数结合不同的回归技术，评估了玉米不同生育期AGB的估算效果； 肖武等[8]利用22种植被指数结合BP神经网络有效地估算采煤沉陷区上方玉米作物的AGB； Qi等[7]使用8种植被指数通过人工神经网络估算花生不同生育期的叶面积指数； 孙诗睿等[9]运用10种植被指数结合随机森林估算冬小麦叶面积指数； Zheng[10]和Brinkhoff等[11]分别通过19种和6种植被指数估算冬小麦和水稻氮素含量，结果发现RDVI和NDRE植被指数估算效果最佳； 陈鹏等[12]和奚雪等[13]分别使用植被指数估算马铃薯和冬小麦叶绿素含量。

以上基于无人机多光谱遥感技术的研究，主要通过单个或多个植被指数估算作物不同生育期的理化参数，没有引入新的模型变量去参与建模。然而，随着生育期的推进，植被指数可能会达到饱和，当理化参数值较大时，会存在低估的现象。为了更全面地估算作物多生育的理化参数，一些学者通过多类别变量或者回归技术提高了估算模型精度，较好地实现作物长势监测[14-19]。但还未见有报道关于从无人机多光谱影像中提取高频信息(high frequency information，HFI)和作物株高(Hdsm)并结合植被指数估算马铃薯不同生育期的AGB。因此，本研究基于光谱信息以及融合结构信息采用偏最小二乘回归(partial least squares regression，PLSR)和岭回归(ridge regression，RR)方法构建马铃薯不同生育期的AGB估算模型，探究融合新的模型因子对AGB估算结果的影响，从而为提高马铃薯AGB监测效果提供新的方法。

1 实验部分

1.1 研究位置与试验设计

马铃薯试验在北京市昌平区小汤山镇国家精准农业研究示范基地进行，该区域位于北纬40°10′32″—40°10′36″，东经116°26′34″—116°26′41″，气候类型为暖温带半湿润大陆性季风气候。试验设计为小区完全随机试验，共设计密度试验、 氮素试验、 钾肥试验3个试验区。试验品种均为早熟的中薯5和中薯3，每个品种进行相同的控制试验，每种试验重复3次，每种重复进行不同程度的密度、 氮素和钾肥处理，其中密度、 氮素和钾肥具体处理详情见文献[14]。小区总计48个，每个小区面积为32.5 m2。为了更精确地获取试验田的数字表面模型，在试验小区周围均匀布控11个地面控制点(ground control point，GCP)，并用差分GPS测定其三维空间位置。

1.2 无人机多光谱数据获取及处理

分别于天空晴朗、 无风无云的2019年4月20日(裸土期)、 2019年5月28(块茎形成期)、 2019年6月10日(块茎增长期)和2019年6月20日(淀粉积累期)进行无人机多光谱遥感数据采集。各次起飞地点固定和飞行航线基本保持一致，操作时间为12:00，飞行高度为30 m，航向重叠度为80%，旁向重叠度为85%。获取影像数据前，首先采集传感器自带的光谱反射率校正板数据，用于多光谱影像像元亮度值的标定。采用八旋翼电动无人机搭载Parrot Sequoia农业遥感专用的4通道多光谱相机，其由光照传感器和多光谱传感器组成，能获取1个1 600万像素的RGB影像和4个120万像素的单波段影像。多光谱相机内置的波段参数信息见表1所示。

表1 多光谱传感器的波段参数Table 1 Band parameters of multispectral sensor

无人机多光谱影像获取后需要预处理，处理前，首先对无人机采集的多光谱影像进行筛选，去除姿态角异常，成像存在问题的影像。将选取的多光谱影像和11个地面控制点导入到以运动结构为核心算法的Agisoft PhotoScan Professional软件中，通过GCP的空间位置对多光谱影像作地形精纠正，各时期的校正误差均小于2 cm。基于密集点云生成试验区域的数字正射影像(digital orthophoto map，DOM)和数字表面模型(digital surfacemap，DSM)。采用ENVI 5.1软件通过校正板获取各波段的反射率影像，并以TIF格式储存。使用ArcGIS软件基于划分的小区绘制出最大面积矢量并对其编号，分别计算统计出所有感兴趣区域在各波段的平均光谱反射率，将其作为各小区马铃薯冠层在该波段下的光谱反射率。

1.3 地面数据采集及处理

地面数据采集时，获取了各小区马铃薯实测AGB和株高数据。马铃薯AGB和株高的测定方法见文献[14]。

1.4 植被指数的选取

根据以往的研究成果，各生育期选取9种在AGB估算方面表现较好的多波段组合植被指数(multi-band combined vegetation indices，MVIs)和本研究提取的4个单波段植被指数(single-band vegetation indices，SVIs)用于构建马铃薯AGB估算模型，具体的模型参数见表2所示。

1.5 高频信息获取

二维离散小波变换的影像分解技术通过包含低通和高通的滤波器组对无人机多光谱影像做两次滤波，分别得到水平方向(HL)、 垂直方向(LH)和对角线方向(LL)3种高频信息和1种低频信息，其中高频信息反映的是影像的大致概貌和轮廓，与影像的真实信息接近。由于红边是植被反射率从红光波段到近红外波段快速升高变化的区域，是区分植被和地物的最显著标志，经常被用来估算作物的理化参数[2]。因此，提取各生育期红边波段的3种高频信息(REGHL，REGLH和REGLL)用来参与模型构建。

1.6 分析方法

马铃薯株高的提取，利用数据管理、 分析和可视化显示为平台的ArcGIS软件，将马铃薯块茎形成期、 块茎增长期和淀粉积累期的DSM分别与裸土期的DSM进行差值运算，可以得到相应的作物高度模型，最后利用ROI工具提取出各生育期试验小区的马铃薯平均高度Hdsm。

表2 模型参数Table 2 Model parameters

构建AGB估算模型时，若输入变量存在严重的共线性问题，则会降低模型的稳定性和准确性。因此，采用2种方法(PLSR和RR)估算马铃薯各生育期AGB。PLSR利用了多元线性回归、 典型相关分析和主成分分析结合为一体，可以提供一种多对多的线性回归建模方法，能够消除自变量之间的相关性，用较少数据来估测因变量[14]。RR是一种专用于共线性数据分析的有偏估计回归方法，实质上是一种改良的最小二乘估计法，通过放弃最小二乘法的无偏性，以损失部分信息、 降低精度为代价获得回归系数更为符合实际、 更可靠的回归方法[2]。

1.7 精度评价

各生育期提取48组数据构成样本数据集，选取32个用来构建模型，剩余16个用来验证模型效果。为了评估不同模型的估算效果和稳定性，选用决定系数(coefficient of determination，R2)和标准均方根误差(normalized root mean square error, NRMSE)作为精度评价指标。

2 结果与讨论

2.1 马铃薯株高的提取

3个生育期共提取144个马铃薯株高，为了评估提取的精度，将地面采集的实测株高数据和Hdsm作比较分析，R2为0.87，NRMSE为14.34%，说明提取的株高Hdsm较为可靠。

2.2 植被指数、 高频信息和Hdsm与AGB的相关性

将表2的13种植被指数、 红边波段的3种高频信息和基于无人机多光谱遥感技术提取的Hdsm分别与马铃薯各生育期的AGB作相关性分析，得到结果见表3所示。由表3可知，各生育期提取的模型参数均与AGB达到0.01显著水平，相关性均从块茎形成期到淀粉积累期先升高后降低。整体上，各模型参数与AGB的相关性由高到低依次为多波段组合植被指数、 高频信息、 单波段植被指数和Hdsm。块茎形成期相关性最高的模型参数为RVI，相关性系数为0.750，块茎增长期和淀粉积累期相关性最高的模型参数均为GNDVI，相关性系数分别为0.762和0.759。

表3 植被指数和Hdsm与马铃薯AGB的相关性系数Table 3 The correlation coefficients of vegetationindices and Hdsm with potato AGB

2.3 基于植被指数估算AGB

为了评估原始4个单波段植被指数(SVIs)和9个多波段组合植被指数(MVIs)估算AGB的能力，分别使用PLSR和RR方法构建各生育期AGB估算模型，其结果见表4所示。由表可知，各生育期使用2种方法基于SVIs和MVIs构建的模型效果均从块茎形成期到淀粉积累期先好后变差。使用同种方法以MVIs构建的模型精度更高，稳定性更强，且均在块茎增长期达到最佳估算效果(PLSR： 建模R2=0.65，NRMSE=17.48%； 验证R2=0.68，NRMSE=16.71%。RR： 建模R2=0.62，NRMSE=18.42%； 验证R2=0.65，NRMSE=17.75%)。另外，从3个生育期的估算模型精度和稳定性来看，使用PLSR方法估算AGB的效果要优于RR方法。

2.4 基于植被指数结合HFI或Hdsm估算AGB

为了探究融合新的模型因子(HFI或Hdsm)对AGB估算结果的影响，将各生育期提取的13种植被指数(vegetation indices，VIs)分别结合红边波段的3种高频信息和Hdsm使用PLSR和RR方法建立各生育期的AGB估算模型，其结果见表5所示。从表5可以看出，各生育期基于植被指数结合HFI或Hdsm利用2种方法估算AGB的效果同表4的结果保持一致，也从块茎形成期到淀粉积累期先好后变差。相较于仅以植被指数构建的AGB估算模型(表4)，结合HFI或Hdsm明显提高了估算模型的拟合性和稳定性(表5)，其中植被指数融合红边波段的高频信息效果较优，2种方法均在块茎增长期达到最高的估算效果(PLSR： 建模R2=0.72，NRMSE=15.44%； 验证R2=0.74，NRMSE=15.33%。RR： 建模R2=0.67，NRMSE=16.34%； 验证R2=0.70，NRMSE=15.62%)。通过分析各生育期的表5建模和验证结果，发现使用PLSR方法基于植被指数结合高频信息或Hdsm估算AGB的效果也优于RR方法。

表4 基于植被指数使用PLSR和RR估算马铃薯AGBTable 4 Estimation of potato AGB using partial least squares regression and ridge regression based on vegetation indices

表5 基于植被指数结合HFI或Hdsm使用PLSR和RR估算马铃薯AGBTable 5 Estimation of potato AGB using partial least squares regression and ridgeregression based on vegetation indices combined with HFI or Hdsm

2.5 基于植被指数结合HFI和Hdsm估算AGB

将植被指数结合HFI和Hdsm作为模型输入参数，同样使用PLSR和RR方法构建各生育期AGB估算模型，其结果见表6所示。从表6可以看出，整个生育期基于VIs+Hdsm+HFI使用2种方法估算AGB效果变化趋势(先好后变差)与表4和表5结果相一致，其中基于融合所有特征为变量得到的估算结果最出色(表6)。2种方法也均在块茎增长期达到最佳估算精度(PLSR： 建模R2=0.73，NRMSE=15.22%； 验证R2=0.75，NRMSE=14.62%。RR： 建模R2=0.69，NRMSE=15.56%； 验证R2=0.71，NRMSE=15.47%)。各生育期基于融合所有特征的PLSR-AGB模型的R2分别提高7.69%和2.94%，7.35%和1.38%，4.83%和3.17，RR-AGB模型的R2分别提高10%和6.45%，6.15%和2.99%，10.34%和4.9%(相较于表4)。分析表6建模和验证结果可知，马铃薯3个生长期使用PLSR方法基于融合特征估算AGB效果也同样优于RR方法。

表6 基于植被指数结合HFI和Hdsm使用PLSR和RR估算马铃薯AGBTable 6 Estimation of potato AGB using partial least squares regression and ridgeregression based on vegetation indices combined with HFI and Hdsm

通过马铃薯3个生育期的无人机多光谱影像，结合GCP生成了试验田的DSM，提取了各生育期株高Hdsm，实测株高和Hdsm拟合的R2为0.87，证实了基于DSM提取的Hdsm效果较优，此结论与他人提取作物株高的研究基本相同[5]。但是，基于无人机多光谱影像提取的马铃薯株高相比于实测株高偏小，这是因为多光谱传感器获取的是马铃薯冠层数据，而地面实测的株高最高点位置在生成3D点云时被去除，导致基于DSM提取的Hdsm较小； 另外，获取的多光谱影像中包含一些裸土像元，在提取Hdsm时也参与了运算，这同样会使提取的株高偏低。

将4个单波段的植被指数、 9个多波段组合的植被指数、 红边波段的3种高频信息和Hdsm分别与AGB作相关性分析，结果表明各生育期的模型参数与AGB的相关性均达到0.01显著水平，这说明提取的各类参数都能够反映作物的长势情况[1]。3个生育期的相关性均从块茎形成期到淀粉积累期先升高后降低，这与马铃薯作物自身的生长状况相关。前期主要表现为马铃薯生殖器官发育，茎节和叶片不断地增长，当进入块茎增长期，马铃薯植株地上各部位的鲜重达到峰值，植被覆盖度为整个生育期的最优时期，此阶段提取的冠层光谱反射率能够真实反映作物AGB的变化情况。而生长后期地上的同化物需要向地下块茎转移，基部叶片自下而上逐渐衰老变黄，马铃薯作物长势变差，植被覆盖度也明显降低，此时提取的光谱信息并不是马铃薯植株冠层真实的反射率，使得上述3种光谱参数与AGB的相关性降低。

以SVIs(x1)，MVIs(x2)，VIs结合Hdsm(x3)，VIs结合HFI(x4)，以及VIs结合Hdsm和HFI(x5)为模型输入变量，使用PLSR和RR方法估算马铃薯各生育期的AGB。结果发现，基于5种变量使用同样的方法构建的模型效果变化趋势，均从块茎形成期到淀粉积累期先好后变差，这与模型参数和AGB的相关性变化趋势相一致。各生育期以不同变量使用同样的方法估算AGB的精度由高到低依次为x5>x4>x3>x2>x1，主要因为融合结构信息(HFI和Hdsm)解决了生育期效应引起的植被指数饱和性问题，提高了模型预测能力，这与刘杨等[14]得到的结论一致。相较于单波段植被指数(x1)估算模型，基于多波段组合植被指数(x2)构建的模型精度和稳定性都较优，主要因为通过多波段组合的植被指数能够去除或者降低背景土壤对马铃薯冠层光谱信息的影响，增强了植被指数与AGB的敏感性(表3)，以此提高了AGB估算精度，如Bispo等[15]研究结果也认为多波段组合的植被指数估算AGB效果较优。当采用波段组合的植被指数估算AGB时，大多通过绿、 红和近红外波段的反射率经过波段运算得到宽波段参数[15]，忽略了红光波段与近红外区域的红边参数，红边是植被特有的光谱特征，在整个生育期内与作物参数敏感性较高，因此红边位置对于研究AGB的动态变化非常重要[18]。然而，仅仅通过植被指数估算作物不同生育期的AGB，随着生育期的推进，植被指数会出现饱和的现象，这会造成估算AGB不准确[4]。因此，将提取红边波段的3种高频信息(HFI)和作物株高(Hdsm)一起融入到植被指数中形成新的模型因子来估算各生育期的AGB，结果表明结合作物光谱信息和结构信息构建的模型精度最高，稳定性最强(表6)，这与Li等[1]和Yang等[3]研究作物理化参数结果一致，都表明融入作物结构信息能够解决植被指数造成的低估现象。

为了减弱模型参数之间的自相关性，使用PLSR和RR方法构建各生育期AGB估算模型，探究了这2种方法估算AGB的效果。结果表明每种变量以PLSR方法构建的模型R2较大，NRMSE较小，说明此方法估算效果要优于RR方法，这与Tao等[17]和Yue等[19]研究冬小麦AGB结果一致，都证明PLSR方法估算效果较优。综上，研究作物不同生育期AGB时，模型参数类别和建模方法的选取对构建AGB估算模型的精度影响不同。因此利用无人机多光谱遥感技术估算作物AGB时，模型参数的适用性和建模方法的筛选还需要进一步研究。此外，还需要获取不同年限和不同地点的马铃薯冠层光谱数据对模型的外推性进行验证，以便得到一个更具有代表性的AGB估算模型。

3 结 论

(1)基于无人机多光谱影像提取的Hdsm与马铃薯实测株高之间有较高的拟合性(R2=0.87，NRMSE=14.34%)，表明提取的Hdsm精度可靠。

(2)选取的4个单波段植被指数、 9个多波段组合的植被指数、 红边波段的3种高频信息和提取的Hdsm与AGB相关性均达到0.01显著水平，各种变量与AGB的相关性变化趋势保持一致，均从块茎形成期到淀粉积累期先升高后降低。块茎形成期，相关系数最高的模型变量为RVI(r=0.750)，块茎增长期和淀粉积累期，相关系数最高的模型变量均为GNDVI(r=0.762和0.759)。

(3)各生育期以5种变量使用相同的回归方法估算马铃薯AGB，从块茎形成期到淀粉积累期估算效果先好后变差，其中估算精度由高到低依次为x5>x4>x3>x2>x1。

(4)各生育期使用PLSR以不同变量估算AGB的效果要优于RR方法，其中在块茎增长期使用PLSR结合x5变量估算马铃薯AGB，精度达到最高(建模R2=0.73，NRMSE=15.22%； 验证R2=0.75，NRMSE=14.62%)。