基于哨兵-2A模拟反射率及其影像的冬小麦收获指数估算

任建强 张宁丹 刘杏认 吴尚蓉

(1.中国农业科学院农业资源与农业区划研究所， 北京 100081；2.中国农业科学院农业环境与可持续发展研究所， 北京 100081)

0 引言

收获指数(Harvest index，HI)是评价作物单产水平和栽培成效的重要生物学参数，也是作物单产进一步提高的重要决定因素之一[1]。对粮食作物来说，一般的收获指数是指成熟期作物籽粒产量占作物地上生物量的百分数，该指标本质反映了作物同化产物在籽粒和营养器官中的分配比例[2]。近年来，为了定量描述籽粒灌浆过程中作物收获指数逐步形成至达到最大值的变化过程，部分学者以冬小麦为例提出了作物动态收获指数(Dynamic harvest index，D-HI)的概念，这不仅丰富了收获指数的概念内涵，也对作物生长过程和产量形成的定量模拟具有重要意义[3]。因此，快速准确地获取成熟期作物收获指数及其动态过程信息对于作物单产准确模拟[4-6]、作物生物量估算[7]、作物品种选育和表型信息获取[8-9]、作物栽培技术优化与效果评价[10-11]具有重要科学意义，同时，对农业管理部门及时掌握农作物长势和作物产量估算信息，有效开展农业生产管理也具有重要指导意义[12-14]。

一般的作物收获指数研究主要基于田块尺度的农学试验层面，大多侧重于作物收获指数的数学模拟、收获指数与相关农学参数关系、作物生长环境及其管理措施对收获指数的影响评价等[15-16]。作物收获指数获取方法主要有直接法和间接法。直接法主要通过田间取样计算获得，该方法虽然准确，但费时费力且很难在大范围内开展，也无法获得收获指数的连续空间分布信息。间接法主要包括2类：①通过模拟灌浆期作物收获指数的逐步增加变化过程与时间之间函数关系实现作物收获指数的准确估算[17-18]。②基于作物生长过程中的植被信息(如各种植被指数、生物量等)、作物生长环境影响因子(如温度、光照、土壤含水率等)与HI建立函数关系实现HI的估算[19]。

近些年来，遥感技术凭借其覆盖范围大、快速和准确获取地表作物参数信息的优势，为准确获取区域作物HI空间信息提供了可靠的技术手段[20-21]。其中，通过宽波段多光谱传感器(如NOAA、MODIS、MERIS等)数据，国内外学者利用发芽-开花和开花-成熟2个阶段能够反映作物长势状况的时序植被遥感信息(如归一化差值植被指数和叶面积指数等)开展了成熟期冬小麦HI遥感估算研究，对利用遥感信息获取区域尺度作物HI具有重要借鉴意义[15,17,22-24]。然而，已有研究只进行了成熟期最终作物收获指数估算，且使用开花后时序NDVI表征籽粒灌浆过程，对作物HI形成过程中动态收获指数信息的考虑不足[25]。针对这一问题，张宁丹等[26]以河北省衡水市深州市冬小麦为研究对象，在构建花后累积地上生物量比例动态参数(Dynamic fG, D-fG)基础上，提出了敏感波段中心构建归一化差值光谱指数(Normalized difference spectral index，NDSI)估算D-fG的作物动态收获指数估测方法，实现了田间冠层尺度D-HI的准确获取。在敏感波段中心基础上，通过波段扩展确定了冬小麦D-fG估算敏感波段最大宽度，实现了最大波宽下田间冠层尺度D-HI的遥感准确获取。上述冠层尺度D-HI估算研究为基于高光谱和多光谱等遥感卫星开展区域尺度作物动态收获指数获取奠定了一定基础。但是，该方法由田间冠层尺度扩展到区域尺度的应用效果还有待进一步评价，在已获得的敏感波段中心和最大波宽基础上如何开展基于遥感卫星的区域冬小麦收获指数遥感估算还有待进一步深入研究，特别是如何充分利用反映作物生长状况的红光、红边、近红外等波段信息进行大范围D-HI遥感估算成为亟待研究的问题。

基于以上分析，考虑到哨兵-2A(Sentinel-2A)是具有较高重访频率的宽幅高分辨率多光谱成像卫星，是目前唯一在红边区域含有3个独特红边波段的光学遥感数据，这对监测植被健康信息和植物生长状况非常有效[27-29]，而且该数据已在大范围作物监测中发挥了重要作用[30-33]。但是，基于Sentinel-2数据的作物HI估测研究尚未见报道。因此，为进一步验证基于花后累积地上生物量比例的D-HI遥感估算方法在大范围获取作物HI空间信息的可行性和有效性，在前期田间冠层尺度作物HI估算研究确定的冬小麦D-fG和HI估算敏感波段中心和最大波段宽度基础上，本研究利用近地高光谱和Sentinel-2A影像数据，开展基于Sentinel-2A遥感模拟数据和真实遥感数据的区域冬小麦HI估算研究，以期为大范围作物HI空间信息卫星遥感准确获取提供技术方法借鉴。

