退耕还林（草）工程对黄土高原地表产沙的影响

管亚兵 ，赵长森 ，杨胜天 ，※，白 娟 ，彭睿文

（1. 北京师范大学水科学研究院，北京100875；2. 遥感科学国家重点实验室，北京师范大学地理科学学部，北京100875；3. 贵州师范大学地理与环境科学学院，贵阳550001）

0 引言

作为陆地生态系统的主要组分，植被是影响土壤侵蚀与水土流失的重要因子。植被具有显著的减水减沙效益，主要体现在：（1）植被冠层通过截留降雨减小降雨侵蚀力，同时增加地表糙度降低径流对地表冲刷能力；（2）植被根系可以改善土壤结构，增加土壤抗冲性能和入渗能力［1-2］。大量研究已经证明植被对流域水沙变化有深刻影响，表征植被长势的植被覆盖度指标与产沙系数之间存在较好的相关性［3-7］。基于以上认识，植被生长状况现已被广泛用来作为土壤侵蚀估算和风险评价的重要指标［7-10］。

黄土高原是世界上水土流失最严重的区域之一，约60%以上的土地已经遭受水土流失［11］，土壤侵蚀模数可达15 000 t/km2/yr。 土壤侵蚀可能导致严重的土地问题，如耕地蚕食、土壤有机质和养分流失［12］、土壤物理性质退化［13］。 进入河道的泥沙又将在下游不断淤积，抬高河床形成地上河，危害水坝和水库安全。因此，黄土高原水土流失不仅限制了当地经济社会的发展，而且威胁到黄河下游人民群众的生命财产安全［14］。

自1999 年退耕还林草政策实施以来，黄土高原植被覆盖率从1999 年的31.6%增加到2013 年的59.6%，增幅达88%，植被覆盖呈明显增加趋势［15］。同期黄河水沙锐减现象日益明显。黄河天然来沙量 16 亿t，天然径流量560 亿m3，但在2003—2013 年降雨略偏丰情况下，潼关实测年均来沙量只有2.6 亿t（其中河潼间偏少85%），花园口站天然径流量只有488 亿m3（其中河潼间偏少40%）［6］。在此背景下，很多学者研究了黄河流域林草植被变化特征［16］、林草植被变化对产流产沙的影响机制［6］、林草植被减沙方程构建［4-5］、林草植被减蚀临界效应［17］。然而，以往黄河流域植被变化与地表产沙关系研究多在试验小区或流域尺度上开展，以不同地貌分区为单元的研究还鲜有报道。该文基于多源卫星遥感数据估算了1980 年、2000 年、2010 年、2016 年黄河流域植被覆盖度，通过分析黄土高原重点研究区植被覆盖度的时空变化特征识别植被变化背景下不同地貌分区侵蚀产沙风险变化，试图回答对黄河减水减沙有重要贡献的潜在区域。

1 研究区概况

中国黄土高原是地球上最集中最大的黄土地区，总面积达 64 万km2。 由于土壤松散、地形破碎、植被稀疏、沟壑纵横，黄土高原是黄河流域主要的泥沙来源，也是世界上水土流失最为严重的和生态环境最脆弱的区域之一［18-19］。已有研究表明，黄土高原约 9.12 万 km2土地土壤侵蚀速率超过 8 000 t/km2/yr［20］。为改善黄土高原水土流失现状，自1999 年起，国家开始提出并实施退耕还林（草）、荒山造林等一系列大规模生态修复工程，黄土高原地区成为政府首批“退耕还林还草”工程试点区域之一［4，21］。

该文以黄土高原中部区域为重点研究区（图1），该区域水土流失严重，是黄河流域的主要产沙区，面积约为20.4 万km2，气候多为干旱半干旱类型。在地貌类型上，该重点研究区覆盖了平原、黄土塬、黄土破碎塬、片沙黄土丘陵、峁状黄土丘陵、梁状黄土丘陵、宽谷长梁丘陵、沙丘、风蚀沙化丘陵、石质山地、黄土低山、陇西丘陵12 种类型。

2 数据来源与方法

2.1 数据来源

