基于RUSLE模型的太行山区土壤侵蚀时空分异特征及影响因子研究*

王 丰,刘金铜,付同刚,高 会,齐 菲

(1.中国科学院遗传与发育生物学研究所农业资源研究中心/中国科学院农业水资源重点实验室/河北省土壤生态学重点实验室 石家庄 050022; 2.中国科学院大学 北京 100049)

土壤侵蚀作为一个重大的全球环境问题,严重地威胁着区域生态环境质量和社会经济发展。土壤侵蚀是一个综合性的研究领域,其研究需要气候、地质、地貌、土壤、水文、生态等相关学科的基本知识。中国是世界上水土流失最严重的国家之一,水土流失面广量大。根据第二次全国土壤侵蚀遥感调查,我国水蚀和风蚀面积为356.9万km,占国土面积的37.6%。因此,土壤侵蚀的研究对生态恢复及可持续发展均具有重要意义。

根据模型估算和预测土壤侵蚀是较为普遍的一种研究手段,依据模型构建方法主要分为经验统计模型和物理成因模型。Wischmeier和Smith最早于1965年建立USLE (Universal Soil Loss Equation)土壤侵蚀模型。1993年,美国农业部正式发布修正的通用土壤侵蚀方程RUSLE,且已经在国内外的土壤侵蚀估测研究中得到了普遍的应用。基于USLE模型,20世纪90年代以来,我国学者也建立了一些经验侵蚀模型,比较有代表性的是江忠善等浅沟裸露地坡面土壤预报模型和刘宝元的中国土壤流失方程CSLE。土壤侵蚀物理成因模型中比较著名的是美国的WEPP、欧洲的EUROSEM、荷兰的LISEM,以及我国学者蔡强国等的次降雨侵蚀产沙模型、段建南等的坡耕地土壤侵蚀过程数学模型SLEMSEP等。但由于此类模型参数复杂和各地区情况千差万别,其推广和应用受到了一定限制,但USLE和RUSLE依然是应用最为广泛的模型。

USLE和RUSLE模型自从引入我国后,在土壤侵蚀较严重的黄土高原地区进行了较为广泛的应用,但对易发生土壤侵蚀的太行山区的研究则相对较少。太行山区是京津冀豫重要的生态屏障区,但由于自然环境恶化以及人类不合理的土地利用,导致该地区土壤侵蚀状况严峻。目前,关于太行山区土壤侵蚀的研究多集中在小流域和某一区域尺度,覆盖整个太行山区尺度的研究相对较少; 此外诸多研究通常只是关注某一年度的土壤侵蚀量,未能很好地反映土壤侵蚀的时间变化情况。因此,本文将RUSLE模型与ArcGIS软件相结合,计算了太行山区2000-2015年4期的土壤侵蚀模数(单位面积年均土壤侵蚀量),分析了土壤侵蚀的时空分异特征,并对土壤侵蚀的影响因子进行了探究,旨在为研究区的水土保持和生态环境治理工作提供一定的科学参考。

1 研究区概况

图1 研究区太行山区地理位置Fig.1 Geographical location of the study area of Taihang Mountain area

2 研究数据与方法

2.1 数据来源

本研究所使用的数据类型及来源如表1所示,各数据均通过ArcGIS栅格重采样工具统一成100 m空间分辨率,投影坐标统一为WGS_1984_UTM_Zone_49N。土壤数据采取FAO-90分类系统; 土地利用数据根据各地类遥感影像进行目视解译和更新,分类综合精度达87.75%。

表1 本研究所用数据及来源Table 1 Dataset and resources used in the study

2.2 研究方法

采用的土壤侵蚀模型是修正的通用土壤流失方程RUSLE模型:

式中:为预测的单位面积的年土壤侵蚀量(t·hm·a);为降雨侵蚀力因子[(MJ·mm)·(hm·h·a)];为土壤可蚀性因子[(t·hm·h)·(MJ·hm·mm)]; LS为坡长坡度因子(无量纲);为植被覆盖与管理因子(无量纲);为水土保持措施因子(无量纲)。