1 研究区与数据

1.1 研究区概况

图1 研究区位置与采样点空间分布Fig.1 Location of study area and spatial distribution of sampling points

研究区为中国北方粮食生产基地黄淮海平原内河北省衡水市(115.17°～116.57°E，37.05°～38.38°N)，典型试验区为深州市(115.35°～115.80°E，37.70°～38.17°N)(图1)。衡水市总面积8.836×107km2，地处河北冲积平原，地势自西南向东北缓慢倾斜，海拔12～30 m。该研究区域属于暖温带半干旱区季风气候，年均降水量约480 mm，年均温度约13.4℃，无霜期200 d左右，研究区为冬小麦-夏玉米一年两熟轮作种植制度。其中，冬小麦种植时间为10月上中旬，返青期为次年3月上中旬，拔节期为3月下旬至4月上中旬，抽穗开花期为4月下旬至5月上旬，灌浆-乳熟期为5月中下旬，成熟期为6月上旬。

1.2 数据获取与处理

1.2.1地面数据采集与处理

数据获取内容主要包括冠层高光谱、作物地上干生物量、动态籽粒产量和GPS定位信息等。本研究分别在冬小麦开花期(2021年5月3日)、灌浆前期(2021年5月15日)、灌浆后期(2021年5月25日)、成熟期(2021年6月5日)共4个关键时期进行数据采集。深州市典型研究区共18个样方，每个样方布设5个采样点，每次试验获取90个样本数据。衡水市区域调查共布设50个调查点，每次调查获取50个样点数据，主要包括作物地上生物量、动态籽粒产量和GPS定位信息等。以冬小麦开花期为时间基准，在深州市共获取了灌浆前期、灌浆后期和成熟期共3个时期270个地面样本数据，根据比例3∶2 将其分为建模数据集(162个)和验证数据集(108个)。衡水市10个县50个调查点在灌浆前期、灌浆后期和成熟期3个时期共获得150个地面样本数据，主要作为后期衡水地区作物收获指数的区域验证数据。

(1)地上干生物量

在冬小麦关键生育期(包括开花期、灌浆前期、灌浆后期和成熟期)，首先根据各个采样点测定的GPS位置信息，在每个采样点分别割取1行长度20 cm的冬小麦地上部分作为样本。在将冬小麦茎叶穗分离基础上，将分离后冬小麦茎叶穗放入105℃干燥箱并进行30 min杀青处理。然后，在85℃条件下连续干燥48 h以上，直至样本质量恒定再进行称量。其次，在将冬小麦茎叶穗干质量进行相加的基础上，根据种植密度和样本干质量换算成单位面积冬小麦地上干生物量。最后，在对各样点冬小麦穗脱粒处理基础上，分别记录各个采样点的籽粒质量。

(2)冬小麦冠层高光谱

利用美国ASD(Analytical Spectral Devices Inc.)公司生产的Field Spec 4型光谱辐射仪(波长350～2 500 nm)采集冬小麦冠层高光谱。其中，波长范围350～1 000 nm采样间隔1.4 nm，波长范围1 000～2 500 nm采样间隔2 nm，重采样后光谱间隔为1 nm。选择天气状况良好、阳光照射充足条件下进行冠层高光谱测定，观测时间为10:00—14:00。高光谱测量前用标准白板校正，测量时保证探头垂直向下，光谱装置探头视角为25°，为了保证冬小麦样本处于探测视场内，同时减少下垫面光谱反射对测定结果的影响，探头距离作物冠层顶部高度约0.5 m。测量时，每个采样点获取10条光谱数据，取其均值作为该采样点的光谱反射率。然后，对每个采样点的光谱曲线利用9点加权移动平均法进行去噪平滑处理。图2为深州市调查样方采样点经过光谱平均和光谱平滑处理后的不同生育期冬小麦冠层高光谱曲线。

图2 深州市典型样方冬小麦不同生育期冠层 高光谱曲线Fig.2 Canopy hyperspectral curves of winter wheat at different growth stages of quadrats in Shenzhou City

1.2.2遥感数据获取与预处理

