砾石含量对土壤可蚀性因子估算的影响*

杨苗苗，杨勤科†，张科利，李玉茹，3，王春梅，庞国伟

（1. 西北大学城市与环境学院，西安 710127；2. 北京师范大学地理科学学部，北京 100875；3. 自然资源部第一地理信息制图院，西安 710054）

土壤可蚀性因子（soil erodibility factor，K）是USLE[1]、RUSLE[2]，以及EPIC[3]、AGNPS[4]等模型的基本参数；可蚀性因子值又是采用USLE/CSLE等模型来进行区域土壤侵蚀调查与制图所必须的重要数据[5-7]。自Wischmeier[1]正式提出土壤可蚀性因子（K）以来，国内外学者对此进行了大量研究，提出了若干种土壤可蚀性因子估算方法。比较成熟的算法有USLE-K算法[1]、RUSLE2-K算法[8]、EPIC-K算法[3]和基于平均几何粒径的算法[9]等四种。这些算法均基于土壤理化性状指标（粒径、有机质、结构和渗透性等）来计算土壤可蚀性K值（式（1）～式（3））[1]。土壤可蚀性因子计算模型中的粒径包括黏粒含量、粉粒含量、砂粒含量，但均不涉及土壤剖面砾石含量和土壤表面砾石覆盖（Rock Fragments，≥2mm颗粒）。

而大量研究表明，土壤剖面中的砾石含量（以下简称剖面砾石）和土壤表面砾石覆盖对土壤水文、土壤侵蚀产沙过程和土地利用等具有重要影响[10-13]。Brakensiek等[14]建立了基于剖面砾石含量推算饱和导水率的方法。Poesen等[15]证明，土壤表层砾石覆盖可以减少土壤可蚀性，并可以用指数衰减函数来定量估算。Zhou和Shao[16-17]认为，当土壤中砾石含量大于40%时，土壤渗透速率下降，土壤饱和导水率增加。Jomaa等[18]室内水槽试验表明，砾石覆盖 可以减缓雨滴打击，并减少地表水流，降低径流输移能力，一定程度上可以减轻土壤侵蚀。Hlaváčiková等[19]发现剖面砾石可以使饱和导水率降低，进而增加地表径流和侵蚀产沙。马晨雷等[20]也印证了砾石覆盖使坡面总侵蚀量呈负指数减少，且砾石覆盖在冲刷流量较小的情况下，对坡面侵蚀量的减少效果更明显。

尽管Römkens等[9]在RUSLE手册中提到了砾石对土壤可蚀性因子的影响，并将其分解为土壤剖面中砾石含量影响和土壤表面砾石覆盖影响两个方面。但是在土壤可蚀性因子计算与制图研究中（如梁音等[21]；王彬[22]），均未考虑剖面砾石和砾石覆盖对土壤可蚀性的影响。Panagos等[23]在欧洲土壤可蚀性因子计算时，对剖面砾石和砾石覆盖的影响进行了分析，并认为砾石校正因子的应用使K-因子平均降低15%。若不考虑剖面砾石的影响，剖面含砾石较多土壤可蚀性因子将被低估；若不考虑砾石覆盖的影响，分布于山地丘陵土壤的可蚀性将被高估。迄今为止，仍然没有全球范围内砾石含量对土壤可蚀性因子计算值的订正方法和结果。

本研究在全球范围内，利用比较容易获取的表层土壤（0～15 cm）质地等级、土壤剖面砾石含量、土壤表面砾石覆盖等分析估算剖面砾石含量和砾石覆盖对土壤可蚀性因子的影响，以期使土壤可蚀性的计算更加完整准确，进而提高区域土壤侵蚀预报的精度，本研究对全球土壤碳收支的模拟研究也具有一定意义。

1 材料与方法

1.1 研究区域及数据

