植被是陆地生态环境的重要组成部分，在一定程度上能代表土地覆盖的变化，对全球能量传输、气候和环境变化具有重要的指示作用^[1-2]。植被覆盖度(fractional vegetation cover, FVC or C_FV)指单位面积内植被冠层垂直投影所占百分比^[3]，是衡量地表植被覆盖及区域内生态环境状况的重要指标^[4]。因此，植被覆盖演变特征研究能从一定程度上反映全球气候变化的趋势，是生态环境领域研究热点，对国家公园保护研究具有重要现实意义和学术价值^[5]。前人基于植被指数监测植被覆盖变化已取得大量成果，邓晨晖等基于归一化差值植被指数(normalized difference vegetation index, NDVI or I_NDV)研究秦岭植被覆盖时空变化特征^[6]；甄计国等基于3种植被指数探索兴隆山国家级自然保护区复杂地形调节下植被指数的年际变化特征^[7]；齐亚霄等基于NDVI研究天山北坡植被覆盖动态变化特征^[2]。

武夷山国家公园(Mount Wuyi National Park, MWYNP)是我国唯一一个既是国家公园体制试点，又是世界生物圈保护区，同时还是列入世界文化和自然双遗产名录的自然保护区，其保护研究价值在中国乃至全球都是不可替代的。众多学者在武夷山植物群落结构及其生态环境^[8-9]、社区居民^[10]、生态旅游^[11-13]、景观空间格局与评价、生态系统服务价值^[14]等领域对武夷山自然保护区和风景名胜区做了大量研究，但自武夷山国家公园成立以来，在全球生态系统和气候变化下植被覆盖演变特征研究尚鲜有报道。因此，以武夷山国家公园为研究对象，对2000年(武夷山加入“世界自然与文化遗产”的第1年)、2005年(2003、2004年出现大旱^[15]，2005年是旱情后的第1年)、2010年(武夷山自然保护区承前启后正式迈进国家公园阶段)、2015年(武夷山被发改委确定为我国首批国家公园体制试点之一)以及2018年[武夷山国家公园正式发布总体规划及专项规划(2017—2025年)初稿] 5个关键时间节点内近18 a植被覆盖变化的时空特征演变规律以及地形因子、生态系统服务价值等因子进行分析，为武夷山国家公园管理和生态恢复提供借鉴。

1 研究区概况

武夷山国家公园位于武夷山脉北部(27°31′20″~27°55′49″N，117°24′12″~117°59′19″E)，北与江西省铅山县交界，南至建阳区黄坑镇，西至光泽县崇仁乡，东至武夷山市武夷街道，包括东部武夷山国家级风景名胜区、中部九曲溪生态保护区、西部武夷山生物多样性保护区、光泽武夷天池国家森林公园及周边公益林、邵武市国有林场龙湖场部分区域以及城村闽越王城遗址等区域，总面积为1 027.8 km²。武夷山国家公园山高坡陡、峡谷纵深，最高主峰黄岗山高2 158 m，地势起伏大，地形复杂，海拔1 500 m以上的山峰有112座，地形坡度一般为0°~70°；相对海拔最高达2 000 m，随着海拔的递增、坡向和坡度变化，气温和降水、生物、土壤垂直变化明显，植被依次分布为常绿阔叶林带、针阔过渡林带、温性针叶林带、中山苔藓矮曲林带和中山草甸带5个垂直带谱；雨量充沛，年均降水量1 684~1 780 mm，年均相对湿度为78%，雾日在100 d以上，年均气温17~19 ℃。武夷山国家公园是全球同纬度带最完整、最典型、面积最大的中亚热带原生性森林生态系统，包含了我国中亚热带所有的植被类型，孕育着多种多样的植被类型^[16]，还有原生性森林植被290 km²，且未受到人为破坏，是我国亚热带东部地区森林植被保存最完好的区域^[17-19]。

2 研究方法 2.1 数据来源和预处理

