孔鋒　王一飛　呂麗莉　閆緒嫻

摘要 采用基于GIS的非累計流量的坡長直接提取算法（NCSL）和空間分析提取法（SAC）對陜北省安塞縣10個樣區(qū)的坡長進行提取，并計算相應的坡長因子。對不同DEM分辨率下的坡長提取結果對比，結果表明：NCSL對坡長的提取精度明顯好于SAC的提取結果，其中DEM在5 、10 m分辨率下的提取精度最好，且二者計算的坡長因子值差異不大，因此，可采用10 m分辨率NCSL方法提取坡長。

關鍵詞 DEM；分辨率；坡長；精度差異；土壤侵蝕模型

中圖分類號 S157 文獻標識碼 A 文章編號 0517-6611（2017）29-0104-05

Comparison of Slope Length Factor Extraction in Hillslope Soil Erosion Model with Different DEM Resolution

KONG Feng1，2，3，4，5， WANG Yifei1， L Lili1，2 et al

（1.China Meteorological Administration Training Center， Beijing 100081；2. China Meteorological Administration Development Research Center， Beijing 100081；3.State Key Laboratory of Earth Surface Processes and Resource Ecology， Beijing Normal University， Beijing 100875；4.Academy of Disaster Reduction and Emergency Management， Ministry of Civil Affairs & Ministry of Education， Beijing 100875；5. College of Management Science and Engineering， Shanxi University of Finance&Economics，Taiyuan， Shanxi 030006 ）

Abstract In this research，noncumulative slope length （NCSL） calculation method and spatial analytical calculation （SAC） method were respectively applied to extract slope length and slopelength factor from 10 sample areas， which were located in Ansai County， north Shaanxi Province. The comparison of computation precision among variable DEM resolutions showed that NCSL was superior to SAC entirely. And the results were best when DEM resolutions was 5m and 10 m；besides，the result of slopelength factor was nearly the same under the same condition. So DEM of 10 m resolution can be used to extract slopelength.

Key words DEM；Resolution；Slope length；Precision differentiation；Soil erosion model

基金項目 中國氣象局氣象軟科學重點項目“基于綜合風險防范視角的中國及周邊國家安全和全球戰(zhàn)略研究”（2017〔21〕）；中國氣象局氣象軟科學自主項目“新常態(tài)下中國自然災害風險時空格局和綜合防災減災工作的現(xiàn)狀、趨勢、挑戰(zhàn)及戰(zhàn)略對策范式研究”（2017〔36〕）；中國氣象局氣象軟科學自主項目“中國氣象災害防御能力評估及政策建議”（2017〔35〕）。

作者簡介 孔鋒（1986—），男，山西臨汾人，助理研究員，博士，從事自然災害與環(huán)境演變研究。

收稿日期 2017-08-02

快速城鎮(zhèn)化背景下人類的不合理生產(chǎn)活動使土壤侵蝕的強度日漸加劇，由此引發(fā)嚴重的水土流失，不僅造成全球范圍內(nèi)的土壤肥力下降，土地生產(chǎn)力喪失，而且導致諸多江河的河床抬高，水庫庫容減少以至報廢，進而加劇洪澇等災害的發(fā)生。土壤侵蝕已經(jīng)成為人類生存與可持續(xù)發(fā)展的世界性重大環(huán)境問題之一。

1 土壤侵蝕模型與坡長提取

1.1 土壤侵蝕模型進展

為了更好地掌握土壤侵蝕發(fā)生的規(guī)律和機理，減小因嚴重的水土流失而帶來的災難和損失，世界各國的科學家相繼在該領域開展了大量的研究工作。建立定量的模型是研究土壤侵蝕的重要手段之一。 Wischmeier等[1]于1959年提出了以下通用土壤流失方程USLE（universal soil loss equation）：

A=R×K×LS×C×P （1）

式中，各因子分別代表降雨侵蝕力、土壤可蝕性、坡長、坡度、覆蓋和管理以及水土保持措施。該公式發(fā)表后，世界各國學者又不斷地修訂并改進上述方程中的各因子，這些改進成果被統(tǒng)稱為RUSLE（revised universal soil loss equation）[2-3]?；谏鲜鏊悸?，世界許多國家相繼建立了更為適合本國實際的土壤侵蝕模型。例如我國學者劉寶元等[4]曾提出中國土壤流失預報方程，將上述公式中的覆蓋和管理因子以及水土保持措施因子變?yōu)樗帘３稚锎胧˙）、工程措施（E）與耕作措施（T）3個因子：

