坡长坡度因子计算工具

符素华，刘宝元†，周贵云，孙中轩，朱小立

(1.北京师范大学 地表过程与资源生态国家重点实验室 地理学与遥感科学学院，100875，北京；2.电子科技大学资源与环境学院，611731，成都)

坡长坡度因子计算工具

符素华1，刘宝元1†，周贵云2，孙中轩2，朱小立1

(1.北京师范大学 地表过程与资源生态国家重点实验室 地理学与遥感科学学院，100875，北京；2.电子科技大学资源与环境学院，611731，成都)

地形(坡长坡度)因子是坡面土壤侵蚀模型USLE(通用土壤流失方程)或CSLE(中国水蚀方程)中的重要参数。本文选择了适合我国土壤侵蚀特点的坡长坡度因子计算公式，基于Visual Studio 2010平台进行了程序编写，开发了LS计算工具。该工具界面友好且计算速度快，在32位计算机上可快速计算1万行×1万列数据区域的坡长坡度因子，在64位计算机上能计算4万行×4万列数据区域的坡长坡度因子。本软件的开发可为区域土壤侵蚀评价以及水土保持措施规划服务。

坡度因子； 坡长因子； 土壤侵蚀； USLE； CSLE

在通用土壤流失方程USLE和中国土壤流失方程CSLE中都用地形因子来反映地形对土壤侵蚀的影响。这2个模型的地形因子定义一样，是坡长坡度因子的统称。其中坡长因子是指降雨、土壤、坡度、地表状况等条件一致时，某种坡长的坡面土壤侵蚀量与22.13 m坡长的坡面土壤侵蚀量比值[1]，该比值反映了土壤侵蚀量与坡长的定量关系。坡度因子是指其他条件一致的情况下，某坡度下的坡面土壤侵蚀量与坡度为5.14°时的坡面土壤侵蚀量比值[1]；这个比值反映了土壤侵蚀量与坡度之间的定量关系。当USLE和CLSE应用于区域时，地形因子在数字高程模型DEM的基础上生成，计算过程较为复杂，不能直接由通用ArcGIS软件生成，限制了USLE和CSLE在区域土壤侵蚀评价中的应用；因此R.Hickey等[2-3]先后用Arc map软件的宏语言开发了基于USLE[4]版本坡长坡度因子计算公式的算法。R.D.VanRemortel等[5]采用Arc map软件宏语言开发了基于RUSLE[6]版本坡长坡度因子计算公式的算法。2004年，R.D.Van Remortel等[7]用ANSI C++TM语言改写了R.D.Van Remortel等[5]2001年的算法。我国杨勤科等[8]在R.D.Van Remortel 等[7]算法的基础上，修改了水流流向的算法，并结合我国的实际情况，坡度因子计算公式中增加了Liu Baoyuan等[9]的陡坡坡度因子公式。

但在这些计算LS因子的版本里，坡长因子的计算公式都采用的是整坡坡长公式。该公式仅适用于均匀坡，不分段的整坡情况。对于用若干栅格来代表的任何一个区域，每一个栅格仅是一面坡上的一段。在此情况下用整坡坡长因子公式来计算每一栅格的坡长因子就会带来误差[10]。此时应采用G.R.Foster等[11]1974年提出的分段坡坡长因子来计算区域上每一栅格的坡长因子[10，12]。到目前为止，国内外还没有采用分段坡坡长因子公式计算坡长因子的软件；因此本研究是综合国内外的最新研究成果，集成一个能运用于我国区域坡面土壤侵蚀量评价的LS因子计算工具，为我国区域土壤侵蚀评价以及水土保持规划服务。

1 区域坡长坡度因子计算原理

根据国内外的研究成果，在不同的坡度范围下分别选用不同的公式来计算坡度因子。10°以下的坡度选用D.K.McCool[13]的公式：

S=10.8sinθ+0.03，θ<5°；

(1)

S=16.8sinθ-0.5，5°≤θ<10°。

(2)

10°以上的坡度选用Liu Baoyuan等[9]的公式：

S=21.91sinθ-0.96，θ≥10°

(3)

