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

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

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

引 言

地上生物量(above ground biomass， AGB)是评估作物营养状况和长势情况的重要理化参数， 与品质和产量有着密切联系， 及时准确地估算AGB有助于提高农田生产管理和作物监测水平[1]。 传统的AGB测量方法需要人工破坏性取样， 虽能达到较高精度， 但过程复杂， 费时费力， 无法满足大面积应用需求。 然而， 遥感技术通过远距离、 非接触式探测目标的电磁波特性， 为无损高通量估算作物AGB提供了有效途径。 无人机高光谱遥感虽传感器价格昂贵， 数据后续处理复杂， 但因具有操作灵活、 图谱合一和光谱分辨率高等优势， 在农情监测和产量预测等方面得到广泛关注[2]。

近年来国内外学者应用高光谱成像技术监测作物AGB， 取得了一定的研究进展。 如Yue等[3]基于冬小麦冠层光谱数据构建多种植被指数， 结合不同的回归技术有效地估算作物多生育AGB。 Tao等[4]利用相关性分析法(correlation analysis method， CAM)筛选出与冬小麦AGB敏感的光亲指数和红边参数， 使用2种回归方法准确地评价作物长势情况。 Kanke等[5]分析了红和红边波段组成的植被指数与水稻AGB的定量关系， 最终确定了较好估算AGB的敏感光谱波段。 刘斌等[6]研究高光谱数据任意两波段构建的NDVI与冬小麦AGB的相关关系， 通过比较模型间的精度， 得出估算AGB最优中心波段和波段宽度。 陶惠林等[7]基于无人机高光谱数据构建的多种植被指数和提取的红边参数， 采用多元线性回归方法估算冬小麦多生育期AGB， 结果表明加入红边参数信息能够较好地改善模型精度。 石雅娇等[8]基于玉米冠层高光谱数据构建敏感植被指数， 结合BP神经网络估算玉米AGB， 建模的决定系数达到0.99， 估测效果极好。 这些研究成果表明， 通过光谱分析技术能够很好地估测冬小麦、 水稻和玉米等作物AGB， 然而马铃薯植株形态结构与上述作物存在显著差异， 且在作物养分吸收、 输送和转移等方面表现也明显不同。 前期地下营养物质往地上供应， 促进植株茎叶生长， 而后期地上生殖器官光合作用积累的有机物向地下块茎转移， 促进果实膨大， 整个生育期AGB值呈现先升高后降低的变化趋势， 故已有的其他作物AGB监测成果无法在马铃薯作物上直接应用。

目前， 马铃薯作物多生育期的长势动态监测， 主要集中于冠层氮素含量、 叶面积指数、 叶绿素和叶片含水量的估测研究。 如Liu等[9]基于马铃薯冠层高光谱影像， 使用偏最小二乘回归(partial leastsquares regression， PLSR)分别结合全谱、 可见光-近红外和短波红外的光谱反射率估算氮素含量， 结果表明利用全谱波段信息构建的模型效果较好。 Luo等[10]从马铃薯冠层高光谱数据中提取不同类型的光谱参数， 分别结合一元回归技术估算不同水分条件的叶面积指数， 结果表明水分饱和条件下使用深度面积比指数估算精度最高。 孙红等[11-12]基于马铃薯冠层高光谱影像， 通过CAM和随机蛙跳算法(random frog algorithm， RFM)分别挑选出与叶绿素和叶片含水量相关的特征波长， 结合PLSR方法有效地实现马铃薯理化参数的反演。

这些研究成果对马铃薯作物监测有一定的参考和应用价值， 但是， 针对马铃薯多生育期AGB的估测研究， 鲜有报道。 因此， 为了探究生长过程中AGB的动态光谱响应机制， 本研究采用无人机成像高光谱技术获取不同生育期的马铃薯冠层影像， 同时为了减少全谱数据的冗余度和增强模型的稳定性， 通过CAM， RFM和高斯过程回归波长分析工具(Gaussian process regression-band analysis tool， GPR-BAT)分别对冠层原始光谱(canopy original spectra， COS)和一阶导数光谱(first derivative spectra， FDS)筛选敏感波长， 结合使用PLSR和GPR方法建立AGB估算模型， 从而确定估算AGB的最佳波长和最优反演模型， 以期为通过无人机高光谱遥感快速无损地监测马铃薯长势情况以及分析AGB含量提供方法依据。

