基于连续投影算法的博斯腾湖西岸湖滨绿洲土壤有机碳含量的高光谱估算

牛芳鹏，李新国*，麦麦提吐尔逊·艾则孜，赵慧

（1.新疆师范大学地理科学与旅游学院，乌鲁木齐 830054；2.新疆干旱区湖泊环境与资源实验室，乌鲁木齐 830054）

土壤是碳在陆地生态系统中最大的储存库，并影响其生产力的发展，其中土壤有机碳（soil organic carbon,SOC）是构成土壤碳库的重要部分[1]。干旱半干旱区是稳定全球生态环境变化的重要环节之一，绿洲及其生态系统在对抗与缓解干旱区气候方面起着重要作用，对土壤碳的研究特别是绿洲土壤有机碳的研究对全球碳循环及生态环境保护具有重要意义[2]。土壤光谱中包含丰富的土壤信息，便携快速、无损与高精度的波谱特性使高光谱遥感技术在估算SOC含量的应用中，从定性到定量越来越成熟[3]。筛选具有代表性的光谱响应波段是构建高光谱估算模型的基础，学者们在进行光谱数据特征波段的优选方面已做了大量研究[4-5]。WANG 等研究表明，利用连续投影算法（successive projection algorithm,SPA）选择的光谱特征波段可以有效地提高偏最小二乘回归（partial least square regression,PLSR）模型的决定系数（coefficient of determination,R2）和残余预测误差（residual prediction deviation,RPD）性能[6]；王涛等采用相关性分析法和SPA 筛选特征波长，构建的最优模型精度R2=0.98，均方根误差（root mean square error,RMSE）=0.30 g/kg，RPD=9.32[7]；章海亮等采用遗传算法结合SPA 挑选出18个特征波段并构建PLSR模型，其预测集的R2=0.83，RMSE=0.20 g/kg，RPD=2.45[8]。然而，利用土壤光谱反射率对SOC含量的估算，大多是通过建立线性方程来进行，这些传统的回归方法受数据本身影响较大；而支持向量机（support vector machine,SVM）的发展，有效解决了样本维数问题[9]。张子鹏等通过构建SVM 模型，比较了不同维度光谱参数对SOC含量估测的准确性[10]；SUN 等对光谱数据进行SG（Savitzky-Golay）平滑与多种数据变换后，利用PLSR 方法构建了复垦土地SOC 含量估算模型，其验证集R2=0.78，RMSE=1.81 g/kg，RPD=2.69[11]；HONG 等利用SVM 算法对华中地区土壤有机质（soil organic matter, SOM）含量进行估算，结果表明，用1.5 阶微分构建的SVM 模型精度较高，R2达0.88，RPD=2.89[12]。

本研究以博斯腾湖西岸湖滨绿洲为研究区，通过测量采集的255 个样本的SOC 含量与高光谱数据，并对光谱反射率进行SG平滑、标准化正态变换（standard normal variate, SNV）与一阶微分（first derivative,1st Der）预处理，运用SPA 与相关系数法（correlation coefficient,CC）从原始光谱和预处理后的光谱中筛选特征波段，对比用PLSR 与SVM 2 种方法构建的SOC含量的高光谱估算模型的差异，为利用可见- 近红外光谱（visible-near infrared spectroscopy,Vis-NIR）技术快速无损地估算湖滨绿洲SOC含量提供方法支持。

1 研究区概况与研究方法

1.1 研究区概况

博斯腾湖西岸湖滨绿洲位于新疆焉耆盆地，行政区隶属于新疆博湖县，地理位置为41°45′—42°10′N，86°15′—86°55′E，海拔1 047～1 048 m，地势由西北向东南倾斜，为典型的山前湖泊绿洲。夏季月平均气温22.8 ℃，干旱少雨，蒸发强烈；冬季月平均气温9.0 ℃，寒冷干燥，无霜期176～200 d，多年平均降水量83.55 mm，年均温8.0～8.6 ℃，季节过渡快，属于大陆性荒漠气候区；主要的自然植被类型有胡杨、柽柳和梭梭等；主要的土壤类型有绿洲潮土、草甸土、荒漠林土和沼泽土等[13]。

1.2 土样采集与制备

综合考虑到研究区范围内的主要土地利用类型、植被覆盖类型和微地形等因素，土样采集按照“S”形线路随机均匀布点，采样时间为2019 年10月，共布设51 个典型土壤剖面，并进一步分割出5个不同深度（0～10、＞10～20、＞20～30、＞30～40和＞40～50 cm）土层进行取样，去除杂物后采用四分法选取200 g土装袋，共计255个样品。带回实验室自然风干后进行研磨和过60目孔筛，一部分用作土壤高光谱数据的测试；另一部分用重铬酸钾-外加热法测定SOC含量[14]。

