李佳蕾,孫然好,熊木齊,楊國成
1 中國科學(xué)院生態(tài)環(huán)境研究中心城市與區(qū)域生態(tài)國家重點實驗室,北京 100085 2 中國科學(xué)院大學(xué),北京 100049
土壤侵蝕由于對經(jīng)濟和環(huán)境帶來嚴重的影響而成為了一個全球性的難題[1]。全球受水蝕影響的土地面積為10.94億hm2,其中751萬hm2受侵蝕嚴重[2]。世界上大部分農(nóng)業(yè)用地的土壤正在以每年13 t/hm2到40 t/hm2的速度流失,每年約有1000萬hm2的農(nóng)田因水土流失而退化[3]。我國是土壤侵蝕最嚴重國家之一,根據(jù)第3次水土流失普查,我國年土壤流失量約50億t[4],水蝕區(qū)平均侵蝕強度約為3800 t km-2a-1[5]。根據(jù)《第一次全國水利普查水土保持情況公報》,2011年,全國水力侵蝕面積達129.32萬km2。水土流失已經(jīng)成為制約我國經(jīng)濟社會可持續(xù)發(fā)展的重要因素之一[4]。
多數(shù)土壤侵蝕研究集中在小尺度區(qū)域,如坡地、坡面和小流域。這些研究多基于小區(qū)實驗[6],也有的研究利用模型計算,Karydas整理歸納了8類、82個水蝕模型[7]。大尺度土壤侵蝕定量研究主要是通過尺度上推和宏觀因子評價的方法進行[8]。定量評估土壤侵蝕的經(jīng)驗?zāi)P椭?最常用的是美國農(nóng)業(yè)部(United States Department of Agriculture, USDA) 提出的通用土壤流失方程(universal soil loss equation, USLE)和修正版的通用土壤流失方程 (revised universal soil loss equation, RUSLE)[9]。它以降雨侵蝕力、土壤可蝕性、地形、植被管理和保護措施為輸入,采用基于因子計算的方法估計每年因侵蝕造成的土壤損失[10]。自中國引入USLE和RUSLE模型以來,基于此模型的土壤侵蝕的研究在小區(qū)域和小流域上有著很多的研究應(yīng)用,而大流域和大尺度的研究相對較少。比如,傅伯杰等[11]利用USLE模型計算了我國黃土高原的水蝕,得出黃土高原每年水土流失量從2000年的3362 t/km2降低到了2008年的 2405 t/km2;還有其他學(xué)者利用USLE或RUSLE模型估算了中國南方丘陵山區(qū)[12]、廣東省[13]、江西省[14]和全中國[15]的年平均土壤侵蝕率。雖然這些也是較大尺度上的研究,但是大部分研究將中國的各個地區(qū)的土壤侵蝕計算方法歸一化,整個研究區(qū)用同種算法計算各個因子,忽略了各個因子因人類活動和地理分異帶來的區(qū)域差異。
本研究基于國際上較新發(fā)表的RUSLE模型中的C因子、R因子和P因子的計算方法,結(jié)合中國的實際情況對各因子進行修正,并應(yīng)用在中國土壤侵蝕的估算中。與以往的研究相比,利用大尺度高精度(30 m分辨率)數(shù)據(jù)、按區(qū)域性質(zhì)對中國不同地域進行分區(qū)運算,提高了土壤侵蝕計算結(jié)果的空間異質(zhì)性與精準度,并對2000年至2015年的中國土壤侵蝕量進行時空動態(tài)分析,探究中國土壤侵蝕強度變化的原因,為中國水土流失的治理與水土保持服務(wù)功能的保護提供理論基礎(chǔ)。
本研究所使用的主要基礎(chǔ)數(shù)據(jù)及其來源見表1。各個數(shù)據(jù)的中國分布情況見圖1。

表1 基礎(chǔ)數(shù)據(jù)

