自2002年4月国务院发布《关于进一步完善退耕还林政策措施的若干意见》后, 我国全面实施了退耕还林工程。据统计, 截止到2022年, 我国先后在25个省(自治区、直辖市)和新疆生产建设兵团实施两轮退耕还林还草工程, 完成退耕还林还草任务14.2万km², 同时完成配套荒山荒地造林和封山育林20.7万km², 共有4100万农户、1.58亿农牧民参与并受益^[1]。该工程在提升农业绿色全要素生产率、改善农村生态环境、增加农民收入等方面发挥了重要作用, 极大地促进了人与自然和谐共生, 加快推进了我国生态文明建设进程。2011年, 中共中央、国务院发布《中国农村扶贫开发纲要(2011—2020年)》, 明确将六盘山区等14个集中连片特困地区确立为扶贫攻坚主战场, 这些地区地势结构复杂, 水土流失和土地沙化严重, 生态环境极为脆弱, 曾是我国扶贫工作和生态保护的重点实施区域。2012—2020年期间, 党和国家相继在这些地区大力推行相关精准扶贫政策, 原集中连片特困地区的经济社会发展水平有显著提高^[2—4], 生态环境得到明显改善^[5—7]。而退耕还林还草作为其中一项重点生态减贫和生态修复工程, 对原集中连片特困地区的脱贫摘帽与经济社会发展产生了重要影响。因此, 探究原集中连片特困地区退耕还林还草工程的实施成效, 对新发展阶段进一步完善退耕还林还草政策、巩固和拓展退耕还林还草成果极具现实意义。

自退耕还林还草工程实施以来, 国内外许多学者聚焦于该生态修复工程的环境效益和经济效益展开了大量研究。有关退耕还林还草的环境效益, 多数研究是以某个特定区县或典型地貌为例, 对其实施退耕还林还草工程后的土壤蓄水能力^[8]、气候变化^[9]、土地利用结构^[10—12]、植被覆盖^[13—14]以及土壤有机碳^[15—20]等生态系统结构和功能进行探究, 并认为该项工程的实施有效消除了土壤侵蚀对生态系统所造成的负面影响^{[18, 21]}, 提升了植被的覆盖面积和总初级生产力^[22], 为保护生态系统服务功能和优化土地资源配置带来了较好的生态效益^[23—26]。而针对退耕还林还草的经济效益, 不少研究是通过对某个地区展开问卷调查的方式分析退耕农户收入变化来反映该效益^[27—32], 并认为退耕还林还草所产生的退耕补贴和劳动力转移是实现农户收入增长的关键动因^[33—37], 该项生态建设工程有效提高了农户的总收入水平^[38]。然而, 对于退耕还林还草工程的这种经济效益, 目前仍存在质疑, 有学者认为由于退耕还林还草工程无法从根本上改变农业生产结构, 因此难以从真正意义上实现农户增收^[39—43]。

总体而言, 当前有关退耕还林还草工程的研究已有较为丰富的成果, 尤其是针对该项工程所带来的生态效益。基于此, 本研究借鉴前人的研究成果, 考察集中连片特困区退耕还林还草的现状, 通过系统分析集中连片特困区的退耕区与未退耕区生态系统年总初级生产力(AGPP)和归一化植被指数(NDVI)的变化差异, 以期能够评估出退耕还林还草工程实施对集中连片特困区所带来的生态效应。为进一步完善集中连片特困区的退耕还林还草政策提供经验依据, 同时为该项生态修复工程的实施提供科学支撑。

1 研究设计 1.1 研究区概况

