李亞平,盧小平,張 航，路澤忠，王舜瑤

(河南理工大學測繪與國土信息工程學院，焦作 454003)

0 引言

土壤侵蝕是世界三大環(huán)境問題之一，會造成土壤有機質流失、泥沙淤積和土地退化等，嚴重威脅區(qū)域的可持續(xù)發(fā)展[1]。隨著地理信息系統(geographic information system,GIS)和遙感(remote sensing,RS)技術的日益成熟，結合RS/GIS技術和修正通用土壤流失方程(revised university soil loss equation,RUSLE)的區(qū)域土壤侵蝕定量評估，顯著提高了監(jiān)測及預測效率和空間可視化能力[2,3]，在世界中小流域或區(qū)域的土壤流失研究中得到了廣泛應用[4]。

坡長L因子是RUSLE模型中的重要影響因子，而實際研究中大尺度的坡長很難通過實測獲取，且利用DEM提取坡長的方法存在許多復雜性和不確定性，使L因子成為模型應用受到限制的主要因素[5,6]。Mahala[7]和Yue Li等學者[8]的眾多研究都以規(guī)則格網的單元格代表坡面的一個坡段，忽略了單元格間的匯流過程，使其形成水文孤島[9]。Desmet等[10]基于匯水面積的L因子提取方法得到了廣泛應用。雍斌等[11]基于多流向算法對累積分配中單位等高線匯流面積進行了改進。然而這些研究都忽略了上坡地表覆蓋/植被對匯流的影響，植被和地形的耦合作用必然會導致匯流變化，進而對侵蝕產生影響[12]。秦偉等[13]考慮了上坡土地利用/覆蓋，但其實驗僅針對黃土高原地區(qū)，且缺乏對比實驗與驗證，對多流向上坡累計匯流算法的介紹也不夠詳盡。針對現有研究中L因子提取算法不科學和對上坡匯流面積計算不準確等問題，本文改進了RUSLE模型中地形因子的提取算法，提出了考慮地表土地利用/植被對匯流影響的基于多流向的上坡匯流累積的L因子提取算法，并以淮河流域的商城縣為例，選擇適合該區(qū)域的其他因子算法獲取土壤侵蝕量及其分布特征，通過對比證明了本文方法的有效性。

1 地形因子提取方法

1.1 通用算法

通用的RUSLE模型中提取地形因子的計算公式為

(1)

L=(λ/22.1)m，

(2)

(3)

式中：L和S分別為坡長和坡度因子;θ為利用數字高程模型(digital elevation model，DEM)數據提取的坡度值，°；m為坡長指數；λ為坡長，m。

通用的RUSLE模型提取的L因子忽略了上坡匯流作用，導致L因子小于真實值，影響土壤侵蝕評定的準確性。本文采用基于上坡匯流的L因子算法對L因子的提取進行了改進。

1.2 本文提出的改進算法

基于上坡匯流的L因子算法首先要計算上坡累積匯流面積。ArcGIS軟件中匯流累積計算采用的是單流向算法，而真實的水流方向并不唯一，多流向算法更加符合真實三維地形中的匯流狀況[14]。本研究采用Quinn等[15]提出的多流向算法，使用3×3的柵格窗口計算中心像元累積的上坡匯流面積，通過下坡方向加權后的有效等高線長度分配到相鄰八像元中下坡方向像元的流量，權重由下坡方向決定，如圖1。

圖1 多流向徑流分配方法Fig.1 Flow partition using the multiple flow direction algorithm

各下坡方向的累積匯流面積公式為

(4)

式中：Ai為第i個相鄰下坡單元的累積匯流面積；A為總上坡面積；tanβi為第i個下坡方向的坡角，βi為坡度；Wi為第i個與流向垂直的有效等高線長，為柵格尺寸與權重的積(主方向上權重為0.5，對角方向為0.354)；k為下坡單元個數。

Griffin等[10]在考慮匯流累積及發(fā)散后，提出匯流計算L因子的算法，即

(5)

坐標(i，j)處L因子可以表示為

(6)

式中ASi,j-out和ASi,j-in分別為坐標(i，j)出水口和入水口單位匯水面積。

多流向的匯流累積算法只反映了地形對匯流的影響，忽略了地表植被對產流的影響，而事實上植被對降雨再分配和阻滯地表徑流等影響顯著。雖然RUSLE中的植被覆蓋與管理因子C反映了植被覆蓋對侵蝕的影響，但C因子僅代表所在單元格的影響，并無法體現上坡土地利用和植被對匯流的作用[13]。因此本文使用不同土地利用類型的匯流面積貢獻率tn更加真實地反映上坡植被和地形的作用。首先，根據淮河流域有關研究成果得到不同土地利用類型的平均產流系數[16]，考慮到RUSLE模型是以休閑耕地為基準建立的，tn以裸地的產流系數為標準換算得到(表1)，得到單元上坡實際匯水面積，即

