基于无人机高光谱遥感的水稻氮营养诊断方法

许童羽 白驹驰 郭忠辉 金忠煜 于丰华

(1.沈阳农业大学信息与电气工程学院，沈阳 110866；2.辽宁省农业信息化工程技术研究中心，沈阳 110866)

0 引言

氮素是对作物生长、发育，以及最终产量影响最明显的元素，其含量变化会对光合作用、蛋白质合成以及碳氮代谢产生影响。氮素的缺乏会抑制作物地上部分和根系的生长，限制繁殖器官的形成和发育，并显著影响作物最终产量以及品质[1-2]。因此，快速、精准、大面积地对水稻田间氮素需求情况进行诊断，并依据诊断结果实现精准施肥，是实现水稻田间精准管理和保证水稻产量的重要手段。传统的水稻氮营养诊断标准主要基于临界氮浓度曲线(Critical N concentration curve)。GREENWOOD等[3]最早提出了临界氮浓度的概念，即作物达到最大生物量所需的最小氮浓度。随后LEMAIRE等[4]提出了基于临界氮浓度理论的氮营养指数(Nitrogen nutrition index，NNI)。氮营养指数作为作物氮营养诊断的重要指标，能够定量描述作物的氮营养丰缺状况[5-6]，为作物氮亏缺量的确定提供了理论依据。吕茹洁等[7]基于临界氮浓度理论，计算出超级杂交稻和常规稻的临界氮浓度曲线，并依据曲线确定了二者的氮亏缺量与适宜施肥量，为水稻氮亏缺量的确定提供了理论支撑。

然而，基于氮营养指数计算氮亏缺量需要通过田间采样获取数据，成本较高、测量周期长，且结果具有滞后性[8]，难以对实际农业生产做出指导。近年来，随着高光谱遥感技术的发展，利用无人机高光谱遥感技术获取水稻氮营养状况信息成为精准农业领域的重要发展方向。高光谱遥感技术能够通过分析作物光谱数据来获得其生长信息，具有快速、无损、准确等优点[9]。对此，国内外已在相关领域取得一定的研究进展[10-11]。这些研究大多数集中在利用高光谱对水稻氮含量进行反演，而单纯的氮含量信息难以反映水稻的氮营养丰缺状况，因此，结合无人机高光谱遥感技术与临界氮浓度曲线理论，通过光谱反演氮营养指数的形式对水稻进行氮营养诊断成为当前的研究热点[12-13]。ZHA等[14]通过无人机遥感获取水稻冠层光谱数据，构建多种植被指数并结合机器学习算法对水稻地上生物量、氮吸收量与氮营养指数进行反演建模，结果显示随机森林算法显著提高了氮营养指数反演精度；QIU等[15]提取水稻RGB光谱信息并构建植被指数，最后比较多种机器学习方法反演水稻氮营养指数的精度，认为随机森林算法在各生育期表现最佳。

前人研究大都以对水稻进行氮营养诊断为主，而单纯的氮营养诊断只能判断水稻的生长状况，需要进一步获取水稻氮亏缺量数据才能以此对精准施肥做出定量指导。因此，本文拟以水稻无人机高光谱数据与氮亏缺量为研究对象，以水稻各个时期氮亏缺量约等于0的冠层光谱为标准，对水稻高光谱数据进行比值、差值与归一化差值变换，然后利用竞争性自适应重加权采样法处理原始光谱与变换光谱，筛选出效果最佳的光谱特征波段，并建立基于多元线性回归、极限学习机与蝙蝠算法优化极限学习机的水稻氮亏缺量反演模型，最后比较模型反演效果，确定基于高光谱数据的水稻氮亏缺量最佳反演方法，以期为快速获取水稻氮亏缺量提供研究方向，也为基于田间水稻氮营养状况的精准施肥提供定量指导。

1 数据获取与反演模型建立

1.1 实验设计