2.2.1 降雨侵蚀力因子()

降雨侵蚀力是降雨侵蚀的潜在能力,它是土壤流失方程中首要的基础因子。降雨是引起土壤侵蚀的主要驱动力,降雨侵蚀力表征了降雨引起土壤发生侵蚀的潜在能力。目前在中国多以年降雨量和月降雨量对因子进行估算,本文采用马志尊的方法计算太行山区:

式中:为年平均降雨量(mm),p为第月降雨量(mm)。根据研究区88个气象站的降雨数据,由以上公式计算每个站点的年降雨侵蚀力,在ArcGIS软件中进行Kriging插值和重采样操作,得到降雨侵蚀力因子的空间分布图(图2)。

图2 2000—2015年太行山区降雨侵蚀力(R)的空间分布Fig.2 Spatial distribution of rainfall erosion factor (R)in Taihang Mountain area from 2000 to 2015

2.2.2 土壤可蚀性因子()

土壤可蚀性因子主要表达土壤的性质,反映土壤被降雨侵蚀力分离、冲蚀和搬运的难易程度。土壤值是土壤抵抗水蚀能力的一个综合性指标,土壤值越大,土壤抵抗水蚀能力越弱,越易受侵蚀;反之亦然。本文基于EPIC模型中土壤质地和土壤有机碳含量计算因子,并根据张科利等的研究对结果进行校正,最终得到因子的空间分布图(图3):

图3 太行山区土壤可蚀性因子(K)的空间分布Fig.3 Spatial distribution of soil erodibility factor (K)in Taihang Mountain area

式中:SAN表示砂粒(粒径0.05～2 mm)质量分数(%); SIL表示粉粒(粒径0.002～0.05 mm)质量分数(%); CLA表示黏粒(粒径＜0.002 mm)质量分数(%);表示有机碳含量(%)。

2.2.3 坡长坡度因子(LS)

坡长是指坡面上从产生地表径流的起点至径流汇集到沟谷的距离,坡长越大,径流的速度就越大,汇聚的流量也越大,其侵蚀力就越强; 坡度表示地表某点的倾斜程度,坡度通过影响水流的速度进而影响渗透率、径流量以及对土壤的冲刷量,是影响区域内土壤侵蚀的最重要因素。在流域层面上,坡长、坡度可通过数字高程DEM提取,坡长因子()的计算采用Wischmeier提出的坡长因子计算公式:

式中:为水平投影坡长(m);为坡长指数,取值方法为:当坡度()≤1º,=0.2; 当1º＜≤3º时,=0.3; 当3º＜≤5º时,=0.4; 当＞5º时,=0.5;利用DEM提取。

RUSLE模型中坡度因子()的计算采用Mccool等建立的公式:

Liu等对陡坡情况下因子的计算做了研究,提出了如下计算公式:

太行山区地形因子LS的空间分布如图4所示。

图4 太行山区坡长坡度因子(LS)的空间分布Fig.4 Spatial distribution of topographical factor (LS)in Taihang Mountain area

2.2.4 植被覆盖与管理因子()

植被覆盖与管理因子是评价植被因素抵抗土壤侵蚀能力的重要指标,指在相同的土壤、坡度和降雨条件下,有特定植被覆盖或田间管理的土地上的土壤流失量与实施清耕、无覆盖裸露休闲地上的土壤流失量之比,其值在0～1之间。归一化植被指数(NDVI)是计算因子最为普遍的数据,本文采用Vander Knijff等提出的利用NDVI确定因子的方法(图5):

图5 2000—2015年太行山区植被覆盖与管理因子(C)的空间分布Fig.5 Spatial distribution of vegetation management and coverage factor (C)in Taihang Mountain area from 2000 to 2015

式中:和为无量纲因子,决定因子与NDVI关系曲线图,经Vander Knijff M发现=2、=1是一个合理的取值,有着良好相关性。

2.2.5 水土保持措施因子()