作为特定时期党和国家脱贫攻坚重点支持区域, 原集中连片特困区曾是我国深度贫困地区, 根据2007—2009年各县的人均GDP、农民人均纯收入等指标, 划定包括六盘山区、秦巴山区、武陵山区、乌蒙山区、滇贵黔石漠化山区、滇西边境山区、大兴安岭南麓山区、燕山-太行山区、吕梁山区、大别山区、罗霄山区等11个集中连片特困区, 再加上已明确实施特殊扶持政策的四省藏区、西藏以及新疆南疆三地州, 总共14个片区, 680个贫困县, 国土面积约为392万km², 涵盖680个县及县级单位, 其中有440个国家扶贫开发工作重点县, 252个革命老区县, 448个地质灾害高发县, 661个县属于地方病病区, 覆盖全国经济发展相对落后的县和贫困人口较为集中的地区。到2009年, 总人口2.36亿人, 其中乡村人口约为2.3亿人, 占全国的17.7%。集中连片特困区主要分布在我国中西部地区, 地区分布密度高、覆盖面积广, 不同连片特困区相互毗邻, 地势结构复杂、气候环境恶劣、地貌类型多种多样, 主要包括山地、丘陵、盆地以及草原等, 土地、耕地以及水资源普遍短缺, 生态环境脆弱；多数贫困县处于山区, 远离区域中心城市和城市群, 经济辐射带动作用也因基础设施建设落后而受限。由于存在少部分的贫困地区名称未显示在地图上, 本研究仅提取到了原集中连片特困区625个贫困县(图 1), 因此以这625个贫困县作为研究对象进行分析。

1.2 数据来源

