基于反地形坡長(zhǎng)地表匯流模擬逼近算法

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

（1西華師范大學(xué)生命科學(xué)學(xué)院，四川南充 637009；2西華師范大學(xué)地理科學(xué)學(xué)院，四川南充 637009；3西華師范大學(xué)四川省干旱河谷土壤侵蝕檢測(cè)與控制工程實(shí)驗(yàn)室，四川南充 637009；4中國(guó)科學(xué)院水利部、成都山地災(zāi)害與環(huán)境研究所，四川成都 610000）

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

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

為此，本研究以黃土高原甘泉和延川小流域?yàn)檠芯繀^(qū)，基于5 m分辨率數(shù)據(jù)，提出在反地形條件下不斷逼近真實(shí)坡長(zhǎng)結(jié)果的算法。該研究不僅降低了原有坡長(zhǎng)計(jì)算方法中由于填洼造成的誤差，也有效減輕了計(jì)算過程中的工作量，為坡長(zhǎng)計(jì)算提供了一定的思路和借鑒。

1 研究區(qū)概況

本研究選取黃土高原侵蝕劇烈、地形起伏大的甘泉和延川2個(gè)典型的丘陵溝壑地貌區(qū)作為研究對(duì)象[14-15]。研究區(qū)海拔為600～1 800 m，甘泉位于陜西省延安市中部，屬于陜西省北部黃土高原黃土丘陵，以黃土梁狀丘陵溝壑地貌為主，溝壑密度約5.6 km/100 hm2。延川則以梁峁?fàn)钋鹆甑孛矠橹鳎横蛊鸱?，殘塬、梁、峁相間分布，溝壑縱橫、河谷深切，溝壑密度在所選樣區(qū)中最大。研究區(qū)內(nèi)以大陸性季風(fēng)氣候?yàn)橹?，大陸性氣候特征明顯，年均降水量約447 mm，多短時(shí)強(qiáng)降雨，土壤侵蝕強(qiáng)烈，如圖1所示。

圖1 研究區(qū)地形

2 研究方法

2.1 研究原理

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

圖2 坡長(zhǎng)示意

本研究在反地形條件下計(jì)算坡長(zhǎng)時(shí)，為確定水流的起點(diǎn)，將原始地形分解為若干個(gè)集水盆地。當(dāng)反轉(zhuǎn)地形時(shí)，原始正向地貌的山地反轉(zhuǎn)成為負(fù)向流域及集水區(qū)系統(tǒng)，相應(yīng)的山頂點(diǎn)也演化為各集水區(qū)中的洼地點(diǎn)，原有正向地形上的洼地成為向上凸起的小丘。反地形水流的起點(diǎn)在原始地形中集水盆地的凹陷點(diǎn)上。該方法雖省略了填洼過程，但并不會(huì)對(duì)原始地形流路和相關(guān)數(shù)據(jù)的提取產(chǎn)生影響。這種計(jì)算過程不僅簡(jiǎn)單易懂，而且可以減小多種數(shù)據(jù)融合后計(jì)算過程中產(chǎn)生的誤差。

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

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

圖3 水文分析

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

利用網(wǎng)格計(jì)算器構(gòu)建地圖表達(dá)式α=arcsin（δh/M），并計(jì)算坡面與水平面之間的包含角α；使用三角函數(shù)cos構(gòu)建函數(shù)表達(dá)式L=M×cos（α），計(jì)算結(jié)果L即坡長(zhǎng)（圖4）。

圖4 坡長(zhǎng)計(jì)算示意

2.2 坡長(zhǎng)逼近算法驗(yàn)證

為驗(yàn)證ASL計(jì)算坡長(zhǎng)的精度，以相對(duì)差值系數(shù)和線性函數(shù)為描述指標(biāo)，對(duì)3種坡長(zhǎng)計(jì)算結(jié)果進(jìn)行分析。

2.2.1 LS_TOOL（LST）LS_TOOL 算法是基于GIS 技術(shù)和多方位地形系數(shù)的估算方法，用于獲取區(qū)域尺度的地形坡長(zhǎng)坡度參數(shù)等數(shù)據(jù)。以LS_TOOL 算法計(jì)算的坡長(zhǎng)（L）作為實(shí)際坡長(zhǎng)，驗(yàn)證本文提出的近似算法的可靠性[18]。

2.2.2 SAGA（SGL）SAGA軟件中的坡長(zhǎng)計(jì)算工具可以計(jì)算特殊情況下的坡長(zhǎng)。首先用Zevenbergen &Thorne 方法計(jì)算每個(gè)單元的斜率，然后使用D8單流向算法和閾值標(biāo)準(zhǔn)，將坡長(zhǎng)值向下累積。

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

為能夠更清楚地比較坡長(zhǎng)算法之間的相關(guān)關(guān)系，將其轉(zhuǎn)換為相對(duì)差異系數(shù)α，表示比較算法和基準(zhǔn)算法之間的偏差，其公式如下：