研究所用的数据主要包括1 980 年前后MSS、2000 年前后TM、2010 年前后TM 和2016 年前后Landsat-8 OLI 多源遥感数据（表1）。这些数据均来源于美国地质勘探局网站（https：//earthexplorer.usgs.gov/），获取时间均为植被生长季，较好代表了一年中植被生长状况。此外，还使用到1∶100 万土壤类型图、四期（1980 年、2000 年、2010 年、2016 年）目视解译的土地利用图和黄土高原地貌分区图，分别辅助植被盖度信息提取和区域分析。

图1 黄土高原重点研究区域空间位置Fig.1 Location of the key study area on the Loess Plateau

表1 遥感数据成像信息Table 1 Information of the data sources

2.2 方法

景可［22］认为无论植被类型、降雨条件和其它下垫面条件如何，当植被覆盖度大于70%，地表侵蚀量都极其微弱。该文运用≥70%植被盖度面积变化指数揭示对地表侵蚀产沙有较强控制作用的林草植被变化情况，公式为：

式（1）中，VC≥70%表示≥70%植被盖度变化指数；2016VC≥70%表示2016 年植被盖度≥70%的面积百分比；1980VC≥70%表示1980 年植被盖度≥70%的面积百分比。VC≥70%越大，表示对地表侵蚀产沙有较强控制作用的林草植被面积增加越多。

对于植被盖度信息提取，由于归一化植被指数（NDVI）与植被长势、生物量、盖度和叶面积指数等有较强相关性，而且能部分消除辐照条件变化对反演参数的影响，该文采用式（2）所示像元二分模型反演植被盖度［23-24］。

式（2）中，VC为像元植被盖度；NDVI为像元归一化植被指数值；NDVIsoil为裸土或无植被覆盖区NDVI 值；NDVIveg为完全被植被覆盖区NDVI 值；NIR为近红外波段反射值；R为红光波段反射值。

为保证多时相遥感数据NDVI 有可比性，原始影像数据经辐射定标、简单大气校正和几何精校正后，生成的NDVI 用式3 所示相邻图幅递进回归分析法消除时相差异［25］。

式（3）中，Image′i为校正后影像数据矩阵，Imagei为原始影像数据矩阵，ki-1和qi-1分别为影像数据矩阵Image′i-1⌒Imagei的回归分析斜率与截距，n为影像图幅数量。

由于受大气效应、地表水分和植被类型等影响，NDVIsoil和NDVIveg会随时间和空间而改变。因此，使用土壤类型图分图斑统计NDVI 的累积频率，选取土种单元内累积频率为5%的NDVI 值作为NDVIsoil；使用土地利用图分类型统计NDVI 的累积频率，选取林地和草地累积频率为95%的NDVI 值作为NDVIveg。

3 结果与讨论

3.1 植被覆盖变化特征

3.1.1 植被盖度均值变化特征

使用像元二分模型计算得到1980 年、2000 年、2010 年和2016 年黄土高原重点研究区植被盖度。经多次野外考察GPS 定位拍摄典型样地照片，在植被盖度无显著变化区域内，用照相法植被覆盖度动态监测系统获取样地植被盖度。其中，“十二五”期间获取126 个样地植被盖度，2015—2017 年4 次野外考察获取78 个样地植被盖度。检验植被盖度遥感反演精度达90.18%，表明遥感反演植被盖度结果能满足研究需求。统计得到黄土高原重点研究区1980 年全区植被盖度均值为34.99%、2000 年为36.24%、2010 年为58.86%、2016 年为63.64%。表明黄土高原重点研究区植被盖度均值呈增加趋势，尤其是2000 年退耕还林（草）工程实施后增加的速度更为显著，且2010—2016 年间增速相对于2000—2010 年间较缓（图 2）。

3.1.2 植被盖度分级变化特征

为了集中反映黄土高原重点研究区植被变化特征，将研究区1980 年、2000 年、2010 年和2016 年植被盖度值从0～100%等距分为10 个等级，用GIS 空间分析方法分别统计4 个时期的植被盖度等级变化，如图3 和图4 所示。低等级植被盖度主要分布在研究区的北部和西北区域，高等级植被盖度主要分布在研究区东南部，呈现出由西北至东南方向的植被状况变好的空间趋势。