水土保持措施因子是采取专门措施后的土壤流失量与采用顺坡种植时土壤流失量的比值。因子的控制措施主要有等高耕作、等高带状种植和修筑梯田,常见的水土保持措施有梯田、造林、种草等,其值为0～1,0代表基本不发生土壤侵蚀,1代表未采取水土保持措施或措施完全失效。

大量研究表明,土地利用信息能够间接反映水土保持措施因子,在大尺度的流域土壤侵蚀研究中,可以根据土地利用类型赋值的方法确定,由此得到太行山区水土保持措施因子分布图(图6),林地、草地、未利用地一般均赋值为1,水体、建筑用地及裸岩赋值为0,对于耕地,通常坡度越大,水土保持措施的作用越突出,因此耕地依据表2坡度范围赋值。

图6 2000—2015年太行山区水土保持措施因子(P)的空间分布Fig.6 Distribution of soil and water conservation factor (P)in Taihang Mountain area from 2000 to 2015

表2 不同坡度范围的耕地的水土保持措施因子 (P)值[30]Table 2 Values of soil and water conservation factors (P)of cultivated land in different slope conditions

3 结果与分析

3.1 太行山区土壤侵蚀的时间变化

3.1.1 土壤侵蚀模数的时间变化

将计算所得的降雨侵蚀力因子、土壤可蚀性因子、坡长坡度因子、植被覆盖与管理因子和水土保持措施因子栅格图层统一为100 m×100 m分辨率后,利用ArcGIS中的栅格计算器进行叠加相乘,2000年、2005年、2010年和2015年太行山区平均土壤侵蚀模数分别为4434.14 t·km·a、2984.65 t·km·a、1761.93 t·km·a和1833.81 t·km·a,呈明显减小的趋势,16年间土壤侵蚀模数减少58.64% (图7),土壤侵蚀模数在2010-2015年间整体趋于稳定。

图7 2000—2015年太行山区土壤侵蚀模数变化Fig.7 Changes of soil erosion modulus in Taihang Mountain area from 2000 to 2015

3.1.2 土壤侵蚀强度的时间变化

根据水利部颁布的土壤侵蚀分类分级标准SL190-2007,根据侵蚀模数的大小,土壤侵蚀强度分为微度(＜500 t·km·a)、轻度(500～2500 t·km·a)、中度(2500～5000 t·km·a)、强烈(5000～8000 t·km·a)、极强烈(8000～15 000 t·km·a)以及剧烈(＞15 000 t·km·a)6个级别。首先,利用ArcGIS重分类工具将侵蚀模数分为6种强度,进而利用相交工具得到侵蚀强度的转移矩阵(表3)。太行山区土壤侵蚀强度的面积变化随时间表现各异,微度侵蚀的面积呈持续增加趋势,轻度、中度及强烈侵蚀面积呈先增加后降低的趋势,而极强烈和剧烈侵蚀面积主要表现为减小趋势。

表3 2000—2015年太行山区土壤侵蚀强度转移矩阵Table 3 Transfer matrix of soil erosion intensity in Taihang Mountain area from 2000 to 2015km 2

在侵蚀强度转移1级的面积统计中,2000-2005年、2005-2010年和2010-2015年3个变化时间段轻度侵蚀向微度侵蚀转移的面积比例分别为14.47%、35.76%和30.48%,其中2000-2005年的转移面积比例最小; 3个变化时间段中度侵蚀向轻度侵蚀转移的面积比例分别为33.45%、68.61%和42.29%,其中2005-2010年的转移面积比例最大。2010-2015年轻度侵蚀面积向中度侵蚀转移的比例达13.05%,但其值仍远小于轻度侵蚀转微度的面积。在侵蚀强度转移两级别及以上的面积统计中,2000-2005年没有极强烈、剧烈侵蚀向微度侵蚀转移的面积,而2005年之后的两个时间段均有向微度侵蚀等级转移的面积,且转移的面积呈现增大趋势。

因此,2000-2015年微度侵蚀的面积一直在增加,极强烈和剧烈侵蚀的面积主要表现为降低趋势,太行山区土壤侵蚀强度整体向较低的等级转移,2005-2010年是侵蚀强度降低变化率最高的时间段,太行山区土壤侵蚀防治成效在此时期最为显著。

