抚仙湖流域植被覆盖的时空演变规律

Temporal and Spatial Evolution of Vegetation Cover in Fuxian Lake Basin, Yunnan Province

摘要: 为揭示上世纪70年代以来抚仙湖流域内植被覆盖的时空演变，以抚仙湖流域1988年~2018年的长序列旱季遥感影像为基础，采用“3S”技术、Sen’s趋势分析、Mann-Kendall非参数检验和小波分析等方法，研究了流域内植被覆盖的时空演变规律，并分析了其影响因素。结果显示：1) 抚仙湖流域植被指数年际间波动较大，以0.012/10a的速率上升，2001年流域内植被覆盖发生了突变，5年为植被指数周期变化的第一主周期。2) 流域内植被覆盖增加区、减少区的面积为615.93 km2、218.85 km2，分别占陆地面积的73.78%、26.22%，增加区集中在山坡而减少区则集中于城镇建设区。3) 抚仙湖流域植被覆盖与海拔间存在正相关关系，1988年~2018年核心保护区、开发区、控制开发区和禁止开发区的变化值分别为0.016、−0.019、0.085和0.035；4) 流域内植被覆盖与地面坡度间表现为正相关关系，1988年至2018年[8˚~15˚)坡度范围内植被覆盖最高且增加了0.033，而[0˚~8˚)坡度范围内覆盖度最低且减少了0.004。本研究表明，抚仙湖流域内植被覆盖逐步提升，但空间变化上存在差异，这是人类活动的结果。

1. 引言

植被是陆地生态系统的重要组成要素，植被覆盖度能够真实客观地反映植物进行光合作用面积的大小，以及植被生长的茂盛程度，表征了植被的不同生长状态和生长态势 [1]。气候变化和人类活动是植被覆盖变化的重要驱动因素，自20世纪80年代以来，通过分析区域植被覆盖度演变进而研究气候变化和人类活动规律成为了重要的研究方法 [2] [3] [4]。植被指数是地表植被覆盖的有效度量指标，在目前已经定义的多种植被指数中，归一化植被指数(NDVI)因其客观性、精确性和连续性而被广泛使用 [5]。

抚仙湖面积216.6 km²、平均水深95.2 m、蓄水量206 × 108 m³，占据我国天然淡水湖泊I类淡水总量的91.4%，是珠江源区第一大湖，也是我国最大的深水型淡水湖泊和重要的战略备用水源地。由此，保护抚仙湖成为了国家政策制定和科学研究的焦点问题。近年来，在遥感和地理信息技术的支持下，分别从局部乃至整个流域的土地利用、植被动态等方面开展了多角度的研究工作。首先，从土地利用方面来看，发现抚仙湖流域主要土地利用类型为耕地、林地、水域和未利用土地，除水域外其余各类土地面积均发生了较大变化，未利用土地不断减少，耕地先减少后增加，其余5种地类呈现逐年增加的态势 [6] [7] [8]。人类活动对土地利用格局的形成影响强烈，人口增加、经济发展和城市化是研究区土地利用变化、退化的重要因素 [9] [10]。其次，从植被动态研究来看，抚仙湖流域自然植被集中分布的坡度区间主要在2˚~25˚之间，旱地集中分布在海拔2100 m以下，不同植被覆盖区面积与坡度间存在密切关系，年际间植被覆盖度有显著差异 [11] [12] [13]。然而，自上世纪70年代加大保护力度、实施整体保护以来，流域内植被覆盖的时空演变规律仍是亟需解决的难题。

自2010年开始，课题组开始关注本流域生态环境演变特征，系统获得了流域内自然地理及社会经济数据，为研究植被覆盖的时空演变奠定了了基础。鉴此，本研究以长序列卫星影像资料为基础，从时空变化角度分析抚仙湖流域植被覆盖动态规律，不仅有助于把握历史时期流域内的植被恢复效果，同时也将为《抚仙湖保护和开发利用总体规划(2018~2035)》的实施及流域内社会经济持续发展提供直接的科学支撑。

2. 研究区与数据

2.1. 研究区概况

抚仙湖流域位于云南高原中部(102˚39'~103˚06'E，24˚13'N~24˚47'N)，面积1084 km²，跨越了玉溪市的澄江、江川、华宁3县市，其中抚仙湖及其上游湖泊星云湖水面面积分别为216.6 km²、34.7 km²，陆地面积为832.7 km²。流域内气候属亚热带高原季风气候，年平均气温15.6℃~23.8℃，极端最高气温33.7℃，极端最低气温−3.9℃。年降雨量960~1250 mm，相对湿度76%，全年日照总时数2200 h，≥10℃的活动积温5200℃。土壤类型以红壤为主，地带性植被类型为亚热带常绿阔叶林，自然植被分布面积较大。

