1 引言

土壤有机质(SOM)含量是衡量土壤质量的一个重要指标,尽管有机质仅占土壤总量的很小一部分,但在促进植物生长发育、改善土壤物理性质等多方面的作用显著^[1-2]。 高光谱遥感能够快速、大范围、无损地获取土壤信息^[3-4],但高光谱包含几百乃至上千个变量,其中有些变量与待测样品信息并无关系,因此对原始光谱进行特征变量或敏感波段选择,研究其是否能够替代全波段,获得较高的预测精度,减少模型工作量及提高模型效率非常有意义。早期如刘焕军等^[5]以松嫩平原的黑土、草甸土、黑钙土等耕层土壤有机质含量为研究对象,采用相关分析法得出土壤有机质含量的敏感波段为520 nm。卢艳丽等^[6]对东北平原的黑土和潮土光谱采用多元线性回归方法提取特征波段,得出土壤有机质含量的敏感波段为550~830 nm。近年来竞争性自适应重加权(CARS)、稳定竞争性自适应重加权(sCARS)、连续投影算法(SPA)、遗传算法(GA)和迭代保留有效信息变量(IRIV)等变量选择方法被用于土壤有机质估算研究^[7]。朱亚星等^[8]利用去除有机质实验结合无信息量消除法(UVE)和CARS变量筛选方法得出土壤有机质的敏感波段为561~721 nm与1920~2280 nm。

筛选土壤有机质含量的光谱敏感波段是简化模型和提高模型预测能力的关键。 Thielebruhn等^[9]对连续小波变换后的可见光-近红外(Vis-NIR)光谱采用CARS进行变量选择,结合偏最小二乘回归(PLSR)模型预测SOM含量,决定系数R²从0.81提高到0.93。林志丹等^[10]对采自涡阳县的130个砂姜黑土土壤样本预测其SOM含量,采用GA提取23个特征变量建立主成分回归模型,决定系数R²达到0.93。以上研究表明基于变量选择算法挑选特征变量建模,使得模型精度和稳健性提高。

除了光谱变量,预测模型的选择对于Vis-NIR光谱的预测精度也至关重要^[11]。 Viscarra Rossel等^[12]通过对比PLSR模型与多种数据挖掘算法如SVM和RF模型,发现数据挖掘算法精度高于PLSR模型。RF模型不仅有较快的拟合速度,而且对异常值和噪声的敏感度更低,稳健性更好,在模型拟合能力方面优于其他算法^[13]。葛翔宇等^[14]基于CARS耦合机器学习预测土壤含水量,相较于线性模型,其决定系数R²从0.617提高到0.918。不同变量选择算法挑选的特征变量并不相同,这会影响模型精度,目前文献中多数为多种变量选择方法结合线性模型,与SVM、RF模型相结合的研究并不多见。因此,需要将多种变量选择方法与多种回归方法相结合对土壤属性进行预测,为设计更高效率、便携性更好的光谱仪器提供理论依据。

本文以青海省湟水流域401个表层土样为研究对象,原始光谱经预处理后,采用sCARS、SPA、GA、IRIV和sCARS-SPA挑选光谱特征波段,并基于全波段和特征波段分别建立PLSR、SVM和RF预测模型,为利用Vis-NIR 光谱分析技术快速无损地估测农田土壤有机质含量提供方法支持。

2 数据与方法

2.1 研究区概况

湟水流域位于青海省东北部,为青藏高原一个特殊自然地理单元和生态屏障。地理位置介于36°02'—37°28'N,100°42'—103°04'E之间,流域面积为1.62×10⁴ km²。流域地势整体西北高,东南低,东西长,南北窄,海拔(DEM)在1655~4860 m之间。气候为高原干旱、半干旱大陆性气候。湟水流域土壤类型以栗钙土、黑钙土、灰钙土、山地草甸土、高山草甸土为主,土壤肥沃,为青海省主要粮食生产基地,主要种植春小麦、油菜、马铃薯、青稞、燕麦和玉米等农作物。

