试验区位于江苏省扬州市广陵区沙头镇试验基地(32° 24'N, 119° 26'E), 属亚热带温润气候, 年平均气温14.8 ℃, 光热条件丰富, 土质为沙壤土, 以小麦和水稻轮作为主。 如图1所示, 该试验区按照四个施氮水平NO、 N14、 N18、 N22(N0: 0 kg· ha^-1, N14: 210 kg· ha^-1, N18: 270 kg· ha^-1, N22: 330 kg· ha^-1), 和24个品种处理(如表1)被划分为96个小区, 小区大小为3 m× 3 m。 每个小区都接受了相同的灌溉和田间管理。

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 冬小麦试验区位置 (a): 试验地块位置; (b): 冬小麦试验田可见光无人机图像Fig.1 Location of the research area of winter wheat (a): Location of the research area; (b): UAV RGB images of the winter wheat research fields

表1

Table 1

表1(Table 1)

表1 本研究所涉及24个小麦品种 Table 1 The 24 types of winter wheat involved in this study 小麦品种 最佳生长习惯 小麦品种 最佳生长习惯 扬麦15 春性 宁麦13 春性 扬麦20 春性 宁麦26 春性 扬麦23 春性 宁麦资126 春性 扬麦25 春性 农麦88 春性 扬麦29 春性 国红6 春性 扬辅麦4 春性 烟麦1 春性 镇麦9 春性 郑麦9023 弱春性 镇麦10 春性 洛麦24 弱春性 镇麦18 春性 烟农19 半春性 镇麦13 春性 济麦22 半春性 镇麦12 春性 徐麦33 半春性 镇麦168 春性 淮麦33 半春性

表1 本研究所涉及24个小麦品种 Table 1 The 24 types of winter wheat involved in this study

根据江苏省冬小麦生长发育期观测资料分析, 本试验的无人机多光谱影像数据采集以及田间实测冠层叶片SPAD值的采集时间选择为2021年1月13日和3月7日, 分别对应冬小麦的越冬期和拔节期。 在无人机数据采集时, 同时进行田间实测SPAD数据的采集。

田间实测叶片冠层SPAD数据的采集使用Konica Minolta SPAD-502plus型手持式叶绿素测量仪, 对试验田内96个小区进行测量。 在每个小区中, 随机选择10张顶层叶片。 在每张选择的叶片上等间距测量3个点的SPAD值, 求得平均值作为该植株的SPAD值, 5株小麦的平均值作为该小区的SPAD值。 每个生育期各有96个样本。

试验用的多光谱无人机遥感平台采用深圳大疆创新公司生产的精灵4多光谱版。 该无人机搭载6个镜头, 包括1个用于可见光成像的彩色传感器和5个用于多光谱成像的单色传感器(蓝(B): 450 nm; 绿(G): 560 nm; 红(R): 650 nm; 红边(RE): 730 nm; 近红外(NIR): 840 nm)。 使用该多光谱无人机能够采集农作物的光谱信息, 并且快速、 精准地判断作物的生长情况及健康水平等情况。 该无人机具有观测面积大的优点, 无论是针对单株植物还是整片农田, 均可高效完成数据采集。 影像的采集选择晴朗无风的时刻, 利用大疆创新公司的GS PRO软件进行常规航线的任务设置, 设置飞行高度为15 m, 飞行速度为3 m· s^-1, 航向和旁向重叠度均为80%。 无人机开始拍摄前, 在获取影像范围内按顺序摆放10组标准灰度梯度板, 以进行后续无人机影像的辐射定标^[9]。

1.3.1 辐射校正

获取的多光谱无人机数据利用大疆创新公司的大疆智图产品进行二维多光谱合成, 然后利用标准灰度梯度板对各波段图像进行辐射校正, 得到反射率数据。 辐射校正的过程是将DN值转化为实际物理意义的大气顶层辐射亮度或反射率, 它的原理是建立数字量化值与对应视场中辐射亮度值之间的关系, 以消除传感器本身产生的误差。 对获取的多光谱无人机影像进行辐射校正的主要过程是, 在视场中按顺序摆放10组梯度灰度板, 即稳定反射辐射特征的像元, 并且可以测得在相同天气状况下的反射率数据, 这一系列灰度梯度板在不同时相下的反射率和遥感图像之间存在的某种线性关系, 就可以通过波段运算对遥感图像进行辐射校正。

1.3.2 植被指数阈值法剔除土壤背景

