植被作为自然生态系统的重要组成部分, 通过光合作用、呼吸作用与生物圈的其他自然要素形成紧密联系, 从而在生态系统的物质循环、能量流动与信息传输等方面发挥着重要作用^[1-2]。在植被生长过程中, 气温、降水、海拔与人类活动等环境因素具有重要影响, 植被对环境变化尤其是气候变化也最为敏感^[3]。因此, 研究植被动态变化及其与自然环境、人类活动因素的响应机制, 可以为全球变化提供重要的理论依据^[4], 对评价区域环境质量与维护生态平衡具有重要的现实意义^[5], 对人类社会经济建设与生态环境保护也具有相当的借鉴意义。

归一化植被指数(Normalized Difference Vegetation Index, NDVI)与植被覆盖度、生物量、叶面积指数以及土地利用等密切相关^[6], 能够精准地反映植被绿度、光合作用强度, 植被代谢强度及时间变化^[7], 在很大程度上展现区域植被覆盖度的高低以及植被生长状态, 被认为是监测陆地植被变化的最佳指示因子^[8], 已广泛应用于生态、环境和农业等领域^[9-12]。利用NDVI时序数据研究植被覆盖变化及对外界干扰响应, 对生态系统及其生态效益具有重要意义^[13], 已经成为全球变化研究的重要方向之一。

国内外学者基于长时间序列的NDVI数据集, 在不同时空尺度上对地表植被时空变化规律进行了深入研究, 主要集中于区域尺度NDVI变化趋势与气候因子之间的相互关系及对气候变化的响应机制^{[5, 14-16]}。部分学者从植物生长机制出发, 作了一些新的研究。如李辉霞等^[17]通过残差分析分离气候要素与人类活动对三江源NDVI的贡献, 定量评估生态保护建设工程成效。邓晨晖等^[18]基于NDVI数据探究了秦岭植被变化特征及对气候变化和人类活动的双重响应机制。徐芝英等^[19]应用小波分析研究了浙江省NDVI与高程、坡度和土地利用强度等因子之间的多尺度相关关系。Peng等^[20]从自然环境角度量化了四川省NDVI时空变化的自然因子驱动力。目前, 鲜有学者对影响植被变化的气候环境因素、非气候环境因素和人类活动因素进行整体驱动机制研究。因此, 本文选取长江上游地区喀斯特与丘陵地貌分布广泛, 地形条件较为复杂, 生态环境相对脆弱的宜宾市为典型案例, 利用2000—2015年MODIS MOD 13Q1数据, 采用Sen氏趋势分析与Manna-Kendall检验、Hurst指数分析宜宾市NDVI在年际尺度、季节尺度上的时间变化特征、变化趋势与未来可持续特征, 并应用空间转移矩阵和重心迁移模型分析其空间变化特征, 最后通过地理探测器模型探测影响NDVI的主导因素, 进一步揭示NDVI的空间分异特征及其驱动因子, 从而全面的分析宜宾市植被动态变化规律, 探讨生态建设工程成效, 为该地区乃至长江上游地区未来植被保护、水土保持与生态环境建设提供科学依据和决策参考。

1 研究区域概况

宜宾市位于四川省东南部(如图 1), 地处云贵川三省结合部, 金沙江、岷江、长江在此交汇, 川南丘陵与云贵高原过渡地带, 地跨27°50′N—29°16′N、103°36′E—105°20′E, 幅员面积13283 km²。境内地形整体呈东北低、西南高态势, 以中低山地和丘陵为主, 喀斯特地貌分布广泛。气候以亚热带湿润季风气候为主, 雨热充足、光照适宜、气候温润、四季分明。植被覆盖以亚热带次生性常绿阔叶林为主, 竹林等其他植被丰富多样。自2000年开始实行第一轮退耕还林工程以来, 境内森林覆盖率稳步增加, 2015年底全市森林覆盖率达44.08%。

2 数据与方法 2.1 数据来源与预处理