A=R×K×LS×B×E×T （2）

RUSLE是目前應用最多的水蝕預報經(jīng)驗模型。近年來，該領域又有新的進展，特別值得一提的是基于土壤侵蝕機理過程（水蝕和風蝕等）的物理模型的出現(xiàn)，開辟了土壤侵蝕定量研究的新紀元。

坡長指從地表徑流原點到坡度減小直至有沉積出現(xiàn)的地方之間的距離，或到一個明顯渠道之間的水平距離，這是關于完整坡的定義[1-3]。然而在絕大多數(shù)情況下，實際的坡面形態(tài)是非常復雜的，故上述完整坡長的定義不再適用。為此，F(xiàn)oster等[3]提出分段處理不規(guī)則坡面的方法，該方法視每一段坡長為上游各個分段坡長的累加，我國學者湯國安等[5-7]據(jù)此將坡長定義為地面上一點沿水流方向到其源點的最大距離在水平面上的投影。該定義是基于DEM坡長自動提取的理論基礎。在上述通用土壤流失方程USLE和改進后的通用土壤流失方程RUSLE中，坡長是以坡長因子的形式出現(xiàn)的，用來描述坡長與土壤侵蝕量的關系。在RUSLE中，坡長因子的計算公式為：

L=（λ/22.13）m （3）

式中，λ為實際坡面的水平投影距離；22.13是標準化坡長的水平投影長度，其單位是m；m為可變的坡長因子指數(shù)，其取值與細溝侵蝕與細溝間侵蝕的比例有關。對于不規(guī)則坡和分段坡，人們作分段處理，將每個分段內(nèi)的地形特征看作是均勻的，則從坡頂開始第i分段對應的坡長因子計算公式為：

Li= λm+1i-λm+1i-1 （λi-λi-1）22.13m （4）

式中，λi為第i段坡底到上游邊界的長度。這種方法考慮坡面的復雜性，為坡長因子量化提供了依據(jù)。如果考慮到水流的二維性，坡面上某點的徑流量和土壤侵蝕量并非完全取決于坡長值，而取決于向該點匯水的單位匯水面積，這一觀點更能體現(xiàn)徑流的分散和匯集，其水文學意義更明顯。接著，Desmet等[8]提出了在二維情況下土壤侵蝕模型中的坡長因子應由坡長長度轉(zhuǎn)為單位匯水面積，例如單位等高線長度的上游匯水面積?；谏鲜鏊悸?，Kinnell[9]對USLE中的地形因子的計算方法進行了改進，得出了一種基于DEM自動提取的新算法。該算法用單位匯水面積A代替坡長λ：

Li= Am+1i，j-out-Am+1i，j-in （Ai，j-out-Ai，j-in）22.13m （5）

式中，Ai，j-in表示流入i，j格網(wǎng)中的上游匯水面積，Ai，j-out代表從該網(wǎng)格流出的匯水面積，后者由于包含網(wǎng)格本身的面積而大于前者，上式又可改為：

Li，j= （Ai，j-in+g2）m+1-Am+1i，j-in g2xi，j（22.13）m （6）

式中，g為格網(wǎng)分辨率，xi，j為該格網(wǎng)相對應的有效等高線長度系數(shù)。對于RUSLE[1-3]，坡長因子計算公式為：L=（m+1）（λ/22.13）m。轉(zhuǎn)化為單位匯水面積形式為：

Li，j=（m+1） 2Ai，j-in+g2 2gxi，j（22.13） m （7）

在實際應用中，人們采用GRASS-GIS進行坡長因子的計算[10-12]。國內(nèi)有人曾經(jīng)使用該算法求得坡長因子，然后帶入公式（1）中反演實際坡長，結果表明這種方法只能整體上反映區(qū)域內(nèi)的坡長分異規(guī)律，而坡長值的提取結果精度較低[13-16]。

1.2 4種坡長提取算法對比

基于GIS的地形因子提取多研究坡度、地形等因子，而對其在坡長提取中的應用研究尚不很成熟。目前，非累計流量的直接計算方法（noncumulative slope length，NCSL）、基于累計流量的單位匯水面積計算法（specific contribution area，SCA）、

