基于RUSLE的引黄入晋北干线沿线地区土壤侵蚀定量研究

郭子萍, 王乃昂, 屈志勇

(1.兰州大学 资源环境学院/干旱区与沙漠研究中心, 甘肃 兰州 730000; 2.山西省气象信息中心, 山西 太原 030000)

土壤侵蚀是当今人类面临的一种最普遍、持续性最强的地质灾害[1-3]，会引起土地退化、土壤肥力下降、水环境恶化，从而影响生态安全和人类的可持续发展。黄土高原是中国水土流失最严重的地区[4]，而土壤侵蚀是黄土高原河流泥沙的主要来源[5]，因此，定量计算土壤侵蚀量、土壤侵蚀强度，评价水土流失时空分布特征，对该地区水土保持及土地的合理利用具有重要的意义。土壤侵蚀模型是土壤侵蚀量定量评价和预报的重要方法。国内外已经研究开发出许多实用的模型，其中修正的通用土壤流失方程(RUSLE)应用最为广泛[6-8]。

自20世纪90年代以来，在GIS和RS技术的支持下，中国学者利用RUSLE模型开展了广泛的土壤侵蚀研究，其中黄土高原作为研究的热点区域受到了广泛的关注，但相关研究大都集中于对侵蚀面积、侵蚀总量、空间分布的分析，对与影响因子的关系以及对人为因素作用日渐增强的地区关注度不够。引黄入晋工程是为促进地区经济发展，满足能源基地供水而建设的水利工程，对煤电基地建设、提高人民生活质量具有重要意义，但也不可避免的加重了水土流失问题。目前对于该工程的研究主要集中在水价分析[9-10]和社会经济影响评价研究[11]等社会经济领域；抗震性分析[12]、泥质膨胀岩工程地质研究[13]等工程地质领域；环境保护设计[14]等环境工程领域，针对土壤侵蚀的研究较少，几乎没有涉及。

本研究拟采用RUSLE模型计算引黄入晋北干线2005，2010，2015年3 a土壤侵蚀模数，分析土壤侵蚀时空变化特征，以期为当地水土保持和生态环境建设提供科学理论依据。

1 研究区概况与数据收集

1.1 研究区概况

山西省万家寨引黄工程位于山西省西北部，北干线主要任务是自万家寨水库引水，向山西北部的大同、平朔等能源基地供水，年引水量2.96×108m3，供水区总供水面积5 273 km2[14]。研究区属温带大陆性季风气候，年均气温5.7～7.9 ℃，年平均降水量380～440 mm,降水季节分配不均，降水变率大。研究区海拔为859～2 405 m，地形起伏变化较大[15]。该工程北干线2009年2月开始实施，2011年9月16日竣工。

1.2 数据来源

本研究所使用的数据包括： ①朔州市市辖区、山阴县、怀仁县、大同市市辖区2005，2010，2015年月降水量数据，来源为山西省气象信息中心。②研究区GDEMV230 M分辨率数字高程模型数据以及研究区2005，2010，2015年LandsatTM4-5遥感影像图，均来源于地理空间数据云。③研究区2005，2010，2015年土地利用图、研究区土壤类型图，来源于中国科学院资源环境科学数据中心。

2 研究方法

2.1 修正的通用土壤流失方程

修正的通用土壤流失方程(revised universal soil loss equation, RUSLE)作为估算年均土壤流失量的方法，需要考虑降水(月降水量)、土壤类型、坡度因子、坡长因子、植被覆盖与管理因子、水土保持措施管理因子6大主要因子，方程定义为[16]：

A=R·K·L·S·C·P

(1)

式中：A——土壤侵蚀模数〔t/(hm2·a)〕;R——降雨侵蚀力因子〔MJ·mm/(hm2·h·a)〕；K——土壤可蚀性因子〔t·hm2·h/(hm2·MJ·mm)〕；L——坡长因子(无量纲)；S——坡度因子(无量纲)；C——植被覆盖与管理因子(无量纲)；P——水土保持措施因子(无量纲)。

