基于Sentinel-2多光谱数据的枸杞叶面积指数反演

余刚冯克鹏

(1.宁夏大学土木与水利工程学院，宁夏 银川 750001；2.宁夏节水灌溉与水资源调控工程技术研究中心,宁夏 银川 750021；3.旱区现代农业水资源高效利用教育部工程研究中心，宁夏 银川 750021)

引言

叶面积指数(Leaf area index，LAI)是模拟陆地生态过程、生物地球化学循环以及水热循环的重要参数[1]，在农业、林业、生态、气候变化等领域已经得到了广泛的应用[2]。LAI能够反映植被长势个体和群体特征，是农作物长势监测的一个关键生态参数[3]，因此快速、准确、可靠地获取农作物LAI具有重要意义。

随着遥感技术的飞速发展，过去的几十年时间里，许多新发射的中分辨率对地观测卫星投入了研究应用，如1999年12月18日、2002年5月4日发射的MODIS卫星，2013年2月11日发射的Landsat-8卫星，2013年4月26日发射的高分1号卫星，2015年6月23日、2017年3月7日由欧空局发射的Sentinel-2A和Sentinel-2B多光谱遥感卫星。Sentinel-2多光谱遥感卫星搭载包括空间分辨率为10m的3个可见光波段以及1个近红外波段，空间分辨率为20m的特有3个红边参数波段和短波红外波段，空间分辨率为60m的近红外波段等，共计13个波段。Sentinel-2A和Sentinel-2B双轨卫星同时运行，时间分辨率缩短至5d，为遥感卫星数据在农业遥感中的应用提供了更多的可能性[4-7]。

Sentinel-2卫星所搭载的红边参数波段在农业遥感中具有很高的有效性[8]。Frampton等[9]基于模拟的Sentinel-2数据提出了用于估算作物叶绿素和叶面积指数LAI的IRECI(inverted red edge chlorophyll index)和 S2REP(Sentinel-2 red-edge position)红边参数植被指数，Korhonen等[10]开展的森林叶面积指数估算，以及苏伟等[11]利用Sentinel-2数据进行的玉米叶面积指数估算，这些研究基于Sentinel-2遥感影像所构建的的各种估算模型，都证实了Sentinel-2多光谱数据用于LAI估算的可行性，但研究中所构建的LAI估算模型均是针对特定研究区域和研究对象，无法直接用于宁夏地区枸杞LAI估算。同时在这些研究中，所构建的LAI估算模型大多是物理模型和经验模型，无法直接应用在新的数据源上。

机器学习算法通过相关数据拟合一个灵活的模型，学习输入变量与输出变量之间的关系，将模型参数调整为在一个独立验证集中预测误差最小，使模型获得最好的泛化能力，其可以同时输入多个与LAI相关性高的植被指数，并具有很好的非线性拟合能力，可提高模型的稳定性及估测精度[12]，是农作物LAI估算的重要手段。近年来，已有众多学者使用机器学习算法进行LAI遥感估算[13-15]。梁栋等[16]基于支持向量机回归算法，将不同植被指数与波段反射率作为输入参数，在冬小麦多个生育期反演LAI具有较好的适用性。

本研究通过大面积田间实验实测LAI，构建基于Sentinel-2多光谱遥感影像提取的波段反射率和计算的植被指数，构建3种不同数据输入模型，运用机器学习算法进行估算，对比分析反演精度，研究结果将进一步丰富枸杞LAI估算方法，提高估算精度。

1 材料与方法

1.1 研究区概况

研究区位于宁夏回族自治区，地理坐标为E104°17′～107°39′，N35°14′～39°23′。海拔1100～1200m，平均气温-9℃，年降水量在150～600mm，属中温带干旱气候带，干旱少雨，光照条件充足，蒸发强烈。

图1 研究区位置

1.2 采样区设置

研究区内枸杞种植广泛，本研究选取包括2年期、4年期和多年期枸杞植株种植园，共6个采样区，分布于宁夏中北部。枸杞采取沟垄相间的种植方式，在每个采样区内选取面积较大且相对均质的枸杞地块，每个地块随机设置1个10m×10m样方。在全部采样区共设置了80个样方，并采用GPS记录样方经纬度坐标。采区设置及枸杞植株物理量如表1所示。

表1 样区设置及枸杞植株物理量

1.3 枸杞LAI测定

