基于GIS和RUSLE的淮河流域土壤侵蝕研究
——以信陽(yáng)市商城縣為例

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

(河南理工大學(xué)測(cè)繪與國(guó)土信息工程學(xué)院，焦作 454003)

0 引言

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

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

1 地形因子提取方法

1.1 通用算法

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

(1)

L=(λ/22.1)m，

(2)

(3)

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

通用的RUSLE模型提取的L因子忽略了上坡匯流作用，導(dǎo)致L因子小于真實(shí)值，影響土壤侵蝕評(píng)定的準(zhǔn)確性。本文采用基于上坡匯流的L因子算法對(duì)L因子的提取進(jìn)行了改進(jìn)。

1.2 本文提出的改進(jìn)算法

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

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

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

(4)

式中：Ai為第i個(gè)相鄰下坡單元的累積匯流面積；A為總上坡面積；tanβi為第i個(gè)下坡方向的坡角，βi為坡度；Wi為第i個(gè)與流向垂直的有效等高線長(zhǎng)，為柵格尺寸與權(quán)重的積(主方向上權(quán)重為0.5，對(duì)角方向?yàn)?.354)；k為下坡單元個(gè)數(shù)。

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

(5)

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

(6)

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

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

表1 不同土地利用匯流面積貢獻(xiàn)率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為所在柵格的有效等值線長(zhǎng)度；θi,j為坡向。將(7)代入(6)中得到

(9)

算法實(shí)現(xiàn)首先要對(duì)DEM進(jìn)行填洼處理；其次要注意坡長(zhǎng)的截?cái)辔恢?，本文選擇水體作為坡長(zhǎng)L因子的計(jì)算終點(diǎn)；最后使用MATLAB編程完成算法。技術(shù)路線如圖2。

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

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

2.1 研究區(qū)概況

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

2.2 數(shù)據(jù)收集

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

2.3 RUSLE模型

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

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分別為坡長(zhǎng)和坡度因子；C為植被覆蓋與管理因子；P為水土保持措施因子。

2.3.1 降雨侵蝕力因子R

降雨侵蝕力是指由降雨導(dǎo)致土壤侵蝕的潛在能力，是侵蝕的主要驅(qū)動(dòng)力[18]。本文R因子的計(jì)算公式為

R=4.1×F-152，

(11)

(12)

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

2.3.2 土壤可蝕性因子K

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

(13)

Sn=1-Sa/100，

(14)

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

2.3.3 植被覆蓋與管理因子C

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

(15)

(16)

(17)

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

2.3.4 水土保持措施因子P

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

表2 不同土地利用類(lèi)型P值Tab.2 P values of different land uses

3 結(jié)果與分析

3.1 對(duì)比驗(yàn)證

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

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

圖3 土壤侵蝕強(qiáng)度空間分布

Fig.3Soilerosiondistributionmap

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

圖4 兩種結(jié)果對(duì)比驗(yàn)證

Fig.4Comparisonandverificationoftworesults

3.2 土壤侵蝕坡度分布特征

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

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

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

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

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

4 結(jié)論

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

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

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