湖南蒸水流域近20年土壤侵蝕時空變化

王海濤， 謝紅霞?， 楊勤科， 周 清， 段良霞

(1.湖南農業大學資源環境學院，410128，長沙；2.西北大學，710127，西安)

土壤侵蝕破壞地表植被，造成土壤養分流失，降低土地生產力，導致土地資源破壞、流域水質污染、水庫泥沙堆積和洪澇災害，危及生態安全，是全球性的生態環境問題之一[1- 3]。USLE和RUSLE模型的大力發展[4]，廣泛應用于土壤侵蝕研究，多在坡面和小流域尺度上進行。當前國內學者對于流域土壤侵蝕時空變化特征以本地小流域的研究較多，且大多局限于短時間尺度的預測，對于大型流域長時間尺度的研究較少。蒸水流域曾經是湖南省水土流失嚴重的地區之一，流域內土壤母質風化程度高，防洪、排澇設施少，洪水來臨導致河水漫灘，造成較大財產損失[5]。為了遏制水土流失，當地相關部門在蒸水流域開展了水土保持和水環境綜合治理工作，如蒸水南堤草橋- 蒸水橋防洪工程建設等。要了解流域土壤侵蝕歷史及現狀，需要在蒸水流域土壤侵蝕狀況進行定量評價與分析。筆者基于土壤侵蝕的基本理論，利用GIS和RS技術，將CSLE模型與長時序地面觀測數據結合，定量評價蒸水流域1995、2000、2005、2010和2015年的土壤侵蝕情況，了解流域內土壤侵蝕環境背景和土壤侵蝕差異，掌握其水土流失整體現狀，研究結果可望為水土流失防治、水土保持規劃和技術管理提供科學參考，以便更好的利用水土資源。

1 研究區概況

蒸水屬洞庭湖水系湘江的一級支流，流域介于E 111°53′～112°37′、N 26°52′～27°10′之間，流域面積3 480 km2，河長194 km，形狀呈“乙”字型，河流坡降0.54‰。北、東、南3面環山，山地丘陵為2 720 km2，約占流域總面積的71%；中、下游地區為平原、低地，面積為750 km2，約占總面積的20%[6]；其余為水域。流域內及流域周邊有雙峰站、邵東站等9個氣象站，神山頭站為水文站(圖1)。

圖1 蒸水流域高程圖Fig.1 Elevation map of Zhengshui River Basin

2 數據與方法

研究采用中國土壤流失方程CSLE模型來估算蒸水流域的土壤侵蝕量，CSLE是Liu等[7]基于我國土壤侵蝕現狀，以USLE和RUSLE為原型，改進提出的適用于我國的預測坡耕地年土壤流失的經驗模型。該模型結構簡單、參數要求低，也適用于溝壑縱橫，山高坡陡地區。模型為

A=RKLSBET。

(1)

式中：A為土壤侵蝕模數，表示多年平均土壤流失速率，t/(hm2·a)；R為降雨侵蝕力因子，MJ·mm/(hm2·h·a)；K為土壤可蝕性因子，t·hm2·h/(hm2·MJ·mm)；L為坡長因子；S為坡度因子；B為植被覆蓋與生物措施因子；E為工程措施因子；T為耕作措施因子。參考水利部辦公廳文件辦水保(2018)189號文件《關于印發區域水土流失動態監測技術規定(試行)的通知》及水利部監測中心頒布的技術規定[8]：R因子采用殷水清等[9]的利用冷暖季日雨量資料計算半月降雨侵蝕力公式，累加半月降雨侵蝕力得到年降雨侵蝕力，經過普通克里金空間插值后，形成蒸水流域降雨侵蝕力因子分布圖。K因子數據由劉寶元教授提供，計算方法采用的國務院第1次全國水利普查計算K值方法[10]。LS因子采用Liu等[7]的CSLE算法。

