基于CSLE和高分辨率航空影像的孤山川流域土壤侵蚀定量评价

王 凯，夏燕秋，马金辉，屈 创

（兰州大学 资源环境学院，兰州730000）

土壤侵蚀，是指地球表面的土壤及其母质受水力、风力、冻融、重力等外力的作用，在各种自然因素和人为因素的影响下，发生破坏、分离、搬迁和沉积的现象［1］。流域内土壤侵蚀可导致土地退化、土壤肥力下降和破坏土地资源等，造成洪涝和泥石流等灾害，引起生态环境恶化。黄河中游多沙粗沙区是我国水土流失最为严重的区域和黄河下游淤积泥沙的主要来源区，为国家级水土流失重点治理区，是黄河流域水土保持生态建设的重中之重，也是水土保持监测的重点区域，对该流域进行定量评价是非常必要的。因而近年来众多学者对该区域开展了许多研究：马丽梅等［2］利用黄河中游主要水文站输沙量及降雨数据，分析了黄河中游输沙与减沙特征，以及降雨与人类活动对输沙量变化的影响；高亚军等［3］探讨总结了研究区粗颗粒泥沙、平均粒径以及中数粒径的变化规律；康玲玲等［4］应用三川河流域所对应的不同降水条件、不同措施质量等级的水土保持蓄水拦沙指标，计算分析了流域的水土保持蓄水拦沙效益。

土壤侵蚀模型的研究一直是一个热门领域。从20世纪80年代至今，随着对土壤侵蚀物理过程研究的不断深入，借助GIS与RS等建立了一大批土壤侵蚀模型，如RUSLE，WEPPE，UROSEM，LISEM 等。目前，我国比较常用的方法是依据水利部颁布的土壤侵蚀分类分级标准［5］进行评价，该方法简单实用、可操作性较强，在实际工作中给予广泛应用。但也存在一些问题，主要包括：不能直接反映气候、土壤的影响，不能全面反映水保措施的影响。2000年，刘宝元等通过分析已有的大量观测小区规格，提出我国标准小区规格，并对分布在全国不同地区小区观测资料的综合分析，提出了适用于中国土壤侵蚀特征的中国土壤 流 失 模 型 方 程 （Chinese Soil Loss Equation，CSLE）［6－7］。该模型的主要优点是根据我国水土保持的实际情况，充分考虑了生物措施、工程措施和耕作措施对土壤侵蚀和水土流失过程与结果的影响，同时考虑到我国地形地貌和美国的差异，对地形因子的算法，尤其是陡坡地形因子的改进，使过程描述和结果表达更能充分反映中国实际。程琳等［8］利用CSLE模型和TM等数据对陕西省土壤侵蚀进行了定量评价分析；孙泽祥等［9］基于GIS和CSLE，利用SPOT等数据对淮河流域沂源县土壤侵蚀监测进行了调查研究。大量研究表明，利用CSLE模型进行土壤侵蚀评价的结果是比较可靠的。

本文选择CSLE模型的另一个优势在于利用高分辨率航空影像可以克服传统方式下（采用中低分辨率遥感影像）依靠人工统计粗略获得水土保持措施数量的缺点，保证数据的时效性和准确性，提高水土保持监测精度，满足水土流失治理及决策需要。

1 研究区概况

孤山川流域位于东经 110°34′—111°04′，北纬39°02′—39°26′，属于黄河干流水系，涉及陕西省府谷县、内蒙古自治区准格尔旗12个乡（镇），其总面积1 272km2，海拔811.3～1 380m。根据黄土高原综合治理分区，孤山川流域为陕蒙沙盖黄土丘陵区，降水集中且多暴雨，77.8%以上的降雨集中在汛期6—9月，最大年降雨量为最小年降雨量的4.3倍。从年际变化分析，该流域属于典型的纬向型暴雨，雨区遍布全流域，暴雨中心雨量为210mm，全流域平均雨量为144mm，平均雨强13.1mm／h，常常山洪暴发，农田摧毁，水土流失严重。

