坡長坡度因子計算工具

符素華，劉寶元?，周貴云，孫中軒，朱小立

(1.北京師范大學 地表過程與資源生態(tài)國家重點實驗室 地理學與遙感科學學院，100875，北京；2.電子科技大學資源與環(huán)境學院，611731，成都)

坡長坡度因子計算工具

符素華1，劉寶元1?，周貴云2，孫中軒2，朱小立1

(1.北京師范大學 地表過程與資源生態(tài)國家重點實驗室 地理學與遙感科學學院，100875，北京；2.電子科技大學資源與環(huán)境學院，611731，成都)

地形(坡長坡度)因子是坡面土壤侵蝕模型USLE(通用土壤流失方程)或CSLE(中國水蝕方程)中的重要參數(shù)。本文選擇了適合我國土壤侵蝕特點的坡長坡度因子計算公式，基于Visual Studio 2010平臺進行了程序編寫，開發(fā)了LS計算工具。該工具界面友好且計算速度快，在32位計算機上可快速計算1萬行×1萬列數(shù)據(jù)區(qū)域的坡長坡度因子，在64位計算機上能計算4萬行×4萬列數(shù)據(jù)區(qū)域的坡長坡度因子。本軟件的開發(fā)可為區(qū)域土壤侵蝕評價以及水土保持措施規(guī)劃服務。

坡度因子； 坡長因子； 土壤侵蝕； USLE； CSLE

在通用土壤流失方程USLE和中國土壤流失方程CSLE中都用地形因子來反映地形對土壤侵蝕的影響。這2個模型的地形因子定義一樣，是坡長坡度因子的統(tǒng)稱。其中坡長因子是指降雨、土壤、坡度、地表狀況等條件一致時，某種坡長的坡面土壤侵蝕量與22.13 m坡長的坡面土壤侵蝕量比值[1]，該比值反映了土壤侵蝕量與坡長的定量關系。坡度因子是指其他條件一致的情況下，某坡度下的坡面土壤侵蝕量與坡度為5.14°時的坡面土壤侵蝕量比值[1]；這個比值反映了土壤侵蝕量與坡度之間的定量關系。當USLE和CLSE應用于區(qū)域時，地形因子在數(shù)字高程模型DEM的基礎上生成，計算過程較為復雜，不能直接由通用ArcGIS軟件生成，限制了USLE和CSLE在區(qū)域土壤侵蝕評價中的應用；因此R.Hickey等[2-3]先后用Arc map軟件的宏語言開發(fā)了基于USLE[4]版本坡長坡度因子計算公式的算法。R.D.VanRemortel等[5]采用Arc map軟件宏語言開發(fā)了基于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因子的版本里，坡長因子的計算公式都采用的是整坡坡長公式。該公式僅適用于均勻坡，不分段的整坡情況。對于用若干柵格來代表的任何一個區(qū)域，每一個柵格僅是一面坡上的一段。在此情況下用整坡坡長因子公式來計算每一柵格的坡長因子就會帶來誤差[10]。此時應采用G.R.Foster等[11]1974年提出的分段坡坡長因子來計算區(qū)域上每一柵格的坡長因子[10，12]。到目前為止，國內(nèi)外還沒有采用分段坡坡長因子公式計算坡長因子的軟件；因此本研究是綜合國內(nèi)外的最新研究成果，集成一個能運用于我國區(qū)域坡面土壤侵蝕量評價的LS因子計算工具，為我國區(qū)域土壤侵蝕評價以及水土保持規(guī)劃服務。

1 區(qū)域坡長坡度因子計算原理

根據(jù)國內(nèi)外的研究成果，在不同的坡度范圍下分別選用不同的公式來計算坡度因子。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年提出的分段坡坡長因子公式來計算區(qū)域上每一柵格的坡長因子

(4)

式中：Li為第i個柵格的坡長因子；λout、λin分別為柵格出口及入口的坡長，m；m為坡長指數(shù)。

根據(jù)劉寶元等[14]的研究結果，坡長指數(shù)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因子計算工具設計及實現(xiàn)

2.1 總體設計及計算框圖

設計LS因子計算工具的目的是實現(xiàn)LS因子的計算，便于用戶使用。數(shù)據(jù)基礎是DEM柵格。計算流程圖如圖1所示。

2.2 計算步驟

1)填洼處理。DEM中一般會存在一些高程值低于周圍的凹陷點，這些點無法進行流向判斷而造成流路中斷。本次計算采用L.W.Martz等[16]提出的掃描窗口法進行填洼處理。其主要思路是先找到每個平地柵格或是洼地底點柵格(入流柵格)，并標記出這個柵格，然后從已標記的柵格中找出潛在的出流點，找到最低的潛在出流點后，比較其和洼地柵格的高程。如果出流點高程高于洼地柵格，那么洼地是一個凹地，否則是一個平地。對于凹地，把洼地集水區(qū)域內(nèi)所有低于出流點的柵格高程升高至出流點高程。這樣，凹地就成為一個平地。對于平地，按照L.W.Martz等[16]中使用起伏平地的算法進行處理。

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

2)計算水流流向。采用最大坡降算法，即柵格坡度的最佳代表值是以之為中心3×3窗口內(nèi)其周圍8個方向坡度最大值。水流方向與最大坡降一致。

