青藏高原農耕區土地利用變化對土壤侵蝕的影響

周 檉，邵東國，顧文權，姚明磊

（武漢大學水資源與水電工程科學國家重點實驗室，武漢 430072）

0 引 言

受地形、凍融、風力、水力、重力等多種侵蝕營力影響，青藏高原土壤侵蝕嚴重［1］，致使土層變薄、生態失衡、環境退化［1-3］。近年來，青藏高原大面積的退耕或開墾等土地利用變化在改變局部微地形的同時，會對土壤侵蝕動力造成影響，減弱或誘發土壤侵蝕［4］。因此，探究土地利用變化對土壤侵蝕的影響及其對土壤減蝕的貢獻，不僅對揭示該地區土地利用對土壤侵蝕的擾動效應具有重要的科學意義，而且對合理開發利用水土資源、防治土壤侵蝕具有實用前景。

青藏高原作為我國重要的生態屏障，其生態環境亟待保護，備受中國政府與學者的關注。國務院于2011年印發了《關于印發青藏高原區域生態建設與環境保護規劃（2011-2030年）的通知》，以加強青藏高原的生態建設與環境保護。多年來，我國學者就青藏高原的土壤侵蝕特征［5，6］、強度［7，8］及其影響因素［9，10］開展了大量研究。胡玲等［11］采用野外沖刷試驗研究了青藏高原公路邊坡在不同植被防護措施下的產流產沙特征，揭示了初始產流時間、流量和植被蓋度的關系，邊坡產流量、產沙量與植被蓋度的關系，以及不同植被措施的邊坡防護效果；康琳琦等［12］基于RUSLE 模型分析了青藏高原的土壤侵蝕強度時空變化特征，發現不同時期內土壤侵蝕變化不大，侵蝕強度由南向北逐漸減弱；以及許多學者在局部地區開展了地質地貌［13］、降水［14］、植被［15］與土壤侵蝕的關系研究，發現拉薩河流域坡度為15°～25°的區域土壤侵蝕最為嚴重；三江源地區在降雨次數和降雨量增多的情況下，降雨侵蝕力在1961-2012年間不斷增大；三江源地區的植被覆蓋度隨物種多樣性增大而減小，土壤侵蝕隨之增大。現有青藏高原土地利用的研究成果，主要聚焦于估算單一年份不同土地利用類型的土壤侵蝕強度［13］，從空間維度分析不同土地利用的侵蝕特征，以此指導土地利用規劃，而時空維度上的分析欠缺量化，即土地利用變化對土壤侵蝕的驅動作用，在青藏高原地區尚缺乏量化研究，需要深入研究同一地塊上年際間土地利用變化對土壤侵蝕的驅動作用機制，揭示土地利用變化對當下土壤侵蝕的影響與減蝕貢獻度，為未來土地利用發展規劃提供科學依據。

本文以青藏高原農耕區為研究對象，基于CSLE 模型，解析了農耕區2015年土壤侵蝕狀態，探討了土壤侵蝕強度變化綜合指數及區域減蝕貢獻度，量化了不同耕地變化模式對自身土壤侵蝕的影響及其對區域土壤減蝕的貢獻，以期為青藏高原各個地區的耕地開發利用和保護、土壤侵蝕防治等提供參考。

1 材料與方法

1.1 研究區概況

青藏高原位于中國西南部地區，涉及六個省級行政區（圖一），南鄰喜馬拉雅山脈南緣，北鄰昆侖山、阿爾金山和祁連山北緣，西鄰帕米爾高原和喀喇昆侖山脈，東鄰秦嶺山脈西段和黃土高原，介于北緯26°00'～39°47'，東經73°19'～104°47'之間，面積達257.24 萬km2。降水季節性顯著，主要分布在5-9月份，區域性顯著，自東南向西北逐漸減少；地域遼闊，水土資源豐富，其中水資源總量占全國22.71%，土地資源主要用于牧業，其次用于林業，耕地稀少，僅占0.9%，主要分布在一江兩河流域與湟水流域，其他區域分布較為破碎，集約化程度較低。

圖1 青藏高原省級行政區劃圖Fig.1 The provincial administrative divisions of the QTP

