轿子山国家级自然保护区气候变化特征及对NDVI的影响

任钇潼，程清平,2,3，金韩宇，张宏英，王 平

(1.西南林业大学地理与生态旅游学院，昆明 650224；2.西南生态文明研究中心，昆明 650224；3.中国科学院西北生态环境与资源研究院玉龙雪山冰冻圈与可持续发展国家(云南省)科学野外观测研究站，兰州 730000；4.云南师范大学地理学部，昆明 650500)

保护区在减少栖息地损失和维持物种的可持续方面发挥着至关重要的作用，气候变化和人为压力已对保护区生物多样性构成巨大威胁。已有研究指出，在气候变化下，江西武夷山、河北大海陀、吉林长白山、宁夏贺兰山、西藏羌塘、西藏雅鲁藏布大峡谷和甘肃祁连山等保护区NDVI对气温、降水响应敏感，除气候影响外，张玉琴等、李娟等、付健新等指出地形地貌因子(坡度、坡向)亦会影响NDVI变化趋势。已有研究主要侧重于平均气温、降水对保护区NDVI的影响。近年来，极端气候事件呈增加、增强、持续时间延长的态势，将可能显著影响区域植被动态变化，保护区气候变化风险与日俱增。因此，势必需要进一步揭示保护区NDVI对极端气候变化的响应。当前，中国气候脆弱性较高的保护区主要分布在西南地区和青藏高原的部分地区，尤其是云南气候变化风险更高。轿子山国家级自然保护区(简称轿子山)位于滇中高原，作为我国少有的高山保护区之一，保存有代表性的滇中高原植被与生境垂直带谱，对气候变化尤为敏感，识别极端气候事件对保护区植被的影响特征，可为保护区气候风险评估和监测预警提供参考依据。

基于此，本研究利用保护区所在的东川区与禄劝彝族苗族自治县(禄劝县)20世纪50—80年代气候要素资料和1961—2019年逐日气候要素数据集、2017—2019年保护区周边自动气象站逐时气温、降水要素数据集、1999—2016年1 km月尺度归一化植被指数集(NDVI)，围绕2个问题展开研究：(1)气候(极端)要素和温度带垂直组合随时间如何变化？(2)保护区NDVI如何响应极端气候变化？

1 研究区概况

轿子山国家级自然保护区(102.80°—102.98°E，26.00°—26.20°N)位于东川区与禄劝县交汇处，分为轿子山片区和普渡河片区，总面积16 456 hm(图1)。山势南北延伸，主峰为雪岭(4 344.1 m)，该峰是青藏高原以东海拔最高的山地，也是北半球该纬度带上最高山地之一。轿子山位于北亚热带半湿润气候区，水分和热量条件垂直分异显著，夏季受热带海洋气团控制多云雨，冬季受热带大陆气团控制少降水。区域属森林生态系统类型自然保护区，主要保护对象为以攀枝花苏铁、须弥红豆杉、林麝等为代表的珍稀濒危野生动植物资源及其栖息环境，分布有我国面积最大的高山柏林和海拔最低的高山松林。

图1 站点分布与研究区概况

2 数据来源与方法

2.1 数据来源