图2 黄土高原重点研究区全区植被盖度均值统计Fig.2 Statistical characteristics of the mean value of vegetation coverage in the whole region

图3 黄土高原重点研究区植被盖度空间分布Fig.3 Spatial distribution of vegetation coverage in key study area on Loess Plateau

统计分析研究区各等级植被盖度面积占比情况由图4 所示，1980 年植被盖度值主要集中在0～50%，共占全区面积比例为80.14%；2000 年植被盖度值主要集中在10%～60%，共占全区面积比例为86.47%；2010 年植被盖度值主要集中在30%～80%，共占全区面积比例为68.46%；2016 年植被盖度值主要集中在40%～90%，共占全区面积比例为74.56%。可以看出从1980—2016 年，峰值植被盖度等级由1980 年和2000 年的20%～30%过渡到2010年的50%～60%，再到2016 年的60%～70%，呈现出由低等级向高等级逐渐过渡的趋势；低值植被盖度区间0～20%占全区面积百分比呈现出持续减小的趋势，高值植被盖度区间50%～100%呈现出增加的趋势，较低值植被盖度区间20%～40%呈现出先增加后锐减的趋势。总的来说，黄土高原重点研究区植被盖度等级向着植被持续恢复的趋势发展。

图4 黄土高原重点研究区全区植被盖度分级统计Fig.4 Statistical characteristics of different vegetation coverage intervals in the whole region

3.2 产沙风险等级变化特征

该文根据野外实地调查确定产沙风险等级划分标准如表2 所示。 黄土高原重点研究区1980 年、2000 年、2010 年和2016 年4 期产沙风险等级如图5 所示。可以看出4 期产沙风险等级图均是产沙高风险和易产沙区主要分布在研究区的北部和西北部，不产沙区和产沙微弱区主要分布在研究区东南部。

表2 产沙风险等级划分标准Table 2 Criteria for classification of sediment yield risk

各产沙风险等级类型面积占比变化趋势如图6 所示。4 个时期均为易产沙区面积比例最高，1980 年和2000 年易产沙区面积比例均远高于产沙微弱区；2010 年和2016 年产沙微弱区面积比例与易产沙区面积比例差距逐渐缩小。1980—2016 年，产沙高风险区面积比例持续减小；易产沙区面积比例先小幅增加，后显著减小；产沙微弱区面积比例先小幅减小，后显著增加；不产沙区面积比例先小幅减小，后增大至稳定。表明由于退耕还林（草）工程实施，植被状况改善，研究区产沙高风险区面积比例显著减小。

图5 产沙风险等级分区图Fig.5 Spatial distribution of different sediment yield risk zones

图6 各产沙风险等级类型面积占比变化趋势Fig.6 The trend of the proportion of the area of each sediment yield risk class

3.3 植被盖度变化的地貌差异

3.3.1 植被盖度变化特征

图7 展示了不同地貌区的平均植被盖度变化情况。4 个时期均为黄土低山地貌区植被盖度均值最高，植被盖度均值均在70%以上；2000 年、2010 年和2016 年均为风蚀沙化丘陵植被盖度均值最低；2010 年和2016 年黄土塬和黄土破碎塬植被盖度均值均在70%以上；2016 年峁状黄土丘陵地貌区植被盖度也达到了70%以上；1980 年和2000 年除黄土低山外的所有地貌类型植被盖度均值均在70%以下；各地貌区均呈现2016 年植被盖度均值最高，2010 年植被盖度均值次之。表明自1999 年实施退耕还林（草）工程以来，不同地貌分区的植被状况均持续好转。

图7 不同地貌区平均植被盖度变化Fig.7 Variations of the mean value of vegetation coverage in different landforms

3.3.2 ≥70%植被盖度变化