其中降雨數(shù)據(jù)是全國839個國家級站點的日數(shù)據(jù)經(jīng)SPSS軟件統(tǒng)計得出年的各個站點的年降水量,以及各氣象站點的經(jīng)緯度坐標和海拔高度。DEM數(shù)據(jù)的原始數(shù)據(jù)是USGS下載的中國分區(qū)數(shù)據(jù),再合成整個中國的30 m分辨率數(shù)據(jù)。中國年度植被指數(shù)(NDVI)空間分布數(shù)據(jù)是基于連續(xù)時間序列的SPOT/VEGETATION NDVI衛(wèi)星遙感數(shù)據(jù),采用最大值合成法生成的。土壤成分數(shù)據(jù)是在ISRIC中下載的2017年全球土壤的黏粒、粉粒、砂粒和有機碳的含量數(shù)據(jù),再提取中國的土壤成分數(shù)據(jù)。降水數(shù)據(jù)、NDVI、土地利用類型數(shù)據(jù)都選取了2000、2005、2010、2015年的數(shù)據(jù)。
1.2.1通用土壤流失方程(RUSLE)
RUSLE表達式為[16]:
A=R×K×LS×C×P
(1)
式中,A是模型預(yù)測的年土壤侵蝕量,單位為t km-2a-1。R是降雨侵蝕力因子,單位為MJ mm hm2h-1a-1,是降水產(chǎn)生的徑流對土壤造成侵蝕的動力指標[17],降雨的強度和持續(xù)時間對侵蝕有重要影響。LS為地形因子,無量綱;L即坡長因子,是指標準化到22.13 m坡長上的土壤侵蝕量;S即坡度因子,是指標準化到5.14°坡度下的土壤侵蝕量[18];一般的小尺度研究直接利用實測地形數(shù)據(jù),大尺度研究會利用DEM數(shù)據(jù)提取信息來計算LS因子。K是土壤可蝕性因子,單位為thm2h hm-2MJ-1mm-1,其反映了土壤對侵蝕營力分離和搬運作用的敏感性[19-20]。C是植被覆蓋與管理因子,無量綱,表示植被覆蓋和管理措施對土壤侵蝕的作用[21]。P因子即水土保持措施因子,無量綱,指在特定水土保持措施的土壤流失與起伏地耕作的相應(yīng)土壤流失之比[16]。C和P因子可反映人為控制土壤侵蝕的作用。
1.2.2基于氣候分區(qū)的R因子修正
Naipal等[22]基于美國降雨站點監(jiān)測值,擬合降雨侵蝕力因子R值和年降水量(P)、高程(Z)、日降水強度(SDII)的回歸方程,得出修正的不同氣候區(qū)降雨侵蝕力R計算方法:
R=f(P,Z,SDII)
(2)
式中,日降水強度SDII等于每日降水量≥1 mm的總降雨量與天數(shù)的比值。不同的氣候區(qū)有不同的計算公式,其氣候區(qū)的劃分標準為Koppen-Geiger氣候分區(qū)法[23]。
本研究將利用Naipal的方法,將中國不同的氣候區(qū)R值分別進行計算。對于精度較差的地區(qū),以及無SDII數(shù)據(jù)的氣候區(qū),則利用Renard and Freimund指數(shù)[16](表2)。
首先利用全國800多個站點的日降水數(shù)據(jù)計算出各個站點的年平均降雨侵蝕力,再根據(jù)李璐等基于地統(tǒng)計學(xué)的降雨侵蝕力插值方法的研究,選用模擬精確度較好的普通克里金插值方法來進行空間插值[24]推測全國的R值。

表2 各個氣候區(qū)的R因子的計算方法
(1)氣候分區(qū)的命名方式為字母縮寫組合,字母縮寫意義為B:干旱區(qū)Arid,C:溫和區(qū)Temperate,D:寒冷區(qū)Cold,W:沙漠區(qū)Desert,S:干旱草原Steppe,k:Cold,f:Withoutdryseason,w:DryWinter,a:HotSummer,b:WarmSummer
1.2.3基于土地覆被的C因子修正方法
利用Borrelli計算C因子的方法[25],并結(jié)合中國實際土地覆蓋類型和農(nóng)作物種類對中國的耕地和非耕地的C因子分別進行調(diào)整與計算(C因子的取值見表3)。