3.2 太行山区土壤侵蚀的空间变化

3.2.1 土壤侵蚀强度的水平空间变化

本文运用ArcGIS空间叠加分析功能,揭示了太行山区2000-2015年土壤侵蚀强度的空间分异特征(图8)。

图8 2000—2015年太行山区土壤侵蚀强度等级变化(括号内数字为面积占比)Fig.8 Changes of soil erosion intensity levels in Taihang Mountain area from 2000 to 2015 (data in parenthesis is area proportion)

2000-2005年土壤侵蚀强度降低的区域主要分布在山西和河北两省交界的较高海拔地区,东经114º附近,例如晋中东部、石家庄西部等; 土壤侵蚀强度增加的区域总体分布在太行山东北部区域,以北京西南部、张家口东南部以及保定北部为主。2005-2010年土壤侵蚀强度降低的区域在整个太行山区分布较为均匀; 而侵蚀强度增加的区域很少,此时间段研究区土壤侵蚀的状况改善明显。2010-2015年土壤侵蚀强度降低的区域主要分布在阳泉、石家庄以南地区; 土壤侵蚀强度增加的区域整体上集中在北纬38º线以北,例如张家口、保定、石家庄等水土流失敏感性较高的西部山区。

综上所述,太行山区土壤侵蚀具有明显的空间异质性。3个变化时期至少有74.77%的地区侵蚀强度并没有表现出等级上的转移,侵蚀强度降低一级的面积始终比增加一级的大。

3.2.2 土壤侵蚀模数的垂向空间变化

根据高会等对太行山区海拔的分类标准,通过GIS分区统计得到土壤侵蚀的垂直空间分布图(图9)。具体来看,在海拔小于200 m的山前平原区,土壤侵蚀模数很小,平均低至244 t·km·a; 在200～1600 m的海拔区间,主要是低山丘陵区和中山区,土壤侵蚀与海拔未呈现出相关性,平均侵蚀模数始终在2000 t·km·a左右,属于轻度侵蚀; 在海拔大于1600 m的亚高山区,土壤侵蚀模数随着海拔升高表现为迅速增大,在2600～2800 m附近达极大值,即13 057 t·km·a,属于极强烈侵蚀等级。

图9 太行山区不同海拔的土壤侵蚀变化Fig.9 Changes of soil erosion in different altitudes in Taihang Mountain area

因此,在中低山区海拔高的地方侵蚀模数并不一定大,一定海拔高度后土壤侵蚀状况有其特殊性,中低山区集中了侵蚀总量(2.5×10t)的86.54%,原因是中低山区的面积约占研究区的91.12%。土壤侵蚀防治工作应着眼于中低山区,对于亚高山区应减少人为因素的干扰,采取更多的封山育林措施。

3.3 太行山区土壤侵蚀的影响因子

3.3.1 土壤侵蚀与坡度的关系

选取2015年土壤侵蚀与坡度的栅格图层进行叠加分析,得到太行山区土壤侵蚀模数与坡度的关系(图10 a),土壤侵蚀模数总体上随坡度增加而增大,侵蚀模数与坡度呈二次多项式拟合关系=259.18+221.16-2.76,拟合系数=0.96,侵蚀模数与坡度的关系极显著(＜0.01)。当坡度在40º左右时侵蚀模数达极大值,约4693 t·km·a; 而当坡度大于40º时,侵蚀模数并未继续增大,而是呈现波动下降的趋势。因此,土壤侵蚀存在临界坡度使得其不会随坡度的增加一直增加,这也与众多的学者研究相似,基于试验方法计算的临界坡度多在30º以下,而基于理论分析得出的结论大多认为临界值在30º以上,临界坡度的大小因试验方法和客观条件的差异而不尽相同。

