楊 波,王全九,郭思琪

(1.咸陽師范學(xué)院資源環(huán)境與歷史文化學(xué)院,712000,陜西咸陽;2.西安理工大學(xué)水利水電學(xué)院,710048,西安;3.長安大學(xué)地球科學(xué)與資源學(xué)院,710054,西安)

黃土高原是中國土壤侵蝕最為嚴(yán)重的地區(qū),總面積64萬km2,土壤侵蝕面積達39萬km2,土壤侵蝕模數(shù)≥1萬5 000 t/(km2·a)的劇烈水蝕面積3.67萬km2,占全國同類面積的89%[1]。陜北風(fēng)沙區(qū)位于毛烏素沙漠和黃土高原的過渡地帶,植被破壞比較嚴(yán)重,水土流失加劇,為了治理水土流失,改善自然環(huán)境,國家從1999年開始試行退耕還林生態(tài)工程,陜北地區(qū)植被指數(shù)(NDVI)有了明顯提高[2]。禿尾河是黃河中游陜北風(fēng)沙區(qū)的典型流域之一,因此該流域的土壤侵蝕時空變化能代表風(fēng)沙區(qū)的總體趨勢,有著重要的意義。土壤侵蝕的時空變化是水土保持綜合治理成果體現(xiàn),傳統(tǒng)的統(tǒng)計方法難以及時掌握水土治理工作的成效,有一定的滯后性。為了及時掌握土壤侵蝕變化的時空分布,筆者解譯3景不同時期禿尾河流域Landsat衛(wèi)星TM影像,與USLE土壤侵蝕模型結(jié)合,估算該流域25年以來的土壤侵蝕變化情況,為陜北風(fēng)沙區(qū)水土流失治理提供參考。

1 研究區(qū)域概況

禿尾河是黃河一級支流,發(fā)源地于陜西省神木縣瑤鎮(zhèn)鄉(xiāng)宮泊海子,禿尾河全長140.0 km,流域面積3 294.0 km,河道平均比降3.87‰,流域面積為3 295 km2,自西北向東南分別流經(jīng)神木、榆林、佳縣等3個縣16個鄉(xiāng)(鎮(zhèn))。流域內(nèi)地貌類型可分為草灘區(qū)、流動風(fēng)沙區(qū)、蓋沙區(qū)和黃土丘陵區(qū)4個區(qū)。大陸性季風(fēng)半干旱氣候區(qū),多年平均降水量394.5 mm,降水總體上由西北向東南遞增,年蒸發(fā)能力1 500 mm[3]。截至2011年9月,流域水土保持措施面積為2 137.56 km2,占流域總面積的64.9%。上游為風(fēng)沙區(qū),地表植被稀疏,主要植物有沙竹(Phyllostachys propinqua)、沙芥(Pugionium cornutum)、沙蓬(Agriophyllum squarrosum)、沙蒿(Artemisia desertorum)、沙柳(Salix cheilophila)、沙棘(Hippophae rhamnoides)等;中下游為黃土丘陵溝壑區(qū),面積為1 121.95 km2,植被側(cè)柏(Platycladus orientalis)、刺槐(Robinia pseudoacacia)、人工喬木和檸條(Caragana korshinskii)、沙棘、草地主要有紫花苜蓿(Medicago sativa)等,天然草以蒿類為主[4]。

圖1 禿尾河流域Fig.1 Map of Tuwei river watershed

2 材料與方法

2.1 數(shù)據(jù)來源

筆者研究使用的數(shù)據(jù)有：1)1988年8月29號,2000年8月24號的landsat 5的TM影像,2013年8月6號landsat 8的TM影像,空間分辨率為30 m,來源于中科院遙感與地球研究所下載地址(http：∥ids.ceode.ac.cn/Index.aspx);2)土壤數(shù)據(jù)來源于陜西省第2次土壤普查數(shù)據(jù)集、陜西土壤和南京土壤所的土壤類型圖(shp格式),土壤屬性數(shù)據(jù)和空間數(shù)據(jù)在ArcGIS中關(guān)聯(lián)后,得到土壤屬性空間數(shù)據(jù)集;3)1∶1萬DEM數(shù)據(jù)用于計算SL因子;4)高家川站、高家堡站、榆林、神木等6個氣象站點日降雨資料來源于黃土高原科學(xué)數(shù)據(jù)中心和中國氣象共享網(wǎng)。2000年和2013年影像在ENVI解譯后,土地利用類型和Googlearth典型地物對比精度符合要求。

