Sentinel-2宽波段光谱指数预测拔节后受冻冬小麦减产率

赵爱萍，马浚诚，武永峰*，胡 新，任德超，李崇瑞

1. 中国农业科学院农业环境与可持续发展研究所，北京 100081 2. 商丘市农林科学院小麦研究所，河南 商丘 476000

引 言

低温引起的作物灾害，是世界性的农业气象灾害。 霜冻是中国发生范围广、 危害作物种类多、 造成经济损失大的一种气象灾害，北至黑龙江，南至广东、 广西均存在霜冻灾害。 晚霜冻害是冬小麦进入拔节期以后所发生的零下低温灾害，此时正是冬小麦幼穗形成的关键期，对低温敏感，霜冻后幼穗最先受害，无法正常结实导致减产[1]。 减产率能较好地表征作物在低温胁迫下受冻害的程度，对灾情评估和后续田间管理具有指导作用。 进行霜冻胁迫后冬小麦减产率早期预测，对灾害监测和生产管理决策具有重要的现实意义。

基于不同遥感载荷平台获取的高光谱和多光谱反射率数据，可分别构建各种组合形式的窄波段与宽波段光谱指数，以增强作物冠层反射光谱的敏感性，从而捕获灾害胁迫状态下的更多作物变化信息。 与多光谱相比，高光谱具有光谱分辨率高(波段宽度<10 nm)、 波段连续性强(在400～2 500 nm范围内有几百个波段)、 光谱信息量大等特点。 Nuttall等[2]通过定制试验霜箱处理手段，分析认为小麦叶片窄波段光化学反射指数(photochemical response index，PRI)和归一化植被指数(normalized difference vegetation index，NDVI)与低温有显著关系，可用以表征低温胁迫的程度。 高光谱技术能在肉眼观测到小麦植株变化前对霜冻害进行早期检测[3]，而基于卫星遥感的多光谱技术则在揭示和评价地物时空变化方面更具优势[4]。 Wang等[5]利用中分辨率成像光谱仪(moderate-resolution imaging spectroradiometer，MODIS)遥感数据计算的植被指数对小麦霜冻敏感性进行研究，发现增强植被指数(enhanced vegetation index，EVI)可以有效反映霜冻胁迫特征。 Alessia等[6]利用Sentinel-2数据评估葡萄园发生晚霜冻事件后的损害情况和恢复时间，研究发现Sentinel-2的红边波段3和EVI等可以有效响应霜冻。 近些年来，多光谱遥感技术在作物监测中得到广泛应用[7]，根据作物特征光谱信息建立光谱指数并构建与作物产量间定量关系的方法常用于农作物产量估算和预测[8]。 陈爱莲等[9]开展了基于Sentinel-2数据的大豆种植分布及产量估算的研究，曹娟等[10]利用Sentinel-2A影像进行了低温冷害下大豆估产以及损失评估的研究。 目前利用Sentinel-2卫星多光谱数据来建立晚霜冻胁迫后冬小麦产量损失预测方法的研究还鲜有报道。

从光谱响应霜冻的机制出发开展Sentinel-2宽波段光谱指数预测拔节后受霜冻胁迫的冬小麦减产率研究，采用人工霜冻模拟试验和光谱重采样技术，通过光谱指数与减产率的线性回归决定系数，候选出响应霜冻敏感性较为显著的Sentinel-2宽波段光谱指数，利用Sentinel-2影像多光谱数据计算候选出的光谱指数，与区域自然霜冻后的冬小麦地面采样点减产率数据做线性拟合，通过精度评价验证所选宽波段光谱指数在空间区域内的适用性，以期为分区域制定冬小麦拔节后霜冻胁迫的最佳应对措施提供参考。

1 实验部分

1.1 试验设计