1.2 數據來源

本研究選用的數據為：①降雨數據，利用由氣象站點日觀測數據而生成的年降水量數據集，獲取R因子的空間分布；②土壤數據，利用世界土壤數據庫（HWSD）中的土壤質地和土壤有機碳，獲取K因子的空間分布；③地形數據：利用DEM 數據生成坡度及坡長數據，獲取LS因子的空間分布；④土地利用及植被數據：利用NDVI 生成植被覆蓋度圖，結合土地利用類型圖，獲取B因子的空間分布；土地利用類型圖，結合地形數據，獲取T因子的空間分布。基礎數據、分辨率及來源如表1所示。

表1 本文使用的基礎數據Tab.1 Basic data used in this paper

1.3 研究方法

1.3.1 土壤侵蝕模數計算

Liu 等［16］建立了中國土壤流失方程（CSLE），該方程相比USLE 模型，充分考慮中國坡面的土壤侵蝕實際特征和土壤侵蝕植物覆蓋與生物防治措施，且已被眾多學者廣泛應用于中國地區的土壤侵蝕定量評價［17，18］，因此本文采用中國土壤流失方程CSLE 定量評價青藏高原農耕區的土壤侵蝕狀況，其基本結構形式如下：

式中，A為單位公頃的年土壤流失量，即土壤侵蝕模數，t/(hm2·a)；R為降雨侵蝕力因子，MJ·mm/(hm2·h·a)；K為土壤可蝕性因子，t·hm2·h/(hm2·MJ·mm)；L為坡長因子，無量綱；S為坡度因子，無量綱；B為植被覆蓋與生物措施因子，無量綱；E為工程措施因子，無量綱；T為耕作措施因子，無量綱。其各項因子的計算方法如下所述。

（1）降雨侵蝕力因子R。R因子的計算方法主要分為經典算法和簡易算法，前者計算復雜、數據要求高，難以推廣至大規模的土壤侵蝕研究。在簡易算法中，雨量模式由于數據的易獲取性而備受青睞。由于現有數據的限制，本文采用年雨量估算模型［19］計算R因子，其計算模型如下：

式中：R為降雨侵蝕力因子，MJ·mm/(hm2·h·a)；P為年降水量，mm。

（2）土壤可蝕性因子K。由于侵蝕-生產力計算模型（EPIC模型）所需數據較易獲取，本文中K因子采用EPIC 模型進行計算，其模型如下［20］：

式中：Sa為砂粒（0.05～2 mm）百分含量，%；Si為粉砂（0.002～0.05）百分含量，%；Ci為黏粒（＜0.002 mm）百分含量，%；C為有機碳含量，% ；EPIC模型計算K值為美制單位，t·arce·h/(100 arce·ft·tf·in)，需將其再乘以0.131 7 方轉化為國際制單位，t·hm2·h/(hm2·MJ·mm)。

（3）坡長因子L和坡度因子S。某一點的坡長在土壤侵蝕計算中是流域頂部至該點的坡長距離，即該點的上坡長，坡長因子L和坡度因子S采用劉寶元算法計算，公式如下：

式中：L為坡長因子；S為坡度因子；λ為坡長值，m；θ為坡度值，°；m為坡長系數。

（4）植被覆蓋與生物措施因子B。由于不同的土地利用類型和植被覆蓋度會致使植被對土壤侵蝕的抑制作用程度不同，故本文參考已有研究成果［21］，通過土地利用類型和植被覆蓋度確定B因子，如表2所示。

表2 B因子賦值Tab.2 The B factor values of different land-use types on the QTP

根據農村土地調查，坡度≤2°的耕地為平耕地，坡度＞2°的耕地為坡耕地。植被覆蓋度由NDVI計算得來，見下式［22］：

式中：C1為植被覆蓋度；NDVI為歸一化植被指數；NDVImax與NDVImin分別為研究區內NDVI的最大值與最小值。

（5）工程措施因子E。E因子依賴于工程措施的收集，但在大尺度研究區上，難以進行收集以致沒有基礎數據，故本文暫未考慮該因子對土壤侵蝕的影響。