NDVI数据提取自美国国家航空航天局网站(http://reverb.echo.nasa.gov/reverb/)发布的MODIS MOD 13Q1陆地3级标准数据产品, 其空间分辨率为250 m, 时间分辨率为16 d。本文选取了2000年3月至2016年2月的MODIS NDVI数据共368期影像为数据源, 利用MRT、ENVI和ArcGIS等软件进行提取、剪裁和投影转换等批处理, 提取2000—2015年NDVI半月数据, 采用最大值合成法(Maximum Value Composite, MVC)消除云层、大气和太阳高度角的干扰^[12], 合成月NDVI数据。基于目视解译发现冬季NDVI数据存在云层较厚, 数据质量较差的情况, 应用Savitzky-Golay滤波法对其进行滤波处理。为了充分反映NDVI变化的年际和年内季节特征, 根据宜宾市植被物候特点, 进行春季(3—5月)、夏季(6—8月)、秋季(9—11月)和冬季(12—次年2月)的划分, 然后提取平均值生成2000—2015年各季度和年度NDVI数据, 以指示各季度和年际变化。

用于精度验证的3景Landsat 8-OLI数据来源于地理空间数据云(http://www.gscloud.cn/), 其轨道号分别为128/040、128/041、129/040, 经过辐射定标、大气校正、镶嵌裁剪等预处理后, 进行NDVI计算及异常值去除, 最终得到宜宾市30 m分辨率NDVI数据。

气候、土壤、DEM数据和2015年人口密度、土地利用类型、GDP数据来源于中国科学院资源环境科学数据中心(http://www.resdc.cn/), 通过ArcGIS软件进行裁剪、投影变换、重采样等处理。其中, DEM数据原始空间分辨率为90 m, 经过重采样处理后, 在250 m的DEM数据基础上提取宜宾市海拔、坡度和坡向数据。

2.2 研究方法 2.2.1 Sen氏趋势分析与MK检验

与传统的最小二乘法进行线性回归趋势分析相比, Sen氏斜率趋势分析可以有效避免时间序列数据缺失和数据分布形态的影响, 并且可以消除异常值对时间序列的干扰^[21]。Sen氏斜率计算公式如下：

式中, 1 < i < j < n, i、j为时间序列数, NDVI_i、NDVI_j分别为i、j时间序列的NDVI数值。当斜率Q大于0时则表示存在上升趋势, 小于0表示存在下降趋势。

Manna-Kendall检验(MK检验)不需要样本数据遵从特定的分布, 能够有效的剔除异常值, 适用于非正态分布的数据^[22]。所以, 本文采用Sen氏趋势分析与MK检验相结合, 分析NDVI在像元尺度上的变化趋势及显著性；在区域尺度上则采用线性回归分析NDVI变化的总体趋势。

2.2.2 Hurst指数

自然界中普遍存在具有长程依赖性的时间序列, 如水文、地质和气候等, 估计Hurst指数是定量描述长程依赖性的主要方法之一^[23-24]。关于Hurst指数H的获取可以选用多种方法, 常用的是R/S分析法^[25]。其基本原理如下^[26]：

对于给定的时间序列{NDVI_(t)}, t=1, 2, …, n, 定义均值序列：

累积离差为：

极差为：

标准差为：

则根据以上公式, 求比值R_(T)/S_(T)≌R/S, 若R/S∝T^H, 说明在分析的时间序列内存在Hurst现象, Hurst指数H通过log(R/S)_n=a+H ×log(n)利用最小二乘法拟合得到^[2]。如果0.5 < H < 1, 表明时间序列内NDVI具有长程正相关性, H越接近1, 持续性越强；如果H = 0.5, 表明NDVI变化趋势在时间序列内是随机序列, 不存在变化持续性；如果0 < H < 0.5, 表明NDVI趋势在时间序列内具有反持续性, 即未来趋势与过去趋势相反, H越接近于0, 反持续性越强^[26]。

2.2.3 空间转移矩阵与重心迁移模型

空间转移矩阵来源于系统分析中对系统状态与状态转移的定量描述^[27], 可以定量识别某一要素不同等级在某一时间间隔的空间格局变化, 除了反映其不同等级面积变化, 还能直观反映各等级面积转入转出情况^[28], 因而广泛应用于土地利用、植被覆盖等方面。

重心迁移模型可以反映某一要素在空间演变过程中的时空聚集和迁移特征。NDVI所表征的植被覆盖空间变化可以通过不同植被覆盖类型重心迁移过程来描述其总体变化趋势和空间变化特征^[29]。计算第t年植被覆盖的重心公式为：

式中, X_t、Y_t分别表示不同类型植被覆盖分布重心的经纬度坐标；C_it表示第i个植被覆盖斑块第t年的面积；X_it、Y_it分别表示第i个植被覆盖斑块第t年重心的经纬度坐标。

2.2.4 地理探测器

地理探测器是基于“因子力”度量指标, 结合GIS空间叠加技术和集合论提出的, 用于探测空间分异及分异机制背后驱动力的计量模型^[30], 包括风险探测、因子探测、生态探测和交互探测4部分。传统的植被覆盖变化驱动力研究一般采用相关分析方法, 计算和检验影响因子与植被覆盖之间的相关关系。当面临多个影响因子时, 其计算较为繁琐, 而地理探测器能够较好的克服这一问题^[31]。其因子探测模型如下：

式中, L为因变量NDVI或影响因子X的分层, 即分类或分区；N_h和σ_h²分别为层h的单元数和方差；N和σ²分别为研究区域整体单元数和方差。q的值域为[0, 1], 值越大表示因子X对NDVI空间分异的解释力越强, 反之则越弱^[32]。

交互探测用于识别不同影响因子X_S之间的交互作用, 即对因子X₁和X₂共同作用时对因变量NDVI解释力的影响进行评估。这是地理探测器区别于其他统计方法的最大优势之一^[33]。生态探测利用F统计量, 分析两影响因子X₁与X₂对因变量NDVI空间分异的影响是否存在显著差异。如果存在显著差异, 记为“Y”, 否则记为“N”, 据此可以判断哪一因子对NDVI更具影响力。风险探测是用于判断两个子区域间的属性均值是否存在显著差异, 均值显著性越大的区域, 其因变量值越大, 据此可以搜寻影响因变量NDVI的因子适宜区域。

2.3 精度验证

基于遥感数据不确定性, 本文通过同期高分辨率影像计算NDVI, 与研究所用的MODIS NDVI数据进行对比, 以验证其精确性。选择研究区域云量较少且植被茂盛, 符合验证需求的2015年7月3景Landsat 8-OLI遥感数据, 应用ArcGIS生成500个随机点, 提取2015年7月MODIS NDVI数据与Landsat 8-OLI数据估算的NDVI数据进行相关分析和曲线拟合。结果如图 2所示, 两组数据相关系数r=0.8006, 在P < 0.001水平上显著强相关。说明MODIS NDVI数据具有较高精度, 符合本文研究要求。

3 结果与分析 3.1 植被覆盖空间格局

通过NDVI表征植被覆盖度, 参考相关研究^[34-35]并结合宜宾市多年平均NDVI实际情况, 利用等间距法将植被状况划分为如表 1的5种类型, 最终得到宜宾市植被覆盖空间格局如图 3。总体上看, 宜宾市植被覆盖整体较好, 受地形与土地利用方式影响呈自南向北逐渐减少的趋势。与土地利用类型图对比可以发现, 在西部、南部达到高植被覆盖的区域也是主要的林地, 占全市面积52.06%；较高植被覆盖度则与耕地、草地相对应, 占全市面积45.83%；较低植被覆盖度主要分布在沿江地带和建成区等以水生植被和人工绿化植被为主的植被稀少区域, 仅占全市面积0.64%；低植被覆盖度主要沿长江干流分布, 水面开阔而反射率高, NDVI值极低。总之, 整体植被覆盖状况与土地利用方式高度相关。

3.2 NDVI时间变化趋势 3.2.1 区域尺度NDVI动态变化

图 4是基于线性回归的宜宾市区域尺度全年NDVI均值逐年变化情况。总体上看, NDVI处于较高的水平, 多年平均值为0.598, 且在2000—2015年间呈波动趋势上升。平均年际变化率达到0.007/a, R²为0.7736, 标准误为0.0186, 表明数据离散度较低, 拟合预测效果较好；全年NDVI标准差为0.0378, 变异系数为0.0632, 表明波动趋势明显。

图 5是区域尺度四季NDVI均值逐年变化情况。具体来看四个季节中, NDVI均值整体较高, 夏季最高为0.751, 秋季、春季次之为0.624、0.577, 冬季最低为0.443。NDVI年际波动均明显, 冬季尤为剧烈, 其标准误、标准差和变异系数最大, 分别为0.0454、0.0571和0.1291, 线性拟合较差, 春季、秋季年际波动相对次之, 夏季年际波动相对稳定。最后, 各季节中平均年际变化率最高的是春季和秋季, 均到达0.0083/a, 冬季次之为0.0073/a, 夏季最低为0.0039/a。其中, 夏季NDVI年际波动比较平稳但增长率较低, 可能是宜宾市植被以亚热带常绿阔叶林为主, 夏季草木繁盛, 增长相对其他季节不够明显；此外, NDVI饱和效应对夏季也存在很大影响。总的来说, 宜宾市NDVI在2000—2015年间逐渐波动增长, 但各季节增长速率与年际波动状况各有不同。

3.2.2 像元尺度NDVI动态变化

如图 6和图 7、8所示, 分别显示了宜宾市2000—2015年在全年及季节尺度上NDVI变化趋势与变化显著性。根据Sen氏趋势分析结果将变化趋势分为5级, 根据MK检验在α=0.05和α=0.01两个显著性水平上的结果将变化显著性分为6级。

全年NDVI变化趋势。总体上看, 全年NDVI在近16年中以极显著增长趋势为主, 极显著增长区域分布广泛, 达全市面积70.18%, 其中, 出现显著增长的区域又存在增长速率差异。快速增长区域占5.61%, 主要分布在宜宾县、翠屏区、高县、兴文县和筠连县山地区域；缓慢增长占据绝对地位为80.99%, 这与实施生态建设工程以来在原有耕地基础上人工造林、封山育林存在很大关系。同时, 随着宜宾市城市化与工业化的发展, 城市面积快速扩张, 人类活动干扰对植被造成了剧烈影响, 在沿江区域和建设用地区域等小部分区域也出现NDVI显著退化, 占全市面积0.22%。最后, 不存在显著变化趋势区域占全市面积14.92%, 主要分布在屏山县、江安县北部及沿江水域, 其NDVI基本保持稳定。屏山县境内变化不显著区域主要是大凉山余绪, 作为主要林区, 其植被茂密, 年际变化不大；此外, NDVI夏季饱和效应对此也存在很大影响。江安县北部是川南丘陵区, 作为宜宾市主要商品粮种植区, 天然植被较少, 耕地耕种收获所造成的植被年际变化也不大。造成以上区域NDVI发生微弱变化的影响因素包括自然与人为等多方原因综合, 难以用线性模型进行准确刻画^[36]。

各季节NDVI变化趋势。从整体变化显著性来看, 春季、秋季都以极显著增长趋势为主, 分布区域分别占全市面积48.82%、43.98%；夏季、冬季则以不显著增长趋势为主, 分布区域分别为46.29%、51.79%。从变化速率上看, 春季、秋季与冬季均以缓慢增长为主, 其分布区域分别达66.39%、64.65%与68.02%；夏季则是缓慢增长与保持稳定两种状态占据绝对地位, 分别达49.49%、48.20%。从空间分布上看, 春季、秋季大部分区域NDVI均出现极显著增长, 广泛分布于西部、南部地势较高区域。冬季NDVI出现极显著增长主要分布在西部、南部等生态建设一线区域, 常绿植被出现明显增长, 不显著增长则主要是东北部地势较低的丘陵区域。夏季则呈破碎化分布, 在东南部原本植被覆盖相对较低的喀斯特地貌区出现明显增长。此外, 夏季高植被覆盖的西北部屏山县、西南部山区变化往往不够显著, 这与NDVI饱和效应关系密切。各季节变化速率的空间分布也与变化显著性大体一致。与全年NDVI均值变化趋势相类似, 出现显著退化的区域也主要分布在沿江区域与县城等建成区。NDVI总体上都呈现海拔地势越高, 林地分布越多, 变化越明显且出现明显增长；海拔地势越低, 耕地分布越多, 变化越不显著且变化相对稳定趋势。

总的来说, 通过全年和季节NDVI变化均可发现, 自2000年实施生态建设工程以来, 执行退耕还林、封山育林的主要区域植被覆盖稳步增加, 生态环境不断优化, 生态建设成效明显。

3.3 基于Hurst指数的NDVI变化可持续性

通过逐像元计算Hurst指数如图 9, 可以看出, H值高于0.7的区域较少且分布比较破碎, 其面积仅占全市面积3%；H值高于0.6的区域分布相对集中于中部、西南部, 占13.76%。H值低于0.5的区域分布广泛, 占53.62%, 其中低于0.4的区域占20.37%, 主要是东南部、南部与西北部等山地区域出现了明显的反持续性特征。这表明宜宾市NDVI变化受自然与人为等难以量化的原因影响, 尤其是人类活动的影响, 其植被生长持续性不强, 需进一步加强生态保护工作, 特别是出现反持续性特征的山地区域。

根据H值是否大于0.5的正反持续性, 将Hurst指数与全年NDVI均值变化趋势通过栅格计算进行叠加后重分类, 得到如图 9的NDVI未来变化趋势分布图。据图可知, 继续保持快速增长的区域较少, 仅占2.56%, 继续保持缓慢增长的区域分布广泛占37.68%, 与Hurst指数高于0.6的分布区域较为一致, 主要分布在中部南溪区、翠屏区、宜宾县南部与西南部高县等区域。同时, 出现波动的区域占53.65%。主要是原本保持稳定的地区出现波动, 占46.42%, 分布在屏山县境内大凉山余绪、兴文县喀斯特地貌区以及江安县北部川南丘陵区。继续出现退化的区域分布极少, 仅占0.43%, 也主要分布在相关建成区。未来生态建设需要精准施策, 对波动、退化区域加强生态保护力度。

3.4 NDVI空间分布变化 3.4.1 植被覆盖空间转移矩阵

根据前文植被覆盖分级, 将各级植被通过ArcGIS空间叠加分析, 得到2000—2005年、2005—2010年、2010—2015年3个时段植被覆盖空间转移矩阵如表 2, 可以发现：

2000—2005年这一时间段较高植被覆盖与高植被覆盖面积占比大幅增加, 较低植被覆盖面积占比略有增加, 主要是由低一级植被类型转化而来。其中, 中植被覆盖与较高植被覆盖之间的相互转化占比较大, 植被覆盖保持稳定不变的区域占88.47%, 植被覆盖有所改善区域占8.79%, 植被覆盖发生退化区域占2.74%。总的来说转化趋势以正向演进为主, 生态建设成效初显, 生态环境发生良好变化。

2005—2010年这一时间段较高植被覆盖与高植被覆盖面积占比有所增加, 也主要是由低一级植被类型转化而来, 同时低植被覆盖面积占比存在微弱增加, 是由较低植被覆盖发生反复退化所致。其中, 植被覆盖保持稳定不变的区域占91.02%, 植被覆盖有所改善区域占6.88%, 发生退化区域占2.1%。其转化趋势也以正向演进为主, 发生在“中植被覆盖→较高植被覆盖→高植被覆盖”这一路径, 生态环境持续优化。

2010—2015年这一时间段高植被覆盖的面积占比大幅增加, 主要是由低一级植被类型转化而来。第一轮退耕还林工程实施十余年, 原本是坡耕地的区域通过植树育苗, 部分树种成熟, 林地立体结构初显, 植被覆盖率大幅增长。在这五年间, 植被覆盖保持稳定不变的区域占54.2%, 植被覆盖有所改善区域占45%, 发生退化区域占0.8%。其转化趋势也以正向演进为主, 显著发生于“较高植被覆盖→高植被覆盖”这一路径, 经过长期生态建设, 生态环境显著优化。

3.4.2 植被覆盖重心迁移过程

如图 10以3年为一个刻度计算不同植被覆盖类型重心迁移过程。总体上看, 不同植被覆盖类型的重心分布与植被覆盖空间分布格局具有较高一致性, 说明不同植被覆盖类型重心迁移过程能够很好的反映植被空间变化特征^[29]。低植被覆盖、较低植被覆盖重心主要稳定在宜宾市中北部沿江区域, 存在东西向曲折迁移过程, 这与两种植被类型空间分布严格一致。较高植被覆盖、高植被覆盖重心主要稳定在中部翠屏区、高县区域, 由南向北迁移, 与两种植被类型空间分布一致, 且与实施退耕还林后北部地区植被覆盖增加这一实际情况相吻合。

具体来看, 低植被覆盖重心在2000—2015年总体上向西南方向迁移, 迁移距离为16.76 km, 其原因主要是建设用地快速增加, 城市建成区尤其是宜宾市市区的快速扩张。较低植被覆盖重心主要在东西方向上迁移, 迁移距离为2.89 km；中植被覆盖重心主要在南北方向上呈S线型迁移, 迁移距离为3.44 km, 变化均不明显。较高植被覆盖重心主要向东北方向迁移, 迁移距离为12.29 km。高植被覆盖重心主要是向南北方向上呈折线型迁移, 迁移距离为16.25 km。在这16年间较高植被覆盖与高植被覆盖存在由南向北大规模扩张的趋势, 其重心迁移方向、迁移距离亦说明实施生态建设工程以来, 生态环境不断优化。

3.5 NDVI空间分布的地理探测

一般来说, 影响植被覆盖空间分布的因素主要包括气候因素、地表因素与人类活动因素^[37]。因此, 本文从以上三方面出发, 选取相应的代理变量, 形成如表 3的宜宾市NDVI影响因子。在ArcGIS中设置2 km规则渔网划分研究区, 再将NDVI及各影响因子分区提取到点, 最后纳入地理探测器模型计算。

3.5.1 探测因子对NDVI的影响力分析

根据因子探测, 计算各因子对NDVI的影响力大小, 结果如表 4。各因子对NDVI影响程度的排序为：海拔(X₅)>年均温(X₁)>土地利用类型(X₁₁)>人口密度(X₁₀)>日照时数(X₃)>坡度(X₆)>土壤类型(X₈)>GDP(X₉)>降水(X₂)>相对湿度(X₄)>坡向(X₇)。从各因子对NDVI的解释力来看, 海拔、年均温、土地利用类型和人口密度对NDVI的解释力都超过25%, 分别为39.9%、31.0%、28.1%和25.4%, 是主要影响因子；日照时数、坡度、土壤类型与GDP对NDVI的解释力都超过20%, 分别为24.7%、23.9%、23.2%与22.7%, 是次要影响因子。其余因子解释力相对微弱, 没有超过20%, 且相对湿度与坡向P值过大, 影响不显著。

表 5为生态探测结果。根据统计检验, 海拔与各气候因素对NDVI空间分布的影响具有显著性差异, 与其他因子无显著性差异；年均温与除海拔外的其他因子无显著差异；土地利用类型与年均温、日照时数、海拔、人口密度之间无显著差异；人口密度与解释力较强的因子之间均无显著差异。这进一步表明, 海拔、年均温对NDVI影响最大；土地利用类型、人口密度对NDVI影响极大, 但这种作用同时受到其他因子影响；日照时数、坡度、土壤类型与GDP对NDVI影响较大, 其他因子影响较小。

综合各因子可以发现, 以海拔、年均温为主的部分自然环境因子是影响宜宾市NDVI空间分布的主要基础因子。人类活动对NDVI空间分布影响很大, 但这种影响作用存在区域差异, 且与海拔、年均温的影响作用之间无显著差异。究其原因, 宜宾市土地利用存在明显的高低指向性, 建设用地、耕地等主要分布在地势低平的沿江及川南丘陵地带, 而地势较高的山地则以林地、草地为主, 是开展退耕还林等生态建设工程的重点区域, 人类活动相对较少。同时, 降水、相对湿度等影响植物生长的重要因素在宜宾市空间差异较小, 因而对NDVI空间分布影响微弱。因此, 可以确定海拔、年均温、土地利用类型与人口密度为影响宜宾市NDVI空间分布的主导因子。

3.5.2 主导因子交互作用分析

根据交互探测, 对NDVI影响因子进行交互作用分析, 选取主导因子交互作用结果如表 6。可以发现, 所有主导因子之间都存在显著的双协同增强作用, 即两个因子叠加大大增强了二者对NDVI空间分布的影响。这进一步佐证了以上因子对NDVI的主导作用, 并且主导作用的发挥通常是呈双协同增强状态。值得注意的是, 所有因子与海拔进行交互探测后, 都存在及其明显的增强效应, 特别是人类活动因素。这表明, 自然环境因素与人类活动因素共同作用下对NDVI空间分布影响更大。

3.5.3 影响因子的适宜范围或类型

根据风险探测, 计算分析各因子对NDVI分布变化的适宜范围或类型。考虑到人类活动因素受自然环境因素影响较大, 故对四个主导因子及次要自然环境因子的适宜范围或类型进行探讨, 其结果如表 7。不考虑人类活动影响, 在宜宾市仅有的中低山与丘陵地形, 当海拔高于1390 m, 坡度大于31.8°时, NDVI值较高, 这些区域山高坡陡, 不太适宜生产生活, 因而植被覆盖度较高。同时, 在年均温为14.8—15.6 ℃, 年日照时数为1196—1273 h时, 具有较好的光温条件, 植被覆盖度也比较高。此外, 暗黄棕壤作为广泛分布于亚热带中山上部区域, 对该地区植被生长也具有积极作用, 侧面反应了海拔对NDVI的主导作用。

总的来说, 通过地理探测器发现, 宜宾市气候要素空间差异较小, NDVI空间分布主要由海拔导致的地表差异, 包括海拔、年均温、土地利用类型、人口密度等构成的影响因素分异所决定。

4 结论与讨论 4.1 讨论

长江上游地区地形复杂破碎, 水土流失严重, 人类活动对植被的影响也极为剧烈。宜宾市作为三江交汇、高原盆地过渡区域, 其地理环境在长江上游流域极具代表性。本文从不同尺度对2000—2015年宜宾市NDVI时空变化进行了分析, 并在中小尺度上探究植被生长与自然环境、人类活动之间的响应机制。与以往关于植被变化和气候要素响应的研究结论有所不同, 研究发现影响植被空间分布的主要是海拔导致的地表差异。植被生长受到自然环境和人类活动的共同影响, 但人类活动也存在明显的环境导向。

通过区域与像元尺度的NDVI时空变化分析, 发现长江上游实施退耕还林、天保工程以来, 宜宾市植被覆盖稳步增加, 生态环境得到优化。这是人类活动对自然环境的正反馈, 说明生态建设成效明显。但是, 由于NDVI存在饱和效应, 尤其是夏季, 使得其生态建设成效评估仅局限于研究区域整体, 不能深入到更小尺度。在湿润环境下以常绿阔叶林为主的地区, 可以考虑使用EVI进行生长季变化研究, 较之NDVI具有更好的表现^[38]。本文主要对植被年际变化进行探讨, 受饱和效应影响相对较小, 故基于NDVI展开研究。此外, 不论是使用NDVI或是EVI进行植被变化研究, 都需注意植被指数对植被覆盖度真实值的偏离。本文使用同期高分辨率影像和MODIS数据验证证明, MODIS NDVI对研究区植被覆盖度具有较好的表征效果, 但MODIS遥感影像仍具有不确定性, 需要进一步检验^[39]。

地理探测器能够精准识别多因子之间的关系及其交互作用, 在驱动力机制分析方面应用广泛^[30]。目前应用地理探测器探究长时间序列NDVI空间分异机制的研究相对较少, 本文用其探究NDVI空间分异与影响因子之间的关系, 得出的主导因子符合长江上游地区植被生长自然规律, 与同类型研究结论也相对一致^{[20, 31]}, 地理探测器在这一领域的应用具有科学性。但是, 地理探测器在应用过程中对影响因子的空间分区没有明确的划分依据, 同时因子选取可能不够完备, 其因子解释力存在微弱的主观性。

4.2 结论

(1) 2000—2015年宜宾市的植被生长状态良好, 整体植被覆盖度较高, 受地形主要是海拔影响呈自南向北逐渐减少趋势。其中, 高植被覆盖在西部、南部山区分布广泛, 占全市面积52.06%。

(2) 从NDVI变化趋势来看, 2000—2015年宜宾市NDVI呈波动趋势上升, 平均年际变化率达到0.007/a, 全年NDVI标准误为0.0186, 标准差为0.0378, 变异系数为0.0632, 波动较为明显。同时, 极显著增长区域面积达70.18%, 广泛分布于各区县, 增长速率上以缓慢增长为主, 说明生态建设成效明显。

(3) 从NDVI变化可持续性来看, 宜宾市植被覆盖总体持续性不强, 许多区域出现波动, 主要是东南部、南部与西北部等山区出现了明显的反持续性特征。在未来可能出现波动的区域占53.65%, 主要分布在西北部、东南部及东北部, 需进一步加强退耕还林、封山育林等生态建设工作。

(4) 从NDVI空间变化来看, 2000—2015年宜宾市植被生长以正向演进为主, 主要是发生在“中植被覆盖→较高植被覆盖→高植被覆盖”这一路径, 生态环境不断得到优化；不同时段存在不同的演进效率和迁移路径, 演进过程与已有关于长江上游地区生态建设工程成效监测的研究结论相互佐证^{[34, 40-41]}。

(5) 地理探测器探究结果表明, 宜宾市气候因素空间差异较小, NDVI空间分布主要由一系列地表差异, 包括海拔、年均温、土地利用类型、人口密度所决定, 其因子解释力均超过25%, 自然环境因素与人类活动共同作用对NDVI影响更加显著。