处于生长早期的冬小麦由于叶片在图像中占比较小, 所以在遥感图像中多以土壤背景为主, 因此在进行冠层反射率提取时要消除土壤背景的影响。 本研究采用植被指数阈值法^[10], 通过计算整景图像的NDVI, 划分出土壤和小麦的阈值。 利用决策树分类和掩膜提取来剔除非冠层部分的信息(如图2)。 然后在ENVI中对各波段经过波段运算的图像进行感兴趣区构建和像元统计, 提取各小区的反射率均值作为样本在该波段的光谱反射率。

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 剔除土壤背景 (a): 原始图像; (b): NVDI阈值划分; (c): 剔除背景Fig.2 Removing soil background (a): Original image; (b): Threshold partition of NDVI; (c): Removing background

对不同波段反射率的进行组合, 可对地表植被冠层进行定性和定量的分析, 从而反映作物长势。 根据各植被指数的适用性, 利用Excel对提取的小区各波段反射率数据, 并选择15种植被指数进行计算, 一共得到20种遥感变量, 如表2。 为充分利用多光谱无人机平台能够获取植被红边信息的特点, 其中包含了6种红边参数的指数。

表2

Table 2

表2(Table 2)

表2 遥感变量及其计算公式 Table 2 Vegetation indices and calculating formulas 遥感变量 计算公式 参考文献 红光波段反射率 Red / 绿光波段反射率 Green / 蓝光波段反射率 Blue / 红边波段反射率 Red Edge / 近红外波段反射率 NIR / 归一化植被指数NDVI (NIR-R)/(NIR+R) [7] 红边重归一化植被指数RERDVI (NIR-Red Edge)/(NIR+Red Edge) [7] 叶绿素吸收比值指数CARI (Red Edge-R)-0.2× (Red Edge+R) [11] 优化土壤调节植被指数OSAVI (NIR-R)/(NIR+R+0.16) [12] 红边简单比SRred edge NIR/Red Edge [12] 核心红边三角植被指数RTVIcore 100× (NIR-Red Edge)-10× (NIR-G) [13] 红边叶绿素指数C I⊕red edge (NIR/Red Edge)-1 [13] 绿色归一化植被指数GNDVI (NIR-G)/(NIR+G) [14] 绿波段优化土壤调节植被指数GOSAVI 1.16× [(NIR-G)/(NIR+G+0.16)] [7] 改进简单比植被指数MSR (NIR/R-1)/(NIR/R+1) [14] 红边优化土壤调节植被指数REOSAVI 1.16× [(NIR-R)/(NIR+R+0.16)] [7] 归一化蓝绿差异指数NGBDI (G-B)/(G+B) [16] 归一化绿红差异指数NDRGI (G-R)/(G+R) [16] 陆地叶绿素指数MTCI (NIR-Red Edge)/(Red Edge+Red) [17] 增强植被指数2 EVI2 2.5× (NIR-Red)/(NIR+2.4× Red+1) [18]

表2 遥感变量及其计算公式 Table 2 Vegetation indices and calculating formulas

BP神经网络是一种后向反馈网络, 相较于多层感知机, 优化了模型中的系数和偏移量, 对非线性数据具有更强的映射、 自适应和泛化能力^[19]。 近年来, 在模式识别, 函数逼近、 分类、 数据挖掘、 数据压缩等领域都有较好的表现。

BP神经网络是一种有监督学习的算法。 它的核心思想是将每次根据训练得到的结果与预想结果进行误差分析, 进而修改权值和阈值, 逐步缩小误差, 得到输出和预想结果相对一致的模型。 以三层神经网络为例, BP神经网络是包含输入、 隐含和输出层三层结构的算法。 输入层接收数据, 输出层输出数据, 前一层神经元连接到下一层神经元, 收集上一层神经元传递来的信息, 经过“ 激活” 把值传递给下一层。 常用的激活函数有Sigmoid、 Tanh、 ReLu等。 其输出结果采用前向传播, 误差采用反向(Back Propagation)传播方式进行。 反向传播的基本思想就是利用梯度下降法原理, 通过计算输出层与期望值之间的误差来调整网络参数, 从而使得误差变小。

在构建BP神经网络过程中, 有数个较为关键的训练参数( 如训练函数、 目标函数、 激活函数种类和学习率等) 成为神经网络结构的超参数。 引入激活函数的目的是在模型中引入非线性。 一般来说, 在神经网络的中间层更加建议使用ReLu函数, 主要有两个原因。 一是ReLu函数计算简单, 可以加快模型速度; 二是避免当层数比较多的时候可能会造成梯度消失的现象, 从而模型无法收敛。

