1980—2020年延河甘谷驛流域土壤侵蝕評(píng)價(jià)與驅(qū)動(dòng)因子分析

陳方磊 王計(jì)平 程復(fù) 謝海燕

摘要：采用日降雨量、DEM、土壤類型、泥沙含量及多期NDVI等數(shù)據(jù)，基于修正通用土壤流失方程（RUSLE）和地理探測(cè)器，研究了國(guó)家生態(tài)退耕還林還草工程實(shí)施前后近41年延河甘谷驛流域土壤侵蝕動(dòng)態(tài)與驅(qū)動(dòng)因子。結(jié)果表明，1980—2020年研究區(qū)土壤侵蝕強(qiáng)度總體呈波動(dòng)變化趨勢(shì)，1980年、1990年、2000年、2010年和2020年平均侵蝕模數(shù)分別為6 746.30、5 740.28、6 389.56、5 450.46、5 480.56 t/（km2·年）。1980—2000年研究區(qū)整體侵蝕強(qiáng)度逐漸增強(qiáng)，強(qiáng)烈及以上等級(jí)侵蝕面積占比逐漸增加，表現(xiàn)為“增蝕升級(jí)”的特點(diǎn)；2000年后研究區(qū)內(nèi)土壤侵蝕強(qiáng)度開始降低，強(qiáng)烈及以上等級(jí)的侵蝕面積減少，總體表現(xiàn)為“減蝕降級(jí)”的特點(diǎn)。研究區(qū)土壤侵蝕強(qiáng)度隨著坡度的升高而加劇，同時(shí)發(fā)現(xiàn)海拔1 000～? ? 1 200 m和1 200～1 400 m是研究區(qū)內(nèi)侵蝕發(fā)生的主要高程帶。2020年土地利用類型因子解釋力最為突出，表明退耕還林還草工程實(shí)施效果顯著，大面積的耕地向林草地轉(zhuǎn)換是使得研究區(qū)2000年后土壤侵蝕強(qiáng)度降低的最主要原因。土壤侵蝕各影響因子的協(xié)同作用明顯強(qiáng)于單一因子的影響。

關(guān)鍵詞：土壤侵蝕；修正通用土壤流失方程（RUSLE）；地理探測(cè)器；驅(qū)動(dòng)因子；延河甘谷驛流域

中圖分類號(hào)：S157.1? ? ? ? ?文獻(xiàn)標(biāo)識(shí)碼：A

文章編號(hào)：0439-8114（2024）06-0027-08

DOI：10.14088/j.cnki.issn0439-8114.2024.06.005 開放科學(xué)（資源服務(wù)）標(biāo)識(shí)碼（OSID）：

Evaluation of soil erosion and analysis of driving factors in the Ganguyi Watershed of Yanhe River from 1980 to 2020

CHEN Fang-lei1， WANG Ji-ping2，3， CHENG Fu4， XIE Hai-yan1

（1.College of Resources and Environment， Xinjiang Agricultural University， Urumqi? 830052， China； 2.Institute of Ecological Protection and Restoration， Chinese Academy of Forestry， Beijing? 100091， China； 3.Research Center of Saline and Alkali Land， National State Forestry and Grassland Administration， Beijing? 100091， China； 4.Water and Soil Conservation Monitoring Center， Ministry of Water Resources of the Peoples Republic of China， Beijing? 100053， China）

Abstract： The daily rainfall data， DEM data， soil type data， sediment content data and multi period NDVI data were used to study the soil erosion dynamics and driving factors in the Ganguyi Watershed of the Yanhe River in the past 41 years before and after the implementation of the national ecological rehabilitation project of returning farmland to forest and grassland based on the Revised General Soil Loss Equation （RUSLE） and geographic detectors. The results showed that， from 1980 to 2020， the overall soil erosion intensity in the study area showed a fluctuating trend， with an average erosion modulus of 6 746.30 t/（km2·a）， 5 740.28 t/（km2·a）， 6 389.56 t/（km2·a）， 5 450.46 t/（km2·a） and 5 480.56 t/（km2·a） in 1980， 1990， 2000， 2010 and 2020， respectively. From 1980 to 2000， the overall erosion intensity in the study area gradually increased， and the proportion of erosion areas at the strong level and above gradually increased， which was characterized by “erosion increase and upgrading”. After 2000， the intensity of soil erosion in the study area began to decrease， and the area of erosion at the strong level and above decreased， which was characterized by “erosion reduction and degradation”. The intensity of soil erosion in the study area increased with the increase of slope. At the same time， it was found that? ? ?1 000～1 200 m and 1 200～1 400 m were the main elevation zones for erosion occurrence in the study area. The explanatory power of land use type factors was most prominent in 2020， indicating that the implementation of the project of returning farmland to forests and grasslands had a significant effect. The conversion of large areas of farmland to forests and grasslands was the main reason for the decrease in soil erosion intensity in the research area after 2000. The synergistic effect of various influencing factors on soil erosion was significantly stronger than that of a single factor.

Key words： soil erosion； Revised Universal Soil Loss Equation （RUSLE）； geographic detector； driving factor； Ganguyi Watershed of Yanhe River

土壤侵蝕是指結(jié)構(gòu)不穩(wěn)定的土壤因流水或降雨而被沖刷、剝蝕、搬運(yùn)，因重力作用而失穩(wěn)移動(dòng)，因風(fēng)力作用而懸浮移動(dòng)、沉積，或因凍融作用而形成泥質(zhì)流體的現(xiàn)象［1］。中國(guó)是世界上土壤侵蝕狀況最為嚴(yán)重的國(guó)家，目前受侵蝕的土壤總面積達(dá)26萬km2，西北黃土高原地區(qū)尤為嚴(yán)重［2］。黃土高原不僅是中國(guó)水土保持工程治理與生態(tài)修復(fù)的重點(diǎn)地區(qū)，同時(shí)也是土壤侵蝕機(jī)理、治理理論、方法及防治技術(shù)研究與實(shí)踐關(guān)注的熱點(diǎn)區(qū)［3］。因此，開展黃土高原小流域范圍內(nèi)水土流失動(dòng)態(tài)變化分析變得尤為重要。

