基于RUSLE模型的中国土壤水蚀时空规律研究

李佳蕾,孙然好,熊木齐,杨国成

1 中国科学院生态环境研究中心城市与区域生态国家重点实验室,北京 100085 2 中国科学院大学,北京 100049

土壤侵蚀由于对经济和环境带来严重的影响而成为了一个全球性的难题[1]。全球受水蚀影响的土地面积为10.94亿hm2,其中751万hm2受侵蚀严重[2]。世界上大部分农业用地的土壤正在以每年13 t/hm2到40 t/hm2的速度流失,每年约有1000万hm2的农田因水土流失而退化[3]。我国是土壤侵蚀最严重国家之一,根据第3次水土流失普查,我国年土壤流失量约50亿t[4],水蚀区平均侵蚀强度约为3800 t km-2a-1[5]。根据《第一次全国水利普查水土保持情况公报》,2011年,全国水力侵蚀面积达129.32万km2。水土流失已经成为制约我国经济社会可持续发展的重要因素之一[4]。

多数土壤侵蚀研究集中在小尺度区域,如坡地、坡面和小流域。这些研究多基于小区实验[6],也有的研究利用模型计算,Karydas整理归纳了8类、82个水蚀模型[7]。大尺度土壤侵蚀定量研究主要是通过尺度上推和宏观因子评价的方法进行[8]。定量评估土壤侵蚀的经验模型中,最常用的是美国农业部(United States Department of Agriculture, USDA) 提出的通用土壤流失方程(universal soil loss equation, USLE)和修正版的通用土壤流失方程 (revised universal soil loss equation, RUSLE)[9]。它以降雨侵蚀力、土壤可蚀性、地形、植被管理和保护措施为输入,采用基于因子计算的方法估计每年因侵蚀造成的土壤损失[10]。自中国引入USLE和RUSLE模型以来,基于此模型的土壤侵蚀的研究在小区域和小流域上有着很多的研究应用,而大流域和大尺度的研究相对较少。比如,傅伯杰等[11]利用USLE模型计算了我国黄土高原的水蚀,得出黄土高原每年水土流失量从2000年的3362 t/km2降低到了2008年的 2405 t/km2；还有其他学者利用USLE或RUSLE模型估算了中国南方丘陵山区[12]、广东省[13]、江西省[14]和全中国[15]的年平均土壤侵蚀率。虽然这些也是较大尺度上的研究,但是大部分研究将中国的各个地区的土壤侵蚀计算方法归一化,整个研究区用同种算法计算各个因子,忽略了各个因子因人类活动和地理分异带来的区域差异。

本研究基于国际上较新发表的RUSLE模型中的C因子、R因子和P因子的计算方法,结合中国的实际情况对各因子进行修正,并应用在中国土壤侵蚀的估算中。与以往的研究相比,利用大尺度高精度(30 m分辨率)数据、按区域性质对中国不同地域进行分区运算,提高了土壤侵蚀计算结果的空间异质性与精准度,并对2000年至2015年的中国土壤侵蚀量进行时空动态分析,探究中国土壤侵蚀强度变化的原因,为中国水土流失的治理与水土保持服务功能的保护提供理论基础。

1 基础数据与研究方法

1.1 基础数据

本研究所使用的主要基础数据及其来源见表1。各个数据的中国分布情况见图1。

表1 基础数据

其中降雨数据是全国839个国家级站点的日数据经SPSS软件统计得出年的各个站点的年降水量,以及各气象站点的经纬度坐标和海拔高度。DEM数据的原始数据是USGS下载的中国分区数据,再合成整个中国的30 m分辨率数据。中国年度植被指数(NDVI)空间分布数据是基于连续时间序列的SPOT/VEGETATION NDVI卫星遥感数据,采用最大值合成法生成的。土壤成分数据是在ISRIC中下载的2017年全球土壤的黏粒、粉粒、砂粒和有机碳的含量数据,再提取中国的土壤成分数据。降水数据、NDVI、土地利用类型数据都选取了2000、2005、2010、2015年的数据。

1.2 研究方法

1.2.1通用土壤流失方程(RUSLE)

RUSLE表达式为[16]：

A=R×K×LS×C×P

(1)

