-
开放科学(资源服务)标识码(OSID):
-
土壤侵蚀是自然界普遍存在的现象，也是我国面临的重大生态安全问题。土壤侵蚀导致土壤退化，土壤质量下降，易引发严重的自然灾害[1]。土壤侵蚀模数是衡量土壤侵蚀强度的重要指标[2]，通常定义为一定时间内单位面积土壤流失的量[t/(km2·a)]。准确估算土壤侵蚀模数对于评估土壤侵蚀风险、制定防治措施以及实现土地资源的可持续利用具有重要意义。在缺乏长时间序列径流观测小区的地区，学者们往往通过遥感半定量评估或模型的方式开展土壤侵蚀模数估算[3-5]。土壤可蚀性K值是现行主流模型中敏感且基础的因子，比如在通用土壤流失方程(Universal Soil Loss Equation，USLE)、修正通用土壤流失方程(Revised Universal Soil Loss Equation，RUSLE)、中国土壤流失方程(Chinese Soil Loss Equation，CSLE)中，土壤可蚀性K值是量化土壤是否容易受到侵蚀的一个重要因子[6]，是土壤侵蚀研究中的焦点之一。基于此，学者们对K值的估算方法开展了大量的研究。张宪奎等[7]、张爱国等[8]、刘宝元等[9]、张科利等[10]学者根据我国的具体情况，提出我国土壤可蚀性K值的测定方法，并将国外的USLE和RUSLE模型在我国的适用性进行了分析。总结当前国内外K值的估算方法，主要有诺莫方程[11]、环境政策集成气候(Environmental Policy Integrated Climate，EPIC)模型[12]、Sharri模型[13]、Torri模型[14]、径流小区实测法[9]等5种方法。其中，EPIC模型在西藏土壤可蚀性K值估算中被学者们认可[15-16]。
K值的大小受多种因素影响，包括土壤的物理性质、化学性质、环境条件等。张科利等[17]发现K值受土壤物理性质的影响而呈现出季节性变化，K值在小雨强范围内变化更大，随雨强的由小变大，K值会逐渐减小；郑海金等[18]通过对不同坡地建立径流小区进行实地观察，发现地表植被覆盖度高的小区土壤可蚀性小；崔晓薇等[19]通过EPIC模型计算K值，分析环境因子对K值的影响，发现坡度、年降水量和归一化植被指数对K值的影响达到了显著水平。梁音等[20]基于第二次土壤普查数据采用加权法计算了全国范围内的K值，并获得了K值在全国的分布图；王小丹等[21]通过多边形连接(Polygon Connection，PC)的方法生成了西藏高原的K值分布图，并分析了K值的空间分布特征；Baskan等[22]使用普通克里格(Ordinary Kriging，OK)的方法，对半干旱地区K值的空间分布进行了分析。Yu等[23]通过考虑气候、地形、遥感变量等环境因素对K值的影响，运用随机森林(Random Forest，RF)模型的方法对西藏东南部K值的空间变异性进行预测。以上研究为区域土壤可蚀性K值的方法选择和估算提供了科学基础，但在刻画区域复杂环境下土壤可蚀性K值的精度上还有待提升。
西藏作为青藏高原的主体组成部分，生态安全战略地位极其重要，但由于其极高的海拔和中低纬度等特殊的自然地理环境，使得高原生态环境脆弱且敏感。近年来，伴随着全球气候变化和人类活动的多重影响，青藏高原面临草场退化、土地沙化等一系列生态环境问题，尤其是水土流失给高原生态环境带来了极大的挑战[24]。然而高原上的水土流失过程和模数等问题并不明晰，这给区域水土流失防治，生态环境保护带来了极大的困难。过去学者们在高原上开展了关于K值的研究工作，取得了大量成果[25-28]。特别是第二次全国土壤普查(后简称“二普”)为K值研究提供了大量基础数据，但大多是基于有限采样点，通过加权、克里金插值等方法来评估土壤可蚀性K值[29-30]，精度和尺度都存在一定问题。随着地理信息和大数据技术的发展，采用大数据模型，考虑环境协变量来描述目标对象的方法被广泛应用[31]，但不同模型之间算法构架存在较大差异，对于不同情景下的适应性需要甄别。此外，第三次全国土壤普查(后简称“三普”)工作已经展开，样本数量较多，覆盖面也大，但其是一种关于土壤全面的普查，并非土壤侵蚀潜在危险性的专题评估。因此，本研究在Yu等[23]的研究基础上增加不同模型的对比和模型的不确定性估计，对区域土壤可蚀性K值进行评估；在崔晓薇等[19]的研究基础上选择更加系统的环境协变量，研究其对区域土壤可蚀性K值的影响；在梁音等[20]的研究基础上选择对土壤可蚀性K值有影响的环境协变量加入到空间插值中，以提高结果的准确性。基于此，本研究选择EPIC模型估算土壤可蚀性K值，并对当前广泛应用的OK、RF、PC 3种模型估算区域土壤可蚀性K值的方法开展精度评估后，以影响成土的五大因子中的相应因素作为环境协变量，评价其对土壤可蚀性K值的重要性，然后对西藏土壤可蚀性K值进行空间估算。研究结果可为西藏地区土壤可蚀性K值的空间估算、基于“三普”成果的土壤可蚀性评价、区域侵蚀危险性评估和防治提供理论参考。
HTML
-
本研究内容中涉及的样点一共有595个，每个样点的数据包含了其土壤机械组成、有机质含量和样点的海拔、坡度、地形起伏度、年平均气温、年平均降水和年归一化植被指数(Normalized Difference Vegetation Index，NDVI)等数据。其中包含“二普”数据242个样点，虽然“二普”工作已过去30多年，但考虑到高原绝大部分区域受人为因素干扰较弱，且土壤理化性质一般都相对稳定，所以选用该数据能很好地弥补现有采样工作的不足，且在“三普”数据未正式公布前，在短时间内获取西藏自治区全部区域的土壤理化性质数据非常困难。鉴于此种情况，本研究加入2023-2024年间实地采样点(353个)数据进行补充，实地采样点数据来源为本团队采样数据，基本涵盖西藏自治区范围(图 1)。
实地采样数据在2023-2024年采集，样点主要通过ArcGIS Pro软件对西藏不同地区土壤类型进行选择，以此为依据，选取西藏具有代表性的土壤类型，基本覆盖西藏范围，均匀采样。由于土壤侵蚀均发生在地表位置，所以采样点的土层深度为0~20 cm。采样过程中，设置1m×1m的样方，采用5点混合采样法，采用铁锹采集散土1 kg，用于分析土壤粒径组成；通过环刀法采集环刀土样3个和原状土1 kg，用于测定土壤的物理性质。土壤机械组成采取吸管法测定；土壤有机质含量采用重铬酸钾外加热法测定[32]。采用Microsoft Excel 2019进行数据整理；IBM SPSS Statistics 24.0软件和GS+软件用于地统计分析；R Studio 4.1.1和ArcGIS Pro软件用于K值的空间预测和制图。
-
本研究采用Sharpley等[12]开发的EPIC模型，根据土壤机械组成和有机碳数据计算每个土壤采样点的K值。公式如下：
式中：土壤粒径分级采用国际制；WSAN为砂粒(0.05~2.0 mm)比例(%)；WSIL为粉粒(0.002~0.05 mm)比例(%)，WCLA为黏粒(＜0.002 mm)比例(%)，WC为有机碳比例(%)；WSN1=1-WSAN/100。
由于“二普”资料中用于土壤粒径数据分析的标准为国际制单位，因此使用EPIC模型计算土壤可蚀性K值之前需将土壤粒径数据转化为国际制，即＞2 mm为石砾、0.1~2 mm为砂粒、0.05~0.1 mm为极细砂、0.002~0.05 mm为粉粒、＜0.002 mm为黏粒。本研究采用样条插值法，用MATLAB软件进行土壤粒径转换。通过EPIC模型计算之后，再将K值统一使用0.131 7转换系数转换为国际制进行数据计算。
-
本研究选取3种模型对K因子的空间分布进行预测，分别是：OK、RF、PC。其中，OK模型[33]是在非采样位置进行空间插值的方法。该模型比较灵活，允许使用研究空间自相关的变量[34]。RF模型是Breiman[35]提出的由多个决策树集成的机器学习算法，目前广泛应用于数字土壤制图，以建立土壤与环境之间复杂的定量关系。PC模型[36]是使用全国土壤发生分类系统展开的一种制图方法。首先获取土壤的基本信息，然后通过发生系统的5个基本等级，利用面积加权平均法进行逐级归并，最后对计算得到的K值进行拼接处理，得到土壤可蚀性K值分布图。
-
使用分位数回归森林(Quantile Regression Forest，QRF)的方法对预测模型的不确定性进行评估，具体方法为：使用第95个百分位预测值和第5个百分位预测值之间的差值，除以第50个百分位预测值来量化预测模型的不确定性[36]。不确定性值越小，表明对模型的预测效果越好。不确定性值越大，表明对模型的预测效果越差。
式中：Kuncertainty表示K值的不确定性；Kupper为第95个百分位的预测值；Klower为第5个百分位的预测值；Kn为第50个百分位的预测值。
-
为了评估模型的精度，采用4/5的数据作为训练集，1/5的数据作为验证集；采用决定系数(R2)、均方根误差(Root Mean Square Error，RMSE)、平均相对误差(Mean Relative Error，MRE)、平均绝对误差(Mean Absolute Error，MAE)4个指标对模型的预测值和观测值进行分析[37]。公式如下：
式中：Zi为K的实测值；Zi′为观测K的平均值；
$ \overline{Z_i}$ 为K的预测值；R2值越大，RMSE、MAE和MRE值越小，表明K的预测误差越小，精度越高。 -
环境协变量基于SCORPAN(土壤、气候、地形、母质物质、年龄、空间)的理论框架进行选择[38]。该框架通过总结前人对土壤预测建模和数据获取的方法，提出具有一定空间自相关误差的数字土壤制图的方法。具体来说，就是通过对前人研究的空间方法进行总结，通过已知土壤性质预测未知土壤属性。其中K值的计算与土壤机械组成和有机碳含量密切相关，因此选择的协变量大多与土壤颗粒的形成、转运和储存有关，从而保证具有事实依据。环境协变量如表 1所示进行了分类。由于数据来源不同，部分数据分辨率难以统一。为确保数据之间的空间匹配性，为后续分析提供可靠基础，本研究将分辨率较低的数据统一进行重采样处理，调整为30 m分辨率进行空间运算。
-
本研究使用Pearson法来分析K值与环境协变量的相关性，使用RF模型中的均方误差(%IncMSE)指标评估特征重要性，该方法通过打乱某一特征的值，观察其模型预测误差变化，从而量化每一个特征值对模型性能的贡献大小[39]。均方误差值越大，表明环境协变量对K值的预测结果越重要。值越小，表明环境协变量对K值的预测结果越不重要。
1.1. 数据采样与分析
1.2. K值的估算
1.3. K值空间分布预测模型
1.4. 模型不确定性估计
1.5. 模型精度验证
1.6. K值空间预测的环境协变量选择
1.7. 环境协变量与K值的关系
-
由表 2可知，K因子观测值和有机碳含量分别为0.010 4~0.060 5 t·hm2·h/(hm2·MJ·mm)和1.160~180.456 g/kg。K值的变异系数(CV)为27.0%，属于中等程度变异，说明土壤可蚀性K值在区域内存在一定的空间变化[40]。K值的偏度和峰度分别为0.689和0.123，表明K值呈右偏分布，数据整体偏向K值较低的一侧，低可蚀性土壤在西藏分布较为广泛。峰度接近0，说明数据整体接近正态分布，结果在研究区域内分布相对均匀，不存在极端集中或离散的情况。有机碳的变异系数为84.7%，表现为强变异性，表明有机碳在空间上变化较大，该结果与区域内土壤类型、土地利用类型以及环境条件等因素密切相关。黏粒、粉粒和砂粒的平均值分别为51.214%、40.123%和10.123%，其范围分别为2.217%~97.393%、0.300%~81.333%和0.115%~76.800%，三者含量在空间上不同位置的分布差异较大，数据的离散程度高。变异系数分别为39.5%、46.0%和82.0%，表明黏粒和粉粒具有较弱程度的变异性，而砂粒则是具有中等偏强程度的变异性。
利用SPSS软件，对研究区595个采样点的K值进行正态分布检验，结果表明，K值的平均值为0.029 1 t·hm2·h/(hm2·MJ·mm)，标准差为0.008 5 t·hm2·h/(hm2·MJ·mm)，K值总体服从正态分布，结果如图 2所示。说明本研究所得样本数值可靠，能够满足地统计分析的要求。
-
表 3和图 3分别展示了3个模型十倍交叉验证的结果，其中，OK模型已被证明在地形相对平坦、土壤类型分布相对均匀的区域时，精度较高。而西藏地区地形复杂，植被类型多样，导致OK模型的准确性最低[41]，R2为0.311，RMSE、MAE和MRE值分别为0.008、0.006和0.003；PC模型的准确度较OK模型有所增加，其R2为0.654，RMSE、MAE和MRE值分别为0.006、0.004和0.002，PC模型的使用依赖于全国土壤发生系统，通过土壤类型面积加权平均法由小到大逐级回归得到K值分布图[42]，由于西藏土壤类型分布相对均匀，与土壤发生分类系统匹配度较高，因此该模型能够较好地发挥作用。与前两种方法相比，RF模型作为一种机器算法，能够建立起土壤与环境变量之间的复杂关系[43]，在西藏复杂的环境下仍能显示出最高的精度。RF的R2为0.781，RMSE、MAE和MRE值分别为0.004、0.003和0.001，表明模型解释了78.1%的K值在西藏地区的空间变异，其中MRE的值接近于0，表明模型几乎无偏差。3种模型中，RF模型精度最高。因此，本研究运用RF模型对西藏自治区土壤可蚀性K值进行空间预测估算。
-
图 4是对K值、土壤质地、有机碳和环境协变量之间进行Pearson相关性分析，并进行可视化的结果，K值与环境协变量的相关性结果显示，K值与DEM呈显著正相关，p＜0.001；K值与S、TRI、MAT、NDVI、NSH、NP、NL均呈显著负相关，p＜0.05。其中K值与TRI的负相关关系最明显，p＜0.001。AS、MAP、PD、RVI和DVI 5项指标与K值之间未呈现出明显的相关关系。其中，AS对太阳辐射、植被生长有着直接影响，阴坡低温潮湿有利于有机质积累，阳坡水分蒸发较多，植被覆盖少，导致土壤侵蚀风险增加，且受流域内的地形影响，迎风坡比背风坡的降雨量大[44]。西藏地区以高山为主，其海拔和坡度的变化对水热条件的影响远大于坡向对水热条件的影响，故坡向与K值没有呈现出相关关系；人类活动主要通过破坏地表植被，改变土壤结构以影响K值，而西藏PD仅为3.6人/km2，且集中于雅鲁藏布江流域，导致广阔的藏北草原PD＜1人/km2。故PD在西藏范围内与K值未能呈现出显著相关。
此外，本研究还计算了DEM、NP、NSH、NDVI、MAP、MAT、TRI之间的相关性。其中，DEM与NP、NSH、MAT和MAP等气候因子之间均呈现显著负相关，其中与MAT的相关性最大，p＜0.001，与MAP的相关性最小，p＜0.01。DEM与NDVI和TRI 2项植被指标也均达到了显著负相关关系，p＜0.001。
由以上结果可知，K值与所选环境协变量中的地形、气候和植被变量之间显著相关。
-
使用均方误差(%IncMSE)计算环境协变量对预测K值的相对重要性(图 5)。其结果表明，所选环境协变量中，对预测K值影响最大的是NDVI，其%IncMSE达到了13.64%；其次是MAT、S和RVI，这3个环境协变量的%IncMSE＞9%；DEM、AS、NSH、TRI、DVI和MAP 6个环境协变量的%IncMSE＞7%。对预测K值影响最小的为PD数据，其%IncMSE为0.85%。地形因子的重要性占35.04%，植被因子的重要性占31.31%，气候因子的重要性占29.88%，人类活动的重要性占3.77%。结果表明，地形、植被和气候因子是影响K值空间变异性的主要环境因素，人类活动对K值的影响较小。
-
使用RF模型，将对其重要的环境数据(图 5)整合到预测的空间网格中。在此基础上，生成研究区预测K值的空间分布图(图 6)。将预测K值利用自然间段点分级法分为6个水平，分别是：高(≥0.036 38)、中高(0.030 15~0.036 37)、中(0.027 01~0.030 14)、中低(0.025 42~0.027 00)、相对低(0.022 28~0.025 41)、低(≤0.022 27)。研究区预测K值为0.016 05~0.048 70 t·hm2·h/(hm2·MJ·mm)。高可蚀性区域约为0.95×105 km2，占西藏地区总面积的7.90%，主要分布在藏西、藏中部分地区；中高可蚀性的面积为2.02×105 km2，占16.83%；中可蚀性的面积为2.87×105 km2，占23.86%；中低可蚀性的面积为2.49×105 km2，占20.75%；相对低可蚀性的面积为2.22×105 km2，占18.43%；低可蚀性的面积为1.47×105 km2，占12.23%。K值总体呈现由西至东降低的趋势。
以上结果表明，西部干旱荒漠区的年均降水量低，植被覆盖稀少，土壤发育不完全，团聚体稳定性较差，因此该地区K值较高。而在东南部湿润的林区范围内，降雨丰富，植被覆盖度高，土壤发育完全，土壤稳定性高，因此该地区K值较低。在中部高原湖盆区，其高原面地势平坦，广泛分布湖积物和冰渍物，兼具一定的保水能力和结构稳定能力。
-
图 7和图 8显示了预测K值在第5百分位和第95百分位的分布情况，两个值的取值范围分别为0.011 15~0.015 49和0.016 02~0.038 09。此外，由图 7和图 8可以看出K值的不确定性在两个百分位时具有一致性，都具有自西向东逐渐减小的趋势，总体西部和中部相对较高，东部、北部和中南部相对较低。图 9显示了不确定性指数的分布情况，其取值范围为0.271~0.816，其分布格局总体为中部、西部及北部地区的不确定性较高，东部、南部地区的不确定性较低。
结果表明，高不确定性区主要分布于西藏的中部、西部及北部地区，由于该区域属于西藏无人区，且本研究仅有12个样点分布在该区域，远少于东部、南部地区，RF模型由于数据不足导致预测不确定性增大；中西部地形过渡带由于海拔的陡然爬升，导致年降水量增加，土壤类型变化巨大，因此模型预测不确定性增大。
2.1. K值与土壤性质观测值的描述性统计分析
2.2. K值预测模型的选择
2.3. 土壤可蚀性K值空间估算与不确定性评估
2.3.1. 环境协变量对K值的影响
2.3.1.1. K值与环境协变量的相关性
2.3.1.2. 环境协变量影响K值的重要性排序
2.3.2. 预测K值的空间分布特征
2.3.3. K值的不确定性估计
-
环境因子对K值的影响评价显示，地形、气候和植被是最主要的影响驱动因素，其中地形的贡献最大。弓小平等[41]、Yu等[23]与本研究的结论基本一致。地形变量中的海拔和坡度因子与K值的空间变化密切相关。其原因是，海拔控制地表水和热量的再分配，而地表水和热量影响着土壤颗粒的形成和转运，因此海拔也控制着K值[45]。此外，Gu等[46]在云南地区的研究发现，当坡度小于15°时，K值随坡度增大而增大，当坡度大于15°时，K值随坡度增大而减小，且减小趋势更为明显。因此，地形对土壤K值的影响值得深入研究。气候和植被一同决定了有机质的形成和储存[45]，影响土壤结构的演化与稳定，从而使得土壤有机质和土壤颗粒呈现出明显的环境变异特征，最终导致土壤可蚀性K值的空间差异性。人类活动对K值的空间变异影响最小，这一结论与西藏地广人稀、人类活动强度较弱的特点一致。尽管局部地区的放牧会对土壤产生一定的影响，但其范围和强度均有限，不足以影响整个K值的空间分布格局。
K值的空间预测结果显示，K值的高和中高侵蚀区、中和中低侵蚀区、相对低和低侵蚀区的面积分别占西藏面积的24.73%、44.61%和30.66%。王小丹等[3]研究结果中的侵蚀面积与本研究有一定的偏差，主要原因可能是本研究对通过RF模型预测K值的计算中考虑了环境协变量对K值的影响；前人研究多聚焦于西藏局部地区或特定地貌类型[23-28]，而本研究覆盖整个西藏范围，采样区域包含了前人未充分采样的藏西北高寒荒漠区、藏北高原等含复杂地貌类型地区，导致西藏整体K值与局部研究结论存在一定的出入。此外，本研究整合了“二普”数据和2023-2024年实地采样数据，而前人的研究多依赖单一的“二普”数据或短期采样数据进行分析。K值空间变异总体上呈现由西南至东南降低的趋势，东南部处于较低水平，中部地区处于中和中低水平。这是由于藏西南地区的气候主要表现出干燥、降雨量低、年均温度低等特点，导致地表物质更加破碎，其土壤存在很多粗骨砾，导致K值较高。藏东地区高可蚀性的原因是该地区有干热河谷的分布，主要表现出高温、低湿度、低降雨量等特点，导致该地区的K值较高。通过本研究，明确西藏高可蚀性区域的分布范围，可以重点部署、精准施策、避免资源浪费，为西藏的土地利用规划、水土流失防治措施制定等方面提供科学依据。
尽管本研究的样本点位数已经有将近600个，但是某些地区由于技术原因无法到达，导致交叉验证的结果可能存在一定偏差。例如，中北部由于没有样点的覆盖而导致RF预测的K值不确定性偏大。未来，针对北部高海拔无人区采样困难的问题，可以采用无人机遥感反演技术和地面验证结合的方法进行，优化采样方法，通过精细化采样、多源数据融合和建模分析，进一步优化西藏土壤可蚀性K值的空间预测精度。
-
1) RF模型在西藏土壤可蚀性K值空间预测中表现最优，其精度显著高于OK模型和PC模型；RF模型预测K值的取值范围为0.016 05~0.048 70 t·hm2·h/(hm2·MJ·mm)。RF的R2达到了0.781，RMSE为0.004。
2) 地形(35.04%)、植被(31.31%)、气候(29.88%)是影响K值的主要环境因子，人类活动贡献度仅3.77%，反映出西藏土壤可蚀性K值主要受到自然条件的影响。
3) K值空间预测结果显示，在空间上呈现“西南、东北部高，中部、东南部低”的分布格局，其高和中高侵蚀区、中和中低侵蚀区、相对低和低侵蚀区的面积分别占西藏面积的24.73%、44.61%和30.66%；K值空间预测不确定性分布格局总体表现为：中部、西部及北部地区的不确定性较高，东部、南部地区的不确定性较低。
DownLoad: