郭 晴，任永康，陳 濤*

（1. 中國地質(zhì)大學（武漢） 地球物理與空間信息學院，湖北 武漢 430074）

土壤侵蝕對土壤分級、農(nóng)業(yè)生產(chǎn)、水文系統(tǒng)和環(huán)境的不利影響，長期以來被認為是影響人類可持續(xù)性發(fā)展的嚴重問題[1]。截至2019 年，我國土壤侵蝕總面積已達到271.08 萬km2，成為世界上土壤侵蝕最嚴峻的國家之一。其中，水力侵蝕面積為113.47 萬km2，占總面積的41.86%；風力侵蝕面積為157.61 萬km2，占總面積的58.14%[2]。土壤侵蝕會引起土壤肥力衰退，進而威脅農(nóng)業(yè)生產(chǎn)安全。徑流泥沙搬移的污染源將影響侵蝕鄰近區(qū)生態(tài)和經(jīng)濟的發(fā)展，且侵蝕泥沙的搬運使得土壤碳、氮、磷等元素的含量與組分產(chǎn)生變化，從而影響全球生源要素的循環(huán)，甚至成為全球氣候變化的驅(qū)動要素之一[3]。

土壤侵蝕研究通常采用實地測算和構(gòu)建土壤侵蝕模型來開展。由于實地觀測難以直接監(jiān)測空間大范圍土壤侵蝕情況，且易受主觀因素和監(jiān)測技術(shù)的限制，因此不能滿足土壤侵蝕治理工作對時效性、快速性與準確性的需求[4]；而土壤侵蝕模型不僅可以解決實地觀測中的矛盾，還保障了土壤侵蝕預測和評價的準確性。根據(jù)研究方法，土壤侵蝕模型可分為物理模型和經(jīng)驗統(tǒng)計模型[5]，物理模型中最具代表性的是水蝕預報模型[6]，雖對土壤侵蝕的描述十分全面，但所需數(shù)據(jù)量大且難以在大區(qū)域中實現(xiàn)[7]；經(jīng)驗統(tǒng)計模型中認可度最高、應用范圍最廣的是Kenneth G R[8]等提出的土壤流失方程（USLE）和修正后的通用土壤流失方程（RUSLE）[6]。隨著土壤侵蝕研究的不斷深入，許多學者結(jié)合研究區(qū)的真實狀況對模型中各因子進行了修正，并與GIS 和遙感技術(shù)進行集成，這樣不僅具有了RUSLE 模型結(jié)構(gòu)簡單、所需數(shù)據(jù)量少的特點[9]，而且具有了GIS 和遙感技術(shù)快速獲取、有效管理并實時分析時空大數(shù)據(jù)的優(yōu)勢，為土壤侵蝕提供了科學、快速、準確的監(jiān)測和評估。例如，劉寶元[10]等根據(jù)我國國情，結(jié)合GIS 和遙感技術(shù)，輔以遙感圖像，修正了RUSLE 中的各個因子；楊佳佳[11]等基于遙感和GIS 技術(shù)，利用修正后的RUSLE模型對黑龍江典型黑土區(qū)土壤侵蝕進行了監(jiān)測；王丹媛[12]等利用遙感數(shù)據(jù)修改了桂南沿海諸河流域USLE 模型中的因子，并評價了土壤侵蝕強度。近年來，眾多學者將RUSLE模型與GIS相結(jié)合以研究大尺度土壤侵蝕和水土保持格局與規(guī)律[13-14]，如陳峰[15]等利用RUSLE 模型研究了滇南山區(qū)土壤侵蝕的時空演變情況；王志杰[16]等運用RUSLE模型探究10 a間貴陽市土壤侵蝕空間演變規(guī)律。

