基于GIS和CSLE的高淳慢城土壤侵蚀评估

邱娅柳，耿 韧，洪静雨，吴七斤，洪大林*

（1.南京市高淳区水务局，江苏 南京 211300；2.南京水利科学研究院水文水资源与水利工程国家重点实验室，江苏 南京 210029；3.南京瑞迪建设科技有限公司，江苏 南京 210009）

土壤侵蚀通过剥蚀、搬运土壤，切割地表，破坏地面完整性，造成土地退化，危及农业生产，侵蚀泥沙淤积河湖库塘，加剧洪涝、干旱等灾害的发生，同时土壤侵蚀物携带营养物质造成水体污染，严重威胁人类的生存和发展，是全球面临的主要环境问题之一[1]。定量计算土壤侵蚀量可为合理利用和管理土地资源提供科学依据[2]，土壤侵蚀模型是定量计算土壤侵蚀量的重要方法之一[3-5]。通用土壤侵蚀模型（USLE）是世界上应用最为广泛的经验模型[3]，近些年随着地理信息系统（GIS）技术的快速发展，GIS结合ULSE已经广泛应用于土壤侵蚀量的预报[1，6]。

ULSE是基于缓坡和年降雨的土壤侵蚀模型，其外延性不是甚好[7]。ULSE在我国的应用过程中，应注意降雨侵蚀力指标、单位转换、标准小区设计、坡长坡度因子计算等问题，否则会导致结果不可比[8]。基于以上建模思想，刘宝元结合我国水土流失的实际情况，提出水土保持措施的三分法分类，建立了适用于我国陡坡农地的陡坡公式，采用不同精度降水资料估算降雨侵蚀力的计算方法，计算了我国主要土壤可蚀性因子值，建立了中国土壤流失方程（CSLE），并应用于第一次全国水利普查水土保持情况普查[9]。本研究选择南京市高淳慢城地区为研究对象，采用CSLE与GIS结合的方法，研究该地区土壤侵蚀量的时间和空间变化，以期精确评价区域水土流失时空变化，为区域水土流失科学治理提供技术指导。

1 材料与方法

1.1 研究区概况

高淳慢城总面积81.60 km2，主要位于南京市高淳区境内（图1）。慢城区域为低山丘陵区，东部为茅山山脉的余脉环绕，青山连绵，山势平缓，以海拔150 m以内低山丘陵为主，呈东北西南走向延伸，是太湖水系与水阳江、青弋江水系的分水岭。局部山体如大山、荆山等坡度稍陡，但仍然属于30°以内缓坡。研究区位于北亚热带和中亚热带过渡地区，受季风环流影响，区域性气候明显，常年四季分明，平均气温16 ℃，最高气温39.8 ℃，最低气温-14.0 ℃。多年平均降水量1190.8 mm，降水量年内分配不均，年际变化较大，最大年降水量1878.6 mm（1991年），最小年降水量569.5 mm（1978年）。土壤以水稻土类为主，主要包括青泥土、白泥土、黄泥土、马肝土、泥骨土，旱地土壤为沙土、黄土、夜潮土等。区域内农用地主要由水田、旱地、水浇地、田坎组成，其中水田面积最大，占51.95%。

图1 研究区地理位置示意图

1.2 预报模型

采用中国土壤流失方程CSLE（Chinese Soil Loss Equation）评价典型水土流失现状和工程建设导致的水土流失变化[6]。该方程为第一次全国水利普查水土保持情况所使用的土壤侵蚀计算方程，是将影响水土流失的7个因子用连乘的形式组成，侵蚀模数（单位面积年均土壤侵蚀量）计算式为：

式中，A为土壤侵蚀模数，t/（hm2·a）；R 为降雨侵蚀力因子，MJ·mm·hm-2·h-1·a-1；K 为土壤可蚀因子，t·hm2·h·MJ-1·mm-1·hm-2；L、S 分别为坡长、坡度因子，无量纲；B为植被盖度与生物措施因子，无量纲；E为工程措施因子，无量纲；T为耕作措施因子，无量纲。

1.2.1 土地利用信息获取

本研究基于不同高分辨率遥感影像资料解译和野外实地调查获取慢城地区不同时段的土地利用数据，遥感影像数据源来自Googleearth，覆盖慢城地区的有2015年、2014年和2007年等年度的0.6 m分辨率的Quickbird和1 m分辨率的IKONOS影像。为了研究工程建设等人类活动对慢城地区水土流失风险的影响，本文选择2015年的影像作为土地利用现状资料，选择2007年的影像作为建设前的土地利用资料，以2014年影像作为建设中期的土地利用资料数据源。结合野外实地调查的土地利用方式及其在高分辨率影像中的特征，采用人工目视解译的方法判读影像中的土地利用方式信息，在ArcGIS中编辑面图层shp文件，勾勒出每种土地利用类型斑块的轮廓并赋予相应的属性，最后计算每种土地利用斑块的面积。

1.2.2 降雨侵蚀力计算

本研究采用日雨量计算降雨侵蚀力，其计算公式如下：