參照《區域水土流失動態監測技術規定(試行)》附錄7的要求，并結合前人對B因子的研究和對華中區林下蓋度的野外調查成果[11- 13]，得到不同用地類型的B因子值，其中林地和草地的B因子利用TM多光譜影像(1995—2010年利用Landsat4-5 TM衛星影像，2015年利用Landsat 8 OLI_TIRS衛星影像)中紅波段和近紅外波段，計算歸一化植被指數( normalized difference vegetation index， NDVI)，基于NDVI數據計算得到30 m空間分辨率的植被覆蓋度(fractional vegetation cover，FVC)。

灌木林地和草地的B值采用計算公式：

(2)

(3)

果園、其他園地、有林地和其他林地的B值采用計算公式：

B=0.444 68e(-3.200 96GD)- 0.040 99e(FVC-FVC×GD)+0.025。

(4)

式中：GD為喬木林的林下蓋度。根據前人對華中區林下蓋度的研究和野外調查[12- 13]，分別擬定有林地、疏林地和其他林地的林下蓋度值，代入到上式計算。耕地、城鎮建設用地、農村建設用地、其他建設用地、水域、裸土地和裸巖的B值分別直接賦值為1.000、0.010、0.025、0.100、0、1.000和0[8]，形成蒸水流域植被覆蓋與生物措施因子圖。

蒸水流域內工程措施主要為土坎水平梯田，參照《區域水土流失動態監測技術規定(試行)》，將水田賦E因子值為0.084，其余用地類型賦值為1.000。蒸水流域所在區域在全國輪作區內均處于2湖平原丘陵水田中三熟二熟區，將耕地T因子賦值為0.312，其他用地類型賦值為1。形成蒸水流域工程措施因子圖和耕作措施因子圖。

由于土壤侵蝕變化受到了自然和人類活動的綜合作用影響，其中自然因素主要考慮于降水變化，人類活動主要考慮土地利用變化和水土保持措施實施的變化。為了區分降水和人類活動對蒸水流域土壤侵蝕的不同影響，利用控制變量方法，分別模擬只考慮降水變化和只考慮土地利用變化情景下2015年的土壤侵蝕狀況[14]。具體步驟如下：1)不考慮其他因素變化，只考慮降水變化情景下的2015年土壤侵蝕模數：A2015R=R2015KLSB1995E1995T1995；2)不考慮其他因素變化，只考慮土地利用變化情景下的2015年土壤侵蝕模數：A2015=R1995KLSB2015E2015T2015；3)同時考慮降水變化和土地利用變化情景下的2015年土壤侵蝕模數：A2015RBET=R2015KLSB2015E2015T2015；4)將計算的1995年實際土壤侵蝕模數與模擬情景下的2015年土壤侵蝕模數進行對比。

3 結果與分析

3.1 土壤侵蝕因子

基于CSLE的因子算法，利用1995—2015年的5個年度數據，獲得了蒸水流域30 m分辨率R、K、LS、B、E、T因子圖。

3.1.1 降雨侵蝕力因子 統計蒸水流域周邊9個氣象站每年侵蝕性降雨的出現時間，繪制柱狀圖(圖2)。1995—2015年各年的平均侵蝕性降雨時間分別為32.6、44.6、44.1、47.9和48.3 d；非侵蝕性降雨時間分別為332.4、321.4、320.9、317.1和316.7 d。1995年侵蝕性降雨占全年比例最小，而2010和2015年侵蝕性降雨占全年比例最大。

圖2 蒸水流域各氣象站5個年度各級侵蝕性降雨全年出現次數Fig.2 Frequency of different erosive rainfall of 5 periods in each climate station of Zhengshui River Basin

如表1與圖3所示，2010年R值最大，均值為5 934 MJ·mm/(hm2·h·a)，地區間年降雨侵蝕力分布差異最大，標準差為1 280 MJ·mm/(hm2·h·a)；1995年R值最小，均值4 275 MJ·mm/(hm2·h·a)，標準差為215 MJ·mm/(hm2·h·a)；2005年R值為5 177 MJ·mm/(hm2·h·a)，降雨侵蝕力空間分布最為均勻，標準差為199 MJ·mm/(hm2·h·a)。

