黄小芳, 林丽蓉,2,3†, 高 超, 丁树文,2,3, 朱新远, 郭靖东
(1.华中农业大学资源与环境学院, 430070, 武汉; 2.华中农业大学农业农村部长江中下游耕地保育重点实验室, 430070, 武汉;3.国家环境保护土壤健康诊断与绿色修复重点实验室, 430070, 武汉; 4.湖北省水土保持监测中心, 430071, 武汉)
我国于20世纪50年代开始土壤侵蚀模型的定量研究,在80年代引入美国通用土壤流失方程(USLE)后,许多学者开始以USLE模型为蓝本,参考径流小区的观测数据,根据区域实际情况对USLE模型进行修正,并建立具有区域特征的土壤侵蚀预报模型[1-8]。吴发启等[5]对黄土高原地区,杨子生[7]对滇东北山区分别进行土壤侵蚀预报模型探索。其中,最典型且应用范围最广的为刘宝元提出的中国土壤流失方程(Chinese soil loss equation,CSLE)[9]。CSLE模型是基于黄土丘陵沟壑区的安塞、离石、延安、子洲等地的径流小区监测数据,并充分考虑我国水土保持措施和当地水土流失实际情况而建立起来的。该模型是SL190—2007《土壤侵蚀分类分级标准》中确定水力侵蚀模数的重要方法之一,曾被用于第一次全国水利普查水土保持情况普查和第四次全国土壤侵蚀普查[10-12]。
CSLE模型的基本形式为:
A=RKLSBET。
(1)
式中:A为年均土壤流失量,t/(hm2·a);R为降雨侵蚀力因子,MJ·mm/(hm2·h·a);K为土壤可蚀性因子,t·hm2·h/( hm2·MJ·mm);L、S分别为坡度、坡长因子或统称为地形因子,量纲为1;B为植被覆盖与生物措施因子,量纲为1;E为工程措施因子,量纲为1;T为耕作措施因子,量纲为1。
CSLE模型将USLE模型中的作物管理覆盖(C)与水土保持措施(P)2大因子扩展为植被覆盖与生物措施(B)、工程措施(E)和耕作措施(T)3大水土保持措施因子,以反映我国的农业生产和水土保持状况[13];此外,CSLE模型中对陡坡的地形因子进行改进,使之更加符合我国地形地貌实际情况。
在CSLE模型中,T因子反映水土保持耕作措施对土壤侵蚀的作用[14]。T因子取值范围为0~1,T因子值越小表明水土保持耕作措施效果越好[15]。目前,国内T因子的研究主要集中在应用CSLE模型进行土壤侵蚀强度评价和消长分析,T因子的时空异质性以及T因子赋值系统等方面[16-19]。目前,T因子尚没有统一的获取方法,水利普查结果查询、区域水土流失动态监测查询、土地利用类型划分公式、径流小区公式、坡度分级及水土保持耕作措施赋值这几种方法均有应用。而T因子获取方法的不统一,在一定程度上影响区域之间的CSLE模型计算结果对比。此外,部分学者在应用CSLE模型时没有充分考虑研究区的气候、土壤、地形、耕作方式等区域差异性,而直接采用已有文献其他地区的研究成果进行T因子赋值[20-21]。研究表明,同一耕作措施水土保持效益因不同地区气候、地形和土壤类型有所不同,T因子值差异较大[15,22]。这无疑会增加CSLE模型的计算误差,无法客观反映研究区的实际土壤侵蚀状况,也不利于区域之间的土壤侵蚀量对比。
笔者基于已有相关研究成果,总结归纳水土保持耕作措施的定义、分类及水土保持耕作措施因子(T)的获取方法,并比较各方法的优缺点及适用条件,以期为今后开展区域土壤侵蚀评价和比较不同水土保持耕作措施的水土保持效益提供参考依据。
水土保持耕作措施,是指通过改变微地形、增加地表覆盖或改善土壤物理性质等,提高土壤抗蚀性,防治土壤侵蚀的耕作栽培措施[23]。水土保持耕作措施分类是保证T因子正确获取的重要基础。USLE模型将水土保持耕作措施的影响归在覆盖与管理因子C和水土保持措施因子P中,而国内学者在水土保持措施一级分类中将耕作措施单独划分出来,以充分体现我国农地土壤侵蚀和水土保持耕作的特征[24]。GB/T 16453.1~6—2008《水土保持综合治理技术规范》将耕作措施类型分为等高耕作、等高沟垅种植、垄作区田、掏钵(穴状)种植、草田轮作、间作与套种、休闲地绿肥、留茬少耕、免耕和轮作等14类[25]。水土保持耕作措施的准确分类,有利于我国T因子值确定的规范化,提高不同区域T因子值的可比性。
T因子值主要根据不同坡度条件下耕作措施减少水土流失的作用来确定的;但不同学者关于T因子值的计算存在不同的理解。有学者认为T因子值等于采取某种耕作措施下的土壤流失量与同等条件下无耕作措施(裸地小区)的土壤流失量之比[26]。也有研究将采取某种耕作措施的土壤流失量与同等条件下传统耕作(一般指顺坡平作或垄作)措施的土壤流失量之比作为T因子值[15,27]。由此可见,不同学者对T因子计算过程中对照小区还存在不同的理解,而对照小区的不同选择(或不统一)势必影响T因子值的最终结果。
不同水土保持耕作措施,尤其是不同农作物的轮作,对土壤流失的影响有明显差异。T因子既可以反映水土保持耕作措施影响,也可以反映轮作措施形成的作物覆盖影响。前者是通过耕作以改变微地形和增加地表覆盖、土壤入渗等,进而提高土壤抗蚀性能,达到保水保土和防治土壤侵蚀的方式[23]。后者(轮作)是指同一田块有顺序地在季节间或年份间轮换种植不同作物或复种组合的种植方式,反映作物生长过程中覆盖的影响,作物种类不同、种植结构的时间和空间差异都会对土壤侵蚀有明显影响[28]。当土地利用类型是耕地时,应考虑轮作的影响。根据我国主要农作物(稻谷、小麦和玉米等10类)生长期的降雨侵蚀力季节变化进行修正可以得到T因子值[15]。基于上述T因子含义,笔者依据相关研究成果,总结以下几种获取T因子值的方法。
2010—2012年,我国开展第一次全国水利普查水土保持情况普查,基于全国各地径流小区的观测数据对T因子值进行修正,并发布了区划县(市)的T因子参考值[29]。水利普查结果查询法,主要是通过查询研究区所属的轮作制度分区类型,进而在水利普查结果中查询对应的T因子值,对研究区进行T因子赋值。齐斐等[30]根据水利普查结果提供的措施因子参考值,研究沂蒙山区不同抽样密度和推算方法对土壤侵蚀状况及其精度的影响。李子轩等[31]基于CSLE模型,通过水利普查结果查询法获取T因子值用以对比全域覆盖计算与抽样单元推算法计算县域土壤侵蚀模数的差异。Duan等[32]利用查询法获得不同作物轮作系统的作物轮作图和T因子值,并使用不同作物轮作系统的实测数据对T因子值进行校正。该方法虽被广泛应用,但需要结合研究区域内的调查单元和径流小区的数据资料对T因子值进行校正。
2018年8月,水利部水土保持监测中心颁布《区域水土流失动态监测技术规定(试行)》,其中全国轮作区分区参照《中国耕作制度70年》的区划县(市)名录。当土地利用为耕地时,学者可通过技术规定中查询研究区所属的轮作区代码及其对应的T因子值,并进行赋值。如陈锐银等[33]参考区域水土流失动态监测技术规定,查找研究区轮作区名称及代码对T因子进行赋值,从而计算四川省省级水土流失重点防治区土壤侵蚀状况。石秀梅[34]在福建省重点治理区水土流失动态监测研究中,查询区域水土流失动态监测技术规定中研究区域内的T因子值。刘桂成[35]从区域尺度下定量评价江西省土壤侵蚀及其空间分布中,经查询法获取T因子值。区域水土流失动态监测查询法与水利普查结果查询法的T因子赋值原理一样,均根据轮作措施进行赋值。该获取方法相对简单,因此被普遍应用于流域及区域尺度研究中,但是不同研究区的水土流失环境及田间具体管理措施不同,T值的准确性难以保证。
为减小查询法获取T值所产生的区域差异性,有学者根据高分辨率遥感影像的土地利用解译和水利普查结果查询法分别获取研究区的土地利用状况和耕作措施因子,并结合野外单元调查结果,通过加权平均方法得到不同土地利用类型的T因子值,计算公式如下:
(2)
式中:T为某一土地利用类型的耕作措施因子值;Ti为某一耕作措施的T值;βi为某一耕作措施占该土地利用类型的面积比例,%;m为耕作措施的类别数。
杨韶洋[36]根据水土保持情况普查培训教材提供的措施因子值表,结合野外调查单元调查结果和计算式(2),确定沂蒙山重点治理区不同土地利用的T因子值。陈美淇等[37]通过水蚀野外调查结合土地利用调查,按照比例平均分配计算T因子值。此外,Zhang等[38]参考水利普查结果的T因子值,通过加权平均法对山东省沂蒙山区沂水县T因子值进行计算。土地利用类型划分公式法充分考虑各水土保持措施面积比例的影响,可应用于小流域、流域、区域几种较大尺度,但这种方法确定的T因子值易受遥感影像中土地利用/覆盖类型解译精度的限制。
基于水土保持耕作措施因子值的定义,利用径流小区观测资料,可获得不同地区相关的T因子值。但由于径流小区坡度不一致,计算各种措施T因子值时,需要将研究区监测点所有径流小区的土壤侵蚀量校正为标准小区(20 m×5 m,坡度15°)的土壤侵蚀量[39],校正公式为:
(3)
式中:Ai15为第i小区校正到标准小区的土壤侵蚀量,kg;Ai为第i小区的土壤侵蚀量,kg;S15为15°小区的坡度因子;Si为第i小区的坡度因子。其中,坡度因子S的计算方法是结合McCool等[40]的缓坡条件和刘宝元等[41]的陡坡条件下的坡度因子公式(式4)。
(4)
式中:S为坡度因子,量纲为1;θ为坡度,(°)。
范建荣等[42]在东北地区利用径流小区公式法对大豆的顺坡种植、横坡种植和地梗植物带的T因子值进行计算和比较。郭继成等[43]利用贵州龙里、遵义、毕节等地的野外径流小区资料,将当地主要农作物(马铃薯和玉米)的生长阶段分为苗期、发育期、成熟期、残茬期4个阶段,并获得各阶段的土壤流失比及T因子值(0.005 2~0.834 1)。虽然用该方法获取T因子值较为准确,但仅适用于小区域、小范围的研究,且需要对径流小区进行长期连续观测。
当因研究区面积过大,学者难以通过遥感影像和实地调查完整准确获得所有的水土保持耕作措施类型时,可以用坡度分级法获取T因子值,即根据不同坡度条件下等高耕作减少的土壤侵蚀量来确定(没有耕作措施的地方的T因子值为1)[44]。2006年,刘宝元[45]通过不同坡度条件下的等高耕作减少土壤流失量来获取T因子值,对西北黄土高原区进行了土壤侵蚀定量评价。此后,众多学者均参考该方法,分别对陕西[12]、广西岑溪市和大化县[46]、山西[47]、福建宁化县谢家小流域[48]等地区进行区域坡度分级,计算T因子值(表1)。通过坡度分级给T因子赋值,可以不依赖标准径流小区观测,简单快速,但是需准确获取研究区域的坡度信息,并对研究区域坡度进行合理的划分。
表1 各地区坡度分级和T因子赋值Tab.1 Slope classification and T factor assignment in different regions
该方法基于郭乾坤等[15]的研究成果,通过收集多年的文献资料和监测资料,对全国水蚀区等高沟垄种植和垄作区田等耕作措施进行统计,以区域差异不显著分析结果为依据,最终得到适用于全国的主要T因子值。周晓莹等[22]通过该方法获取我国6个水土保持区的6种耕作措施的T因子值。该方法修订和计算了我国主要T因子值(表2),保证了第一次全国水力侵蚀普查的顺利实施,但对其他耕作措施需要进一步观测和研究。
表2 我国水蚀区主要耕作措施的T因子值Tab.2 T factor values of main tillage measures in the water erosion area of China
T因子反映了轮作制度形成的作物覆盖和水土保持耕作措施对土壤侵蚀的影响。本文主要总结了6种T因子值的有关获取方法,不同方法的考虑范围、优缺点和应用范围均有所不同(表3)。水利普查结果查询法,因使用便捷,被广泛应用,但需结合研究区域内调查单元、径流小区数据,对T因子值进行校正;区域水土流失动态监测查询法,简单易获取,但只针对我国主要10种农作物,种植其他作物的T值将难以确定;土地利用类型划分公式法便于计算,且可应用于小流域、流域、区域尺度,但需要综合遥感解译和查询法获取土地利用和各种耕作措施的具体因子值,随着尺度增大需要花费较多的时间、精力、费用等;径流小区公式法是最基础的方法,其用于获取T因子值较为准确,但对照小区和有措施小区需要定期维护管理和长期监测小区泥沙数据。对不同尺度的流域来说,径流、土壤侵蚀和产沙量等都明显不同,比较不同耕作措施对土壤侵蚀的影响很困难;坡度分级法可操作性强,不依赖标准径流小区,但需要对坡度进行准确提取;水土保持耕作措施赋值法操作简单,便于查找,但缺乏其他主要耕作措施的校正,在不同研究区的准确性与适用性有待考察。
表3 T因子估算方法优缺点及适用范围比较Tab.3 Comparison of advantages and disadvantages of T factor estimation method and its application scope
随着国内学者应用CSLE模型开展水土流失研究的不断深入,水土保持措施耕作因子作为CSLE模型较难准确获取的重要参数之一,其相关确定方法及计算应用必须全面梳理。笔者较系统地回顾了前人有关T因子的研究成果,总结6种获取T因子值的不同方法,并分析和比较了这些方法的适用条件及优缺点。国内学者针对T因子开展了大量研究,主要集中在:1)T因子值被应用于CSLE模型中,用以评价和分析研究区土壤侵蚀强度和消长趋势等。2)T因子值的时空异质性研究,不同学者研究表明,同一水土保持耕作措施因不同地区的气候、地形和土壤类型不同,T因子值差异较大[15,22];3)T因子赋值系统完善,我国学者对东北黑土区、北方土石山区、西北黄土高原区和南方红壤丘陵区等水土保持区均有进行研究,积累了大量的经验与数据库,并根据全国轮作制度区划确定了主要作物的轮作因子值。
T因子值估算研究中,还需对以下4方面继续深入研究。1)尽管全国轮作制度区划及轮作措施均按照各地区主要10种农作物进行划分赋值,但由于农民耕作方式的自主性,土地利用方式并不统一,建议各地要结合本地区的轮作制度开展较长期的径流小区不同耕作措施的水土流失动态监测,以获取能表征当地实际水土流失状况的T因子值;2)当研究区为大尺度范围时,需要结合全国轮作制度区划对T因子进行赋值,并综合空间信息技术手段和野外实地调查对研究区域其他耕作措施的T因子进行计算和查缺补漏,从而准确获取水土保持耕作措施效益;3)相关研究表明不少学者在水土流失相关研究中直接参考前人的T因子值进行赋值,而没有对其适用性与准确性进行验证。在今后的研究中建议充分考虑研究区域的特殊性以及所选择方法本身的适用范围,确保获得的T因子值更准确;4)不同学者通过径流小区公式法计算T因子值时,因耕作措施的对照小区选择不统一,导致实际获取T因子值有差异,这方面需要进一步探讨。