基于GIS和USLE的土壤侵蚀定量分析研究——以四川省洪雅县为例

周湘山，孙保平，李锦荣，赵 岩，钟晓娟，王引乾，冯磊，邱一丹

(北京林业大学水土保持学院水土保持与荒漠化防治教育部重点实验室，北京100083)

土壤侵蚀已成为全球性的环境问题之一，引起了世界各国的广泛关注。我国是世界上土壤侵蚀问题最严重的国家之一，其范围遍及全国各地，主要发生在黄河中上游黄土高原地区、东北黑土地区和长江中上游丘陵地区。土壤侵蚀是土地退化、河道和湖泊淤积的主要原因，也是导致生态环境恶化的重要因素。对区域土壤侵蚀类型和过程研究是一项基础性的工作。然而，传统的土壤侵蚀研究方法难以进行大尺度区域的土壤侵蚀状况的调查与分析。20世纪80年代以来，GIS开始被引入与通用土壤流失方程(USLE)结合进行流域土壤侵蚀量的预测和估算。90年代初至今，二者的集成被广泛地应用于土壤侵蚀的风险评价、关键源区的识别、不同水土保持方案的模拟评价以及流域的管理规划等，有效地提高了土壤侵蚀量的预测效率和结果的显示度[1]。USLE作为使用最广泛的土壤侵蚀模型，已经发展成为一个相对标准化的方法，在我国得到广泛的推广应用。例如，蔡崇法等[2]依据实地调查资料，建立了典型小流域地理数据库，并在地理信息系统IDRISI支持下，采用USLE土壤侵蚀预测模型预测了小流域土壤侵蚀量;洪华生等人[3]基于GIS和 USLE对下庄小流域开展土壤侵蚀量的预测研究，实现了土壤侵蚀量估算和土壤侵蚀关键源区的标识。杨娟[4]等人基于GIS和USLE的卧龙地区小流域土壤侵蚀预报，分析小流域的土壤侵蚀强度的空间分布格局，以及土地利用/覆盖对流域土壤侵蚀的影响。基于此，以四川省洪雅县为研究区域，选取合适的因子算法，定量确定了土壤侵蚀因子指标。在ArcGIS的支持下，建立空间数据库，根据土壤侵蚀预测模型对数据库进行图形运算，实现了对洪雅县土壤侵蚀强度的分等定级，并根据洪雅县的土地利用情况与经济发展规划模拟了水土保持的治理措施，为开展县域水土保持规划工作提供科学依据。

1 研究区概况

洪雅县位于四川盆地西南边缘，地理位置跨东经102°49′-103°32′，北纬 29°24′-30°00′。斜置于峨眉山、大相岭和总岗山之间，形如哑铃。全县幅员面积1 896.49 km2，其中耕地 27 987 hm2，林地 139 117.27 hm2，牧草地4 196.08 hm2。洪雅县位于中亚热带湿润季风气候北缘的山地平原交错带，地理环境复杂，生态条件多样。气候、土地、生物和水等农业自然资源丰富，年均温 16.8℃，无霜期 304 d，≥0℃积温6 149.2℃，≥10℃积温 5 309.6℃，多年平均降水量1 493.8 mm。洪雅地处四川盆地和凉山山原的过渡地带，具有平坝、丘陵和山地等地貌类型，土地资源以林地为主(占总面积的68.8%，人均5.05 hm2)。该县农耕地以水稻土、紫色土为主，林地以山地黄棕壤为主。

2 数据来源和研究方法

2.1 数据来源

本研究采用的图形数据库包括:土地利用类型图、土壤类型图、地形图、流域数字高程模型(DEM)等。其中土地利用图根据2009年的Landsat/TM影像，TM影像为 LANDSAT7波段遥感数据，其地面分辨率为30 m×30 m，根据不同土地利用类型的光谱反映特征建立解译标志，采用目视解译法识别影像的特征属性，并结合野外资料对影像进行监督分类，得到遥感分类图。土壤数据来自全国第二次土壤普查图件及各地土壤志。其精度为1∶100万，包括土壤类型矢量数据以及黏粒(%)、粉砂(%)、砂粒(%)、有机质(%)等数据。属性数据包括:洪雅县1998-2007年的气象数据。

运用ArcGIS的空间数据管理和分析功能，建立流域的数字高程模型(DEM)、土地利用现状图、土壤类型图等矢量图，对其属性数据进行相应的数据编码操作，再借助ArcGIS将上述各图栅格化(Grid格式)，考虑到县域尺度和计算精度，栅格大小为10 m×10 m，获得USLE方程的因子图，再通过ArcGIS将各因子图连乘，获得洪雅县的土壤侵蚀强度等级图。

