武夷山国家自然保护区NDVI对坡度的响应

张玉琴 ，李丽纯 *，彭继达 ，廖廓 ，刘桐恺

1. 福建省灾害天气重点实验室，福建 福州 350001；2. 福建省气象科学研究所，福建 福州 350001；3. 福建省气象台，福建 福州 350001

自然保护区是指对有代表性的自然生态系统、珍稀濒危野生动植物物种的天然集中分布和有特殊意义的自然遗迹等保护对象所在的陆地、陆地水域或海域，依法划出一定面积予以特殊保护和管理的区域（祝萍等，2018）。各类自然保护区是保护生物多样性及其生境最有效的途径之一，已成为全球自然资源保护基石，对促进国家的国民经济持续发展和科技文化事业发展具有十分重大的意义（杨绘婷等，2020），但自然因素和人类活动常常对生态系统造成威胁。武夷山国家自然保护区气候环境复杂，是西风带系统东进南下和东风波系统西进北上的必经之路，主要灾害性天气为暴雨，地势高，起伏大，坡度陡峭，短时强降水常导致山体滑坡和泥石流，影响植被生态质量。坡度是描述地形与地貌形态的重要指标，会影响气候条件的空间分布，并在一定程度上限制人类活动，也决定了地表植被的立地条件（如土壤结构、土壤水分和土壤厚度等），进而影响植被生长恢复（仙巍等，2007；梁志权等，2015）。研究坡度与植被的定量关系，揭示在植被长期生长适应过程中坡度如何影响其时空分布，可提高植被参数的建模能力，对保护区植被保护与管理具有指导意义。

现有坡度对植被影响的研究主要集中在植被覆盖的地形分异（潘霞等，2019；马士彬等，2019）和坡度对退耕和水土流失的影响（卢爱刚等，2010；王超等，2014）等，也有研究定性分析了坡度对植被变化的影响（王钊等，2018；王瑾等，2019）。武夷山国家自然保护区作为中亚热带森林生态系统的典型代表，区内坡度差异大，许多学者针对其森林生态系统结构和功能（王邵军等，2009；王国兵等，2011）、物种多样性（封磊等，2004；王同亮等，2015）和植被分布（汪荣等，2007；邱炳文等，2009）等展开研究，但对保护区植被覆盖与坡度关系的定量研究未见报道。

归一化植被指数（Normal Difference Vegetation Index，NDVI）是表征植被状况的常用指标（Choudhury et al.，1994；Gillies et al.，1995），已有很多学者基于NDVI分析植被状况，以及自然因素（地形、气候条件等）、人类活动等要素变化对植被的影响（邹伟成等，2015；张佳琦等，2019；罗新兰等，2020；刘海等，2020）。本研究采用长时间序列Landsat卫星遥感资料，提取归一化植被指数（NDVI），揭示2004—2018年武夷山国家自然保护区平坡、缓坡、斜坡、陡坡、急坡、险坡 6种坡度下NDVI的分布特征和变化趋势，以及不同坡度下植被的改善和退化情况，阐明坡度对植被生长的长期影响，并指出影响各坡度植被生长的主要气象因子，以深化对保护区的植被生长变化情况的了解，为保护区制定生态保护规划、实施针对性地保护和应对气候变化影响提供参考。

1 研究区概况与数据来源

1.1 研究区概况

武夷山国家自然保护区（117°27′—117°51′E，27°33′—27°54′N）位于福建省北部和武夷山脉北端，地跨武夷山、建阳、光泽和邵武四市（县），总面积56527.3 hm2（图1）。该保护区成立于1979年，被联合国教科文组织接纳为国际生物圈保护区，是全球同纬度带保存最为完整、面积最大的典型中亚热带原生性森林生态系统。区内植被类型丰富，森林覆盖率超过 95%。保护区具有地势高，起伏大，多垭口的地貌特征。平均海拔 1200 m，最高处达2158 m，最低处仅300 m，相对高差达1858 m，高差极为悬殊；坡度一般为20°—30°，最陡达70°。