东川区和禄劝县气象数据资料均由东川区和禄劝县气象局提供，所搜集气象数据资料见文献[18-21]。1 km月尺度归一化植被指数(NDVI，1999—2016年)来自于中国科学院资源环境科学与数据中心(https://www.resdc.cn)，30 m数字高程模型(DEM)来自地理空间数据云(http://www.gscloud.cn/search)。

2.2 研究方法

采用一元线性回归、Pettitt非参数法分析气候要素、极端气候指数、NDVI趋势特征。通过Arcgis 10.5反距离权重插值20世纪50—80年代、2017—2019年平均气温、年降水量的栅格并利用等值线工具分别提取等温线与等降水量线。对DEM(30 m)进行栅格重分类，划分间隔为500 m，将保护区划分为4个海拔带，统计不同海拔带内栅格值，用以研究不同海拔带NDVI变化特征。此外，为探究植被对极端气候事件的响应，将NDVI与极端气候指数进行Spearman相关分析。

2.2.1 极端气候指数 结合保护区气候特征，从ETCCDI定义的极端指数中选取14个表征极端气温指数和9个表征极端降水指数，定义见http://etccdi.pacificclimate.org/list_27_indices.shtml。采用R ClimDex软件计算东川、禄劝站极端气温/降水指数。

2.2.2 温度带划分 参考郑景云等和王宇温度带划分标准，探讨全球变暖下轿子山垂直方向温度带的动态变化，各指标分类见表1。以日平均气温稳定≥10 ℃积温日数作为主要指标，能够准确辨别山地地区热量地域分异。落雪站(可代表轿子山高海拔地区气候状况)历史1月平均气温为1.3 ℃，7月平均气温为11.8 ℃，由此可知，1月平均气温和7月平均气温并不是轿子山地区主要限制因子。故本文以日平气均温≥10 ℃持续日为主要指标，以≥10 ℃积温为参考指标，1月平均气温与7月平均气温为辅助指标划分轿子山温度带。

表1 云南温度带划分标准[26-27]

对于垂直温度带的划分，2017—2019年研究区有5个自动气象站和2个国家站的气象数据可用，其中东川和禄劝站在2个时段均有，且最低点均是东川站，尽管前后站点不同，但在空间上均比较接近，因而可用于揭示2017—2019年相对于20世纪50—80年代垂直温度带的变化。依据表2推算轿子山每隔100 m高度下的≥10 ℃积温日数、≥10 ℃积温、1月平均气温以及7月平均气温，并以此进行温度带识别划分。

表2 轿子山国家级自然保护区≥10 ℃积温日数、≥10 ℃积温、1月平均气温、7月平均气温垂直递减率

3 结果与分析

3.1 区域气候

3.1.1 气候要素时间变化趋势 从图2可以看出，轿子山年平均气温显著升高，年日照时间、年平均风速、年雾日呈显著降低态势，年降水量、年平均相对湿度趋于缓慢降低。从显著突变点来看，年雾日、年日照时间突变年相对提早，年平均气温、年平均风速较为滞后，年平均相对湿度突变年最为滞后。

图2 1961-2019年轿子山国家级自然保护区气候要素年际变化

3.1.2 平均气温空间变化特征 由图3a～图3f可知，20世纪50—80年代轿子山年平均气温、1月平均气温、7月平均气温总体分布呈南高北低、西高东低趋势。轿子山西部、北部等温线分布较为密集，表明温差大，可能与轿子山地形关系密切。相较于20世纪50—80年代，2017—2019年年平均气温、1月平均气温、7月平均气温表现为北高南低、东高西低趋势。等温线分布特征均表现为东密西疏，可能是由于2017—2019所使用的部分站点与20世纪50—80年代不同而对平均气温空间分布格局产生影响。总体来看，轿子山平均气温变化具有海拔依赖性的特点，且各局部地区平均气温趋于上升。从平均气温变幅来看，轿子山西部(地势较高)气温增幅小于轿子山东部(地势较低)。

3.1.3 年降水量空间变化特征 由图3g～图3h可知，轿子山年降水量空间分布特征为西多东少。20世纪50—80年代，轿子山西部等降水量线分布较为密集，向东等降水量线分布逐渐稀疏，即年降水量的空间差异向东逐渐减小。由于20世纪50—80年代部分站点终止使用，2017—2019年所使用的站点数量以及站点地理参数上的不同可能对等降水量线的分布产生影响，但不会造成最终结论偏差，因此仍有研究意义。基于此，相较于20世纪50—80年代，2017—2019年等降水量线表现为东密西疏，即年降水量的空间差异东大西小，轿子山除西南角、东北角年降水量趋于减少外，大部分地区呈增加态势，西南区域变幅最大。

图3 20世纪50-80年代与2017-2019年轿子山国家级自然保护区平均气温、年降水量分布

3.2 温度带划分

从图4可以看出，20世纪50—80年代，轿子山垂直温度带依次为南亚热带、中亚热带、北亚热带、暖温带、中温带、山地寒温带和山地亚寒带。相较于20世纪50—90年代，2017—2019年轿子山基带温度带为北热带，南亚热带下界随之抬升200 m，中亚热带上界降低约100 m，北亚热带、暖温带、中温带、山地寒温带、山地亚寒带下界下降约100～300 m。

图4 20世纪50-80年代与2017-2019年轿子山国家级自然保护区温度带

3.3 极端气温与降水事件分析

3.3.1 极端气温 由图5可知，东川站暖指数(除TXx)趋于增长，冷指数(除TXn、TNn)趋于下降，且仅TN10p通过0.05显著检验。突变点通过0.05显著性检验的极端气温指数中，暖指数多于冷指数，TXx突变相对提早(1974年)，TNx、SU25、TX90p、TN10p突变相对同步(2008—2010年)。禄劝站暖指数均趋于上升，以TXx、TNx、SU25、TR20、TX90p、TN90p、WSDI最为显著，冷指数(除TNn)呈下降态势，以TNn、FD0、TN10p、CSDI最为显著。突变点通过0.05显著性检验的极端气温指数中，暖指数多于冷指数，冷指数(1989—2000年)突变相对同步性强于暖指数(1991—2005年)。

图5 1961—2019年轿子山国家级自然保护区东川、禄劝站极端暖、冷指数年际变化

3.3.2 极端降水 由图6可知，东川站强度指数、极端降水量指数(R95p除外)趋于不显著增长，极端降水日指数(R10 mm除外)趋于不显著降低，极端干旱日指数趋于轻微下降，未有极端降水指数突变点通过显著性检验。禄劝站强度指数趋于不显著增长，极端降水量指数(除R99p)、极端降水日指数趋于不显著降低，极端干旱日指数趋于轻微上升。仅RX1与R99p存在显著突变且均为1984年。综上，东川、禄劝站降雨强度均有所增强，极端降水日趋于减少，东川站极端降水量趋于缓慢增加，禄劝站趋于缓慢减少。

图6 1961—2019年轿子山国家级自然保护区东川、禄劝站极端降水指数年际变化

3.4 NDVI时空特征以及对极端气候指数的响应

3.4.1 NDVI时空特征 由图7可知，不同海拔范围下NDVI均趋于显著增长，NDVI增速随海拔升高减小。1 736～3 236m海拔带NDVI值为0.6～0.8，3 236～3 736 m海拔带内NDVI值范围为0.4～1.0。2003年以前，海拔1 736～2 236，2 236～2 736 m区域内NDVI值为0.6～0.8，并且稳定不变，2 736～3 236，3 236～3 736 m海拔带内分别约有11.86%，10.77%的栅格NDVI值降低。2003年—2010年，4个海拔带(从低到高)分别有约33.33%，16.67%，11.86%，4.62%的栅格NDVI值增长，2010—2012年又存在不同程度的NDVI值下降。2012—2016年，各海拔带分别有约33.33%，50.00%，53.34%，28.41%的栅格NDVI值趋于增长。

图7 1999-2016年轿子山国家级自然保护区不同海拔带NDVI趋势变化以及值结构

3.4.2 NDVI对极端气候事件的响应 图8为东川站和禄劝站极端气候指数与轿子山年平均NDVI(NDVI)、年最大NDVI(NDVI)、年最小NDVI(NDVI)的Spearman相关。NDVI、NDVI与暖/冷指数多呈正/负相关，且NDVI与暖指数中的TXx(东川)、SU25(禄劝)、TR20呈显著正相关，与冷指数中的TX10p(东川)呈显著负相关。NDVI与暖/冷指数多呈负/正相关，与其中的WSDI(东川)、TNx呈显著负相关。

注：*表通过0.05显著性检验。

总体来看，3个NDVI指数与极端降水指数相关性偏弱，仅NDVI与RX1(禄劝)，NDVI与R25 mm(东川)呈显著正相关，即降雨强度、极端降雨日数分别对NDVI、NDVI升高存在促进作用。

4 讨 论

4.1 气候变化特征

从表3可以看出，多数山地保护区仅研究气温、降水2个气候要素：气温趋于升高或降低，降水除轿子山外均趋于增加。本研究表明，轿子山年日照时间、年平均风速、年雾日趋于显著降低。此外，平均气温与降水量等值线分布未考虑地形影响，其插值结果可能存在一定的偏差(图3)。轿子山暖指数、降水强度指数趋于增长，冷指数、极端降水日指数趋于下降，与Shao等的研究结论一致。此外，张善红等和Guan等分别发现，秦岭地区和青藏高原温度带(气候带)正向高海拔攀升，本研究表明，轿子山同样存在该现象，即基带温度带由南亚热带转为北热带，原南亚热带整体高度抬升。

表3 国家级山地自然保护区气候要素趋势比较

4.2 NDVI趋势的海拔差异

轿子山海拔1 736～2 236，2 236～2 736，2 736～3 236，3 236～3 736 m区域内NDVI斜率分别为7.28×10，5.85×10，5.28×10，3.75×10a，表明随海拔升高，NDVI增速减缓。Tao等研究表明，西南地区不同海拔带NDVI均趋于增长，且增速随海拔升高而上升，以海拔1 200～3 500 m范围内NDVI增速的垂直变化率最高，每1 000 m升高0.06%/a；Jiang等研究发现，西南地区海拔<3 400 m，NDVI显著增长，3 400 m以上NDVI先增后降不一致；而Li等研究表明，雅鲁藏布江大峡谷保护区NDVI增速随海拔升高显著降低，海拔≤3 000 m、NDVI增速为0.05/10～0.08/10 a、海拔>3 000 m，NDVI增速为0.01/10～0.03/10 a。上述差异可能与研究区域大小和海拔带划分标准不同有关。

除海拔作用外，张玉琴等、李娟等、付建新等、魏建瑛等研究指出，坡度、坡向不同，区域内NDVI趋势变化也会出现差异。因此，NDVI趋势变化受海拔作用外，坡度、坡向可能也有所影响。

4.3 NDVI与极端气候指数相关

气温通过影响植物的碳素同化和消耗过程，从而控制植被生长，同时也通过调节区域蒸散过程控制土壤水分影响植被生长；降水则通过调控土壤水分间接影响植被变化。本研究发现，轿子山极端气温指数与NDVI、NDVI、NDVI的相关性整体强于极端降水指数，与Shao等的研究结果西南地区极端气温指数对NDVI、LAI、GPP的贡献率大于极端降水指数一致。进一步来看，NDVI对极端暖指数TXx、TR20、SU25呈显著正相关，而NDVI与暖指数TNx、WSDI呈显著负相关，表明NDVI对夜间和白昼气温响应存在显著差异，这与Peng等的研究结论昼夜升温存在对NDVI的非同步性作用一致。极端降水指数中，仅强度指数中的RX1、极端降水日指数中R25 mm分别与NDVI、NDVI呈显著正相关，表明降水强度、大雨日数对NDVI影响较大。然而，NDVI与极端降水指数均未通过显著性检验，表明极端降雨指数对NDVI的影响不明显。

Jiang等、Yin等、Wei等分别揭示了西南地区、横断山区、澜沧江流域不同植被类型NDVI在不同海拔带对平均气温/极端气温、降水/极端降水存在显著不同的海拔依赖性变化。因此，若轿子山不同海拔下可用气象站更多(如落雪、雪山等站持续监测或现自动站监测气象序列更长)，可进一步揭示轿子山不同海拔带不同植被类型对不同海拔带极端气候指数的响应。

气候变化已对中国西南气候脆弱性地区的生物多样性构成巨大威胁，保护区气候风险阈值的精确识别是气候风险评估与管理的基础。目前，我国自然保护区建设与管理的主要技术标准尚未纳入气候变化风险调查与管控的相关内容，缺乏成熟的极端气候变化风险评价和监测预警等技术，未来可设立自然保护区适应气候变化的专项资金，以加强自然保护区应对气候变化风险的管控能力建设。

5 结 论

(1)轿子山年平均气温(0.13 ℃/10 a)趋于显著升高，年日照时间(-58.73 h/10 a)、年平均风速(-0.21(m·s)/10 a)和年雾日(-4.8 d/10 a)均显著下降。

(2)20世纪50—80年代、2017—2019年轿子山基带温度带分别为南亚热带、北热带，相较于20世纪50—80年代，2017—2019年南亚热带攀升。

(3)1999—2016年，海拔1 736～2 236(0.07/10 a)，2 236～2 736(0.06/10 a)，2 736～3 236(0.05/10 a)，3 236～3 736 m(0.04/10 a)的区域NDVI均趋于显著增长，但增幅随海拔升高而减小。

(4)NDVI对暖指数、降水强度指数、极端降水日指数响应更为敏感，白昼气温升高促进/抑制NDVI/NDVI增长，夜间气温升高对NDVI产生较大影响。