植被覆蓋和降雨因子變化及對東北黑土區土壤侵蝕的影響

王 碩, 方海燕, 和繼軍

(1.首都師范大學 城市環境過程和數字模擬國家重點試驗室培育基地, 北京100048; 2.中國科學院地理科學與資源研究所 陸地水循環及地表過程重點試驗室, 北京 100101; 3.中國科學院大學 資源與環境學院， 北京 100049)

土壤侵蝕是世界性的環境問題之一[1]，其不僅破壞土地資源，降低土壤肥力，還會誘發洪澇災害，嚴重影響生態環境和可持續發展[2]。因此，探索土壤侵蝕變化的原因，成為控制水土流失的關鍵。20世紀Wischmeier[3-4]建立了通用土壤流失方程(USLE)，該方程將諸多影響土壤侵蝕速率的物理和管理因素歸納為6個主要因子，分別為降雨侵蝕力因子R，土壤可蝕性因子K，坡度坡長因子LS，植被覆蓋與管理因子C和水土保持措施因子P。近些年來，全球氣候變化開始對土壤侵蝕產生影響，同時大規模的人為活動在加重土壤侵蝕的產生[5]。國內外都開展了這方面的工作，例如Teng[6]和Chakrabortty[7]曾在青藏高原和印度東部研究氣候變化影響下的土壤侵蝕，吳昌廣[8]和余新曉[9]也曾結合降雨侵蝕力與植被覆蓋因子在三峽與黃土區進行土壤侵蝕研究，最終表明植被覆蓋和降雨侵蝕力因子已成為影響土壤侵蝕變化的兩個關鍵因素。

東北黑土區是中國重要的商品糧基地[10]，素有“北大荒”之稱。然而由于黑土區獨特的地理環境[11]，不合理的開墾耕作[12]以及近些年氣候變暖的影響，使得黑土區內植被覆蓋率和降雨發生了巨大變化[13]，最終導致其水土流失現象日益突出。許多研究者[14-16]曾采用RUSLE模型，在東北黑土區的諸多區域開展過工作。然而以往開展的工作中，單因子變化對侵蝕造成影響的研究相對較少，特別是基于植被覆蓋和降雨因子變化對整個黑土區侵蝕影響的研究更是未見報道。近20 a，該區域植被覆蓋率和降雨量變化較大，因此有必要針對這些變化的影響程度進行探討，有利于揭示東北黑土區土壤侵蝕的空間變化特征，對于深化該區的土壤侵蝕研究更具深遠意義。因此，本文利用氣象、遙感等資料，采用RUSLE模型估算土壤侵蝕，并將東北黑土區劃分為6個區域，從時空角度結合降雨侵蝕力因子、植被覆蓋與管理因子變化來揭示土壤侵蝕空間差異及其對降雨侵蝕力和植被覆蓋變化的敏感性，進而為水土保持規劃提供數據支持。

1 研究區與研究方法

1.1 研究區概況

東北黑土區(38°42′—53°36′，115°24′—135°12′)包括黑龍江、吉林、遼寧及內蒙古東北部地區，總面積約為1.24×106km2，海拔在0～2 667 m之間。大興安嶺、小興安嶺和長白山分別坐落于研究區的西部、東北和東南，遼河平原、松嫩平原、三江平原三大平原位于三大山脈之間并且由南至北排開，遼河、嫩江、松花江流淌其中。研究區屬于溫帶季風氣候，四季分明。近40 a年均降雨量由東南部的900 mm衰減至西北部的300 mm左右，近70%降水來自夏季。研究區主要有黑土、黑鈣土、潮土、白漿土等土壤類型。山區以森林為主，平原地區以耕地為主。為了方便研究區進行描述，本文依據松嫩平原區域、嫩江、松花江及遼河等地勢與水系為邊界，將研究區分為6個區域，分別為遼東半島及長白山南部(Ⅰ區)、三江平原及長白山北部(Ⅱ區)、小興安嶺(Ⅲ區)、蒙古北部和大興安嶺北部(Ⅳ區)、遼河平原及大興安嶺中南部(Ⅴ區)及松嫩平原(Ⅵ區)。

1.2 數據來源

本文2000，2005，2010，2015和2018年0.5°×0.5°的日降雨數據來自國家氣象科學中心。土壤數據來自聯合國糧食與農業組織The Harmonized World Soil Database (HWSD)。數字高程影像(DEM； 90 m分辨率)數據來源于中國科學院地理空間數據云平臺。2000，2005，2010，2015和2018年歸一化植被指數(NDVI)與土地利用數據(100 m分辨率)均來自中國科學院資源環境科學數據中心。

