郭子萍, 王乃昂, 屈志勇

(1.蘭州大學(xué) 資源環(huán)境學(xué)院/干旱區(qū)與沙漠研究中心, 甘肅 蘭州 730000; 2.山西省氣象信息中心, 山西 太原 030000)

土壤侵蝕是當(dāng)今人類面臨的一種最普遍、持續(xù)性最強的地質(zhì)災(zāi)害[1-3]，會引起土地退化、土壤肥力下降、水環(huán)境惡化，從而影響生態(tài)安全和人類的可持續(xù)發(fā)展。黃土高原是中國水土流失最嚴重的地區(qū)[4]，而土壤侵蝕是黃土高原河流泥沙的主要來源[5]，因此，定量計算土壤侵蝕量、土壤侵蝕強度，評價水土流失時空分布特征，對該地區(qū)水土保持及土地的合理利用具有重要的意義。土壤侵蝕模型是土壤侵蝕量定量評價和預(yù)報的重要方法。國內(nèi)外已經(jīng)研究開發(fā)出許多實用的模型，其中修正的通用土壤流失方程(RUSLE)應(yīng)用最為廣泛[6-8]。

自20世紀90年代以來，在GIS和RS技術(shù)的支持下，中國學(xué)者利用RUSLE模型開展了廣泛的土壤侵蝕研究，其中黃土高原作為研究的熱點區(qū)域受到了廣泛的關(guān)注，但相關(guān)研究大都集中于對侵蝕面積、侵蝕總量、空間分布的分析，對與影響因子的關(guān)系以及對人為因素作用日漸增強的地區(qū)關(guān)注度不夠。引黃入晉工程是為促進地區(qū)經(jīng)濟發(fā)展，滿足能源基地供水而建設(shè)的水利工程，對煤電基地建設(shè)、提高人民生活質(zhì)量具有重要意義，但也不可避免的加重了水土流失問題。目前對于該工程的研究主要集中在水價分析[9-10]和社會經(jīng)濟影響評價研究[11]等社會經(jīng)濟領(lǐng)域；抗震性分析[12]、泥質(zhì)膨脹巖工程地質(zhì)研究[13]等工程地質(zhì)領(lǐng)域；環(huán)境保護設(shè)計[14]等環(huán)境工程領(lǐng)域，針對土壤侵蝕的研究較少，幾乎沒有涉及。

本研究擬采用RUSLE模型計算引黃入晉北干線2005，2010，2015年3 a土壤侵蝕模數(shù)，分析土壤侵蝕時空變化特征，以期為當(dāng)?shù)厮帘３趾蜕鷳B(tài)環(huán)境建設(shè)提供科學(xué)理論依據(jù)。

1 研究區(qū)概況與數(shù)據(jù)收集

1.1 研究區(qū)概況

山西省萬家寨引黃工程位于山西省西北部，北干線主要任務(wù)是自萬家寨水庫引水，向山西北部的大同、平朔等能源基地供水，年引水量2.96×108m3，供水區(qū)總供水面積5 273 km2[14]。研究區(qū)屬溫帶大陸性季風(fēng)氣候，年均氣溫5.7～7.9 ℃，年平均降水量380～440 mm,降水季節(jié)分配不均，降水變率大。研究區(qū)海拔為859～2 405 m，地形起伏變化較大[15]。該工程北干線2009年2月開始實施，2011年9月16日竣工。

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

本研究所使用的數(shù)據(jù)包括： ①朔州市市轄區(qū)、山陰縣、懷仁縣、大同市市轄區(qū)2005，2010，2015年月降水量數(shù)據(jù)，來源為山西省氣象信息中心。②研究區(qū)GDEMV230 M分辨率數(shù)字高程模型數(shù)據(jù)以及研究區(qū)2005，2010，2015年LandsatTM4-5遙感影像圖，均來源于地理空間數(shù)據(jù)云。③研究區(qū)2005，2010，2015年土地利用圖、研究區(qū)土壤類型圖，來源于中國科學(xué)院資源環(huán)境科學(xué)數(shù)據(jù)中心。