基于水流強度指數(shù)的間接計算法（stream power index）、空間提取算法是比較常用的、也是相對較成熟的基于GIS的坡長自動提取算法，被廣泛應用于水利、國土、交通、環(huán)保和農(nóng)業(yè)等部門。

1.2.1 非累計流量的直接計算方法（noncumulative slope length，NCSL）。該算法主要基于“水流方向與最大坡降的方向一致”這一假設。其主要原理是利用“八聯(lián)通”法則，水流由源頭（如山脊線，山峰等）像素出發(fā)，朝著最大坡降方向像素流動，依次類推， 得出水流路徑，進而求得坡長。該算法的特點是算法簡單，效率較高，尤其是對洼地和平坦區(qū)域有較強的處理能力。但是該算法認為水流方向與最大坡降一致，實地試驗表明這與現(xiàn)實的地形和水流情況并非完全吻合[15]。

1.2.2

基于累計流量的單位匯水面積計算法（specific contribution area，SCA）。該算法認為水流分布具有分散性質(zhì)，其基本思想是用多流向算法得到單位匯水面積，進而取代方程中的坡長因子。由于匯水面積的計算是通過流量累積的方式得到，故稱之為基于流量累積的坡長計算方法。該算法解決了方法1中的水流流向問題，但是算法設計極為復雜，并且要考慮諸多特殊情況，因此雖然該算法能夠精確地算出匯水面積，但是程序復雜且運行可靠性不高而應用受限[14]。

1.2.3

基于水流強度指數(shù)的間接計算法（stream power index）。該算法是基于流量指數(shù)的坡長，非直接物理意義上的坡長，而是計算坡長坡度的合成因子（USLE中的LS因子），并認為LS因子是地表徑流輸沙能力的度量，從而將代表地面曲面形態(tài)的LS因子的計算解譯為流量和坡度呈非線性函數(shù)關系的無量綱輸沙能力的指數(shù)的計算。值得注意的是，基于水流強度指數(shù)的間接計算法得到的坡長其意義已經(jīng)不是地理學中的物理坡長。該算法較少應用于地理學領域，而在水土保持和水文學方面利用較多[16]。

1.2.4

空間提取算法。除上述需編程實現(xiàn)的算法之外，GIS軟件的水文分析功能有直接的坡長提取算法，我們稱之為空間分析提取法，其大致原理是如果坡地近似滿足坡面水流方向與山脊線垂直，則只要計算出每個點沿垂直方向到山脊線的水平距離可以作為該點坡長的近似值，相對以上3種提取方法而言，該算法較為簡單，其大致原理是先采用負地形法提取出研究區(qū)域的山脊線，提取出山脊線以后，求取格網(wǎng)單元到最近鄰的山脊線的垂直距離作為近似坡長。由于該方法操作簡便，計算快捷，又稱為快速計算方法?？焖儆嬎惴椒ǖ膬?yōu)點是計算快捷，操作簡單。缺點是適用范圍比較窄——僅僅局限在等高線與山脊大致平行的地貌（例如川和梁等）或者坡面較均勻、坡度較緩的坡地地貌，計算精度也相對不是很高[6，15]。

綜上可知，坡度、坡長和坡向是坡面的3個最基本的地形要素。然而相對于諸如坡度、坡向、溝壑密度、坡面曲率等其他地形要素，人們對坡長的定義仍然有很大爭議；有關坡長及坡長因子的研究還不夠深入。已有的不同提取算法也有各自的優(yōu)劣，而對于這些優(yōu)缺點和方法的適應性并未給出明確的限定和詳細的討論，生成的坡長往往誤差較大，精度也不好控制。此外，對于基于GIS的提取算法，在DEM（數(shù)字高程模型）高分辨率（<20 m）條件下尚無研究，因DEM分辨率不同而造成的提取結果精度的差異。以足夠大比例尺地形圖為源數(shù)據(jù)并數(shù)字化，利用功能強大的GIS軟件可以很方便地提取到足夠精確的DEM。目前已知從1 ∶10 000比例尺地形圖中提取得到的DEM的理想分辨率為5 m[5-7]。在源數(shù)據(jù)精度得以保證的情況下，DEM的分辨率越高，與真實地形的吻合程度也越高，相應地形要素的提取效果也越好。但是利用程序從高分辨率DEM中提取坡長意味著高耗時、低效率，這一缺陷在處理海量數(shù)據(jù)的時候體現(xiàn)得尤為明顯。在DEM高分辨率條件下，其高低對坡長的提取和坡長因子的計算精度有無顯著性影響目前尚無定論。