式中,A是模型预测的年土壤侵蚀量,单位为t km-2a-1。R是降雨侵蚀力因子,单位为MJ mm hm2h-1a-1,是降水产生的径流对土壤造成侵蚀的动力指标[17],降雨的强度和持续时间对侵蚀有重要影响。LS为地形因子,无量纲；L即坡长因子,是指标准化到22.13 m坡长上的土壤侵蚀量；S即坡度因子,是指标准化到5.14°坡度下的土壤侵蚀量[18]；一般的小尺度研究直接利用实测地形数据,大尺度研究会利用DEM数据提取信息来计算LS因子。K是土壤可蚀性因子,单位为thm2h hm-2MJ-1mm-1,其反映了土壤对侵蚀营力分离和搬运作用的敏感性[19-20]。C是植被覆盖与管理因子,无量纲,表示植被覆盖和管理措施对土壤侵蚀的作用[21]。P因子即水土保持措施因子,无量纲,指在特定水土保持措施的土壤流失与起伏地耕作的相应土壤流失之比[16]。C和P因子可反映人为控制土壤侵蚀的作用。

1.2.2基于气候分区的R因子修正

Naipal等[22]基于美国降雨站点监测值,拟合降雨侵蚀力因子R值和年降水量(P)、高程(Z)、日降水强度(SDII)的回归方程,得出修正的不同气候区降雨侵蚀力R计算方法：

R=f(P,Z,SDII)

(2)

式中,日降水强度SDII等于每日降水量≥1 mm的总降雨量与天数的比值。不同的气候区有不同的计算公式,其气候区的划分标准为Koppen-Geiger气候分区法[23]。

本研究将利用Naipal的方法,将中国不同的气候区R值分别进行计算。对于精度较差的地区,以及无SDII数据的气候区,则利用Renard and Freimund指数[16](表2)。

首先利用全国800多个站点的日降水数据计算出各个站点的年平均降雨侵蚀力,再根据李璐等基于地统计学的降雨侵蚀力插值方法的研究,选用模拟精确度较好的普通克里金插值方法来进行空间插值[24]推测全国的R值。

表2 各个气候区的R因子的计算方法

(1)气候分区的命名方式为字母缩写组合,字母缩写意义为B：干旱区Arid,C：温和区Temperate,D：寒冷区Cold,W：沙漠区Desert,S：干旱草原Steppe,k：Cold,f：Withoutdryseason,w：DryWinter,a：HotSummer,b：WarmSummer

1.2.3基于土地覆被的C因子修正方法

利用Borrelli计算C因子的方法[25],并结合中国实际土地覆盖类型和农作物种类对中国的耕地和非耕地的C因子分别进行调整与计算(C因子的取值见表3)。

表3 不同土地覆盖类型的C值

植被覆盖因子C和农作物种类紧密相关[26-27],所以本研究中耕地采用国家统计局发布的各省耕地主要作物种类和播种面积,将已发布的作物归为10类,通过以下公式计算过去全国耕地的C值：

(3)

式中,Ccropn是作物n的C值,%RegionCropn是作物n的播种面积占各省总面积的比例。

植被覆盖因子C取决于植被覆盖度,本研究结合已有文献中非耕地的各种植被覆盖类型的C值[25, 28],利用土地利用数据和NDVI数据,计算非耕地的C值。

CNonArable= Min(CNA) + Range(CNA) × (1-Fcover)

(4)

Fcover=VFC = (NDVI- NDVImin)/(NDVImax- NDVImin)

(5)

式中,Min(CNA)是CNA取值范围的最小值,Range(CNA)是CNA的最大值与最小值的差,Fcover是植被覆盖度。

1.2.4基于地形特征的P因子修正方法

在USLE/RUSLE支持下的P因子在大尺度的土壤侵蚀风险建模中很少被考虑,因为很难对大面积进行估算[29]。Xiong等依据大量文献总结[10, 30-31],归纳不同坡度耕地及不同的水土保持措施的P值差异。针对本研究,将Xiong的P因子赋值方法结合历年土地利用类型,将水平稻田区域赋值为0.2,其他农田区域根据坡度为P因子赋值。其中坡度小于等于10°的农田P取值为0.5,坡度大于10°小于等于25° 的P取值为0.6,坡度大于25°小于等于45°的P取值为0.8,坡度大于45°的P取值为1。

1.2.5其他因子的计算方法

DEM的单个栅格大小会影响计算精确度,栅格越大坡度减小,估算的土壤侵蚀量随之减小[32-33]。基于数据可靠性和计算可行性,本研究利用30 m分辨率的中国DEM数据,选择如下公式[34- 36]计算L因子：

L=(γ/22.3)m

(6)

m=β(1+β)

(7)

