云南山地胶园土壤有机质高光谱估算

陈桂良，刘忠妹，许木果，黎小清，丁华平，杨春霞

（云南省热带作物科学研究所，云南景洪 666100）

0 引言

云南是中国种植面积最大、产胶最多、单产最高的优质天然橡胶生产基地[1]。云南橡胶园为典型的山地胶园，地形地貌复杂，小气候复杂多变，成土母质丰富，土壤养分空间差异明显[2-3]，植胶初期土壤肥力较高，但在长期植胶后，土壤退化严重[4-5]。土壤有机质是反映橡胶园土壤养分状况的重要指标，其含量变化受到广泛关注[6-8]。传统化学分析方法存在周期长、成本高、污染环境等不足，而高光谱技术凭借快速、简便、无污染等优势，为土壤有机质含量的快速测定开辟了新的途径[9-10]。土壤有机质含量快速、高精度估算模型的建立，可以更好地指导云南山地胶园精细化生产管理，对云南天然橡胶产业的高质量发展意义重大。

国内外学者普遍认为土壤光谱信息中存在有机质敏感波段，利用土壤有机质敏感波长建立土壤有机质含量估算模型，可以简化模型，保持甚至提升模型精度[11-13]。为突出土壤有机质光谱信号，通常会对原始光谱反射率进行适当变换，然后进行特征波长优选，再建立相应的预测模型，以提高模型的稳定性和精准度[14-17]。有关光谱反演土壤有机质含量的研究中，模型的建立大多数以经验统计方法为主，由于不同区域土壤类型多样且成因复杂，使得高光谱特征和估算模型的差异很大，难以建立统一的有机质估算模型[18]。为了提高模型的有效性，针对特定区域建立区域化的有机质光谱预测模型是常用手段[19]。很多学者对国内典型区域进行了土壤有机质高光谱估算研究，如东北黑土[20]、干旱区典型绿洲[21]、红壤区[22]、华南地区[23]、黄河三角洲[24]等。目前还未见利用高光谱技术估算云南山地胶园土壤有机质的研究报道。

笔者以云南典型山地胶园东风农场为研究区，基于获取的土壤有机质和光谱反射率数据，在优选光谱变换模式和特征波长的基础上，采用多元线性回归、偏最小二乘回归和支持向量回归，构建云南山地胶园土壤有机质含量的高光谱最优估算模型，以期实现云南山地胶园土壤有机质含量快速检测。

1 材料与方法

1.1 研究区与土壤样品采集

云南省景洪市东风农场位于景洪市西南部，在21°30′00″—21°46′12″N、100°34′48″—100°49′48″E 之间，海拔600～900 m 之间，成土母质主要是花岗岩、千枚岩、砂页岩和老冲积物，土壤主要为砖红壤和砖红壤性红壤，下辖6 个分场，现有天然橡胶1.02 万hm2。2019 年9 月，在东风农场选取了103 个GPS 定位采样点，根据地形的差异，在每个GPS定位采样点采集1～3个橡胶树保护带土壤样品(0～20 cm)，共采集到土壤样品225个。图1为研究区位置及采样点分布图。

图1 研究区位置及采样点分布图

1.2 土壤有机质含量与光谱反射率测定

土壤样品带回实验室后，将每个土壤样品混合均匀，剔除树根、石块等杂质，经室内自然风干后研磨并过0.15 mm 孔径筛，用于土壤有机质和光谱反射率测定。采用重铬酸钾容量法测定土壤有机质[25]。采用FieldSpec4光谱仪（美国ASD公司产）进行土壤光谱反射率测定，整个测定工作在一个同等于暗室的实验室进行。用培养皿（半径10 cm、深1.5 cm）装满土壤样品，并将表面刮平，放在反射率近似为0的黑色橡胶垫上。光源为50 W卤素灯，光源入射角度为45°，距土壤样品表面70 cm。传感器探头采用25°视场角，位于土壤样品表面垂直上方15 cm处。光谱测定之前先进行白板校正，测定期间每隔15 min 进行白板校正，每个土壤样品采集10条光谱曲线，算术平均后得到该土壤样品的光谱反射率数据。

