基于GIS的小流域尺度水土流失敏感性评价——以通渭县牛谷河项目区为例

徐剑春，王 博，罗进选

(甘肃省水土保持科学研究所，甘肃 兰州 730020)

水土流失是当今全球共同面临的一个严峻的环境和灾害问题，严重的水土流失会造成生态环境急剧恶化，导致土壤表层营养成分流失、土地生产力下降等，同时农业用地的化肥、农药残留物随水土流失进入河湖水库，也是水环境恶化的主要污染源之一[1]。水土流失敏感性是指区域生态系统水土流失生态过程发生的潜在可能性及其程度，是评价区域生态环境质量、人口负荷、土地利用合理程度的指标之一，也是实施区域生态环境治理与管理的重要基础依据[2]。关于水土流失敏感性的研究成果很多，各类评价方法也较为成熟，但对小尺度、生态治理前后对比的水土流失敏感性评价的研究尚不多见。本研究以甘肃省通渭县牛谷河流域为研究区。该流域在2008—2013年实施了国家水土保持重点建设工程，至2018年底各项治理措施均已稳定发挥效益。本研究采用遥感数据结合实地测量的方法，以降雨侵蚀力、土壤可蚀性、地形起伏度和地表植被覆盖等为主要影响因子，开展小流域尺度治理前后的水土流失敏感性评价，探索水土流失敏感性监测体系及方法，希望能为后续甘肃省水土保持生态工程的监测、规划提供技术参考。

1 研究区概况

1.1 区域概况

甘肃省是全国水土流失最严重的省份之一，水土流失面积大、范围广、类型多、强度高、危害重，严重制约着区域经济社会发展。甘肃省通渭县牛谷河流域属于国家水土保持重点治理区，位于渭河一级支流牛谷河流域上中游、黄土丘陵沟壑区第三副区，地理位置介于104°54′41″～105°17′20″E、35°10′08″～35°23′34″N，海拔1 774～2 521 m，最大相对高差747 m；区内包含马营、水岔、李家大河、朱家营滩、段家峡、万家岔和蒋家川等7条完整的小流域，共涉及马营、华岭、北城、平襄、三铺和陇阳6个乡镇42个行政村；地貌特点以黄土梁峁沟壑为主，沟壑密度为1.43 km/km2，沟壑面积占项目区总面积的9.0%；干沟长37.5 km，主支沟长64.6 km，其他支毛沟长491.91 km，干沟平均比降2.6%，主支沟平均比降3.7%，干沟、主支沟多为U形沟，小支毛沟多为V形沟；属中温带大陆性气候区，年均气温6.6 ℃，≥10 ℃年活动积温1 226 ℃，年均太阳辐射总量129 kJ/(cm2·a)，年均日照时数2 239 h，年均无霜期131 d，年均大风日数32 d，人口密度166人/km2，人均耕地面积0.37 hm2。

1.2 建设概况

研究区在2008—2013年治理期内共完成各类水土保持措施面积148 km2，其中修建梯田4 007 hm2、栽植水土保持林2 085 hm2、种草1 664 hm2、实施封禁治理7 017 hm2，新建土谷坊210座、水窖96眼，实施沟头防护4.24 km。措施布设主要是在坡度<15°且距村庄较近的坡耕地建设水平梯田；在坡度15°～25°且距村庄较远、立地条件和水分条件较好的坡耕地退耕种植紫花苜蓿；在坡度>25°的陡坡耕地及15°～25°植被稀疏的荒山荒沟进行水平阶整地，种植侧柏及沙棘；在陡坡耕地及荒坡进行鱼鳞坑整地，种植山杏、柠条、沙棘；在坡度>35°且水土流失严重的荒坡荒沟和不宜布设其他防治措施的荒地采用封禁措施；此外，根据治理区域具体地形条件和沟道特征选择适宜的地点，全面建设谷坊及沟头防护措施。

2 评价方法

2.1 评价因子及模型

