陈 浩 ，张晓萍，权 伟，郭晋伟，袁 阅张思瑞

(1.延安大学 生命科学学院，陕西 延安 716000；2.陕西省区域生物资源保育与利用工程技术研究中心，陕西 延安 716000；3.西北农林科技大学 水土保持研究所 黄土高原土壤侵蚀与旱地农业国家重点实验室，陕西 杨凌 712100；4.乾县农业农村局，陕西 乾县 713300）

土壤侵蚀是全球面临的突出生态环境问题之一［1］,防治水土流失是黄河流域生态治理的根本措施［2］。1999 年开始实施退耕还林(草)工程以来,在气候变化不显著的背景下,黄河中游地区土壤侵蚀与入黄泥沙状况发生了显著变化［3-4］。土壤侵蚀是各种因子综合作用的结果,研究土壤侵蚀时空变化及其与影响因子的关系对防治土壤侵蚀具有重要意义［5］。北洛河流域是黄河中游粗泥沙主要来源区之一,也是全国实施水土流失治理和退耕还林(草)工程的重点地区,笔者采用通用土壤流失方程RUSLE,结合GIS 和RS 技术,定量估算北洛河上游流域实施退耕还林(草)工程前后不同时段土壤侵蚀模数和土壤侵蚀量,分析土壤侵蚀变化与土地利用、植被覆盖、地形等的关系,以期为该区域水土保持生态建设及相关决策提供依据。

1 数据与方法

1.1 研究区概况

北洛河流经陕西、甘肃两省的5 个地(市)18 个县(区),于陕西省大荔县注入渭河,属黄河二级支流。北洛河流域属温带大陆性季风气候区,多年平均降水量459.1 mm,5—9 月降水量占全年的74.5%,降水空间分布不均,东南部降水较多、西北部降水较少。吴旗水文站控制的北洛河上游流域属黄土丘陵沟壑区,是黄河中游粗泥沙集中来源区之一,流域面积3 408 km2,土壤以黄绵土为主,土层深厚、土质疏松、抗蚀性较差,沟谷密度为3.0～4.5 km/km2［6］,梁峁起伏、沟壑纵横,水力侵蚀极为强烈,年均土壤侵蚀模数7 006 t/(km2·a),河源区一带土壤侵蚀模数可达10 000 t/(km2·a)。自1999 年开始实施退耕还林(草)工程以来,流域内的林草覆盖面积显著增加,形成了以落叶阔叶及灌木草丛为主的次生植被类型［7］,因沟道比降较大而淤地坝工程较少且规模较小,有效拦沙库容非常有限,目前大都基本淤满［8］,因此人工植被建设是北洛河上游流域水土保持生态建设的主要措施。

1.2 数据来源

本研究采用的主要数据及其来源:①研究区域及周边分布的吴起、盐池、定边、靖边、志丹、安塞、华池、环县等8 个雨量站1981—2010 年逐日降水量数据,来源于国家气象信息中心；②吴旗水文站1981—2010 年逐日径流泥沙数据,来源于黄河水文年鉴及中国科学院水土保持研究所馆藏数据；③1 ∶50 万土壤属性数据,下载于黄土高原科学数据中心(http://loess.data.ac.cn/),不同土壤类型(黄绵土、黑垆土、红土、新积土)的理化性质数据从《陕西土壤》［9］中查得；④地形数据(分辨率为30 m 的ASTER GDEM),下载于地理空间数据云(http://www.gscloud.cn/)；⑤1990 年、2000 年、2010 年Landsat TM5 遥感影像数据(分辨率为30 m),来源于地理空间数据云；⑥1990 年和2000年土地利用(土地利用类型分为耕地、草地、林地、水域、建设用地)数据下载于黄土高原科学数据中心,2010 年土地利用数据通过解译Landsat TM5 遥感影像并经实地调查验证后得到［10］。

1.3 通用土壤流失方程

在众多土壤侵蚀模型中,通用土壤流失方程(RUSLE)因结构简单、计算简便、精度较高而得到广泛应用［11-13］,其形式为

式中:A为土壤侵蚀模数,t/(hm2·a)；R为降雨侵蚀力因子,MJ·mm/(hm2·h·a)；K为土壤可蚀性因子,t·hm2·h/(hm2·MJ·mm)；L为坡长因子、S为坡度因子,通常综合考虑坡长、坡度对土壤侵蚀模数的影响,把LS称为地形因子；C为植被覆盖因子；P为水土保持措施因子。

(1)降雨侵蚀力因子R值的确定。采用逐日降雨资料估算年降雨侵蚀力的简易算法［14］,得到研究区8个雨量站1981—2010 年降雨侵蚀力,然后通过克里金法进行空间插值,得到研究区逐年降雨侵蚀力。根据研究区1981—2010 年降雨侵蚀力与吴旗水文站年输沙量双累积曲线斜率变化情况,可知该曲线在2000 年出现较明显拐点［15］。考虑到降雨侵蚀力年际波动较大的情况,将研究期1981—2010 年分为1981—1990年、1991—2000 年和2001—2010 年3 个时段(编号分别为I、Ⅱ、Ⅲ,把1990 年、2000 年、2010 年作为3 个时段的代表年)。3 个时段研究区年均降雨侵蚀力分别为1 093.72、1 390.20、1 260.74 MJ·mm/(hm2·h·a),表明1981—2010 年研究区降雨侵蚀力呈现波动上升趋势。

(2)土壤可蚀性因子K值的确定。根据第一次全国水利普查所得水土保持专项调查成果(土壤有机碳、土壤颗粒组成数据等)［16］,采用Williams 等［17］提出的公式进行计算,并基于黄土高原径流小区资料进行修正后得到研究区土壤可蚀性因子K值为0.039 8～0.055 2 t·hm2·h/(hm2·MJ·mm),平均值为0.051 7 t·hm2·h/(hm2·MJ·mm)。

(3)地形因子LS值的确定。根据研究区DEM 数据,采用张宏鸣等［18］基于GIS 提取区域坡长坡度因子的算法进行计算,并经过必要的尺度变换修正［19］,得到研究区LS值范围为0.01～71.61,平均值为12.72。

(4)植被覆盖因子C值的确定。根据研究区土地利用数据和基于归一化植被指数NDVI分布图,在前人研究的基础上,得到研究区不同土地利用类型及其植被覆盖度的C值,其中:耕地C值为0.44(根据研究区常见旱地农作物谷子、玉米、大豆、马铃薯C值分别为0.53、0.28、0.51、0.47［20］,其种植面积比例为1 ∶2 ∶1 ∶6,经加权平均计算得到)；林地、草地的C值按照江忠善等［21］建立的黄土丘陵区人工林、草地植被覆盖度与C值的关系式进行计算(非固定值)；建设用地和水域C值为0。利用ArcGIS 软件得到研究区C值的分布。1990 年、2000 年、2010 年研究区平均C值分别为0.548、0.402、0.222,呈现持续减小趋势,原因是研究区植被覆盖度持续提高从而使C值持续减小。

(5)水土保持措施因子P值的确定。通过Landsat TM5 遥感影像难以提取水土保持措施,故利用研究区土地利用类型数据,通过赋值的方法对水土保持措施因子进行赋值。参考前人研究成果［10］,对耕地、林地、草地P值分别赋值为0.31、0.05、0.16,对建设用地及水域P值均赋值为1,然后利用ArcGIS 软件获得研究区P值的分布。1990 年、2000 年、2010 年研究区P值分别为0.220 7、0.220 5、0.157 0,相比较而言,2010 年较2000 年明显减小,原因是2000 年以后研究区耕地被大面积转化为林地、草地。

2 结果与分析

2.1 土壤侵蚀时空变化情况

按前述RUSLE 分别计算出研究区时段Ⅰ、Ⅱ、Ⅲ土壤侵蚀模数［单位换算为t/(km2·a)］,依据《土壤侵蚀分类分级标准》(SL190—2007),利用ArcGIS 软件获得研究区3 个时段各级土壤侵蚀空间分布情况(见图1),3 个时段各级土壤侵蚀面积及占比统计见表1。3 个时段平均土壤侵蚀模数分别为8 612.62、8 078.91、3 060.71 t/(km2·a),呈持续减小趋势。由表1 可知,3 个时段中度土壤侵蚀面积及占比差别相对较小,而强烈及以上土壤侵蚀面积持续减小、占比持续降低(时段I 面积为2 052.90 km2、占比为60.24%,时段Ⅱ面积为1 925.59 km2、占比降为56.50%,时段Ⅲ面积减小到645.25 km2、占比降至18.93%),轻度及微度土壤侵蚀面积持续增大、占比持续提高。时段Ⅲ强烈及以上土壤侵蚀面积大幅减小、轻度及微度土壤侵蚀面积显著增大,这与1999 年开始实施的退耕还林(草)工程密不可分,退耕还林(草)工程的实施使研究区植被覆盖度大幅提升,尤其研究区东南部的吴起县植被恢复最为明显,因而土壤侵蚀模数在空间上呈现西北高、东南低的空间分异特征。

表1 研究区3 个时段各级土壤侵蚀面积及占比统计

图1 研究区各时段土壤侵蚀强度分布情况

2.2 不同土地利用类型的土壤侵蚀变化情况

将各时段土壤侵蚀强度分布图与土地利用图进行叠加,可得到各时段不同土地利用类型的土壤侵蚀状况。鉴于3 个时段水域、建设用地面积占比均较小(分别为0.08%～0.51%、0.06%～0.52%),且本研究在对水域、建设用地进行土壤侵蚀判读时认定为不发生土壤侵蚀,因此只统计了研究区3 个时段耕地、林地、草地的面积占比及土壤侵蚀模数,见表2 。

表2 研究区各时段主要土地利用类型土壤侵蚀模数

由表2 可知,时段I 和时段Ⅱ土地利用类型以耕地和草地为主,二者约占流域面积的97%；时段Ⅲ土地利用格局发生根本性变化,呈现以草地和林地为主的格局,林地、草地面积占比分别提高至20.95%和63.43%,而耕地面积占比由42.62%～42.69%锐减至14.59%,这种变化是实施退耕还林(草)工程的结果。

研究区3 种主要土地利用类型土壤侵蚀模数大小为耕地＞草地＞林地,从土壤侵蚀模数变化情况来看,从时段Ⅰ到时段Ⅱ耕地土壤侵蚀模数小幅提高、草地土壤侵蚀模数小幅降低、林地土壤侵蚀模数明显降低,从时段Ⅱ到时段Ⅲ耕地、草地、林地土壤侵蚀模数均明显降低。

2.3 不同植被覆盖度的土壤侵蚀变化情况

把植被覆盖度分为低覆盖度(小于10%)、较低覆盖度(10%～30%)、中覆盖度(30%～50%)、较高覆盖度(50%～70%)、高覆盖度(大于70%)5 个级别［22］,各时段土壤侵蚀强度分布图与各级植被覆盖度分布图叠加,可得3 个时段各级覆盖度的土壤侵蚀模数。由表3 可知,时段I 植被覆盖度以低覆盖度和较低覆盖度为主(二者面积占比达89.83%),时段Ⅱ植被覆盖度以较低覆盖度和中覆盖度为主(二者面积占比为86.86%),时段Ⅲ低覆盖度和较低覆盖度面积占比明显降低、较高覆盖度和高覆盖度面积占比明显提升,植被覆盖度的变化与实施退耕还林(草)工程后生态环境质量持续向好相吻合。随着植被覆盖度的提高,土壤侵蚀模数明显降低,原因是植被覆盖度越高抵抗侵蚀的作用越强。从土壤侵蚀模数变化情况来看,从时段I 至时段Ⅱ,在降雨侵蚀力增大的背景下同一植被覆盖度的土壤侵蚀模数有所提高；从时段Ⅱ至时段Ⅲ,相同植被覆盖度的土壤侵蚀模数明显降低,这一方面与气候变化引起降雨侵蚀力有所减弱有关,另一方面主要与植被覆盖度上升有关。

表3 研究区3 个时段各级植被覆盖度土壤侵蚀模数

2.4 不同地形的土壤侵蚀变化情况

2.4.1 不同高程的土壤侵蚀变化情况

结合流域实际情况对高程进行分级,将高程分级图与各时段土壤侵蚀强度分布图进行叠加,统计得到各时段不同高程的土壤侵蚀模数(见表4)。从各级高程面积占比来看,1 350～1 450、1 450～1 550 m 这两级高程的面积占比达76.78%,为研究区主要高程带。从各级高程的土壤侵蚀模数来看,1 350～1 550 m 高程带是研究区土壤侵蚀的主要高程带(土壤侵蚀量约占研究区土壤侵蚀总量的80%)。从各级高程土壤侵蚀模数的变化来看,实施退耕还林(草)工程后的时段Ⅲ土壤侵蚀模数较时段Ⅰ和时段Ⅱ均大幅度降低(平均降幅为62.26%)。

表4 研究区各时段不同高程的土壤侵蚀模数

2.4.2 不同坡角的土壤侵蚀变化情况

把研究区坡角分级(0°～5°、5°～8°、8°～15°、15°～25°、25°～35°、＞35°)图与各时段土壤侵蚀强度分布图进行叠加,统计得到各时段不同坡角的土壤侵蚀模数(见表5)。土壤侵蚀模数随坡角增大而提高,即坡面越陡土壤侵蚀模数越大。从各级坡角坡面面积占比和土壤侵蚀模数综合来看,15°以上坡面是研究区土壤侵蚀的主要区域(土壤侵蚀量约占研究区土壤侵蚀总量的72%)。从各时段土壤侵蚀模数的变化情况来看,时段Ⅲ各级坡角的土壤侵蚀模数较时段Ⅰ和时段Ⅱ显著降低(平均降幅为60.53%)。

表5 研究区各时段不同坡角的土壤侵蚀模数

2.4.3 不同坡向的土壤侵蚀变化情况

按水土流失常规调查中的四分法将坡向分为阴坡、半阴坡、阳坡、半阳坡,把坡向分布图与各时段土壤侵蚀强度分布图进行叠加,统计得到各时段不同坡向的土壤侵蚀模数(见表6)。各时段不同坡向的土壤侵蚀模数均呈现阳坡＞半阳坡＞半阴坡＞阴坡,原因是坡向影响水、热分配,从而影响植被生长状况,进而导致不同坡向土壤侵蚀强度存在差异,阴坡相对阳坡来讲,土壤水分含量大,植被生长迅速,地表覆盖度高,因而土壤侵蚀模数相对较小。从各时段土壤侵蚀变化来看,时段Ⅲ各坡向的土壤侵蚀模数较时段Ⅰ和时段Ⅱ均显著降低(平均降幅为60.49%)。

表6 研究区各时段不同坡向的土壤侵蚀模数

2.5 土壤侵蚀量变化情况

依据上述各时段主要土地利用类型面积占比和土壤侵蚀模数计算的年均土壤侵蚀量见表7。从各时段年均土壤侵蚀量变化情况来看,从时段I 到时段Ⅱ,耕地的年均土壤侵蚀量变化很小(约增大0.2%),占土壤侵蚀总量的比例由57%提高到61%；草地的年均土壤侵蚀量减小14%,其占土壤侵蚀总量的比例由42%下降到39%；林地的年均土壤侵蚀量减小40%,其占土壤侵蚀总量的比例由0.6%下降到0.4%。从时段Ⅱ到时段Ⅲ,耕地的年均土壤侵蚀量减小77%,其占土壤侵蚀总量的比例由61%下降到37%；草地的年均土壤侵蚀量减小40%,其占土壤侵蚀总量的比例由39%提高到62%；林地的年均土壤侵蚀量增大58%,其占土壤侵蚀总量的比例由0.4%提高到1.6%。时段Ⅲ耕地土壤侵蚀量及其占比大幅度下降、林草地土壤侵蚀量及其占比有所提高的原因是,实施退耕还林(草)工程后耕地面积显著减小、林草地面积显著增大。

表7 各时段主要土地利用类型年均土壤侵蚀量 万t/a

2.6 讨论

研究结果显示,时段Ⅰ、时段Ⅱ、时段Ⅲ研究区平均土壤侵蚀模数分别为8 612.62、8 078.91、3 060.71 t/(km2·a),可以看出,实施退耕还林(草)工程前从时段Ⅰ到时段Ⅱ土壤侵蚀模数稍有降低,而实施退耕还林(草)工程后的时段Ⅲ土壤侵蚀模数大幅度下降,表明退耕还林(草)工程的实施在北洛河上游流域取得了显著的水土保持效益,这与刘文超等［23］对陕北地区退耕还林还草工程保土效应的研究结果相一致。退耕还林(草)工程实施以后,北洛河上游流域土壤侵蚀强度在空间上呈现西北高、东南低的特征,原因是东南部地势较低,气候较为湿润,植被恢复效果明显,而西北部多为河源区,山高坡陡,水热条件较差,植被恢复难度大甚至在局部地区植被略有退化［24］,因此植被恢复效果较好的东南部土壤侵蚀模数明显低于西北部土壤侵蚀模数。

研究区主要土地利用类型土壤侵蚀模数遵循耕地＞草地＞林地的规律,退耕还林(草)工程的实施,使土地利用格局与植被覆盖度发生显著变化(耕地面积占比大幅度下降、林草地面积占比显著提高),土壤侵蚀模数随植被覆盖度的提高而降低,这与有关学者［10,25］的研究结果一致。值得注意的是,耕地是土壤侵蚀最为严重的土地利用类型,其土壤侵蚀模数在退耕还林(草)工程实施前的时段Ⅰ为11 481.84 t/(km2·a)、时段Ⅱ为11 528.87 t/(km2·a),时段Ⅱ较时段I 高的原因是降雨侵蚀力增大；在退耕还林(草)工程实施后的时段Ⅲ,土壤侵蚀模数大幅度降低到7 745.10 t/(km2·a),其原因除时段Ⅲ降雨侵蚀力略有减小外,主要是随着退耕还林(草)工程的实施耕地的平均坡度减缓,即土壤侵蚀严重的陡坡耕地转化为林草地,而保留的耕地为土壤侵蚀强度较低的缓坡地。

地形是影响土壤侵蚀的主要因素之一,不同时段土壤侵蚀模数与各地形因子的关系表现出一致性,研究区土壤侵蚀主要发生在1 350～1 550 m 高程带、坡角＞15°的陡坡地、阳坡和半阳坡,这些地带是研究区今后土壤侵蚀防治的重点区域。

土壤侵蚀模数的变化主要受气候条件、土壤、地形、植被覆盖度等因子的影响。土壤和地形为相对静态因子,决定了土壤侵蚀模数的空间分异特征；气候、植被覆盖度为相对动态因子,影响土壤侵蚀模数随时间的动态变化［26］。统计分析表明,1981—2010 年北洛河上游流域降雨侵蚀力年际变化较大,3 a 滑动平均处理结果显示,降雨侵蚀力呈增大趋势,线性倾向率为7.3 MJ·mm/(hm2·h·a)［15］,而土壤侵蚀模数呈下降趋势,尤其是时段Ⅲ土壤侵蚀模数急剧下降,这主要归因于退耕还林(草)工程的实施使植被覆盖度显著提高、土壤侵蚀得到有效控制。

计算的时段Ⅰ、时段Ⅱ、时段Ⅲ年均土壤侵蚀量分别为2 946.4 万、2 763.0 万、1 047.6 万t/a,吴旗水文站实测的时段Ⅰ、时段Ⅱ、时段Ⅲ年均输沙量分别为2 690.4万、4 559.0 万、1 033.8 万t/a,二者相关性较好(相关系数为0.64),但存在一定差别,二者在时段I 和时段Ⅲ比较接近,而时段Ⅱ计算的土壤侵蚀量远小于实测输沙量,其原因:一是采用RUSLE 计算的土壤侵蚀量是坡面侵蚀量,没有考虑沟蚀、重力侵蚀以及泥沙输移过程中的沉积等；二是时段Ⅱ期间的1994 年8 月31 日北洛河流域出现了超过50 a 一遇的特大暴雨,产生了高含沙洪水,该年输沙量是多年平均输沙量的7.3倍［27］,具有明显的峰值特征,因此时段Ⅱ实测年均输沙量远大于计算的年均土壤侵蚀量。

总体而言,采用RUSLE 模型估算的研究区土壤侵蚀模数较为合理。然而,受所采用数据来源及精度的制约,估算结果可能与实际情况存在较大差异(如时段Ⅱ尤其是1994 年)。受Landsat TM5 遥感影像分辨率的制约,本研究在无法获取研究区水土保持措施空间分布信息的情况下,采用土地利用数据间接反映水土保持措施对土壤侵蚀的抑制作用,其与实际情况可能存在差异。相关研究中,大中流域尺度地形因子大都是基于开放免费的低精度DEM 数据提取的［28］,这类数据对地形的表达有限甚至存在地形描述误差,据此提取地形因子存在坡度衰减、坡长增长的问题［29］,并影响土壤侵蚀模数估算精度。因此,今后应加强相关研究,以提高土壤侵蚀模型的精度。

3 结论

(1)北洛河上游流域土壤侵蚀模数在空间上呈现西北高、东南低的分异特征,1981—2010 年在降雨侵蚀力呈现波动上升趋势的情况下,流域平均土壤侵蚀模数呈持续下降趋势,尤其退耕还林(草)工程实施以后土壤侵蚀模数显著下降,时段I、时段Ⅱ、时段Ⅲ土壤侵蚀模数分别为8 612.62、8 078.91、3 060.71 t/(km2·a)。

(2)实施退耕还林(草)工程以后,研究区耕地面积锐减,林草地面积显著增加、植被覆盖度大幅提升,主要土地利用类型土壤侵蚀模数遵循耕地＞草地＞林地的规律,即土壤侵蚀模数随着植被覆盖度的提高呈现下降趋势。

(3)3 个时段土壤侵蚀随地形的分布特征具有一致性,1 350～1 550 m 是研究区土壤侵蚀分布的主要高程带,其侵蚀量约占研究区土壤侵蚀总量的80%,研究区72%的土壤侵蚀量分布在坡角＞15°的坡面,各时段不同坡向的平均土壤侵蚀模数均遵循阳坡＞半阳坡＞半阴坡＞阴坡的规律。