沙棘造林对吴起县土壤水蚀的影响

郭建英，何京丽，李锦荣，殷丽强，刘铁军，荣 浩

（1.水利部牧区水利科学研究所，内蒙古 呼和浩特 010020；2.水利部水土保持植物开发管理中心，北京 100083）

土壤侵蚀是制约人类生存和社会可持续发展的重大环境问题，是我国各种生态问题的集中反映，又是导致生态进一步恶化和贫困的根源［1］。黄土高原土壤侵蚀是目前世界瞩目的环境问题，严重的水土流失使土地资源遭到破坏、土壤肥力和质量下降，泥沙淤积河道引发洪水灾害，破坏地表植被，恶化生态环境，加剧土地和小气候的干旱程度以及其他自然灾害的发生［2-3］。因此有效地控制水土流失，恢复和保护环境是社会经济可持续发展的需要，为此国家在黄土高原地区实施了大规模的退耕还林还草工程，对于治理水土流失取得了显著的成效［4-5］。而选择什么样的树种才能有效地恢复该地区恶化的生态环境是重中之重，沙棘因其具有抗逆性强、适应性广、生态效益好、开发价值大的诸多特点，在吴起县实施退耕还林工程建设中大面积的种植，占吴起县退耕还林地总面积的82.3%。目前关于沙棘造林对于吴起县土壤侵蚀的影响程度，尚缺乏定量分析和研究。定量评价区域土壤侵蚀量、土壤侵蚀强度及其空间分布特征，对于采取相应措施减少水土流失、保护和合理利用土地资源具有重要的意义［7-9］。随着科学技术的不断发展以及数据源的多样化和完整化，RS和GIS技术在宏观把握水力侵蚀现状和预测其发展趋势方面得到了广泛应用［10］。本研究基于遥感和GIS技术，应用修改后的通用土壤流失模型（RUSLE）对吴起县沙棘造林前后的土壤侵蚀状况进行定量计算。试图通过对吴起县沙棘造林以来土壤侵蚀量、土壤侵蚀强度及其空间分布特征的研究，揭示沙棘造林对县域土壤侵蚀的控制效果，以期为沙棘造林水土保持效益监测与评价技术路线的确定提供参考。

1 研究区概况

吴起县位于延安市西北部，东经107°38′37″～108°32′49″，北纬36°33′33″～37°24′27″，地处毛乌素沙地南缘农牧过渡地带，属黄土高原梁峁沟壑区，县境内有无定河、北洛河两大流域，白于山、子午岭两大山系。全县面积3791.5km2，海拔高度1233～1809m。年平均气温7.8℃，无霜期96～146d，年平均降雨量478.3mm，且64%以上集中在7～9月份，其它季节多为无效降雨，旱灾、雹灾、冻灾、风灾等自然灾害频繁，多年平均年陆面蒸发量为400～450mm，属典型干旱半干旱地区［11-12］。吴起县1998年在全县实行封山禁牧，于1999年实行一次性退耕还林还草，2000年被列为国家退耕还林还草科技示范点，截止到2010年，吴起县经国家确认退耕还林面积116846.66hm2，占吴起县总面积的30.82%，1999年一次性的退耕地还林面积占总退耕还林面积80.24%；退耕还林模式主要有沙棘×山杏混交林、刺槐×沙棘混交林、山杏、沙棘、油松纯林，其中沙棘纯林面积最大，生长较好，种植面积占吴起县退耕还林总面积的82.3%［13］。

2 研究方法

2.1 数据来源与处理

遥感数据与地理信息系统相结合为定量化研究土壤侵蚀的时空变化提供了技术支持。本研究的基础数据包括：遥感影像包括多波段的Landsat5 TM（p128r34）,成像时间为1998年7月30日和2010年8月3日（空间分辨率30m）；吴起县行政边界,依据国家基础地理信息系统网站提供1∶400万的全国县界电子图进行确定；DEM数据及由其产生的坡度、坡长等数据，采用国家测绘局提供的1∶5万电子数字化图生成；土壤数据来自陕西省第2次土壤普查数据［14］，并依据当地的土壤分布图（1∶50万）和土地利用方式对各种类型土壤于2009年7月进行了实地取样，野外采样时取0～20cm的表层土壤，每一采样点周围取3个点，混合土样，四分法取样,采样的同时，利用GPS获取土壤采样点的地理坐标共采集土样165个，带回实验室对土壤有机质含量进行重铬酸钾法测定及土壤机械组成进行吸管法测定。