实验于2021年6—9月在辽宁省鞍山市海城市耿庄镇沈阳农业大学精准农业航空科研基地(北纬40°58′45.39″，东经122°43′47.01″)进行，为避免阴雨多云天气对遥感数据采集造成误差，数据采集过程选择天气晴好的日期，如遇云量超过20%或不利于遥感数据采集天气则顺延采集。水稻品种为北粳1705。实验区域如图1所示，实验田分为2个大区，其中实验区1设立5个氮肥梯度，施氮量分别为N0(0 kg/hm2)、N1(75 kg/hm2)、N2(150 kg/hm2)、N3(225 kg/hm2)、N4(300 kg/hm2)，实验区2也设计5个氮肥梯度，施氮量分别为：N0(0 kg/hm2)、N1(50 kg/hm2)、N2(100 kg/hm2)、N3(150 kg/hm2)、N4(200 kg/hm2)，氮肥基追比为5∶3∶2。实验区1对每个梯度设立3个重复小区，共设立3×5=15个小区。每个小区面积为5 m×8 m=40 m2，实验区2每个小区面积为660 m2，除氮肥梯度外两组实验区田间管理一致，磷钾肥施用量采用当地标准施用量进行，其中磷肥标准施用量为144 kg/hm2，钾肥标准施用量为192 kg/hm2，基追比为1∶1，其余田间管理同常规高产管理。田间采样由分蘖期至抽穗期，采样间隔9 d，每次采样在各实验小区选取具有代表性的3穴水稻进行冠层光谱、水稻氮浓度与地上干物质量的获取，最后结果取平均值作为该小区水稻氮浓度与地上干物质量。实验区1共采集小区数据120组，实验区2共采集小区数据88组，共计208组样本。

图1 实验小区分布图

1.2 数据获取

1.2.1水稻光谱参量获取

无人机高光谱平台采用深圳大疆创新公司M600 PRO型六旋翼无人机，高光谱成像仪选用四川双利合谱公司GaiaSky-mini内置推扫式机载高光谱成像系统，高光谱波段范围为400～1 000 nm，分辨率为3 nm，有效波段数为253个。无人机高光谱遥感平台数据采集时间为选择太阳光强较为稳定的时段11:00—12:00，无人机飞行高度为150 m。运用 ENVI 5.3+IDL工具软件对获取的高光谱遥感影像进行小区高光谱数据提取，运用波谱角填图方法去除干扰地物光谱的影响，对每个小区的感兴趣区计算平均光谱，再利用Matlab软件对平均光谱进行重采样处理，使光谱分辨率降为1 nm，最后利用高斯滤波器对重采样后光谱进行去噪处理，结果作为每个实验小区的高光谱信息。

1.2.2水稻农学参量获取

测定样本氮浓度与地上干物质量时，首先对每个小区的水稻进行破坏性采样，将样本带至实验室，然后将样本放入干燥箱中105℃杀青30 min后，80℃干燥至恒质量；之后测量干燥后样本的干物质量，并根据种植密度换算成地上部干物质量；最后将干燥后的样本磨碎，对研磨后的样本通过凯氏定氮法测量植株氮浓度。

1.3 基于临界氮浓度曲线的氮亏缺量计算方法

临界氮浓度是作物达到最大生物量所需的最小氮浓度，因此以临界氮浓度为基准可以推导出临界氮积累量方程，进而建立氮亏缺量方程，并确定水稻氮亏缺量[7]。

根据JUSTES等[16]提出的临界氮浓度曲线构建方法，以所测水稻氮浓度与地上干物质量为基础构建东北水稻的临界氮浓度曲线，具体方法如下：对比不同氮肥梯度下测定的氮浓度与干物质量，通过方差分析将样本分为受氮营养限制组与不受氮营养限制组；对受氮营养限制组的数据，将地上干物质量与氮浓度进行线性拟合；对不受氮营养限制组，对同一时期样本地上干物质量取平均值代表该时期最大地上干物质量；取各时期拟合曲线于最大地上干物质中位置的截点为该时期的理论临界氮浓度点；将各时期临界氮浓度点进行幂函数拟合，构建作物临界氮浓度曲线

Nc=aM-b

(1)

式中Nc——水稻临界氮质量比，g/g

M——地上干物质量，t/hm2

a、b——曲线参数

根据临界氮浓度曲线推导出临界氮积累量方程与氮亏缺量方程，推导过程为

Ncna=10aM1-b

(2)

Nand=Ncna-Nna

(3)

式中Ncna——临界氮浓度条件下植株氮积累量，kg/hm2