2.2 USLE模型

通用土壤流失方程(USLE)由Wischemier和Smith[5]于1958年提出，为提高各因子计算的通用性，美国农业部又提出了修正方程RUSLE。本研究采用这一模型，其数学表达式为:

式中:A——单位面积上年平均土壤流失量;R——降雨侵蚀力因子;K——土壤可蚀性因子;L——坡长因子;S——坡度因子;C——植被覆盖和管理因子;P——土壤保持措施因子。

2.2.1 R值计算 降雨侵蚀力因子是评价降雨引起土壤分离和搬运的潜在动力指标，虽然一些基于小区实验的回归方程给出了降水因子组合算法，但是这些相对精准算法的数据往往难以获得，所以国内外的很多研究人员都偏向于采用一些降水侵蚀因子R的简易算法，其中应用比较多的就是Wischmeier提出的直接利用多年各月平均降量的经验公式。Wischmeier经验公式，既考虑了年总降雨量，又考虑了降雨的年内分布，比较准确地反映出区域降雨对土壤侵蚀的贡献率。本研究R值的计算采用Wischmeier经验公式[6]:

式中:Pi——第i月的平均降雨量(mm);P——全年或多年平均降雨量(mm);i=1，2，…，12。Pi和P 来自洪雅县气象站的降雨资料。

2.2.2 K值计算 依据Ei-Swaify等的研究，K值的大小与土壤质地有较高的相关性。目前土壤可蚀因子K值的求法主要有:直接测定法、诺谟图法、公式法[7]。直接测定K值法被认为最符合田间实际土壤对侵蚀力的敏感尺度，但是直接测定K值所需的时间较长，经费较多;诺谟图法不仅需要较多参数，特别是土壤结构级别和土壤渗透级别很难准确获得;公式法则比较快捷，也较为准确。美国通用土壤流失方程USLE及其修正模型RUSLE中K值的确定方法，与诺谟图法一样，需要的指标不容易确定。因此，本研究，运用Williams等人在EPIC模型中发展的土壤可蚀性因子K的估算方法进行计算[8](图1)。

式中:Sa——砂粒(0.1～2 mm)含量(%);Si——粉粒(0.002～0.1 mm)含量(%);Ci——黏粒(＜0.002 mm)含量(%);C——有机碳含量(%);Sn=1-Sa/100。

图1 洪雅县土壤可蚀性K值

2.2.3 LS因子计算 在小流域和区域尺度上，坡度和坡长很难通过实地测量获得，因此，长期以来缺少坡度和坡长等地形指标的估算方法。国内外研究学者对LS因子计算原理和关键技术的研究仍显薄弱。目前运用较为广泛的就是通过数字高程模型(Digital Elevation Model DEM)对LS因子进行提取。

(1)L坡长因子。运用ArcGIS建立洪雅县数字高程模型，利用 ArcGIS的 GRID模块进行基于DEM的地形特征分析，提取坡度坡长图(图2)。用通用土壤流失方程中的坡长指数L的计算方法计算坡长因子。

式中:λ——坡长(m);m——坡长效应指数，当坡度＜l%时 ，m=0.2;坡度 1%～ 3%，m=0.3;坡度 3%～5%，m=0.4;坡度＞5%，m=0.5[9]。

(2)S坡度因子的计算。缓坡上选用McCool等(1987年)研究的坡度公式;在陡坡上采用刘宝元的坡度公式[10]:

式中:S——坡度因子;θ——坡度。

2.2.4 C作物管理及植被盖度因子的确定

(1)植被盖度的提取。植被盖度的计算现在大多数是利用遥感提取植被指数，利用植被指数与外业调查的植被盖度建立相关关系，最终得到植被盖度图。而最常用的就是利用植被NDVI指数，来估算植被盖度。

图2 洪雅县坡度坡长因子LS值

根据安培浚等[11]在西北干旱区利用NDVI所建立的植被指数与植被覆盖度的关系来确定植被覆盖度c。

其中:NDVI——所求像元的归一化植被指数。NDVImin，NDVImax——研究区 NDVI的最小值和最大值。

(2)C因子的确定。C值主要受植被覆盖度和土地利用现状的制约，蔡崇法等[12]通过建立坡面产沙量与植被覆盖度的数学模型得到了C因子值与植被覆盖度之间的关系。本次研究参考Wischmeier的原始C因子取值对照表的定义及相关研究文献，结合洪雅县的实际情况，利用式(7)计算各地类年均C因子(图3)。

图3 洪雅县作物管理及植被盖度因子C值