图8 和图9 展示了不同产沙风险等级在不同地貌区的历年绝对面积变化情况和不同地貌分区≥70%植被盖度变化指数。对于不产沙区（VC＞90），黄土低山该类型面积比例较高，2010 年和2016 年面积占比均达到了70%；其余所有地貌区内该类型面积比例均较低，在30%以下。对于产沙微弱区（VC ≥70），黄土低山、黄土塬、黄土破碎塬、峁状黄土丘陵、梁状黄土丘陵、陇西丘陵该类型面积比例显著增加，至2016 年均达到了40%以上。对于易产沙区（VC＜70），除黄土低山外，其余地貌区1980 年和2000 年该类型面积比例均在80%以上，说明这两个时期研究区整体产沙风险很高；至2016 年，各地貌类型区该类型面积比例均呈现下降的趋势，其中黄土塬、黄土破碎塬、峁状黄土丘陵、陇西丘陵降幅比较大，均在40%以上。对于产沙高风险区（VC＜30），1980 年和2000 年平原区、片沙黄土丘陵、沙丘、风蚀沙化丘陵该类型占比比较高，约在60%以上；各地貌类型区植被状况均有所改善，至2016 年，该类型面积比例降幅比较大的有平原区、片沙黄土丘陵、宽谷长梁丘陵、沙丘和风蚀沙化丘陵，约在50%以上。

图9 可知，1980—2016 年，≥70%植被盖度变化指数最高的为黄土塬，其次为黄土破碎塬，最低的为沙丘。≥70%植被盖度变化指数小于30%的地貌区排序为：沙丘＜风蚀沙化丘陵＜平原＜黄土低山＜片沙黄土丘陵。由于黄土低山4 期植被盖度均值均在70%以上，因此沙丘、风蚀沙化丘陵、平原、片沙黄土丘陵地貌区是未来水土保持工作需要关注的区域，植被改善工作仍需加强。

图8 各地貌类型区不同产沙风险等级面积变化趋势Fig.8 Trends of variations of area of different sediment yield risk classes in each landform

图9 1980—2016 年不同地貌区≥70%植被盖度变化指数（VC≥70%）Fig.9 Variations of VC≥70% index in different landforms during the study period

4 结论与展望

基于多源遥感数据利用像元二分模型可实现大区域尺度长时间序列的植被覆盖动态变化监测。通过监测分析黄土高原重点研究区1980 年、2000 年、2010 年及2016 年的植被覆盖度状况，得到结论如下。

（1）黄土高原重点研究区全区植被盖度均值呈增加趋势，由1980 年的34.99%提升至2016 年的63.64%，且2000 年退耕还林（草）工程实施后增加的速度更为显著，2010—2016 年间增速相对于2000—2010 年间较缓。在空间上，植被状况呈现出由西北至东南方向逐渐变好的趋势。

（2）4 期产沙风险等级图均是产沙高风险和易产沙区主要分布在研究区的北部和西北部，不产沙区和产沙微弱区主要分布在研究区东南部。由于退耕还林（草）工程实施，植被状况改善，研究区产沙高风险区面积比例显著减小。

（3）自1999 年实施退耕还林（草）工程以来，不同地貌分区的植被状况均持续好转。沙丘、风蚀沙化丘陵、平原、片沙黄土丘陵地貌区是未来水土保持工作需要关注的区域，植被改善工作仍需加强。

相关研究表明黄土高原不同地貌类型区遏制侵蚀产沙的林草植被覆盖率存在不同阈值［17］。此外，王万忠基于水文站年输沙量数据，通过对比分析近50 年来黄土高原治理前、一期治理、二期治理和退耕后（ 2000—2009 年）4 个阶段侵蚀产沙的时间变化特征，表明受坡改梯、退耕还林（草）、淤地坝建设等多种水土保持措施的综合影响，侵蚀产沙量和减水减沙效益均表现出不同的区域分布特征［26］。该文采用了不同于上述研究的空间尺度，受研究观测资料的限制，进行了概化处理，主要包括：（1）在整个黄土高原重点研究区均按照景可等人提出的70%的植被覆盖度标准来划分黄土高原重点研究区侵蚀产沙风险等级；（2）仅考虑了植被覆盖度这一指标，分析了植被恢复背景下黄土高原重点研究区侵蚀产沙风险的时空变化。这些概化处理势必会对黄土高原重点研究区侵蚀产沙风险的区域时空分布特征造成一定影响，在未来的研究中需要构建考虑植被减沙阈值差异性、多水土保持措施综合影响下的区域侵蚀产沙风险评价方法，全面评估今后黄土高原治理的重点区域。