2.1.1 降雨侵蚀力因子R降雨侵蚀力因子是预测预报土壤侵蚀模型中的重要因子，本研究采用Wischmeier等[17]提出并经Arnoldus修改的比较适用于黄土高原地区且广泛应用的经验公式(2)来计算研究区2005，2010，2015年的R值。

(2)

式中：P——年降水量(mm)；Pi——第i月降水量(mm)[17]。

利用2005，2010，2015年北干线沿线完整且均匀分布的4个气象站点的气象数据，根据公式(1)计算得到研究区3 a的R值，利用ArcGIS工具连接功能将R值赋给各站点，并利用空间插值对R值进行空间化处理。研究采用普通克里格方法，半变异模型选用球形函数模型，进行研究区R值空间插值，生成研究区3 a的R值空间分布。

2.1.2 土壤可蚀性因子K土壤可蚀性因子K(soil erodibility)是土壤抵抗侵蚀能力的综合体现，影响因素有：土壤结构、土壤质地、土壤渗透率、紧实度、含水率以及黏土矿物性质等[18]，根据中国科学院资源环境科学数据中心公布的土壤类型数据，获得研究区土壤类型，以及《山西土壤》[19]中记载的关于各种土壤类型典型剖面的理化性质，采用Wischmeier等[7]于1990年在EPIC模型中建立的计算公式(3)计算研究区各土壤类型的K值。

(3)

式中：SAN——砂粒含量(%); SIL——粉砂含量(%); CLA——黏粒含量(%);C——有机碳含量(%)，SN1=1-SAN/100。

张科利等[20]的研究指出，与美国相比，中国土壤可蚀性K值普遍偏小，用EPIC模型所建立的公式计算的可蚀性K值，通过公式(4)进行修正转换，方可用于中国的土壤可蚀性估算。

K=0.013 83+0.515 75

(4)

利用公式(3),(4)计算的研究区K值，并在ArcGIS中生成K值分布。

2.1.3 坡度因子S通用土壤流失方程中允许计算的最大坡度为18%(10°)，研究区位于黄土高原，这样的条件不适用于研究区丘陵沟壑中的陡坡地形(>10°)，因此，在利用DEM提取坡度以后，要分段对坡度因子S进行计算，缓坡采用McCool坡度公式，陡坡采用刘宝元[21]提出的黄土高原S因子计算公式(5)。

(5)

式中：S——坡度因子；θ——坡度。

2.1.4 坡长因子L采用Hickey提出的非累计流量直接计算法估算坡长，坡长因子L的提取采用Wischmeier提出的坡长因子计算公式(6)：

L=(λ/22.13)m

m=β/(1-β)

(6)

β=(sinθ/0.089 6)/〔3.0(sinθ)0.8+0.56〕

式中：λ——水平坡长(m)；m——坡长指数； 22.13——标准径流小区坡长(m)；β——细沟侵蚀与面蚀的比率；θ——利用DEM提取的坡度。ArcGIS计算的坡度以度为单位，要转换为弧度表示，即乘以π/180。

2.1.5 植被覆盖与管理因子CNDVI是估算C因子最普遍应用的数据[22]。以经过辐射定标和大气校正的TM影像为数据，在ENVI 5.1中计算归一化植被指数NDVI，在提取NDVI之后，根据像元二分模型原理[23]，计算研究区的植被覆盖度F，公式为：

F=(NDVI-NDVIsoil)/(NDVIveg-NDVIsoil)

(7)

式中：F——有植被覆盖的面积比例，即像元的植被覆盖度； NDVIveg——纯植被覆盖像元的NDVI值； NDVIsoil——无植被覆盖像元的NDVI值。在求解植被覆盖度F后，参考有关研究成果，赋予了研究区不同土地利用类型和不同植被覆盖度下的C值(见表1)。

表1 研究区不同土地利用类型和不同植被盖度下的C值

2.1.6 水土保持措施因子Pp值介于0～1之间，0代表根本不会发生土壤侵蚀的地区，1代表未采取任何水土保持措施的地区。国内学者根据区域特点对不同的土地利用方式进行赋值，本研究区位于黄土高原，参照前人有关研究成果[24-25]赋值P，将水体、林地、草地、建筑用地、未利用地赋值为1，对于耕地，通常坡度越大，水土保持措施的作用越突出，因此耕地依据表2按照坡度范围进行赋值。

