土壤侵蝕時空分布規(guī)律及侵蝕熱點分析
——以陜西省王益區(qū)為例

陳四賓，霍艾迪，2，陳思名，趙志欣，陳 建

(1.長安大學 水利與環(huán)境學院，西安 710054； 2.長安大學 旱區(qū)地下水文與生態(tài)效應(yīng)教育部重點實驗室，西安 710054)

1 研究背景

土壤侵蝕是全球性的環(huán)境問題之一[1-2]，它會減少表層土壤的養(yǎng)分含量，破壞土壤肥力造成土地退化等問題[3-6]。嚴重的土壤侵蝕還會影響區(qū)域農(nóng)業(yè)經(jīng)濟發(fā)展[7-8]。黃土臺塬是寶貴的宜農(nóng)、宜耕、宜居的土地資源，其塬面地勢平坦，土地肥沃，易于耕作。然而受土壤侵蝕問題的困擾，黃土塬周圍溝道叢生，溝頭不斷前進導致塬面萎縮[9]，嚴重影響當?shù)氐慕?jīng)濟發(fā)展。自20世紀50年代以來，黃河流域開展的小流域治理、植樹造林、退耕還林等治理工程成效顯著，但是在治理程度方面仍然存在區(qū)域差異[10]。定量分析黃土高塬溝壑區(qū)的土壤侵蝕狀況及其時空變化特征，對當?shù)赝寥狼治g防治和生態(tài)建設(shè)具有重要意義。

一般而言，受研究區(qū)域客觀條件限制，采用實地觀測的土壤侵蝕評價方法往往存在較大困難[11]，實際應(yīng)用過程中局限性較大，因此常采用土壤侵蝕預(yù)報模型進行土壤侵蝕的定量評價。從20世紀50年代開始，我國學者逐步結(jié)合中國土壤侵蝕的特點開展土壤侵蝕模型的研究，有關(guān)黃土高原地區(qū)侵蝕模型的研究也得到了快速發(fā)展[12]。劉寶元等[13]基于通用土壤流失方程(Universal Soil Loss Equation，USLE)利用黃土丘陵溝壑區(qū)實測資料建立了中國土壤流失預(yù)報方程(Chinese Soil Loss Equation，CSLE)；Zheng等[14]根據(jù)黃土高原不同小流域的實測資料建立了預(yù)測次洪產(chǎn)沙量的比例函數(shù)，流域模擬產(chǎn)沙結(jié)果比較準確。物理模型方面，王建勛等[15]研究發(fā)現(xiàn)水蝕預(yù)報模型(Water Erosion Prediction Project，WEPP)坡面版在黃土丘陵溝壑區(qū)具有良好的適用性；此外，Huo等[16]基于水土評價工具(Soil and Water Assessment Tool，SWAT)模型模擬了不同固溝保塬工程下隴東黃土高原典型流域的徑流泥沙變化。隨著現(xiàn)代技術(shù)的快速發(fā)展，采用遙感和地理信息系統(tǒng)(GIS)技術(shù)進行土壤侵蝕的評價分析，可以快速、準確地獲取具有一定認可程度的研究結(jié)果[17]。

在生態(tài)保護與城市建設(shè)發(fā)展的雙重影響背景下，如何科學防治土壤侵蝕，協(xié)調(diào)好人地關(guān)系成為當前城市建設(shè)和生態(tài)發(fā)展工作中亟須解決的問題。本文以典型的黃土高塬溝壑區(qū)——銅川市王益區(qū)為研究對象，分析了其2000—2019年的土地利用和植被覆蓋度變化，依據(jù)水利部頒布的《土壤侵蝕分類分級標準》(SL 190—2007)[18](以下簡稱“標準”)估算分析王益區(qū)土壤侵蝕強度變化及其時空分布特征?；谕寥狼治g強度風險指數(shù)(Soil Erosion Dangerous Index，SEDI)和探索性空間數(shù)據(jù)分析(Exploratory Spatial Data Analysis-GIS，ESDA-GIS)方法進一步分析了區(qū)域內(nèi)的土壤侵蝕強度風險狀況和土壤侵蝕冷熱點分布(土壤侵蝕高強度集聚區(qū)和低強度集聚區(qū))情況，以期為王益區(qū)水土流失重點防治區(qū)的劃分、固溝保塬工程實施和生態(tài)環(huán)境保護提供一定的科學依據(jù)。

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