1.3 研究方法

1.3.1 光谱预处理及建模集划分 剔除噪声较大的350～399 nm波段，剩余的400～2500 nm光谱反射率数据用于本研究的建模研究。将光谱反射率重采样为5 nm间隔(400、405、410、…、2495、2500 nm)。 采用Kennard-Stone(KS)算法对225 个土壤样本进行优选，选出150个作为校正集，剩余75个作为验证集，校正集和验证集的土壤样品有机质含量统计见表1。

表1 用于模型校正和验证的土壤样品有机质含量统计

1.3.2 光谱变换为了消除光谱数据无关信息和噪声，提高建模精度，首先采用3 种方法对重采样后的光谱反射率R进行光谱变换处理，即倒数对数变换[log(1/R)]、多元散射校正(MSC)、标准正态变换(SNV)，然后对光谱反射率R以及3 种变换形式光谱数据进行Savitzky-Golay(SG)平滑或导数变换，以得到最佳的光谱变换模式。

1.3.3 特征波长筛选方法

（1）竞争适应重加权采样。竞争适应重加权采样(competitive adaptive reweighted sampling，CARS)，通过优选偏最小二乘回归模型回归系数绝对值相对较大的波长变量，剔除模型回归系数绝对值相对较小的波长变量，获得模型交叉验证均方根误差最小的优选变量子集[26-27]。

（2）连续投影算法。连续投影算法(successive projections algorithm，SPA)，是在光谱矩阵中应用变量投影操作寻找含有冗余信息最低、共线性最少的光谱特征变量组，最大程度地避免光谱信息重叠、简化模型结构、提高建模的速度和效率[28]。

（3）CARS-SPA。将CARS与SPA 2种方法结合用于特征波长筛选，首先用CARS 方法获得优选的变量子集，然后再应用SPA 方法对CARS 方法获得的变量子集进一步优选。

1.3.4 建模方法

（1）多元线性回归。多元线性回归(multiple linear regression，MLR)是由多个自变量的最优组合通过回归预测因变量的一种回归方法，在已知多组自变量和所对应的单组因变量时，可以对它们之间的关系进行很好的拟合和预测[29]。

（2）偏最小二乘回归。偏最小二乘回归(partial least squares regression，PLSR)广泛用于高光谱数据回归建模，是建立稳健线性光谱定量校正模型的通用方法。PLSR 集主成分分析、线性回归分析和典型相关分析的优点于一身，适于处理自变量存在多重共线性的回归问题[30]。

（3）支持向量回归。支持向量回归(support vector regression，SVR) 是支持向量机(support vector machine，SVM)的重要应用分支，用于回归模型的构建。SVM 通过核函数将低维非线性问题转换成高维的线性问题，通过结构风险最小化原则提高泛化能力，在保证最小化样本的同时缩小了模型泛化误差的上界，广泛应用于分类与回归[31-33]。

1.3.5 模型验证模型验证采用决定系数(R2)、均方根误差(RMSE)和相对分析误差(RPD)，R2越接近1，RMSE越小，RPD越大，说明模型预测效果越好。1.5＜RPD＜2 表明模型只能对样品进行粗略估算，2.0＜RPD＜2.5表明模型具有较好的估算能力，2.5＜RPD＜3.0 表明模型具有很好的估算能力，RPD＞3.0 表明模型具有极好的估算能力[18]。

2 结果与分析

2.1 土壤有机质含量与光谱反射率

