黄土高原典型区土壤保持服务效应研究

周日平

(中煤地质集团北京大地高科地质勘查有限公司，北京 100040)

0 引言

土壤侵蚀在陆地生态系统普遍存在，是全球最严重的环境问题之一[1]。减缓或控制土壤侵蚀一直是世界性的课题[2]。黄土高原地处于半湿润、半干旱和干旱区的过渡带[3]，是我国乃至全球水土流失最严重[4]、土壤侵蚀规律最复杂的地区之一[5]。为了遏制严重的水土流失，中国政府采取了调整土地利用结构、恢复植被、改进耕作方式、在坡面修建梯田以及在沟道修建淤地坝等一系列水土保持措施[6-8]，其中20世纪90年代以来，植被恢复措施和大范围的植树造林使得植被覆盖显著增加。截止2010年，造林面积约9.68万km2，使黄土高原土地覆被由水土流失较为严重的耕地，尤其是坡耕地等，陆续转变为水土保持功能强的草地和林地，区域土壤侵蚀状况得到了显著改善[9-11]，进而使得黄河下游的径流量和输沙量显著减少[12-14]。按照2010年12月国务院发布的全国主体功能区划，黄土高原中部地区是我国“两屏三带”生态安全建设的重要组成部分，定量分析黄土高原生态屏障区退耕还林还草背景下的土壤侵蚀时空动态以及对流域泥沙的影响，不仅有助于揭示黄土高原侵蚀环境的演变规律，也可为黄土高原生态建设和黄土泥沙的治理提供理论依据。

运用模型的方法开展定量测度是土壤侵蚀研究的常用手段[15]。20世纪60年代，美国学者Wischmeier等[16]最早建立了通用土壤流失方程(universal soil loss equation，USLE)。随着农业开发的需要，1993年美国农业部对USLE的局限性进行修正，研发了RUSLE(revised universal soil loss equation)模型[17]。Prasannakumar等[18]将RUSLE模型与地理信息系统(geographic information system，GIS)技术相结合对印度的Siruvani流域进行了土壤侵蚀风险评估； Asis等[19]将RUSLE模型与线性光谱混合分析方法相结合对菲律宾拉梅萨流域进行了土壤侵蚀评估； 卜兆宏等[20]基于RUSLE模型的结构，引入遥感数据开发出了与RUSLE模型相媲美的水土流失定量遥感方法； Liu等[21]以USLE/RUSLE为基础，通过研究坡面侵蚀量预报经验模型后，建立了适用于全中国土壤流失预报方程； 江忠善等[22]将沟间地与沟谷地区别对待，分别建立了与之对应的侵蚀模型。20世纪90年代相继出现了很多基于土壤侵蚀过程的物理模型，其中以美国的WEPP(water erosion prediction project)[23]模型最具代表性，分布式模型则以SHE(system hydrologic european)[24]最为典型，但由于土壤侵蚀过程的物理模型和分布式模型需要众多参数，它们的实用性均有所受限。而RUSLE模型从与土壤侵蚀密切相关因子出发(如降雨因子、植被覆盖和地形因子等)研究土壤侵蚀与这些因子的关系。该模型因操作简单，计算方便，加之随着RS和GIS的迅猛发展，是目前应用最广泛、具有较好实用性的土壤流失遥感定量模型[25-27]，并且在我国黄土高原地区应用广泛。

本文以黄土高原生态屏障区为研究对象，采用RUSLE模型，定量评估黄土高原典型区土壤侵蚀时空变化特征以及不同植被覆盖区和退耕还林还草区域土壤侵蚀状况，并且利用点面结合的方法，探讨区域生态建设的生态效益，以期对黄土高原山水林田湖草综合治理提供一定的科学依据。

1 研究区概况及数据来源

1.1 研究区概况

我国黄土高原生态屏障区地理位置为E105°1′～112°21′，N34°01′～38°13′，面积约11.67万km2，占黄土高原总面积的19.5%。行政区域共涉及山西省、陕西省、甘肃省和宁夏回族自治区等4个省区的60个区县，如图1所示。

图1 黄土高原生态屏障区研究范围

