黄 麟,曹 巍,祝 萍
中国科学院地理科学与资源研究所中国科学院陆地表层格局与模拟重点实验室, 北京 100101
退耕还林还草工程是我国乃至世界上投资最大、政策性最强、涉及面最广、群众参与程度最高的一项重大生态工程[1- 3]。工程始于1999年,目标是将水土流失、沙化、盐碱化、石漠化严重及生态地位重要、粮食产量低而不稳的耕地,逐步停止耕种并因地制宜地造林种草、恢复植被[4]。据统计,截至2017年,退耕还林还草工程累计投入资金4500多亿元,完成退耕还林还草面积0.3亿hm2。国内外学者在退耕还林还草工程的生态、经济、社会效应方面开展了大量研究[5-7],通过分析土地利用变化、景观结构、植被覆盖状况、水土流失等指标变化评价生态效应[8-10],退耕农户家庭收入等分析经济效应[11],退耕农户适龄儿童入学、就医、家庭文化消费支出等分析社会效应[12-14]。评价方法从定性描述转向定量评估,评价内容从单一因素和单一目标评价向多因素、多功能、多指标转变。
国家层面形成了退耕还林工程建设效益监测评价国家标准[15]和国家报告[16],基于固定样地、森林生态站、工程效益监测点,以及全国退耕还林工程生态连清数据集,开展工程生态效益年度监测,测算了长江、黄河流域中上游及北方沙化地区不同植被恢复类型和不同林种类型的7项生态系统服务物质量与价值量。在局地尺度,国内外学者分别开展了退耕还林的工程模式、配置方式、植物种类、时间尺度对水土保持[17-20]、径流水文[21-22]、土壤质量[23-24]、抗蚀性[25-27]、碳储量与碳汇[28-30]、植被与动物多样性[31-33]、区域气候[34]等方面的作用评价。研究认为黄土高原、喀斯特山区、三峡库区、东北平原等区域退耕还林还草工程在增加植被覆盖、蓄水保土、增加基流、削弱侵蚀、减少泥沙输移等方面发挥了生态效益。然而,有研究认为退耕还林还草引起耕地大量减少从而对粮食安全存在不利影响,干旱半干旱区大规模人工造林导致严重的水资源短缺,导致黄土高原大规模土壤水分下降,影响水沙关系,部分区域出现覆盖度降低、植物种类多样性减少等负面作用[35-36]。
因此,准确、客观、科学地反映退耕还林还草工程的生态效应,对退耕还林还草工程的可持续发展具有极其重要的意义。对于如此大规模的工程,是否达到了预期的效果?退耕还林还草工程规划的“陡坡耕地基本退耕还林,严重沙化耕地基本得到治理”是否实现?实施退耕还林还草工程的县域生态状况是否得到较大改善?退耕还林还草工程是否对县域生态状况变化产生了正面效应?基于站点对比观测、统计数据分析、生态模型模拟等方法,可反映局地退耕还林还草的生态效应,但是由于缺乏工程的空间分布,难以确定工程区生态状况变化源于气候变化还是生态工程,或源于某个具体工程,难以在大时空尺度反映退耕还林还草工程的生态效应。针对这些问题,本文基于多源遥感获取的耕地转林地、耕地转草地时空分布信息确定退耕还林还草工程空间范围,根据工程规划目标选择退耕还林还草面积、植被覆盖度、土壤侵蚀量作为指标,分析退耕还林还草工程实施15年来目标实现情况以及工程的宏观生态效应。本研究将为退耕还林还草工程的巩固与滚动实施提供政策建议,为大尺度生态工程生态效应评估提供技术方法。
第一轮退耕还林还草工程建设范围覆盖25个省域的1897个县域,涉及国土面积达742万km2、3200多万农户、1.24亿农民。工程区划分为西南高山峡谷区、川渝鄂湘山地丘陵区、长江中下游低山丘陵区、云贵高原区、琼桂丘陵山地区、长江黄河源头高寒草原草甸区、新疆干旱荒漠区、黄土丘陵沟壑区、华北干旱半干旱区、东北山地及沙地区等10个类型区(图1)。工程规划目标:完成退耕地造林1467万hm2,宜林荒山荒地造林1733万hm2,陡坡耕地基本退耕还林,严重沙化耕地基本得到治理,工程区林草覆盖率增加4.5个百分点,工程治理地区的生态状况得到较大改善。根据林业统计年鉴,第一轮工程累计完成退耕地造林0.1亿hm2,退耕还草0.02亿hm2,实施荒山荒地造林和封山育林0.21亿hm2,退耕还林面积保存率达到98.9%。
图1 2000年中国耕地分布与退耕还林还草工程实施县域的分区分布Fig.1 The distribution of cropland in the year of 2000 and regionalization of the Grain for Green Project in China
基于中国科学院的中国土地利用与土地覆被变化(LUCC)数据库,提取2000年耕地分布与2000—2015年耕地转换林地、草地的动态变化数据。该数据集以陆地卫星TM/ETM+遥感图像为信息源,结合中巴资源卫星、环境小卫星图像数据,经图像精校正和拉伸处理后,通过人工解译获得土地利用类型及其变化空间分布,包括6个一级类型和25个二级类型[37]。利用野外调查资料进行统一质量检查[38],按10%县数比例开展精度验证,土地利用一级类型综合评价精度达到94.3%,二级类型达到91.2%[37]。2000—2015年耕地转换林地、草地的动态变化数据是对1 km大小的矢量栅格进行切割,得到每个栅格内耕地转为林地、林地转为草地的动态变化面积百分比及其比例。
收集2000—2015年全国范围1 km空间分辨率、16天时间分辨率的MODIS归一化植被指数(NDVI)数据(MOD13Q1),通过格式转换、重投影、拼接、重采样和S—G滤波处理,采用最大合成法(MVC)得到连续时间序列的半月NDVI数据。根据像元二分模型理论,认为一个像元的NDVI值是由绿色植被部分贡献的信息与无植被覆盖部分贡献的信息组合而成,利用NDVI根据如下公式计算半月尺度的植被覆盖度:
(1)
其中,VFC为植被覆盖度,NDVIveg是纯植被像元的NDVI值,NDVIsoil是完全无植被覆盖像元的NDVI值,依据1 km栅格百分比土地利用数据确定纯植被和无植被覆盖的像元。
采用最小二乘法分析植被覆盖度的年际变化趋势,计算公式为:
(2)
式中,i为2000年到2015年的年序号,i=1,2,3,…,n,某栅格像元的趋势线是这个像元的VFC值用一元线性回归模拟出来的一个总变化趋势,slo即这条趋势线的斜率,斜率为正,说明此像元植被覆盖度在该时间段的变化趋势增加,反之则减少。
本文中土壤侵蚀主要考虑水蚀和风蚀,土壤水蚀模数与风蚀模数的估算分别采用修正的通用水土流失方程(RUSLE)和修正的风蚀方程(RWEQ)。RUSLE表示为:
Mwater=R·K·L·S·C·P
(3)
式中,R为降雨侵蚀力,基于日降雨量资料的半月降雨侵蚀力模型[39]估算;K为土壤可蚀性因子,采用Nomo图法计算,公式如下:
K=[2.1×10-4(12-OM)M1.14+3.25(S-2)+2.5(P-3)]/100×0.1317
(4)
式中,OM为土壤有机质含量百分比(%),M为土壤颗粒级配参数,即粉粒、极细砂与粘粒百分比之积,S为土壤结构系数,P为渗透等级。土壤属性数据来源于寒区旱区科学数据中心(http://westdc.westgis.ac.cn/)的中国土壤特征数据集,该数据集来源于第二次土壤普查的1∶100万中国土壤图。
坡长因子(L)和坡度因子(S)基于McCool等[40]和Liu等[41]的方法。计算L时,把生态系统类型边界、道路、河流、沟塘湖泊等地表要素作为径流的阻隔因素,改进传统算法中通过相邻栅格间的坡向以及坡度变化率确定坡长终止点的方法,避免坡长因子的高估。
(5)
(6)
(7)
式中,λ为坡长(m),β取决于坡度百分比值θ,坡度s单位为弧度。DEM数据(SRTM3 V4.1)空间分辨率为90 m。
C值根据蔡崇法等[42]的方法利用植被覆盖度获得,
(8)
RWEQ表示为[43-44]:
(9)
Qmax=109.8(WF·EF·SCF·K′·COG)
(10)
式中,Mwind表示土壤风蚀模数,x表示地块长度,Qx表示x处的沙通量(kg/m),Qmax表示风力的最大输沙能力(kg/m),s表示关键地块长度(m)。WF表示气象因子,EF表示土壤可蚀性成分因子,SCF表示土壤结皮因子,K′表示土壤糙度因子,COG 表示植被因子。K′反映地表对风速减弱作用及对风沙活动的影响,即粗糙度反映了地表抗风蚀的能力,其大小取决于地表粗糙元的性质,通过文献参数整理得到草地、沙地、农田等类型的土壤糙度因子值。COG为枯萎植被、直立残茬和生长植被覆盖的土壤流失比率的乘积。
(11)
式中,WS2为2 m处风速(m/s),WSt为2 m处临界风速;N为风速观测次数(一般500次);Nd为试验天数;ρ为空气密度(kg/m3),g为重力加速度(m/s2);SW为无量纲的土壤湿度因子;SD为雪覆盖因子,即计算时段内积雪覆盖深度大于25.4 mm的概率。
(12)
(13)
式中,Sa为土壤砂粒含量,Si为土壤粉砂含量,Cl为土壤粘土含量,OM为有机质含量,CaCO3为碳酸钙含量。
退耕还林还草工程的资金和人力投入皆以县域为基本核算单元,因此本文对其生态效应的评估亦聚焦到县域及县域内耕地转林地或草地区域。首先,评估县域退耕还林还草工程与规划目标相比的面积完成率,进而通过退耕还林还草面积占县域耕地面积的比例、占县域15—25度以上耕地面积的比例,评价退耕还林还草工程实施的合理性。然后,基于工程规划目标,分析县域退耕还林还草区域的植被覆盖度、土壤侵蚀模数变化趋势,将此趋势与整个县域植被覆盖度、土壤侵蚀模数变化趋势相比,用于评价退耕还林还草工程提高植被覆盖度、增加土壤抗蚀效应在县域生态变化中的作用。第三,由于退耕还林还草一期工程涉及我国1897个县域,区域覆盖我国各种气候类型、地貌类型,加之各省域在具体实施过程中采用不同方式,因此,进一步分析退耕还林还草工程产生生态效应的区域差异。
据遥感数据估算(表1),2000—2015年,退耕还林还草工程县域的耕地转林地面积12.75万km2,其中2000—2010年11.27万km2,2010—2015年1.47万km2,面积增加较多的县域主要分布在东北山地及沙地区、川渝鄂湘山地丘陵区、云贵高原区、长江中下游低山丘陵区等;耕地转草地面积9.43万km2,其中2000—2010年8.75万km2,2010—2015年0.68万km2,面积增加较多的县域主要分布在华北干旱半干旱区、黄土丘陵沟壑区、云贵高原区等(图2)。
表1 2000-2015年不同区域耕地转林地或草地的面积统计
图2 2000—2015年退耕还林还草工程县域的耕地转林地、耕地转草地时空分布Fig.2 The spatial distribution of land cover from cropland to forest or grassland at the county scale from 2000 to 2015
2000—2015年,退耕还林还草工程区平均植被覆盖度年增加0.17%,增加较多的县域主要分布在黄土丘陵沟壑区、华北干旱半干旱区,特别是黄土丘陵沟壑区年增加0.56%(表2)。其中,耕地转林地区域的平均植被覆盖度年增加0.32%,增加较多的县域主要分布在黄土丘陵沟壑区、新疆干旱荒漠区、华北干旱半干旱区、琼桂丘陵山地区等;耕地转草地区域的平均植被覆盖度年增加0.43%,增加较多的县域主要分布在黄土丘陵沟壑区、新疆干旱荒漠区、川渝鄂湘山地丘陵区、东北山地及沙地区等。耕地转林地或草地区域与全县域植被覆盖度变化相比(图3),差值为正的区域说明退耕还林还草工程是该县域植被覆盖度增加的主要驱动力,特别是黄土丘陵沟壑区和东北山地及沙地区的县域;差值为负的区域说明退耕还林还草工程仅是该县域植被覆盖度增加的驱动力之一,比如西南高山峡谷区的县域。
表2 2000—2015年不同退耕还林工程区的植被覆盖度变化趋势
图3 2000—2015年退耕还林还草工程县域耕地转林地区、耕地转草地区、全县域的植被覆盖度变化,以及耕地转林/草区与全县域的植被覆盖度变化差值Fig.3 The variation trends of vegetation coverage in regions of cropland to forest, cropland to grassland, whole county, and the differences between regions of cropland to forest and grassland and the whole county at the county scale from 2000 to 2015
2000—2015年,退耕还林还草工程区平均土壤侵蚀模数减少较多的县域主要分布在黄土丘陵沟壑区、西南高山峡谷区和云贵高原区(图4)。其中,土壤水蚀模数年均减少0.13 t/hm2,特别是黄土丘陵沟壑区、川渝鄂湘山地丘陵区、长江黄河源头高寒草原草甸区(表3);土壤风蚀模数年均减少0.68 t/hm2,特别是新疆干旱荒漠区和华北干旱半干旱区。耕地转林地区的土壤水蚀模数年均减少0.43 t/hm2,土壤风蚀模数年均减少0.21 t/hm2,减少较多的县域主要分布在黄土丘陵沟壑区、川渝鄂湘山地丘陵区、西南高山峡谷区和云贵高原区;耕地转草地区的土壤水蚀模数年均减少0.55 t/hm2,土壤风蚀模数年均减少0.94 t/hm2,减少较多的县域主要分布在华北干旱半干旱区、黄土丘陵沟壑区、川渝鄂湘山地丘陵区和云贵高原区。耕地转林地/草地区与全县域土壤侵蚀模数变化相比,差值为正的区域说明退耕还林还草工程是该县域土壤侵蚀减弱的主要驱动力,特别是黄土丘陵沟壑区、云贵高原区的县域;差值为负的区域说明退耕还林还草工程仅是该县域土壤侵蚀减弱的驱动力之一,比如东北山地及沙地区、西南高山峡谷区等的县域。
图4 2000—2015年退耕还林还草工程县域耕地转林地区、耕地转草地区、全县域的土壤侵蚀模数变化,以及耕地转林/草区与全县域的土壤侵蚀模数变化差值Fig.4 The variation trends of soil erosion modulus in regions of cropland to forest, cropland to grassland, whole county, and the differences between regions of cropland to forest and grassland and the whole county at the county scale from 2000 to 2015
表3 2000—2015年不同退耕还林工程区的土壤侵蚀模数变化统计
基于工程规划目标的完成率与几个生态指标变化评价退耕还林还草工程的生态效应,可以看出:(1)与工程规划目标“完成退耕地造林0.15亿hm2,宜林荒山荒地造林0.17亿hm2”相比,遥感估算得到退耕还林(耕地转林地面积0.13亿hm2)的面积完成率为87%。(2)川渝鄂湘山地丘陵区、云贵高原区、黄土丘陵沟壑区、华北干旱半干旱区的县域退耕还林还草比例超过县域耕地面积的5%(图5)。根据遥感估算得到县域15度及以上坡耕地退耕比例,黄土丘陵沟壑区、川渝鄂湘山地丘陵区、长江中下游低山丘陵区、云贵高原区,特别是华北干旱半干旱区、东北山地及沙地区的农牧交错带、三江平原,基本达到 “陡坡耕地基本退耕还林”的规划目标。(3)遥感估算得到的工程区全县域平均植被覆盖度增加了2.6%、耕地转林地区植被覆盖度增加了4.8%、耕地转草地区植被覆盖度增加了6.5%,达到了“工程区林草覆盖率增加4.5个百分点”的规划目标。(4)植被覆盖度与土壤侵蚀变化反映“工程治理地区的生态状况得到较大改善”在部分区域显现成效,工程区年均土壤水蚀、风蚀模数分别减少0.13 t/hm2和0.68 t/hm2,其中耕地转林地区分别减少0.43 t/hm2和0.21 t/hm2,耕地转草地区分别减少0.55 t/hm2和0.94 t/hm2。耕地转林/草区与全县域的指标变化相比,可以看出退耕还林还草工程提高植被覆盖度与土壤抗蚀效应在县域生态恢复中的正面作用。
图5 2000—2015年退耕还林还草面积占县域耕地面积、15度及以上坡耕地退耕比例Fig.5 The proportions of cropland to forest and grassland accounted for the total cropland area in the county, and it accounted for the cropland area with slopes over 15 degrees
本文根据耕地转林地、耕地转草地时空分布信息确定退耕还林还草工程空间范围,针对工程规划目标,选择退耕还林还草面积、植被覆盖度、土壤侵蚀模数作为指标,评价了退耕还林还草工程实施15年来规划目标实现情况,工程宏观生态效应及其地域分异特征。得到主要结论如下:根据遥感估算结果统计,2000—2015年,退耕还林还草工程区耕地转林地面积12.75万km2,耕地转草地面积9.43万km2;工程区平均植被覆盖度年增加0.17%,其中耕地转林地区年增加0.32%,耕地转草地区年增加0.43%;工程区土壤水蚀和风蚀模数分别减少0.13 t/hm2和0.68 t/hm2,其中耕地转林地区分别减少0.43 t/hm2和0.21 t/hm2,耕地转草地区分别减少0.55 t/hm2和0.94 t/hm2。与工程规划目标相比,遥感估算得到退耕还林的面积完成率达到87%,工程区林草覆盖率增加4.8%—6.5%,大部分县域15度及以上坡耕地退耕比例超过50%。
退耕还林还草工程使得县域植被覆盖度得以逐年提高,生态环境持续改善。但是,适宜退耕还草、退耕还灌、退耕还林还是封育自然恢复,与工程区水分条件有直接关系,干旱半干旱区大规模人工造林可诱发严重的水资源短缺、产生对生物多样性的负面影响。因此,该工程在川渝鄂湘山地丘陵、长江中下游低山丘陵、云贵高原、琼桂丘陵山地等区域明显具有可持续性。然而,其生态效应在长江黄河源头高寒草原草甸、新疆干旱荒漠、黄土丘陵沟壑、华北干旱半干旱、东北山地及沙地等区域是否会缩水,值得进一步商榷。
本文结果的不确定性主要来源于数据精度和评估指标:首先,由于遥感数据本身时空分辨率的局限,基于多源遥感数据解译得到的耕地转林地、耕地转草地信息存在误差,比如初期造林难识别、农田作物与草地植被易混淆;其次,可以明确识别的退耕还林还草工程措施只有耕地转林地、草地,荒山造林和封育等是各类重大生态工程皆采用的普适措施,难以识别是否属于退耕还林还草工程,因此本文未涉及;第三,土壤侵蚀方程中的参数需要深入本地化,并考虑是否可以采用其他易行的替代方法,比如Nomo图法需要土壤结构系数和渗透级别资料等较多输入参数,而EPIC模型的K值估算只需要土壤有机碳和土壤颗粒组成即可;此外,本文仅利用面积完成率、植被覆盖度变化、土壤侵蚀量变化作为评价工程宏观生态效应的指标,可以回答工程主要规划目标是否实现的问题,但是未涉及工程在水源涵养、碳源汇、粮食供给等其他方面的生态效应,结果不够全面。下一步还需要深入分析退耕还林还草对县域核心生态系统服务权衡与协同的影响,以及退耕还林还草的可持续性问题。