退耕还林(草)以来陕北白于山区植被覆盖与土壤侵蚀强度变化

黑 哲, 王 飞, 韩剑桥, 史尚渝, 戈文艳

(1.中国科学院 水利部 水土保持研究所, 陕西 杨凌 712100; 2.西北农林科技大学水土保持研究所, 陕西 杨凌 712100; 3.中国科学院大学, 北京 100049)

植被是陆地生态系统中的重要组成部分，在全球物质循环和能量流动中有重要作用[1]。研究地表植被盖度及其变化，对于揭示植被生长状况，空间分布规律，评价区域生态环境质量具有重要意义[2-3]。同时植被在防治土壤侵蚀过程中有重要作用[4]。一方面，植物冠层能够拦截降雨，植物枯枝落叶层能够增加入渗来减少径流[5]。另一方面，植被能够阻挡泥沙运移，减少水土流失[6-8]。在植被与土壤侵蚀的关系研究中，植被覆盖是影响土壤侵蚀的重要指标。传统的植被覆盖测量虽然在不断改进，但仍然有野外操作不便、成本较高、难以快速大范围动态监测等缺点，而遥感技术可以有效解决这一问题[9-11]。归一化植被指数(Normalized Difference Vegetation Index，NDVI)能够相当精确地反映植被盖度、生长状况、生物量以及光合作用强度等，是目前最常用的植被指数[12-14]。为了进一步认识植被与土壤侵蚀的关系，有必要研究植被的时空变化对土壤侵蚀风险性的影响。

土壤侵蚀风险性研究主要依据土壤侵蚀预报模型，国际上先后开发了通用土壤流失方程(Universal Soil Loss Equation,USLE)、修正的通用土壤流失方程(Revised Universal Soil Loss Equation,RUSLE)、水蚀预报模型(Water Erosion Prediction Project,WEPP)、欧洲水蚀预报模型(European Soil Erosion Model,EUROSEM)等。此外一些学者基于USLE建立了若干侵蚀预报模型[15]，但是以上土壤侵蚀预报模型都需要获取降雨侵蚀力、土壤可蚀性、坡度坡长、作物覆盖和水土保持措施因子基础数据，且需要多年长期监测数据，在基础设施薄弱的地区并不适用。有学者基于植被盖度、地形和土地利用进行水土流失风险评价，使得大范围长时间尺度的土壤侵蚀研究成为可能[16-18]。遥感技术能够对地形、植被覆盖、土地利用方式和地表裸露程度等土壤侵蚀的重要标志进行规则重复的观测，是进行动态环境监测的重要技术手段，目前被广泛应用于土壤侵蚀监测[19-20]。彭文甫等[21]利用Landsat5/8影像研究了岷江汶川—都江堰段植被覆盖动态变化；彭双云等[22]利用RUSLE模型对滇池流域1999—2014年土壤侵蚀时空演变进行了分析；陈学兄[23]利用坡度和植被盖度的指标对中国水土流失进行评价。

白于山区是陕西省水土流失最为严重的地区之一，该地区气候干旱，降雨集中，地形破碎，加之退耕实施不彻底，经济发展落后等因素的影响，导致区域水土流失现象并未明显缓解，土壤侵蚀风险一直较高。但目前对其土壤侵蚀的评估还较少涉及，因此快速评估白于山区植被状况和土壤侵蚀风险，对水土保持措施决策提供依据具有非常重要的意义。为此，本文利用MODIS-NDVI和Landsat数据，结合GIS手段，分析白于山区退耕还林(草)以来植被覆盖和土壤侵蚀强度变化特征，为该地区生态环境建设提供科技支持。

1 材料与方法

1.1 研究区概况

