2000-2019年青海省同仁市NDVI时空动态变化

徐 佳, 高云飞, 刘伟伟, 党恬敏, 曾建青

(1.黄河水利委员会 黄河上中游管理局, 陕西 西安 710021；2.青海省水土保持中心, 青海 西宁 810000； 3.互助县水利局， 青海 海东 810500)

植被作为陆地生态系统的重要组成部分，受到土壤、气候和人类活动的综合影响，是生态环境质量的综合体现。2000年以后，大规模的退耕还林还草、封禁、坡改梯等生态工程的实施，使得黄土高原的植被覆盖显著提高，生态环境明显改善[1-4]。一方面，黄土高原植被建设取得了举世瞩目的成效，另一方面，局部地区的植被退化仍然存在。研究表明，黄土高原中部、东南部地区的植被恢复明显，而西北部地区的植被恢复却受到一定程度的抑制[5-6]。黄河源区高寒草地的退化趋势没有得到有效遏制[7]。城镇化建设导致大量农田和林地被侵占，植被面积减少[8]。由于地形地貌、水热条件、地理位置和植被类型的不同，植被的变化趋势可能存在很大的差异[9-10]。

归一化植被指数(normalized difference vegetation index, NDVI)是利用植物叶片在近红外波段具有高反射值，叶绿素在红光波段具有强吸收的特征，在多光谱遥感影像中，用近红外(IR)/红波段(R)，通过反映植被绿度变化来表征植被生长状况的指标，是目前公认的表征陆地植被覆盖度和生长状况的有效指标[11-13]，同时NDVI具有时空的连续性，因此在区域植被变化研究中得到了广泛应用[6,14-15]。目前最常用的NDVI 数据包括SPOTVGT NDVI[8]、AVHRR NDVI[16]和MODIS NDVI[17]。其中，MODIS-NDVI数据空间分辨率相对较高，具有时间序列较长、数据源较稳定的优势，适用于市、县或小流域尺度的植被变化研究[18-19]。目前的生态治理工程大多以县级行政区为单元进行实施，以县为单元研究植被变化情况，可以更为细致地分析区域生态环境治理效果，为下一步治理决策提供切实可行的依据。

青海省同仁市位于黄土高原最西端，区域内自然环境差异较大，生态环境脆弱，植被对人类活动和气候变化的响应敏感。作为隆务河省级水土流失重点治理区，近年来，同仁市实施了大面积的生态工程。本研究采用时间连续的MODIS遥感影像，对2000—2019年同仁市NDVI动态变化进行研究，分析NDVI与土地利用、海拔高度的关系，并对NDVI未来变化趋势进行预测，旨在客观评价同仁市20 a来的植被变化情况，为下一步生态治理提供参考，为黄土高原农牧交错区水土流失防治提供依据，助推黄河流域生态环境高质量发展。

1 数据与方法

1.1 研究区概况

同仁市地处青海省东南部(35°00′—35°47′N， 101°37′—102°27′E)，南北长85 km，东西宽75 km，全市面积3 465 km2。同仁市位于黄土高原最西端，与青藏高原接壤，属于农、牧业交错区，大陆性高原凉温、冷温半干旱气候，温度垂直变化明显，气温日差较大，光照充足，年均降水量401 mm，降水变率大，雨热同季，年平均气温5.6 ℃。

黄河一级支流隆务河纵贯全市南北，形成东西部山区和中部河谷地区。境内山峦起伏，沟谷相间，地势南高北低，地形可分为河谷川地、低山丘陵、中高山脑山、高山4个区域，海拔2 181～4 954 m。境内草场资源丰富，约占全市面积的70%。耕地沿隆务河河谷呈带状分布，多年来实施了大规模的梯田、封禁、造林、小流域综合治理等生态工程，生态环境得到显著改善。同仁市主要存在水力侵蚀和冻融侵蚀两种土壤侵蚀类型，水力侵蚀主要分布在中部隆务河谷地区；冻融侵蚀主要分布在东西部山区的高海拔地区。

1.2 数据来源

本研究采用的MODIS-NDVI数据通过NASA网站(https:∥search.earthdata.nasa.gov/)下载，分辨率250m。已有的研究表明，研究区NDVI在每年的7—8月达到峰值[20]，因此本研究采用最大合成法对2000—2019年7—8月NDVI数据进行最大化处理，消除云、大气等因素的干扰。海拔高度数据采用DEM数据提取，DEM数据通过地理空间数据云(http:∥www.gscloud.cn/)下载，分辨率为30m。土地利用数据采用2018年青海省水土流失动态监测成果，分辨率为2 m。

1.3 分析方法