（6）耕作措施因子T。現有研究中［23］，一般是根據不同的地形坡度下，等高耕作抑制土壤侵蝕的效果得到耕作措施因子T，參考相關研究成果［24］，耕地中與該因子的關系見表3。

表3 青藏高原耕地在不同坡度條件下的T因子Tab.3 The T factor values of the different slope of farmland on the QTP

1.3.2 土壤侵蝕強度變化綜合指數及減蝕貢獻度計算

穆天龍等［25］構建了土壤侵蝕強度綜合指數，用以研究不同土地利用類型土壤侵蝕程度，得到了廣泛認可和應用［26，27］。但其公式采用分級值計算，將所有處于某一級侵蝕的地塊按照同一侵蝕模數處理，與真實計算結果不符，易產生較大的系統誤差。因此，在穆天龍公式的基礎上，以單元值代替分級值，提出了土壤侵蝕強度變化綜合指數和區域土壤侵蝕強度變化綜合指數的概念及計算公式，以此評估土地利用變化對自身和農耕區土壤侵蝕的影響。為評估不同土地利用變化對區域水土流失防治的貢獻，提出了區域土壤減蝕貢獻度的概念及計算公式，以量化土地利用變化對區域土壤減蝕的貢獻。

（1）土壤侵蝕強度變化綜合指數，指某種土地利用變化類型引起的土壤侵蝕強度變化加權值，反映不同土地利用變化的自身土壤侵蝕的變化程度。其計算式為：

式中：j為土地利用變化類型；i為土地利用變化單元；n為土地利用變化單元數；ΔIj為j類土地利用變化下，土壤侵蝕強度變化綜合指數，t/(hm2·a)；ΔAi為第i個單元的土壤侵蝕模數變化，t/(hm2·a)；Sji為第i個單元的j類土地利用變化面積，hm2；Sj為j類土地利用變化總面積，hm2。

其中，ΔAi的計算如下：

通過CSLE 模型可知，在短時間內，土地利用變化實質上是通過引起植被覆蓋與生物措施因子B、工程措施因子E和耕作措施因子T的變化，進而引起土壤侵蝕的變化。因此，本文將2015年的降雨侵蝕力因子R代入2000年的CSLE 計算式［式（9）］，獲得在2000-2015年之間土地利用類型不變的假定情景下，2015年的土壤侵蝕模數情景估算值，將情景估算值與2015年的實際值作比較，差值即為土地利用變化引起的土壤侵蝕模數變化（式（10））。

式中：A估為土地利用不變情景下，2015年的土壤侵蝕模數情景估算值，t/(hm2·a)；R2015為2015年的R因子；B2000、T2000為2000年的B和T因子；其余因子為常數，同式（1）。

式中：ΔA為土地利用變化引起的土壤侵蝕模數變化，t/(hm2·a)；A2015為2015年的土壤侵蝕模數實際值，t/(hm2·a)。

（2）區域土壤侵蝕強度變化綜合指數，指區域內的土地利用變化引起的土壤侵蝕強度變化加權值，反映出不同土地利用變化類型引起的區域土壤侵蝕變化，以及該區域因土地利用變化形成的土壤侵蝕變化程度。其計算式為：

式中：m為土地利用變化類型數；ΔI為區域土壤侵蝕強度變化綜合指數，t/(hm2·a)；S為區域總面積，hm2。

（3）區域土壤減蝕貢獻度，指某種土地利用變化所引起的土壤侵蝕變化，占區域內各類土地利用變化引起的土壤減蝕總和的比例，能反映不同土地利用變化類型對研究區土壤減蝕的貢獻程度。計算式為：

式中：k為土地利用變化類型，k∈[0，j]；Cj為j類土地利用變化類型對區域土壤減蝕的貢獻度；ΔIk為k類土地利用變化下，土壤侵蝕強度變化綜合指數，其中ΔIk＜0，t/(hm2·a)；Sk為k類土地利用變化總面積，hm2。

2 結果與分析

2.1 土壤侵蝕空間分布