表2 不同坡度范围内耕地的p值

2.2 土壤侵蚀量估算

将上述各因子图层栅格单元统一为1 000 m×1 000 m,并标准化为相同的投影坐标系统WGS1 984_UTM_zone_49 N，地理坐标系GCS_WGS_1984,在ArcGIS软件的栅格计算器工具中将每个网格(Grid)中各因子连乘，得到研究区3期土壤侵蚀模数。RUSLE计算的土壤流失量默认单位是美制单位〔ton/(area·a)〕，要将结果乘以224.2得到中国常用的公制单位t/(km2·a)。

3 结果与分析

3.1 土壤侵蚀强度分级

在ArcGIS中根据水利部《土壤侵蚀分类分级标准(SL190-2007)》，对研究区土壤侵蚀模数进行重分类，得到3期土壤侵蚀强度图如附图2所示。

3.2 土壤侵蚀强度时空动态变化

由表3可知，从土壤侵蚀模数看，研究区2005，2015年的土壤侵蚀模数均在500～2 500 t/(km2·a)之间，属于轻度侵蚀，而2010年土壤侵蚀模数为3 340.41 t/(km2·a)，属于中度侵蚀。

表3 研究区土壤侵蚀变化统计

从2005—2010年，土壤侵蚀模数增长了80.91%；从2010—2015年下降了47.87%。3 a的土壤侵蚀模数呈先增后减，总体减小的趋势。但3 a中土壤侵蚀模数均大于水利部公布的黄土高原土壤允许流失量1 000 t/(km2·a)的标准。这说明研究区土壤侵蚀程度严重，应重点给予关注。引黄入晋北干线2009年2月实施，2011年9月竣工，2010年土壤侵蚀模数骤增与该工程的实施关系密切。工程建成后，为了区域的生态安全和经济的可持续发展，实施了一系列生物措施和工程措施，产生了良好的水土保持效果，因此2015年土壤侵蚀模数减小。

从图1中可以看出，2005,2015年均以微度侵蚀所占的面积最大，其次为轻度、中度、强度、极强度、剧烈侵蚀，土壤侵蚀面积随着侵蚀强度的增加而逐渐减小。与2005年相比，2010年由微度侵蚀转移到中度以上侵蚀的面积高达11.75%。中度以上侵蚀强度所占面积增大，无疑会加重土壤侵蚀。2015与2010年相比，微度、轻度侵蚀面积增加，中度以上侵蚀所占面积减小12.63%，侵蚀状况明显好转，且各侵蚀强度所占的面积与2005年基本相同。

图1 引黄入晋北干线沿线地区不同年份各侵蚀强度所占面积比

由表4可知，从空间上看，研究区土壤侵蚀呈块状侵蚀和点状侵蚀相间分布的特点。3 a中土壤侵蚀面积均为朔州市市辖区最大，块状侵蚀分布较多，其次是大同市辖区、山阴、怀仁，且3 a中土壤侵蚀量也均为朔州市辖区最大。从土壤侵蚀模数看，4个地区土壤侵蚀模数从2005—2015年均呈先增后减的趋势，在2010年处达到峰值，因此每个地区的土壤侵蚀模数与总的土壤侵蚀模数分布特点一致。

由附图2可以看出，以北干线为界，线上方土壤侵蚀明显高于线下方，两侧侵蚀状况差异显著，这种趋势在2010年尤为明显。土壤侵蚀的空间分布与研究区的DEM也有较高的吻合度，分析研究区DEM可以明显看出，研究区分为2个地貌单元，分别为侵蚀中低山区和冲洪积平原区[15]，侵蚀中低山区土壤侵蚀明显比冲洪积平原侵蚀严重。

3.3 土壤侵蚀主控因子分析