2.2 研究方法

隨著GIS技術(shù)的發(fā)展,USLE通用土壤侵蝕模型逐漸和GIS、RS融合,學(xué)者們在黃土高原地區(qū)利用該模型在小流域[5-6]和大流域[7-9]對土壤侵蝕進行了估算。由于黃土高原地區(qū)細(xì)溝侵蝕比較嚴(yán)重,該模型估算出來的侵蝕模數(shù)比真實值偏小,在延河流域的估算值在82%～92%之間[7];但該方法仍是區(qū)域尺度下土壤定量侵蝕估算的最佳選擇之一。USLE土壤侵蝕模型如下：

式中：A為土壤流失量,t/(km2·a);R為降雨侵蝕力因子,MJ/(mm·hm2·h);L為坡長因子;S為坡度因子;K 為土壤可蝕性因子,t·hm2·h/(hm2·MJ·mm);C為植被覆蓋和管理因子;P為水土保持措施因子。

2.2.1 降雨侵蝕力R因子的確定 日降雨量超過12 mm會產(chǎn)生土壤侵蝕,被稱為侵蝕性降雨[9]。運用公式計算降雨侵蝕力

式中：Mi為第i個半月的降雨侵蝕力[MJ/(mm·hm2·h)];Dj是半個月內(nèi)第j天降雨量,mm。這里的降雨量為侵蝕性降雨的降雨量,謝云等[10]在對降雨標(biāo)準(zhǔn)的研究中,要求使Dj≥12 mm,若＜12 mm計為0。k是研究期內(nèi)半月的時段數(shù),將每個月前15 d視為一個半月時段,剩下為另一個半月時段,全年共計24個時段;α和β為模型參數(shù)[11],其計算公式為：

式中：Pd12為日降雨量≥12 mm的日平均雨量,mm;Py12為年平均降雨量,mm。

2.2.2 坡長L因子與坡度S因子的確定 LS因子是地形對土壤侵蝕的影響,在區(qū)域尺度下,一般SL因子利用DEM計算。汪幫穩(wěn)等[12]利用修正算法在黃土高原地區(qū)通過DEM提取LS因子值,并對其準(zhǔn)確性進行了驗證,因此,本文將采用該算法。

坡度S因子修正后的算法為：

式中：S為坡度因子;[slope]為坡度的弧度形式。

坡長L因子的算法為：

式中：L為坡長因子;m為坡長指數(shù);λ為水平坡長,m;β為細(xì)溝和細(xì)溝間侵蝕的比率。

2.2.3 土壤可蝕性K因子的確定 土壤可蝕性是土壤自身的抗沖蝕能力[13]。經(jīng)過研究發(fā)現(xiàn),幾何平

式中：Dg為平均幾何粒徑,mm;fi為第i粒徑等級的比例;mi為該粒徑等級中粒徑最大值和粒徑最小值的算數(shù)平均值,mm;OM為有機質(zhì)比例。

2.2.4 植被覆蓋與管理C因子的確定 植被覆蓋與管理C因子被定義為在土地植被、作物覆蓋的特定條件下土壤的侵蝕量與土地?zé)o作物、連續(xù)休憩條件下土壤的侵蝕量的比值,反映植被、作物等覆蓋與管理措施對土壤侵蝕的影響。C因子量綱為1,取值范圍為[0,1]。歸一化植被指數(shù)(NDVI)是能夠很好的反映植被、作物覆蓋度的指數(shù)。本文采用的算法是參考蔡崇法等[14]研究提出的計算C因子的算法,即與植被、作物覆蓋度有關(guān)：

式中：NDVI為歸一化植被指數(shù);NDVImin和NDVImax分別為研究區(qū)域NDVI的最小值和最大值。