Sentinel-2A数据主要从欧空局哥白尼数据中心(https:∥scihub.copernicus.eu/)获取，卫星数据产品级别为Level-1C(L1C)。Sentinel-2A遥感影像有13个光谱波段，包括可见光、红边、近红外、水汽、卷云以及短波红外波段等。其中，L1C级数据为已经过辐射定标和几何校正的大气顶反射率数据，因此，只需要进行大气校正处理便可得到地表反射率数据。利用Sen2cor、SNAP(Sentinel application platform)和ENVI软件对卫星影像分别进行大气校正、格式转换和镶嵌裁剪等预处理工作。为了保证遥感数据各波段空间分辨率的一致性，在SNAP平台中使用最邻近内插法将遥感数据各波段重新采样至10 m。结合冬小麦物候信息以及地面观测试验时间，根据选择影像云量较低原则(小于10%)，获得了衡水市2021年5月17日、2021年5月27日、2021年6月6日共3个时期Sentinel-2A影像数据。考虑到前期冠层高光谱的作物收获指数估算敏感波段中心和最大波段宽度仅对应Sentinel-2A的蓝光、绿光、红光、红边-1、红边-2、红边-3、近红外等8个波段，因此，仅使用了相应的B2波段至B8a波段信息，具体Sentinel-2A多光谱影像波段信息如表1所示[34-36]。

1.2.3其他辅助数据

(1)光谱响应函数

由于ASD冠层高光谱与Sentinel-2A卫星影像数据的波段宽度不同，为了将前期冠层尺度ASD高光谱确定的冬小麦D-fG估算敏感波段中心以及最大波段宽度应用到宽波段多光谱Sentinel-2A卫星影像数据，需要利用卫星的光谱响应函数将近地面高光谱数据转换为卫星波段反射率，从而将近地面高光谱数据与Sentinel-2A卫星遥感各个波段反射率对应。Sentinel-2A光谱响应函数如图3所示。

表1 本研究选取的Sentinel-2A卫星遥感主要 波段参数Tab.1 Main band information of Sentinel-2A selected in this study

图3 Sentinel-2A光谱响应函数Fig.3 Curves of spectral response function of Sentinel-2A

(2)冬小麦空间分布信息

为了获得研究区的冬小麦空间分布信息，采用被广泛使用的监督分类中支持向量机(Support vector machine，SVM)方法进行研究区2021年冬小麦空间分布提取，利用的遥感影像为Sentinel-2A影像。支持向量机是一种建立在统计学基础上的机器学习方法[37]，即在有限的分类样本信息前提下，在模型的复杂性和学习能力之间寻求最佳结果，保证得到的极值解是全局的最优解[38]。

样点数据主要包括训练样点数据和验证样点数据。其中，训练样点用于支持向量机分类进行作物分布提取，验证样点用于对冬小麦空间分布结果的精度验证。根据冬小麦实际调查点和Google Earth上选取的冬小麦和非冬小麦样点，在衡水研究区最终获得冬小麦样点和非冬小麦样点共计4 845个。其中，选择冬小麦样点870个、非冬小麦样点1 175个作为训练样本，其余1 300个冬小麦样点和1 500个非冬小麦样点作为验证样本对冬小麦空间分布结果进行精度验证，具体地面样点数据分布信息如图4所示。最终，本研究基于2021年4月19日的Sentinel-2A影像和SVM分类方法实现冬小麦空间分布提取。通过验证，冬小麦空间分布提取总体精度为91.77%，Kappa系数为0.835 6，精度达到较高水平，可以满足本研究区域冬小麦收获指数遥感估算所需作物空间分布的精度要求，衡水市冬小麦空间分布结果如图5所示。

图4 衡水市2021年冬小麦提取地面样点数据分布Fig.4 Distribution of ground sampling points for winter wheat mapping in Hengshui City in 2021

图5 衡水市2021年冬小麦空间分布结果Fig.5 Spatial distribution results of winter wheat in Hengshui City in 2021

2 研究方法

2.1 总体技术路线

在筛选出的敏感波段中心及其最大波段宽度基础上，首先，确定对应Sentinel-2A遥感数据的波段范围，并利用遥感数据的光谱响应函数和地面高光谱数据进行宽波段遥感反射率模拟；其次，根据模拟的宽波段反射率数据构建归一化差值光谱指数(NDSI)进行D-fG和D-HI的估算和精度验证。在此基础上，筛选出冬小麦D-fG和D-HI估算最优波段组合。最后，在最优波段组合信息筛选结果下，基于Sentinel-2A真实卫星遥感数据进行区域冬小麦D-fG和D-HI的估算与验证。主要技术路线如图6所示。

图6 技术路线Fig.6 Flowchart of the research

2.2 主要参数构建与计算

2.2.1D-fG参数构建