2.2. 数据来源与处理

本研究主要数据源为美国地质勘探局(https://glovis.usgs.gov/)网站获取的Landsat4/5 TM、Landsat7 ETM+和Landsat8 OLI的30 m分辨率卫星数字产品，以云量 < 5%为标准筛选出1988~2018年旱季的遥感影像序列。在ENVI/IDL 5.0软件支持下完成影像的图像增强、辐射定标和大气校正，并进行植被指数NDVI值计算。然后在ArcInfo Desktop10.3中裁剪出研究区范围，提取NDVI属性值并分析时空变化。

流域内澄江气象站年均气温、年降水量数据来源于国家气象信息中心(http://data.cma.cn)，使用前对数据进行了统一的质量检验，包括界限值检验、气候极值检验和内部一致性分析。

2.3. 研究方法

时间序列分析是时空演变分析中的重点，本研究中植被覆盖的时间演变趋势分析采用了Sen’s斜率法，该方法具有高效剔除干扰噪声的优良特性而被广泛采用。计算过程中，对于植被覆盖时间序列 $x_i=\left(x_1,x_2,\cdots ,x_n\right)$，Sen’s斜率的计算公式为 [14]：

$\beta =\text{Median}\left(\frac{x_j-x_i}{j-i}\right)$， $\forall j>i$ (1)

式(1)中： $\beta$ 为Sen’s斜率，正值代表分析序列呈上升趋势，反之为减少趋势， $\beta \ast 10$ 为植被覆盖数据每10a的变化量；Median为中值函数。

由于Sen’s斜率法无法完成植被覆盖序列变化趋势的显著性检验和突变分析，故变化趋势的显著性判断、突变年份确定采用了非参数Mann-Kendall秩次检验法 [15]。构造时间序列X，具n个样本量：

$S_k={\displaystyle {\sum }_{i=0}^k{r}_i}$， $r_i=\{\begin{array}{l}1,x_i>x_j,j=1,2,\cdots ,i\\ 0,\text{else}\end{array}$ (2)

秩序列S_k是第i时刻数值大于j时刻数值个数的累计数。时间序列随机独立的情况下，统计量UF_k为：

$U{F}_k=\frac{S_k-E\left(S_k\right)}{\sqrt{Var\left(S_k\right)}}$， $k=1,2,\cdots ,n$ (3)

式中，E(S_k)为均值、Var(S_k)为方差。

$X_1,X_2,\cdots ,X_n$ 相互独立，相同连续分布，计算公式如下：

$E\left(S_k\right)=\frac{n\left(n-1\right)}{4}$， $Var\left(S_k\right)=\frac{n\left(n-1\right)\left(2n+5\right)}{72}$ (4)

UF和UB是按照标准正态分布的算法，将时间序列X分布按照顺序和逆序排列之后得到的统计量。若UF值 > 0，呈上升趋势，反之呈下降趋势。当UF超过的临界值±1.96时，则通过了95%的置信度检验，呈显著上升或下降趋势，超过临界线的范围即为突变时间区域。若UF和UB在临界线之间出现交点，交点所在时刻即为突变时刻。该点即为突变点。

而植被指数的多时间尺度特征则采用了复值Morlet小波分析方法 [16]。该方法的基本数学思想是用一簇函数去表示或逼近一个信号或函数，它通过一个基本小波函数的平移和伸缩得到。对于连续函数 $f\left(x\right)$ $\left(f\left(x\right)\in L^2\left(R\right)\right)$，其连续小波变换公式为：

$W_{f\left(a,b\right)}=\langle f\left(x\right),{\psi }_{\left(a,b\right)}\left(x\right)\rangle ={\left|a\right|}^{-1/2}{\displaystyle {\int }_{R}f\left(x\right)\bar{\psi }}\left(\frac{t-b}{a}\right)\text{d}x$ (5)

式中， $W_{f\left(a,b\right)}$ 为连续小波变换系数， $\bar{\psi }\left(x\right)$ 为 $\psi \left(x\right)$ 的复共轭函数，a反映的是函数的尺度，为尺度因子， $a\in R$ 且 $a\ne 0$。它不仅影响窗口在频率轴上的位置，也影响窗口的形状，而b反映的是它沿时间轴的平移位置以及相位延迟，为时间平移因子，表示与频率相关的伸缩。

3. 结果与分析

3.1. 植被指数时间演变

自20世纪70年代以来，抚仙湖流域的保护治理逐步得到强化，尤其是“九五”以来，对抚仙湖的保护治理成为了国家和地方的重要任务。图1(a)为抚仙湖流域1988~2018年干季植被NDVI的变化趋势。从图中可以看出，植被指数在年际间波动较大，1998年达0.315而1995年仅为0.161，极差为0.154，变差系数Cv值为0.177。Sen’s斜率的 值为0.0012，表明流域内植被盖度呈现出上升的趋势，其速率为0.012/10a，这是强化流域整体保护的效果。

从抚仙湖流域植被NDVI的Mann-Kendall检验的统计值来看(图1(b))，研究时段内正向UF曲线和逆向UB曲线波动变化，两者间存在9个明显的交叉点，表明流域内植被覆盖可能存在9次突变过程。此外，从图1(b)中还可以看出，在2001年前两者虽然有交叉，但是UB曲线主要位于UF曲线之上，而2001年后则UB曲线主要位于UF之上，加之2001年的交叉点位于临界值−1.96~1.96之间，故可以确定2001年为流域内植被覆盖发生突变的年份。

图2(a)为抚仙湖流域植被NDVI序列的复值Morlet小波系数实部。从图中可以看出，小波系数实部在不同时间尺度上的震荡差异显著，在5年和20年时间尺度上周期性波动规律较为明显，这在小波系数模值(图2(b))中也可以得出同样的结论。在植被NDVI时序的小波方差图中(图略)存在2个较为明显的峰值，依次对应着5年和20年的时间尺度。其中，最大峰值对应着5年的时间尺度，表明5年左右的周期变化最强，为植被变化的第一主周期，20年时间尺度为第二主周期。

Figure 1. NDVI change trend of vegetation in Fuxian Lake Basin (a) and Mann Kendall test (b)

图1. 抚仙湖流域植被NDVI变化趋势(a)及Mann-Kendall检验(b)

Figure 2. Real part (a) and modulus (b) of Morlet wavelet coefficients of NDVI series

图2. 抚仙湖流域植被NDVI序列的Morlet小波系数实部(a)及模(b)

3.2. 植被覆盖空间变化

抚仙湖流域植被空间覆盖度差异较大，从图3(a)中可以看出，以抚仙湖、星云湖水面为中心，植被的分布具有环状分布特征，流域北侧梁王山、西侧谷堆山、南侧狮子山和东侧大黑山等山脉上植被覆盖度较高，而湖面周边，尤其是澄江市、江川市驻地，植被覆盖度在整个流域中最低。从1988~2018年流域植被覆盖的空间变化来看(图3(b))，植被覆盖增加区、减少区的面积为615.93 km²、218.85 km²，分别占流域陆地面积的73.78%、26.22%。此外，从图3(b)也可以看出，研究时段内，植被覆盖增加的区域为自山脚至山岭的过渡区，而减少的区域则主要发生在城镇建设区，即澄江、江川市政府驻地及澄江市的龙街镇一带。

对于植被覆盖与海拔的关系，从图4(a)中可以看出，两者间存在正相关关系，即随着海拔的升高，植被覆盖总体呈现上升的趋势。参照《抚仙湖流域禁止开发控制区规划(2006~2020)》，将抚仙湖流域环湖110 m范围内即从湖面至海拔1735 m划为核心保护区，1735~1770 m划为开发区，1770~1980 m和1980~2820 m分别划为控制开发区和禁止开发区。通过提取1988年~2018年各海拔范围内植被NDVI变化值后发现，核心保护区、开发区、控制开发区和禁止开发区的变化值分别为0.016、−0.019、0.085和0.035，即除了开发区外，其余陆地范围的植被盖度都呈现出增加的趋势，其中控制开发区范围内植被NDVI增加最为明显。

Figure 3. Spatial distribution (a) and change (b) of NDVI in Fuxian Lake Basin

图3. 抚仙湖流域植被NDVI空间分布(a)及变化(b)

Figure 4. Relationship between NDVI and elevation (a) and slope (b) in Fuxian Lake Basin

图4. 抚仙湖流域植被NDVI与海拔(a)及坡度(b)的关系

从植被覆盖与地面坡度的关系来看，在图4(b)中可以看出，整体上两者间同样表现为正相关关系，随着坡度的增加，植被覆盖呈现上升的趋势。为分析不同坡度植被覆盖的变化特征，将抚仙湖流域陆地区域划分为自[0˚~8˚)至[35˚~90˚)的5个等级，然后提取不同坡度植被覆盖均值及变化情况，结果见表1。

Table 1. NDVI characteristics of vegetation with different slopes in Fuxian lake basin from 1998 to 2018

表1. 抚仙湖流域1998~2018年不同坡度植被NDVI变化特征

从表1中可以清晰地看出，抚仙湖流域坡度在[0˚~8˚)的陆地面积达327.87 km²，占总陆地面积的39.37%，而坡度在[35˚~90˚)的陆地面积19.38 km²，仅占流域内陆地面积的2.33%。不同坡度陆域植被覆盖的差别不明显，无论是1988年还是2018年，均为[8˚~15˚)坡度范围内植被覆盖最高而[0˚~8˚)范围内覆盖度最低。统计各等级坡度内植被覆盖的变化后发现，自1988年至2018年，[8˚~15˚)坡度范围内植被覆盖增加了0.033，而[0˚~8˚)坡度范围内则减少了0.004。

4. 讨论

流域植被覆盖的时空变化影响因素众多，可以梳理为自然和人为影响方面。短时间尺度上自然因素特别是气温与降水变化是导致植被覆盖变化的重要因子，而长时间尺度上人类活动尤其是土地利用转型对植被覆盖的影响却更明显 [17]。

4.1. 气候因子与植被时间演变

已有研究表明，气候对植被覆盖的影响主要表现为对植被生长期年内韵律的控制方面 [18]。图5(a)为抚仙湖流域控制站澄江1988~2018年年均温变化图，由图可知，近30年来流域年均气温为16.8℃，波动于16.1℃~17.9℃之间，研究时段内气温变化趋势不明显。年均降雨量为944.5 mm，自1988年以来，降水量呈现出微弱减少的趋势，其中最大值为1994年达1362.2 mm，最小值为2011年仅为550.6 mm (图5(b))。流域植被覆盖与年均气温、年降水量的变化趋势存在明显差异，表明植被覆盖长期的增加趋势并非气候因素直接导致。

Figure 5. Variation characteristics of temperature (a) and precipitation (b) at Chengjiang station

图5. 抚仙湖流域澄江站气温(a)及降水(b)变化特征

此外，利用Morlet小波进行澄江站1988~2018年平均气温、年降水量数据变换，提取小波系数实部后可知：研究时段内年均气温在5、25年时间尺度上的周期性最显著，振荡具有全域性特征(图6(a))，而年降雨量在5、20年时间尺度上最强烈，具有全域性振荡特征(图6(b))。与同期流域植被覆盖的周期性进行对比后可以确定，气候因素尤其是年降水量的变化是植被覆盖产生周期性波动的原因。

Figure 6. Real part of Morlet wavelet coefficients of temperature (a) and precipitation (b) in Fuxian Lake Basin

图6. 抚仙湖流域气温(a)与降水(b)的Morlet小波系数实部

4.2. 人类活动与植被植被空间分布

人类活动影响日益剧烈的今天，植被覆盖变化深刻地反映了人类活动的特征。人类活动对植被的影响分为正面作用(如江河源石漠化区治理等生态工程)和负面影响(如城市扩展、森林破坏等)，其中退耕还林还草工程对区域生态环境的恢复和改善起到了巨大作用 [19]。对于抚仙湖流域，早在1993年云南省就出台了《云南省抚仙湖管理条例》，随后于2002年全面启动实施国家退耕还林工程，2007年《云南省抚仙湖保护条例》正式施行，抚仙湖也由依法管理步入了依法保护阶段。2012年，澄江县被列为国家石漠化综合治理试点县。这一系列水源地管理保护、石漠化综合治理工程的系统实施，极大促进了流域植被的恢复。至今，澄江县已累计完成人工造林31.3 km²、封山育林78.4 km²。正是这些工程的实施，使流域植被得以恢复，植被覆盖逐年增加。

5. 结论与展望

本研究探究了抚仙湖流域植被覆盖的时空演变规律及其影响因素。研究发现：

1) 1988~2018年抚仙湖流域植被指数以0.012/10a的速率上升，2001年为流域内植被覆盖发生了突变，植被指数存在5年和20年时间尺度上的周期性波动，其中5年为时间尺度变化的第一主周期。

2) 植被覆盖增加区、减少区的面积分别为615.93 km²、218.85 km²，增加区集中在山坡，而减少的区域则主要发生在城镇建设区。

3) 植被覆盖与海拔、坡度间存在正相关关系，随着海拔的升高植被覆盖总体呈现上升的趋势。

4) 植被覆盖周期性变化主要是降水变化的控制，而趋势变化则是人为活动的结果。

本研究成果将为《抚仙湖保护和开发利用总体规划(2018~2035)》的实施及珠江源区水源地管理提供直接的科学支撑。由于数据局限，本研究对植被覆盖时空变化的影响因素探讨局限在定性分析对比，下一步可以通过建立量化模型，定量揭示不同要素的贡献率。

基金项目

国家自然科学基金项目(41761105)资助Supported by the National Natural Science Foundation of China (41761105)。

参考文献

NOTES

^*通讯作者。