植被覆盖不仅是生态系统的主体部分^[1—3], 更是生态环境变化的重要指标^{[1, 4—5]}, 其变化是气候变化和人类活动共同作用的结果。其中, 人类活动对植被覆盖的影响成为研究全球变化的热点问题。20世纪中叶以来, 随着人类活动对自然环境干扰的不断加剧, 土地退化、水土流失等一系列生态环境问题频发, 植被覆盖出现一定程度退化。为有效解决区域的生态环境问题, 自1999年开始中国政府实施一系列生态恢复工程^[6]。人类通过保护现有植被、封山育林、退耕还林、人工造林等方式修复和重建被破坏的森林等自然生态系统。如今, 美国国家航空航天局(National Aeronautics and Space Administration, NASA)最新卫星资料显示：地球比20年前更绿了, 中国的植被增加量占过去近20年全球植被总增加量的25%^[7—8]。可以说, 植被覆盖变化是人类活动的深刻反映。

已有学者从国家或区域尺度出发^[9—11], 围绕人类活动和植被覆盖变化的关系做了大量的研究。赵安周等^[6]利用Hurst指数和非参数化趋势度(Sen)方法计算退耕还林(草)工程实施前后黄土高原植被覆盖的变化情况。程东亚等^[12]利用像元二分模型和相关分析等方法研究贵州沅江流域植被覆盖的时空变化特征, 并做了地形效应和人口效应的分析。李昊等^[13]在贵州省毕节地区, 建立NDVI—气候响应模型, 利用残差法分析退耕还林工程为主的人为因素在当地植被覆盖恢复中的作用。成佩昆等^[14]基于PCSE修正的面板数据模型, 研究了陕西省退耕还林工程对植被恢复的影响效应。在当前研究中, 大多着眼于黄土高原、青藏高原等热点研究区^[15—17], 对于西南地区植被覆盖变化的研究相对较少。虽然人类实施了多种植被覆盖恢复工程, 但目前研究多将目光集中在退耕还林方面, 对于封禁育林等其他工程对植被恢复作用的研究目前仍停留在定性的认识层面, 尚未见到相关定量研究的报道。

我国西南山地地形复杂, 生态环境脆弱, 植被恢复任务艰巨。重庆市作为成渝双城经济圈中的重要区域, 以及我国第三大林区(西南林区)的重要地区之一, 同时面对快速发展和生态保护的双重需求, 更兼具三峡库区重要生态战略地位。近20年, 随着城镇化高速发展、产业结构调整、人口流动等, 以及退耕还林、封禁育林等一系列生态修复工程的实施, 植被覆盖发生深刻演变。可以说, 重庆市是开展解析人类活动对于植被覆盖变化作用研究的典型试验区。因此, 本文基于重庆市2000—2020年归一化植被指数、气象数据、土地利用数据、封禁育林统计数据, 采用slope趋势分析和M-K检验分析重庆市植被覆盖时空变化特征; 利用GIS技术, 整合残差分析、土地利用转置矩阵等方法, 分析人类活动对植被覆盖变化的贡献程度, 定量评估封禁育林生态工程对植被恢复作用的大小, 为西南山地生态恢复工程政策提供科学依据。

1 研究区概况与数据来源 1.1 研究区概况

重庆市作为西南地区的门户城市, 地跨105°11′—110°11′E、28°10′—32°13′N, 幅员82402.95 km²。辖区内多山脉, 具有西南山地的显著特点。重庆市属亚热带季风性湿润气候, 年均温差小, 雨量丰沛, 太阳辐射弱, 日照时间短。重庆市是长江上游最大的经济中心、西南工商业重镇和水陆交通枢纽, 突出的经济实力和城市现代化进程对周边人口的吸引明显, 整体人口流动强度大。土地覆被类型以林地、农田和灌丛为主, 土地覆被类型以转换为林地和建设用地为主, 其中, 林地以常绿针叶林、常绿阔叶林和落叶阔叶林为主, 农田以耕地(水田和旱地)为主, 灌丛以阔叶灌丛为主^[18](图 1)。

