基于反地形坡長地表匯流模擬逼近算法

桑云云 羅明良 楊 慧 白雷超 姜佳艷 劉維明

（1西華師范大學生命科學學院，四川南充 637009；2西華師范大學地理科學學院，四川南充 637009；3西華師范大學四川省干旱河谷土壤侵蝕檢測與控制工程實驗室，四川南充 637009；4中國科學院水利部、成都山地災害與環境研究所，四川成都 610000）

我國土壤侵蝕分布范圍較廣，是土壤侵蝕較為嚴重的國家之一[1-2]。坡長是土壤侵蝕領域研究中的重要因子之一[3]，當其他條件相同時，土壤侵蝕的強度由坡面的長度來決定，坡面越長，匯聚的流量越大，侵蝕力和沖刷力就越強[4-5]。目前，關于坡度因子的計算已經相對成熟，而坡長因子的精確、快速提取是目前研究的重點之一[6-7]。

對于坡長的計算，學者提出了多種方法，例如：Hickey 等[8]在水文分析的基礎上提出了非累積流量的直接算法；Desmet 等[9]提出了基于累積流量的單位匯水面積算法；楊昕等[10]基于Arcview空間分析提出坡長快速近似算法。隨著計算機和GIS 研究方法的不斷創新，在DEM 數據和GIS 空間分析方法的支持下，坡長提取的方法得以不斷改進。其中，非累積流量的直接算法可以通過按照不同的流向進行賦值來體現坡長的準確特征[11]，但此方法需要大量的計算和復雜地形因子的提取；快速近似算法雖然計算過程簡便，但由于地形類型復雜，局部高程點無法準確識別，導致計算結果誤差較大[8]。基于網格DEM 的坡長計算可用來模擬地表的水流路徑，但要求水流在DEM 模擬的地形面上自由流動，因此地表的自然凹陷對坡長的計算結果有很大影響，需要在確定水流方向之前對DEM凹陷處進行填充和平整。此外，坡長的計算還需要考慮很多因素，不僅包括地形特征，還需要考慮具體地塊的土地利用類型[12-13]。坡長作為土壤侵蝕計算中不可忽視的因子，迫切需要一個高度精確的計算方法。

為此，本研究以黃土高原甘泉和延川小流域為研究區，基于5 m分辨率數據，提出在反地形條件下不斷逼近真實坡長結果的算法。該研究不僅降低了原有坡長計算方法中由于填洼造成的誤差，也有效減輕了計算過程中的工作量，為坡長計算提供了一定的思路和借鑒。

1 研究區概況

本研究選取黃土高原侵蝕劇烈、地形起伏大的甘泉和延川2個典型的丘陵溝壑地貌區作為研究對象[14-15]。研究區海拔為600～1 800 m，甘泉位于陜西省延安市中部，屬于陜西省北部黃土高原黃土丘陵，以黃土梁狀丘陵溝壑地貌為主，溝壑密度約5.6 km/100 hm2。延川則以梁峁狀丘陵地貌為主，梁峁起伏，殘塬、梁、峁相間分布，溝壑縱橫、河谷深切，溝壑密度在所選樣區中最大。研究區內以大陸性季風氣候為主，大陸性氣候特征明顯，年均降水量約447 mm，多短時強降雨，土壤侵蝕強烈，如圖1所示。

圖1 研究區地形

2 研究方法

2.1 研究原理

坡長一般定義為從坡面徑流的起點到徑流被攔截點或流路中斷點的水平距離[8，16]。如圖2 中A′B′是徑流線在水平面上的投影，即坡長。

圖2 坡長示意

本研究在反地形條件下計算坡長時，為確定水流的起點，將原始地形分解為若干個集水盆地。當反轉地形時，原始正向地貌的山地反轉成為負向流域及集水區系統，相應的山頂點也演化為各集水區中的洼地點，原有正向地形上的洼地成為向上凸起的小丘。反地形水流的起點在原始地形中集水盆地的凹陷點上。該方法雖省略了填洼過程，但并不會對原始地形流路和相關數據的提取產生影響。這種計算過程不僅簡單易懂，而且可以減小多種數據融合后計算過程中產生的誤差。

