刘艳,聂磊,杨耘
1.中国气象局乌鲁木齐沙漠气象研究所,新疆 乌鲁木齐 830002;2. 中亚大气科学研究中心,新疆乌 鲁木齐 830002;3.武汉大学测绘遥感信息工程国家重点实验室,湖北 武汉 430079;4. 长安大学地质工程与测绘学院,陕西 西安 710054
草地是新疆畜牧业发展的基础,新疆现有草地面积5.12×107hm2,占全国草地面积的19.52%(刘新平等,2009)。新疆天山南北坡地理地貌、气候、生态差异巨大,牧草品种较多,作为新疆畜牧业的重要基地,承担着畜牧生产、生态平衡可持续发展的重任(许鹏,1993)。因此,定量评估天山南北坡草地覆盖变化对区域生态环境建设和牧业生产管理具有重要意义。随着遥感技术的进步,遥感监测可以提供植被生长相对即时连续的信息。其中,植被指数(Normalized Difference Vegetation Index,NDVI)可反映植被生长季和年际动态及气候条件变动对草地生态系统影响评价信息(Muukkonen et al.,2007;Piao et al.,2007;严建武等,2008;王新欣等,2008;Jin et al.,2011;Leisher et al.,2012;Ouyang et al.,2012)。NOAA/NDVI植被指数产品在空间和光谱分辨率上劣于MODIS/NDVI,相对获取不便,特别是近年来的植被指数产品(严建武等,2008)。MODIS/NDVI空间分辨率有250 m、500 m和1 km 3种,有单日、16日和月合成三类数据,产品时空分辨率能很好地反映草地植被时空变化和估算草地生物量 ,在森林及草地资源监测中有一定的有效性(Piao et al.,2007;王新欣等,2008;Jin et al.,2011;Ouyang et al.,2012;Leisher et al.,2012)。
目前,新疆大部分天然草地出现了不同程度的退化,草地生态环境破坏、畜牧业经济可持续发展受到影响。草地生态变化研究重点应集中在山区等生态脆弱区和自然地理过渡地带(夏浩铭等,2015)。但是,现有研究主要集中在天山北麓和伊犁河谷,主要针对植被时空格局及变化显著性水平开展(谢国辉等,2007;闫俊杰等,2013;刘芳等,2014;曹孟磊等,2016),很难了解整个天山地区及其典型草地空间分布特征和长时间序列变化规律及草地水热影响强弱的空间分布特征。因此,本文利用天山地区2001—2015年MODIS MOD13Q1数据集中的NDVI数据,定量分析天山地区不同草地类型 NDVI时空演变特征及其气候驱动时空特征,探讨不同类型草地水热影响强弱的空间特征,为天山地区草地生态环境评价保护提供科学依据。
根据1∶400万《中华人民共和国植被图》(中国科学院植物研究所,1979),将天山地区分为无植被、栽培植被和自然植被地段。其中,自然植被包含草丛、草原、草甸、荒漠、沼泽、高山稀疏植被、高山植被、灌从、阔叶林、针叶林和针叶阔叶混交林。研究区覆盖天山北坡西段伊犁河谷草原畜牧业区、天山南坡中段高山盆地草地限牧恢复区、天山北坡中段山地草原限牧恢复区、天山北坡东段山间盆地草原限牧恢复区和天山南坡东段-吐哈盆地草地禁牧保护区,涵盖高寒草原、草甸、温性草原化荒漠、沼泽等12种典型草地类型(图1)。研究区草地类型具有一定的典型性和完整性,空间分布上具有一定的连续性,是开展天山地区草地资源时空变化及水热影响研究的典型区域。
2.1.1 植被指数数据
植被指数数据为 MOD13Q1 NDVI(空间分辨率250 m,时间分辨率16 d),数据格式EOS-HDF,在全球正弦曲线投影SIN(Sinusoidal Projection)中编号h23v04、h23v05、h24v04和h24v05,时间为2000年6月—2015年12月,共1430景。利用MRT(MODIS Reprojection Tools)批量完成HDF-TIF格式转换、定义投影(WGS84)和影像拼接、最大值合成(Maximum Value Composition,MVC)2000—2015年NDVI年最大值以表征天山地区当年植被最好长势和生成7月NDVI均值。转换投影过程中,将原数据 ISIN(Integerized Sinusoidal)投影转换为 Albers,坐标系选取WGS84,像素大小250 m,校正类型采用最邻近像元法,输出类型为TIFF格式。
2.1.2 草地产草量数据
草地产草量地面调查集中在 2009—2015年 7—8月进行。在草地类型均一的典型区域设计了235块调查样地,样地基本特征调查中主要记录样地所隶属行政区、草地类型、地形、季节利用方式和利用状况等。一个样地内布设3个样方:草本及矮小灌木草原布设样方为1 m×1 m正方形,当样地草地分布呈斑块状或者较为稀疏时将样方扩大到 2~4 m2;具有灌木及高大草本植物且数量较多或分布较为均匀的样地,则布设1个10 m×10 m的正方形样方,也可为20 m×5 m的长方形样方,不做重复;齐地割取地上全部植株测量地上生物量,即为草地产草量。最终获取草地产草量采样点793个(图2)。利用GPS测定样方经纬度和海拔,同时在样方内采用常规植被调查法测定植物种数、植被盖度、群落平均高度等指标,获取植被盖度数据294个(图2)。
图1 研究区草地利用类型分区、典型草地类型及其位置示意图Fig. 1 Partitioning of the grassland use types, typical grassland types and their locations in the study area
图2 2009—2015年7月底—8月初草地总产量实测位置分布图Fig. 2 Spatial distribution of the in-situlocations of total herbage yield from the end of July to the beginning of August during the years 2009—2015
2.1.3 气象数据
气象数据为研究区内51个地面气象站2000—2015年日平均气温(℃)和降水量(mm)数据。已有研究指出植被生长和温度及降水量响应通常具有1~3个月的时间延滞(Piao et al.,2003;Chen et al.,2005)。因此,以生长季7月NDVI均值与前1个月、2个月及3个月≥10 ℃积温和累积降水量做空间相关分析,了解其是否具有时间延滞效应并以决定系数(r)来评估回归效果。利用ARCGIS-IDW(Inverse Distance Weighted)法将上述气象要素插值为分辨率、投影与NDVI一致的栅格数据,通过天山地区边界数据掩膜获得水热因子序列栅格数据。
2.2.1 植被指数处理
选取一年中草地最茂盛时期的NDVI值定量分析不同牧草类型变化趋势。由于气候影响,不同年份草地长势最茂盛的时间是不确定的。因此,对一年23景NDVI影像进行逐像元最大值提取,计算公式如下:
式中,MAX(NDVI,i)为第i年该像元最大NDVI,即是一年内该像元上草地植被最茂盛时期的NDVI;NDVI(i,j)为第i年第j景的NDVI值。
2.2.2 Sen and Mann-Kendall趋势分析
运用Sen and Man-Kendall趋势分析法模拟每个栅格单元的变化趋势,来综合反映不同草地类型时空格局演变(王佃来等,2013)。Sen趋势度β的计算方法为:
式中,n为时间序列长度(通常大于 10);Qi为时间序列中i时刻值。此处 n=15,Qi为对应年最大NDVI。。当β>0时,时间序列呈上升趋势,反之呈下降趋势。使用Man_Kendall趋势检验方法中的标准正态检验统计量 ZS对β趋势度进行显著性检验,计算方法为:
式中,S是Man_Kendall检验统计量;m是序列中重复出现数据组的个数;ti是第i组重复数据组中重复数据个数。当β>0且|Z|>1.96,序列呈显著上升趋势;当β>0且|Z|≤1.96,序列呈上升但不显著趋势;同理当β<0且|Z|>1.96时,序列呈显著下降趋势;当β<0且|Z|≤1.96序列呈下降但不显著趋势。
2.2.3 相关及显著性分析
对NDVI与气象栅格数据进行相关分析和F显著性检验,显著性仅代表趋势性变化可置信程度的高低,与因变量变化快慢无关。计算公式为:
式中,Rxy分别为 x、y两变量的相关系数;xi为第i年(月)的x要素值;yi为第i年(月)的y要素值;n为样本数。为要素x的平均值,为要素y的平均值。
式中,U为误差平方和;Q为回归平方和;yi为第i年(月)的实际要素值;Yi为其拟合回归值;为多年(月)平均值;n为样本数。
利用与准同期草地总产量及盖度估算数据拟合的方法对 MOD13Q1的 NDVI数据精度进行验证,以相关系数(r)来评估将其应用到天山地区草地资源时空演变的有效性。按照冷季放牧和暖季放牧,294个实测点比对结果显示(图3),冷季放牧形式下草地盖度和NDVI的相关系数r达0.897,暖季放牧形式下草地盖度和 NDVI相关系数 r达0.882,草地总产量与 NDVI的拟合决定系数R2分别为0.54和0.63,说明将NDVI数据应用到天山地区草地资源时空演变研究中是有效的。
为了消除季节变化带来的不确定性,对每年23景 NDVI逐像元最大值合成数据进行多年平均,生成2000—2015年NDVI年均值代表研究区统计年内草地平均覆盖状况。总体而言,NDVI≤0.2的地区占47.9%,0.7≤NDVI<1的地区占13%左右,0.4≤NDVI<0.7的地区占18%左右(图4)。天山南北坡东西跨度大,山体特征和在大气环流中所处的位置不同,承受西来水汽程度差别较大,草地植被的垂直带结构在西、中、东三段存在着明显的差异。整体而言,天山地区 NDVI表现为南部高于北部,西部高于东部,呈现由西南—中部—东北部逐渐降低的特征。按照天山南北坡看,天山北坡东段山间盆地草原限牧恢复区和天山南坡东段—吐哈盆地草地禁牧保护区整体NDVI偏低,NDVI集中在(0, 0.2],天山北坡中段山地草原限牧恢复区整体NDVI偏高,NDVI集中在(0.2, 0.7],天山北坡西段伊犁河谷草原畜牧业区NDVI偏高,集中在(0.7, 1],草地覆盖程度较高。天山南坡中段高山盆地草地限牧恢复区 NDVI偏高,集中在(0.2, 1](图 5)。
图3 实测草地盖度、总产草量和NDVI相关性分析Fig. 3 Actual measured grassland coverage, total herbage yield and NDVI correlation analyze
图4 天山地区2000—2015年年均NDVI等级图及统计面积(单位:km2)Fig. 4 Spatial distribution map and area statistics of annual averaged NDVI in Tianshan Mountain Area
图5 草地利用类型分区2000—2015年年均NDVI分类面积统计信息Fig. 5 Area statistics of the annual averaged NDVI of different grassland use types in 2000—2015 years
统计天山地区2001—2015年年最大NDVI多年均值和标准差,前者反映当年NDVI平均覆盖情况,后者反映相应的空间变异程度。由图6可知,NDVI均值维持在 0.32~0.35,表明天山地区草地覆盖多年来较稳定,NDVI均值和标准差随时间序列具有较大的波动,直观上无法进行趋势判断。NDVI均值序列曼-肯德尔(Mann-Kendal1)非参数统计检验(曹孟磊等,2016)结果也表明,统计年内区域NDVI增减趋势不显著,NDVI均值趋势斜率为正,呈现波浪式缓慢发展趋势,说明总体略呈增加态势。
图6 NDVI年均值和标准差趋势及NDVI年均值MK检验结果Fig. 6 The trend of annual mean of NDVI and standard deviation, and MK test results of averaged NDVI annually
结合研究区草地利用类型数据(引自新疆维吾尔自治区草原研究所)统计天山地区 12种典型草地面积和不同NDVI等级的面积比。结果表明,天山地区12种典型草地NDVI值域分布差异显著。温性草原化荒漠类、温性荒漠类NDVI均值较低,集中在[0.2, 0.4],属于较低覆盖度;低地草甸、高寒草甸类、高寒草原、温性草原类和沼泽NDVI均值较大,集中在[0.4, 0.7],属于中覆盖度;山地草甸、改良草地和热性灌草丛类NDVI均值较大,集中在[0.7, 1.0],属于高覆盖度;温性荒漠草原类和温性草甸草原类[0.4, 0.7],[0.7, 1.0]区间内像元基本均等,属于较高覆盖度(图7)。
从空间分布上看,43.82%的区域NDVI呈增加趋势,其中8.3%的区域显著增加,集中分布在天山北坡中段山地草原限牧恢复区;10.6%的区域NDVI呈显著下降趋势,主要分布在天山北坡西段—伊犁河谷草原畜牧业区。借用ARCGIS对其进行进一步量化分析以探究增长和退化具体出现在哪一种草地类型中。结果表明,山地草甸、温性草甸草原类和改良草地及热性灌草丛类下降趋势显著,温性草原化荒漠类、温性荒漠类、低地草甸和沼泽增加趋势显著,其他类型变化趋势不显著(图8)。
图7 基于NDVI年均状况分类图的不同草地NDVI分类面积比例统计图Fig. 7 A statistical chart of NDVI of area proportion of different herbage types based on category map of annually averaged NDVI
图8 NDVI年均状况发展趋势空间分布及统计面积和典型草地类型增减趋势面积统计Fig. 8 Spatial distribution map and area statistics of annual averaged NDVI and area statistics of development trend of typical grassland types
以7月NDVI均值为例,将5月、6月、7月及5—7月累积降水量与NDVI进行空间相关性分析,结果显示,NDVI与三者均具有一定的相关性,随着累积时段的增加,相关面积趋于增加,表明降水对植被生长的影响具有一定的滞后性,研究时段越长,越有利于准确辨识降水变化对植被生长的作用(图9)。以5—7月累积降水量为例,45.94%的区域NDVI与降水呈不显著相关,其中 26.8%的呈不显著正相关,46.65%的区域NDVI与累积降水呈显著正相关,其中只有3.2%的区域呈极显著正相关关系(图 9)。从草地类型看,温性荒漠类、温性草原化荒漠类、温性草原类和温性荒漠草原类与降水呈正相关,山地草甸类、改良草地与降水呈负相关,其他类型草地与降水相关性不显著(图10)。
以 3—8月≥10 ℃积温表征整个植被生长季节的热量状况,逐像元计算 2001—2015年 7月NDVI值与3—8月≥10 ℃积温线性相关系数。结果显示,约61.41%区域NDVI与 3—8月≥10 ℃积温呈不显著相关,其中24.42%的区域呈不显著正相关,30.84%的区域呈显著负相关,7.75%的区域呈显著正相关(图 11)。从草地类型看,温性荒漠类、温性草原化荒漠类、高寒草原和温性荒漠草原类与生长季热量状况呈正相关,山地草甸类、低地草甸、改良草地与生长季热量状况呈负相关,其他类型草地与生长季热量状况相关性不显著(图11)。
图9 不同累积时段降水量与NDVI空间相关分类面积比信息Fig. 9 Categorized area statistics of spatial correlation between NDVI and precipitation in different accumulated periods
图10 2001—2015年7月NDVI均值与不同时段累积降水相关性分布图Fig. 10 Correlation map between accumulated precipitation and NDVI mean of July in different periods during 2001—2015 years
图11 2001—2015年7月NDVI均值与同期3—8月≥10 ℃积温相关性分布图及典型草地相关性等级统计面积Fig. 11 Correlation distribution between the mean of NDVI in July and accumulated temperature more than 10 degree from March to August and area statistics of correlation level of typical herbage
以上研究结果与Heumann et al.(2005)、Seaquist et al.(2009)、杜加强等(2016)的研究结论一致,干旱和半干旱地区的植被物候因降水量的变动而发生改变,植被生长旺盛时期主要受到水分限制,区域尺度NDVI与降水量显著正相关。
本文利用MOD13Q1NDVI数据定量分析了天山地区典型草地2001—2015年时空特征和演变规律,对影响NDVI变化的主要水热因子进行量化分析。结果表明,天山地区不同草地年均 NDVI维持在0.32~0.35,呈波浪式缓慢发展趋势,天山地区草地状况总体发展稳定略呈增加趋势,南部高于北部,西部高于东部,呈现由西南—中部—东北部逐渐降低的特征,天山南北坡存在显著差异。影响NDVI变化主要水热因子量化分析显示,约 40.77%的区域与降水量呈显著正相关,30.84%的区域与3—8月≥10 ℃积温呈显著负相关,说明草地生态系统受气候环境因子的影响程度很大。
CHEN X Q, HU B, YU R. 2005. Spatial and temporal of phenological growing season and climate change impacts in temperate Eastern China[J]. Global Change Biology, 11(7): 1118-1130.
HEUMANN B W, SEAQUIST J W, EKLUNDH L, et al. 2007. AVHRR derived phenological change in the Sahel and Soudan, Africa, 1982—2005 [J]. Remote Sensing of Environment, 108(4): 385-392.
JIN Y X, XU B, YANG X C, et al. 2011. Remote sensing dynamic estimation of grass production in Xilinguole, Inner Mongolia [J].Scientia Sinica, 41(12): 1185-1195.
LEISHER C, HESS S, BOUCHER T M, et al. 2012. Measuring the impacts of community-based grasslands management in Mongolia's Gobi [J].Plos One, 7(2): e30991.
MUUKKONEN P, HEISKANEN J. 2007. Biomass estimation over a large area based on standwise forest inventory data and ASTER and MODIS satellite data: A possibility to verify carbon inventories [J]. Remote Sensing of Environment, 107(4): 617-624.
OUYANG W, HAO F H, SKIDMORE A K, et al. 2012. Integration of multi-sensor data to assess grassland dynamics in a Yellow River sub-watershed [J]. Ecological Indicators, 18: 163-170.
PIAO S L, FANG J Y, ZHOU L M, et al. 2003. Interannual varia-tions of monthly and seasonal normalized difference vegetation index (NDVI)in China from 1982 to 1999 [J]. Journal of Geophysical Research,108(D14): 4401.
PIAO S L, FANG J Y, ZHOU L M, et al. 2007. Changes in biomass carbon stocks in China's grasslands between 1982 and 1999 [J]. Global Biogeochemical Cycles, 21(2): B2002(1-10).
SEAQUIST J W, HICKLER T, EKLUNDH L. 2009. Disentangling the effects of climate and people on Sahel vegetation dynamics [J].Biogeosciences, 6(3): 469-477.
曹孟磊, 肖继东, 陈爱京, 等. 2016. 伊犁地区不同草地类型植被指数与气候因子的关系[J]. 沙漠与绿洲气象, 10(6): 73-80.
杜加强, 赵晨曦, 贾尔恒·阿哈提, 等. 2016. 近30 a新疆月NDVI动态变化及其驱动因子分析[J]. 农业工程学报, 32(5): 172-181.
刘芳, 张红旗, 董光龙. 2014. 伊犁河谷草地植被NDVI变化及其降水敏感性特征[J]. 资源科学, 36(8): 1724-1731.
刘新平, 吕晓. 2009. 新疆牧草地资源利用动态变化及其绩效分析[J].干旱区地理, 32(1): 81-85.
王佃来, 刘文萍, 黄心渊. 2013. 基于 Sen+Mann-Kendall的北京植被变化趋势分析[J]. 计算机工程与应用, 49(5): 13-17.
王新欣, 朱进忠, 范燕敏, 等. 2008. 利用EOS/MODIS植被指数建立草地估产模型的研究[J]. 新疆农业科学, 45(5): 843-846.
夏浩铭, 李爱农, 赵伟, 等. 2015. 2001—2010年秦岭森林物候时空变化遥感监测[J]. 地理科学进展, 34(10): 1297-1305.
谢国辉, 李晓东, 周立平, 等. 2007. 气候因子影响天山北坡植被指数时空分布研究[J]. 地球科学进展, 22(6): 618-624.
许鹏. 1993. 新疆草地资源及其利用[M]. 乌鲁木齐: 科技卫生出版社.
闫俊杰, 乔木, 周宏飞, 等. 2013. 基于MODIS/NDVI的新疆伊犁河谷植被变化[J]. 干旱区地理(汉文版), 36(3): 512-519.
严建武, 李春娥, 袁雷, 等. 2008. EOS-MODIS数据在草地资源监测中的应用进展综述[J]. 草业科学, 25(4): 1-9.