参照GB/T 29570—2013 附录C 胶园土壤养分含量正常指标[34]，将所有样本按有机质含量高低划分为高(＞25 g/kg)、中(20～25 g/kg)、低(＜20 g/kg)3 种有机质含量水平，并计算每种有机质含量水平土壤样品的平均反射率光谱，发现400～2500 nm 波段范围内反射率与有机质含量负相关，且不同有机质含量水平土壤的反射率光谱曲线有相似的反射特征，总体呈现先增加后降低的抛物线型（图2）。在可见光波段，受土壤发色团和有机质本身黑色的影响，存在较宽的吸收波段，反射率总体较低，随着波长的增加而急剧上升；在900 nm 附近出现典型的Fe3+的吸收谷；在1400、1900、2200 nm 附近有3 个明显的吸收峰，深度略有差别，这可能是黏土矿物中含有的水分子和羟基的吸收带[35]。

图2 3种有机质含量水平土壤的平均原始反射率光谱

2.2 光谱变换模式筛选

对重采样后的光谱反射率R进行光谱变换处理后，对光谱反射率R以及3 种变换形式光谱数据进行SG 平滑或导数变换模式优选。优选方法如下：基于SG 平滑或导数变换后的全波段光谱数据与有机质含量数据，采用PLSR 及留一交叉验证方法建模，其中，SG 平滑或导数变换的滤波窗口在3～101 的奇数中筛选，多项式次数1～9 中筛选，导数阶数在0、1、2 中优选，PLSR主成分数在1～20中筛选，按照交叉验证均方根误差(root mean square error of cross-validation，RMSECV)选择最佳的光谱变换模式。

log(1/R)、MSC 和SNV 变换在Unscrambler 9.7 中实现，SG 平滑或导数变换和PLSR 模型的建立通过MATLAB 2015 编程实现。研究发现MSC 和SNV 变换对模型影响不大，log(1/R)变换能够显著提升模型性能，log(1/R)变换后再进行SG平滑或导数变换，可以进一步提升模型预测能力，但不同的导数阶数、滤波窗口大小对模型的预测能力影响较大（图3）。log(1/R)结合SG平滑变换是最佳的光谱变换模式，其中，SG平滑变换模式为导数阶数0，SG 滤波窗口5，多项式次数2 或3（表2）。

表2 不同光谱变换和不同导数阶数下最优模型预测效果

图3 不同光谱变换下各阶导数不同滤波窗口大小对应的最小交叉验证均方根误差

采用PLSR 模型和留一交叉验证法，分别基于全波段的光谱反射率R和最佳光谱变换光谱数据，建立土壤有机质高光谱估算模型，利用验证集对模型进行验证，相比光谱反射率R，基于最佳光谱变换光谱数据建立的模型精度显著提高，RPD从2.433 增加到2.795（图4）。表明应用本研究的最佳光谱变换模式对光谱反射率R进行光谱变换处理，能够有效突出土壤有机质光谱信号，有利于提高模型的稳定性和精准度。

图4 基于全波段的光谱反射率R(a)和最佳光谱变换光谱数据(b)建立的PLSR模型预测效果

2.3 特征波长筛选

为进一步提高模型的稳定性和精准度，基于筛选得到的最佳光谱变换光谱数据与土壤有机质含量数据，选择CARS[36]、SPA[37]、CARS-SPA 3 种特征波长筛选方法，提取得到3个特征波长组合（表3）。

表3 不同特征波长筛选方法提取的特征波长

2.4 模型构建与验证

基于校正集的土壤有机质含量与最佳光谱变换模式下的光谱数据，分别选用表3所示的3种特征波长组合对应的光谱特征变量作为自变量，土壤有机质含量作为因变量，并采用MLR、PLSR和SVR 3种方法构建土壤有机质高光谱估算模型。其中，MLR 和PLSR 模型在Unscrambler 9.7 中实现，结合留一交叉验证方法进行建模。SVR 模型利用LIBSVM 软件[38]实现，选取SVR 类型为epsilon-SVR，核函数类型为radial basis function，通过网格搜索和留一交叉验证进行参数优选。利用验证集对所建模型进行精度验证，模型的预测效果如图5 所示。从模型预测效果来看，CARSSVR最优。

