趙宇豪,李崇明,何海珊,王清濤,吳健生
(1.北京大學深圳研究生院 城市規(guī)劃與設計學院 城市人居環(huán)境科學與技術重點實驗室,廣東 深圳 518055;2.北京大學城市與環(huán)境學院 地表過程分析與模擬教育部重點實驗室,北京 100871;3.河北工程大學 園林與生態(tài)工程學院,河北 邯鄲 056038)
土地可為人類提供多種生態(tài)服務,但目前全球約30%的土地發(fā)生了退化,預計到2050年95%的地球表面將受到土地退化的影響,聯(lián)合國已將土地退化防治納入可持續(xù)發(fā)展目標[1-3]。 土壤侵蝕是導致土地退化的主要原因之一,大面積的土壤侵蝕使土地嚴重退化、生產(chǎn)力降低、糧食減產(chǎn),并造成荒漠化、泥石流、河道及水庫淤積、水質(zhì)惡化等一系列后果[4-8]。 隨著全球氣候進一步變暖,土壤侵蝕可能還將加劇并威脅人類的生存與可持續(xù)發(fā)展[9-11]。 因此,土壤侵蝕的相關研究受到各國學者的廣泛重視,目前研究的熱點內(nèi)容有土壤侵蝕過程與機理、土壤侵蝕模型修正、土壤侵蝕時空分異特征、土壤侵蝕與景觀格局的關系、土壤侵蝕驅動因素等[12-14]。
陜西省是我國水土流失嚴重省份和黃河的主要產(chǎn)沙區(qū)[15-17],許多學者對陜西省土壤侵蝕進行了相關研究,如:鐘莉娜等[18]探討了陜北黃土丘陵溝壑區(qū)降雨和土地利用格局對土壤侵蝕的影響,張雪才等[19]對陜西省渭河流域的水土流失風險進行了評估,王濤等[20]探究了陜北無定河流域土壤侵蝕時空演變特征,程金文等[21]分析了陜南地區(qū)1960—2014年降雨侵蝕力的變化特征,火紅等[22]分析了延安市實施退耕還林前后各土地利用類型分布的差異及轉移特征。 綜上所述,對陜西省土壤侵蝕的研究主要集中在部分地市和黃河支流,鮮有以全省為尺度進行研究。
土壤侵蝕驅動因子識別是因地制宜制定治理策略的基礎,以往多依據(jù)RUSLE 模型計算的土壤侵蝕量與外部影響因子相關聯(lián)進而識別土壤侵蝕驅動因子[23-26]。 本研究提出一種通過比較RUSLE 模型中各土壤侵蝕因子的變化率來識別土壤侵蝕變化驅動因子的方法(稱為最大變化率法),并以陜西省為例分析2000—2015年土壤侵蝕時空變化的主要驅動因子,以期為陜西省及黃河流域土壤侵蝕治理提供理論依據(jù)。
陜西省位于東經(jīng)105°29′—111°15′、北緯31°42′—39°35′,總面積約20.53 萬km2,轄10 個地級市,境內(nèi)有高原、山地、平原和盆地等多種地貌類型,土壤類型多樣,由北向南依次穿越中溫帶、暖溫帶和北亞熱帶等3個氣候帶,涵蓋半干旱區(qū)、半濕潤區(qū)和濕潤區(qū)[27]。 北山和秦嶺將全省分為陜北黃土高原區(qū)、關中平原區(qū)和陜南秦巴山區(qū)等三大區(qū)域(分別簡稱陜北、關中、陜南,見圖1),其中:陜北包括延安、榆林兩市,面積約8.00 萬km2,位于黃土高原中部,地勢西北高東南低,西北部為風沙區(qū),中南部大都為丘陵溝壑區(qū),處于暖溫帶大陸性季風半濕潤氣候向溫帶半干旱氣候的過渡區(qū),年均氣溫8~12 ℃,年均降水量350~600 mm[28-29];關中包括銅川、渭南、咸陽、寶雞、西安等五市,面積約5.53 萬km2,地勢西高東低,包含秦嶺山地、渭河平原以及渭北黃土臺塬,屬于大陸性季風氣候區(qū),年均氣溫6~13 ℃,年均降水量500~800 mm[30];陜南包括漢中、安康、商洛三市,面積約7.00 萬km2,北靠秦嶺、南依巴山,具有“兩山夾一川”的地貌結構,西部屬于北亞熱帶季風氣候區(qū),東部為北亞熱帶氣候與暖溫帶氣候過渡區(qū),年均氣溫12~15 ℃,年均降水量700~1300 mm[31-32]。

