黃土高原典型流域土壤侵蝕對退耕還林土地利用變化的響應(yīng)

汪濱，張志強

（1. 北京林業(yè)大學水土保持學院，水土保持與荒漠化防治教育部重點實驗室，北京 100083；2. 山西吉縣森林生態(tài)系統(tǒng)國家野外科學觀測研究站，臨汾 042200）

0 引言

由于土壤抗蝕能力較弱以及人類長期的不合理土地開發(fā)利用，黃土高原屬于世界上土壤侵蝕最嚴重的地區(qū)[1-2]，為了從根本上有效控制土壤侵蝕，中國于1999年起在黃土高原實施了大規(guī)模的退耕還林工程。該工程以改變不合理的土地利用、遏制土壤侵蝕、減少入黃泥沙并兼顧生態(tài)、經(jīng)濟和社會效益為主要目的，以坡耕地退耕還林還草、荒山造林、封山育林和梯田建設(shè)等為主要措施[3-4]，實施后土地利用結(jié)構(gòu)明顯優(yōu)化，取得了顯著的綜合效益尤其是水土保持效益[4-6]。退耕還林工程歷時較長，不同階段實施了不同的主要措施，直接驅(qū)動土地利用持續(xù)發(fā)生變化，深入研究這一變化對土壤侵蝕的影響，對于促進退耕還林成果鞏固和推動退耕還林高質(zhì)量發(fā)展具有重要意義。

土壤侵蝕的發(fā)生是降雨、土地利用、地形、土壤等各種因素相互影響與制約的綜合結(jié)果[7]，其中降雨和土地利用是2個主要影響因素[8-9]。黃土高原有關(guān)試驗與觀測研究表明，降雨量、降雨強度和土地利用都對土壤侵蝕產(chǎn)生了明顯影響[10-12]。運用模型情景模擬（如RUSLE（revised universal soil loss equation）[13]、SWAT（soil and water assessment tool）[14]、MMF（Morgan,Morgan and Finney）[15]）和數(shù)理統(tǒng)計（如雙累積曲線[16]、線性回歸[17]、彈性系數(shù)[18]）等方法，黃土高原已開展了降雨和土地利用對退耕還林期間土壤侵蝕變化影響程度方面的研究，得出的共同結(jié)論是土地利用變化在一定程度上減輕了土壤侵蝕，但在影響程度方面存在著2 種完全不同的結(jié)論，一些研究觀點認為土地利用變化的影響程度相對較大[16,18-20]，另一些研究的結(jié)論則相反[13,17,21]。退耕還林工程的具體實施對象是相關(guān)的土地利用類型（以下簡稱地類），實施的結(jié)果是這些地類中部分地區(qū)（或地塊）發(fā)生了地類的轉(zhuǎn)換（即類別變更），如坡耕地通過退耕還林轉(zhuǎn)至林地、草地通過工程整地轉(zhuǎn)至耕地[14,22]，與此同時這些地類內(nèi)未發(fā)生轉(zhuǎn)換的地區(qū)（或地塊）也在發(fā)生變化，如坡耕地實施了梯改、林草地實施了封禁，且從已有的研究來看未發(fā)生轉(zhuǎn)換的面積較大[14,22-23]。土地利用是人類對土地自然屬性的利用方式和狀況[24]，地類轉(zhuǎn)換地區(qū)所發(fā)生的變化和未轉(zhuǎn)換地區(qū)所發(fā)生的變化通過改變侵蝕環(huán)境，如耕地轉(zhuǎn)草地后植物由作物轉(zhuǎn)變?yōu)椴輩?、林草地的封禁提高了植被覆蓋度，都會對土壤侵蝕產(chǎn)生直接影響，故同時研究這2 種變化對土壤侵蝕的影響尤其是剔除降雨變化后的影響，可完整地揭示土地利用變化對土壤侵蝕的影響過程。目前黃土高原已開展了大量有關(guān)這方面的研究，其中以對地類或區(qū)域土壤侵蝕整體變化影響的研究為主[25-27]，也有一些研究側(cè)重轉(zhuǎn)換對侵蝕的影響[28-30]，但前者沒有區(qū)分這2 種變化，后者沒有關(guān)注未轉(zhuǎn)換地區(qū)的變化，且均疊加了降雨變化的影響，不能直接反映土地利用變化對土壤侵蝕的影響，關(guān)于剔除降雨變化后這2 種變化影響方面的研究少有報道。實施退耕還林工程后，一些地類對土壤侵蝕的影響最大，黃土高原已有的研究表明影響最大的地類主要有旱地和中低覆蓋草地[31]、林地和草地[32]等，有的地區(qū)建設(shè)用地也發(fā)揮了重要作用[33]，這些研究主要依據(jù)降雨變化下地類整體或其轉(zhuǎn)換地區(qū)侵蝕變化的重要性，關(guān)于剔除降雨變化后地類轉(zhuǎn)換地區(qū)和未轉(zhuǎn)換地區(qū)侵蝕變化重要性方面的研究也少有報道。