地面数据采集于2019年8月15—18日，使用LAI-2200C冠层分析仪测量样方内每行枸杞LAI值，取样方内所有行平均值为该样方LAI测量结果。LAI-2200C植物冠层分析仪是利用“鱼眼”光学传感器测量植被冠层5个不同天顶角方向的散射天空辐射衰减，通过5个天顶角方位的透光度结果计算出叶面积指数、空隙比等树冠结构参数。LAI测量时应尽量避免在直射光背景下测量，测量模式采用AAAABBBB。第1个A值采用白色散射帽，测量人员背对太阳，使阳光直射散射帽，保持探杆水平测量；第2个A值将探杆移动至测量人员头部阴影下测量；第3个A值将白色散射帽移除，保持在测量者头部阴影下测量；第4个A值使用45°散射帽，正常测量。B值在样方内部采用两垄之间对角线测量，每条对角线上采集4个B值，第1个B值取在垄上，第2个B值取在两垄之间1/4处，第3个B值取在两垄中间，第4个B值取离垄3/4处。实测LAI统计如表2所示。

表2 样区实测LAI统计

1.4 Sentinel-2遥感影像

遥感影像使用Sentinel-2多光谱卫星影像，通过欧洲航空局的数据共享网站(https://scihub. copernicus.eu/dhus/#/home)下载，获取与地面观测时间同步的影像数据，空间分辨率10m，成像时间分别为2019年8月15日和18日。下载的遥感影像均为已经过辐射校正和几何校正处理的Level-1C级大气上层表观反射率。遥感影像大气校正采用SNAP-Sen2Cor软件处理，经大气校正后的遥感影像通过最近邻插值法重采样至10m后用于提取单波段反射率和计算植被指数。Sentinel-2多光谱波段信息见表3。

表3 Sentinel-2多光谱传感器参数

1.5 植被指数

利用Sentinel-2多光谱卫星可见光及近红外波段计算多种植被指数，主要包括土壤校正型植被指数、红边参数植被指数、叶绿素含量植被指数以及传统的近红外植被指数，共计21种植被指数[17-21]。各植被指数基于Sentinel-2数据各波段计算公式如表4所示。本研究所采用的植被指数使用欧洲航空局SNAP软件计算得到。

表4 植被指数及基于Sentinel-2数据的计算方法

1.6 机器学习算法

多元线性回归模型(Multivariable Linear Regression，MLR)通常用来描述变量y和x之间的随机线性关系，用来研究一个应变量依赖多个自变量的变化关系，特点是形式简单、易于建模[33]。

逐步多元线性回归模型(Stepwise Multiple Linear Regression，SMLR)是将变量逐个引入模型，每引入一个解释变量后都要进行F检验，并对已经选入的解释变量逐个进行t检验。这是一个反复的过程，直到既没有显著的解释变量选入回归方程，也没有不显著的解释变量从回归方程中剔除为止，以保证最后所得到的解释变量集是最优、最简单的[34]。

支持向量机(Support Vector Machine，SVM)是一种二分类模型，其目的是寻找一个超平面对样本进行分割，分割的原则是间隔最大化，最终转换为一个凸二次规划问题来求解[35]。

高斯过程回归模型(Gaussian Process Regression，GPR)是使用高斯过程先验对数据进行回归分析的非参数模型[36]。GPR有着严格的统计学习理论基础,对处理高维数、小样木、非线性等复杂的问题具有很好的适应性，且泛化能力强[37,38]。

树的集成模型(Ensembles of Trees)是使用基本模型创建组合模型的学习方法。主要包括梯度提升树和随机森林2种，这2个算法都使用决策树作为基础模型，两者都是学习集成树的算法，但是训练过程不同[31]。

回归树模型(Regression Trees)由Breiman等人在1984年提出，是应用广泛的决策树学习方法。由特征选择、树的生成及剪枝组成，既可以用于分类也可以用于回归。回归树算法流程类似于分类树，分枝时穷举每一个特征的每一个阈值，来寻找最优切分特征j和最优切分点s，衡量的方法是平方误差最小化，分枝达到预设的终止条件(如叶子个数上限)就停止[32]。

1.7 评价指标

选用决定系数(R2)、均方根误差(RMSE)进行分析检验，其中R2表示模拟值与实测值的拟合优度，其值越接近1，表示模型的参考价值越高；RMSE主要用于模型验证，反映了模拟值与真实值的偏离程度，其值越小，模型精度越高。指标的计算公式：

2 结果与分析

2.1 遥感影像反射率特征

对Sentinel-2多光谱遥感影像各波段反射率在不同LAI值时的变化特征进行分析，LAI值选取本研究中实测LAI最大值(LAI=2.62)、最小值(LAI=0.22)和平均值(LAI=1.01)的样方所对应的Sentienl-2多光谱反射率进行对比，如图2所示。3种不同LAI值情况下，光谱曲线均在400～700nm的可见光区域冠层反射率均呈现先上升再下降趋势，在蓝光波长450nm和黄光波长650nm左右有1个吸收谷，在绿光波长550nm左右有1个反射峰，这是因为植物叶片中叶绿素对蓝光和红光辐射吸收较强，对绿光辐射反射较强。在700～850nm的近红外区间内，反射率随着波长增大而增大，此谱段急剧上升，具有陡而近于直线的形态，在红边波段850nm左右出现最大值，形成“红边”现象。主要原因是由叶片细胞结构造成的，在这一区域，叶片结构导致对光谱的反射和透射达到最大，吸收最小，所以出现反射峰。随着LAI值增大，反射率数值上开始逐渐增大，光谱曲线的位置增高，这是由于随着叶片大小、数量和层数的增加，光线在叶片的不同层面发生多次反射，光线除了在上层叶片形成反射以外，还可以透过上层叶片在下层叶片发生二次反射，导致植株总体反射率急剧上升。在850～2200nm区间，反射率随着波长增大而减小，这是因为在这一波长范围内，光线被植物叶片中的水分和二氧化碳吸收，导致反射率减小。

图2 不同LAI值反射率变化规律图

2.2 波段反射率及植被指数与LAI相关关系分析

Sentinel-2多光谱遥感影像的波段反射率及各类植被指数与枸杞LAI的皮尔逊相关性分析见表5。

表5 波段反射率及植被指数与实测LAI相关关系

对比各单波段反射率与LAI的相关性，可见光波段B1～B4和B11、B12波段与LAI呈负相关；红边参数波段B5～B7以及近红外波段B8和B8a与LAI均在(P<0.01)水平呈极显著相关且相关系数均大于0.6，其中近红外波段B8相关性表现最佳，这与已有研究证实的LAI反演中红边参数波段及近红外波段为敏感波段，以及LAI与这些波段的反射率呈极显著线性相关关系的结论一致[39]。对比各类植被指数的表现发现，本研究所选植被指数与LAI均在(P<0.01)水平呈极显著相关，其中红边参数植被指数IRECI与LAI相关性最好，具有最大的相关系数。

2.3 LAI估算模型构建

随机将样方分为训练集和测试集，其中训练集60组，测试集20组。基于Sentinel-2多光谱遥感影像提取的波段反射率和计算的植被指数，本研究共建立3种数据输入模式。模式1以波段反射率为自变量，对应实测LAI值为因变量；模式2以植被指数为自变量，对应实测LAI值为因变量；模式3以波段反射率和植被指数为自变量，对应实测LAI值为因变量。将训练集数据驱动MLR、SML、GPR、SVM、Rtree和Etree算法，进行多次计算，分别统计各种算法决定系数(R2)最高的10次计算结果，并求算术平均数。结果见表6和图3。

从表6、图3可知，训练集数据中，在3种数据输入模式中，GPR算法均表现了较强的预测能力，决定系数(R2)分别为0.803、0.668、0.699，均高于同模型中其它算法；对应均方根误差(RMSE)为0.402、0.453、0.441，均低于同模型中其它算法。

图3 训练集中各输入模式不同算法R2和RMSE箱线图

表6 训练集中各输入模式不同算法R2和RMSE统计

将训练好的模型带入测试集中，统计3种数据输入模式的决定系数(R2)和均方根误差(RMSE)，结果如表7，箱线图如4所示。可知基于最优的GPR算法，3种数据输入模式中，基于波段反射率的模式1具有最好的反演精度，决定系数(R2)高于其它模式，均方根误差(RMSE)低于其它模式。

表7 测试集中GPR算法各模式R2和RMSE统计

图4 测试集中各输入模式R2和RMSE箱线图

3 结论

输入模式的选择对农作物生理参数遥感反演的精度影响较大，各种机器学习算法预测精度差距较大，建立不同数据输入模式和定量评价不同算法反演精度对农作物生长遥感监测意义重大。本研究利用Sentinel-2多光谱遥感卫星数据及枸杞实测LAI数据，探讨了Sentinel-2多光谱卫星数据用于枸杞LAI估算的特点和精度。主要得到以下结论。

分析基于Sentinel-2多光谱数据提取的波段反射率和实测LAI值，可以得出在400～700nm可见光区域冠层反射率呈现先上升再下降趋势，在700～850nm的近红外区间内，反射率随着波长增大而增大，在红边波段850nm左右出现最大值，形成“红边”现象，在850～2200nm区间，反射率随着波长增大而减小。

基于Sentinel-2多光谱卫星数据提取的波段反射率和以此计算的植被指数，可见光波段B1～B4和B9～B12波段与LAI呈负相关；红边参数波段B5～B7以及近红外波段B8和B8a与LAI呈极显著相关，其中近红外波段B8相关性表现最佳。对比各类植被指数的表现发现，植被指数与LAI均呈极显著相关，其中红边参数植被指数IRECI与LAI相关性最好，具有最大的相关系数。

基于Sentinel-2多光谱卫星数据提取的波段反射率和以此计算的植被指数，在划分的训练集中，3种数据输入模式下，GPR算法均表现较强的预测能力，具有最高的决定系数(R2)和最低的均方根误差(RMSE)。将训练好的模型带入测试集中，以波段反射率为数据输入模式具有最好的反演精度。