圖1 陜西省分區(qū)及主要氣象站分布
研究所需數(shù)據(jù)及其來源見表1。

表1 數(shù)據(jù)來源
氣象數(shù)據(jù)包括省內(nèi)21 個氣象站(分布情況見圖1)1999—2001年、2004—2006年、2009—2011年、2014—2016年的逐月降水量(其中1999—2001年只有20 個氣象站的數(shù)據(jù)),把上述4 個時段的數(shù)據(jù)平均值分別作為2000年、2005年、2010年、2015年的數(shù)據(jù)(用當年及其前后各一年數(shù)據(jù)的平均值以避免研究數(shù)據(jù)受極端降水的影響)。 土壤屬性數(shù)據(jù)包括沙粒、粉粒、黏粒含量和有機碳含量。 歸一化植被指數(shù)NDVI和土地利用類型包括2000年、2005年、2010年、2015年逐月數(shù)據(jù),年度NDVI值取當年3月、6月、9月、12月的平均值,土地利用類型分為耕地、林地、草地、水域、建設用地、未利用地等6 類。 對所有柵格數(shù)據(jù)的分辨率,采用最鄰近法重采樣至100 m。
1.3.1 RUSLE 模型及各因子的計算
RUSLE 模型[23]也稱修正的通用土壤流失方程,形式為

式中:A為土壤侵蝕模數(shù), t/(hm2·a);R為降雨侵蝕力,MJ · mm/(hm2·h·a);K為土壤可蝕性,t·hm2·h/(hm2·MJ·mm);LS為地形因子(其中L為坡長因子、S為坡度因子);C為植被覆蓋因子;P為水土保持措施因子。
(1)降雨侵蝕力R。 降雨侵蝕力R反映降雨對土壤剝離、搬運和沖刷能力的大小,采用基于月降雨量的經(jīng)驗公式[33]來計算:

式中:ri為年內(nèi)第i個月的降雨量,mm;r為年降雨量,mm。
區(qū)域(全省或分區(qū))年降雨侵蝕力的計算:分別計算研究區(qū)各氣象站各年份的R值,采用協(xié)同克里金法以高程為協(xié)變量進行插值,得到研究區(qū)每年的R值。
(2)土壤可蝕性K。 土壤可蝕性K反映土壤被外部營力分離、搬運的難易程度,其值越大表示土壤越容易被侵蝕,采用EPIC 模型[34]來計算:

其中

式中:SAN、SIL、CLA和Carbon分別為土壤中沙粒、粉粒、黏粒含量和有機碳含量,%。
(3)地形因子LS。 地形因子LS反映坡長和坡度對土壤侵蝕的影響,依據(jù)DEM 數(shù)據(jù),采用如下公式[35]計算:

其中

式中:LSi為柵格i(i為柵格編號)的地形因子;Si為柵格i的坡度因子;Ai-in為柵格i入口處的匯流面積,m2;D為柵格邊長,m;ai為柵格i的坡向;m為坡長指數(shù)。
(4)植被覆蓋因子C。 植被覆蓋因子C反映植被對土壤侵蝕的抑制作用,取值范圍為0~1,數(shù)值越大說明植被的抑制作用越小、土壤侵蝕越嚴重,依據(jù)研究區(qū)歸一化植被指數(shù)NDVI,參考蔡崇法等[36]的研究對C進行分段賦值:

式中:f為植被覆蓋度,%;NDVIsoil、NDVIveg分別為純裸土柵格、純植被柵格的NDVI值。
(5)水土保持措施因子P。 水土保持措施因子P用來衡量特定水土保持措施對土壤侵蝕的影響,取值范圍為0~1,即采取水土保持措施后的土壤侵蝕量與無水土保持措施的順坡耕作時的土壤侵蝕量的比值。參考有關學者的研究[8,32,37],本研究設定耕地、林地、草地、水域、建設用地、未利用地的P值分別為0.4、1.0、1.0、0.0、0.0、1.0。
1.3.2 驅動因子識別
依據(jù)RUSLE 模型中各因子的變化率來識別土壤侵蝕變化的主要驅動因子。 用變量X表示RUSLE 模型因子R、K、LS、C、P中的任意一個,采用如下公式計算各土壤侵蝕因子時空變化率,進而識別土壤侵蝕變化的主要驅動因子。土壤侵蝕因子時間變化率計算公式為