表3 不同土地覆蓋類型的C值
植被覆蓋因子C和農(nóng)作物種類緊密相關(guān)[26-27],所以本研究中耕地采用國家統(tǒng)計局發(fā)布的各省耕地主要作物種類和播種面積,將已發(fā)布的作物歸為10類,通過以下公式計算過去全國耕地的C值:
(3)
式中,Ccropn是作物n的C值,%RegionCropn是作物n的播種面積占各省總面積的比例。
植被覆蓋因子C取決于植被覆蓋度,本研究結(jié)合已有文獻中非耕地的各種植被覆蓋類型的C值[25, 28],利用土地利用數(shù)據(jù)和NDVI數(shù)據(jù),計算非耕地的C值。
CNonArable= Min(CNA) + Range(CNA) × (1-Fcover)
(4)
Fcover=VFC = (NDVI- NDVImin)/(NDVImax- NDVImin)
(5)
式中,Min(CNA)是CNA取值范圍的最小值,Range(CNA)是CNA的最大值與最小值的差,Fcover是植被覆蓋度。
1.2.4基于地形特征的P因子修正方法
在USLE/RUSLE支持下的P因子在大尺度的土壤侵蝕風險建模中很少被考慮,因為很難對大面積進行估算[29]。Xiong等依據(jù)大量文獻總結(jié)[10, 30-31],歸納不同坡度耕地及不同的水土保持措施的P值差異。針對本研究,將Xiong的P因子賦值方法結(jié)合歷年土地利用類型,將水平稻田區(qū)域賦值為0.2,其他農(nóng)田區(qū)域根據(jù)坡度為P因子賦值。其中坡度小于等于10°的農(nóng)田P取值為0.5,坡度大于10°小于等于25° 的P取值為0.6,坡度大于25°小于等于45°的P取值為0.8,坡度大于45°的P取值為1。
1.2.5其他因子的計算方法
DEM的單個柵格大小會影響計算精確度,柵格越大坡度減小,估算的土壤侵蝕量隨之減小[32-33]。基于數(shù)據(jù)可靠性和計算可行性,本研究利用30 m分辨率的中國DEM數(shù)據(jù),選擇如下公式[34- 36]計算L因子:
L=(γ/22.3)m
(6)
m=β(1+β)
(7)
(8)
式中,γ為坡面的水平投影長度,為方便計算直接用柵格長度(30 m)代替;m為坡長因子指數(shù),β為細溝侵蝕和面蝕的比值[37],θ為坡度。
對于S因子,則利用劉寶元提出的CSLE模型中計算S因子的方法[38],對不同坡度段分別進行計算:

(9)
RUSLE中的K值是對具體土壤的由試驗確定的定量數(shù)值,劉寶元提出了適用于我國的土壤可蝕性測定方法,即在15°坡度、20 m坡長、清耕休閑地上,單位降雨侵蝕力所引起的土壤流失量[20]。但本研究尺度較大,因此土壤可蝕性因子K利用EPIC模型[39]計算。其公式為:

(10)
(11)
式中,SAN為砂粒含量百分數(shù),SIL為粉粒含量百分數(shù),CLA為黏粒含量百分數(shù),C為有機碳含量百分數(shù)。
將年均土壤侵蝕量的柵格圖層按照時間序列進行一元線性回歸分析,采用最小二乘法計算2000—2015年土壤侵蝕變化斜率,以得出2000—2015年中國土壤侵蝕強度的時間變化趨勢,斜率為正則表示土壤侵蝕有增加的趨勢,斜率為負則表示土壤侵蝕有減小的趨勢。斜率的計算公式如下:
(12)
式中,θslope是回歸方程的斜率,n為所計算的年數(shù),Ai為第i年的平均土壤侵蝕量。變化趨勢的顯著性采取F檢驗,以檢驗中國土壤侵蝕的時間變化規(guī)律是否顯著。該統(tǒng)計量的計算構(gòu)造為:
(13)
利用上述計算各因子的方法得出的各因子的空間分布圖如圖2。其中,C因子和R因子的年際變化較明顯,歷年的土壤侵蝕的變化即體現(xiàn)在降雨、植被覆蓋及土地利用中。由于本研究所計算的R因子、C因子和P因子利用了較新的算法,較之前的國內(nèi)RUSLE模型的應(yīng)用研究有了一定程度的細化。
2000年、2005年、2010年和2015年中國水蝕區(qū)C因子的年平均計算結(jié)果為0.1464、0.1482、0.1498和0.1538,2015年較2000年提高了5.05%。而計算后的2000年、2005年、2010年和2015年中國水蝕區(qū)年平均R因子分別為2903.79、3034.62、3515.44、3444.77 MJ mm hm2h-1a-1,2015年比2000年增加了18.63%,說明中國降水帶來的侵蝕潛力有增強的趨勢。


