湖南蒸水流域近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值数据，感谢西北农林科技大学张宏鸣教授对地形因子提取的指导。