岭回归是一种线性模型, 它以最小二乘法为理论基础进行改良, 通过降低精度和损失部分信息为代价, 从而获得更可靠回归系数。 梯度提升树是以决策树为基函数的集成学习算法, 它的思想是弱学习器的合并, 通过提升模型的复杂度, 来使模型有更多机会纠正去训练集上的错误, 具有占用的内存小, 预测速度快的特点。

采用决定系数R²(coefficient of determination)和均方根误差RMSE(root mean squared error)来综合评价冬小麦SPAD值的反演模型精度。 R²越接近于1, RMSE越小, 说明模型精度越高, 反演估算能力越好。 计算公式如式(1)和式(2)。

R2=∑i=1n(xi-x¯)2(yi-y¯)2∑i=1n(xi-x¯)2∑i=1n(yi-y¯)2(1)

RMSE=∑i=1n(xi-yi)2n(2)

其中, x_i是冬小麦SPAD的实测值, x̅是SPAD的实测值的平均值; y_i为估算模型的SPAD估算值, y̅为SPAD估算值的平均值; n是样本的数目。

本研究基于Python3.8, 利用数学分析包Numpy和Pandas库进行数据分析, 计算出了各变量与实测SPAD值之间的Pearson相关性系数(表3)。 结果表明, 在越冬期, 相关性较低的遥感变量较多, 其中蓝波段反射率与SPAD的相关性最低, 红边重归一化植被指数(RERDVI)的相关性最高, 达到0.81。 在拔节期, 其中大部分相关性绝对值都在0.6以上, 相关性较好, 核心红边三角植被指数(RTVI core)、 近红外波段和SPAD值的相关性最大, 相关性绝对值在0.84以上。 整体上遥感变量和SPAD值相关性呈极显著关系, 拔节期的光谱指数与SPAD值的相关性整体上优于越冬期。

表3

Table 3

表3(Table 3)

表3 冬小麦SPAD值与遥感变量的相关系数 Table 3 Correlation between spectral indices and SPAD in different growth periods of winter wheat 遥感变量 相关系数 越冬期 拔节期 Red 0.53* -0.75* Green 0.49 -0.55* Blue 0.005 -0.73* NIR 0.71* 0.84* * NDVI 0.47 0.8* * OSAVI 0.58* 0.81* * GNDVI 0.31 0.83* * GOSAVI 0.47 0.83* * MSR 0.47 0.8* * REOSAVI 0.58* 0.81* * NGBDI 0.41 0.57* NDRGI 0.27 0.74* EVI2 0.63* 0.81* * Red Edge 0.67* 0.77* RERDVI 0.81* * 0.83* * CARI 0.47 0.83* * SRrededge 0.8* * 0.83* * RTVIcore 0.78* 0.84* * CIrededge 0.8* * 0.83* * MTCI 0.78* 0.83* * 注: * 和* * 分别表示差异达显著(p< 0.05)和极显著(p< 0.01)水平 Note: Note: * and * * represent that the correlations are significant at p< 0.05 and p< 0.01 levels, respectively

表3 冬小麦SPAD值与遥感变量的相关系数 Table 3 Correlation between spectral indices and SPAD in different growth periods of winter wheat

对冬小麦的越冬期、 拔节期的20个遥感变量分别进行特征重要性排序, 为模型构建之前的特征筛选提供参考依据。 如图3所示, 按照特征重要性的阈值对排序中的特征进行筛选, 确定0.05为重要性阈值。 在越冬期, 得到7个最优的遥感变量, 分别是RTVIcore、 NDRGI、 Red、 SRrededge、 NGBDI、 NIR和MTCI; 在拔节期, 得到了6个最优的遥感变量, 分别是RTVIcore、 NIR、 GOSAVI、 CIrededge、 RERDVI和MTCI。 在两个生育期, RTVIcore(核心红边三角植被指数)均是与SPAD值最重要的特征, 体现了利用多光谱无人机监测小麦长势的有效性。

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 特征重要性排序 (a): 越冬期; (b): 拔节期Fig.3 Feature importance sort (a): Overwinter stage; (b): Jointing stage

