植被作为陆地生态系统的重要组成部分, 具有维持区域生态平衡, 促进全球物质循环和能量流动等作用, 其覆盖变化对区域生态环境状况具有重要指示意义^[1—3]。而归一化植被指数(Normal Difference Vegetation Index, NDVI)作为监测区域植被和生态环境变化的最佳指示因子^[4—5], 其数值可定量反映植被生长状况、植被生物量等信息, 因此在植被覆盖变化研究中得到广泛应用^[6—7]。重庆市地处我国西南山地, 具有差异明显的环境梯度和高度异质化的自然生境^[8], 是植被对环境变化极为敏感的区域。同时, 重庆市大部位于三峡库区腹心地带, 覆盖三峡库区85.6%的区域, 是长江上游重要的生态屏障。因此, 全面评价重庆市植被覆盖时空变化特征, 研究其与气候、地形、人类活动等因素各因子之间的作用机制, 对掌握近年来重庆市乃至整个西南山地植被生长状况, 认识区域生态环境变化, 进而揭示陆地生态系统演变规律等具有重要意义。

目前, 已有学者基于不同区域数据, 采用不同方法对植被覆盖变化及其驱动因素进行了大量研究。其中, 国外学者Park等^[9]对1997—2006年北纬50°以南的东亚地区植被覆盖进行了研究, 认为降水是影响植被覆盖变化的主要因子；Tucker等^[10]通过对1982—1999年北半球中纬度地区植被覆盖的研究, 指出植被覆盖变化与气温密切相关；Piao等^[2]通过对整个欧亚大陆长时间序列植被覆盖变化的研究, 发现气温和降水对植被覆盖变化具有交叉作用。国内学者栾金凯等^[11]利用多元回归方法分析了陕西榆林地区植被覆盖变化, 认为气象因素对植被生长与演化具有较强促进作用；张亮等^[12]基于偏相关分析法研究了长江流域植被覆盖与气象因子的相关性, 发现长江流域植被覆盖与气温的相关性较高；王建邦^[7]等运用趋势分析法研究了人为因素对中国植被覆盖的影响, 认为中国植被覆盖与人类活动整体呈负相关性。具体到重庆地区, 李学梅^[13]、朱林富^[14]、李惠敏^[15]等同国内其他学者类似, 主要通过回归分析、趋势分析和相关分析等方法, 探讨了气温、降水、地形等因素对重庆植被覆盖变化的影响。

综上, 目前研究植被覆盖变化及其驱动因素的成果已较为丰富^[16—18], 但主要集中在分析气候、地形等自然因子对植被覆盖变化的影响上, 探讨人类活动强度等人为因子对植被覆盖变化影响的则较少, 而综合研究气候、地形、人类活动等因素对植被覆盖变化驱动机制的则更少。鉴于此, 本研究基于MODIS-NDVI数据、气候数据、地形数据和人类活动数据, 采用趋势分析等方法, 应用地理探测器模型^[19—20], 全面分析2000—2019年重庆市植被覆盖时空变化特征, 并探究各地理因子及其交互作用对重庆市植被覆盖的影响程度与作用机制, 以期为重庆市及西南山地区域植被恢复、生态系统保护、生态文明建设等提供一定参考。

1 数据与方法 1.1 研究区概况

重庆市地处青藏高原与长江中下游平原的过渡地带, 地理范围为105°11′—110°11′E、28°10′—32°13′N, 幅员面积8.24×10⁴km²(图 1)。地形复杂多样, 以丘陵、山地为主, 地势起伏大, 有大巴山、巫山、武陵山、大娄山等山脉绵延境内。气候地域差异较大, 且立体气候明显, 年均温13—19℃, 年降水量914—1844mm, 年日照时数1000—1400h, 年均相对湿度高达70%, 具有冬暖春早, 夏热秋迟, 无霜期长等特点。境内江河纵横, 地表水资源丰富, 主要有长江、嘉陵江、乌江等, 其中长江干流自西向东横贯全境, 流程长达665km。土壤类型多样, 以紫色土、冲积土和水稻土为主。植被资源丰富, 类型多样, 2019年全市森林覆盖率达50.1%。

1.2 数据来源

