2.1 模型选择
美国农业部于1997年在通用水土流失方程USLE基础上改进的修正通用土壤流失方程RSULE模型,因其需要参数少,形式简单灵活,参数获取方便并且评价土壤侵蚀的各种影响因素较为全面而得到世界各地的普遍使用。与USLE模型相比,RUSLE模型在众多方面得到了改进,例如:填补了之前数据的空白,数据来源不再单一,能够模拟不同的系统等[13]。RUSLE模型在我国西南山区已有较长时间应用,适用性已得到有效验证[14],因此本文选取RUSLE模型进行研究区土壤侵蚀的量化,公式为:
A=R·K·L·S·C·P(1)
式中:A为由于降雨及其径流作用于坡面使其出现细沟或细沟间发生侵蚀所产生的土壤年平均侵蚀量[t/(hm2·a)]; R为降雨侵蚀力因子[(MJ·mm)/(hm2·h·a)]; K为土壤可蚀性因子[(t·hm2·h)/(hm2·MJ·mm)]; LS为坡长和坡度因子,无量纲; C为植被覆盖和管理因子,无量纲; P为水土保持措施因子,无量纲?厱
2.2 数据源及预处理
论文分析所用数据包括1996年、2006年、2012年、2018年日降雨量数据(来源于云南省气象台); 土壤属性和土壤类型数据通过HWSD查询砂粒含量、粉粒含量、黏粒含量和有机碳含量得到; 由地理空间数据云平台下载得到DEM数据和贡山县1996年、2006年、2012年、2018年4个时期的Landsat TM遥感影像,对获取的DEM数据根据贡山县矢量边界进行裁剪,Landsat TM遥感影像空间分辨率为30 m×30 m,影像时间为10—12月,每景影像的含云量均小于5%,数据质量较好,并利用ENVI 5.3软件对遥感影像进行辐射定标、FLAASH大气校正、图像增强、图像镶嵌和裁剪等处理。NDVI由遥感影像反演获得,土地利用类型数据由遥感影像目视解译提取得到。
2.3 RUSLE各因子的确定
2.3.1 降雨侵蚀力因子
降雨产生土壤侵蚀主要是通过雨滴溅落对土壤颗粒的冲击和地表径流的搬运作用而产生的,它是引起土壤侵蚀最重要的外部驱动因素。降雨侵蚀力因子(R)是评估区域降雨引起土壤分离和搬运的动力指标[15]。章文波等[16]通过大量的对比分析发现基于日雨量计算多年平均降雨侵蚀力的精度最高。因此,本文根据日降雨量数据采用半月逐日雨量模型计算降雨侵蚀力,模型表达式如下:
Rj=∑ki=1(Pi)β(2)
β=0.8363+(18.144)/(Pd12)+(24.455)/(Py12)(3)
α=21.586β-7.1891(4)
式中:Rj为第j个半月时段的降雨侵蚀力值[(MJ·mm)/(hm2·h·a)]; Pi为半月时段内第i天的侵蚀性日降雨量(mm),模型中要求Pi≥12,否则Pi值计入0; K为该半月时段内的天数; α,β均为模型待定参数。Pd12为日降雨量≥12 mm时的日均降雨量; Py12为日降雨量≥12 mm时的年均降雨量。虎雄岗等[17]通过对5种空间数据插值方法进行对比与分析发现,协同克里金法考虑了地形高程对降雨量的影响,因而更适合贡山县高山峡谷地区降雨数据的空间插值,具有更好的插值效果。因此,本文将各雨量站点收集到的贡山县日降雨数据代入公式(2)计算出各站点的降雨侵蚀力R,然后使用协同克里金插值法得到贡山县1996年、2006年、2012年、2018年4个时期降水侵蚀力R的空间分布图(图1)。
图1 1996-2018年贡山县降雨侵蚀力因子空间分布
2.3.2 土壤可蚀性因子
土壤可侵蚀因子(K)可以反映土壤被降雨侵蚀力剥蚀和搬运作用的敏感程度,也可以评价土壤是否易受侵蚀营力的破坏,是土壤性质中的一个重要因子[18]。土壤可侵蚀因子K值越大,表明土壤更易受到侵蚀[19]。K值的估算方法较多,最常用的有通过天然小区直接测定、作物模型(EPIC模型)、中值粒径法。本文因研究区域较小,且已获取到研究区的砂砾、粉粒、黏粒和有机碳的含量,故采用被广泛用来进行K值计算的Williams等[20]在EPIC模型中的估算方法,该方法由蔡崇法等[21]在小流域进行应用研究,其公式为:
K={0.2+0.3exp[-0.0256Wd(1-(Wi)/(100))]}×((Wi)/(Wi+Wt))0.3×
[1-(0.25Wc)/(Wc+exp(3.72-2.95Wc))]×
[1-(0.7Wn)/(Wn+exp(-5.51+22.9Wn))](5)
式中:Wd为土壤砂粒含量(%); Wi为土壤粉粒含量(%); Wt为土壤黏粒含量(%); Wc为土壤有机碳含量(%),其中:
Wn=1-(Wd)/(100)(6)
通过公式(5)—(6)计算不同土壤类型的K值,其单位为(t·hm2·h)/(MJ·hm2·mm),然后利用GIS空间属性赋值功能给各种土壤类型赋予相应的K值,得到研究区土壤可蚀性因子K的空间分布图(图2A)。
2.3.3 坡长坡度因子
地形因子(LS)包含斜坡长度因子L和坡度因子S,是影响土壤侵蚀强度的重要因素,其中坡度的影响最大。我国西南土石山区山高坡陡、地形复杂,大部分区域存在25°以上的坡地,甚至有些坡耕地的坡度已达到45°[22-23]。因此根据贡山县坡度大于25°的区域占比面积较大的特点,借鉴刘斌涛等[24]提出的适用于西南土石山区土壤流失方程坡度因子计算公式:
S={10.8sinθ+0.03 θ≤5°
16.8sinθ-0.50 5°<θ≤10°
20.204sinθ-1.2404 10°<θ≤25°
29.585sinθ-5.6079 θ>25°(7)
式中:S为坡度因子; θ为坡度(°)。
采用Wischmeier等[25]提出的经验公式估算坡长因子L值,计算公式为:
L=(λ/(22.13))m(8)
式中:λ为水平投影坡长; 22.13为标准小区的坡长; m为可变的坡长因子指数,不同的坡度选取不同的m值,其中θ≤1°时m=0.2,1°<θ≤3°时m=0.3,3°<θ≤5°时m=0.4; θ>5°时m=0.5。由公式(7)—(8)得到LS分布图(图2B)。
2.3.4 植被覆盖与管理因子
植被覆盖与管理因子(C)是指在一定植被覆盖和管理措施条件下土壤流失量与同等条件下适时翻耕、连续休闲对照地上土壤流失量之比[26],属于无量纲数,其值为0~1。植物根系对土壤有固着作用,因此,植被覆盖度越高的地区,越不容易发生土壤侵蚀[27]。目前主要有3种方法来确定C值:小区试验法、经验法和基于植被覆盖度的遥感数据定量估算法。本文采用蔡崇法等[21]的算法来计算植被覆盖与管理因子C,其计算公式如下:
C={1
0.6508-0.3436lgf
0 f=0
0<f≤78.3%
f>78.3%(9)
f=(NDVI-NDVImin)/(NDVImax-NDVImin)(10)
式中:C为植被覆盖因子; f为植被覆盖度; NDVI为归一化植被指数。
2.3.5 水土保持措施因子
水土保持措施因子(P)是指在一定水土保持措施的作用下产生的土壤侵蚀量与对应未采取保持措施的土壤侵蚀量之比,其值为0~1,0代表基本不发生土壤侵蚀; 1表示未采取水土保持措施或措施完全失效,土壤侵蚀十分严重。本文基于彭建[28]、张岩[29]等对山区土壤侵蚀研究基础上,参照其他学者的研究成果[30-31],对研究区不同土地利用类型的水土保持措施因子P进行赋值(表1)。在ArcGIS中计算得到研究区CP因子的空间分布图(图3)。
图3 1996-2018年贡山县植被覆盖与水土保持措施因子(CP)空间分布