圖2 中國RUSLE(1)模型各因子空間分布(2)Fig.2 Distribution of RUSLE factors in China(1)RUSLE:修正版的通用土壤流失方程,revised universal soil loss equation。(2)圖例中,降雨侵蝕力因子單位為MJ mm hm2 h-1 a-1,土壤可蝕性因子單位為t hm2 h hm-2MJ-1mm-1
在將所有的RUSLE模型中的因子利用重采樣統(tǒng)一為100 m分辨率后,將各因子相乘計算出各年土壤侵蝕量。其中,2000年、2005年、2010年、2015年中國的年平均土壤侵蝕量分別為38.63、37.35、49.03、47.84 t km-2a-1。根據(jù)中華人民共和國水利部發(fā)布的《土壤侵蝕分類分級標準》(SL190—2007)[40],將計算后的中國土壤侵蝕度分成6個等級,即微度、輕度、中度、強烈、極強烈和劇烈,其空間分布見圖3,分級統(tǒng)計情況如表4。

圖3 中國2000—2015年土壤水力侵蝕程度空間分布Fig.3 Distribution map of water erosion in China

表4 中國2000—2015年土壤侵蝕強度統(tǒng)計
從表4可以看出,中國水蝕區(qū)大多(60%以上)的空間為微度侵蝕,微度侵蝕的變化幅度不大,但2010年的微度侵蝕的面積占比是4年中最小的;劇烈侵蝕的空間占比在2005年后有所增加,其中2010年劇烈侵蝕的面積比最高,比2000年增加了2.00%。2000年以來,全國土壤侵蝕總面積在減少。具體來講,侵蝕輕微區(qū)在增加、侵蝕劇烈區(qū)也在增加,說明了生態(tài)修復(fù)的效果體現(xiàn)在了侵蝕總面積和侵蝕改善區(qū)的面積增加;而由于部分區(qū)域的耕地增加及降水增加,也會導(dǎo)致個別區(qū)域的侵蝕強度增加,是后續(xù)需要進一步加強侵蝕管理和控制的區(qū)域。
中國水蝕熱點分布在中國南方,呈點狀分布分散于南方各地;北方雖有水蝕較強的區(qū)域,但水蝕劇烈的地方相對較少。在南方,云貴高原與四川盆地之間的區(qū)域水蝕劇烈,尤其分布于貴州省。通過研究其地形和土地覆被,分析其原因可能是其處于喀斯特地貌區(qū),海拔較高(1000 m以上)、坡度較大,土地覆被類型大多為耕地,人口密集、人為干擾較多,土壤保持力較為低下[41]。中國東南部的丘陵區(qū)土壤侵蝕劇烈的分布也比較廣,但大多為點狀分布,一般分布于坡耕地上,這類地區(qū)降水豐沛、坡度較大、人為干擾較多,水土流失極易發(fā)生。
北方土壤侵蝕熱點集中于黃土高原、山東丘陵區(qū)、東北大興安嶺和西北的昆侖山與塔里木盆地的交界處。黃土高原由于其特殊的地貌一直受到很多學(xué)者的關(guān)注[42],其不穩(wěn)定的土壤特性和稀少的植被覆蓋一直是其成為水土流失熱點的原因。大興安嶺和山東丘陵的海拔都較周圍高,土壤侵蝕的敏感性也較強。而中國西北部的昆侖山與塔里木盆地之間的區(qū)域,其坡度較大,有較多的冰雪消融和降雨帶來的徑流,植被以草地為主且較為稀疏,城市分布較多人為影響較大,其水蝕也比較劇烈。
將年均土壤侵蝕量按照時間序列進行一元線性回歸分析后,得出2000年到2015年中國土壤侵蝕強度的時空變化格局(圖4),其中有28.72%的區(qū)域基本穩(wěn)定;37.22%的區(qū)域土壤侵蝕有增加趨勢,34.04%的區(qū)域為有減小的趨勢。根據(jù)1.3的方法將中國年均土壤侵蝕演變顯著程度分為5類(圖4)并進行統(tǒng)計,其中土壤侵蝕有極顯著減小趨勢的區(qū)域面積為3417 km2,占總有效統(tǒng)計面積的0.74%;有顯著減小趨勢的區(qū)域面積為132593 km2,占2.88%;有極顯著增加趨勢的區(qū)域面積為20823 km2,占0.45%;有顯著增加趨勢的區(qū)域面積為82783 km2,占1.80%。土壤侵蝕顯著增加的區(qū)域分布在中國長江中下游平原、江南丘陵、云貴高原、四川盆地、東北平原和黃土高原的西北部;顯著改善的區(qū)域有黃土高原南部、秦嶺地區(qū)和東南沿海地區(qū)。