(8)

式中,γ为坡面的水平投影长度,为方便计算直接用栅格长度(30 m)代替；m为坡长因子指数,β为细沟侵蚀和面蚀的比值[37],θ为坡度。

对于S因子,则利用刘宝元提出的CSLE模型中计算S因子的方法[38],对不同坡度段分别进行计算：

(9)

RUSLE中的K值是对具体土壤的由试验确定的定量数值,刘宝元提出了适用于我国的土壤可蚀性测定方法,即在15°坡度、20 m坡长、清耕休闲地上,单位降雨侵蚀力所引起的土壤流失量[20]。但本研究尺度较大,因此土壤可蚀性因子K利用EPIC模型[39]计算。其公式为：

(10)

(11)

式中,SAN为砂粒含量百分数,SIL为粉粒含量百分数,CLA为黏粒含量百分数,C为有机碳含量百分数。

1.3 中国土壤侵蚀年际变化分析

将年均土壤侵蚀量的栅格图层按照时间序列进行一元线性回归分析,采用最小二乘法计算2000—2015年土壤侵蚀变化斜率,以得出2000—2015年中国土壤侵蚀强度的时间变化趋势,斜率为正则表示土壤侵蚀有增加的趋势,斜率为负则表示土壤侵蚀有减小的趋势。斜率的计算公式如下：

(12)

式中,θslope是回归方程的斜率,n为所计算的年数,Ai为第i年的平均土壤侵蚀量。变化趋势的显著性采取F检验,以检验中国土壤侵蚀的时间变化规律是否显著。该统计量的计算构造为：

(13)

2 结果

2.1 RUSLE各因子的空间分异

利用上述计算各因子的方法得出的各因子的空间分布图如图2。其中,C因子和R因子的年际变化较明显,历年的土壤侵蚀的变化即体现在降雨、植被覆盖及土地利用中。由于本研究所计算的R因子、C因子和P因子利用了较新的算法,较之前的国内RUSLE模型的应用研究有了一定程度的细化。

2000年、2005年、2010年和2015年中国水蚀区C因子的年平均计算结果为0.1464、0.1482、0.1498和0.1538,2015年较2000年提高了5.05%。而计算后的2000年、2005年、2010年和2015年中国水蚀区年平均R因子分别为2903.79、3034.62、3515.44、3444.77 MJ mm hm2h-1a-1,2015年比2000年增加了18.63%,说明中国降水带来的侵蚀潜力有增强的趋势。

图2 中国RUSLE(1)模型各因子空间分布(2)Fig.2 Distribution of RUSLE factors in China(1)RUSLE：修正版的通用土壤流失方程,revised universal soil loss equation。(2)图例中,降雨侵蚀力因子单位为MJ mm hm2 h-1 a-1,土壤可蚀性因子单位为t hm2 h hm-2MJ-1mm-1

2.2 中国土壤侵蚀的空间热点区域

在将所有的RUSLE模型中的因子利用重采样统一为100 m分辨率后,将各因子相乘计算出各年土壤侵蚀量。其中,2000年、2005年、2010年、2015年中国的年平均土壤侵蚀量分别为38.63、37.35、49.03、47.84 t km-2a-1。根据中华人民共和国水利部发布的《土壤侵蚀分类分级标准》(SL190—2007)[40],将计算后的中国土壤侵蚀度分成6个等级,即微度、轻度、中度、强烈、极强烈和剧烈,其空间分布见图3,分级统计情况如表4。

图3 中国2000—2015年土壤水力侵蚀程度空间分布Fig.3 Distribution map of water erosion in China

表4 中国2000—2015年土壤侵蚀强度统计

从表4可以看出,中国水蚀区大多(60%以上)的空间为微度侵蚀,微度侵蚀的变化幅度不大,但2010年的微度侵蚀的面积占比是4年中最小的；剧烈侵蚀的空间占比在2005年后有所增加,其中2010年剧烈侵蚀的面积比最高,比2000年增加了2.00%。2000年以来,全国土壤侵蚀总面积在减少。具体来讲,侵蚀轻微区在增加、侵蚀剧烈区也在增加,说明了生态修复的效果体现在了侵蚀总面积和侵蚀改善区的面积增加；而由于部分区域的耕地增加及降水增加,也会导致个别区域的侵蚀强度增加,是后续需要进一步加强侵蚀管理和控制的区域。