3)溝道提取。如果柵格流入?yún)R水面積大于匯水面積閾值，則認為當前柵格為溝道。柵格流入?yún)R水面積為流入當前柵格的所有上游柵格面積之和匯流面積閾值由用戶輸入，其值的確定原則是讓由DEM生成的溝道盡可能地與實際情況相符。

4)計算非累計坡長L/m。即單個柵格對后續(xù)流入及流出坡長的貢獻，計算規(guī)則如下：

(1)如果沒有水流流入當前單元格，非累計坡長L為柵格大小Δx/m；

(2)如果入流方向與出流方向同為東北、西北、東南和西南方向，則L= 1.414 2Δx；

(3)如果入流方向與出流方向同為東、西、南和北方向，則L=Δx；

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

5)柵格入口/出口坡長計算。依據(jù)水流流向和中斷因子計算當前柵格的入口坡長。柵格入口坡長為流入柵格中出口坡長最大值；柵格出口坡長為柵格入口坡長加上與入流及出流方向對應的非累計坡長。根據(jù)坡長的定義，在坡長計算過程中用以下2個控制條件來確定流出坡長終點。

(1)溝道：當柵格是溝道時，則認為坡長終止，該柵格的流出坡長為0。

(2)中斷因子：中斷因子用于確定坡面坡度由陡坡到緩坡變化時，由于出現(xiàn)泥沙淤積而導致的坡長中斷。中斷因子定義為當前柵格的坡度值與其流入方向柵格坡度比值的臨界值，這一比值反映了當前柵格坡度相對于流入方向柵格坡度的減小程度：若當前柵格的坡度小于中斷因子與流入柵格坡度的乘積時，會出現(xiàn)泥沙淤積，坡長中斷，此時當前柵格的流出坡長重新設置為柵格大小。本次計算時，緩坡(坡度<3°)的中斷因子設為0.7，陡坡(坡度≥3°)的中斷因子設為0.5。

6)計算坡長坡度因子。依據(jù)坡度的不同分別采用式(1)～式(3)計算坡度因子；根據(jù)用戶輸入的坡長因子計算方法采用式(4)或式(5)計算坡長因子。

2.3 技術實現(xiàn)

軟件的開發(fā)工具為Visual Studio 2010，軟件算法由C++語言實現(xiàn)，軟件界面使用C#語言實現(xiàn)，文件讀寫用開源庫GDAL來實現(xiàn)，圖形顯示不依賴任何第三方圖形庫。

軟件安裝時需要Windows XP或以上操作系統(tǒng)，Microsoft.NET Framework 2.0或以上環(huán)境；LS計算工具支持32位和64位2個操作系統(tǒng)的計算機；32位版本能夠支持的數(shù)據(jù)一般不超過1萬行×1萬列的大小，我國大部分縣計算一次即完成LS因子的計算；64位版本能支持的數(shù)據(jù)大小由系統(tǒng)的物理內(nèi)存來決定。一個8 GB內(nèi)存的64位系統(tǒng)能夠處理2萬行×2萬列的浮點型柵格數(shù)據(jù)。一個32 GB內(nèi)存的64位系統(tǒng)可以處理4萬行×4萬列的浮點型柵格數(shù)據(jù)。

軟件設計考慮了用戶使用的方便性，界面友好、清楚，參數(shù)輸入界面如圖2所示。

圖2 LS計算工具軟件用戶界面Fig.2 User interface of the software for calculating LS factors

需要輸入的DEM文件格式為柵格文件*.aux、*.xml或*.tif。程序主要輸出以下變量(表1)。

表1 LS計算工具輸出文件

表1(續(xù))

3 應用案例

以青?，敹嗫h為試驗區(qū)域，計算了該地區(qū)的坡長坡度因子?，敹嗫h行政面積為2.5萬km2，海拔在3 915～5 262 m之間，平均海拔為4 392 m。平均坡度為8.2°，58%的區(qū)域坡度小于8°。計算時運用了30 m分辨率的DEM，該縣共有7 500×6 080個柵格，DEM文件大小為43 MB，僅運用不到2 min的時間完成計算，說明該程序具有很快的運算速度。計算的坡度、坡長、以及坡長坡度因子分別如圖3所示。統(tǒng)計參數(shù)如表2所示。

表2 主要輸出圖層統(tǒng)計參數(shù)

圖3 瑪多縣坡度、坡長以及坡度坡長因子分布圖Fig.3 Distribution of slope steepness, slopes length, slope steepness factor and slope length factor in Maduo County

4 結論

本文應用Visual Studio 2010平臺開發(fā)了LS因子計算工具，該軟件界面簡潔、清晰，便于用戶操作。計算結果可在軟件中直接查看，也可在Arc map軟件中進行查看。運用在32位計算機上可進行一般縣級區(qū)域的LS因子計算，在64位計算機上可運算更大區(qū)域的LS因子值。應用案例說明，該工具軟件具有很快的運算速度。該工具軟件的開發(fā)，對于區(qū)域土壤侵蝕評價和水土保持規(guī)劃提供了極為有效的工具。

[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

項目名稱：中央高?；究蒲袠I(yè)務費專項資金資助

符素華(1973—)，博士，教授，博士生導師。主要研究方向：水土流失與水土資源管理。E-mail: suhua@bnu.edu.cn

?通信作者簡介：劉寶元(1958—)，博士，教授，博士生導師。主要研究方向：土壤侵蝕。E-mail: baoyuan@bnu.edu.cn

S157.1

A

1672-3007(2015)05-0105-06