人工霜冻模拟试验地点位于河南省商丘市农林科学院小麦基地(34°30′21″N，115°41′13″E)，海拔高度52 m，气候类型为大陆性季风气候。 试验时间为2017年-2018年和2018年-2019年冬小麦生长季, 冬小麦品种为周麦22。 分别于2017年10月8日和2018年10月5日盆栽种植，每盆种植小麦11株。 试验用盆为直径约25 cm、 高度约35 cm的两头均开口的圆柱型空心管，播种前埋于试验地大田并使盆栽顶部与地表持平，这样尽量使冬小麦生长的盆栽条件与大田环境一致。 从黄淮麦区晚霜冻发生时间来看，亦以药隔后冻害较为常见[11-12]。 本研究将霜冻处理时间安排在幼穗发育的药隔形成阶段(对应的雌蕊分化阶段为小凹期至柱头伸长期)。 2018年试验设置9个低温梯度(-6～-14 ℃)，分别在幼穗发育进入小凹期和柱头伸长期进行低温处理； 2019年试验也设置9个低温梯度(-5～-13 ℃)，在小凹期进行低温处理； 对照组(CK)不做低温处理，正常生长； 每个低温处理及CK均有3盆重复。 盆栽小麦置于低温气候箱(型号为PU-4K； ESPEC公司，日本)中进行低温处理，处理时间为6 h，最开始的3 h内降至最低设定温度后并保持2 h，之后1 h内升至周围环境温度。 处理后再搬回田间，以继续保持大田生长环境，期间的田间管理依据《农作物品种(小麦)区域试验技术规程》(NY/T 1301-2007)[13]进行。 为避免2018年4月6日-7日发生在商丘地区的自然霜冻事件对人工霜冻试验的影响，在降温当晚采取了塑料大棚覆盖的方式以保护盆栽小麦不受霜冻侵袭。

1.2 数据来源

1.2.1 高光谱数据

采用野外便携式ASD地物光谱仪(FieldSpec®3; ASD 公司，美国)测定冬小麦冠层反射率，所有冠层光谱测定工作均于降温处理结束后3 d内，在晴朗、 无风的中午时段(11:00-14:00)完成。 该光谱仪视场角为25°，波段范围为350～2 500 nm，其中350～1 050 nm光谱采样间隔为1.4 nm，光谱分辨率为3 nm； 1 050～2 500 nm光谱采样间隔为2 nm，光谱分辨率为10 nm。 光谱测定工作前，先对光谱仪预热30 min以上，采集目标光谱前利用白板校正，对同一被测目标采集30条光谱，以其平均值表示被测目标光谱反射率值。 去除1 340～1 450，1 800～2 050和2 350～2 500 nm的异常波段后，在Origin 2018软件中，采用Savitzky-Golay卷积平滑方法进行平滑，窗口点数为20，以消除噪声干扰。 2018年试验的光谱采集时间分别为3月30日和4月7日，2019年试验的光谱采集时间为4月3日。

1.2.2 Sentinel-2数据

2018年4月6日-7日，商丘地区出现了一次明显的降温过程。 商丘各县气象观测台站的最低草面温度记录如下，睢县、 民权县、 商丘市市辖区、 虞城县、 柘城县、 宁陵县、 夏邑县、 永城市分别为-3.1，-8.1，-0.5，-7.2，-5.0，-5.8，-7.9和-8.0 ℃。 在本次自然霜冻事件发生之日至4月底期间，选取商丘地区具备晴天条件的4月20日作为研究日期。 从ESA scihub网站下载经过几何校正的Sentinel-2 L1C大气顶反射率影像，然后通过ESA官网提供的Sen2Cor插件对该影像进行大气校正，得到L2A大气底反射率影像。 提取出单波段影像，使用最近邻插值法将影像的空间分辨率采样至10 m。 基于真彩色合成影像，采用面向对象分类方法提取商丘地区冬小麦分布信息，其空间分布及34个随机采样点位置如图1所示。

图1 冬小麦种植区空间分布和随机采样点分布Fig.1 Spatial distribution of winter wheat planting area and the distribution of random sampling points

1.2.3 小麦测产数据及减产率计算