2.1 研究區(qū)概況

王益區(qū)位于陜西省銅川市中部，地處渭北黃土高原中部南緣地帶(圖1)，地理坐標為108°57′E—109°09′E，34°57′N—35°12′N，全區(qū)面積162.2 km2。經(jīng)過多年的地殼運動，現(xiàn)整體地勢西北高、東南低，區(qū)域內(nèi)梁峁川塬山丘交錯相間，溝壑縱橫。西北部是以黃土梁為主的丘陵山地，中部的黃土川塬地帶面積大，黃土厚度深，東南部為黃土丘陵，海拔1 000～1 400 m，多溝壑殘塬，地形復雜。該地區(qū)屬于暖溫帶大陸性季風氣候，四季分明，年平均氣溫8.7～12.1 ℃，年均降水量589.2 mm，年蒸發(fā)量837.1 mm，年均無霜期182 d。

王益區(qū)地處半干旱地帶，屬于典型的黃土高塬溝壑區(qū)，土壤侵蝕以水力侵蝕[19]為主，黃土塬塬面受土壤侵蝕影響嚴重，塬面萎縮，農(nóng)耕地肥力下降。區(qū)域內(nèi)土壤侵蝕風險可能進一步加劇。

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

本文采用王益區(qū)數(shù)字高程模型(Digital Elevation Models，DEM)數(shù)據(jù)、土地利用數(shù)據(jù)以及植被覆蓋數(shù)據(jù)。DEM數(shù)據(jù)來自地理空間數(shù)據(jù)云網(wǎng)站(www.gscloud.cn)，分辨率為30 m?；诘乩砜臻g數(shù)據(jù)云網(wǎng)站(www.gscloud.cn)下載2000年的Landsat ETM+影像和2019年的Landsat OLI影像，提取王益區(qū)植被覆蓋度變化，兩期影像成像時間分別為9月份和8月份，軌道號127/36，影像成像時間接近且質(zhì)量較好。在1.96 m分辨率遙感影像支持下，結(jié)合Landsat遙感影像的光譜特征，采用人機交互解譯的方法，判別地物的幾何形狀，顏色紋理特征和空間分布情況，目視勾繪出地類邊界，解譯得到所需土地利用數(shù)據(jù)。

2.3 指標因子提取

本文基于遙感影像，采用像元二分模型[20]進行植被覆蓋度的估算，李苗苗等[21]基于歸一化植被指數(shù)(NDVI)構(gòu)建的像元二分模型在植被覆蓋度提取方面可操作性強，精度較高。相比較于直接利用NDVI反映植被覆蓋度，該方法在削減植被類型、土壤和大氣的影響方面更具有優(yōu)勢[22]。其計算公式為

(1)

式中：FVC表示植被覆蓋度(%)；NDVIsoil表示裸土或無植被覆蓋區(qū)域的NDVI值，其變化范圍在-0.1～0.2之間[21，23]；NDVIveg表示完全由植被覆蓋區(qū)域的NDVI值。參照NDVI的頻率統(tǒng)計表，分別取累計頻率5%和95%的NDVI值作為NDVIsoil和NDVIveg的取值。

基于DEM數(shù)據(jù)，利用ArcGIS軟件進行研究區(qū)坡度的提取，并根據(jù)相應(yīng)標準進行重分類，由于研究區(qū)坡度>35°的區(qū)域面積較小，故將(25°，35°]和>35°區(qū)間統(tǒng)一合并為>25°一個區(qū)間。依據(jù)相應(yīng)標準劃分為(0，5°]、(5°，8°]、(8°，15°]、(15°，25°]、>25°五個坡度區(qū)間。

2.4 土壤侵蝕量估算

采用水利部頒布的《土壤侵蝕分類分級標準》(SL 190—2007)，可以基于植被覆蓋度、坡度和土地利用較為快速準確地獲取研究區(qū)土壤侵蝕強度分布狀況[24]。