为了估算作物收获指数，KEMANIAN等[39]提出了fG参数，该参数为作物开花期-成熟期累积地上生物量与成熟期地上生物量的比值，但一般fG参数只应用于成熟期fG计算，而未考虑开花期至成熟期之间fG参数的动态变化过程，这在一定程度上会降低利用fG参数进行作物收获指数估算的精度。为了提高作物收获指数估算模型的精度，增强模型的稳定性和可靠性，降低成熟期fG参数不稳定对收获指数估算精度的不利影响，在已有研究基础上，考虑了开花期-成熟期fG参数的动态变化过程，采用动态参数D-fG进行作物动态收获指数估算，即作物开花期至采样时期累积的地上生物量与对应采样时期地上生物量间比值[26]。D-fG指标计算公式为

(1)

式中DfG——花后累积地上生物量比例动态参数(D-fG)指标

Wpost——冬小麦开花期至采样日期累积地上生物量，kg/hm2

Wwhole——冬小麦播种至采样日期累积地上生物量，kg/hm2

Wt——播种至采样日期t累积地上干物质量，kg/hm2

Wa——开花期地上干物质量，kg/hm2

2.2.2动态收获指数构建

一般粮食作物(如小麦、玉米、水稻等)开花至成熟期间，随着籽粒不断灌浆，作物籽粒产量占作物地上部干物质量百分比呈逐步增加状态，直到成熟期收获指数达到最大值，其中，用于定量描述作物灌浆过程中籽粒产量占作物地上部干物质量百分比逐步增加至最大值的收获指数变化过程指标，称为动态收获指数[3]。因此，本文在已有研究基础上，利用冬小麦籽粒灌浆过程中各个地面观测时间的收获指数动态变化信息获得动态收获指数信息，即在获得各采样点冬小麦茎叶穗干质量基础上，分别对各采样点小麦穗进行脱粒处理，并记录各采样点的籽粒质量。最后，计算冬小麦灌浆至成熟期期间各个地面观测时间的冬小麦动态收获指数(D-HI)。D-HI计算公式为

(2)

式中DHI——作物动态收获指数(D-HI)

WZ,t——灌浆至成熟期期间采样日期t的冬小麦籽粒干质量，kg/hm2

WA,t——采样日期t的冬小麦地上干生物量，kg/hm2

2.2.3NDSI计算

光谱指数(Spectral index，SI)可以通过某些特定波段的组合来指示绿色植被内部的色素含量、水分变化和营养状态等[40]。为更好地利用多光谱遥感各个波长所包含的信息，本文将冠层高光谱模拟的多光谱遥感反射率及其真实遥感反射率分别进行任意两两波段组合，从而构建归一化差值光谱指数(Normalized difference spectral index，NDSI)，构建形式为

(3)

式中NDSI(λ1,λ2)——波长λ1、λ2计算的NDSI指数，值域为[-1,1]

Rλ1——波长λ1所对应光谱反射率

Rλ2——波长λ2所对应光谱反射率

考虑到作物光谱在波段1 350～1 415 nm和波段1 800～1 950 nm受大气和水蒸气影响较大[41]，因此，在可见光-近红外波段350～1 000 nm内进行D-fG估算敏感波段筛选和动态作物收获指数遥感估算，即式(3)中λ1、λ2表示在350～1 000 nm内的任意波长λ1和波长λ2。

2.3 地面高光谱模拟宽波段遥感数据方法

为了利用模拟遥感数据筛选作物D-fG和D-HI估算的最优波段，需利用灌浆前期、灌浆后期和成熟期的模拟Sentinel-2A多光谱数据开展冬小麦关键生育期动态收获指数估算。在已有近地面冠层高光谱确定的冬小麦D-fG和D-HI估算敏感波段中心和最大波段宽度研究结果基础上[26]，在对应最大波段宽度范围内[26](表2)，借助Sentinel-2A光谱响应函数，将实测的近地面高光谱窄波段反射率转换为Sentinel-2A多光谱波段反射率模拟数据，波段反射率转换公式为

(4)

式中Rrs——模拟卫星波段反射率

λm——传感器光谱探测起始波长

λn——传感器光谱探测终止波长

S(λ)——传感器在波长λ处的光谱响应函数值

R(λ)——冠层光谱在波长λ处光谱反射率

表2 D-fG和D-HI估算敏感波段中心最大波宽对应的 地面高光谱波段范围Tab.2 Wavelength range of ground hyperspectral bands corresponding to center of sensitive band and its maximum bandwidth for estimation of D-fG and D-HI

2.4 归一化光谱指数NDSI与D-fG间模型构建

基于模拟Sentinel-2A反射率数据和真实Sentinel-2A卫星反射率数据构建的归一化光谱指数NDSI分别进行D-fG遥感估算，为冬小麦动态收获指数获取奠定基础。NDSI与D-fG间的线性模型为

DfG=aNDSI(λ1,λ2)+b