2.2.5 P保护措施因子的确定 水土保持措施因子P是采用专门措施后土壤流失量与顺坡种植时的土壤流失量的比值，反映植被的管理措施差异引起的土壤流失量差别，其范围在0～1之间。本次研究主要参照美国农业部手册703号和洪雅县2009年TM遥感卫片解译出的土地利用现状图及相关的文献来进行P值的确定[13](图4)。

图4 洪雅县保护措施因子P值

3 结果分析

3.1 计算结果和数据分析

3.1.1 RUSLE各要素的计算结果 运用Matlab数学软件对洪雅县土壤粒径进行3次样条的插值，用苏联制为基础进行插值，转换得出美国制、国际制的主要有效粒径对应的土壤颗粒累积百分含量。再用ArcGIS对土壤属性数据处理得到洪雅县土壤可蚀性K值图;利用遥感提取植被指数，算出植被覆盖度c。用ArcGIS对植被覆盖度c进行计算得出洪雅县作物管理及植被盖度因子C值图;运用ArcGIS建立洪雅县数字高程模型，计算得出洪雅县坡度坡长因子LS值图;运用遥感采用目视解译法识别影像的特征属性并结合野外资料对影像进行监督分类，得到遥感分类图。用ArcGIS对遥感分类图进行赋值运算得到洪雅县保护措施因子P值图。考虑到县域尺度和计算精度，栅格大小为10 m×10 m，将得出的USLE方程的因子图，通过ArcGIS将各因子图连乘，获得洪雅县的土壤侵蚀强度等级分布图(图5)。

图5 洪雅县2009年土壤侵蚀强度分布

3.1.2 土壤侵蚀强度统计分析 通过图5统计得到表1，可以得出洪雅县土壤侵蚀强度以中度以下等级为主，其中中度以下侵蚀强度面积共170 115.15 hm2，占土地总面积的89.7%;中强度侵蚀面积共19 533.85 hm2，占土地总面积的10.3%。

表1 洪雅县2009年土壤侵蚀强度分级面积

3.1.3 土地利用类型统计分析 运用遥感采用目视解译法识别影像的特征属性并结合野外资料对影像进行监督分类，得到遥感分类图，运用ArcGIS对遥感分类图进行处理得出洪雅县2009年土地利用类型图(图6)。

图6 洪雅县2009年土地利用类型

表2 洪雅县2009年土地利用类型分类面积统计

通过图6和表2可以得出洪雅县土地利用类型中耕地面积 27 726.68 hm2，占土地总面积的14.62%;灌木林地面积17 523.57 hm2，占土地总面积的9.24%;有林地面积123 404.6 hm2，占土地总面积的65.07%;水域面积7 301.49 hm2，占土地总面积的3.85%;城镇村及工矿用地面积9 463.49 hm2，占土地总面积的4.99%;草地面积4 229.17 hm2，占土地总面积的2.23%。洪雅县土地利用类型主要以耕地和林地为主，耕地和林地一共168 654.85 hm2，占全县土地总面积的88.93%。

表3 洪雅县2009年土地利用类型的侵蚀强度的面积百分比 %

根据表3可以得出洪雅县各土地利用类型的侵蚀强度分布中:耕地和草地主要以中度以下侵蚀为主，其中无明显侵蚀面积分别占20.95%和32.34%，微度侵蚀面积分别占21.59%和24.53%，轻度侵蚀面积分别占48.79%和36.26%，中强度侵蚀面积分别共占8.67%和6.87%。灌木林地和有林地中无明显侵蚀均占有较大比例，分别为50.67%和57.35%，其次是微度和轻度侵蚀面积，分别共占有41.29%和36.65%，中强度侵蚀面积最小，分别共占8.04%和6%。水域和城镇村及工矿用地中度侵蚀和轻度侵蚀占有较大面积，分别共占84.01%和82.31%，相比其他土地利用类型，水域和城镇村及工矿用地中的强度侵蚀面积占有比例相对过高，分别为 4.39%和9.81%，微度侵蚀和无明显侵蚀面积占有比例显著减少，分别共占11.6%和7.88%。洪雅县土地利用单元的侵蚀强度存在显著差异，不同的土地利用类型跟土壤侵蚀强度高低有着较为密切的联系。

3.2 土壤侵蚀强度分区

结合图5的洪雅县2009年土壤侵蚀强度分布情况，图6的洪雅县2009年土地利用类型状况和洪雅县乡镇边界确定洪雅县土壤侵蚀强度三区分区边界，得到洪雅县土壤侵蚀强度分区分布图(图7)。

图7 洪雅县土壤侵蚀强度分区分布