王益區(qū)土壤侵蝕以水力侵蝕為主[19]，根據(jù)面蝕分級指標(表1)，將土壤侵蝕強度劃分為“微度”“輕度”“中度”“強烈”“極強烈”5類，侵蝕模數(shù)分別為<1 000、[1 000，2 500)、[2 500，5 000)、[5 000，8 000)、[8 000，15 000) t/(km2·a)。對于王益區(qū)城鄉(xiāng)居民用地和工礦用地，參考城市化土壤侵蝕等級分級評分標準中以斜坡高差、植被覆蓋度等為評價項目進行分級的方法[25-27]，以及侵蝕強度調(diào)查中，對非耕地和未利用地室內(nèi)測算采用坡度和地面物質(zhì)為主要指標的分類方法[28]。對大面積地面硬化條件較好的建成城區(qū)，土壤侵蝕強度以微度侵蝕為主[29]；其余工礦裸露用地主要考慮坡度、下墊面條件及區(qū)域植被覆蓋情況進行侵蝕強度劃分。此外，基于2019年的高分辨率影像及野外調(diào)查發(fā)現(xiàn)，王益區(qū)梯田面積較大，在使用“標準”進行土壤侵蝕強度劃分時，為區(qū)別于坡耕地，參考通用土壤流失方程模型中P因子(采取水土保持措施與不采取任何措施的地塊土壤流失量的比率)的取值[30]，判斷不同坡度下梯田水土保持作用的折減效果，并對區(qū)域內(nèi)梯田的土壤侵蝕強度進行劃分。而對于坡度>24°的梯田，水土保持效果不明顯[31]，仍將其劃分為極強烈侵蝕。

表1 土壤侵蝕分類分級標準

王益區(qū)屬于西北黃土高原區(qū)，容許土壤流失量為1 000 t/(km2·a)，即該地區(qū)允許存在微度侵蝕。研究過程中，取各個土壤侵蝕模數(shù)區(qū)間的中值作為平均侵蝕模數(shù)進行王益區(qū)土壤侵蝕量的估算。

2.5 土壤侵蝕強度風險分析

為進一步分析王益區(qū)土壤侵蝕風險狀況，采取SEDI分析研究區(qū)內(nèi)的土壤侵蝕嚴重程度。SEDI可以綜合反映出不同土壤侵蝕強度對區(qū)域內(nèi)侵蝕狀況的貢獻程度大小，評判土壤侵蝕嚴重情況，計算公式[32-33]為

(2)

式中；Fi表示第i級侵蝕強度的量化分級值；Ai表示評價區(qū)域第i級侵蝕強度的面積；A表示評價區(qū)域的總面積；n為土壤侵蝕強度級別數(shù)目；100為擴大系數(shù)。微度侵蝕、輕度侵蝕、中度侵蝕、強烈侵蝕、極強烈侵蝕對應(yīng)的量化分級值分別為0、2、4、6、8。相應(yīng)的土壤侵蝕強度風險等級分為輕微(SEDI<50)、一般(50200)。

2.6 空間數(shù)據(jù)分析

全局空間自相關(guān)分析是利用全局Moran’ I指數(shù)統(tǒng)計不同土壤侵蝕強度之間的平均關(guān)聯(lián)程度和顯著性，描述土壤侵蝕強度空間上整體的關(guān)聯(lián)程度和差異性[34-35]，計算公式為

(3)

(4)

式中：Z為正態(tài)檢驗統(tǒng)計量；E(I)和VAR(I)分別為期望值和方差，若其通過5%的顯著性水平檢驗且Z>1.96，則說明具有顯著的空間自相關(guān)。

局部空間自相關(guān)可以反映區(qū)域內(nèi)不同土壤侵蝕強度之間的空間相關(guān)程度，彌補Moran’ I不能反映局部相關(guān)類型和聚集區(qū)的不足[34-35]，計算公式為

(5)

式中Ii表示空間關(guān)聯(lián)局域指標，取值在-1～1之間，其指標值判定與Moran’ I指數(shù)一致。

(6)

(7)

3 結(jié)果與分析