1.3 研究方法

由于本研究區較大，RUSLE模型相比于CSLE模型而言其參數更易獲取，且該模型已經在其他大尺度地區成功開展過工作，C，P因子能有所借鑒。因此本文采用RUSLE模型計算土壤侵蝕強度，表達式為：

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

(1)

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

1.3.1 降雨侵蝕力R因子R值反映了降雨對土壤的潛在侵蝕能力。章文波[17]指出我國東北地區利用逐日降雨數據計算R值精度較高。因而，本研究采用日降雨數據計算R值：

(2)

式中：Mi為第i個半月時段的侵蝕力值〔(MJ·mm)/(hm2·h·a)〕;K為該半月時段內的天數；Dj為半月時段內第j天的侵蝕性日雨量，要求日雨量大于12 mm，否則以0計算;α，β為模型的待定系數。

(3)

α=21 568×β-7.189

(4)

式中：Pd12為≥12 mm的日平均雨量(mm)；Py12為≥12 mm的年平均雨量(mm)。通過以上方法計算R并對其進行插值[18]，最終得到R值分布。

1.3.2 土壤可蝕性K因子 張科利[19]表明修正后的EPIC模型在我國東北地區適用性較好。因此本研究通過Sharpley[20]和Williams[21]提出的EPIC模型方法，再結合張科利修正公式進行計算：

(5)

K=0.515 75×Kepic-0.013 83

(6)

式中：Kepic為修正前的K因子；CLA，SIL，SAN分別為土壤黏粒,粉粒和砂粒含量(%);C為有機碳含量(%); SN1=1-SAN/100。

1.3.3 坡度坡長LS因子LS反映地貌特征對土壤侵蝕的影響。考慮到我國東北地勢起伏，而RUSLE是美國農業部以美國緩坡農用地為對象建立的。因此，本研究采用Liu[22-23]提出基于我國陡坡的LS因子修正方法：

(7)

(8)

式中：λ為坡長(m);θ為坡度(°)。

1.3.4 植被覆蓋與管理C因子 由于研究區較大，植被覆蓋度變化于0%～76%，為避免出現斑塊化結果和準確獲取C值，經參考黑土區相關文獻[24-25]后采用蔡崇法[26]建立的C值模型計算C因子值。計算方法為：

(9)

(10)

式中：c為植被覆蓋度; NDVI為歸一化植被指數; NDVImax為研究區全植被覆蓋的NDVI值; NDVImin為全裸地或無植被覆蓋的NDVI值。

1.3.5 水土保持措施P因子 根據《黑土區水土流失綜合防治技術標準》，規定<3°坡耕地實施等高改壟；3°～5°坡耕地采取地埂植物帶；>5°修筑梯田；研究區內的坡耕地進行了大規模治理[28]。根據相關文獻[29]賦值方法，最終將水域區P值取0，梯田P值為0.03，等高改壟與地埂植物帶的旱田P值取0.352，順耕耕作的農田及其他自然植被區P值取1。

1.4 敏感系數

本文采用敏感系數法來衡量侵蝕模數對因子變化的敏感程度。計算公式為：

(11)

式中：Q為敏感度,表示侵蝕模數在影響因子變動下的變化快慢,數值越大敏感性越強; ΔA為侵蝕模數年均變化量; ΔX為影響因子的年均變化量。

(12)

(13)

式中：A1為研究期第一年的侵蝕模數;A2為研究期最后一年的侵蝕模數;X1為研究期第一年的平均因子數值;X2為研究期最后一年的平均因子數值;T為研究期時間跨度(a)。

為了消除數量級的影響，將敏感度進行標準化處理：

(14)

式中：M為標準化敏感度，即敏感系數；Qmax為區域內敏感度最大值；Qmin區域敏感度最小值。敏感系數分級[30]詳見表1。

表1 土壤侵蝕敏感系數的分級

1.5 數據處理

本文借助ArcGIS 10.2軟件，利用2000，2005，2010，2015和2018年5期C因子與其他各因子的多年平均值進行運算，得到C因子變化下土壤侵蝕模數；同理，利用2000，2005，2010，2015和2018年5期R因子數據得到R因子變化下的土壤侵蝕模數；再結合5期對應年份的各個因子數據，得到研究區內各個年份的土壤侵蝕模數。最后根據水利部《土壤侵蝕分類分級標準(SL190-2007)》，將土壤侵蝕模數進行重分類，得到C因子變化下土壤侵蝕模數時空分布、降雨侵蝕力因子變化下的侵蝕模數分布和各個年份的土壤侵蝕模數分布。