(5)

式中a——一次项系数b——常数项

当利用模拟的宽波段反射率数据进行D-fG估算时，λ1、λ2分别为模拟多光谱遥感数据的波长，NDSI为模拟多光谱遥感波段λ1、λ2构建的NDSI光谱指数。当基于多光谱遥感卫星数据进行区域D-fG估算时，NDSI为真实宽波段遥感数据构建的光谱指数。

2.5 冬小麦动态收获指数估算模型

在KEMANIAN等[39]提出的基于成熟期实测fG的作物收获指数估算方法基础上，提出了基于D-fG遥感信息的D-HI遥感估算方法。D-fG和动态收获指数D-HI间统计关系模型为

DHI=HI0+kDfG

(6)

式中HI0——截距，即在作物开花期之后生物量不发生变化情况下的动态收获指数

k——斜率常数

2.6 精度评价

精度评价指标包括决定系数(Coefficient of determination，R2)、均方根误差(Root mean square error，RMSE)、归一化均方根误差(Normalized root mean square error，NRMSE)和平均相对误差(Mean relative error，MRE)。其中，基于模拟Sentinel-2A的冬小麦收获指数估算中，主要使用深州市灌浆前期、灌浆后期和成熟期共3个时期270个地面样本数据，将其分为建模数据集(n=162)和验证数据集(n=108)。基于真实Sentinel-2A卫星遥感波段的冬小麦收获指数的估算中，共获取2021年5月17日、2021年5月27日、2021年6月6日3个时期遥感影像，综合使用深州市270个地面样本数据和衡水区域调查的150个验证数据，最终得到建模数据集(n=162)和验证数据集(n=258)。

3 结果与分析

3.1 遥感数据间波段对应关系及NDSI计算

3.1.1Sentinel-2A模拟数据波段范围确定

在综合考虑敏感波段中心最大宽度对应的地面高光谱波段范围和Sentinel-2A真实遥感数据波段范围基础上，进一步确定了每个敏感波段中心最大波段宽度所对应的模拟Sentinel-2A波段范围。由于Sentinel-2A卫星遥感宽波段数据的波段宽度比地面高光谱波段数据宽，因此，构建NDSI的两个敏感波段在同一个宽波段内的可能性比较大，筛选出的波段组合相对较少。考虑到筛选出的地面高光谱波段范围与真实Sentinel-2A波段范围并非完全一致，为了便于利用光谱响应函数模拟Sentinel-2A波段反射率，当地面高光谱波段范围超出真实Sentinel-2A波段范围时，应用Sentinel-2A波段范围内的数据；当地面高光谱波段在Sentinel-2A波段范围内时，应用地面高光谱波段内数据[42]。地面高光谱波段范围与Sentinel-2A卫星遥感及其模拟数据之间波段对应关系如表3所示。其中，当D-fG敏感波段中心在(366 nm，489 nm)时，对应的地面高光谱波段范围(351～381 nm)不在Sentinel-2A真实波段范围内，因此，无法构建相应的NDSI；同理，当D-fG敏感波段中心在(443 nm，495 nm)时，地面高光谱波段(409～477 nm和461～529 nm)所对应的Sentinel-2A真实波段均处于蓝光波段内，也无法构建相应的NDSI，故本研究对上述两种情况进行了舍弃处理。

表3 地面高光谱波段范围与Sentinel-2A卫星遥感及其模拟数据之间对应关系Tab.3 Correspondence between wavelength ranges of the ground hyperspectrum and bands of Sentinel-2A satellite remote sensing images and their simulated date

3.1.2NDSI计算

根据表3中高光谱敏感波段对应的模拟Sentinel-2A波段数据和真实Sentinel-2A波段数据，计算了冬小麦灌浆前期、灌浆后期、成熟期Sentinel-2A模拟数据及其真实遥感影像各个波段组合构建的NDSI。限于篇幅，本文仅展示地面高光谱波段λ1(672～680 nm)和λ2(855～875 nm)对应Sentinel-2A红光波段和窄近红外波段构建的NDSI的空间分布图，冬小麦灌浆前期、灌浆后期和成熟期NDSI空间信息计算结果如图7所示。通过分析可知，冬小麦灌浆前期NDSI范围在0.40～0.95之间，研究区冬小麦NDSI平均值为0.88；冬小麦灌浆后期NDSI范围在0.20～0.80之间，研究区冬小麦NDSI平均值为0.58；冬小麦成熟期NDSI范围在0.10～0.50之间，研究区冬小麦NDSI平均值为0.35。

图7 基于Sentinel-2A卫星红光和窄近红外波段构建的冬小麦NDSI空间分布图Fig.7 NDSI spatial distribution map constructed by red band and NIR2 band of Sentinel-2A

