基于MODIS-Landsat时空融合的陕北黄土高原植被覆盖变化研究

刘咏梅，马 黎,黄 昌,凯 楠

(1.西北大学 城市与环境学院, 陕西 西安 710127；2.水利部 黄土高原水土流失过程与控制重点实验室，河南 郑州 450003)

植被作为气候和土壤之间的联系在全球气候变化和物质能源流动的交换中发挥着重要作用[1],植被覆盖变化监测已成为生态环境研究的重心[2-5]。陕北黄土高原地表破碎、沟壑纵横,生态环境脆弱[6],植被恢复对该区域的水土保持起着重要作用[7]。自1999年退耕还林生态建设工程实施以来,该区域的植被状况发生了很大变化,因此,开展陕北黄土高原植被覆盖的动态监测对黄土高原生态环境建设和退耕还林工程效果评价具有重要意义。

多位学者对陕北黄土高原植被覆盖变化及影响因子开展了系列研究,李登科等[6]利用GIMMS，SPOTVEGETATION两种数据对陕北黄土丘陵沟壑区不同时期的植被覆盖变化及其与气候的关系进行了分析;宋富强等[8]运用MODIS数据,从不同土地利用类型和不同坡度NDVI植被指数变化两方面分析了退耕还林对陕北黄土高原植被变化的影响;湛青青等[9]运用Landsat TM数据分别计算1990，2000，2010年的植被覆盖度,分析该区植被覆盖的空间变化情况。上述研究开展时必须在时空尺度选择上进行取舍:MODIS，SPOTVEGETATION等数据可以较好地表达植被覆盖的时序变化,但250～1 000m的低空间分辨率对于陕北黄土高原的复杂地形来说，在空间细节的表达上不够精确[10];Landsat等传感器获取的30m多光谱数据空间分辨率较高但有较长的重访周期,使得植被监测缺少时间上的连续性。时空数据融合综合了Landsat数据高空间分辨率和MODIS数据高时间分辨率的优点[11],Zhu等[12]学者提出了一种考虑像元反射率时间变化特征的增强型自适应时空融合模型ESTARFM(enhanced spatial and temporal adaptive reflectance fusion model),模型提高了对非均匀异质景观的预测精度,在复杂地形条件ESTARFM模型的预测结果与真实影像间的相关性最高[13]。

本文基于ESTARFM模型,采用Landsat和MODIS数据融合获得2008—2016年6～8月高时空分辨率的NDVI时序数据,并利用像元二分模型估算植被覆盖度,结合土地利用分类数据和气温、降水等气象数据,分析陕北黄土高原植被覆盖时空变化规律及影响因素,探讨时空分辨率融合影像在植被动态监测中的适用性,以期为该区域水土保持、环境保护等相关应用提供基础数据和科学参考。

1 研究区与数据

1.1 研究区概况

陕北黄土高原地处于陕西省北部,地理位置在35°21′～39°34′N,107°15′～111°14′E之间,海拔高度为800m～1 800m。东隔黄河与山西省毗邻,南边与陕西省渭南市、铜川市、咸阳市相接,西以子午岭为界与宁夏回族自治区、甘肃省接壤,北接内蒙古自治区[14]，如图1所示。气候上处在暖温带大陆性季风半湿润气候向温带半干旱气候的过渡地带[8],四季变化交替明显,年平均温度介于7.5～12.3℃之间,年降水量为350～660mm[15]。陕北黄土高原北部属于风沙高原区,南部属于黄土高原区,地表包括农田、森林、草地、荒漠等类型,地形复杂破碎,沟壑纵横[8]。区内植被类型多样,且分布差异较明显,主要分布有温带针叶林、落叶阔叶林、温带落叶灌丛、温带中生草甸、温带草原和暖温带草甸草原等植被类型,并呈现由东南向西北的落叶阔叶林向温带草原过渡的特征[16]。

图1 研究区地理位置示意图Fig.1 The geographical location of the study area

1.2 数据源与预处理