表1 蒸水流域降雨侵蝕力統計特征值Tab.1 Statistics of rainfall erosivity in Zhengshui River Basin MJ·mm/(hm2·h·a)

圖3 蒸水流域5個年度降雨侵蝕力因子圖Fig.3 Rainfall erosivity factor map of 5 periods in Zhengshui River Basin

3.1.2 土壤可蝕性因子 蒸水流域表層土壤可蝕性因子K值主要介于0～0.006 748 t·hm2·h/(hm2·MJ·mm)之間，均值為0.004 333 t·hm2·h/(hm2·MJ·mm)。高值區主要分布在流域南部、東南部以及中部區域，低值區主要分布在流域西北部地區。造成這種空間異質性的原因主要取決于土壤類型的異質性； 高值區的主要土壤類型為水稻土，低值區的主要土壤類型為石灰土(圖4)。

圖4 蒸水流域土壤可蝕性因子K值圖Fig.4 Soil erodibility factor K map in Zhengshui River Basin

3.1.3 坡長坡度因子 利用DEM數據提取得到蒸水流域坡度分布圖和地形因子圖(圖5和6)蒸水流域LS值均值為4.45。流域東北部、北部、西南山區LS值較高，中部、南部與東南部LS值較低，呈馬蹄形分布。

3.1.4 植被覆蓋與生物措施因子和土地利用情況 1995、2000、2005、2010和2015年蒸水流域的年均B因子分別為0.451 1、0.450 7、0.449 6、0.448 5和0.445 6。B因子逐漸減小，年際差異主要來源于土地利用變化情況。西南、北部區域B值較低，東南部區域B值較高(圖7)。

圖5 蒸水流域坡度分布圖Fig.5 Slope map of Zhengshui River Basin

圖6 蒸水流域坡度坡長因子圖Fig.6 Slope and slope length factor LS map of Zhengshui River Basin

圖7 蒸水流域5個年度植被覆蓋與生物措施因子圖Fig.7 Vegetation cover and biological practice factor B map during 5 periods in Zhengshui River Basin

1995年蒸水流域內以耕地和林地為主，其面積分別占總面積的42.79%和54.16%，呈現以農、林地為主的土地利用結構特征。此后隨著退耕和建設占用，耕地面積逐年下降。到了2015年，耕地和林地為主要用地類型的格局沒有變化，分別占流域面積的42.22%和53.89%。相比于1995年，水田面積占比減少0.37%，旱地面積占比減少0.19%，林地面積比例減少0.27%，而城鎮用地、農村居民點等建設用地面積比例增加0.80%，這種變化是流域內城市化發展的結果。

3.1.5 工程措施因子與耕作措施因子 1995—2015年蒸水流域的年均E因子分別為0.710 9、0.711 5、0.711 9、0.712 1和0.714 3，西南、北部區域E值較高，中部與東南部區域E值較低。E因子的空間異質性來源于梯田，而西南及北部地區為非水田區域，工程措施因子都為1(圖8)。1995—2015年蒸水流域的年均T因子分別為0.705 6、0.706 1、0.706 8、0.707 1和0.709 5。西南、北部區域T值較高，中部與東南部區域T值較低。T因子的空間異質性來源于2湖平原丘陵輪作區的耕地，西南及北部地區用地類型為非耕地，耕作措施因子都為1(圖9)。

圖8 圖8 蒸水流域2015年工程措施因子E值圖Fig.8 Engineering practice factor E map of Zhengshui River Basin in 2015

圖9 蒸水流域2015年耕作措施因子T值圖Fig.9 Tillage practice factor T map of Zhengshui River Basin in 2015

3.2 土壤侵蝕強度

運用CSLE模型計算得到蒸水流域1995、2000、2005、2010和2015年的土壤侵蝕模數(圖10)。各年的平均土壤侵蝕模數分別為412、520、479、530和528 t/(km2·a)，1995與2005年屬于微度侵蝕等級，2000、2010和2015年屬于輕度侵蝕等級。根據水利部行業標準SL190—2007《土壤侵蝕分類分級標準》規定，蒸水流域所在的南方紅壤丘陵區土壤侵蝕模數所推薦的容許土壤流失量為500 t/(km2·a)，這5年的平均侵蝕模數為494 t/(km2·a)，低于容許土壤流失量，可以判定蒸水流域水土流失情況較輕。