根据土壤侵蚀发生的动力条件，甘肃省水土流失类型主要为水力侵蚀和风力侵蚀，局部地区有冻融侵蚀及重力侵蚀。本研究主要是对以水力侵蚀为主的水土流失敏感性进行评价。采用通用土壤流失方程(USLE)为理论指导，选取降雨侵蚀力、土壤可蚀性、地形起伏度和地表植被覆盖等因子，将反映各因子对水土流失敏感性的单因子评价数据用ArcGIS转化为栅格数据后，再利用栅格计算器运算得到水土流失敏感性指数。水土流失敏感性指数评价模型为

(1)

式中：SSi为第i个空间单元的水土流失敏感性指数；Ri为第i个空间单元的降雨侵蚀力因子指数；Ki为第i个空间单元的土壤可蚀性因子指数；LSi为第i个空间单元的地形起伏度因子指数；Ci为第i个空间单元的地表植被覆盖因子指数。

2.2 数据来源

根据上述评价模型，收集本研究水土流失敏感性评价所需的气象、土壤、高程和遥感影像等，并对遥感数据尺度较大的采用现场实测数据加密后进行插值计算。

(1)降雨侵蚀力因子(R)。选取研究区附近的2个气象观测站——华家岭气象站(纬度35.38°、经度104.83°、海拔2 450.6 m)、通渭气象站(纬度35.22°、经度105.23°、海拔1 768.2 m)降雨系列资料中的年均降雨量、年均最大60 min降雨量、年均最大24 h降雨量，采用式(2)计算R因子值，利用ArcGIS软件插值并绘制研究区R因子分布栅格图(图1)。R值计算公式[3]为

图1 研究区R因子分布

(2)

(2)土壤可蚀性因子(K)。水土流失发生的主体是土壤甚至是其母质，所以土壤类型从根本上决定了土壤可蚀性[4]。K因子是指土壤颗粒被水力分离和搬运的难易程度，主要与土壤质地、有机质含量、土体结构、渗透性等土壤理化性质有关。采用应用比较广泛的EPIC模型[5]计算得出修正前的土壤可蚀性因子(KEPIC)，计算结果采用张利科等[6]的研究成果进行修正，并乘以0.131 7转化为国际单位制[7]，其计算公式分别为

KEPIC={0.2+0.3exp[-0.025 6ms(1-msilt/100)]}×

[msilt/(mc+msilt)]0.3×

{1-0.25orgC/[orgC+exp(3.72-2.95orgC)]}×

{1-0.7(1-ms/100)/[1-ms/100+

exp(-5.51+22.9-22.9ms)]}

(3)

K=(-0.013 83+0.515 75KEPIC)×0.131 7

(4)

上式中：K为修正后的土壤可蚀性因子,t·hm2·h/(hm2·MJ·mm)；KEPIC为修正前的土壤可蚀性因子，t·hm2·h/(hm2·MJ·mm)；mc为黏粒(粒径<0.002 mm)所占的百分比，%；msilt为粉粒(粒径0.002～0.05 mm)所占的百分比，%；ms为砂粒(粒径0.05～2 mm)所占的百分比，%；orgC为有机碳的百分比含量，%。

根据土地利用类型、地形、地表覆盖情况等分别进行取样，取样深度10 cm，检测、计算获取mc、msilt、ms、orgC数值，在Excel表格中利用式(3)和(4)计算K值，并根据现场调查绘制区域土壤类型图，最后以土壤类型图为工作底图，在ArcGIS中将K值依据取样经纬度插入后进行插值计算，转换形成研究区K因子分布栅格图(图2)。

图2 研究区K因子分布

(3)地形起伏度因子(LS)。从土壤侵蚀与地形坡度的关系看，基本表现为随地形坡度增大土壤侵蚀加剧的特征[8]。LS因子主要利用反映地形对土壤侵蚀影响的坡长因子(L)和坡度因子(S)进行计算，计算参数可通过研究区DEM数据获取。本研究中利用地面一定距离范围内最大地形高差，作为区域地形起伏度因子。利用研究区治理前原始地貌DEM数据，选择高程数据集，利用ArcGIS设置各像元的最大值和最小值，即得到高程数据集的最大值和最小值，在栅格计算器中利用公式(最大值-最小值)计算得出治理前(2008年，下同)LS因子栅格图(图3)。研究区主要的水土保持措施为梯田工程，其对于地形起伏度因子影响较大。通过治理后无人机测绘高程数据或实地GPS调查高程数据，形成措施实施后的DEM数据，并计算得出治理后(2018年，下同)的LS因子栅格图(图4)。