式中：S为坡度因子；θ为坡度，(°)。本软件中坡长因子设计了2种算法。一是采用G.R.Foster 等[10]于1974年提出的分段坡坡长因子公式来计算区域上每一栅格的坡长因子

(4)

式中：Li为第i个栅格的坡长因子；λout、λin分别为栅格出口及入口的坡长，m；m为坡长指数。

根据刘宝元等[14]的研究结果，坡长指数m取值如下：

m=0.2，θ<0.5°；

m=0.3，0.5°≤θ<1.5°；

m=0.4，1.5°≤θ<3°；

m=0.5，θ≥3°。

二是采用P.J.J.Desmet等[15]于1996年提出的考虑汇流对坡长因子影响的公式

(5)

式中：Aout、Ain分别为栅格出口及入口的汇流面积，m2；Δx为栅格分辨率，m；L为与栅格入口、出口水流方向相关的非累计坡长，m。

2 LS因子计算工具设计及实现

2.1 总体设计及计算框图

设计LS因子计算工具的目的是实现LS因子的计算，便于用户使用。数据基础是DEM栅格。计算流程图如图1所示。

2.2 计算步骤

1)填洼处理。DEM中一般会存在一些高程值低于周围的凹陷点，这些点无法进行流向判断而造成流路中断。本次计算采用L.W.Martz等[16]提出的扫描窗口法进行填洼处理。其主要思路是先找到每个平地栅格或是洼地底点栅格(入流栅格)，并标记出这个栅格，然后从已标记的栅格中找出潜在的出流点，找到最低的潜在出流点后，比较其和洼地栅格的高程。如果出流点高程高于洼地栅格，那么洼地是一个凹地，否则是一个平地。对于凹地，把洼地集水区域内所有低于出流点的栅格高程升高至出流点高程。这样，凹地就成为一个平地。对于平地，按照L.W.Martz等[16]中使用起伏平地的算法进行处理。

图1 LS因子计算流程Fig.1 Flow chart illustrating the process of calculating LS factors

2)计算水流流向。采用最大坡降算法，即栅格坡度的最佳代表值是以之为中心3×3窗口内其周围8个方向坡度最大值。水流方向与最大坡降一致。

3)沟道提取。如果栅格流入汇水面积大于汇水面积阈值，则认为当前栅格为沟道。栅格流入汇水面积为流入当前栅格的所有上游栅格面积之和汇流面积阈值由用户输入，其值的确定原则是让由DEM生成的沟道尽可能地与实际情况相符。

4)计算非累计坡长L/m。即单个栅格对后续流入及流出坡长的贡献，计算规则如下：

(1)如果没有水流流入当前单元格，非累计坡长L为栅格大小Δx/m；

(2)如果入流方向与出流方向同为东北、西北、东南和西南方向，则L= 1.414 2Δx；

(3)如果入流方向与出流方向同为东、西、南和北方向，则L=Δx；

(4)其他情况，L=1.207Δx。

5)栅格入口/出口坡长计算。依据水流流向和中断因子计算当前栅格的入口坡长。栅格入口坡长为流入栅格中出口坡长最大值；栅格出口坡长为栅格入口坡长加上与入流及出流方向对应的非累计坡长。根据坡长的定义，在坡长计算过程中用以下2个控制条件来确定流出坡长终点。

(1)沟道：当栅格是沟道时，则认为坡长终止，该栅格的流出坡长为0。

(2)中断因子：中断因子用于确定坡面坡度由陡坡到缓坡变化时，由于出现泥沙淤积而导致的坡长中断。中断因子定义为当前栅格的坡度值与其流入方向栅格坡度比值的临界值，这一比值反映了当前栅格坡度相对于流入方向栅格坡度的减小程度：若当前栅格的坡度小于中断因子与流入栅格坡度的乘积时，会出现泥沙淤积，坡长中断，此时当前栅格的流出坡长重新设置为栅格大小。本次计算时，缓坡(坡度<3°)的中断因子设为0.7，陡坡(坡度≥3°)的中断因子设为0.5。

6)计算坡长坡度因子。依据坡度的不同分别采用式(1)～式(3)计算坡度因子；根据用户输入的坡长因子计算方法采用式(4)或式(5)计算坡长因子。