图1 武夷山国家自然保护区区划及坡度图Fig. 1 Geographical location and slope map of Wuyishan National Nature Reserve

1.2 数据来源与预处理

综合考虑保护区数据的图像分辨率、可获取性、云覆盖程度和影响成像日期等，选取 2004—2012年 Landsat-7 ETM+以及 2013—2018年Landsat-8 OLI遥感影像为基础数据源（https://earthexplorer.usgs.gov/），原始遥感影像已做过系统级辐射校正和几何粗校正，空间分辨率为30 m，时间分辨率为16 d，轨道号为120/41，时间跨度为2004—2018年，共60景遥感影像。由于Landsat-7卫星机载扫描行校正器在2003年5月出现故障，因此在对 Landsat-7 ETM+影像进行预处理前，先使用ENVI 5.3的去条带补丁landsat_gapfill对其进行去条带处理，再结合ENVI 5.3和ArcGIS 10.2对数据先后进行辐射定标、大气校正、投影转换和剪裁等预处理。

数字高程模型（DEM）数据为中国范围SRTMDEM数据，数据来源于中国科学院计算机网络信息中心地理空间数据云平台（http://www.gscloud.cn/），空间分辨率为90 m。借助ArcGIS 10.2空间分析功能计算得到保护区坡度分布图。

气象观测数据为国家基本气象站武夷山、建阳、光泽和邵武站的气温、降水月值数据，时间跨度为 2004—2018年，数据来源于福建省气象局。将4站平均后的气温、降水月值作为保护区的气温、降水月值，再基于月数据形成保护区年均温、年累积降水量逐年时间序列。

2 研究方法

2.1 遥感数据预处理

基于同年四季多景影像计算 NDVI值，采用最大值合成法得到该年NDVI均值（田庆久等，1998；Lovell et al.，2001）。采用时间序列谐波分析（HANTS）对数据进行滤波，重构NDVI序列，解决福建云雨影响导致的 NDVI数据波动较大，无法进行趋势分析的问题。HANTS是基于傅立叶分析改进的一种有效的算法，它将时间序列数据分解为几个最为显著的频率，采用最小二乘法拟合各频率的余弦或正弦分量，其中各频率的振幅和相位都是通过不停的迭代得到，其结果能够有效真实地反映曲线周期性变化（Julien et al.，2006；潘学鹏等，2015）。

2.2 样本选择及坡度分级

在保护区范围内随机均匀选择样点 2000个（图2），控制各样本互不重叠，并将2000个样本纳入回归模型。统计各样本区域NDVI变化趋势及平均坡度，建立二者回归模型，分析二者关系。同时，将坡度分为6个等级（汪荣等，2007），即平坡（0°—5°）、缓坡（5°—15°）、斜坡（15°—25°），陡坡（25°—35°）、急坡（35°—45°）、险坡（＞45°），分析不同坡度NDVI变化趋势。

2.3 趋势线分析

图2 武夷山国家自然保护区坡度和研究样本分布图Fig. 2 Slope and research sample distribution map of Wuyishan National Nature Reserve

对NDVI进行一元线性回归处理，以逐像元的时间变化规律来反映保护区整体的空间变化规律。回归方程斜率反映了植被覆盖度在研究时段的变化趋势，若斜率小于0说明该时段内植被呈减少趋势，反之则为增加趋势（Li et al.，2016）。计算公式如下：

式中，yNDVI为NDVI的变化趋势（即斜率）；n为研究时段的年数，即15年；xj为第j年保护区NDVI值。本文参考谭学玲等（2018）对植被覆盖变化的分级标准，结合武夷山国家自然保护区实际情况，将植被覆盖变化分为4个等级，即明显改善（yNDVI≥0.01），轻微改善（0＜yNDVI≤0.01），轻微退化（-0.01＜yNDVI≤0.01），严重退化（yNDVI≤-0.01）。

2.4 地形面积差异修正系数

地形面积差异修正系数（k）用于消除各地形因子条件下因地形绝对面积不同而引起的植被恢复评价的不确定性（谭学玲等，2018）。计算公式如下：