2.2.5 水土保持措施P因子的確定 水土保持措施P因子,是指采用特殊侵蝕控制措施后的土壤侵蝕量與采用順坡種植時土壤侵蝕量的比值,取值在[0,1]之間[15],文中P因子的取值主要參照傅伯杰等[8]對黃土高原地區(qū)的研究成果,并在對黃土高原地區(qū)調(diào)查的基礎(chǔ)上,結(jié)合坡度確定P值因子。坡耕地 P 因子的賦值為：0°～5°是 0.100、5°～10°是0.221、10°～15°是 0.305、15°～20°是 0.575、20°～25°是0.705,在 ＞25°是 0.800。

3 結(jié)果與分析

3.1 主要因子結(jié)果

3.1.1 降雨侵蝕力動態(tài)變化分析 用禿尾河流域氣象站日降雨數(shù)據(jù),計算1988—2013年降雨侵蝕力,如圖2所示：25年來,降雨侵蝕力最大和最小出現(xiàn)在2012和1991年,分別是2 849.53和254.46 MJ·mm/(hm2·h)。2013年和1988年降雨侵蝕力比較接近,分別為 2 050.31 和 2 047.00 MJ·mm/(hm2·h)。1988—2000年,12年來多年平均降雨侵蝕力為779.55 MJ·mm/(hm2·h),2001—2013 年,13 年來多年平均降雨侵蝕力為1 479.70 MJ·mm/(hm2·h),1988—2000年和2001—2013年相比,降雨侵蝕力明顯增加。均粒徑結(jié)合有機質(zhì)模型最接近陜西地區(qū)真實K值。計算公式為：

圖2 1988—2013年降雨侵蝕力R因子統(tǒng)計圖Fig.2 Statistics of R factor from 1988 to 2013

3.1.2 植被覆蓋因子分析 1999—2010年神木縣累計實施退耕還林還草造林78.52萬畝[16](1 hm2=15畝)。1988、2000和 2013年平均NDVI值分別為0.09、0.11和0.15。退耕還林的水土保持效益在USLE模型中直接表現(xiàn)在C因子上。利用蔡崇法模型計算的C因子1988、2000和2013年的值分別為0.4、0.34和0.11。2013比1988年C因子的值減少0.29。

3.2 土壤侵蝕時間動態(tài)變化分析

按照水利部土壤侵蝕分類標(biāo)準(zhǔn)將土壤侵蝕強度劃分為微度(0＜A≤500)、輕度(500＜A≤2 500)、中度(2 500＜A≤5 000)、強烈(5 000＜A≤8 000)、極強烈(8 000＜A≤15 000)和劇烈(A＞15 000)6級[17],單位為 t/(km2·a)。 1988、2000 和 2013 年不同等級侵蝕強度占流域面積比例見(表1)。2013年8—9月野外采集樣本35個,微度6個、輕度8個、中度9個、強烈及以上12個,2013年計算結(jié)果進行對比檢驗,估算正確個數(shù)為：微度5個,輕度7個,中度7個,強烈及以上8個,正確率為77.14%。

表1 1988、2000、2013年土壤侵蝕統(tǒng)計Tab.1 Statistics of soil erosion categories in 1988,2000 and 2013

1988和2013年降雨侵蝕力分別2 050.31和2 047.00 MJ·mm/(hm2·h),降雨侵蝕力相似,土壤侵蝕模數(shù)分別為12 434.47 t/(km2·a)和3 721.08 t/(km2·a),相差 8 713.39 t/(km2·a)。 2000 年的降雨侵蝕力是25年最低的,僅有332.77 MJ·mm/(hm2·h),土壤侵蝕模數(shù)為 1 776.70 t/(km2·a)。 從侵蝕的空間分布(圖3)來看,1988年中度及以上侵蝕占70.47%約2 479.84 km2,在整個流域內(nèi)呈現(xiàn)面狀和帶狀連續(xù)分布,輕度和微度侵蝕呈零星分布,強烈及以上侵蝕在整個流域的上中下游普遍分布。2000年,全流域中度以及上侵蝕占21.13%,約694.36 km2。2013年,中度及以上侵蝕占到42.14%約1 483.53 km2,中度及以下侵蝕占總面積的76.71%,約2 700.34 km2。在流域的分布特征為,中上游地區(qū)以中度及以下侵蝕為主,成面狀分布,強烈及以上侵蝕為主要分在中下游地區(qū),以零星或線狀分布。其余各等級侵蝕呈現(xiàn)零星交錯分布?？傮w相比1988年等級減弱,中度及以上侵蝕減少28.33%,約996.31 km2?？偳治g量1988年和2013年分別為4 376萬t和1 310萬t,退耕還林14年以來該流域土壤侵蝕無論從面積還是強度上都有大幅下降。將1988年和2013年的土壤侵蝕強度等級圖做差運算后。土壤侵蝕模數(shù)減少8 713.39 t/(km2·a),總侵蝕量減少3 066萬t,2013年相比1988年總侵蝕量減少68.33%。整個流域土壤侵蝕模數(shù)大幅減少。尤其是中度以上侵蝕由1988年的面狀分布,變成2013年的線狀和星狀分布。

