2000-2016年黄河源区植被NDVI变化趋势及影响因素

刘启兴, 董国涛, 景海涛, 周俊利, 党素珍, 贾培培

(1.河南理工大学 测绘与国土信息工程学院, 河南 焦作 454000；2.黄河水利委员会 黄河水利科学研究院 水利部 黄土高原水土保持流失过程与控制重点实验室, 郑州 450003)

植被作为土地覆被系统中的主要组分，是陆地生态系统存在的基础条件，也是连接土壤、大气、水分和人类土地利用的自然“纽带”[1]。植被是陆地表面能量交换过程、生物地球化学循环过程和水文循环过程中重要的下垫层，在土地利用、覆被变化、全球变化研究中起着“指示器”的作用。归一化植被指数(NDVI)被认为是反映植被生长状态及植被覆盖程度的最佳指示因子，是监测区域或全球植被和生态环境变化最有效指标[2-4]，监测植被动态变化以及分析其与气候关系的响应情况己经成为全球变化研究中的重要应用方向之一[5-6]。黄河源区位于青藏高原腹地，是黄河流域上游最重作要的产流区和水源涵养区，也是我国生态屏障的重要保障。因此，研究黄河源区植被覆盖的变化情况，对于整个黄河流域和青藏高原环境生态系统的演变具有重大意义[7-8]。

国内外学者利用NDVI数据来监测植被覆盖变化以及探寻植被覆盖与气象因子之间存在的响应关系。对NDVI影响因素的研究有以下3个思路：将提取的研究区域植被指数面平均值与气象要素进行简单相关分析或回归分析[9]；在像元尺度上进行简单相关分析[10]；选取研究区气象站点周围的NDVI值与对应气象站点气象数据进行简单相关分析和多元回归分析[11]。研究发现，植被覆盖与气象因子之间存在的响应关系受地理因素的影响[12]，不同植被类型NDVI与气温、降水的响应程度存在一定差异，但所有植被类型与气温、降水都呈现为正相关关系[13-14]，并且这种响应关系存在一定的滞后性[15-16]。刘宪峰等基于MODIS NDVI数据分析了2000—2011年黄河源区植被覆盖的时空变化特征及驱动因素，结果表明黄河源区植被呈增加趋势，植被覆盖的增加主要与气候的暖湿化及生态保护工程的实施有关[17]。李辉霞等探讨了三江源地区不同植被类型对气候变化的响应，通过分离气候变化与人类活动对植被NDVI的贡献定量评估了生态工程对植被变化的影响，认为气候要素对该区NDVI的贡献大于人类活动[18]。目前尚未有研究对黄河源区NDVI进行像元尺度的多元回归分析，缺乏定量化区分气象因素与人类活动对黄河源区内植被指数的影响的手段。开展NDVI像元尺度的分析，可以准确分析研究区域内每个地点的时空变化及气象因素对小范围的影响，空间连续性和异质性能够更好地体现和分析；利用多元回归分析，同时考虑多个气象因素对植被指数的影响，能准确分析气象因素对植被指数的影响及贡献，实现气象因素与人类活动对植被指数影响的定量区分，比较准确地基于像元尺度预测未来植被覆盖状况[19-21]。

本文利用M-K检验及趋势分析在像元尺度上对黄河源区2000—2016年植被NDVI进行时空变化分析，并利用相关分析法及复直线回归分析方法，定量分析气象因素和人为因素对植被动态变化的影响。

1 研究区概况

黄河源区(31.5°—36.5°N，95.5°—103.5°E)位于青藏高原东北部，面积约12.2万km2,横跨青海、四川、甘肃三省(图1)，绝大部分在海拔3 000～5 000 m以上，区域内分布有高山、盆地、峡谷、草甸、湖泊和冰川、冻土等地貌。由于地处高海拔寒冷地区，黄河源是对气候变化响应最敏感和生态环境最脆弱的地区。黄河源区自然环境类型多样，高寒植被分布广泛，有高寒草甸、高寒草原、高寒沼泽、高山稀疏植被、高寒灌丛、常绿针阔叶林等植被类型，其中主要以高寒草甸和高寒草原为主，其约占黄河源区总面积的70%以上。