圖4 2000—2015年中國土壤侵蝕的時空變化格局Fig.4 Spatio-temporal variation pattern of soil erosion in China, 2000—2015
由于年際土壤侵蝕的變化體現(xiàn)在C因子和R因子中,可在通過對C因子和R因子進行一元線性回歸分析后得出變化斜率(圖4),分析出有明顯土壤侵蝕強度變化區(qū)域的變化因素。通過圖4,可以看出中國C因子有明顯增大的區(qū)域有東北平原、新疆的準噶爾盆地和吐魯番盆地、南方丘陵區(qū),以及黃土高原東部。分析其原因,東北和西北地區(qū),主要是自然植被的退化和農(nóng)業(yè)逐漸發(fā)達。隨著南方經(jīng)濟的發(fā)展,我國農(nóng)業(yè)逐漸向新疆和東北地區(qū)轉(zhuǎn)移[43],天然植被被耕地取代,耕地面積逐漸擴大,糧食作物種類與耕種面積變化,植被因子也隨之增大,土壤侵蝕也逐漸增強。而部分地區(qū),如黃土高原南部、秦嶺一帶,其土壤侵蝕顯著改善,原因是陜西省、山西省、河南省實施了退耕還林和修建梯田等水土保持措施,其植被覆蓋增加、C因子減小,土壤侵蝕的情況有所好轉(zhuǎn);并且部分區(qū)域伴隨著R因子的減少,土壤侵蝕力也有所下降。南方東南沿海區(qū)散狀分布的C因子的改變,推測與部分農(nóng)田的荒廢或退耕還林的措施有關(guān)。而R因子明顯增大的區(qū)域則在中國長江中下游平原和兩廣丘陵與云貴高原之間,這些地區(qū)年降水量在逐漸增大,這也是其水蝕強度增大的主要原因。
本研究計算的2000年、2005年、2010年、2015年中國的平均土壤侵蝕率分別為38.63、37.35、49.03、47.84 t km-2a-1,其中有2.25%的土地土壤侵蝕情況明顯惡化。整體上中國的平均土壤侵蝕量有所增加,但全國土壤侵蝕總面積在減少,其原因主要為C因子和R因子的變動,其表現(xiàn)是部分地區(qū)降水量的增減、植被覆蓋的變化和耕地糧食作物種類和播種面積的變化。中國社會在2000年到2015年處于快速發(fā)展的時期,人口增多、城市化加速,一方面對建筑用地的需求增多,另一方面對糧食的需求也在不斷增多,面臨著糧食安全的問題[44]。為了達到供求平衡,我國的糧食生產(chǎn)中心逐漸北移[43],農(nóng)田的開墾和土壤肥力的消耗使得東北和西北地區(qū)的土壤侵蝕強度增大。為了改善中國的水土流失,我國實施了一系列修復(fù)措施,有34.04%的土地的土壤侵蝕強度有所減小。比如我國水土流失最嚴重的黃土高原區(qū),在1999年啟動退耕還林還草工程后[45],土壤侵蝕現(xiàn)象有了一定的改善。
大尺度土壤侵蝕的定量計算與分析在實現(xiàn)空間準確性和真實性上有一定的困難。雖然RUSLE模型設(shè)計的初衷是用于坡面土壤侵蝕的計算,但由于其簡潔性和參數(shù)的可獲得性在大尺度土壤侵蝕估算中有了較為廣泛的應(yīng)用。然而大尺度的計算會提高數(shù)據(jù)粗糙度,對算法進行合理的修正可以彌補數(shù)據(jù)不足帶來的誤差。本研究所利用的C因子、R因子和P因子的算法,將不同土地覆蓋類型和氣候帶的參數(shù)進行區(qū)分,提高了計算結(jié)果的空間異質(zhì)性與區(qū)域準確度。對于R因子,本研究所用的Naipal法對不同氣候區(qū)進行了分異,并針對不同的氣候區(qū)利用了不同的算法,同時也結(jié)合了海拔和降雨集中性的指標,使得計算結(jié)果在中國這個大尺度研究區(qū)中更能體現(xiàn)出異質(zhì)性。目前在國內(nèi)廣泛應(yīng)用的方法是章文波利用日降水量計算R因子的方法[46]。將易玲[47]利用章文波法計算的中國年平均R因子,與本研究應(yīng)用的Naipal法的計算結(jié)果相比,章文波法計算出的中國平均年降雨侵蝕量比本研究的結(jié)果偏大。其原因可理解為丘陵地區(qū)的R因子在調(diào)整后比未調(diào)整的R因子有所降低[22]。而目前國內(nèi)廣泛應(yīng)用的C因子的計算方法一般是兩種:一種是利用NDVI確定C因子,而另一種是根據(jù)土地覆蓋類型數(shù)據(jù),對不同的土地覆蓋類型賦予不同的C值。而本研究利用Panagos計算歐洲C因子的算法,區(qū)分耕地和非耕地并將土地利用類型、NDVI和土地覆蓋度數(shù)據(jù)融合,得出更加精準的中國C因子分布。P因子則根據(jù)農(nóng)田類型和坡度賦值,體現(xiàn)了人類土壤保持措施的作用。高分辨率的地形因子(30m分辨率)使得坡度因子的精確度達到提升。
雖然數(shù)據(jù)的空間分辨率的提升和算法的改進使得結(jié)果有了更高的準確度,但由于空間尺度大,土壤侵蝕具有空間尺度效應(yīng)[6],大尺度模型估算與小尺度的土壤調(diào)查數(shù)據(jù)會一定程度的偏差。大尺度的模型計算的原數(shù)據(jù)精度與計算方法的不同也會導(dǎo)致不同研究的結(jié)果不同。例如江西省的年土壤侵蝕率,本研究計算的2015年平均土壤侵蝕率為9060 t/km2,而周夏飛等[14]利用USLE模型在江西省尺度下定量評估的結(jié)果為6375 t/km2,Teng等[15]利用RUSLE模型在全國尺度下計算結(jié)果為 92 t/km2,陳思旭等[12]基于RUSLE模型在南方丘陵的尺度下計算的結(jié)果為354 t/km2。不同研究的研究尺度不同、計算方法不同、數(shù)據(jù)來源不同,結(jié)果差異較大,但分布規(guī)律與變化規(guī)律大致相似。
中國土壤侵蝕研究應(yīng)進一步優(yōu)化各個因子的計算方法、提高數(shù)據(jù)精度,同時修正風力侵蝕和凍融侵蝕區(qū)的土壤侵蝕估算方法,提高其他類型土壤侵蝕量估算的空間準確性。這將為土壤侵蝕的大尺度估算提供更準確的地域信息,也可以為中國不同區(qū)域制定針對性的水土保持措施提供理論依據(jù)。
(1)2000年、2005年、2010年、2015年中國的平均土壤侵蝕率為38.63、37.35、49.03、47.84,總體上中國的土壤侵蝕量有所增加,但近期國家對水土保持的重視以及實施了一系列保護措施,土壤侵蝕的強度有所緩解。
(2)土壤侵蝕量較高的地區(qū)集中在中國南方長江中下游平原區(qū)、云貴高原、黃土高原區(qū)、昆侖山山麓,占統(tǒng)計總面積的9.65%。而土壤侵蝕量明顯增大的區(qū)域有西北農(nóng)田區(qū)、四川盆地、兩廣丘陵西部、長江中下游平原和東北平原,面積達10.36×104km2。其原因有自然的降雨增多,也有開墾農(nóng)田、農(nóng)作物更替等人為行為干擾。
(3)農(nóng)業(yè)活動增強了西北地區(qū)和東北地區(qū)的土壤侵蝕;南方城鎮(zhèn)化破壞了天然植被,導(dǎo)致土壤侵蝕呈現(xiàn)分散狀的增加;但黃土高原南部和東部的水土修復(fù)工程的實施使其土壤侵蝕有所緩解。