（1）各因子的空間分布。CSLE 各因子的空間分布如圖2，年降雨侵蝕力因子R如圖2（a）所示，R的均值為1 137.87 MJ·mm/(hm2·h·a)，最小的R因子(R= 41.63MJ·mm/(hm2·h·a))位于農耕區的西北部地區，主要是由于喜馬拉雅山和喀喇昆侖山阻擋了暖濕氣流，致使全年降雨極少［28］；最大的R因子(R=2 737.75 MJ·mm/(hm2·h·a)）位于農耕區的東南部地區，這是因為雅魯藏布江下游的峽谷向南展開的喇叭口，便于孟加拉灣暖濕空氣進入青藏高原東南部，形成降水豐沛的“舌狀地帶”［28］。

土壤可蝕性因子K如圖2（b）所示，K的均值為0.034 8 t·hm2·h/(hm2·MJ·mm)，可蝕性最低的土壤［K=0.025 51 t·hm2·h/(hm2·MJ·mm)]位于柴達木盆地，富含碳酸鈣、黏膠狀物質，使得土壤顆粒膠結在一起，不易分離。最容易被侵蝕的土壤(K=0.039 18 t·hm2·h/[hm2·MJ·mm)]位于阿里山地，冷漠土廣布，植被稀疏，常見裸地，易被侵蝕。K的分布和Teng 等人［29］、劉斌濤等人［20］的成果相似。

坡度坡長因子LS如圖2（c）所示，LS的均值為11.53，LS的最小值(LS= 0.56)位于柴達木盆地，海拔低且地勢平坦。LS的最大值(LS= 29.20)位于川西藏東高山深谷及云貴高原，處于中國地勢第二級階梯上，呈波狀起伏，山高谷深，落差大。

植被覆蓋與生物措施因子B如圖2（d）所示，B的最小值(B= 0.122)位于東喜馬拉雅南翼的熱帶雨林中，植被茂密，此外，在川西高山和云南高原里的熱帶雨林，B因子較小。B的最大值(B= 0.484)位于昆侖山北翼，植被覆蓋度低至15.60%～20.73%之間。

耕作措施因子T如圖2（e）所示，反映的是人類耕作措施對耕地侵蝕的減少，集中分布在農田，在一定范圍內，地勢越平坦，耕作措施效果越好，T因子越小。

圖2 CSLE各因子空間分布Fig.2 Distribution of the CSLE factors of the QTP

（2）土壤侵蝕空間分布特征。將CSLE 各因子的圖層相互乘積，計算得到2015年的青藏高原農耕區的土壤侵蝕模數，根據《土壤侵蝕分類分級標準》（SL190-2007）（表5），得到土壤侵蝕強度等級圖（圖3）。

圖3 2015年土壤侵蝕強度等級空間分布Fig.3 Spatial distribution of soil erosion intensity grades in 2015

結果表明，2015年，青藏高原農耕區的平均土壤侵蝕模數為195.61t/(hm2·a)，潛在土壤流失量達4.465 1 億t/a。農耕區土壤侵蝕呈明顯的空間分異特征，東南部地區侵蝕相對強烈，西北部地區侵蝕相對較輕，這與降雨、地形因子的空間分布規律相似，根據地理探測器的核心思想，即當自變量對于因變量的影響至關重要時，自變量與因變量的空間分布理應存在相似性，因此，在青藏高原農耕區，降雨、地形因素是土壤侵蝕的重要影響因素。

未利用地的土壤侵蝕模數最高（表4），但由于面積有限，造成的土壤侵蝕量并不大。耕地的土壤侵蝕模數相對較高，侵蝕量占總量的99.52%。草地土壤侵蝕模數相對較小，林地的土壤侵蝕模數最低，這與其他地區的相關研究結果相似［4，30］，說明土地利用對土壤侵蝕的影響規律在不同地區具有一定可復用性。

表4 不同土地利用類型的土壤侵蝕情況Tab.4 Soil erosion of different land use types

微度侵蝕的面積最大，占49.04%（表5），主要分布在西北部，強烈侵蝕面積最小，占4.07%，主要分布在中部，劇烈侵蝕分布在東南部，占19.12%，土壤侵蝕強度由西北向東南逐漸增大，這很可能是受降雨和地形因子的影響。

表5 研究區的各土壤侵蝕等級面積及占比Tab.5 Areas and proportions of different soil erosion grades