图1 研究区位置及地形

2 数据源及分析方法

2.1 数据源

2.1.1 NDVI数据来源及处理 本文采用遥感数据MOD13Q1为NASA网站(https:∥ladsweb.modaps.eosdis.nasa.gov/search/)免费下载，时间范围为2000—2016年。MOD13Q1遥感影像数据是由美国对地观测计划EOS/Terra卫星携带的中分辨率成像光谱仪MODIS(Moderate Resolution Imaging spectro radiometer)获取的采用正弦曲线投影(SIN，sinusoidal Projection)方式的3级网格数据产品，具有1 km的空间分辨率和16 d的时间分辨率，单景影像覆盖面积为1 200×1 200 km2，数据采用格式为HDF—EOS。黄河源区地理范围覆盖全球正弦曲线投影SIN系统中编号为h25v05和h26v05的两个影像分带。

原始影像通过MRT(MODIS Reprojection Tools)工具统一进行投影拼接转换，转换为UTM投影，坐标系为WGS84坐标，并利用黄河源区矢量边界对处理后的数据进行裁剪。为了消除云层、大气与太阳高度角等干扰影响，将每年NDVI数据采用最大值合成法(MVC，maximum Value Composites)进行合成，获得从2000—2016年的17幅影像图。

2.1.2 气象数据来源及处理 气温和降水数据来源于中国气象科学数据共享服务网(http:∥data.cma.cn/)下载的全国气象站点数据。利用研究区内站点站号，通过IDL编程提取出黄河源区间内9个气象站点的月降水量和气温数据。由于气象因子存在空间不均匀性，同时为了逐像元地分析气象因素和人类活动等不确定因素对植被覆盖的影响，采用反距离权重插值法将研究区9个气象站点实测降水、气温数据插值为与MODIS NDVI数据集相同的空间分辨率栅格影像，获得每一个像元的降水和气温数据。采用反距离权重插值法将站点实测降水、气温数据插值为与MODIS NDVI数据集相同的空间分辨率栅格影像。

2.2 分析方法

2.2.1 Mann-Kendall检验法 Mann-Kendall检验法被广泛应用在植被等的趋势分析中。其优点在于检测范围宽、人为影响小、定量化程度高。该方法对于变化要素从一个相对稳定状态变化到另一个状态的变化检验非常有效[22]。利用Mann-Kendall检验法对黄河源区MODIS NDVI随时间序列变化的特征进行分析，研究结果不仅从统计学意义上揭示植被随时间的变化趋势，还揭示植被发生突变的时间段。在Mann-Kendall突变检验中，对于具有n个样本量的时间序列x，构造一个秩序列：

(1)

其中

(2)

可见，秩序列sk是第i时刻数值大于j时刻数值个数的累计数。

在时间序列随机独立的假定下，定义统计量：

(3)

式中：UF1=0；E(sk)，var(sk)是累计数sk的均值和方差，在x1,x2,…,xn相互独立，且有相同连续分布时，它们可由下式算出：

(4)

(5)

给定显著性水平α,将UFk和UBk两个统计量曲线和显著性水平线绘在同一个图上，若UFk和UBk的值大于0，则表明序列呈上升趋势，小于0则呈下降趋势。当超过临界直线时，表明上升或下降趋势显著，超过临界线的范围确定为突变的时间区域。如果UFk和UBk两条曲线出现交点，那么交点对应的时刻便是突变开始的时间。

2.2.2 趋势分析 趋势分析是通过对一组随时间变化的变量进行线性回归分析，从而预测其变化趋势的方法。该方法可计算研究区域监测时段内每个栅格点NDVI的变化趋势，反映该区域植被NDVI变化的方向和速率。线性回归方法是研究植被变化趋势的常用方法，其优点在于通过利用各个年份的数据值进行拟合，可以消除研究时段内偶发性异常因素对植被长势的影响，更真实地反映长时间序列植被变化趋势。

