劉歡歡, 剛成誠, 溫仲明,, 陳同德, 劉 悅, 陳 印

(1.西北農林科技大學 草業(yè)與草原學院, 陜西 楊凌 712100; 2.西北農林科技大學 水土保持研究所,陜西 楊凌 712100; 3.中國科學院 水利部 水土保持研究所, 陜西 楊凌 712100)

土壤侵蝕是一個全球性的環(huán)境問題，不僅降低土壤質量、破壞土地資源，而且淤積河道、加劇洪澇災害，對農業(yè)生產和生態(tài)環(huán)境構成了極大的威脅，嚴重制約了社會經濟的可持續(xù)發(fā)展[1-2]。據估算，全球土壤侵蝕的面積約為1.64×107km2，約占陸地總面積的10.95%[1]。水利部2018年全國水土流失動態(tài)監(jiān)測結果顯示，2018年全國發(fā)生土壤侵蝕的面積約為2.73×106km2，占國土總面積的28.6%[3]。我國每年因土壤侵蝕造成的經濟損失高達百億元以上[1,4]。因此，科學計算流域土壤侵蝕動態(tài)變化是深入認識土壤侵蝕發(fā)生過程與機理的重要步驟，對于評估區(qū)域水土保持措施有效性及合理調整具有重要的理論與現(xiàn)實意義。

土壤侵蝕模型的發(fā)展為計算不同時空尺度的土壤侵蝕提供了有力的技術手段[5-7]。目前在我國應用比較廣泛的土壤侵蝕模型包括通用土壤流失方程(Universal Soil Loss Equation, USLE)、修正的通用土壤流失方程(Revised Universal Soil Loss Equation, RUSLE)和中國土壤流失方程(Chinese Soil Loss Equation, CSLE)。如王濤等[8]利用RUSLE模型計算了2000—2014年陜北無定河流域土壤侵蝕的時空變化；Teng等則利用RUSLE模型預測了未來青藏高原的土壤侵蝕[9]；USLE被用于計算祁連山[10]、湄公河流域[11]及長株潭城市群等[12]地區(qū)的土壤侵蝕和水土流失特征。CSLE也被用于不同縣域、省(市)區(qū)及流域尺度的土壤侵蝕研究[13-15]。

黃土高原由于其特殊地理位置及氣候條件，加之長期以來不合理的人類活動，使其成為中國乃至全球土壤侵蝕最為強烈的地區(qū)之一，生態(tài)環(huán)境極其脆弱[16-17]。自上個世紀末，國家先后實施了“退耕還林還草”等一系列生態(tài)工程，黃土高原的生態(tài)環(huán)境得到了有效的改善，植被覆蓋度增加，多年平均入黃泥沙量從16億t銳減至2～3億t[18-19]，可見植被在抵抗土壤侵蝕中的重要作用。現(xiàn)有土壤侵蝕模型中一般采用給不同土地利用類型賦值或者基于植被覆蓋度來表征植被在水土保持中的作用，即植被覆蓋管理因子(C)[6-7,15,20]。但由于植被覆蓋度僅能體現(xiàn)植被的投影蓋度，即植被的水平結構信息，無法體現(xiàn)植被的層次結構，即垂直結構信息，因此其不能完全真實的反映植被的水土保持功能。據此，雷婉寧和溫仲明等在野外大量觀測數據的基礎上，提出了結構化植被指數(Cs)的概念[21]。與傳統(tǒng)的歸一化植被指數(NDVI)相比，Cs在高植被覆蓋區(qū)，能夠有效避免植被指數飽和現(xiàn)象，而且能夠將不同層次的植被信息相耦合，以表達植被群落更為完整的水土保持功能[21-22]。

為進一步了解植被結構對土壤侵蝕的影響，本研究分別采用Cs和NDVI計算植被覆蓋管理因子，在此基礎上，利用RUSLE模型模擬2000—2018年延河流域逐年土壤侵蝕量，分析其時空動態(tài)特征，比較NDVI和Cs在流域水土流失評價中的可靠性。研究結果一方面為土壤侵蝕模型的植被因子參數優(yōu)化提供科學依據，另一方面為黃土高原植被結構優(yōu)化調控和管理提供科學參考。