圖3 1988和2013年土壤侵蝕等級分布Fig.3 Map of soil erosion categories in 1988 and 2013

為了進一步分析禿尾河流域上中下游的土壤侵蝕變化情況,在ArcGIS水文分析模塊中,將禿尾河進一步細(xì)分為17個子流域(圖4),1988年編號為1,5,11的子流域是強烈侵蝕,編號為6的子流域是劇烈侵蝕,其余13個子流域均為極強烈侵蝕,2000年由于降雨偏少,流域內(nèi)均為輕度侵蝕,2013年流域內(nèi)編號為1,2,3,4,5,7的子流域均為輕度侵蝕,編號為6,7,8,9,10,11,12,13,14,15的子流域均為中度侵蝕,編號為16,17的子流域為強烈侵蝕。禿尾河流域土壤侵蝕強度從上游向下游逐漸增加的特征。

利用1988年和2013年土壤侵蝕強度等級圖,繪制土壤侵蝕等級轉(zhuǎn)移矩陣(表2),25年來微度侵蝕25.94%的面積增加為輕度侵蝕;輕度侵蝕35.79%轉(zhuǎn)變成微度侵蝕;中度侵蝕中68.09% 轉(zhuǎn)為輕度及以下侵蝕;強烈侵蝕中82.33% 轉(zhuǎn)為中度及以下侵蝕;極強烈侵蝕中88.87% 轉(zhuǎn)為強烈及以下侵蝕;劇烈侵蝕中89.06% 轉(zhuǎn)變?yōu)闃O強烈及以下侵蝕。1988—2013年禿尾河流域微度、輕度、中度、強度侵、極強和劇烈侵蝕的穩(wěn)定率分別為57.07%、37.65%、16.10%、8.33%、7.60%和10.94%。中度及以上侵蝕的穩(wěn)定率較低,是因為該等級土壤侵蝕強度大部分轉(zhuǎn)換為更低一級侵蝕等級,從土壤侵蝕轉(zhuǎn)移矩陣來看,強度侵、極強和劇烈侵蝕的面積在不斷下降。

圖4 1988、2000和2013子流域土壤侵蝕等級圖Fig.4 Map of sub-catchment soil erosion categories in 1988,2000 and 2013

表2 1988、2013年土壤侵蝕強度轉(zhuǎn)移矩陣Tab.2 Transformation matrix of soil erosion intensity from 1988 to 2013 %

3.3 不同坡度、高程等級下土壤侵蝕變化特征

木獨兔、邊老楞、瑤鎮(zhèn)、大保當(dāng)處于流域上游地區(qū),地貌以沙漠為主平均坡度2.65°,高佳堡鎮(zhèn),大河塔,安崖鎮(zhèn),喬岔灘、上高寨鄉(xiāng),劉國具鄉(xiāng)處于中下游地區(qū),地貌以黃土溝壑為主,平均坡度7.82°。將1988和2013年土壤侵蝕強度等級和坡度等級(0°～ 5°,5°～ 8°,8°～ 15°,15°～ 20°,20°～ 25°,＞25°)做疊置分析(表3)。2013比1988年坡度0°～5°強烈及以上減少 38.81%,5°～8°和 8°～15°范圍內(nèi),極強烈侵蝕減少26.87%和23.36%;15°～20°、20°～25°, ＞ 25°劇烈侵蝕分別減少 32.34%、32.52%和31.46%。在相同坡度條件下,中度及以上侵蝕在逐漸向中度及以下侵蝕轉(zhuǎn)變趨勢。