3.2 动态收获指数估算模型建立

根据深州市建模数据集(n=162)中的地上生物量数据和灌浆过程中籽粒产量动态数据，计算162个冬小麦样本点的D-fG和动态收获指数D-HI。在此基础上，对实测D-fG和实测动态收获指数D-HI间的线性关系进行拟合，得到D-fG参数和动态收获指数D-HI间拟合方程

DHI=0.104 3+0.763 1DfG

(7)

其中，基于实测D-fG的动态收获指数D-HI线性估算模型R2为0.940 9(图8)，这为开展基于D-fG遥感参数信息的动态收获指数遥感估算奠定了基础。

图8 基于D-fG参数的D-HI估算模型Fig.8 D-HI estimation model based on D-fG parameter

3.3 基于模拟Sentinel-2A波段数据的冬小麦收获指数估算

3.3.1冬小麦D-fG参数估算和精度验证

利用Sentinel-2A光谱响应函数模拟地面高光谱波段所对应的卫星遥感反射率数据，并以模拟反射率结果构建相应的NDSI；然后，构建模拟NDSI与实测D-fG之间的估算模型，并估算冬小麦D-fG。最后，得到基于模拟Sentinel-2A波段构建NDSI与D-fG间统计关系及D-fG估算精度，结果如表4所示。从表4可以看出，9个模拟Sentinel-2A波段构建的NDSI拟合D-fG在P<0.01水平上均达到极显著水平，估算模型R2在0.663 6～0.836 2之间。通过验证数据集(n=108)进行精度检验可知，其精度均达到了较高水平。其中，精度验证R2在0.908 9～0.925 9之间，RMSE在0.038 6～0.048 2之间，NRMSE在11.18%～13.96%之间，MRE在10.07%～12.31%之间。可以看出，相较于红、绿、蓝波段的组合，红光波段、红边波段与近红外波段间组合的HI估算精度普遍较高，这主要由于红边波段是指示绿色植物生长状况的敏感波段，能够有效反映植被养分状况、健康状态和生理生化参数等信息。其中，λ1(672～680 nm)和λ2(855～875 nm)模拟Sentinel-2A反射率构建的NDSI进行冬小麦D-fG估算的精度最高，其RMSE、NRMSE和MRE分别为0.038 6、11.18%、10.07%，这主要是由于近红外与红光波段交界处快速变化的区域能够对植被冠层结构和叶绿素含量等微小变化和植被生长状况进行有效反映，同时，红光和近红外波段反射率有明显的反差。因此，红光和处于近红外区域的窄近红外波段组合构成的NDSI估算D-fG精度更高，效果更为理想。

3.3.2冬小麦动态收获指数估算和精度验证

在模拟波段构建的NDSI估算D-fG的条件下，根据式(6)分别计算冬小麦D-HI估算结果，并进行精度验证。冬小麦动态收获指数总体精度验证结果如表5所示。从表5可知，估算结果验证均达到了高精度水平，其精度验证R2在0.889 9～0.909 7之间，RMSE在0.040 4～0.051 5之间，NRMSE在10.83%～13.81%之间，MRE在9.56%～12.38%之间。其中，红光波段、红边波段与近红外波段间组合进行D-HI估测精度均较高，且基于波段λ1(672～680 nm)和λ2(855～875 nm)模拟Sentinel-2A反射率数据估算D-fG参数的D-HI估测结果精度最高，RMSE、NRMSE和MRE分别为0.040 4、10.83%、9.56%。由于红光和近红外波段的反射率具有明显的反差，故由上述红光波段和窄近红外波段构成的NDSI对不同植被光谱的变化更为敏感，能很好地反映作物长势、生长状况以及D-fG，加之D-fG和D-HI间具有较好的正相关关系，因此，基于波段λ1(672～680 nm)和λ2(855～875 nm)获得的D-HI估算精度最高。上述基于模拟Sentinel-2A遥感数据波段组合进行D-HI估测的精度结果，对基于真实遥感数据的D-HI估测中遥感数据源选择、波段组合优选具有重要指导意义。

表4 基于模拟Sentinel-2A波段构建NDSI的D-fG估算精度验证Tab.4 Verification of estimation accuracy of D-fG based on NDSI constructed by simulated bands of Sentinel-2A

表5 基于模拟Sentinel-2A波段的D-HI估算模型总体 精度验证Tab.5 Verification of overall accuracy of D-HI estimation model based on simulated bands of Sentinel-2A

3.4 基于Sentinel-2A卫星遥感数据的冬小麦收获指数估算

3.4.1冬小麦D-fG空间信息估算和精度验证