待冬小麦成熟后，分别开展基于盆栽小麦和大田小麦的产量测定工作。

(2)对样本库中数据进行PCA方法降维。对上浆浓度、上浆流量、清水加入量、浆厚度、洗涤水温度、纸浆硬度进行主元分析，得到载荷矩阵和特征值矢量。寻找第一个载荷矢量中绝对值最大的系数对应的变量，并剔除。再次进行主元分析。

(1)盆栽小麦产量测定。 将盆栽种植的小麦随机选取5株，将小麦植株整株拔起装袋带回实验室，待晾晒风干后分别统计低温处理和对照组小麦植株的粒重。

(2)大田小麦产量测定。 在整个商丘地区的小麦种植区，随机采样34个点(图1)。 每个样点的取样方式为利用手持式GPS记录经纬度信息，采集1 m×1 m范围内4行全部冬小麦植株，将整株拔起装袋并带回实验室，待晾晒风干后统计全部植株的粒重，作为实际产量。 由于大田自然霜冻条件下小麦对照产量难以直接获得，因此本研究通过估算正常生长条件下冬小麦理论产量的办法，来获得对照产量。 具体计算步骤为： 统计每个样点采集的所有小麦植株(包括不抽穗的茎秆)，将其作为正常生长条件下的小麦穗数； 挑选20个无冻害的麦穗作为正常条件下的麦穗，统计其粒数和粒重，进而计算穗粒数和千粒重； 利用以上穗数、 穗粒数和千粒重估算冬小麦的正常产量作为对照产量。

基于以上测产结果，构建人工霜冻模拟试验和大田自然霜冻条件下的冬小麦减产率(PYLF，%)，计算如式(1)所示

(1)

式(1)中，GWnormal为正常生长条件下冬小麦的粒重(即对照)，g； GWfrosted为模拟霜冻胁迫下盆栽冬小麦或自然霜冻胁迫下大田冬小麦的粒重，g。

1.3 光谱重采样

光谱重采样是将地面实测光谱与已知传感器的光谱进行匹配。 本研究将ASD地物光谱仪实测的冬小麦冠层高光谱数据重采样为Sentinel-2卫星携带的多光谱仪传感器所采集的光谱数据。 一般情况下，用于遥感的传感器各波段的光谱响应通常呈中间高、 两边低，光谱响应函数与高斯函数趋势一致，有较高的拟合度，因此根据房秀凤等研究采用高斯函数拟合的方法模拟光谱响应函数[14]。 根据模拟的光谱响应函数进行重采样，重采样计算如式(2)所示

(2)

式(2)中，ρ(λ)为重采样后的Sentinel-2卫星传感器对应波段范围的宽波段反射率；ρASD(λ)为ASD地物光谱仪实测的冬小麦冠层光谱的窄波段反射率；f(λ)为Sentinel-2卫星传感器相应波段范围的光谱响应函数值；λmax和λmin分别为Sentinel-2卫星传感器中各个波段范围的上下限。

1.4 光谱指数选建

1.4.1 已有光谱指数选择

表1 本研究选用的已有光谱指数及其表达式Table 1 Published spectral indices and corresponding expressions selected in this study

Note：RNIR,RR,RG,RBandRSWIRare the refmectance in near-infrared, red, green, blue and short-wave infrared wavebands, respectively;Rlrepresent the spectral reflectance of different wavelength

1.4.2 基于Sentinel-2波段的光谱指数构建

Sentinel-2卫星携带的多光谱仪传感器设置有13个光谱波段，覆盖可见光、 近红外和短波红外光谱范围。 因Sentinel-2第10波段的光谱范围是1 365～1 385 nm，位于高光谱冠层反射率被除去的水汽波段内，因此不参与光谱指数构建。 所选取的Sentinel-2波段基本参数如表2所示，将选取的Sentinel-2 12个波段进行光谱指数构建，按照简单比值(simple ratio，SR)、 简单差值(simple difference，SD)和归一化(normalized difference，ND)3种形式组合构建宽波段光谱指数，计算如式(3)—式(5)所示