白于山区(107°16′17″—109°1′55″E，36°49′11″—37°32′55″N)位于陕北榆林市与延安市的交界处，北邻毛乌素沙地，南为黄土丘陵沟壑区。全区包括定边县、靖边县的南部，吴起县的大部，总面积约8 706.48 km2，最高海拔1 913 m，最低海拔1 178 m。白于山梁长沟深，沟坡陡峻，由主梁分出的次一级长梁分别向东、南、北方向延伸，呈波状阶梯式降低，为无定河、洛河、延河、清涧河等陕北主要河流的发源地。该区属于温带半干旱大陆性季风气候，年均气温为7.8～10℃，年降雨量320～490 mm，主要集中在5—9月。土壤类型多为细沙黄绵土，质地均一，土质疏松，抗侵蚀能力差，加之坡面流水侵蚀、重力侵蚀活跃，退耕还林实施面积小，坡耕地较多，水土流失严重。

1.2 数据来源及处理

本研究使用的遥感数据来源于美国国家航天航空局(National Aeronautics and Space Administration,NASA)EOS/MODIS卫星的MOD13Q1 v006数据产品数据集(https:∥lpdaac.usgs.gov/)，空间分辨率为250 m，时间分辨率为16 d，2000—2018年内共有434期影像。该产品属于MODIS三级数据产品，进行过大气校正和几何纠正。研究中使用MRT(MODIS Reprojection Tool)工具对数据进行波段提取、投影转换、格式变换等处理并用最大值合成法(maximum Value Compositess,MVC)获得月栅格图像和年栅格图像。为了更精确地研究植被盖度的地形分布特征和土壤侵蚀特征，本研究利用30 m分辨率的Landsat数据与相同分辨率的DEM数据和土地利用数据进行叠加。Landsat数据来源于地理空间数据云(http:∥www.gscloud.cn/)，空间分辨率为30 m，选取了2010年、2017年和2018年7月份无云覆盖的影像。利用ENVI 5.1进行辐射定标和大气校正并用像元二分模型计算植被盖度。数字高程模型(DEM)是空间分辨率为30 m的ASTER GDEM。土地利用数据来源于清华大学Finer Resolution Observation and Monitoring of Global Land Cover数据库(http:∥data.ess.tsinghua.edu.cn/)，空间分辨率为30 m。

1.3 研究方法

1.3.1 最大值(MVC)合成法 最大值合成法是将每个月的最大值假设为天气晴朗时获取的影像，通过最大值合成方法可以较好地去除云和大气的干扰。基本表达式如下：

NDVIi=(NDVIj)max

(1)

式中：NDVIi为第i月NDVI影像(i=1，2，…，12)，NDVIj为当月两景NDVI影像(j=1，2)。

1.3.2 NDVI趋势分析 为了对NDVI变化逐像元分析，采用一元线性回归分析法可计算每个栅格的NDVI值在研究时段中的变化趋势，并估计变化幅度。本文使用该方法模拟2000—2018年19 a平均的变化趋势，其计算公式如下：

(2)

表1 植被覆盖变化趋势分类标准

1.3.3 土壤侵蚀分级 研究中使用ArcGIS将研究区的土地利用类型、坡度和植被盖度数据叠加分析，并参照《土壤侵蚀分类分级标准》(SL190—2007)，将土壤侵蚀强度分为5个等级[24](表2)。

1.3.4 像元二分法 为了评估研究区土壤侵蚀强度，需要利用像元二分模型将植被指数转换为植被盖度(Fractional vegetation coverage，FVC)，再对植被指数进行等级划分。植被覆盖度计算公式如下：

(3)

式中：FVC为植被盖度值；NDVI为各像元的NDVI值；NDVIsoil，NDVIveg分别为纯土壤覆盖像元和纯植被覆盖像元的NDVI值。根据李苗苗等[25]提出的估算NDVIsoil和NDVIveg的方法，将研究区内NDVI累计频率为0.5%，99.5%时的值作为全土壤覆盖和全植被覆盖下的NDVI值，然后利用上式计算每像元植被覆盖度。

表2 土壤侵蚀分级指标

2 结果与分析

2.1 植被盖度空间变化特征

2.1.1 NDVI随高程变化特征 白于山区高程分布如图1所示，大部分区域高程范围在1 300～1 700 m，高程最高的区域位于定边县南部白湾子镇，从中心到四周高程环形递减。随着山区内高程的增加(图2)，植被盖度均呈现出“增加—减少”的趋势，不同年份均在1 250～1 400 m范围内达到峰值。