為了開展土壤侵蝕的定量評(píng)估等相關(guān)工作，國(guó)內(nèi)外水土保持研究相關(guān)學(xué)者提出了眾多土壤侵蝕相關(guān)模型，這些模型可分為物理過程性模型、產(chǎn)量模型和經(jīng)驗(yàn)?zāi)Ｐ?。物理過程性模型主要有WEPP模型等［4］；產(chǎn)量模型有EPIC［5］和PI［6］；原理簡(jiǎn)單和適用性強(qiáng)是經(jīng)驗(yàn)?zāi)Ｐ偷膬?yōu)勢(shì)所在，主要有土壤流失方程（USLE）［4］、修正通用土壤流失方程（RUSLE）［7］和中國(guó)土壤流失方程（CSLE）等［8］。黃土高原土壤侵蝕的相關(guān)研究表明，RUSLE可以更為準(zhǔn)確地模擬黃土高原土壤侵蝕狀況［9-11］。

延河流域是黃土高原水土流失最為嚴(yán)重的區(qū)域［12］，整個(gè)流域受水土流失影響，土壤肥力下降，水庫(kù)淤積，河床上升，旱災(zāi)頻發(fā)，農(nóng)業(yè)生產(chǎn)條件落后，生態(tài)環(huán)境脆弱，嚴(yán)重制約經(jīng)濟(jì)社會(huì)高質(zhì)量發(fā)展［13］。以小流域?yàn)閱卧乃亮魇ЬC合治理是中國(guó)目前防治水土流失的主要方式［14］。本研究以延河流域甘谷驛水文站控制區(qū)為研究區(qū)，運(yùn)用RUSLE模型，分析該地區(qū)1980—2020年長(zhǎng)時(shí)間研究序列的土壤侵蝕時(shí)空變化特征，闡明研究區(qū)近41年土壤侵蝕變化特征及驅(qū)動(dòng)因素，以期為該區(qū)未來制定土壤侵蝕應(yīng)對(duì)策略提供理論依據(jù)。

1 數(shù)據(jù)與方法

1.1 研究區(qū)概況

延河流域甘谷驛水文站控制區(qū)（36°22′N—37°19′、108°39′—109°48′E，圖1）位于黃河流域黃土高原中部，控制區(qū)總面積為5 872 km2。延河流域地表破碎，黃土侵蝕極為劇烈，溝間地以長(zhǎng)梁、斜梁、梁峁和峁為主，同時(shí)也有少量殘塬存在，部分地區(qū)有石質(zhì)丘陵與裸露基巖，整個(gè)流域土壤穩(wěn)定性較差，因此長(zhǎng)期遭受強(qiáng)烈剝蝕［15］。

1.2 數(shù)據(jù)來源

數(shù)據(jù)包括逐日降雨量數(shù)據(jù)、土地利用類型數(shù)據(jù)、土壤類型數(shù)據(jù)、DEM數(shù)據(jù)、歸一化植被指數(shù)（Normalized difference vegetation index，NDVI）數(shù)據(jù)、泥沙含量數(shù)據(jù)和人口數(shù)據(jù)。降雨量數(shù)據(jù)下載于國(guó)家氣象科學(xué)數(shù)據(jù)中心-中國(guó)氣象數(shù)據(jù)網(wǎng)（http：//data.cma.cn） ，選用“中國(guó)地面氣候資料日數(shù)據(jù)集（V3.0）”中延河流域甘谷驛水文站控制區(qū)內(nèi)及周邊共13個(gè)氣象站點(diǎn)（大寧、靖邊、清澗、橫山、甘泉、志丹、延長(zhǎng)、延川、延安、宜川、安塞、子長(zhǎng)、吳起）1980—2020年逐日降雨量數(shù)據(jù)；水文站輸沙模數(shù)數(shù)據(jù)來源于水利部黃河水利委員會(huì)編制的黃河泥沙公報(bào)（www.yrcc.gov.cn）。土壤類型和質(zhì)地?cái)?shù)據(jù)來源于世界土壤數(shù)據(jù)庫(kù)中國(guó)土壤數(shù)據(jù)集（Harmonized world soil database version 1.1，HWSD）。地形數(shù)據(jù)為30 m分辨率的DEM數(shù)據(jù)，來源于地理空間數(shù)據(jù)云（http：//www.gscloud.cn/）。歸一化植被指數(shù)數(shù)據(jù)來自Google Earth Engine平臺(tái)，由于研究?jī)?nèi)容時(shí)序跨度大，目前沒有單一數(shù)據(jù)產(chǎn)品可以滿足。本研究采用250 m分辨率MODIS（Moderate-resolution imaging spectroradiometer）-NDVI產(chǎn)品和8 km分辨率GIMMS（Global inventory modelling and mapping studies）-NDVI產(chǎn)品作為數(shù)據(jù)源，通過合成處理得到1981—2020年（1980年數(shù)據(jù)缺失，1980—1981年未發(fā)生重大變化）NDVI數(shù)據(jù)。1980年、1990年、2000年、2010年和2020年5期土地利用數(shù)據(jù)均來源于中國(guó)科學(xué)院地理科學(xué)與資源研究所/地理國(guó)情監(jiān)測(cè)云平臺(tái)（http：//www.dsac.cn）。人口數(shù)據(jù)來自LandScan數(shù)據(jù)庫(kù)。

1.3 土壤侵蝕模數(shù)計(jì)算

計(jì)算土壤侵蝕模數(shù)選用修正后的通用土壤流失方程（RUSLE），表達(dá)式如下。

A=R×L×S×K×C×P? ? ? ? ? ? ? ? ? ? ? ? ? （1）

式中，A為土壤侵蝕模數(shù)；R為降雨侵蝕力因子；L、S分別為坡長(zhǎng)、坡度因子；K為土壤可蝕性因子；C為植被覆蓋與管理因子；P為水土保持措施因子。