1.3.1 趋势分析 趋势分析采用一元线性回归分析方法，对2000—2019年的20个 NDVI图层进行逐像元拟合，得到线性回归方程的斜率，即NDVI多年变化趋势Slope图层，表示多年植被覆盖变化的方向和速率[21]。若Slope>0，则NDVI呈现上升趋势；若Slope<0，则NDVI呈现下降趋势；Slope绝对值越大，表示NDVI值变化越快。Slope的计算公式为：

(1)

式中：Slope为NDVI与时间变量拟合的一元线性回归方程的斜率；i为时间变量，等于1到n的整数；n为研究时段的年数，取值20； NDVIi为第i年的NDVI值。

1.3.2F检验F检验[22]又叫做联合假设检验，假设回归方程很好地符合数据集要求，用来检验回归方程中被解释变量与解释变量之间线性关系在总体上是否显著。本研究采用F检验对NDVI变化趋势进行显著性分析，用于检验NDVI变化趋势置信度的高低。F检验法只代表趋势变化的置信度，与变化快慢程度无关。

1.3.3 高海拔和低海拔区域界定 野外调查发现，同仁市高海拔地区存在冻融侵蚀造成的植被退化现象，为了进一步分析冻融侵蚀对植被的影响，本研究结合冻融侵蚀区边界，将NDVI变化趋势图和冻融侵蚀区相结合，分析不同海拔区域NDVI变化情况。本研究采用第一次全国水利普查青藏高原区冻融侵蚀区边界的确定方法(公式2)[23]，利用经纬度计算冻融侵蚀区海拔下界，即该界限以上为冻融侵蚀发生的主要区域，本文将该界限以上定义为高海拔地区，该界限以下定义为低海拔地区。

(2)

式中：H为冻融侵蚀区下界的海拔(在此界限以上的区域，地表处于冻结或冻融交替的状态，为冻融侵蚀区,m)；X1为纬度(°)；X2为经度(°)。

1.3.4 Hurst指数 Hurst指数是基于重标极差(R/S)的分析方法，用来定量描述植被覆盖的长程依赖性[20,24]，基本原理是给定一个时间序列{ξ(t)},t=1,2,…，对于任意正整数τ=1，定义均值系列:

(3)

式中：{ξ(t)}为一个给定的时间序列，其中，t=1,2,…。τ为任意正整数，τ=1,2,…。

(4)

式中：X(t,τ)为累积离差，1≤t≤τ。

(5)

式中：R(τ)为极差，τ=1,2,…

(6)

式中：S(τ)为标准差，τ=1,2,…

(7)

式中：H为Hurst指数；yi为第i年的lnR/S值；n为研究时序。

若存在R/S∝τH，则说明时间序列{ξ(t)},t=1,2,…存在Hurst现象，H值称为Hurst指数，主要有3种形式(详见表1)。H值可在双对数坐标系(lnτ,lnR/S)中用最小二乘法拟合得到。

表1 Hurst指数存在形式

2 结果与分析

2.1 2000-2019年同仁市NDVI时间动态变化

图1为2000—2019年同仁市NDVI平均值年际变化情况，可以看出NDVI总体呈波动上升趋势，平均上升速率为0.027/10 a。2000年为研究时段内的最低值0.649 5，2019年上升到0.711 7。整个研究时段内NDVI最大值出现在2018年，为0.749 4。 2000—2019年NDVI变化大致可以分为4个阶段：①2000—2005年，NDVI值持续上升，由2000年的0.649 5，上升到2005年的0.718 7；②2005—2010年，NDVI先下降，然后波动上升，2010年上升至0.738 2；③2010—2015年，NDVI值先下降，后又上升至0.737 0，随后波动下降，之后缓慢持续下降，至2015年降至0.694 0；④2015—2019年，NDVI波动变化幅度较大，2018年NDVI为研究时段内最高值0.749 4。

图1 2000-2019年青海省同仁市NDVI平均值年际变化

近年来黄土高原NDVI整体呈现的波动上升趋势[1,5,25]与本研究结果基本一致，随着国家和整个社会对生态环境的日益重视，退耕、退牧、封禁等植被保护措施的广泛实施，黄土高原植被整体改善。由于黄土高原半干旱的气候特征，年际降水量的变化是造成NDVI波动的主要原因，降雨量偏少的年份NDVI值可能出现明显下降。

2.2 2000-2019年同仁市NDVI空间动态变化