孤山川流域主要土壤类型有黄土性土壤和风沙土两大类，其次还有轻垆土、黑垆土、淤土、潮土、紫色土等。黄土性土壤结构疏松，抗蚀力弱，透水性良好，易形成地表径流，另外流域内植被稀少，林草覆盖率平均在20%左右，水土保持效益不佳。

2 研究方法

本研究采用的是刘宝元参考USLE的有关思想针对中国实际提出适用于全国范围的中国土壤流失方程 CSLE（Chinese Soil Loss Equation）模型［10］，模型公式如下：

式中：A——年均单位面积的土壤侵蚀模数［t／（hm2·a）］；R——降雨侵蚀力因子［（MJ·mm）／（hm2·h·a）］；K——土壤可蚀性因子［（t·h）／（MJ·mm）］；L——坡长因子；S——坡度因子；B——生物措施因子；E——工程措施因子；T——耕作措施因子；其中L，S，B，E，T因子无量纲。

模型所采用的数据及数据来源：① 孤山川流域周边200km范围内的13个气象站点，每个站点包含1952—2013年（62a）期间的年、月降雨数据。②2012年9月拍摄的0.38m高分辨率航空影像和5 m分辨率的DEM数据（均来源于黄河中游多沙粗沙区孤山川等重点支流水土保持监测项目）。③ 联合国粮农组织（FAO）、国际应用系统分析研究所（ⅡASA）、荷兰ISRIC－World Soil Information、中国科学院南京土壤研究所（ISSCAS）、欧洲委员会联合研究中心（JRC）于2009年3月共同发布的1∶500万比例尺的1km格网世界土壤数据库（简称HWSD）土壤类型数据。依据DEM分辨率大小，将上述各因子均重采样成5m栅格图进行叠加分析得到土壤侵蚀空间分布结果图（5m分辨率）。

2.1 降雨侵蚀力因子R

降雨侵蚀力因子反映了由降雨引起土壤分离和搬运的动力。比较经典的算法是采用降雨动能和最大30min雨强的乘积来度量降雨侵蚀力，但是由于长时间序列的降雨过程资料难以获得，该算法难以实施。因而本研究中采用Wischmeier提出的基于多年的月平均降雨量和多年的年平均降雨量经验公式［11］：

式中：pm——多年的月平均降雨量（mm）；pmi——多年的年均降雨量（mm）。

本研究中对河曲、五寨、兴县、伊金霍洛旗、榆林等13个气象站62a（1952—2013年）的降雨数据，利用上述公式计算出各个站点的降雨侵蚀力因子值（见表1），然后考虑地形因素采用反距离内插方法计算出整个孤山川流域的降雨侵蚀力因子，生成流域降雨侵蚀力R值图。

表1 孤山川流域周边13站点降雨侵蚀力因子（MJ·mm）／（hm2·h·a）

2.2 土壤可蚀性因子K

土壤可蚀性因子K是一项评估土壤遭受侵蚀敏感程度的指标。以世界土壤数据库（简称HWSD）中的孤山川流域土壤类型图为数据基础，并对该土壤数据中的砂粒、粉粒、黏粒和有机碳含量进行野外实地抽样调查与验证。运用Williams等在EPIC模型中发展的土壤可蚀性因子K值估算方法［12］，得到孤山川流域不同类型土壤的可蚀性K值，生成流域土壤可蚀性K值图（图3b）。K值计算模型公式：

式中：SAN，SIL，CLA，C——砂粒、粉粒、黏粒和有机碳含量（%），其中SN1＝1－SAN／100。

K值单位为美制单位，（t·acre·h）／（100·acre·ft·tonf·in），但国际上通常用公制单位（t·hm2·h）／（hm2·MJ·mm），两者的换算关系：美制×0.1317＝公制单位。根据美国土壤类型数据分类，研究区内主要包括的土壤类型及其成分含量如表2所示。

表2 孤山川流域各土壤类型的SAN，SIL，CLA和C含量

2.3 地形因子（LS）

地形因子是指坡度因子（S）与坡长因子（L）。本文在对于坡度因子（S）的提取中，当坡度θ小于5°时，采用 McCool等［13］提出的坡度因子公式（式4），当坡度θ≥5°时，采用的是刘宝元［14］在黄土高原建立的坡度坡长因子计算方法（式5—6），对于坡长因子（L）的计算，采用的是 W.H.维希迈尔和D.D.史密斯［11］提出的经验公式（式7）。