1 实验部分

1.1 试验设计

试验地点位于北京市昌平区小汤山镇国家精准农业研究示范基地， 地处北纬40°10′35″， 东经116°26′40″。 马铃薯种植试验为小区完全随机试验设计， 共设计密度试验(N区)、 氮素试验(S区)、 钾肥试验(K区)3个试验区。 试验品种采用早熟的中薯5(Z5)和中薯3(Z3)， 每个品种进行相同的控制试验， 每种试验重复3次， 每种重复进行不同程度的密度(T1， T2和T3)、 氮素和钾肥处理， 其中密度、 氮素和钾肥具体处理详情见文献[1]。 小区总计48个， 每个小区面积为32.5 m2。 为了更精准地获取试验田的位置， 在试验小区周围均匀布控11个地面控制点(G1—G11)， 并用差分GPS测定其三维空间位置， 详细的试验设计见图1所示。

图1 试验区域和试验设计Fig.1 Experimental area and experimental design

1.2 无人机高光谱影像采集及处理

采用八旋翼电动无人机搭载德国Cubert公司生产的UHD185成像光谱仪(波长范围为450～950 nm， 光谱分辨率4 nm， 共有125个光谱通道， 光谱采集时保证仪器垂直向下)， 分别在2019年5月28日(块茎形成期)、 2019年6月10日(块茎增长期)和2019年6月20日(淀粉积累期)进行高光谱遥感作业。 各时期飞行时间为12:00—13:00， 此时天空晴朗， 无风无云， 飞行高度为50 m， 获得的影像空间分辨率为13 cm。 每次无人机采集数据前， 需要在地面利用黑白板进行高光谱影像辐射校正。

高光谱数据预处理过程主要包括两个部分： (1)影像拼接和几何地形纠正， 首先使用Cuber公司生产的Cuber-Pilot软件融合航带内的灰度影像和高光谱影像， 进而形成新的融合后的高光谱影像； 其次通过PhotoScan软件利用地面控制点的位置信息进行影像地形校正， 各时期的校正误差均小于2 cm； 最后基于高密度点云数据完成影像拼接， 生成马铃薯3个生育期的数字正射影像(DOM)和数字表面模型(DSM)。 (2)提取冠层光谱反射率， 在ArcGIS软件中根据各小区的最大面积矢量数据， 通过IDL语言提取出各小区平均光谱， 将平均光谱反射率作为不同小区冠层光谱反射率。

1.3 地面数据采集及处理

地面数据采集与无人机高光谱遥感作业同步进行， 每个生育期共获取48组实测AGB数据， 具体采集过程见文献[2]。

1.4 特征波段筛选方法

为了降低光谱数据的冗余度， 提高模型的稳定性和预测能力， 使用相关性分析法(CAM)、 随机蛙跳算法(RFM)和高斯过程回归波长分析工具(GPR-BAT)分别对马铃薯3个生育期的冠层原始光谱(COS)和一阶导数光谱(FDS)进行特征波长筛选。 CAM是数据分析过程中最为常用的敏感波长筛选方法， 根据相关系数的大小决定建模的个数和最优模型参数。 RFM是一种类似于可逆跳转马尔科夫链蒙特卡洛的变量选择方法， 通过迭代的方式计算变量在每次迭代过程中被选择的概率， 根据概率值的高低来评价变量的重要性[12]。 GPR-BAT通过光谱数据依次迭代， 去除对模型贡献最小的波长， 采用10折交叉验证法进行模型内部验证， 基于交叉验证均方根误差RMSECV的最小值来确定最优特征波长数。 CAM在Excel2019软件中进行， RFM和GPR-BAT在Matlab R2020b软件中进行。

1.5 建模方法