1.3.1 降雨侵蝕力因子 采用章文波等［16］提出的基于逐日降雨量的降雨侵蝕力計(jì)算方法，其表達(dá)式如下。

[R半月=αk=1n（Pk）β]? ? ?（2）

[R年=i=124R半月i]? ? ? ? ? ? ? ? ? （3）

[β=0.836 3+18.144/Pd12+24.455/Py12]? ? （4）

[α=21.586β-7.189 1]? ? ? ? ? ?（5）

相關(guān)學(xué)者根據(jù)黃河流域徑流實(shí)驗(yàn)站降雨和徑流觀測(cè)資料的分析結(jié)果擬定侵蝕性降雨量標(biāo)準(zhǔn)定為12 mm［17］。式中，[R半月]為半月降雨侵蝕力；[R年]為年降雨侵蝕力；k為某半月內(nèi)侵蝕性降雨日數(shù)；[Pk]為半月內(nèi)第k天的侵蝕性降雨日雨量；[Pd12]為一年內(nèi)侵蝕性降雨日雨量均值；[Py12]為侵蝕性降雨年總量的多年均值；[α]、[β]為該模型的兩個(gè)參數(shù)。

由于降雨具有偶然性，為了減少單個(gè)年份降雨存在的偶然性，故將1980—2020年分為5個(gè)時(shí)段計(jì)算多年平均降雨侵蝕力，第1時(shí)段為1980—1985年（時(shí)段Ⅰ），第2時(shí)段為1986—1995年（時(shí)段Ⅱ），第3時(shí)段為1996—2005年（時(shí)段Ⅲ），第4時(shí)段為2006—2015年（時(shí)段Ⅳ），第5時(shí)段為2016—2020年（時(shí)段Ⅴ）。最后采用反距離權(quán)重法［18］進(jìn)行空間插值，得到研究區(qū)各個(gè)時(shí)段降雨侵蝕力的空間分布。

1.3.2 土壤可蝕性因子 以世界土壤數(shù)據(jù)庫(kù)中國(guó)土壤數(shù)據(jù)集為基礎(chǔ)，運(yùn)用EPIC模型計(jì)算研究區(qū)內(nèi)各類型土壤的K，計(jì)算式如下。

[K=0.2+0.3exp-0.025 6SAN1-SIL100×SILCLA+SIL0.3×1-0.25CC+exp3.72-2.95C×1-0.7SN1SN1+exp-5.51+22.9SN1] ? （6）

式中，SAN為沙粒含量；SIL為粉粒含量；CLA為黏粒含量；C為有機(jī)碳含量；SN1=1-SAN/100。

1.3.3 坡長(zhǎng)、坡度因子 坡長(zhǎng)、坡度因子是評(píng)估土壤侵蝕的重要參數(shù)，本研究采用符素華等［19］提出的坡長(zhǎng)、坡度因子計(jì)算方法，計(jì)算式如下。

[S=10.8sin θ+0.03? ? ? ? ? ? ? ? θ≤5°16.8sin θ-0.05? ? ? ? ? ? ? ? 5°<θ≤10°21.9sin θ-0.96? ? ? ? ? ? ? ?θ>10°] （7）

[L=λ/22.1m， m0.2? ? ? ? ? ? ? ? ? ? θ≤1° 0.3? ? ? ? ? ? ? ? ? ? 1°<θ≤3°0.4? ? ? ? ? ? ? ? ? ? 3°<θ≤5°0.5? ? ? ? ? ? ? ? ? ? θ>5°] （8）

式中，λ表示坡長(zhǎng)；m為坡長(zhǎng)系數(shù)；[θ]為坡度。

1.3.4 植被覆蓋與管理因子 為滿足1980—2020年長(zhǎng)時(shí)間序列的研究需求，需選取2種數(shù)據(jù)源相同時(shí)序（2002—2006年）NDVI月最大合成數(shù)據(jù)，將MODIS-NDVI（250 m）平均聚合至8 km空間分辨率，得出逐月空間分布系數(shù)，并將MODIS-NDVI（8 km）與GIMMS-NDVI（8 km）進(jìn)行回歸統(tǒng)計(jì)，得到二者的回歸方程；利用所得的回歸擬合方程對(duì)GIMMS-NDVI（8 km）數(shù)據(jù)進(jìn)行修正，隨后將修正后的結(jié)果與MODIS-NDVI（250 m）逐月數(shù)據(jù)的空間分布系數(shù)相乘，最終得到空間分辨率為250 m的GIMMS-NDVI數(shù)據(jù)［20］。由于1980年GIMMS-NDVI數(shù)據(jù)缺失，該地區(qū)歷史相關(guān)資料表明研究區(qū)在1980—1981年并未發(fā)生大規(guī)模土地利用類型轉(zhuǎn)變，同時(shí)未發(fā)生重大氣候?yàn)?zāi)害，故用1981年NDVI數(shù)據(jù)作為1980年缺失數(shù)據(jù)代入公式計(jì)算植被覆蓋與管理因子。數(shù)據(jù)降尺度前后對(duì)比如圖2所示。

土壤侵蝕強(qiáng)度在不同的植被覆蓋類型和管理方式下呈現(xiàn)出差異。植被覆蓋與管理因子在土壤侵蝕計(jì)算模型中反映植被覆蓋和管理措施對(duì)土壤侵蝕的影響，取值范圍為0～1。本研究參考相關(guān)研究［21］成果對(duì)研究區(qū)C進(jìn)行賦值和計(jì)算，耕地賦值為0.44，水域和生產(chǎn)建設(shè)用地的植被覆蓋度較低，土壤侵蝕強(qiáng)度較小，故將其C賦值為0，未利用地賦值為1。根據(jù)江忠善等［22］提出的方法，草地和林地的植被覆蓋與管理因子計(jì)算式如下。

[CG=1? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? fvc≤5%exp-0.041 8fvc-5? ? ? fvc>5%] （9）