2.3 技术实现

软件的开发工具为Visual Studio 2010，软件算法由C++语言实现，软件界面使用C#语言实现，文件读写用开源库GDAL来实现，图形显示不依赖任何第三方图形库。

软件安装时需要Windows XP或以上操作系统，Microsoft.NET Framework 2.0或以上环境；LS计算工具支持32位和64位2个操作系统的计算机；32位版本能够支持的数据一般不超过1万行×1万列的大小，我国大部分县计算一次即完成LS因子的计算；64位版本能支持的数据大小由系统的物理内存来决定。一个8 GB内存的64位系统能够处理2万行×2万列的浮点型栅格数据。一个32 GB内存的64位系统可以处理4万行×4万列的浮点型栅格数据。

软件设计考虑了用户使用的方便性，界面友好、清楚，参数输入界面如图2所示。

图2 LS计算工具软件用户界面Fig.2 User interface of the software for calculating LS factors

需要输入的DEM文件格式为栅格文件*.aux、*.xml或*.tif。程序主要输出以下变量(表1)。

表1 LS计算工具输出文件

表1(续)

3 应用案例

以青海玛多县为试验区域，计算了该地区的坡长坡度因子。玛多县行政面积为2.5万km2，海拔在3 915～5 262 m之间，平均海拔为4 392 m。平均坡度为8.2°，58%的区域坡度小于8°。计算时运用了30 m分辨率的DEM，该县共有7 500×6 080个栅格，DEM文件大小为43 MB，仅运用不到2 min的时间完成计算，说明该程序具有很快的运算速度。计算的坡度、坡长、以及坡长坡度因子分别如图3所示。统计参数如表2所示。

表2 主要输出图层统计参数

图3 玛多县坡度、坡长以及坡度坡长因子分布图Fig.3 Distribution of slope steepness, slopes length, slope steepness factor and slope length factor in Maduo County

4 结论

本文应用Visual Studio 2010平台开发了LS因子计算工具，该软件界面简洁、清晰，便于用户操作。计算结果可在软件中直接查看，也可在Arc map软件中进行查看。运用在32位计算机上可进行一般县级区域的LS因子计算，在64位计算机上可运算更大区域的LS因子值。应用案例说明，该工具软件具有很快的运算速度。该工具软件的开发，对于区域土壤侵蚀评价和水土保持规划提供了极为有效的工具。

[1] Wischmeier W H, Smith D D. Rainfall erosion losses from cropland east of the rocky mountains. Guide for selection of practices for soil and water conservation[M]. Agriculture Handbook 282, 1965：8-9

[2] Hickey R, Smith A, Jankowski P. Slope length calculations from a DEM within ARC/INFO GRID[J]. Computers,Environment and Urban Systems, 1994, 18(5): 365-380

[3] Hickey R. Slope angle and slope length solutions for GIS[J]. Cartography, 2000, 29(1): 1-8

[4] Wischmeier W H, Smith D D. Predicting rainfall erosion losses[M]. USDA Agricultural Handbook 537, 1978：12-15

[5] Van Remortel R D, Hamilton M E, Hickey R J. Estimating the LS factor for RUSLE through iterative slope length processing of digital elevation data within Arclnfo grid[J]. Cartography, 2001, 30(1): 27-35

[6] Renard K G, Foster G R, Weesies G A, et al. RUSLE a guide to conservation planning with the revised universal soil loss equation[M]. USDA Agricultural Handbook 703,1997：105-117

[7] Van Remortel R D, Maichle R W, Hickey R J. Computing the LS factor for the Revised Universal Soil Loss Equation through array-based slope processing of digital elevation data using a C++executable[J]. Computers &Geosciences, 2004, 30(9): 1043-1053

[8] 杨勤科,郭伟玲,张宏鸣,等.基于 DEM 的流域坡度坡长因子计算方法研究初报[J].水土保持通报, 2010,30(2): 203-206

[9] Liu Baoyuan, 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

[10] Fu Suhua, Wu Zhiping, Liu Baoyuan, et al. Comparison of the effects of the different methods for computing the slope length factor at a watershed scale[J]. International Soil and Water Conservation Research, 2013, 1(2): 64-71