本研究以山西省吉縣清水河流域為對象，應(yīng)用RUSLE模型分析流域2000—2020 年土壤侵蝕強度的變化特征，通過情景模擬判別土地利用變化和降雨變化對土壤侵蝕的影響程度，將土地利用變化分解為土地利用轉(zhuǎn)換和土地利用改造2 種形式，并在剔除降雨變化影響的基礎(chǔ)上，分析退耕還林工程實施過程中這2 種變化形式對土壤侵蝕的影響及其在影響中的相互關(guān)系和各自作用，識別影響流域土壤侵蝕變化的主要地類，以完整地反映土地利用變化對土壤侵蝕的影響程度和過程，為黃土高原退耕還林成果鞏固及高質(zhì)量發(fā)展有效措施的制定提供科學依據(jù)。

1 研究流域及其退耕還林工程實施概況

清水河發(fā)源于山西省吉縣高天山西麓，經(jīng)州川河向西南匯入黃河，是黃河中游河口鎮(zhèn)至龍門區(qū)間左岸的一級支流。吉縣水文站控制的清水河流域地理位置為110°36′47″E～110°56′0″E，36°2′18″N～36°16′23″N（圖1）。

圖1 清水河流域區(qū)位及水文站和雨量站分布圖Fig.1 Location of the Qingshuihe watershed and distributions of hydrologic and rainfall stations

跨吉縣、鄉(xiāng)寧兩縣，總土地面積436 km2，其中吉縣境內(nèi)占94.67%。為黃土殘塬溝壑區(qū)典型流域，土壤類型以褐土為主。屬暖溫帶半濕潤大陸性季風氣候，多年平均氣溫10.2 ℃，多年平均降水量541.5 mm，年內(nèi)降水不均，根據(jù)吉縣水文站2000—2020 年的觀測結(jié)果，汛期5—10 月降雨量占全年降水量的86.58%。

清水河流域從2000 年起屬于黃土高原退耕還林工程的重點實施地區(qū)，經(jīng)歷了不同的發(fā)展時期：2000 年下半年開始實施，2005 年完成退耕和造林任務(wù)；從2008 年起進入退耕還林成果鞏固期，流域所在的吉縣在繼續(xù)實施荒山造林和封山育林的同時，大力開展坡改梯與基本農(nóng)田建設(shè)；2010 年起吉縣發(fā)揮地方優(yōu)勢，將發(fā)展蘋果產(chǎn)業(yè)作為成果鞏固的重要舉措；2014 年下半年起開始實施新一輪工程。為了完整地反映流域退耕還林工程的實施過程，本研究以2000、2005、2010、2014 和2020 年5個標志年份為節(jié)點，將流域已進行的退耕還林實施時期（2000—2020 年）分為2000—2005、2005—2010、2010—2014 和2014—2020 年4個時段，并依據(jù)工程實施期間主要措施的不同，將其對應(yīng)的實施階段分別分為第1 階段（2000—2005 年）、第2 階段（2006—2009 年）、第3階段（2010—2014 年）和第4 階段（2015—2020 年）。

2 數(shù)據(jù)與方法

2.1 數(shù)據(jù)來源與處理

2.1.1 遙感影像與土地利用分類

本研究選取2000 年5 月21 日Landsat 7、2005 年6月12 日Landsat 5、2011 年6 月29 日Landsat 5、2014 年5 月20 日Landsat 8 和2020 年6 月5 日Landsat 8 影像分別作為5個標志年份的代表影像（因2009 年G22 臨吉高速公路的開工建設(shè)對2010 年沿線地區(qū)土地利用的影響較大，故用2011 年影像替代，相應(yīng)的時段分別調(diào)整為2005—2011 和2011—2014 年，對應(yīng)的實施階段分別調(diào)整為第2 階段（2006—2010 年）和第3 階段（2011—2014 年）），這些影像均源自中國科學院計算機網(wǎng)絡(luò)信息中心地理空間數(shù)據(jù)云（http://www.gscloud.cn/）。以中國土地利用現(xiàn)狀遙感監(jiān)測數(shù)據(jù)庫中的土地利用分類系統(tǒng)為基礎(chǔ)，并參照土地利用現(xiàn)狀分類國家標準（GB/T 21010—2017），將土地利用類型劃分為耕地、園地、有林地、疏林地、灌木林地、草地、建設(shè)用地和水域8 類。在ERDAS 9.2 中經(jīng)監(jiān)督分類并借助Google Earth 和91衛(wèi)圖的高空間分辨率影像，進行人工目視判讀與修正[34]后得到5 期土地利用空間分布圖。利用野外實地采點數(shù)據(jù)、地塊利用歷史實地訪問調(diào)查資料等對5 期分類結(jié)果進行精度檢驗，總體解譯精度均達到90%以上。

2.1.2 降雨數(shù)據(jù)與侵蝕性降雨標準

流域土壤侵蝕主要發(fā)生在汛期，流域內(nèi)及周邊共10個站點的汛期逐日降雨數(shù)據(jù)來源于《中華人民共和國水文年鑒黃河流域水文資料》。依據(jù)謝云等[35]基于黃土高原試驗小區(qū)實測數(shù)據(jù)的研究結(jié)果，將12 mm/d 作為侵蝕性降雨標準。

2.1.3 土壤類型的分布和屬性數(shù)據(jù)

流域土壤類型的分布數(shù)據(jù)來源于《山西省土壤圖集》，屬性數(shù)據(jù)來源于《山西土種志》。