1 材料與方法

1.1 研究區(qū)概況

延河是黃河中游的一級支流，發(fā)源于靖邊縣境內，由西北向東南流經志丹、安塞和延安，在延長縣南河溝鄉(xiāng)涼水岸附近匯入黃河，全長286.9 km[7]。主要支流包括杏子河，西川、潘龍川和南川等，整個流域面積為7 725 km2，位于陜北黃土高原中部丘陵溝壑區(qū)(36°21′—37°19′N, 108°38′—110°29′E)，地勢西北高東南低(圖1)，地形破碎復雜，黃綿土占流域總面積的85%以上，抗侵蝕能力差[7,23]。

該流域屬于溫帶大陸性半干旱氣候區(qū)，年均溫約為9℃，年降水量500 mm，大部分地區(qū)干旱少雨，降雨季節(jié)分配不均，夏季多暴雨，降水變化率大。植被類型以次生植被和中、低覆蓋度草地為主，天然林集中部分在流域的西南部。常見的植物種從南向北有遼東櫟、刺槐、油松。延安至安塞之間為檸條和白羊草，安塞以北為百里香、長芒草等[23-24]。

圖1 延河流域地形

1.2 數據來源

本研究所用的數據包括：(1) DEM數據，來源于中國科學院資源環(huán)境科學與數據中心，空間分辨率為30 m；(2) 氣象數據：流域內及周邊氣象站數據來源于中國氣候數據共享網；(3) 土地利用數據：來源于中國科學院資源環(huán)境科學與數據中心，空間分辨率為30 m；(4) 遙感影像數據：2000—2018年延河流域6—9，11—2月的Landsat 7 ETM和Landsat 8 OLI遙感影像(影像編號127034，127035，126035)，來源于美國USGS網站，利用ENVI軟件進行缺失條帶插補、輻射定標和大氣校正等操作，提取遙感指數應用于RUSLE模型中。

1.3 土壤侵蝕模型

本研究采用修正的通用土壤流失方程RUSLE來計算延河流域的土壤侵蝕模數：

A=R·K·LS·C·P

(1)

式中:A為年均土壤侵蝕模數,是指單位面積單位時間和空間的平均土壤流失量〔t/(km2·a)〕;R為降雨侵蝕力因子〔MJ·mm/(km2·h·a)〕;K為土壤可蝕性因子〔t·hm2·h /(km2·MJ·mm)〕;L和S分別為坡長和坡度因子;C為植被覆蓋與管理因子;P為水土保持措施因子。

1.3.1 降雨侵蝕力因子R利用逐月雨量估算延河流域的降雨侵蝕力，其計算方法為：

(2)

R=αFβ

(3)

式中:pi為逐月降雨量(mm);p為年均降水量(mm);n為產生侵蝕性降雨(10 mm)的月份數[25];F為修正參數;R為年均降雨侵蝕力〔MJ·mm/(km2·h·a)〕;α和β為模型參數,依據現(xiàn)有研究[7,26],其值分別取0.183 3,1.995 7。將延河流域內及周邊氣象站的月降水數據利用薄板樣條插值方法插值為空間分辨率為30 m的空間圖,再利用上述公式計算出降雨侵蝕力的空間分布。

1.3.2 土壤可蝕性因子K利用EPIC模型計算延河流域不同土壤類型的K值[27-28]，其計算公式為：

K={0.2+0.3exp〔0.0256SAN(1-SIL/100)〕}×

(4)

式中:K為土壤可蝕性因子;SAN,SIL和CLA分別是砂粒、粉粒和黏粒的有機碳含量(%),其中SN1=1-SAN/100。

1.3.3 坡度和坡長因子LS在流域尺度，坡度和坡長可以用過DEM提取。本研究采用劉寶元等在黃土高原建立的坡度坡長因子計算方法[29-30]：

L=(λ/22.1)α

(5)

(6)

α=β/(1+β)