Nand——氮亏缺量，kg/hm2

Nna——不同施氮量下植株实际氮积累量，kg/hm2

1.4 光谱数据处理

1.4.1高光谱数据变换

由于临界氮浓度状态下作物地上部氮营养状况与冠层叶片结构达到最佳状态，而这些都是影响水稻冠层光谱反射率的直接因素[17-18]，因此，相比非临界氮浓度状态下的作物，临界氮浓度下水稻同其余水稻无论在光谱层面还是氮营养状况层面都存在一定的差异。为了在光谱层面突出临界氮浓度下水稻与其余水稻的差距，本研究选取各个时期Nand≈0 kg/hm2样本对应的光谱反射率，将其与同时期其余样本光谱反射率分别做比值、差值和归一化差值变换，计算方法为

(4)

RD=RNC-Ri

(5)

(6)

式中RSV——比值光谱反射率，%

RD——差值光谱反射率，%

RND——归一化差值光谱反射率，%

RNC——临界氮浓度状态下样本光谱反射率，%

Ri——同时期其余光谱反射率，%

1.4.2高光谱特征波段提取

由于全波段光谱相邻波段数据相似度较高，且存在大量与所求变量无关的冗余信息，基于全波段光谱的建模往往存在运行速度较慢，模型误差较高，反演精度较低等缺点[19-21]。为了减少输入变量个数，降低数据冗余性，提高建模速度与精度，本文采用竞争性自适应重加权采样法(Competitive adaptive reweighted sampling，CARS)对光谱反射率数据进行特征提取。CARS算法基于达尔文进化论简单有效的生存原则“适者生存”，对不适应的波长变量进行逐步淘汰。首先，利用蒙特卡罗采样法(Monte Carlo simulation，MCS)采样N次，再以偏最小二乘回归(Partial least squares regression，PLSR)系数绝对值作为衡量标准，保留回归系数绝对值大的波长变量，去除回归系数绝对值小的波长变量，从而获得一系列的波长变量子集；然后对比每次产生的PLSR模型的交互验证均方差值(Root mean square error of cross validation，RMSECV)，RMSECV 最小的那个模型所对应的变量子集被选为最优变量子集。

1.5 反演建模方法

选用多元线性回归(Multivariable linear regression，MLR)、极限学习机(Extreme learning machine, ELM)与蝙蝠算法优化极限学习机(Bat algorithm optimized extreme learning machine, BA-ELM)3种算法进行建模，并通过3种模型对训练集与测试集预测结果的决定系数R2与均方根误差(Root mean square error，RMSE)来判断模型的反演精度与鲁棒性。

多元线性回归是一种解释单因变量与多自变量间线性关系的传统回归建模方法，模型结果唯一，并能够更好地体现各个自变量与因变量间的相关程度。

极限学习机是一种基于单层前馈神经网络理论的改进算法，具有学习速率快、泛化能力强、训练精度较高等优点[22-23]。但由于极限学习机每次运行中输入层与隐含层之间的连接权值和隐含层神经元的阈值都是随机生成的，其在训练稳定性上存在不足，且容易陷入局部最优解。因此，本研究采用蝙蝠算法对极限学习机隐含层初始权重进行优化。

蝙蝠算法是一种新的元启发式优化算法。该算法模拟自然界中蝙蝠利用回声定位来探测猎物和规避障碍物的原理，控制蝙蝠发出不同频率的搜索脉冲来寻找最优解。相较于经典的粒子群算法，蝙蝠算法引入了局部搜索机制，在每轮迭代中最优蝙蝠个体会采用随机搜索的方式进行局部寻优，有助于算法跳出局部最优解，增强算法的全局搜索性[24-26]。

2 结果与分析

2.1 植株氮亏缺量确定

2.1.1植株氮浓度与地上部干物质量统计分析

东北水稻全生育期植株氮浓度与地上部干物质量的基本信息如表1所示，由于实验包含全生育期采样数据，实验区1与实验区2中的植株氮浓度与地上部干物质量离散程度较高，最大值与最小值差别较大。其中实验区1由于施氮梯度较高，其整体氮浓度与地上部干物质量高于实验区2。实验区2由于在分蘖期与拔节期采样密度较高，其平均氮浓度较高，平均地上部干物质量较低。