(6)

(7)

(8)

式中：xi表示研究时段内第i个年份；b为变化趋势线的斜率；a为截距；NDVIg为各个像元的NDVI；R为单个像元NDVIg变化率；x1和xn分别为研究时段的起始和终止年份。上述公式被广泛应用在NDVIg的时间序列分析，具有较好的稳定性[20-22]。

2.2.3 相关性分析 采用person相关性分析法与t检验对NDVI与同期气象数据进行相关分析和显著性检验，计算公式如下：

(9)

(10)

2.2.4 复直线回归分析 本文对NDVI与温度、降水之间进行复直线回归分析，回归方程为：

z=a+bx+cy

(11)

参数计算公式为：

(12)

(13)

(14)

NDVI模拟残差计算公式：

residual=NDVI-NDVIp

(15)

式中：residual为NDVI复直线回归方程的残差；NDVI为NDVI时间数据集；NDVIp是根据复直线回归模型预测的NDVI值。

3 结果与分析

3.1 黄河源区植被时空变化特征

3.1.1 黄河源区植被NDVI在时间上的动态变化分析 图2为利用matlab对黄河源区2000—2016年的MODIS NDVI进行MK趋势检验图。

图2 2000-2016年黄河源区NDVI非参数检验

黄河源区2000—2016年的植被NDVI趋势分析表明：UF统计量在2004—2016年大于零，说明黄河源区NDVI整体呈增加趋势；虽然在2002—2004年UF统计量小于零，并且UF线与UB线在2004年出现交点，说明在这期间NDVI的变化趋势发生较大变化，但交点统计量并未超过显著性α<0.05的临界值1.96，大体可以判断出黄河源区NDVI在2000—2016年期间整体比较平稳，没有特别明显的突变年份；从2010年起UF统计量超过了显著性α<0.05的临界值1.96，说明黄河源区NDVI在2010年后增长趋势显著。

3.1.2 黄河源区植被NDVI在空间上的动态变化分析 图3为源区植被NDVI年际变化趋势空间分布图，由图可以看到大部分区域呈浅色，说明2000—2016年黄河源区大部分区域NDVI是增加的，大部分区域的变化趋势为0～0.004/a，变化趋势最大达到0.213/a；NDVI在2000—2016年在减少的区域主要分布在扎陵湖和鄂陵湖西部与南部、黄河源区的中南部以及若尔盖湿地的北部。变化趋势最小值为-0.031/a。源区整体植被NDVI变化趋势均值为0.001 2/a。

图3 2000-2016年黄河源区NDVI年际变化趋势空间

表1 黄河源区NDVI一元线性回归分析变化

表1为黄河源区NDVI不同级别变化趋势面积统计表。植被NDVI年变化趋势<0的面积为3.54万km2，占黄河源区总面积的29.54%，其中变化趋势在-0.001～0之间的面积为1.68万km2，占黄河源区总面积的14.02%；植被NDVI年变化趋势>0的面积为10.31万km2，占黄河源区总面积的70.46%，其中变化趋势在0～0.004/a之间的面积为70 189.20 km2，占黄河源区总面积的58.55%。结果表明黄河源区在2000—2016年期间NDVI整体呈增加趋势，且变化趋势都在0～0.004/a，黄河源区部分地区植被覆盖情况正逐年好转。

3.1.3 2000年和2016年植被NDVI变化比较 黄河源区2016年植被NDVI年最大值影像减去2000年植被NDVI年最大值影像，得到两个年份NDVI差值空间分布图(图4)，并且做出不同差值范围面积比例统计表(表2)。2016年植被NDVI同2000年相比，大部分区域植被NDVI差值在-0.2～0，0～0.2这两个范围。说明近17 a来，黄河源区植被覆盖并未出现特别显著的变化。NDVI差值分布图中，NDVI减少地区域分布和NDVI年际变化趋势空间分布中变化趋势<0的区域基本相同。