该地区地处黄河中游丘陵沟壑地带，境内有黄河穿流而过，支流交错，沟壑纵横，地形破碎，平均海拔为1 400 m，年均降雨量为400～776 mm，年内降雨分布不均，60%～70%的降雨量集中在6—9月份，且以暴雨为主； 土地类型以黄绵土为主，土质疏松，易于侵蚀，是我国水土流失的严重地区之一，也是退耕还林还草的主要地区。该区地带性植被类型是森林草原带，分布有黄龙山和子午岭等森林，是我国土壤保持重要功能区。

1.2 数据来源及预处理

采用的气象数据是来自中国国家气象信息中心(https: //data.cma.cn/)的2000—2010年全中国776个气象站点年、月平均降雨量，该数据用于计算降雨侵蚀力因子； 土壤类型数据为来自国家土壤信息服务平台(http: //www.soil.csdb.cn)的中国1∶ 100万土壤数据集，该数据用于计算土壤可蚀性因子； 地形数据是来自中国科学院数据云(http: //www.csdb.cn/n)的中国30 m SRTM-DEM数据产品，该数据用于提取坡长和坡度因子； 采用的遥感影像数据为2000—2010年间MODIS月最大合成数据产品，空间分辨率为250 m，本文对该数据进行了几何精纠正、辐射校正和大气校正等预处理，主要用于提取植被覆盖因子； 土地覆被数据为2000年、2005年和2010年3期30 m空间分辨率栅格数据，来自中国生态系统评估与生态安全格局数据库(http: //www.ecosystem.csdb.cn/ecogj/index.jsp)。采用的数据投影为Albers投影，地球椭球体采用Krasovsky 1940。为了进一步提高评估精度，本研究采用最高精度原则，将所有栅格数据重采样到30 m空间分辨率。

2 研究方法

RUSLE是目前应用广泛、具有较好实用性的土壤流失定量评价模型[28]，并且在我国黄土高原地区应用广泛。土壤侵蚀量是考虑任何植被覆盖和水土保持措施条件下的土壤侵蚀量，其计算公式为

A1=R·K·LS·C·P，

(1)

A2=1-R·K·LS·C·P，

(2)

式中：A1和A2分别为土壤侵蚀量和保持量，t · hm2· a-1；R为降雨侵蚀力因子，MJ·mm/(hm2·h·a)；K为土壤可蚀性因子，t·hm2·h/(MJ·hm·mm)；L和S分别为坡长和坡度因子，无量纲；C为植被覆盖因子，无量纲；P为水土保持措施因子，无量纲。

2.1 降雨侵蚀力因子

降雨侵蚀力因子采用Wischmeier等[16]提出的利用各月降雨量推求的经验公式计算。经验证明该模型在黄土高原地区具有较好的适宜性，其计算公式为

(3)

式中：p为年降雨量，mm；pi为月降雨量，mm；i为月份。

2.2 土壤可蚀性因子

RUSLE模型中土壤可蚀性因子是经过试验获得的定量数值，通常从标准小区单位降雨侵蚀力所造成的土壤流失量获得，因此计算的最好方法是通过天然小区直接测定，而在研究区大规模布设天然小区很难做到，故采用Williams等[29]在EPIC模型中土壤可蚀性因子的计算方法。该方法由蔡崇法等[30]在小流域进行应用研究，利用该公式计算的结果与实测值基本无差异，且方法更简便，其计算公式为

(4)

SNI=1-SAN/100，

(5)

式中：SAN，SIL和CLA分别为砂粒、粉粒和黏粒含量，%；Q为土壤有机碳含量，%。

2.3 坡长和坡度因子

在区域尺度上基于GIS提取坡长和坡度因子是应用RUSLE模型的难点问题之一[31]。本研究采用Van Remortel开发的模型RUSLE_PC.AML提取坡长，该模型坡度的提取采用栅格累积法[21]，即

(6)

(7)

式中：λ为坡长；m为坡长指数；θ为坡度值。

2.4 植被覆盖因子

植被覆盖因子的计算目前主要有3种方法： 人工赋值法、小区实验法和基于植被覆盖度的遥感数据定量估算法[32]。通过建立植被指数与植被覆盖因子的关系式可以更精确地估算植被覆盖度[33]。本文采用蔡崇法等[30]的植被覆盖因子计算公式，即

(8)

(9)

式中：f为植被覆盖度；NDVI为归一化差值植被指数；NDVImax和NDVImin分别为NDVI的最大值和最小值。

2.5 水土保持措施因子

由于黄土高原主要的管理因子是淤地坝和梯田，遥感数据的尺度分析很难反映出来[33]，故在计算中不予考虑，取值为1。

3 结果与分析

3.1 土地覆盖变化分析

3.1.1 土地利用变化分析

2000—2010年间黄土高原生态屏障区土地类型发生了显著变化。其中草地和林地面积增加幅度较大，面积占比分别由2000年23.56%和40.69%增加到2010年的26.46%和41.16%，其面积分别增加了3 342.19 km2和896.98 km2； 耕地呈现下降态势，面积由2000年的39 934.27 km2降低到2010年的35 575.34 km2，10 a间减少了4 358.93 km2，所占比例由2000年的34.21%降低到2010年的30.48%，降低了3.73%。黄土高原屏障区是重点退耕还林还草区域，按照退耕还林还草方式，可分为退耕还草、退耕还林以及草地转林地3种类型，经统计分析，在2000—2010年间，研究区退耕还林还草总面积为4 238.07 km2，占研究区总面积的3.63%，其中退耕还草面积为3 287.01 km2，退耕还林面积为211.22 km2和草地转林地面积为739.31 km2，分别占研究区总面积的2.82%，0.18%和0.63%。

3.1.2 植被覆盖度变化分析

2000—2010年间，黄土高原生态屏障区植被覆盖度呈上升趋势，由2000年的55.05%增加至2010年的69.26%，增速为1.29%/a，如图2所示。按照等值划分，将植被覆盖度划分为5个等级，即[0，20)，[20，40)，[40，60)，[60，80)和[80，100]，分别为低、中低、中、中高和高覆盖度。结果表明，2000—2010年间，植被在中低覆盖度区间的覆盖面积明显减少，低覆盖度和中低覆盖度面积所占比例分别由2000年的1%和28.79%降低到2010年的0.1%和1.8%，主要集中在退耕还林还草区域； 中等覆盖度比例基本保持平衡； 而中高覆盖度和高覆盖度比例均在增加，其面积比例分别由2000年的17.74%和19.67%增加到36.60%和28.52%。研究区植被覆盖整体在好转。

图2 2000—2010年植被覆盖度变化

3.2 生态屏障区土壤侵蚀特征

3.2.1 土壤侵蚀时间变化特征

2000—2010年间黄土高原生态屏障区土壤侵蚀模数减少趋势明显，2000年土壤侵蚀模数为6 579.55 t·km-2·a-1，2010年土壤侵蚀模数降低到1 986.66 t·km-2·a-1，降低了69.8%。但2003年，土壤侵蚀模数呈现较大幅度的增长，这与降雨密切相关。降雨是影响土壤侵蚀动态变化的主要因素之一[15]。20世纪90年代开始，黄土高原广泛开展了退耕还林还草工程，研究区植被覆盖度得到明显改善，尽管黄土高原降雨量呈现明显增长态势，年平均降雨量由2000年的455.41 mm增加到2010年的573.1 mm，降雨侵蚀力在增强，但由于良好的植被覆盖，土壤侵蚀出现明显降低，可见植被起到了良好的防治水土流失的作用。但在极端年份，土壤侵蚀仍然非常严重。如2003年研究区降雨量为704.14 mm，高出研究时段的平均降雨量182 mm，土壤侵蚀模数高达5 633.82 t·km-2·a-1。根据相关分析表明[27]，植被覆盖度的提高，可有效地阻止84%的水土流失，而降雨对土壤侵蚀的影响仅占21.8%。经统计，2000—2010年间黄土高原生态屏障区土壤侵蚀强度变化情况如表1所示。

表1 2000—2010年黄土高原生态屏障区土壤侵蚀强度面积百分比变化

从表1可以看出，2010年黄土高原生态屏障区微度土壤侵蚀占研究区总面积的72.54%，其次是中度、轻度、强度和极强侵蚀，分别占8.77%，6.62%，5.08%和4.64%； 而剧烈侵蚀仅占2.35%。2010年与2000年相比较，土壤侵蚀状况得到了明显改善，其中微度和轻度面积比例分别提高了43.44%和5.75%，但中度、强度、极强和剧烈侵蚀类型的面积分别降低了13.08%，35.29%，57.20%和83.67%，并且以剧烈侵蚀强度减小幅度最大。

3.2.2 土壤侵蚀类型变化特征

2000—2005年和2005—2010年间黄土高原生态屏障区土壤侵蚀强度转移矩阵分别如表2—3所示。

表2 2000—2005年黄土高原生态屏障区土壤侵蚀强度转移矩阵

表3 2005—2010年黄土高原生态屏障区土壤侵蚀强度转移矩阵

由表2可知，2000—2005年间，土壤侵蚀剧烈、极强、强度、中度和轻度向微度转化明显，其面积比例分别为32.33%，40.45%，47.01%，54.21%和60.84%，其次是各类型向中度转移幅度较大，其中剧烈、极强和强度分别转移了19.10%，19.12%和17.52%。但不可忽视的是，各较低侵蚀类型均有不同程度向更强一级侵蚀转移的现象，说明黄土高原生态屏障区水土流失治理任重道远。由表3可知，2005—2010年间，土壤侵蚀在中度和轻度土壤侵蚀等级上的转化比较明显。轻度和中度分别有47.31%和29.14%转化为微度侵蚀，强度侵蚀中有27.54%转化为中度侵蚀、极强侵蚀中有22.81%转化为强度侵蚀，剧烈侵蚀中31.90%转为极强侵蚀。研究区内土壤侵蚀是由剧烈侵蚀向微度侵蚀转变，侵蚀等级在逐渐降低。

3.3 植被梯度土壤侵蚀特征

分区统计2000—2010年间黄土高原生态屏障区不同植被覆盖度级别的土壤侵蚀强度变化百分比如图3所示。

图3 2000—2010年黄土高原生态屏障区不同植被覆盖度土壤侵蚀强度变化百分比

从图3中可以看出，在不同植被覆盖度下，微度侵蚀面积在逐渐扩大，其他各侵蚀类型面积在减少。剧烈侵蚀在低、中低、中、中高和高覆盖度的土壤侵蚀面积分别由2000年的77.93 km2，2 001.53 km2，7.06 km2和0.27 km2降低到2010年的3.19 km2，98.63 km2，335.47 km2和0.01 km2，且高覆盖度区域内不存在剧烈土壤侵蚀。2000年，在低、中低、中、中高和高覆盖度区域内，微度侵蚀的面积分别为21.52 km2，611.88 km2，2 402.24 km2，2 989.17 km2和3 384.42 km2； 2010年，其对应的微度侵蚀面积分别变化为3.50 km2，57.37 km2，2 009.74 km2，6 138.44 km2和5 306.62 km2。经过10 a的生态恢复，研究区内低覆盖度-剧烈土壤侵蚀面积在大幅度降低，而高覆盖度-微度土壤侵蚀类型面积在大幅度提升。

3.4 退耕还林还草区域土壤侵蚀特征

退耕还林还草区域土壤侵蚀强度面积统计结果如表4所示。

表4 2000—2010年黄土高原退耕还林还草区域土壤侵蚀强度面积统计

从表4可以看出，在2000—2010年间，黄土高原生态屏障区退耕还林还草区域土壤侵蚀基本呈现直线下降态势。在退耕还草区域，剧烈、极强和重度侵蚀面积均在下降，下降速率分别为12.14 km2/a，5.39 km2/a和1.75 km2/a； 而中度、轻度和微度侵蚀面积均在扩大，扩张速率分别为0.51 km2/a，1.85 km2/a和17.02 km2/a； 退耕还林区域和草地转林地区域各侵蚀面积变化基本相似，除微度侵蚀面积增加外，其余各类型面积均在下降。

3.5 区域土壤保持与观测站点泥沙关联分析

根据地理学第一定律及时空邻近度等概念[34-35]，区域土壤侵蚀与相邻站点泥沙观测数量应存在一定的关联性。本研究收集了距离研究区最近的华县和龙门站点泥沙数据，其中华县观测站点接收了泾河和渭河水量，土壤侵蚀主要集中在泾河和渭河流经的六盘山黄土丘陵沟壑地带； 而关中平原地形平坦，随着河流比降的下降，泥沙呈现大幅度的沉积，土壤侵蚀微弱，可忽略不计。龙门站点主要接收了来自窟野河、秃尾河、延河以及汾河流域水量，晋陕峡谷和白于山的黄土丘陵沟壑地带土壤侵蚀严重，是黄河泥沙(特别是粗泥沙)的主要来源，据统计，黄河干流56%的泥沙主要来源于此区间。同时本研究收集了各子流域2000—2010年间河流泥沙监测站点数据，并将该监测数据和对应的子流域土壤侵蚀量和保持量进行相关分析，发现二者之间存在较弱的相关性，主要是因为该区域地形平坦，土壤侵蚀微乎其微； 相反，将站点数据和研究区土壤侵蚀强度进行Persion相关分析，得出该区域站点监测数据和研究区土壤侵蚀强度具有较高的相关性。

分区统计研究区泾河—渭河流域和晋陕黄河峡谷流域土壤侵蚀量，并将其分别与华县和龙门站点的含沙量及输沙量进行相关分析，结果如图4—5所示。

(a) 泾河—渭河流域(b) 晋陕黄河峡谷流域

(a) 泾河—渭河流域 (b) 晋陕黄河峡谷流域

分析结果表明： 研究区泾河—渭河流域和晋陕黄河峡谷流域土壤侵蚀变化趋势基本一致。相比较，晋陕黄河峡谷流域土壤侵蚀量明显大于泾河—渭河流域。2000年华县和龙门站点含沙量分别为41.9 kg/m3和13.9 kg/m3，由于降雨量增加，2003年华县和龙门站点含沙量分别增加至89.6 kg/m3和21.4 kg/m3。2000—2010年华县和龙门站点含沙量随土壤侵蚀减弱而降低，10 a间分别下降了41.5%和73%。华县和龙门2个站点输沙量变化和含沙量变化规律基本相同。2000—2010年间华县和龙门站点输沙量分别降低了70.1%和64.47%。利用SPSS软件统计流域土壤侵蚀量与观测站点含沙量和输沙量之间的Person相关系数，结果表明，子流域土壤侵蚀强度与相关站点含沙量和输沙量均呈正相关，泾河—渭河流域和晋陕黄河峡谷流域土壤侵蚀量与站点含沙量的相关系数分别为0.34和0.60。泾河—渭河流域和晋陕黄河峡谷流域土壤侵蚀量与站点输沙量的相关系数分别为0.51和0.58。

4 结论与讨论

本研究通过应用修正的通用土壤流失方程，以黄土高原典型生态屏障区域为研究区，定量评估了2000—2010年间实施退耕还林还草生态工程后研究区的土壤侵蚀空间分异特征，得到以下结论：

1)尽管黄土高原降雨量呈现明显的增长态势，降雨侵蚀力在增强，但由于良好的植被覆盖，2000—2010年间黄土高原生态屏障区土壤侵蚀模数由2000年的6 579.55 t·km-2·a-1降低到2010年的1 986.66 t·km-2·a-1，降低了69.8%，并且低覆盖度-剧烈土壤侵蚀面积在大幅度降低，而高覆盖度-微度土壤侵蚀类型面积在大幅度提升。可见植被的恢复起到了良好的防治水土流失的作用。

2)土壤侵蚀状况发生了明显的改善，其中微度侵蚀面积提高了43.44%，而中度、强度、极强和剧烈侵蚀类型的面积比例在降低，其中剧烈侵蚀强度减小幅度最大，降低了83.67%。研究区内土壤侵蚀是由剧烈侵蚀向微度侵蚀转变，侵蚀等级在逐渐降低。

3)研究区以退耕还草为主，退耕还草面积为3 287.01 km2，占研究区总面积的2.82%。相比较，退耕还草区域土壤侵蚀变化较大，剧烈、极强和重度侵蚀以12.14 km2/a，5.39 km2/a和1.75 km2/a速率下降； 而中、轻度和微度侵蚀面积均在扩大。按照汇水面统计，晋陕黄河峡谷流域土壤侵蚀强度大于泾河—渭河流域，并且晋陕黄河峡谷流域与相关站点含沙量和输沙量Person相关系数均大于泾河—渭河流域。

由于本研究所用模型中地形因子是作为不变因子的，但实际上地形因子一直处在微小的变化中，而这种变化利用宏观遥感数据很难探测出来，故研究结果存在一定的偏差。在今后的研究中，需要进一步考虑模型的反演精度和在地形变化环境下的土壤侵蚀状况。