2.2 土壤样品采集

分别于2015、2016年10—11月进行土壤采样,此时庄稼已收割完毕,共采集401个湟水流域表层(0~20 cm)土壤样品。考虑到土壤有机质易受到坡度、坡向、土壤母质等自然因素,耕种、施肥灌溉管理措施等人为因素的影响,河谷区域的水浇地及地形平坦区域采用“梅花状”5点采样,坡耕地采用“S”型7点采样,采样点分布见图1。将土样装入密封袋中编号,并于室内自然风干、研磨、过100目筛(150 μm),将过筛后的土壤样品分为两份,分别用于有机质含量测试分析和土壤光谱数据采集。有机质含量采用重铬酸钾-外加热法测定,并应用浓度梯度法确定建模校准集与验证集样本,即将401个土样有机质含量从高到低排序,按2∶1比例划分校准集和验证集样本。土壤有机质含量(质量分数)最大值(Max)、最小值(Min)、平均值(Mean)和标准差值(SD)统计见表1。

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

Fig. 1. Location of the study area and distribution of soil sampling sites

下载图片 查看所有图片

2.3 土壤光谱数据采集及预处理

土壤样品的Vis-NIR光谱数据采集设备为美国ASD FieldSpec 4光谱仪,光谱范围为350~2500 nm,室内土壤光谱测量参照文献[ 15]。去除边缘噪声较大的350~400 nm和2401~2500 nm,并参考文献[ 16]将土壤有机质含量分为高、中、低和非常低4类,分别为大于30 g·kg^-1、23~30 g·kg^-1、12~22 g·kg^-1和小于12 g·kg^-1。土壤样品原始光谱曲线如图2(a)所示,在401~1100 nm波长区域,光谱反射率值随波长变化增加比较明显,且不同等级有机质含量的平均光谱反射率值差异较大;而在1400~2400 nm波长区域,光谱反射率值随波长变化增加平缓,且不同等级有机质含量的反射率曲线差别不明显。

光谱采集过程中易受仪器噪声、土壤颗粒分布不均匀及测量随机误差等因素干扰,使所测样品光谱中含有光谱噪声,影响预测模型精度^[17-18],故采用多元散射校正(MSC)、中值滤波(MF)和一阶微分(1st Derivative)对原始光谱依次进行预处理。图2(b)为经MSC-MF-1st Derivative预处理后的光谱反射曲线。从图中可以看出,光谱预处理后400~900 nm区域不同等级有机质含量光谱曲线之间差异性减小,光谱细节特征更加明显。

图 2. 土壤样品反射率曲线。(a) 原始光谱;(b) MSC-MF-1st Derivative预处理光谱

Fig. 2. Spectral reflectance curves of soil samples. (a) Raw spectra; (b) spectra after MSC-MF-1st Derivative pre-processing

下载图片 查看所有图片

2.4 特征变量选择方法

2.4.1 sCARS变量选择

sCARS以变量的稳定性作为变量选择衡量指标,增强了变量选择的稳定性,并延续CARS变量筛选流程^[19]。具体步骤如下:

1) 计算每个波长变量的稳定性值C_i,C_i定义为

Ci=b¯is(bi) i=1,2,…,P,(1)

式中:C_i为M次蒙特卡罗采样中第i个变量的稳定性值, b¯i和s(b_i)分别为M次采样中第i个变量的均值和标准偏差,P为变量数。显然, b¯i值越大,s(b_i)值越小,第i个变量的重要性越强;

2) 使用强制波长选择和自适应性重加权采样方法(ARS)筛选出一组稳定性较好的变量子集,筛选出的变量数占全波段的比率由指数衰减函数(EDF)计算;

3) 循环步骤1)和步骤2),最终得到K个变量子集,并基于这些变量建立PLSR模型,对应的均方根误差(R_MSE)值最小的变量子集为最后的特征变量,其中K为sCARS的循环次数。sCARS通过在MATLAB 2010b中编写代码实现。