表1 不同土地利用匯流面積貢獻率Tab.1 Contribution rate of confluence area for different land uses

(7)

Di,j=R·(sinθi,j+cosθi,j)，

(8)

式中：Ai,j-out和Ai,j-in分別為(i，j)處出水口和入水口的匯水面積；R為柵格寬度；Di,j為所在柵格的有效等值線長度；θi,j為坡向。將(7)代入(6)中得到

(9)

算法實現首先要對DEM進行填洼處理；其次要注意坡長的截斷位置，本文選擇水體作為坡長L因子的計算終點；最后使用MATLAB編程完成算法。技術路線如圖2。

圖2 L因子提取流程圖Fig.2 Flow chart for extraction of L factor

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

2.1 研究區(qū)概況

淮河流域大別山區(qū)是集山區(qū)、貧困區(qū)、水土流失嚴重區(qū)于一體的特殊地區(qū)。水土流失致使區(qū)域生態(tài)環(huán)境惡化，嚴重阻礙了經濟社會發(fā)展[17]。本文選擇淮河流域大別山區(qū)內的信陽市商城縣作為研究區(qū)。該區(qū)域位于大別山北麓，地理位置E115°06′～115°37′，N31°23′～32°05′，土地總面積為2 109.66 km2，流域面積占全縣總面積的60%，海拔高程44.5～1 584 m，地勢由南向北傾斜，形成中低山、低山丘陵和丘陵壟崗3大自然區(qū)。

2.2 數據收集

研究使用的數據包括：①Landsat衛(wèi)星2010年空間分辨率為30 m的TM遙感影像；②2010年空間分辨率為2.5 m的SPOT5數據；③空間分辨率為10 m的DEM數據；④淮河流域內及周邊共58個氣象站點1981—2010年間的月降雨數據及年降雨數據，來源于中國氣象局氣象中心；⑤中國土壤數據集來源于寒區(qū)旱區(qū)科學數據中心。

2.3 RUSLE模型

RUSLE模型充分考慮了影響土壤侵蝕的自然要素，具有較強的實用性。其表達式為

A=R·K·LS·C·P，

(10)

式中：A為年平均土壤侵蝕量，t·hm-2·a-1；R為降雨侵蝕力因子，MJ·mm·hm-2·h-1·a-1；K為土壤可蝕性因子，t·hm2·h·hm-2·MJ-1·mm-1;LS為地形因子，L和S分別為坡長和坡度因子；C為植被覆蓋與管理因子；P為水土保持措施因子。

2.3.1 降雨侵蝕力因子R

降雨侵蝕力是指由降雨導致土壤侵蝕的潛在能力，是侵蝕的主要驅動力[18]。本文R因子的計算公式為

R=4.1×F-152，

(11)

(12)

式中:F為修正參數；Pj為第j月的降雨量，mm；P為年平均降雨量，mm。將計算所得的各氣象站點的年降雨侵蝕力在ArcGIS中進行克呂格空間插值，得到研究區(qū)年降雨侵蝕力圖。

2.3.2 土壤可蝕性因子K

土壤可蝕性是反映土壤抵抗雨滴擊濺或徑流沖刷作用的能力。計算公式為

(13)

Sn=1-Sa/100，

(14)

式中：Sa為砂粒的含量，%；Si為粉粒的含量，%；Cl為粘粒的含量，%；C0為有機碳的含量，%。

2.3.3 植被覆蓋與管理因子C

植被覆蓋與管理因子反映了植被覆蓋或作物管理措施對土壤侵蝕量的影響[19]。計算公式為

(15)

(16)

(17)

式中：NDVI為歸一化植被指數；ρNIR和ρR分別為近紅外波段和紅波段反射率；f為植被覆蓋度，%；NDVImax和NDVImin分別為NDVI的最大和最小值。

2.3.4 水土保持措施因子P

運用研究區(qū)遙感影像在ENVI軟件平臺支持下建立解譯標志，采用人機交互解譯，將研究區(qū)分為耕地、林地、草地、建設用地、水體和濕地6類，對分類結果進行外業(yè)驗證與內業(yè)修正，直到滿足精度要求，得到研究區(qū)土地利用分類數據。結合土地利用數據和由DEM生成的坡度信息，參照方廣玲等[20]和黃金良等[21]的研究成果對P因子賦值，見表2。

表2 不同土地利用類型P值Tab.2 P values of different land uses

3 結果與分析

