基于USLE原理和3S技术的水土流失定量监测方法及其应用研究*

马 力 卜兆宏† 梁文广 夏立忠 彭桂兰 徐姗姗 王冬梅 钟 敏

（1 中国科学院南京土壤研究所，南京 210008）

（2 江苏省水利科学研究院湖泊研究所/江苏省水利遥感工程研究中心，南京 210017）

（3 江苏省南京市水利局，南京 210036）

水土流失是造成区域生态环境恶化的重要影响因素。水土流失的产生是由自然因素决定并受到人为因素的影响。包括气候、植被、土壤、地形等在内的自然因素，决定着一个地区水土流失的产生和发展，而人为因素是水土流失的外在因素，它既有积极利用、改造自然的控制流失行为，也有过度利用、违背自然的加剧流失行为[1-2］。河湖水库水质富营养化的污染源，很重要一部分来自于由降雨引发水土流失进入水域的泥沙、养分等面源污染[3］。水土流失的监测原理与其成因密切相关，任何一种监测水土流失的方法或技术均必须获取该区域自然和人为各项因素的准确参数，才有可能准确解析该地区水土流失的实际状况。近一个多世纪以来，对于土壤侵蚀和面源污染定量监测的基本规律和方法，国内外学者已进行了大量卓有成效的研究，如通用土壤流失方程（Universal Soil Loss Equation, USLE/RUSLE）,水土资源评价工具SWAT（Soil and Water Assessment Tool），农业非点源模型AGNPS（Agricultural Nonpoint Source）等，经过数十年的修订完善，均是目前应用广泛的水土流失和面源污染预测预报模型[4-10］。我国学者近20年来在土壤可蚀性及相关模型因子的估算[11-12］，基于水蚀预报模型（Water Erosion Prodiction Project, WEPP)的应用[13］，以及基于GIS与RS技术在土壤侵蚀模型应用等方面做了大量工作[14］。刘宝元等[15-16］以RUSLE为基础，利用西部黄土高原区域的土壤侵蚀实测数据，建立了中国土壤流失模型CSLE (Chinese Soil Loss Equation)，并在我国不同区域得到了很好的应用。就不同土壤侵蚀监测预报模型的应用而言，因我国地域辽阔，农业用地的地形、地貌条件和管理方式等与国外具有较大差异，相关模型计算和参数的设定需要根据不同地区的实测数据来进行不断修订和完善，如果直接应用现有参数，则可能使部分区域的监测结果难以客观和准确。因此，基于大量实测数据，进一步开发完善适用于我国南北方实际并可大范围应用的水土流失和面源污染定量监测方法具有重要意义[17-19］。

本文阐述的水土流失定量监测方法于1983—2001年期间提出，受四项国家基金项目支持并在全国各地应用，2001年以后又进行深入研究、改进和完善，被改进应用的我国各省市或流域的累计总面积已超过30万km2[17-22］。该系统方法基于USLE的基本原理，利用我国南北方各水蚀区水土保持试验站径流小区监测和标桩法开展的土壤流失量及其影响因素周年实测获得的大量成套观测数据，建立了类似RUSLE又适宜我国南北方的以3S（GIS、GPS、RS）技术为支持的分布式水土流失定量监测模型，该模型监测年流失量精度接近水文站实测精度，目前已在我国山东全省山区、福建泉州市、太湖流域、云南龙川江流域、江西UNDP试区小流域、江苏全省（两次：1997—2000年和2001—2007年）和涉及北京河北的密云水库流域（两次：2001—2002年和2001—2005年）等地得到实际应用[6,18-24］。近年来，通过对土地利用、土壤养分状况、土壤流失对径流污染负荷形成影响的研究，又建立了分布式污染物负荷监测模型，形成了改进后的土壤与养分流失综合监测方法[20］。该方法具有操作方便、省时、能快速更新和监测精度高的优点。它不仅能获得流域出口或区域排出的固体和养分污染总量，而且能显示其固体泥沙和养分污染来源地。

该定量监测方法自提出以来，发表了一系列文章针对各模型因子算式算法和不同监测区域的应用进行了深入阐述[17-29］，但尚缺对该方法进行一次总结性的系统论述。从在我国不同地区的应用效果来看，在1997—2015年期间，对江苏南京全市区域内各区县和各小流域不同年份水土流失分布面积、程度及总量等时空变化状况进行了定量系统监测，该方法在南京市的应用时间最长并最具代表性[30-31］。因此，本文系统总结该方法的组成、结构和实现原理，旨在进一步总结完善该定量监测方法体系，为我国水土流失和面源污染一体化定量监测提供技术支撑；另一方面，研究结果也为今后不同区域的水土流失防治和生态环境保护工作提供科学决策依据。

1 材料与方法

1.1 研究区概况

本监测方法在我国山东、福建、云南、江苏、太湖流域、密云水库流域和江西UNDP试区小流域等多地得到实际应用[6,18-22］，其中在南京全市区域的监测应用时间最长，最具代表性[30-31］。南京市位于31°14′～32°36′N，118°32′～119°14′E，包括玄武、鼓楼、建邺、秦淮、雨花台、栖霞、江宁、六合、浦口、溧水和高淳共11个市辖区。全市总面积约6 590 km2，其中丘陵岗地面积占全市总面积的58.4%，平原洼地占41.6%（江河湖泊等水面占10.2%）。该市属于北亚热带季风气候区，年均气温17.8 ℃，日照时数1 687 h，四季明显，降水丰沛，年均降雨量1 034 mm，6—9月降雨量占全年的60%～70%。南京市主要水域为长江、秦淮河、水阳江、滁河、固城湖和石臼湖等，总蓄水量约4.28亿m3[30-31］。

1.2 数据资料获取与分析

本项监测工作选择的雨量代表站是南京下关、六合和天生桥站，境内外共39个雨量一般站。根据各代表站获得降雨特征数据计算出各次降雨的侵蚀动能E、降雨侵蚀力R值和I30B、汛期月份，进而获得代表站和一般站的R值。监测所使用的行政区界图、DEM图、土壤图、土地利用图等基础数据或纸图均由南京市水利局提供。卫星遥感数据图分别使用了2000年4月接收的南京地区TM遥感数据、2008年4月接收的中巴资源卫星CCD遥感数据和2014年5月1日接受的我国环境卫星HJ遥感数据。水土流失强度分级采用中华人民共和国水利部批准的《土壤侵蚀分类分级标准》（SL190-2007）。获得的所有数据资料通过系统软件进行汇总分析。

1.3 定量监测方法原理、组成与实现

水土流失定量监测方法所采用的模型因子算式算法较能反映水土流失和面源污染的成因，它不仅具备3S技术特征，而且具有更为快速、准确和一体化监测水土流失和面源污染的功能。该水土流失的定量新方法基于三条基础原理：（1）以实测数据建立监测模型与模型因子算式的原理；（2）以遥感（RS）与非RS数据图的复合配置配准实现流失量监测的原理；（3）以与水文站实测数据对监测结果检验较为客观公正的原理。该方法主要由水土流失定量监测预报模型、面源污染监测模型、模型因子算式算法原则及其软件、DEM精度评价新方法、GPS快速实测更新GIS的实用方法和监测结果建库查询法等几项内容组成[17-18,20,25,28］。

（1）土壤年流失量A监测模型

式中，A为土壤年流失量，单位t·km-2·a-1；R为降雨侵蚀力因子，单位MJ·mm·hm-2·h-1·a-1；K为土壤可蚀性因子，单位t·hm2·h·hm-2·MJ-1·mm-1；LS为地形的坡长因子L与坡度因子S之积，无量纲单位；CP为植被、作物覆盖因子C与水保措施因子P之积，无量纲单位；f为转换系数100，将土壤流失量A的单位t·hm-2·a-1转换为我国常用单位t·km-2·a-1（由于国际制单位的R值太大和K值太小，故在实际应用中R和K皆采用美国常用单位，即1美国常用单位的R等于17.02国际制单位的R，1美国常用单位的K等于0.1317国际制单位的K，而f为224.2）；i为栅格像元号。该模型形式与美国的USLE相同，但模型因子算式算法系根据我国各水蚀区大量实测数据所建，并且是基于栅格像元运算，故其算式算法和软件与USLE有显著不同[11,25］。主要的变化有以下几个方面：

R因子的算法：R因子是不同地区或同一地区不同年代的侵蚀流失量大小差异的主要决定因素。R因子计算需收集代表站的降雨自记纸，并按要求摘取符合侵蚀雨标准的各次降雨时段的ΔTj及其雨量Pj和连续30 min最大雨量P30，再按式（2）和式（3）计算各次降雨的动能E和R值及各月和年的R值[22, 26-27］。

在确定侵蚀汛期月份后，按式（4）计算该代表站的I30B年代表值。

按式（5）计算出代表站和一般站的R值。

式中，Pf为该区域汛期各月的降雨总量，mm，I30B为该区代表站的连续30 min最大降雨强度的年代表值，cm·h-1。代表站按式（3）和式（5）获得的两个R值应超过85%的一致性，是评价降雨数据摘取是否合适和准确计算R因子的关键[26］。

K因子的算法：K因子主要由土壤的机械组成粒级含量（N1，N2）、土壤有机质含量（OM）、土壤结构级别（S）、土壤渗透级别（P）等土壤性状决定。K因子计算需利用收集土壤理化分析数据。中美两国土壤的K值变化范围几乎相同，因而可采用美国的算式[2,17］。在土壤粒径分析中，一般缺乏与0.1 mm有关的数据，依据已有土壤机械分析数据建立2～0.1 mm含量算式和准确建立结构级别、渗透级别的确定原则是准确计算K因子的关键[18, 28］。

LS因子的算法：LS是侵蚀动力的加速因子。我国坡耕地的坡度一般较美国大，而且运用坡度因子计算需考虑0°～90°变化范围，故不能完全采用USLE的坡度因子算法。待求像元i的坡度θi算式为：

式中，hi为待求像元高程，hj为其邻像元高程，j有1～8个方向；Di为其与邻像元间距，为像元边长或像元边长。

依据大量成套实测数据按照不同的像元坡度有θi0°～1.4°、1.4°～3.5°、3.5°～5.14°、5.14°～8°、8°～40.72°和40.72°～90°的坡度因子。Si算式为式（8），像元坡度＜5.14°的Si与USLE基本相同（但因我国谷地平原水田或水浇地皆有地埂，故较其小1/3或1/2）。

其他三种坡度形式均为：Si=Α·Bθi或Si=A·B ·sinθi。其中，坡度θi为5.14°～8°的像元坡度因子算式：

像元坡长算式：

式中，m值随坡度而变，当坡度θi≥5.1 4°时m=0.5；当5.14°＞θi≥1.72°时m=0.4；当1.72°＞θi≥0.57°时m=0.3；当θi＜0.57°时m=0.2。

本方法中像元地形因子LSi计算已由开发的专用软件完成。它需要利用数字高程模型DEM为输入数据，输出像元坡度θi图和像元坡度因子与像元坡长因子乘积LSi图。该DEM在输入前必须评价其适用性和高程精度，必要时重新构建与评价，因为它是水系、水文雨量站位R值、土壤类型图、土地利用图和遥感数据图像等能套合配置配准的几何校正的依据基础。具有专用软件和获得适用的DEM，是准确计算LSi因子的关键[18,28］。

CP因子的算法：CP因子是侵蚀动力的抑制因子，起着水土保持作用。C、P单个因子值小于或等于1。CPi像元因子积的算法是直接利用遥感数据排序特征的赋值算法。在IDRISI系统软件中利用不同波谱段（蓝430～520 nm，绿520～600 nm，红630～690 nm和近红外750～900 nm）遥感数据按值选取几类水体和十几类植被训练样区，并建立各类在各波谱段的标识文件和运行最大似然分类法获得遥感分类图；将土地利用图归并为三大类（一大类含水田、水浇地、有埂水平梯地、道路、城镇居民地、水域等流失量很小区域，二大类含旱作地、工矿用地、果园、茶园、菜地等土壤扰动流失量较大区域，三大类含林地、草地、灌木地、裸土、裸岩等荒山自然植被区域）后，分别以三大地

像元坡长因子算式：类分割遥感分类图；依据实测数据编制各大类CP赋值文件（一大地类赋值0.0001～0.2，二大地类赋值0.0001～0.6，三大类赋值0.0001～0.45），并运行开发的专用赋值软件分别运算获得三个地类的遥感分类图上每类植被覆盖度的CP值；最后用IDRISI系统的逻辑加将三个赋了CP值的遥感分类图合并为全监测区的植被覆盖与保土措施CPi像元因子图。在IDRISI系统中用0.63～0.69 μm波谱段选准样区并建好赋值文件，是准确运算CPi像元因子图的关键[18］。

（2）防治强度Y预报模型

式中，CP’为确保在降雨侵蚀力周期均值R时土壤年流失量≤500 t·km-2所要求的植被与保土措施因子之积值，数值为500/(fRKSL)。当Yi为负值时，表明为无需治理区；反之，则为预防监督区、需要或急需治理区[18］。

（3）土壤养分面源污染Wsp监测模型

土壤养分面源污染(年增量)WSp监测模型为[25-26］：

式中，WSp的单位为kg·km-2;Ai为所监测i像元被侵蚀进入水域的土壤年流失量，单位为 t·km-2；Spj为土壤j种(分别有全氮、全磷、铵态氮、有效磷)养分在所监测土壤利用地类中的含量，全氮和全磷的单位为 g·kg-1，铵态氮和有效磷的单位为mg·kg-1。

（4）监测方法的软件系统与实现

监测水土流失的定量方法建立了更符合流失实际的模型因子算式算法和原则，如LS、K中土壤粒级制转换等算式、DEM构建法、R算法和汛期月份的确定原则，均进行了改进，使得监测结果更符合当地实况和流失实际。模型算式算法采用了当前通用的软件系统，主要有ArcView、Arc/Info或ArcGIS和IDRISI，并自行研制了它们之间互相转换数据接口的专用软件。通过运行专用软件，使建立一个区域（流域或省）定量监测系统的速度加快、精度更高，也为所建系统实现和完成每年监测的时间缩短为1～2个月。同时，可以完成遥感和非遥感的GIS、GPS数据的几何校正、信息采编、矢栅转换和按模型要求进行运算，最终获得定量监测结果和建立信息查询系统。本方法采用的GPS仪由ArcPad、PDA和GPS接收器自行组装，实现GPS与GIS的集成，可直接由4颗以上GPS卫星解算WGS84 (World Geodetic System，1984)经纬度并转换为北京54坐标的新的GIS图（如土地利用图），可实现GIS快速实测更新。

本方法的软件系统和实现原理如图1所示。方法同时建立了一套定量监测结果的精度评价法，所有水土流失和面源污染的监测结果，均采用与水文站点测流速、流量和采水样分析后，测算的下排泥沙量、养分污染量，进行比较评价，以其相对一致性作为监测精度，这样评价较客观。此外，还有建立了监测结果的建库查询方法，以便于治理应用[23,31］。

图1 水土流失定量监测方法软件系统与实现原理Fig. 1 Software system and principle used for realization of quantitative monitoring of soil and water losses

（5）模型监测结果的精度评价方法

为评价模型监测结果的精度，除了野外实地核对监测结果图是否符合实际外，主要依据上下游水文站点水文观测计算的输出泥沙量（Ngc上Ngc下）、养分量（Yfc上Yfc下）进行评价。首先，利用上游汇水界线和含上游的下游汇水界线分别统计出监测图的上游流失总量 A上总、各养分总量Yf上总和含上下游的流失总量A上下游总、各养分Yf上下游总；然后，按式（13）和式（14）计算泥沙、养分下排系数和监测结果图的精度。

2 结 果

2.1 2001—2015年南京市水土流失状况变化

采用水土流失定量监测方法分析计算了1997—2015年期间研究区域包括降雨侵蚀因子R、植被覆盖和水保措施因子CP、地形因子LS和土壤可蚀因子K值，并获得了全市区域的栅格像元图（图2），再进一步按照模型公式运算获得了全市区域年水土流失面积和总量数据，并利用ArcGIS软件，结合水系、行政区划等数据，完成了全市水土流失定量监测图，这里给出2001—2010年和2011—2015年两个时间段南京全市水土流失总量平均值的定量监测图（图3）。

图2 南京市年降雨侵蚀因子R（a）、植被覆盖和水保措施因子CP（b）、地形因子LS（c）和土壤可蚀因子K（d）值像元图Fig. 2 Pixel map of rainfall erosivity factor R (a), vegetation cover and soil conservation factor CP (b), slope length factor LS (c) and soil erosibility factor K (d) in Nanjing area

图3 南京市水土流失总量平均值（a. 2001—2010和b. 2011—2015）的监测图Fig. 3 Monitoring map of mean total soil/water loss (2001—2010 and 2011—2015) in Nanjing area

由南京全市1997—2015年的连续19年的监测结果（图4和图5）得出，全市19年平均轻度以上年水土流失面积863.7 km2，占全市总面积的13.0%；平均微度水土流失面积5091.5 km2，占全市总面积的76.8%；水域面积671.02 km2，占全市总面积的10.2%。全市19年期间的平均年水土流失总量为276.5×104t，其中轻度以上流失量为229.2×104t，占全市流失总量的8 2.9%；微度水土流失量45.3×104t，占全市总流失量的16.4%。南京市轻度以上平均水土流失面积所占比例较小，但水土流失量占总流失量80%以上，而大于80%的区域属于微度水土流失，其水土流失量则不足全市流失总量的20%。从水土保持角度来看，轻度以上水土流失治理仍是水土流失治理中的重点方向。年际间轻度以上水土流失面积与流失量变化较大，但总体存在下降的趋势。由图4和图5可以看出，在19年期间，1998年轻度以上水土流失面积和流失量为19年中最大（1 509.9 km2和414.3×104t），2013年轻度以上水土流失面积和流失量为19年中最小（367.2 km2和102.3×104t），水土流失面积、流失量的大小，与年降雨量大小有关，特别是汛期雨量为决定因素。结果显示，19年期间平均极强烈与剧烈水土流失面积分别为39.0 km2和20.2 km2，分别占轻度以上水土流失面积4.5%和2.3%，但水土流失量分别为41.8×104t和50.8×104t，分别占轻度以上水土流失量的18.2%和22.2%，其侵蚀模数为分级侵蚀量中最大。换句话说仅6.8%的极强烈和剧烈水土流失面积，却占有轻度以上水土流失量的40%以上。这说明了当降雨集中和雨量大的时期，如果生产建设活动未采取有效水保措施或措施不到位，在全市范围水土流失等级高和生态环境较差的区域将会产量严重的水土流失，进而加剧土壤侵蚀和生态环境恶化。

R因子均值的变化是全市不同年份的侵蚀流失量大小差异的主要决定因素，由图6看出，南京市19年期间的R均值有明显的增大趋势，这与该地区气象要素变化，如降雨量和降雨强度增大有密切关系，说明在此期间降雨侵蚀力有所增强，而图4和图5结果则反映了水土流失面积和流失总量的降低，这说明了当地水土保持措施产生了明显效果，在降雨侵蚀力增加的条件下，有效减轻了该地区的水土流失状况。而由图6也发现，全市不同年份的土壤侵蚀模数也有随时间变化增大的趋势，产生此现象原因是南京市在此19年期间的快速建设发展，使土地利用方式发生了较大变化，在今后的水土流失定量监测工作中，应当注重监测区域土地利用图的及时更新，以满足未来新的监测工作需要。

图4 南京市1997—2015年不同等级水土流失面积变化Fig. 4 Variation of soil erosion areas of different erosion levels from 1997 to 2015 in Nanjing area

图5 南京市1997—2015年不同等级水土流失总量变化Fig. 5 Variation of total soil/water loss from 1997 to 2015 in Nanjing area relative to erosion level

2.2 1997—2015年南京市水土流失治理成效

自1997年开始至今，本项监测方法在南京市进行了三次监测成果的分析（1997—2000，2001—2010和2011—2015）。2011—2015年5年平均流失量与前两次的结果相比，轻度以上水土流失面积和流失量均明显下降（图7）。与1997—2000年监测结果相比，轻度以上水土流失面积由1288.7 km2降至515.0 km2，减少了60.0%；年土壤流失量由335.6×104t降至166.8×104t，减少了50.3%；结合前文对1997—2015年南京市R均值变化的分析，这说明近年来，全市水土流失治理措施和水土保持监督管理措施加大，成效非常明显。南京市水土保持不同流失等级监测结果表明，微度水土流失面积由4630.3 km2增加至5403.1 km2，增加了16.7 %，说明了相当一部分轻度以上级别的水土流失区域的土壤侵蚀强度降到微度等级，治理卓有成效；极强烈和剧烈水土流失面积分别由57.6和33.4 km2减少至34.4和15.1 km2，降幅分别为40.2%和53.6%，极强烈水土流失总量分别由62.0和84.4×104t减少至 36.5和37.3×104t，降幅为41.1%和55.8%，说明部分流失等级高和生态环境较差区域的土壤侵蚀也得到了有效控制，这体现了水保措施的实施和植被覆盖面积增加带来的效果。通过实地调查验证了监测结果的可靠性，但也发现南京市范围一些区域的开发建设和采石、采矿业对周边生态环境造成的严重破坏，尚未采取有效的根治措施，使得个别区域的土壤侵蚀和环境破坏情况加剧，需要引起足够的重视。

图6 南京市1997—2015年R均值（a）和土壤侵蚀模数（b）变化Fig. 6 Variation of mean R value (a) and soil erosion modulus (b) from 1997 to 2015 in Nanjing area

2.3 水土流失定量监测结果的可靠性

在运用本方法开展南京市水土流失定量监测过程中，保证了GIS、RS和GPS三种技术的有效结合，通过实测采样点、水文站或段面位置，订正和配准已有的图件，使3S三种数据不仅能套合一致，且界线和内容更符合当地实际。基于固体泥沙和氮磷养分的面源污染来自流失的土壤，而植被覆盖率是减少降雨侵蚀土壤的重要因素，故以遥感数据像元作为成果运算的基础，不仅充分发挥了遥感数据反映植被信息的优势，而且确保了监测结果的精细详尽。模型方法的科学性以及重视从野外采样到室内化验和数据处理中各环节的质量查验，确保了所获监测成果有较高的可信度。对于监测结果，采用与上、下游水文站观测、分析计算的泥沙与水质数据比较的方法评价其精度，较其他评价方法更加客观公正和可信。在监测结果分析中进行了野外实地核对，并与已有前期的监测数据进行对比分析，最后与水文实测泥沙数据做出分析。实地核查表明，凡是中度以上水土流失等级的监测数据均与实地基本一致。与水文实测数据的比较分析表明，水土流失的平均监测精度达到了83.36%，土壤全氮流失的平均监测精度为78.1%，土壤全磷流失的平均监测精度为87.8%。因此，在采用本系统方法获得的监测结果是合理和可靠的。

3 结 论

本项水土流失定量监测方法是基于通用土壤流失方程（USLE）的基本原理，利用我国南北方各水蚀区的土壤流失量及其相关生态因子的大量实测数据，建立适宜我国南北方的以3S技术为支持的定量监测水土流失方法。该方法建立了更符合流失实际的模型因子算式算法和原则，监测精度和可靠性较高，已在我国山东、福建、江苏、云南、北京、河北等地累计应用面积超过30万km2，应用效果得到了较好的验证。掌握水土流失定量监测方法中各模型因子运算的关键操作，是准确快速完成新老监测区任务的必要条件。该定量监测方法在南京市的应用时间最长并最具代表性，在1997—2015年的19年期间，南京全市平均年水土流失总量为276.5×104t，平均轻度以上年水土流失面积863.7 km2，占全市总面积的13.0%。轻度以上水土流失治理仍是南京市今后水土流失治理和生态环境保护的重点方向。与2000年之前的监测结果相比，轻度以上，尤其是极强烈和剧烈等级的水土流失面积和年土壤流失总量均有大幅的降低，微度水土流失面积明显增加，说明近年来南京全市水土流失治理措施和生态保护监督管理措施有所加大，治理成效非常明显。

图7 南京市1997—2015年三次定量监测不同流失强度下的流失面积(a)和流失总量(b)Fig. 7 Comparison of the three rounds of quantitative monitoring during the period from 1997 to 2015 in Nanjing area in erosion area(a) and total soil/water loss (b) relative to erosion intensity