[CF=1? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? fvc≤5%exp-0.008 5fvc-51.5? ? ? fvc>5%] （10）

式中，CG為草地植被覆蓋與管理因子；CF為林地植被覆蓋與管理因子，fvc表示植被覆蓋度。根據(jù)研究區(qū)土地利用數(shù)據(jù)和基于NDVI數(shù)據(jù)計(jì)算的植被覆蓋度，得到研究區(qū)5個(gè)時(shí)期的植被覆蓋與管理因子。

1.3.5 水土保持措施因子 土地利用類型可以反映出水土保持措施因子，本研究參考文獻(xiàn)［23］，采用經(jīng)驗(yàn)賦值法得到研究區(qū)不同時(shí)期的水土保持措施因子P分布，林地、草地和耕地分別賦值0.05、0.16和0.31，水域和未利用地均賦值1.0，生產(chǎn)建設(shè)用地賦值0。

1.4 地理探測(cè)器

地理探測(cè)器是探測(cè)空間分異性并揭示其背后驅(qū)動(dòng)因子的統(tǒng)計(jì)學(xué)方法［24］。采用該工具中的單因子探測(cè)器來定量評(píng)估延河流域甘谷驛水文站控制區(qū)土壤強(qiáng)度變化的影響因子，定量分析各因子對(duì)土壤侵蝕強(qiáng)度的影響程度。在地理探測(cè)器中，運(yùn)用因子探測(cè)器來探測(cè)因變量的空間分異性，用q表示自變量對(duì)因變量的解釋度，表達(dá)式如下。

[q=1-?=1LN? σ2?Nσ2=1-SSWSST]? ? ?（11）

[SSW=?=1LN?σ2?]? ? ? ? ? ? （12）

[SST=Nσ2]? ? ? ? ? ? ? ? ? ? ?（13）

式中，L為自變量或因變量的分層或分區(qū)；Nh和N分別為層h和全區(qū)的單元數(shù)；[σ2?]和[σ2]分別為層h和全區(qū)Y值的方差；SSW和SST分別為層內(nèi)方差之和與全區(qū)總方差；q的值域?yàn)椋?，1］［25］。

交互作用探測(cè)器通過對(duì)比單因子q和雙因子q來確定2個(gè)因子間的交互作用方向和方式，交互作用的判定依據(jù)見表1。

2 結(jié)果與分析

2.1 流域土壤侵蝕強(qiáng)度及其時(shí)空變化總特征

按照上述各因子方法逐個(gè)計(jì)算1980年、1990年、2000年、2010年和2020年共5期延河流域甘谷驛水文站控制區(qū)單位面積土壤侵蝕模數(shù)，根據(jù)《土壤侵蝕分類分級(jí)標(biāo)準(zhǔn)》對(duì)研究區(qū)土壤侵蝕強(qiáng)度計(jì)算結(jié)果進(jìn)行分級(jí)，利用GIS生成延河流域甘谷驛水文站控制區(qū)土壤侵蝕強(qiáng)度分布（圖3）。查閱黃河水資源和泥沙公報(bào)中延河流域甘谷驛水文站年輸沙模數(shù)數(shù)據(jù)，對(duì)比發(fā)現(xiàn)本研究土壤侵蝕計(jì)算結(jié)果變化趨勢(shì)與延河甘谷驛水文站實(shí)測(cè)輸沙模數(shù)數(shù)據(jù)變化趨勢(shì)一致。同時(shí)發(fā)現(xiàn)1987—2010年年均輸沙模數(shù)為5 510 t/（km2·年），計(jì)算出的1980—2010年年均土壤侵蝕模數(shù)約為6 082 t/（km2·年），計(jì)算結(jié)果比輸沙模數(shù)高出10.4%，呈現(xiàn)出這樣結(jié)果的主要原因是由于RULSE是基于坡面來計(jì)算坡面侵蝕的模型，在工程措施以外同時(shí)存在泥沙沉淀和重力侵蝕，因此通過模型計(jì)算所得的侵蝕模數(shù)與水文站觀測(cè)到的輸沙模數(shù)有一定的出入。

由土壤侵蝕強(qiáng)度分級(jí)可以看出，與1980年相比，1990年土壤侵蝕強(qiáng)度空間分布特征變化較小，中度及以下的侵蝕主要分布在河谷周圍和研究區(qū)南部植被覆蓋度較高的地帶；2000年與前2期結(jié)果不同，研究區(qū)中部地帶土壤侵蝕強(qiáng)度明顯增大；2010年較2000年侵蝕情況變化最為明顯，沿河道向外兩側(cè)侵蝕強(qiáng)度顯著降低；2020年流域西北部高海拔地區(qū)侵蝕強(qiáng)度較高，研究區(qū)中部和南部區(qū)域侵蝕強(qiáng)度較低，基本上以中度、輕度和微度侵蝕為主。

土壤侵蝕模數(shù)計(jì)算結(jié)果表明，1980—2020年延河流域甘谷驛水文站控制區(qū)年均土壤侵蝕強(qiáng)度呈波動(dòng)變化的發(fā)展趨勢(shì)，1980年平均土壤侵蝕模數(shù)為6 746.30 t/（km2·年）；1990年平均侵蝕模數(shù)減小到5 740.28 t/（km2·年），比1980年減少14.9%；2000年平均侵蝕模數(shù)增加到 6 389.56 t/（km2·年），比1990年增加11.3%；2000年后研究區(qū)平均侵蝕模數(shù)開始大幅降低，2010年和2020年平均侵蝕模數(shù)分別為5 450.46、5 480.56 t/（km2·年），相較于2000年分別減少14.7%和14.2%。

表2為研究區(qū)1980—2020年各侵蝕強(qiáng)度等級(jí)的面積及占比情況。通過分析5期土壤侵蝕強(qiáng)度的面積分布可知，1980年、1990年、2000年、2010年和2020年微度和輕度侵蝕面積占比較大，分別為48.42%、50.58%、46.31%、56.02%和62.85%，強(qiáng)烈及以上等級(jí)侵蝕面積占比分別為40.06%、37.36%、42.36%、31.99%和27.90%，均小于各時(shí)期微度和輕度侵蝕所占面積。