2 研究方法

2.1 修正的通用土壤流失方程

修正的通用土壤流失方程(revised universal soil loss equation, RUSLE)作為估算年均土壤流失量的方法，需要考慮降水(月降水量)、土壤類型、坡度因子、坡長因子、植被覆蓋與管理因子、水土保持措施管理因子6大主要因子，方程定義為[16]：

A=R·K·L·S·C·P

(1)

式中：A——土壤侵蝕模數(shù)〔t/(hm2·a)〕;R——降雨侵蝕力因子〔MJ·mm/(hm2·h·a)〕；K——土壤可蝕性因子〔t·hm2·h/(hm2·MJ·mm)〕；L——坡長因子(無量綱)；S——坡度因子(無量綱)；C——植被覆蓋與管理因子(無量綱)；P——水土保持措施因子(無量綱)。

2.1.1 降雨侵蝕力因子R降雨侵蝕力因子是預(yù)測預(yù)報土壤侵蝕模型中的重要因子，本研究采用Wischmeier等[17]提出并經(jīng)Arnoldus修改的比較適用于黃土高原地區(qū)且廣泛應(yīng)用的經(jīng)驗公式(2)來計算研究區(qū)2005，2010，2015年的R值。

(2)

式中：P——年降水量(mm)；Pi——第i月降水量(mm)[17]。

利用2005，2010，2015年北干線沿線完整且均勻分布的4個氣象站點的氣象數(shù)據(jù)，根據(jù)公式(1)計算得到研究區(qū)3 a的R值，利用ArcGIS工具連接功能將R值賦給各站點，并利用空間插值對R值進行空間化處理。研究采用普通克里格方法，半變異模型選用球形函數(shù)模型，進行研究區(qū)R值空間插值，生成研究區(qū)3 a的R值空間分布。

2.1.2 土壤可蝕性因子K土壤可蝕性因子K(soil erodibility)是土壤抵抗侵蝕能力的綜合體現(xiàn)，影響因素有：土壤結(jié)構(gòu)、土壤質(zhì)地、土壤滲透率、緊實度、含水率以及黏土礦物性質(zhì)等[18]，根據(jù)中國科學(xué)院資源環(huán)境科學(xué)數(shù)據(jù)中心公布的土壤類型數(shù)據(jù)，獲得研究區(qū)土壤類型，以及《山西土壤》[19]中記載的關(guān)于各種土壤類型典型剖面的理化性質(zhì)，采用Wischmeier等[7]于1990年在EPIC模型中建立的計算公式(3)計算研究區(qū)各土壤類型的K值。

(3)

式中：SAN——砂粒含量(%); SIL——粉砂含量(%); CLA——黏粒含量(%);C——有機碳含量(%)，SN1=1-SAN/100。

張科利等[20]的研究指出，與美國相比，中國土壤可蝕性K值普遍偏小，用EPIC模型所建立的公式計算的可蝕性K值，通過公式(4)進行修正轉(zhuǎn)換，方可用于中國的土壤可蝕性估算。

K=0.013 83+0.515 75

(4)

利用公式(3),(4)計算的研究區(qū)K值，并在ArcGIS中生成K值分布。

2.1.3 坡度因子S通用土壤流失方程中允許計算的最大坡度為18%(10°)，研究區(qū)位于黃土高原，這樣的條件不適用于研究區(qū)丘陵溝壑中的陡坡地形(>10°)，因此，在利用DEM提取坡度以后，要分段對坡度因子S進行計算，緩坡采用McCool坡度公式，陡坡采用劉寶元[21]提出的黃土高原S因子計算公式(5)。

(5)

式中：S——坡度因子；θ——坡度。

2.1.4 坡長因子L采用Hickey提出的非累計流量直接計算法估算坡長，坡長因子L的提取采用Wischmeier提出的坡長因子計算公式(6)：

L=(λ/22.13)m

m=β/(1-β)

(6)

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