根據坡長的定義，如果將坡上的水流路徑（M）設為直角三角形的斜邊，則坡長可定義為坡上流線的長度乘以cos（α），其中α為流向與水平方向的夾角[17-19]。為計算α，坡長逼近算法（ASL）首先計算從水流流動起點到流路上某一點的高差，定義為δh，然后用反正弦函數求解α，即α=arcsin（δh/M）。按照這種思路，坡長逼近算法ASL的設計步驟如下。

反地形和相關水文分析（圖3），如流向、匯流累積量和分水嶺等[20]。其中，默認的流向算法為D8流向算法，即單元梯度的最佳代表值是以其為中心的3×3窗口中周圍8個方向的最大梯度，水流的流動方向與最大梯度相同。

圖3 水文分析

分區統計后利用柵格計算器計算前一步分區統計結果與原始DEM的差值，利用重分類工具將上一步柵格計算的結果進行賦值顯示后，復制為1 的區域即為區域最高點。接著使用水文分析中的水流長度工具基于所得到的流向數據計算出反地形下的水流長度（上坡長度：計算沿水流路徑從每個像元到分水嶺頂部的最長上坡距離），即為地形坡面長度。

利用網格計算器構建地圖表達式α=arcsin（δh/M），并計算坡面與水平面之間的包含角α；使用三角函數cos構建函數表達式L=M×cos（α），計算結果L即坡長（圖4）。

圖4 坡長計算示意

2.2 坡長逼近算法驗證

為驗證ASL計算坡長的精度，以相對差值系數和線性函數為描述指標，對3種坡長計算結果進行分析。

2.2.1 LS_TOOL（LST）LS_TOOL 算法是基于GIS 技術和多方位地形系數的估算方法，用于獲取區域尺度的地形坡長坡度參數等數據。以LS_TOOL 算法計算的坡長（L）作為實際坡長，驗證本文提出的近似算法的可靠性[18]。

2.2.2 SAGA（SGL）SAGA軟件中的坡長計算工具可以計算特殊情況下的坡長。首先用Zevenbergen &Thorne 方法計算每個單元的斜率，然后使用D8單流向算法和閾值標準，將坡長值向下累積。

2.2.3 相關性分析 利用線性回歸方法擬合出算法間的一元線性函數，當2 種算法之間的坡長值完全一致時，則函數關系為X=Y，XY 散點分布圖是一條斜率為1的直線；當2種算法之間的坡長值存在差異時，另一種函數關系為Y=a+bx，表明2種算法之間存在系統差異。接著計算相對差異系數，相對差異系數通過逐個柵格單元進行對比，可分區域分析算法之間的相似性。當2 種算法進行比較時，選擇其中一個作為基準算法，另外一個則作為比較算法，它們之間的差值θ如公式 （1）所示，可以通過減去相應位置的結果而得到（是基準算法在網格單元i的坡長為整個DEM的平均坡長，是比較算法在相應位置的坡長。）

為能夠更清楚地比較坡長算法之間的相關關系，將其轉換為相對差異系數α，表示比較算法和基準算法之間的偏差，其公式如下：

當α=1 時，2 種計算方法的結果完全一致，沒有差異；α值越小，則算法之間的差異越大；當α<0時，則2種算法基本沒有可比性。

3 結果與分析

3.1 不同區域3種坡長算法

ASL、LST和SGL 3種計算方法得到的甘泉和延川樣區的坡長結果如圖5～6 所示。從圖5～6 可以看出，同一樣區3種算法的坡長變化大致相同，但坡長結果均有差異，主要表現在3 種坡長算法的最大值上，ASL算法的最大值分別為564和468，LST算法的最大值為493 和353，SGL 算法的最大值為399 和389，ASL 算法相對于其他2 種算法較大的坡長值主要集中在山脊處。

圖6 延川樣區坡長（aASL、bLST、cSGL）

經過抽稀提取2 個研究區20%的坡長結果，并將結果劃分為3個部分繪制累積坡長（圖7）。從圖7可以看出，3種算法得到的坡長在0～300 m，占提取坡長總數的98%左右；而大于300 m的坡長較少，占2%左右。ASL受流向和網格分辨率的影響，在甘泉樣區有超過500 m的離群值，但其累計數量不大，僅占坡長總數的0.001 7%；而延川樣區超過400 m 的離群值占總數的0.002 4%，這些離群值均出現在山脊處。