提取冬小麦两个生育期的多光谱反射率数据与之对应的小区地面实测数据, 作为训练样本数据集, 每个波段共得到96组数据, 随机选取80%的数据作为训练集, 构建冬小麦SPAD值反演模型, 20%作为测试集, 用于模型评价。

为了提升回归预测模型的表现, 通常采用减少特征数量的方法, 称为特征变量筛选或降维。 这里采用基于交叉验证的递归特征消除法(recursive feature elimination CV)。 递归特征消除法是一种寻找最优特征子集的贪心算法^[20]。 评估器选择随机森林(Random Forest), 根据特征重要性属性排序, 以特征重要性平均值0.05作为阈值, 选择此时期与因变量(SPAD值)的特征重要性在此阈值之上的遥感变量作为自变量, 进行建模。

在特征筛选之后, 选择岭回归(ridge), 梯度提升树(gradient boosting decision tree)和BP神经网络(back propagation)等3种机器学习算法, 建立冬小麦SPAD回归预测模型。

在模型构建时, 通过调节模型参数可以提升模型的泛化能力。 在岭回归算法中, 影响模型表现最大的因子是正则化系数alpha, 绘制岭迹图可以确定最佳正则化系数。 梯度提升树是迭代的决策树算法, 在树模型中, 最重要的是对n_estimators和max_depth参数的搜寻, 采用交叉验证和网格搜索的方式, 来绘制学习曲线来确定最优参数。

BP算法采用3层神经网络结构, 以各时期筛选的遥感变量作为输入层节点, 小麦叶片SPAD值为输出层, 建立包括输入层、 隐藏层和输出层的 BP 网络。 输入向量首先经过了归一化处理, 初始权重和阈值为任意值。 BP 网络的具体参数见表4。

表4

Table 4

表4(Table 4)

表4 BP神经网络超参数设置 Table 4 Parameters of BP algorithm 参数 设置 优化器 Adam 输入层到隐层激活函数 ReLu 隐层到输出层激活函数 Linear Epoch 200 学习率 0.01 隐含层层数 2

表4 BP神经网络超参数设置 Table 4 Parameters of BP algorithm

将两生长期筛选的光谱信息通过三种机器学习算法建模, 并在测试集上验证模型表现, 结果如表5、 表6和图4。 结果表明, 适应在训练集上建立的三种估算模型, R²和RMSE的结果基本都在合理范围内, 说明这些模型对于监测冬小麦SPAD值是可行的。 其中, 基于随机森林算法的递归特征消除筛选得到的遥感变量, 再经BP神经网络训练的反演模型(RF-RFE-BP)在两个生育期都表现出最佳的学习和预测能力。

表5

Table 5

表5(Table 5)

表5 冬小麦冠层叶片SPAD值估算反演模型比较 Table 5 Comparison of inversion model of SPAD in winter wheat 生育期 筛选 变量 模型 训练集 测试集 R2 RMSE 样本 R2 RMSE 样本 越冬期 RFE 7 Ridge 0.550 2.256 76 0.518 2.687 20 GBD 0.618 2.080 0.584 2.724 BP 0.813 1.512 0.806 1.861 拔节期 6 Ridge 0.757 0.543 0.653 0.576 GBD 0.675 0.636 0.488 0.749 BP 0.835 0.499 0.827 0.507

表5 冬小麦冠层叶片SPAD值估算反演模型比较 Table 5 Comparison of inversion model of SPAD in winter wheat

表6

Table 6

表6(Table 6)

表6 模型验证结果 Table 6 Model validation results 时期 算法 模型验证回归关系 (x为实测值, y为预测值) 样本数 Ridge y=0.466 3x+20.368 越冬期 GBD y=0.466 4x+21.271 BP y=0.717 9x+10.085 20 Ridge y=0.758 5x+8.224 5 拔节期 GBD y=0.7111 x+9.8546 BP y=0.892 6x+3.255 2

表6 模型验证结果 Table 6 Model validation results

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 验证模型分析散点图 (a): 越冬期小麦叶片SPAD实测值与模型预测值关系; (b): 拔节期小麦叶片SPAD实测值与模型预测值关系Fig.4 Scatter plot used for model verification (a): Relationship between measured SPAD and prediction values for wheat leaves in overwinter stage; (b): Relationship between measured SPAD and predicted values for wheat leaves in jointing stage