從5期土壤侵蝕強(qiáng)度面積變化來看，1980—1990年，土壤侵蝕強(qiáng)度只有劇烈侵蝕的面積減小，面積占比也降低3.29個(gè)百分點(diǎn)，可見劇烈侵蝕面積減少是該段時(shí)間內(nèi)流域土壤侵蝕量減小的主要原因；1990—2000年，強(qiáng)烈及以上等級(jí)侵蝕面積占比累計(jì)增加了5.00個(gè)百分點(diǎn)，微度、輕度、中度侵蝕面積占比分別減少了1.84、2.43、0.73個(gè)百分點(diǎn)，可見該時(shí)間段內(nèi)土壤侵蝕量的增加與強(qiáng)烈及以上等級(jí)侵蝕面積增加有關(guān)；2000—2010年，強(qiáng)烈及以上等級(jí)侵蝕面積占比共計(jì)減少了10.37個(gè)百分點(diǎn)，可見該時(shí)間段內(nèi)土壤侵蝕強(qiáng)度的大幅減弱與強(qiáng)烈及以上等級(jí)侵蝕面積的減少有關(guān)；2010—2020年，劇烈侵蝕面積增加0.64個(gè)百分點(diǎn)，土壤侵蝕量的細(xì)微增加可能與劇烈侵蝕面積的增加有關(guān)。

2.2 流域土壤侵蝕強(qiáng)度的空間分布變化

為進(jìn)一步分析延河流域甘谷驛水文站控制區(qū)退耕還林還草措施實(shí)施前后土壤侵蝕的空間分布變化，利用GIS對(duì)侵蝕結(jié)果進(jìn)行疊加分析，得出土壤侵蝕強(qiáng)度空間轉(zhuǎn)化，結(jié)果如圖4所示。由圖4可知，1980—2000年，土壤侵蝕強(qiáng)度空間分布變化明顯，整體上呈侵蝕強(qiáng)度升級(jí)的趨勢(shì)，流域中部地區(qū)坪橋鎮(zhèn)、建華鎮(zhèn)、化子坪鎮(zhèn)、真武洞鎮(zhèn)尤為顯著；2000—2020年土壤侵蝕強(qiáng)度減弱面積明顯大于強(qiáng)度增大的面積，土壤侵蝕強(qiáng)度減弱的面積占變化總面積的67.5%，土壤侵蝕量在這個(gè)時(shí)段下降明顯，侵蝕強(qiáng)度減弱主要發(fā)生在化子坪鎮(zhèn)、招安鎮(zhèn)、真武洞鎮(zhèn)和建華鎮(zhèn)等流域中部地區(qū)的鄉(xiāng)鎮(zhèn)。

通過分析1980—2020年土壤侵蝕強(qiáng)度空間分布變化結(jié)果，得到1980—2000年、2000—2020年2個(gè)時(shí)段研究區(qū)的土壤侵蝕強(qiáng)度等級(jí)轉(zhuǎn)移矩陣，見表3和表4。1980—2000年，流域侵蝕強(qiáng)度整體呈增大趨勢(shì)，土壤侵蝕強(qiáng)度增大的區(qū)域面積占比為16.1%，說明在這期間研究區(qū)土壤侵蝕強(qiáng)度加劇。2000—2020年，47.7%的區(qū)域土壤侵蝕強(qiáng)度發(fā)生變化，其中高侵蝕等級(jí)向低侵蝕等級(jí)轉(zhuǎn)換的面積占比為35.5%，這期間土壤侵蝕強(qiáng)度降低明顯。

2.3 流域土壤侵蝕強(qiáng)度的地形分布變化

根據(jù)研究區(qū)的海拔范圍，分為<1 000 m、1 000～1 200 m、1 200～1 400 m、1 400～1 600 m和1 600～? ? 1 800 m共5個(gè)海拔段，基于GIS進(jìn)行疊加分析和統(tǒng)計(jì)，得出研究區(qū)1980—2020年不同海拔段的土壤侵蝕參數(shù)，結(jié)果見表5?？梢钥闯觯芯繀^(qū)內(nèi)土壤侵蝕強(qiáng)度與海拔間密切相關(guān)，隨著海拔的升高，土壤侵蝕強(qiáng)度先明顯增大后減小。在<1 000 m高程帶，1980—2020年土壤侵蝕面積占比較低且強(qiáng)度變化不明顯，微度侵蝕面積占比最高。從侵蝕面積占比來看，1 000～1 200 m和1 200～1 400 m是研究區(qū)內(nèi)侵蝕發(fā)生的主要高程帶，在這2個(gè)高程帶內(nèi)，2000年之前各侵蝕強(qiáng)度所占比例變化并不明顯，2000年后微度和輕度侵蝕等級(jí)面積占比明顯增大，強(qiáng)烈及以上等級(jí)侵蝕面積減少，表明導(dǎo)致2000年后研究區(qū)土壤侵蝕狀況向好發(fā)展的主要改善區(qū)域在該高程帶內(nèi)。研究區(qū)1 400～1 600 m高程帶內(nèi)各等級(jí)土壤侵蝕強(qiáng)度面積占比變化不大，強(qiáng)烈和極強(qiáng)烈侵蝕面積占比呈減小趨勢(shì)。1 600～1 800 m高程帶內(nèi)各時(shí)期土壤侵蝕面積占比幾乎不發(fā)生變化，可見該高程帶并非近年來改善水土保持措施實(shí)施的重點(diǎn)區(qū)域。

坡度是影響坡面土壤侵蝕模型計(jì)算結(jié)果中的重要影響因子［19］。利用流域高程數(shù)據(jù)提取坡度信息，將坡度按照≤5°、5°～8°、8°～15°、15°～25°和>25°進(jìn)行分級(jí)，隨后將坡度分級(jí)結(jié)果與5期土壤侵蝕強(qiáng)度結(jié)果帶入ArcGIS軟件進(jìn)行疊加分析，得到5期不同坡度的土壤侵蝕分布狀況，結(jié)果見表6?？傮w來看，研究區(qū)土壤侵蝕強(qiáng)度隨著坡度的增大而加劇。坡度