图3 治理前LS因子分布

图4 治理后LS因子分布

(4)地表植被覆盖因子(C)。植被覆盖度对生态敏感性具有重要影响，一般来说，植被覆盖度越高，植物层次就越丰富，抗破坏能力就越强，生态敏感性就会降低[9]。在对光谱信号进行分析的基础上，通过建立归一化植被指数与植被覆盖度的转换信息，直接提取植被覆盖信息。C因子计算公式为

C=(NDVI-NDVIsoil)/(NDVIveg-NDVIsoil)

(5)

式中：NDVIveg为植被完全覆盖地表所贡献的信息；NDVIsoil为无植被覆盖地表所贡献的信息。

本研究采用覆盖全国的MODIS NDVI数据，来源于美国国家航空航天局的EOS/MODIS数据产品(http://e4ft101.cr.usgs.gov)，空间分辨率为250 m×250 m，时间分辨率为16 d。由于大部分植被覆盖类型是不同植被类型的混合体，所以不能采用固定的NDVIsoil和NDVIveg，通常根据NDVI的频率统计表，计算NDVI的频率累积值，即累积频率为2%的NDVI值为NDVIsoil，累积频率为98%的NDVI值为NDVIveg，再使用栅格计算器依据公式(5)，计算植被覆盖因子。根据治理前后的植被NDVI影像数据，形成治理前C因子分布栅格图(图5)及治理后C因子分布栅格图(图6)。

图5 治理前C因子分布

图6 治理后C因子分布

2.3 水土流失敏感性评价

(1)敏感性等级值。综合采用自然分界法与相关文献确定分级赋值标准参考值，根据研究区4个影响因子的分布区间，以突出区域特点和评价区域治理前后对比性为出发点，修正后确定各类影响因子的分级赋值标准和不同评价指标对应的敏感性等级值，见表1。

表1 水土流失敏感性的评价指标及分级

(2)模型运算。将4类评价因子栅格图统一形成5 m分辨率的栅格数据，并按照表1进行敏感性指数分级赋值，以赋值结果在ArcGIS中利用栅格计算器根据式(1)计算得到水土流失敏感性指数。研究区治理前后水土流失敏感性指数分布见图7、8。

图7 治理前水土流失敏感性指数分布

图8 治理后水土流失敏感性指数分布

3 评价结果

研究区治理措施实施前后评价因子中R、K因子变化不大，因此本研究主要分析研究区LS、C因子及水土流失敏感性指数的变化情况。利用计算得出的R、K、LS、C因子值，按照表1中的敏感性指数分级赋值表进行赋值，依据面积加权计算各因子敏感性指数，再利用4类评价因子赋值后的敏感性指数，根据式(1)计算得到研究区治理前后的水土流失敏感性指数，见表2。从表2可以看出，水土流失综合治理工程实施后通渭县牛谷河流域水土流失敏感性指数从2.285下降至2.226，LS因子敏感性指数从2.968下降至2.854，C因子敏感性指数从2.035下降至1.831，说明研究区内重度及中度敏感区面积均有显著减少。

表2 研究区治理前后LS、C因子和水土流失敏感性变化分析

4 结 语

基于RS、GIS、物联网监测技术进行生态项目监测，是水土保持监测行业未来发展的重要方向。目前，个别区域的小流域尺度的水土保持监测要么流于形式，要么主要依托国家重点监测项目，采用的监测方法较为传统，以调查访问、 资料收集和样地定位监测为主，难以准确反映治理前后的水土保持效益。因此，探索小流域尺度的GIS监测评价体系及方法，研究现场实测结合遥感分析的数据获取及处理技术，具有很强的现实指导意义。本研究尝试利用GIS技术针对具体治理项目开展了小流域治理前后的水土流失敏感性评价，尚存在遥感数据精度不高、现场气象监测等数据来源较少、通用计算模型与区域实际不完全匹配等一些问题，还需进一步展开相关研究，逐步完善水土流失敏感性评价的监测体系和方法。