當(dāng)α=1 時(shí)，2 種計(jì)算方法的結(jié)果完全一致，沒有差異；α值越小，則算法之間的差異越大；當(dāng)α<0時(shí)，則2種算法基本沒有可比性。

3 結(jié)果與分析

3.1 不同區(qū)域3種坡長(zhǎng)算法

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

圖6 延川樣區(qū)坡長(zhǎng)（aASL、bLST、cSGL）

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

圖7 坡長(zhǎng)累積頻率分布（a甘泉坡長(zhǎng)結(jié)果頻率分布，b延川坡長(zhǎng)結(jié)果頻率分布）

3.2 相關(guān)性分析

為檢驗(yàn)ASL算法的可靠性，利用ArcGIS 10.2中的抽稀工具提取了3 種算法坡長(zhǎng)結(jié)果20%的數(shù)據(jù)，并分別計(jì)算ASL 與LST、ASL 與SGL、SGL 與LST 算法之間的相關(guān)關(guān)系。

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

表1 3種算法的相關(guān)性分析結(jié)果

算法之間的線性函數(shù)的R2值均在0.900 以上，其中最大的R2值出現(xiàn)在延川樣區(qū)的ASL 與LST 之間，為0.993；最小R2值出現(xiàn)在延川樣區(qū)的SGL 與LST之間，為0.937，表明3種算法的坡長(zhǎng)結(jié)果之間存在強(qiáng)相關(guān)性。

3.2.2 相對(duì)差異系數(shù) 利用抽稀得到的2個(gè)樣區(qū)20%坡長(zhǎng)數(shù)據(jù)計(jì)算其相對(duì)差異系數(shù)（a，b），結(jié)果見表2。

表2 3種算法的相對(duì)差異系數(shù)

從表2 可以看出，算法間的相對(duì)差異系數(shù)均在0.80 以上，說明每2 個(gè)樣區(qū)的坡長(zhǎng)結(jié)果高度相關(guān)。其中，甘泉區(qū)ASL 與LST 的相對(duì)差異系數(shù)最高，為0.96；延川樣區(qū)ASL與SGL的相對(duì)差異系數(shù)最低，為0.80。說明整體上這2 個(gè)樣區(qū)的坡長(zhǎng)結(jié)果相似度很高，但不同樣區(qū)算法間的相對(duì)差異系數(shù)不完全相同，表示算法之間存在一定的差異。

4 討論和結(jié)論

4.1 討論

影響坡長(zhǎng)計(jì)算結(jié)果并產(chǎn)生差異的原因，首先是地形表達(dá)的有限分辨率限制了地形高度的表達(dá)精度[21-22]。受網(wǎng)格大小的限制，一些山峰或山脊表達(dá)為平坦地形，導(dǎo)致集水區(qū)的劃定結(jié)果出現(xiàn)一定誤差。由圖8 中的等高線可以看出，地形變化比較明顯，頂部有大面積的平坦區(qū)域。結(jié)合谷歌圖片，可以看出坡長(zhǎng)出現(xiàn)較大值時(shí)存在大量的梯田、黃土脊?fàn)钌角鸬绕教沟孛差愋停匦蜗鄬?duì)破碎，這對(duì)DEM反映真實(shí)地形帶來了較大挑戰(zhàn)，進(jìn)而影響了坡長(zhǎng)的計(jì)算。如樣本區(qū)有較大的平坦區(qū)域，則會(huì)影響反地形下的分割精度，從而導(dǎo)致坡長(zhǎng)計(jì)算結(jié)果誤差較大[23]。本研究通過反地形計(jì)算坡長(zhǎng)，去除了在原始地形中填充洼地帶來的誤差，同時(shí)利用函數(shù)方法計(jì)算坡長(zhǎng)提高了坡長(zhǎng)計(jì)算的準(zhǔn)確度，降低了試驗(yàn)過程中可能帶來的誤差，同時(shí)也減輕了計(jì)算過程中的工作量。

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

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

4.2 結(jié)論

本研究基于5 m 分辨率高程數(shù)據(jù)、利用ASL、LST、SGL 3種坡長(zhǎng)計(jì)算方法對(duì)黃土高原丘陵溝壑區(qū)2 個(gè)樣區(qū)的坡長(zhǎng)進(jìn)行計(jì)算，并利用線性回歸和相對(duì)差異系數(shù)對(duì)計(jì)算結(jié)果進(jìn)行了相似性分析，得出以下結(jié)論。

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

坡長(zhǎng)因子作為估計(jì)土壤侵蝕狀況的重要參數(shù)，尋求一個(gè)不斷逼近真實(shí)坡長(zhǎng)結(jié)果的算法對(duì)提高土壤侵蝕治理水平十分重要。本研究提出的坡長(zhǎng)逼近算法為坡長(zhǎng)計(jì)算提供了更多的思路。