式中：θ——由孤山川流域分辨率为5m的数字高程模型（DEM）提取坡度图获得；γ——坡长；m——可变的坡长指数，当θ＜0.57°时，取0.2；0.57≤θ＜1.72°时，取0.3；1.72°≤θ＜5.14°时，取0.4；θ≥5.14°时，取0.5［15］。然后根据公式（4—7）分别计算出坡度和坡长因子，将二者相乘便得到LS值分布图。

2.4 水保措施因子（BET）

CSLE模型将水保措施因子区分为生物措施因子（B）、工程措施因子（E）以及耕作措施因子（T）。

（1）生物措施因子（B）。利用高分辨率（0.38m）航空影像解译的土地利用类型图和基于归一化植被指数NDVI反演出的植被盖度图，参考谢红霞等［16］在延河流域的研究，赋予不同土地利用类型和不同植被盖度下的B值（表3）。根据B值，结合土地利用图生成生物措施因子。

（2）工程措施因子（E）。E因子计算也采用谢红霞等在延河流域研究的算法公式：

式中：St——梯田面积；Sd——淤地坝控制面积；S——土地面积；α，β——梯田和淤地坝的减沙系数。由航空影像解译出的结果中，梯田面积St为17.95 km2，淤地坝的控制面积Sd为6.55km2，流域总面积S为1 272km2。α，β分别为0.763和1.0。经过计算得出孤山川流域的工程措施因子E为0.984。

（3）耕作措施因子（T）。根据高分辨率航空影像解译出的水土保持措施，结合前人研究结果［17－22］，对水土保持措施的耕作因子赋值，没有任何耕作措施的类型赋值为1。

经过地图代数计算方法，将上述生物措施因子、工程措施因子和耕作措施因子相乘，便可得到水土保持措施因子BET图。

表3 孤山川流域不同土地利用类型不同植被覆盖生物措施因子

3 结果与分析

3.1 土地利用和植被覆盖现状

通过航空影像，结合土地利用类型和外业实地调绘成果，将研究区土地利用现状解译为8个一级类（耕地、园地、林地、草地、居民点及工矿交通用地、水域及水利设施用地、未利用地和其他土地），27个二级类，33个三级类（表4），同时生成植被覆盖度现状图。

表4 不同水保措施耕作因子值

3.2 计算结果

根据上述研究方法，分别计算出孤山川流域各因子值的栅格图（图1）。

利用上述计算的降雨侵蚀力因子、土壤可蚀性因子、地形因子、水土保持措施因子结果和CSLE模型，运用地图代数得到孤山川流域2012年土壤侵蚀图（图2）。

孤山川流域侵蚀类型主要是水力侵蚀，该流域2012年侵蚀总量达274.89万t，年均侵蚀模数为2 161.06t／（km2·a）。从图中可以看出，该流域南部较北部侵蚀严重。依据土壤侵蚀分类分级标准对该流域进行土壤侵蚀强度分级，其中微度侵蚀年侵蚀总量27.86万t，年均侵蚀模数为343.84t／（km2·a）；轻度侵蚀年侵蚀总量36.55万t，年均侵蚀模数为1 574.41 t／（km2·a）；中度侵蚀年侵蚀总量37.91万t，年均侵蚀模数为3 531.11t／（km2·a）；强烈侵蚀年侵蚀总量30.64万t，年均侵蚀模数为6 290.29t／（km2·a）；极强烈侵蚀年侵蚀总量44.15万t，年均侵蚀模数为10 726.43t／（km2·a）；剧烈侵蚀年侵蚀总量97.78万t，年均侵蚀模数为30 206.96t／（km2·a）。整体而言，流域以微度侵蚀和轻度侵蚀为主（表5）。

图1 小流域R，K，LS，BET因子值图

图2 孤山川流域2012年土壤侵蚀现状

3.3 土壤侵蚀结果分析