式中：λ——水平坡長(m)；m——坡長指數(shù)； 22.13——標準徑流小區(qū)坡長(m)；β——細溝侵蝕與面蝕的比率；θ——利用DEM提取的坡度。ArcGIS計算的坡度以度為單位，要轉(zhuǎn)換為弧度表示，即乘以π/180。

2.1.5 植被覆蓋與管理因子CNDVI是估算C因子最普遍應(yīng)用的數(shù)據(jù)[22]。以經(jīng)過輻射定標和大氣校正的TM影像為數(shù)據(jù)，在ENVI 5.1中計算歸一化植被指數(shù)NDVI，在提取NDVI之后，根據(jù)像元二分模型原理[23]，計算研究區(qū)的植被覆蓋度F，公式為：

F=(NDVI-NDVIsoil)/(NDVIveg-NDVIsoil)

(7)

式中：F——有植被覆蓋的面積比例，即像元的植被覆蓋度； NDVIveg——純植被覆蓋像元的NDVI值； NDVIsoil——無植被覆蓋像元的NDVI值。在求解植被覆蓋度F后，參考有關(guān)研究成果，賦予了研究區(qū)不同土地利用類型和不同植被覆蓋度下的C值(見表1)。

表1 研究區(qū)不同土地利用類型和不同植被蓋度下的C值

2.1.6 水土保持措施因子Pp值介于0～1之間，0代表根本不會發(fā)生土壤侵蝕的地區(qū)，1代表未采取任何水土保持措施的地區(qū)。國內(nèi)學(xué)者根據(jù)區(qū)域特點對不同的土地利用方式進行賦值，本研究區(qū)位于黃土高原，參照前人有關(guān)研究成果[24-25]賦值P，將水體、林地、草地、建筑用地、未利用地賦值為1，對于耕地，通常坡度越大，水土保持措施的作用越突出，因此耕地依據(jù)表2按照坡度范圍進行賦值。

表2 不同坡度范圍內(nèi)耕地的p值

2.2 土壤侵蝕量估算

將上述各因子圖層?xùn)鸥駟卧y(tǒng)一為1 000 m×1 000 m,并標準化為相同的投影坐標系統(tǒng)WGS1 984_UTM_zone_49 N，地理坐標系GCS_WGS_1984,在ArcGIS軟件的柵格計算器工具中將每個網(wǎng)格(Grid)中各因子連乘，得到研究區(qū)3期土壤侵蝕模數(shù)。RUSLE計算的土壤流失量默認單位是美制單位〔ton/(area·a)〕，要將結(jié)果乘以224.2得到中國常用的公制單位t/(km2·a)。

3 結(jié)果與分析

3.1 土壤侵蝕強度分級

在ArcGIS中根據(jù)水利部《土壤侵蝕分類分級標準(SL190-2007)》，對研究區(qū)土壤侵蝕模數(shù)進行重分類，得到3期土壤侵蝕強度圖如附圖2所示。

3.2 土壤侵蝕強度時空動態(tài)變化

由表3可知，從土壤侵蝕模數(shù)看，研究區(qū)2005，2015年的土壤侵蝕模數(shù)均在500～2 500 t/(km2·a)之間，屬于輕度侵蝕，而2010年土壤侵蝕模數(shù)為3 340.41 t/(km2·a)，屬于中度侵蝕。

表3 研究區(qū)土壤侵蝕變化統(tǒng)計

從2005—2010年，土壤侵蝕模數(shù)增長了80.91%；從2010—2015年下降了47.87%。3 a的土壤侵蝕模數(shù)呈先增后減，總體減小的趨勢。但3 a中土壤侵蝕模數(shù)均大于水利部公布的黃土高原土壤允許流失量1 000 t/(km2·a)的標準。這說明研究區(qū)土壤侵蝕程度嚴重，應(yīng)重點給予關(guān)注。引黃入晉北干線2009年2月實施，2011年9月竣工，2010年土壤侵蝕模數(shù)驟增與該工程的實施關(guān)系密切。工程建成后，為了區(qū)域的生態(tài)安全和經(jīng)濟的可持續(xù)發(fā)展，實施了一系列生物措施和工程措施，產(chǎn)生了良好的水土保持效果，因此2015年土壤侵蝕模數(shù)減小。