式中，ΔVi/Ti为特定地形条件下某植被变化类型的面积比；ΔVi为某植被变化类型（明显改善、轻微改善、轻微退化、严重退化）在特定坡度条件（平坡、缓坡、斜坡、陡坡、急坡、险坡）下所占的面积；Ti为特定坡度条件的面积；ΔV/S为保护区某植被变化类型的总面积；ΔV为某植被变化类型的总面积，S为武夷山国家自然保护区总面积。k＞1表明在该特定坡度条件下保护区内此植被变化类型为增强趋势，k＜1则相反，k=1表明此植被变化类型在该特定坡度条件下的分布与保护区内该植被变化类型的分布比例相同，其变化趋于平稳。

2.5 滑动平均

滑动平均是趋势拟合的基础方法，用确定时间序列的平滑值来显现出变化趋势，滑动平均序列计算公式如下（魏凤英，2007）：

式中，k为滑动长度，本文将2004—2018保护区年均温、年累积降水量和年NDVI序列进行3年滑动平均计算，有助于减小序列的波动性，显示出整体变化的趋势，便于进行相关性分析。

3 结果与分析

3.1 NDVI的变化特征

计算2004—2018年逐年不同坡度的NDVI均值（图3），可见，2004—2018年不同坡度NDVI年际变化都呈上升趋势，其中，险坡（＞45°）NDVI上升趋势最为明显（sig.＜0.05），趋势倾向率为0.044/10 a，急坡（35°—45°）次之（sig.＜0.1），趋势倾向率为0.034/10 a，其余坡度NDVI上升趋势均不明显，趋势倾向率为0.018—0.025/10 a。15年里有2个上升高峰期，分别为2006—2009年和2012—2015年。各坡度年均NDVI也存在差异，整体而言，陡坡（25°—35°）NDVI值最大，平坡（0°—5°）NDVI值最小，最低值为2006年，最高值为2009年。2006年保护区NDVI主要集中于0.7—0.9，平均值为0.81；2009年NDVI主要集中于0.8—1.0，平均值为0.91，说明武夷山国家自然保护区植被状况良好，总体维持高植被覆盖状况（图4）。

图3 2004—2018年武夷山国家自然保护区不同坡度NDVI均值Fig. 3 Annual average NDVIs on different slopes in Wuyishan National Nature Reserve from 2004 to 2018

图4 武夷山国家自然保护区2006（a）和2009（b）年NDVI空间分布Fig. 4 Spatial distribution map of NDVI in Wuyishan National Nature Reserve in 2006 (a) and 2009 (b)

图5 2004—2018年武夷山国家自然保护区NDVI年际变化趋势空间分布Fig. 5 Spatial distribution map of NDVI interannual variation trend in Wuyishan National Nature Reserve from 2004 to 2018

逐像元计算NDVI年际变化趋势率（图5），结合坡度分析发现，2004—2018年武夷山国家自然保护区植被覆盖总体呈现趋好态势，以东部和西部高坡度地区增加最为明显，中部和北部低坡度区呈现减少趋势。其中，2.2%的区域NDVI值呈现明显改善趋势，57.2%区域呈现轻微改善趋势，40.4%的区域NDVI值呈现轻微退化趋势，0.2%区域呈现严重退化趋势。武夷山国家自然保护区自建立起，各级人民政府相继出台了多项政策和管理办法，如严格限制外来人口迁入、控制人口进出等，减轻人类活动对自然资源的压力，使得植被覆盖度总体呈上升趋势，在生态环境保护方面取得良好成效。

3.2 NDVI变化坡度分异

为了分析不同坡度 NDVI变化趋势的差异，计算了各坡度上NDVI变化趋势的统计量（表1）。从坡度范围NDVI变化趋势均值来看，NDVI变化趋势最明显的坡度范围为险坡（＞45°），其次为急坡（35°—45°），0°—25°的区域 NDVI变化趋势都较不明显。从标准差来看，在0°—25°范围内，标准差随坡度增加而减小，说明 NDVI变化趋势分布随坡度增加呈集中趋势，即坡度较小区域 NDVI变化趋势波动较大；在＞25°范围内，标准差随坡度增加而增加，说明 NDVI变化趋势分布随坡度增加呈分散趋势，即坡度较大区域NDVI变化趋势波动较大。