(7)

β=(sinθ/0.089)/〔3.0(sinθ)0.8+0.56〕

(8)

式中:L為坡長因子;λ為由DEM提取的坡長(m);22.1為標準小區(qū)坡長;α為坡度坡長指數;S為坡度因子；θ為由DEM提取的坡度值。

1.3.4 植被覆蓋管理因子C目前，已有大量關于植被指數和C因子的關系式，本研究采用下列公式進行計算：

(9)

式中:α和β的取值分別為2,1。由于水體和建筑用地的NDVI較低,接近于0,但他們的土壤侵蝕強度很小,因此將水體和建筑用地的C值賦值為0[7]。

NDVI的計算公式為：

(10)

式中:NIR和R分別為Landsat 7,8影像的近紅外和紅光波段。

在計算流域尺度Cs時，土壤調整植被指數(MSAVI)能夠有效的去除土壤背景值的影響，用來反映植被綠色覆蓋狀況。在提取枯枝落葉層信息時，歸一化耕作指數(NDTI)與Cs的相關性要高于其他植被黃度指數[21-22]。因此，利用夏季和冬季的遙感影像提取MSAVI和NDTI來綜合表征Cs，其公式為：

(11)

(12)

Cs=0.175MSAVI+0.522NDTI

(13)

式中:NDTI為歸一化耕作指數;SWIR1和SWIR2分別為Landsat 7,8影像的短波紅外波段1,2;MSAVI為修正的土壤調整植被指數。關于Cs的驗證過程詳見參考文獻[21-22]。

1.3.5 水土保持措施因子P依據現(xiàn)有研究[7,20]，采用給不同土地利用類型賦值的方法確定P值。水體、建筑用地、林地和草地一般未采取水土保持措施，因此賦值為1。對于耕地而言，通常坡度越大，水土保持措施的作用越重要，因此，依據表1按照坡度范圍對耕地進行賦值。

2 結果與分析

2.1 2000-2018年延河流域土壤侵蝕強度的時間動態(tài)

2000—2018年，基于Cs的延河流域平均土壤侵蝕模數為8 354.62 t/(km2·a)，高于基于NDVI得到的6 421.13 t/(km2·a)(圖2)?；贑s和NDVI的土壤侵蝕模數的變化趨勢相同，在2000—2006年下降明顯，之后趨于平穩(wěn)。最高值均出現(xiàn)在2001年，分別為13 916.34，15 867.67 t/(km2·a)，之后逐漸下降，在2006年分別下降為4 092.22，5 962.21 t/(km2·a)。2014年土壤侵蝕模數達到最低值，之后又呈現(xiàn)上升的趨勢。基于Cs和NDVI的2018年土壤侵蝕模數分別為7 102.62，5 053.78 t/(km2·a)，與2000年相比，分別下降了44.15%和56.41%。

表1 不同坡度范圍耕地的P值

依據水利部頒發(fā)的《土壤侵蝕分類標準》，將土壤侵蝕強度分為微度侵蝕〔0～500 t/(km2·a)〕、輕度侵蝕〔500～2 500 t/(km2·a)〕、中度侵蝕〔2 500～5 000 t/(km2·a)〕、強烈侵蝕〔5 000～8 000 t/(km2·a)〕、極強烈侵蝕〔800～15 000 t/(km2·a)〕和劇烈侵蝕〔>15 000 t/(km2·a)〕。統(tǒng)計得到2000—2018年延河流域不同土壤侵蝕等級面積百分比(圖3)。2000—2018年，基于NDVI的土壤侵蝕強度為強烈、極強烈和劇烈侵蝕的面積分別下降了29.66%，44.48%和62.75%，而基于Cs的土壤侵蝕強度為強烈、極強烈和劇烈侵蝕的面積分別下降了12.17%，28.43%和51.21%。輕度和微度侵蝕的面積都呈現(xiàn)上升的趨勢。基于Cs和NDVI的結果表明，在2000—2018年，延河流域的土壤侵蝕強度明顯下降。在2018年，基于Cs的強烈、極強烈和劇烈侵蝕的面積占流域總面積的36.84%，而基于NDVI的三者面積占比為27.22%。

