植被作为陆地生态系统的主体, 在物质循环和能量流动中起着重要作用, 是联结大气、水分和土壤关系的自然纽带, 被认为是全球变化研究中的“指示器”^[1-3], 其动态变化受土壤、气候、人类活动等多重因素的共同影响, 气候因子中以温度和降水的影响最为直接^[4-7]。归一化植被指数(Normalized Difference Vegetation Index, NDVI)能较好的反映区域植被覆盖与植被生长状况, 在全球范围内得到广泛应用^[8-11]。

近百年来, 气候系统正经历以全球变暖为主要特征的显著变化, 尤其是在北半球高海拔、高纬度(如蒙古高原、北极)等自然过渡带与生境脆弱区, 温度升高显著高于全球温度升高的平均水平^[12-13], 由这种变化所引起的陆地生态系统效应已受到人类的高度重视。山地因在较小的空间具有浓缩的环境梯度和高度异质化的生境, 对气候变化的敏感性仅次于极地, 亦成为全球变化背景下研究植被动态对气候响应的热点区域^[14-16]。

秦巴山区(102°24′—112°40′E和30°43′—35°29′N)西接青藏高原东缘, 东临华北平原, 汉江以北为秦岭, 南为大巴山, 中部分布有汉中盆地和安康盆地。该区是中国地理格局的重要枢纽、大尺度东西向的生态廊道, 又是最高级别的生态安全保障区、扶贫开发重点区, 以复杂、多样、独特的自然环境孕育了丰富的天然动植物资源, 对其生态环境变化的研究具有重要意义。作为气候变化的敏感区和生境脆弱区, 秦巴山区植被覆盖变化及其驱动力研究是全球变化研究领域的重点内容之一, 国内学者在这方面已做了大量工作且取得一定成果^{[14, 17-19]}。雒新萍^[20]研究秦巴山区1982—2003年植被年均NDVI与气候因子的相关性, 认为秦巴山区植被覆盖与温度的相关性高于降水, 这与刘宪锋等^[21]的研究结果相反。崔晓临等^[22-23]分析了2000—2009年秦岭山地植被覆盖与温度的响应, 发现研究时段内秦岭植被覆盖整体呈增加趋势, 高海拔植被NDVI对温度变化更为敏感。邓晨辉等^[24]则指出2000—2015年间秦岭覆盖植被变化主要受人类活动影响, 其次是降水与温度。张清雨等^[25]对秦岭山地林地类型时空变化格局进行研究, 发现1978—2006年研究区林地面积增加但破碎化日益严重, 经济发展、适耕地面积等是林地类型转换的主要驱动力。然而, 目前关于秦巴山区植被覆盖变化及其驱动力研究主要集中于陕西省境内^{[21-22, 26-28]}, 对西秦岭、伏牛山和大巴山研究较少, 且前人对秦巴山区植被覆盖变化的主要气候驱动因子争议较大, 大多仅用一种植被数据对植被覆盖变化进行探究, 全面、准确评价整个秦巴山区植被覆盖时空变化特征及其对气候因子的响应机制、整合碎片化、局域化资料显得尤为重要与紧迫。鉴于此, 本文使用GIMMS_3g(1982—2015年)、SPOT VEG(1998—2015年)、MODIS(2000—2017年)3种不同时空分辨率的NDVI数据集, 对秦巴山区1982—2017年植被覆盖的时空变化特征及其与主要气候因子的相关性进行研究, 并以秦岭山地(陕西境内)为重点, 分析了植被覆盖的海拔梯度差异及其与气候因子的垂直响应模式, 以期在全球气候变化的背景下为山地植被保护与恢复提供科学依据。

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