2 結果與分析

2.1 模擬結果的驗證

2.1.1C因子結果驗證 本文采用的C值模型結果與張憲奎[27]和張雪花[16]的C因子賦值標準基本一致，因此該C值模型在東北黑土區可適用(表2)。

表2 C因子結果驗證對比

2.1.2 侵蝕模數結果驗證 在以往的研究中，東北黑土區采用137Cs示蹤、徑流小區監測與模型估算等方法開展了侵蝕研究的相關工作。例如，方華軍[31]、李國強[32]與張克新[33]同樣采用137Cs示蹤法分別在松花江鎮、拜泉縣與遼東灣的得到了220，245與1 739～3 892 t/(km2·a)的結果，與本研究所得的187.62，237.46與2 568.57～4 323.36 t/(km2·a)的估算值相差不大；劉寶元[34]還在鶴山農場利用徑流小區監測得到坡耕地的侵蝕模數為845～1 157 t/(km2·a)，對比所用數據為鶴山縣的土壤侵蝕模數，因此小于徑流小區結果。還有盛美玲[35]利用WATEM/SEDAM模型在得到的2005年土壤侵蝕模數為351.2 t/(km2·a)，顧治家[36]利用CSLE模型得到的結果為677 t/(km2·a)，Fang[29]利用RUSLE模型在東北黑土區尺度下計算2010年的侵蝕模數為943.7 t/(km2·a)，本研究結果均與以上模型估算值相近(表3)。本研究中土壤侵蝕模數偏小的原因可能來自兩個方面，其一是當前所能獲取的NDVI與土地利用數據分辨率過低平滑了數值[38]，降低了土壤侵蝕模數；其二則是研究時間和區域、以及賦值精細度不同所造成。

表3 土壤侵蝕模數結果驗證對比

2.2 RUSLE各因子分布及變化特征

2.2.1 研究區R因子變化分布及變化特征 2000—2018年R因子呈現先增加后降低再增加的趨勢。空間上，R因子由西南向東北部減少，區域上呈現：Ⅰ>Ⅱ>Ⅵ>Ⅲ>Ⅴ>Ⅳ的分布特點(圖1)。從局部上看，Ⅰ區是R值最大區并且以丹東市為中心向四周降低。遼東半島東側與西側受到長白山脈阻隔造成降雨量差距明顯，其中2018年東西側R值相差近2 000 (MJ·mm)/(hm2·h·a)。

圖1 研究區各年份降雨侵蝕力因子分布

在Ⅱ，Ⅲ和Ⅵ區R因子由765.42 (MJ·mm)/(hm2·h·a)增至1 417.79 (MJ·mm)/(hm2·h·a)。研究區內最小R值出現在Ⅳ和Ⅴ區，該區R值穩定升高但變化幅度較小，呈現大興安嶺中部大于南北部的分布特點。

2.2.2 研究區K因子與LS因子變化分布及變化特征 研究區內K因子最高值分布在大興安嶺南部(Ⅴ)和遼河平原北部(Ⅵ)內；最低值則出現在大興安嶺北部(Ⅳ)、遼河平原西南部(Ⅴ)和松嫩平原(Ⅵ)。研究區內LS因子的高值主要分布在三大山脈，研究區內的三大平原均由低值構成。

2.2.3 研究區C因子變化分布及變化特征 2000—2018年C值降低了42%，呈現中部與西南部偏高，東部與北部偏低的特征(圖2)。C因子的變化主要受夏季平均氣溫影響[39]。數據表明[40]，2000—2011年夏季平均氣溫持續升高，而2011年后夏季平均氣溫開始穩定，這趨勢與C值變化曲線相吻合。C因子降低受政策因素的作用也不容忽視。隨著三北防護林、退耕還林等工程的大規模實施，東北全區植被覆蓋面積出現了穩固的增加[41]。從局部看，Ⅵ區受到政策因素影響C值呈下降趨勢，18 a間C值小于0.1的面積占比由27.29%增至69.62%，為植被覆蓋率增加最快的區域。Ⅰ，Ⅱ，Ⅲ和Ⅳ區的C值處于持續降低的態勢，該范圍曾經存在大量的高值區域(C>0.1)，至2010年高值區域面積占比由31.62%下降至12.23%，這也與夏季平均氣溫的逐年升高息息相關。Ⅴ區在研究期間C值低于0.1的面積占比總體增加了34.12%。該區的C值與海拔高度的氣溫呈現負相關[42]，在高海拔地區常年氣溫較低，因此造成高海拔地區C值并未發生明顯改變。