中国水蚀热点分布在中国南方,呈点状分布分散于南方各地；北方虽有水蚀较强的区域,但水蚀剧烈的地方相对较少。在南方,云贵高原与四川盆地之间的区域水蚀剧烈,尤其分布于贵州省。通过研究其地形和土地覆被,分析其原因可能是其处于喀斯特地貌区,海拔较高(1000 m以上)、坡度较大,土地覆被类型大多为耕地,人口密集、人为干扰较多,土壤保持力较为低下[41]。中国东南部的丘陵区土壤侵蚀剧烈的分布也比较广,但大多为点状分布,一般分布于坡耕地上,这类地区降水丰沛、坡度较大、人为干扰较多,水土流失极易发生。

北方土壤侵蚀热点集中于黄土高原、山东丘陵区、东北大兴安岭和西北的昆仑山与塔里木盆地的交界处。黄土高原由于其特殊的地貌一直受到很多学者的关注[42],其不稳定的土壤特性和稀少的植被覆盖一直是其成为水土流失热点的原因。大兴安岭和山东丘陵的海拔都较周围高,土壤侵蚀的敏感性也较强。而中国西北部的昆仑山与塔里木盆地之间的区域,其坡度较大,有较多的冰雪消融和降雨带来的径流,植被以草地为主且较为稀疏,城市分布较多人为影响较大,其水蚀也比较剧烈。

2.3 中国土壤侵蚀的时间变化趋势

将年均土壤侵蚀量按照时间序列进行一元线性回归分析后,得出2000年到2015年中国土壤侵蚀强度的时空变化格局(图4),其中有28.72%的区域基本稳定；37.22%的区域土壤侵蚀有增加趋势,34.04%的区域为有减小的趋势。根据1.3的方法将中国年均土壤侵蚀演变显著程度分为5类(图4)并进行统计,其中土壤侵蚀有极显著减小趋势的区域面积为3417 km2,占总有效统计面积的0.74%；有显著减小趋势的区域面积为132593 km2,占2.88%；有极显著增加趋势的区域面积为20823 km2,占0.45%；有显著增加趋势的区域面积为82783 km2,占1.80%。土壤侵蚀显著增加的区域分布在中国长江中下游平原、江南丘陵、云贵高原、四川盆地、东北平原和黄土高原的西北部；显著改善的区域有黄土高原南部、秦岭地区和东南沿海地区。

图4 2000—2015年中国土壤侵蚀的时空变化格局Fig.4 Spatio-temporal variation pattern of soil erosion in China, 2000—2015

由于年际土壤侵蚀的变化体现在C因子和R因子中,可在通过对C因子和R因子进行一元线性回归分析后得出变化斜率(图4),分析出有明显土壤侵蚀强度变化区域的变化因素。通过图4,可以看出中国C因子有明显增大的区域有东北平原、新疆的准噶尔盆地和吐鲁番盆地、南方丘陵区,以及黄土高原东部。分析其原因,东北和西北地区,主要是自然植被的退化和农业逐渐发达。随着南方经济的发展,我国农业逐渐向新疆和东北地区转移[43],天然植被被耕地取代,耕地面积逐渐扩大,粮食作物种类与耕种面积变化,植被因子也随之增大,土壤侵蚀也逐渐增强。而部分地区,如黄土高原南部、秦岭一带,其土壤侵蚀显著改善,原因是陕西省、山西省、河南省实施了退耕还林和修建梯田等水土保持措施,其植被覆盖增加、C因子减小,土壤侵蚀的情况有所好转；并且部分区域伴随着R因子的减少,土壤侵蚀力也有所下降。南方东南沿海区散状分布的C因子的改变,推测与部分农田的荒废或退耕还林的措施有关。而R因子明显增大的区域则在中国长江中下游平原和两广丘陵与云贵高原之间,这些地区年降水量在逐渐增大,这也是其水蚀强度增大的主要原因。

3 讨论

本研究计算的2000年、2005年、2010年、2015年中国的平均土壤侵蚀率分别为38.63、37.35、49.03、47.84 t km-2a-1,其中有2.25%的土地土壤侵蚀情况明显恶化。整体上中国的平均土壤侵蚀量有所增加,但全国土壤侵蚀总面积在减少,其原因主要为C因子和R因子的变动,其表现是部分地区降水量的增减、植被覆盖的变化和耕地粮食作物种类和播种面积的变化。中国社会在2000年到2015年处于快速发展的时期,人口增多、城市化加速,一方面对建筑用地的需求增多,另一方面对粮食的需求也在不断增多,面临着粮食安全的问题[44]。为了达到供求平衡,我国的粮食生产中心逐渐北移[43],农田的开垦和土壤肥力的消耗使得东北和西北地区的土壤侵蚀强度增大。为了改善中国的水土流失,我国实施了一系列修复措施,有34.04%的土地的土壤侵蚀强度有所减小。比如我国水土流失最严重的黄土高原区,在1999年启动退耕还林还草工程后[45],土壤侵蚀现象有了一定的改善。