MODIS NDVI数据来源于美国地质调查局(United States Geological Survey, USGS)提供的MOD13A3产品(http://glovis.usgs.gov/), 空间分辨率为1km×1km, 时间分辨率为30d。该产品经过几何纠正、辐射校正、大气校正等预处理；然后使用MRT(MODIS Reprojection Tools)工具进行投影和拼接, 经裁剪得到重庆市NDVI数据集；为进一步消除异常值的影响, 采用最大值合成法(Maximum Value Composites, MVC)^[21]合成月NDVI数据, 最终得到2000—2019年重庆市年NDVI数据。

气象数据包括气温、降水量和太阳辐射, 均来源于国家气象信息中心(http://data.cma.cn/)。其中气温、降水量通过逐站点统计重庆34个气象站点的日平均气温和日降水量, 最终得到年平均气温和年降水量数据；太阳辐射数据由于重庆地区的太阳辐射实测站点较少, 为提高数据精度, 利用目前模拟精度较高且应用广泛的Angstrom模型及日照时数数据, 模拟得到34个气象站点的年太阳辐射数据^[22]。最后, 利用ANUSPLIN气象插值软件的平滑样条函数对气温、降水量和太阳辐射数据进行空间插值^[23], 得到重庆市各气象因子栅格数据。

地形数据包括高程、坡度和坡向, 其中高程数据来源于地理空间数据云(http://www.gscloud.cn/), 而坡度、坡向数据则基于高程数据通过ArcGIS10.6的Spatial Analyst工具计算得到。

人类活动强度指数(Human Activity Intensity of Land Surface, HAILS)数据借鉴徐勇^[24]等学者构建的数学模型, 通过重庆市土地利用数据计算得到, 其中重庆市土地利用数据来源于中国科学院资源环境科学数据中心(http://www.resdc.cn/)。本研究中栅格数据空间分辨率均为1km, 其投影方式与NDVI数据一致。

1.3 研究方法 1.3.1 趋势分析

在植被覆盖时空变化分析中, 为研究重庆市2000—2019年植被覆盖变化趋势, 采用一元线性回归模型, 将NDVI与时间序列做回归分析^[25—26], 计算公式如下：

(1)

式中, θ为NDVI在多年时间序列上变化趋势的斜率, n为时间序列的年数, NDVI_i为第i年的NDVI值。若θ为正则说明植被覆盖随时间序列呈增加趋势, 反之呈减少趋势；并且θ数值的绝对值越大, 表明植被覆盖状况变化越明显。

1.3.2 偏差分析

为进一步研究重庆市植被覆盖在2000—2019年间相对变化(波动)程度, 通过偏差分析计算各年份NDVI偏离多年NDVI均值的大小, 以反映研究时段内植被覆盖的稳定性^[27], 计算公式如下：

(2)

式中, D_i表示第i年NDVI的偏差值, NDVI_i为第i年NDVI值, 为n年NDVI平均值。若D_i为正则说明第i年植被覆盖较多年均值高, 反之则低；并且D_i数值的绝对值越大, 说明偏离平均水平程度越大；同时, 通过分析连续多年NDVI偏差值, 可在一定程度上反映区域生态系统的稳定度。

1.3.3 地理探测器

在植被覆盖变化驱动因素分析中, 为定量研究气候、地形、人类活动等不同因素各地理因子(表 1)及其交互作用对植被覆盖变化的影响程度与作用机制, 本研究应用地理探测器模型^[19—20]中的因子探测模块和交互探测模块, 对各单因子作用强度及双因子交互作用强度进行探测。

其中, 因子探测主要用来探测NDVI的空间分异性和不同地理因子对NDVI空间分异性的解释程度, 其解释力大小用q值衡量, 表达式为：

(3)

式中, q为地理因子对NDVI的解释力, 值域为[0, 1], 值越大表示地理因子对NDVI解释力越强, 反之则越弱；h=1, 2, …, L, L为地理因子的分层数, N_h和N分别为地理因子层h和整个重庆市的样本数；σ_h²和σ²分别是地理因子层h和整个重庆市的NDVI值的方差。

交互探测用于识别不同地理因子之间的交互作用, 即评估地理因子共同作用时是否会增加或减弱对因变量NDVI的解释力, 或这些因子对NDVI的影响是否相互独立。通过对比单因子q值及交互作用q值即可判断交互作用的方式和方向。具体情况如下^[19—20]：q(X1∩X2) < Min(q(X1), q(X2)), 则交互后非线性减弱；Min(q(X1), q(X2)) < q(X1∩X2) < Max(q(X1), q(X2)), 则交互后单因子非线性减弱；q(X1∩X2)>Max(q(X1), q(X2)), 则交互后双因子增强；q(X1∩X2) =q(X1)+q(X2), 则交互后独立；q(X1∩X2)>q(X1)+q(X2), 则交互后非线性增强, 其中, X1和X2代表不同的地理因子。

地理探测器要求输入变量为分类变量, 即需对输入的连续型数据进行离散化处理。基于王劲峰等提出的数据离散化方法^[19], 本研究采用自然断点法对气象、地形、人类活动等因素各地理因子的连续型数据进行离散化处理, 将其分为10类。使用ArcGIS10.6中的Create Fishnet功能生成间隔为2km的渔网点数据, 提取栅格数据作为样本, 剔除异常值后, 共有19679条数据参与计算。

2 结果分析 2.1 植被覆盖变化特征 2.1.1 时间变化特征

趋势分析显示(图 2), 2000—2019年重庆市植被覆盖整体呈波动上升趋势, 且存在明显阶段性特征。2000—2019年重庆市植被覆盖状况正以4.4%/10a的增速逐步改善, 以2002年和2011年为节点, 植被覆盖变化明显表现为3个阶段：①前期(2000—2002年), NDVI年均值呈下降趋势, 下降速率为10.5%/10a, 作为对环境变化极为敏感的生态脆弱区, 受2000年和2001年连续特大干旱影响^[28], 重庆市植被覆盖状况出现明显下降；②中期(2002—2011年), NDVI年均值呈小幅增长趋势, 增速为2.9%/10a, 期间受三峡库区建成蓄水^[29]和极端高温事件^[25]等影响而出现轻微震荡, 但植被覆盖状况整体稳定；③后期(2011—2019年), NDVI年均值呈快速增长趋势, 增速高达7.8%/10a, 其间由于持续性的区域生态修复和生态文明建设等因素, 虽然不时受异常气候等因素的扰动^[26], 但植被覆盖增加明显。

偏差分析表明(图 2), 2000—2019年重庆市植被覆盖波动显著, NDVI偏差值呈现先减小后增加趋势。其中, 2000—2012年, 除2010年外其余年份NDVI偏差均为负值, 植被覆盖低于多年平均水平, 且2000年、2001年和2002年偏差值较大, 偏离多年平均水平程度较严重；2013—2019年, NDVI偏差均为正值, 植被覆盖高于多年平均水平, 且偏差值整体较大, 但各年份偏差数值差异明显。说明虽然近年来一系列生态工程的不断推进对改善区域植被覆盖状况起到了一定积极作用, 但生态系统依然脆弱, 还需进一步增强区域水源涵养、水土保持等生态功能, 逐步提升区域生态系统稳定性。

2.1.2 空间变化特征

从空间分布看(图 3), 2000—2019年重庆市植被覆盖空间异质性显著, 整体呈“东高西低”空间分布格局。植被覆盖多年均值介于0.26—0.88之间, 其中, 高值区(0.82—0.88)主要分布于渝东北的城口、巫溪、巫山等大巴山脉绵延区, 以及渝东南的武隆、丰都、石柱等武陵山边缘地区；低值区(0.26—0.58)则主要分部于受蓄水淹没影响的三峡库区消落带, 以及人口高度集聚、城镇化快速发展的渝西和主城地区。

从格局演化看(图 3), 2000—2019年重庆市植被覆盖整体呈改善趋势, 但高低集聚格局基本稳定。通过ArcGIS10.6软件中趋势面分析工具对重庆市植被覆盖格局演化进一步分析发现, 2000—2019年重庆市植被覆盖在东西方向上由“一字型”向“斜线型”演化, 东部地区植被覆盖增速大于西部地区, 东西部植被覆盖差距相对增大；南北方向上由“浅U型”向“深U型”演化, 北部和南部地区植被覆盖增速大于中部地区, 且高值区面积不断扩大, 低值区面积相对减少, 但植被覆盖高低集聚格局基本稳定。

2.2 驱动因素分析 2.2.1 因子探测分析

选取表征气候、地形和人类活动3类因素7个地理因子, 分别探测2000年、2005年、2010年、2015年、2019年各因子对植被覆盖的影响大小。同时, 以2000—2019年平均值从静态角度进一步分析各因子对植被覆盖的影响。从静态角度看(表 2), 不同地理因子对植被覆盖空间分异性解释力存在显著差异。各地理因子解释力大小依次为：年均温(0.3459)>高程(0.3281)>年均降水量(0.2305)>人类活动强度指数(0.1747)>坡度(0.1008)>总辐射(0.0552)>坡向(0.0034), 且均通过显著性检验(P < 0.05)。其中, 年均温、高程、年均降水量对植被覆盖空间分异性解释力较大, 是影响植被覆盖变化的主要因子, 人类活动强度指数、坡度次之, 总辐射、坡向最小。整体来看, 气候因素对植被覆盖空间分异性解释力最大, 解释力高达0.6316, 且通过气温、降水等因子直接影响植被覆盖；其次是地形因素, 解释力为0.4323, 通过在一定程度上影响降水、太阳辐射等气候因素的空间再分配过程及人类活动因素, 进而影响植被覆盖空间分布；人类活动因素则最小, 解释力为0.1747。

从动态角度看(图 4, 表 2), 各地理因子对植被覆盖空间分异性解释力表现出不同变化趋势, 但主导因子相对一致。2000—2019年年均温、坡向、人类活动强度指数解释力总体呈增加趋势, 年均降水量、总辐射、高程、坡度解释力总体呈减小趋势, 且年均温为第一主导因子, 高程和年均降水量分别为第二、第三主导因子。其中, 2000—2005年年均降水量、坡向解释力呈增加趋势, 其余各因子解释力呈减小趋势, 第一、第二、第三主导因子分别为年均温、高程、年均降水量；2005—2010年除坡向解释力呈减小趋势外, 其余各因子解释力均呈增加趋势, 且高程上升为第一主导因子, 年均温下降为第二主导因子, 年均降水量仍为第三主导因子；2010-2015年坡向、人类活动强度指数解释力呈增加趋势, 其余各因子解释力呈减小趋势, 人类活动强度指数上升为第三主导因子, 年均温和高程分别为第一、第二主导因子；2015—2019年高程、坡度、坡向、人类活动强度指数解释力呈增加趋势, 其余各因子解释力呈减小趋势, 高程、年均温、年均降水量分别为第一、第二、第三主导因子。

2.2.2 交互探测分析

同样, 在分别探测2000年、2005年、2010年、2015年、2019年各地理因子交互作用的同时, 进一步探测各地理因子平均值对植被覆盖变化的交互作用。交互探测结果表明(表 3), 2000—2019年重庆市各地理因子对植被覆盖变化存在交互作用, 且呈双因子增强和非线性增强两种类型, 说明各地理因子交互作用时会增强对植被覆盖变化的解释力, 而不存在相互独立作用或减弱对植被覆盖变化解释力的交互因子。其中, X₁∩X₂、X₃∩X₄、X₁∩X₄、X₄∩X₅、X₂∩X₃、X₄∩X₇、X₂∩X₄、X₂∩X₅、X₅∩X₇、X₂∩X₇在2000—2019年均为双因子增强型；X₃∩X₆、X₁∩X₆、X₄∩X₆、X₅∩X₆在2000—2019年均为非线性增强型；X₁∩X₅、X₁∩X₇、X₂∩X₆在2000—2019年主要为双因子增强型, 个别年份为非线性增强型；X₁∩X₃、X₃∩X₅、X₃∩X₇、X₆∩X₇在2000—2019年主要为非线性增强型, 个别年份为双因子增强型。

同时, 交互探测结果显示(表 3), 2000—2019年重庆市植被覆盖变化第一主导交互因子为年均温∩人类活动强度指数, 第二主导交互因子为高程∩人类活动强度指数, 第三主导交互因子为年均温∩总辐射。根据相关研究, 气温对植被生长具有重要影响作用, 气温低于或高于植被生长所能忍受的范围时, 均不利于植被生长, 甚至出现死亡^[20]；高程在一定程度上会影响人类活动因素对植被的作用大小, 同时重庆市作为典型的山地区域, 受气温垂直递减率影响, 在一定范围内随海拔升高而降雨量增加, 通过改变区域水热条件进而影响植被覆盖变化^[30]；人类活动对植被覆盖变化同样具有重要影响, 一方面通过退耕还林还草、生态保护修复等措施可促进植被覆盖, 另一方面通过城镇化建设、资源过度开发等行为可破坏植被覆盖^[30]；而太阳辐射则通过光合作用影响植被生长。总之, 2000—2019年重庆市各地理因子交互作用对植被覆盖变化的解释力不是独立作用或减弱的, 也不是简单的叠加过程, 而是解释力显著增强的双因子增强或非线性增强。

3 结论与讨论

本研究基于MODIS-NDVI数据、气候数据、地形数据和人类活动数据, 采用趋势分析等方法, 应用地理探测器模型, 全面分析2000—2019年重庆市植被覆盖时空变化特征, 并探究各地理因子及其交互作用对重庆市植被覆盖的影响程度与作用机制。主要结论如下：

(1) 时间变化分析显示, 2000—2019年重庆市植被覆盖整体呈波动上升趋势, NDVI偏差值呈先减小后增加趋势, 植被覆盖变化特征明显且表现为3个阶段。2000—2019年重庆市植被覆盖增长率为4.4%/10a, 前期(2000—2002年), NDVI年均值呈下降趋势；中期(2002—2011年), NDVI年均值呈小幅增长趋势；后期(2011—2019年), NDVI年均值呈快速增长趋势。

(2) 空间变化分析显示, 2000—2019年重庆市植被覆盖空间异质性显著, 但高低集聚格局基本稳定。2000—2019年重庆市植被覆盖呈“东高西低”空间分布格局, 高值区主要分布于渝东北的大巴山脉绵延区和渝东南的武陵山边缘地区；低值区则主要分部于三峡库区消落带及渝西、主城地区。植被覆盖空间格局在东西方向上由“一字型”向“斜线型”演化, 南北方向上由“浅U型”向“深U型”演化, 但植被覆盖高低集聚格局基本稳定。

(3) 因子探测发现, 不同地理因子对植被覆盖空间分异性解释力存在显著差异, 且变化趋势不同, 但主导因子相对一致。2000—2019年重庆市植被覆盖各地理因子解释力大小依次为：年均温(0.3459)>高程(0.3281)>年均降水量(0.2305)>人类活动强度指数(0.1747)>坡度(0.1008)>总辐射(0.0552)>坡向(0.0034)。其中, 年均温、坡向、人类活动强度指数解释力总体呈增加趋势, 年均降水量、总辐射、高程、坡度解释力总体呈减小趋势。

(4) 交互探测发现, 2000—2019年重庆市各地理因子对植被覆盖变化存在交互作用, 且呈双因子增强和非线性增强两种类型, 而不存在相互独立作用或对植被覆盖变化解释力减弱的交互因子。2000—2019年重庆市植被覆盖变化第一主导交互因子为年均温∩人类活动强度指数, 第二主导交互因子为高程∩人类活动强度指数, 第三主导交互因子为年均温∩总辐射。

本研究通过对2000—2019年重庆市植被覆盖时空变化特征及其驱动因素的分析, 在一定程度上探究了气象、地形和人类活动等因素各地理因子及其交互作用对重庆市植被覆盖的影响程度与作用机制, 但植被覆盖变化受多因素综合作用, 是一个非常复杂的动态过程, 因此, 在后续的研究中还需进一步完善指标因子, 同时考虑建立更长时间尺度的分析模型。此外, 今后的研究可在深入探究各因素对植被覆盖作用机制与驱动规律的基础上, 进行植被覆盖变化模型的构建和区域未来植被覆盖情况的预测。