本研究采用的陕北黄土高原Landsat数据分别为2008—2012年的Landsat5 TM和2013—2016年的Landsat8 OLI影像(http://www.gscloud.cn)。MODIS数据采用500m空间分辨率8天合成的MOD09A1数据产品(http://www.nasa.gov),数据获取时间与Landsat数据相对照且间隔不超过3天。土地利用数据采用基于Landsat TM/OLI影像解译的2015年全国土地利用数据。气象数据来自寒区旱区科学数据中心(http://westdc.westgis.ac.cn)的SWAT模型中国大气同化驱动数据集(CMADS V1.0),空间分辨率为1/3°×1/3°的逐天格网数据,时间尺度为2008—2016年。

利用ENVI5.1软件对Landsat影像进行辐射定标、大气校正和几何校正。运用MRT(MODIS Reprojection Tool)对MOD09A1数据进行投影转换,重采样到Landsat的30m空间分辨率。对两种影像分别拼接后，用研究区矢量边界裁剪。由于影像数据来自不同的传感器,需调整影像的波段顺序相一致(表1)。为确保数据大小、影像栅格行列号等完全一致，裁剪时需在统一的图像处理软件下裁剪。

表1 Landsat和MODIS影像波段对照表Tab.1 Corresponding bands between Landsat and MODIS

2 研究方法

2.1 ESTARFM时空融合模型

ESTARFM时空融合模型旨在将Landsat影像与MODIS影像融合,从而得到具有高空间高时间分辨率特征的数据。对于指定区域,假设在同一日期由不同传感器获得的遥感数据可以相互比较,并在数据预处理后相互间有关联。ESTARFM的主要思想是利用关联来融合多源数据,同时最小化系统偏差。模型以预测像元为中心设置一定大小的滑动窗口,对窗口内像元利用权重函数进行卷积运算,得到中心像元预测值[11]。滑动窗口在整幅影像上逐一移动,从而得到需要预测的影像(式(1))[12]。

RM(xi,yi,tt,B))

(1)

式中:RL和RM分别为Landsat和MODIS影像数据;w为滑动窗口大小;(xw/2,yw/2)为待测像元的中心位置;(xiyi)为第i个相似像元的位置;tt为影像获取时间;B为影像波段;N为相似像元的数量;Wi是由空间、时间、光谱的距离共同决定的权重大小;Vi为转换系数[12]。选择ta和tb时的两对MODIS和Landsat数据分别通过公式(1)结合tp时的MODIS数据预测tp时的Landsat数据,结果记为RLa(xw/2,yw/2,tp,B)和RLb(xw/2,yw/2,tp,B)。结合这两种预测结果,预测的中心像元反射率会更准确。以越靠近预测时间则权重更大为准则,该权重计算为式(2)[12]。最后中心像元反射率的预测值为公式(3)[12]。

βt=

(t=a,b)

(2)

(3)

ESTARFM模型可以通过BI(Blend-then-Index)和IB(Index-then-Blend)两种方式实现[17],BI即将原始遥感影像输入到模型中,融合得到影像后再进行后续指数的计算[17];IB通过原始遥感影像计算研究所需的指数,再将计算好的指数输入到模型中,以模拟预测日期的指数[17]。由于运用ESTARFM模型时IB方式具有更高的准确性且计算时间较短[17],所以本文通过ESTARFM-IB方式进行归一化植被指数的模拟。即采用经过ta和tb时刻的两对MODIS NDVI和Landsat NDVI数据,结合tp时刻的MODIS NDVI数据预测同一时刻的高空间分辨率的Landsat NDVI数据。

2.2 植被覆盖度计算

经前期处理共获得陕北黄土高原2008—2016年6～8月共27期高时空分辨率NDVI数据,采用像元二分模型估算植被覆盖度,其计算公式如式(4)，(5)[18]。

NDVI=(NIR-R)/(NIR+R)，

(4)

Fc=(NDVI-NDVIsoil)/(NDVIveg-NDVIsoil)

(5)

式中:NIR为近红外波段反射率,R为红波段反射率;NDVIveg为纯植被覆盖时的NDVI值,NDVIsoil为纯土壤地表的NDVI值[19]。根据整幅影像的NDVI灰度分布频率累计表,选择5%和95%的置信度区间内上下限累计频率的值分别为NDVIsoil和NDVIveg。