2 研究區(qū)概況與方法

2.1 研究區(qū)概況

選擇的試驗區(qū)位于陜西省延安市安塞縣境內(nèi)（圖1），該縣地處陜北黃土高原核心區(qū)域，鄂爾多斯盆地南緣，黃河一級支流的延河流域，縣域面積2 950 km2。區(qū)域地貌類型復雜多樣，溝壑密度達4.7 km/km2，平均海拔1 371.9 m，整個區(qū)域均被黃土覆蓋，厚度大，基本為全新世黃土。氣候?qū)僦袦貛Т箨懶园敫珊导撅L氣候，年均氣溫9.0 ℃，年均降水量約500 mm。水土流失嚴重，長期觀測數(shù)據(jù)表明，該區(qū)域年平均侵蝕輸沙模數(shù)高達9 370 t/（km2·a），自然條件惡劣，生態(tài)環(huán)境惡化。

2.2 數(shù)據(jù)來源

數(shù)據(jù)主要涉及從1 ∶10 000大比例尺地形圖上進行試驗小區(qū)的選擇和DEM的生成。試驗小區(qū)選擇的好壞直接關系到研究結果的可靠性。基于科學性原則、代表性原則和信息全面原則，同時考慮該研究的實際情況，選取典型流域小區(qū)共10塊，這些地塊的覆蓋面積介于0.5～1.0 km2，以這些小區(qū)的1 ∶10 000比例尺的矢量地形圖作為源數(shù)據(jù)（表1）。以1 ∶10 000地形圖為源數(shù)據(jù)所生成的DEM的理想分辨率為5 m，是分辨率的上限。楊勤科等[17]在榆林、延安和秦嶺地區(qū)的3個樣區(qū)基于10、50 m分辨率DEM提取的坡長表明，低分辨率下的提取結果是高分辨率下提取結果的近4倍，由此可見，對于坡長提取尤其是小流域的坡長提取來說，50 m分辨率設置實在太低，提取結果必然誤差很大。為探討DEM分辨率差異帶來的計算結果的差異，該研究以5 m分辨率為上限，20 m分辨率為下限，每隔5 m選取一個分辨率值，共選擇4組。對每一塊試驗小區(qū)均生成4組由高到低不同分辨率的DEM，然后利用算法提取坡長，分析結果。

2.3 研究方法

利用2種不同的既有坡長提取算法，以若干具有代表性的小流域試驗區(qū)的不同分辨率DEM為源數(shù)據(jù)，計算相應小區(qū)的相應分辨率的坡長和坡長因子，然后采用數(shù)理統(tǒng)計分析法進行結果數(shù)據(jù)處理，得出相應指標，確定DEM分辨率與提取精度的關系，并分析產(chǎn)生誤差的原因?？紤]到算法實現(xiàn)的復雜性和論文篇幅限制，采用非累計流量的直接計算方法（noncumulative slope length，NCSL）和空間分析提取法（快速計算方法）提取坡長。如果DEM分辨率高低對提取結果的精度沒有顯著性影響，則在保證精度的情況下可以通過降低DEM分辨率來提高計算效率，這在處理海量數(shù)據(jù)時是有一定意義的。采取的統(tǒng)計方法主要是常規(guī)數(shù)理統(tǒng)計分析法。主要包括以下4種：①絕對誤差和相對誤差：提取值和標準坡長值的差值的絕對值為絕對誤差；絕對誤差與標準坡長值的比值乘以100%為相對誤差。②合格率：設30%為相對誤差上限，當提取得到的坡長值與標準值的誤差小于30%為合格點，反之為不合格。以合格點數(shù)目與總數(shù)的比值乘以100%即得到合格率。③相關系數(shù)。④確定性系數(shù)：該系數(shù)由Nash和Sutcliffe于1970年提出，用以分析提取值與實際值的相關程度。該系數(shù)反映的是計算值與實測值1 ∶1線的接近程度，其取值范圍是（-∞，1）。當確定性系數(shù)大于0時，表明提取值與標準值相當吻合；當確定性系數(shù)越接近1，提取值越精確；當確定性系數(shù)為0時，表示計算結果與實測值之差的平方和恰好等于實測值與其均值的平方和，表明用實測值的平均值來估算與計算值的效果一樣。此外，采用單因素方差分析法以確定DEM分辨率是否對坡長因子計算結果有顯著性影響。