图2为2000年和2019年同仁市NDVI空间分布图。与2000年相比，2019年NDVI分布图深色区域明显增加，浅色区域明显减少，即高覆盖度植被面积增加，低覆盖度面积减少。以NDVI值介于0.1～0.5的区域为例，2000年该区域为661.54 km2，占全市面积的19.09%；2019年该区域为350.21 km2，占全市面积10.11%。NDVI值介于0.8～1的区域2000年为388.63 km2，占全市面积的11.22%；2019年增长到1 066.92 km2，占全市面积的30.79%(表2)。从2000年到2019年，低覆盖度植被向高覆盖度转化明显。隆务河河谷地区NDVI增加主要归因于造林和小流域综合治理的广泛实施。东西部山区主要是在实施退牧和封禁措施后，植被自然恢复形成，东西部山区是高覆盖度植被的主要分布区。

表2 2000年和2019年同仁市NDVI分布区间面积和比例

图2 2000-2019年青海省同仁市NDVI空间分布特征

图3为2000—2019年同仁市NDVI变化趋势Slope的空间分布图。Slope为正值的区域，即NDVI上升的区域，为2 635.82 km2，占全市面积的84.42%(图3,表3)，广泛分布在隆务河谷地区和周围山区，其中林地626.13 km2，草地2 094.11 km2。

Slope为负值的区域，即NDVI下降的区域，为495.63 km2，占全市面积的15.58%。NDVI下降的区域有一小部分在隆务河河谷地区，是人类活动密集的区域，主要由于城镇建设造成；大部分则分布在西部和南部山区。NDVI下降的区域中林地105.58 km2，草地393.10 km2(表3)。由于草地是同仁市的主要植被类型，NDVI下降的区域植被类型也主要是草地，NDVI下降的草地占全市草地总面积的15.80%。建设用地中有24.21 km2的区域NDVI上升，而10.07 km2的区域NDVI下降。

图3 2000-2019年青海省同仁市NDVI变化趋势空间分布

表3 不同土地利用类型NDVI变化区域面积和比例

总体来看，同仁市NDVI在2000—2019年显著上升，Slope值主要分布在0～0.01之间(图4)，这部分区域占全市面积的81.94%，Slope值介于0.01～0.02的区域占全市面积的2.47%，介于0.02～0.03的区域占0.02%。

图4 2000-2019年青海省同仁市Slope值比例

采用F检验法进一步对同仁市NDVI变化趋势的显著性进行分析，结果表明，2000—2019年，同仁市NDVI极显著上升的区域占46.68%，显著上升的区域占16.20%；极显著下降的区域占2.42%，显著下降的区域占2.64%(图5)。可以看出，同仁市植被整体上显著改善，但仍然有部分区域植被存在退化趋势。

图5 2000-2019年青海省同仁市NDVI变化显著水平空间分布

2.3 NDVI变化和海拔高度关系

同仁市海拔高度在2 081～4 954 m之间，根据公式(2)计算得到同仁市冻融侵蚀下限海拔为3 583 m，因此本研究将同仁市海拔2 181～3 583 m的地区定为“低海拔地区”，将海拔3 583～4 954 m的地区定为“高海拔地区”(图6)，低海拔地区面积2 005.62 km2，高海拔地区面积1 459.38 km2。

将Slope栅格图和海拔高度图进行叠加分析(图7)，可以看出，在低海拔地区，NDVI整体呈上升趋势，NDVI上升的区域占92.16%，NDVI下降的区域仅占7.84%；在高海拔地区，NDVI下降趋势明显，NDVI上升的区域占73.71%，下降的区域占26.29%。从全市范围来看，NDVI上升的区域主要分布在低海拔地区，而NDVI下降的区域主要分布在高海拔地区。

图6 青海省同仁市低海拔和高海拔地区分布特征

图7 青海省同仁市NDVI变化和海拔关系分布特征

2.4 NDVI未来变化趋势预测

基于同仁市2000—2019年NDVI数据，采用Hurst指数对同仁市未来NDVI变化趋势进行预测(表4)。结果表明，同仁市NDVI的Hurst指数介于0.259 4～0.893 5之间，平均值为0.647 5。其中，Hurst指数>0.5的区域，即未来NDVI变化呈现正向演化特征的区域占92.28%，Hurst指数<0.5的区域，即未来NDVI变化呈现负向演化特征的区域占7.72%。

表4 2000-2019年青海省同仁市不同海拔地区NDVI变化面积和比例

将NDVI变化趋势Slope栅格图(图3)与Hurst指数栅格图进行叠加，得出NDVI未来变化趋势(图8)。可以看出，同仁市NDVI未来持续上升的区域占79.17%，持续下降的区域占13.13%，由下降转为上升的区域占2.45%，由上升转为下降的区域占5.27%。同仁市NDVI未来整体上持续上升，但仍有部分区域存在下降趋势。未来持续下降区域主要分布在西部和南部的高海拔地区，少部分分布在隆务河谷的城镇周边。

图8 青海省同仁市NDVI未来变化趋势

3 讨 论