（3）結果合理性分析。將本文計算得到的2015年土壤侵蝕估算值與國家青藏高原科學數據中心的數據集［31］作比較，不同土壤侵蝕等級面積的重合率為90.91%，不重合面積來源于本文對劇烈等級面積的估算較大，其原因是本文采用年降雨量模型進行計算，未區分降雨和非降雨時段，一定程度上放大了土壤侵蝕強度，劇烈等級面積的差異比例為7.42%，但對劇烈及以下侵蝕等級的估算結果良好，其土壤侵蝕量相對誤差為0.48%。說明模型估算結果是合理的。

2.2 不同土地利用變化對土壤侵蝕的影響

基于公式（8）～（10），得到農耕區不同土地利用變化下的土壤侵蝕強度變化綜合指數，計算結果如圖4 所示。土壤侵蝕強度變化綜合指數大于0，表明該類土地利用變化引起增蝕，數值越大，增蝕作用越強；小于0，引起減蝕，數值越小，減蝕作用越強；絕對值越大，土壤侵蝕強度變化越大，說明該土地利用變化對土壤侵蝕的影響越大。

圖4 不同土地利用變化的土壤侵蝕強度變化綜合指數Fig.4 Comprehensive index of soil erosion intensity change for different land-use changes

就研究區而言，綜合指數大小排序為：退耕還林＜退耕還草＜未利用地轉耕地＜草地轉耕地＜林地轉耕地＜耕地轉未利用地，反映出退耕還林、退耕還草和未利用地轉耕地減弱了土壤侵蝕，減蝕作用依次降低；草地、林地開墾為耕地和耕地退化為未利用地增大了土壤侵蝕，增蝕作用依次增大，該規律在不同地域基本表現一致，說明土地利用是影響土壤侵蝕的關鍵因素。

綜合指數絕對值大小排序為：退耕還林＞耕地轉未利用地＞退耕還草＞林地轉耕地＞草地轉耕地＞未利用地轉耕地，其中，退耕還林的綜合指數絕對值最大，原因一是土地利用類型之間的自身土壤保持能力差異，林地最高，草地次之，耕地和未利用地分列三四，其中林地和耕地之間的差別最大；另一個原因是降水、地形的區域性顯著，退耕還林區有66.35%的面積分布在四川和云南，領先其他變化類型，而這些區域降雨量最大且地形有利于降水的沖刷，因此，當下墊面有了更強的降雨植物截留能力，能更大程度減輕降雨對土地的沖刷作用，從而引起了更大的減蝕。其他土地利用變化的綜合指數絕對值排序的原因亦然，主要是源于自身土壤保持能力差異和降水、地形的地區性差異。

與此同時，降水、地形及植被覆蓋的地區性差異使得土地利用變化對土壤侵蝕的影響具有一定的區域性特征：

（1）不同土地利用變化在不同區域的表現稍有差異，如在青海和西藏，退耕還草的綜合指數絕對值比退耕還林更大。

（2）相同土地利用變化在不同區域的指數不同，如，退耕還林在各省的綜合指數都不相等。

（3）甚至局部地區的土壤侵蝕變化方向與研究區不一致，如青海海西蒙古族藏族自治州和西藏阿里地區的草地轉耕地區，其綜合指數均小于0，減弱了土壤侵蝕，主要原因是原有草地的植被覆蓋度僅31.37%和20.48%，草地抗侵蝕能力大大下降，轉為耕地后，由于坡度較緩，僅1.2°和2.2°，耕作措施起到了水土保持作用。

2.3 不同土地利用變化對區域土壤侵蝕的影響及減蝕貢獻

基于2.2 節的土壤侵蝕強度變化綜合指數，按公式（11）計算農耕區不同土地利用變化類型與區域的土壤侵蝕侵蝕強度變化綜合指數，計算結果如圖5 所示。土壤侵蝕強度變化區域綜合指數大于0，表明該類土地利用變化引起區域的增蝕，數值越大，增蝕作用越強；小于0，引起減蝕，數值越小，減蝕作用越強。

圖5 不同土地利用變化的區域土壤侵蝕強度變化綜合指數Fig.5 Comprehensive index of regional soil erosion intensity change for different land-use changes

