郭 晴，任永康，陳 濤*

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

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

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

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

1 研究區(qū)概況

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

圖1 密云水庫(kù)概要圖

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

2.1 數(shù)據(jù)來(lái)源與預(yù)處理

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

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

2.2 研究方法

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

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

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

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

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

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

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

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

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

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

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

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

2.3 轉(zhuǎn)移矩陣

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

2.4 相關(guān)性分析

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

3 研究結(jié)果分析

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

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

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

由面積轉(zhuǎn)移矩陣可知（表2）， 2001—2020 年研究區(qū)土壤侵蝕以微度侵蝕少量轉(zhuǎn)出、其他侵蝕等級(jí)不同程度轉(zhuǎn)入為主要特征；侵蝕強(qiáng)度等級(jí)的轉(zhuǎn)出方向以向低等級(jí)侵蝕轉(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)向低等級(jí)侵蝕；但極強(qiáng)烈侵蝕大部分面積轉(zhuǎn)入為更高等級(jí)，強(qiáng)烈侵蝕轉(zhuǎn)向低等級(jí)侵蝕的面積占比為44.7%，因此流域內(nèi)20 a 間土壤侵蝕呈“整體好轉(zhuǎn)、局部惡化”的趨勢(shì)。

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

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

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

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

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

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

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

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

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

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

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

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

3.3 影響因素分析

3.3.1 自然因素

本文對(duì)密云水庫(kù)流域研究期間的平均土壤侵蝕模數(shù)與平均植被覆蓋變化、平均降雨侵蝕模數(shù)進(jìn)行回歸分析，結(jié)果表明，平均土壤侵蝕模數(shù)與平均植被覆蓋變化存在明顯的負(fù)相關(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)系散點(diǎn)圖

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

3.3.2 人為因素

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

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

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

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

4 結(jié) 語(yǔ)

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