利用ArcGIS10.0以及Erdas9.1等地理信息系统软件，对遥感数据进行影像校正、影像镶嵌、影像裁剪、影像增强和融合等预处理，将预处理之后的遥感影像，用于解译土地利用/覆盖信息；在获得植被、地形、土地利用等基础数据后，建立包括矢量和栅格两种形式的地理信息数据库，对其属性数据进行相应的数据编码，再借助ArcGIS10.0将上述各类型图层栅格化（Grid格式），考虑到县域尺度和计算精度，栅格大小为10m×10m，获得通用土壤流失方程（USLE）的因子图层，再通过Arc-GIS10.0将各因子图连乘，获得研究区不同时期的土壤侵蚀强度等级及其分布图；根据水利部颁发的土壤侵蚀强度分级标准（SL190-2007）确定土壤侵蚀强度，将本研究区内土壤侵蚀分为微度、轻度、中度、强烈和极强烈侵蚀5种强度级别［15-16］。按照遥感影像判读精度，对照影像预判结果，选择具有代表性的点和区域，结合土壤、地质、气候和实地GPS定位调查结果等因子对解译结果进行综合分析及纠正［17］。

2.2 土壤侵蚀模型选择

目前用来评估土壤侵蚀量的模型有很多，其中最早建立且较为成熟的侵蚀预报模型是Wishmeier等人1965年基于大量小区观测资料和人工模拟降雨试验资料建立的的通用土壤流失方程USLE（Universal Soil Loss Equation）［18］，我国自20世纪80年代以来，开始引进通用土壤流失方程，根据实地观测参数，对通用土壤流失方程进行修订，该模型在我国得到了广泛的应用［3,7-11,19-23］。其方程表达式如下：

式中：A-侵蚀模数［t/km2·a］；R-降雨和径流因子［MJ·m m/km2·h·a］；K-土壤可蚀性因子［t·h/MJ·m m]；L S-坡长坡度因子，无量纲；C-植被与经营管理因子，无量纲；P-水土保持因子，无量纲。

2.3 土壤侵蚀因子计算

2.3.1 降雨侵蚀力R因子

根据研究区气象局降雨资料的分析，该地区水土流失季节为5～10月份，当降雨量≥10mm，最大30min雨强≥4mm时将会引起土壤流失。侵蚀降雨多集中在7、8、9三个月，占5～10月降雨量63.5%～94.6%，占年降水总量51.1%～80.3%，与孙保平、孙立达等人研究降雨侵蚀力的区域基本相似，因此采用其建立的5～10月份降雨量（h）与年侵蚀力回归方程：R=1.77h-133.03（相关系数r=0.85）［24-25］。按以上公式计算，吴旗县1999年R=537.5MJ·mm/hm2·h·a，2007年R=570.2MJ·mm/hm2·h·a。

2.3.2 土壤可蚀性K因子

K因子反映了土壤对侵蚀的敏感性及降水所产生的径流量与径流速率的大小。影响K因子的因素很多，但一般来说主要与土壤的粒径和有机质有关。

依据Ei-Swaify等的研究，K值的大小与土壤质地有较高的相关性。本文采用吕喜玺等人对Williams等人在EPIC棋型中土壤可蚀性因子K的改进公式［26］。

式中：Sa是砂粒（0.1～2mm）含量%，Si是粉粒（0.002～0.1mm）含量%，Cl是黏粒（＜0.002mm）含量%，C是有机碳含量%，Sn=1-Sa100。

