退耕還林前后北洛河上游流域土壤侵蝕時空變化分析

陳 浩 ，張曉萍，權 偉，郭晉偉，袁 閱張思瑞

(1.延安大學 生命科學學院，陜西 延安 716000；2.陜西省區域生物資源保育與利用工程技術研究中心，陜西 延安 716000；3.西北農林科技大學 水土保持研究所 黃土高原土壤侵蝕與旱地農業國家重點實驗室，陜西 楊凌 712100；4.乾縣農業農村局，陜西 乾縣 713300）

土壤侵蝕是全球面臨的突出生態環境問題之一［1］,防治水土流失是黃河流域生態治理的根本措施［2］。1999 年開始實施退耕還林(草)工程以來,在氣候變化不顯著的背景下,黃河中游地區土壤侵蝕與入黃泥沙狀況發生了顯著變化［3-4］。土壤侵蝕是各種因子綜合作用的結果,研究土壤侵蝕時空變化及其與影響因子的關系對防治土壤侵蝕具有重要意義［5］。北洛河流域是黃河中游粗泥沙主要來源區之一,也是全國實施水土流失治理和退耕還林(草)工程的重點地區,筆者采用通用土壤流失方程RUSLE,結合GIS 和RS 技術,定量估算北洛河上游流域實施退耕還林(草)工程前后不同時段土壤侵蝕模數和土壤侵蝕量,分析土壤侵蝕變化與土地利用、植被覆蓋、地形等的關系,以期為該區域水土保持生態建設及相關決策提供依據。

1 數據與方法

1.1 研究區概況

北洛河流經陜西、甘肅兩省的5 個地(市)18 個縣(區),于陜西省大荔縣注入渭河,屬黃河二級支流。北洛河流域屬溫帶大陸性季風氣候區,多年平均降水量459.1 mm,5—9 月降水量占全年的74.5%,降水空間分布不均,東南部降水較多、西北部降水較少。吳旗水文站控制的北洛河上游流域屬黃土丘陵溝壑區,是黃河中游粗泥沙集中來源區之一,流域面積3 408 km2,土壤以黃綿土為主,土層深厚、土質疏松、抗蝕性較差,溝谷密度為3.0～4.5 km/km2［6］,梁峁起伏、溝壑縱橫,水力侵蝕極為強烈,年均土壤侵蝕模數7 006 t/(km2·a),河源區一帶土壤侵蝕模數可達10 000 t/(km2·a)。自1999 年開始實施退耕還林(草)工程以來,流域內的林草覆蓋面積顯著增加,形成了以落葉闊葉及灌木草叢為主的次生植被類型［7］,因溝道比降較大而淤地壩工程較少且規模較小,有效攔沙庫容非常有限,目前大都基本淤滿［8］,因此人工植被建設是北洛河上游流域水土保持生態建設的主要措施。

1.2 數據來源

本研究采用的主要數據及其來源:①研究區域及周邊分布的吳起、鹽池、定邊、靖邊、志丹、安塞、華池、環縣等8 個雨量站1981—2010 年逐日降水量數據,來源于國家氣象信息中心；②吳旗水文站1981—2010 年逐日徑流泥沙數據,來源于黃河水文年鑒及中國科學院水土保持研究所館藏數據；③1 ∶50 萬土壤屬性數據,下載于黃土高原科學數據中心(http://loess.data.ac.cn/),不同土壤類型(黃綿土、黑壚土、紅土、新積土)的理化性質數據從《陜西土壤》［9］中查得；④地形數據(分辨率為30 m 的ASTER GDEM),下載于地理空間數據云(http://www.gscloud.cn/)；⑤1990 年、2000 年、2010 年Landsat TM5 遙感影像數據(分辨率為30 m),來源于地理空間數據云；⑥1990 年和2000年土地利用(土地利用類型分為耕地、草地、林地、水域、建設用地)數據下載于黃土高原科學數據中心,2010 年土地利用數據通過解譯Landsat TM5 遙感影像并經實地調查驗證后得到［10］。

1.3 通用土壤流失方程

在眾多土壤侵蝕模型中,通用土壤流失方程(RUSLE)因結構簡單、計算簡便、精度較高而得到廣泛應用［11-13］,其形式為

式中:A為土壤侵蝕模數,t/(hm2·a)；R為降雨侵蝕力因子,MJ·mm/(hm2·h·a)；K為土壤可蝕性因子,t·hm2·h/(hm2·MJ·mm)；L為坡長因子、S為坡度因子,通常綜合考慮坡長、坡度對土壤侵蝕模數的影響,把LS稱為地形因子；C為植被覆蓋因子；P為水土保持措施因子。

(1)降雨侵蝕力因子R值的確定。采用逐日降雨資料估算年降雨侵蝕力的簡易算法［14］,得到研究區8個雨量站1981—2010 年降雨侵蝕力,然后通過克里金法進行空間插值,得到研究區逐年降雨侵蝕力。根據研究區1981—2010 年降雨侵蝕力與吳旗水文站年輸沙量雙累積曲線斜率變化情況,可知該曲線在2000 年出現較明顯拐點［15］?？紤]到降雨侵蝕力年際波動較大的情況,將研究期1981—2010 年分為1981—1990年、1991—2000 年和2001—2010 年3 個時段(編號分別為I、Ⅱ、Ⅲ,把1990 年、2000 年、2010 年作為3 個時段的代表年)。3 個時段研究區年均降雨侵蝕力分別為1 093.72、1 390.20、1 260.74 MJ·mm/(hm2·h·a),表明1981—2010 年研究區降雨侵蝕力呈現波動上升趨勢。

(2)土壤可蝕性因子K值的確定。根據第一次全國水利普查所得水土保持專項調查成果(土壤有機碳、土壤顆粒組成數據等)［16］,采用Williams 等［17］提出的公式進行計算,并基于黃土高原徑流小區資料進行修正后得到研究區土壤可蝕性因子K值為0.039 8～0.055 2 t·hm2·h/(hm2·MJ·mm),平均值為0.051 7 t·hm2·h/(hm2·MJ·mm)。

(3)地形因子LS值的確定。根據研究區DEM 數據,采用張宏鳴等［18］基于GIS 提取區域坡長坡度因子的算法進行計算,并經過必要的尺度變換修正［19］,得到研究區LS值范圍為0.01～71.61,平均值為12.72。

(4)植被覆蓋因子C值的確定。根據研究區土地利用數據和基于歸一化植被指數NDVI分布圖,在前人研究的基礎上,得到研究區不同土地利用類型及其植被覆蓋度的C值,其中:耕地C值為0.44(根據研究區常見旱地農作物谷子、玉米、大豆、馬鈴薯C值分別為0.53、0.28、0.51、0.47［20］,其種植面積比例為1 ∶2 ∶1 ∶6,經加權平均計算得到)；林地、草地的C值按照江忠善等［21］建立的黃土丘陵區人工林、草地植被覆蓋度與C值的關系式進行計算(非固定值)；建設用地和水域C值為0。利用ArcGIS 軟件得到研究區C值的分布。1990 年、2000 年、2010 年研究區平均C值分別為0.548、0.402、0.222,呈現持續減小趨勢,原因是研究區植被覆蓋度持續提高從而使C值持續減小。

(5)水土保持措施因子P值的確定。通過Landsat TM5 遙感影像難以提取水土保持措施,故利用研究區土地利用類型數據,通過賦值的方法對水土保持措施因子進行賦值。參考前人研究成果［10］,對耕地、林地、草地P值分別賦值為0.31、0.05、0.16,對建設用地及水域P值均賦值為1,然后利用ArcGIS 軟件獲得研究區P值的分布。1990 年、2000 年、2010 年研究區P值分別為0.220 7、0.220 5、0.157 0,相比較而言,2010 年較2000 年明顯減小,原因是2000 年以后研究區耕地被大面積轉化為林地、草地。

2 結果與分析

2.1 土壤侵蝕時空變化情況

按前述RUSLE 分別計算出研究區時段Ⅰ、Ⅱ、Ⅲ土壤侵蝕模數［單位換算為t/(km2·a)］,依據《土壤侵蝕分類分級標準》(SL190—2007),利用ArcGIS 軟件獲得研究區3 個時段各級土壤侵蝕空間分布情況(見圖1),3 個時段各級土壤侵蝕面積及占比統計見表1。3 個時段平均土壤侵蝕模數分別為8 612.62、8 078.91、3 060.71 t/(km2·a),呈持續減小趨勢。由表1 可知,3 個時段中度土壤侵蝕面積及占比差別相對較小,而強烈及以上土壤侵蝕面積持續減小、占比持續降低(時段I 面積為2 052.90 km2、占比為60.24%,時段Ⅱ面積為1 925.59 km2、占比降為56.50%,時段Ⅲ面積減小到645.25 km2、占比降至18.93%),輕度及微度土壤侵蝕面積持續增大、占比持續提高。時段Ⅲ強烈及以上土壤侵蝕面積大幅減小、輕度及微度土壤侵蝕面積顯著增大,這與1999 年開始實施的退耕還林(草)工程密不可分,退耕還林(草)工程的實施使研究區植被覆蓋度大幅提升,尤其研究區東南部的吳起縣植被恢復最為明顯,因而土壤侵蝕模數在空間上呈現西北高、東南低的空間分異特征。

表1 研究區3 個時段各級土壤侵蝕面積及占比統計

圖1 研究區各時段土壤侵蝕強度分布情況

2.2 不同土地利用類型的土壤侵蝕變化情況

將各時段土壤侵蝕強度分布圖與土地利用圖進行疊加,可得到各時段不同土地利用類型的土壤侵蝕狀況。鑒于3 個時段水域、建設用地面積占比均較小(分別為0.08%～0.51%、0.06%～0.52%),且本研究在對水域、建設用地進行土壤侵蝕判讀時認定為不發生土壤侵蝕,因此只統計了研究區3 個時段耕地、林地、草地的面積占比及土壤侵蝕模數,見表2 。

表2 研究區各時段主要土地利用類型土壤侵蝕模數

由表2 可知,時段I 和時段Ⅱ土地利用類型以耕地和草地為主,二者約占流域面積的97%；時段Ⅲ土地利用格局發生根本性變化,呈現以草地和林地為主的格局,林地、草地面積占比分別提高至20.95%和63.43%,而耕地面積占比由42.62%～42.69%銳減至14.59%,這種變化是實施退耕還林(草)工程的結果。

研究區3 種主要土地利用類型土壤侵蝕模數大小為耕地＞草地＞林地,從土壤侵蝕模數變化情況來看,從時段Ⅰ到時段Ⅱ耕地土壤侵蝕模數小幅提高、草地土壤侵蝕模數小幅降低、林地土壤侵蝕模數明顯降低,從時段Ⅱ到時段Ⅲ耕地、草地、林地土壤侵蝕模數均明顯降低。

2.3 不同植被覆蓋度的土壤侵蝕變化情況

把植被覆蓋度分為低覆蓋度(小于10%)、較低覆蓋度(10%～30%)、中覆蓋度(30%～50%)、較高覆蓋度(50%～70%)、高覆蓋度(大于70%)5 個級別［22］,各時段土壤侵蝕強度分布圖與各級植被覆蓋度分布圖疊加,可得3 個時段各級覆蓋度的土壤侵蝕模數。由表3 可知,時段I 植被覆蓋度以低覆蓋度和較低覆蓋度為主(二者面積占比達89.83%),時段Ⅱ植被覆蓋度以較低覆蓋度和中覆蓋度為主(二者面積占比為86.86%),時段Ⅲ低覆蓋度和較低覆蓋度面積占比明顯降低、較高覆蓋度和高覆蓋度面積占比明顯提升,植被覆蓋度的變化與實施退耕還林(草)工程后生態環境質量持續向好相吻合。隨著植被覆蓋度的提高,土壤侵蝕模數明顯降低,原因是植被覆蓋度越高抵抗侵蝕的作用越強。從土壤侵蝕模數變化情況來看,從時段I 至時段Ⅱ,在降雨侵蝕力增大的背景下同一植被覆蓋度的土壤侵蝕模數有所提高；從時段Ⅱ至時段Ⅲ,相同植被覆蓋度的土壤侵蝕模數明顯降低,這一方面與氣候變化引起降雨侵蝕力有所減弱有關,另一方面主要與植被覆蓋度上升有關。

表3 研究區3 個時段各級植被覆蓋度土壤侵蝕模數

2.4 不同地形的土壤侵蝕變化情況

2.4.1 不同高程的土壤侵蝕變化情況

結合流域實際情況對高程進行分級,將高程分級圖與各時段土壤侵蝕強度分布圖進行疊加,統計得到各時段不同高程的土壤侵蝕模數(見表4)。從各級高程面積占比來看,1 350～1 450、1 450～1 550 m 這兩級高程的面積占比達76.78%,為研究區主要高程帶。從各級高程的土壤侵蝕模數來看,1 350～1 550 m 高程帶是研究區土壤侵蝕的主要高程帶(土壤侵蝕量約占研究區土壤侵蝕總量的80%)。從各級高程土壤侵蝕模數的變化來看,實施退耕還林(草)工程后的時段Ⅲ土壤侵蝕模數較時段Ⅰ和時段Ⅱ均大幅度降低(平均降幅為62.26%)。

表4 研究區各時段不同高程的土壤侵蝕模數

2.4.2 不同坡角的土壤侵蝕變化情況

把研究區坡角分級(0°～5°、5°～8°、8°～15°、15°～25°、25°～35°、＞35°)圖與各時段土壤侵蝕強度分布圖進行疊加,統計得到各時段不同坡角的土壤侵蝕模數(見表5)。土壤侵蝕模數隨坡角增大而提高,即坡面越陡土壤侵蝕模數越大。從各級坡角坡面面積占比和土壤侵蝕模數綜合來看,15°以上坡面是研究區土壤侵蝕的主要區域(土壤侵蝕量約占研究區土壤侵蝕總量的72%)。從各時段土壤侵蝕模數的變化情況來看,時段Ⅲ各級坡角的土壤侵蝕模數較時段Ⅰ和時段Ⅱ顯著降低(平均降幅為60.53%)。

表5 研究區各時段不同坡角的土壤侵蝕模數

2.4.3 不同坡向的土壤侵蝕變化情況

按水土流失常規調查中的四分法將坡向分為陰坡、半陰坡、陽坡、半陽坡,把坡向分布圖與各時段土壤侵蝕強度分布圖進行疊加,統計得到各時段不同坡向的土壤侵蝕模數(見表6)。各時段不同坡向的土壤侵蝕模數均呈現陽坡＞半陽坡＞半陰坡＞陰坡,原因是坡向影響水、熱分配,從而影響植被生長狀況,進而導致不同坡向土壤侵蝕強度存在差異,陰坡相對陽坡來講,土壤水分含量大,植被生長迅速,地表覆蓋度高,因而土壤侵蝕模數相對較小。從各時段土壤侵蝕變化來看,時段Ⅲ各坡向的土壤侵蝕模數較時段Ⅰ和時段Ⅱ均顯著降低(平均降幅為60.49%)。

表6 研究區各時段不同坡向的土壤侵蝕模數

2.5 土壤侵蝕量變化情況

依據上述各時段主要土地利用類型面積占比和土壤侵蝕模數計算的年均土壤侵蝕量見表7。從各時段年均土壤侵蝕量變化情況來看,從時段I 到時段Ⅱ,耕地的年均土壤侵蝕量變化很小(約增大0.2%),占土壤侵蝕總量的比例由57%提高到61%；草地的年均土壤侵蝕量減小14%,其占土壤侵蝕總量的比例由42%下降到39%；林地的年均土壤侵蝕量減小40%,其占土壤侵蝕總量的比例由0.6%下降到0.4%。從時段Ⅱ到時段Ⅲ,耕地的年均土壤侵蝕量減小77%,其占土壤侵蝕總量的比例由61%下降到37%；草地的年均土壤侵蝕量減小40%,其占土壤侵蝕總量的比例由39%提高到62%；林地的年均土壤侵蝕量增大58%,其占土壤侵蝕總量的比例由0.4%提高到1.6%。時段Ⅲ耕地土壤侵蝕量及其占比大幅度下降、林草地土壤侵蝕量及其占比有所提高的原因是,實施退耕還林(草)工程后耕地面積顯著減小、林草地面積顯著增大。

表7 各時段主要土地利用類型年均土壤侵蝕量 萬t/a

2.6 討論

研究結果顯示,時段Ⅰ、時段Ⅱ、時段Ⅲ研究區平均土壤侵蝕模數分別為8 612.62、8 078.91、3 060.71 t/(km2·a),可以看出,實施退耕還林(草)工程前從時段Ⅰ到時段Ⅱ土壤侵蝕模數稍有降低,而實施退耕還林(草)工程后的時段Ⅲ土壤侵蝕模數大幅度下降,表明退耕還林(草)工程的實施在北洛河上游流域取得了顯著的水土保持效益,這與劉文超等［23］對陜北地區退耕還林還草工程保土效應的研究結果相一致。退耕還林(草)工程實施以后,北洛河上游流域土壤侵蝕強度在空間上呈現西北高、東南低的特征,原因是東南部地勢較低,氣候較為濕潤,植被恢復效果明顯,而西北部多為河源區,山高坡陡,水熱條件較差,植被恢復難度大甚至在局部地區植被略有退化［24］,因此植被恢復效果較好的東南部土壤侵蝕模數明顯低于西北部土壤侵蝕模數。

研究區主要土地利用類型土壤侵蝕模數遵循耕地＞草地＞林地的規律,退耕還林(草)工程的實施,使土地利用格局與植被覆蓋度發生顯著變化(耕地面積占比大幅度下降、林草地面積占比顯著提高),土壤侵蝕模數隨植被覆蓋度的提高而降低,這與有關學者［10,25］的研究結果一致。值得注意的是,耕地是土壤侵蝕最為嚴重的土地利用類型,其土壤侵蝕模數在退耕還林(草)工程實施前的時段Ⅰ為11 481.84 t/(km2·a)、時段Ⅱ為11 528.87 t/(km2·a),時段Ⅱ較時段I 高的原因是降雨侵蝕力增大；在退耕還林(草)工程實施后的時段Ⅲ,土壤侵蝕模數大幅度降低到7 745.10 t/(km2·a),其原因除時段Ⅲ降雨侵蝕力略有減小外,主要是隨著退耕還林(草)工程的實施耕地的平均坡度減緩,即土壤侵蝕嚴重的陡坡耕地轉化為林草地,而保留的耕地為土壤侵蝕強度較低的緩坡地。

地形是影響土壤侵蝕的主要因素之一,不同時段土壤侵蝕模數與各地形因子的關系表現出一致性,研究區土壤侵蝕主要發生在1 350～1 550 m 高程帶、坡角＞15°的陡坡地、陽坡和半陽坡,這些地帶是研究區今后土壤侵蝕防治的重點區域。

土壤侵蝕模數的變化主要受氣候條件、土壤、地形、植被覆蓋度等因子的影響。土壤和地形為相對靜態因子,決定了土壤侵蝕模數的空間分異特征；氣候、植被覆蓋度為相對動態因子,影響土壤侵蝕模數隨時間的動態變化［26］。統計分析表明,1981—2010 年北洛河上游流域降雨侵蝕力年際變化較大,3 a 滑動平均處理結果顯示,降雨侵蝕力呈增大趨勢,線性傾向率為7.3 MJ·mm/(hm2·h·a)［15］,而土壤侵蝕模數呈下降趨勢,尤其是時段Ⅲ土壤侵蝕模數急劇下降,這主要歸因于退耕還林(草)工程的實施使植被覆蓋度顯著提高、土壤侵蝕得到有效控制。

計算的時段Ⅰ、時段Ⅱ、時段Ⅲ年均土壤侵蝕量分別為2 946.4 萬、2 763.0 萬、1 047.6 萬t/a,吳旗水文站實測的時段Ⅰ、時段Ⅱ、時段Ⅲ年均輸沙量分別為2 690.4萬、4 559.0 萬、1 033.8 萬t/a,二者相關性較好(相關系數為0.64),但存在一定差別,二者在時段I 和時段Ⅲ比較接近,而時段Ⅱ計算的土壤侵蝕量遠小于實測輸沙量,其原因:一是采用RUSLE 計算的土壤侵蝕量是坡面侵蝕量,沒有考慮溝蝕、重力侵蝕以及泥沙輸移過程中的沉積等；二是時段Ⅱ期間的1994 年8 月31 日北洛河流域出現了超過50 a 一遇的特大暴雨,產生了高含沙洪水,該年輸沙量是多年平均輸沙量的7.3倍［27］,具有明顯的峰值特征,因此時段Ⅱ實測年均輸沙量遠大于計算的年均土壤侵蝕量。

總體而言,采用RUSLE 模型估算的研究區土壤侵蝕模數較為合理。然而,受所采用數據來源及精度的制約,估算結果可能與實際情況存在較大差異(如時段Ⅱ尤其是1994 年)。受Landsat TM5 遙感影像分辨率的制約,本研究在無法獲取研究區水土保持措施空間分布信息的情況下,采用土地利用數據間接反映水土保持措施對土壤侵蝕的抑制作用,其與實際情況可能存在差異。相關研究中,大中流域尺度地形因子大都是基于開放免費的低精度DEM 數據提取的［28］,這類數據對地形的表達有限甚至存在地形描述誤差,據此提取地形因子存在坡度衰減、坡長增長的問題［29］,并影響土壤侵蝕模數估算精度。因此,今后應加強相關研究,以提高土壤侵蝕模型的精度。

3 結論

(1)北洛河上游流域土壤侵蝕模數在空間上呈現西北高、東南低的分異特征,1981—2010 年在降雨侵蝕力呈現波動上升趨勢的情況下,流域平均土壤侵蝕模數呈持續下降趨勢,尤其退耕還林(草)工程實施以后土壤侵蝕模數顯著下降,時段I、時段Ⅱ、時段Ⅲ土壤侵蝕模數分別為8 612.62、8 078.91、3 060.71 t/(km2·a)。

(2)實施退耕還林(草)工程以后,研究區耕地面積銳減,林草地面積顯著增加、植被覆蓋度大幅提升,主要土地利用類型土壤侵蝕模數遵循耕地＞草地＞林地的規律,即土壤侵蝕模數隨著植被覆蓋度的提高呈現下降趨勢。

(3)3 個時段土壤侵蝕隨地形的分布特征具有一致性,1 350～1 550 m 是研究區土壤侵蝕分布的主要高程帶,其侵蝕量約占研究區土壤侵蝕總量的80%,研究區72%的土壤侵蝕量分布在坡角＞15°的坡面,各時段不同坡向的平均土壤侵蝕模數均遵循陽坡＞半陽坡＞半陰坡＞陰坡的規律。