從圖1中可以看出，2005,2015年均以微度侵蝕所占的面積最大，其次為輕度、中度、強度、極強度、劇烈侵蝕，土壤侵蝕面積隨著侵蝕強度的增加而逐漸減小。與2005年相比，2010年由微度侵蝕轉(zhuǎn)移到中度以上侵蝕的面積高達11.75%。中度以上侵蝕強度所占面積增大，無疑會加重土壤侵蝕。2015與2010年相比，微度、輕度侵蝕面積增加，中度以上侵蝕所占面積減小12.63%，侵蝕狀況明顯好轉(zhuǎn)，且各侵蝕強度所占的面積與2005年基本相同。

圖1 引黃入晉北干線沿線地區(qū)不同年份各侵蝕強度所占面積比

由表4可知，從空間上看，研究區(qū)土壤侵蝕呈塊狀侵蝕和點狀侵蝕相間分布的特點。3 a中土壤侵蝕面積均為朔州市市轄區(qū)最大，塊狀侵蝕分布較多，其次是大同市轄區(qū)、山陰、懷仁，且3 a中土壤侵蝕量也均為朔州市轄區(qū)最大。從土壤侵蝕模數(shù)看，4個地區(qū)土壤侵蝕模數(shù)從2005—2015年均呈先增后減的趨勢，在2010年處達到峰值，因此每個地區(qū)的土壤侵蝕模數(shù)與總的土壤侵蝕模數(shù)分布特點一致。

由附圖2可以看出，以北干線為界，線上方土壤侵蝕明顯高于線下方，兩側(cè)侵蝕狀況差異顯著，這種趨勢在2010年尤為明顯。土壤侵蝕的空間分布與研究區(qū)的DEM也有較高的吻合度，分析研究區(qū)DEM可以明顯看出，研究區(qū)分為2個地貌單元，分別為侵蝕中低山區(qū)和沖洪積平原區(qū)[15]，侵蝕中低山區(qū)土壤侵蝕明顯比沖洪積平原侵蝕嚴重。

3.3 土壤侵蝕主控因子分析

3.3.1 土壤侵蝕與降雨 從圖2可以看出，研究區(qū)不同地域單位R因子均在2010年處達到最大值，而折線圖顯示，土壤侵蝕模數(shù)也在2010年處達到峰值。在其他因素不變的情況下，降雨侵蝕力因子R是土壤侵蝕的加速因子。因此，當(dāng)R因子的值增大時，土壤侵蝕模數(shù)也會相應(yīng)的增大。分析所獲得的降水?dāng)?shù)據(jù)，研究區(qū)2005年年降水量為367.73 mm，2010年為423.5 mm，2015年為380 mm，由此可以看出2010年降水量顯著增大是導(dǎo)致R因子值增大的直接原因，R因子與土壤侵蝕模數(shù)呈明顯正相關(guān)。

3.3.2 土壤侵蝕與土壤類型 由表5可以看出，栗褐土所占的土壤侵蝕面積最大，為4 199 km2，栗褐土的土壤侵蝕量也最高。從研究區(qū)土壤類型分布圖中，我們可以清楚的看到，研究區(qū)栗褐土的面積為46.36%，栗褐土區(qū)為半干旱一年一熟旱作區(qū)，分布有大量的耕地，但是耕地的質(zhì)量較差，絕大部分為坡耕地，且坡度大的耕地占很大比例，耕地的輸出少，只用不養(yǎng)，導(dǎo)致農(nóng)田生態(tài)惡化，土壤干旱、貧瘠，侵蝕嚴重，再生能力低。因此在栗褐土分布區(qū)，應(yīng)重點防范水土流失。

圖2 引黃入晉北干線沿線地區(qū)各年份R因子值與土壤侵蝕模數(shù)

表4 引黃入晉北干線沿線地區(qū)各年份不同地域土壤侵蝕變化統(tǒng)計