利用陕西省1∶50万土壤分类图，将野外采样点所得的土壤粒径组成和有机质含量数据输入上式计算各采样点的K值，采用克里格插值法得到吴起县土壤可蚀性K值分布图。

2.3.3 坡度坡长因子的计算

①坡长因子L的计算

利用DEM栅格数据，利用Arcgis提取坡长λ。利用L和坡长的关系求取坡长因子。

式中λ——坡长（m）；m——坡长效应指数，当坡度＜l%时，m=0.2；坡度1%～3%，m=0.3；坡度3%～5%，m=0.4；坡度＞5%，m=0.5。

②坡度因子S的计算

缓坡上选用McCool等（1987年）研究的坡度公式［27］，见式［2-1］和式［2-2］；在陡坡上采用刘宝元（1994年）的坡度公式［37］，见式［2-3］。

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

按照上述方法和计算公式，计算得到吴起县坡长坡度因子LS图。

2.3.4 植被与经营管理C因子的计算

C因子反映的是所有有关覆盖和管理变量对土壤侵蚀的综合作用，其值范围在0～1之间。对于C而言，基本没有土壤侵蚀危险的地区被赋予0；1值被赋给那些最容易受到侵蚀的地区，如裸地；其值不仅取决于具体的作物覆盖、轮作顺序及管理措施的综合作用等，还取决与作物不同生长期侵蚀性降雨的多寡，特别是在作物覆盖与田间管理措施所提供的水土保持作用最差时段内侵蚀性降雨的数量，C值的取值主要与植被覆盖和土地利用类型有关［28-29］。

①植被盖度的提取

根据安培浚等在西北干旱区利用NDVI所建立的植被指数与植被覆盖度的关系来确定植被覆盖度c［39］。利用NDVI计算植被覆盖度c的公式为：

其中：NDVI为所求像元的归一化植被指数。NDVImin和NDVImax分别为研究区NDVI的最小值和最大值。

②C因子的确定

C值主要受植被覆盖度和土地利用现状的制约，本次研究将参考国内外的研究成果［21,31-32］，结合吴起县的实际情况利用下式计算得出各地类年均C因子。利用ARCGIS中的raster calculator，采用下面公式计算得到研究区退耕还林前后C因子图。

2.3.5 保护措施因子P的确定

P因子指采用专门措施后的土壤流失量与采用顺坡种植时的土壤流失量的比值。国内P值的获得基本上是根据区域特点对土地的不同利用方式赋值。最终参照齐述华等人的研究成果［21,32-33］，结合遥感解译分类得到的土地利用现状图进行P因子赋值。坡耕地在实施沙棘造林工程后，农业用地基本上属于缓坡等高耕作和梯田，侵蚀量比较大而且没有任何保护措施的坡耕地基本都实施了沙棘造林，部分转化为草地。所以可以据此对不同土地利用下的P因子赋值：退耕前林地、高覆盖度草地和中覆盖度草地的均赋值为1，建筑用地和水域的赋值为0，而农田的则赋值为0.3；退耕后考虑到吴起县从1999年的退耕还林工程的实施，＞25°的坡耕地全部退耕地为林草地，转化之后的有林地和灌木林地均有侵蚀防治措施，有林地和灌木林地进行反坡整地、鱼鳞坑等。所以，在2010年的土地利用新增有林地和灌木林地的P值取0.8，未变化有林地和灌木林地取1，退耕后耕地基本为平坦的旱地和梯田，P值取0.35，其他土地利用类型P值沙棘造林前后无变化（见表1）。

表1 吴起县沙棘造林前后P值变化

3 结果与分析

3.1 沙棘造林前后土壤侵蚀变化分析

利用Arcgis10.0的空间分析模块，执行图形叠加运算功能，生成土壤侵蚀量图层，在栅格土壤流失量图的基础上，根据水利部颁发的土壤侵蚀强度等级划分标准（表2），对栅格土壤流失量进行分类，将侵蚀量在同一侵蚀等级的栅格进行合并，得到吴起县沙棘造林工程建设前后的土壤侵蚀强度分布图（图1、图2），由于其单位为英制，需要进行单位转换，乘以系数224.2，即可转换为t/km2·a的公制单位，得到各像元的年土壤流失量。