≤5°的侵蝕強(qiáng)度較弱，不發(fā)生強(qiáng)烈及以上等級(jí)的侵蝕，15°～25°和>25°是研究區(qū)內(nèi)土壤侵蝕強(qiáng)度最高的地區(qū)。1980—2000年，8°～15°坡度帶中度侵蝕面積占比增加，這也是導(dǎo)致在此期間土壤侵蝕強(qiáng)度整體上升的原因之一。2000年以后15°～25°和>25°坡度帶內(nèi)強(qiáng)烈、極強(qiáng)烈和劇烈侵蝕面積明顯減少，這與當(dāng)?shù)卦谠摃r(shí)段所實(shí)施的生態(tài)恢復(fù)相關(guān)工程有密切關(guān)聯(lián)。

2.4 流域土壤侵蝕驅(qū)動(dòng)因子分析

為進(jìn)一步正確認(rèn)識(shí)在黃土高原地區(qū)開展的“退耕還林”工程背景下，工程后研究區(qū)內(nèi)土壤侵蝕的自然與社會(huì)因素對(duì)土壤侵蝕的影響，選取研究區(qū)平均降雨量（X1）、植被覆蓋度（X2）、海拔（X3）、坡度（X4）、土地利用類型（X5）、土壤類型（X6）和人口密度（X7）7個(gè)因子作為代入運(yùn)算的自變量，以2020年土壤侵蝕強(qiáng)度作為因變量Y代入地理探測(cè)器中進(jìn)行運(yùn)算。地理探測(cè)器中要求輸入的自變量為類型數(shù)據(jù)，本研究采用王勁峰等［26］提出的數(shù)據(jù)離散化方法，使用自然斷點(diǎn)法將時(shí)段內(nèi)降雨量、高程和人口密度數(shù)據(jù)離散化為6類，土壤類型數(shù)據(jù)和土地利用按各自類別分為6類，植被覆蓋度數(shù)據(jù)分為≤0.3、0.3～0.4、0.4～0.5、0.5～0.6、0.6～0.7、0.7～0.8、0.8～0.9、0.9～1.0共8類，坡度分為≤5°、5°～10°、10°～15°、15°～20°、20°～25°、>25°共6類，采用GIS軟件的漁網(wǎng)功能，將研究區(qū)劃分成1 km×1 km格網(wǎng)，共提取5 870個(gè)采樣點(diǎn)代入地理探測(cè)器運(yùn)行，影響土壤侵蝕結(jié)果的7個(gè)因子q由高到低依次為土地利用類型（0.380 0）>坡度（0.145 0）>海拔（0.026 1）>植被覆蓋度（0.026 0）>土壤類型（0.024 0）>平均降雨量（0.011 0）>人口密度（0.009 0），結(jié)果顯示所有因子均通過顯著性檢驗(yàn)（P<0.000 1）。

因子探測(cè)器的應(yīng)用結(jié)果表明，退耕還林還草工程實(shí)施后，不同因子對(duì)土壤侵蝕強(qiáng)度的解釋力存在差異，土地利用類型因子的解釋力最強(qiáng)，q為0.380 0，是影響研究區(qū)土壤侵蝕空間分布的主導(dǎo)因子。延河流域地理位置和環(huán)境特殊，流域內(nèi)土壤穩(wěn)定性較差，生態(tài)環(huán)境脆弱。自1999年開始的退耕還林還草工程實(shí)施以來，黃土高原土地利用類型發(fā)生根本性改變［27，28］，根據(jù)計(jì)算分析得到研究區(qū)2000—2020年耕地向林地和草地轉(zhuǎn)移面積分別為204.75 km2和502.68 km2，植被覆蓋度顯著增加，生態(tài)環(huán)境得到明顯改善，根據(jù)探測(cè)結(jié)果，研究區(qū)內(nèi)坡度因子也是影響土壤侵蝕的另一主導(dǎo)因子。上述結(jié)果表明，自然和人為因素共同影響著延河流域甘谷驛水文站控制區(qū)土壤侵蝕空間分布格局。

利用交互探測(cè)器研究了延河流域甘谷驛水文站控制區(qū)土壤侵蝕和各驅(qū)動(dòng)因素間的交互作用，發(fā)現(xiàn)各因素間并不存在相互獨(dú)立的作用，而是表現(xiàn)為非線性增強(qiáng)和雙因子增強(qiáng)2種交互作用，結(jié)果如表7所示。平均降雨量與土壤類型、平均降雨量與人口密度、海拔與坡度、海拔與土壤類型、海拔與人口密度、坡度與土壤類型間的交互作用以雙因子增強(qiáng)的形式影響土壤侵蝕分布格局變化，土壤類型與人口密度交互探測(cè)結(jié)果為單因子非線性增強(qiáng)，其他因子交互探測(cè)結(jié)果均為非線性增強(qiáng)。土地利用類型與植被覆蓋度、海拔、坡度、土壤類型的交互影響力均在0.40以上，其中土地利用類型與坡度的交互影響力為0.600 0，表明這2個(gè)因子組合對(duì)研究區(qū)土壤侵蝕分布格局影響最大。

表7 2020年土壤侵蝕各因子交互探測(cè)結(jié)果