式中:Xi、Xj分別為基準年份i、目標年份j的土壤侵蝕因子值;t為基準年份i至目標年份j的時段長(即ji),a;TX為土壤侵蝕因子從基準年份i到目標年份j的年均變化率,%,TX為正值表示土壤侵蝕因子數(shù)值增大,TX為負值表示土壤侵蝕因子數(shù)值減小。
土壤侵蝕因子空間變化率計算公式為

式中:XC為省內(nèi)某地市或某區(qū)域土壤侵蝕因子值;XP為土壤侵蝕因子的全省平均值;SX為土壤侵蝕因子空間變化率,%,當某地市或某區(qū)域土壤侵蝕因子值大于全省平均值時SX為正值,反之則為負值。
在土壤侵蝕各因子中,把時間變化率、空間變化率絕對值最大的因子分別作為土壤侵蝕變化的主要時間驅動因子、主要空間驅動因子。 若多個因子的變化率絕對值與主要驅動因子變化率絕對值相差不超過10%時,則認為土壤侵蝕變化由多因子共同驅動。
(1)陜西省各級土壤侵蝕空間分布格局。 依據(jù)土壤侵蝕模數(shù),按照水利部頒布的《土壤侵蝕分類分級標準》(SL 190—2007)[38],把土壤侵蝕強度分為微度、輕度、中度、強烈、極強烈、劇烈等6 個等級,2000年、2005年、2010年、2015年陜西省各級土壤侵蝕分布情況見圖2。 由圖2可以看出:陜西省土壤侵蝕以微度和輕度為主,劇烈、極強烈和強烈侵蝕面積相對較小;從時間上看,2000—2015年土壤侵蝕強度呈逐漸降低的趨勢;從空間上看土壤侵蝕等級空間分布格局變化不大,陜北以劇烈和極強烈侵蝕為主,關中以微度和輕度侵蝕為主、部分區(qū)域存在劇烈和極強烈侵蝕,陜南以微度和輕度侵蝕為主。

圖2 典型年份陜西省各級土壤侵蝕分布情況
(2)陜西省各級土壤侵蝕面積占比及轉化情況。2000年、2005年、2010年、2015年各級土壤侵蝕面積占比及轉化情況見圖3(圖中百分數(shù)為各級土壤侵蝕面積占比)。

圖3 陜西省各級土壤侵蝕面積占比及轉化情況
從2000—2015年總體來看:高強度土壤侵蝕面積占比下降、中度及以下土壤侵蝕面積占比上升,劇烈和極強烈侵蝕面積占比分別下降了4.70、0.79 個百分點(這兩級土壤侵蝕面積及其占比較小,雖然從絕對數(shù)值上看下降不多,但是下降比例即減少面積占原面積的比例分別達到了34.81%和9.71%),強烈侵蝕面積占比基本未變,中度、輕度、微度侵蝕面積占比分別上升0.54、1.71、3.20 個百分點,表明雖然土壤侵蝕現(xiàn)象依然存在、土壤侵蝕總面積基本沒變化,但是土壤侵蝕強度整體上趨于降低。
從3 個時段來看:2000—2005年,微度和劇烈侵蝕面積占比下降(分別下降3.78、0.38 個百分點),其他強度等級的侵蝕面積占比上升,總體上土壤侵蝕強度呈現(xiàn)提高趨勢,具體來說,劇烈侵蝕面積有12.72%向極強烈侵蝕轉化,各級土壤侵蝕面積均存在向相鄰等級轉化的情況,除輕度侵蝕有9.15%向微度侵蝕轉化外,鮮有其他強度向微度侵蝕轉化的情況,微度侵蝕面積有84.56%侵蝕強度保持不變、有11.77%向輕度侵蝕轉化、有3.67%向中度侵蝕轉化;2005—2010年,各級土壤侵蝕面積的轉化與上一時段的轉化明顯不同,除微度侵蝕面積占比上升外,其他強度等級的土壤侵蝕面積占比均呈現(xiàn)下降趨勢,并且高強度等級向低強度等級轉化的數(shù)量明顯增加,尤其是輕度侵蝕向微度侵蝕的轉化(有28.89% 轉化為微度侵蝕),表明該時段土壤侵蝕強度等級整體呈現(xiàn)降低趨勢;2010—2015年,微度侵蝕面積占比略微上升,輕度和中度侵蝕面積占比呈小幅上升趨勢,而劇烈、極強烈和強烈侵蝕面積占比呈下降趨勢(其中劇烈侵蝕面積占比下降較為顯著),整體上土壤侵蝕強度等級進一步降低。不同時段各級土壤侵蝕面積轉化情況較為復雜,表明大范圍的土壤侵蝕控制是一項復雜的系統(tǒng)工程,在控制高強度等級土壤侵蝕的同時,也應防止低強度等級的土壤侵蝕向高強度等級轉化。
(3)各地區(qū)各級土壤侵蝕面積占比變化情況。 陜北、關中和陜南各級土壤侵蝕面積占比變化情況見圖4。 總體上看,2000—2015年陜西省土壤侵蝕強度等級均呈現(xiàn)從高強度向低強度轉化的趨勢,但各地區(qū)轉化方式不盡相同。 陜北2000—2015年除劇烈和極強烈侵蝕面積占比呈下降趨勢外,其他強度等級的侵蝕面積占比均呈上升趨勢,其中:劇烈侵蝕面積占比一直呈下降趨勢,極強烈侵蝕面積占比在2000—2005年小幅上升后保持下降趨勢,強烈侵蝕面積占比在2000—2005年小幅上升后較為穩(wěn)定,中度、輕度、微度侵蝕面積占比在研究時段整體上升。 關中2000—2015年除輕度和微度侵蝕面積占比上升外其他強度等級侵蝕面積占比均呈下降趨勢,其中:劇烈侵蝕面積占比在2000—2010年持續(xù)上升而后在2010—2015年迅速下降,極強烈和強烈侵蝕面積占比在2000—2005年有所上升之后呈下降趨勢,中度和輕度侵蝕面積占比在3個時段呈上升、下降、上升的波動狀態(tài),微度侵蝕面積占比在2000—2005年有所下降后一直保持上升趨勢。陜南劇烈和極強烈侵蝕面積占比2000—2010年持續(xù)上升、2010—2015年下降,微度侵蝕面積占比呈現(xiàn)先降后升又降的波動變化,其他強度等級侵蝕面積占比變化情況與關中的類似。

