程介虹, 陈争光*, 张庆华

(1.黑龙江八一农垦大学电气与信息学院, 黑龙江 大庆 163319 2.大庆技师学院计算机工程系, 黑龙江 大庆 163524)

土壤作为地球表面五大典型地物(水域、农田、山脉、城市、土壤)之一，其可见光-近红外反射光谱特性研究对当前精准农业、土壤数字制图、土壤资源遥感调查等工作有着重要作用, 备受众多科研人员的关注[1]。土壤有机质(soil organic matter,SOM)泛指土壤中来源于生命的物质，是土壤固态部分的重要组成成分，能促进植物的生长发育，改善土壤的物理性质。研究表明，近红外光谱对土壤生物量碳和生物量氮、全碳和全氮以及有机质含量均有较好的评估效果[2]。由于近红外光谱检测技术具有快速、成本低、无污染等优点，所以被用作一种“快速无损”的技术来检测SOM含量，可以避免大量的土壤取样。

近年来，王昶等[2]结合偏最小二乘法对长江中下游粮食生产区水稻土的土壤质量(全碳、全氮、碳氮比和土壤pH)进行评估，决定系数均大于0.9，效果较好，可进行快速测定；李伟等[3]用偏最小二乘和人工神经网络方法对近红外光谱建模，可测定土壤碱解氮、速效磷和速效钾含量；吴金卓等[4]结合偏最小二乘和主成分分析法，建立土壤全氮和碱解氮含量测定的定标模型，在野外快速测定土壤全氮和碱解氮含量。徐夕博等[5]通过主成分波段信息判别提取6个主成分，建立MLR和BPN预测模型，发现BPN模型预测精度优于MLR模型。单海斌等[6]通过对土壤光谱进行倒数、对数、一阶微分、倒数的一阶微分、对数的一阶微分变换，运用单相关分析法提取土壤光谱特征波段，建立多元回归方程。Ge等[7]通过外部参数正交化建立偏最小二乘回归模型，对不同含水率土壤样品的粘土和有机碳含量进行了预测。Kirshnan等[8]采用逐步多元线性回归对土壤反射率数据进行分析，确定土壤有机质含量预测的最佳波长为564和624 nm。 Brown等[9]以PLSR和BRT回归证实大量特征波段都可用来对SOM进行估算。Wang等[10]利用不同预处理方式改善光谱与土壤性质间相关性。Bilgili等[11]使用PLSR和MARS对土壤有机质含量进行建模研究，取得良好的预测精度。

上述研究多基于全谱或特征波段建模来估测土壤的养分信息，但近红外光谱所采集的样本变量数目庞大，影响建模速度，并且近红外波段是合频和倍频的吸收区域，光谱信息重叠严重。若选用波段信息，变量中会出现干扰、共线性变量，所以仅仅建立一部分波段与土壤某个属性之间的线性关系无法达到最优估算结果[5]。

为了解决这一问题，波长选择成为近红外光谱分析中非常重要的一步。特征波长选择是从全谱数据中提取部分涵盖有用信息的光谱，去除噪声光谱及无用信息，建立一个更为简约、稳定的近红外光谱模型。波长选择可大大减少变量数目，加快模型的计算效率，提高模型的稳健性[12]。本文以土壤为研究对象，Rinnan等[13]已证明近红外光谱与SOM含量密切相关，相关系数超过0.9。为提高其预测精度，本文选取连续投影算法(successive projections algorithm,SPA)消除共线性信息变量，偏最小二乘法(interval partial least squares,IPLS)选择特征波长区间间隔，竞争自适应重加权采样法滤除冗余信息变量(competitive adaptive reweighted sampling, CARS),从三个不同方面有代表性算法进行特征波长选择，建立模型，比较其预测性能。

1 材料与方法

1.1 样本数据

土壤样本数据来自于Quality & Technology[13]，是Havstrom等在瑞典北部阿比斯科(68°21′N，18°49′E)进行长期野外试验时，模拟了亚北极气候的预测变化对SOM的影响而得到的样本数据。该数据为一组公开的数据集，基于此数据的方法改进和所建模型具有一定的可比性。样本数据包含两部分，样本的NIR光谱及其化学性质。共计108个样本。