3 结果与分析

3.1 ESTARFM模型融合结果与分析

由于Landsat数据存在云量大和数据缺失等问题,依据相关研究成果[10-13]和实际检验，本研究中采取预测日期前或后的两对MODIS和Landsat数据进行融合。例如2010年6月27日和2011年7月17日的两组MODIS NDVI和Landsat NDVI数据作为基础数据,通过ESTARFM-IB融合方式得到的2011年6月17日的ESTARFM NDVI数据,用2011年6月17日真实的Landsat NDVI数据作为参考为例,验证融合数据结果的准确和可靠性。

图2(a)～(b)显示整体上融合影像与真实影像的相似程度较高,NDVI像元值在空间分布上基本一致,清晰再现了30m空间分辨率的空间细节信息,数据纹理结构明显。相关性分析表明(图2(c)),ESTARFM NDVI与真实Landsat NDVI的散点大部分分布在y=x附近,斜率K为0.874 3,两者的相关性系数达到0.972 8(P<0.01)。因此融合影像能够较好的反映同期Landsat NDVI的影像光谱信息,可以作为后续计算研究区植被覆盖度的基础数据。

图2 预测NDVI与真实NDVI对比分析Fig.2 Comparison between predicted NDVI and real NDVI

3.2 陕北黄土高原植被覆盖时空变化特征

3.2.1 植被覆盖年际变化特征 计算陕北黄土高原2008—2016年每年6～8月的平均覆盖度并绘制其变化曲线图。图3显示,2008—2016陕北黄土高原整体植被覆盖变化呈显著增加的趋势,平均植被覆盖度由2008年的0.402增加到2016年的0.639,增速为0.031 2/a,增幅达58.9%,线性增长趋势呈显著水平。研究区植被覆盖曲线变化呈现波动变化上升的状态,但其年际间波动较大。表明研究区植被覆盖状况得到了明显改善,也显示出当地生态环境的脆弱性。

图3 2008—2016陕北黄土高原植被覆盖变化曲线Fig.3 The variation of annual FVC of Loess Plateau in Northern Shaanxi Province during 2008—2016

3.2.2 植被覆盖空间变化特征 根据陕北黄土高原的植被生长状况将植被覆盖度分为极低覆盖度(0～20%)、低覆盖度(20～40%)、中覆盖度(40～60%)、中高覆盖度(60～80%)、高覆盖度(80%～100%)5个等级[9]。2008，2012和2016年该区域的植被覆盖度空间分布如图4所示,呈现较为明显的区域差异,由东南向西北递减,高值区主要分布在陕北黄土高原东南部地区,低值区分布在西北部地区。2008年北部定边、靖边、横山、榆阳、神木、子洲、绥德等县区主要为极低、低、中覆盖度植被区域,南部甘泉、富县、洛川、黄陵、黄龙、宜川等县区为中高、高植被覆盖度。陕北黄土高原北部,植被稀少,稀树灌木丛草原较多;南部特别是桥山、黄龙山林区保有大量天然次生林,植被繁茂[2]。2012年全区极低植被覆盖度区域面积占比缩小,中覆盖度面积扩大,约占整个研究区的60%。到2016年中高覆盖度面积占比明显增多,主要增多的区域位于研究区中部、中东部和东北部各县区。陕北黄土高原植被覆盖度等级结构有所好转。

为研究近10年来陕北黄土高原植被覆盖变化趋势的空间差异,对每个像元的植被覆盖度与时间进行一元线性回归分析。图5显示陕北黄土高原大部分区域的植被覆盖度均表现增加趋势,东北部府谷、神木、佳县和中部米脂、绥德到中南部的清涧、延长等县的植被覆盖度多呈中度和明显增加的趋势,该区域经过近10年退耕还林等生态工程的持续实施,林地、草地等面积的快速增加使得植被恢复效果较好。南部的黄龙、黄陵、富县和中北部的定边、榆阳等县区的部分区域多呈现基本不变或略微减少的趋势。南部地区由于天然次生林较多,植被生长状况较为稳定,植被变化不明显;中北部部分地区属于长城沿线风沙区,在退耕还林工程实施中多采用封禁措施[8]，加之降水稀少使得植被的恢复效果相对较低。对比图4,植被覆盖增加趋势较为明显的区域多为中、中高植被覆盖度,基本不变或降低的区域多为高覆盖度。运用自然间断点分级法将变化斜率分为7类(表2),由统计结果看来2008—2016年研究区78%的土地面积植被覆盖度在逐年增加。