表5 引黃入晉北干線沿線地區(qū)各年份不同土壤類型土壤侵蝕變化統(tǒng)計

3.3.3 土壤侵蝕與坡度 以全國《土壤侵蝕強度劃分標準(SL190-2007)》為坡度分類標準，將坡度重分類為6類：<5°，5°～8°，8°～15°，15°～25°，25°～35°，>35°。據(jù)相關(guān)文獻可知，地表類型按坡度可以分為平坡(0～5°)、緩坡(5°～15°)、斜坡(15°～25°)、陡坡(25°～35°)、急坡(35°～45°)、險坡(≥45°)[26]，根據(jù)面積分布情況，可以看出，研究區(qū)以平坡和緩坡為主，斜坡為輔，急坡和險坡占很小比重。為了進一步探討坡度與土壤侵蝕量之間的關(guān)系，利用空間疊加工具將坡度圖層與侵蝕模數(shù)圖層疊加，計算研究區(qū)2005，2010，2015年6類坡度下的土壤侵蝕量模數(shù)及土壤侵蝕量比如表6所示。

表6 引黃入晉北干線沿線地區(qū)各年份不同坡度土壤侵蝕變化統(tǒng)計

從表6中可以看出，2005，2010，2015年在不同坡度上的土壤侵蝕格局變化不大，不同坡度與土壤侵蝕關(guān)系密切。

從侵蝕面積來看，8°～15°坡度帶最大，這主要與研究區(qū)緩坡分布面積大有關(guān)。從平均侵蝕模數(shù)看，當(dāng)坡度大于25°時土壤侵蝕模數(shù)較高，盡管此帶的分布面積小，但侵蝕量大，因此侵蝕模數(shù)高，侵蝕級別2005，2015年達到強度侵蝕，2010年為極強度侵蝕，主要是因為在北方，25°以下的區(qū)域?qū)儆诰徠聟^(qū)域，退耕還林工程的實施使得該區(qū)域成為人工林和梯田改造后的農(nóng)田，土壤侵蝕狀況稍有改善，因此低于25°以上的區(qū)域的土壤侵蝕模數(shù)。從各坡度帶的土壤侵蝕模數(shù)看，在3 a中各坡度帶的值均為2010年最大，并且，在每1 a中，坡度小于35°時，隨著坡度的增加，土壤侵蝕模數(shù)增大。總體而言，25°以上的坡度帶應(yīng)成為重點預(yù)防和治理的區(qū)域。

3.3.4 土壤侵蝕與土地利用類型 從土地利用類型可以看出，研究區(qū)大部分地區(qū)分布著耕地、林地和草地，將土地利用類型與土壤侵蝕強度進行疊加可以得出如表7所示結(jié)果。

表7 引黃入晉北干線沿線地區(qū)各年份不同土地利用類型土壤侵蝕變化統(tǒng)計

(1) 不同土地利用類型之間，土壤侵蝕狀況差異顯著，3 a中土壤侵蝕面積均為耕地最大，其次為草地、林地、建筑用地、水體、未利用土地。耕地以輕度侵蝕為主，這是因為陡坡的耕地面積大，在黃土高原地區(qū)，耕地一般沿著等高線帶狀耕作[4]，耕作方式比較粗放，缺少水土保持的措施，加之人類活動頻繁干擾，導(dǎo)致耕地侵蝕。

(2) 林地和草地以輕度和中度侵蝕為主，這是因為林地和草地大都分布在海拔1 000 m以上的山區(qū)，林地被人為毀林開荒，草地基本為自然草地，當(dāng)降雨發(fā)生時，這些未采取水保措施的區(qū)域，在降雨和陡坡的作用下，會不同程度的產(chǎn)生土壤侵蝕。

(3) 建筑用地的水土保持良好，而水體基本無土可蝕，未利用土地所占面積雖小，但土壤侵蝕強度高，在山陰縣南部和懷仁縣北部，有部分區(qū)域地表裸露，無植被覆蓋且坡度大，土壤侵蝕嚴重。