根据2007年、2014年和2015年高淳站的逐日降水资料，由以上方程式，计算可得2007年、2014年和2015年慢城的降雨侵蚀力因子分别为2412.1MJ·mm·hm-2·h-1·a-1、7109.48MJ·mm·hm-2·h-1·a-1和 6549.55MJ·mm·hm-2·h-1·a-1。由于慢城地区的降雨量分布比较均匀，整个地区只采用一个降雨侵蚀力因子值。在慢城边界polygon文件的属性列表中添加相应字段并分别赋2007年、2014年和2015年降雨侵蚀力值。在Arcgis中生成2007年、2014年和2015年的慢城地区的降雨侵蚀力的栅格文件。

1.2.3 土壤可蚀性计算

王彬对比分析了USLE、RUSLE2、EPIC和Dg-OM模型的土壤可蚀性计算方程[10]，结果表明Dg-OM模型计算我国土壤可蚀性的效果更好。因此，本研究采用Dg-OM来计算慢城地区的土壤可蚀性，其计算公式如下：

式中，fi为土壤中粒径百分比（%）；mi为小于该粒径算术平均值（mm）；OM为土壤有机质含量（g kg-1）；Dg为平均几何粒径（mm）。

查阅我国第二次土壤普查资料，可知慢城地区小于0.002 mm土壤粒径的含量为23.3%，粒径组成0.002～0.02 mm的含量为30.3%，粒径组成0.02～0.05 mm的含量为44.6%，粒径组成0.05～2 mm的含量为1.4%，土壤有机质含量为1.69g kg-1。将该地区的土壤粒径数据和土壤有机质数据，代入式（6）可得慢城地区的土壤可蚀性因子值为 0.025726 t·hm2·h·MJ-1·mm-1·hm-2。在慢城地区polygon文件的属性列表中添加K字段，字段类型为double，赋值为0.025726，运用Arcgis进行矢量向栅格数据的转换。

1.2.4 LS因子计算

根据已获得的30 m数字高程栅格数据DEM（ASTER GDEM），来计算慢城地区坡长因子L和坡度因子S。其计算公式如下：

式中，L为坡长因子，无量纲；λ为坡长（m）；S为坡度因子，无量纲；θ为坡度值（°）。

根据以上方程，用区域LS因子计算工具（LS-TOOL）计算出慢城地区坡度和坡长因子。

1.2.5 BET因子计算

植被覆盖与生物措施因子B采用直接赋值法，根据前人研究成果并结合本研究区土地利用实际情况确定慢城地区不同土地利用的B因子值[6]。8种土地利用方式具体B值赋值结果：居民地为0；水体为0；林地为0.05；旱地为0.1943；建筑用地为0；草地为0.045[6]；水田为0.1；道路为0。由于慢城地区的水土保持工程措施和耕作措施较少，因此本研究不考虑水土保持工程措施和耕作措施对土壤侵蚀的影响，即工程措施因子和耕作措施因子统一按1赋值处理。将2007年、2014年和2015年的土地利用polygon文件中属性列表中分别添加字段，并按照上文确定的B值进行赋值，再运用Arcgis10.2生成2007年、2014年和2015年的B因子栅格文件。

1.2.6 侵蚀模型计算

得到各个地块的各因子栅格图层后，根据中国土壤流失方程（CSLE）的土壤侵蚀模数计算公式，用Arcgis10.2来计算慢城地区2007年、2014年和2015年的土壤侵蚀模数M。然后根据水利部颁布的《土壤侵蚀分类分级标准》（SL 190-2007）来确定土壤侵蚀分级指标，依据此指标用Arcgis10.2栅格再分类模块（Reclassify）生成研究区土壤侵蚀强度等级图。

2 结果与分析

2.1 土壤侵蚀的时间变化

随着慢城区域开发过程，区域土壤侵蚀模数呈现先增加后下降的变化趋势。2014年、2015年慢城区域平均土壤侵蚀模数分别为12.53 t/（km2·a）、11.54 t（/km2·a），分别为2007年（4.38 t（/km2·a））的2.86倍、2.63倍，2015年慢城区域平均土壤侵蚀模数较2014年下降了7.9%（图2）。2014年、2015年慢城区域土壤侵蚀模数极大值分别为 995.04 t/（km2·a）、916.67 t/（km2·a），分别为2007年（301.43 t/（km2·a））的3.30 倍、3.04倍，2015年慢城区域平均土壤侵蚀模数较2014年下降了7.9%（图2）。

图2 2007年、2014年和2015年慢城地区土壤侵蚀模数

总体而言，2007年慢城区域土壤侵蚀模数小于2014年和2015年土壤侵蚀模数，2015年土壤侵蚀模数较2014年小幅下降。这是由于2010年后伴随慢城区域旅游开发、交通建设等过程，区域土壤侵蚀不断加剧，使得土壤侵蚀模数大幅增加。在后期由于建设强度的降低，同时对于区域生态环境保护意识的加强，区域土壤侵蚀得到了有效控制，土壤侵蚀模数逐渐降低。此外，不同年份降雨量的差异也是导致区域土壤侵蚀模数变化的重要因素。