图4 2008，2012，2016年陕北黄土高原植被盖度空间分布图Fig.4 Spatial distribution of vegetation coverage of Loess Plateau in Northern Shaanxi Province in 2008, 2012 and 2016

图5 2008—2016陕北黄土高原植被覆盖度变化斜率分级图Fig.5 The classified changing slope of FVC of Loess Plateau in Northern Shaanxi Province(2008—2016)

Tab.2 Statistical results of vegetation cover variation trend of Loess Plateau in Northern Shaanxi Province(2008—2016)

植被覆盖变化斜率范围变化程度面积/km2百分比/% slope<-0.045明显减少1760.22 -0.045≤slope<-0.027中度减少3600.45 -0.027≤slope<-0.009略微减少1 6802.10 -0.009≤slope<0.009基本不变15 28019.1 0.009≤slope<0.027略微增加33 49541.8 0.027≤slope<0.045中度增加20 26425.3 0.045≤slope明显增加8 72010.9

图6 不同土地利用类型植被盖变化曲线Fig.6 The change of vegetation coverage trend of different types of land use types

3.2.3 不同土地利用类型植被覆盖变化 目前针对陕北黄土高原地区不同植被类型覆盖度变化特征的研究相对较少。本文对该区域的主要植被类型:有林地、灌木林、疏林地、其他林地、高覆盖度草地、中覆盖度草地和低覆盖度草地的植被覆盖变化进行分析可知(图6),2008—2016年陕北黄土高原7种主要植被类型6～8月的盖度均值介于0.372～0.744之间,依次为有林地>灌木林>高覆盖度草地>疏林地>其他林地>中覆盖度草地>低覆盖度草地,其年均数值分别为0.744 9，0.695 5，0.695 2，0.643 6，0.502 0，0.473 3，0.372 8。各类型平均盖度的曲线变化特点基本一致,均呈波动增加的趋势。计算各类型年均植被盖度与年份的相关系数可知灌木林和疏林地与年份的相关性最高(P<0.05),表明这两种类型的植被状况恢复较好,其余各类型与年份的相关性分析均没有通过显著性检验。有林地和高覆盖度草地的标准差较小,分别为0.028和0.021,表明在9a间这两种植被类型的覆盖变化较为稳定和平衡,没有明显的增加或减低的趋势;疏林地和其他林地的标准差较高,约为0.05,这两种类型属于植被分布稀疏的区域,主要由各种园地、苗圃构成,受人类活动影响而植被盖度年际波动较大。

3.3 气温和降水对植被覆盖变化的影响

植被的变化受自然和人为因素的共同影响,其中自然因素中气温和降水是对植被变化影响较大的两个因素[8]。对研究区71个站气象资料的分析表明(图7),2008—2016年该区6～8月平均气温呈微弱的下降趋势(slope=-0.021℃/a),6～8月降水量呈增加趋势(slope=12.85mm/a),同时气温和降水的年际变化较大。

图7 2008—2016年6～8月陕北黄土高原降水、气温的变化趋势Fig.7 Trends of precipitation and temperature of loess plateau in Northern Shaanxi Province from june to august(2008—2016)

本文从空间分布的角度进一步分析该区域6～8月平均植被盖度与同期气温、降水的关系,以像元为单位计算植被覆盖度与气温和降水的相关系数(图8,表3)。植被盖度与降水间呈显著正相关的区域主要分布在陕北黄土高原中部的定边至吴堡和东北部的佳县、府谷等县区(P<0.05),占全区总面积的10.1%。中部个别县区的部分地区相关系数呈极显著相关(P<0.01),占全区总面积的6.35%。中北部区域的主要植被类型为草地[8],在降水丰沛的时段长势较好,地表植被覆盖度就会增加,反之亦然。与降水相关系数为负的区域主要分布在研究区南部的黄陵、黄龙、宜川各县,其中0.33%的区域呈极显著负相关、1.79%呈显著负相关。这些区域主要是针阔混交林[8],南部林区有大面积的水源涵养林,这些森林植被短期内受降水的影响较小。植被盖度与气温呈显著和极显著正相关的区域主要分布在南部各县区,仅占全区总面积1.36%。植被盖度与气温呈显著和极显著负相关的区域极少,仅占0.78%,零星分布于北部地区。相比研究区植被盖度与同时期温度变化的关系,研究区植被盖度对降水因子的响应更敏感。