利用ArcGIS軟件統計功能，得到各土壤侵蝕強度等級的侵蝕面積與比例。按照水利侵蝕強度分級對土壤侵蝕模數進行劃分，蒸水流域1995—2015年發生的主要侵蝕類型為微度侵蝕，均占到全流域面積的70%以上；輕度侵蝕區域均占全流域面積的20%以上；中度侵蝕區域占全流域面積的2%以內；強烈及以上等級侵蝕區域占全流域面積的1%以內(表2)。將土壤侵蝕強度空間分布圖(圖10)與流域坡度分布圖(圖5)比較，發現土壤侵蝕強度在空間上的分布與坡度區間在空間上的分布相似，坡度較大的區域，土壤侵蝕強度較高。

2015年與1995年相比，時間跨度較長，R值與B、E和T值差異較大，A值差異也較大。通過模擬考慮降水變化和考慮土地利用變化情景下2015年的土壤侵蝕狀況，結果發現：1)只考慮降雨變化情景下的2015年單位面積土壤侵蝕量均值為526 t/(km2·a)；2)只考慮土地利用變化情景下的2015年單位面積土壤侵蝕量均值為412 t/(km2·a) ；3)同時考慮降水和土地利用變化情景下的2015年單位面積土壤侵蝕量均值為528 t/(km2·a)。通過對比可以看出：從1995年到2015年，降水變化加劇了蒸水流域的土壤侵蝕，土地利用變化對土壤侵蝕變化不明顯，其中降水變化使流域侵蝕量增加27.67%，土地利用變化使流域侵蝕量增加0.18%，降水和土地利用變化共同作用使侵蝕量增加28.16%。

圖10 蒸水流域5個年度土壤侵蝕模數空間分布Fig.10 Spatial distribution of soil erosion modulus during 5 periods in Zhengshui River Basin

表2 蒸水流域5個年度各土壤侵蝕強度的面積Tab.2 Area of each soil erosion intensity during 5 periods in Zhengshui River Basin km2

分區統計不同用地類型的土壤侵蝕情況(表3)。流域內面積大小占比靠前的土地利用類型分別為有林地(35.94%)、水田(31.24%)、疏林地(15.20%)、旱地(11.00%)，以上4類用地占到研究區總面積的93.38%。土壤侵蝕強度從大到小分別為旱地>有林地>疏林地>水田，旱地土壤侵蝕模數大于有林地土壤侵蝕模數，主要原因是旱地較有林地的植被覆蓋度小，其B值遠大于有林地的B值；有林地土壤侵蝕模數大于疏林地土壤侵蝕模數，主要原因是有林地的平均海拔高、坡度大，其LS值遠大于疏林地的LS值；疏林地土壤侵蝕模數略大于水田土壤侵蝕模數，主要原因是水田在工程措施與耕作措施的共同影響下，保土效果略優于植被覆蓋與生物措施對疏林地的影響。

表3 蒸水流域5個年度主要用地類型與土壤侵蝕模數Tab.3 Soil erosion modulus under different land use type during 5 periods in Zhengshui River Basin

3.3 各縣區的侵蝕強度分析

根據模擬估算結果表明，土壤侵蝕現象在蒸水流域并不常見。蒸水流域涉及的行政區面積前3的區縣為衡陽縣、邵東縣和衡南縣，面積占比分別為64.98%、16.47%和11.84%，其余各縣區總面積比例不足7%(表4)。南方紅壤丘陵區所推薦的容許土壤流失量為500 t/(km2·a)，衡陽縣僅在2010年達到輕度侵蝕等級，平均侵蝕模數為462 t/(km2·a)；衡南縣在2000年、2005年和2015年達到輕度侵蝕等級，平均侵蝕模數為517 t/(km2·a)；邵東縣在5個年度均為輕度侵蝕，平均侵蝕模數為591 t/(km2·a)。流域內土壤侵蝕模數>500 t/(km2·a)的地區集中分布在邵東縣東部，衡陽縣西部、北部、東北部，衡南縣西北部，南岳區，以上地區是未來水土保持的重點地區。