2.4.2 SPA变量选择

SPA的优点是从光谱矩阵中选择最小信息冗余的波长组合,降低波长之间共线性的影响,从而降低模型的复杂度,提高模型的稳定性和准确性^[20]。设光谱矩阵为X(A×B),m为提取的特征变量个数,n为迭代次数。

1) 任意选择X中的一列j,把建模集的第j列赋值给x_j,记为变量x_k(0);

2) 未被选入的列向量位置集S为

S={j,1≤j≤B,j∉{k(0),…,k(n⁃1)}},(2)

分别计算x_k(n-1)对未被选中列变量的投影,记为 Kxj,K为投影算子:

Kxj=xj(xjsxk(n⁃1))xk(n1)(xk(n⁃1)sxk(n⁃1))-1;(3) 3)k(n⁃1)=arg(maxKxj),j∈S,n=n+1,(4)

如果n<N,则令x_k(n-1)选入列向量,返回2)再计算;

4) 采用多元线性回归方法对新选择的变量集{x_k(0),x_k(1),…,x_k(m-1)}进行评估,最后得到m个最优波长子集{x_k(0),x_k(1),…,x_k(m-1)}。

2.4.3 GA变量选择

GA通过模拟自然进化过程搜索最优解,具有较高的自适应和全局优化能力^[21]。GA具体步骤为:

1) 波长编码:对波长进行0/1二进制编码,0表示未选中该波长,1表示选取该波长,0和1随机组合生成一条染色体;

2) 选择初始群体:记初始群体个数为N,每一个染色体长度(波长数)为q,随机产生N个q位的0/1随机组合二进制群体;

3) 适应度函数:采用PLSR模型中交叉检验均方根误差(R_MSECV)为适应度函数;

4) 遗传操作:包括选择、交叉和变异3种方式,一般而言,选择轮盘赌法,变异概率为0.01,交叉概率为0.5;

5) 终止条件:重复4),若达到设定的最大繁殖代数,则进化过程中得到的具有最大适应度的个体作为最优解输出,计算终止。

2.4.4 IRIV变量选择

IRIV由随机子集生成、子集模型建立、模型参数分析3个环节构成,是一种基于模型集群分析策略的波长选择算法^[22]。相对于一般的波长选择算法,IRIV具有在变量选择时呈现软收缩的特点,因此在保留有效变量方面更为稳妥。需注意的是IRIV虽可以更好地保留变量间的协同效应,但在迭代运算过程中需要建立大量子模型,这使得该算法的计算量较大。

2.4.5 sCARS结合SPA变量选择

sCARS的优点是速度快,最终选出的特征变量的化学意义也比较容易解释,SPA挑选的特征变量冗余度最低,共线性最小。但SPA计算量较大,且选中的波长子集中很可能会纳入一些无关信息甚至是干扰信息,sCARS和SPA的联合使用不仅可以最大程度地降低光谱信息冗余,还可以降低无效波长对SPA计算过程的干扰^[23]。

2.5 模型精度评价

采用均方根误差(R_MSE)、 决定系数(R²)、验证集标准差(SD)与验证集相对分析误差(R_PD)对预测模型进行精度评价。验证集均方根误差(R_MSEVAL)越小,R²、R_PD越大,模型预测效果越好。当R_PD≥2时,表明模型预测性能较好,可对样本进行有效估测;当1.4≤R_PD<2时,表明模型可对样本进行粗略估测;当R_PD<1.4时,表明模型不能对样本进行估测。

3 结果

3.1 特征变量筛选

