1980—2020年延河甘谷驛流域土壤侵蝕評價與驅動因子分析

陳方磊 王計平 程復 謝海燕

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

關鍵詞：土壤侵蝕；修正通用土壤流失方程（RUSLE）；地理探測器；驅動因子；延河甘谷驛流域

中圖分類號：S157.1? ? ? ? ?文獻標識碼：A

文章編號：0439-8114（2024）06-0027-08

DOI：10.14088/j.cnki.issn0439-8114.2024.06.005 開放科學（資源服務）標識碼（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

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

為了開展土壤侵蝕的定量評估等相關工作，國內外水土保持研究相關學者提出了眾多土壤侵蝕相關模型，這些模型可分為物理過程性模型、產量模型和經驗模型。物理過程性模型主要有WEPP模型等［4］；產量模型有EPIC［5］和PI［6］；原理簡單和適用性強是經驗模型的優(yōu)勢所在，主要有土壤流失方程（USLE）［4］、修正通用土壤流失方程（RUSLE）［7］和中國土壤流失方程（CSLE）等［8］。黃土高原土壤侵蝕的相關研究表明，RUSLE可以更為準確地模擬黃土高原土壤侵蝕狀況［9-11］。

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

1 數據與方法

1.1 研究區(qū)概況

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

1.2 數據來源

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

1.3 土壤侵蝕模數計算

計算土壤侵蝕模數選用修正后的通用土壤流失方程（RUSLE），表達式如下。

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

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

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

[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）

相關學者根據黃河流域徑流實驗站降雨和徑流觀測資料的分析結果擬定侵蝕性降雨量標準定為12 mm［17］。式中，[R半月]為半月降雨侵蝕力；[R年]為年降雨侵蝕力；k為某半月內侵蝕性降雨日數；[Pk]為半月內第k天的侵蝕性降雨日雨量；[Pd12]為一年內侵蝕性降雨日雨量均值；[Py12]為侵蝕性降雨年總量的多年均值；[α]、[β]為該模型的兩個參數。

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

1.3.2 土壤可蝕性因子 以世界土壤數據庫中國土壤數據集為基礎，運用EPIC模型計算研究區(qū)內各類型土壤的K，計算式如下。

[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為有機碳含量；SN1=1-SAN/100。

1.3.3 坡長、坡度因子 坡長、坡度因子是評估土壤侵蝕的重要參數，本研究采用符素華等［19］提出的坡長、坡度因子計算方法，計算式如下。

[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）

式中，λ表示坡長；m為坡長系數；[θ]為坡度。

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

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

[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表示植被覆蓋度。根據研究區(qū)土地利用數據和基于NDVI數據計算的植被覆蓋度，得到研究區(qū)5個時期的植被覆蓋與管理因子。

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

1.4 地理探測器

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

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

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

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

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

交互作用探測器通過對比單因子q和雙因子q來確定2個因子間的交互作用方向和方式，交互作用的判定依據見表1。

2 結果與分析

2.1 流域土壤侵蝕強度及其時空變化總特征

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

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

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

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

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

2.2 流域土壤侵蝕強度的空間分布變化

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

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

2.3 流域土壤侵蝕強度的地形分布變化

根據研究區(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個海拔段，基于GIS進行疊加分析和統(tǒng)計，得出研究區(qū)1980—2020年不同海拔段的土壤侵蝕參數，結果見表5。可以看出，研究區(qū)內土壤侵蝕強度與海拔間密切相關，隨著海拔的升高，土壤侵蝕強度先明顯增大后減小。在<1 000 m高程帶，1980—2020年土壤侵蝕面積占比較低且強度變化不明顯，微度侵蝕面積占比最高。從侵蝕面積占比來看，1 000～1 200 m和1 200～1 400 m是研究區(qū)內侵蝕發(fā)生的主要高程帶，在這2個高程帶內，2000年之前各侵蝕強度所占比例變化并不明顯，2000年后微度和輕度侵蝕等級面積占比明顯增大，強烈及以上等級侵蝕面積減少，表明導致2000年后研究區(qū)土壤侵蝕狀況向好發(fā)展的主要改善區(qū)域在該高程帶內。研究區(qū)1 400～1 600 m高程帶內各等級土壤侵蝕強度面積占比變化不大，強烈和極強烈侵蝕面積占比呈減小趨勢。1 600～1 800 m高程帶內各時期土壤侵蝕面積占比幾乎不發(fā)生變化，可見該高程帶并非近年來改善水土保持措施實施的重點區(qū)域。

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

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

2.4 流域土壤侵蝕驅動因子分析

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

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

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

表7 2020年土壤侵蝕各因子交互探測結果

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

3 小結

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

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

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

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

參考文獻：

［1］ 洪宇辰，姬鑫慧，戴曉愛.基于GIS/RS和USLE模型的土壤侵蝕研究及主要侵蝕因子識別［J］.測繪， 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］ 火 紅.基于地形因子的延安市土壤侵蝕及其地表驅動因素的時空分異特征研究［D］.西安：長安大學， 2021.

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

［5］ 方興義.基于EPIC模型的農業(yè)旱災風險模糊評估方法［J］.災害學，2020， 35（3）： 55-58.

［6］ 段興武，謝 云，張玉平，等. PI模型在東北松嫩黑土區(qū)土壤生產力評價中的應用［J］.中國農學通報， 2010， 26（8）： 179-188.

［7］ RENARD K G，FERREIRA V A. RUSLE model description and database sensitivity［J］. Journal of environmental quality，1993，22（3）：458-466.

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

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

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

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

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

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

［14］ 王福嶺，吳 鵬，郭子琦.以小流域為單元的水土流失強度評估［J］.中國水土保持， 2023（1）： 44-47.

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

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

［17］ 謝 云，劉寶元，章文波.侵蝕性降雨標準研究［J］.水土保持學報，2000，14（4）： 6-11.

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

［19］ 符素華， 劉寶元， 周貴云， 等. 坡長坡度因子計算工具［J］. 中國水土保持科學， 2015， 13（5）： 105-110.

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

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

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

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

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

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

［26］ 王勁峰，徐成東.地理探測器：原理與展望［J］.地理學報，2017，72（1）：116-134.

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

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