基于NDVI与EVI的作物长势监测研究
2019-10-10白燕英高聚林张宝林
白燕英 高聚林 张宝林
(1.内蒙古农业大学水利与土木建筑工程学院, 呼和浩特 010018; 2.内蒙古农业大学农学院, 呼和浩特 010019;3.内蒙古师范大学化学与环境科学学院, 呼和浩特 010018)
0 引言
植被是陆地最重要的生态系统,它覆盖了地球大部分陆地表面并强烈地影响地球的生态环境[1]。研究区域尺度植被的时空动态变化具有重要科学和应用价值。遥感影像由于其覆盖面积广,能频繁获取地表信息,已被广泛应用于区域尺度植被覆盖变化的动态监测研究[2]。植被指数是遥感领域中用来表征地表植被覆盖的简单有效度量参数[3],被广泛应用于植被类型监测[4-7]、长势监测[8]、植被覆盖度监测[9-12],及生物量估算[13]、作物估产[14]等研究。
归一化植被指数NDVI已经成为监测植被覆盖情况和生长状况的最佳遥感指数[2],但基于NIR和Red比值的NDVI算式是以其易饱和为代价来减少大气的影响,主要表现在对大气干扰的处理有限,在低植被覆盖区易受土壤背景和植被冠层的影响而偏高,而在植被高覆盖区则容易饱和[15-19]。植被指数饱和是指在密集植被区,当植被越来越茂密时,植被指数不再随生物量的增加而增长的现象[19]。增强型植被指数EVI采用“抗大气植被指数”和“土壤调节植被指数”削弱了气溶胶和土壤背景的影响,并克服了植被饱和现象,能敏感地监测稀疏植被和茂密植被生长与衰落状况,可以弥补NDVI的不足[20-23]。
然而关于植被指数NDVI与EVI的差异研究主要基于MODIS和NOAA/AVHRR等数据产品,该数据时间分辨率高,空间分辨率低,适用于大中尺度区域的研究。基于中高分辨率影像的NDVI和EVI监测植被差异的研究鲜见报道,对二者定量化监测植被的比较研究更少。本文基于2015年大气校正后的Landsat8影像,辅以Landsat7影像,以像元为基本计算单元,定量分析NDVI和EVI随植被覆盖度增加的变化规律及二者关系比较;分析低、中、高植被覆盖度下,NDVI与EVI对作物长势变化的监测效果及其差异性,为作物类型及长势监测选取适宜的植被指数提供科学依据。
1 研究区概况及数据处理
1.1 研究区概况
土默特右旗位于内蒙古包头市东南部,北倚阴山山脉,南邻黄河,本文以土默特右旗阴山山脉以南的平原区为研究区,地理位置范围110°13′~111°8′E,40°14′~40°39′N,面积1 667.5 km2(图1)。研究区属于典型温带大陆性气候,降水量少,蒸发量大,年平均降水量357.6 mm,年平均蒸发量1 947.86 mm;主要依赖黄河水灌溉,土地平坦,水源充足,适宜多种农作物种植。
图1 研究区地理位置分布Fig.1 Geographical distribution of study area
1.2 遥感影像数据及处理
选用Landsat8 OLI影像,辅以Landsat7 ETM+影像为遥感数据源,共筛选出2015年研究区云量较少的覆盖作物生育期的12景影像(表1),影像尺寸178 km×178 km,多光谱、全色波段空间分辨率分别为30、15 m。过境周期16 d,研究区所在行列号为127-32,过境时间为北京时间11:18左右。影像数据预处理主要包括辐射定标、大气校正、图像裁剪、去云处理和对ETM+的去条带处理。选取6S模型对影像进行大气校正,采用SPLINE线性插值法进行去云和去条带处理。
1.3 研究区作物物候
土默特右旗主要种植小麦、玉米、豆类、马铃薯等粮食作物和葵花、西葫芦、西瓜、甜菜及蔬菜等经济作物。2015年小麦、玉米、葵花和西葫芦种植面积约占农作物总种植面积的89.2%,其中玉米种植面积最大,约占总播种面积的70%。2015年Landsat8影像过境时间对应的主要作物生育期如下:①4月28日小麦处于苗期。②5月30日小麦处于拔节期,玉米和西葫芦处于苗期。③6月15日小麦处于生育高峰期(孕穗期),玉米处于拔节期,西葫芦处于初花期。④7月1日小麦开始乳熟,长势开始衰减,玉米进入大喇叭口期,西葫芦进入生育高峰期(结果期),葵花处于苗期。⑤7月25日小麦已经全部收割,玉米处于生育高峰期(抽雄期),西葫芦果实逐渐成熟,长势开始衰落,葵花处于拔节期。⑥8月2日玉米进入灌浆期,长势开始衰落。西葫芦果实开始成熟,葵花进入生育高峰期(开花期)。⑦8月26日玉米进入乳熟期,西葫芦已经收获,葵花正值生育高峰期(开花期),长势茂密。⑧9月3日玉米开始成熟,葵花进入成熟期,长势开始衰落。⑨9月19日玉米和葵花都进入收获期。整个生育期内林地和荒草地从4月28日到8月2日长势一直呈现逐渐增加趋势,从8月2日到10月5日呈逐渐下降趋势。
表1 2015年Landsat遥感影像Tab.1 Landsat images of study area in 2015
2 植被指数
对绿色植物强吸收的可见光红波段和对绿色植物高反射的近红外波段对同一植被的光谱响应截然相反,二者形成明显的反差,这种反差随着叶冠结构、植被覆盖度而变化。红色和近红外波段的反差是对植物量很敏感的度量,常用这两个波段反射率的比值、差分、线性组合等多种组合构建各种植被指数来增强或揭示隐含的植物信息[24]。常见植被指数有比值植被指数RVI、差值植被指数DVI、绿度植被指数GVI、归一化植被指数NDVI、抗大气植被指数ARVI、土壤调节植被指数SAVI、增强型植被指数EVI等。目前应用最为广泛的是归一化植被指数NDVI,其计算公式为
(1)
式中ρNIR——近红外波段反射率
ρred——红波段反射率
为了克服NDVI指数存在的缺陷,文献[24]引入了背景调节参数L和大气修正参数C1、C2,同时减少背景和大气噪声的影响,建立了增强型植被指数EVI,其计算公式为
(2)
式中ρblue——蓝波段反射率
L——土壤调节参数,取1
C1——大气修正红光校正参数,取6.0
C2——大气修正蓝光校正参数,取7.5
土壤的红波段反射率ρred与近红外波段反射率ρNIR接近,随着地表植被覆盖度的增加,其红波段反射率ρred减小,而近红外反射率ρNIR大幅度增加,植被覆盖度越高,红光反射ρred越小,近红外反射ρNIR越大,因此比值植被指数ρNIR/ρred能灵敏指示绿色植物生物量[19,24]。当植被越来越茂密,红波段反射率ρred越小,近红外波段反射率ρNIR越高,导致ρNIR/ρred快速增长[19]。
3 植被指数关系分析
3.1 NDVI与EVI与ρNIR/ρred关系分析
通过分析植被指数NDVI和EVI随ρNIR/ρred的变化规律,揭示二者随植被覆盖度增加的变化规律及差异。土壤的ρNIR/ρred接近1,随着地表开始出现植被,ρNIR/ρred开始逐渐增加,当地表的光谱特征开始表现为植被光谱特征时,ρNIR/ρred达到2[25]。
将ENVI软件计算得到的NDVI和EVI图像以像元为单位输出文本格式,每个像元输出一个X、Y坐标及植被指数。用Matlab软件读入文本数据,绘制NDVI、EVI与ρNIR/ρred关系图,见图2a;绘制ρNIR/ρred与NDVI增量、EVI拟合曲线增量关系图,见图2b。
图2a表明NDVI随着ρNIR/ρred的增加而增加,这是由式(1)本身决定的。EVI引入调节土壤背景参数和大气修正参数,随着ρNIR/ρred的增加呈幂函数趋势增加,图中用其拟合曲线说明其变化情况。当地表为裸土,ρNIR/ρred为1时,NDVI与EVI均为0;当地表刚出现植被时,NDVI和EVI的增加速度最快,当ρNIR/ρred达到2时, NDVI升到0.333,EVI升到0.234,NDVI的增加速度明显大于EVI的增加速度。图2b表明随着ρNIR/ρred的增加,NDVI和EVI持续增加,但增加值逐渐减小,当ρNIR/ρred达到10以后,NDVI与EVI增加极为缓慢。从地表开始出现植被到植被长势茂密,NDVI的增加幅度一直高于EVI的增加幅度,这是植被无论在生长初期还是生育高峰期NDVI值都高于EVI值的根本原因。
图2 植被指数与ρNIR/ρred的关系(2015年7月1日影像)Fig.2 Relationship between ρNIR/ρred and VI(Remote sensing image on July 1, 2015)
3.2 NDVI、EVI与植被覆盖度关系分析
采用二分模型法计算植被覆盖度,计算公式为
(3)
(4)
式中fcNDVI——NDVI值计算的植被覆盖度
fcEVI——EVI值计算的植被覆盖度
NDVIsoil——裸土的NDVI值
EVIsoil——裸土的EVI值
NDVIveg——植被的NDVI值
EVIveg——植被的EVI值
3.2.1纯裸土和纯植被的植被指数的确定
2015年5月14日遥感影像的裸土像元较多,用于确定纯裸土像元的NDVI和EVI值比较理想。经过统计分析,像元分布频率高于1%裸土的NDVI值范围为0.08~0.15,EVI值范围为0.06~0.1,该范围的裸土像元数分别占总裸土总像元数的72.3%和62.4%。裸土的NDVI和EVI值主要受土壤背景明暗程度的影响,土壤背景越暗,NDVI和EVI越高,土壤背景越亮,NDVI和EVI越低。本文取纯裸土像元NDVI值为0.15,EVI值取0.1。2015年7月1日遥感影像高植被覆盖像元较多,用于确定纯植被像元的NDVI和EVI值比较理想,因此根据图2a确定纯植被像元,选取NDVI和EVI的饱和值作为纯植被像元的NDVI值和EVI值,分别为0.95和0.7。
3.2.2植被覆盖度的计算
将研究区植被覆盖度分为4个等级:低植被覆盖(0 经过统计分析,低植被覆盖地表NDVI和EVI值范围分别为0~0.4和0~0.3,中低植被覆盖地表NDVI和EVI值范围分别为0.4~0.55和0.3~0.4,中高植被覆盖度地表的NDVI和EVI值范围分别为0.55~0.8和0.4~0.6,高植被覆盖地表的NDVI和EVI值范围分别为0.8~1.0和0.6~1.0。 表2 植被指数对应的植被覆盖度Tab.2 Vegetation coverage corresponding to VI NDVI和EVI的计算方法决定了NDVI始终高于EVI,二者差值随植被长势的增加而增加,当植被覆盖度达到30%时,NDVI与EVI的差值为0.1,植被覆盖度达到50%时,NDVI与EVI的差值达到0.2,植被覆盖度达到80%时,NDVI与EVI的差值增加到0.27,植被覆盖度达到100%时,NDVI与EVI的差值减小到0.26。 根据NDVI与EVI散点图,分析研究区作物生育前期(5月30日)、生育高峰期(7月1日)及生育末期(9月19日)3个时期的NDVI和EVI关系,见图3。图3表明3个时期的散点图分布规律基本一致,呈正比例散点图分布,且主要分布在1∶1直线上方,说明绝大部分像元的NDVI值均高于EVI值。3个时期的散点图分布略有区别,作物生育前期(5月30日),高植被覆盖地表少,因此在高值区的散点少。7月1日的高植被覆盖地表较多,植被指数较高,高值区的散点多而密集。9月19日很多作物已经开始收获,植被长势明显枯萎,在高值区的散点明显减少。 图3 EVI与NDVI散点关系图及拟合曲线Fig.3 Scatter diagrams and fitting curves of EVI and NDVI EVI与NDVI散点图的拟合曲线形式为二次函数,根据拟合曲线斜率将EVI与NDVI散点图的拟合曲线分为3段,EVI为0~0.3、0.3~0.6、0.6~1.0。拟合曲线在EVI小于0.3范围内斜率大于1∶1直线,说明低植被覆盖下随着植被覆盖度的增加,NDVI增加速度高于EVI值增加速度。拟合曲线在EVI为0.3~0.6范围内斜率接近1∶1直线,说明在中等植被覆盖区随着植被覆盖度的增加,NDVI增加速度与EVI的增加速度基本一致。拟合曲线在EVI为 0.6~1.0范围内斜率小于1∶1直线,说明在高植被覆盖区随着植被覆盖度的增加,NDVI的增加速度低于EVI增加速度。 不同NDVI(EVI)值对应的像元数与研究区总像元数比值即NDVI(EVI)值的分布频率,植被指数对应的分布频率越高,表明该植被指数对应的像元数越多,植被指数对应的分布频率越低,表明该植被指数对应的像元数越少。作物不同生育时期遥感影像的NDVI和EVI分布频率曲线见图4和图5。 图4 不同时期植被指数分布频率曲线Fig.4 VI distribution frequency curves in different periods 图5 同一时期NDVI与EVI分布频率曲线比较Fig.5 Comparison of NDVI and EVI distribution frequency curves in the same time period 图4表明不同时期的NDVI(EVI)分布频率曲线差异较大,5月30日NDVI和EVI分布频率曲线在低值区呈尖峰状,说明NDVI和EVI的分布比较集中,表明此时地表大部分为裸土和低植被覆盖。6月15日曲线峰值明显下降且右移,但曲线峰值对应的植被指数仍属于低植被覆盖范围,说明地表仍以低植被覆盖为主,但地表的植被覆盖明显增加。从7月1日开始,曲线明显下移且分布范围扩大,表明高植被覆盖地表明显增加,此时大部分作物处于生育高峰期,长势茂密。9月3日分布频率曲线从右侧收缩,表明长势茂盛的作物开始衰落,植被指数下降。9月19日分布曲线从右侧继续收缩,右侧峰值左移,说明作物长势继续衰落,植被指数继续下降。 比较不同时期的NDVI(EVI)分布频率曲线目的在于比较不同时期地表植被覆盖度的变化,结果表明:不同时期NDVI(EVI)分布频率曲线对比能明显反映地表植被覆盖度的变化。 图5表明同一时期NDVI与EVI分布频率曲线形状相似,由于NDVI值域范围宽为0~0.95,EVI值域范围窄为0~0.7,因此EVI频率曲线分布比NDVI频率曲线分布集中。由于同一像元NDVI值高于EVI,且随着地表覆盖度的增加,NDVI与EVI的差值逐渐增加,在图上表现为EVI曲线较NDVI曲线左移,随着植被覆盖度的增加,左移距离明显增加。NDVI与EVI分布曲线差异表现在: 6月15日NDVI曲线峰值明显低于EVI曲线峰值,9月3日和9月19日NDVI曲线左侧峰值略低于EVI曲线左侧峰值,说明NDVI对地表开始生长植被反映敏感,估计值偏高,因此刚出现植被的地表像元NDVI值增加较EVI快,分布范围较EVI广,导致NDVI峰值低于EVI峰值。 图6 同一作物NDVI与EVI时序曲线比较Fig.6 Comparison of the same crop NDVI and EVI time series curves 比较同一时期NDVI与EVI分布频率曲线的目的在于揭示NDVI和EVI指示研究区地物不同植被覆盖度的分布结果。分析结果表明:同一时期的NDVI与EVI的频率分布曲线形状相似,说明NDVI和EVI描述不同植被覆盖度分布的结果基本一致。 3.5.1同一作物NDVI与EVI时序曲线比较 在不同时相NDVI图像和EVI图像上提取作物采样点的NDVI值和EVI值,构建作物时间序列NDVI曲线和EVI曲线,见图6。图6表明每种作物的NDVI时序曲线和EVI时序曲线趋势基本一致,都能清晰反映作物的生长和衰落变化。作物的EVI曲线始终低于NDVI曲线,二者差值随长势增加而增大,随长势衰落减小。在作物生长初期或中低植被覆盖下,两种曲线变化趋势一致,对植被描述能力相似。 在高植被覆盖下,两种曲线描述作物生长或衰落有所差异。主要表现在作物继续生长时,NDVI增长趋势不如EVI增长趋势明显,如玉米从7月1日大喇叭口期到7月25日抽雄期,小麦从5月22日拔节期到6月15日孕穗期。高植被覆盖作物开始成熟时NDVI呈继续增加趋势而EVI则呈下降趋势,如玉米从7月25日抽雄期到8月2日灌浆期,小麦从6月15日孕穗期到7月1日乳熟期。说明由于NDVI的饱和性,高植被覆盖度作物继续生长时,NDVI增加不如EVI增加明显;高植被覆盖作物衰落时,EVI呈现下降趋势而NDVI表现为继续增加。 NDVI时序曲线在作物长势明显增加或衰落时呈突增或突降趋势,EVI曲线变化较NDVI曲线平缓,主要因为NDVI仅与ρNIR/ρred有关,当植被长势变化明显时,其ρNIR/ρred变化明显导致NDVI突增或突减。EVI由于引入调节土壤背景参数和大气修正参数,计算结果较为稳定,不会因为ρNIR/ρred的较大变化而呈突增或突降。 3.5.2不同作物间NDVI及EVI时序曲线比较 分别将上述研究的6种植被(作物)典型采样点的NDVI、EVI时序曲线放在一起分析其差异,比较结果见图7。 图7 不同作物植被指数时间序列曲线比较Fig.7 Comparison of VI time series curves of different crops 图7表明不同植被(作物)NDVI和EVI时序曲线的差异规律基本一致,都能清晰反映植被(作物)在生育期内随长势变化的规律和不同植被(作物)在同一时期的长势差异。 两条时序曲线的明显差异表现在植被(作物)从6月15日到7月1日的变化趋势,小麦的NDVI曲线表现为增加趋势,EVI曲线显示为降低趋势,原因在于6月15日小麦处于生育高峰期孕穗期,由营养生长变为生殖生长,此时长势最旺盛,此后开始逐渐乳熟,长势逐渐衰落。7月1日玉米、西葫芦、小麦和林地的NDVI值接近,难以区分,但它们的EVI值差异明显,西葫芦最高,玉米和林地接近,小麦略低。上述现象表明NDVI由于其易饱和性对高植被覆盖作物长势变化反应不敏感,而EVI对高植被覆盖作物长势变化反应敏感。 (1)图2植被指数与ρNIR/ρred关系和图3b EVI与NDVI散点图都显示植被茂盛期7月7日NDVI的取值范围在0~0.95,EVI的取值范围在0~0.95,但图5频率分布曲线显示7月1日植被生育高峰期EVI超过0.7的像元分布频率接近0,达到0.7的像元分布频率为0.02%,大于0.7的像元数共有4 684个像元,累积值仅占研究区面积0.25%,因此EVI最大值取0.7。NDVI超过0.95的像元分布频率接近0,NDVI最大值取0.95。 (2)在植被生长初期或低植被覆盖下,NDVI和EVI增加都很快,NDVI的增加速度明显高于EVI的增加速度,导致作物生长初期NDVI和EVI都过高地估计植被覆盖度,NDVI估计值略高于EVI的估计值。在植被生育中期或中等植被覆盖下,NDVI和EVI随植被长势增加速度接近,对植被描述能力相似。在植被生育高峰期或高植被覆盖下,EVI的增加速度高于NDVI,植被继续生长或开始衰落时,NDVI值变化不如EVI变化明显,甚至有的情况下当植被开始衰落时,EVI值显示下降,NDVI值仍显示增加,说明NDVI由于其公式决定的高植被覆盖下易饱和性,不能敏感监测植被长势变化,而EVI则能敏感监测高植被覆盖植被长势变化。 (1)定量分析了NDVI与EVI随植被覆盖度变化关系。地表刚出现植被时,NDVI和EVI的增加速度最快,随着地表植被覆盖度的增加,NDVI与EVI的增加速度减缓,高植被覆盖时,二者增加速度极为缓慢。NDVI和EVI的公式决定了NDVI始终高于EVI,二者差值随植被长势的增加而增加,当植被覆盖度达到30%时,NDVI与EVI的差值为0.1,植被覆盖度达到50%时,NDVI与EVI的差值达到0.2,植被覆盖度达到80%时,NDVI与EVI的差值增加到0.27,植被覆盖度达到100%时,NDVI与EVI的差值减小到0.26。 (2)量化了低、中、高植被覆盖度下NDVI和EVI的取值范围。低植被覆盖地表NDVI和EVI值范围分别为0~0.4和0~0.3,中低植被覆盖地表的NDVI和EVI值范围分别为0.4~0.55和0.3~0.4,中高植被覆盖地表的NDVI和EVI值范围分别为0.55~0.8和0.4~0.6,高植被覆盖地表的NDVI和EVI值范围分别为0.8~1.0和0.6~1.0。 (3)比较不同时期的NDVI(EVI)分布频率曲线目的在于比较不同时期地表植被覆盖度的变化,结果表明:不同时期NDVI(EVI)分布频率曲线对比能明显反映地表植被覆盖度的变化。比较同一时期NDVI与EVI分布频率曲线的目的在于揭示NDVI和EVI指示研究区地物不同植被覆盖度的分布结果。分析结果表明:同一时期的NDVI与EVI的频率分布曲线形状相似,说明NDVI和EVI描述不同植被覆盖度分布的结果基本一致。 (4)通过比较几种作物NDVI与EVI时间序列曲线发现,NDVI时序曲线在作物长势明显增加或衰落时呈突增或突降趋势,EVI由于引入调节土壤背景参数和大气修正参数,计算结果较为稳定,不会因为ρNIR/ρred的较大变化而呈突增或突降,更符合作物生长发育特点。在作物生长初期或低植被覆盖下,NDVI、EVI都过高估计植被覆盖度,NDVI估计值略高于EVI的估计值。在作物与生育中期或中等植被覆盖下,对作物描述能力相似。在作物生育高峰期或高植被覆盖下,监测作物长势变化EVI比NDVI更敏感。3.3 NDVI与EVI散点图及拟合曲线关系分析
3.4 NDVI与EVI分布频率比较
3.5 作物NDVI与EVI时序曲线比较
4 讨论
5 结论