图3为采用sCARS挑选特征变量过程。从图3(a)中可以看出,随着sCARS迭代次数的增加,所保留的波长数量逐渐减少,且减少速度由快到慢;图3(b)中十折交叉验证均方根误差(R_MSECV)值呈现由大到小再由小到大的变化趋势;图3(c)为波长变量稳定性轨迹图,图中每条曲线为每个变量的稳定度随迭代次数的变化趋势,将R_MSECV值最小时的最优变量子集数量用星号标记。sCARS挑选特征变量过程中具有“粗选”和“精选”两个阶段,当sCARS运行次数为27次时, R_MSECV值最小,此时共选择51个特征变量,仅占总变量数的2.55%。

图 3. sCARS变量筛选流程。(a)变量变化趋势;(b)十折交叉均方根误差;(c)变量回归系数

Fig. 3. Variable selection process by sCARS. (a) Changing trend of variables; (b) 10-fold RMSECV values; (c) regression coefficients of variables

下载图片 查看所有图片

结合光谱数据, SPA运行时设定特征变量最小波段数为5,最大波段数为100,当模型均方根误差最小时对应的变量数为最佳波段数,提取的变量为最佳波段。本研究中利用SPA共选择出5个最优特征变量,分别为1361、1758、1909、2049、2213 nm。

图4为所有变量被选择的频率,图中3条横线表示选择频率阈值,频率阈值越大,选择的变量数越少。横线以上表示保留的特征变量,用于模型构建,横线以下是未被选中变量,不用于建模分析。本研究中采用GA,从原始光谱中共选取186个特征变量,占Vis-NIR光谱全部变量的9.3%。

图 4. GA特征变量筛选过程

Fig. 4. Characteristic variable selection process by GA

下载图片 查看所有图片

IRIV能够将把单一指标进行波长硬性删除的传统策略转化为较为柔性的波长保留策略,进行波长选择时呈现软收缩的特点,因此在保留有效波长方面更为稳妥。本研究中IRIV设定光谱矩阵为X(500×2000),利用五折交叉验证建立PLSR模型选择特征变量,最大主因子数为10,最终 IRIV共选择特征变量63个,占全部变量的3.15%。

在sCARS处理后,再采用SPA挑选特征变量,光谱变量数由51个减少到17个,如图5所示。由图可知随着筛选变量数量的增加,R_MSE值先迅速下降,当变量数为17时,R_MSE值趋于稳定状态。图5(a)中空心小方块表示采用SPA得到的最优变量子集个数,图5(b)为采用SPA挑选的最优特征变量点在一条光谱曲线上的分布,分别为440、444、458、460、493、561、938、976、1905、1916、2025、2253、2305、2314、2318、2390、2398 nm。

图 5. 预处理光谱sCARS-SPA特征变量筛选过程。(a)模型变量数;(b)变量指数

Fig. 5. Characteristic variable selection process by sCARS-SPA from the pre-processing spectrum. (a) Number of variables in the model; (b) variable index

下载图片 查看所有图片

图6为5种变量筛选方法挑选特征变量分布图。由图可知,5种变量筛选方法挑选的特征变量主要分布在近红外光谱区域。这可能是由于近红外光谱记录的主要是土壤中有机分子基团(O—H、N—H、C—H及C=C等)在不同波长近红外光的倍频与和频吸收信息,而含氢基团是大多数有机物分子结构的基本单元。张娟娟等^[24]系统分析了我国东部、中部5种不同土壤类型近红外光谱反射率,认为与SOM含量相关性较好的波段范围主要集中在1350~1420 nm、1860~1890 nm、2000~2350 nm 3个区域。Krishnan等^[25]基于1130、1350、1398、2210 nm 4个波段建立了SOM估测模型。Bendor等^[26]的研究也表明1000~2500 nm的近红外波段区域存在SOM相关的吸收波段。

图 6. 不同变量筛选方法挑选特征变量分布

Fig. 6. Distribution of characteristic variables with different variable selection methods

下载图片 查看所有图片

3.2 PLSR模型