图1 白于山区高程分布

图2 2000年、2010年、2018年白于山区NDVI随高程的变化

2.1.2 NDVI随坡度变化特征 以2°为间距来估算研究区NDVI与坡度关系表明(图3)，2000年、2010年、2018年NDVI随坡度的变化趋势均为“增加—减少”，2000年植被盖度整体偏低，在坡度≤50°范围内NDVI呈现缓慢增长趋势，>50°陡坡处达到最高值0.34。2010年NDVI均值为2000年的1.38倍，2018年NDVI相比2010年增加了45%，其中在坡度0°～25°范围内快速增长，在坡度为56°时达到最高值0.65，在56°以上，各年份植被盖度均显著减少。

图3 2000年、2010年、2018年白于山区NDVI随坡度的变化

2.1.3 NDVI随坡向变化特征 将研究区坡向按照平地(-1°)、阴坡(0°—45°,315°—360°)、半阳坡(45°—135°)、阳坡(135°—225°)、半阴坡(225°—315°)分为5个坡向，按照年份分别统计研究区NDVI随坡向的变化特征研究，结果表明植被覆盖在坡向上的差异较小，2000年、2010年总体植被盖度随坡向变化趋势相同，从小到大依次为：阳坡<半阳坡<平地<半阴坡<阴坡，2018年平均植被盖度随坡向变化由小到大依次为：平地<阳坡<阴坡<半阳坡<半阴坡，表明坡面植被得以较快恢复。

2.1.4 NDVI空间分布特征 白于山区2000—2018年平均NDVI空间分异比较明显(图4)，其中植被覆盖最好的地区位于吴起县，NDVI最高值可达0.72，白于山区的NDVI多年平均值主要分布在0.2～0.3，占整个区域73.63%，整个研究区NDVI平均值仅有0.25，整体上植被覆盖偏低(表3 )。

图4 2000-2018年白于山区平均NDVI空间分布

2.2 植被盖度时间变化特征

2.2.1 NDVI年际变化特征 通过最大值合成法(MVC)处理得到2000—2018年研究区各像元年最大NDVI值并在整个研究区进行平均，以年平均值代表研究区植被覆盖最佳状况。图5是2000—2018年研究区植被覆盖年变化趋势，植被盖度在2000—2018年呈极显著上升趋势(p<0.001)，其增速为0.065/10 a。总体来说研究区植被盖度的年际变化大致可以分为3个阶段，第1阶段为2000—2003年，是大规模生态建设初期，植被覆盖随年份呈稳定增加趋势，其增速达到0.36/10 a(p<0.001)。第2阶段为2004—2015年，该阶段植被覆盖处于相对稳定状态，NDVI最低值(0.41)出现在2015年。第3阶段为2016—2018年，植被呈快速增长趋势，增速达到0.59/10 a(p<0.001)。

表3 2000-2018年白于山区平均NDVI空间分布

注：研究区平均NDVI值为0.248。

图5 2000-2018年白于山区NDVI年际变化趋势

为了进一步分析2000—2018年白于山区植被覆盖的年际空间变化趋势，经过Slope分析并通过Mann-Kendall检验计算出置信度水平。根据计算结果将趋势变化划分为严重退化到明显改善，将M-K检验计算结果划分为无显著变化、弱显著变化(p<0.1)、显著变化(p<0.05)和极显著变化(p<0.01),结果见图6—7和表4，空间格局的主要变化特征有：(1) 呈退化趋势的区域比重占3.18%，分布在定边县白于山南部丘陵沟壑区和靖边县中部的零星区域；其中定边县东部零星区域呈现极显著(p<0.01)退化趋势，主要是因为工程建设导致植被退化。(2) 基本不变区域比重为32.62%，连片分布在定边县中部地带即白于山南缘和靖边县中南部区域；(3) 轻微改善区域占研究区面积的45.66%，连片分布于吴起县东部、靖边县北部和南部区域、定边县北部区域。而且这些区域的改善趋势都通过了0.001的信度检验，为极显著改善；(4) 中度改善的区域比重为16.38%，主要是吴起县的东部区域、定边县白于山北部和靖边县的部分区域，改善趋势为极显著(p<0.01)；(5) 明显改善的区域较小，比重为2.16%，零星分布在吴起县的东部和靖边县的南部。