作為北京市最大的飲用水源供應地，密云水庫擔負著北京生產(chǎn)生活用水的重要任務(wù)。然而，密云水庫流域生態(tài)環(huán)境較脆弱，土壤侵蝕不利于流域內(nèi)的水循環(huán)，進而將威脅北京市民的飲用水安全[17]。因此，研究密云水庫的土壤侵蝕情況，為改進水庫周邊水土保持工作提供科學性的指導，對于保護其生態(tài)環(huán)境具有重要意義。目前，針對密云水庫流域土壤侵蝕的研究主要為單一年份，缺乏時空動態(tài)變化的分析，在一定程度上不利于決策和管理，阻礙了流域的可持續(xù)發(fā)展。本文基于遙感和GIS技術(shù)，利用RUSLE模型對密云水庫流域2001—2020年土壤侵蝕情況進行了動態(tài)監(jiān)測，分析了時空動態(tài)變化特征及其影響因素，為密云水庫流域的水土保持提供了科學合理的依據(jù)。

1 研究區(qū)概況

密云水庫流域位于北京市城北 13 km2（41°14′~41°05′N、116°07′~117°30′E），面積為180 km2， 包括北京市的密云、懷柔、延慶以及河北省的灤平、豐寧、赤城等10個區(qū)縣（圖1）；坐落在潮河以及白河中游偏下，流域面積為15 788 km2。水庫流域年均氣溫為9～10 ℃，年均降水量為488.9 mm，主要集中在汛期。流域內(nèi)土壤主要包括褐土、棕壤、草甸土和粟鈣土，植被以闊葉混交雜木林、油松、刺槐和落葉松為主。

圖1 密云水庫概要圖

2 數(shù)據(jù)來源與研究方法

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

本文采用的遙感數(shù)據(jù)包括2001 年、2005 年、2009 年 、 2011 年 、 2015 年 和 2020 年 6 期 共 24 景Landsat8/TM 遙感影像。由于部分年份遙感影像存在云量較多、噪聲明顯和數(shù)據(jù)條帶丟失等問題，因此在等間隔的基礎(chǔ)上從研究年份中選擇數(shù)據(jù)精度較高的6 期數(shù)據(jù)作為研究對象；由于研究區(qū)遙感影像數(shù)據(jù)具有多時相、多光譜、多傳感器等特點，因此利用ENVI 5.3對原始影像進行輻射定標、大氣校正等預處理，形成標準化的多源遙感數(shù)據(jù)集。本文采用的地形數(shù)據(jù)為ASTER GDEM 30M產(chǎn)品；降水量數(shù)據(jù)包括2001—2020年流域內(nèi)豐寧、赤城、沽源、崇禮、承德、懷來、灤平、懷柔、密云、興隆、延慶和張家口地區(qū)的年均降水量；土地利用數(shù)據(jù)包括2001年、2005年、2009年、2011 年、2015 年和 2020 年 6 期 1 km 分辨率的中國土地利用數(shù)據(jù)（包含耕地、林地、草地、水域、建設(shè)用地和未利用地6個一級地類）；土壤數(shù)據(jù)為土壤科學數(shù)據(jù)庫中的中國第二次普查數(shù)據(jù)，包括土壤質(zhì)地砂礫、粉粒和有機碳含量等數(shù)據(jù)；各區(qū)縣生產(chǎn)總值來源于當?shù)卣俜骄W(wǎng)站數(shù)據(jù)；人口數(shù)據(jù)來源于中國統(tǒng)計信息網(wǎng)（表1）。

表1 數(shù)據(jù)列表

2.2 研究方法

基于GIS和遙感技術(shù)，本文利用RUSLE模型計算得到土壤侵蝕模數(shù)，并利用轉(zhuǎn)換矩陣分析了2001—2020年土壤侵蝕強度的空間動態(tài)變化情況。RUSLE模型的計算公式為：

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

1）降雨侵蝕力因子表示降雨引起土壤侵蝕的潛在能力。本文基于監(jiān)測區(qū)域代表站點連續(xù)30 min的最大降雨強度和每次降雨的降雨量，利用基于日降雨量的擬合模型[18-19]計算降雨侵蝕因子。為了得到流域內(nèi)降雨侵蝕力因子的連續(xù)數(shù)據(jù)，在各測站降雨侵蝕力值和空間分布的基礎(chǔ)上，利用GIS 的樣條插值實現(xiàn)。其計算公式為：

式中，Pf為汛期雨量；I30B為侵蝕性降雨最大30 min雨強的年代表值，單位為mm；i 為侵蝕性降雨次數(shù)；Pt為汛期月份的月降雨總量。