研究中的DEM高程数据来源于地理空间数据云(https://www.gscloud.cn/), 其分辨率为1km×1km。矢量数据来自于国家基础地理信息中心1 ∶ 400万数据库(http://www.ngcc.cn/ngcc)。2000年和2020年两期土地利用数据来源于中国科学院地理科学与资源研究所(https://www.resdc.cn), 其数据分辨率为1km×1km。土地利用类型采用三级分类系统包括一级分为6类, 二级分为25类。根据研究需要, 采用ArcGIS10.7对一级分类中的耕地、林地、草地、湿地、未利用土地以及城乡建设用地等六种用地类型进行分析, 以观察原集中连片特困地区土地利用结构在退耕还林还草工程实施前后的变动情况。陆地生态系统年总初级生产力(AGPP)是粮食生产和陆地生态系统固碳的基础, 集中连片特困区2000—2020年AGPP数据来源于范雪等^[44]基于ChinaFLUX的长期联网观测数据和公开数据集, 利用随机森林回归树计算获得, 其分辨率为30arcsecond。归一化植被指数(NDVI)是反映植被覆盖情况的指标, 其数据来源于Landsat5/7/8(https://search.earthdata.nasa.gov/search), 分辨率为1km×1km。考虑当前针对生态系统功能的性状主要聚焦于AGPP和NDVI, 因而选取二者来评估集中连片特困区退耕还林还草的生态效应, 通过对比分析研究区内退耕区像元与未退耕区像元的AGPP和NDVI变化差异来反映该项生态修复工程所带来的生态效益。同时, 考虑虽然退耕还林还草工程和气候变化均会对AGPP和NDVI变化产生影响, 但退耕区与未退耕区在空间上相互交错, 彼此相邻, 自然因素对两个区块像元所产生的影响在统计特征上无显著差异, 因而可认为两个区块的AGPP和NDVI的变化差异主要是由退耕还林还草所导致^[21]。

1.3 研究方法

为探究集中连片特困区AGPP和NDVI的年际变化趋势, 采用Sen斜率估计^[45]对其进行计算, 该方法是一种稳健的非参数统计的趋势计算方法, 对测量误差和离群数据不敏感, 可以对长时间序列数据进行趋势分析。同时借助M-K统计检验法对AGPP和NDVI变化趋势的显著性进行检验, 该方法不需要测量数据服从正态分布, 不受缺失值和异常值的影响, 可以对长时间序列数据进行趋势显著性检验。Sen趋势分析的具体公式如下：

(1)

式中, β表示AGPP或NDVI的变化趋势, 当β>0时, 表示AGPP或NDVI存在增加趋势, 当β < 0时, 表示AGPP或NDVI存在减少趋势；i和j分别为时间序数；x_i和x_j分别为第i、j时间的AGPP或NDVI值。

进一步地, 结合AGPP和NDVI的变化趋势(线性回归的斜率)和二者多年的平均值, 计算AGPP和NDVI年际相对变化率, 从而得到每个像元AGPP和NDVI的相对自身平均水平的变化速率, 具体公式如下：

(2)

(3)

式中, δ_AGPP、δ_NDVI分别表示AGPP和NDVI的年际相对变化率, 和分别表示AGPP和NDVI在2000—2020年的多年均值。

2 结果分析 2.1 2000—2020年集中连片特困区土地利用变化

林地主要集中于西南部连片特困区, 如滇桂黔石漠化区、滇西边境山区、乌蒙山区、武陵山区、罗霄山区以及秦巴山区等连片特困区被大面积林地覆盖, 西藏区和四省藏区主要以草地为主, 新疆南疆三地州则以未利用土地为主(图 2)。同时, 还发现研究期内西藏区存在大面积土地由草地向未利用地转化。相比2000年, 集中连片特困区耕地、湿地以及建设用地面积分别减少了2.80%、5.61%、72.99%, 林地、草地以及未利用土地面积分别增加了30.54%、16.88%、27.27%, 退耕还林还草总面积约为178554km²(表 1)。退耕还林还草主要集中发生在滇西边境山区、滇桂黔石漠化区、乌蒙山区、秦巴山区、武陵山区、六盘山区以及燕山-太行山区, 大别山区、罗霄山区、新疆南疆三地州、四省藏区以及大兴安岭南麓山区均有少部分耕地向林地和草地转化, 且相对分散(图 3)。此外, 大兴安岭南麓山区、大别山区、六盘山区、秦巴山区的未退耕土地也相对集中, 尤其是大兴安岭南麓山区和大别山区最为突出, 其余地区的未退耕土地相对分散且分布范围较广(图 3)。总体来看, 集中连片特困区发生耕地转林地和草地的像元总数为178554, 占2000年整个区域耕地总面积的44.71%, 耕地未发生变化的像元总数为192679, 所占比例为48.24%, 这意味着退耕还林还草工程可能在集中连片特困区并未实现全覆盖, 还存在部分区县受地形、气候等自然因素的制约导致生态修复工程实施较为困难。

2.2 2000—2020年集中连片特困区AGPP和NDVI年均值变化及趋势变化分析

2000—2020年集中连片特困区AGPP和NDVI的多年平均值分别为376.87g/m²、0.4948, 且二者的空间分布格局呈现出“西南高、西北低”的特征(图 4)。AGPP和NDVI年均值的低值区间集中于新疆南疆三地州、西藏区以及四省藏区, 三者毗邻海拔高、地形复杂多样, 日温差较大, 以草地为主, 植被覆盖率不高, 因而AGPP和NDVI的数值相对较小。中高值区域则集中分布于西南部连片特困区, 包括滇西边境山区、滇桂黔石漠化区、武陵山区、乌蒙山区、秦巴山区以及六盘山区等, 这些地区降水充足、多气候带、动植物资源丰富, 因而具有较高的初级生产力。大兴安岭南麓山区、燕山-太行山区以及吕梁山区等东部偏北地区的AGPP多年的年均值处于中低值区域, AGPP相对较低的可能原因是受气温、日照等条件的影响, 由于偏北部地区, 气温相对较低, 日照时间相对较短, 导致植被的光合作用不强, NDVI年均值则处于中高值区域。

在研究期内, 集中连片特困区AGPP和NDVI的变化整体呈现出增加趋势, 其中AGPP和NDVI表现为极显著和显著增加趋势的像元个数总和占比分别为41.91%、48.03%(表 2), 且这些像元集中分布在滇西边境山区、滇桂黔石漠化区、武陵山区、乌蒙山区、吕梁山区、六盘山区、燕山-太行山区以及大兴安岭南麓山区(图 5)。AGPP和NDVI表现为微显著增加和不显著增加的像元个数总和占比分别为25.53%、36.47%, 这些像元主要分布在西藏区、四省藏区以及新疆南疆三地州, AGPP无变化的区域占比较大为22.06%, 集中于新疆南疆三地州和四省藏区, 而NDVI无变化的区域占比较小仅为2.26%。集中连片特困区AGPP和NDVI具有减少趋势的像元较少, 其中AGPP和NDVI表现出微显著减少和不显著减少的像元个数总和占比分别为9.06%、12.18%, 这些像元分布相对分散, 而表现为显著减少和极显著减少的像元个数总和占比分别仅为1.45%、1.06%。综合以上分析可知, 研究期内集中连片特困区的AGPP和NDVI均有明显提高, 生态环境有进一步改善的趋势。

2.3 集中连片特困区退耕还林还草对AGPP年际变化的影响

对比所有退耕像元与未退耕像元AGPP的年均值、变化趋势以及相对变化率, 可以发现退耕区与未退耕区的空间分布情况大体一致, 这说明两个分区受自然条件影响并无明显差别, 从而可以排除自然因素对AGPP变化影响的差异性。进一步具体分析, 研究期内退耕区像元AGPP年均值的平均值(1136.67g/m²)高于整个集中连片特困区(376.87g/m²)和未退耕区(988.26 g/m²)。这可能是由于在自然条件相近的前提下, 退耕区林地和草地面积扩张, 植被覆盖增加, 使生态系统的初级生产力水平提升。同时, 还可以看出退耕区和未退耕区大部分像元AGPP年均值处于中值区域(图 6), 六盘山区、吕梁山区、燕山-太行山区以及大兴安岭南麓山区的AGPP年均值则主要处于低值区域。

退耕区与未退耕区AGPP整体表现出显著增长趋势, 其中退耕区AGPP具有极显著和显著增加趋势的像元个数总和占比为65.30%, 未退耕区则为60.04%(表 3), 主要位于偏北部的连片特困区, 包括六盘山区、吕梁山区、燕山-太行山区以及大兴安岭南麓山区(图 7)；退耕区与未退耕区AGPP增加趋势不明显的像元个数总和占比分别为19.11%、26.24%, 主要位于西南部的连片特困区, 分散在秦巴山区、武陵山区、乌蒙山区、滇桂黔石漠化区以及滇西边境山区。退耕区与未退耕区无变化的像元个数极少, 占比分别为0.70%、0.74%。相比增加趋势, 退耕区与未退耕区所表现出的减少趋势较不明显, 其中退耕区与未退耕区AGPP减少趋势不显著和微显著的像元个数总和占比分别为7.40%、10.67%, 显著减少和极显著减少的像元个数总和占比分别为1.24%、2.30%。

对比研究期内退耕区与未退耕区AGPP的相对变化率, 可以发现, 退耕区和未退耕区像元AGPP的相对变化率均值分别为0.40%和0.39%, 两个分区大部分像元AGPP的相对变化率处于低值区域(图 8), 且主要位于西南部地区, 而处于中值区域的像元则集中在六盘山区、吕梁山区、燕山-太行山区以及大兴安岭南麓山区, 高值区域不明显。综合分析可知, 相比未退耕区, 退耕区AGPP的总体水平更高, 增长趋势更为明显, 相对变化率更快。这表明实施退耕还林还草工程为集中连片特困区的生态环境带来了较好的生态效益。

2.4 集中连片特困区退耕还林还草对NDVI的年际变化影响

对比分析所有退耕像元和未退耕像元NDVI的年均值、变化趋势以及相对变化率, 发现两个分区的空间分布情况也大致相同, 也即再次验证了自然因素对退耕区和未退耕区植被变化影响无显著区别。进一步具体分析, 整个研究区像元NDVI多年的平均值(8.8116)要远大于未退耕区(0.6913)和退耕区(0.7111), 同时退耕区NDVI年均值要高于未退耕区。两个分区NDVI年均值处于高值区域的像元主要集中于西南部和东北部连片特困区(图 9), 包括滇西边境山区、黔桂石漠化区、武陵山区、乌蒙山区、秦巴山区、大别山区以及大兴安岭南麓山区, 中值区域则分布在六盘山区、吕梁山区以及燕山-太行山区, 低值区域不明显。

研究期内退耕区和未退耕区的NDVI整体具有增长的态势, 其中退耕区与未退耕区NDVI增呈现出显著和极显著加趋势的像元个数总和占比分别为82.17%、67.93%(表 4), 且连片特困区中的大部分退耕区和未退耕区均表现出极显著增加趋势(图 10)。而增加趋势表现出微显著和不显著的像元个数总和占比分别为17.27%、22.81%, 其像元分布相对分散。无变化的像元个数占比极少, 均为0.01%。两个区块NDVI的总体减少趋势不明显, 其中退耕区与未退耕区NDVI表现出显著和极显著减少趋势的像元个数总和占比分别为0.69%、2.03%, 表现出微显著和不显著减少趋势的像元个数总和占比分别为4.20%、8.23%。可以看出退耕区和未退耕区的NDVI增长趋势明显, 特别是退耕区大部分像元表现出极显著的增长趋势, 而两个分区呈现显著减少趋势的像元均较少。

退耕区NDVI的相对变化率均值(4.76%)要高于未退耕区NDVI的相对变化率均值(4.20%), 但二者均小于整个集中连片特困区NDVI相对变化率均值(8.81%)。退耕区和未退耕区NDVI相对变化率集中处于中高值区域, 低值区域像元所占比例极小且分散(图 11)。这表明两个分区的整体NDVI变化具有较高增速, 尤其是退耕区的NDVI增长速度更快。总体来看, 相比未退耕区, 退耕区植被覆盖度和相对变化率更高, NDVI具有更明显的增长趋势, 这也进一步验证了退耕还林还草工程为集中连片特困区的生态环境改善作出了重要贡献。

2.5 集中连片特困区退耕还林还草对AGPP和NDVI年际变化的影响

对比2000—2020年所有退耕像元、未退耕像元以及整个研究区所有像元的AGPP和NDVI年际变化趋势, 可以发现退耕区多年AGPP年均值和NDVI年均值的中位数(分别为1164.41g/m²和0.7078)均要高于整个集中连片特困区(中位数分别为813.83g/m²和0.4935)和未退耕区(中位数分别为1017.01g/m²和0.6988)。研究期内退耕区AGPP和NDVI的年均值均集中分布在较大值区域；未退耕区AGPP的年均值集中在较大值区域, 而NDVI均值则集中在较小值区域；整个集中连片特困区的AGPP和NDVI年均值的分布比较均匀(图 12)。

进一步地, 从整个研究区、未退耕区以及退耕区AGPP和NDVI在2000—2020年的年均值变化趋势可以发现(图 13), 从整体来看, 无论是AGPP还是NDVI, 在区域尺度和像元尺度二者的时间序列变化趋势大体一致, 均表现出递增态势, 且退耕区的增长趋势始终保持在最上方, 整个研究区的AGPP和NDVI的年均值变化趋势相对小于退耕区和未退耕区, 这可能是由于部分连片特困区尤其是西藏区、新疆南疆三地州以及四省藏区, 占地面积广阔且以未利用土地为主, 因地理位置、气候以及土壤等因素条件的制约, 导致植被稀少, 植被覆盖面积小, 从而影响了整个研究区域AGPP和NDVI的数值大小。综合以上分析, 再次验证了退耕还林还草工程在集中连片特困区AGPP和NDVI的增长过程中发挥了重要作用。

3 结论与讨论 3.1 结论

(1) 2000—2020年集中连片特困区耕地、湿地、建设用地面积分别减少了2.80%、5.61%、72.99%, 林地、草地以及未利用土地面积分别增加了30.54%、16.88%、27.27%, 其中耕地转化为林地和草地的总面积为178554km², 约占2000年耕地总面积的44.71%, 耕地面积未发生变化的比例为48.24%。

(2) 2000—2020年集中连片特困区AGPP和NDVI多年平均值分别为376.87g/m²、0.4948, 二者的中高值区域集中于西南部和东北部的连片特困区, 而新疆南疆三地州、西藏区以及四省藏区的AGPP和NDVI年均值则主要处于低值区域。集中连片特困区AGPP和NDVI整体变化具有增长态势, 且二者表现出增加趋势的像元总和占比分别为68.44%、84.50%。

(3) 2000—2020年退耕区和未退耕区像元多年AGPP的平均值分别为1136.67g/m²、988.26g/m², 且两个区域AGPP呈现增加趋势的像元总和占比分别为65.30%、60.04%, 而呈现减少趋势的像元占比分别为8.98%、8.64%。此外, 退耕区与未退耕区像元的AGPP平均相对变化率分别为0.40%和0.39%。总体来看, 相比未退耕区, 退耕区AGPP的增长更为明显。

(4) 2000—2020年退耕区与未退耕区像元多年NDVI的平均值分别为0.7111、0.6913, 退耕区与未退耕区NDVI表现出增加趋势的像元总和占比分别为82.17%、67.93%, 而呈现减少趋势的像元占比分别为4.89%、10.26%。此外, 退耕区与未退耕区像元的NDVI平均相对变化率分别为4.76%、4.20%。总体而言, 相较于未退耕区, 退耕区的植被覆盖度有更为显著的提升。

(5) 相比于未退耕区, 退耕区的AGPP和NDVI年均值更集中于较大值区域。同时, 在研究期内整个集中连片特困区、未退耕区以及退耕区AGPP和NDVI的年际变化趋势大体一致, 均表现出递增态势, 且退耕区的AGPP和NDVI增长变化始终保持在最上方。

3.2 讨论

通过研究发现, 退耕还林还草工程对原集中连片特困区AGPP和NDVI产生了显著影响, 耕地转林地和草地后AGPP和NDVI均有明显提高, 这与以往对退耕还林还草工程的生态效益研究结果一致。如姚楠等^[46]和衣鹏慧等^[47]对黄土高原退耕还林还草工程的固碳效益进行评估后, 发现退耕后黄土高原的土壤有机碳储量有显著提升, 尤南山等则发现黄土高原的退耕区总初级生产力的增长趋势明显高于未退耕区^[22]。相比于以往以黄土高原作为研究对象, 本研究选择从原集中连片特困区这一生态环境极为脆弱区域入手, 考察退耕还林还草工程实施后相关连片特困区AGPP和NDVI的变化, 以评估退耕还林还草所带来的生态效应。可以发现, 不同连片特困区退耕后的AGPP和NDVI呈现出不同程度的增长态势, 其中AGPP和NDVI呈现极显著增长趋势的面积占退耕区总面积比例分别为55.72%和68.95%, 与以往发现黄土高原退耕区总初级生产力(GPP)和NDVI上升区域面积占全区域面积的67.3%和53.34%一致^{[13, 22]}。同时, 还可以发现以中部和偏北部的连片特困区退耕后的生态效益最为良好, 如六盘山区、吕梁山区、燕山-太行山区以及大兴安岭南麓山区等地区, 退耕后的AGPP和NDVI均表现出极显著的增长趋势, 存在极少部分相对分散的退耕区AGPP和NDVI增长趋势较不明显, 表明中部和偏北部的连片特困区对退耕还林还草工程的响应较为敏感, 这也说明针对不同连片特困区的生态环境状况, 需因地制宜, 遵循“宜草则草、宜林则林”的原则, 以便最大限度地发挥出退耕还林还草工程的生态效应, 从而进一步降低集中连片特困区的生态脆弱程度。

本研究基于2000和2020年两期土地利用数据将原集中连片特困区中的退耕像元和未退耕像元分离出来, 并计算出2000—2020年所有退耕像元、未退耕像元以及整个研究区像元的AGPP和NDVI年均值、变化趋势以及相对变化率, 试图从AGPP和NDVI两个视角考察原集中连片特困区实施退耕还林还草工程对生态环境所带来的生态效益, 但自然生态系统不仅受降水、气温等自然因素的影响, 同时还因退耕还林还草、植树造林、自然植被恢复等生态修复工程的实施而发生明显变化。因而植被的生长会因气候条件、人类活动等因素的影响, 导致AGPP和NDVI的提取值与真实值之间存在偏差, 但从已有研究来看, 由于退耕区和未退耕区相互毗邻, 气候变化对两个区域植被的影响无明显区别^[22], 因此在差异分析中并未讨论气候对AGPP和NDVI变化的干扰。在后续研究中, 应当综合考虑自然因素和退耕还林还草等生态建设工程的共同作用, 更加精准分离这些因素对生态系统服务功能的影响, 以更科学准确地评估出研究区退耕还林还草的生态效应。此外, 鉴于本研究仅分析了AGPP和NDVI的变化差异, 这无法全面反映退耕还林还草工程实施所带来的生态效应, 因此在之后的研究中也需进一步考虑SOC、NPP、NEP等生态指标的变化, 才能更加精准评价原集中连片特困区退耕还林还草的环境效益。