作为全球土壤侵蚀制图研究计划的一部分，本研究在全球范围内展开。考虑到南极大陆和北极的格陵兰地区数据资料的缺失，同时这些地区多被冰雪覆盖，因此本研究不包括南极大陆和格陵兰等地。为了便于全面分析土壤剖面砾石的影响及其空间分布，研究区包括了水蚀地区、风蚀地区和高寒地区。研究区总面积约1.32亿km2，约占地球表面面积（1.48亿km2）的89.2%，与全球土壤属性数据产品（SoilGrids）覆盖的范围基本一致。

本研究所用基本数据为 ISRIC 发布的SoilGrids[24]，包括土壤类型图和土壤性状图（表层0～15 cm和土壤表面）等（表1）。该数据最新版的原始分辨率为7.5弧秒，分析过程中重采样成30弧秒。

1.2 土壤可蚀性因子计算

考虑到USLE模型中的土壤可蚀性因子（K）算法较为全面地考虑了各种土壤理化性状，也在区域土壤侵蚀调查制图中得到了比较广泛的应用（如Lu等[25]；梁音等[21]；Panagos等[23]；Borrelli等[5]），本研究以此为例讨论土壤粗颗粒物质对土壤可蚀性因子的影响。USLE的土壤可蚀性因子计算公式如下：

式中，Mc为黏粒含量，Msilt为粉砂粒含量，Mvfs为极细砂粒含量，OM为有机质含量（%），s为土壤结构等级，Pf为渗透性等级，Kp为土壤渗透指数对可蚀性的影响，不考虑剖面砾石影响的计算结果为Kpf，考虑剖面砾石影响的计算结果为Kpc。

1.3 剖面砾石含量影响的计算

本研究根据土壤孔隙度和饱和导水率关系，及Brakensiek等[14]的方法，通过分析剖面砾石对土壤渗透等级的影响来估算剖面砾石对土壤可蚀性因子的影响。具体计算过程如下：用土壤质地等级和表 2规则[26-28]，推算未考虑剖面砾石条件下的渗透等级指数Pf；用Pf据表2反推饱和导水率Ks；然后推算有剖面砾石（也叫田间土壤）和没有剖面砾石的饱和导水率（分别为Kb和Ks）的比例（式（4）），并根据式（5）反推包含了剖面砾石的土壤饱和导水率Kb（式（5））。

式中，Kb为含有剖面砾石的土壤饱和导水率，Ks为没有剖面砾石的土壤饱和导水率，CRFVOL为土壤剖面（0～15 cm）砾石含量。

根据表2中Pf与最后一列拟合成一个简单模型 （式（6）），再将Kb代入，计算田间土壤（包含剖面砾石）的渗透性等级指数Pc。

表2 基于土壤质地等级/饱和导水率和渗透系数（P）关系表 Table 2 Relationship between classes of soil texture/saturated hydraulic conductivity and permeability coefficient（P）

最后，将Pf和Pc代入式（7）和式（8）计算未考虑剖面砾石的土壤可蚀性因子（Kf）和有剖面砾石影响的土壤可蚀性因子（Kc）。

1.4 砾石覆盖影响的计算

野外调查和观察表明，丘陵山区的土壤也广泛存在砾石覆盖[10，29]，但目前尚无砾石覆盖度的数据产品。本研究将土壤表面砾石含量近似作为砾石覆盖度，采用Poesen的算法（式（9））[15]计算砾石覆 盖下的可蚀性衰减系数St。

式中，St为土壤可蚀性衰减系数，Rc为砾石覆盖率（%）。

1.5 剖面砾石含量和砾石覆盖的总影响计算

将St与式（7）和式（8）计算结果结合，土壤属性数据代入式(1)，利用本项目开发的土壤可蚀性因子计算工具（K_tool），计算得到质地分量kt、有机质分量k0、土壤结构性分量ks，通过不考虑剖面砾石和考虑剖面砾石的土壤渗透性分量（Kpf和Kpc）推算得到不考虑剖面砾石和考虑剖面砾石的土壤可蚀性因子（Kf,Kc）。将砾石覆盖影响因子（St）与Kc相乘，得到考虑剖面砾石和砾石覆盖影响的土壤可蚀性因子Kcs，将砾石覆盖影响因子（St）与Kf相乘，得到未考虑剖面砾石和有砾石覆盖影响的土壤可蚀性因子Kfs。