其中, 从单一生育期来看, BP神经网络算法的预测结果最优, 在越冬期, 在测试集上的效果R²达到了0.806, 相较于比岭回归提升了35%, 比梯度提升树提升了24%; RMSE为1.861, 比岭回归降低了30%比梯度提升树降低了32%; 在拔节期, 测试集的R²达到了0.827, 相较于比岭回归提升了21%, 比梯度提升树提升了40%; RMSE为0.507, 比岭回归降低了13%, 比梯度提升树降低了32%。 实测值和预测值之间的回归曲线与零误差线的偏角也是最小的。 整体来看, 三种机器学习算法构建的模型, 在拔节期都比越冬期的模型效果要更好, 预测的准确率更高。 由表4可知, 拔节期输入的参数量比越冬期少, 但反演模型的精度和稳定性略有提高, 说明, 输入模型训练的特征并不是越多越好, 经过了最优遥感变量筛选减少了冗余的特征, 保留了信息含量大且重要性较高的特征参数。

综上, BP神经网络构建的估算模型更适合用于多光谱无人机遥感估测冬小麦冠层SPAD值。 相较于传统的对冬小麦SPAD值的反演研究, 只涉及到单一或很少的品种^{[21, 22, 23]}, 不适合用于品种筛选实验。 RF-RFE-BP模型对不同施氮水平和不同品种均取得了较好的估测效果, 因为模型的训练是在品种筛选的前提下进行的, 此模型可以用于所有涉及的小麦品种。

在遥感变量选择方面, 应用到的方法有逐步回归法, 或直接根据自变量和SPAD值之间的相关性系数高低进行筛选^{[6, 7]}。 逐步回归是最常用也是最简单的特征选择方法, 它能够比较直接地删除回归过程中的冗余特征, 但是往往通过逐步回归分析得到的反演精度并不是很高, 回归方程也不是最优解。 若直接根据自变量与因变量之间的相关性进行特征选择, 就需要考虑变量之间是否存在多重共线性。 在通过遥感手段获取特征参数的过程中, 特征数量越多, 影响预测精度。 前人在对森林生物量的估算研究中提出一种后向迭代的特征选择方法, 这与本研究的变量筛选原理一致。 首先从数据特征预处理的角度出发, 利用递归特征消除法(RFE), 选择随机森林(RF)作为评估器, 通过计算特征重要度, 采用后向迭代的方法逐步简化特征数量, 对20种遥感变量进行特征筛选。 进行特征变量筛选的目的, 一是为了剔除相关性不高的冗余变量, 二是为了提升模型表现和建模速度。 光谱信息过少的模型容易导致模型受到背景因素的干扰而缺乏稳定性, 光谱信息过多的模型会导致模型复杂且易用性低, 因此需要平衡两者之间的关系^[23], 因为所选择计算的植被指数具有一定的主观性, 并不是所计算的遥感变量与SPAD值之间都具有非常重要的指示作用。

现有的研究中, 对作物冠层叶片的SPAD值进行遥感反演的算法有很多, 如随机森林、 偏最小二乘、 多元线性回归等^{[4, 24, 25]}, 均得出了相应的结论。 结果表明, BP在两个生育期对于冬小麦SPAD值的估测能力同其他两种算法来说, 具有很好的优势。 经过对比发现, 在受土壤背景影响较大的越冬期, 该模型仍然具有估算精度高, 容错能力强的特点, 对于越冬期之后的一段生育时期的SPAD值预测, 也具有很好的可用性。

其中, 岭回归在两期数据集上的表现能力都较为稳定, 但线性模型对于该问题的预测, 准确率较低^[26], 它的应用范围有限, 难以处理复杂的数据之间的映射关系。 以往基于小麦近地多光谱图像对山农15和泰农18两个品种的越冬期数据进行SPAD估测, 结果显示在越冬期构建的线性模型的决定系数为0.705 3。 以集成学习为代表的梯度提升树算法, 在预测准确度上, 相比于本研究所选用的线性模型岭回归, 集成学习模型的预测表现有所提升。 这与马文君^[27]等在棉花冠层SPAD值预测研究中的结果一致。 在两个生育期, BP都达到了最佳的效果, 其原因可能是, 其一, 经过特征筛选后, 遥感变量数量较少, 其他两种模型的稳定性发生了显著变化。 其二, 这得益于神经网络算法对于非线性问题的映射能力更强的优点。 BP神经网络构建的SPAD值估算模型在各生育期上的效果都是最好的, 当特征变量数较少时, 模型仍然具有很强的泛化能力, 在训练集和测试集上的效果相差不大。