表2为全波段和特征波段PLSR模型建模结果,其中 Rcal2和 R2val分别为校准集和验证集决定系数,R_MSECAL和R_MSEVAL分别为校准集和验证集均方根误差。由表2可知,基于特征波段的PLSR模型的预测效果均优于全波段PLSR模型,特别是基于sCARS挑选的特征变量建模,建模变量数为51个,仅占全部变量数的2.55%,验证集 R2val为0.883,获得的R_PD为2.9,模型预测效果最佳。图7为sCARS-PLSR模型校准集和验证集样本实测值和预测值的散点图,正方形代表校准集样本点,圆形代表验证集样本点。从图中可以看出,除验证集C17、C18、C21、C30、C32、C45、C98及C259点和校准集的V8、V108点外,其他样本点均分布在1∶1直线两侧,这可能是由于这些样本点存在光谱异常或测量值异常。

图 7. sCARS-PLSR模型预测SOM含量散点图

Fig. 7. Scatter plot for the measured and predicted value by sCARS-PLSR model

下载图片 查看所有图片

3.3 SVM模型

表3为全波段和特征波段SVM模型建模结果。首先,从表中可知,基于全波段建立的SVM模型,验证集R²为0.74,获得的R_PD为1.9,可粗略对样本进行预测;而其校准集R²为0.91,显著高于验证集R²,这可能是由于光谱变量中存在干扰变量或无效变量,导致模型过于复杂,模型出现过拟合。其次,基于sCARS、SPA、GA、IRIV和sCARS-SPA挑选特征波长变量后建立预测模型,验证集R²均显著高于全波段模型,且与校准集的R²较为接近,说明所建模型稳定性较好,同时获得的验证R_PD均大于2,可较好地对样本进行预测,表明建模之前对全波段进行特征变量筛选不仅能够有效提高建模效率,还能提高模型的预测精度。采用SPA挑选特征变量建模,建模变量数为5个,验证集R²为0.889,获得的R_PD为2.9,模型预测效果最佳。图8为SPA-SVM模型校准集和验证集样本实测值和预测值的散点图,从图中可以看出,除验证集C17、C18、C21、C30、C45、C98及C259点和验证集V8、V108点外,其他样本点均离1∶1直线较近。

图 8. SPA-SVM模型预测SOM含量散点图

Fig. 8. Scatter plot for the measured and predicted value by SPA-SVM model

下载图片 查看所有图片

3.4 RF模型

RF模型是一种分层非参数方法,用于估计独立变量和因变量之间复杂的非线性关系。RF模型为了提高模型预测精度并避免过拟合,引入两个随机性因素,即森林中单棵树的分类强度和森林中树与树之间的相关强度。这两个随机性因素的引入,使得RF模型不容易因变量数远大于建模样本数陷入过拟合,且具有很好的抗噪声能力^[27]。RF模型通过MATLAB 2010b编程实现,分类树数目n_tree=500,节点数m_try为 p,p为自变量数目。表4为全波段和特征波段RF建模结果。从表中可知,全波段和特征波段RF模型验证集R²值无明显差异,这说明采用特征波段构建RF模型对模型预测精度的改善不明显,但其构建模型的变量数却显著减少,可提高建模效率。这可能与RF模型具有很好的抗噪声能力,不易因变量数大于建模样本数而陷入过度拟合有关。其中IRIV-RF模型预测效果略好,校准集和验证集R²分别为0.941和0.96,验证集R_PD为4.8,可较好地预测SOM含量。图9为IRIV-RF模型校准集和验证集样本实测值和预测值的散点图。从图中可以看出,IRIV-RF模型校准集和验证集数据点除C17、C18、C21、C98、C117、C132、C150和V8、V10外其他数据点均靠近1∶1直线,且均较PLSR、SVM模型数据点更靠近1∶1直线的两侧。

图 9. IRIV-RF模型预测SOM含量散点图

Fig. 9. Scatter plot for the measured and predicted value by IRIV-RF model

下载图片 查看所有图片

4 讨论

