在气候变化和人类开发等多种因素的压力之下, 水土流失、植被退化和生物多样性丧失等生态问题逐渐凸显, 生态退化已经成为一个不容忽视的问题^[1]。在这种情况下, 生态恢复被认为是应对生态退化、改善环境质量的有效手段, 并得到广泛应用^[2]。在中国的生态恢复进程中, 黄土高原地区一直是恢复的热点和重点区域, 以退耕还林(草)工程为代表的大规模生态恢复项目则是恢复的主要手段。经过十余年的发展, 退耕还林(草)工程在植被恢复方面产生了显著的生态效益, 这已经在许多相关研究中得到体现^[3-5]。但是, 目前的研究多是基于一个特定时间段来分析植被变化的空间格局, 对于植被随时间的变化过程和特征关注相对较少。因此, 在黄土高原地区, 对植被随恢复年限延长而出现的时间及空间演化进行定量评估, 有助于深入了解植被变化机制。

另一方面, 生态系统服务, 即人类从生态系统中获得的利益^[6], 近年来已经成为学术界的研究热点^[7-8]。特别是2012年生物多样性和生态系统服务政府间科学政策平台(Intergovernmental Platform on Biodiversity and Ecosystem Services, IPBES)的成立进一步提高了学术界以及决策者对生态系统服务的关注。植被作为生态系统服务产生和传递的基础, 其恢复势必对区域生态系统服务产生重要影响。因此, 有必要在明确植被恢复态势的基础上, 进一步对区域主要生态系统服务变化进行定量模拟, 明确植被变化对生态系统服务影响的范围和趋势。

选取黄土高原中部的4个地级市, 即陕西的延安市、榆林市和山西的临汾市、吕梁市作为典型区, 对植被恢复范围及过程进行定量评估, 明确其时空演变特征及格局; 在此基础上, 对主要生态系统服务类型(土壤保持、水文调节、碳固定)变化进行模拟, 分析植被恢复对生态系统服务的影响, 以期为未来生态恢复工程的科学实施和可持续管理提供借鉴与参考。

1 研究区概况

由陕西省延安市、榆林市以及山西省吕梁市、临汾市组成的黄土高原中部典型研究区(35°25′~39°35′ N, 107°15′~112°34′ E), 总面积为121 391.7 km², 总人口约为1 376万人(2013年), 平均海拔约为1 190 m。该区域属于半暖温带向中温带的过渡区域, 大陆性气候特征明显, 气候干燥, 蒸发强烈, 四季分明; 降水时空分布不均, 呈自东南向西北递减趋势, 集中在7—9月。草原植被广泛分布, 主要物种包括茭蒿、铁杆蒿和冷蒿等; 林地在山地连片分布, 以次生林为主, 主要物种包括桦木、油松和辽东栎等^[9]。研究区生态环境较脆弱, 水土流失等环境问题较突出。为改善区域生态环境质量, 退耕还林(草)工程于1999年在黄土高原率先实施, 上述典型区属于项目试验示范和重点实施区域。

2 数据与方法 2.1 数据来源