2.1.4 DEM 和NDVI 數(shù)據(jù)

DEM 數(shù)據(jù)為空間分辨率30 m 的ASTER GDEM 影像，源自中國科學院計算機網(wǎng)絡(luò)信息中心地理空間數(shù)據(jù)云（http://www.gscloud.cn/）。NDVI 數(shù)據(jù)中2000、2005、2011 和2014 年數(shù)據(jù)為5—10 月共24 期分辨率500 m 的MODIS NDVI 月合成影像，同樣源自地理空間數(shù)據(jù)云（http://www.gscloud.cn/）。因地理空間數(shù)據(jù)云無2020 年分辨率500 m 的月合成影像，2020 年數(shù)據(jù)選取美國航空航天局（NASA）地球觀測數(shù)據(jù)與信息系統(tǒng)（EOSDIS）（https://search.earthdata.nasa.gov/）4 月22 日至10 月30日共12 期分辨率500 m 的MODIS NDVI 16 日合成影像。

2.2 研究方法

2.2.1 RUSLE 模型的應(yīng)用

黃土高原溝蝕問題突出[36]，其中淺溝侵蝕分布較廣[37]，模擬試驗結(jié)果表明淺溝侵蝕量占坡面侵蝕總量的26.6%～59.7%[38]。流域?qū)儆跉堒珳羡謪^(qū)，溝蝕嚴重，參考黃土高原地區(qū)的有關(guān)研究[39-41]，其侵蝕強度的估算采用增加了淺溝侵蝕影響因子G的改進RUSLE 模型[42]：