[q1 q2 q1+q2 q1∩q2 結(jié)果 X1=0.011 0 X2=0.026 0 0.037 0 0.045 0 非線性增強(qiáng) X1=0.011 0 X3=0.026 1 0.037 1 0.043 0 非線性增強(qiáng) X1=0.011 0 X4=0.145 0 0.156 0 0.166 0 非線性增強(qiáng) X1=0.011 0 X5=0.380 0 0.391 0 0.398 0 非線性增強(qiáng) X1=0.011 0 X6=0.024 0 0.035 0 0.034 0 雙因子增強(qiáng) X1=0.011 0 X7=0.009 0 0.020 0 0.016 0 雙因子增強(qiáng) X2=0.026 0 X3=0.026 1 0.052 1 0.056 0 非線性增強(qiáng) X2=0.026 0 X4=0.145 0 0.171 0 0.198 0 非線性增強(qiáng) X2=0.026 0 X5=0.380 0 0.406 0 0.427 0 非線性增強(qiáng) X2=0.026 0 X6=0.024 0 0.050 0 0.060 0 非線性增強(qiáng) X2=0.026 0 X7=0.009 0 0.035 0 0.047 0 非線性增強(qiáng) X3=0.026 1 X4=0.145 0 0.171 1 0.167 0 雙因子增強(qiáng) X3=0.026 1 X5=0.380 0 0.406 1 0.410 0 非線性增強(qiáng) X3=0.026 1 X6=0.024 0 0.050 1 0.043 0 雙因子增強(qiáng) X3=0.026 1 X7=0.009 0 0.035 1 0.034 0 雙因子增強(qiáng) X4=0.145 0 X5=0.380 0 0.525 0 0.600 0 非線性增強(qiáng) X4=0.145 0 X6=0.024 0 0.169 0 0.166 0 雙因子增強(qiáng) X4=0.145 0 X7=0.009 0 0.154 0 0.161 0 非線性增強(qiáng) X5=0.380 0 X6=0.024 0 0.404 0 0.420 0 非線性增強(qiáng) X5=0.380 0 X7=0.009 0 0.389 0 0.393 0 非線性增強(qiáng) X6=0.024 0 X7=0.009 0 0.033 0 0.013 2 單因子非線性增強(qiáng) ]

3 小結(jié)

1）利用RUSLE模型對(duì)延河流域甘谷驛水文站控制區(qū)1980—2020年土壤侵蝕強(qiáng)度進(jìn)行計(jì)算，結(jié)果表明，1980年平均土壤侵蝕模數(shù)為6 746.30 t/（km2·年）；1990年平均侵蝕模數(shù)減小到5 740.28 t/（km2·年），比1980年減少14.9%；2000年平均侵蝕模數(shù)增加到 6 389.56 t/（km2·年），比1990年增加11.3%；2000年后研究區(qū)平均侵蝕模數(shù)開始大幅降低，2010年和2020年平均侵蝕模數(shù)相較于2000年分別減少14.7%和14.2%。

2）延河流域甘谷驛水文站控制區(qū)1980—2000年土壤侵蝕強(qiáng)度逐漸增強(qiáng)，強(qiáng)烈及以上等級(jí)侵蝕面積占比逐漸增加，表現(xiàn)為“增蝕升級(jí)”的特點(diǎn)。2000年后研究區(qū)內(nèi)土壤侵蝕強(qiáng)度開始降低，強(qiáng)烈及以上等級(jí)的侵蝕面積減少，輕度和微度侵蝕面積增大，總體表現(xiàn)為“減蝕降級(jí)”的特點(diǎn)。結(jié)合DEM數(shù)據(jù)分析得出研究區(qū)土壤侵蝕強(qiáng)度隨著坡度的升高而加劇，坡度≤5°的侵蝕強(qiáng)度較弱，不發(fā)生強(qiáng)烈及以上等級(jí)的侵蝕，15°～25°和>25°是研究區(qū)發(fā)生土壤侵蝕最為嚴(yán)重的地區(qū)，同一坡度區(qū)間內(nèi)各時(shí)段土壤侵蝕強(qiáng)度變化特征同區(qū)域整體變化特征一致。同時(shí)發(fā)現(xiàn)1 000～1 200 m和1 200～1 400 m是研究區(qū)內(nèi)侵蝕發(fā)生的主要高程帶。

3）通過地理探測(cè)器對(duì)研究區(qū)退耕還林還草工程實(shí)施背景下土壤侵蝕強(qiáng)度的影響因素分析發(fā)現(xiàn)，土地利用類型因子解釋力較為突出，表明退耕還林還草工程實(shí)施后水土流失治理效果顯著，大面積的耕地向林草地轉(zhuǎn)換是研究區(qū)2000年后土壤侵蝕強(qiáng)度降低的最主要原因。交互探測(cè)結(jié)果表明，各影響因子的協(xié)同作用明顯強(qiáng)于單一因子的影響。

本研究以RUSLE模型計(jì)算結(jié)果作為基礎(chǔ)，對(duì)研究區(qū)1980—2020年土壤侵蝕的時(shí)空變化規(guī)律進(jìn)行探索，使用地理探測(cè)器對(duì)研究區(qū)7個(gè)影響因子進(jìn)行探測(cè)分析，通過對(duì)研究結(jié)果分析可知，研究中采用的方法合理，預(yù)期研究任務(wù)基本完成。相較以往研究，本研究著重突出近41年長(zhǎng)時(shí)間序列數(shù)據(jù)的變化，可以較好地降低研究時(shí)段內(nèi)數(shù)據(jù)突變的影響，同時(shí)能更加全面地掌握研究區(qū)土壤侵蝕的變化特點(diǎn)和歸因，使研究結(jié)果更有價(jià)值。本研究仍存在不足之處，在計(jì)算結(jié)果驗(yàn)證中，如果可以獲取研究區(qū)水土保持實(shí)地監(jiān)測(cè)數(shù)據(jù)，模型計(jì)算結(jié)果會(huì)更具說服力。針對(duì)地理探測(cè)器的使用，如果在今后研究中可以加入GDP和種植作物類型數(shù)據(jù)等社會(huì)經(jīng)濟(jì)因子，研究結(jié)果會(huì)更有價(jià)值。

參考文獻(xiàn)：

［1］ 洪宇辰，姬鑫慧，戴曉愛.基于GIS/RS和USLE模型的土壤侵蝕研究及主要侵蝕因子識(shí)別［J］.測(cè)繪， 2017， 40（6）： 277-283.