图8 植被覆盖度与降水、气温的相关性分析图Fig.8 Correlation coefficient maps of vegetation coverage with precipitation and temperature

Tab.3 The vegetation coverage correlation significance percentage with temperature and precipitation of Loess Plateau in Northern Shaanxi Province(2008—2016)

相关性植被覆盖度与降水/%植被覆盖度与气温/% 极显著负相关0.330.15 显著负相关1.790.63 不显著负相关25.3047.83 不显著正相关56.1450.04 显著正相关10.101.18 极显著正相关6.350.18

4 结 论

1)ESTARFM模型得到的黄土高原地区NDVI时序数据整体和局部细节上与真实数据相似性程度很高,获得了较高精度的光谱与空间信息,表明通过ESTARFM模型得到的高时空分辨率NDVI数据能够更好地应用于陕北黄土高原地形复杂地区的植被覆盖变化监测研究。

2)2008—2016年陕北黄土高原地区植被覆盖度整体呈显著增加的趋势,空间上呈现东南高西北低的特点。植被覆盖结构整体上好转,趋势分析表明近10年来研究区植被盖度整体呈现稳定的增加状态。该区域7种主要植被类型的盖度均呈波动增加趋势,灌木林和疏林地的增长趋势最为显著;有林地和高覆盖度草地的波动变化较为稳定,疏林地和其他林地的植被覆盖变化年际波动大。

3)2008—2016年陕北黄土高原6～8月的降水量呈上升趋势,而气温则呈现微弱的下降趋势。空间相关性分析表明,整体上6～8月的降水和气温对植被覆盖度的影响区域差异较为明显。植被覆盖度与降水呈显著相关的面积为18.57%,与气温呈显著相关的面积仅为2.04%,整体上植被盖度变化与气候因素的相关性较弱。黄土高原实施退耕还林还草工程以来,植被恢复状况较好,也使得植被覆盖对气温、降水等自然因子的影响敏感度减弱。

5 讨 论

在黄土高原地区植被覆盖变化的相关研究中, 多采用MODIS NDVI 250m/16d， SPOT VGT NDVI 1km/10d， GIMMS NDVI 8km/15d等低空间分辨率的时序数据。 其中MODIS， SPOT数据多用在2000年以后的黄土高原相关区域植被覆盖变化的研究中[2,8,23,25], 而GIMMS数据多用在较长时间尺度的研究[1,6,24,26]中。 上述研究数据在空间精度和时间连续性上并不能达到统一, 因而研究结论多偏向区域整体植被覆盖的变化或在研究时段的选择上有所限制。 而对于陕北黄土高原特殊复杂的地形特点来说,基于ESTARFM模型获取的高时空分辨率的遥感影像数据相对较好的解决了以上问题, 这种数据同时具备高度的空间细节表达力和快速时序变化能力, 为在陕北黄土高原开展高精度的植被动态监测研究提供了支撑。 同时研究发现连续的时空分辨率合成Landsat数据提高了MODIS等时序数据的空间分辨率, 对于陕北黄土高原复杂的地形特点来说, 可以更加精确的得到每一个栅格像元对应的NDVI数值, 确保了植被盖度趋势演变过程的可靠性。 ESTARFM模型为综合多源遥感数据在黄土丘陵沟壑地形复杂地区开展高精度的植被覆盖遥感动态监测提供了较优的数据源。本文针对研究区植被覆盖的整体时空变化、各等级植被覆盖变化和不同土地利用类型植被覆盖变化进行了探讨,今后在对高时空分辨率数据的挖掘和细化方面需要做进一步研究。