3.1 對比驗證

依照上述方法得到研究區(qū)各因子空間分布圖，分別與通用算法和本文改進算法提取的LS因子疊加運算，得到兩種算法的土壤侵蝕模數圖，并依據水利部頒布的《土壤侵蝕分類分級標準》進行分級，得到研究區(qū)土壤侵蝕強度空間分布圖(圖3)。結果顯示，通用算法研究區(qū)年平均侵蝕模數為12.2 t·hm-2·a-1，屬輕度侵蝕；改進算法的年平均侵蝕模數為28.16 t·hm-2·a-1，屬中度侵蝕，與信陽市水土保持監(jiān)測站所提供數據結果吻合[22]。將2種算法的結果進行疊加分析，可知2種算法的侵蝕等級相同的面積居多，但整體上改進算法得到的結果侵蝕更嚴重。從2種結果異同區(qū)域均勻選取122個特征點結合SPOT5影像進行人工目視解譯和野外核查驗證，可知本文方法有118處與解譯結果一致，驗證結果混淆表如表3所示，說明本文提出的算法計算結果更為真實有效。如圖4所示區(qū)域為采礦用地、坡地，無植被，斑塊為亮白色，與周圍地物特征差異明顯，判定為劇烈侵蝕。本文方法結果為劇烈侵蝕，通用方法結果卻為強烈侵蝕。其下坡方向的區(qū)域，因其帶來的雨水沖刷，侵蝕程度也應較強，本文方法結果為極強烈或劇烈侵蝕，通用方法結果卻為中度或強烈侵蝕，不及本文方法準確可靠。

(a)通用算法 (b)本文算法

圖3 土壤侵蝕強度空間分布

Fig.3Soilerosiondistributionmap

表3 本文方法試驗的混淆矩陣Tab.3 Confusion matrix of the method test in this paper

圖4 兩種結果對比驗證

Fig.4Comparisonandverificationoftworesults

3.2 土壤侵蝕坡度分布特征

統計分析本文算法的土壤侵蝕分布結果，可知研究區(qū)總侵蝕面積達905.95 km2，占區(qū)域總面積的42.94%，以輕度侵蝕為主。將侵蝕等級與坡度等級疊加統計分析(見表4)可知隨著坡度的增加，平均侵蝕模數顯著增加，且包含的主要的侵蝕強度等級越高；極強烈和劇烈主要分布在[15°，35°)坡度帶，該區(qū)域是受人類活動影響較大的易遭受侵蝕的坡耕地和林地的主要分布區(qū)域，所以侵蝕明顯，應作為治理的重點區(qū)域，對于坡耕地最好坡改梯或等高耕種，坡度≥25°的要堅決退耕還林；侵蝕面積并沒有隨著坡度的增加持續(xù)增加，在≥35°區(qū)域反而減少了，這與胡世雄等[23]、王星等[24]等的研究結果一致。

表4 不同坡度的土壤侵蝕分布Tab.4 Soil erosion distribution of different slope degrees

3.3 土壤侵蝕土地利用分布特征

將土壤侵蝕結果與土地利用分布圖疊加分析(見表5)，可知林地和草地年平均侵蝕模數較高；極強和劇烈侵蝕主要分布在林地。研究區(qū)森林覆蓋率占49.46%，使得植被覆蓋整體較高而侵蝕強度相對較低。但研究區(qū)林地多為實施退耕還林后的人工林、順坡種植且無任何水土保持措施的經濟林和疏林地，致使土壤侵蝕嚴重。草地包含封山育林形成的低覆蓋度荒草地和人工種植草地，面積比例不高但侵蝕明顯。

表5 不同土地利用的土壤侵蝕分布Tab.5 Soil erosion distribution of different land uses

4 結論

1)本文基于多流向匯流累積算法提取L因子，并考慮了地表土地利用/覆蓋的影響加以改進。綜合利用RS/GIS技術，基于RUSLE模型將通用算法和改進后算法的土壤侵蝕成果進行對比驗證，實驗證明了本文算法的有效性。

2)研究區(qū)年平均土壤侵蝕模數為28.16 t·hm-2·a-1，屬于中度侵蝕，總侵蝕面積達905.95 km2，以輕度侵蝕為主。受地形影響，侵蝕主要分布在中南部地勢較陡區(qū)域。

3)隨著坡度升高，平均侵蝕模數和包含的主要土壤侵蝕等級也升高，極強烈和劇烈侵蝕主要分布在[15°，35°)坡度區(qū)域，應重點治理。加強經濟林和疏林地改造，保護高密度草地和濕地，防止演化為水土易流失的其他土地利用/覆蓋類型是研究區(qū)長期的任務。