［2］ CHEN R Y，YAN D C，WEN A B， et al. The regional difference in engineering-control and tillage factors of Chinese Soil Loss Equation［J］. Journal of mountain science， 2021， 18（3）： 658-670.

［3］ 火 紅.基于地形因子的延安市土壤侵蝕及其地表驅(qū)動(dòng)因素的時(shí)空分異特征研究［D］.西安：長(zhǎng)安大學(xué)， 2021.

［4］ 陳朝良，趙廣舉，穆興民，等.基于RUSLE模型的湟水流域土壤侵蝕時(shí)空變化［J］.水土保持學(xué)報(bào)， 2021， 35（4）： 73-79.

［5］ 方興義.基于EPIC模型的農(nóng)業(yè)旱災(zāi)風(fēng)險(xiǎn)模糊評(píng)估方法［J］.災(zāi)害學(xué)，2020， 35（3）： 55-58.

［6］ 段興武，謝 云，張玉平，等. PI模型在東北松嫩黑土區(qū)土壤生產(chǎn)力評(píng)價(jià)中的應(yīng)用［J］.中國(guó)農(nóng)學(xué)通報(bào)， 2010， 26（8）： 179-188.

［7］ RENARD K G，F(xiàn)ERREIRA V A. RUSLE model description and database sensitivity［J］. Journal of environmental quality，1993，22（3）：458-466.

［8］ 周正朝，上官周平.土壤侵蝕模型研究綜述［J］.中國(guó)水土保持科學(xué)， 2004，2（1）： 52-56.

［9］ 郭 達(dá)，宋小寧，董 震，等.基于RUSLE與GIS的黃土高原水土流失評(píng)價(jià)研究——以寧夏中衛(wèi)地區(qū)為例［J］.泥沙研究， 2020， 45（5）： 55-60.

［10］ 王澤宇，陳旭陽(yáng)，馬彩詩(shī)，等.陜北榆林市退耕還林前后土壤侵蝕及生態(tài)服務(wù)價(jià)值變化［J］.西北林學(xué)院學(xué)報(bào)， 2021， 36（3）： 59-67.

［11］ 張 素，熊東紅，吳 漢，等.基于RUSLE模型的孫水河流域土壤侵蝕空間分異特征［J］.水土保持學(xué)報(bào)， 2021， 35（5）： 24-30.

［12］ 鐘 雪.延河流域社會(huì)水文耦合規(guī)律研究［D］.西安：西北大學(xué)， 2020.

［13］ 張永峰，王明華.延安市水土保持高質(zhì)量發(fā)展先行區(qū)示范創(chuàng)建［J］.中國(guó)水土保持， 2022（9）： 37-40.

［14］ 王福嶺，吳 鵬，郭子琦.以小流域?yàn)閱卧乃亮魇?qiáng)度評(píng)估［J］.中國(guó)水土保持， 2023（1）： 44-47.

［15］ 孫 虎.陜西延河流域地貌組合類型的模糊聚類劃分［J］.陜西師范大學(xué)學(xué)報(bào)（自然科學(xué)版）， 1996，24（4）： 83-88.

［16］ 章文波，謝 云，劉寶元.利用日雨量計(jì)算降雨侵蝕力的方法研究［J］.地理科學(xué)， 2002，22（6）： 705-711.

［17］ 謝 云，劉寶元，章文波.侵蝕性降雨標(biāo)準(zhǔn)研究［J］.水土保持學(xué)報(bào)，2000，14（4）： 6-11.

［18］ 鄒 強(qiáng).流域降雨空間插值方法比較［J］.節(jié)水灌溉， 2009（7）： 12-14.

［19］ 符素華， 劉寶元， 周貴云， 等. 坡長(zhǎng)坡度因子計(jì)算工具［J］. 中國(guó)水土保持科學(xué)， 2015， 13（5）： 105-110.

［20］ 林錦闊.河西地區(qū)土壤侵蝕時(shí)空分異及其驅(qū)動(dòng)因素［D］.蘭州：蘭州大學(xué)， 2020.

［21］ 秦 偉，朱清科，張 巖.基于GIS和RUSLE的黃土高原小流域土壤侵蝕評(píng)估［J］.農(nóng)業(yè)工程學(xué)報(bào)， 2009， 25（8）：157-163，4.

［22］ 江忠善，鄭粉莉，武 敏.中國(guó)坡面水蝕預(yù)報(bào)模型研究［J］.泥沙研究， 2005（4）： 1-6.

［23］ 閆 瑞，張曉萍，李夠霞，等.基于RUSLE的北洛河上游流域侵蝕產(chǎn)沙模擬研究［J］.水土保持學(xué)報(bào)， 2017， 31（4）： 32-37.

［24］ 宋麗潔，戴昭鑫，劉 新. 2000—2020年?yáng)|營(yíng)市人口分布的時(shí)空特征演變規(guī)律及影響因素分析［J］.北京測(cè)繪， 2022， 36（3）： 260-265.

［25］ 姜旭海，韓 玲，白宗璠，等.內(nèi)蒙古自治區(qū)沙漠化敏感性時(shí)空演變格局和趨勢(shì)分析［J］.生態(tài)學(xué)報(bào)， 2023， 43（1）： 364-378.

［26］ 王勁峰，徐成東.地理探測(cè)器：原理與展望［J］.地理學(xué)報(bào)，2017，72（1）：116-134.

［27］ 龔直文，姚順波，黨晶晶.黃土高原退耕前后土地利用/覆被變化及驅(qū)動(dòng)力分析——以志丹縣為例［J］.林業(yè)經(jīng)濟(jì)，2017，39（2）：54-58， 66.

［28］ 楊 燦，魏天興，李亦然，等.黃土高原水蝕風(fēng)蝕交錯(cuò)區(qū)退耕還林工程前后NDVI時(shí)空變化特征［J］.北京林業(yè)大學(xué)學(xué)報(bào)， 2021， 43（6）： 83-91.