表3 不同坡度下土壤侵蝕強度比例疊置統(tǒng)計Tab.3 Overlay statistics of soil erosion intensity under different slope %

1988、2000、2013年侵蝕模數(shù)和坡度等級進行疊置分析(圖5),得到不同坡度下的侵蝕模數(shù)變化。侵蝕模數(shù)隨著坡度的增加而迅速增長。2013年相比1988年,不同坡度下的土壤侵蝕模數(shù)迅速減少,0°～ 5°減少 6 712.31 t/(km2·a);5°～ 8°減少7668.51 t/(km2·a),8°～15°減少 8 885.27 t/(km2·a);在20°～25°減少2 萬3 357.46 t/(km2·a);15°～20°減少 1萬 5 375.43 t/(km2·a);20°～25°減少2萬1 262.71 t/(km2·a);＞25°減少 2 萬 3 357.46 t/(km2·a)。 ＞25°的坡耕地退耕還林以后,先種植了沙棘,隨后幾年種適應(yīng)性強的油松和樟子松,隨著樹苗的生長,根系加固了松散的地表土壤;地表覆蓋的枯枝落葉層避免了降雨對地面的直接沖刷,有效的減少的坡地的土壤侵蝕量,同時植被也起到了涵養(yǎng)水源的作用,減少地表徑流。

圖5 不同坡度土壤侵蝕模數(shù)分布Fig.5 Distribution of soil erosion modulus under different slope

為分析禿尾河流域不同高程下的土壤侵蝕特征,DEM按照100 m的間隔分成7個等級。1988年,2000年和2013年土壤侵蝕模數(shù)與高程等級進行疊置分析(圖6)。高程在677～900 m土壤侵蝕最為嚴(yán)重,對應(yīng)的地貌是黃土高原溝壑地區(qū),1988年土壤侵蝕模數(shù)達到2萬2 146.6 t/(km2·a),2013年降低為9 201.26 t/(km2·a)。降低58%。900～1 100 m區(qū)間土壤侵蝕模數(shù),1988和2013年土壤侵蝕模數(shù)為1萬2 114.33 t/(km2·a)和4 589.55 t/(km2·a)。減少7 524.78,減少62.12%。退耕還林在這一地區(qū)取得了明顯的水土保持效益。1 100～1 385 m高程區(qū)間,地貌逐漸向風(fēng)沙區(qū)過渡,地面起伏不大,2013年土壤侵蝕模數(shù)降低到5 000 t/(km2·a)以下。其余高程區(qū)間內(nèi),劇烈、極強烈、強烈侵蝕面積均減少,侵蝕等級向低一級變化。

3.4 不同土地利用下土壤侵蝕變化特征

禿尾河流域1988—2013年土地利用類型發(fā)生了比較大的變化(表4)。耕地由總面積的27.72%減少到19.52%,減少面積為260.29 km2,林地由總面積的4.02%增加到7.58%,增加127.17 km2;草地由總面積的33.96%增加到50.42%,約552.56 km2;未利用土地由34.31%減少到14.22%,約649.05 km2。

表4 1988—2013年不同土地利用類型下的土壤侵蝕變化Tab.4 Statistics of soil loss change under different land use type from 1988 to 2013

圖6 不同高程土壤侵蝕模數(shù)分布Fig.6 Distribution of soil erosion modulus under different elevation