表2 土壤侵蚀强度分级标准

图1 沙棘造林前土壤侵蚀强度分布图

图2 沙棘造林后土壤侵蚀强度分布图

表3 吴起县沙棘造林前后土壤侵蚀强度面积百分比分级表 单位：t/km2·a

从图1、图2及表3上可以明显看到，与退耕前的1998年相比，土壤侵蚀侵蚀强度总体上降低，较高强度土壤侵蚀的土地面积在大幅度减少。极强烈、剧烈的高等级土壤侵蚀面积2010年比1998年分别减少了74.58%和78.93%，使得极强烈、剧烈的高等级土壤侵蚀面积由1998年占国土总面积的56.67%下降到2010年的13.38%；依据表3的统计结果可知：吴起县2010年的全县平均土壤侵蚀为5285t/km2·a，属于强烈侵蚀，其研究结果与张岩等人利用中国土壤侵蚀模型（CSLE）对吴起县2009年度县域平均土壤侵蚀强度5504t/km2·a的研究结果基本相似［43］。总体来看，吴起县沙棘造林前以极强烈侵蚀为主，占吴起县土地总面积的42.76%，沙棘造林后的2010年以中度侵蚀为主，占吴起县土地总面的43.62%。由此可见吴起县沙棘造林后水土流失得到了有效遏制，生态建设成果突出，但土壤侵蚀状况仍然较为严重，需要继续加强巩固退耕还林工程建设成果，充分发挥沙棘造林工程建设的水土保育效应。

3.2 土壤侵蚀强度的空间分布

对所获得的吴起县沙棘造林前后土壤侵蚀分级的结果进行统计分析，如图3，研究区沙棘造林后土壤侵蚀分布以中度侵蚀最多，其面积占吴起县国土总面积的43.62%；其次为强烈侵蚀（27.35%）；而沙棘造林前吴起县土壤侵蚀以极强烈侵蚀为主，占吴起县国土总面积的42.76%,其次为强烈侵蚀，占吴起县国土总面积的31.48%,由此可见，沙棘造林工程建设对于控制地区水土流失起到积极的促进作用，改善效果明显。

图3 吴起县沙棘造林前后各侵蚀级别面积统计

对照土壤侵蚀强度分布图1、图2及研究区坡度图，可以看出，各土壤侵蚀级别的空间分布具有一定的规律性。其中，无明显侵蚀的分布集中分布在塬面、梁顶等坡度较小、地势较平缓的地段。轻度侵蚀多分布于无明显侵蚀地段的周围的塬坡或梁坡上，这些地段多在0°～5°之间。另外沿沟道平缓地也有轻度侵蚀分布，目前多种植以沙棘+紫花苜蓿为主的灌草复合林，植被覆盖度一般在45%以上。中度侵蚀面积大，成片分布在沟坡、沟边、塬坡或梁坡上，坡度一般在15°～20°之间，目前土地利用以草地为主，植被覆盖度在40%～50%之间。强度侵蚀区域的坡度更陡，一般在20°～25°之间，主要分布于沟谷和塬坡或梁坡的陡坡上，以及较大的沟内陡坡和沟床的斜坡面上，坡度多大于30°，植被稀疏，甚至完全无植被覆盖，这是泻溜、崩滑重力侵蚀现象的多发区域。产生这种现象的原因是由于土壤侵蚀严重的地区植被覆盖度、水土保持管理措施并不十分理想。这说明坡度、植被覆盖度是影响吴起县土壤侵蚀的主要因素。

4 结论与讨论

沙棘造林工程是吴旗县实施范围最广、投资量最大的生态建设工程，其改变了不合理的土地利用方式，改善了区域内的环境条件，防治或减轻自然灾害，控制水土流失，形成有利于人类与动植物生存的生态环境，保障了该地区的生态安全体系建设。本文基于3S技术，利用修改后的通用土壤流失模型（RUSLE），结合野外实地调查结果，客观地反映了吴起县沙棘造林生态建设前后土壤侵蚀的变化过程，结果表明：吴起县实施沙棘造林10年来，植被恢复明显，水土流失得到了有效遏制，生态建设成果突出。