1.2 数据来源与处理

本文所需的基础数据包括：遥感数据、土地利用数据以及气象数据等(表 1)。遥感数据来源于NASA EOS/MODIS, 选用2000—2020年归一化植被指数(Normalized Difference Vegetation Index, NDVI)500 m分辨率数据集, 使用最大值合成法生成逐月最大NDVI数据。2000—2020年降雨、气温、土地利用数据来源于国家地球系统科学数据中心, 精度分别为1000 m与30 m。数字高程模型(Digital Elevation Model, DEM)数据来源于中国地理空间数据云, 空间分辨率为90 m。封禁育林面积由重庆市各区县统计年鉴收集获取。

2 研究方法 2.1 趋势分析

利用线性回归对NDVI在时间序列上作趋势分析。公式如下^[19]：

(1)

式中, C_i为第i年的NDVI实际观测值, n为总年数。当slope>0时, 说明植被覆盖呈增加趋势; 当slope < 0时, 说明植被覆盖呈下降趋势; 当slope=0时, 说明植被覆盖不变。

线性回归的趋势分析显著性检验采用t检验法, 计算方法如下所示：

(2)

(3)

式中, 代表线性回归系数, 即所求slope; 其中y_i为第i年的NDVI实际值, 为第i年的NDVI的回归值; x_i为时间自变量; a表示显著性水平。

2.2 Mann-Kendall分析

Mann-Kendall方法是一种非参数统计检验的方法, 具有适用范围广、人为性少、定量化程度高的优点^[20], 在长时间序列的趋势检验和分析中得到广泛应用。假定时间序列(x₁, x₂, …, x_n) 是n个独立的、随机变量同分布的样本; 备选假设H₁是双边检验。对于所有的i＜j≤n, 且i≠j, X_i和X_j的分布是不相同的。定义检验统计量S:

(4)

式中, sign()是符号函数, 主要对括号中X_i-X_j的值进行判别, 设有以下规定：

(5)

当S为正态分布, 其均值为0, 其方差为：

(6)

M-K统计量公式在S>0, S=0, S <0是Z对应的值分别计算为：

(7)

式中, Z为正值代表增加趋势, 负值代表减少趋势。Z的绝对值大于等于1.28、1.96、2.32分别代表通过了可信度90%、95%、99%的显著性检验。

2.3 残差分析

残差分析的思路是通过剔除气候变化对植被覆盖的影响, 得到人类活动对植被覆盖变化影响的部分。本文对每个像元建立年均NDVI与气候因素(降雨和气温)的回归模型, 得到像元尺度上的年均NDVI预测值, 进而计算实测值与预测值之间的差值, 用来表征人类活动对植被覆盖变化的影响程度。其表达式如下：

3 结果分析 3.1 植被覆盖变化特征分析

2000—2020年重庆市植被覆盖呈现显著上升趋势, 上升速率为0.06/10a;M-K检验Z值为2.386, 置信水平2.32, 通过99%的显著性检验(图 2)。重庆市植被覆盖变化呈现出显著的时空分异特征, 大部分地区植被覆盖呈增加趋势, 最大增加趋势达到0.032, 主要集中于渝东北地区; 而最大减少趋势值达到-0.0214, 主要集中于渝中地区。重庆市近20年90%以上地区植被覆盖有所增加, 植被恢复状况良好。

根据表 2可知：植被覆盖显著增加区域所占面积最大(45.18%); 其次为变化不显著区域和极显著增加区域, 分别占30.53%和23.24%;极显著减少区域与显著减少区域最少, 仅占整个研究区的0.66%与0.39%。总体来看, 显著和极显著增加区域达到研究区面积的68.4%, 显著和极显著减小区域仅为1.05%, 研究区2000—2020年间植被覆盖增加趋势明显(图 3)。