样本光谱的波长范围为400～2 500 nm，采样间隔为2 nm。样本的化学特性为土壤有机质含量。本文采用SOM含量作为因变量进行近红外光谱数据建模预测分析。

1.2 变量选择算法

1.2.1SPA算法步骤 SPA是由Araujo等[14]提出的一种消除变量间共线性的波长选择算法，该方法通过提取全谱的几个特征波长消除原始光谱中冗余的信息[15]。其主要原理是利用向量的投影分析，假设已给出初始波长k(0)和所需提取波长数目N，算法步骤[14,16]如下。

Step 0:在第一次迭代(n=1)之前，将校正集Xcal的第j列光谱数据赋值给xj,j=1,…,J，J为总波长数。

Step 1:没有被选择的列向量的集合记为S,S={j,1≤j≤J,j∉{k(0),…,k(n-1)}}

Step 3:令k(n)=arg(max‖Pxj‖),j∈S

Step 4:令xj=Pxj,j∈S

Step 5:令n=n+1，如果n

End:最后得到的波长为{k(n);n=0,…,N-1}

其中，k(0)和N的选择是很关键的一个步骤。k(0)取值在1～J之间变化，而N的变化范围为1≤N≤Mcal，Mcal为校正集样本数，这是可以由SPA选择的最大波长数。

1.2.2CARS算法步骤 CARS是由Li等[17]在2009年提出来的一种结合蒙特卡洛(Monte Carlo,MC)采样与偏最小二乘(partial least squares,PLS)模型回归系数的一种波长选择方法，通过消除变量来建立一个高性能的校准模型。该方法可以对无信息变量进行有效去除，最终选择出对预测目标起关键作用的变量。

具体步骤为：首先，随机抽取80%的样本做为校正集建立PLS模型，计算回归系数的绝对值和每个波长点对应的权重，然后，利用指数衰减函数(exponentially decreasing function，EDF)和自适应重加权采样法(adaptive reweighted sampling，ARS)对变量进行选择，通过交叉验证的方法计算交叉验证均方根误差(RMSECV)。N次MC采样后选择N个子集，得到N个RMSECV，选择RMSECV最小的波长子集，该子集所包含的变量即为最优变量组合[18]。

1.2.3IPLS算法步骤 IPLS是由Norgaard等[19]于2000年提出的一种特征波长区间选择方法，它的主要作用是在不同的光谱分支中提供有关信息的整体图像，从而聚焦于重要的光谱区域，并消除来自其他区域的干扰[19]。具体步骤为：将全谱分成若干相同宽度的光谱子区间，然后在每个子区间建立偏最小二乘回归模型，通过比较区间模型和全谱模型的RMSECV来比较预测性能，RMSECV最低时的区间模型即为最优波段。

1.3 模型建立与评价

偏最小二乘回归(partial least squares regression,PLSR)建模方法在近红外光谱分析中使用较为广泛，它考虑光谱数据与性质之间的内在联系，模型更加稳健，但模型建立过程复杂、抽象，计算速度较慢且繁琐。当各变量高度线性相关时，用PLSR法非常有效[20]。而多元线性回归(multiple linear regression,MLR)的特点适用于线性较好的数据，不考虑参数之间的相互干扰、计算简单、公式含义清楚，产生的模型比PLSR更简单，更容易解释。故本文基于所选取的特征波长，通过建立MLR模型分析其预测能力。

MLR是一种常规的校正方法，直观简单，且具有良好的统计特性，应用非常普遍。MLR采用最小二乘法进行回归计算，模型一般为

近红外光谱模型的预测能力主要通过变量数(n)、校正相关系数(Rc)、校正均方根误差(RMSEC)、预测相关系数(Rp)、预测均方根误差(RMSEP)指标来评价。其中，R取值越接近1，RMSEC和RMSEP越接近0，模型的拟合性越好，预测精度越高[4]。

1.4 数据分析工具

所用软件包括MATLAB R2015b和The Unscrambler X 10.3 (64-bit)。The Unscrambler软件具有广泛的数据预处理方案,确保数据都适合进行多因素分析，并且具有先进的回归、分类和预测建模工具。本文光谱数据的预处理、建模分析在Unscrambler软件中实现，各变量选择方法、图形的绘制在MATLAB中实现。