式中A為年平均土壤侵蝕模數(shù)，t/(hm2·a)；R為降雨侵蝕力因子，MJ·mm/(hm2·h·a）；K為土壤可蝕性因子，t·hm2·h/(hm2·MJ·mm）；L和S分別為坡長和坡度因子，無量綱；G為淺溝侵蝕影響因子，無量綱；C為覆蓋管理因子，無量綱；P為水土保持措施因子，無量綱。本研究根據(jù)數(shù)據(jù)的可獲取情況，選用了各因子的適當算法。

1）降雨侵蝕力因子R。降雨侵蝕力反映了降雨引起土壤侵蝕的潛在能力[43]，采用章文波等[44]基于日降雨資料提出的算法：

式中Ri為第i個半月的降雨侵蝕力值，MJ·mm/(hm2·h）；Pj為該半月內(nèi)第j日的侵蝕性降雨量，mm；k為該半月內(nèi)的侵蝕性降雨日數(shù)，d；α和β為參數(shù)；Pd和Py分別為研究時段內(nèi)侵蝕性日降雨量和侵蝕性年降雨量的多年平均值，mm。采用泰森多邊形法為各雨量站所控制的區(qū)域賦值。

2）土壤可蝕性因子K。土壤可蝕性反映了土壤性質(zhì)對侵蝕的敏感程度[45]，采用EPIC 模型[46]中的算法以及張科利等[47]基于中國土壤屬性提出的因子值修正算法：

式中Sa為砂粒質(zhì)量分數(shù)，%；Si為粉粒質(zhì)量分數(shù)，%；Cl為黏粒質(zhì)量分數(shù)，%；Co為有機碳質(zhì)量分數(shù)，%；Sn=1—Sa/100；K為基于中國土壤屬性對KEPIC的修正值。

3）坡長因子L和坡度因子S。坡長和坡度是誘發(fā)土壤侵蝕的重要地形因素，L因子和10°以下坡度的S因子直接應(yīng)用RUSLE 用戶手冊[48]中的算法，10°以上的S因子采用LIU等[49]基于實測數(shù)據(jù)提出的算法：

式中λ為水平坡長，m；m為坡長指數(shù)；γ為細溝侵蝕和細溝間侵蝕的比值；θ為坡度，（°）。

4）淺溝侵蝕影響因子G。在無植被覆蓋的黃土陡坡條件下，淺溝侵蝕發(fā)生的臨界坡度為15°，采用江忠善等[42]基于實測數(shù)據(jù)提出的算法：

5）覆蓋管理因子C。植被是土壤侵蝕的重要抑制因素，采用蔡崇法等[50]基于實測數(shù)據(jù)提出的算法：

式中f為植被覆蓋度，%。通過采用像元二分模型[51]對NDVI 數(shù)據(jù)進行反演獲取。

6）水土保持措施因子P。水土保持措施也是土壤侵蝕的重要抑制因素。采用謝怡凡等[52]關(guān)于耕地、建設(shè)用地和水域的賦值方法，耕地（不包括梯田）按坡度分級賦值，坡度≤5°、＞5°～10°、＞10°～15°、＞15°～20°、＞20°～25°和＞25°分別賦值0.100、0.221、0.305、0.575、0.705 和0.800，建設(shè)用地和水域按不發(fā)生侵蝕賦值0。采用秦偉等[53]關(guān)于梯田、園地、林地和草地的賦值方法，梯田賦值0.084，園地參照魚鱗坑造林地一并賦值0.187，未實施水土保持措施的林地和草地賦值1。

2.2.2 土地利用變化形式的分解及數(shù)據(jù)提取

在人類活動驅(qū)動下土地利用與土地覆被變化關(guān)系的分析中，TURNER Ⅱ等[54]分解并定義了不同土地利用下土地覆被的2 種變化形式：一是土地覆被的轉(zhuǎn)換，即一種覆被類型完全轉(zhuǎn)至另一種覆被類型，如森林皆伐改種牧草后轉(zhuǎn)至草地，二是土地覆被的改造，即土地覆被屬性的改變，如草地過牧后成為退化草地，并且無論是轉(zhuǎn)換還是改造，均需對土地覆被進行維護，如梯田維修、森林撫育等。本研究應(yīng)用這一土地覆被變化形式分解原理，將土地利用變化分解為土地利用轉(zhuǎn)換和土地利用改造，兩者分別指同一空間位置上土地利用類別發(fā)生變更和未發(fā)生變更的變化。轉(zhuǎn)換區(qū)和改造區(qū)的面積數(shù)據(jù)提取基于土地利用轉(zhuǎn)移矩陣。

2.2.3 土地利用和降雨變化對土壤侵蝕影響程度的判別

為了判別土地利用變化和降雨變化各自對土壤侵蝕的影響程度，本研究采用RUSLE 模型情景模擬方法，提出了依據(jù)侵蝕模數(shù)升降變化量差值的算法：首先以土地利用與降雨共同變化實際情景下研究時段初的侵蝕模數(shù)為參照，分別計算研究時段末較時段初土地利用與降雨共同變化實際情景下侵蝕模數(shù)的升降變化量（ΔA）、僅土地利用變化模擬情景下侵蝕模數(shù)的升降變化量（ΔAl）和僅降雨變化模擬情景下侵蝕模數(shù)的升降變化量（ΔAr），然后計算ΔAl和ΔAr分別與ΔA的差值（D），最后將這一差值在總差值（均取絕對值）中的比重作為影響程度（E，%）的判別依據(jù)。僅土地利用變化和僅降雨變化對土壤侵蝕影響程度（分別為El和Er）的計算方法如下：

式中Dl和Dr分別為僅土地利用變化模擬情景和僅降雨變化模擬情景下侵蝕模數(shù)升降變化量的差值（t/(hm2·a)），Ab和Ae分別為研究時段初和時段末土地利用與降雨共同變化實際情景下的侵蝕模數(shù)，Al和Ar分別為研究時段末僅土地利用變化模擬情景和僅降雨變化模擬情景下的侵蝕模數(shù)，t/(hm2·a)。

在RUSLE 模型因子中，K、L、S和G因子短時期內(nèi)相對穩(wěn)定，R和C因子逐年發(fā)生變化[17]，因水土保持措施貫穿于整個退耕還林工程實施過程之中，P因子也會在短時期內(nèi)發(fā)生變化，參考黃土高原地區(qū)的有關(guān)研究，將C和P作為反映土地利用變化的因子[19-20]。Ab、Ae、Al和Ar的計算如下：

以上各因子的計算方法與式（1）相同。僅土地利用變化模擬情景下各地類及流域轉(zhuǎn)換區(qū)和改造區(qū)時段初、末的平均侵蝕模數(shù)均分別依據(jù)式（20）和式（22）計算。

2.2.4 影響土壤侵蝕變化的主要土地利用類型識別

為了兼顧地類侵蝕的強度和面積，避免地類侵蝕強度變化大而面積變化小或侵蝕強度變化小而面積變化大2 種極端情況并剔除降雨變化的影響，本研究依據(jù)僅土地利用變化情景的模擬結(jié)果，將侵蝕量變化的大小作為主要影響地類的識別指標。由于各地類轉(zhuǎn)換區(qū)與改造區(qū)侵蝕量的變化可能有增有減，且這樣的增和減均是地類變化影響土壤侵蝕的結(jié)果，本研究提出了依據(jù)侵蝕量增減變化絕對值來識別主要影響地類的算法，計算式如下：

式中Wi為地類i侵蝕量合計變化量的貢獻率，%；Mi為地類i侵蝕量的合計變化量，萬t；Mci和Mmi分別為地類i轉(zhuǎn)換區(qū)和改造區(qū)侵蝕量的增減變化量，萬t。

3 結(jié)果與分析

3.1 土壤侵蝕變化基本特征

根據(jù)水力侵蝕強度分級國家行業(yè)標準[55]，將流域土壤侵蝕強度劃分為微度、輕度、中度、強烈、極強烈和劇烈6個級別。在ArcGIS 9.3 中計算RUSLE 模型，得到流域2000、2005、2011、2014 和2020 年的土壤侵蝕強度圖，對侵蝕強度進行分級，得到土壤侵蝕強度分級空間分布圖（圖2），疊加分析土壤侵蝕強度圖與土地利用圖，得到流域各地類的平均土壤侵蝕模數(shù)（表1）。

表1 清水河流域2000—2020 年土地利用類型的平均土壤侵蝕模數(shù)Table 1 Average soil erosion moduli of land use types in the Qingshuihe watershed from 2000 to 2020

圖2 清水河流域 2000—2020 年土壤侵蝕強度分級空間分布Fig.2 Spatial distribution of soil erosion intensity grades in the Qingshuihe watershed from 2000 to 2020

由圖2 可知，微度侵蝕主要分布在高天山和人祖山山區(qū)、東部和中部大部殘塬地區(qū)以及清水河河谷地帶，與2000 年相比較，2005、2011、2014 和2020 年其分布范圍依次明顯擴大；輕度、中度和強烈侵蝕主要分布在山地周邊、部分殘塬和河谷兩側(cè)地區(qū)，2005、2011 和2014 年其分布范圍穩(wěn)定縮小，2020 年明顯縮??；極強烈和劇烈侵蝕主要分布在溝壑地區(qū)，2005、2011、2014 和2020 年其分布范圍總體上呈依次縮小趨勢。由表1 可知，2000 年流域平均侵蝕模數(shù)36.21 t/(hm2·a)，侵蝕強度屬于中度級別，2005 年升至41.02 t/(hm2·a)，接近強烈級別，2011 年大幅降至24.93 t/(hm2·a)，進入輕度級別，2014 年繼續(xù)降至23.72 t/(hm2·a)，仍保持輕度級別，2020年大幅降至8.24t/(hm2·a)，進入微度級別。

由表1 還可以看出：2000 年草地平均侵蝕模數(shù)高達53.04 t/(hm2·a)，屬于強烈級別，而園地只有5.43 t/(hm2·a)，屬于微度級別，表明除建設(shè)用地和水域外，退耕還林初各地類之間的平均侵蝕模數(shù)差異較大；與前一年份相比較，2005 年除草地上升外，其他地類均為下降，而2014 年除有林地下降外，其他地類均為上升，2011 和2020 年各地類則均為下降，2005、2011 和2014 年各地類平均侵蝕模數(shù)自高至低的排序依次為草地、疏林地、灌木林地、耕地、有林地和園地，2020 年的排序則依次為草地、疏林地、耕地、灌木林地、園地和有林地，表明通過不同階段退耕還林措施的實施，各地類的平均侵蝕模數(shù)總體處于下降趨勢，但地類之間仍存在明顯差異。

3.2 土地利用變化對土壤侵蝕的影響程度

在土地利用變化和降雨變化對土壤侵蝕影響程度的情景模擬過程中，將侵蝕性降雨量作為降雨變化的計算依據(jù)。流域1995—2020 年侵蝕性降雨量的變化見圖3，其中退耕還林實施前（1995—1999 年）與4個階段（分別為 2000—2005、2006—2010、2011—2014 和2015—2020 年）的年均侵蝕性降雨量分別為269.0、333.6、317.9、434.6 和331.3 mm。TFPW—MK 趨勢性檢驗結(jié)果為：顯著性檢驗統(tǒng)計量Z=1.76＞1.64，表明在顯著性水平0.10 下1995—2020 年的侵蝕性降雨量呈顯著增加趨勢（P=0.078＜0.10，年均變化率3.7 mm）。

圖3 清水河流域1995—2020 年侵蝕性降雨量的變化Fig.3 Change of annual erosive rainfall in the Qingshuihe watershed from 1995 to 2020

依據(jù)式（13）～（23）模擬了流域土地利用和降雨變化各自對土壤侵蝕的影響程度，結(jié)果見表2。

表2 清水河流域土地利用和降雨變化對土壤侵蝕影響程度的模擬結(jié)果Table 2 Simulations of effect degrees of land use and rainfall change on soil erosion in the Qingshuihe watershed

由表2 可知第1 階段、第3 階段末較階段初僅土地利用變化情景下的流域平均侵蝕模數(shù)分別下降5.54 和7.34 t/(hm2·a)，僅降雨變化情景下的流域平均侵蝕模數(shù)分別上升14.94 和8.67 t/(hm2·a)，第1 階段僅土地利用變化和僅降雨變化對土壤侵蝕的影響程度分別為49.46%和50.54%，差距較小，表明土地利用變化和降雨變化在流域土壤侵蝕的變化中起了同等重要的作用，第3 階段的影響程度分別為61.71%和38.29%，差距明顯，表明土地利用變化起了主導作用。

第2 階段、第4 階段末較階段初僅土地利用變化情景和僅降雨變化情景下的流域平均侵蝕模數(shù)均為下降，第2 階段分別下降11.90 和6.04 t/(hm2·a)，第4 階段分別下降13.92 和12.95 t/(hm2·a)，這2個階段僅土地利用變化的影響程度分別為70.58%和61.86%，明顯大于僅降雨變化的影響程度（分別為29.42%和38.14%），表明土地利用變化在流域土壤侵蝕的變化中都起了主導作用。

綜上可知，僅土地利用變化情景下的流域侵蝕模數(shù)階段平均下降9.67 t/(hm2·a)，僅降雨變化情景下侵蝕模數(shù)階段平均上升1.15 t/(hm2·a)，兩者對土壤侵蝕的影響程度分別為75.23%和24.77%，前者明顯大于后者，表明土地利用變化在流域土壤侵蝕的變化中總體上起了主導作用。

3.3 土地利用變化對土壤侵蝕的影響過程

3.3.1 土地利用轉(zhuǎn)換和改造過程對土壤侵蝕的影響

為了分析流域退耕還林不同階段土地利用轉(zhuǎn)換和改造過程對土壤侵蝕的影響，依據(jù)了土地利用轉(zhuǎn)移矩陣和僅土地利用變化情景下階段初和階段末各地類轉(zhuǎn)換區(qū)與改造區(qū)平均侵蝕模數(shù)的模擬結(jié)果。第1 階段退耕還林工程實施的主要措施為坡耕地退耕、荒山造林（包括補植，配套了魚鱗坑整地等）和封山育林，未退耕地和園地加強了相關(guān)水土保持措施，如坡耕地梯改和水保耕作、園地整治等，并且這些措施實施后草地、灌木林地、疏林地和有林地內(nèi)部及之間植被的恢復與更新演替開始見效，由此驅(qū)動的土地利用轉(zhuǎn)換和改造過程影響了各地類轉(zhuǎn)換區(qū)和改造區(qū)土壤侵蝕強度的變化。耕地轉(zhuǎn)草地、草地轉(zhuǎn)3 類林地是此階段的主要轉(zhuǎn)換過程，地類類別的變更（改變了轉(zhuǎn)換區(qū)的土壤侵蝕環(huán)境）及所實施主要相關(guān)措施的差異直接影響流域土地利用轉(zhuǎn)換區(qū)侵蝕強度的變化。以坡耕地退耕為例，地類類別從耕地主要轉(zhuǎn)至草地，所實施的主要相關(guān)措施由水保耕作轉(zhuǎn)為草地植被封禁，由于地類的侵蝕強度本身已經(jīng)蘊涵了地類類別和實施措施的影響，因而這一過程主要屬于侵蝕強度較低地類向侵蝕強度較高地類轉(zhuǎn)換的過程，耕地轉(zhuǎn)換區(qū)的平均侵蝕模數(shù)由17.03 t/(hm2·a)升至21.48 t/(hm2·a)（圖4a）。耕地、灌木林地和草地的改造是此階段的主要改造過程，所實施的主要相關(guān)措施直接影響流域土地利用改造區(qū)侵蝕強度的變化。以灌木林地為例，補植和封禁實施后改造區(qū)植被的恢復與更新演替開始見效，平均侵蝕模數(shù)由16.58 t/(hm2·a)降至6.90 t/(hm2·a)（圖4b）。

圖4 清水河流域僅土地利用變化情景下土地利用類型轉(zhuǎn)換區(qū)和改造區(qū)平均土壤侵蝕模數(shù)的升降變化Fig.4 Variations of average soil erosion moduli in conversion regions and modification regions of land use types under scenario of only land use change in the Qingshuihe watershed

第2 階段為退耕還林成果鞏固期，在繼續(xù)實施第1階段造林和封禁措施的同時，地方政府開展了新果園建設(shè)和荒坡墾殖，以提高退耕后農(nóng)民的生活水平并解決糧食問題，期間草地和3 類林地內(nèi)部及之間植被的恢復與更新演替效果明顯。耕地轉(zhuǎn)園地、疏林地和灌木林地轉(zhuǎn)有林地、草地轉(zhuǎn)耕地和3 類林地是此階段的主要轉(zhuǎn)換過程，以疏林地為例，在主要轉(zhuǎn)至有林地的過程中，主要實施措施均為補植和封禁，但地類類別發(fā)生變更，轉(zhuǎn)換過程主要屬于侵蝕強度較高地類向侵蝕強度較低地類的轉(zhuǎn)換，轉(zhuǎn)換區(qū)的平均侵蝕模數(shù)由階段初的27.11 t/(hm2·a)降至階段末的20.41 t/(hm2·a)（圖4c）。有林地、灌木林地和草地的改造是此階段的主要改造過程，以有林地為例，補植和封禁后植被的恢復與更新演替明顯見效，改造區(qū)的平均侵蝕模數(shù)由8.51 t/(hm2·a)降至6.61 t/(hm2·a)（圖4d）。

第3 階段在延續(xù)第2 階段主要措施的同時，坡耕地梯改持續(xù)推進，地方政府為了發(fā)展蘋果產(chǎn)業(yè)，進一步加大了果園整治和新果園建設(shè)力度，且期間草地和3 類林地內(nèi)部及之間植被的恢復與更新演替已開始顯著見效。此階段的主要轉(zhuǎn)換過程為耕地轉(zhuǎn)園地、灌木林地轉(zhuǎn)有林地、草地轉(zhuǎn)園地和3 類林地，以灌木林地為例，主要轉(zhuǎn)至有林地，同上階段的疏林地一樣也屬于侵蝕強度較高地類向侵蝕強度較低地類的轉(zhuǎn)換，轉(zhuǎn)換區(qū)的平均侵蝕模數(shù)由階段初的10.95 t/(hm2·a)降至階段末的6.33 t/(hm2·a)（圖4e）。此階段的主要改造過程為園地、有林地、灌木林地和草地的改造，以草地為例，通過嚴格封禁，植被的恢復與更新演替使得其改造區(qū)平均侵蝕模數(shù)由55.37 t/(hm2·a)降至50.25 t/(hm2·a)（圖4f）。

第4 階段為新一輪退耕還林工程實施期，地方政府圍繞生態(tài)林業(yè)、民生林業(yè)兩大主體，突出了第1 輪退耕還林成果鞏固、植被恢復和生態(tài)果園建設(shè)等工作，由此果園整治和新果園建設(shè)穩(wěn)步推進，草地和3 類林地內(nèi)部及之間植被的恢復與更新演替已顯著見效。此階段的轉(zhuǎn)換過程主要為疏林地轉(zhuǎn)有林地、灌木林地轉(zhuǎn)有林地和疏林地、草地轉(zhuǎn)園地和3 類林地，以草地為例，因其自身的侵蝕強度在各地類中最高，故在地類類別變更及所實施主要措施（草地實施的措施主要為封禁）差異的影響下，轉(zhuǎn)至任何其他地類都屬于侵蝕強度較高地類向侵蝕強度較低地類的轉(zhuǎn)換，轉(zhuǎn)換區(qū)的平均侵蝕模數(shù)由階段初的42.11 t/(hm2·a)降至階段末的12.61 t/(hm2·a)（圖4g）。此階段的主要改造過程仍為園地、有林地、灌木林地和草地的改造，以園地為例，整治措施的實施直接降低了侵蝕強度，改造區(qū)的平均侵蝕模數(shù)由4.83 t/(hm2·a)降至3.23 t/(hm2·a)（圖4h）。

3.3.2 土地利用轉(zhuǎn)換和改造在影響土壤侵蝕中的相互關(guān)系

由于各地類內(nèi)的侵蝕強度在空間上并非均勻分布，故轉(zhuǎn)換區(qū)和改造區(qū)的起始侵蝕強度存在差異：如果某地類內(nèi)侵蝕強度相對較高的部分地區(qū)（或地塊）轉(zhuǎn)至其他地類，其轉(zhuǎn)換區(qū)侵蝕強度變化的起點就相應(yīng)較高（反之較低），此時剩余地區(qū)（或地塊）的侵蝕強度也就相對較低，改造區(qū)侵蝕強度變化的起點也就相應(yīng)較低（反之較高），兩者之間存在著此消彼長的關(guān)系。

為了分析這一關(guān)系，對各階段初僅土地利用變化情景下流域土地利用轉(zhuǎn)換區(qū)和改造區(qū)的平均侵蝕模數(shù)（表3）進行了比較：流域4個階段轉(zhuǎn)換區(qū)的起始侵蝕模數(shù)分別為27.07、24.02、21.84 和23.01 t/(hm2·a)，分別較改造區(qū)的起始侵蝕模數(shù)（37.59、45.21、26.35 和23.90 t/(hm2·a)）低27.99%、46.87%、17.12%和3.72%，改造區(qū)的階段平均起始侵蝕模數(shù)為33.63 t/(hm2·a)，較轉(zhuǎn)換區(qū)的階段平均起始侵蝕模數(shù)（23.44 t/(hm2·a)）高43.47%。這一結(jié)果表明，各階段土地利用轉(zhuǎn)換主要是從侵蝕強度相對較輕的地區(qū)開始的，而土地利用改造則正好相反，改造區(qū)水土流失綜合治理的難度總體上大于轉(zhuǎn)換區(qū)。

表3 清水河流域僅土地利用變化情景下土地利用轉(zhuǎn)換區(qū)和改造區(qū)各時段初與末的平均土壤侵蝕模數(shù)Table 3 Average soil erosion moduli in land use conversion regions and modification regions under scenario of only land use change at the beginning and the end of each period in the Qingshuihe watershed

3.3.3 土地利用轉(zhuǎn)換和改造在影響土壤侵蝕中的各自作用

土地利用轉(zhuǎn)換是在不同地類之間進行的，而改造是在同一地類內(nèi)進行的，在起始侵蝕強度差異的基礎(chǔ)上，轉(zhuǎn)換和改造過程的不同直接導致轉(zhuǎn)換區(qū)和改造區(qū)之間侵蝕強度的升降變化量出現(xiàn)差異。為了分析這一差異，利用表3 對僅土地利用變化情景下流域土地利用轉(zhuǎn)換區(qū)和改造區(qū)各階段末較階段初平均侵蝕模數(shù)的升降變化量進行了比較：流域4個階段轉(zhuǎn)換區(qū)侵蝕模數(shù)的下降量分別為10.11、12.05、14.15 和14.73 t/(hm2·a)，分別較改造區(qū)的下降量（4.86、11.86、4.21 和13.71 t/(hm2·a)）高108.02%、1.60%、236.10%和7.44%，轉(zhuǎn)換區(qū)侵蝕模數(shù)的階段平均下降量為13.18 t/(hm2·a)，較改造區(qū)的階段平均下降量（8.74 t/(hm2·a)）高50.80%。這一結(jié)果表明土地利用轉(zhuǎn)換擁有較強的侵蝕模數(shù)調(diào)節(jié)能力，在降低其實施地區(qū)的侵蝕強度方面發(fā)揮了重要作用。

土地利用轉(zhuǎn)換和改造在土壤侵蝕變化中的作用還與各自面積變化的大小有關(guān)。從表3 可以看出，流域4個階段轉(zhuǎn)換面積均明顯小于改造面積，這樣的面積差異與上述各階段末較階段初侵蝕強度升降變化量差異的結(jié)合，直接導致這2 種變化區(qū)之間侵蝕量的增減變化量也隨之出現(xiàn)差異。為了分析這一差異，利用表3進行了侵蝕量增減變化量的計算，結(jié)果顯示：流域4個階段轉(zhuǎn)換區(qū)侵蝕量的減少量分別為5.76、10.39、19.40 和13.12 萬t，改造區(qū)侵蝕量的減少量分別為18.42、41.49、12.58 和47.57 萬t，除第3 階段轉(zhuǎn)換區(qū)是改造區(qū)的1.54 倍外，其他3個階段改造區(qū)分別是轉(zhuǎn)換區(qū)的3.20、3.99 和3.63 倍，4個階段改造區(qū)侵蝕量的合計減少量120.06 萬t，是轉(zhuǎn)換區(qū)階段合計減少量（48.67 萬t）的2.47 倍，占流域階段合計減少量（168.73 萬t）的71.16%。這一結(jié)果表明，由于實施面積較大，土地利用改造擁有較強的侵蝕量調(diào)節(jié)能力，在減少流域侵蝕總量方面發(fā)揮了重要作用。

3.3.4 對流域土壤侵蝕變化影響最大的土地利用類型

依據(jù)僅土地利用變化情景下各地類轉(zhuǎn)換區(qū)和改造區(qū)侵蝕量增減變化的模擬結(jié)果，利用式（24）～（25）計算各地類侵蝕量合計變化量的貢獻率，結(jié)果見表4。

表4 清水河流域僅土地利用變化情景下各土地利用類型土壤侵蝕量合計變化量的貢獻率Table 4 Contribution rates of total variation of soil erosion amount of land use types under scenario of only land use change in the Qingshuihe watershed

由表4 可知：草地侵蝕量合計變化量（12.12、47.68、22.68 和 39.01萬t）及貢獻率（46.37%、91.31%、69.32%和63.71%）在4個階段均位居各地類之首，4個階段合計變化量為121.49 萬t，對流域侵蝕量階段合計變化量（172.31 萬t）的貢獻率為70.51%，在4個階段草地合計變化量中，改造區(qū)的階段合計變化量（81.90萬t）占67.41%。這一結(jié)果表明在各地類中草地對流域土壤侵蝕變化的影響最大，并主要源自其改造的影響。

4 結(jié)論與建議

2000、2005、2011、2014 和2020 年清水河流域平均侵蝕模數(shù)分別為36.21、41.02、24.93、23.72 和8.24 t/(hm2·a)，侵蝕強度從中度降至微度級別，地類之間侵蝕強度差異明顯，土地利用和降雨變化對土壤侵蝕的階段平均影響程度分別為75.23%和24.77%，退耕還林土地利用變化是土壤侵蝕變化的主導影響因素。通過將土地利用變化分解為土地利用轉(zhuǎn)換和改造2 種形式并通過情景模擬剔除降雨變化的影響，得到如下結(jié)論：

1）流域土地利用轉(zhuǎn)換區(qū)侵蝕強度的變化直接受轉(zhuǎn)換過程中地類類別變更及所實施主要相關(guān)措施差異的影響，而改造區(qū)侵蝕強度的變化直接受改造過程中所實施主要相關(guān)措施的影響。

2）流域土地利用改造區(qū)的階段平均起始侵蝕模數(shù)為33.63 t/(hm2·a)，較轉(zhuǎn)換區(qū)的階段平均起始侵蝕模數(shù)（23.44 t/(hm2·a)）高43.47%，其水土流失綜合治理的難度總體上大于轉(zhuǎn)換區(qū)。

3）流域土地利用轉(zhuǎn)換區(qū)的侵蝕模數(shù)階段平均下降13.18 t/(hm2·a)，較改造區(qū)的階段平均下降量（8.74 t/(hm2·a)）高50.80%，土地利用轉(zhuǎn)換在降低其實施地區(qū)侵蝕強度方面發(fā)揮了重要作用，而改造區(qū)的侵蝕量階段合計減少120.06 萬t，占流域階段合計減少量（168.73 萬t）的71.16%，土地利用改造因?qū)嵤┟娣e較大在減少流域侵蝕總量方面發(fā)揮了重要作用。

4）草地轉(zhuǎn)換區(qū)和改造區(qū)侵蝕量的階段合計變化量為121.49 萬t，占流域階段合計變化量（172.31 萬t）的70.51%，在草地階段合計變化量中，改造區(qū)的階段合計變化量（81.90 萬t）占67.41%，草地是對流域土壤侵蝕變化影響最大的地類，并主要源自其改造的影響。

流域總體侵蝕強度已降至微度級別，但草地和部分林地的侵蝕仍較嚴重，建議后續(xù)退耕還林工程的實施及成果維護應(yīng)將提高林草生態(tài)系統(tǒng)的質(zhì)量作為重點，尤其應(yīng)注重草地改造措施的實施，以期獲得更大的綜合效益。本研究中土壤侵蝕強度的估算主要基于RUSLE 模型，降雨變化的情景模擬主要針對降雨侵蝕力R因子，沒有剔除降雨變化對覆蓋管理C因子的影響，建議深入開展這方面的相關(guān)研究。

致謝：非常感謝北卡羅來納大學教堂山分校宋從和教授與自然資源部國土衛(wèi)星遙感應(yīng)用中心李憲文研究員的指導以及吉縣林業(yè)局在實地考察中所提供的幫助。