3.1 王益區(qū)土地利用時空變化特征

3.1.1 土地利用時間變化特征

王益區(qū)現(xiàn)有土地利用類型中，耕地(坡耕地+梯田)面積最大(圖2)，2000年和2019年分別為107.12 km2和79.85 km2，其中2019年耕地面積中梯田占比約47.63%，為38.03 km2。研究期間，水體面積占比極小且變化不明顯，除坡耕地外，其他土地利用類型面積均不同幅度的增加。坡耕地作為其他土地利用類型轉(zhuǎn)入量的最大來源，2019年其面積占比減少至25.78%。與此同時，林草地和城鄉(xiāng)居民用地呈現(xiàn)較大幅度增長，前者增幅最為顯著(39%)，高達17.05 km2。城鄉(xiāng)居民用地增加了10.15 km2。耕地、林草地、城鄉(xiāng)居民用地的變化，側(cè)面反映出研究區(qū)內(nèi)退耕還林的成效以及城市化的進程。

研究期間土地利用變化主要集中在林草地、耕地和城鄉(xiāng)居民用地、梯田的變化上。因此土地利用轉(zhuǎn)移的空間分布主要考慮6種土地利用變化類型(圖2)。

3.1.2 土地利用空間變化特征

王益區(qū)溝壑縱橫，地勢起伏變化大，(8°，25°]坡度區(qū)間面積為94.24 km2，占全區(qū)面積的58.10%。空間分布上，隨著坡度增加，各土地利用類型的面積總體呈現(xiàn)出先增后減的趨勢(圖3)，不同時期林草地多集中在坡度>8°區(qū)域內(nèi)，在(15°，25°]坡度區(qū)間內(nèi)達到峰值，該坡度區(qū)間內(nèi)的林草地面積分別占林草地總面積的38.11%和38.54%。2000年(0°，25°]坡度區(qū)間內(nèi)的優(yōu)勢土地利用類型均為耕地，到2019年，(15°，25°]區(qū)間內(nèi)的優(yōu)勢土地利用類型則轉(zhuǎn)化為林草地。城鄉(xiāng)居民工礦用地的總體面積占比較小且主要均勻分布在(0°，15°]坡度區(qū)間內(nèi)。

總體而言，各土地利用類型雖呈現(xiàn)出不同幅度的變化，但是除(15°，25°]坡度區(qū)間外，其他坡度區(qū)間內(nèi)的優(yōu)勢土地利用類型未發(fā)生變化。

3.2 林草地植被覆蓋度變化

根據(jù)相應(yīng)標準將植被覆蓋度等級劃分為5類，分別為(0，30%]，(30%，45%]，(45%，60%]，(60%，75%]，>75%。研究期間王益區(qū)植被覆蓋度變化如圖4所示。

植被覆蓋度FVC>60%的林草地面積大幅增長，增加了18.93 km2，約占2000年林草地總面積的42.94%；其中植被覆蓋度>75%的林草地面積增幅最大，增加了約0.7倍。此外，植被覆蓋度<60%的林草地面積小幅下降，其面積占比由25.5%減少到15.4%。林草地面積增加的同時伴隨著高植被覆蓋度面積的增長，研究區(qū)內(nèi)的植被覆蓋整體上明顯好轉(zhuǎn)。

結(jié)合研究區(qū)DEM(圖1)，土地利用及其空間分布(圖2、圖3)，王益區(qū)林草地多分布于>8°的高坡度區(qū)域，高植被覆蓋度比例的增長，說明王益區(qū)封山禁牧取得了一定的成效，對退耕還林成果的鞏固較好。

3.3 2000—2019年王益區(qū)土壤侵蝕時空變化特征

3.3.1 土壤侵蝕時間變化特征

根據(jù)相應(yīng)標準，2000—2019年王益區(qū)土壤侵蝕面積變化如圖5所示：2000—2019年王益區(qū)土壤侵蝕面積由115.04 km2減少至90.58km2，侵蝕比例由70.92%減少至55.84%。侵蝕強度結(jié)構(gòu)由中度侵蝕>強烈侵蝕>輕度侵蝕>極強烈侵蝕變化為輕度侵蝕>中度侵蝕>強烈侵蝕>極強烈侵蝕。極強烈侵蝕面積減少雖小(5.27 km2)，但減幅較大(48.75%)，強烈侵蝕面積減少最大(22.45 km2)，減幅也最大(74.26%)。一方面由于陡坡耕地的減少，另一方面也得益于梯田的修建以及高植被覆蓋度林草地面積的增加。