2 结果与分析

2.1 光谱数据特征

各样本400～2 500 nm区间原始近红外光谱变化趋势相同，无明显异常样本。为了消除基线漂移、消除背景影响、提高光谱分辨率，使用窗口为5的二阶Savitzky-Golay(S-G)求导法对光谱进行预处理，结果如图1所示。

图1 原始光谱及预处理后光谱Fig.1 Original spectrum and pre-processed spectrum

将108个样本通过SPXY(sample set portioning based on joint x-y distance)算法分为75%训练集和25%预测集，建模集包含81个样本，预测集包含27个样本， 土壤有机质含量统计数据结果如表1所示。划分后建模集的SOM含量范围涵盖预测集的SOM含量，证明建模集具有代表性。

图1 原始光谱及预处理后光谱Fig.1 Original spectrum and pre-processed spectrum

分别对原始光谱和预处理后的光谱建立PLSR模型，结果如表2所示。经S-G求导预处理后模型的预测性能有所提高，故将预处理后的数据作为后续波长选择的数据。

表2 原始光谱和预处理后的光谱PLSR模型参数对比Table 2 Comparison of PLSR models of the raw spectrum and pretreated spectrum

2.2 特征波长选取

2.2.1SPA变量选择结果 由SPA方法选择特征波长建立MLR模型，选取RMSEP最小值对应的波长个数为最终的特征波长个数。图2A中的正方形标记所示为SPA多元线性回归模型选择的变量数，RMSEP的最小值为1.214 4，之后RMSEP值基本达到稳定，此时选择特征波长数为6个。图2B中方块所对应的6个点即为SPA选择的最佳特征波长,分别为636、688、1 934、2 224、2 330、2 426 nm。

2.2.2CARS变量选择结果 利用CARS 算法对光谱数据进行特征波长选择，MC抽样的次数设置为100次。图3A所示为变量变化的趋势，由于指数衰减函数的作用，在刚开始进行MC采样时进行快速选择，变量减少的速度非常快。然后在精细选择中，使用自适应重加权采样法来选择基于回归系数的关键变量，变量减少的速度较慢，直到得到最优子集。图3B所示为十折交叉验证RMSECV的变化趋势，随着MC采样次数的递增，PLS模型的十折交叉验证RMSECV的值逐渐减少，可知剔除了大量无信息变量。但当达到最小值后(虚线处)又逐渐增大，表明光谱数据中一些与预测相关的信息被剔除，导致模型的性能变差。图3C所示为回归系数路径，星号垂线处是为第55次抽样，此时的RMSECV值最小为2.032 6，有34个回归系数波长为非零值，即特征波长共计34个变量。

A:模型包含的变量数和模型RMSEP的关系;B：特征波长点。A: Relationship between the number of variables and RMSEP of the model; B: Characteristic wavelength points.图2 SPA的特征波长选择Fig.2 Characteristic wavelength selection of SPA

A:变量变化趋势；B:十折交叉验证RMSECV的变化趋势；C:回归系数路径。A: Variation trend of variables; B: Trend of of RMSECV under 10-fold cross validation; C:Regression coefficients path.图3 CARS的特征波长选择Fig.3 Characteristic wavelengths selection of CARS

2.2.3IPLS变量选择结果 利用IPLS算法对光谱数据进行特征波长选择，将预处理后的光谱数据依次等分成20、25、30、35、40个区间，然后在每个子区间建立PLS模型，比较这些区间模型和全谱模型的RMSECV，选择效果最佳的子区间。IPLS选择结果如表3所示。

从表3可以看出，当光谱区间划分为30时，对应的RMSECV最小。因此,将光谱数据等分成30个区间进行特征波长选择。由图4A可知，当全谱的主成分数为3时，RMSECV基本平稳。由图4B可知,建立在第15区间的PLS模型的RMSECV最低。因此选择1 392～1 460 nm区间，共计35个波长点作为波长选择结果。

表3 IPLS子区间优选结果Table 3 Results of IPLS subinterval optimization