表1 植株氮浓度与地上部生物量统计

2.1.2临界氮浓度曲线与氮亏缺量计算

根据1.3节提出的临界氮浓度曲线的构建方法，本文将每个采样日获得的氮浓度与干物质量进行回归拟合，计算出每个采样日的东北水稻临界氮浓度，之后根据各个临界氮浓度与对应的干物质量构建东北水稻临界氮浓度曲线(图2)。曲线方程为

图2 临界氮质量比拟合曲线

Nc=2.026M-0.460 3

(10)

决定系数R2为0.874 7，RMSE为0.382 5 g/g。曲线参数a、b分别为2.026、-0.460 3，与SONG等[27]计算出的东北水稻临界氮浓度曲线参数a=1.99、b=-0.44非常接近，可进一步用于东北水稻的氮亏缺量计算。

根据1.3节中的计算方法，本文结合样本数据与临界氮浓度曲线，计算各样本的氮亏缺量。各个氮肥梯度的平均氮亏缺量如图3所示，随着氮肥梯度的增高，氮亏缺量整体呈下降趋势，综合实验区1与实验区2氮亏缺量的下降趋势可知，水稻最佳施肥量在150～200 kg/hm2之间。

图3 各氮肥梯度平均氮亏缺量

2.2 光谱数据处理

2.2.1光谱数据变换结果

根据1.4.1节提出的方法，本研究选取各个时期Nand≈0 kg/hm2的样本对应光谱反射率，将其与其余样本光谱反射率分别作比值、差值、归一化差值变换，变换结果与相关性分析结果如图4、5所示，对比可得，不同施氮水平下原始光谱全波段反射率变换趋势基本一致，在可见光波段(波长400～700 nm)光谱反射率均呈“先升后降”的变化规律。同时，在绿波段(波长550 nm附近)出现明显的反射峰，在红波段(波长680 nm附近)出现明显的吸收谷。在波长680～770 nm内不同施氮水平的光谱反射率曲线变化基本一致，反射率急剧升高，在近红外波段(770～1 000 nm)形成较高的反射平台；不同施氮水平下比值、差值与归一化差值光谱反射率有一定差别，但在绿波段与红波段差别较为明显;在与氮亏缺量的相关性方面，4种光谱反射率与氮亏缺量的相关性图像趋势基本一致，在红波段与绿波段形成2个波峰，其中归一化差值光谱在2个波峰处相关性最佳，其次分别是比值光谱、差值光谱与原始光谱。

图4 原始与变换光谱反射率

2.2.2特征波长选择

利用CARS算法对水稻原始光谱与归一化差值光谱反射率进行特征波长的筛选，根据交叉验证集RMSECV确定光谱特征波长后，删去相邻的特征波长，并通过相关性分析选取与Nand相关性最好的波长作为最佳特征波长，以原始光谱选择过程为例：由图6a可以看出，随着运行次数的增加，原始光谱被选出的波段数逐渐减少，图6b中RMSECV整体呈下降趋势，说明筛选过程中剔除的变量与SOM去除量无关，而72次迭代以后，RMSECV呈回升趋势，表明反射率光谱中与SOM无关的大量信息或噪声被添加，从而导致RMSECV上升。图6c为所有变量在每次采样过程中的回归系数路径图，表示随着运行次数的增加各波段变量回归系数的变化趋势。结合图6b分析发现当运行次数为第72次时，RMSECV最小即所选择的光谱变量子集最优。原始光谱与3种变换光谱的特征选择结果如表2所示。

图5 原始与变换光谱相关性

图6 原始光谱CARS降维结果

表2 CARS提取后原始光谱与变换后光谱特征波长与运行结果

2.3 氮亏缺量反演模型构建

为比较不同光谱变换方法与反演建模方法的性能，分别以经CARS算法筛选出的原始光谱、比值光谱、差值光谱、归一化差值光谱特征波长为输入变量，水稻氮亏缺量为输出变量，分别构建基于MLR、ELM与BA-ELM的水稻氮亏缺量反演模型。

2.3.1基于MLR的氮亏缺量反演