3.3.1 土壤侵蚀与降雨 从图2可以看出，研究区不同地域单位R因子均在2010年处达到最大值，而折线图显示，土壤侵蚀模数也在2010年处达到峰值。在其他因素不变的情况下，降雨侵蚀力因子R是土壤侵蚀的加速因子。因此，当R因子的值增大时，土壤侵蚀模数也会相应的增大。分析所获得的降水数据，研究区2005年年降水量为367.73 mm，2010年为423.5 mm，2015年为380 mm，由此可以看出2010年降水量显著增大是导致R因子值增大的直接原因，R因子与土壤侵蚀模数呈明显正相关。

3.3.2 土壤侵蚀与土壤类型 由表5可以看出，栗褐土所占的土壤侵蚀面积最大，为4 199 km2，栗褐土的土壤侵蚀量也最高。从研究区土壤类型分布图中，我们可以清楚的看到，研究区栗褐土的面积为46.36%，栗褐土区为半干旱一年一熟旱作区，分布有大量的耕地，但是耕地的质量较差，绝大部分为坡耕地，且坡度大的耕地占很大比例，耕地的输出少，只用不养，导致农田生态恶化，土壤干旱、贫瘠，侵蚀严重，再生能力低。因此在栗褐土分布区，应重点防范水土流失。

图2 引黄入晋北干线沿线地区各年份R因子值与土壤侵蚀模数

表4 引黄入晋北干线沿线地区各年份不同地域土壤侵蚀变化统计

表5 引黄入晋北干线沿线地区各年份不同土壤类型土壤侵蚀变化统计

3.3.3 土壤侵蚀与坡度 以全国《土壤侵蚀强度划分标准(SL190-2007)》为坡度分类标准，将坡度重分类为6类：<5°，5°～8°，8°～15°，15°～25°，25°～35°，>35°。据相关文献可知，地表类型按坡度可以分为平坡(0～5°)、缓坡(5°～15°)、斜坡(15°～25°)、陡坡(25°～35°)、急坡(35°～45°)、险坡(≥45°)[26]，根据面积分布情况，可以看出，研究区以平坡和缓坡为主，斜坡为辅，急坡和险坡占很小比重。为了进一步探讨坡度与土壤侵蚀量之间的关系，利用空间叠加工具将坡度图层与侵蚀模数图层叠加，计算研究区2005，2010，2015年6类坡度下的土壤侵蚀量模数及土壤侵蚀量比如表6所示。

表6 引黄入晋北干线沿线地区各年份不同坡度土壤侵蚀变化统计

从表6中可以看出，2005，2010，2015年在不同坡度上的土壤侵蚀格局变化不大，不同坡度与土壤侵蚀关系密切。

从侵蚀面积来看，8°～15°坡度带最大，这主要与研究区缓坡分布面积大有关。从平均侵蚀模数看，当坡度大于25°时土壤侵蚀模数较高，尽管此带的分布面积小，但侵蚀量大，因此侵蚀模数高，侵蚀级别2005，2015年达到强度侵蚀，2010年为极强度侵蚀，主要是因为在北方，25°以下的区域属于缓坡区域，退耕还林工程的实施使得该区域成为人工林和梯田改造后的农田，土壤侵蚀状况稍有改善，因此低于25°以上的区域的土壤侵蚀模数。从各坡度带的土壤侵蚀模数看，在3 a中各坡度带的值均为2010年最大，并且，在每1 a中，坡度小于35°时，随着坡度的增加，土壤侵蚀模数增大。总体而言，25°以上的坡度带应成为重点预防和治理的区域。

3.3.4 土壤侵蚀与土地利用类型 从土地利用类型可以看出，研究区大部分地区分布着耕地、林地和草地，将土地利用类型与土壤侵蚀强度进行叠加可以得出如表7所示结果。

表7 引黄入晋北干线沿线地区各年份不同土地利用类型土壤侵蚀变化统计

(1) 不同土地利用类型之间，土壤侵蚀状况差异显著，3 a中土壤侵蚀面积均为耕地最大，其次为草地、林地、建筑用地、水体、未利用土地。耕地以轻度侵蚀为主，这是因为陡坡的耕地面积大，在黄土高原地区，耕地一般沿着等高线带状耕作[4]，耕作方式比较粗放，缺少水土保持的措施，加之人类活动频繁干扰，导致耕地侵蚀。