1.3 土壤光谱的测定及预处理

光谱数据通过ASD FieldSpec3 地物光谱仪于室外采集，光谱波长为350～2 500 nm。选择云量小于5%、无风或风力低于3 级的晴朗天气，采集时间为12：00—14：00；光谱采集前先预热仪器30 min并进行白板校正；光纤探头视场角不超过25°，垂直于土样表面15 cm处，使用五点梅花采样法采集，在每个位置测量3条光谱曲线，共15条光谱曲线记录存档，每测完一组剖面土样采集一次暗电流，同时进行白板优化校正，减小误差[15]。由于环境不可控因素，统一去除噪声较大的尾部波段（2 451～2 500 nm）以及受环境水汽影响的1 300～1 450和1 800～1 950 nm 波段；此外，为减少高频噪声的影响，提升光谱数据信噪比，对原始光谱曲线进行窗口数为5的二次多项式SG平滑处理，并联合使用SNV与1st Der 进行光谱数据预处理[16]。图1 为预处理后的土壤光谱反射率曲线图。

图1 SG 平滑及SG-SNV-1st Der 预处理后的土壤光谱反射率曲线Fig.1 Soil spectral reflectance curve after SG smoothing and SG-SNV-1st Der pretreatments

1.4 连续投影算法的特征变量选择

连续投影算法（SPA）是一种前向变量选择算法，运用向量投影分析选择最大向量，最终通过校正模型提取几个特征波长[17]。其优点是从光谱矩阵中选择最小共线性的变量组合，从而降低模型的冗余度，提高模型的稳定性和准确性。SPA 的具体步骤如下：

记xk(0)和N分别为初始迭代向量与需要提取的变量个数，光谱矩阵为J列。

1）任选光谱矩阵的1列（第j列），把建模集的第j列赋值给xj，记为xk(0)。2）将未选入的列向量位置的集合记为s，

3）用xj分别对剩余的列向量投影进行计算，

4）选取向量投影最大的光谱波长，

5）令xj=Px，j∈s。

6）n=n+1，如果n＜N，则按公式（1）循环计算。

最后，提取出的变量为{xk(n)=0，…，N-1}，分别构建关于每一次循环中k(0)和N的多元线性回归（multiple linear regression, MLR）模型，得到建模集交互验证的RMSE，以及对应的不同候选子集，其中最小的RMSE值对应的k(0)和N就是最优值。

1.5 集合划分与模型构建

采用基于联合X-Y距离（sample set partitioning based on jointX-Ydistances, SPXY）方法进行建模样本的划分，该方法充分考虑了X和Y的空间可变性，计算样本间的距离时能同时考虑SOC含量的影响与土壤光谱信息特征[18]。SPA、PLSR模型与SVM模型构建运用MatLab R2020a 软件完成。PLSR 模型中采用“留一法”交叉验证来决定最佳主成分数；SVM 模型选择Poly 核函数进行分析。模型精度选用R2、RMSE、RPD 和四分位数间距性能比（ratio of performance to interquartile distance, RPIQ）进行评价，其中：R2的取值范围为0～1.00，R2越大，RMSE越接近0，模型拟合度越高。当RPD＜1.40时，模型预测能力较差；当1.40≤RPD＜2.00 时，模型预测效果一般；当RPD≥2.00 时，模型预测能力较好[19]。RPIQ为样本观测值第三、四分位数Q3和第一、四分位数Q1的差与RMSE的比值，RPIQ越大，表示模拟结果分布的偏态程度越小[20]。

2 结果与分析

2.1 土壤有机碳含量描述

对255个样本光谱反射率数据进行主成分分析后，剔除11 个异常样本，运用SPXY 方法对剩余样本进行训练集与验证集的划分。由表1 可知：SOC质量分数变化范围为0.75～48.13 g/kg，均值为13.31 g/kg，呈中等变异性，变异系数为63.19%；训练集与验证集的SOC质量分数平均值分别为12.67和12.16 g/kg，分别为右偏平顶峰与右偏尖顶峰，呈中等变异性，变异系数分别为52.33%和58.56%。

表1 土壤有机碳含量描述性统计Table 1 Descriptive statistics of soil organic carbon(SOC)content

2.2 土壤有机碳含量光谱曲线特征