圖2 研究區各年份植被覆蓋與管理因子分布

2.2.4 研究區P因子變化分布及變化特征 2000—2018年平均P值變化率不足1%，且變化均出現于黑龍江的北部(Ⅳ)和大興安嶺南部(Ⅴ)。相比而言，C與R因子變化率都超過了50%。因此，本文僅考慮C因子與R因子對于東北黑土區土壤侵蝕模數變化的影響。

2.3 單因子變化下的土壤侵蝕時空特征

2.3.1 研究區C因子變化下土壤侵蝕特征 在C因子變化下，2000—2018年研究區土壤侵蝕模數持續降低，呈現南部侵蝕嚴重、北部侵蝕輕微的特征(圖3和表4)。局部上看，研究區內中部和北部主要以微輕度侵蝕為主。Ⅱ，Ⅲ和Ⅳ區侵蝕模數與C因子的年均變化量都較小，敏感系數均小于0.3，屬不敏感區域。Ⅵ區C因子年均變化量達0.021，高于全區的平均水平，該因子大幅度下降勢必導致該區侵蝕模數隨之大幅度下降。根據敏感度計算公式，Ⅵ區敏感度為0，屬不敏感區域，表明該區在單位植被覆蓋度增加的情況下，對侵蝕的緩解程度不明顯。以上不敏感區域大部分處于地勢平坦且土壤可蝕性較低的地區，所以低值的K與LS因子削弱了植被覆蓋對侵蝕變化的影響。

表4 研究區不同分區土壤侵蝕模數變化量與敏感系數

圖3 研究區C變化下的土壤侵蝕模數分布

土壤侵蝕最嚴重的區域為Ⅰ與Ⅴ區，研究期間強度以上侵蝕向中度以下侵蝕轉變的面積近63 000 km2。這兩個區域的C因子年均變化量僅為Ⅵ區的1/2，但其年均侵蝕變化量卻為Ⅵ區的3～4倍，分別達到了-214.33與-271.31 t/(km2·a)，同時Ⅰ與Ⅴ區的敏感系數分別為0.95和1.00，是研究區內的兩個強度敏感區域。Ⅰ區和Ⅴ區分別坐落于長白山中南部與大興安嶺中南部，地勢陡峭，LS值最大；并且Ⅰ區瀕臨渤海，降雨充沛，R值最大；而Ⅴ區是K，P值最高的區域。在其他因子值較高的情況下，C因子若發生輕微變化，勢必會造成侵蝕模數發生較大的改變。因此應優先提高這兩個區域的植被覆蓋程度，能夠更為高效的減少土壤侵蝕。

2.3.2R因子變化下土壤侵蝕時空特征 2000—2018年，受R因子變化的影響，土壤侵蝕模數呈現先增加后減少再增加的趨勢，空間上呈現西南部侵蝕嚴重，東北部侵蝕輕微的特點(圖4和表5)。局部來看，Ⅰ區年均侵蝕增加量雖為全區最大的69.81 t/(km2·a)，但該區的敏感系數僅為0.08，屬不敏感區域。Ⅱ，Ⅲ和Ⅵ區侵蝕模數變化小，敏感系數均在0.1以下，是不敏感區域。以上4個區域不敏感的原因在于C因子較小，較高的植被覆蓋可有效攔截降雨，減少降雨擊濺侵蝕，緩解強降雨對侵蝕造成的影響。在R因子偏低且浮動較小的Ⅳ和Ⅴ區，出現較嚴重的侵蝕和侵蝕強度變化，敏感系數分別為1.00和0.45，屬中度敏感和強度敏感。Ⅳ和Ⅴ區分布著大興安嶺與蒙古高原，地勢陡峭并且包含大量的K，LS，C和P因子的高值區域，便造成這兩個分區敏感。因此該區域應優先減少坡耕地，實行免耕等耕作方式，降低降雨對侵蝕所造成的影響。

圖4 研究區R變化下的土壤侵蝕模數分布

表5 研究區內區域變化量與敏感系數

2.4 不同年份土壤侵蝕時空特征

研究期間東北黑土區年均侵蝕模數為893.53 t/(km2·a)，遠超過該區的容許土壤流失量200 t/(km2·a)。Ⅰ和Ⅴ區侵蝕較為嚴重，其余地區屬于微度侵蝕和輕度侵蝕且變化較小。從土壤侵蝕面積來看，從2000—2015年研究區微度侵蝕面積占比增加了3.73%，而中度以上侵蝕的面積占比縮幅為3.84%，說明研究區出現中度以上侵蝕向輕微度侵蝕轉移的跡象。從2015—2018年，微度侵蝕面積縮小了55 989.67 km2，輕度侵蝕面積增長了50 462.60 km2，研究期后3 a出現微度向輕度以上侵蝕轉移的情況(圖5和表6)。