由以上分析可知，區(qū)域內(nèi)耕地的侵蝕面積大，林地、草地、未利用土地的侵蝕級別高，因此應(yīng)加強以上土地利用類型的合理利用和管理規(guī)劃，制定相應(yīng)的水土保持措施，改善區(qū)域內(nèi)的環(huán)境狀況。

3.4 土壤侵蝕的人為因素作用

2005—2010年土壤侵蝕模數(shù)增長將近1倍主要是由于工程建設(shè)時期的一系列施工活動會對地表進行開挖、擾動、再塑，還會損毀水土保持設(shè)施，造成大面積的土地裸露，加重水土流失。具體到各因子中主要表現(xiàn)為： ①施工過程中，對地表的擾動再塑，會形成大量的人工坡面和懸空面，在短期內(nèi)坡度、坡長等地形因子會發(fā)生劇烈的變化； ②項目建設(shè)過程中會產(chǎn)生大量的棄土棄渣，這些物質(zhì)穩(wěn)定性差，在無植被覆蓋的情況下，極易產(chǎn)生風(fēng)蝕和水蝕； ③對地表覆蓋物(林地、草地)的清除或砍伐，會使得防風(fēng)固沙的能力下降，極易引起風(fēng)力侵蝕，并且大量的土地直接裸露，土壤直接遭受雨滴擊濺和沖刷，極易引起水力侵蝕。

2010—2015年土壤侵蝕模數(shù)下降主要是由于工程建成后，實行了一系列的工程措施和生物措施來恢復(fù)該地的生態(tài)環(huán)境。如坡面綠化、庫岸治理護岸、道路護坡固土、整地措施、抗旱造林措施、攔截措施等有效的減緩了水土流失[15]。

4 討論與結(jié)論

4.1 討 論

本研究仍存在一定的不足：本文的研究結(jié)果仍需水文站實測的泥沙數(shù)據(jù)進行驗證，各因子的數(shù)據(jù)精度不一，難免在柵格計算過程中產(chǎn)生誤差，例如降雨侵蝕力的計算只采用了4個雨量站點的數(shù)據(jù)，如果采用更多的雨量站點數(shù)據(jù)，可以進一步提高計算精度。P因子的賦值不確定性因素較高，如有定量計算P因子的方法，計算的結(jié)果將更為客觀。未來的研究將根據(jù)研究區(qū)背景特征修正模型中因子的參數(shù)，結(jié)合土壤允許流失量的閾值來對研究區(qū)土壤侵蝕風(fēng)險進行評估。

4.2 結(jié) 論

(1) 2005—2015年，引黃入晉北干線所經(jīng)地區(qū)土壤侵蝕模數(shù)呈先增后減，總體減小的趨勢。2005,2015年屬于輕度侵蝕，而2010年屬于中度侵蝕，大于水利部公布的黃土高原土壤允許流失量1 000 t/(km2·a)的標準。引黃入晉北干線的實施是2010年土壤侵蝕模數(shù)增長近1倍的重要原因，表明人類活動的影響在土壤侵蝕的研究中不可忽視。

(2) 研究區(qū)土壤侵蝕呈塊狀侵蝕和點狀侵蝕相間分布的特點。以北干線為界界線兩側(cè)土壤侵蝕差異顯著。土壤侵蝕的空間分布與研究區(qū)的地貌單元關(guān)系密切，侵蝕中低山區(qū)高，沖洪積平原低。朔州市轄區(qū)的土壤侵蝕面積和土壤侵蝕量最大，朔州市轄區(qū)、山陰、懷仁、大同市轄區(qū)的土壤侵蝕模數(shù)均與總的土壤侵蝕模數(shù)分布特點一致。

(3) 各因子中，R因子與土壤侵蝕模數(shù)呈明顯正相關(guān)，當(dāng)土壤類型為栗褐土、坡度為8°～15°、土地利用類型為耕地時土壤侵蝕面積最大，當(dāng)坡度大于25°、土地利用類型為草地和林地時土壤侵蝕量最高。