图5 不同模型的预测效果

3 结论

本研究尝试对土壤光谱反射率R进行多种光谱变换后，再结合SG平滑或导数变换模式筛选，得到最佳的光谱变换模式为log(1/R)结合SG 平滑变换，其中，SG 平滑变换模式为导数阶数0、SG 滤波窗口5、多项式次数2 或3。基于最佳光谱变换光谱数据与土壤有机质含量数据，选择CARS、SPA、CARS-SPA 3 种方法进行特征波长筛选，并采用MLR、PLSR和SVR 3种方法构建土壤有机质高光谱估算模型。结果表明，MLR、PLSR、SVR 3 种模型类型中，最佳模型分别为CARS-MLR、SPA-PLSR、CARS-SVR，RPD分别为2.745、2.617、2.947。CARS-SVR 模型估算精度最高，R2、RMSE、RPD分别为0.897、3.990 g/kg、2.947，该模型RPD位于2.5～3.0之间，具有很好的估算能力。本研究应用高光谱技术，对云南省景洪市东风农场典型山地胶园土壤有机质含量进行估算，为云南山地胶园土壤有机质含量的快速估算提供了参考。

4 讨论

本研究采用3种方法[log(1/R)、MSC、SNV]对5 nm间隔光谱反射率R(400、405、…、2495、2500 nm)进行光谱变换处理，并对光谱反射率R以及3 种变换形式光谱数据进行SG 平滑或导数变换模式优选，得到最佳的光谱变换模式为log(1/R)结合SG 平滑变换，其中，SG 平滑变换模式为导数阶数0、SG 滤波窗口5、多项式次数2 或3。研究发现，log(1/R)变换能够明显提升土壤有机质估算模型性能，这与方少文等[22]的研究一致。log(1/R)变换后再进行SG 平滑或导数变换，可以进一步提升模型预测能力，但不同的导数阶数、滤波窗口大小对模型的预测能力影响较大。

基于全波段的最佳光谱变换光谱数据建立土壤有机质PLSR模型，并利用验证集对模型进行验证，RPD达2.795。基于筛选得到的最佳光谱变换光谱数据与土壤有机质含量数据，选择CARS、SPA、CARS-SPA 3种方法进行特征波长筛选，并采用MLR、PLSR 和SVR 3 种方法构建土壤有机质高光谱估算模型，利用验证集对各模型进行验证，结果表明，MLR、PLSR、SVR 3 种模型类型中，最佳模型分别为CARS-MLR、SPA-PLSR、CARS-SVR，RPD分别为2.745、2.617、2.947。研究发现，PLSR 作为一种建立稳健线性光谱定量校正模型的通用方法，在全波段模型中具有一定优势，但在特征波段模型中，不如SVR，甚至不如MLR，与王涛等[39]的研究结果一致，这可能也与特征波长提取方法有关。SVR作为一种非线性建模方法，特别是基于特征波长建模，在土壤有机质高光谱估算研究中具有较大潜力。

为了尽可能地提高光谱数据与土壤有机质含量之间的相关性，有效突出土壤有机质光谱信号，本研究尝试对土壤光谱反射率R进行多种光谱变换后，再结合SG平滑或导数变换模式筛选，这为土壤有机质含量高光谱估算模型的建立提供了新思路。本研究采集的土壤样品包括花岗岩、千枚岩、砂页岩和老冲积物等成土母质，也是云南山地胶园主要母质类型[5]。由于建立的有机质估算模型主要针对云南山地胶园，为了模型的区域普适性，并未按母质分类建模，而是将全部土壤样品混合建模，但仍然取得了较好的效果。本研究仅用同一农场的土壤样本数据对模型进行了验证，在后续研究中，还需在云南山地胶园开展更为广泛的模型验证及优化。