冷杉叶片氮含量高光谱反演技术研究

宋雪莲,王志伟,张文，2，张威，丁磊磊，柳嘉佳，阮玺睿,王普昶

(1.贵州省农业科学院草业研究所，贵州 贵阳 550006；2.贵州阳光草业科技有限责任公司，贵州 贵阳 550006；3.贵州省水利水电勘测设计研究院有限公司,贵州 贵阳 550002；4.贵州省农业科学院,贵州 贵阳 550006)

氮是植物生长的基本营养元素，是陆地生态系统植物生长的主要限制因子，与植物的光合作用和细胞生长分类等重要生理活动相关[1]。缺氮会影响植物的光合作用能力，导致叶片生理及形态结构的相应变化，因而引起叶片光谱反射特性的变化[2]。传统的化学方法能够精确地监测出植物氮含量，但具有高损耗、复杂、时滞性等缺点，不能进行全面、大范围、快速的营养诊断。高光谱光谱分辨率高，能够获得连续的光谱信息，已成为监测植物氮含量的有效手段。

目前有不少学者采用敏感波段筛选和植被指数进行高光谱植物氮含量检测的相关研究。朱艳等[4]研究表明，单波段光谱在610 nm和680 nm处的水稻冠层反射率与叶片氮含量具有较好的相关性，提出采用回归系数来提高水稻叶片氮含量估测的准确性。刘冰峰等[5]指出720 nm处的反射光谱以及550，720 nm处的一阶光谱导数是夏玉米氮含量的敏感参数。Tarpley等[6]研究认为，利用红边位置和短波近红外波段光谱反射率比值可预测棉花叶片氮浓度。Pacheco-Labrador等[7]采用归一化指数NDIS和三波段指数TBIs有效估算Holm oak叶片氮含量；Abdel-Rahman等[8]采用一阶导数的SR指数(波段743，1 316 nm)(波段743，1 317 nm)(波段741，1 323 nm)估算的叶片氮含量的决定系数在0.75左右；Yao等[9]采用归一化指数NDSI和比值指数RSI一阶导数对叶片氮的累积量进行了估算，其决定系数高达0.81；Ullah 等[10]利用MERIS数据采用植被指数(NDVI，NBDI，SAVI，TSAVI，REIP，MTCI)以及波段深度参数对草地生物量及氮含量进行了估算和对比。Sanches等[11]采用反射率，吸光度及其衍生形式和偏最小二乘方法对植物氮磷钾进行估算，结果表明一阶导数形式的估算精度最高。另有学者着重从算法角度进行了植物氮含量高光谱反演技术研究，Zhang等[12]利用单变量线性回归、逐步多元线性回归、偏最小二乘对湿地植物芦苇氮含量进行估算，并采用留一交叉验证比较3种模型的精度，认为氮的敏感波段集中在红波段和绿波段，单变量线性回归对665 nm和680 nm处的归一化植被指数效果最好，3种方法中逐步线性回归的精度最高。Yao等[13]采用原始反射率、一阶导数、植被指数等变量利用SMLR、PLSR、ANNS、SVMs等反演了冠层叶片氮含量，结果表明一阶导数的支持向量机方法精度更高。Wang 等[3]验证了多核支持向量机在估算小麦叶片氮含量的有效性，并与多元线性回归，偏最小二乘，人工神经网络，单核支持向量机进行了对比。总体来说，现有的关于高光谱植物氮含量监测的大部分研究集中在植被指数反演氮含量以及敏感波段筛选方面，关于建模算法的研究较少[3]。且已有的研究或集中于建模之前敏感特征的筛选或集中于建模算法的比较，本文尝试将不同特征提取方法与建模方法相结合，通过不同方法间的组合，筛选出一套完整的能有效估算叶片氮含量的流程。

1 材料和方法

1.1 光谱数据与氮含量数据获取