表4 各区土壤侵蚀强度分级面积及平均侵蚀模数

根据洪雅县的土壤侵蚀强度分布情况、土地利用类型分布情况和乡镇分布情况，大致将可以将洪雅县分为以下3个区域(图7)。

(1)北部轻度侵蚀区。本区包括汉王、中山、洪川 、余坪 、中保 、止戈 、将军 、三宝 、东岳等乡镇 ，面积约665.54 km2，占全县的 35.09%，人口 24.8万人，占全县总人口的74.97%。该区土地利用类型以耕地、水域、城镇村及工矿用地为主。该区地貌以沿江河谷冲积平原、台地、丘陵为主，水、热、气条件好，人口密度大，人均耕地少，能源、交通条件较好，又是全县农、工综合开发区，经济最发达。该区平均土壤侵蚀模数为2 118.20 t/(km2◦a)，土壤侵蚀强度以轻度为主。

(2)中部微度侵蚀区。本区包括花溪、柳江、高庙、桃源等乡镇，面积约 491.69 km2，占全县的25.93%。人口5.2万人，占全县总人口 15.72%。该区土地利用类型以耕地、林地和水域为主。该区地貌以中山为主兼有低山，气候呈立体分布人口稀少，山地宜林宜牧面积大，饲草充足，矿物贮量较丰富，是全县林、药、土特产的主要生产基地。该区平均土壤侵蚀模数为1 055.30 t/(km2◦a)，土壤侵蚀强度以微度为主。

(3)西南部无明显侵蚀区。本区包括瓦屋山、张村等乡镇，全区幅员面积约739.25 km2，占全县总面积的 38.98%。人口 3.08万人，占全县总人口9.31%。该区土地利用类型主要包括:林地、草地和水域。区内地处县境西南部，绝大多数为国有林地，有国家级森林公园瓦屋山，是人烟稀少的天然林区。由于山峰高耸，河谷深切，气候寒冷，人迹罕至，所以本区生态环境尚属天然或原始状态。加之位于青衣江流域上游，对水源涵养和生态平衡具有重要作用。该区平均土壤侵蚀模数为498.82 t/(km2◦a)，土壤侵蚀强度以无明显侵蚀为主，瓦屋山水库等局部地区伴有中强轻度侵蚀。

3.3 土壤侵蚀强度分区治理措施

根据洪雅县的土壤侵蚀强度分布情况、土地利用类型分布情况和乡镇分布情况， 可以知道洪雅县土壤侵蚀强度三区各区的生产方式和社会经济功能， 从而界定三区的水土流失主要方式和原因，做出针对性和准确性的分区治理措施。

(1)北部轻度侵蚀区。该区域包括洪雅县城区及其它7个主要乡镇，地势较为平坦，人口密集，交通、能源条件好，是全县的农工综合开发区和粮食作物耕作区。造成该区域水土流失的主要因素包括城镇、工矿等项目的开发建设，不合理的耕作方式和沿江流域的水力侵蚀。

该区应加强城市水土流失预防监督工作，认真贯彻水土保持相关法律法规，督促开发建设项目业主认真履行水土流失防治义务，按相关法律法规明确水土流失防治责任，确定水土流失防治方案和措施。工业园区、商业贸易区、经济开发区、文化娱乐区、生活居住区等，根据不同的开发建设目标，采取与之相适应的水土保持措施，同时城区的水土保持措施与城市水系布局、河道整治、排洪工程结合，提高水系绿化指数和城市雨洪调蓄能力，减轻因水土流失而产生的洪涝灾害和城市水源污染[14]。在人口相对密集的平坝、丘陵等农业生产区域，结合农业结构调整，加强控制化肥农药的使用，发展与水源保护相适应的生态农业、观光农业、休闲农业，积极发展城郊型农业服务基地，减少面源污染。加强小型水利水保设施建设，改善生产条件，因地制宜在村庄及旅游景点等人类活动和聚居区加强农村污水垃圾处理，改善人居环境。

(2)中部微度侵蚀区。该区地貌以中山为主兼有低山，土地利用类型包括耕地、林地和水域。其中耕地以坡耕地为主，坡耕地改造工程已成为水土保持生态建设的重中之重。水土保持法规定25°以上的坡耕地退耕还林还草，25°以下的坡耕地要采取相应的水土保持措施方可耕作。坡耕地治理措施的主要通过改变坡面小地形的方法防止坡地水土流失，将雨水及融雪水就地拦蓄，使其渗入地下，减少或防止形成坡面径流，保持土壤养分，增加农作物可利用的土壤水分，改善农业生产基础条件，提高土地生产能力。天然沟壑及其两侧山坡地形成的小流域，需因地制宜地采取工程措施、植物措施和保护性耕作等措施，进行坡耕地和沟道水土流失综合治理。河道两侧的水土流失治理应保育植被，恢复湿地，清理河道垃圾、障碍物，恢复景观生态，有效发挥灌木和水生植物的水质净化功能，维护河道及湖库周边生态平衡。