（1）不同坡度的土壤侵蚀分析。对该流域坡度栅格图按照坡度分级标准进行分级，再与土壤侵蚀结果进行叠置分析，得出每个坡度级别的土壤侵蚀情况。结果表明，坡度在0～25°范围带内年均侵蚀模数小于该流域年均侵蚀模数，这是由于25°以下属缓坡区域，且经过退耕还林工程后主要为人工林地和保留的基本农田，植被在缓坡易于生长，林草覆盖度较高，保留的基本农田大多进行了梯田改造，种植作物通常采用等高耕作，因此，缓坡区域侵蚀强度较弱。随着坡度的增加，年均侵蚀模数增加；在坡度大于25°坡度带范围内的年均侵蚀模数达3 766.01 t／（km2·a），其土壤侵蚀总量占流域总侵蚀量的70.70%（表6）。所以陡坡区域应是今后该流域水土流失的重点治理对象，着重对坡度大于25°的未利用地等土地利用类型进行综合治理，以减少其土壤侵蚀。

（2）不同植被覆盖度的土壤侵蚀分析。对该流域植被覆盖度图与土壤侵蚀结果进行叠加，得到各个植被覆盖度等级的土壤侵蚀情况（表7）：其中高覆盖度土壤侵蚀总量仅占流域总侵蚀量的0.06%，中高覆盖度土壤侵蚀总量占流域总侵蚀量的1.03%，中覆盖度土壤侵蚀总量占流域总侵蚀量的10.10%，中低覆盖度土壤侵蚀总量占流域总侵蚀量的15.73%，高覆盖度土壤侵蚀总量占流域总侵蚀量的25.21%，其他（包括耕地等非植被）占流域总侵蚀量的47.87%。其中高覆盖年均侵蚀模数仅127.53t／（km2·a），而低覆盖达2 804.57t／（km2·a）。

表5 孤山川流域土壤侵蚀强度分级

表6 孤山川流域不同坡度等级土壤侵蚀情况

由表7得出，随着植被覆盖度的降低，年均侵蚀模数显著升高；由此可见植树造林与退耕还林工程建设对于控制土壤侵蚀的效果极为明显，建议该流域在今后的水土流失治理中要特别注重植树造林与退耕还林方面的水保效益，应积极巩固植树造林与退耕还林工程建设，加强植树造林与退耕还林工程管理，使其长效发挥其生态功能。

（3）不同土地利用类型的土壤侵蚀分析。将该流域土地利用现状结果图与土壤侵蚀图进行叠加分析，得到各个土地利用类型的土壤侵蚀情况（表8）。结果表明：没有任何水保措施的裸地年均侵蚀模数最高达12 264.69t／（km2·a），其次侵蚀比较严重的依次有居民点及交通工矿用地、荒草地、草地和坡耕地等。

表7 孤山川流域不同植被覆盖土壤侵蚀情况

表8 孤山川流域不同水土保持类型土壤侵蚀情况

由表8可以得出，该区域不同土地利用类型水保效益顺序依次：苗圃＞水浇地＞水域及水利设施用地＞鱼鳞坑＞梯田＞其他旱地＞果园＞乔木林＞灌木林＞疏林地＞其他林地＞坡耕地＞草地＞荒草地＞居民点及交通工矿用地＞裸地。据此可为该区域今后的水保措施建设方面提供一定的依据，例如该区域应积极将坡耕地采取“坡改梯”的措施等。

4 结 论

本文采用CSLE和高分辨率航空影像结合对流域进行土壤侵蚀定量分析，不仅克服了传统工程措施数据收集的困难，更重要的是可以精确地得到该流域植被覆盖、工程措施（主要是梯田和坝地）和耕作措施分布情况，提高了流域土壤侵蚀定量计算的精度。

但是在评价过程中还存在一些问题：① 本模型的降雨侵蚀力R因子应采用EI30来计算，但是由于长时间序列的降雨过程资料难以获得，故本研究采用的是简易的月均数据计算公式而未考虑雨强的影响。② 关于K因子的计算，由于缺乏该区域土壤类型数据，采用的是世界土壤数据库中的土壤类型图数据，土壤类型比较粗糙。故计算出的降雨侵蚀力和土壤可蚀性因子的精度还有待进一步提高。