A:RMSECV与全谱模型的主成分数；B:各子区间模型和全谱模型RMSECV对比，各柱状图上斜体数字代表子区间的主成分数。A:RMSECV versus PLS components for global model；B:Comparison of RMSECV between each interval model and full spectrum model, and the italic number on each rectangle represents the number of principal components of the subinterval.图4 IPLS的特征波长选择 Fig.4 Characteristic wavelengths selection of IPLS

2.3 模型建立与比较

本研究拟利用MLR对全谱以及经三种特征波长选择算法后得到的数据进行建模，比较预测精度，但由于全谱含有1 050个波长点，数据过于庞大，而MLR只适用于对波长点少于样本数的数据进行建模，所以对全谱进行PLS回归；而使用IPLS方法进行特征波长选择，所选择的波长点是连续的，存在共线性，所以对IPLS选择的35个特征波长点进行PLS回归；对SPA选择的6个特征波长点和CARS选择的34个特征波长点依然建立MLR模型，得到的模型的校正、预测相关系数和校正、预测均方根误差的值如表4所示。

表4 不同波长选择方法所得模型对比Table 4 Comparison of models obtained by different wavelength selection methods

由表4可以发现，经过特征波长提取后的预测精度均高于原始光谱，进一步证明了对全谱进行特征波长选择的重要性。通过对比三种特征波长选择算法的预测效果，发现SPA仅挑选全谱的6个波长点，RMSEP为1.214 4，预测效果最好。SPA不仅降低了建模的复杂度，还大大提高了模型的预测精度，这种方法可以有效地应用于特征波长的提取中，简化模型复杂度，提高模型计算效率。

图5分别为对预测集样本不同模型的土壤有机质实测值和预测值相关图。其中，SPA-MLR模型得到的方程如下。

图5 不同模型下土壤有机质的实测值和预测值相关性Fig.5 Correlation between measured and predicted values of soil organic matter under different models

y=-12.00-5 410.73x(636 nm)-375.76x(688 nm)-15 957.55x(1 934 nm)+23 717.22x(2 224 nm)-50 740.52x(2 330 nm)+45 713.68x(2 426 nm)

3 讨论

王岩[21]利用此数据，通过IPLSR、PLSR对土壤有机质含量建立模型，预测模型的RMSEP分别为2.787 3和2.709 8，而本文使用三种变量选择方法建立模型的预测能力均得到不同程度的改善，进一步证明了对全谱进行特征波长选择的重要性。文中SPA预测能力最优，CARS次之。基于SPA选择的特征波长MLR模型预测集的相关系数为0.970 2，RMSEP为1.214 4。并且仅选取全谱的6个波长点，大大降低模型复杂度，提高模型效率。这与吴龙国等[22]应用不同波长选择方法对土壤含水率进行建模研究时得到的结论一致，相比于CARS，SPA-MLR模型的预测精度更优。Chen等[23]利用IPLS和SPA结合MLR建模，对谷物的蛋白质、总碳水化合物和粗脂肪的含量进行预测，发现SPA-MLR的预测精度优于IPLS，这也与本文结果一致，说明SPA-MLR在某些光谱模型中具有一定的有效性。

在SPA方法中，通过不断提高MLR模型的预测能力来选择变量。与PLS方法相比，SPA投影只用于选择变量，不会修改原始数据，因此保留了光谱数据和样本化学性质间的关系。SPA在进行波长选择前，对各变量进行标准化处理，各波长点向量模相同，因此进行投影时选择最大投影向量的过程其实就是选择相互正交的向量的过程，所以SPA可以将变量之间的共线性降到最低。而MLR对线性较好的数据回归效果好，所以可以解释SPA-MLR的预测性能最好。

土壤的NIR与SOM之间存在较好的相关性，其预测集的相关系数大于0.94，更加证明了土壤的光谱特性可以对SOM含量进行预测。本研究通过三种不同特征波长提取算法对土壤有机质含量进行建模，取得了较满意的成果。现有研究大多对土壤有机质、重金属及氮磷钾的含量进行测定，土壤的微生物含量对农业生产也有着影响，今后应对土壤的微生物含量的预测进行进一步研究。