PLSR是解决模型参数共线性问题最为常用的方法， 在数据拟合过程中通过对光谱反射率矩阵和AGB含量矩阵同时分解， 提取最佳主成分后， 将二者进行关联， 建立线性的回归模型， 从而达到估算AGB的目的[12]。 GPR是一种非参数概率统计模型， 基于贝叶斯定理来学习自变量和因变量之间的关系， 利用均值和协方差函数根据最大似然估计法来训练样本， 可以提供预测及其相关的置信区间， 这能够评估预测结果的可靠性。 与常规机器学习法相比， 参数优化更简单， 更适合训练小样本数据， 最大的优势在于能够通过GPR-BAT工具箱自动识别最佳光谱特征[13]。

1.6 精度评价

选取67%样本数据(32个， 重复一和重复二)作为建模集， 33%样本数据(16个， 重复三)作为验证集， 利用PLSR和GPR方法构建各生育期马铃薯AGB估算模型。 为了评估不同模型的拟合效果和稳定性， 采用决定系数R2、 均方根误差RMSE和标准均方根误差NRMSE作为精度评价指标。R2越高， RMSE和NRMSE越低， 构建的模型精度就越高。

2 结果与讨论

2.1 马铃薯冠层高光谱反射率曲线

通过标准正态变量(standard normal variable， SNV)校正方法对马铃薯各生育期获取的冠层高光谱数据进行修正， 以此减少背景噪声、 地物纹理等多种因素对光谱反射率的影响[11]。 以块茎增长期为例， n13—n15， s17—s20， k02和k05小区经过SNV处理后的冠层原始光谱曲线如图2所示。 由图可知， 马铃薯冠层反射率曲线符合绿色植被特征， 在450和650 nm附近存在吸收谷， 而在550 nm附近存在小的反射峰。 由于叶肉细胞结构的影响， 在680～750 nm范围内， 光谱反射率急剧增加， 形成植物特有的红边特征。 对处理后的光谱反射率进行一阶微分， 以及利用CAM， RFM和GPR-BAT筛选与AGB相关的敏感波长。

图2 马铃薯块茎增长期经SNV处理后的光谱反射率曲线Fig.2 Spectral reflectance curves treated with SNV during potato tubergrow period

2.2 CAM筛选特征波长

将马铃薯3个生育期的COS和FDS分别与AGB作相关性分析， 得到结果如图3所示。 为了避免边界波长震荡效应的影响， 只对458～946 nm波段范围作进一步研究。 以各生育期的相关系数绝对值由大到小依次排列， 基于COS和FDS块茎形成期分别筛选出相关系数大于0.7的特征波长总数为28个和12个， 块茎增长期分别筛选出相关系数大于0.7的特征波长总数为32个和18个， 淀粉积累期基于COS筛选出相关系数大于0.7的特征波长总数为30个， 基于FDS筛选出相关系数大于0.6的特征波长总数为21个。 各生育期使用CAM得到的具体模型参数见表1所示。

表1 各生育期基于FDS和COS使用CAM筛选的敏感波长Table 1 The sensitive wavelengths selected by CAM based on FDS and COS during each growth period

图3 马铃薯各生育期COS和FDS与AGB的相关系数图Fig.3 Correlation coefficient diagrams of COS and FDS with AGB during each growth period of potato

2.3 RFM筛选特征波长

在MatlabR2020b软件中， 设置RFM初始变量个数为2个， 迭代次数10 000次， 运行结果如图4所示。 马铃薯3个生育期基于COS和FDS选择概率值大于0.2的变量为特征波长， 块茎形成期选择的特征波长数分别为12个和23个， 块茎增长期对应为8个和28个， 淀粉积累期对应为15个和33个。 各生育期使用RFM得到的具体模型参数见表2所示。

图4 马铃薯各生育期COS和FDS的变量选择概率Fig.4 Variable selection probabilities of COS and FDS during each growth period of potato

表2 各生育期基于FDS和COS使用RFM筛选的敏感波长Table 3 The sensitive wavelengths were selected by RFM based on FDS and COS during each growth period

2.4 GPR-BAT筛选特征波长

