2025, Vol. 46
2. 新疆干旱区湖泊环境与资源重点实验室，乌鲁木齐 830017
2. Xinjiang Key Laboratory of Lake Environment and Resources in Arid Zone, Urumqi 830017, China
土壤有机碳（soil organic carbon，SOC）是影响土壤性质和结构，反映土壤质量和健康状况的关键指标[1]. 作为陆地生态系统中最大的碳库，在1 m深度范围内的土壤有机碳含量约为生物圈碳库总量的3倍，大气碳库总量的2倍[2]，准确估测土壤有机碳的含量及其变化对全球碳循环和气候变化以及实现碳达峰和碳中和目标起着极其重要的作用[3，4]. 随着遥感技术和机器学习技术的发展，数字土壤制图技术被广泛地应用到土壤有机碳含量的估测研究[5，6].
利用遥感技术可快捷方便地以较低成本获取丰富的土壤属性信息，同时遥感影像的易获取性和时效性也为土壤有机碳含量的准确估测提供了丰富的数据来源. 相较于单幅遥感影像来说，利用时间序列所得到的合成影像能够包含更加丰富的影像信息，并且一定程度上能够减少土壤水分或其他干扰因子对土壤有机碳含量估测带来的影响，提高土壤有机碳含量的估测精度[7]. Yang等[8]利用NDVI时间序列提取的物候参数对土壤有机碳含量进行估测并将预测精度提高了50%. Wang等[9]利用合成影像与单一影像对比发现多时相合成影像能够提高遥感影像波段和土壤有机碳含量的相关性与建模精度. Shi等[10]对中国东北和比利时两个区域进行的多时相合成影像SOC含量模型效果始终比单期影像模型有更好的估测精度. 利用遥感影像获取的地物光谱信息虽然能够间接反映土壤信息，但土壤有机碳含量还受到气候、地形、母质和生物等多种因素的影响[11]. 祝元丽等[12]将多个环境变量应用到土壤有机碳含量的估测研究，发现气温和降水等因素是影响土壤有机碳含量的重要因素. 闫坤等[13]通过研究表明高程和地形湿度指数等地形因子在土壤有机碳含量估测中发挥着重要作用. 为了提高估测模型的精度和效率，变量筛选算法也是必要的. 作为一种简单实用的变量筛选算法，由Boruta算法筛选得到的特征变量构建的模型精度更高[14]. 土壤有机碳含量与建模变量之间存在着较为复杂的非线性关系，利用机器学习算法可以有效建立土壤有机碳含量与变量之间的联系. Keskin等[15]通过对比随机森林（random forest，RF）、支持向量机以及偏最小二乘等8个模型发现，随机森林模型对土壤有机碳含量具有最高的估测精度和稳定的估测效果. 由于土壤有机碳在空间上存在较大的空间异质性和区域特征，导致当气候、地形和生物等环境条件发生改变时，土壤有机碳含量也相应发生变化. 因此，利用多时相合成影像提取的光谱变量与环境变量构建土壤有机碳含量的估测模型仍需进行深入探讨.
位于中国西北干旱区的新疆博斯腾湖流域，在气候变化和人类活动影响下，其土壤环境要素较为敏感. 目前不同合成影像应用于干旱区土壤有机碳含量估测的精度和效果仍需深入研究，且不同变量组合对土壤有机碳含量的估测精度还有待进一步提高. 基于以上分析，本文以新疆博斯腾湖流域土壤有机碳含量为研究对象，按照哨兵2号（Sentinel-2）多时相卫星影像合成方法生成光谱变量并结合气候和地形等环境变量，通过Boruta算法筛选最优变量组合，利用随机森林模型估测研究区耕层土壤有机碳含量，以期为博斯腾湖流域陆地生态系统碳储量估算提供依据.
1 材料与方法 1.1 研究区概况博斯腾湖流域（简称博湖流域）位于天山中段南麓，在行政区划上属于巴音郭楞蒙古自治州，流域主要包括和静县、和硕县、博湖县和焉耆回族自治县等市县的部分区域. 博湖流域的海拔高度介于1 040.28~4 812.65 m之间，地势西高东低，整体上自西北向东南方向倾斜，地形主要以山地和盆地为主（见图1），流域内景观随着海拔高低从上到下依次为积雪带、草甸草原带、绿洲及荒漠带等[16]，土地利用类型主要由耕地、草地、荒漠和水域构成. 由于深居内陆且受高山阻挡使得水汽无法到达，流域内总体呈现出温差较大且降水稀少的特征，年平均气温在8~8.6℃之间，年降水量在70~600 mm之间，属于温带大陆性干旱气候[17]. 流域内主要以种植业、畜牧业和林业为经济发展的支柱性产业，主要粮食作物包括玉米、小麦、番茄和辣椒等，此外还种植有棉花、甜菜和向日葵等多种经济作物.
|
图 1 博斯腾湖流域土壤样点分布示意 Fig. 1 Soil sampling sites distribution in the Bosten Lake Basin |
Sentinel-2中的L2A级产品凭借其涵盖多个波段、最高10 m的空间分辨率、5 d重访周期的高时间分辨率、免费获取以及无需进行辐射定标及大气校正等多个优点被广泛应用到土壤有机碳含量估测研究中，且取得了较好的效果[18]. 利用谷歌地球引擎（Google Earth Engine，GEE）平台选择与研究区内相同采样月份的T44和T45投影网格区域，瓦片编号涵盖TPN、TQN、TUH、TVH、TVG、TWH和TWG的共63景Sentinel-L2A遥感影像. 分别利用最小值合成（min）、中值合成（median）以及均值合成（mean）这3种影像合成方式，通过 QA60和SCL波段对存在云的区域进行识别和掩膜去除，并分别进行影像拼接得到高质量的合成遥感影像数据. 由于存在云的区域主要集中在研究区中部及东北山地，人类活动痕迹少，所以不存在地物变化等情况. 因此对于存在云的区域，通过选择相同年份相邻月份以及相邻年份相同月份的遥感影像数据对云量较高的区域进行填充以达到去云的效果. 对经过预处理后的影像选择空间分辨率10 m的蓝（B2）、绿（B3）、红（B4）和近红外（B8）波段，以及20 m分辨率的红边1（B5）、红边2（B6）、近红外窄波段1（B7）、近红外窄波段2（B8A）、短波红外1（B11）和短波红外2（B12）波段. 为了统一尺度，使用最邻近插值法将20 m空间分辨率的6个波段重采样到10 m分辨率. 利用影像B2~B8和B8A、B11以及B12等波段数据以及不同波段运算得到的光谱指数，包括归一化植被指数（normalized difference vegetation index，NDVI）[19]、增强型植被指数（enhanced vegetation index，EVI）[19]、比值植被指数（ratio vegetation index，RVI）[20]、差值植被指数（difference vegetation index，DVI）[21]、土壤调整总植被指数（soil-adjusted total vegetation index，SATVI）[22]、转换植被指数（transformed vegetation index，TVI）[23]、核归一化植被指数（kernel normalized difference vegetation index，KNDVI） [24]、绿光归一化植被指数（green normalized difference vegetation index，GNDVI）[25]、绿红植被指数（green-red vegetation index，GRVI）[26]、比例差植被指数（ratio difference vegetation index，RDVI）[27]和土壤调整植被指数（soil adjusted vegetation index，SAVI）[28]共11个植被指数，以及归一化盐分指数（normalized differential salinity index，NDSI）[28]、水分胁迫指数（moisture stress index，MSI）[29]和地表水指数（land surface water index，LSWI）[30]3个土壤及水分指数，共14个光谱指数对土壤表层有机碳含量进行估测. 光谱指数的具体计算公式见表1.
|
表 1 光谱指数计算公式1） Table 1 Spectral index calculation formula |
1.2.2 环境数据获取
除了生物因素外，气候和地形对土壤有机碳的空间分布影响较大. 气候通常能够决定SOC空间分布的格局，其他环境因素则导致SOC含量在小范围内产生变化[31]. 本研究中收集了包括气象因子和地形因子在内的环境变量作为SOC含量估测的建模变量. 其中气候变量主要包括气温、降水以及近地表气压等气象因子. 地形变量主要包括高程、坡度、坡向和河谷深度等，所有的环境变量及其缩写见表2. 选择由欧洲中期天气预报中心开发的ERA5数据集作为本研究中的气候变量. ERA5数据集具有长时间序列以及丰富的气候参数，可详细记录气温、降水、地表压力和地表径流等特点. 本研究中选择2 m高度平均温度、土壤温度、总降水量、总蒸散量、裸土蒸散量、地表压力和径流量总和这7个变量，数据来源于GEE平台. 由于ERA5数据的空间分辨率较低，因此在ArcGIS 10.7版本中，根据ERA5数据的空间分辨率创建对应的渔网点数量并提取气候变量点值，利用普通克里金插值法得到10 m空间分辨率的气象因子数据. 研究中的地形数据来源于哥白尼数字高程模型（Copernicus DEM，COP-DEM），该数据是由欧洲航天局（European Space Agency，ESA）所发布，具有精度高、地形细节表现准确、数据现势性强以及开源免费获取等优势，被认为是最佳全球性开源DEM数据之一[32]. 在GEE平台获取COP-DEM作为输入数据，利用SAGA-GIS软件，计算得到坡度、坡向、谷地深度和地形湿度指数等共15个地形变量. 由于地形变量和影像光谱数据具有不同的空间分辨率，因此，在GEE中对所收集的地形因子按照最邻近采样法重采样到10 m.
|
|
表 2 环境变量及缩写 Table 2 Environmental variables and abbreviations |
1.2.3 土壤有机碳含量数据获取
在干旱区，植被可以有效缓解土壤养分的流失，对维持土壤肥力具有重要作用. 因此本研究选择植被生长旺盛期的2023年7月中旬，在博斯腾湖流域进行土壤样品采集，以反映土壤有机碳含量在高植被覆盖时的状况. 由于正处植物生长期，样方和影像中并不完全是纯土壤像元. 采集土壤样品时，在避免人类和动物活动干扰的前提下，选择10 m×10 m的区域，根据梅花土壤采样法采集5个耕层（深度为0~20 cm）土壤样本，通过混合得到最终样本. 在采样的同时，使用手持式GPS定位仪记录每个采样点的地理坐标，获得实际土壤样点空间分布（图1）. 由于土地利用类型的差异将直接影响SOC含量，因此采样时分别采集了58个耕地以及34个荒漠林草地样品，共92个土壤样品. 对采集到的所有样品进行自然风干、去除杂质、粉碎以及过2 mm筛，并对处理好的土壤样品按照重铬酸钾容量-外加热法测定SOC含量.
由于流域范围较大且地形复杂，部分区域难以进行现场采样. 为了全面地对研究区域进行较为准确的SOC含量估测，需要对存在数据空白的区域进行填充插补. 本研究分别利用国际土壤参考和资料中心SoilGrids土壤数据库与OpenLandMap土壤数据库以及光谱和环境变量，与实测SOC含量进行多元逐步线性回归拟合插补. 计算实测样点与对应拟合样点的绝对误差和相对误差，结果如图2所示. 从中可以看出，利用土壤数据库进行插补拟合得到的样点误差（包括绝对误差和相对误差）通常要大于用光谱变量结合环境变量（即非土壤数据库）拟合得到的样点误差. 因此为了更准确地模拟流域SOC含量空间分布，在研究中根据遥感光谱变量、环境变量和实测SOC含量数据之间的关系建立多元逐步线性回归模型，按照渔网分割的方法对流域内采样困难区域插补样点29个，剔除实测样本点中的一个异常数值，共得到120个样点数据.
|
图 2 拟合样点对比误差 Fig. 2 Error of fitting sampling sites |
处理多个复杂的光谱和环境变量需选择合适的变量筛选方法，本研究采用Boruta算法进行变量筛选. Boruta算法通过创建所有输入变量的影子特征来增加变量集合的随机性，并在包括影子特征的新变量集合的基础上通过运行随机森林分类器对每个原始输入变量进行迭代. 根据模型结果准确率的升降来评估每个原始变量的重要性，如果原始输入变量的得分高于影子特征的最大得分，表明此变量在模型估测中发挥着重要作用并被保留，否则变量被认为对模型估测贡献较小从而被删除[33]. 这一过程持续循环直至所有原始变量都被保留或删除，或者达到模型所设定的最大迭代运行次数时，Boruta算法终止运行. 在本研究中，利用RStudio 4.2.2版本实现Boruta算法，运行时将doTrace设置为1，maxRuns设置为30 000以便能够获得较为详细的变量重要性相关信息.
1.4 样本划分方法为了对训练集和验证集进行更好地划分，发掘数据本身的真实可用信息，选择采用SPXY（sample set partitioning based on joint x⁃y distance）算法进行样本划分. SPXY算法以K⁃S（kennard-stone）算法为基础并对其进行了优化改善，使得能够选取更具有代表性的样本集. K⁃S算法以选取样本中欧氏距离最远的两个作为训练集的初始样本开始运行，接着依次迭代运算剩余样本与当前训练集中每个样本之间的最大和最小距离并进行选择，当达到规定样本数量时，算法终止运行. 在K⁃S算法的基础上，SPXY算法通过进一步探究样本特征变量与目标变量的关系，使得样本划分更加合理和全面[34]. 在Matlab软件中，按照训练集占总样本数量的70%，验证集占总样本数量的30%，对收集到的120个土壤样品数据按照SPXY方法进行划分.
1.5 建模反演方法随机森林（RF）模型是一种基于多决策树的机器学习模型，凭借其处理高维非线性关系的能力好、对噪声具有较强鲁棒性、参数设置简单以及模型效率高且效果较好等优点，已经被广泛地应用于SOC含量空间反演研究中[35~37]. RF模型通过多个决策树（ntree）和树中每个节点所选择的变量数量（mtry）来决定模型精度. 模型中的每棵决策树都依赖随机且相互独立的样本进行构建，通过将这些决策树整合到一起对所要预测的对象给出单独的预测，最终的预测结果是每棵树预测结果的平均值. 本研究使用RStudio 4.2.2版本构建RF模型，并利用RF包中的tuneRF函数选择袋外误差最小的mtry值. 通过绘制模型训练时误差率随ntree增加的变化曲线，找到使误差率趋于稳定的mtry的值. 通过上述方法发现在本研究中mtry值设置为输入变量的2/3，mtry值设置为100能够较好地表现模型的性能.
1.6 模型评价方法为了验证研究中所使用模型的估测效果，使用决定系数（R2）和均方根误差（root mean square error，RMSE）以及相对分析误差（relative percent deviation，RPD）这3个指标对模型估测精度进行评价. 其中R2和RPD值越大，RMSE值越小，所建立模型的估测效果就越好. 反之，模型估测效果越差. 当模型的RPD<1.4时，表示所建立的模型不适合进行估测；当1.4≤RPD<2时，表示该模型的估测能力一般，可以用来进行简单的估测；当RPD≥2时，表示模型的估测能力很好，可以用来进行较为准确的估测[38].
| (1) |
| (2) |
| (3) |
| (4) |
式中，n表示样本总数；
按照SPXY样本划分方法建立训练集和验证集，并分别对总样本以及划分后的2个样本进行SOC含量统计特征分析，结果如表3所示. 博斯腾湖流域样本点的ω（SOC）在0.974~54.471 g·kg-1之间变化，表现出较大的SOC含量差异. 总样本ω（SOC）平均值为20.226 g·kg-1，训练集的平均值为21.949 g·kg-1，验证集的平均值为16.205 g·kg-1，3个样本集之间的SOC含量平均值差异不大. 各样本集ω（SOC）的标准差在15.523~16.882 g·kg-1之间变化，表明每个样本集内部均存在较大的差异. 变异系数用来衡量样本内的空间异质性，结果显示3种样本集的变异系数在76.9%~95.8%区间，进一步说明博斯腾湖流域各土壤样本SOC含量的离散程度较高，在整个流域内呈现出较为明显的空间异质性.
|
|
表 3 SOC含量统计 Table 3 Soil organic carbon content statistics |
2.2 相关性分析
通过相关性分析筛选掉与SOC含量无显著相关的光谱和环境变量以提高模型的准确性和建模效率，3种光谱变量集以及环境变量与SOC含量的相关关系如图3所示. 从中可以看出，光谱变量B6、B7、B8、B8A、B11、RVI、GRVI、KNDVI、LSWI和MSI在3种光谱数据集中都通过了相关显著性检验（P<0.01）. 而DVI、EVI和SATVI仅在中值和均值数据集中与SOC含量存在显著相关性. 分析最小值合成影像、中值合成和均值合成影像可以看出，3种光谱数据集与SOC含量的相关性排序依次为：中值>均值>最小值. 在环境变量中，所有的气候变量都通过了相关显著性检验，而地形变量中仅有ah、cd、cnbl、ele、slope和vd这6个变量与SOC含量存在显著相关性（P<0.01）. 其中，气候变量中的soilevap与SOC含量的相关系数为0.96，高于其他建模变量. 除了径流量和降水以外，其余5个气候变量与SOC含量的相关性均大于0.9，说明气候因素对SOC含量具有重要影响. 在地形变量中，高程（ele）和通道网络基础层（cnbl）的相关性明显高于其他地形变量，且比相关性排名第三的分析式山体阴影（ah）高出了0.56.
|
图 3 SOC含量与建模变量相关关系 Fig. 3 Correlation between SOC content and modeling variables |
在利用相关性分析剔除掉不显著相关（P>0.01）的变量后，为充分发掘模型潜力并提高精度，同时提高建模效率，采用Boruta算法对输入变量进行特征筛选，并与经过相关性分析后未进行筛选的全变量集进行对比. 分别对3种光谱和环境变量集和进行筛选，得到的结果如图4所示. 可以看出对于最小值合成的变量集合，保留了包括全部气候变量和大多数光谱变量以及两个地形变量在内的共17个变量，筛选掉了4个环境变量以及B8和B8A两个光谱变量. 对于均值合成的变量集合，保留了包括全部气候变量、大多数光谱变量以及两个地形变量在内的共18个变量，筛选掉了4个环境变量以及B7、B8、B8A和GRVI这4个光谱变量. 对于中值合成的变量集合，保留了包括全部气候变量、大多数光谱变量以及两个地形变量在内的共18个变量，筛选掉了4个环境变量以及B6、B8、B8A和GRVI这4个光谱变量. 结果表明，气候和地形变量在3个自变量集合中的作用都尤为重要.
|
1.shadowMin，是Boruta算法运行过程中在所有变量基础上生成的阴影变量的最小值，表示所有影子特征重要性的最小值，用于进行变量初次筛选；2.ah，3.shadowMean，是Boruta算法运行过程中在所有变量基础上生成的阴影变量的均值，表示所有影子特征重要性的平均值，用于进行变量第二次筛选；4.cd，5.vd，6.B8A，7.slope，8.GRVI，9.B7，10.B8，11.shadowMax，是Boruta算法运行过程中在所有变量基础上生成的阴影变量的最大值，表示所有影子特征重要性的最大值，用于进行变量最终筛选，若原始变量重要性高于shadowMax则被归类为确定变量；12.B6，13.MSI，14.LSWI，15.B11，16.DVI，17.EVI，18.KNDVI，19.RVI，20.SATVI，21.evap，22.cnbl，23.ele，24.temp，25.precip，26.runoff，27.pressure，28.soiltemp，29.soilevap 图 4 Boruta算法变量筛选结果 Fig. 4 Boruta algorithm variable screening results |
对3个变量集分别构建RF模型，并利用R2、RMSE和RPD来评价模型的效果，不同模型精度如表4所示. 从中可以看出，通过在全光谱变量集合上增加气候和地形等环境变量可以有效提升全光谱变量模型的估测精度，其中R2提高最多的为最小值合成影像模型，相比Set Ⅰ提高了13.83%. 利用Boruta算法筛选后的数据集建模效果较Set Ⅱ数据集在不同程度上都得到了提高. 其中R2提高最多的为最小值合成影像模型变量集合，筛选后的模型相对Set Ⅱ模型提高了1.15%. RMSE提高最多的为中值合成影像模型变量集合，筛选后的模型相对Set Ⅱ模型提高了13.23%. RPD提高最多的是最小值合成影像模型变量集合，Boruta算法筛选后的模型相对Set Ⅱ模型提高了14.63%. 表明Boruta算法在一定程度上能够提高SOC含量反演模型的估测精度. 在所有建模模型中，模型精度最好的为mean-Set Ⅲ模型，R2为0.970，RMSE为2.919 g·kg-1，RPD为5.319，其效果在所有模型中均为最优.
|
|
表 4 估测精度对比1） Table 4 Comparison of estimation accuracy |
结合图5来看，对于随机森林模型可以看出，验证集样本预测值主要集中在15 g·kg-1以下的低值和40 g·kg-1附近的高值区域，与实测SOC含量样本值的大小和分布相一致. 对于3种变量数据集，可以看出全光谱变量集合的预测区间都大于全变量集合与Boruta变量集合，Boruta变量集合的效果要优于全变量集合的效果. 对比3种影像合成方法，可以看出均值合成影像数据集建模的整体效果优于中值合成和最小值合成影像数据集. 对于模型精度最好的为mean-Set Ⅲ模型，其线性拟合趋势线、预测区间以及R2相对于其他模型效果都更好，因此可以选择mean-Set Ⅲ模型进行博斯腾湖流域土壤有机碳含量的空间反演.
|
min表示按最小值合成的影像，median表示按中值合成的影像，mean表示按均值合成的影像；Set Ⅰ表示全光谱变量集合，Set Ⅱ表示包括环境变量和光谱变量在内的全变量集合，Set Ⅲ表示Set Ⅱ经过Boruta算法筛选得到的特征变量集合；min-Set Ⅰ表示最小值合成影像的全光谱变量集合模型，其余组合以此类推；横坐标的实测值表示实际采样得到的SOC含量，纵坐标的预测值表示所建立模型对实测土壤有机碳含量进行预测得到的SOC含量 图 5 RF模型的不同变量组合及其估测结果对比 Fig. 5 Comparison of different variable combinations and prediction results of the RF model |
基于表4和图5可以看出，使用Boruta算法筛选得到的均值合成光谱及环境变量集合所建立的RF模型在所有模型中的估测效果最好. 按照该组合对博斯腾湖流域SOC含量的空间分布进行估测，结果如图6所示. 从中可以看出，博湖流域内部ω（SOC）介于1.907~46.521 g·kg-1之间，且流域内部SOC含量存在较大差异. 在空间分布上，SOC含量较高的区域主要位于流域上游，流域下游的SOC含量相对较低，西部SOC含量显著高于东部，北部SOC含量高于南部. SOC含量最高的区域集中在上游水分充足、温度较为适宜的草地，表明研究区内水分和温度与SOC含量关系密切，植被在土壤固碳方面发挥着积极作用. SOC含量较低的区域主要位于流域东南部的荒地和盐渍地. 整体来看，SOC含量呈现出由高海拔地区到低海拔地区逐渐减少的趋势.
|
基于研究区主要土地利用类型，选择耕地（A）、耕地和荒漠（B）以及草地（C）分别对比分析Set Ⅰ、Set Ⅱ和Set Ⅲ中最佳均值合成影像模型以及min-Set Ⅲ和median-Set Ⅲ模型的SOC含量空间反演差异 图 6 SOC含量空间分布 Fig. 6 Spatial distribution of SOC content |
不同变量组合以及不同光谱合成方式结合RF模型所得到的SOC含量反演结果也各不相同. 从图7中可以看出，A、B和C这3个区域相对于其他组合所得到的结果中，基于光谱指数得到的SOC含量反演结果在范围区域内容易存在一些高值或者低值的异常值，从而影响模型的表现效果，并且对不同地块以及不同地形部位之间的SOC含量差异表现得不太明显. 对比mean-Set Ⅱ 和mean-Set Ⅲ 可以看出虽然前者得到的反演结果能够较好区分耕地与荒地以及不同地形部位之间的差异，但由于不同地块的耕作管理方式和种植作物情况等不同，不同地块SOC含量也存在不同情况的差异. 根据mean-Set Ⅱ 反演结果可以看出，其对不同耕地地块内部的差异区分不是特别明显，对同一地块内部的估测也相对破碎.
|
A~C对应图6的位置 图 7 RF模型的不同变量组合及其空间估测结果对比 Fig. 7 Comparison of different combinations of variables and spatial estimation results in the RF model |
从区域A和区域B可以看出，min-Set Ⅲ 模型反演结果不仅同一地块内部相对破碎，并且不同地块之间的区分也没有其他估测结果相对较准确. 对于median-Set Ⅲ 模型反演结果，虽然其能够准确区分出不同地块及其之间的SOC含量差异，但同一地块内部SOC含量预测结果相对破碎，对小地形部位的SOC含量区分能力相对较弱. 对于本研究中最优模型的mean-Set Ⅲ，可以看出不同地块之间的区分最为明显，同一地块内部的差异也较小，能够较为真实地表现出不同地块与同一地块内部的SOC含量，对于地块与不同地形部位以及不同土地利用类型之间的差异也相对明显. 对于这3种影像合成方式的最优模型来说，草地SOC含量预测结果之间的差异除了预测值的大小以及不同地形部位之间有细小差异外，对于整体SOC含量变化趋势都相对不如耕地明显.
3 讨论 3.1 影像合成方法及光谱变量分析光学遥感影像容易受到例如云、雾和雨等天气因素的影响从而导致影像中存在大量无法使用的像元，影响影像的实际使用效果. 多时相影像合成方法能够弥补单期影像因云层覆盖导致的使用缺陷，满足实际的研究需要. 并且多时相影像不仅可以反映同一地区在不同时间的光谱特征差异，减少环境因素的影响，还能提供更多的光谱信息，而单期影像所包含的光谱信息相对较少且易受到环境因素的干扰[9]. 通过比较本研究中影像光谱反射率发现，月末光谱反射率相比月初最多变化了114.44%，合成影像光谱反射率相比月初和月末最多变化了103.43%. 这表明光谱特征在单月内存在较大变化，合成影像相对单景影像能够包含更丰富信息. 使用多时相遥感卫星影像合成方法有助于提高SOC含量估测的模型精度[10，39]. 对比多种不同的影像合成方式，Gasmi等[40]利用均值影像合成方法创建的裸土图像对土壤黏土含量预测最为准确. Luo等[41]通过研究发现中值合成影像对于土壤有机碳含量的建模效果要优于其他影像合成方式. 在本研究中，中值合成的全光谱变量集合和全变量集合建模效果以及经过变量筛选后的变量集合相对于均值合成影像集合的建模效果较差. 因此，对于不同多时相影像合成方法的建模效果仍需进行更多的研究对比.
光学遥感影像无法直接观测到植被覆盖地区下的土壤属性和状况，但通过不同波段之间进行比值或者归一化计算得到的光谱指数可以较好地对SOC和进行响应[42]. 根据表4可以看出利用多个遥感光谱变量结合随机森林模型也能够对SOC含量进行较好地估测. 在全光谱变量集中，B6、B7、B8、B8A以及B11这5个波段在3个光谱变量集中都通过了相关性检验，主要是由于红边波段、近红外波段对SOC具有独特的响应特征[43]. 红波和绿波段在本研究中没有通过相关性检验，但作为SOC吸收和反射的特征波段，通过对红波和绿波段与其他波段进行比值或归一化运算得到的光谱指数，例如GRVI、RVI、DVI、KNDVI以及SATVI等植被指数通过了相关性检验并在模型估测反演中发挥着重要作用. 使用全光谱变量进行建模的结果表明：对于非裸土像素，利用遥感影像波段和合适的光谱指数，可以较好地建立SOC含量的反演模型.
3.2 环境变量的重要性分析不同建模变量对土壤有机碳含量估测的效果各不相同，通过在光谱变量的基础上增加气象因子和地形等环境变量，提高了SOC含量的估测精度. 结合图4的Boruta变量筛选结果，可以看出在所有的环境变量中，气候变量尤其是土壤蒸散量、土壤表面温度、地表压力、近地面气温、降水以及径流量这6个变量在SOC含量建模中的重要性都大于8，显著高于其他变量的重要性，这表明水热状况在本研究中对土壤有机碳含量具有重要的影响. 温度与土壤有机碳的分解速度呈负相关，较高的温度会使土壤有机碳分解速度加快，而较低的温度会导致有机碳分解减少，从而增加有机碳积累[44]，因此土壤表面温度、近地面气压以及近地面气温与SOC含量都呈显著的负相关. 在干旱半干旱地区，土壤和空气湿度的增加以及降水能够给微生物和植物提供良好的生存发展条件，增强土壤团聚体的稳定性，为土壤有机碳含量的积累创造良好的条件[45]. 因此裸土蒸散量和地表总蒸散量与SOC含量呈显著负相关，降水量和径流总和与SOC含量呈正相关. 地形等环境变量主要是通过对降水和温度等水热状况进行重新分配进而影响SOC含量的空间差异[46]. 对于地形差异较大的区域，高程是影响SOC空间分布及含量大小差异的主要因素，影响其他地形和气象因子的变化，对土壤有机碳的积累有显著影响[47]. 本研究中高程与SOC含量相关系数为0.91，具有显著的正相关性（P<0.01）. 在本研究中，环境变量在SOC含量的估测中起着重要作用. 然而由于环境变量空间分辨率较低，仅仅依赖环境变量进行建模反演无法全面捕捉土壤有机碳含量的复杂空间变化. 特别是对于表现土壤有机碳含量空间分布及局部细节变化的效果较差. 将环境变量结合较高空间分辨率的光谱变量进行SOC含量空间反演，可以在充分发挥环境变量作用的基础上更好地优化模型估测效果，表现空间细节.
3.3 样本数量对SOC含量反演效果的影响土壤数据库通过综合利用不同地区和国家的土壤数据，提供了大范围区域甚至全球范围内的土壤数据描述，但由于数据集样点数量以及空间差异的限制，数据集只能对部分地区的土壤属性进行概括性描述，具有一定的质量限制[48]. 在大范围区域有限样点的情况下，如果不填充区域空间上的数据空白，而采用机器学习算法进行外推反演会产生较线性模型更差的估测效果[49]. 在本研究中，通过对拟合样点插补前后最优模型的SOC含量反演结果对比发现，利用实测样点得到的反演结果在数据空白区域中出现了明显的条带和界限，而插补后的反演结果更能凸显不同地形和土地利用类型之间SOC含量的差异. 因此在研究中通过插补样点来填补部分区域的数据空白有利于反映流域真实的SOC含量状况. 且通过单独提取并计算模型中实测样点的R2、RMSE和RPD发现，估测精度会略微降低，但仍然保持和本研究中相同的结论，即利用Boruta算法选择合适的光谱变量和环境变量能够有效提高模型的估测精度，加入环境变量后模型的估测效果明显要优于单独使用光谱变量的模拟效果. 因此，足量的样本数据和合适的建模变量组合能够有效提高SOC含量估测结果的准确性.
4 结论（1）引入地形变量能够明显提高土壤有机碳含量的估测精度，气候因素对博斯腾湖流域土壤有机碳含量的估测具有重要作用.
（2）通过对比全变量集和筛选变量集的建模效果，发现利用Boruta算法筛选得到的变量集合能够提高随机森林模型的估测精度.
（3）不同影像合成变量集的估测结果表明，利用均值合成的光谱变量结合环境变量的建模效果要优于最小值以及中值合成影像的估测效果.
| [1] |
李健明, 康雨欣, 蒋福祯, 等. 基于Meta分析的煤矿区植被恢复对土壤有机碳储量的影响[J]. 环境科学, 2024, 45(3): 1629-1643. Li J M, Kang Y X, Jiang F Z, et al. Effect of vegetation restoration on soil organic carbon storage in coal mining areas based on meta-analysis[J]. Environmental Science, 2024, 45(3): 1629-1643. |
| [2] |
张蕊, 赵学勇, 王瑞雄, 等. 干旱半干旱区土地利用/覆被类型变化对土壤有机碳动态影响的研究进展[J]. 生态学杂志, 2024, 43(7): 2222-2230. Zhang R, Zhao X Y, Wang R X, et al. Review of the effects of land use/cover change on soil organic carbon dynamic in arid and semi-arid areas[J]. Chinese Journal of Ecology, 2024, 43(7): 2222-2230. |
| [3] | Penman D E, Rugenstein J K C, Ibarra D E, et al. Silicate weathering as a feedback and forcing in Earth's climate and carbon cycle[J]. Earth-Science Reviews, 2020, 209. DOI:10.1016/j.earscirev.2020.103298 |
| [4] |
邓旭哲, 韩晨, 薛利祥, 等. 增温施肥对稻麦农田土壤有机碳及其活性组分的影响[J]. 环境科学, 2023, 44(3): 1553-1561. Deng X Z, Han C, Xue L X, et al. Effects of warming and fertilization on soil organic carbon and its labile components in rice-wheat rotation[J]. Environmental Science, 2023, 44(3): 1553-1561. |
| [5] |
任必武, 陈瀚阅, 张黎明, 等. 机器学习用于耕地土壤有机碳空间预测对比研究——以亚热带复杂地貌区为例[J]. 中国生态农业学报(中英文), 2021, 29(6): 1042-1050. Ren B W, Chen H Y, Zhang L M, et al. Comparison of machine learning for predicting and mapping soil organic carbon in cultivated land in a subtropical complex geomorphic region[J]. Chinese Journal of Eco-Agriculture, 2021, 29(6): 1042-1050. |
| [6] | Dharumarajan S, Kalaiselvi B, Suputhra A, et al. Digital soil mapping of soil organic carbon stocks in western Ghats, South India[J]. Geoderma Regional, 2021, 25. DOI:10.1016/j.geodrs.2021.e00387 |
| [7] | Dou X, Wang X, Liu H J, et al. Prediction of soil organic matter using multi-temporal satellite images in the Songnen Plain, China[J]. Geoderma, 2019, 356. DOI:10.1016/j.geoderma.2019.113896 |
| [8] | Yang L, He X L, Shen F X, et al. Improving prediction of soil organic carbon content in croplands using phenological parameters extracted from NDVI time series data[J]. Soil and Tillage Research, 2020, 196. DOI:10.1016/j.still.2019.104465 |
| [9] | Wang X, Wang L P, Li S J, et al. Remote estimates of soil organic carbon using multi-temporal synthetic images and the probability hybrid model[J]. Geoderma, 2022, 425. DOI:10.1016/j.geoderma.2022.116066 |
| [10] | Shi P, Six J, Sila A, et al. Towards spatially continuous mapping of soil organic carbon in croplands using multitemporal Sentinel-2 remote sensing[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2022, 193: 187-199. DOI:10.1016/j.isprsjprs.2022.09.013 |
| [11] | Pastur G M, Acuña M C A, Silveira E M O, et al. Mapping soil organic carbon content in Patagonian forests based on climate, topography and vegetation metrics from satellite imagery[J]. Remote Sensing, 2022, 14(22). DOI:10.3390/rs14225702 |
| [12] |
祝元丽, 冯向阳, 闫庆武, 等. 基于GBDT的望奎县农田土壤有机碳主控因子研究[J]. 中国环境科学, 2024, 44(3): 1407-1417. Zhu Y L, Feng X Y, Yan Q W, et al. Spatial distribution and main controlling factors of soil organic carbon under cultivated land based on GBDT model in black soil region of Northeast China[J]. China Environmental Science, 2024, 44(3): 1407-1417. |
| [13] |
闫坤, 杨慧敏, 王德彩. 基于NDVI时间序列特征的森林土壤有机碳数字制图——以济源南山林场为例[J]. 土壤通报, 2024, 55(4): 921-931. Yan K, Yang H M, Wang D C. Digital mapping of forest soil organic carbon based on NDVI time series features: a case study of Jiyuan Nanshan forest farm[J]. Chinese Journal of Soil Science, 2024, 55(4): 921-931. |
| [14] | Zhou Q, Ding J L, Ge X Y, et al. Estimation of soil organic matter in the Ogan-Kuqa River Oasis, Northwest China, based on visible and near-infrared spectroscopy and machine learning[J]. Journal of Arid Land, 2023, 15(2): 191-204. DOI:10.1007/s40333-023-0094-4 |
| [15] | Keskin H, Grunwald S, Harris W G. Digital mapping of soil carbon fractions with machine learning[J]. Geoderma, 2019, 339: 40-58. DOI:10.1016/j.geoderma.2018.12.037 |
| [16] |
张发, 玉素甫江·如素力. 基于LUCC追踪分析的生态系统服务价值时空变化研究——以博斯腾湖流域为例[J]. 北京林业大学学报, 2021, 43(7): 88-99. Zhang F, Rusuli Y. Spatio-temporal variation of ecosystem service value based on LUCC trajectories: a case study of Bosten Lake Watershed[J]. Journal of Beijing Forestry University, 2021, 43(7): 88-99. |
| [17] | Yang X, Ji G X, Wang C, et al. Modeling nitrogen and phosphorus export with InVEST model in Bosten Lake Basin of Northwest China[J]. PLoS One, 2019, 14(7). DOI:10.1371/journal.pone.0220299 |
| [18] | Castaldi F, Koparan M H, Wetterlind J, et al. Assessing the capability of Sentinel-2 time-series to estimate soil organic carbon and clay content at local scale in croplands[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2023, 199: 40-60. DOI:10.1016/j.isprsjprs.2023.03.016 |
| [19] |
白燕英, 高聚林, 张宝林. 基于NDVI与EVI的作物长势监测研究[J]. 农业机械学报, 2019, 50(9): 153-161. Bai Y Y, Gao J L, Zhang B L. Monitoring of crops growth based on NDVI and EVI[J]. Transactions of the Chinese Society for Agricultural Machinery, 2019, 50(9): 153-161. |
| [20] | Xu S Z, Xu X G, Blacker C, et al. Estimation of leaf nitrogen content in rice using vegetation indices and feature variable optimization with information fusion of multiple-sensor images from UAV[J]. Remote Sensing, 2023, 15(3). DOI:10.3390/rs15030854 |
| [21] |
李顿, 王雪梅, 李坤玉, 等. 基于变量筛选与机器学习算法的渭-库绿洲土壤有机质含量估测研究[J]. 地球与环境, 2024, 52(3): 375-385. Li D, Wang X M, Li K Y, et al. Estimation of soil organic matter content in Wei-Ku oasis based on variables screening and machine learning algorithms[J]. Earth and Environment, 2024, 52(3): 375-385. |
| [22] | Yang Y Y, Shang K, Xiao C C, et al. Spectral index for mapping topsoil organic matter content based on ZY1-02D satellite hyperspectral data in Jiangsu Province, China[J]. ISPRS International Journal of Geo-Information, 2022, 11(2). DOI:10.3390/ijgi11020111 |
| [23] | Ma T Y, Hu Y, Wang J, et al. A novel vegetation index approach using sentinel-2 data and random forest algorithm for estimating forest stock volume in the Helan Mountains, Ningxia, China[J]. Remote Sensing, 2023, 15(7). DOI:10.3390/rs15071853 |
| [24] | Camps-Valls G, Campos-Taberner M, Moreno-Martínez Á, et al. A unified vegetation index for quantifying the terrestrial biosphere[J]. Science Advances, 2024, 7(9). DOI:10.1126/sciadv.abc7447 |
| [25] | Gu H B, Mills C, Ritchie G L, et al. Water stress assessment of cotton cultivars using unmanned aerial system images[J]. Remote Sensing, 2024, 16(14). DOI:10.3390/rs16142609 |
| [26] |
郑舒元, 海燕, 何孟琦, 等. 高区分度的高分六号影像可见光波段植被指数构建[J]. 光谱学与光谱分析, 2023, 43(11): 3509-3517. Zheng S Y, Hai Y, He M Q, et al. Construction of vegetation index in visible light band of GF-6 image with higher discrimination[J]. Spectroscopy and Spectral Analysis, 2023, 43(11): 3509-3517. |
| [27] | Lin C, Zhu A X, Wang Z F, et al. The refined spatiotemporal representation of soil organic matter based on remote images fusion of Sentinel-2 and Sentinel-3[J]. International Journal of Applied Earth Observation and Geoinformation, 2020, 89. DOI:10.1016/j.jag.2020.102094 |
| [28] |
傅楷翔, 贾国栋, 余新晓, 等. 基于改进遥感生态指数的吐鲁番-哈密地区生态环境质量评价及驱动机制分析[J]. 生态学报, 2024, 44(9): 3911-3923. Fu K X, Jia G D, Yu X X, et al. Evaluation of ecological environment quality and analysis of drivingmechanism in Tulufan-Hami region based on improved remote sensing ecological indices[J]. Acta Ecologica Sinica, 2024, 44(9): 3911-3923. |
| [29] | Liu Y, Qian J X, Yue H. Combined Sentinel-1A with Sentinel-2A to estimate soil moisture in farmland[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2021, 14: 1292-1310. DOI:10.1109/JSTARS.2020.3043628 |
| [30] |
薛朝辉, 钱思羽. 融合Landsat 8与Sentinel-2数据的红树林物候信息提取与分类[J]. 遥感学报, 2022, 26(6): 1121-1142. Xue Z H, Qian S Y. Fusion of Landsat 8 and Sentinel-2 data for mangrove phenology information extraction and classification[J]. National Remote Sensing Bulletin, 2022, 26(6): 1121-1142. |
| [31] | Sleutel S, De Neve S, Hofman G. Assessing causes of recent organic carbon losses from cropland soils by means of regional-scaled input balances for the case of Flanders (Belgium)[J]. Nutrient Cycling in Agroecosystems, 2007, 78(3): 265-278. DOI:10.1007/s10705-007-9090-x |
| [32] | Guth P L, Geoffroy T M. LiDAR point cloud and ICESat‐2 evaluation of 1 second global digital elevation models: Copernicus wins[J]. Transactions in GIS, 2021, 25(5): 2245-2261. DOI:10.1111/tgis.12825 |
| [33] |
毛继华, 赵恒谦, 金倩, 等. 河北铅锌尾矿库区土壤重金属含量高光谱反演方法对比[J]. 农业工程学报, 2023, 39(22): 144-156. Mao J H, Zhao H Q, Jin Q, et al. Comparative study on the hyperspectral inversion methods for soil heavy metal contents in Hebei lead-zinc tailings reservoir areas[J]. Transactions of the Chinese Society of Agricultural Engineering, 2023, 39(22): 144-156. DOI:10.11975/j.issn.1002-6819.202307092 |
| [34] | Tian H, Zhang L N, Li M, et al. Weighted SPXY method for calibration set selection for composition analysis based on near-infrared spectroscopy[J]. Infrared Physics & Technology, 2018, 95: 88-92. |
| [35] |
孟鑫鑫, 于雷, 周勇, 等. 基于可见近红外和中红外近地面光谱数据融合的土壤有机碳含量反演[J]. 土壤通报, 2022, 53(2): 301-307. Meng X X, Yu L, Zhou Y, et al. Predicting organic carbon using data fusion of visible near-infrared and middle infrared spectra by proximal soil sensing[J]. Chinese Journal of Soil Science, 2022, 53(2): 301-307. |
| [36] | Gomes L C, Faria R M, de Souza E, et al. Modelling and mapping soil organic carbon stocks in Brazil[J]. Geoderma, 2019, 340: 337-350. DOI:10.1016/j.geoderma.2019.01.007 |
| [37] | Azizi K, Ayoubi S, Demattê J A M. Controlling factors in the variability of soil magnetic measures by machine learning and variable importance analysis[J]. Journal of Applied Geophysics, 2023, 210. DOI:10.1016/j.jappgeo.2023.104944 |
| [38] |
安柏耸, 王雪梅, 黄晓宇, 等. 基于连续小波变换的土壤重金属镉含量的高光谱估测[J]. 地球与环境, 2023, 51(2): 246-253. An B S, Wang X M, Huang X Y, et al. Hyperspectral estimation of heavy metal cadmium content in soil based on continuous wavelet transform[J]. Earth and Environment, 2023, 51(2): 246-253. |
| [39] | Luo C, Wang Y A, Zhang X L, et al. Spatial prediction of soil organic matter content using multiyear synthetic images and partitioning algorithms[J]. CATENA, 2022, 211. DOI:10.1016/j.catena.2022.106023 |
| [40] | Gasmi A, Gomez C, Lagacherie P, et al. Mean spectral reflectance from bare soil pixels along a Landsat-TM time series to increase both the prediction accuracy of soil clay content and mapping coverage[J]. Geoderma, 2021, 388. DOI:10.1016/j.geoderma.2020.114864 |
| [41] | Luo C, Zhang X L, Meng X T, et al. Regional mapping of soil organic matter content using multitemporal synthetic Landsat 8 images in Google earth engine[J]. CATENA, 2022, 209. DOI:10.1016/j.catena.2021.105842 |
| [42] | Zhang Z P, Ding J L, Wang J Z, et al. Prediction of soil organic matter in Northwestern China using fractional-order derivative spectroscopy and modified normalized difference indices[J]. CATENA, 2020, 185. DOI:10.1016/j.catena.2019.104257 |
| [43] | Xie B Q, Ding J L, Ge X Y, et al. Estimation of soil organic carbon content in the Ebinur Lake wetland, Xinjiang, China, based on multisource remote sensing data and ensemble learning algorithms[J]. Sensors, 2022, 22(7). DOI:10.3390/s22072685 |
| [44] | Georgiou K, Koven C D, Wieder W R, et al. Emergent temperature sensitivity of soil organic carbon driven by mineral associations[J]. Nature Geoscience, 2024, 17(3): 205-212. DOI:10.1038/s41561-024-01384-7 |
| [45] | Denef K, Stewart C E, Brenner J, et al. Does long-term center-pivot irrigation increase soil carbon stocks in semi-arid agro-ecosystems?[J]. Geoderma, 2008, 145(1-2): 121-129. DOI:10.1016/j.geoderma.2008.03.002 |
| [46] | Zhu M, Feng Q, Qin Y Y, et al. The role of topography in shaping the spatial patterns of soil organic carbon[J]. CATENA, 2019, 176: 296-305. DOI:10.1016/j.catena.2019.01.029 |
| [47] |
霍瑢彦町, 刘京, 张锐, 等. 黄土丘陵区土壤有机碳含量影响因素分析及空间估测[J]. 环境科学, 2025, 46(4): 2301-2312. Huo R Y T, Liu J, Zhang R, et al. Influencing factors analysis and spatial estimation of soil organic carbon content in hilly areas of Loess Plateau[J]. Environmental Science, 2025, 46(4): 2301-2312. |
| [48] | Guevara M, Arroyo C, Brunsell N, et al. Soil organic carbon across Mexico and the conterminous United States (1991-2010)[J]. Global Biogeochemical Cycles, 2020, 34(3). DOI:10.1029/2019GB006219 |
| [49] | Hengl T, Mendes de Jesus J, Heuvelink G B M, et al. SoilGrids250m: global gridded soil information based on machine learning[J]. PLoS One, 2017, 12(2). DOI:10.1371/journal.pone.0169748 |