2 结 果

2.1 剖面砾石含量对土壤可蚀性的影响

比较未考虑剖面砾石和考虑剖面砾石时土壤可蚀性因子Kf和Kc发现，考虑剖面砾石含量（图1a）使土壤可蚀性因子值增加，全球范围内土壤可蚀性因子值增加0.001123（t·hm2·h）·（hm–2·MJ–1·mm–1），平均增加为5.43%。究其原因是因为土壤剖面中砾石的增加会降低土壤孔隙度，进而减少入渗，增加地表径流和土壤侵蚀产沙[9，23]。本研究计算结果表明，剖面砾石使土壤饱和导水率平均减少11.57%，使土壤渗透性等级增加5.68%（土壤入渗速率降低）。土壤可蚀性因子增加最多的地方在剖面砾石含量较高（土壤平均几何粒径较大）的撒哈拉地区（图1b和图1c），增幅大于20%；减小最少的地方见于北美大平原、亚马逊平原、拉普拉塔平原、东欧平原、西西伯利亚平原和华北平原，减幅<1%（图1d）。一般而言，影响较大的区域为山地和高原，包括：南北美洲的科迪勒拉山系（南美安第斯山、北美落基山脉）、加拿大东部拉布拉多高原、美国东部阿巴拉契亚山脉、东非高原和北部的撒哈拉沙漠、欧洲西北部斯堪的纳维亚山脉、西亚地区的伊朗高原、阿拉伯高原、青藏高原和天水、东西伯利亚山地、澳洲大陆的大部分地区（图1d）。这些地方的地形起伏度多大于200 m。与剖面砾石含量（图1a）对比表明，剖面砾石含量越高的地方影响越明显，特别是剖面砾石含量大于20%的地区。影响比较小的区域多见于平原，如美国中西部大平原、南美平原、印度平原、东欧平原、中国的东北和华北平原、及西伯利亚平原等地区，这些地方的地形起伏度多小于60 m。

将Kc与Kf的差值（图1c）与土壤类型图对比可见，当土壤类型为砂性土（Arenosols）、低活性淋溶土（Lixisols）、灰壤（Podzols）、薄层土（Leptosols）、疏松岩性土（Regosols）和铁铝土（Ferralsols）时，对土壤可蚀性因子影响较大。因为这类土壤发育较弱，剖面砾石含量和砂粒含量较高，平均值分别为10.511%和61.211%，黏粒含量和有机碳含量较低，均值分别为19.28%和23.14%，且大多分布在干旱地区或呈酸性，土壤比较松散[29]，所以较容易受到侵蚀。当土壤类型为硅胶结土（Durisols）、滞水潜育土（Stagnosols）、黏磐土（Planosols）、漂白淋溶土（Albeluvisols）和有机土（Histosols）时，对土壤可蚀性因子影响较小。这是因为这类土壤多分布在森林地带，剖面砾石含量和砂粒含量较低，平均值分别为5.185%和41.72%，黏粒含量和有机碳含量较高，平均值分别为24.44%和71.60%。

2.2 砾石覆盖对土壤可蚀性的影响

全球范围内而言，砾石覆盖使土壤可蚀性因子衰减系数平均为0.911，山地高原荒漠地区可蚀性衰减系数St的均值分别为0.641和0.681，使土壤可蚀性因子降低约8.9%和18.7%。砾石覆盖使K值和土壤侵蚀减小的原理在于，砾石覆盖减缓了雨滴打击，减少了地表水流速度，降低了径流输移能力，进而减轻土壤侵蚀[17-18，20]。由砾石覆盖度和砾石衰减系数图（图2）可以看出，砾石覆盖较大的区域主要见于山地和高原，如南北美洲的太平洋沿岸山地（科迪勒拉山系）、加拿大东部拉布拉多高原、西亚地区的伊朗高原、青藏高原及澳大利亚零星地区。与前述剖面砾石影响比较，范围较小。这些地方土壤表面砾石覆盖度较大，平均值为21.16%。而其他地区（平原、波状平原和俄罗斯远东地区的缓坡丘陵等）则影响轻微，砾石衰减系数一般为0.9左右。就土壤类型看，影响比较大的有暗色土（Umbrisols）、薄层土（Leptosols）、疏松岩性土（Regosols）和寒冻土（Cryosols）；St平均值较低，衰减系统介于0.57～0.63，平均为0.60（表3）。这组土壤发育程度较低，粗颗粒成分和地表砾石覆盖度较大[29]，平均砾石覆盖度为16.13%，所以对St的影响更大。而其他地区（平原、缓坡丘陵等）则影响轻微，土壤类型为聚铁网纹土（Plinthosols）、铁铝土（Ferralsols）、黏绨土（Nitisols）和黏磐土（Planosols）等土壤，衰减系统介于0.84～0.92，平均为0.88；这组土壤发育相对较好、砾石覆盖度较低，平均为4.31%。

表3 几种土壤类型的St均值 Table 3 Mean of St（soil erodibility reduction factor）relative to soil types

2.3 剖面砾石含量和砾石覆盖的总影响

前述分析可见，剖面砾石和砾石覆盖对土壤可蚀性因子均具有影响，但方向相反，这与前期研究的认识一致[15，30]。对不考虑剖面砾石和砾石覆盖、只考虑剖面砾石、只考虑砾石覆盖、同时考虑了剖面砾石和砾石覆盖影响的土壤可蚀性因子计算（Kf、Kc、Kfs、Kcs）的统计（图3）表明，剖面砾石使土壤可蚀性因子增加（增幅0.001 1（t·hm2·h）·（hm–2·MJ–1·mm–1））、砾石覆盖使土壤可蚀性因子减小（减幅0.0024（t·hm2·h）·（hm–2·MJ–1·mm–1））、当剖面砾石和砾石覆盖均考虑时总体上使土壤可蚀性因子减小（减幅0.001（t·hm2·h）·（hm–2·MJ–1·mm–1））。土壤可蚀性因子降低区域占总面积62.7%，增加区域占总面积31.1%，减少和增加的均值分别为0.006 4（t·hm2·h）·（hm–2·MJ–1·mm–1）和0.000 9（t·hm2·h）· （hm–2·MJ–1·mm–1），其中明显增减的面积分别约占20%。这说明从总体上看，砾石覆盖的作用更大，如果不考虑剖面砾石和砾石覆盖，则总体使土壤可蚀性因子被高估。不考虑剖面砾石和砾石覆盖的土壤可蚀性因子和考虑了砾石含量的土壤可蚀性因子差值图（图4）可以清晰地表达粗颗粒物质的综合影响。砾石覆盖影响为主的情形（相对于Kf，土壤可蚀性因子值减小，图面上的紫色区域）见于山地、高原和戈壁，地理区域详见3.1，这些区域土壤可蚀性因子减小的平均值为0.009 1（t·hm2·h）· （hm–2·MJ–1·mm–1）；主 要 的 土 壤 类 型 为 薄 层 土（Leptosols）、高活性淋溶土（Luvisols）、钙积土（Calcisols）、灰壤（Podzols）等。剖面砾石影响为主的情形（相对于Kf，土壤可蚀性因子值增加，图面上的绿色区域）见于土壤剖面砾石含量比较高的地区，但这些地区零星分布，主要见于北美的东北部、南美的东部和西北部、非洲的东部、西欧以及澳洲大陆的北部；这些区域土壤可蚀性因子增加的平均值为0.001 9（t·hm2·h）·（hm–2·MJ–1·mm–1）。主要的土壤类型为砂性土（Arenosols）、铁铝土（Ferralsols）、雏形土（Cambisols）和寒冻土（Cryosols）等。基本未影响的情形（相对于Kf，土壤可蚀性因子值变幅绝对值小于0.000 1，图面上的灰色区域）见于平原，如美国中西部大平原、南美平原、印度平原、东欧平原、中国的东北和华北平原、及西伯利亚平原等地区，面积占6.18%。土壤类型主要有有机土（Histosols）、低活性强酸土（Acrisols）、漂白淋溶土（Albeluvisols）等。这些地方基本上为平原，土壤剖面砾石含量（12.3%）、和砾石覆盖（12.7%）含量小，对土壤可蚀性因子没有明显的影响。

2.4 砾石含量对土壤侵蚀评价的影响

6个样区的分析表明（表4），以砾石覆盖为主的情况如下（SS1, SS3和SS6样区），土壤剖面砾石的作用使土壤水蚀流失速率略增加（幅度为2.3%～4.7%，平均为2.6%），土壤表面砾石覆盖和剖面砾石共同作用使流失速率减少17%～28.4%、平均为22.8%（土壤水蚀流失速率减小181 t·(km–2·a–1)）；由此推算单独砾石覆盖使之减少19.5%～31.6%，平均为27.3%。以剖面砾石作用为主的样区（SS2, SS5），剖面砾石的作用使流失速率略有增加，幅度为3.7%～10.5%、平均为7.0%，砾石覆盖和剖面砾石共同作用下流失速率增加2.1%～7.6%，平均为4.85%（土壤水蚀流失速率增加7 t·(km–2·a–1)），单独砾石覆盖使之减少1.6%～2.9%，平均为2.3%。两种作用均有的一个样区（SS4）位于澳大利亚，剖面砾石的作用使流失速率增加17.1%，共同作用减小8.4%（土壤水蚀流失速率减小2 t·（km–2·a–1）），推算砾石覆盖单独减少25.5%，剖面砾石的作用略大。全球范围内，剖面砾石作用使土壤水蚀流失速率增加4.8%，砾石覆盖和剖面砾石共同作用使流失速率减少14.7%（土壤水蚀流失速率减小54 t·(km–2·a–1）），由此推算单独砾石覆盖使之减少19.5%。

表4 典型样区砾石含量对土壤水蚀流失速率的影响 Table 4 Effects of rock fragment on soil erosion rate in the sample areas

3 讨 论

3.1 砾石覆盖

文献报道（如Miller和Guthrie[10]；Weil 和 Brady[28]）、野外观察和1 cm分辨率高分无人机影像（图5a）均表明，丘陵山区土壤广泛存在砾石覆盖。在戈壁地区地表砾石覆盖度可达到50%[10，31-32]，黄土地区碳酸盐结核（俗名料姜石）出露地表，也属于类似情况[33]。美国西部大约有17%的土系砾石含量大于35%，所占面积超过15%[10]。所以，尽管土壤侵蚀模型研究者注意到了这个问题（如毕小刚[34]；Römkens[9]），但由于缺少覆盖大区域的砾石覆盖度数据（类似于 LUCAS database，Panagos[23]），目前最高分辨率的卫星影像也无法提取这些信息（图5b），所以长期以来的K值计算大多数未考虑砾石覆盖的影响。

作为一种初步的尝试和无奈的研究方案，本文用土壤表面砾石含量代替砾石覆盖，虽然实现了对砾石覆盖影响的分析，但因为土壤学中所用数据和砾石覆盖的概念不完全相同[10，27]，因此本文分析的土壤表面砾石覆盖，未考虑大石块（cm级甚至更大）的覆盖，理论上讲是不严密的。与Panagos等[23]的分析相比（图6），一方面其空间格局相似、并且总影响与欧洲地区的接近（土壤侵蚀速率减少16%左右），一方面本文的结果偏大。今后可通过改善数据（如典型地区野外调查、查阅土壤普查报告等），分析土壤普查和制图数据中对于砾石相关的土相的记载[27]，逐步改善砾石覆盖含量对土壤可蚀性影响的全面估算。

3.2 关于相关K因子算法的修正

本文引言部分提到，土壤可蚀性因子的估算有多种算法，如USLE-K算法[1]、RUSLE2-K算法[8]、EPIC-K算法[3]和基于平均几何粒径的算法[9]等。本文作为一种方法的探讨，以USLE-K算法为例进行分析，形成的计算方法对其他算法基本适用。砾石覆盖计算，只需要砾石覆盖度数据，用式（9）完成砾石覆盖下的可蚀性衰减系数St计算，然后用式（11）可完成对任一可蚀性因子计算结果的砾石覆盖的修正，这一过程与土壤可蚀性因子算法无关。但是对于剖面粗砂影响的修正，因为需要通过土壤饱和导水率—土壤渗透性（式（6）和式（8））来推算K值的增减，所以式（11）对EPIC算法、基于几何平均粒径的算法还不能直接应用，需要进一步研究。

3.3 砾石含量对K的影响机理

砾石覆盖对土壤产流、侵蚀产沙的影响比较复杂。一些研究与本文的计算结果并不一致，如Zavala等[35]的研究表明，砾石存在会增加地表积水、阻滞产流、增加入渗、降低土壤流失速率。王蕙等[13]的研究表明，土壤中砾石的存在将对整个水文循环产生影响，嵌套于土壤中的砾石可以增加地表入渗，阻碍坡面侵蚀。同时，还有试验研究表明，砾石（包括土壤剖面中的和土壤表面覆盖的）影响，需要考虑的因素较多，包括砾石的位置和形状、土壤孔隙类型、土壤表面坡度和侵蚀的微形态特征等[11，30，36]。 这些研究均表明，砾石覆盖对土壤侵蚀的影响比较复杂，有待从侵蚀产沙物理过程出发集成更多的试验研究数据，全面评估砾石覆盖对入渗、产流和侵蚀产沙的影响，进而准确评价其对土壤可蚀性因子的影响。

3.4 砾石含量对风蚀的可蚀性影响

研究表明砾石覆盖对风蚀速率有重要影响[32，37-38]。一些风蚀预报模型，如风蚀方程（Wind Erosion Equation，WEQ）[39]、修正风蚀方程（Revised Wind Erosion Equation，RWEQ）[40]、风蚀预报系统（Wind Erosion Prediction System，WEPS）等[41]，均需要土壤可蚀性因子（soil erodibility factor）。尽管风蚀模型中的土壤可蚀性因子与USLE的K因子的估算方法有所不同，但砾石覆盖、剖面砾石含量均对其有影响[38]。因此，如果将本研究的式（8）、 式（10）和式（11）做必要的修改，有可能用于风蚀预报方程WEP的土壤可蚀性因子计算、提高风蚀预报的精度。

4 结 论

基于全球1 km分辨率土壤属性数据SoilGrids的计算表明，土壤剖面中砾石的存在，使土壤渗透性等级增加5.68%，土壤饱和导水率降低11.57%，全球土壤可蚀性因子增加4.43%；砾石覆盖通过保护土壤免遭雨滴打击和径流冲刷减少土壤侵蚀，在山地、荒漠（沙漠和戈壁）地区，这一影响会使土壤可蚀性值减小约18.7%；土壤剖面砾石含量和土壤表面砾石覆盖的综合影响，使土壤可蚀性因子降低5.52%。以砾石覆盖影响为主的地区，土壤可蚀性因子可降低0.0091（t·hm2·h）·（hm–2·MJ–1·mm–1）；以剖面砾石影响为主的地区，土壤可蚀性因子增加0.0019（t·hm2·h）·（hm–2·MJ–1·mm–1）；初步计算砾石含量可使全球土壤流失速率减少约54 t·km–2·a–1（14.7%）。在土壤可蚀性因子（K）计算和制图中充分考虑砾石含量的影响，计算所得土壤可蚀性因子会更加准确；同时，本研究对风蚀预报也具有参考价值。

致 谢中山大学上官微副教授对本文数据应用提供了指导意见，特此致谢！