圖2 2000-2018年延河流域土壤侵蝕模數

圖3 基于NDVI(A)和Cs(B)計算的2000-2018年延河流域不同侵蝕強度比例變化

2.2 2000-2018年延河流域土壤侵蝕的空間動態(tài)

2000年、2009年、2018年延河流域土壤侵蝕的空間分布見圖4。在2000年，土壤侵蝕模數大于5 000 t/(km2·a)的區(qū)域廣泛分布于延河流域的上中下游地區(qū)，輕度和微度侵蝕的區(qū)域主要分布在延河沿岸地帶及流域的西南部；在2009年，輕度和微度侵蝕的面積增加，主要分布在延河流域的中游地區(qū)，而且基于NDVI的增加面積大于基于Cs的面積；在2018年，延河流域中游和下游地區(qū)微度和輕度侵蝕的面積進一步增加，土壤侵蝕模數大于5 000 t/(km2·a)的區(qū)域主要集中在延河流域的上游地區(qū)。

根據2.1中的研究結果，為進一步揭示土壤侵蝕強度變化的空間分布特征，分析了延河流域2000—2007年和2007—2018年土壤侵蝕強度空間轉化(圖5)?；贜DVI的2000—2007年，土壤侵蝕強度等級不變的面積占流域總面積的36.13%，主要分布在延河的中游和沿岸地區(qū)；土壤侵蝕等級降低的區(qū)域則占62.42%，其中下降一級的面積占34.17%，其廣泛分布于整個流域；在2007—2018年，土壤侵蝕強度等級不變的面積增加到45.20%，侵蝕等級下降的面積降低為35.38%，主要分布在流域的中游和下游，侵蝕強度增加的區(qū)域主要分布在流域的上游。在2000—2018年，土壤侵蝕強度下降的面積占流域總面積的62.45%。

在2000—2007年，基于Cs的延河流域內土壤侵蝕強度無變化的區(qū)域占流域總面積的54.84%，侵蝕強度下降的面積占42.14%；在2007—2018年，流域內土壤侵蝕強度無變化和下降的面積分別下降到46.44%和36.71%。從整個研究時間范圍來看，土壤侵蝕強度以下降趨勢為主，面積占53.70%，其主要分布在延河流域的中游和下游地區(qū)。在2000—2018年，土壤侵蝕強度下降地區(qū)的面積占流域總面積的53.70%，低于基于NDVI的結果?？傮w而言，基于Cs和NDVI的結果均表明延河流域土壤侵蝕狀況得到明顯改善。

圖4 基于NDVI(A, B, C)和CS(D, E, F)的延河流域2000年、2009年、2018年土壤侵蝕模數變化

圖5 基于NDVI(A, B, C)和Cs(D,E,F)的延河流域2000-2007年、2007-2018年、2000-2018年土壤侵蝕強度的變化

3 討 論

植被和降雨是影響土壤侵蝕動態(tài)的兩個主要因素。本研究發(fā)現(xiàn)，在2000—2018年，延河流域的NDVI和Cs均呈現(xiàn)上升的趨勢，增加速率分別為0.012/a和0.014/a，尤其是在2000—2006年，其增加速率均高達0.05/a(圖6)。自20世紀90年代初開始，延河流域開展了大規(guī)模的退耕還林還草、封山育林等植被恢復工程，造林種草的實施力度增強，截至2006年，延河流域林草植被種植面積高達2 370.05 km2[31-32]。因此，延河流域土壤侵蝕模數下降的主要原因是黃土高原大范圍開展的植被恢復措施，這與其他研究的結果一致[7,33]。但由于某些年份的降雨量較大，導致土壤侵蝕模數仍較高，如2011年、2013年、2017年，降雨侵蝕力分別達到了2 371.14，2 950.10，2 378.16 MJ·mm/(km2·h·a)，其對應的土壤侵蝕模數高于6 000 t/(km2·a)。盡管如此，土壤侵蝕強度仍低于2001年、2003年，進一步說明了生態(tài)工程措施的有效性。