根据水利部关于土壤侵蚀坡度等级划分的标准,将太行山区土壤侵蚀栅格图按照坡度划分为6个等级,太行山区土壤侵蚀强度与坡度的关系如图10b所示,在不同坡度带不同土壤侵蚀强度的面积占比变化明显,微度侵蚀的面积占比随坡度的增加而逐渐减小,而中度及以上侵蚀强度的面积占比则逐渐增大。因此,随着坡度的增加,较高等级强度的侵蚀发生的概率逐渐增大。太行山区通过工程措施进行微地形改造,通过调整坡耕地的坡度,有助于控制土壤侵蚀的发生。

图10 太行山区土壤侵蚀模数随坡度变化(a)和不同坡度土壤侵蚀强度面积百分比(b)Fig.10 Changes of soil erosion modulus (a)and percentages of different soil erosion intensity levels (b)under different slope conditions in Taihang Mountain area

3.3.2 土壤侵蚀与土地利用的关系

太行山区土壤侵蚀模数与土地利用的关系如图11所示,耕地、林地和草地是主要的土地利用类型,面积占比分别为35.87%、29.22%和25.58%,三者平均土壤侵蚀模数分别为501.72 t·km·a、2475.46 t·km·a和3605.73 t·km·a。

图11 太行山区不同土地利用类型的侵蚀模数Fig.11 Soil erosion modulus of different land use types in Taihang Mountain area

3种主要的土地利用类型中,耕地的侵蚀模数最小,几乎达到微度侵蚀等级。可能的原因有以下3点:1)太行山区耕地的坡度主要分布在15º以下,统计分析(表4)表明耕地的平均坡度为4.90º,较低坡度使得土壤侵蚀模数相对较小。2)21世纪以来,退耕还林生态工程主要集中在大于15º的耕地上,2000-2015年这些地区平均NDVI植被指数从0.56变为0.68,相应的平均侵蚀模数从5536.12 t·km·a降低为2474.29 t·km·a。3)林草地通常未采取相应的水土保持措施,耕地的水土保持措施因坡度而存在差异。因此,太行山区耕地的土壤侵蚀模数较低。

表4 太行山区主要土地利用类型特征Table 4 Characteristics of main land use types in Taihang Mountain area

与耕地相比,林地和草地主要分布在太行山区较大的坡度范围内,较高土壤侵蚀发生的概率增加。草地是侵蚀模数最高的土地利用类型,原因是一方面草地的NDVI指数远小于林地; 另一方面林地往往具有更稳定的群落结构,冠层对降水具有截留作用,这削弱了雨水到达地表的雨滴动能。

太行山区土壤侵蚀防治工作应着眼于耕地、林地和草地管理。对陡坡或交通不便的耕地适时地推进退耕还林还草,提高生态恢复植被的成活率; 对荒山和草地应加强保护力度,预防草地的生态退化。

3.3.3 土壤侵蚀与NDVI的关系

太行山区土壤侵蚀与NDVI植被指数的关系如图12所示,4个年度的平均NDVI分别为0.58、0.61、0.67和0.66,总体上随时间呈现逐渐增大趋势。侵蚀模数与NDVI的二次多项式拟合关系=141 296-414 662+308 084,拟合系数=0.99,侵蚀模数与NDVI的关系极显著(=0.022)。植被是影响土壤侵蚀的重要因子,关系到作物覆盖与管理因子的大小,由此可见,太行山区植被恢复整体趋好,这对防治水土流失具有一定的积极作用。

图12 太行山区侵蚀模数和NDVI的关系Fig.12 Relationship between soil erosion modulus and NDVI in Taihang Mountain area

4 讨论与结论

4.1 讨论

RUSLE模型的各因子在计算过程中会产生误差,本研究中因子的计算仅使用了月度的降雨数据,未考虑太行山区海拔因素对降雨的影响,总体上计算精度还不够,今后研究中考虑使用Anusplin插值方法得到日降雨数据降雨侵蚀力因子;因子还是根据传统的方法对土地利用类型进行赋值,没有考虑等高耕作和修梯田的影响,今后研究中可以考虑采用高分辨率遥感数据对这些典型区域进行识别。