通过ARTMO软件中MLRA模块的GPR-BAT对马铃薯3个生育期的COS和FDS作敏感性分析， 各生育期进行内部交叉验证的平均RMSECV、 标准差SD以及min-max极值范围见图5所示。 根据RMSECV的最小值确定估算AGB的特征波长， 块茎形成期， 基于COS和FDS筛选的敏感波长个数分别为6个和10个， 块茎增长期分别为2个和4个， 淀粉积累期分别为3个和5个。 各生育期使用GPR-BAT得到的具体模型参数见表3所示。

表3 各生育期基于FDS和COS使用GPR-BAT筛选的敏感波长Table 3 The sensitive wavelengths were selected by GPR-BAT based on FDS and COS during each growth period

图5 马铃薯各生育期使用COS(上)和FDS(下)的交叉验证图Fig.5 Cross validation maps of using COS (upper) and FDS (lower) during each growth period of potato

2.5 构建AGB估算模型

将通过CAM， RFM和GPR-BAT基于COS和FDS筛选的敏感波长作为自变量， 马铃薯AGB作为因变量， 分别使用PLSR和GPR方法构建各生育期的AGB估算模型， 各模型的精度评价指标见表4和表5所示。 由表可知， 基于不同冠层光谱数据， 各生育期通过3种方法筛选的敏感波长， 利用PLSR和GPR建立的模型效果表现一致， 均从块茎形成期到淀粉积累期由好变差， 其中各生育期基于FDS得到的模型变量估算AGB精度更高， 模型较为稳定。 基于同种冠层光谱数据， 通过GPR-BAT筛选的敏感波长使用2种方法估算AGB的效果最优， 其次为RFM， 而通过CAM筛选的特征波长估算效果最差。 基于不同方法对COS和FDS筛选的敏感波长， 各生育期以相同变量使用PLSR构建的AGB估算模型效果优于相应地GPR-AGB估算模型。 基于COS通过GPR-BAT筛选的特征波长， 使用PLSR方法估算AGB，建模R2从块茎形成期到淀粉积累期的变化范围为0.62～0.71， RMSE变化范围为211.03～307.31 kg·hm-2， NRMSE变化范围为16.37%～21.18%， 验证结果与建模结果保持一致，R2越大， RMSE和NRMSE值越小； 而使用GPR方法建模R2相对较小， 变化范围为0.61～0.67， RMSE和NRMSE值相对较大， 变化范围分别为228.30～321.82 kg·hm-2和17.71%～22.18%， 验证效果与建模结果一致，R2越小， RMSE和NRMSE值越大。 基于FDS通过GPR-BAT筛选的特征波长， 使用PLSR方法估算AGB效果最佳， 块茎形成期到淀粉积累期建模R2从0.65变化到0.73， RMSE从203.07 kg·hm-2变化到301.95 kg·hm-2， NRMSE从15.84%～20.81%， 验证R2也是先增大后减小， RMSE和NRMSE值都是先减小后增大； 使用GPR方法估算AGB的建模R2相应地从0.62变化到0.69， RMSE从216.95 kg·hm-2变化到306.30 kg·hm-2， NRMSE从16.83%变化到21.11%， 验证R2， RMSE和NRMSE的变化趋势与建模结果相同，R2先增大后减小， RMSE和NRMSE先减小后增大。

表4 基于COS使用PLSR和GPR估算AGB的建模和验证精度Table 4 Modeling and verification accuracies of estimating AGB based on COS by PLSR and GPR

表5 基于FDS使用PLSR和GPR估算AGB的建模和验证精度Table 5 Modeling and verification accuracies of estimating AGB based on FDS by PLSR and GPR