2000—2019年同仁市NDVI整体显著上升，和2000年相比，2019年NDVI高值区明显增加，主要分布在东西部山区。NDVI上升的面积占全市的84.42%，广泛分布在隆务河谷地区和周围山区。近年来，同仁市实施的一系列生态工程促进了植被恢复，统计数据显示，截止2015年，同仁市实施草地禁牧面积693 km2，退牧、轮牧等草场保护措施促进了草地资源的恢复。在生产建设活动中，人们的生态环境保护意识也在不断提高，对于城镇区域的绿化也更加重视，同仁市建设用地中有24.21 km2的区域植被得到改善。人类生产、生活方式的改变促进了植被的恢复。

与此同时，同仁市NDVI呈下降趋势的面积占全市的15.58%，少部分分布在隆务河河谷地区，主要由于城镇化建设造成，大部分则分布在西部和南部山区，而且主要分布在高海拔地区，高海拔地区NDVI下降的面积占全市NDVI下降总面积的70.93%。作为高海拔地区水土流失加剧的重要原因，冻融侵蚀改变和破坏着土壤的物理性质，降低了农牧业生产能力[26]。研究人员针对青海省兴海盆地冻融侵蚀区植被的研究表明，冻融侵蚀会造成草被层呈斑状、鳞片状分布，草皮层在水平方向上撕裂，在垂直方向上与下面土层分离[27]。高海拔地区受季节性冻融作用的影响，土壤流失量的50%以上都发生在冻土层解冻时期[28]。调查发现，同仁市高海拔地区的坡面存在冻融侵蚀造成的草皮撕裂现象，草甸沿坡面向下滑动，造成土壤或基岩裸露，原本覆盖度很高的坡面大片变为裸地。对于高海拔地区植被退化还需要结合广泛的外业调查来进一步研究。

已有研究显示，多年冻土在黄土高原西北部缓慢退化，导致表层土壤水分减少，对植被产生不利影响[6]。寒冷干旱环境下，低温和干旱双重胁迫，导致植被对气候变化的响应更加复杂[29-30]。青藏高原及其附近的高山区是冻融侵蚀分布最集中且侵蚀最强烈的区域[31]。降雨的增加一方面可以改善植被，另一方面会加剧冻融侵蚀，造成植被破坏[27]。

同仁市NDVI未来变化趋势整体上持续上升，但仍有部分地区存在下降趋势。同仁市处于黄土高原和青藏高原的过渡地带，生态环境脆弱，植被对生态环境变化的响应敏感，受人类活动影响较大。NDVI存在下降趋势的区域是在科学研究和生态环境治理中需要重点关注的区域。气候变暖可能造成部分区域植被改善，也可能造成高海拔地区表层土壤水分减少、冻融侵蚀加剧，最终造成植被退化。目前对于冻融侵蚀的野外原地貌观测研究很少[32]。在未来的研究中，建议加强对高海拔地区植被变化和土壤侵蚀的野外定位观测，对NDVI下降区域进行深入分析，为制定合理的植被保护和恢复措施提供依据。

4 结 论

本文从时间尺度和空间尺度分析了青海省同仁市2000—2019年NDVI动态变化，初步分析了植被与土地利用、海拔高度的关系，对植被未来变化趋势进行了预测。

(1) 2000—2019年同仁市NDVI整体呈波动上升趋势，平均上升速率为0.027/10 a，2000年NDVI平均值为0.6495，2019年达到0.7117。NDVI高值区域增加明显，NDVI值介于0.8～1的区域由2000年的388.63 km2增加到2019年的1 066.92 km2，主要分布在东西部山区，主要为草地。

(2) 从全市范围来看，2000—2019年NDVI显著上升，NDVI上升区域为2 925.21 km2，占全市面积的84.42%，其中NDVI极显著上升的区域占46.68%，显著上升的区域占16.20%。NDVI上升的区域广泛分布在隆务河河谷地区和周围山区，林地NDVI上升的面积为626.13 km2，草地NDVI上升的面积为2 094.11 km2。

(3) 同仁市NDVI下降区域为539.79 km2，占全市面积的15.58%。NDVI下降的区域有少部分分布在隆务河河谷地区，主要为建设用地；大部分分布在西部和南部山区。林地和草地NDVI下降的面积分别为105.58 km2和393.10 km2。NDVI下降区域在海拔3 583 m以上和以下范围内的面积分别为383.62 km2和157.19 km2，分别占全市NDVI下降区域总面积的70.93%和29.07%，NDVI下降区域主要分布在高海拔地区。

(4) 同仁市NDVI未来整体上将持续上升，但仍有部分区域存在下降趋势。持续上升的区域占全市面积的79.17%，持续下降的区域占13.13%。由下降转为上升的区域占2.45%，由上升转为下降的区域占5.27%。