耕地土壤侵蝕模數(shù)由1988年的4 019.02 t/(km2·a)減少為 2013 年的 1 219.35 t/(km2·a),總流失量減少78.41%,284萬t,占總流失量的比例有27.72%減少到19.52%,減少8.2%。林地土壤侵蝕模數(shù)由1988年的6 948.62 t/(km2·a)減少為2013年的2 784.82 t/(km2·a),總流失量減少21%,19.1萬 t,占總流失量的比例由4.02%增加到7.85%,增加3.83%。草地土壤侵蝕模數(shù)由1988年的 12 675.15 t/(km2·a)減少為 2013年的4 677.85 t/(km2·a),總流失量減少 44.65%,626萬t,占總流失量的比例由33.96%增加到50.24%,增加16.46%。未利用土地中95%為沙地,土壤侵蝕模數(shù)由1988年的1.58萬t/(km2·a)到2013年變?yōu)榈? 021.07 t/(km2·a),總流失量減少89.36%,1 580萬t,占總流失量的比例由34.31%減少到14.22%,減少20.09%,在該流域的防風(fēng)固沙措施取得明顯的效果。在水土流失的結(jié)構(gòu)變化中,耕地所占比例迅速下降,而林地和草地上升,造成這一變化的是土壤侵蝕在大量減少過程中,相對比重變化造成的。禿尾河流域沙漠近25年來,荒漠化土地表現(xiàn)出向西北方向遷移趨勢[18],荒漠化程度減輕。

4 結(jié)論與討論

禿尾河流域25年來土壤侵蝕從強度和面積上都有明顯改善,流域中下游侵蝕強度高于上游地區(qū)。土壤侵蝕模數(shù)隨高程增加呈現(xiàn)遞減趨勢,流域上游的木獨兔鎮(zhèn)、大保當(dāng)鎮(zhèn)侵蝕強度明顯小于下游大河塔鎮(zhèn)、安崖鎮(zhèn)等地區(qū)。高程在677～900 m土壤侵蝕最為嚴(yán)重,1988年土壤侵蝕模數(shù)達到2.21萬t/(km2·a),2013 年變?yōu)闉?9 201.26 t/(km2·a),降低了58%。各個高程區(qū)間內(nèi)劇烈、極強烈、強烈侵蝕面積均減少,侵蝕等級向低一級變化。高程和土壤侵蝕模數(shù)關(guān)系體現(xiàn)了地貌對土壤侵蝕的影響。在降雨侵蝕力相似的1988年和2013年,土壤侵蝕模數(shù)分別為 1.24 萬 t/(km2·a)和 3 721.08 t/(km2·a),總侵蝕量分別為0.44億 t和0.13億 t,減少了68.33%。25年以來該流域NDVI均值從0.09增加到0.15,不同土地利用類型下土壤侵蝕強度大小依次是未利用土地、草地、林地和耕地。在水土流失治理后,未利用土地和草地的水土保持效益最明顯。退耕還林后14年風(fēng)沙區(qū)典型流域的水土流失治理取得了明顯效益。

在退耕還林前,修建淤地壩是主要的水土流失治理方式,20世紀(jì)90年代流域內(nèi)開荒川的淤地壩攔截效益為達到了52.9%[19]。目前禿尾河流域現(xiàn)有淤地壩430座(其中風(fēng)蝕區(qū)38座、水蝕區(qū)392座),治理水土流失面積2 137.56 km2。占流域面積的64.9%[4]。由于未能收集到淤地壩數(shù)據(jù),故未做其攔截效益分析,但魏艷紅等[20]評估禿尾河臨近皇甫川流域淤地壩泥沙攔截效益,認(rèn)為2000年后皇甫川流域輸沙量減幅超過85%,淤地壩攔沙已不是輸沙量減少的主要原因,而植被恢復(fù)等因素的作用已明顯得以發(fā)揮。土壤侵蝕的減少是氣候和植被變化共同影響的結(jié)果,1966—2010年禿尾河流域年降雨在400 mm左右,不同時段降雨量在流域內(nèi)部自北向南表現(xiàn)出遞增特征[21]。筆者計算流域內(nèi)本世紀(jì)前13年和上世紀(jì)后12年相比,降雨侵蝕力有所增加。氣候變化對土壤侵蝕的貢獻率在不同地區(qū)存在明顯差異[22],退耕還林對榆林市整體水土保持貢獻率約為40%～50%左右[22]。但是黃土高原地區(qū)近年來局部出現(xiàn)土壤干化趨勢[23],徑流明顯減少[24],在未干旱條件下可能威脅到草木生長。在后退耕還林時代的水土保持工作,要統(tǒng)籌考慮該地區(qū)的水資源承載力,在今后的生態(tài)治理中更加合理規(guī)劃樹草的種植品種和密度,是今后研究的重要方向。