研究期內(nèi)，王益區(qū)土壤侵蝕量減少幅度明顯(39.36%)，由53.81×104t減少至32.63×104t。年均土壤侵蝕模數(shù)由0.332×104t/(km2·a)變化為0.201×104t/(km2·a)。與此同時，王益區(qū)土壤侵蝕強度風險指數(shù)(SEDI)由2000年的310減少至198，區(qū)域內(nèi)的土壤侵蝕狀況整體呈好轉(zhuǎn)趨勢。

3.3.2 土壤侵蝕空間變化特征

通過全局空間自相關(guān)(Global Moran’ I)分析不同土壤侵蝕強度之間的平均關(guān)聯(lián)程度和顯著性，結(jié)果顯示2000年和2019年王益區(qū)的全局Moran’ I指數(shù)分別為0.53，0.51，相應(yīng)的Z得分分別為20.78、22.43，且通過了5%的顯著性水平檢驗，說明王益區(qū)土壤侵蝕存在明顯的空間正相關(guān)性。

基于ArcGIS的聚類和異常值分析(Anselin Local Moran I)計算得到王益區(qū)土壤侵蝕強度高值HH(High-High)、低值LL(Low-Low)的空間聚類和空間異常值HL(High-Low)，LH(Low-High)的分布圖(LISA圖)。其能夠直觀反映不同土壤侵蝕強度的空間聚集及相異性特征，彌補全局空間自相關(guān)在空間聚集位置表達上的不足。結(jié)果(圖6)顯示，2000年和2019年的土壤侵蝕強度空間聚集狀態(tài)明顯，2000年的HH聚集區(qū)面積為32.85 km2，占全部聚集區(qū)面積的54.32%，在王益區(qū)內(nèi)呈連片集中分布。到2019年LL聚集區(qū)面積增幅明顯(49.60%)，其增幅與研究期間微度侵蝕面積增幅(51.86%)相近。與此同時HH聚集區(qū)面積減少了8.4 km2，減幅高達25.57%。區(qū)域內(nèi)LL聚集程度的增強和HH聚集程度的降低與全區(qū)侵蝕強度變化趨勢相符。

總體上，研究期間土壤侵蝕強度的空間演化較為明顯，LL聚集區(qū)面積的增加與HH聚集區(qū)面積減少，也反映土壤侵蝕狀況得到了有效緩解。

疊加分析研究區(qū)侵蝕熱點分布與坡度分級。由表2可知，2000年和2019年土壤侵蝕熱點主要集中在(15°，25°]區(qū)間內(nèi)，分別約占熱點區(qū)總面積的39.95%和40.34%；其次為(8°，15°]和>25°區(qū)域。此外，(8°，15°]和(15°，25°]區(qū)域內(nèi)的熱點區(qū)面積減少較明顯且極為接近。這與前述林草地主要分布在>8°區(qū)域且(15°，25°]坡度區(qū)域內(nèi)的優(yōu)勢土地利用類型轉(zhuǎn)化為林草地關(guān)系密切。

表2 不同坡度區(qū)間土壤侵蝕熱點面積分布

研究期間侵蝕熱點區(qū)面積大幅減少，土壤侵蝕強度風險指數(shù)降低，進一步說明王益區(qū)土壤侵蝕狀況明顯好轉(zhuǎn)。究其原因，一方面由于王益區(qū)封山禁牧政策的實施；另一方面，是由于區(qū)域塬面澇池的修建以及部分坡改梯工程的實施，對區(qū)域內(nèi)土壤侵蝕防治起到一定的促進作用。

3.3 討 論