数据集采用Accelerated canopy chemistry program(ACCP)(https://doi.org/10.3334/ORNLDAAC)[14]，该数据集旨在研究不同生态系统中植被冠层氮含量及木质素含量遥感反演的理论基础，包含野外样品的实验室化学分析数据、光谱数据、小型冠层试验的化学分析数据及冠层建模数据。本研究选择ACCP中91组道格拉斯冷杉幼苗期新鲜叶片光谱数据及相应的实验室化学分析数据，采集于1992年11月。光谱数据为400～2 498 nm处的吸光度[log(1/反射率)]，光谱间隔为2 nm，分辨率为10 nm。在采集光谱数据后的6 h，采集叶片测量叶片氮及叶绿素含量。采用常规的实验室方法，对每棵树幼苗期的叶片样本进行干燥、研磨和全氮、全叶绿素(a+b)分析。样品经硫酸-氧化汞催化剂(Perstorp分析)在块状消化器中消化后，用Alpkem连续流自动分析仪测定总氮,氮含量的单位为单位干重百分比。

1.2 数据分析

利用Matlab对试验数据进行了异常数据剔除、计算了原始反射率、一阶导数R′、反射率倒数的对数log(1/R)及其一阶导数log(1/R)′，以及这些参数与氮含量的相关系数，采用连续投影算法、LASSO、相关系数筛选以上4种参数形式的特征波段，并采用随机森林、偏最小二乘以及支持向量机建立氮含量反演模型。其中，相关系数、连续投影算法筛选波段组合，偏最小二乘建模在Matlab中完成，LASSO算法筛选变量组合、随机森林和支持向量机建模在R中完成。

1.2.1 相关系数法 相关系数法是最简单直接的能够帮助理解特征与相应变量之间关系的方法，衡量的是变量之间的线性相关性，其值绝对值越大表示相关性越强。很多研究通过选择最大相关性的变量作为特征变量进行建模研究。λ1-λ2，λ1/λ2表示两变量间的差异性，其与因变量间的相关系数表示两变量组合对因变量的解释量，选取相关系数最大的两个波长即为特征波长。

1.2.2 连续投影算法 连续投影算法(SPA)是一种矢量空间共线性最小化的前向变量选择算法，能够有效从全波段中提取特征波段，消除原始光谱矩阵中的冗余信息，降低模型的复杂度，在波长选取中取得了较好的效果[15]。连续投影算法是一种前向循环选择算法，设样本集M和波段数K组成一个光谱矩阵XMXK，分别记Xk(0)和N为初始的迭代向量和需要提取的波段个数。从一个波长开始，每次循环计算它在未选入的波长上的投影，将投影向量最大的波长引入波长组合，直到循环N次，每一个新入选的波段，都与前一个线性关系最小。通过循环会得到N×K对波段组合，对每一对XK(0)和N所决定的组合分别建立多元回归模型，并用预测均方根误差RMSE来决定所建模型的优劣，最小的RMSE对应的Xk(0)为最佳的波段组合。

1.2.4 随机森林 随机森林是一种并行式集成学习法，以决策树为基础，并在决策树训练过程中引入随机属性选择。它通过自助法重采样技术，从原始训练样本集N中有放回地重复随机抽取K个样本生成新的训练样本集合，然后根据自助样本集生成K个分类树组成随机森林[17-18]。

1.2.5 支持向量机 支持向量机是一种基于统计学习理论的机器学习算法，它建立在VC(Vapnik-Chervonenkis Dimension)维理论和结构风险最小原理基础上，能较好地解决小样本、非线性、高维数和局部极小点等实际问题，同时能获得较好的泛化能力[19]。

1.2.6 偏最小二乘 偏最小二乘回归是一种新型的多变量回归分析方法，可以实现回归建模、简化数据结果和分析两组变量见得相关性，给多元数据统计分析带来极大便利，它能在回归建模过程中采用数据将为、信息综合和筛选技术，提取对系统最佳解释能力的新综合成分[20]。

1.3 建模方法与精度评价

分别采用偏最小二乘，随机森林与支持向量机建立叶片氮含量反演模型，并采用决定系数R2、均方根误差来评价各个模型的精度。R2表示相关密切程度，RMSE用来衡量估测值与真实值之间的偏差程度。R2越大，RMSE的绝对值越小，表示模拟结果精度更高。

2 结果与分析

2.1 不同光谱形式与氮含量的相关性分析

原始反射率R与氮含量呈负相关关系，在500～730 nm波段的相关系数超过-0.6，在550、702 nm左右相关性达到最大，且差异显著(P<0.05)，在550 nm左右相关性最大是由于叶绿素对绿光的反射作用，702 nm位于红边区域，有相关研究发现红边波段与植物叶绿素含量、生物量等参数存在显著相关性；氮含量与原始光谱的相关性在近红波段较小(表1)。

一阶导数R′能够有效消除叶面积变化的影响，反射率的一阶微分与氮含量的相关性变化较为复杂，在490～544、556～710和736～800 nm超过0.6，其中R′在530、574、626，694 nm达到最大，近红波段的相关性较小，但在1 200、1 672、2 164、2 310 nm相关性分别达到峰值0.6以上(图1)。

根据工程的复杂程度，对设计资质提出具体要求，如对提灌站、压力输水管道、100 m以上的深井、规模较大的建筑物等技术含量较高的工程必须聘请有相应资质的单位进行设计，从源头抓好工程质量。

图1 R，R′，log(1/R)，log(1/R)′与氮含量的相关性

Log(1/R)相关性曲线与原始反射率的相关性曲线线型一致，但相关性相反。log(1/R)′与氮含量的相关性变化也较为复杂，在490～546，620～670，672～694，702～804 nm相关性在0.6以上，在红外波段1 204，1 672，2 158 nm的相关性形成峰值，达到0.6以上。

应用λ1-λ2，λ1/λ2与氮含量的相关系数筛选出的波段组合如表1所示，R与log(1/R)筛选出的波段集中在可见光波段，两种导数形式筛选出的波段集中在近红波段；原始波段反射率及其变换形式筛选出的波段组合都包含522，526 nm；导数形式筛选出的波段组合均包含2 158，2 044和2 070 nm 3个波段，虽然这几个单独的波段与氮含量的相关性并不强，但其波段组合与氮含量的相关性达到0.86以上，说明该波段组合包含了氮含量的大部分信息。

表1 相关系数法筛选的波段组合

2.2 连续投影算法筛选敏感波段特征筛选

ACCP叶片数据集共91组鲜叶片数据，随机选取其中65组数据进行连续投影算法。连续投影方法将65组数据分为训练集和校正集，以校正集的预测均方根误差来确定最佳的光谱变量总数，均方根误差和决定系数来评价最佳的波段变量组合。对以上4种形式的光谱参数分别采用连续投影算法，确定氮含量反演的最佳光谱变量(表2)。

表2 连续投影法筛选的波段组合

原始反射率所选的波段在波段范围内分布较为均匀(图2)，位于可见光的波段有6个，2个波段位于红边区域，位于波峰或波谷处的波段有6个。一阶导数只选择了3个变量，其中一个位于可见光波段。Log(1/R)投影算法选择了两个红光波段，其余均为近红波段。log(1/R)′仅选择了3个红外波段。校正数据集RMSE最低的是Log(1/R)形式，RMSE最高的是一阶导数形式。

图2 连续投影算法筛选出的波段位置

2.3 LASSO算法

采用Lasso算法对4种光谱形式进行变量筛选，各光谱形式筛选出的变量数均比连续投影筛选出的变量数多。原始光谱反射率筛选出22个波段，其中有12个波段位于可见光区。Log(1/R)筛选出的变量数为16个，其中9个可见光波段。两种导数形式均筛选出23个变量数，一阶导数筛选出的可见光波段8个，Log(1/R)′筛选出的可见光波段6个(表3)。

表3 LASSO筛选的波段组合

2.4 模型预测结果与对比

分别将相关系数筛选法、连续投影算法、LASSO算法筛选出的波段作为偏最小二乘、随机森林、支持向量机的输入来估算叶片氮含量，以65组数据为建模数据，26组数据为验证数据。不同方法组合的反演结果如表4及图3所示。

表4 不同方法组合的反演结果

3种变量筛选方法和3种建模方法的组合中，RMSE在0.19～0.38，R2在0.6～0.89，其中，随机森林算法对相关系数法和LASS0算法筛选出的Log(1./R)′形式变量的反演误差最小，RMSE在0.19～0.20，能够解释氮含量的变化，但该算法对筛选出的其他形式的变量反演误差较大，表现并不稳定;不论采用何种建模方法，对于相关系数法，其log(1/R)′形式的变量反演效果较好；对于连续投影法，其筛选出的R′变量，反演误差较小；对于LASSO算法，其筛选出的R'和log(1/R)′形式的变量反演效果要明显优于其余两种变量形式。同时，分别采用偏最小二乘，随机森林、支持向量机进行全波段的建模，反演结果分别为R2=0.605、RMSE=0.4542，R2=0.659、RMSE=0.3352，R2=0.730 4、RMSE=0.348 4。表明，相关系数&Log(1./R)′，LASSO&Log(1./R)′，SPA&R′能够有效筛选出冷杉叶片氮含量的敏感波段组合，且前两种方法组合效果更优。

3 讨论

3.1 不同变量筛选方法和建模方法的组合效果比较

从图3可以看出两种导数形式估算出的氮含量与实测氮含量1:1的对应关系更好，散点基本上沿着y=x线分布，具有较好的线性关系。对于R及log(1/R)形式的反演，三种建模方法的结果与实测值相比普遍偏高，大部分点分布在y=x线以上。两种导数形式的反演结果较为均匀紧致地分布在y=x线两侧，反演结果与实测数据更接近。从实测数据与反演结果的拟合图可以看出，相关系数法&log(1/R)′以及LASSO&log(1/R)′中实测数据与反演结果的拟合函数具有较大的斜率，较小的截距，且R2较大，RMSE较小，具有较高的反演精度(图3)。

图3 实测氮含量与不同方法反演氮含量散点图

3.2 3种筛选方法的敏感波段分析

4 结论

采用相关系数法、连续投影法、LASSO算法对4种形式的光谱变量R，R′，log(1/R)，log(1/R)′进行敏感波段筛选，对筛选后的敏感波段分别采用偏最小二乘、随机森林、支持向量机对氮含量进行建模反演，得出以下结论：

1)两种导数形式变量的反演误差最小。

2)相关系数&Log(1/R)′，LASSO&Log(1/R)′能够有效筛选出冷杉叶片氮含量的敏感波段组合，无论采取何种建模方法，其估算效果在几种筛选方法组合中最好，R2>0.84，RMSE在0.19～0.24，估算效果明显优于采用3种建模方法进行全波段建模。

3)随机森林算法对相关系数法和LASS0算法筛选出的Log(1/R)′形式变量的反演误差最小，但对其他形式筛选出的变量反演结果的误差变化范围较大，表现并不稳定。

4)3种变量筛选方法筛选出的R形式的变量与前人研究相符，连续投影算法能筛选出更多与叶片其他化学含量相关的波段，证明了3种变量筛选方法的有效性。