2016年NDVI值较2000年减少的区域面积为57 362 km2，占黄河源区总面积的47.84%，其中NDVI差值在-0.2～0的区域面积为57 073 km2，占总面积的47.6%。NDVI差值大于0的区域面积为62 517 km2，占黄河源区总面积的52.15%；其中NDVI差值在0～0.2的区域面积为61 931 km2，占总面积的51.66%。结果表明2016年黄河源区植被覆盖较2000年基本保持不变有轻微改善。

图4 2016年与2000年NDVI差值分布

表2 2016年与2000年NDVI差值分布

3.2 黄河源区NDVI动态变化驱动因素定量分析

已有研究表明，植被NDVI对降水、气温的响应关系存在一定的滞后性。为了较为准确地定量分析植被NDVI对降水和气温的响应关系，本文在获得每个像元的降水和气温数据之后，分别求出黄河源区植被NDVI与历年降水、气温数据的6月、7月、8月、6月和7月均值、7月和8月均值及6—8月3个月均值的相关系数，并进行t分布检验。结果表明黄河源区NDVI与7月份的降水总量和6月份的平均气温相关系数最大，且均通过了α>0.05的置信度检验。所以本文选择了2000—2016年9个站点每年7月份的降水总量和6月份的平均气温作为降水和气温的研究序列。

图5 黄河源区2016年植被NDVI气候因素贡献分布

在复直线回归分析中，利用NDVI的2000—2016年时间序列数据和降水气温数据求出a，b，c参数，接着利用确定的NDVI与降水、气温之间的回归方程以及相应的降水和气温数据得出NDVI预测值。这个值即为气候因素对NDVI的贡献值，用实际NDVI减去NDVI预测值得到的残差即为人为因素和其他不确定因素对NDVI的影响。

图5为通过复直线回归分析得到的2016年黄河源区NDVI预测值即降雨和气温对NDVI贡献值空间分布，黄河源区东部和东南部区域NDVI预测值多在0.5以上，而在北部和西北部区域NDVI预测值较小。这是由于黄河源区属于大陆性高原气候，气温降水总体上呈现为由东南向西北逐渐递减的空间分布特征且北部和西北部地区靠近柴达木盆地，土壤水分不足，植被自我更新和自我生长能力较差，植被生长受降水、气温等气候条件限制更大。

表3是气温、降水对黄河源区2016年植被NDVI贡献值。99%以上的区域NDVI气候因素贡献值为正值，其中贡献值在0.5～0.6的面积为17 054 km2，占总面积的14.23%；贡献值在0.6～0.7的面积为21 651 km2，占总面积的18.06%；贡献值在0.7～0.8的面积为34 565 km2，占总面积的28.83%，三者总计面积为73 271 km2，占总面积的61.12%。这说明仅在气候因素的作用下，2016年黄河源区NDVI大部分在0.5～0.8的水平，而2016年NDVI实际观测值在0.5～0.6的面积为16 733 km2，占黄河源区总面积的13.96%；观测值在0.6～0.7的面积为21 449，占总面积的17.89%；观测值在0.7～0.8面积为34 443 km2，占总面积的28.73%，观测值和实际值的差别不大，可以看出目前气候因素仍是黄河源区植被覆盖的主要驱动因素。

为了定量化分析人类活动对黄河源区植被覆盖的影响，将实际观测的NDVI减去复直线回归分析预测的NDVI，即计算出残差，该残差即可认为是人类活动对植被NDVI影响的部分(图6和表4)。人类活动对植被NDVI起削减作用的区域面积为53 265 km2，占黄河源区总面积的44.43%；主要分布在黄河源区西北部的扎陵湖和鄂陵湖周边以及东南部的若尔盖盆地南部等地区。人类活动对黄河源区植被NDVI起促进作用的区域面积为66 613 km2，占黄河源区总面积的55.57%，其中贡献值在0～0.02的区域面积为64 190 km2，占总面积的53.54%。说明这些区域人类活动等因素促进了植被覆盖度增加。