2）土壤可蝕性因子是指土壤對侵蝕的敏感性，與降雨、土壤自身的理化特性等相關(guān)。采用Wisch?meier W H[20]等提出的公式，根據(jù)土壤質(zhì)地、土壤中的有機質(zhì)含量等主要土壤性質(zhì)數(shù)值指標估算土壤可蝕性因子。其計算公式為：

式中， M 為優(yōu)勢粒徑組成的乘積；a 為有機質(zhì)含量百分比；b 為土壤結(jié)構(gòu)等級；c 為土壤滲透等級。

3）坡度坡長因子。坡長因子是指任意坡長的單位面積土壤流失量與標準坡長單位面積土壤流失量之比；坡度因子是指任意坡度下的單位面積土壤流失量與標準坡度下單位面積土壤流失量之比。由于坡度從20%增加到40%和60%時，坡長指數(shù)無變化[6]，因此當坡度<21%時，采用式（5），當坡度＞21%時，采用式（6）[21]。

4）植被覆蓋與管理因子是指在土壤、坡度和降雨條件相同的情況下，具有特定植被覆蓋的土地的土壤流失量與裸地的土壤流失量的比值。本文采用蔡崇法[22]等提出的利用植被覆蓋度與定量估算管理因子的回歸方程，即

式中， f 為植被覆蓋度；C 為定量估算管理因子

植被覆蓋度計算采用李苗苗[23]等提出的改進的像元二分模型，則有：

式中， NDVI 為歸一化植被指數(shù)值，無量綱；NDVImax、NDVImin為遙感圖像中給定置信度置信區(qū)間內(nèi)的最大值與最小值，本文采用的置信區(qū)間為[5%，95%]。

5）水土保持措施因子是指采取一定的水保措施后的土壤流失量與標準狀況下土壤流失量的比值。水土保持措施因子目前主要通過布設(shè)天然小區(qū)試驗得到。由于難以獲得水土保持措施因子，在考慮財力物力人力等綜合因素后，結(jié)合前人的研究結(jié)果[24-25]，本文將水土保持措施因子設(shè)定為1。

2.3 轉(zhuǎn)移矩陣

土壤侵蝕轉(zhuǎn)移矩陣反映了研究區(qū)在某一時段期初和期末不同侵蝕等級面積之間的轉(zhuǎn)化信息，包括期初各類面積的轉(zhuǎn)出量和期末各類面積的轉(zhuǎn)入量[26]。本文計算得到2001 年和2020 年不同侵蝕等級土壤侵蝕模數(shù)的轉(zhuǎn)移矩陣，并探究了侵蝕等級的變化特征。

2.4 相關(guān)性分析

相關(guān)性分析能反映因變量與自變量之間相互關(guān)系的密切程度。本文將土壤侵蝕模數(shù)分別與年均降雨量與平均植被覆蓋度進行統(tǒng)計回歸分析[27]，以探究年均降雨量、平均植被覆蓋度與土壤侵蝕模數(shù)的相關(guān)性。

3 研究結(jié)果分析

3.1 土壤侵蝕模數(shù)的時間變化特征

密云水庫流域年均土壤侵蝕模數(shù)呈“快速減小、小幅增長、又迅速回落”的趨勢，20 a間研究區(qū)土壤侵蝕狀況整體得到明顯控制，2009—2015年土壤侵蝕呈現(xiàn)惡化趨勢。根據(jù)國家水利部發(fā)布的SL190-2007《土壤侵蝕分類分級標準》[28]，本文將研究區(qū)土壤侵蝕分為微度、輕度、中度、強烈、極強烈和劇烈侵蝕6個等級，并統(tǒng)計分析不同時期密云水庫流域土壤侵蝕強度的面積比例。20 a間研究區(qū)土壤侵蝕以輕度和微度侵蝕為主，二者占總面積的70%以上，且呈波浪式變化；其次為中度侵蝕；除強烈侵蝕和極強烈侵蝕面積以外，其他侵蝕強度面積呈總體下降趨勢。2001—2020年各侵蝕強度等級的年變化率（降幅）均在5%~10%，但2001—2020年強烈和極強烈侵蝕增幅明顯（圖2）。