(2) 林地和草地以轻度和中度侵蚀为主，这是因为林地和草地大都分布在海拔1 000 m以上的山区，林地被人为毁林开荒，草地基本为自然草地，当降雨发生时，这些未采取水保措施的区域，在降雨和陡坡的作用下，会不同程度的产生土壤侵蚀。

(3) 建筑用地的水土保持良好，而水体基本无土可蚀，未利用土地所占面积虽小，但土壤侵蚀强度高，在山阴县南部和怀仁县北部，有部分区域地表裸露，无植被覆盖且坡度大，土壤侵蚀严重。

由以上分析可知，区域内耕地的侵蚀面积大，林地、草地、未利用土地的侵蚀级别高，因此应加强以上土地利用类型的合理利用和管理规划，制定相应的水土保持措施，改善区域内的环境状况。

3.4 土壤侵蚀的人为因素作用

2005—2010年土壤侵蚀模数增长将近1倍主要是由于工程建设时期的一系列施工活动会对地表进行开挖、扰动、再塑，还会损毁水土保持设施，造成大面积的土地裸露，加重水土流失。具体到各因子中主要表现为： ①施工过程中，对地表的扰动再塑，会形成大量的人工坡面和悬空面，在短期内坡度、坡长等地形因子会发生剧烈的变化； ②项目建设过程中会产生大量的弃土弃渣，这些物质稳定性差，在无植被覆盖的情况下，极易产生风蚀和水蚀； ③对地表覆盖物(林地、草地)的清除或砍伐，会使得防风固沙的能力下降，极易引起风力侵蚀，并且大量的土地直接裸露，土壤直接遭受雨滴击溅和冲刷，极易引起水力侵蚀。

2010—2015年土壤侵蚀模数下降主要是由于工程建成后，实行了一系列的工程措施和生物措施来恢复该地的生态环境。如坡面绿化、库岸治理护岸、道路护坡固土、整地措施、抗旱造林措施、拦截措施等有效的减缓了水土流失[15]。

4 讨论与结论

4.1 讨 论

本研究仍存在一定的不足：本文的研究结果仍需水文站实测的泥沙数据进行验证，各因子的数据精度不一，难免在栅格计算过程中产生误差，例如降雨侵蚀力的计算只采用了4个雨量站点的数据，如果采用更多的雨量站点数据，可以进一步提高计算精度。P因子的赋值不确定性因素较高，如有定量计算P因子的方法，计算的结果将更为客观。未来的研究将根据研究区背景特征修正模型中因子的参数，结合土壤允许流失量的阈值来对研究区土壤侵蚀风险进行评估。

4.2 结 论

(1) 2005—2015年，引黄入晋北干线所经地区土壤侵蚀模数呈先增后减，总体减小的趋势。2005,2015年属于轻度侵蚀，而2010年属于中度侵蚀，大于水利部公布的黄土高原土壤允许流失量1 000 t/(km2·a)的标准。引黄入晋北干线的实施是2010年土壤侵蚀模数增长近1倍的重要原因，表明人类活动的影响在土壤侵蚀的研究中不可忽视。

(2) 研究区土壤侵蚀呈块状侵蚀和点状侵蚀相间分布的特点。以北干线为界界线两侧土壤侵蚀差异显著。土壤侵蚀的空间分布与研究区的地貌单元关系密切，侵蚀中低山区高，冲洪积平原低。朔州市辖区的土壤侵蚀面积和土壤侵蚀量最大，朔州市辖区、山阴、怀仁、大同市辖区的土壤侵蚀模数均与总的土壤侵蚀模数分布特点一致。

(3) 各因子中，R因子与土壤侵蚀模数呈明显正相关，当土壤类型为栗褐土、坡度为8°～15°、土地利用类型为耕地时土壤侵蚀面积最大，当坡度大于25°、土地利用类型为草地和林地时土壤侵蚀量最高。