秦巴山区(图 1)地跨甘、川、陕、渝、鄂和豫等6个省市, 东西长约1000 km, 南北宽约300 km, 总面积约为30万km²。以秦岭-大巴山为主体的秦巴山区是长江-黄河流域的分水岭, 也是我国南北过渡带的主体区。过渡带的地理位置使得该区植被区系复杂, 秦岭北坡主要为暖温带落叶阔叶林, 南坡主要为混有常绿树种的落叶阔叶林, 大巴山北坡与秦岭南坡植被相似, 南坡则为典型的常绿阔叶林, 随海拔升高, 气候、植被等呈明显的垂直地带性分布^[29-30]。该地区年平均气温介于12—16℃, 年平均降水量介于709—1500 mm。因其地理环境的高度复杂性、生物多样性和气候敏感性, 该区已成为我国地学、生物学研究的天然实验室。

1.2 数据来源与预处理

为了延长对区域植被覆盖变化的研究时间, 避免遥感数据本身的机械误差, 本研究采用了三套不同时空分辨率的植被遥感数据集。其中GIMMS_3g V1.0(1982—2015年)数据来源于美国全球检测与模型组(Global Inventor Modeling and Mapping Studies, GIMMS) 15 d合成的空间分辨率为8 km(https://ecocast.arc.nsds.gov)的NDVI数据集；SPOT VEG(1998—2015年)数据来源于中国科学院资源环境科学数据中心提供的逐月NDVI数据集(http://www.resdc.cn), 空间分辨率为1 km；MODIS(2000—2017年)数据来源于美国国家航空航天局(NASA)16 d合成的空间分辨率为250 m MOD13Q1产品。利用MRT(MODIS Reprojection Tool)工具对MODIS数据集进行拼接、剪裁等预处理, 提取NDVI影像。为减少云、气溶胶、视角等影响, GIMMS_3g和MODIS数据均采用最大值合成法MVC(Maximum Value Composites)^[31]得到月NDVI。

气象数据来源于中国科学院资源环境科学数据中心(http://www.resdc.cn)提供的中国年降水量、年平均气温空间插值数据集, 空间分辨率为1 km。

DEM数据(空间分辨率为30 m)来源于中国科学院计算机网络信息中心国际科学数据镜像网站(http://datamirror.csdb.cn)。

1.3 研究方法

采用均值法得到年均NDVI和多年累计平均NDVI, 为了避免遥感数据本身带来的误差, 将年均NDVI进行z-score标准化处理, 利用线性趋势分析法探讨研究区植被覆盖年际变化特征, 并基于像元尺度的线性趋势分析法对多年年均NDVI进行处理, 以回归方程斜率表征不同区域植被覆盖空间变化特征, 斜率为正时, 植被呈改善趋势(NDVI随时间推移而增加), 斜率为负时, 植被呈退化趋势(NDVI随时间推移而降低), 且斜率绝对值越大, 植被覆盖变化趋势越明显；最后采用相关分析法来探明区域植被覆盖变化对气候因子的响应状况。另外, 水热组合、土壤养分等条件随海拔变化导致植被覆盖有很大的差异^[32-35], 本文以秦岭山地(陕西境内)为重点区, 参考《秦岭山地自然地理》^[36], 根据秦岭植被分布的垂直带性将海拔划分为6个梯度(< 500 m, 500—1000 m, 1000—1600 m, 1600—2700 m, 2700—3200 m, >3200 m)以观察植被覆盖的垂直变化特征。

由于本文采用了3种不同时空分辨率的植被遥感数据集, 根据研究需要, 在分析秦巴山区植被覆盖长时间序列(1982—2017年)年际变化特征时采用3种数据, 这样既能对各种遥感数据质量互相验证又延长了对区域植被变化的研究时段；分析植被覆盖空间分布、空间变化(包括不同海拔带)特征时采用高分辨率的MODIS数据；分析植被覆盖变化对气候因素的响应时采用和气象数据具有相同空间分辨率的SPOT VEG数据。

2 结果与分析 2.1 植被NDVI(GIMMS_{3g 1982—2017}、SPOT VEG_1998—2015和MODIS_2000—2017)年际变化特征

将秦巴山区1982—2017年(GIMMS_3g、SPOT VEG和MODIS数据)年均NDVI进行标准化处理, 用其年际变化趋势代表研究区植被覆盖年际变化趋势(图 2)。

结果表明, 基于不同分辨率的3种遥感数据在重叠年份的年际变化趋势一致, 充分证明了该区域遥感数据的可靠性, 适用于分析秦巴山区植被覆盖的时空变化格局。2000—2017年间, 三种遥感数据NDVI均呈现稳定增加趋势, 且年际波动表现为较好的一致性, GIMMS_{3g 2000—2015}、SPOT VEG_2000—2015和MODIS_2000—2017增速分别为2.2%/10a(P < 0.001)、4.1%/10a(P=0.001)和4.5%/10a(P < 0.001), 由于传感器自身差异, GIMMS_3g数据增速略小于SPOT VEG和MODIS数据增速, 但整体上可以准确显示研究区植被覆盖变化趋势。由GIMMS_3g数据可知, 秦巴山区1982—2015年间植被覆盖整体为增加趋势, 增速为1.4%/10a, 其中1982—2000年增速为1.6%/10a(P=0.098), 植被覆盖变化主要受气候因子影响, NDVI水平较低且年际波动较大, 仅1990、1991、1995、1997和1998年NDVI高于平均值, 其余年份均低于平均水平, 在1998—2000年间NDVI迅速下降(与SPOT VEG数据表现一致), 该时间段为研究区实施退耕还林工程初期, 耕地面积急剧减少^[37], 导致植被覆盖度降低；2000年之后, 生态保护工程成效显著, 在人为作用的影响下, 植被覆盖呈稳定增加趋势(1.7%/10a, P < 0.001), 除2005、2011、2012年有较大下降, 其余年份NDVI值均处于较高水平且稳定上升。计算GIMMS_{3g 1982—2000} NDVI年际变异系数(CV)为0.04, 2000—2015年为0.02, SPOT VEG_2000—2015和MODIS_2000—2017 CV值分别为0.04和0.05, 对比不同时期GIMMS_3g数据的CV值可知, 2000年之前, 秦巴山区植被覆盖变化受气候因素影响年际波动较大, 而2000年之后, 在人为因素(如植树造林工程)影响下稳定增加。

2.2 植被NDVI(MODIS_2000—2017)空间分布特征

由于植被类型、土地利用、气候因素等不同, 秦巴山区2000—2017年植被覆盖呈现明显的区域差异(图 3a), NDVI介于0.27—0.74之间, 平均值为0.59。对比地形图(图 1)可以发现, 研究区植被覆盖最好的区域主要为太白山、米仓山、雪宝顶和神农架等山地, NDVI达到0.60—0.74, 这些区域植被覆盖以森林为主, 植被组成主要有常绿阔叶林、落叶阔叶林以及暖温带针阔叶林等；研究区中部的汉中盆地土地利用类型主要为农耕地及人工园林, 与林地相比, NDVI有所下降, 集中分布于0.4—0.5之间；研究区东北部海拔较低, 人口密集, 植被覆盖受人为干扰严重相对较低, NDVI多分布于0.30—0.40之间, 另外, 西秦岭(研究区西北、西南边缘的岷县、武都县以及西和县)以及研究区海拔最高的西南部等地区, 因受水热条件制约, 植被类型以亚高山落叶阔叶灌丛、禾草、杂草草甸为主, 植被覆盖度相对较差, NDVI小于0.35。

由图 3可以看出, 秦岭山地随着海拔的升高, NDVI值呈先增大后减小的变化特征, 海拔 < 500 m时, NDVI相对较小且分布分散, 平均值为0.52, 主要介于0.50—0.60之间, 占该海拔带总面积的55.50%, 其次为0.40—0.50和0.60—0.70, 所占比例分别为23.62%、13.52%, 小于0.30和0.30—0.40之间分布较少；500—1000 m, NDVI有所上升且分布较为集中, 平均值为0.59, 主要介于0.50—0.6与0.6—0.7之间, 所占比例分别为41.69%、49.21%, 0.30—0.40分布较少；1000—1600 m, NDVI继续上升, 平均值为0.63, 介于0.60—0.70之间的比例为71.13%；1600—2700 m, 平均值为0.67, 植被覆盖状况最好且分布最为集中, NDVI主要介于0.60—0.70之间, 所占比例为71.47%, 大于0.7的分布比例次之；2700—3200 m, NDVI略有下降, 平均值为0.59, 0.6—0.7之间所占比例为39.00%；>3200 m时, NDVI最小, 平均值为0.34, 小于0.4的区域占78.11%。

植被分布受多种因素的影响, 其中人为因素与气候因素最为重要^[38-41]。海拔小于500 m时, 植被类型主要为农耕栽培植被, 500—1600 m海拔带植被覆盖有所升高, 但由于受人类活动影响, 不少山坡已被开垦, 原始森林植被遭到破坏, 植被类型主要为人工林地, 植被组成单一^[42], 且山麓气候干燥, 故植被覆盖较差；海拔升高至1600—2700 m时, 人为干扰较小, 森林植被保存完好, 水热条件适宜, 植被覆盖状况最好, 植被组成主要为针阔叶混交林^[36]；2700—3200 m为亚高山植被带, 植被成分单纯, NDVI值迅速下降, 海拔高于3200 m时, 受温度限制, 植被主要由灌丛和草甸组成, NDVI值最小。

2.3 植被NDVI(MODIS_2000—2017)空间变化特征

利用线性趋势分析法对2000—2017年各年均NDVI值与其相对应的时间序列进行逐项元回归分析, 将其回归斜率划分为五类, 用以表征该时间段内植被覆盖变化趋势, 划分标准及结果见图 4、表 1。结果表明, 18 a间秦巴山区植被覆盖改善的区域面积约20.07万km², 占总面积的96.90%, 退化的区域面积约0.93万km², 所占比例为3.10%, 即研究区植被覆盖整体呈增加趋势。其中轻度改善的面积最大, 为13.79万km²(45.97%), 中度改善的面积次之, 为13.42万km²(44.73%), 明显改善的区域面积为1.86万km²(6.19%), 之后为轻度退化和中度退化, 面积分别为0.71万km²(2.35%)和0.14万km²(0.47%), 明显退化的面积最小, 为0.09万km²(2.9%)。从空间分布来看, 植被明显改善的区域主要集中在西秦岭(甘南的礼县、武都县、西和县), 另外, 安康中部和十堰北部等低海拔区域植被覆盖也以明显改善为主；中度改善的区域较为零散, 遍布整个研究区；轻度改善的区域主要集中在研究区北部高海拔区域, 太白山、米仓山和神农架等山地也多表现为轻度改善；植被退化的面积较少, 主要集中在研究区西南部的少数地区, 如成都北部和德阳北部的交界处, 汉中盆地及南阳西部边缘区域也有零星分布。由上述分析可得, 近18 a秦巴山区植被覆盖整体呈改善趋势, 以轻度改善和中度改善为主体, 变化斜率较大的区域主要集中在海拔相对较低、人类活动较为频繁的地区, 高海拔区域主要受气候因素影响, 植被覆盖变化主要表现为基本不变(轻度退化、轻度改善)。

为进一步探究海拔对植被覆盖变化的影响, 统计秦岭地区不同海拔带多年来植被覆盖的变化情况, 结果显示, 不同海拔带植被覆盖均呈增加趋势, 但增速有所不同(图 4)。海拔 < 500 m、500—1000 m、1000—1600 m时, 植被类型以农作物与人工园林为主, 植被覆盖变化主要受人类活动(耕作制度、灌溉及施肥方式、退耕还林政策等)影响, NDVI增速较大, 分别为4.8%/10a、5.7%/10a、5.1%/10a, 且均通过了置信水平为99%的显著性检验；海拔在1600—2700 m与2700—3200 m时, 植被覆盖变化受人为干扰较小, 主要受温度和降水等气候因子影响, 在全球变暖的大背景下, 该区植被覆盖有所增加, 但与低海拔相比增速减小, 分别为3.2%/10a(P < 0.001)、2.6%/10a(P < 0.05)；随着海拔的继续升高(>3200 m), NDVI变化不显著, 增速为1.2%/10a(P>0.05)。进一步说明人为因素对植被覆盖影响的速度大于气候因素对植被覆盖影响的速度, 人类活动可在短时间内改变区域植被覆盖格局。

2.4 植被NDVI(SPOT VEG_1982—2015)变化的气候因子驱动力分析

温度和降水是影响植被覆盖最直接、重要的气候因子^[7], 在全球气候变化的大背景下, 1982—2015年秦巴山区以增温效应为主, 温度增加速率为0.49℃/10a(R=0.826, P < 0.001, n=34), 增湿效应不显著。在这种情况下, 本文以年为单位, 逐项元分别计算2000—2015年植被NDVI与年均温和年降水的相关系数, 得到其空间分布图, 用以分析研究区植被覆盖对气候因子的响应模式, 并对相关系数进行统计分类(表 2)。结果显示, 研究区植被覆盖变化与温度主要为正相关, 与降水为正负相关并存, 与温度的相关性较降水高。

2000—2015年秦巴山区植被覆盖变化与温度的平均相关系数为0.22, 区域正相关系数和负相关系数所占比例分别为88.81%和11.19%(图 5, 表 2), 即在研究时段内, 研究区植被覆盖变化与温度主要呈正相关。在正相关区域中, 低度正相关的区域占总面积的57.12%, 分布于整个研究区；中度正相关和高度正相关比例分别为10.61%、0.58%, 主要分布在太白山、米仓山和神农架等森林面积较大的高海拔区域, 这些地区降水充沛但温度较低, 热量为该区域植被生长的主要限制性因子, 温度升高有利于植被覆盖的增加。负相关的区域面积较小, 主要分布在四川北部的松藩县、平武县和茂县等地区, 另外, 豫西地区人口密集、人为干扰较大的洛宁县、栾川县、内乡县一带也与温度呈负相关, 但相关性均不显著。

相比之下, 植被覆盖变化与降水的正负相关性较为平衡, 平均相关系数为0.03, 正相关与负相关所占比例分别为41.53%和58.47%。其中呈低度正相关的区域比例为49.14%, 主要分布在研究区北部和南部边缘地区, 中度和高度正相关的区域比例分别为7.86%和1.66%, 主要位于研究区西北部的临潭县、岷县和温县等降雨较少的地区, 这些地区温度较高, 土壤水分蒸发严重, 降水增加将促进植被生长；呈低度负相关和中度负相关的区域比例为34.98%和6.55%, 主要集中在研究区东部城镇化水平较高的大面积低海拔区域, 另外, 西南部的高海拔区域也呈负相关, 这些区域降水量相对较多, 植物生长所需的水分足够, 降水增加将导致光照强度下降、光合作用减弱, 从而抑制植被的生长^[43]。

不同海拔带植被覆盖变化对气候因子的响应具有明显的差异性(图 5)。结果表明, 随着海拔的升高, 植被覆盖与温度的相关性出现由负到正的转变。海拔小于1600 m时, 植被NDVI与温度主要呈负相关, 但均未通过显著性检验, 1000—1600 m海拔带与温度相关系数最大；海拔在1600—3200 m时, 相关性转为正相关, 2700—3200 m时, 相关系数略有上升, >3200 m时, 相关系数最大为0.35, 且通过了置信水平为90%的显著性检验, 即高海拔区域植被对温度变化更为敏感；与温度相反, 随着海拔的升高, NDVI与降水的相关性出现由正到负的转变, 海拔在低于2700 m时两者为负相关, 高于2700 m时两者为正相关, 低海拔区域相关系数更大, 但不同海拔带相关性均未通过显著性检验。造成这种现象的原因主要为：在低海拔区域, 植被覆盖类型主要为农作物和人工园林, 植被覆盖与人类活动的相关性高于与气候因子的相关性, 且栽培植被对水分要求较高, 因此与降水呈正相关, 随着海拔的升高, 降水增加而温度降低, 植被生长主要受热量条件限制, 因此与温度呈正相关。

3 结论

(1) 时间变化上, 1982—2017年3种NDVI遥感数据在秦巴山区均呈显著上升趋势。1982—2015年NDVI-GIMMS_3g增速为1.4%/10a, 其中1982—2000年NDVI-GIMMS_3g整体较低且波动幅度较大, 增速为1.6%/10a, 年际间变异系数(CV)为0.04, 2000年后NDVI较高且呈稳定增加状态, 增速为1.7%/10a, CV为0.02；2000—2015年NDVI-SPOT VEG增速为4.1%/10a, CV为0.04；2000—2017年NDVI-MODIS增速为4.5%/10a, CV为0.05。

(2) 空间分布上, 2000—2017年(MODIS数据)秦巴山区高植被覆盖区主要为太白山、米仓山、神农架等山地, 低植被覆盖区主要为研究区西部的高海拔区和东部受人为干扰严重的低海拔区；随海拔升高, 秦岭山地NDVI呈先增加后减少的变化趋势。

(3) 空间变化上, 2000—2017年(MODIS数据)秦巴山区植被覆盖增加的区域大于减少的区域, 增加和减少的区域分别占总面积的96.90%和3.10%, 低海拔区域变化显著, 高海拔区域基本不变；秦岭山地不同海拔梯度带植被覆盖均呈增加趋势, 低海拔增速较高海拔大。

(4) 1982—2015年(SPOT VEG数据)秦巴山区增温效应显著, 增速为0.49℃/10a, 增湿效应不显著；2000—2015年研究区植被覆盖与温度以正相关为主, 主要分布在山地等降水较多的地区；与降水为正负相关性并存, 正相关主要分布在研究区西北部降水较少的半干旱区, 负相关主要分布在东部城镇化水平较高的低海拔区域；秦岭山地随着海拔的升高, 植被覆盖与温度的相关系数呈由负到正的转变, 高海拔区域对温度更敏感, 与降水的相关系数呈由正到负的转变, 低海拔区对降水更敏感。

研究结果表明, 1982—2017年秦巴山区植被覆盖呈显著增加趋势, 1982—2000年植被NDVI较低且年际波动大, 2000年后植被NDVI较高且稳定增加。植被动态变化受气候、人为和历史等多种因素共同影响, 在研究时段内, 研究区温度和降水均呈增加趋势, 增温效应尤为明显, 但对2000—2015年植被覆盖与气温和降水进行相关分析, 发现植被覆盖与温度以正相关为主, 与降水为正负相关并存, 且高海拔区域植被对温度变化更为敏感, 低海拔区域对降水更为敏感, 这与以往的研究结果相似^{[13, 20, 22, 27, 36]}, 但无论是气温还是降水, 呈明显相关的区域面积并不是很大^[44], 造成这种现象的可能性主要有以下几点：首先, 气候虽为植被覆盖变化的决定性因子, 但两者的相关性受研究时段和空间尺度影响^[5], 研究的时空尺度不同, 相关系数将存在较大的差异性；其次, 植被覆盖变化对气候的响应模式存在显著的空间差异性, 并可能存在一定的滞后性^{[14, 45-47]}, 本文仅研究了年均温、年降水量对植被NDVI的影响, 未对其滞后性进行分析；另外, 有研究表明, 近年来植被覆盖变化的主要驱动力已由气候因素转为人为因素^[48], 人类活动能在短时间内改变区域环境, 退耕还林、植树造林等生态保护工程可快速增加区域植被覆盖^{[21, 49-50]}, 自2000年后, 研究区全面实施退耕还林政策, 林地增加面积为2.54×10⁶ hm², 约占全国退耕还林总面积的10%^[37], 植被覆盖变化受人为因素干扰较大。因此, 秦巴山区植被覆盖变化对不同时段气候因子的响应以及定量分析人为因素(社会经济、城市扩张等)对植被覆盖的影响是以后工作的重点内容, 以期更为全面的揭示秦巴山区植被覆盖对气候因素与人为因素的响应机制, 为山地植被保护和恢复提供科学依据。