将4种光谱特征波长作为输入变量，分别构建基于MLR的反演模型，反演结果如图7所示：其中基于归一化差值光谱的模型精度最高，训练集与测试集的R2分别为0.654 0、0.715 8，RMSE分别为0.955 7、1.204 8 kg/hm2；其次是基于差值光谱反演模型，其R2分别为0.600 5、0.524 2，RMSE分别为1.103 1、1.184 2 kg/hm2；然后是基于原始光谱的反演模型，R2分别为0.595 6、0.582 6，RMSE分别为1.189 4、1.083 9 kg/hm2。基于比值光谱的模型精度精度R2分别为0.597 3、0.517 5，RMSE分别为1.034 0、1.064 6 kg/hm2；基于归一化差值光谱的反演模型反演效果最好。

图7 基于MLR的不同反演模型结果

2.3.2基于ELM的氮亏缺量反演

将4种光谱特征波长作为输入变量，分别构建基于ELM的反演模型，映射函数为Sigmoid，隐含层个数分别为13、14、12、13，训练结果如图8所示。其中基于归一化差值光谱的模型精度最高，训练集与测试集的R2分别为0.721 1、0.683 2，RMSE分别为0.857 8、1.151 2 kg/hm2；其次是基于比值光谱的反演模型，R2分别为0.685 0、0.634 4，RMSE分别为0.914 5、0.958 4 kg/hm2；然后是基于差值光谱反演模型，其R2分别为0.632 2、0599 4，RMSE分别为1.057 0、1.069 9 kg/hm2；最后是原始光谱的模型精度最低，R2分别为0.619 2、0.547 7，RMSE分别为1.154 6、1.175 3 kg/hm2。基于归一化差值光谱的反演模型R2与RMSE均优于其余3种模型。

图8 基于ELM的反演模型结果

2.3.3基于BA-ELM的氮亏缺量反演

将4种光谱特征波长作为输入变量，分别构建基于BA-ELM的反演模型，其中蝙蝠算法的最大迭代次数为80，种群大小为50，适应度函数为验证集RMSE，隐含层个数分别为50、60、57、65，结果如图9所示：基于归一化差值光谱的模型精度最高，训练集与验证集的R2分别为0.816 0、0.830 6，RMSE为0.696 8、0.814 1 kg/hm2；基于比值光谱的模型精度R2分别为0.759 3、0.770 1，RMSE分别为0.799 4、0.730 9 kg/hm2；基于差值光谱的模型精度R2分别为0.733 9、0.659 0，RMSE分别为0.900 3、0.997 1 kg/hm2；原始光谱的模型精度，训练集与验证集的R2分别为0.692 2、0.680 8，RMSE为1.037 5、0.952 8 kg/hm2。经过模型结果对比可知：基于BA-ELM的反演模型精度均优于基于MLR与ELM的反演模型；其中基于归一化差值光谱的反演模型在R2与RMSE上均优于其余3种光谱模型；基于比值光谱在MLR模型中表现较差，但在ELM与BA-ELM模型中表现仅次于归一化差值光谱，这可能是因为比值光谱与氮亏缺量之间存在较为复杂的非线性关系，MLR等线性模型对它的解释能力较差。在3种模型中，差值光谱与原始光谱表现相近，可见对光谱进行单纯的差值处理难以突出其与氮亏缺量的关系。因此采用基于归一化差值光谱的BA-ELM模型在反演水稻氮亏缺量上最有优势。

图9 基于BA-ELM的反演模型结果

2.4 针对不同施肥量的模型预测能力

在实验区1与实验区2中，以反演效果最优模型的氮亏缺量估测值与对应梯度的氮亏缺量实测值对比，结果如图10所示，实验区1与实验区2中，整体上估测值与实测值较为相近，对水稻氮亏缺量的估测能力较强；随着施氮量梯度的提升，对氮亏缺量的估测精度呈下降趋势，在临界氮浓度状态附近估测能力较差，这可能是因为临界氮浓度状态下与非临界氮浓度状态下水稻光谱变化规律不一致，模型同时估测2种状态的水稻光谱预测难度较高所致。

图10 不同梯度氮亏缺量估测效果

3 讨论