從區域角度分析，由研究區的區域土壤侵蝕強度變化綜合指數可知，農耕區在2000-2015年間的土地利用變化，使得研究區的土壤侵蝕強度減少245.54 t/(hm2·a)，較土地利用不變時下降55.66%，且各省農耕區的區域土壤侵蝕強度變化綜合指數均小于0，說明不論是在整體上還是局部地區，近年來農耕區土地利用規劃有助于水土流失的減少和生態環境的保護，其中四川省的區域變化綜合指數最小，這是因為退耕還林和退耕還草的土壤侵蝕強度變化綜合指數最小，即減蝕效果最佳，且四川省擁有45.74%的退耕還林和29.12%的退耕還草，因此四川省的減蝕幅度最大。

基于區域土壤侵蝕強度變化綜合指數，按公式（12）計算得不同土地利用變化的研究區減蝕貢獻度如表6所示。從土地利用變化角度分析，對研究區而言，不同土地利用變化對農耕區土壤減蝕的貢獻排序為：退耕還草＞退耕還林＞未利用地轉耕地＞林地轉耕地＞耕地轉未利用地＞草地轉耕地，該次序主要是受土壤侵蝕強度變化綜合指數和變化類型的面積的雙重影響，如草地相比于林地，缺少繁茂的樹冠、樹干、枝葉和發達的植物根系，抗侵蝕能力稍遜一籌，因此在同面積下，退耕還草的減蝕效果弱于退耕還林，因此退耕還草的土壤侵蝕強度變化綜合指數絕對值小于退耕還林，但退耕還草的面積是退耕還林的4.88倍，綜合計算，退耕還草貢獻了比退耕還林更多的土壤減蝕量。

表6 不同土地利用變化的研究區減蝕貢獻度Tab.6 Contribution of regional erosion reduction of different land-use changes

由前述可知，不論是對自身還是農耕區而言，退耕還林和退耕還草都起到了較大的減蝕作用，林地、草地開墾為耕地和耕地轉未利用地引起了6.53%土壤減蝕的流失，因此，為減少農耕區的水土流失，應側重于開展退耕還林、退耕還草，重點防止耕地退化，盡量避免將林地、草地開墾為耕地，而應多關注未利用地，科學論證未利用地開發為農用地的可行性。

3 結 論

本文研究了年際間土地利用變化對土壤侵蝕的驅動作用，與單一年份土地利用的土壤侵蝕特征相比［14，17］，揭示了土地利用變化對當下土壤侵蝕的影響，提出了土壤侵蝕強度變化綜合指數和區域減蝕貢獻度的概念和計算公式，量化分析了2000-2015年間土地利用變化對土壤侵蝕的驅動效果以及對區域土壤減蝕的貢獻度。主要結論如下：

（1）2015年，青藏高原農耕區的平均土壤侵蝕模數為195.61 t/(hm2·a)，潛在土壤流失量達4.460 1 億t/a。農耕區以微度侵蝕為主，受降雨、地形因子的影響，土壤侵蝕呈明顯的空間分異特征，由西北向東南逐漸增大。

（2）土地利用變化是引起土壤侵蝕變化的重要因素，退耕還林、退耕還草和未利用地轉耕地的土壤侵蝕強度變化綜合指數分別為-378.24、-118.54 和-5.44 t/(hm2·a)，減蝕作用依次降低；草地、林地開墾為耕地和耕地退化為未利用地的土壤侵蝕強度變化綜合指數分別為16.53、32.34和200.61 t/(hm2·a)，增蝕作用依次增大，效果的差異可能主要來源于土地利用自身的土壤保持能力差異和降水、地形的地區性差異，并且這些差異使得土地利用變化對土壤侵蝕的影響具有一定的區域性特征。

（3）青藏高原農耕區近年來的土地利用規劃有助于水土流失的防治，使得青藏高原農耕區的土壤侵蝕強度減少245.54 t/(hm2·a)，較土地利用不變時下降55.66%。不同土地利用變化對農耕區土壤減蝕的貢獻排序為：退耕還草＞退耕還林＞未利用地轉耕地＞林地轉耕地＞耕地轉未利用地＞草地轉耕地。□