1998-2018年藏北申扎县植被NDVI时空变化及其影响因素
2021-12-09张晓克孟宏志
张晓克， 孟宏志
(1.河海大学公共管理学院，江苏 南京 211100；2.河海大学环境与社会研究中心，江苏 南京 211100)
增进民生福祉是我国经济社会发展的根本目标，而民生福祉与生态系统服务之间是一种相互依存的密切关系[1]，植被作为陆地自然生态系统的一个重要组成部分，其在交换能量、稳定气候和调节碳平衡等方面发挥了关键作用[2]，在生态系统建设中扮演着尤为重要的角色。归一化植被指数(Normalized difference vegetation index,NDVI)作为一个衡量地表植被覆盖度和初级生产力的强有力的良好指标，被广泛应用于检测区域植被动态变化及其对气候的反馈等研究之中[3]。在青藏高原的NDVI时空变化方面，国内学者已做了大量研究，卓嘎等[4]的研究结果表明2000—2016年间青藏高原NDVI值在不同研究时段呈增加趋势；张江等[5]通过对近30年来青藏高原高寒草地的NDVI值变化的研究，得出青藏高原高寒草地生长季的NDVI整体上是呈增长、改善的趋势，空间上有局部区域呈退化趋势。
目前国内外学者关于影响NDVI变化的驱动因素方面的研究，集中于NDVI与地形、气候和人类活动等驱动因素的相关性[6-10]。植被作为陆地自然生态系统的重要组成部分，其形成与变化在一定程度上与其所处的区域和地理条件有着密切的联系，地形条件是直接影响植被生境的关键因素，通过对高程、坡度、坡向等地形影响因子的研究可以掌握植被的空间变化规律，已成为揭示植被变化规律的重要研究方法[11-13]。刘梁美子等[14]通过对黔桂区域喀斯特山区年均NDVI与地形效应的研究，发现研究区年均NDVI随高程和坡度的增加呈单峰曲线状的变化趋势。刘晓婉等[15]在对雅鲁藏布江流域NDVI与高程因子的关联性进行研究时，发现研究区NDVI整体上呈现出随高程的增加而线性减小的变化趋势。研究生态系统对气候响应的时空模式尤为重要，植被动态变化及其对气候的响应是确定陆地生态系统气候变化机制的关键[16]。韩炳宏等[17]研究发现青藏高原NDVI与气候因子呈显著正相关，其中与气温的相关性是最高的。此外，许多学者的研究表明人类活动对植被的影响也十分深刻[6,18-19]，拉巴等[19]应用主成分分析的方法，找出1982—2014年藏北地区植被NDVI变化的主要的驱动因子分别是年均降水量、年均气温及年末牲畜存栏量这3个因子。
青藏高原是我国重要的自然资源库区和国家生态安全战略屏障，拥有丰富的水资源、草地资源和多样性生物资源，具有重要的生态系统服务功能，但在气候变化和人类压力下，生态系统变得较为脆弱[20-21]。近年来，国内外学者们大多基于全球、全国或较大尺度区域范围展开对植被生长变化的研究，其研究结果无法细化到各小范围区域内植被生长的具体变化趋势，也无法针对性地提供生态治理与改善的具体依据。本文选取位于我国藏北高原腹地南部的申扎县作为研究对象，申扎县人口密度较高、经济相对发达，具有较强的代表性[22]。基于1998—2018年申扎县NDVI、地形和气象等数据，本文采用线性回归、残差分析等研究方法分析申扎县NDVI年际时空变化特征及其与当地的地形、气候、人类活动的关联性，以期为申扎县及生态环境与其相似的周边县域的生态治理与改善提供科学依据，为评价植被变化状况及和其它生态环境相关领域的研究提供参考。
1 材料与方法
1.1 研究区概况
西藏自治区那曲市申扎县位于西藏中部，地理坐标为东经87°40′～89°48′，北纬30°10′～32°10′，平均海拔4 750 m，总面积25 546 km2。地处藏北高原腹地南部、冈底斯山和色林错之间，四周与双湖县、班戈县、尼玛县、念青唐古拉山接壤，县内包含盆地、丘陵和高山等多种地貌，设有色林错国家级自然保护区及申扎高寒草原与湿地生态系统观测试验站。申扎县属于高原亚寒带半干旱季风气候地区，其气候特征是寒冷干燥，年均气温在0℃左右，全年平均气温最低的月份是1月，1月的日平均气温为—9℃，年平均霜日数长达280 d，年八级以上大风日数多达104.3 d。高寒草原和高寒草甸是申扎县的主要植被类型，高寒草原以紫花针茅(Stipapurpurea)为主，草层高度为10～20 cm，植被盖度30%～40%，主要分布于海拔4 900 m以下区域；高寒草甸以矮生嵩草(Kobresiahumilis)为主，草层高度20～30 cm，植被盖度为70～90%，主要分布于海拔较高的山地阴坡；高山冰缘植被以风毛菊(Saussureajaponica(Thunb.) DC)、红景天(RhodiolaroseaL.)等为主，草层高度5～15 cm，植被盖度小于10%，主要分布于海拔5 300 m以上区域。全县总人口2万，其中乡村从业人员达1.2万，是一个纯牧业县，以饲养牦牛、绵羊和山羊为主[22-24]。
1.2 数据来源与预处理
本文所使用的1998—2018年申扎县年度和生长季NDVI数据来源于中国科学院资源环境科学数据中心(http://www.resdc.cn)发布的基于SPOT/VEG(1998—2018年)数据，采用最大值合成法得到的中国年度和生长季植被NDVI空间分布数据集，空间分辨率是1 km。另外，考虑到申扎县植被的生长情况，本研究将申扎县植被生长季规定为5—9月[23]。
高程数据来源于资源环境科学数据中心(http://www.resdc.cn)的30 m ASTER GDEM数据产品，本文对DEM数据进行校正和裁剪后得到申扎县高程分布数据，再利用ArcGIS 10.7中的“Slope”和“Aspect”函数得到坡度和坡向数据。将申扎县NDVI生长季均值按不同的地形因子进行平均，可得到不同高程、坡度和坡向的生长季平均NDVI变化。
气象数据来源于气象数据网(http://data.cma.cn)提供的应用空间插值法得到的全国年平均气温、年降水量数据集，其已对气象站点观测数据进行对比、拟合和修正，本文对其进行裁剪后得到申扎县气象数据[25]。
申扎县乡村从业人口数量和年末牲畜存栏量数据来自于西藏统计年鉴(2000—2019)各县(市、区)主要统计指标。其中，牲畜统计数据按照1头大牲畜等于5只羊单位的比例，统一换算为标准羊单位[26]。
1.3 研究方法
1.3.1线性倾向率 为定量分析1998—2018年申扎县NDVI值年际变化趋势，本文采用一元线性回归分析法，基于像元尺度，计算每个像元的NDVI与时间的回归斜率[27]。其计算公式为:
(1)
式(1)中：n为研究年数，本文中n=21；i为21年中的第i年；NDVIi为第i年的NDVI值；θslope为研究时间段内NDVI与时间回归趋势线的斜率，若θslope>0，说明该区域NDVI在这21年间总体上是呈增加的变化趋势，反之，则表示呈减少的变化倾向。
1.3.2相关系数 本文利用相关系数来研究NDVI变化与气候这一影响因子之间的相关程度，计算公式为：
(2)
1.3.3残差分析法 残差分析是一种通过计算实测值与预测值之间的误差来判断数据的可靠性，进而分析其他干扰因素的数据分析方法，它是由Evans和Geerken提出的[28-29]。本文将NDVI与气候因子构建的线性回归模型模拟出的NDVI值与实际的NDVI值进行比较，从时间序列中去除气候因子的影响，得到表征人类活动对NDVI变化的影响的残差值，通过残差值的大小来确定趋势的大小、方向和显著性[30-31]，如果残差趋势不显著，NDVI变化可以用气候变量解释，否则表明NDVI变化可能受到了人类活动的影响。计算公式为：
NDVI模拟=ax1+bx2+c残差=NDVI实测-NDVI模拟
(3)
式(3)中：NDVI模拟表示利用气候变量进行回归后预测出的NDVI模拟值，x1和x2分别表示气温和降水这两项气候因子。若残差>0，表明人类活动对植被生长变化产生正面影响，植被条件的改善可能归因于人类的保护和恢复工作；若残差<0，则表明该地区的人类活动对植被生长变化产生了负面影响。
当新赛季开打前一个多月，吴庆龙接过凯撒的教鞭时，整个夏训已经结束，球队的外援组合、内援引进没有丝毫变化余地。名义上出资人、俱乐部和主教练都换了，但球队残缺不全的人员结构和糟糕的备战状态却是从前任那里全盘接收过来的，西王和吴庆龙注定要为高速和凯撒还债。
1.3.4自然间断点分级法 自然间断点分级法是对一组数据内的相似数据值进行恰当分类，在差异相对较大的数据值处设置边界，使各类之间的差异最大化的一种分类方式[32]。
2 结果与分析
2.1 NDVI年际时空变化特征
基于1998—2018年申扎县年均NDVI数据，运用一元线性回归趋势分析法得到1998—2018年申扎县平均NDVI随时间变化的趋势图(图1)，由拟合出的方程式中斜率为0.004 1可知，1998—2018年申扎县NDVI整体呈波动性上升趋势，从NDVI值最小(0.239)的1998年起，21年间增长了37.23%，年均增长率为1.595%，最大值(0.343)出现在2017年。由r=0.828，P<0.01可知，1998—2018年申扎县NDVI与年份呈显著性正相关，说明在1998—2018年这一时段内申扎县NDVI值增加显著，生态环境质量得到有效改善。
图1 1998—2018申扎县年均NDVI值随时间变化趋势Fig.1 Trend of average NDVI in Xainza over time from 1998 to 2018
基于像元尺度，通过利用线性倾向率拟合出1998—2018年申扎县NDVI空间变化，根据自然间断点分级法和实际数据分布状况，确定|0.01|和|0.002|为分界点，将计算结果分为5个等级：显著减少(slope≤—0.01)、轻微减少(—0.01
2.2 NDVI随地形变化特征分析
2.2.1不同高程带的生长季NDVI变化特征分析 按一定的梯度将高程分为6个高程带(图3)，申扎县的高程总体上大致呈现出北低南高的分布趋势，观察图3a可知，申扎县高程范围在4 350～6 429 m之间，最低的高程带4 350～4 600 m主要分布在北部。由图3b可知，大部分区域(35.4%)处于4 600～4 900 m这一高程带上，处于5 800～6 429 m范围的区域面积最小，只占0.4%。4 900～5 200 m这一高程带的NDVI值最高(0.266 2)，NDVI最小值(0.004 6)则出现在5 800 m以上的高程带。在低于5 200 m的高程带上NDVI均值总体上是随高程的增加而增加的，在5 200 m以上的高程带上，NDVI均值则随高程的增加而减少，其中5 500 m以上的高程的降低趋势较为急剧。
图2 1998—2018年植被NDVI变化空间分布图Fig.2 Spatial distribution of NDVI changing trend from 1998 to 2018
图3 高程带分布(a)及不同高程带的面积占比及生长季平均NDVI变化(b)Fig.3 Elevation distribution(a) and different elevation range area and average NDVI of growing season(b)
2.2.2不同坡度的生长季NDVI变化特征分析 按一定的梯度将坡度分为6个坡度带(图4)，由图4a可以看出，申扎县的坡度范围在0°～62.4°之间，高低坡度相间分布，但总体来看，北部坡度较低，东西来看，西部的坡度要高于东部区域。观察图4b可知，有近65%的区域坡度低于10°，因此可以说大部分区域地势较缓。NDVI随坡度的增加呈现单峰曲线状的变化趋势，10°～15°这一坡度带的NDVI均值(0.232 1)最高，坡度大于25°的坡度带NDVI均值(0.116 4)最低，在坡度低于15°的坡度带上NDVI均值总体上大致是随坡度的增加而增加的，而在坡度大于15°的坡度带上，NDVI均值则随坡度的增加而减少，尤其在20°以上的坡度带的NDVI值随着坡度的增加迅速下降。
图4 坡度分布(a)及不同坡度带面积占比及生长季平均NDVI变化(b)Fig.4 Slope distribution(a) and different slope area and average NDVI of growing season(b)
2.2.3不同坡向的生长季NDVI变化特征分析 根据张建亮等[33]的研究，将申扎县坡向按一定度数划分为四类(图5)。由图5b知，阴坡的面积占比最大(36.71%)，半阳坡的面积占比最小(19.49%)。阳坡因光照资源较为丰富，较其它坡向而言能够更好地促进植物生长，因此其NDVI均值(0.220 6)最大，其次为半阳坡(0.214 4)、半阴坡(0.205 7)，阴坡的光照资源最少，其NDVI均值(0.177 1)最小。但总体来说，各坡向上的NDVI值差异不大，特别是半阴坡、半阳坡和阳坡。
图5 坡向分布(a)及不同坡向面积占比及生长季平均NDVI变化(b)Fig.5 Aspect distribution(a) and different aspect area and average NDVI of growing season(b)
2.3 NDVI年际变化影响机制
2.3.1NDVI随气候因子变化特征分析 申扎县在1998—2018年间的年均降水量范围在358.41～679.01 mm之间，年均气温在—3.70℃～1.39℃之间，通过分析研究区内的1998—2018年NDVI与降水、气温的年际相关关系(图6)及其显著性P检验，得出研究时段内NDVI变化趋势与降水变化呈正相关，但两者的相关性不显著(r=0.271，P=0.234)，与气温变化呈正向显著相关(r=0.556，P=0.009)。这说明申扎县NDVI均受降水、气温的正相关影响，且与气温的相关性较降水高，这与韩炳宏等[17]对2001—2018年青藏高原植被生长状况与气候因子相关性的研究结果相似，从图7知，气温每升高1℃，NDVI增加0.291；降水量每增加100 mm，NDVI增加0.224。
图6 年际NDVI与年均降水(a)和年均气温(b)对比分析Fig.6 Comparison analysis of interannual NDVI with precipitation(a) and temperature(b)
图7 NDVI与气象因子的关系Fig.7 Relationship between NDVI and climatic factors
2.3.2NDVI对人类活动的响应分析 本文采用多元线性回归的方法对NDVI预测值进行估算，以反映气候因素(包括气温和降水)的影响，从实际NDVI值和预测NDVI值的差异中获得残余量，再利用趋势分析对年际残差进行分析。通过对1998—2018年申扎县NDVI实际值与模拟值的对比得出研究时期内NDVI的残差值(图8)，从残差值的变化情况来看，2001—2002年残差转为正值后便持续9年为负值，直至2012年转为正值后便持续为正，并于2017年出现了残差最大值。总体而言，1998—2018年间NDVI残差值呈增长趋势，与时间变化呈正向趋势，两者呈显著相关(r=0.540，P<0.05)，这表明2012年后人类活动对植被生长产生了正向作用。由图9可知，年末牲畜存栏量自2011年起便开始总体呈现出下降的趋势，从2011年的102.61万只下降到2018年的88.95万只。
图8 申扎县1998—2018年NDVI残差随时间变化特征Fig.8 Characteristics of NDVI residual variation in Xainza from 1998 to 2018
为区分人为影响和气候控制，根据由气温和降水拟合的多元线性回归模型计算出申扎县NDVI残差趋势(图10a)，再根据自然间断点分级法，将残差趋势分为显著负残差趋势、不显著残差趋势和显著正残差趋势等3类(图10b)，间断值为—0.037 2和0.038 3。有显著残差变化趋势的区域意味着拟合的多元线性回归模型不能很好地解释这些地区的植被变化，这有可能归因于人类活动[34]，根据像元统计后得出，申扎县55.86%的区域的残差变化趋势不显著，NDVI与人类活动成正、负显著相关的区域占比分别为24.13%，20.01%，说明空间上，申扎县的人类活动的正作用总体上大于负作用。由图11可知，呈显著负残差的区域(20.01%)中超一半是位于4 900 m以下(10.68%)的较低海拔地域和坡度低于15°(13.63%)的缓坡区域，人类活动正作用明显的区域则主要分布在5 200 m以下的高程(22.47%)及坡度低于10°(18.31%)的区域。
图9 1999—2018年申扎县乡村从业人员数量和年末 牲畜存栏量Fig.9 Number of rural workers and year-end livestock stocks in Xainza from 1999 to 2018
图10 申扎县1998—2018年NDVI残差变化空间分布Fig.10 Spatial distribution of residual trend of NDVI in Xainza from 1998 to 2018
图11 各高程段(a)、坡度段(b)不同残差变化趋势占总面积的比例Fig.11 Proportions of different residual trends in the total area of each elevation(a) and slope section(b)
3 讨论
3.1 NDVI时空变化及其与地形、气候因素的关系
通过对1998—2018年申扎县NDVI时空数据分析，表明申扎县NDVI总体呈增长的趋势，生态环境质量得到有效改善。空间上，总体呈增长趋势的区域面积明显大于呈下降趋势区域面积，这与之前学者对青藏高原植被变化的研究结果相一致[4-5]。
申扎县NDVI随地形变化的特征与刘梁美子等[14]的研究结果相似，年均NDVI值是随高程和坡度的增加呈现出单峰曲线状的变化趋势，这种单峰曲线状变化趋势可能是由于较低高程带和缓坡地带的人类活动干扰较多，而海拔较高的地带和陡坡的生长环境较为恶劣，由此导致了中等高程带和坡度带NDVI值最高。而与刘晓婉等[15]的研究结果不同的原因可能是研究区域植被生长环境、研究时段、研究区面积等方面的差异。
研究时段内申扎县植被生长均受降水、气温的正相关影响，但与气温的相关性较降水高，这与以往学者的研究结果也是较为相似的[17,35]，说明申扎县植被生长对气温的敏感性较高。但需注意的是自上世纪五十年代以来，青藏高原一直呈快速增暖趋势，这可能会引起一系列水文与地质灾害[36]，今后可基于气候变化数据进一步探讨气候对陆地生态系统的综合影响，为统筹“山水林田湖草”系统建设提供指导。
3.2 人类活动对NDVI的影响
植被的年际生长变化不仅受气候因素的影响，而且还会受到人为活动等其他驱动因素的影响，通过利用残差趋势分析法对人类活动进行分析，结果得出1998—2018年间申扎县NDVI残差值呈增长趋势，与时间变化呈正向显著相关，从2012年开始，人类活动总体上对植被生长产生了正作用。
通过分析研究区内的1999—2018年NDVI与年末牲畜存栏量的年际相关关系及其显著性P检验，得出NDVI变化趋势与年末牲畜存栏量呈显著负相关(r=-0.531，P<0.05)，说明牲畜量的增加总体上对植被生长产生负作用。而人口数量主要通过退耕还草、禁牧轮牧休牧或过度放牧等方式间接影响草地NDVI，区分其不同行为方式产生的影响还需更加详尽的社会经济统计数据或野外跟踪调查数据，这也是当前间接验证残差分析结果的模型还有待进一步研究的原因[37]。申扎县NDVI与人类活动成负显著相关的区域主要是位于4 900 m以下的较低海拔地域和坡度低于15°的缓坡区域(图11)，这些区域可能是乡村从业人员的主要活动区或放牧集中区，随着申扎县地区经济的发展，1999—2018年间申扎县的乡村从业人员数量总体上呈现增长的趋势(图9)，由此带来不断加大的人类活动干扰也可能会对部分区域的植被生长产生负面影响。24.13%的区域NDVI与人类活动成正显著相关，主要分布在5 200 m以下的高程及坡度低于10°的区域(图11)，而由图9可知，申扎县的年末牲畜存栏量自2011年起便开始总体呈现出下降的趋势，这与近年来西藏地区“禁牧工程”和“草原生态保护补助奖励”等多个生态保护重点工程的有效实施密不可分[38]。
4 结论
1998—2018年申扎县植被NDVI总体上呈现出波动上升的趋势，年均增长率为1.595%，呈增长趋势的区域面积占比达71.89%。申扎县年均NDVI值随高程和坡度的升高呈现单峰曲线状的变化趋势，峰值分别位于4 900～5 200 m和10°～15°；阳坡的NDVI值最大，但总体上各坡向的NDVI差异不大。气候变化对申扎县NDVI整体起正作用，NDVI与气温的正相关性及显著性较降水高。从2012年起，申扎县人类活动主要对植被生长产生正面积极影响。整体而言，1998—2018年申扎县NDVI呈增长趋势，这与其暖湿化的气候变化和地方生态工程的有效实施密不可分。