根据研究区SOC含量状况，并结合前人研究成果[21]，采用K-均值（K-means）聚类分析方法将SOC质量分数划分为≤0.75、＞0.75～17.40、＞17.40～32.72、＞32.72 g/kg 4 类，图2 为4 种不同SOC 含量的平均光谱曲线进行基线校正后的效果图。从中可见，4 类光谱曲线形状变化基本一致。在350～1 800 nm之间同一波段的不同SOC含量光谱反射率差异较明显，SOC含量越高，土壤光谱反射率越低；在1 950～2 450 nm之间，光谱反射率差异较小。当SOC质量分数≤0.75 g/kg时，其反射率均值为0.38；当SOC质量分数＞0.75～17.40 g/kg时，其反射率均值为0.37；当SOC 质量分数＞17.40～32.72 g/kg 时，其反射率均值为0.35；当SOC质量分数＞32.72 g/kg时，其反射率均值为0.29。在350～569 nm 之间光谱反射率增幅较大，在570～1 299 nm 之间光谱反射率变化趋势趋于平缓，在1 451～2 450 nm之间光谱反射率的波动较大；在923～955、1 109～1 141、2 148～2 240与2 333～2 358 nm波段内存在4个较明显的吸收谷。

图2 土壤有机碳含量的光谱反射率曲线特征Fig.2 Characteristics of spectral reflectance curve of soil organic carbon content

2.3 土壤有机碳含量与光谱曲线的相关性分析

由图3 可知：SOC 含量与原始光谱反射率表现为负相关性，－0.62＜相关系数（r）＜－0.07，有1 160 个波段通过极显著性检验（P＜0.01），主要集中在524～1 299、1 469～1 790 与1 973～2 056 nm之间，在661 nm 波长处相关性最高，相关系数∣r∣为0.62。经SG-SNV-1st Der预处理后，光谱反射率曲线呈正负波动，放大了原始光谱曲线的细微变化；通过极显著性检验（P＜0.01）的波段数量缩减到414个，相关性较高波段主要集中在487～575、725～998 和1 464～1 514 nm 范围内，在788、800 与1 768 nm波长处相关性最高，r均大于0.80。

图3 土壤有机碳含量与光谱反射率的相关系数曲线Fig.3 Correlation coefficient curve between soil organic carbon content and spectral reflectance

2.4 基于SPA 对特征变量的筛选结果

由图4 可知：利用SPA 对预处理后的光谱数据进行特征波长的筛选时，随着筛选变量数量的增加，RMSE 先是迅速下降，当变量数为14 时，RMSE趋于稳定状态，其值为5.010 1 g/kg。通过SPA运算后得到14 个特征波长，仅占全光谱数据的0.78%，分别为399、1 011、1 046、1 061、1 073、1 596、1 632、1 667、1 749、2 012、2 103、2 268、2 305、2 341 nm，极大地缩减了光谱信息中的冗余变量。同理，运用SPA对原始光谱数据进行筛选，最终得到19个最优特征变量，占全波段光谱数据的1.06%。

图4 利用SPA筛选预处理光谱特征波长Fig.4 Screening of preprocessing spectral characteristic wavelength by SPA

2.5 基于PLSR 的土壤有机碳含量估算模型

由表2 可知，光谱预处理后，利用PLSR 模型获得的RPD 均大于1.40，可以较好地估算SOC 含量。原始光谱通过SPA 构建的PLSR 模型验证集R2为0.75，RMSE 为3.98 g/kg，获得的RPD 为1.79，RPIQ为2.01，模型精度高于全波段（full-band）建模（R2=0.65）；通过相关系数法筛选出相关系数最高的5个波段进行建模，验证集R2为0.70，RMSE 为6.14 g/kg，RPD 为1.16，RPIQ 为1.30，建模效果略好于全波段。光谱预处理后SPA模型精度最高，训练集R2为0.79，RMSE为5.73 g/kg，验证集R2为0.79，RMSE为3.58 g/kg，RPD为1.99，RPIQ为2.23。

表2 2种光谱模式的PLSR建模结果Table 2 PLSR modeling results of the two spectral modes

综合比较2 种光谱模式下3 种变量的PLSR 模型精度，其模型估算能力表现为SPA＞相关系数法（CC）＞全波段法。由图5可知：光谱预处理后构建的PLSR 模型较原始光谱PLSR 模型数据点更靠近1∶1 线；通过SPA 构建的PLSR 模型样本的SOC 含量预测值小于实测值，大多较均匀地分布在1∶1 线下方。

图5 2种光谱模式下PLSR估算模型结果图Fig.5 Scatter plot of PLSR model of full-band and characteristic bands under the two spectral modes

2.6 基于SVM 的土壤有机碳含量估算模型