表1 2004—2018年武夷山国家自然保护区不同坡度NDVI变化趋势统计量Table 1 Statistics of NDVI variation trends from 2004 to 2018 on different slopes in Wuyishan National Nature Reserve 10-4

为了消除不同坡度绝对面积不同而引起的植被变化评价的不确定性，还计算了不同坡度植被变化类型比和修正系数k值（表2、3），发现坡度对植被退化或改善的影响程度差异较大。在所有坡度范围内，植被变化类型均以轻微改善为主（面积比超过 50%），但在 0°—15°范围内，植被退化类型面积比随坡度增加而增加，说明在此坡度范围内坡度增大不利于植被生长；而在＞15°范围内，随坡度增加改善类型面积比增加，说明在此坡度范围内坡度越大越有利于植被生长。同时，随坡度增大，植被退化类型的k值逐渐减小，植被轻微改善类型的k值逐渐增大，植被明显改善类型的k值先减小后增大；其中，0°—25°内植被严重退化、轻微退化类型的k＞1，表明在0°—25°范围内利于植被退化，而在＞35°范围内植被轻微改善、明显改善类型的k＞1，表明在＞35°范围内利于植被改善，以0°—5°范围内植被严重退化类型和＞45°范围内植被明显改善类型趋于增强最为明显。上述分析发现，坡度明显地影响着武夷山国家自然保护区植被覆盖变化，植被生长恢复时要充分考虑坡度因子。

表2 2004—2018年武夷山国家自然保护区不同坡度植被变化类型面积比Table 2 Area ratio of vegetation variation types on different slopes in Wuyishan National Nature Reserve from 2004 to 2018

表3 2004—2018年武夷山国家自然保护区不同坡度植被变化类型修正系数k值Table 3 Correction coefficient k of vegetation variation types on different slopes in Wuyishan National Nature Reserve from 2004 to 2018

3.3 NDVI变化趋势与坡度的回归分析

统计2000个样本NDVI变化趋势及坡度，制作散点图（图6），发现NDVI变化趋势与坡度呈二次函数关系，以保护区2004—2018年NDVI变化趋势为因变量，坡度为自变量，建立回归模型。该模型F值为28.86，在0.001水平上显著，该回归模型揭示了坡度与NDVI变化趋势的关系，坡度在0.001水平上显著影响 NDVI的变化趋势。关系式如下：

式中，y为NDVI变化趋势，x为坡度。当坡度约17°时，NDVI变化趋势达到最小值；当坡度小于17°时，NDVI变化趋势随着坡度的增加而缓慢减少；当坡度大于17°时，NDVI变化趋势随着坡度的增加而增加，与基于地形面积差异系数的植被变化情况分析结果相似，高坡度区域植被改善程度最好。

3.4 气象因子对不同坡度NDVI的影响

图6 武夷山国家自然保护区各样本2004—2018年NDVI变化趋势和坡度的散点图Fig. 6 Scatter plot of NDVI variation trend from 2004 to 2018 and slope for each sample in Wuyishan National Nature Reserve

气候变化是引发地表覆盖变化的主要因子，地表植被类型的时空变化是对长期气候的适应，也是对短期气候波动的响应，影响植被变化的气象因子主要是气温和降水（Wang et al.，2011；韦振锋等，2013）。将2004—2018年年均温、年累积降水量和不同坡度年 NDVI值进行 3年滑动平均，得到2005—2017年滑动平均值。不同坡度年 NDVI滑动平均值与年均温、年累积降水量的Pearson相关分析表明（表4），全区及不同坡度NDVI与年均温均无显著相关，除险坡（＞45°）外，均为负相关；全区及＞25°范围 NDVI与年累积降水量呈显著正相关，0°—25°内呈不显著正相关。表明在保护区内，尤其在＞35°范围内，影响植被生长的主要气象因子为降水，2004—2018年保护区不同坡度降水量均呈明显上升趋势（sig.＜0.01），有利于植被覆盖度增加。