[11] Foster G R, Wischmeier W H. Evaluating irregular slopes for soil loss prediction[J]. Trans ASAE Gen Ed Am Soc Agric Eng, 1974, 17: 305-309

[12] Fu Suhua, Cao Longxi, Liu Baoyuan, et al. Effects of DEM grid size on predicting soil loss from small watersheds in China[J]. Environmental Earth Sciences,2014, 73(1)：2141-2151

[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-American Society of Agricultural Engineers (USA), 1987, 30(5): 1387-1396

[14] 刘宝元，毕小刚，符素华，等.北京土壤流失方程[M].北京：科学出版社,2010: 60

[15] Desmet P J J, Govers G. A GIS procedure for automatically calculating the USLE LS factor on topographically complex landscape units[J]. Journal of Soil and Water Conservation, 1996, 51(5): 427-433

[16] Martz L W, Garbrecht J. Numerical definition of drainage network and subcatchment areas from digital elevation models[J]. Computers & Geosciences, 1992, 18(6): 747-761

(责任编辑：程 云 郭雪芳)

Calculation tool of topographic factors

Fu Suhua1，Liu Baoyuan1，Zhou Guiyun2，Sun Zhongxuan2，Zhu Xiaoli1

(1.State Key Laboratory of Earth Surface Processes and Resource Ecology, School of Geography, Beijing Normal University, 100875，Beijing, China；2.School of Resources and Environment, University of Electronic Science and Technology of China, 611731，Chengdu, China)

Topographic (slope length and slope gradient) factors (LS) are important parameters in the soil erosion model, for example, universal soil loss equation (USLE) and Chinese soil loss equation (CSLE ). TheLSfactor was usually computed using digital elevation models (DEM) for basin-wide application of the USLE and CSLE. The calculating process is very complicated and is difficult to be directly calculated using common GIS software such as ArcMap software. In this paper, anLStool software is developed on the platform of Visual Studio 2010 software. Source codes are written using C++language. The C++language is used to obtain the window of the software. This tool is easy to use with a friendly interface. To extend the tool suitability, the slope gradient factor equation at steep slope is added in the algorithms. An equation considering segmented slope situation is used to calculate the value of slope length factor. The calculation progress includes the following six steps: 1) filling topographical depression, 2) calculating flow direction, 3) extracting gully net, 4) calculating non-cumulative slope length (NCSL) of each grid cell, 5) calculating the cumulating slope length of each grid cell and 6) calculatingLSfactors of each grid cell. The above six steps are described in detail in this paper. The cutoff slope factor and gully net are used to stop the cumulating slope length. The input file is a DEM file with *.aux, *.xml or *.tif format. The outputs include slope gradient, slope gradient factor, gully, slope length and slope length factor etc.. The parameters of threshold values including slope gradient, slope length, slope cutoff factor, channel initiation and gully length are optional of user input. Maduo county located in Qinghai Province was used as an example to test the application of the software. The area of the county is 25 000 km2. The DEM with 30 m resolution includes 7 500 rows and 6 080 columns. The run time was less than 2 min on the computer with the 32-bit operating system to finish calculating. The application results show that the software has a high calculation capacity and runs efficiently. For the 32-bit operating system, the software can be used to calculate theLSfactors of a region with 10 000 rows and 10 000 columns; for a 64-bit operating system, it can be used for a region with 40 000 rows and 40 000 columns. This tool can be used as a sub-model to evaluate soil loss and to plan the soil conservation practice at a region scale.

slope steepness factor; slope length factor; soil erosion; USLE; CSLE

2015-03-18

2015-08-07

项目名称：中央高校基本科研业务费专项资金资助

符素华(1973—)，博士，教授，博士生导师。主要研究方向：水土流失与水土资源管理。E-mail: suhua@bnu.edu.cn

†通信作者简介：刘宝元(1958—)，博士，教授，博士生导师。主要研究方向：土壤侵蚀。E-mail: baoyuan@bnu.edu.cn

S157.1

A

1672-3007(2015)05-0105-06