（1）吴起县实施沙棘造林后10年来，水土流失控制效果明显，较高强度土壤侵蚀的土地面积在大幅度减少。研究区沙棘造林10年后土壤侵蚀分布以中度侵蚀最多，其面积占吴起县国土总面积的43.62%，而沙棘造林前吴起县土壤侵蚀以极强烈侵蚀为主，占吴起县国土总面积的42.76%；吴起县沙棘造林10年后的全县平均土壤侵蚀为5285t/km2·a，属于强烈侵蚀，其沙棘造林前1998年的平均土壤侵蚀为9779t/km2·a，为极强烈侵蚀，为沙棘造林后的1.85倍，每年约减少土壤侵蚀量1402万t。由此可见，沙棘造林对于吴起县控制水土流失和改善生态环境有着重要作用。

（2）研究区土壤侵蚀主要与坡度和土地利用/覆被有关，微度侵蚀的分布集中分布在塬面、梁顶等地势较平缓的造林地段；轻度侵蚀多分布于微度侵蚀地段的周围的塬坡或梁坡上；中度侵蚀分布在沟坡、沟边、塬坡或梁坡上，坡度一般在10°～15°之间，土地利用以草地和灌木林地为主，植被覆盖度在40%～60%之间；强烈侵蚀区域的坡度更陡，一般在20°～25°之间，主要分布于沟谷和塬坡或梁坡的陡坡上，以及较大的沟内陡坡和沟床的斜坡面上；极强烈多为未利用土地，坡度多大于30°，植被稀疏，甚至完全无植被覆盖，为泻溜、崩滑重力侵蚀现象的多发区域，研究区还尚需加强对中度以上侵蚀区的治理，才能更加有效地控制该地区的水土流失。

（3）应用先进的3S技术与修改后的通用土壤流失模型（RUSLE）相结合，可以快速评估沙棘造林及其林业生态工程建设对于县域水土流失的治理成效，可预测出每个栅格的土壤侵蚀量,便于管理者对较为严重的土壤侵蚀区有针对性提出最佳治理措施，节约大量人力、财力和物力，具有较强的客观性。

［1］李锐.中国土壤侵蚀研究的机遇与挑战［J］.自然杂志,2007,30（1）：6.

［2］于国强,李占斌,裴亮.不同植被类型下坡面径流侵蚀产沙差异性［J］.水土保持学报,2012,26（1）：1-6.

［3］王晓慧,陈永富,陈尔学.基于遥感和GIS的黄土高原中阳县土壤侵蚀评价［J］.山地学报,2011,29（4）：442-448.

［4］许智超,张岩,刘宪春,等.半干旱黄土区退耕还林十年植被恢复变化分析［J］.生态环境学报,2011,20（1）：91-96.

［5］郭碧云,王光谦,张正峰.内蒙古农牧交错区土地利用/覆被变化及土壤侵蚀-以太仆寺旗为例［J］.应用基础与工程科学学报,2011,19（5）：9-16.

［6］周湘山,孙保平,李锦荣.等.基于GIS和USLE的土壤侵蚀定量分析研究［J］.水土保持研究,2011,18（4）：5-9.

［7］李忠锋,王一谋,冯毓荪,等.遥感与GIS在准格尔旗土壤侵蚀监测中的应用［J］.遥 感 技 术 与 应用,2002,17（6）：322-326.

［8］王晓慧,陈永富,陈尔学.基于遥感和GIS的黄土高原中阳县土壤侵蚀评价［J］.山地学报,2011,29（4）：442-448.

［9］宫奎方,范建容,张定容,等.RS和GIS技术支持下的昌都县土壤侵蚀评估［J］.中国水土保持,2011,5：51-53.

［10］李萍,朱清科,王露露,等.吴起县耕地面积变化对其产业结构影响［J］.农业工程学报,2012,28（1）：250-255.