3 結果與分析

3.1 2種提取方法比較

由于在自動選擇100 m見方的中心樣點時有很大的隨機性，故可能會出現(xiàn)樣點落在山脊線或山頂上，以及有些點落在溝谷的底部的情況。首先應該對參與分析的樣點進行篩選，剔除那些不宜與坡長真值作比較的樣點，最終確定參與比較的樣點的數(shù)目為218個（原始的樣點數(shù)目為292個）。統(tǒng)計結果（表2、3）顯示，2種算法的提取精度均隨著DEM分辨率的降低而降低，說明DEM分辨率越高，坡長的提取效果越好；非累計流量的直接計算法的坡長提取效果遠好于空間分析提取法。

3.2 DEM分辨率對坡長因子的影響

通過單因素方差分析確定DEM分辨率對坡長因子的計算結果是否有顯著影響。當DEM分辨率為5、10 m時，非累計流量的直接計算法提取的坡長精度較為理想，均為80%以上，而且二者精度相差不明顯；但當DEM分辨率為15、20 m時，提取到的精度已經(jīng)不夠理想。坡長決定坡長因子，坡長因子的計算精度依賴于坡長的提取精度，故當DEM的分辨率大于10 m時，坡長因子的提取精度必然也不高。因此，只需要確定在DEM分辨率為5、10 m這2種條件下對坡長因子的提取有無顯著影響。由于在選擇樣區(qū)和樣點時均考慮隨機性和典型性，可以將參與統(tǒng)計分析的218個樣點的坡長因子值看作是來自正態(tài)總體的一個樣本?？紤]DEM分辨率分別為5、10 m的2個水平，取顯著性水平α=0.05，檢驗在這2種水平下坡長因子的均值有無顯著性差異，進而可以確定當DEM的分辨率

高于10 m時分辨率對坡長因子的計算效果有無顯著影響。

建立給定問題的原假設和備擇假設：H0：μ1=μ2，H1：μ1≠μ2（μ1、μ2分別為2組數(shù)據(jù)的期望值）。F（0.02）α（0.05）（表4）。故接受H0，認為2組坡長計算結果沒有顯著性差異。經(jīng)過上述分析可以得出，從5 m分辨率的DEM和10 m分辨率的DEM中提取到的坡長因子沒有明顯差別，也就是說當DEM分辨率高于10 m時分辨率對坡長因子的計算效果沒有顯著性影響。

4 結論與討論

4.1 結論

（1）非累計流量的直接計算法提取出的坡長的效果整體上遠好于空間分析法的提取效果。之所以如此，是因為非累計流量算法有明確的物理意義，其設計思路有嚴密的水文地貌學依據(jù)，并且依照分段坡的計算方法計算坡長，與真實的地形特征較為接近。此外，當DEM的分辨率足夠大（高于10 m）時，非累計流量直接計算法的提取精度已經(jīng)相當理想。所以，非累計流量直接計算法已經(jīng)足夠成熟，能夠在坡長提取中推廣應用。

（2）當DEM的分辨率分別為5、10 m時，利用非累計流量的直接計算法提取的坡長和坡長因子的精度均沒有顯著性差異。可見，在提取坡長和坡長因子時，一般情況下將DEM的分辨率設置為10 m即可，沒有必要設置為5 m。降低DEM的分辨率很大程度上意味著提高程序的計算效率。

（3）空間分析計算法的整體精度較差?？臻g分析法提取精確的前提條件是山脊線與等高線基本平行，這顯然與一般的實際地形不符甚至差別很大；還有該算法在提取山脊線時閾值的設置沒有公認的可靠標準，存在盲目性和隨意性，極易導致山脊線的誤判和漏提， 故造成計算精度很低，尤其是在以低分辨率DEM為源數(shù)據(jù)時計算結果誤差極大，無實際價值，表明該算法在坡長提取方面的應用尚不成熟，亟待改進。