根据2000—2018年武夷山区域气象站降水量、日照时间、极端温度等气象资料；《武夷山国家公园总体规划及专项规划(2017—2025年) 》；中国科学院计算机网络信息中心地理空间数据云平台(http://www.gscloud.cn)获取的Landsat影像图(2000、2005、2010、2015、2018年5个时期，影像轨道号均为120/041，成像时间均为植被主要生长季节4—5月且云量低于5%遥感影像)，对原始图像进行波段组合、几何校正、辐射校正、图像融合、图像裁剪和投影变换等预处理，使其具有时空可比性。

2015、2016年调查星村镇至桐木关定位信息，2017—2018年实地调查武夷山风景名胜区大王峰以及星村镇周围植被覆盖状况，共收集近500个GPS定位控制点，导入Google Map软件对2015年遥感影像选取均匀分布的30个控制点进行几何校正，均方根误差小于1个像元。以2015年Landsat-8 OLI/TIRS影像为基准，对其他年份影像配准校正，保证了后期植被覆盖度变化研究时，影像能保持一致。

2.2 植被指数

由于红光波段反射率对复杂地形比较敏感，而传统植被指数忽略山区像元级光谱信号受到植被环境的特殊性和差异性等影响。因此，针对高植被覆盖且地形复杂的武夷山国家公园，采用吴志杰等^[20]提出的归一化差值山地植被指数(normalized difference mountain vegetation index, NDMVI or I_NDMV)、LIU et al^[21]提出的增强植被指数(enhanced vegetation index, EVI or I_EV)以及传统的NDVI等植被指数，然后通过像元二分法模型，来对比估算武夷山国家公园的植被覆盖度，具体植被指数如下。

2.2.1 归一化差值植被指数

NDVI对绿色植被的生长势和生长量表现非常敏感，与植被密度呈正相关，能较好地反映地表植被的茂盛程度，是监测植被生长状况和植被覆盖变化的指示因子，其公式为：

$ I_{\rm NDV}= \frac{{f_{\rm NIR}-f_{\rm RED}}}{{ f_{\rm NIR}+f_{\rm RED} }} $

(1)

式中：f_NIR和f_RED分别代表红光波段和近红外波段的地表反射率。I_NDV∈ [-1, 1]，由于f_RED在高植被覆盖区容易饱和，这时只有f_NIR反应植被变化，I_NDV≈1。

2.2.2 归一化差值山地植被指数

为了对复杂山地的地形影响进行校正，吴志杰等^[20]通过红光波段和近红外波段地表反射率的最小值调整NDVI中红光波段和近红外波段的地表反射率，使得具有相同植被覆盖状况的坡面在植被指数估算值上接近，从而在不需要数字高程模型辅助下，构建适合山地复杂地形的植被指数模型，其公式为：

$ I_{\rm NDMV}= \frac{{ ( f_{\rm NIR}-n ) - ( f_{\rm RED}-r ) }}{{ ( f_{\rm NIR}-n ) + ( f_{\rm RED}-r ) }} $

(2)

式中：n为近红外波段波段调整参数，n=f_{NIR, min}；r为RED波段调整参数，r=f_{RED, min}；2018年影像数据n=11；r=29。

2.2.3 增强植被指数

为了弥补NDVI植被指数的缺陷，LIU et al^[21]提出增强植被指数EVI，该指数是通过引进蓝光波段校正红光波段中的大气干扰以增强植被监测效果，利用土壤调整参数修正土壤背景影响，减少植被冠层背景噪声影响且不易产生过饱和现象^[22]，从而达到优化高植物量区域植被信号，公式为：

$ I_{\rm EV}= G \frac{{f_{\rm NIR}-f_{\rm RED}}}{{ f_{\rm NIR}+C_{\rm 1}f_{\rm RED}-C_{\rm 2}f_{\rm BLUE}+L }} $

(3)

式中：f_BLUE代表蓝光波段的地表反射率；G是增益系数；C₁，C₂是气溶胶阻抗系数，使用蓝光波段校正红光波段中的大气干扰；L是土壤调整参数；根据JIANG et al^[23]研究各参数取值分别为：G=2.5、C₁=6、C₂=7.5、L=1。

2.3 植被覆盖度估算方法

采用FVC像元二分模型中适合高植被覆盖区的GUTMAN et al^[24]模型，由公式(4)表示:

$ C_{\rm FV}= \frac{{S-S_{\rm soil}}}{{ S_{\rm veg}-S_{\rm soil}}} $

(4)

式中：S表征像元的植被指数信息，S_soil表征影像中完全无植被覆盖部分(裸地)的植被指数值，S_veg表征完全由植被覆盖部分的植被指数值。

由于武夷山国家公园是经过近30 a生态保护工作的中亚热带原生性森林生态系统，很难选取完全无植被覆盖的天然裸土地，所以参考前人研究基础上^[25-26]，采用在置信区间0.5%~99.5%内获取最大值和最小值，将NDVI、NDMVI和EVI的计算结果分别代入(4)中，以EVI为例可得：

$ C_{\rm FV,EVI}= \frac{{ I_{\rm EV}-I_{\rm EV,\min}}}{{I_{\rm EV,\max}-I_{\rm EV,\min}}} $

(5)

式中：C_{FV, EVI}代表EVI的植被覆盖度；I_{EV, min}为置信区间0.5%~99.5%内最小值，表征纯裸土或无植被区的I_EV；I_{EV, max}为表征高植被覆盖区的I_EV；C_{FV, NDVI}，C_{FV, NDMVI}用同样的方法计算。

为了使变化结果更加直观，将植被覆盖度以20%为间隔分为5个等级制作植被覆盖等级图，分别为：[0.0, 0.2)(低植被覆盖区)、[ 0.2, 0.4)(中低植被覆盖区)、[ 0.4, 0.6)(中植被覆盖区)、[ 0.6, 0.8)(中高植被覆盖区)和[ 0.8, 1.0] (高植被覆盖区)。

2.4 植被覆盖指数突变分析

Mann-Kendall (M-K)是由Mann和Kendall所提出的一种非参数检验法，常用于分析具有时间序列特征的数据，样本不需要遵循特定分布，且很少受到异常值干扰^[27]。对具有n个样本量的时间序列(x₁, x₂, x₃，…，x_n)，构造秩序列S_k：

$ S_{ k}= \sum\limits^k_{ i=1} R_{ i}，(k=1,2，3，…，n) $

(6)

式中：k是数据集的长度(年份)；R_i为约束条件x_i>x_j (1≤j≤i)的累计数。

$ R_{ i}= \left\{ \begin{array}{l} 1，x_{ i}＞x_{ j} \\ 0，x_{ i}≤x_{ j} \end{array} \right. ,(j=1,2，3，…，i) $

(7)

将随机独立的时间序列数据构造统计量U_Fk：

$ U_{ F_{k}}= \frac{{ S_{ k}-E_{ S_{k}}}}{{ \sqrt{V_{ S_{k}} )}}} ，(k=1,2，3，…，n) $

(8)

式中：当k=1时，U_F1=0。E_Sk表示均值，V_Sk表示方差，其公式为：

$ E_{ S_{k}}= \frac{{ n ( n-1 )}}{{ 4}} $

(9)

$ V_{ S_{k}}= \frac{{ n ( n-1 ) ( 2n+5 )}}{{72}} $

(10)

同时按将时间序列逆序后构造x_n ′，代入公式(8)计算逆序统计量U_Bk，且U_Bk=-U_Fk，当构造的两条曲线出现交点时，且交点在95%置信区间内，则说明交点对应的是突变时刻。

3 结果与分析 3.1 3种植被指数估算植被覆盖度

选取最临近的2018年Landsat遥感影像，根据公式(1) ~ (5)分别计算武夷山国家公园NDVI、NDMVI和EVI植被指数，选取置信区间0.5%~99.5%内的最小值代表纯裸土像元和最大值代表高植被覆盖像元，分别计算对应的植被覆盖度C_{FV, NDVI}、C_{FV, NDMVI}、C_{FV, EVI}，通过ENVI 5.3分别统计3种植被指数估算植被覆盖度标准差，排序为NDVI (0.12)＜EVI (0.13)＜NDMVI (0.18)，其实3种植被指数估算方法都能达到0.1，说明估算效果都很好，排序结果只是说明NDVI估算结果更集中，而NDMVI估算结果波动幅度更大。

通过统计每个像元的植被覆盖度数量图进一步比较3者优劣，从图 1(a) 可知，NDVI在植被覆盖度达到0.796后收敛，说明这时RED波段已经饱和，估算结果趋近于1，所以NDVI估算标准差最小；同时图 1(b) NDMVI估算结果也出现了同样现象，但NDMVI通过波段调整参数使得具有相同植被覆盖状况的坡面在植被指数估算值上接近，所以NDMVI在植被覆盖度达到0.699才开始收敛；图 1(c) EVI增强植被指数从图像上看服从正态分布，没有出现急剧收敛的现象，能够较好的表达出从0.6~1.0范围内的植被覆盖，估算效果能较好的优化高植被覆盖区植被信号的作用。因此选取EVI增强植被指数来估算武夷山国家公园近18 a的植被覆盖演化状况。

3.2 不同时期植被覆盖时间格局演化特征

根据公式(3)和公式(5)分别计算EVI和FVC_EVI，得到各年份武夷山国家公园植被覆盖变化趋势图[图 2(a) ]，并通过MATLAB软件实现基于Mann-Kendall (M-K)的植被指数突变分析[图 2(b) ]。为了便于分析2000—2018年的降水量、日照时间以及极端高温的气象数据，对数据做标准化处理后得到武夷山国家公园气候条件变化图[图 2(c)]。

由图 2(a) 可知2000、2005、2010、2015和2018年间，各年份植被覆盖度均值分别为为94.97%，90.73%，93.55%，98.10%，98.02%。武夷山国家公园平均植被覆盖度均是0.8以上的高覆盖植被区域，而且从2005年起植被覆盖变化呈显著上升趋势，增速为0.29% · a^-1 (P＜0.01)，说明近18 a来武夷山国家公园在不断改善总体植被覆盖状况。

由图 2(c) 可知，2000—2005年间，在2003、2004年间降水量达到最低点，同时日照时间和极端最高温都出现了最高点，金保明^[15]研究发现在2003、2004年南平市出现高温持续时间长、降雨量大幅度减少、蒸发量大的极端气候，武夷山市达到大旱标准，导致武夷山国家公园内大量植被死亡，所以植被覆盖度从2000年的为94.97%下降到2005年的90.73%；同时，可发现2011年武夷山国家公园降雨量达到次低点，2012年出现了极端低温现象。所以，从图 2(b) Mann-Kendall (M-K)突变检验可知，近18 a来武夷山国家公园的年均植被覆盖度突变发生于2012年，且有升高的趋势，说明武夷山国家公园生态修复工作从2005年开始，在2010年国家公园发展进入新阶段后，最终在2012年成效趋于稳定且越来越好。

3.3 不同时期植被覆盖空间格局演化特征

通过EVI植被指数估算武夷山国家公园植被覆盖度(图 3)，2000—2018年，武夷山国家公园内植被覆盖的空间格局变化相似，表现为西北高、东南低的空间格局。从图 3可知，2000、2005、2010、2015、2018年高植被覆盖区面积分别为923.16、855.51、911.12、987.24、985.33 km²，分别占整个国家公园面积比例为90.22%、83.63%、89.04%、96.54%、96.30%，高植被覆盖区主要分布在西北区域，大部分为原国家级自然保护区区域，一直以来都为生态红线区，是珍稀、濒危动植物的集中分布区域，严禁砍伐，植被保存完好，没有人为干扰，是园区内最核心部分，其覆盖面积之广在武夷山国家公园内占据绝对优势，因此生态环境得到很好的保护。

从表 1可以看出，2000—2005年武夷山国家公园高植被覆盖区减少了67.65 km²，相比2000年下降了7.32%，而中高植被覆盖区、中植被覆盖区以及中低植被覆盖区分别增加了24.86、18.72和18.45 km²。从图 3(b) 可知，植被覆盖等级下降主要由于发生在星村镇沿着九曲溪一带以及部分的风景名胜区，说明在这个期间内植被覆盖结构发生较大变化，才呈现出东南区域植被覆盖低的空间格局。

但是从2005年开始植被覆盖度逐年上升，从图 3可知，2005—2018年，武夷山国家公园内高植被覆盖区增加129.82 km²，而中高植被覆盖区、中植被覆盖区以及中低植被覆盖区分别减少57.49、40.66和25.64 km²，主要是九曲溪靠近星村镇附近沿着溪流两岸植被覆盖度不断提高，得益于持续开展水源涵养、河流源头区的保护以及星村镇周边地区的违规茶山的生态恢复工作。

3.4 不同时期功能分区植被覆盖空间格局演化

由于武夷山自然保护区与风景名胜区分属于两个不同管理单位，当武夷山国家公园成立时，将二者归并为武夷山国家公园管理局，故需要对园区范围内重新进行功能分区。《武夷山国家公园总体规划及专项规划(2017—2025年) 》结合武夷山自然保护区和风景名胜区规划，按照武夷山自然生态系统原真性和完整性、物种栖息地连通、保护管理统一的原则，统筹生态保护和利用现状，合理归并后确定了武夷山国家公园的管控和功能分区，实行差别化管理，优化功能分区的边界范围。选取其中核心保护区、生态修复区和传统利用区3个功能分区进行植被空间变化研究(图 4)。其中，核心保护区面积为423.85 km²，占总面积41.42%，主要是武夷山自然保护区内常绿阔叶林带的禁止开发区域；生态修复区面积为358.85 km²，占总面积35.07%，主要以武夷山风景名胜区、部分九曲溪上游保护带等为主的区域，通过将违规开垦的茶山恢复为针阔叶林；传统利用区面积为240.60 km²，占总面积23.51%，包括以星村镇、桐木村、坳头村为主的传统村庄区域。

为了进一步分析武夷山国家公园内不同功能分区下的植被覆盖演变特征，采用ArcGIS软件计算2000—2018年不同功能分区的植被覆盖率(表 2)，并通过影像差值法来揭示武夷山国家公园植被覆盖度的空间演化(图 5)。从表 2和图 5中可知，核心保护区植被覆盖度总体保持在97%以上，但是在2000—2005年期间，在自然保护区南部诸母岗以及北部黄岗山山顶一带出现植被覆盖度下降情况；而从2005年开始，核心保护区的植被覆盖度就稳步上升，直至2015—2018年保持在99%稳定状态，这是武夷山国家公园在该时期内保护工作的最终体现。

传统利用区和生态修复区植被覆盖度主要变化均发生在2000—2005年，图 5(a) 可看到生态修复区大面积下降，这也充分反映当时大旱气候条件以及毁林种茶严重破坏武夷山的生态环境。自2005年开始，高植被覆盖度的区域面积不断增加，中低植被覆盖度(图 3黄色及红色区域)逐渐减少，说明在2015—2018年武夷山成为国家公园体制试点期间，武夷山国家公园管理局在生态恢复方面做了大量工作，园区内生态环境稳步优化。

4 讨论与结论

运用植被指数方法估算武夷山国家公园在2000、2005、2010、2015和2018年5个关键时期的植被覆盖度，揭示了不同发展阶段下武夷山国家公园植被覆盖演变时空特征以及不同功能分区下的演变规律，结果表明：

(1) 以2018年Landsat遥感影像为例，EVI结果图像服从正态分布，能够较好表达植被覆盖，估算效果能较好的优化高植物量区域植被信号的作用。

(2) 植被覆盖总体演变时空格局分析结果表明，2000—2018年间，武夷山国家公园内植被覆盖空间格局极为相似，表现为西北高、东南低。

(3) 2000—2005年植被覆盖度总体水平下降，主要发生在星村镇沿着九曲溪一带以及部分的风景名胜区，根据《武夷山风景名胜区总体规划(2000—2010年) 》规划，该时段武夷山存在问题主要有：由于2003—2004年南平市气候变化迅速，武夷山市达到大旱标准，导致武夷山国家公园内大量植被死亡；景区内出现毁林种茶现象，破坏生态环境，人为干扰活动严重；九曲溪源头是九曲溪水赖以存在关键，但九曲溪上游保护地带的森林砍伐现象屡禁不止，已出现九曲溪水位下降的趋势。这些问题不仅说明村民没有意识到生态环境对九曲溪以及武夷山自然景观的重要性，且只注重眼前利益，还验证了本研究符合当时的发展规律，从而造成武夷山国家公园东南区域植被覆盖低的空间格局。

但2005—2018年期间内，植被覆盖度逐年上升；到2018年，武夷山国家公园高植被覆盖区增加了129.82 km²，而中高植被覆盖区、中植被覆盖区以及中低植被覆盖区分别减少了57.49、40.66和25.64 km²，主要是九曲溪靠近星村镇附近沿着溪流两岸植被覆盖度不断提高，这得益于持续开展水源涵养、河流源头区的保护以及星村镇周边地区的违规茶山的生态恢复工作。

(4) 武夷山国家公园核心保护区的植被覆盖度常年保持在97%以上，2000—2005年期间，自然保护区南部诸母岗以及北部黄岗山山顶一带出现植被覆盖度降低情况，邱堋星^[28]认为黄岗山山顶草甸形成是由于人为干扰或自然火灾作用，破坏了原生性植被，而在山顶恶劣气候环境下，只能生长高耐力、匍匐的草本和灌木植被；任引^[19]研究表明武夷山常绿阔叶林占全区森林面积的1/4，且乔木层的物种多样性和均匀度均高于草本、灌木层。因此，黄岗山山顶的物种多样性相对海拔1 600 m以下的乔木层是较差的，较难抵抗大旱极端气候，所以该期间内植被覆盖度大幅度下降；而从2005年开始，核心保护区的植被覆盖度就稳步上升，直至2015—2018年期间保持在99%稳定状态，说明武夷山国家公园越来越重视核心保护区的生态维护工作。

总而言之，有效利用遥感信息技术对国家公园生态资源实时监测与评价，加强针对性的生态修复，开展不损害生态系统的社区生活生产设施改造和公众游憩、科教娱乐项目，为我国东部人口密度较高的森林生态系统保护和发展，打造国家级“智能化、智慧化”示范园，使具有中国特色的国家公园规模和管理达到世界先进水平，真正实现对山水林田湖草的完整保护。