在基于模拟Sentinel-2A卫星遥感数据获得的D-fG和D-HI估算最高精度波段优选结果基础上，利用Sentinel-2A真实卫星遥感数据进行冬小麦收获指数估算。在利用宽波段Sentinel-2A卫星遥感数据估算冬小麦D-fG过程中，首先根据地面采样点的GPS定位信息对预处理后的遥感影像进行反射率提取，构建相应的NDSI与实测D-fG间的模型；然后，根据地面高光谱敏感波段所对应的Sentinel-2A波段，直接运用波段计算获取相应NDSI的空间分布，在此基础上实现区域冬小麦D-fG的空间提取并进行精度验证。基于Sentinel-2A的红光波段和窄近红外波段反射率构建的NDSI与D-fG间拟合模型及D-fG估算精度如图9、10所示。其中，红光波段和窄近红外波段构建NDSI进行冬小麦D-fG空间信息估算的R2为0.907 1，RMSE、NRMSE、MRE分别为0.044 3、13.13%、11.90%。本研究基于上述最高精度的Sentinel-2A影像宽波段构建的NDSI实现了衡水地区冬小麦D-fG空间信息分布获取，其中，灌浆前期、灌浆后期、成熟期D-fG空间估算结果如图11所示。通过分析可知，灌浆前期D-fG估算值在0.19～0.51之间，研究区D-fG平均值为0.23；灌浆后期D-fG估算值在0.27～0.62之间，研究区D-fG平均值为0.40；成熟期 D-fG估算值在0.44～0.68之间，研究区D-fG平均值为0.53。

图9 基于Sentinel-2A卫星红光和窄近红外波段构建 NDSI的D-fG估算模型Fig.9 D-fG estimation model based on NDSI constructed by red band and NIR2 band of Sentinel-2A

图10 基于Sentinel-2A红光和窄近红外波段卫星 遥感数据的冬小麦D-fG精度验证结果Fig.10 Verification of D-fG of winter wheat based on red band and NIR2 band of Sentinel-2A

图11 基于Sentinel-2A卫星遥感红光和窄近红外波段的冬小麦D-fG空间估算结果Fig.11 Spatial estimation results of winter wheat D-fG based on red band and NIR2 band of Sentinel-2A

3.4.2冬小麦收获指数空间信息估算和总体精度验证

在冬小麦D-fG空间信息获取的基础上，利用Sentinel-2A的红光波段和窄近红外波段构建的NDSI实现了衡水地区冬小麦D-HI空间信息估算，并进行总体精度验证。不同时期D-HI空间信息估算结果如图12所示，D-HI估算总体精度验证结果如图13所示。由图12可知，冬小麦灌浆前期D-HI估算值在0.25～0.49之间，研究区D-HI平均值为0.28；冬小麦灌浆后期D-HI估算值在0.31～0.58之间，研究区D-HI平均值为0.41；冬小麦成熟期D-HI估算值在0.44～0.62之间，研究区D-HI平均值为0.51。通过验证，基于Sentinel-2A卫星遥感红光波段和窄近红外波段估算D-fG参数的D-HI估算结果精度验证R2为0.879 8，RMSE为0.050 2，NRMSE为13.81%，MRE为12.00%。上述基于红光波段和处于近红外区域的窄近红外波段D-HI估算结果对宽波段多光谱卫星遥感数据选择和传感器波段设置具有重要指导意义。

图12 基于Sentinel-2A卫星遥感红光和窄近红外波段的冬小麦D-HI空间信息估算结果Fig.12 Spatial estimation results of D-HI based on red band and NIR2 band of Sentinel-2A

图13 基于Sentinel-2A卫星遥感红光和窄近红外 波段的冬小麦D-HI估算精度验证结果Fig.13 Accuracy verification of D-HI estimation model based on red band and NIR2 band of Sentinel-2A

4 讨论

与以往利用遥感信息获取收获指数精度指标RMSE一般在0.01～0.06之间相比[15,17,22,24]，本研究基于Sentinel-2A模拟遥感数据的作物收获指数估算取得了较高的估算精度，且估算D-HI最高精度R2为0.908 3，RMSE为0.040 4，NRMSE为10.83%，而基于真实卫星遥感数据和相关最优波段的D-HI估算中，总体验证精度R2为0.879 8，RMSE为0.050 2，NRMSE为13.81%。上述结果一方面说明本研究所提方法在田间冠层尺度扩展到区域尺度范围内理论上是可行的，另一方面也说明在基于卫星遥感的区域收获指数估算研究过程中，由于受大气影响造成卫星遥感反射率数据出现一定误差或波动，从而在一定程度上降低了作物收获指数估算模型的应用精度。虽然基于真实卫星遥感数据估算D-HI的精度略有下降，但仍然能够满足区域农业遥感监测中收获指数获取的精度要求，对于大范围区域农作物收获指数空间信息获取仍具有一定指导意义。同时，与以往仅实现成熟期收获指数估算比较[15,17,22,24]，本研究所提基于花后累积地上生物量比例的区域冬小麦收获指数遥感估算方法不仅可以实现成熟期收获指数的准确估算，而且可以实现不同灌浆阶段作物收获指数动态变化空间信息获取。从应用前景看，本研究所提收获指数获取技术方法不仅对提高精准农业农田作物管理能力具有重要意义，对开展作物生长过程定量模拟、作物品种选育和表型信息获取、作物生长环境评价、作物灾害损失定量评价、农业对气候变化的响应等研究和应用也具有重要实用价值[43]。此外，本研究目前只针对冬小麦开展了收获指数研究，但对于其他作物，特别是禾本科粮食作物(如水稻、玉米、高粱等)收获指数遥感估算也具有重要的参考价值。