(1)采用《土壤侵蝕分級分類標準》(SL 190—2007)進行土壤侵蝕量的估算，可以基于土地利用、坡度、植被覆蓋度快速地獲取具有一定認可程度的區(qū)域土壤侵蝕強度變化，在土壤侵蝕相關(guān)研究[4，30，36-37]中應(yīng)用較多。探索性空間數(shù)據(jù)分析方法在土壤侵蝕空間聚集的分布演化方面更具優(yōu)勢，但不能反映土壤侵蝕驅(qū)動因素，后續(xù)研究中應(yīng)綜合考慮降雨、坡長、水土保持措施等多個影響因子在內(nèi)的土壤侵蝕變化情況，進一步分析各個影響因子與研究區(qū)域土壤侵蝕的相關(guān)性，或結(jié)合地理探測器等方法進一步揭示研究區(qū)土壤侵蝕可能的驅(qū)動因素。

(2)已有研究發(fā)現(xiàn)[38]，黃土塬區(qū)耕地面積減少主要由于建設(shè)用地的增加以及其他諸多因素的影響；這與研究期間王益區(qū)土地利用變化趨勢相符，此外還得益于退耕還林工作取得的進展。本文估算的2000年、2019年王益區(qū)侵蝕模數(shù)分別為0.332×104t/(km2·a)和0.201×104t/(km2·a)，與郭力宇等[4]估算的渭北高原溝壑區(qū)年平均土壤侵蝕模數(shù)在0.3×104t/(km2·a)上下浮動較為接近。王益區(qū)溝壑面積大，溝頭延伸，溝岸擴張問題嚴重，一方面可以進一步推進塬面保護工程，溝坡排水工程的實施，栽植水保林和經(jīng)濟林，防治溝道侵蝕、溝頭前進；另一方面，結(jié)合當?shù)胤馍浇琳?，保護好現(xiàn)有耕地及林草地，避免出現(xiàn)過度開墾和濫牧行為，完善水土保持治理體系，實現(xiàn)蓄水保土，固溝保塬和區(qū)域生態(tài)良性循環(huán)。

(3)針對土壤侵蝕較為嚴重的侵蝕溝，積極開展區(qū)域內(nèi)小尺度的土壤侵蝕模擬研究，一方面可評判固溝保塬等水土保持措施的效益，為其提供科學合理的改善建議；另一方面推進溝頭前進主要致災(zāi)因子的研究，制訂更加科學合理的土壤侵蝕防治體系。

4 結(jié) 論

本文基于植被覆蓋度及土地利用等研究了王益區(qū)土壤侵蝕變化，結(jié)論表明：

(1)王益區(qū)耕地面積大幅減少，耕地面積減少了27.75 km2，林草地和城鄉(xiāng)居民用地增長明顯，林草地面積增長最為顯著(17.06 km2)?？臻g上，林草地面積大幅增加，取代耕地成為(15°，25°]區(qū)域內(nèi)的優(yōu)勢土地利用類型。

(2)王益區(qū)植被覆蓋度明顯好轉(zhuǎn)。植被覆蓋度>75%的林草地面積增加了0.7倍，>60%植被覆蓋度的林草地面積增幅高達42.94%，封山禁牧的實施取得了一定的效果，對早期退耕還林工作成果保護較好。

(3)王益區(qū)土壤侵蝕狀況明顯好轉(zhuǎn)。侵蝕強度結(jié)構(gòu)由中度侵蝕和強烈侵蝕為主轉(zhuǎn)化為中度侵蝕和輕度侵蝕為主，2000年和2019年王益區(qū)土壤侵蝕面積分別為115.04 km2和90.58 km2，占全區(qū)總面積的70.92%和55.84%。SEDI由310減少到198，土壤侵蝕量減少了21.18×104t，年均土壤侵蝕模數(shù)由0.332×104t/(km2·a)變?yōu)?.201×104t/(km2·a)。

(4)王益區(qū)土壤侵蝕強度分布呈明顯的空間聚集效應(yīng)。全局Moran’ I指數(shù)由0.53變?yōu)?.51，侵蝕熱點由大面積連續(xù)分布轉(zhuǎn)化為局部集中的零散分布，其空間變化與陡坡林草地變化關(guān)系密切，這也表明退耕還林工作取得了很好的效果。