圖2 2001—2020年不同等級平均土壤侵蝕模數(shù)折線圖

由面積轉(zhuǎn)移矩陣可知（表2）， 2001—2020 年研究區(qū)土壤侵蝕以微度侵蝕少量轉(zhuǎn)出、其他侵蝕等級不同程度轉(zhuǎn)入為主要特征；侵蝕強度等級的轉(zhuǎn)出方向以向低等級侵蝕轉(zhuǎn)化為主要特征，微度侵蝕轉(zhuǎn)入了14.69 km2，占總微度侵蝕轉(zhuǎn)出面積的67.35%，輕度侵蝕轉(zhuǎn)入微度侵蝕的面積為2.91 km2，占輕度侵蝕轉(zhuǎn)出面積的56.7%，中度侵蝕轉(zhuǎn)入微度和輕度侵蝕的面積為0.43 km2，占總中度侵蝕轉(zhuǎn)出面積的67.2%，劇烈侵蝕幾乎全部轉(zhuǎn)向低等級侵蝕；但極強烈侵蝕大部分面積轉(zhuǎn)入為更高等級，強烈侵蝕轉(zhuǎn)向低等級侵蝕的面積占比為44.7%，因此流域內(nèi)20 a 間土壤侵蝕呈“整體好轉(zhuǎn)、局部惡化”的趨勢。

表2 2001—2020年研究區(qū)土壤侵蝕強度面積轉(zhuǎn)移矩陣/km2

3.2 土壤侵蝕模數(shù)的空間變化特征

3.2.1 研究區(qū)不同行政區(qū)內(nèi)土壤侵蝕變化分析

本文疊置分析了密云水庫流域研究期土壤侵蝕強度分布圖與研究區(qū)行政區(qū)劃圖，獲得流域內(nèi)研究期間各縣區(qū)土壤侵蝕強度的面積比例（表3），可以看出，流域內(nèi)土壤侵蝕最嚴重的區(qū)域為密云區(qū)和灤平縣，均值都在1 500 t/(km2·a)以上；密云、崇禮和灤平縣的年均侵蝕模數(shù)降幅最大，分別為53.9%、63.5%和42.5%；崇禮縣水土保持情況最好，平均土壤侵蝕模數(shù)接近 700 t/(km2·a)，且降幅較大。

表3 密云水庫流域2001—2020年各縣區(qū)侵蝕強度面積/km2

3.2.2 區(qū)域內(nèi)不同海拔下土壤侵蝕變化分析

本文疊置分析了土壤侵蝕分布圖與數(shù)字高程模型分級圖，得到密云水庫流域2001—2020年不同海拔下土壤侵蝕變化情況（表4、5），可以看出，隨著海拔的升高，研究區(qū)土壤侵蝕模數(shù)整體呈“高—低—高”的特點，海拔500 m 以下以及海拔1 500 m 以上的地區(qū)侵蝕相對嚴重。在低海拔地帶，由于地形平坦有利于開發(fā)，導致人口分布密集、城市化程度較高，進而對地表植被造成嚴重破壞；而在1 500~3 000 m高海拔地區(qū)，土地類型以高山裸土區(qū)和荒草地為主，地勢陡峻，高強度的降雨易產(chǎn)生滑坡、泥石流等自然災害。

表4 密云水庫流域2001—2009年不同海拔下土壤侵蝕變化分析

表5 密云水庫流域2011—2020年不同海拔下土壤侵蝕變化分析

3.2.3 區(qū)域內(nèi)不同土地利用類型下土壤侵蝕變化分析

本文統(tǒng)計分析了不同土地利用類型的侵蝕狀況，結(jié)果如表6所示。各種土地利用類型中，耕地的平均土壤侵蝕強度普遍較高，而城鄉(xiāng)居民用地和水域的平均土壤侵蝕強度較低。2001年、2009年耕地和林地的侵蝕強度接近于未利用地，處于較高水平；2015年耕地、林地和草地的侵蝕強度均大于1 500 t/（km2·a），林地的侵蝕強度最大，未利用地的侵蝕強度未達到1 500 t/（km2·a）；2020年耕地、林地和草地的侵蝕強度均有所回落，與逐步治理耕地以及大力推行退耕還林政策密不可分。