2007年、2014年和2015年区域土壤侵蚀皆以微度和轻度侵蚀为主（图3），2007年、2014年和2015年微度和轻度侵蚀面积百分比之和分别为97.3%、83.3%和83.3%。2007年微度和轻度侵蚀面积百分比之和分别较2014年、2015年增加14.0%和14.0%（图3）。2014年和2015年各级侵蚀强度面积百分比较为接近，无显著差异（图3）。

总体而言，2007年慢城地区土壤侵蚀强度小于2014年和2015年的土壤侵蚀强度，2014年和2015年慢城地区土壤侵蚀强度差异不大。这一方面可以解释为2007年的年降雨量较2014年和2015年小，而2014年和2015年的年降雨量差异不大；另一方面也可以解释为2007年生产建设项目较少，对植被和地表的扰动较小，因而土壤流失较少。而2014年和2015年，有大量生产建设项目的实施，造成该地区地表的扰动和植被的破坏，因而土壤流失较为严重。

图3 不同年限慢城地区侵蚀强度面积百分比对比图

表1统计了不同土地利用类型的平均土壤侵蚀模数，由于居民点、水体、建筑用地和道路不产生土壤侵蚀，因此对以上土地利用不做分析。4种土地利用方式的平均土壤侵蚀强度均随时间进程呈现先增加后降低的变化趋势。2014年和2015年慢城区域林地、旱地、草地和水田平均土壤侵蚀模数约为2007年平均土壤侵蚀模数的3倍，而2014年和2015年4种土地利用方式平均土壤侵蚀模数差异不大。从表1中可以看出，不同年份不同土地利用类型的平均土壤侵蚀模数均为旱地＞林地＞草地＞水田。

表1 不同土地利用类型土壤侵蚀

2.2 土壤侵蚀的空间变化

从图2中可以看出，不同年份土壤侵蚀模数较大值主要分布在中南部和东北部低山丘陵区，这是由于该区域地形坡度较大，土壤侵蚀较为强烈。西部龙墩河水库附近侵蚀也较为强烈，但分布面积较小，龙墩河水库附近畜禽养殖导致了植被覆盖度较低和地表扰动，因而引起较大的土壤侵蚀。西部和中东部土壤侵蚀模数较小，这些地方地形起伏较为平缓，植被覆盖度高，分布较多水体和建筑物。从2.1部分对慢城地区土壤侵蚀的时间变化的分析和图3中可以看出慢城地区以微度侵蚀和轻度侵蚀为主，与第一次全国水利普查水土保持情况普查中江苏省土壤侵蚀强度状况一致，也与林杰在南京江宁地区的研究结果一致[11]。

从表1中可以看出，旱地在4种平均土壤侵蚀模数中土地利用最高值，这与旱地较为频繁的农事活动造成表层土壤扰动有关[12-13]。慢城地区旱地较大的平均土壤侵蚀模数与许月卿在贵州猫跳河流域[14]、张雪华在东北黑土区[15]的研究结果一致。水田的平均土壤侵蚀模数是4种土地利用中最低的，旱地的平均土壤侵蚀模数为水田平均土壤侵蚀模数的3.7～3.9倍。

2.3 土壤侵蚀防治措施

由以上分析可知，旱地的平均土壤侵蚀模数是4种土地利用的最大值，应加强对旱地的治理，缓坡修建水平梯田，陡坡退耕还林还草。慢城地区中南部和东北部低山丘陵区土壤侵蚀较为强烈，亦应加强这些部位的土壤侵蚀治理工作。对于分水岭应选择耐旱、抗风、根系发达、保土能力强的树种，实行乔、灌、草结合。对于坡地，应营造坡地水保林草，并辅以相应的工程措施，如鱼鳞坑、水平梯田等。对于侵蚀沟，应结合固沟、拦沙工程，如修建谷坊（柳谷坊、柳土坊）和拦沙坝，配套营造沟底（包括沟口）防冲林。慢城地区龙墩河水库地区的土壤侵蚀也较为严重，对于这些部位土壤侵蚀的治理，应在入库径流和泥沙的主要沟道内结合谷坊，设置挂淤林。水库边岸上，林地条件较好，可营造速生优质的杉木、马尾松、刺槐、杨树等。岸边周围，地形平缓，可因地制宜开发经济林。

3 结论

（1）随着慢城区域开发过程，区域土壤侵蚀模数呈现先增加后下降的变化趋势。慢城地区2007年土壤侵蚀模数明显小于2014年和2015年，2015年土壤侵蚀模数较2014年小幅下降。

（2）慢城地区土壤侵蚀以微度侵蚀和轻度侵蚀为主，不同年份土壤侵蚀模数较大值分主要分布在中南部和东北部低山丘陵区，西部龙墩河水库附近侵蚀也较为强烈。旱地平均土壤侵蚀模数为4种土地利用的最高值。

（3）高淳慢城区域水土保持工作重点为加强对旱地的治理、缓坡修建水平梯田、陡坡退耕还林还草，慢城地区中南部、东北部低山丘陵区和龙墩河水库地区为区域水土流失治理的重点区域。