图6 白于山区多年平均NDVI时空演变

图7 白于山区NDVI变化趋势置信度水平

表4 2000-2015年白于山区植被覆盖演变趋势

2.2.2 NDVI年内变化特征 为了研究白于山区季节尺度上NDVI变化特征，采用气象学标准进行季节划分：春季3—5月，夏季6—8月，秋季9—11月，冬季12月—次年2月。图8为白于山区4个季节平均NDVI的年际变化，可以看到冬季的NDVI值最低，其次为春季、秋季，夏季最高，这反映了植被与物候的关系。统计得出，4个季节均呈现显著(p<0.05)上升趋势，但是增速有所不同，植被增速由快向慢依次为夏季(0.069/10 a)>秋季(0.044/10 a)>春季(0.039/10 a)>冬季(0.017/10 a)。

图8 白于山区季度NDVI变化趋势

2.2.3 NDVI月变化特征 由图9可知，多年NDVI的月变化均随年份呈现不同程度的增长趋势，其中1月、2月呈显著(p<0.05)上升趋势，其余10个月均呈极显著(p<0.001)上升趋势。增速最快的为7月(0.082/10 a)，最慢为2月(0.015/10 a)。研究区植被覆盖随月份变化显著，NDVI为0.15～0.46，其中2月最低，8月最高，平均值为0.25，其中6—9月都在平均值以上,植被盖度在年内的变化曲线整体呈单峰型波动。

2.3 土壤侵蚀变化特征

2.3.1 土壤侵蚀等级结构变化特征 根据水利部标准，利用植被盖度、坡度和耕地数据对研究区2010年和2017年土壤侵蚀状况进行评价，并统计出土壤侵蚀面积(图10，表5)。2010年和2017年土壤侵蚀等级分布面积最大的等级均为中度侵蚀，2010年土壤侵蚀面积为中度侵蚀>轻度侵蚀>强烈侵蚀>极强烈侵蚀>剧烈侵蚀，2017年土壤侵蚀面积为中度侵蚀>强烈侵蚀>轻度侵蚀>极强烈侵蚀>剧烈侵蚀。2010—2017年期间，轻度侵蚀面积显著减少，由2010年的3 437.87 km2减少至2017年的1 728.22 km2，轻度侵蚀面积占比减少了19.64%。中度侵蚀的面积变化较小，减少了5.23%。强烈侵蚀的面积变化较大，增加面积为15.24%。极强烈侵蚀和剧烈侵蚀面积均为小幅度增加，分别增加了7.64%，1.98%。

对比分析2010年和2018年研究区土壤侵蚀等级结构变化，计算得到两年土壤侵蚀强度转移矩阵(表6)。2010—2017年期间，轻度侵蚀面积减少19.64%，其中45.64%转化为中度侵蚀，5.68%转化为强烈侵蚀。中度侵蚀面积减少5.23%，其中1.43%的面积好转为轻度侵蚀，44.14%的面积转化为强烈侵蚀，11.80%的面积转化为剧烈侵蚀。强烈侵蚀的面积增加了15.24%，其中31.43%来源于极强烈侵蚀，7.92%来源于剧烈侵蚀，说明土壤侵蚀状况略有好转。极强烈侵蚀的面积增加7.64%，其中有20%来自剧烈侵蚀，11.04%来源于中度侵蚀。剧烈侵蚀的面积增加较少，仅有1.98%，其中4.69%来自于强烈侵蚀，15.56%来自于极强烈侵蚀，说明极强烈侵蚀的土壤在2010—2017年改善不明显。

图9 白于山区月均NDVI变化趋势

图10 2010年和2017年土壤侵蚀强度

表5 2010年、2017年白于山区土壤侵蚀面积