3.2 人类活动对植被覆盖影响的时空特征

一般认为, 植被覆盖变化受到气候条件和人类活动的共同驱动作用^[21—23]。通过计算发现, 研究区2000—2020年年均NDVI与年均气温和年降雨量的相关系数, 分别为0.48和0.63, 均达到显著相关水平。本文采用残差分析法剔除气候条件作用, 得到人类活动对植被覆盖变化的影响程度(残差)。在栅格单元上采用降雨、气温因素预测NDVI值, 然后利用遥感观测的NDVI真实值减去NDVI预测值得到残差, 最后对残差做趋势分析, 得到重庆市残差年际变化趋势空间分布(图 4)。通过统计分析发现, 2000—2020年平均NDVI残差为0.005, 说明人类活动对植被覆盖的影响以正向为主。从时间变化上来看, NDVI残差呈逐年增加的趋势, 增长率为0.043/10a, 最大值出现在2020年(0.09)。从空间分布上看, 2000—2020年研究区人类活动对植被覆盖同时存在正向和负向两个方面的影响, 以正向影响为主导。人类活动对植被覆盖正向影响区域主要集中在研究区东南部和东北部, 梁平、忠县、万州和垫江等区县最为明显。人类活动对植被覆盖负向影响区域主要集中于人类活动密集、城市化进程较快的重庆市主城区以及各区县的城区(如潼南、永川、荣昌和大足等)。以上区域受到城市扩张影响, 大量耕地、草地和林地转变为建设用地造成植被覆盖的减少。

3.3 不同土地利用变化对植被覆盖变化的贡献程度

残差分析法采用残差表征每个栅格单元上人类活动对植被覆盖变化的影响程度, 但并未反映出人类活动对植被覆盖变化影响的具体方式。本文采用2000年和2020年土地利用变化、年均NDVI变化以及残差变化进行空间叠加分析, 得到不同土地利用类型变化区域栅格单元上实测NDVI变化值和残差变化值的总和, 并计算不同土地利用变化类型对植被覆盖变化贡献率的大小(表 3)。

据统计, 对植被覆盖变化贡献程度最大的为林地未变区(53.13%)、耕地未变区(31.05%), 这两个区域占植被覆盖变化总量的84.18%, 其次为建筑用地未变区、草地未变区和水域未变区。土地利用类型发生变化的区域对植被覆盖变化的贡献程度均很低, 其中最大的耕地转林地(退耕还林)贡献仅为1.24%。可以看出, 退耕还林工程对植被恢复的作用虽备受关注, 但对于研究区整体的植被覆盖变化而言, 其贡献程度并不显著。耕地是人类活动的产物, 其植被覆盖变化是人类种植行为作用的结果^[23—25], 因此耕地未变区的植被覆盖变化可以认为完全由人类活动造成。林地未变区植被覆盖变化显著(53.13%), 同时受到气候条件和人类活动的共同作用, 其中受到人类活动的作用程度(残差)达到47.33%, 值得深入探讨。

3.4 林地未变化区人类活动对植被恢复的作用

研究区林地未变区多集中在渝东北和渝东南地区, 这些地区多属于海拔较高的山地, 植被覆盖较高。提取林地未变区的年均残差, 可以看出：人类活动对林地未变区植被覆盖变化存在正向和负向两方面的作用, 以正向为主导。渝东北地区人类活动对林地植被普遍起到恢复作用; 渝东南及南部地区人类活动对林地植被覆盖的正向和负向作用并存, 并呈现相间分布格局。从区县尺度来看, 人类活动对林地植被覆盖负面影响在城口、丰都、彭水、武隆、石柱、酉阳和秀山均有分布; 人类活动对植被覆盖正向影响显著的区县包括巫山、奉节、巫溪、开州和江津等(图 5)。