本研究在利用前期作物冠层高光谱确定的冬小麦D-fG敏感波段中心以及最大波段宽度结果基础上，基于近地高光谱模拟Sentinel-2A数据及其真实卫星遥感数据开展冬小麦收获指数遥感估算研究取得了较好的结果，但研究过程中忽略了一些外部因素的影响，且部分问题还需要进一步深入研究。如筛选出的敏感波段最大波宽所对应的波段范围与遥感数据的波段不能完全对应，这在一定程度上影响了最大波段宽度的实际效果；其次，多光谱遥感数据构建的光谱指数与敏感波段构建的光谱指数在数量以及位置上难以一一对应，这使得部分光谱指数组合形式在作物收获指数估算过程中无法进一步应用；再者，研究中Sentinel-2A遥感数据过境时间与地面实际采样时间存在一定时间差，考虑到较短时间差内作物生长变化较小，因此，本研究尚未考虑遥感数据过境时间与地面采样时间不完全一致对收获指数估算精度的影响，但今后需进一步开展卫星遥感和地面观测同步的作物收获指数估算精度验证研究。此外，本研究在冠层尺度冬小麦收获指数敏感波段中心及其最大波宽优选结果基础上，系统开展了Sentinel-2A模拟数据不同波段组合下的D-HI估测，对今后基于真实遥感数据的D-HI估测遥感数据选择、波段组合效果评价、卫星传感器波段应用与设置(如红边波段)等具有一定指导意义，但本文仅利用Sentinel-2A真实遥感波段数据最优波段组合获取了高精度的区域收获指数空间信息，然而基于真实卫星数据其他波段组合的区域作物收获指数空间信息估测精度还有待进一步研究和评价。同时，除了Sentinel-2以外，在考虑作物物候信息基础上，如何进一步利用具有红边的其他卫星(如RapidEye、国产GF-6等)进行更大范围的作物收获指数估算也将是今后的重要研究内容。

5 结论

(1)在作物收获指数估算敏感波段中心和最大波段宽度基础上，基于Sentinel-2A模拟遥感反射率数据，通过构建不同波段组合NDSI与实测D-fG间统计模型，实现了不同波段组合条件下冬小麦D-fG参数的遥感准确估算。在此基础上，基于D-fG与D-HI间统计关系模型和D-fG参数遥感信息，实现了波段组合条件下灌浆期不同阶段和成熟期冬小麦D-HI的准确估测和最优波段组合的筛选。其中，红光波段、红边波段与近红外波段间组合进行D-HI估测精度均较高，且基于模拟反射率波段λ1(672～680 nm)和λ2(855～875 nm)估算D-HI的精度最高，RMSE、NRMSE和MRE分别为0.040 4、10.83%、9.56%。上述模拟遥感数据不同波段组合的D-HI估测结果说明，Sentinel-2A卫星遥感数据在基于花后累积地上生物量比例的冬小麦收获指数估算中具有一定的应用潜力。

(2)在基于模拟Sentinel-2A多光谱遥感数据进行作物动态收获指数估算筛选出的最高精度波段组合信息基础上，利用真实Sentinel-2A卫星遥感数据开展了花后累积地上生物量比例的作物D-HI遥感估测应用，实现了区域冬小麦D-fG参数遥感获取和区域D-HI准确估测。其中，利用Sentinel-2A的红光波段(R)和窄近红外波段(NIR2)进行区域D-HI估测的RMSE、NRMSE和MRE分别为0.050 2、13.81%、12.00%。上述区域冬小麦动态收获指数估算精度能够满足农业遥感监测中对收获指数动态信息估测精度的要求，证明了基于Sentinel-2A卫星遥感数据和花后累积地上生物量比例的作物动态收获指数估算方法在大范围作物D-HI空间信息获取中的可行性和有效性。