大尺度土壤侵蚀的定量计算与分析在实现空间准确性和真实性上有一定的困难。虽然RUSLE模型设计的初衷是用于坡面土壤侵蚀的计算,但由于其简洁性和参数的可获得性在大尺度土壤侵蚀估算中有了较为广泛的应用。然而大尺度的计算会提高数据粗糙度,对算法进行合理的修正可以弥补数据不足带来的误差。本研究所利用的C因子、R因子和P因子的算法,将不同土地覆盖类型和气候带的参数进行区分,提高了计算结果的空间异质性与区域准确度。对于R因子,本研究所用的Naipal法对不同气候区进行了分异,并针对不同的气候区利用了不同的算法,同时也结合了海拔和降雨集中性的指标,使得计算结果在中国这个大尺度研究区中更能体现出异质性。目前在国内广泛应用的方法是章文波利用日降水量计算R因子的方法[46]。将易玲[47]利用章文波法计算的中国年平均R因子,与本研究应用的Naipal法的计算结果相比,章文波法计算出的中国平均年降雨侵蚀量比本研究的结果偏大。其原因可理解为丘陵地区的R因子在调整后比未调整的R因子有所降低[22]。而目前国内广泛应用的C因子的计算方法一般是两种：一种是利用NDVI确定C因子,而另一种是根据土地覆盖类型数据,对不同的土地覆盖类型赋予不同的C值。而本研究利用Panagos计算欧洲C因子的算法,区分耕地和非耕地并将土地利用类型、NDVI和土地覆盖度数据融合,得出更加精准的中国C因子分布。P因子则根据农田类型和坡度赋值,体现了人类土壤保持措施的作用。高分辨率的地形因子(30m分辨率)使得坡度因子的精确度达到提升。

虽然数据的空间分辨率的提升和算法的改进使得结果有了更高的准确度,但由于空间尺度大,土壤侵蚀具有空间尺度效应[6],大尺度模型估算与小尺度的土壤调查数据会一定程度的偏差。大尺度的模型计算的原数据精度与计算方法的不同也会导致不同研究的结果不同。例如江西省的年土壤侵蚀率,本研究计算的2015年平均土壤侵蚀率为9060 t/km2,而周夏飞等[14]利用USLE模型在江西省尺度下定量评估的结果为6375 t/km2,Teng等[15]利用RUSLE模型在全国尺度下计算结果为 92 t/km2,陈思旭等[12]基于RUSLE模型在南方丘陵的尺度下计算的结果为354 t/km2。不同研究的研究尺度不同、计算方法不同、数据来源不同,结果差异较大,但分布规律与变化规律大致相似。

中国土壤侵蚀研究应进一步优化各个因子的计算方法、提高数据精度,同时修正风力侵蚀和冻融侵蚀区的土壤侵蚀估算方法,提高其他类型土壤侵蚀量估算的空间准确性。这将为土壤侵蚀的大尺度估算提供更准确的地域信息,也可以为中国不同区域制定针对性的水土保持措施提供理论依据。

4 结论

(1)2000年、2005年、2010年、2015年中国的平均土壤侵蚀率为38.63、37.35、49.03、47.84,总体上中国的土壤侵蚀量有所增加,但近期国家对水土保持的重视以及实施了一系列保护措施,土壤侵蚀的强度有所缓解。

(2)土壤侵蚀量较高的地区集中在中国南方长江中下游平原区、云贵高原、黄土高原区、昆仑山山麓,占统计总面积的9.65%。而土壤侵蚀量明显增大的区域有西北农田区、四川盆地、两广丘陵西部、长江中下游平原和东北平原,面积达10.36×104km2。其原因有自然的降雨增多,也有开垦农田、农作物更替等人为行为干扰。

(3)农业活动增强了西北地区和东北地区的土壤侵蚀；南方城镇化破坏了天然植被,导致土壤侵蚀呈现分散状的增加；但黄土高原南部和东部的水土修复工程的实施使其土壤侵蚀有所缓解。