(3)西南部无明显侵蚀区。该区属青衣江流域上游，是四川省西南山地重要的生物基因库，也是国家宝贵的生物种源基地。区内有瓦屋山国家级森林公园和瓦屋山水库，生态环境尚属天然或原始状态，对水源涵养和生态平衡具有重要作用。

该区应加强预防监督管护，封山禁牧、禁伐、禁垦，保天然原始森林存量和保护次生林恢复森林植被，充分发挥植被特别是灌草植被的生态功能，在流域上游形成一道生态防线，以达到涵养水源、保持水土的目的。加强小型水库等水资源保护与开发建设，为下游生产和生活蓄积和提供水源。由于该区地势海拔较高，容易形成泥石流、滑坡等灾害，当地应采取监测、径流排导、削坡减载、支挡固坡、修建拦挡工程等措施，建立监测、预报、预警体系。

4 结语

研究利用ArcGIS强大的空间数据管理能力，进行USLE方程各因子参数的输入和运算，实现对洪雅县土壤侵蚀强度的初步定量分析，给水土保持规划治理的提供了较为可靠的依据。相对于传统的集中式方法进行土壤侵蚀强度预测，这种分布式方法由于运用GIS的栅格数据分析功能，可预测出每个栅格的土壤侵蚀量，便于管理者找出较为严重的土壤侵蚀区，并针对性提出最佳治理措施。不过目前的预测模型和侵蚀因子的算式算法大多是借鉴原有的技术和国内外相关的方法与成果来确定的，在以后的工作中，还应根据研究区水土流失实测结果，进行校准、改进各因子的算式算法，深入了解区域土壤、地形、覆被等影响因素的变化，加强对未来各侵蚀因子变化的定量分析，以提高预测的可靠性、准确性和突出区域特色[15]。

[1] 耿豪鹏，潘保田，王超，等.基于GIS与 USLE的榆中县土壤侵蚀[J]兰州大学学报:自然科学版，2009，45(6):9-13.

[2] 蔡崇法，丁树文，史志华，等.应用 USLE模型与地理信息系统IDRISI预测小流域土壤侵蚀量的研究[J].水土保持学报，2000，14(2):19-24.

[3] 洪华生，杨远，黄金良.基于GIS和USLE的下庄小流域土壤侵蚀量预测研究[J].厦门大学学报:自然科学版，2005，44(5):676-679.

[4] 杨娟，葛剑平，李庆斌.基于GIS和USLE的卧龙地区小流域土壤侵蚀预报[J].清华大学学报:自然科学版，2006，46(9):1526-1529.

[5] Wischmeier W H，Smith D D.Agricultural Handbook No.537[S].Science and Education Administration，U-nited States Department of Agriculture，1978.

[6] Wischmeier W H.A soil erodibility nomograph for farm land and construction sites[J].Journal of Soil and Water Conservation ，1971，26:189-193.

[7] 张金池，李海东，林杰，等.基于小流域尺度的土壤可蚀性K值空间变异[J].生态学报，2008，28(5):2019-2206.

[8] 吕喜玺，沈荣明.土壤可蚀性因子K值的初步研究[J].水土保持学报，1992，9(1):63-70.

[9] 蔡永明，张科利，李双才.不同粒径制间土壤质地资料的转换问题研究[J].土壤学报，2003，40(4):511-517.

[10] 刘宝元，史培军.WEPP水蚀预报流域模型[J].水土保持通 报 ，1998 ，8(5):7-13.

[11] 安培浚，刘树林，颉耀文，等.植被指数遥感定量研究:以民勤绿洲为例[J].遥感技术与应用，2005，20(6):574-580.

[12] 江忠善，王志强，刘志.黄土丘陵区小流域土壤侵蚀空间变化定量研究[J].土壤侵蚀与水土保持学报，1996，2(1):1-9.

[13] 游松财，李文卿.GIS支持下的土壤侵蚀量估算[J].自然资源学报，1999，14(1):62-68.

[14] 张永新.吉林省水土流失现状分析及防治对策[D].长春:吉林大学，2004.

[15] 徐清艳，周跃.基于GIS和 USLE的土壤侵蚀预测研究:以云南新平大红山铁矿为例[J].水土保持研究，2009，16(4):72-75.