2.3.2 土壤侵蚀空间变化特征 为了探究2010—2017年土壤侵蚀等级变化的空间异质性，将白于山区2010年和2017年土壤侵蚀强度图叠加、相减，并分区得到土壤侵蚀等级变化分布图(图11)。其中4～-4代表侵蚀状况由显著改善到显著退化。对比2010年和2017年土壤侵蚀等级，其中，51.30%的面积土壤侵蚀有所改善，42.21%的面积土壤侵蚀无明显变化，6.49%的面积土壤侵蚀状况加剧，即研究区内大部分面积的土壤侵蚀有所改善。其中定边县土壤侵蚀加剧的面积占比为8.06%，稍大于吴起县的6.08%和靖边县的5.97%，整体退化趋势不显著。定边县和靖边县内侵蚀无变化的区域最大，面积占比分别为48.33%，41.47%，说明定边县和靖边县从2010—2017年土壤侵蚀状况未有明显好转，吴起县内土壤侵蚀轻微改善的面积最大，为41.08%。土壤侵蚀中度改善和明显改善的面积为吴起县>靖边县>定边县，即吴起县土壤侵蚀状况有较好改善。

3 讨 论

植被覆盖和土地利用方式是影响土壤侵蚀的重要指标，目前已经受到广大学者的关注[26-27]。本研究基于植被盖度、坡度和土地利用3个因子评价区域土壤侵蚀风险性，与传统的土壤流失方程相比，本文重点在于快速便捷地评价区域内土壤侵蚀风险性[28-29]。同时利用长时间序列的遥感数据对研究区平均NDVI时空变化进行了分析，关于植被与土壤侵蚀风险性的研究表明，随着植被盖度的增加，土壤侵蚀风险逐渐降低。这是由于植被具有较高的水土保持功能，其冠层能够防止雨滴直接击打地面，根系能够吸收降水，共同阻止土壤侵蚀的发生。如张雪峰等[30]对内蒙古锡林河流域草地生态系统的研究表明，提高植被覆盖可以显著提升草地生态系统的土壤保持能力。孙文义等[31]对黄土高原不同生态系统水土保持服务功能评价研究中表明植被抗土壤侵蚀能力随植被覆盖度的增加而增加，认为植被盖度大于70%时，地表土壤侵蚀很微弱。本研究与上述结论基本一致。除植被盖度之外，植被类型对土壤侵蚀的影响也有所不同[32]。例如实施“退耕还林(草)”等生态修复工程时植树造林的植被与农耕地分布有较大差异[33]。之后可以进一步分析白于山区不同植被类型对土壤侵蚀风险的影响，降雨侵蚀力与不同类型植被NDVI的一致性。

表6 2010年和2017年白于山区土壤侵蚀强度转移矩阵 km2

图11 白于山区土壤侵蚀等级变化分布

4 结 论

(1) 研究时段内，年均NDVI呈现显著上升趋势，增长速率为0.065/10 a。

(2) NDVI空间分布上，吴起县NDVI较高，而定边县西部和靖边县中部NDVI较低。随着海拔和坡度的增加,NDVI变化趋势均为先增大后减小，其中在海拔1 250 ～1 400 m达到峰值,在坡度50°～56°达到峰值；NDVI的坡向差异性较小，基本趋势为阴坡植被覆盖大于阳坡。

(3) 研究区植被退化面积为276.88 km2，主要集中在白于山南部丘陵沟壑区，植被保持稳定的面积为2 839.96 km2，主要集中在定边县中部，植被改善的面积为5 589.67 km2，主要集中于白于山北部和靖边县北部区域。

(4) 2010年和2017年二期土壤侵蚀状况表明，轻度侵蚀和中度侵蚀的面积减少，强烈侵蚀、极强烈侵蚀和剧烈侵蚀的面积增加。

(5) 2010年和2017年研究区土壤侵蚀加重的面积占比仅为6.49%，主要集中在定边县的北部和靖边县的中西部区域；土壤侵蚀改善的面积占比为51.30%，主要集中在吴起县。