封禁育林是重庆市实施的重要生态工程, 是林地未变区的标志性人类活动。基于此搜集重庆市各区县2000—2020年封禁育林面积数据, 选取其中数据获取较为全面的15个区县, 进行封禁育林面积、植被覆盖变化以及残差的相关性分析(图 6)。发现：林地未变区年均封禁育林面积与年均植被覆盖变化存在极显著正相关关系(R²=0.82)。说明, 封禁育林生态工程在促进植被覆盖恢复方面起到了重要作用。年均封禁育林面积与残差的极显著相关分析结果表明(R²=0.72), 残差能够较好的体现封禁育林生态工程对植被覆盖变化的作用, 在一定程度上印证了残差分析的可行性。

4 讨论

本文以重庆市为例分析了2000—2020年西南山地植被覆盖时空变化特征, 发现研究区植被覆盖整体处于上升趋势。这与刘灿等^[26]对重庆市2000—2011年植被覆盖变化的研究结论基本保持一致。近二十年, 西南山地植被覆盖呈现出较为明显的恢复态势, 人类活动作为导致全球植被覆盖变化的主要驱动力开始受到关注^[27—29]。土地利用变化是人类活动改变植被覆盖的重要表现形式^[30—32], 本文将残差和土地利用变化结合分析, 明确了不同土地利用变化类型对NDVI变化以及残差变化的相对贡献率。发现：土地利用变化对植被覆盖以及对残差的贡献率具有高度一致性, 并且林地未变区与耕地未变区对植被覆盖变化总量的影响程度最大, 达到84.18%。虽然已有研究已经从定量角度论证退耕还林、人工造林工程对植被覆盖恢复的显著效益^[33—34], 但在西南山地, 有必要对承载更多植被覆盖变化的土地利用未变区开展深入研究。耕地未变区植被覆盖变化可认为主要受到人类活动的驱动作用, 林地未变区则受到人类活动和气候条件的共同影响^[35—37]。本文基于封禁育林面积数据, 在区县尺度上建立年均封禁育林面积和林地未变区植被覆盖变化以及残差变化的回归方程。结果表明, 封禁育林对林地未变区植被恢复起到显著推动作用, 同时发现封禁育林与残差具有较高的拟合结果。因此, 可以认为林地未变区人类活动对植被恢复的作用(残差)与封禁育林工程密切相关, 残差变化的物理含义可以得到较好的阐释。本文整合多种分析方法对研究区植被覆盖变化的人类活动驱动作用进行了较为深入的解析, 首次定量揭示了封禁育林工程对植被覆盖恢复的重要作用, 对丰富区域植被覆盖变化驱动研究具有一定参考价值。

5 结论

本文以重庆市为例, 基于2000—2020年的气象数据、NDVI数据、封禁育林数据以及2000年与2020年土地利用数据, 利用slope趋势分析、回归分析、Mann-Kendall、残差分析和土地利用转置矩阵等方法, 分析了研究区植被覆盖变化趋势以及定量解析了人类活动作用。

具体结论如下：

(1) 研究区2000—2020年植被覆盖整体呈上升趋势, 植被覆盖变化空间差异显著。渝东南和渝东北地区植被覆盖较高, 且植被覆盖增加趋势明显; 渝西南中部地区植被覆盖较低, 且是植被覆盖减少的集中区域。

(2) 残差分析结果表明, 2000—2020年研究区人类活动对植被覆盖同时存在正向和负向两个方面的影响, 以正向影响为主导。人类活动对植被覆盖正向影响区域主要集中在研究区东南部和东北部; 人类活动对植被覆盖负向影响区域主要集中于重庆主城区, 以及各区县的城区范围。

(3) 不同土地利用类型转换中, 林地未变区对研究区植被覆盖变化总量贡献最大(53.13%); 其次为耕地未变区(31.05%), 耕地转林地(退耕还林)贡献仅为(1.24%)。封禁育林面积与林地未变区植被覆盖变化存在极显著正相关关系(R²=0.82), 封禁育林生态工程对植被恢复具有积极且重要的作用。