以往学者研究焦点往往集中在如何提高获取水稻田间氮浓度或氮营养指数信息的准确性[28-29]，然而这些信息只能描述水稻的氮营养状况，难以对水稻氮亏缺状况进行精确的定量描述，进而指导精准施肥。本文首先通过田间采样获取水稻氮浓度、地上生物量与冠层光谱数据，根据前人提出的方法构建东北水稻临界氮浓度曲线，并基于临界氮浓度曲线计算水稻氮亏缺量。由于临界氮浓度表示作物达到最大干物质量所需的最小氮浓度，以临界氮浓度作为判断水稻氮亏缺量的标准更能反映田间水稻的氮营养状况与实际氮需求情况。同时本研究还采用了无人机高光谱遥感技术与数据驱动的反演建模方法，通过水稻冠层光谱数据对氮亏缺量进行反演。无人机高光谱遥感技术能够快速、大面积获取水稻冠层光谱信息，进而对田间水稻氮营养状况进行快速诊断，数据驱动的反演建模方法虽然可解释性与泛化能力较差，但其能够对光谱与氮亏缺量之间的复杂非线性关系进行准确描述，将二者结合能够同时满足追肥决策所需的精确性与时效性，为定量追肥决策的实现提供了解决思路。

由于临界氮浓度状态下作物地上部氮营养状况与冠层叶片结构达到最佳状态，而这些都是影响水稻冠层光谱反射率的直接因素[30]，因此，相比非临界氮浓度状态下的作物，临界氮浓度下作物无论在光谱层面还是氮营养状况层面都存在一定差异。为了在光谱层面上放大临界氮浓度下水稻与其余水稻的差距，去除与氮亏缺量无关的冗余信息，本研究以各时期接近临界氮浓度状态下光谱为标准光谱，分别对原始光谱进行比值、差值、归一化差值变换，之后运用CARS算法提取变换光谱的特征波长。由结果可知3种光谱变换均提升了其与氮亏缺量之间的相关性，其中归一化差值变换相关性与反演精度最高，这可能是由于归一化差值计算既通过差值计算放大了临界氮浓度光谱与其余冠层光谱的差异，又统一了差值光谱数据的量纲，还对光谱数据的统计分布性进行了归纳。

在建模过程中，本文选择蝙蝠算法优化ELM算法的初始权重，并将其与原始ELM进行对比建模。建模结果中BA-ELM模型的反演精度显著高于ELM模型，这可能是由于ELM作为一种基于前馈神经网络的机器学习算法，具有学习速度快，收敛能力强等优点，但同时也存在容易陷入局部最优解的问题；蝙蝠算法作为一种由粒子群算法衍生的元启发式搜索算法，在迭代过程中会在最优解附近再次搜索局部最优解，因此具有优秀的全局搜索能力，通过蝙蝠算法优化极限学习机的初始参数会让极限学习机在最优解附近开始学习，提高了算法的反演精度与全局搜索能力。

虽然本研究利用无人机高光谱数据对水稻氮亏缺量反演取得了较好的效果[31]，但对于实现基于无人机高光谱遥感的水稻精准追肥仍存在以下几点问题：①本文采用的归一化差值处理方法是基于各个采样期临界氮浓度状态的水稻光谱计算的，而这些光谱如果用于其他水稻光谱的计算时，可能会因为采样时间、水稻品种与发育状况等因素造成误差，进而影响氮亏缺量反演效果。②蝙蝠算法优化后的ELM算法也存在收敛精度不高，容易陷入局部最小值等问题。后续的研究重点应集中在如何在无人机高光谱数据中实时提取临界氮浓度状态下的水稻光谱与对反演算法的精度与鲁棒性进行优化这两方面上。

4 结论

(1)构建了东北水稻临界氮浓度曲线，可表示为Nc=2.026M-0.460 3。

(2)通过CARS算法提取的原始光谱、比值光谱、差值光谱与归一化差值光谱特征波长分别为400、411、422、447、570、677 nm，533、695、715、738、774 nm，413、424、444、497、533、600、676、691 nm，484、560、644、672、695、724 nm。

(3)对比基于4种光谱特征波长与ELM、BA-ELM算法的8种建模方法反演精度可知，基于归一化差值的BA-ELM模型反演结果最好，训练集与测试集的R2分别为0.816 0、0.830 6，RMSE为0.696 8、0.814 1 kg/hm2。