圖4 陜西省各區(qū)域各級土壤侵蝕面積占比變化情況
在土壤侵蝕各因子中,土壤可蝕性K、地形因子LS發(fā)生變化所需的時間較長,可以認為在研究時段內(nèi)相對穩(wěn)定,因此本研究主要依據(jù)降雨侵蝕力R、植被覆蓋因子C、水土保持措施因子P的時間變化率來識別土壤侵蝕強度時間變化的驅動因子(見表2)。 由表2可知:2000—2005年,R年均變化率為2.70%(表明該時段降雨侵蝕力增大),C年均變化率為-2.50%(表明該時段植被覆蓋度提高),P年均變化率僅為0.08%,R與C年均變化率絕對值的差別小于10%,因此R和C共同驅動了該時段土壤侵蝕強度的變化;2005—2010年,R大幅增大、C大幅減小、P沒有變化,R年均變化率的絕對值遠大于C的,因此R為該時段土壤侵蝕強度變化的主導驅動因子;2010—2015年,R大幅減小、C大 幅 增 大、P變 化 較 小,年 均 變 化 率 分 別 為-11.57%、2.62%、-0.10%,R仍是土壤侵蝕強度變化的主要驅動因子。 從2000—2015年長時段整體看,長時間序列平抑了R的劇烈波動,C成為土壤侵蝕強度變化的主要驅動因子。 上述情況表明,采用最大變化率法對土壤侵蝕強度時間變化驅動力進行識別時對降雨侵蝕力極端值較為敏感,應盡量擴大時間步長以消除其影響。

表2 土壤侵蝕各因子的年均變化率 %
為了消除土壤侵蝕因子年際波動對空間變化驅動因子識別的影響,把各市及全省各因子2000年、2005年、2010年、2015年這4 個年份的平均值作為基礎數(shù)據(jù),分析土壤侵蝕強度空間變化的主要驅動因子(其結果可能與最新年份存在差異)。 依據(jù)式(10)計算的各市土壤侵蝕不同因子空間變化率(即相對于全省平均值的變化率)SX見表3,各市空間變化率絕對值最大的因子為空間變化(即各市相對于全省均值的差異)主要驅動因子,若變化率為負則表示該因子是土壤侵蝕強度降低的因子(即負向貢獻因子),反之若變化率為正則表示該因子是土壤侵蝕強度提高的因子(即正向貢獻因子)。
由表3可知:陜北榆林市土壤侵蝕強度空間變化的主要驅動因子是C(C空間變化率為正且極高,表明相對于全省來說榆林市植被覆蓋度較低);延安市R和LS對土壤侵蝕強度變化貢獻較大且均為負向貢獻,其中R為主要驅動因子。 關中西安、銅川、寶雞三市主要驅動因子為C,而渭南、咸陽兩市主要驅動因子為LS(可能與這兩市地形較其他市平坦有關);渭南、咸陽、西安三市水土保持措施因子P有較大負向貢獻,原因可能是這三市建設用地和耕地面積占比較大。 陜南三市R、LS、P為正向貢獻因子,C為負向貢獻因子且為主要驅動因子,商洛和漢中市K為負向貢獻因子而安康市K為正向貢獻因子。