表3 气温降水对黄河源区2016年NDVI贡献值

表4 人类活动对黄河源区NDVI贡献值

图6 人类活动对黄河源区2016年NDVI的贡献值空间分布

黄河源区地区政府自2003年起开展退牧还草和减畜工程，2015年平均牲畜存栏量减少至1 959.8万头，减幅达到8.60%[23-24]。同时2000年国家批准成立三江源国家级自然保护区和一期工程，对黄河源区等地的生态环境产生了积极影响[25]。但仍有接近50%的区域残差小于零，说明在全球气候变化和人类活动的影响下，对全球变化高度敏感的高寒生态系统已经开始加速退化，导致高寒草原覆盖度不断降低，景观破碎化加剧，草地水源涵养能力降低，水土流失加剧，这也与众多研究结果相吻合[26-28]。

4 讨论与结论

黄河源区超过70%的区域NDVI是增加趋势，增加速率在0～0.004/a。2016年同2000年相比，黄河源区较大部分地区的植被NDVI是增加的，但在扎陵湖、鄂陵湖西部和南部、黄河源区的中南部以及若尔盖湿地的北部等相当大的一部分地区，植被NDVI减少了0～0.2。

气候因素对植被生长变化起到了主要的促进作用，西部和西北地区贡献较小，西部和西北部地区由于临近柴达木盆地，土壤水分不足，植被自我更新和自我生长能力较差。植被覆盖与气候因子之间的响应作用是一个复杂的交互系统，本文目前只考虑了气候因素中的降水和气温，忽略了日照、风速、湿度等气候因子对植被覆盖的影响，深刻揭示其他各种相关因素对植被生长作用机理的研究还有待进行。

从复直线回归分析结果可以看出：人类活动对黄河源区整体植被覆盖产生一定显著的影响，源区高寒草地的退化趋势未得到有效遏制。黄河源区草场资源丰富，但可利用耕地面积少，是青藏高原地区重要的畜牧业基地之一。黄河源区在历史上很长的一段时间内一直处于严重超载过牧的状态。过牧是引起该区草地生态系统退化的主要因素。公路建设期间大量土石方的搬运工作造成的沿线50 m以内植被的破坏，仍然没有恢复，而且随着历次公路的改扩建活动，破坏的面积和范围在进一步扩大，从公路开始修建到历次的改扩建活动，由于其间隔时间很短，加之这一地区降雨、温度较低和植被的生长期短，自我恢复能力差，使得其沿线植被得不到很好的恢复，进而造成了水土流失、泥石流和滑坡等地质灾害的频发。因此要加强黄河源区生态保护力度，国家已经成立三江源自然保护区，并于2005年颁布《青海三江源自然保护区生态保护和建设总体规划》，采取了牧区产业结构调整、转变经济发展方式、牧民进城以及禁牧补贴等一系列有效措施，提高了群众的环保意识，大幅降低人为因素对长江黄河源区高寒草地破坏，延缓草地的退化。

农作物由于人为管理存在，导致NDVI值较大，存在对黄河源区NDVI分析结果产生影响的可能。通过分析2000年、2005年以及2010年三期黄河源区土地覆被数据发现，黄河源区耕地面积仅占黄河源区总面积的0.4%，因此，农作物对分析黄河源区植被NDVI的影响可以忽略不计。

通过植被NDVI动态分析可以有效地监测黄河源区封禁、退牧还草等政策落实和实施效果；且在现阶段人类活动变化幅度不大的情况下，可以根据各种气候模式预测的未来气温、降水数据预测黄河源区未来植被覆盖状况的理论值，以此为黄河源区生态保护和生态修复提供参考。