SR=Ri/Rj

(3)

SD=Ri-Rj

(4)

ND=(Ri-Rj)/(Ri+Rj)

(5)

式中，Ri和Rj分别为i和j波段对应的光谱反射率。 最后，通过Visual Basic宏语言执行编程任务构建光谱指数。

表2 本研究选取的Sentinel-2波段基本参数Table 2 Essential parameters of the selectedSentinel-2 wavebands in this study

1.5 模型构建及精度验证

首先基于人工模拟霜冻试验，利用光谱重采样手段将ASD FieldSpec 3光谱辐射计获取的冠层反射率模拟为Sentinel-2传感器对应的波段反射率，依据波段反射率计算1.4节选建的光谱指数，构建地面光谱指数与冬小麦减产率间的线性回归模型并筛选出模型决定系数排名前三的光谱指数作为候选指数。 其次基于区域自然霜冻，利用Sentinel-2影像数据计算候选指数，通过地面采样点的实测减产率验证候选指数精度。

选用决定系数(coefficient of determination，R2)和均方根误差(root mean squared error，RMSE)综合评价模型预测精度。R2越大，说明模型拟合程度越高； RMSE越小，说明模型精度越高。

2 结果与讨论

2.1 冬小麦冠层光谱反射率在不同低温胁迫下的差异

冬小麦受到低温胁迫，冠层反射率产生相应的变化。 在2018年的冬小麦小凹期和柱头伸长期2个幼穗发育期设置了低温试验，将对照(CK)和不同低温处理下由ASD地物光谱仪测定的冬小麦冠层反射率光谱数据求平均，绘制随波长变化的反射率曲线如图2(a,b)所示。 由图2可知，随处理温度降低，小凹期和柱头伸长期的冬小麦冠层光谱反射率呈现明显的规律性变化，在近红外区域(760～1 300 nm)主要表现为降低趋势，而在可见光(380～760 nm)和短波红外(1 300～2 350 nm)区域表现为升高趋势，且红边(680～760 nm)呈“蓝移”的现象。 近红外区域反射率还表现出两个明显的特征，一是近红外肩峰(750～900 nm)变得越来越平缓，二是以960 nm为中心的水分吸收波段，其水分吸收特征呈消

图2 低温处理下小凹期(a)和柱头伸长期(b)的冬小麦冠层光谱反射率曲线注: CK 表示对照，不做低温处理，正常生长； 1 340～1 450和1 800～2 050 nm波段由于对水蒸气的强吸收，被去除Fig.2 Reflectance curves of winter wheat canopy during the dimple occurrence (a) andthe stigma elongation (b) under low-temperature treatments

失趋势，这一规律和Wei等[12]的研究结果相似。 在不同低温胁迫下，冬小麦冠层光谱反射率的响应特性，可归因于受冻小麦冠层的叶片失水、 色素含量降低以及细胞结构破坏等[2]。 研究结果显示，冠层光谱反射率能反映不同低温胁迫对冬小麦的影响，为进一步分析冬小麦受冻后产量损失与光谱反射率间的关系提供了理论依据。

2.2 已有光谱指数、 构建宽波段光谱指数与减产率的拟合模型

利用光谱重采样方法将ASD光谱仪测定的冬小麦冠层窄波段反射率模拟为Sentinel-2卫星携带的多光谱仪传感器对应的宽波段反射率，通过重采样前后的反射率计算已有光谱指数，得到已有光谱指数与冬小麦减产率的相关关系(表3)。 由表3可知，基于重采样前的反射率计算的19个已有光谱指数中，除已有光谱指数MCARI外，其余18个已有光谱指数都与冬小麦减产率极显著相关(p<0.001)。 其中，相关性最显著的是对冠层结构敏感的已有光谱指数MTVI2，相关系数达到-0.862。 基于重采样后的Sentinel-2宽波段反射率计算的19个已有光谱指数中，有16个指数达到极显著相关(p<0.001)，相关性最好的也是已有光谱指数MTVI2，相关系数达到-0.850。 结果表明，基于模拟Sentinel-2宽波段光谱反射率计算的光谱指数仍然与霜冻后冬小麦减产率保持了高度相关性。