由表3 可知，基于全波段光谱构建的2 种SVM模型，其验证集R2分别为0.68和0.76，RMSE分别为4.21 和4.26 g/kg，RPD 分别为1.69 和1.67，RPIQ 分别为1.90 和1.88，可粗略对样本进行预测。通过相关系数法（CC）与SPA 进行特征变量筛选后构建的SVM 模型的RPD 均大于2.00，拟合程度较好，模型效果有明显提高。原始光谱基于相关系数法构建的SVM模型训练集和验证集R2分别为0.69和0.70，RPD 为2.00，RPIQ 为2.25；光谱预处理后基于相关系数法构建的SVM 模型训练集和验证集R2分别为0.80 和0.77，RPD 为2.13，RPIQ 为2.39。原始光谱与预处理光谱基于SPA 构建的模型验证集R2分别为0.73 和0.81，RMSE 分别为3.35 和3.16 g/kg，RPD分别为2.13和2.25，RPIQ分别为2.39和2.53。综合比较3种不同变量构建的SVM模型效果，对SOC含量估算的效果依次为SPA＞相关系数法＞全波段法。

表3 2种光谱模式的SVM建模结果Table 3 SVM modeling results of the two spectral modes

由图6 可知：光谱预处理后基于SPA 与相关系数法构建的SOC 含量的SVM 模型估算结果较好，样本均匀地接近1∶1线；全波段SVM模型及原始光谱SPA 模型验证样本逐渐偏离1∶1 线，验证样本数据点分布比较分散，模型预测效果较差，与前文分析一致。

图6 2种光谱模式下SVM估算模型结果图Fig.6 Scatter plot of SVM model of full-band and characteristic bands under the two spectral modes

3 讨论

关于室内与室外土壤光谱反射率的采集方法始终是土壤研究的主要课题之一。室内测试是为了研究土壤中的某个因子对光谱反射率特性的影响，室外测试能较好地反映自然景观的真实性，描述表面反射特征，以便为航空和航天传感器定标[22-23]。徐彬彬等认为，在土壤光谱反射特性研究中，应当注重野外实测，尽管在野外测试中受到当前所用仪器的限制和环境条件的影响，但它还是能较好地反映自然界的部分真实情况[24]；马利芳等于野外采集高光谱数据后，研究了新疆阜康市土壤盐分离子的高光谱特征，为区域尺度的土壤盐分主要离子含量估算提供了良好的支撑[25]。本研究采用室外光谱采集的方法，使获得的数据能更好地接近研究区自然环境条件下土壤有机碳高光谱信息。

原始光谱经SG-SNV-1st Der 预处理后，通过全波段、相关系数法与SPA 构建的PLSR 模型验证集R2分别提高了6.15%、5.71%和5.33%；SVM 模型验证集R2分别提高了11.76%、10.00%和10.96%。这与张子鹏等[16]和李冠稳等[26]的研究结果基本一致。在光谱建模之前对光谱数据进行预处理，既能突出光谱的特征波段，还可以提高模型的拟合效果。本研究运用相关系数法与SPA分别筛选出5和19个特征波段，仅占全波段数据的0.28%和1.06%，验证集R2高于全波段建模1.31%和6.58%。这与韩建等[27]和VISCARRA ROSSEL等[28]的研究结果相一致。SPA 可以有效地消除波段之间的共线性影响，剔除不相关变量，降低数据冗余度，提高模型的建模精度。光谱预处理后，基于3 种变量方法构建的SVM 模型较PLSR 模型的验证集R2分别提高了10.14%、4.05%和2.53%，说明SVM模型能在一定程度上弥补PLSR模型在解决非线性关系问题上的缺陷。这与杨爱霞等[29]、曾胤等[30]的研究结果相一致。本文构建的SOC 含量高光谱估算模型尚未考虑土地利用类型、土壤结构和土壤水分等因素，其对模型的影响还需进一步验证。

4 结论

1）研究区SOC 质量分数变化范围为0.75～48.13 g/kg，平均值为13.31 g/kg，呈中等变异性，变异系数为63.19%；同一波段内，随着土层深度的增加，光谱反射率越高，且SOC含量越高，土壤光谱反射率越低。当SOC 质量分数≤0.75 g/kg 时，其反射率均值为0.38，当SOC 质量分数＞32.72 g/kg 时，其反射率均值为0.29。

2）通过光谱变换可以明显提高相关系数，SOC含量与原始光谱反射率呈负相关性，－0.62＜相关系数（r）＜－0.07；经SG-SNV-1st Der 预处理后，通过极显著性检验（P＜0.01）的波段数达到414个，主要集中在487～575、725～998和1 464～1 514 nm范围内，在788、800 与1 768 nm 波长处的相关系数均大于0.80。

3）SPA 的降维效果优于相关系数法，光谱建模时SNV 非线性方法的效果优于PLSR 线性回归模型。光谱经SG-SNV-1st Der预处理后，运用SPA结合SVM 模型能很好地估算研究区SOC 含量，其训练集与验证集R2分别为0.79 和0.81，RMSE 分别为5.61和3.16 g/kg，RPD为2.25，RPIQ为2.53。