在结果验证方面也存在不足,本文只是参照了已有的相关文献,太行山区的土壤侵蚀研究结果不尽相同,这可能是数据源和参照标准的不同引起的,但一致的结论是研究区土壤侵蚀状况逐渐好转,自20世纪80年代开始,太行山区相继实施了“太行山绿化工程” “退耕还林工程”等多项生态工程,对研究区生态恢复具有重要作用。王娇等对河北太行山区土壤侵蚀敏感性的研究发现,水土流失敏感性中度偏重,由于不同的研究方法和评价尺度,其与本研究结果存在一定差异。

一般而言,耕地的土壤侵蚀模数大于林地和草地,但本研究的结论则不同。首先,太行山区林地的平均坡度是耕地的3倍以上(表4); 其次,耕地的平均NDVI值与草地的几乎相同; 最后,统计分析表明耕地的降雨侵蚀力因子是三者最小的,这可能是研究区耕地侵蚀模数较小的原因。何莎莎在太行山区小流域尺度的研究也得到相似的结论。因此,太行山区水土流失防治工作应遵循客观规律,宜林则林,宜耕则耕,增加陡坡地(＞15º)林草地的覆盖比例,在宜耕的缓坡地上,避免大规模将耕地向林草地转变,以保持现状为主,以微地形改造等农田生态工程措施保持土壤。

本研究中的土壤侵蚀模数存在年平均NDVI的阈值,其他相关研究也表明植被生长状况的好转对土壤侵蚀的作用并不是无限增大的。土壤侵蚀受植被、土壤、降雨、土地利用、地形、人类活动等一系列因子的影响,根据“木桶原理”,往往一两种因子决定了土壤侵蚀防控的上限,合理的生态工程对太行山区生态环境的改善是全方位的,植被的恢复仅仅是生态工程发挥作用的一种体现。因此,太行山区应继续坚持以生态工程措施改善环境质量,优化调整各要素结构,进一步提升对土壤侵蚀的防控能力。

4.2 结论

本文将RUSLE模型与ArcGIS软件相结合,研究了太行山区2000-2015年4期土壤侵蚀的时空分异特征,分析了坡度、土地利用和NDVI与土壤侵蚀的关系,结论如下:

1)2000-2015年,太行山区4期平均土壤侵蚀模数分别为4434.14 t·km·a、2984.65 t·km·a、1761.93 t·km·a和1833.81 t·km·a,呈现明显的减小趋势。

2)土壤侵蚀强度面积转换变化明显,侵蚀强度降低1级的面积始终大于增加1级的面积,微度侵蚀的面积持续增加,而极强烈和剧烈侵蚀的面积主要表现为减小趋势。研究区土壤侵蚀强度整体呈向较低等级转移的趋势,但大部分地区并没有表现出等级上的变化,2005-2010年是土壤侵蚀强度降低变化率最高的时间段。

3)太行山区土壤侵蚀具有明显的空间异质性,2000-2005年土壤侵蚀强度降低的区域主要分布在山西和河北两省交界的较高海拔地区,2005-2010年土壤侵蚀强度降低的区域在整个太行山区分布较为均匀,2010-2015年土壤侵蚀强度降低的区域主要分布在阳泉、石家庄以南地区。

4)中低山区的面积约占太行山区的91.12%,两个区域发生的土壤侵蚀量占总量的86.54%,亚高山区的侵蚀模数相对偏大,可能的原因是亚高山区的地形较为陡峭,中低山区是太行山区人类活动的主要区域,也是土壤侵蚀防治工作的重点区。

5)土壤侵蚀与坡度具有正相关关系,40º左右是土壤侵蚀的极值坡度,随着坡度的增加,微度侵蚀所占的面积逐渐降低,较高等级强度的侵蚀发生的概率逐渐增大,坡度大的区域具有更大的侵蚀风险。太行山区耕地平均坡度为4.90º,远小于林草地,因此其侵蚀模数相对较小; 草地是侵蚀模数最高的土地利用类别,平均侵蚀模数达3605.73 t·km·a。土壤侵蚀随着NDVI增加逐渐降低,年平均NDVI指数达到0.66左右时侵蚀模数降低幅度变缓。