利用模拟后的Sentinel-2 12个宽波段反射率数据，计算以上19个已有光谱指数以及按ND，SD和SR三种形式任意组合构建的Sentinel-2宽波段光谱指数，然后构建光谱指数与冬小麦减产率间的线性回归方程，再将所建线性回归方程代入验证集中，以检验光谱指数的精度。 将4种形式下光谱指数校正集和验证集的R2由高到低排序，基于排序结果的交集部分分别选择排名前三的R2值，共得到12个光谱指数作为候选光谱指数，结果如表4所示。 候选光谱指数的线性回归方程校正集和验证集的R2均大于0.630，精度较好。 在校正集中，本研究构建的简单差值指数B8-B12精度最好，决定系数R2=0.776，均方根误差RMSE=15.437%； 在验证集中，本研究构建的归一化指数(B7-B11)/(B7+B11)精度最好，决定系数R2=0.851，均方根误差RMSE=8.700%。 此外发现基于ND，SD和SR三种组合形式构建的所有光谱指数中得到9个候选光谱指数，其中有6个候选光谱指数包含短波红外波段(B11和B12)，3个候选光谱指数包含近红外水的吸收波段(B9)，表明与水含量变化密切相关的光谱区域对低温胁迫后冬小麦减产率变化敏感。

表3 已有光谱指数与冬小麦减产率的相关性分析

表4 已有和本研究构建光谱指数预测冬小麦减产率的线性回归方程精度Table 4 Accuracy of linear regression equations of published and constructed spectral indicesin this study on predicting the yield reduction rates of winter wheat

2.3 宽波段光谱指数的区域适用性分析

为进一步检验候选出的Sentinel-2宽波段光谱指数在区域尺度上预测受冻冬小麦减产率的能力，以2018年商丘地区自然霜冻事件为契机，利用Sentinel-2影像波段数据计算候选光谱指数，建立候选光谱指数与地面采样点冬小麦减产率的线性回归，得到精度结果如表5所示。 为减少GPS测量造成的定位误差，本研究针对每个地面采样点提取3×3矩阵像素值并求平均，以平均值作为采样点的反射率值。 结果显示，12个宽波段光谱指数中，有9个光谱指数与采样点冬小麦减产率的线性回归通过p<0.01的显著性检验，表明基于人工霜冻模拟试验候选出的部分Sentinel-2宽波段光谱指数，适用于区域尺度内自然霜冻导致的冬小麦减产率预测。 其中，简单差值指数B8a-B12的线性回归精度最高，R2达0.543，RMSE为8.510%； 其次是简单差值指数B8-B12和已有光谱指数EVI，R2分别为0.492和0.486。 指数B8a-B12和B8-B12相比线性回归精度有所提高，由此可见，在区域尺度上Sentinel-2的窄近红外波段(B8a)提高了低温处理下光谱指数对减产率的敏感性，原因是B8a波段包含的信息量丰富，在植被监测中具有较好的效果。 对比分析ND，SD以及SR三种组合形式构建的Sentinel-2宽波段光谱指数线性回归精度，发现SD组合形式构建光谱指数对冬小麦减产率预测精度高于ND和SR组合形式。 研究发现与基于人工霜冻模拟试验的冬小麦减产率预测精度相比，自然霜冻胁迫后的大田冬小麦减产率预测精度明显偏低，其原因可能是Sentinel-2卫星传感器在接收地物反射率过程中，会受到大气条件、 地形起伏以及传感器本身等影响，导致地物光谱反射率与实际情况有差异[6, 18]。 同时，为直观展示冬小麦实测减产率与预测减产率的线性关系，选取线性回归精度较优的候选光谱指数B8a-B12，B8-B12和EVI制作散点图，如图3所示。 由图3(a,b,c)可知随着实测减产率增大，预测减产率逐渐分布在1∶1线以下，表明预测结果有偏低倾向。 其中，候选光谱指数B8a-B12的散点分布更接近于1∶1线[见图3(a)]，离散程度较小，对冬小麦减产率的预测能力更强。