4.2 討論

首先，受時間等客觀條件的制約，試驗小區(qū)的選擇數(shù)量偏少。要想得到更為可靠和令人信服的結論，按照以往研究的慣例，僅安塞縣境內(nèi)就需要選擇試驗小區(qū)近100個。在試驗小區(qū)數(shù)量少的情況下，很難反映出該區(qū)域的地域分異規(guī)律。其次，由于該研究水平十分有限，僅從已有的4種常用算法中選擇了2種較為簡單的進行提取效果的對比，未能開闊思路，發(fā)掘出精度和效率更高的、更適合地理學坡長提取的算法。此外，沒有系統(tǒng)地分析誤差，未能給出造成計算誤差的原因，更沒有提出減少和校正誤差的可行的辦法。這些都是很值得進一步研究的內(nèi)容。最后，算法的精度評價是一個精細而又復雜的過程，尤其對坡長這一原本就有爭議的地形要素來說，其精度評價必然存在著很多尚未解決的問題，涉及這些問題的難度遠遠超出了該研究的實際水平。此外，整個研究中從地形圖的數(shù)字化到坡長真值的量取，均不可避免地存在誤差；設置的評價標準也較為單一。在這樣的條件下得出的結論是否真正可靠有待進一步研究。

參考文獻

[1] WISCHMEIER W H.A rainfall erosion index for a universal soil-loss equation[J].Soil science society of America journal，1959，23（3）：246-249.

[2] ONORI F，BONIS P D，GRAUSO S.Soil erosion prediction at the basin scale using the revised universal soil loss equation （RUSLE） in a catchment of Sicily （southern Italy）[J].Environmental geology，2006，50（8）：1129-1140.

[3] FOSTER G R，MCCOOL D K，RENARD K G，et al.Conversion of the universal soil loss equation to SI metric units[J].Journal of soil & water conservation，1981，36（6）：355-359.

[4] 劉寶元，謝云，張科利，等.土壤侵蝕預報模型[M].北京：中國科學技術出版社，2001.

[5] 湯國安，陳正江，趙牡丹，等.Arc View地理信息系統(tǒng)空間分析方法[M].北京：科學出版社，2002.

[6] 湯國安，劉學軍，房亮，等.DEM及數(shù)字地形分析中尺度問題研究綜述[J].武漢大學學報（信息科學版），2006，31（12）：1059-1066.

[7] 湯國安，楊昕.ArcGIS地理信息系統(tǒng)空間分析實驗教程[M].北京：科學出版社，2006.

[8] DESMET P J J，GOVERS G.GISbased simulation of erosion and deposition patterns in an agricultural landscape：A comparison of model results with soil map information[J].Catena，1995，25（1/2/3/4）：389-401.

[9] KINNELL P I A.The effect of slope length on sediment concentrations associated with side-slope erosion[J].Soil science society of America journal，2000，64（3）：1004-1008.

[10] 李俊.基于DEM的黃土高原坡長的自動提取和分析[D].西安：西北大學，2007.

[11] HOYOS N.Spatial modeling of soil erosion potential in a tropical watershed of the Colombian Andes [J].Catena，2005，63（1）：85-108.

[12] DE VENTE J，POESEN J，VERSTRAETEN G，et al.Spatially distributed modelling of soil erosion and sediment yield at regional scales in Spain [J].Global and planetary change，2008，60（3/4）：393-415.

[13] 周貴云，劉瑜，鄔倫.基于數(shù)字高程模型的水系提取算法[J].地理學與國土研究，2000，16（4）：77-81.

[14] 蔡強國，劉紀根.關于我國土壤侵蝕模型研究進展[J].地理科學進展，2003，22（3）：242-250.

[15] 曹龍熹，符素華.基于DEM的坡長計算方法比較分析[J].水土保持通報，2007，27（5）：58-62.

[16] 張啟旺，安俊珍，王霞，等.中國土壤侵蝕相關模型研究進展[J].中國水土保持，2014（1）：43-46.

[17] 楊勤科，趙牡丹，劉詠梅，等.DEM與區(qū)域土壤侵蝕地形因子研究[J].地理信息世界，2009，7（1）：25-31.