以往研究通过无人机成像高光谱技术估算作物AGB， 大多利用全谱数据构建光谱指数并结合不同的回归技术实现AGB含量监测， 但是这样会增加模型的复杂性和降低运算效率。 为解决这一问题， 本研究通过CAM， RFM和GPR-BAT3种方法对马铃薯各生育期的COS和FDS分别作敏感性分析， 从而筛选特征波长结合PLSR和GPR构建AGB反演模型。 结果显示， 3个生育期构建的模型中， 块茎增长期以不同变量使用同种方法建立的估算模型精度均高于其他生育期， 这是因为马铃薯植株从块茎形成期开始， 由原来的营养生长转变为生殖生长和物质积累， 地上茎叶逐渐发育完善， 到了块茎增长期， 茎叶生长速度、 叶面积指数和茎叶鲜重达到峰值， 此时植被覆盖度最大， 提取冠层光谱反射率时不易受到地面土壤的干扰， 随后由于地上同化的有机物向地下块茎输送， 造成地上茎叶因营养匮乏而枯黄脱落， 此时马铃薯作物长势变差， 间接提取的光谱反射率包括大量的裸土像元， 因此参与建模的变量不能真实反映AGB的实际情况。

各生育期分别用基于COS和FDS通过RFM方法筛选的特征波长来估算AGB的效果要优于相应地CAM方法， 这是因为RFM方法筛选的特征波长间隔大， 跨度广， 包含信息量丰富， 模型变量自相关性较弱， 因此构建的模型精度较高、 稳健性较强， 这与孙红等[11-12]估算马铃薯叶片叶绿素和含水量结论一致， 但本研究利用RFM和CAM方法筛选的敏感波长来估算AGB的精度较低， 主要原因其一可能是估算的理化参数不同， 其二是传感器类型和数据获取场景不同， 本研究通过无人机平台搭载UHD185成像光谱仪在田间获取数据， 而孙红等研究在封闭实验室载物台通过Gaia高光谱成像系统获取数据， 其三是本研究对马铃薯冠层群体植株进行AGB监测， 而孙红等仅以单个叶片为目标实现作物参数估算。 为了验证通过利用GPR-BAT筛选敏感波长来估算AGB的效果， 同样对COS和FDS作敏感性分析， 结果发现， 各生育期基于GPR-BAT筛选波长运行效率较低， 但得到的特征波长个数少， 使用PLSR和GPR建模方法估算AGB效果更优， 这一结论与Fu等[13]研究冬小麦氮素情况相一致， 也表明GPR-BAT筛选的敏感波长与作物理化参数联系更紧密。 基于同种变量， 各生育期利用PLSR方法构建的模型精度更高， 这与Tao等[4]估算冬小麦多生育期AGB结果一致， 这主要与其自身处理光谱数据能力相关， 可以较好地解决变量间的多重共线性问题， 使得能够准确地估算作物理化参数[1-2]。 综上表明， GPR-BAT筛选的敏感波长在符合AGB含量与光谱反射率间变化的规律前提下， 可结合PLSR方法使用较少波长来预测各生育期AGB， 以达到准确监测的目的。 但是， 本研究并未考虑不同水分灌溉下， 基于马铃薯冠层光谱筛选的敏感波长对AGB监测结果的影响， 需要在未来的研究中进行深入探究。

3 结 论

(1)基于COS和FDS使用CAM， RFM和GPR-BAT方法筛选的特征波长个数在块茎形成期分别为28， 12， 6个和12， 23， 10个， 在块茎增长期分别为32， 8， 2个和18， 28， 4个， 在淀粉积累期分别为30， 15， 3个和21， 33， 5个。

(2)在相同条件下， 各生育期基于FDS筛选的敏感波长相比于基于COS筛选的敏感波长更能准确估算AGB。

(3)各生育期通过GPR-BAT筛选的特征波长估算AGB效果最优， 其次为RFM方法， 而CAM方法筛选的特征波长估算效果最差。

(4)各生育期基于FDS通过GPR-BAT筛选敏感波长， 并结合PLSR方法建立的模型精度更高， 块茎形成期建模R2， RMSE和NRMSE分别为0.67， 203.07 kg·hm-2和16.63%， 块茎增长期建模R2， RMSE和NRMSE分别为0.73， 204.19 kg·hm-2和15.84%， 淀粉积累期建模R2， RMSE和NRMSE分别为0.65， 301.95 kg·hm-2和20.81%。