［11］王文田,董哲,葛根巴图.吴起县不同退耕还林模式生态效益研究［J］.林业调查规划,2010,35（6）：34-38.

［12］郭建英.吴起县退耕还林工程效益的监测与评价研究［D］.北京：北京林业大学,2002.

［13］郭兆元,黄自立,冯立孝.山西土壤［M］.北京：科学出版社,1992.

［14］吴亚宁,张虎林,杨丽萍.黄河流域甘肃片土壤侵蚀遥感普查的初步分析［J］.遥感技术与应用,2002,17（6）：398-401.

［15］何兴元,胡志斌,李月辉,等.GIS支持下岷江上游土壤侵蚀动态研究［J］.应用生态学报,2005,16（12）：2271-2278.

［16］李扬,刘霞,姚孝友,等.沂蒙山区典型小流域土壤侵蚀与环境影响因子耦合关系［J］.水土保持研究,2011,18（2）：66-70.

［17］Wischmeier W H,Smith D D.Predicting rainfall erosion losses from cropland east of the Rocky Mountains：A guide for selection of practices for soil and water conservation［J］.Depatment of Agriculture,1965,21（3）：282.

［18］李婷,张世熔,廖明辉.基于3S和USLE的沱江流域中游土壤侵蚀定量评价［J］.水土保持研究,2011,18（1）：24-27.

［19］史志刚.安徽省土壤侵蚀区域分布研究［J］.灌溉排水学报,2012,31（6）：64-66.

［20］齐述华,蒋梅鑫,于秀波.基于遥感和ULSE模型评价1995-2005年江西土壤侵蚀［J］.中国环境科学,2011,31（7）：1197-1203.

［21］傅世锋,查轩.基于GIS和USLE的东圳库区土壤侵蚀量预测研究［J］.地球信息科学,2008,10（3）：390-395.

［22］周为峰,吴炳方.基于遥感和GIS的密云水库上游土壤侵蚀定量估算［J］.农业工程学报,2005,21（10）：46-50.

［23］孙立达,孙保平.西吉县黄土丘陵沟壑区小流域土壤流失量预报方程［J］.自然资源学报,1988,3（2）：141-153.

［24］孙保平,赵廷宁,齐实.在西吉县黄土丘陵沟壑区的应用［C］.中国科学院水利部西北水土保持研究所集刊.1990,12：50-58.

［25］吕喜玺,沈荣明.土壤可蚀性因子K值的初步研究［J］.水土保持学报,1992,6（1）：63-70

［26］B.Y.Liu,M.A.Nearing,P.J.Shi,Z.W.Jia.Slope Length Effects on Soil Loss for Steep Slopes［J］.Soil Science of America Journal,2000（64）：1759-1763

［27］Wischmeier W H,Smith D D.Predicting rainfall erosion losses.A guide to conservation planning.U.S. ［J］.Department of Agriculture,1987,20（6）：537.

［28］Wischmeier W H,Smith D D.Predicting rainfall erosion losses from cropland east of the Rocky Mountains：A guide for selection of practices for soil and water conservation［J］.Department of Agriculture,1965,21（3）：282.

［29］江忠善,王志强,刘志.黄土丘陵区小流域土壤侵蚀空间变化定量研究［J］.土壤侵蚀与水土保持学报,1996,1（2）：1-9.

［30］安培浚,刘树林,领耀文,等.植被指数遥感定量研究—以民勤绿洲为例［J］.遥感技术与应用，2005,20（6）：576-581.

［31］刘得俊,李润杰,王文卿,等.基于地理信息系统的西宁市土壤侵蚀监测的实现［J］.水土保持研究,2006,13（5）：111-113.

［32］傅世锋,查轩.基于GIS和USLE的东圳库区土壤侵蚀量预测研究［J］.地球信息科学,2008,10（3）：390-395.

［33］张岩,刘宪春,李智广,等.利用侵蚀模型普查黄土高原土壤侵蚀状况［J］.农业工程学报,2012,28（10）：165-171.