表5 Sentinel-2宽波段光谱指数预测冬小麦减产率的线性回归精度

图3 本研究构建光谱指数B8a-B12 (a)，B8-B12 (b)和已有光谱指数EVI(c)的冬小麦实测减产率与预测减产率的散点图Fig.3 Scatter plots showing the measured and predicted yield reduction rates of winter wheat based on the constructedspectral indices B8a-B12 (a), B8-B12 (b) in this study and the published spectral index EVI (c)

2.4 受冻冬小麦预测减产率的空间分布

上述分析表明，本研究构建的光谱指数B8a-B12在区域尺度上预测受冻冬小麦减产率的线性回归精度最高。 基于指数B8a-B12估算的商丘地区冬小麦预测减产率的空间分布，如图4所示。 其空间分布统计结果表明，全区域预测冬小麦减产率的平均值为12.93%。 冬小麦减产率小于10%的分布面积占比44.29%，在10%～40%之间的占比54.31%，大于40%的分布面积最小，占比仅为1.4%。 此外，从研究区冬小麦预测减产率的空间分布可以看出，减产率大于10%的区域集中分布于民权县、 宁陵县、 虞城县、 夏邑县和永城市(图4)，受灾相对较重，由这5个区域的气象观测台站记录可知最低草面温度均低于-5.8 ℃。 其中，预测减产率分布面积最广的为民权县，受灾最为严重，其气象观测台站所测最低草面温度也是商丘地区最低为-8.1 ℃。 而睢县、 柘城和商丘市市辖区3个区域，气象台站记录的最低草面温度都在-5.0 ℃以上，明显看出这3个区域减产率大于10%分布面积较少，受灾相对较轻。 表明本研究构建光谱指数B8a-B12预测冬小麦减产率的空间分布与气象观测台站最低草面温度的高低呈现一致趋势。

图4 本研究构建光谱指数B8a-B12预测 冬小麦减产率的空间分布Fig.4 Spatial distribution of the winter wheat predicted yield reduction rates by the constructed spectral index B8a-B12 in this study

3 结 论

以冬小麦为研究对象，利用人工模拟霜冻试验和大田自然霜冻试验，探讨基于Sentinel-2宽波段光谱指数预测拔节后受冻冬小麦减产率的可行性和在区域尺度上的适用能力，得到主要结论如下：

(1)在人工模拟霜冻条件下，利用ASD地物光谱仪实测的高光谱数据模拟Sentinel-2宽波段反射率数据，计算已有光谱指数、 以及以简单比值(simple ratio，SR)、 简单差值(simple difference，SD)和归一化(normalized difference，ND)3种形式组合构建的光谱指数中候选出12个精度较优的光谱指数。 其中，在校正集和验证集中综合表现最好的是本研究构建的简单差值指数(B8-B12)以及归一化指数[(B7-B11)/(B7+B11)]，其验证集决定系数R2分别为0.717和0.851，表明Sentinel-2宽波段候选光谱指数可用于预测受冻冬小麦的减产率。

(2)在自然霜冻条件下，利用Sentinel-2多光谱影像数据计算的12个候选宽波段光谱指数其中有9个光谱指数的预测精度达到0.01显著性水平。 预测精度较高的是基于简单差值(simple difference，SD)形式构建的2个Sentinel-2宽波段光谱指数，即B8a-B12和B8-B12，其线性方程的决定系数R2分别为0.543和0.492，均方根误差RMSE分别为8.510%和8.971%。 表明基于Sentinel-2影像数据计算的光谱指数预测具有可行性，候选出的Sentinel-2宽波段光谱指数可用于区域尺度内自然霜冻导致的冬小麦减产率预测。