不同土壤侵蝕模型及因子的不同計算方法均會導致結果的差異，如李天宏等利用RUSLE模型計算延河流域2006年的土壤侵蝕模數約為3 000 t/(km2·a)[7]，而謝紅霞等利用CSLE模型計算的同年結果為5 009 t/(km2·a)[23]。本研究基于NDVI和Cs得到的2006年土壤侵蝕模數分別為4 092.22，5 962.21 t/(km2·a)。本研究中，水土保持措施因子P采用的是給不同土地利用類型及不同坡度耕地賦值的方法，而未采用Lufafa等[34]利用坡度計算P的方法，該方法在黃土高原亦被廣泛采用[8,24,35]。該方法僅按照坡度計算p值，不區(qū)分土地利用類型，即同一坡度的農田和自然植被有相同的水土保持措施因子。若本文采用此方法計算P，進而得到的土壤侵蝕模數將低于現(xiàn)研究結果。此外，本研究未將淤堤壩和梯田等工程措施計算在內，也會對結果造成一定的誤差。

圖6 2000-2018年延河流域NDVI，Cs及年降雨侵蝕力變化

植被覆蓋管理因子的準確估算對于土壤侵蝕預測和植被水土保持效益評價具有重要的意義。喬木林中若無灌木或草本植被，林下植被結構的不完整可導致植被群落的水土保持功能下降，其水土流失可能仍然較為嚴重。灌草層和枯枝落葉層能夠有效的攔蓄徑流，減少雨滴動能，增加土壤入滲，從而降低土壤侵蝕，在控制土壤侵蝕過程中發(fā)揮著重要作用[36]。傳統(tǒng)的NDVI指數由于只考慮了植被蓋度，無法反映林下植被結構狀況，因此，會存在低估土壤侵蝕強度的現(xiàn)象。而結構化植被指數Cs不僅考慮了植物群落的水平結構特征，還考慮了包括植被喬木層、灌木層、草本層及枯枝落葉層不同垂直結構特征在水土保持作用中的差異，并且較好了克服了NDVI在高植被覆蓋區(qū)易飽和的現(xiàn)象[21]?；谶b感影像，以MSAVI為綠度指數和NDTI為黃度指數提取的Cs能夠在流域尺度更好的反映植被群落的水土保持作用與效益。本文研究結果表明2018年延河流域上游地區(qū)的土壤侵蝕模數仍較高，且Cs的土壤侵蝕強度要大于基于NDVI的結果，說明這些地區(qū)未來水土保持工作仍需進一步加強或者調整。

4 結 論

(1) 在2000—2018年，基于Cs和NDVI得到的延河流域土壤侵蝕模數分別下降了44.15%和56.41%；基于Cs計算的強烈侵蝕、極強烈和劇烈侵蝕的面積分別下降了12.17%，28.43%和51.21%，而基于NDVI計算的比例分別下降了29.66%，44.48%和62.75%；

(2) 在2000—2007年，土壤侵蝕強度下降的面積占流域總面積的42.14%，廣泛分布于整個延河流域；在2007—2018年，侵蝕強度持續(xù)下降的面積占流域總面積的36.71%，主要分布在延河的中游和下游；在整個研究時間范圍內，延河流域土壤侵蝕強度下降的面積占流域總面積的53.70%；

(3) Cs的變化趨勢表明生態(tài)工程的實施促進了植被的恢復，有效降低了土壤侵蝕強度。延河流域的上游土壤侵蝕強度仍較高，未來應進一步加強實施水土保持措施。

(4) 基于Cs的植被覆蓋管理因子不僅能夠反映生態(tài)工程實施以來植被的改善，而且由于其包含植被垂直結構特征，能夠更好的描述植被在控制土壤侵蝕中的作用，為區(qū)域尺度水土保持措施中植被結構優(yōu)化和調控提供了參考。未來可考慮利用無人機遙感或更高時空分辨率的遙感影像在不同的區(qū)域尺度開展進一步研究。