圖7 坡長累積頻率分布（a甘泉坡長結果頻率分布，b延川坡長結果頻率分布）

3.2 相關性分析

為檢驗ASL算法的可靠性，利用ArcGIS 10.2中的抽稀工具提取了3 種算法坡長結果20%的數據，并分別計算ASL 與LST、ASL 與SGL、SGL 與LST 算法之間的相關關系。

3.2.1 線性回歸分析 每2種坡長算法間的一元線性函數以及R2值見表1。

表1 3種算法的相關性分析結果

算法之間的線性函數的R2值均在0.900 以上，其中最大的R2值出現在延川樣區的ASL 與LST 之間，為0.993；最小R2值出現在延川樣區的SGL 與LST之間，為0.937，表明3種算法的坡長結果之間存在強相關性。

3.2.2 相對差異系數 利用抽稀得到的2個樣區20%坡長數據計算其相對差異系數（a，b），結果見表2。

表2 3種算法的相對差異系數

從表2 可以看出，算法間的相對差異系數均在0.80 以上，說明每2 個樣區的坡長結果高度相關。其中，甘泉區ASL 與LST 的相對差異系數最高，為0.96；延川樣區ASL與SGL的相對差異系數最低，為0.80。說明整體上這2 個樣區的坡長結果相似度很高，但不同樣區算法間的相對差異系數不完全相同，表示算法之間存在一定的差異。

4 討論和結論

4.1 討論

影響坡長計算結果并產生差異的原因，首先是地形表達的有限分辨率限制了地形高度的表達精度[21-22]。受網格大小的限制，一些山峰或山脊表達為平坦地形，導致集水區的劃定結果出現一定誤差。由圖8 中的等高線可以看出，地形變化比較明顯，頂部有大面積的平坦區域。結合谷歌圖片，可以看出坡長出現較大值時存在大量的梯田、黃土脊狀山丘等平坦地貌類型，地形相對破碎，這對DEM反映真實地形帶來了較大挑戰，進而影響了坡長的計算。如樣本區有較大的平坦區域，則會影響反地形下的分割精度，從而導致坡長計算結果誤差較大[23]。本研究通過反地形計算坡長，去除了在原始地形中填充洼地帶來的誤差，同時利用函數方法計算坡長提高了坡長計算的準確度，降低了試驗過程中可能帶來的誤差，同時也減輕了計算過程中的工作量。

圖8 研究區坡度和等高線（a甘泉，b延川）

其次，D8流向算法的局限性可能會影響計算結果。流向是8 個可能的網格方向之一，它們之間的間隔為45°，流向是確定的且不連續的。因此，D8算法具有高度敏感性，即較小的高度誤差也會導致水流方向的改變，進而導致坡長計算結果出現偏差。其他2種坡長算法同樣采用D8流向算法，因此不同流向算法對坡長的具體影響，還需要在后續的研究中做進一步的分析。

4.2 結論

本研究基于5 m 分辨率高程數據、利用ASL、LST、SGL 3種坡長計算方法對黃土高原丘陵溝壑區2 個樣區的坡長進行計算，并利用線性回歸和相對差異系數對計算結果進行了相似性分析，得出以下結論。

相同精度的DEM數據得到2個樣區內不同計算方法結果具有差異。本研究提出的逼近算法其坡長結果98%的數據主要集中在0～300 m，并且線性回歸以及相對差異系數分析結果表明3種算法計算得到的坡長結果非常相似。3個坡長計算方法通過回歸函數得到的R2值均在0.9 以上，同時2 個樣區的3個算法間的相對差異系數也主要集中在0.9附近，表明3 種坡長計算結果擬合效果均較好。因此，為獲取更高精度的坡長結果，可將3 種算法的坡長平均值作為最逼近真實的坡長用于分析地區的土壤侵蝕狀況。

坡長因子作為估計土壤侵蝕狀況的重要參數，尋求一個不斷逼近真實坡長結果的算法對提高土壤侵蝕治理水平十分重要。本研究提出的坡長逼近算法為坡長計算提供了更多的思路。