杨苗苗, 杨勤科, 张科利, 王春梅, 庞国伟, 李玉茹,3
(1.西北大学 城市与环境学院, 西安 710127;2.北京师范大学 地理科学学部, 北京 100875; 3.自然资源部 第一地理信息制图院, 西安 710054)
土壤可蚀性因子(K)是土壤属性对土壤侵蚀影响的度量指标[1],该指标是采用USLE、RUSLE和CSLE等模型来进行区域土壤侵蚀调查与制图所必需的重要数据[1-2]。有研究表明,土壤可蚀性因子除受到土壤颗粒组成等影响外,也受到土壤中砾石(直径>2 mm)含量的影响[3-4]。
砾石是由于成土过程不充分,尚未发育为成熟的土壤而形成的[5]。砾石也可从其他地方搬运而来[6],如地表的滑坡过程、耕作使土壤中的砾石向上移动、水选择性冲刷/剥蚀土壤细颗粒等[7]都有可能产生砾石。研究表明,砾石影响地表径流过程和侵蚀产沙过程[8]。对于砾石含量与环境因子关系方面,Poesen等[9]的研究表明,集约化耕作丘陵区砾石覆盖模式与总地形曲率呈正相关关系,凸坡砾石覆盖率大于凹坡;Simanton[10]和朱元骏[11]等发现,山坡坡度影响砾石覆盖层的特征,砾石覆盖率的差异与山坡坡度密切相关;Nyssen[12]和李燕[13]等解释了耕作动态筛分过程导致了砾石在土壤表面的分布,砾石通过耕作分选被迅速带到土壤表面;Li等[14]的研究表明,砾石覆盖率倾向于随着坡度的增加而增大,南坡和西坡的总覆盖率往往略高于北坡和东坡;王小燕[15]和王慧芳[16]等认为,碎石含量随着土层深度的增加而增加;Meng等[17]研究了不同森林石质土壤岩石碎块特征,发现在0—10 cm土层中,混交林中的砾石含量明显少于纯林;Liu等[18]认为,岩性、构造运动和水文过程的差异可能会使土壤表面产生不同的地貌特征,并在土壤中留下不同的砾石覆盖。尽管前人对于砾石含量与环境要素关系中涉及了土壤可蚀性因子砾石效应(以下简称砾石效应)的影响因素分析,但全球尺度下砾石效应受到哪些因素影响,依然有待认识。
杨苗苗等[4]的研究表明,全球范围内,剖面(0—15 cm)砾石的存在会增加土壤可蚀性,而表面覆盖砾石会减少土壤可蚀性,两者综合影响使土壤侵蚀降低。本研究在对砾石效应(SEK)研究基础上,综合考虑SEK与高程、坡度、坡向、土地利用类型、植被覆盖度、多年平均温度和多年平均降雨量等多个环境要素的关系,分析影响SEK的主控因素,旨在为全面认识砾石效应(SEK)以及应用K因子和土壤侵蚀调查提供参考。
本研究在全球范围内展开,所用基本数据为作者前期研究积累的数据[4],包括砾石覆盖(土壤表面)影响(St),剖面(表层,0—15cm)砾石影响(Kcf)、砾石覆盖和剖面砾石综合影响(Kf-cs);同时有通过ISRIC(https:∥files.isric.org/soilgrids/data/aggregated/),NASA(https:∥earthdata.nasa.gov/esds/competitive-programs/measures/nasadem),GEE(https:∥developers.google.cn/earth-engine/datasets),IPCC(https:∥www.ipcc.ch/)网站下载相关的数据,包括砾石覆盖率数据、剖面砾石含量数据、高程数据、土地利用类型数据、多年平均植被覆盖度NDVI数据、全球气候类型分区图、全球多年平均温度和多年平均年降雨量数据,以及由高程派生的坡度数据、坡向数据、地形起伏度数据等(表1)。
表1 基础数据的来源及其用途
本研究根据土壤发生学、土壤侵蚀和土壤利用的原理,以前期研究形成的认识和计算得到的3个表示砾石效应的参数SEK(砾石覆盖影响St、剖面砾石影响Kcf、砾石覆盖和剖面砾石综合影响Kf-cs)为基础,通过相关分析认识属性相关特征,通过断面分析认识其空间格局及其关系,通过主控因子分析形成对砾石效应(SEK)的影响因素的系统认识。
采用Poesen等[20]的算法计算砾石覆盖下的可蚀性衰减系数St。
St=1-e-0.04(Rc-10)
(1)
式中:St为土壤可蚀性衰减系数;Rc为砾石覆盖率(%)。St越大说明砾石覆盖对土壤侵蚀的削弱作用越明显。
未考虑砾石的土壤可蚀性因子[Kf,(t·hm2·h)/(hm2·MJ·mm)]和有剖面砾石影响的土壤可蚀性因子[Kc,(t·hm2·h)/(hm2·MJ·mm)]的计算公式如下:
(2)
(3)
M=(Msilt+Mvfs)×(100-Mc)
(4)
式中:Mc为黏粒含量(%);Msilt为粉砂粒含量(%);Mvfs为极细砂粒含量(%);OM为有机质含量(%);s为土壤结构等级;Kpf为不考虑剖面砾石影响的土壤渗透性分量;Kpc为考虑剖面砾石影响的土壤渗透性分量。
通过St与Kc相乘可计算考虑剖面砾石和砾石覆盖影响的土壤可蚀性因子Kcs[(t·hm2·h)/(hm2·MJ·mm)]
Kcs=(1-St)×Kc
(5)
将有剖面砾石影响的土壤可蚀性因子(Kc)与未考虑砾石的土壤可蚀性因子(Kf)求差值,得到剖面砾石影响系数[Kcf,(t·hm2·h)/(hm2·MJ·mm)]:
Kcf=Kc-Kf
(6)
Kcf越大,说明剖面砾石对土壤侵蚀的增强作用越明显。
将未考虑砾石的土壤可蚀性因子(Kf)与考虑剖面砾石和砾石覆盖影响的土壤可蚀性因子(Kcs)求差值,得到剖面砾石和表面砾石覆盖综合影响系数[Kf-cs,(t·hm2·h)/(hm2·MJ·mm)]:
Kf-cs=Kf-Kcs
(7)
Kf-cs越大,说明剖面砾石和砾石覆盖对土壤侵蚀的影响作用越明显。
(1) 相关分析。在全球范围内,经纬度每间隔2.5°生成一个采样点,在其上采集SEK的参数(St,Kcf,Kf-cs)和环境要素[高程、坡度、植被覆盖度(NDVI)、多年平均降雨量和多年平均温度]数据,分析SEK因子与环境要素的关系,并进行显著性检验。全球范围内共采集样点数据2 642个,样点均匀分布,具有一定代表性。
(2) 格局分析。借助生态学研究中梯度分析的方法[21],通过布设断面分析不同高程和坡度梯度方向上砾石效应(SEK)的变化。同时在全球范围内,统计不同坡向、不同土地利用类型和不同气候类型区的砾石效应均值,分析其值高低的原因。
(3) 主控因子分析。随机森林(random forest)模型是由Breiman[22]在2001年提出来的一种基于分类树的机器学习算法,该模型是利用bootsrap重抽样方法从原始样本中抽取多个样本,对每个bootsrap样本进行决策树建模,然后组合多棵决策树的预测,通过投票得出最终预测结果。本研究在Python语言平台上进行随机森林回归。
对于随机森林回归的精度检验,本研究采用的评价指标包括决定系数(R2)和平均绝对误差(MAE)。
(8)
(9)
2.1.1 砾石效应与高程关系 砾石含量与砾石效应计算结果如图1所示,砾石效应与高程的关系见图2,图2和表2表明,随着高程增加,St,Kcf和Kf-cs均呈上升趋势,且高程与St,Kcf和Kf-cs之间均存在极显著的相关性(p<0.01)(表2),其间的相关系数分别为0.843,0.252,0.807,其中Kcf与高程的相关性较低,St和Kf-cs与高程的相关性较高。由图1可见,影响较大的范围主要分布在青藏高原、喜马拉雅山脉和科迪拉山系等山脉地区。
根据不同的高程梯度,分别在中国、南美洲、蒙古俄罗斯和美国布设高程断面,如图3A所示(白色断面线),归一化后得到不同断面的高程、St和Kf-cs的剖面图(图4)。由图3可见,St和Kf-cs的剖面图与高程剖面图波动趋势大致相同。
图1 砾石效应的参数
图2 砾石效应与高程关系
表2 环境因子与砾石效应相关性结果
2.1.2 砾石效应与坡度关系 图5显示,St,Kf-cs与坡度之间存在显著的相关性(p<0.01)(表2),相关系数分别为0.715,0.662,而Kcf与坡度之间的相关性较弱,相关系数仅为0.09。随着坡度的增加,St,Kcf,Kf-cs均增加。从全球看(图1),砾石覆盖和砾石综合影响比较明显的地方多位于荒漠和山地高原边缘,主要有俄罗斯远东地区中西伯利亚高原、蒙古高原、中国的塔克拉玛干沙漠、青藏高原、喜马拉雅山脉和美洲的科迪拉山系、北美沙漠、阿塔卡马沙漠等。
根据不同的坡度梯度,分别在中国、西欧、北美以及南美部分区域布设坡度断面,见图3B(白色断面线),归一化后得到不同断面的坡度、St和Kf-cs的剖面图(图6)。由图6可知,St和Kf-cs剖面图波动规律与高程剖面图的相似度较高。
2.1.3 砾石效应与坡向关系 在北回归线以北的地区,南坡为阳坡;在南回归线以南的地区,北坡为阳坡;在南北回归线之间,阳坡需要根据太阳直射点的纬度位置决定。分析砾石效应与坡向的关系时,统计北回归线以北和南回归线以南的地区,可以发现,在北半球,南坡和西坡的砾石效应略大于北坡和东坡,即阳坡的砾石效应大于阴坡,南半球正好相反,与Li等[14]的研究结果一致。
图3 断面布设
图4 高程断面
图5 砾石效应与坡度关系
2.2.1 砾石效应与土地利用类型的关系 从图7可以看出,分析全球范围内,草地、林地、耕地、荒漠和裸地等覆盖类型的St,Kcf和Kf-cs的均值,可以发现,耕地和林地的砾石效应较小,草地、荒漠和裸地的砾石效应较明显。统计不同土地利用类型的坡度和地形起伏度均值(表4),发现耕地通常在比较平的地形部位。
图6 坡度断面
表3 不同坡向砾石效应统计
图7 砾石效应与土地利用统计
表4 不同土地利用类型的坡度和地形起伏度均值
2.2.2 砾石效应与植被覆盖的关系 对比砾石效应与植被覆盖度NDVI的关系,St,Kcf,Kf-cs与NDVI之间均呈现显著性(表2),相关系数分别为0.582,0.049,0.164。图8显示,随着NDVI的增大,St,Kf-cs总体上呈下降趋势,其中NDVI对Kcf,Kf-cs的影响较弱。
2.3.1 砾石效应与多年平均降雨量的关系 对比砾石效应与多年平均降雨量的关系,St,Kcf和Kf-cs与多年平均降雨量之间均呈现显著性(表2),相关系数分别为0.546,0.568,0.523。如图9所示,随着多年平均降雨量的增加,St和Kf-cs呈先增后减趋势,而Kcf在雨量<3 000mm时呈下降趋势,然后略有上升。
2.3.2 砾石效应与多年平均温度的关系 对比砾石效应与多年平均温度的关系,St,Kcf和Kf-cs与温度之间均呈现显著性(表2),相关系数分别为0.597,0.468,0.477,说明砾石效应与多年平均温度的相关性较高。图10显示,St和Kf-cs与多年平均温度呈负相关,Kcf与多年平均温度呈正相关。
图8 砾石效应与植被覆盖度关系
图9 砾石效应响与多年平均年降雨量关系
图10 砾石效应与多年平均温度关系
2.3.3 典型气候区的砾石效应分析 根据Peel等[19]得到的全球气候类型分区图的分类标准,将全球分为热带赤道地区、较温暖湿润地区、北方寒冷地区、干旱地区和极地地区(图11),并分别统计不同气候类型区的砾石效应(表5)。可以发现,对于砾石覆盖影响St和砾石综合影响Kf-cs而言,热带赤道地区影响最小,干旱地区影响最大,即温度高降雨量大的区域砾石效应小,温度低降雨量小的区域砾石效应较大。
但较温暖湿润地区比北方寒冷地区砾石覆盖影响和砾石综合影响比北方寒冷地区大,原因是,在北方寒冷地区,主要分布有西西伯利亚平原、东欧平原和北美中部大平原等,平均海拔为548.77m,平均坡度为17.25°;而在较温暖湿润地区,主要分布有巴拉纳高原、阿巴拉契亚山脉、加丹加高原、隆达高原、伊朗高原等,平均海拔为735.48m,平均坡度为24.03°。在这两个地区St和Kf-cs主要受地形影响,受温度和降雨影响较小,所以北方寒冷地区的砾石效应要小一些。对于剖面砾石影响Kcf而言,影响从大到小的地区依次是干旱地区、热带赤道地区、温暖湿润地区和北方寒冷地区,即温度越高剖面砾石效应越大,与前文分析规律一致。
将全部影响因子作为自变量,对SEK(St,Kcf和Kf-cs)进行随机森林回归,并分别计算经过随机森林回归方法得到的SEK预测值和SEK真实值的R2和MAE(图12)。由图12可知,St,Kcf和Kf-cs的R2分别为0.593,0.267,0.718,说明7个环境因子(高程、坡度、坡向、NDVI、土地利用类型、多年平均降雨和多年平均温度)对Kf-cs的解释度最高,其次是St,对Kcf的解释度较弱;St,Kcf和Kf-cs的MAE分别为0.096,0.000 7,0.003,与St,Kcf和Kf-cs的真实值取值范围相比都较小,说明随机森林对St,Kcf和Kf-cs的预测值偏离误差都较小。
图11 全球气候类型分区
由图13可知,高程和坡度对St的影响较大,其次是多年平均温度和植被覆盖度,坡向、多年平均降雨和土地利用影响较小;多年平均温度和多年平均降雨对Kcf的影响较大,其余各因子影响都较为微弱;高程和坡度对Kf-cs的影响较大,其次是多年平均降雨和坡向,土地利用类型、多年平均温度和植被覆盖度的影响较小;说明高程和坡度是St和Kf-cs的主控因子,多年平均降雨和多年平均温度是Kcf的主控因子。原因是St和Kf-cs主要是砾石覆盖的影响作用,在土壤表面,更容易受到地形等外界因素的影响,而Kcf是土壤剖面中砾石的影响,较为稳定,受地形等外部因素的影响仅为间接作用。
表5 不同气候区砾石效应统计
图12 随机森林回归砾石效应验证散点图
图13 影响因子重要性排序
土壤由基岩风化过程、生物化学过程、地表物质运动过程、人为作用过程等作用下形成。土壤形成比较充分则土壤比较黏重,否则粗颗粒较多。砾石正是由于成土过程不充分,尚未发育为成熟的土壤而形成的[5]。另一方面,后期的土壤侵蚀也有可能产生砾石,如荒漠是过于干旱、生物活动微弱、风大所致;高山(冰雪带以下)是温度过低、或水蚀过强所致;砾石也可从其他地方搬运而来[6],如地表的滑坡过程、耕作使土壤中的砾石向上移动、水选择性冲刷/剥蚀土壤细颗粒等[10]都有可能产生砾石。有研究表明,砾石主要分布在坡度大(图3)[11]、高程高(图3)[5]、阳坡[14]、凸坡[9]、较为寒冷干旱的山地和荒漠(图11)[12]等地区,本研究的结论虽与一些认识一致,但是本研究量化分析更能说明砾石效应(SEK)的影响因素。
由于缺少全球范围内的砾石覆盖率的实测数据,而且目前最高分辨率的卫星影像也无法提取这些信息,作为一种初步的尝试,本文用Hengl等[23]通过机器学习得到的Soilgrid土壤表面砾石含量数据代替砾石覆盖率。一些学者尝试通过构建模型来预测砾石覆盖率,如Simanton等[10]的研究中提到,知道土壤剖面50mm内的砾石含量,使用基于坡度的方程可以更好地估计地表岩石碎块覆盖率,计算公式如下
RFC=68.62-80.62/RFPSL0.5
(10)
式中:RFP为土壤剖面砾石含量(%);SL为坡度(%);RFC为土壤表面砾石覆盖率(%)。
基于公式(10)计算结果(图14),基本可以表达丘陵山区砾石覆盖的分布大势,与Soilgrid数据(图15)整体趋势相似,相关性达到0.518。但是均值差异较大(分别为13.40%,9.19%),在较平坦地区也有所不同(如澳大利亚、非洲、中亚等地区)。存在差异的原因应该是上述公式只考虑坡度,没有海拔、土地覆盖类型、气候要素等参数,荒漠地区的估计偏低。然而更全面、精确地分析有待于通过野外抽样调查、无人机影像判读等方法,改善砾石覆盖率数据的精度。
注:图14是基于公式(10)的计算结果。
3.2.1 地形对SEK影响的解释
(1) 高程对SEK影响。高程可以被认为是降水和温度的综合影响的一个综合表现[5]。高程最高的地方在青藏高原、喜马拉雅山脉和科迪拉山系等山脉地区;这种高、旱、寒环境下,以物理风化为主,成土过程弱,因此土壤表面砾石含量高,砾石效应也较大。而且不同高程的水热条件对土壤理化性质的也有间接影响,高程越高,温度较低,水分较少,砾石效应较大。
图15 SoilGrids数据(均值9.19%)
(2) 坡度对SEK影响。坡度的差异引起的砾石在坡面上的分选[15],从而影响不同坡面的侵蚀强度[11]。坡度越大,侵蚀发生越剧烈,砾石在该处发生侵蚀性富集,覆盖度越大,砾石覆盖和砾石综合影响也就越大。随着坡度的增加,砾石覆盖率的增长速度逐渐变慢。山坡坡度下砾石覆盖率的普遍增加是由于植被稀疏的山坡上发生的水蚀和沉积过程造成的。
(3) 坡向对SEK影响。坡向影响给定景观的太阳能吸收率。在北半球,朝南的斜坡比相对应朝北的斜坡更垂直于太阳光线,通常要更加温暖,从而一般湿度较低。因此,朝南的斜坡风化程度不深,植被覆盖率较低,雨水和径流侵蚀率较高。朝南的山坡段经历了强烈的层间和细沟侵蚀,因此具有较高的砾石覆盖率,砾石覆盖影响也更大,南半球则正好相反。
3.2.2 土地利用与覆盖对SEK影响的解释 耕地通常在比较平的地形部位,耕地的坡度和地形起伏度较小,草地、荒漠和裸地地区通常较干旱,生物活动微弱,风大导致有较多的砾石覆盖,而高山(冰雪带以下)地区温度过低,或由于水蚀过强导致砾石含量较多,从而有较大的砾石效应。
对于植被覆盖NDVI而言,随着NDVI的增大,St和Kf-cs整体上呈下降趋势,而NDVI对剖面砾石影响Kcf的影响较弱。原因是植被覆盖度越大的区域,土壤条件更好,砾石含量较低,所以砾石效应也会较小。
3.2.3 水热条件对SEK影响的解释 由于降雨的影响,砾石受到侵蚀分选作用的影响,会缺乏地表富集,而且降雨量较大的区域主要分布在亚洲东南部、非洲中部、北美洲东南部和南美洲北部地区,这些地区地势较低,多植被覆盖,土壤类型主要为冲积土、黑土、黑钙土等,这类土壤土层较深,营养较为丰富,所以砾石含量相对较小,影响也较小。
气温通过影响土壤中的微生物群落及总量使其土壤结构发生变化,对土壤理化性质有较大影响,最终影响土壤可蚀性[24],气温越高的地区,水热条件好,地势平缓,土壤条件更好,土壤表面砾石覆盖较少,砾石覆盖的影响会降低。而当温度较高时,土壤中生物和化学反应速率会增快,土壤剖面中的砾石效应会有一定的增强。
3.2.4 影响SEK的其他要素 已有研究表明,耕作对砾石含量的分布也有一定的影响[25],进而影响土壤侵蚀。这是由于耕作一方面对不同粒径的颗粒具有筛选作用,导致大的颗粒向上运动以及小的颗粒向下运动;另一方面耕作侵蚀使得表层土壤颗粒被径流带走。这两方面原因增加土壤表面砾石覆盖率,客观上对土壤也起到一定的保护作用。但由于全球范围内,耕作的影响无法定量估算,本文暂未考虑。
也有研究表明,砾石分布还与林种[17]、母质的物理和化学性质[15]、动植物扰动、碎屑扰乱(树木连根拔起)、冷冻扰动(冻融)、变硬扰动(黏土的收缩和溶胀))、重力沉降、晶体生长和浪费和地震作用(地震运动)[6]等有关,但在全球范围内,这些数据难以获取和定量估算。今后在小流域范围内研究砾石效应的影响因素,或许可将这些因素也予以考虑,也可利用高精度(如5 m)的实测高程、坡度、砾石覆盖率等数据进行更细致的分析。
砾石效应可在K因子中考虑[3],也可在(R)USLE-C[1]中考虑;如果不予考虑,则砾石含量高的地区(山地、高原和荒漠等)的土壤侵蚀速率会被高估。目前在土壤可蚀性因子计算时,大部分未考虑砾石的影响,导致了评价结果的偏差,如Borrelli等[26]的图上山地被明显高估,Teng等[27]文章中秦岭、川西地区被高估,Teng等[25]青藏高原东部被高估。
以科迪勒拉山系、阿拉伯高原、伊朗高原、阿巴拉契亚山脉、喜马拉雅和喀喇昆仑山区为例,这些地区海拔高、地形起伏度较大,地表以物理风化为主,地表砾石覆盖率很高,土壤侵蚀比较微弱。但是Borrelli等[26]的土壤侵蚀评价结果中,这些地区土壤侵蚀速率明显高于相邻地区。我们猜想可能性的原因是没有充分考虑砾石效应。所以,应该对重点地区(海拔和坡度大且土壤水热条件不够好的山地、降水量较小的荒漠地区),考虑砾石的影响,并对土壤侵蚀评价结果的影响做出系统分析。
生产实践和有关研究中,有认识到砾石的覆盖具有抵抗土壤侵蚀、保护土壤的作用,如我国西北的砂田是一种在土壤表面覆盖一层砾石的农田,主要集中分布在甘肃省中部地区,以及青海、新疆和宁夏的部分地区[28]。在世界上其他降水稀少的地方也有砂田分布,如美国南部干旱地区的德克萨斯州等地[29],瑞士的沙莫松以及西班牙的加那利群岛等地[30]。砂田能有效地减少土壤蒸发和径流,提高土壤水分入渗效率和土壤温度,阻止水土流失和次生盐渍化,增加砂田作物产量、水分生产效率和土壤微生物数量[28]。Nyssen等[12]的研究指出,农民通常不愿意将较小的砾石从他们的土地上拿走,因为他们认为它们对土壤水分保持和表土保护免受侵蚀有积极影响。在RUSLE手册[1]中也有提到,土壤中的砾石会影响径流,进而影响细沟间的侵蚀,使计算得到的土壤可蚀性因子值不够准确。但这些认识尚没有在区域土壤侵蚀调查与评价制图中考虑。
(1) 土壤可蚀性因子(K)的估算与砾石含量有关,砾石效应(SEK)受到高程、坡度、坡向、土地利用类型、植被覆盖、多年平均降雨量和多年平均温度等环境要素的影响;其中,高程和坡度与SEK(St,Kcf和Kf-cs)呈正相关;植被覆盖度、多年平均降雨量和多年平均温度与St和Kf-cs呈负相关;多年平均降雨与Kcf呈负相关;多年平均温度与Kcf呈正相关;阳坡的砾石效应大于阴坡;耕地和林地的砾石效应较小,草地、荒漠和裸地的砾石效应较大。
(2) 砾石覆盖影响St和砾石综合影响Kf-cs的主控因子为高程和坡度,原因是St和Kf-cs主要是砾石覆盖的影响作用,在土壤表面,更容易受到地形等外界因素的影响。Kcf的主控因子为多年平均降雨和多年平均温度,原因是Kcf为土壤剖面中砾石的影响,较为稳定,受地形等外部因素的影响仅为间接作用。
(3) 在大区域乃至全球性土壤侵蚀评价中,对于砾石含量较大、海拔较高、坡度较大的区域(如科迪勒拉山系、阿拉伯高原、伊朗高原、阿巴拉契亚山脉、青藏高原和天山等地区),均应给予特别关注,如果不加考虑,则有可能使土壤侵蚀评价结果出现偏差。
致谢:澳大利亚联邦科工组织水土资源研究所David Jupp先生协助订正英文摘要,特此致谢。