通过本次对孤山川流域土壤侵蚀的定量研究分析，对孤山川流域的土壤侵蚀特点、分布规律及形成原因等有了一定认识，有利于掌握黄河中游多沙粗沙区孤山川等重点支流水土流失的发展趋势和水土保持措施的防治效果，可为制定有关政策和调整防治方略提供依据。

［1］ 张有全，宫辉力，赵文吉.基于遥感与GIS的土壤侵蚀强度快速估测方法［J］.生态环境，2007，16（1）：102－107.

［2］ 马丽梅，王万忠，焦菊英，等.黄河中游输沙与减沙的时空分异特征［J］.水土保持研究，2010，17（4）：67－72.

［3］ 高亚军，王玉明，赫晓慧.黄河中游严重水土流失区土壤粒径分布规律研究［J］.水土保持研究，2006，13（5）：27－29.

［4］ 康玲玲，李莉，刘立斌，等.水土保持措施蓄水拦沙指标在三川河流域的应用［J］.水土保持研究，2006，12（6）：86－88.

［5］ （SL190—2007）SL.中华人民共和国水利行业标准：土壤侵蚀分类分级标准［S］.2008.

［6］ Liu Baoyuan，Zhang Keli，Xie Yun.An Empirical Soil Loss Equation［C］∥ Proceedings 12th International Soil Conservation Organization Conference：Ⅱ.Beijing：Tsinghua University Press，2002.

［7］ 刘宝元.西北黄土高原区土壤侵蚀预报模型开发项目研究成果报告［R］.北京：水利部水土保持监测中心，2006.

［8］ 程琳，杨勤科，谢红霞，等.基于GIS和CSLE的陕西省土壤侵蚀定量评价方法研究［J］.水土保持学报，2009，23（5）：61－66.

［9］ 孙泽祥，孙希华，马倩倩.基于CSLE的淮河流域沂源县土壤侵蚀监测调查研究［J］.山东师范大学学报：自然科学版，2012，27（1）：94－97.

［10］ 刘宝元.西北黄土高原区土壤侵蚀预报模型开发项目研究成果报告［R］.北京：水利部水土保持监测中心，2006.

［11］ Wischmeier W H，Smith D D.Predicting rainfall erosion losses：aguide to conservation planning.Washington：United States Department of Agriculture，1978.58p［R］.Agriculture handbook，537.

［12］ Peel T C.The relation of certain physical characteristics to the erodibility of soils［J］.Soil Science Society Proceedings，1937，2：79－84.

［13］ McCool D K，Brown L C，Foster G R，et al.Revised slope steepness factor for the Universal Soil Loss Equation［J］.Transactions of the ASAE，1987，30：1387－1396.

［14］ Liu B Y，Nearing M A，Risse L M.Slope gradient effects on soil loss for steep slopes［J］.Transactions of the ASAE，1994，37（6）：1835－1840.

［15］ 郭建英，李锦荣，刘铁军，等.退耕还林工程建设对吴起县土地利用／覆被变化及其土壤侵蚀的影响［J］.水土保持研究，2013，（5）.

［16］ 谢红霞.延河流域土壤侵蚀时空变化及水土保持环境效应评价研究［D］.西安：陕西师范大学，2008.

［17］ 李勋贵，魏霞.区域土壤侵蚀演化风险分析：以泾河流域为例［J］.地理研究，2011，30（8）：1361－1369.

［18］ 吴发启，张玉斌，王健.黄土高原水平梯田的蓄水保土效益分析［J］.中国水土保持科学，2004，2（1）：34－37.

［19］ 许月卿，邵晓梅.基于GIS和RUSLE的土壤侵蚀量计算：以贵州省猫跳河流域为例［J］.北京林业大学学报，2006，28（4）：67－71.

［20］ 兰敏.基于GIS的秦巴山区土壤侵蚀空间特征分析：以陕西省宁强县为例［J］.干旱地区农业研究，2012，30（2）：196－200.

［21］ 石生新，蒋定生.几种水土保持措施对强化降水入渗和减沙的影响试验研究［J］.水土保持研究，1994，1（1）：66－69.

［22］ 李斌兵，郑粉莉，龙栋材，等.基于GIS纸坊沟小流域土壤侵蚀强度空间分布［J］.地理科学，2009，29（1）：105－110.