表3 各地市土壤侵蝕因子空間變化率 %
綜上所述,陜西省絕大多數(shù)地市C和LS空間變化率絕對值較大,因此這2 個因子是土壤侵蝕強度空間變化的主要驅動因子。 因子C對土壤侵蝕強度變化的貢獻,在陜北為正而在關中和陜南為負,原因可能是陜北地區(qū)植被覆蓋度明顯比關中和陜南地區(qū)的低。因子LS對土壤侵蝕強度變化的貢獻在陜南三市及關中西安和寶雞兩市為正,而在其他市為負,原因可能是陜南和關中部分區(qū)域山地較多。
上述分析表明,2000—2015年陜西省土壤侵蝕強度整體上降低,這與有關學者[3,39-42]的研究結果一致。從分區(qū)來看,高強度的土壤侵蝕面積占比在陜北一直呈下降趨勢,而在關中和陜南2000—2010年持續(xù)上升、2010—2015年趨于下降,其原因:首先,陜北黃土高原作為全省及國家的水土流失重點治理區(qū)一直受到重視,而對關中和陜南水土流失治理的重視程度相對不足;其次,可能受有關水土流失治理新政策的影響,例如《丹江口庫區(qū)及上游水污染防治和水土保持規(guī)劃》于2006年經(jīng)國務院批復后實施,使丹江口上游地區(qū)水土流失得到有效治理,因此2010—2015年陜南高強度土壤侵蝕面積占比降低。
在陜西省土壤侵蝕強度時間變化的驅動因子方面,2000—2015年C為主要驅動因子,這與王娟等[43]的研究結果相似;2000—2005年R和C為共同驅動因子,2005—2015年R為主要驅動因子,這與Wang等[41]的研究結果相似。 在三江源地區(qū),1997—2012年R為主要驅動因子[8],與本研究所得結果存在差異,原因可能是三江源地區(qū)人為干擾極少,與降雨相比其植被變化不明顯,而1999年開始實施的大規(guī)模退耕還林還草使陜西省植被覆蓋率明顯提高,在較長時間尺度下極端降雨(因子R)對土壤侵蝕的影響被平抑后,C成為主導因子。 土壤侵蝕強度空間變化的主要驅動因子為C和LS,這與鄒雅婧等[44]的研究結果相似。
本研究存在如下不足:受限于所收集數(shù)據(jù)的精度,土壤侵蝕部分因子的計算結果可能不夠準確;受限于所收集數(shù)據(jù)序列的長度,僅分析到2015年,沒有探究更長時間尺度和更近期的土壤侵蝕強度變化情況及其驅動因子;只考慮了水力侵蝕,對其他主要侵蝕類型如水蝕風蝕共同分布的情況則沒有考慮。 此外,各土壤侵蝕等級相互轉化的歸因非常復雜,本研究對其關注不夠。
(1)陜西省土壤侵蝕分布異質(zhì)性較強,劇烈和極強烈侵蝕主要分布在陜北尤其集中分布在榆林市東南部和延安市北部,在關中北部區(qū)域有少量分布,在陜南商洛市西南部和安康市東部有少量分布。
(2)2000—2015年陜西省土壤侵蝕強度總體呈現(xiàn)下降趨勢,其中陜北侵蝕強度一直呈下降趨勢、關中和陜南呈現(xiàn)先升后降的變化趨勢。
(3)陜西省土壤侵蝕強度時間變化的主要驅動因子,2000—2005年為降雨侵蝕力和植被覆蓋,2005—2015年為降雨侵蝕力,2000—2015年長時段總體上為植被覆蓋。
(4)陜西省土壤侵蝕強度空間變化的主要驅動因子為植被覆蓋和地形,其中榆林、銅川、寶雞、西安、商洛、安康和漢中七市的主要驅動因子為植被覆蓋,渭南、咸陽兩市的主要驅動因子為地形,延安市的主要驅動因子為降雨侵蝕力。
(5)最大變化率法可用于識別土壤侵蝕變化的主要時空驅動因子,但在識別主要時間驅動因子時對降雨侵蝕力極端值較為敏感,應盡量擴大時間步長以消除其影響。