表6 密云水庫流域2001—2020年不同土地利用類型下土壤侵蝕變化/（t/（km2·a））

3.3 影響因素分析

3.3.1 自然因素

本文對密云水庫流域研究期間的平均土壤侵蝕模數(shù)與平均植被覆蓋變化、平均降雨侵蝕模數(shù)進行回歸分析，結(jié)果表明，平均土壤侵蝕模數(shù)與平均植被覆蓋變化存在明顯的負相關(guān)性，相關(guān)系數(shù)為-0.77（圖3）；而與平均降雨侵蝕模數(shù)呈正相關(guān)性，相關(guān)系數(shù)為0.71（圖4）。由于平均土壤侵蝕模數(shù)與植被覆蓋度的決定系數(shù)大于其與降雨侵蝕模數(shù)的決定系數(shù)，因此平均土壤侵蝕模數(shù)與植被覆蓋度的相關(guān)程度更高。

圖3 植被覆蓋度與土壤侵蝕模數(shù)關(guān)系散點圖

圖4 平均降雨侵蝕模數(shù)和土壤侵蝕模數(shù)關(guān)系散點圖

3.3.2 人為因素

分析研究區(qū)各縣市人口數(shù)據(jù)（圖5）和土壤侵蝕強度的面積比例統(tǒng)計表可知，研究期內(nèi)密云區(qū)、興隆縣、灤平縣和懷柔區(qū)的平均侵蝕面積均在1 000 km2以上，且人口逐年遞增，其中密云區(qū)和懷柔區(qū)的人口基數(shù)大、增速也較快，人口的迅速增長促使農(nóng)民大量開荒種糧，導致土地墾荒率越來越高，加劇了水土流失；后期（2015—2020年）崇禮縣、沽源縣和延慶區(qū)等地區(qū)人口出現(xiàn)零增長甚至負增長的情況，人口對環(huán)境的壓力下降了，水土流失也有所緩解。

圖5 研究區(qū)各縣區(qū)常駐人口變化圖

分析各區(qū)縣GDP 變化情況（圖6）和土壤侵蝕強度的面積比例統(tǒng)計表可知，土壤侵蝕與社會經(jīng)濟發(fā)展具有相關(guān)性，研究區(qū)初期和中期（2001—2011年）各行政區(qū)的GDP增速較大，如懷柔區(qū)、密云區(qū)和灤平縣，其土壤侵蝕面積相對較大，原因在于其地理位置離北京主城區(qū)較近，經(jīng)濟社會活動較頻繁，經(jīng)濟的發(fā)展伴隨著開發(fā)區(qū)建設(shè)等土建工程的急劇增加，使得大量草地、林地遭到了破環(huán)，加劇了水土流失；后期各行政區(qū)GDP增速逐漸放緩，政府踐行綠色發(fā)展理念，加大對水土保持工作的投入，使得土壤侵蝕面積較前期有所回落。

圖6 研究區(qū)各縣區(qū)GDP變化圖

4 結(jié) 語

本文利用RUSLE 模型對密云水庫流域2001—2020 年的土壤侵蝕情況進行了定量研究。結(jié)果表明，2001—2020年密云水庫流域土壤侵蝕模數(shù)大致經(jīng)歷了先下降后回升又大幅下降的變化過程，最小值出現(xiàn)在2009 年，2020 年比2001 年的土壤侵蝕面積有所減少，說明整體情況在好轉(zhuǎn)；流域內(nèi)侵蝕等級主要由輕度轉(zhuǎn)為微度，但強烈侵蝕和劇烈侵蝕主要集中在灤平縣和密云區(qū)；隨著海拔的上升，流域內(nèi)侵蝕模數(shù)先升高后降低又升高；土壤侵蝕模數(shù)與土地利用類型相關(guān)性程度較高；植被覆蓋度升高，土壤侵蝕模數(shù)下降，降雨侵蝕加劇，土壤侵蝕模數(shù)降低，但前者的影響程度較大；人口與社會經(jīng)濟的發(fā)展也會對土壤侵蝕模數(shù)造成影響。