表4 武夷山国家自然保护区不同坡度植被变化和气象要素的相关分析Table 4 Correlation analysis between NDVI on different slopes and meteorological factors in Wuyishan National Nature Reserve

4 讨论

研究采用地形面积差异修正系数分析了6种坡度下植被的退化和改善程度，在 0°—25°范围内利于植被退化，而在＞35°范围内利于植被改善。有研究也得到坡度较高地区相较缓坡地区NDVI变化趋势较高的相似结论（聂斌斌等，2016；李薇等，2017；刘梁美子等，2018）。同时，建立了NDVI变化趋势与坡度的关系模型，进一步细化不同坡度下植被的改善程度，以便为有关部门开展针对性植被保护提供参考。中坡度区域改善程度最低，高坡度和低坡度区域植被改善程度较高，可能与中坡度区域植被指数高，容易达到饱和有关（李红军等，2007；陈燕丽等，2014），改善程度有限，致使植被指数变化趋势不明显。本研究从数学角度分析了保护区NDVI与坡度的关系，缺乏对生态系统过程的生理生态认知，该结果只能初步反映坡度对植被变化的影响，要获取较为准确的信息还需要通过开展生理生态试验和野外实地观测等进一步明确。

影响保护区内植被生长的主要气象因子为降水，与现有研究结果（顾婷婷等，2009；于泉洲等，2015；叶尔纳尔·胡马尔汗等，2019）相同。武夷山国家自然保护区以森林植被为主，现有研究对森林植被与气象因子的关系研究仍存在较多争议，有研究认为植被覆盖与气温因子的相关性要比降水因子大（徐凯健等，2015；周梦遥等，2018），需要更多高密度气象数据和植被生理生态的相关研究进一步明确。武夷山国家自然保护区海拔梯度明显，可能存在小气候差异，但保护区内无国家基本气象站，研究使用保护区周围武夷山、建阳、光泽和邵武4个国家站气象观测数据平均值作为保护区气象数据进行使用，可能无法准确反映小气候差异。近年来开始建设高密度、多梯度区域自动站，但时间序列仍不够长，随着时间的推进，未来可考虑长时间序列区域站数据进行研究分析。

2004—2018年武夷山国家自然保护区植被覆盖总体呈现趋好态势，与杨绘婷等（2020）研究结果一致，多归因于当地各级人民政府出台的多项生态保护政策和管理办法，以及近 15年降水增多，有利于植被覆盖度增加。但保护区仍存有耕地、居民点和生产经营活动，尤其低坡度区人类活动相对较多，中部和北部低坡度地区植被呈变差趋势，建议适当控制人类活动强度，做好长期规划。同时，NDVI在一定程度上能反映植被的健康状况，但要评估自然保护区的植被状况还需建立更为全面的区域植被生态质量评估模型，加入更多评估指标（植被指数、气温、降水、海拔、坡度、坡向、大气湿度和地表温度等），准确认清区域植被生态质量状况，这是未来研究工作中的重点所在。

5 结论

本研究基于2004—2018年Landsat数据、数字高程模型（DEM）数据和气象观测数据，对武夷山国家自然保护区NDVI的变化及其对坡度和气象因子的响应进行了分析。

（1）2004—2018年武夷山国家自然保护区植被覆盖总体呈现趋好态势，以东部和西部高坡度地区增加最为明显，中部和北部低坡度地区植被呈变差趋势。

（2）在所有坡度范围内，植被变化类型均以轻微改善为主（面积比超过 50%），在 0°—25°范围内利于植被退化，而在＞35°范围内利于植被改善。

（3）NDVI变化趋势与坡度呈二次函数关系，当坡度约17°时，NDVI变化趋势达到最小值；当坡度小于17°时，NDVI变化趋势随着坡度的增加而缓慢减少；当坡度大于17°时，NDVI变化趋势随着坡度的增加而增加。

（4）影响保护区植被覆盖主要气象因子为降水，近15年降水增多有利于植被覆盖度增加。