表4 蒸水流域5個年度各縣區面積與土壤侵蝕模數Tab.4 Soil erosion modulus and area of each county or district under different region during 5 periods in Zhengshui River Basin

3.4 土壤侵蝕量與輸沙量對比

神頭山站是蒸水流域的水文站，其控制面積為2 857 km2，約占全流域面積的82.10%，利用ArcGIS分區統計工具來計算神頭山水文站控制面積的侵蝕模數，并將計算侵蝕量與神頭山站觀測的輸沙量進行對比(表5)。

表5 神山頭站控制范圍內外各年統計特征值Tab.5 Statistical characteristic values of each year within and outside the control area of Shenshantou Station

結果表明1995—2015年神山頭水文站控制范圍侵蝕量呈波動變化，輸沙量逐漸減少，由1995年的26.4萬t減少到2015年的10.2萬t，各年侵蝕量遠遠大于該年輸沙量，泥沙輸移比逐年下降，由0.22下降到了0.07。神山頭站控制范圍內輸沙量遠小于侵蝕量，實測輸沙量逐年降低。究其原因：一是因為本研究的計算侵蝕量是基于坡面模型的每個單元侵蝕量的總和，沒有考慮每個單元侵蝕量在輸移過程中的沉積[15]，神頭山水文站的控制面積占全流域的80%以上，侵蝕土壤需要通過長距離運輸，輸移過程中有大量侵蝕土壤發生沉積，絕大部分侵蝕土壤因就地沉積或沿途沉積等因素最終未能抵達神頭山站；二是由于流域內林草覆蓋度達到50%以上，植被覆蓋除減少坡面侵蝕且攔截侵蝕土壤的作用顯著；三是流域除梯田外還有一些其他水土保持工程措施，在本次評價中未考慮，他們也能攔截侵蝕土壤[16]。

4 討論與結論

1)蒸水流域土壤侵蝕在空間上呈現由西北至東南逐漸減弱的特征。5個年度的侵蝕模數分別為412、520、479、530和528 t/(km2·a)，呈現波動上升的趨勢。5個年度的平均侵蝕模數為494 t/(km2·a)，低于容許土壤流失量，屬于微度侵蝕等級，全流域平均水土流失情況較輕，但存在局部侵蝕嚴重的情況。

2)侵蝕性降雨發生頻次較高，侵蝕性降雨變化是蒸水流域近20年侵蝕時空變化最主要的驅動因素，降水變化使流域侵蝕量增加了26.67%。

3)蒸水流域主要用地類型的土壤侵蝕情況為旱地>有林地>疏林地>水田。

4)神山頭水文站控制范圍內侵蝕計算結果遠大于實測輸沙量，輸沙量逐年降低；流域內高覆蓋度的林草，水土保持工程參與攔截侵蝕土壤過程，使絕大部分侵蝕土壤就地或沿途發生沉積，最終未能到達水文站泥沙出口，使泥沙輸移比很小，而且逐年下降。

5)流域內土壤侵蝕模數>500 t/(km2·a)的地區集中分布在邵東縣東部，衡陽縣西部、北部、東北部，衡南縣西北部，南岳區，以上地區是未來水土保持的重點地區。

本次研究原基于研究區土壤粒徑組成以及土壤有機碳含量，采用Williams的EPIC模型法[17]對蒸水流域的土壤可蝕性K值進行計算，但直接利用EPIC公式估算K值與實測值存在較大差異，最大可達10倍，與其他學者研究結果類似[18]，因此直接使用了現有的土壤可蝕性研究成果評價了蒸水流域的土壤侵蝕。

感謝北京師范大學劉寶元教授提供的湖南省土壤可蝕性因子K值數據，感謝西北農林科技大學張宏鳴教授對地形因子提取的指導。