表6 研究區各年份不同土壤侵蝕強度分級面積 t/(km2·a)

圖5 研究區各年份的土壤侵蝕模數分布

從空間局部上看，研究期間Ⅰ區侵蝕模數呈現減弱趨勢。極強度以上侵蝕面積由9 914.18 km2減少到1 793.61 km2。該區侵蝕模數對C因子呈強度敏感，研究期C因子由0.53下降至0.17，使該區土壤侵蝕強度降低。但其地處長白山南部，地勢陡峭，瀕臨渤海，常年雨量大，也使其侵蝕模數屬于全區偏大的水平。

遼河平原(Ⅴ)是極強度以上侵蝕的主要分布區，侵蝕模數呈現先減后增的趨勢。該區侵蝕對C因子與R因子皆為強度敏感。在2000—2010年，該區C因子由0.55驟降至0.19，使得侵蝕得到極大改善。在2010—2015年，該區R值由2 634.83 (MJ·mm)/(hm2·h·a)降至813.98 (MJ·mm)/(hm2·h·a)，同樣改善了侵蝕情況。2015至2018年間，C因子出現30.44%的回升，導致該區侵蝕模數的少量增加。

大興安嶺(Ⅴ)地區屬于對R與C因子敏感地區，呈現先增后減再增的趨勢。2000年強度以上侵蝕面積為27 848.86 km2，后10 a侵蝕開始由大興安嶺南部蔓延至大興安嶺西南部，強度以上侵蝕面積達到34 475.68 km2，占該年強度以上侵蝕總面積的73.72%。前10 a中R因子增幅達42.85%，導致該區侵蝕嚴重化；2010—2015年，該區C與R值同時降低，改善了該區的侵蝕情況。2015—2018年，該區R和C值同時增加，強度以上侵蝕面積再次出現增加，達到了24 753.39 km2；這些年份盡管大興安嶺中部的侵蝕情況有所緩解，但大興安嶺南部與內蒙古高原接壤地帶仍是高侵蝕地帶，這與該區陡峭的地形有關。

Ⅰ區和Ⅴ區盡管侵蝕情況得到了相當大程度的緩解，但大部分區域仍然處于中度侵蝕的水平。而該地區屬于對C因子強度敏感且地勢陡峭的地區，因此應當制定相應政策減少坡耕地，增加植被覆蓋的程度，能夠更高效的減少當地的土壤侵蝕。Ⅳ區雖一直處于輕中度侵蝕的水平，但研究期間侵蝕模數出現了明顯的增加，僅此更應該盡早的制定農田免耕，梯田等水土保持措施來緩解降雨的影響。

3 結 論

(1) 整體上，東北黑土區平均土壤侵蝕量有所下降，其原因主要為C因子和R因子的變動，其表現是部分地區降雨量的增減、植被覆蓋的變化。

(2) 2000—2018年C因子呈持續降低的趨勢；該因子呈現Ⅵ和Ⅴ區偏高，Ⅰ，Ⅱ，Ⅲ和Ⅳ偏小的分布特征。在該因子影響下，侵蝕模數從1 551.07 t/(km2·a)降低至665.39 t/(km2·a)。研究區Ⅰ和Ⅴ的敏感系數分別為0.95和1.00，是強度敏感的地區；Ⅱ，Ⅲ，Ⅳ和Ⅵ區敏感系數最低，屬于不敏感地區。

(3) 研究期間R因子呈現先增加后降低再增加的趨勢，空間上呈現：區域Ⅰ>Ⅱ>Ⅵ>Ⅲ>Ⅴ>Ⅳ特征。僅在R因子影響下，侵蝕模數呈先增加后降低再增加的趨勢，Ⅴ與Ⅳ區敏感系數分別為1.00和0.45，分別屬于強度敏感和中度敏感地區；Ⅰ，Ⅱ，Ⅲ和Ⅵ區的敏感系數均較低，屬于不敏感地區。

(4) 多年來，東北黑土區區土壤侵蝕主要以微度和輕度侵蝕為主，兩類侵蝕面積占比在92.3%～96.2%之間波動。強度侵蝕以上的地區基本分布在Ⅰ和Ⅴ區，侵蝕強度從西南到東北逐漸遞減。研究期侵蝕模數從1 175.20 t/(km2·a)變化到822.07 t/(km2·a)，且2015年還出現過628.73 t/(km2·a)的最低值。