归一化植被指数(normalized difference vegetation index, NDVI, I_NDV)数据是获得植被覆盖度(fractional vegetation cover, FVC, C_FV)数据的基础, 采用的是MODIS-NDVI产品, 空间分辨率为250 m, 时间分辨率为16 d。土地利用数据包括2000和2010年2期, 土地利用类型包括林地、草地、湿地、耕地、人工表面和其他共6种一级类型, 由中国科学院遥感与数字地球研究所提供。气象数据来自中国气象科学数据共享网(http://www.escience.gov.cn/metdata/page/index.html), 包括温度、降水和日照时间等。土壤属性数据由中国土壤特征数据集(V 1.1) 获得, 来自寒区旱区科学数据中心(http://westdc.westgis.ac.cn)。DEM数据来自地理空间数据云(http://www.gscloud.cn/), 分辨率为90 m。

2.2 分析方法 2.2.1 植被恢复趋势分析

(1) 评估指标

为了明确退耕还林(草)工程实施以来黄土高原典型区植被恢复的趋势, 采用C_FV为主要指标对典型区植被演变进行定量评估。C_FV是基于NDVI数据、采用像元二分法计算获取的^[10-11]。通过对各月16-day的C_FV取均值, 获取对应逐月C_FV, 再将年内各月C_FV均值作为年均C_FV, 从而获取2000—2014年黄土高原典型区年均C_FV数据序列。C_FV计算公式如下:

$ {C_{{\rm{FV}}}} = \left( {{I_{{\rm{NDV}}}}-{I_{{\rm{NDV, soil}}}}} \right)/\left( {{I_{{\rm{NDV, veg}}}}-{I_{{\rm{NDV, soil}}}}} \right)。$

(1)

式(1) 中, I_NDV为MODIS-NDVI数据产品经Savitzky-Golay滤波处理后的结果^[12]; I_NDV，soil和I_NDV，veg分别为纯裸土和纯植被对应的I_NDV值, 根据土地利用及野外实地监测数据^[13]确定。

(2) 趋势分析

线性回归方法是已经得到广泛认可及应用的趋势变化分析方法^[14-15]。为了确定典型区植被随恢复年限的时空演变趋势特征, 采用线性回归方法, 以2005和2010年为时间节点, 对2000—2005、2000—2010和2000—2014年C_FV变化情况进行逐像元定量分析。线性回归具体方法如下:

$ y = a \times t + b + \varepsilon 。$

(2)

式(2) 中, y为像元对应的C_FV; a为线性趋势系数; t为时间序列(值域为[2000, 2014]); b为截距; ε为随机误差。当a>0时, 表明像元表征区域的植被呈恢复态势; 反之则呈退化态势。恢复与退化的显著性水平设为α=0.05。同时也将该方法应用于典型区净初级生产力(net primary productivity, NPP)趋势分析。

2.2.2 生态系统服务评估

(1) 土壤保持服务

土壤保持服务采用土壤侵蚀和土壤保持量作为指标。土壤侵蚀数据基于通用土壤流失方程(universal soil loss equation, USLE)获得, 参考SL 190—2007《土壤侵蚀分类分级标准》, 根据侵蚀程度不同分为微度、轻度、中度、强烈、极强烈和剧烈6级。土壤保持量定义为没有植被覆盖和土壤管理措施情况下的土壤侵蚀量(潜在土壤侵蚀)减去当前土地利用/覆被格局下的土壤侵蚀量(实际土壤侵蚀)^[16], 即:

$ A = {A_{\rm{p}}}-{A_{\rm{v}}} = R \times K \times L \times S \times \left( {1-C \times P} \right), $

(3)

$ {A_{\rm{v}}} = R \times K \times L \times S \times C \times P。$

(4)

式(3)~(4) 中, A、A_p、A_v分别为土壤保持量、潜在土壤侵蚀和实际土壤侵蚀, t·hm^-2·a^-1; R为降水侵蚀因子, MJ·mm·hm^-2·h^-1·a^-1; K为土壤侵蚀因子, t·h·MJ^-1·mm^-1; L为坡长因子(无量纲); S为坡度因子(无量纲); C为植被覆盖因子(无量纲); P为管理因子(无量纲)。各项因子计算公式^[17-21]如下:

$ R = \sum\limits_{i = 1}^{12} {1.735 \times {{10}^{1.5 \times \lg }}\left( {P_i^2/{P_{\rm{a}}}} \right)-0.8188}, $

(5)

$ \begin{array}{l} K = \left[{0.2 + 0.3 \times {{\rm{e}}^{-0.0256 \times {S_{\rm{a}}} \times \left( {1-{S_{\rm{i}}}/100} \right)}}} \right] \times \\ {\left( {\frac{{{S_{\rm{i}}}}}{{{C_{\rm{n}}} + {S_{\rm{i}}}}}} \right)^{0.3}} \times \left[{1.0-\frac{{0.25 \times C}}{{C + {{\rm{e}}^{8.72 + 2.95 \times C}}}}} \right] \times \\ \left[{1.0-\frac{{0.7 \times \left( {1-{S_{\rm{a}}}/100} \right)}}{{1-{S_{\rm{a}}}/100 + {{\rm{e}}^{17.40 - {S_{\rm{a}}}/100}}}}} \right], \end{array} $

(6)

$ L = {\left( {\frac{\lambda }{{22.13}}} \right)^m}, \left\{ \begin{array}{l} m = 0.5, \;\;\;\theta \ge 9\\ m = 0.4, \;\;\;3 \le \theta < 9\\ m = 0.3, \;\;\;1 \le \theta < 3\\ m = 0.2, \;\;\;\theta < 1 \end{array} \right., $

(7)

$ S = {\left( {\sin \theta /0.0896} \right)^{0.6}}, $

(8)

$ C = \left\{ \begin{array}{l} 1, \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;f = 0\\ 0.6508-0.3436 \times \ln f, \;\;\;\;\;0 < f \le 78.3\% \\ 0, \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;f > 78.3\% \end{array} \right.。$

(9)

式(5)~(9) 中, P_a为年降水总量, mm; P_i为各月降水总量, mm; S_a为土壤砂粒质量含量, %; S_i为粉砂质量含量, %; C_n为黏粒质量含量, %; C为有机碳质量含量, %; λ为坡长, m; θ为坡度, %; f为地表植被覆盖度, %。

管理因子P在林地、草地、人工表面、水体和裸地均取值1, 在耕地林地交错带和耕地草地交错带取值0.8, P在耕地的计算方法^[22-24]如下:

$ P = 0.2 + 0.03 \times \theta 。$

(10)

式(10) 中, θ为坡度, %。

(2) 水文调节服务

地表植被的蒸散(evapotranspiration, ET)是植被参与水文过程、影响区域水平衡的重要环节^[25-26]。选择植被实际蒸散指标来表征典型区水文调节服务水平。定义植被蒸散为逐月地表实际蒸散量之和, 各月地表实际蒸散根据ZHANG等^[27]提出的经验模型求得, 如式(11) 所示:

$ {T_{\rm{E}}} = \frac{{Q \times \left( {1 + \omega \times \frac{{{T_{{\rm{PE}}}}}}{Q}} \right)}}{{1 + \omega \times \frac{{{T_{{\rm{PE}}}}}}{Q} + {{\left( {\frac{{{T_{{\rm{PE}}}}}}{Q}} \right)}^{-1}}}}。$

(11)

式(11) 中, T_E为地表实际蒸散, mm; Q为降水量, mm, 由典型区及周边地区气象站的降水监测结果经空间插值获取; ω为无量纲的植物水利用系数, 根据土地利用类型确定, 系数在林地取值为2.0, 灌木取值为1.0, 草地和耕地取值为0.5, 人工表面和裸地取值为0.1^[27-28]; T_PE为地表潜在蒸散量, mm。

T_PE根据HAMON^[29]提出的模型计算得到:

$ {T_{{\rm{PE}}}} = 0.1651 \times {d_{\rm{L}}} \times {T_{{\rm{RHOSA}}}}, $

(12)

$ {T_{{\rm{RHOSA}}}} = 216.7 \times {T_{{\rm{ESA}}}}/\left( {T + 273.3} \right), $

(13)

$ \begin{array}{l} {T_{{\rm{ESA}}}} = 6.108 \times \exp \left[{17.26939 \times } \right.\\ \;\;\;\;\;\;\;\;\;\;\left. {T/\left( {T + 237.3} \right)} \right]。\end{array} $

(14)

式(12)~(14) 中, d_L为逐月日照时数, 单位为12 h; T_RHOSA为月平均气温(t, ℃)条件下的饱和蒸汽密度, g·m^-3; T_ESA为特定温度条件下的饱和蒸汽压, mb^[30-31]。

(3) 碳固定服务

碳固定服务采用的指标为植被NPP^[32]。在土地利用、NDVI以及气象监测数据的支持下, 采用CASA(Carnegie-Ames-Stanford approach)模型来模拟植被吸收的光合有效辐射以及对光合有效辐射的转化比例, 从而实现大空间尺度上NPP的估算。CASA模型如下:

$ \begin{array}{l} {P_{{\rm{NP}}}}\left( {x, t} \right) = {R_{{\rm{APA}}}}\left( {x, t} \right) \times \varepsilon \left( {x, t} \right) = \left[{{L_{{\rm{SO}}}}\left( {x, t} \right) \times } \right.\\ \left. {0.5 \times {R_{{\rm{FPA}}}}\left( {x, t} \right)} \right] \times \left[{{T_\varepsilon }\left( {x, t} \right) \times {W_\varepsilon }\left( {x, t} \right) \times {\varepsilon _{\max }}} \right]。\end{array} $

(15)

式(15) 中, P_NP为像元x处在t时间的净初级生产力(以C计), g·m^-2·a^-1; R_APA为植被吸收的光合有效辐射, MJ·m^-2, 由太阳总辐射量(L_SO)、植被对光合有效辐射的吸收比例(R_FPA)估算获得; ε为光合有效辐射转化为有机碳的效率(以C计), g·MJ^-2, 由最大光能利用率(ε_max, 取0.389 g·MJ^-1)、温度胁迫(T_ε)和水分胁迫(W_ε)估算获得。年际NPP为各月NPP的累加之和。

(4) 植被恢复对生态系统服务的影响

Welch t检验适用于样本量不等、方差不等条件下的数据分析。采用Welch t检验来分析植被恢复对生态系统服务的影响^[33-34]。在黄土高原典型区范围内生成200个随机点, 作为分析样本。各样本点植被变化根据2.2.1节中趋势分析方法确定。对于生态系统服务变化, 计算研究区恢复前后生态系统服务的差值并将其标准化, 将标准化后结果在各缓冲区内的平均值定义为对应样本点的生态系统服务变化。具体检验方法详见文献[34]。

3 结果与分析 3.1 土地利用变化

基于2000和2010年土地利用数据, 结合相关研究成果^[4], 对典型区土地利用变化进行分析。结果显示, 退耕还林(草)工程实施前后典型区土地利用发生了较为明显的变化。工程实施前研究区土地利用以草地和耕地为主。2000—2010年典型区林地和草地面积较恢复前分别增加1.2%和9.8%, 耕地面积则减少16.3%。耕地向林地和草地的转化是耕地减少的主要原因, 两者合计占转出耕地总面积的94.2%, 表明退耕还林(草)工程实施成效显著。通过空间分析获取退耕还林区域的分布格局, 可以看出, 增加的林地主要分布于榆林市东部、延安市北部及中部, 增加的草地则主要位于吕梁市西南部和临汾市西北部(图 1)。至2010年, 林地、草地和耕地面积分别占典型区总面积的24.2%、44.0%和23.8%。

3.2 植被恢复趋势

利用ArcGIS 10.0软件对黄土高原典型区历年FVC数据进行回归分析, 可以定量获取2000—2005、2000—2010和2000—2014年典型区植被恢复分布格局, 确定植被演变随恢复年限增加而表现出的趋势特征(图 2)。

2000—2005年, 黄土高原典型区范围内, 93.6%的区域植被没有发生显著变化, 表明典型区植被覆盖度总体较为稳定。该时期植被显著退化(a < 0, P < 0.05) 区域比例为0.5%, 植被显著恢复(a>0, P < 0.05) 区域比例为5.8%, 区域总体处于植被净增加状态, 增加速率为1 070.1 km²·a^-1。显著恢复植被面积为7 079.9 km², 在分布格局上表现出明显的空间聚集特征, 主要分布于延安市东北部和西北部, 同时在临汾市西部也较为集中, 在典型区北部的榆林市和吕梁市则相对较少。另外, 值得注意的是, 该时期典型区范围内植被表现为不显著恢复的区域面积较大, 达到71 834.5 km², 占典型区总面积的59.2%。此类区域的植被恢复趋势尚不明显, 但是具有较大的恢复潜力。

2000—2010年, 黄土高原典型区范围内, 植被未出现显著变化的区域比例为50.0%, 较前一时段大幅降低, 植被覆盖度变化趋于剧烈。该时期植被显著退化的区域比例相对稳定(0.9%), 而植被显著恢复区域比例增至49.1%, 使得该区域植被净增加态势明显加强, 增加速率提高至5 317 km²·a^-1。在空间上, 与2000—2005年相比, 该时期植被显著恢复区域呈南北向扩展, 覆盖范围包括榆林市东部, 延安市中部及北部, 吕梁市西南部以及临汾市西北部地区。陕西省2市的植被恢复情况相对优于山西省2市。在2000—2005年表现为在植被不显著恢复区域有47 971.4 km²在2000—2010年转变为植被显著恢复区域, 转变比例达66.8%, 证实了植被不显著恢复区域所具有的恢复潜力。

2000—2014年, 黄土高原典型区范围内仅余20.4%的区域植被未发生显著变化, 主要为榆林市西部的毛乌素沙漠以及延安市南部黄龙、黄陵等地的大型林场。在退耕还林(草)工程实施15 a后, 典型区范围内79.0%的区域表现为植被显著恢复趋势, 植被净增加趋势进一步明显, 增加速率达到6 338.6 km²·a^-1, 其分布范围接近覆盖典型区全境。将2000—2005年的趋势与2000—2014年的趋势进行空间分析, 结果显示2000—2005年表现为在植被不显著恢复区域有91.3%在2000—2014年已转为植被显著恢复区域。

综合植被变化趋势和退耕还林区(即图 1所示增加林地和草地, 其余区域视为非退耕还林区)范围进行空间分析, 结果显示, 退耕还林区在2000—2010和2000—2014年表现出植被显著恢复比例为83.1%和96.2%, 明显高于非退耕还林区(47.5%和78.1%), 体现了退耕还林工程对植被恢复进程起到的显著促进作用。此外, 对比3个时段植被变化趋势可以发现, 植被覆盖度的增加速率并不均匀, 而是表现为随着恢复年限的延长而持续加速的过程。这一结果体现了植被对生态恢复响应的非线性特征, 同时表明植被恢复是一个长期过程, 生态恢复措施所产生的效果和效益需要随着时间的推移而逐渐显现。

3.3 生态系统服务变化 3.3.1 土壤保持服务

基于USLE模型对黄土高原典型区2000—2014年土壤侵蚀情况进行模拟分析(图 3~4)。结果(图 3)显示, 15 a间典型区土壤侵蚀以微度、轻度和中度侵蚀3个级别为主。在空间分布上, 各年土壤侵蚀分布格局相对一致, 微度和轻度侵蚀主要分布于典型区西部和南部, 而中度侵蚀地区则集中于榆林市东南部、延安市北部、吕梁市西南部以及临汾市西北部。2000年黄土高原典型区平均侵蚀速率为1 162.6 t·km^-2·a^-1, 土壤侵蚀以轻度侵蚀为主, 面积为63 735.6 km², 占典型区总面积的52.5%, 其次为微度和中度侵蚀, 面积分别为42 109.1和14 831.4 km², 所占比例分别为34.7%和12.2%。2014年典型区平均侵蚀速率为959.6 t·km^-2·a^-1, 中度侵蚀区域面积降低至6 866.4 km², 比2000年下降53.7%, 微度和轻度侵蚀区域面积则相应有所增加, 增至44 736.4和69 475.7 km², 占典型区的36.9%和57.2%。结合退耕还林区范围分析发现, 退耕还林区中有88.7%的区域土壤侵蚀有所减轻, 而在非退耕还林区比例为67.6%。此外, 还分析了2000—2014年黄土高原典型区土壤保持率(土壤保持量占潜在土壤侵蚀的比例)的变化趋势。结果(图 4)表明, 2000—2014年典型区土壤保持率维持在84%以上, 历年数据虽然有所波动, 但总体呈增加趋势, 特别是2010年之后土壤保持率均超过90%。上述结果表明, 2000年之后黄土高原典型区土壤侵蚀得到一定程度的治理, 土壤流失现象减轻, 区域土壤保持服务增强。

3.3.2 水文调节服务

采用地表植被的蒸散(ET)来表征生态系统水文调节服务水平, 对2000和2010年黄土高原典型区水文调节服务进行评估。结果(图 5)显示, 2000年典型区ET范围为135~426 mm, 2010年ET范围为170~404 mm。在空间上, 2000与2010年植被ET分布趋势基本一致, 以2010年为例, ET较高的地区主要分布于延安市南部、吕梁市东部、临汾市东侧边缘及西南部。土地利用数据显示, ET高值区土地利用类型以林地为主, 表明较其他用地类型而言林地的水文调节能力更强。综合2000和2010年的评估结果可以发现, 2000—2010年ET增加区域面积为48 094.1 km², 占全区面积的39.6%。ET增加区域覆盖榆林市大部分地区、延安市北部, 另外, 在吕梁市东南部也有少量分布, 但增幅较小。典型区内退耕还林区和非退耕还林区中ET增加区域所占比例分别为35.3%和42.2%。ET的变化表明随着典型区退耕还林(草)工程的开展和深入, 原有耕地和撂荒地等用地类型向林地和草地转化, 能够对区域地表植被蒸腾起到促进作用, 进而增强了水文调节服务。但是, 退耕还林区和非退耕还林区之间的差异表明这种增强作用可能需要更长的时间才能够体现。此外, 上述增强作用具有较显著的空间异质性, 在陕西省的效果比山西省更明显。

3.3.3 碳固定服务

NPP可以表征生态系统的年固碳量, 可用作典型区碳固定服务的主要指标^[31]。通过对2000—2014年黄土高原典型区NPP变化趋势进行分析, 结果显示, 2000—2014年黄土高原典型区NPP总体处于增加态势, NPP显著增加(P < 0.05) 区域面积为73 204.3 km², 占全区总面积的60.3%(图 6)。在空间上, NPP显著增加区域覆盖较广且连片分布, 其中, 在陕西省境内面积为52 120.4 km², 基本覆盖榆林市全境和延安市中部及北部; 在山西省境内面积为21 083.9 km², 较陕西省境内要小, 主要沿黄河南北向分布于吕梁市西部和临汾市西北部。但是, 同时NPP在延安市南部、吕梁市东部和临汾市南部等地有所下降。结合退耕还林区范围分析发现, 退耕还林区NPP增加区域比例达到98.9%, 非退耕还林区也达到87.3%。随着典型区NPP显著增加区域面积的扩大, 区域NPP总量(以C计)同样明显提高, 由2000年的39.2 Tg增加至2014年的57.0 Tg, 增幅达45.4%。

上述结果表明, 2000年以来黄土高原典型区生态系统固碳数量和规模不断增大, 区域碳固定服务明显提高。

综合上述生态系统服务变化评估结果可以看出, 2000年以来典型区各项生态系统服务均不同程度得到增强。为了进一步分析典型区植被恢复与生态系统服务变化之间的关系, 采用Welch t检验来分析植被得到显著恢复区域的生态系统服务变化(μ_A)是否高于植被未得到显著恢复的区域(μ_B)。检验方法如下:

$ {H_0}:{\mu _{\rm{A}}} = {\mu _{\rm{B}}};{H_1}:{\mu _{\rm{A}}} > {\mu _{\rm{B}}}, \alpha = 0.05。$

结果显示, 土壤保持服务、水文调节服务和碳固定服务的t值分别为4.66、12.54和2.98, 大于各自对应临界值(1.99、1.98和1.97), 拒绝原假设H₀。植被显著恢复区域的生态系统服务变化显著高于植被未显著恢复区域, 表明植被恢复对于生态系统服务具有促进作用。检验结果表明, 退耕还林(草)政策驱动下的植被大规模增加, 巩固并加强了生态系统服务产生的基础, 从而促进典型区生态系统服务的提高。此外, 土壤保持、水文调节和碳固定服务变化趋势间存在差异, 如在延安市西北部、吕梁市西部等区域, 碳固定服务表现出显著增强, 但水文调节服务则有所降低, 即不同生态系统服务间存在权衡关系^[35]。生态系统服务权衡的程度和范围等信息对于生态恢复及其可持续管理具有重要意义。

4 结论

通过基于遥感的植被变化趋势监测和基于模型的关键生态系统服务评估, 揭示了退耕还林(草)工程实施背景下, 黄土高原中部典型区的土地利用变化情况和植被恢复时空变化特征, 以及土壤保持、水文调节和植被固碳服务效应。结果表明, 随着恢复时段的延长, 植被恢复效果加速显现, 带动3类生态系统服务得到提升, 但是也可能引发3类服务之间的权衡问题, 比如ET增加会减少地表产水, 可能会影响可利用水资源量。所以, 大尺度植被恢复的时空过程及生态系统服务效应对于陆地生态系统恢复的规划、评价和可持续管理决策具有重要借鉴意义, 需要给予关注。