前文采用sCARS、SPA、IRIV、GA和sCARS-SPA模型挑选特征变量。基于特征变量的PLSR和SVM模型精度均高于全波段模型精度,基于特征变量的RF模型精度较全波段模型提高不明显,但建模效率大大提升,进一步说明对全波段进行有效变量选择的重要性。

图10为不同变量筛选算法的PLSR、SVM和RF模型建模结果。由图可知,RF模型的预测效果最佳,优于SVM模型和PLSR模型。SOM含量与光谱之间关系较为复杂,PLSR是一种线性方法,在解决非线性问题时表现不佳,而SVM和RF可较好地解决独立变量和因变量之间复杂的非线性关系;但SVM模型易因较高的频谱噪声引起严重的偏差估计,模型精度降低;RF模型融合了随机特征选择和Bagging算法两大机器学习技术,与传统的分类器算法相比,不但能较好地容忍异常值和噪声,使建立的模型精度较高且具有较好的稳健性,而且能同时处理连续型和离散型数据^[28]。如Viscarra等^[29]采用多种建模方法预测土壤有机碳、黏土含量和pH值,结果显示,RF模型精度要优于PLSR、MARS、ANN和SVM模型。同样Douglas等^[30]基于Vis-NIR光谱分析技术,通过比较PLSR和RF模型预测土壤中总石油烃(TPH)含量,发现RF模型在处理光谱中非线性因素时,与线性PLSR模型相比,提供了更高的预测精度。

图 10. 不同变量筛选方法PLSR、SVM、RF模型建模结果

Fig. 10. Results of PLSR, SVM and RF models with different variable selection methods

下载图片 查看所有图片

异常样本的存在对模型的性能产生一定的影响,表5为剔除C17、C18、C21、C30、C45、C98、C259、V8和V108异常点后的模型精度。从表中可知,在剔除这些异常值之后sCARS-PLSR、sCARS-SVM和sCARS-RF模型精度均有提高,特别是sCARS-RF模型验证集R²提高到0.988,R_PD达到7.8,预测SOM含量性能最佳。因此,在以后的工作中将深入研究异常值剔除对模型精度的影响,如高洪智等^[31]基于随机抽取一致性算法(RANSAC)建立有机质含量估算模型,相关系数r从0.891提高到0.963。图11为剔除异常值前后sCARS-RF模型散点图。从图中可知,剔除异常值后,sCARS-RF模型校准集和验证集数据点均较为均匀地分布在1∶1直线的两侧。

图 11. 人工剔除异常值前后sCARS-RF模型散点图。(a)剔除异常值前;(b)剔除异常值后

Fig. 11. Scatter plots for the measured and predicted value by sCARS-RF model before and after artificially eliminating outliers. (a) Contain outliers; (b) eliminate outliers

下载图片 查看所有图片

5 结论

筛选土壤有机质的光谱响应波段是简化模型和提高模型预测能力的关键。本研究以青海省湟水流域表层土样为研究对象,原始光谱经预处理后,采用sCARS、SPA、GA、IRIV和sCARS-SPA 5种算法挑选光谱特征波段,并引入PLSR、SVM和RF模型对优选的特征变量构建SOM预测模型。PLSR模型和SVM模型精度较全波段均有一定提高,而RF模型R²没有明显变化,对模型精度的提高帮助不大,但模型变量数却显著减少,降低了模型复杂度,提高了建模效率。 同时RF模型精度明显优于SVM模型和PLSR模型,其在处理多维特征数据和抗噪声能力方面展现出独特的优势。5种变量选择算法中,基于IRIV挑选的特征变量建立的RF模型反演SOM含量精度和拟合效果最优,较全波段PLSR模型,预测集R²从0.835提高到0.96,R_PD从2.5提高到4.8,实现了对SOM含量的精准预测。变量选择算法耦合回归模型,在保证模型精度的同时,大大降低了模型的复杂度,提高了建模效率,为采用光谱分析技术快速无损估测耕地SOM含量提供了技术支撑,应用前景广阔。