拉曼光谱数据由奥谱天成(厦门)光电有限公司制造的便携式拉曼光谱仪进行采集, 仪器激发波长为785 nm, 测量范围在200~3 400 cm^-1区间, 分辨率为6.58 cm^-1。 光谱强度稳定性低于5%, 符合稳定性标准要求。 数据处理软件采用MathWorks 公司的Matlab R2017b 实现。

实验选用黑龙江省2020年主要种植的4种高蛋白大豆种子为研究对象, 品种分别是黑农88、 黑农98、 绥农71以及绥农76, 由黑龙江省农业科学院提供。 从每个品种中随机选取40 粒成熟、 无损伤的种子作为该品种样本集, 4个品种共计160个样本用于光谱数据采集。 拉曼光谱仪采集数据时, 激光功率设为500 mW, 积分时间设为5 000 ms。 对于每个样本连续进行3次采集, 将平均光谱作为该样本的原始光谱数据。

1.3.1 随机蛙跳算法

随机蛙跳(random frog, RF)算法^[9]借鉴可逆跳跃马尔可夫链蒙特卡罗(reversible jump Markov Chain Monte Carlo, RJMCMC)算法的思想, 建立了一种数学上简单并且计算效率高的特征波长选择算法。 RF通过在模型空间中实现固定维度和跨维度移动搜索获得一个伪MCMC链, 用于计算每个特征变量的被选概率, 作为该特征变量选择的标准, 主要包括以下三个步骤^[9]:

(1) 参数初始化: K₀为初始变量集 F₀的个数, 在1~P之间随机设置(P为全部变量个数); 迭代次数N为10 000; θ 为正态分布控制参数, 取默认值0.3; ω 为大于1的候选变量个数控制参数, 取默认值3; η 为接受概率, 取默认值0.1。

(2) 候选变量子集构造: 通过N(K₀, 0.3K₀)产生一个随机数, 就近取整确定为候选子集变量个数K^* , 进而构造一个包含K^* 个变量的候选子集F^* , F为包含全部p个变量的集合。

① 若K^* =K₀, 则F^* =F₀。

② 若K^* < K₀, 则先用F₀建立PLS模型, 并从F₀中剔除回归系数最小的K₀-K^* 个变量, 则剩余变量构成候选子集F^* 。

③ 若K^* > K₀, 则先从F-F₀中随机选取3(K^* -K₀)个变量, 并与F₀构成新变量子集T后建立PLS-DA模型, 保留模型中回归系数绝对值最大的K^* 个变量构成候选子集F^* 。

候选变量子集F^* 初选后, 判断其是否被接受进行下一次迭代。 分别用F₀和F^* 构建偏最小二乘分析模型, Err₀和Err^* 为两个模型的预测误差。 如果Err^* ≤ Err₀, 则接受F^* 作为F₁, 否则以0.1Err₀/Err^* 的概率接受F^* 作为F₁。 用F₁代替F₀, 进行下一次迭代计算, 直到结束N次循环。

(3) 变量选择概率计算: 经过N次迭代计算共获取N个变量子集。 N_j表示第j个变量在N次迭代中被选择的次数。 则每个变量被选择概率根据式(1)进行计算

probj=NjN, j=1, 2, …, p(1)

由式(1)可知, 变量越重要被随机蛙跳算法选择的概率就越大。 因此, 该被选概率可以作为衡量变量重要性的指标。 通过设定不同阈值, 选取被选概率大于阈值的变量作为特征波长变量, 对比所建模型性能从而确定最优特征波长。

1.3.2 改进随机蛙跳算法

在采用RF算法提取特征变量时, 由于初始变量集F₀是从原始数据集中随机抽取产生的, 可能会引入无用或干扰波长变量, 使得算法的预测能力和收敛速度降低。 因此, 为了提高初始变量集的有效性, 提出一种基于最小绝对收敛和选择算子(least absolute shrinkage and selection operator, LASSO)的改进随机蛙跳算法。 LASSO^[10]是一种基于普通偏最小二乘法的特征波长选择算法, 通过增加L1范数惩罚, 在预测残差平方和最小准则下, 将模型的回归系数的绝对值之和限定在一个常数范围内, 从而将一些不重要变量的回归系数严格收缩为0。 多元线性回归分析模型如式(2)

Y=Xβ+ε(2)

式(2)中, X=[x₁, x₂, …, x_n]^T∈ R^{n× m}为大豆样本组分变量, Y=[y₁, y₂, …, y_n]^T∈ R^{n× 1}为样本对应品种属性, 其中, m为样本的特征变量数, n为大豆样本数量。 β =[β ₁, β ₂, …, β _m]^T∈ R^{m× 1}是回归系数向量。 LASSO准则的目标函数如式(3)

L(β)=argminβ{‖Y-Xβ‖2+λ1|β|1}(3)

式(3)中, λ ₁为L1正则化系数, 当λ ₁=0时为常规最小二乘法, 当λ ₁逐渐增大时L1正则化影响也越大, 越来越多的系数将接近于0或者等于0, λ ₁的取值将根据RMSECV取最小值时确定。

LASSO算法可以筛选与属性变量最相关特征变量, 采用该方法提取大豆重要特征波长点作为RF初始变量集F₀, 可以有效消除初始变量的随机性, 在此基础上开始迭代计算, 可以改善RF算法所需迭代次数大、 算法收敛慢等问题。 同时, RF算法通常选取概率排序1~10的变量, 或者选取概率大于某一阈值的变量, 因此提取的特征波长往往具有不确定性。 改进方法如下: 首先去除选择概率为0的变量, 对于排序后变量以10个波长点为间隔, 每次增加1个间隔构建波长变量与大豆品种属性的偏最小二乘分析模型。 然后计算每个模型的RMSECV找到最优特征波长以提高模型的预测性能。

1.3.3 ElaticNet特征选择方法

ElaticNet方法结合了LASSO回归和岭回归的思想, 在普通最小二乘回归基础上增加L1惩罚和L2惩罚, 可以同时对变量系数的绝对值和平方项进行压缩。 LASSO回归对所有波长变量的回归系数进行同等程度的压缩, 则可能由于过度压缩导致模型的预测准确度下降。 岭回归^[11]是通过增加L2惩罚, 收缩不重要变量回归系数以达到提高预测准确度的目的, 但不具有变量选择的能力。 Elastic Net结合这两种惩罚函数的特点, 在高维数据的情况下剔除冗余变量, 且一般不会过度压缩回归系数, 从而达到提高模型预测精度的目的。 Elastic Net准则的目标函数如式(4)

L(β)=argminβ{‖Y-Xβ‖2+λ1|β|1+λ2‖β‖22}(4)

令α =λ ₁/(λ ₁+λ ₂), λ =λ ₁+λ ₂, 则式(4)为

L(β)=argminβ{‖Y-Xβ‖2+λ[α|β|1+(1-α)‖β‖22]}(5)

从式(5)中可以看出, 当α =0时为岭回归, 当α =1为LASSO回归, 当α ∈ (0, 1)时为Elastic Net回归, 它是L1和L2惩罚的特殊线性组合, 采用α 的大小来调节它们之间的权重。 当λ =0时为普通最小二乘法, 当λ 逐渐增大时正则化的影响也将越大, λ 的取值将根据RMSECV取最小值时确定。

1.3.4 模型构建方法与评价

PLS是在主成分回归基础上提出的多元校正方法, 在计算主成分时, 考虑主成分方差尽量最大以提取有用信息, 同时还保证主成分与样本属性变量尽可能地相关, 使得所建模型具有较高的预测精度和较强的稳健性。 因此以改进随机蛙跳法选择后的特征波长为输入变量, 使用PLS-DA方法建立大豆品种鉴别模型。 并采用均方根误差(root mean square error, RMSE)、 决定系数(coefficient of determination, R²)和准确率(accuracy)三个指标对模型效果进行评价。 RMSE值越小, 样本预测值与实际值之间的差异越小, R²越接近1, 样本预测值与实际值之间的相关度越高, 表明模型的预测性能越好。 对于识别正确率, 分别将黑农88、 黑农98、 绥农71以及绥农76四个大豆品种样本赋值 1、 2、 3、 4, 根据文献[12], 将PLS-DA鉴别模型的误差阈值设置为± 0.5, 结果分别在1± 0.5、 2± 0.5、 3± 0.5和4± 0.5之间时, 大豆样本识别正确, 否则识别错误。

大豆样本的原始拉曼光谱如图1(a)所示。 标准正态变换(standard normal variate transformation, SNV)^[13]是一种有效的预处理方法, 可以校正样本不均匀带来的散射效应对光谱数据的影响, 预处理效果如图1(b)所示。 再利用KS (Kennard-Stone)算法按照3∶ 1的比例划分大豆样本集, 则校正集样本为120个、 预测集样本为40个。 采用PCA方法对校正集和预测集的共160个大豆拉曼光谱数据进行降维处理, 第一、 第二和第三主成分的贡献率分别为88.56%、 4.65%和3.45%, 前三个主成分的累积贡献率达到96.66%, 说明可以反映原始样本所提供的绝大多数信息。 采用前三个主成分绘制大豆样本校正集和预测集的得分图如图2所示, 可知校正集和预测集的样本分布比较均匀, 校正集能够覆盖预测集, 说明该样本划分方法有效并可以用于建模分析。

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 大豆原始及预处理后拉曼光谱图 (a): 原始拉曼光谱; (b): SNV预处理后拉曼光谱Fig.1 Raw and SNV pre-processed Raman spectra of soybean (a): Raw Raman spectra; (b): SNV pre-processed Raman spectra

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 大豆样本前三主成分得分图Fig.2 Score plot of the top three PCA for soybean samples

2.2.1 RF特征波长优选

对RF算法进行参数初始化: 迭贷次数N设为10 000, 初始变量集个数K₀设为50, θ 、 ω 和η 为默认值。 迭代计算结果如图3所示, 波长变量被选概率处于0~0.747 4之间。 通过人为设定不同阈值, 分别建立PLS-DA鉴别模型, 结果如表1所示。 当阈值设定为0.1时, 则概率大于0.1的波长入选, 共提取198个特征波长点, 此时RMEP和 Rp2分别为0.276 9和0.938 7, 识别准确率为95%。 随着阈值降低, 入选有效波长数增加, 模型预测性能得到提升。 当阈值设定为0.07时, 共有284个特征波长点入选, RMEP和 Rp2分别为0.262 1和0.945 0, 识别准确率为100%, 此时模型预测效果最好。 随着阈值继续增加, RMEP和 Rp2参数性能略有提升, 但识别准确率下降, 表明引入了无效或者干扰信息变量。 由此可知, 应选取0.07为最终设定阈值。

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 RF算法的波长被选概率Fig.3 Selected probability of each wavelength by RF

表1

Table 1

表1(Table 1)

表1 不同阈值的优选波长 Table 1 Optimized wavelengths with different thresholds 序号 阈值 选择波 长点数 RMSEC Rc2 RMSEP Rp2 准确率 /% 1 0.1 198 0.065 4 0.996 6 0.276 9 0.938 7 95 2 0.09 223 0.062 7 0.996 9 0.270 2 0.941 6 95 3 0.08 254 0.117 3 0.989 0 0.272 3 0.940 7 97.5 4 0.07 284 0.072 8 0.995 8 0.262 1 0.945 0 100 5 0.06 318 0.044 6 0.998 4 0.252 6 0.949 0 97.5 6 0.05 380 0.033 8 0.999 1 0.244 0 0.952 4 97.5

表1 不同阈值的优选波长 Table 1 Optimized wavelengths with different thresholds

2.2.2 MRF特征波长优选

采用LASSO算法对大豆全光谱数据进行初选。 正则化参数λ ₁通过10折交叉验证法确定, 首先将预处理后大豆样本数据划分为10份, 每次选取其中9份作为校正集, 剩余1份作为预测集, 采用LASSO方法选择特征波长建模并计算每个模型均方根误差, 经过10次迭代计算后取其平均值作为最终RSMECV值。 当λ ₁=0.003 0时, RMSECV最小为0.195 8, 此时共提取了106个与大豆品种属性最相关的特征变量, 将其作为RF算法的初选变量集F₀, 开始进行迭代计算。 将MRF算法的迭代次数分别设置为1 000、 1 500、 2 000和2 500次, 运算结果如表2所示。

表2

Table 2

表2(Table 2)

表2 不同迭代次数的优选波长 Table 2 Optimized wavelengths with different iteration times 序号 迭代次数 选择波长点数 RMSECVmin 1 1 000 170 0.134 0 2 1 500 250 0.121 7 3 2 000 300 0.111 0 4 2 500 270 0.134 2

表2 不同迭代次数的优选波长 Table 2 Optimized wavelengths with different iteration times

对于不同迭代次数选取的特征波长, 首先去除被选概率为0的波长点, 对剩余波长进行排序后以10个波长点为间隔, 每次增加1个间隔建立PLS-DA模型, 选择RMSECV取最小值时的波长组合为优选特征波长。 由表2可知, 当迭代次数为2 000次时, 不同波长组合的RMSECV计算结果如图4(a)所示, RMSECV值最小时为0.111 0, 最终选取300个有效特征波长点, 如图4(b)所示。

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 MRF特征波长优选结果 (a): 特征波长组合的RMSECV值; (b): 优选特征波长分布Fig.4 Results of MRF characteristic wavelength optimization (a): RMSECV of characteristic wavelength combinations; (b): Distributionof optimized characteristic wavelength

2.2.3 ElasticNet特征波长优选

由式(5)可知, ElasticNet的优化函数L(β )包含参数α (0< α < 1)和正则化系数λ (λ > 0)。 通过10折交叉验证法确定调节参数α , 由表3可知, 当参数α 取0.016时, RMSECV值最小为0.145 6, 此时最优正则化系数为λ =0.031 7。 基于所选最优参数(α =0.016, λ =0.031 7), ElasticNet法通过坐标下降法(coordinate descent)迭代计算3201个大豆拉曼光谱数据的稀疏系数, 非零稀疏系数对应成分变量即为所选特征波长点, 通过ElasticNet法共选择了1 345个大豆特征波长点, 波长分布如图5所示。

表3

Table 3

表3(Table 3)

表3 调节参数α 选择 Table 3 Selection of adjustment parameter α 序号 α RMSECV 选择波 长点数 序号 α RMSECV 选择波 长点数 1 0.01 0.154 7 1 661 6 0.016 0.145 6 1 345 2 0.011 0.157 1 1 600 7 0.017 0.156 0 1 313 3 0.012 0.150 6 1 546 8 0.018 0.156 4 1 288 4 0.013 0.156 2 1 501 9 0.019 0.176 0 1 257 5 0.014 0.155 3 1 431 10 0.02 0.166 9 1 212 6 0.015 0.153 4 1 393

表3 调节参数α 选择 Table 3 Selection of adjustment parameter α

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 ElasticNet选择的特征波长分布Fig.5 Characteristic wavelengths distribution by ElasticNet

将预处理后大豆光谱数据, 分别采用RF、 MRF、 LASSO和ElasticNet算法进行特征波长提取, 将这4种方法得到的特征波长点以及全光谱数据作为预测模型输入变量, 大豆品种属性作为输出变量, 建立大豆品种鉴别的PLS-DA模型, 各模型鉴别结果如表4所示。

表4

Table 4

表4(Table 4)

表4 不同波长选择方法的模型对比 Table 4 Model comparisons of different wavelength selection methods 模型 主成分数 选择波长点数 RMSEC Rc2 RMSEP Rp2 准确率/% Full-PLS-DA 8 3 201 0.121 5 0.988 2 0.256 6 0.947 3 97.5 RF-PLS-DA 8 284 0.072 8 0.995 8 0.262 1 0.945 0 100 MRF-PLS-DA 6 300 0.129 1 0.986 7 0.246 9 0.951 2 100 LASSO-PLS-DA 13 106 0.052 3 0.997 8 0.381 4 0.883 6 87 ElasticNet-PLS-DA 14 1 345 0.012 4 0.999 9 0.253 7 0.948 5 97.5 注: 五种建模方法校正集的鉴别准确率均为100% Note: The identification accuracies of the calibration set of five models were 100%

表4 不同波长选择方法的模型对比 Table 4 Model comparisons of different wavelength selection methods

由建模结果可知, Full-PLS-DA模型的预测结果中RMSEP和 Rp2分别为0.256 6和0.947 3, 识别准确率为97%, 模型预测性能优于LASSO方法。 原因在于LASSO是对所有波长变量的回归系数实行同一尺度的压缩, 导致部分有用波长变量的丢失, 而用于建模的特征波长点仅有106个过于稀疏, 在各种波长优选方法中预测性能最差。 ElasticNet方法是在LASSO基础上增加L2惩罚, 可以同时实现变量回归系数的绝对值和平方惩罚, 从而在特征变量的稀疏性和显著性方面得到平衡。 预测结果中RMSEP和 Rp2分别为加了0.253 7和0.948 5, 识别准确率为97.0%, 优于全光谱和LASSO方法。 相比较而言, RF算法不需要严格的数学公式, 从随机变量开始进行迭代计算, 无需像常规RJMCMC法指定先验分布, 更容易实现并且计算效率高, 对于高维光谱数据是一种非常高效的变量选择方法。 当初始变量集个数K₀设为50、 迭代10 000次、 阈值设为0.07时提取有效波长点284个, 仅占全波长的8.87%, RMSEP和 Rp2分别为0.262 1和0.945 0, 这两个参数性能略差于全光谱和ElasticNet方法, 优于LASSO方法。 而且RF算法的预测准确率高达100%, 优于全光谱、 LASSO方法和ElasticNet方法。 针对RF算法初始变量集的随机性, 提出了改进随机蛙跳算法, 从LASSO变量提取的与大豆品种属性最相关的106个变量开始迭代计算, 在迭代2 000次时算法达到收敛, 有效地减少了无用迭代次数, 显著降低了模型计算复杂度。 并且通过交叉验证法选择了300个有效特征波长点, 改善了RF算法根据阈值选取的不确定性, 预测结果中RMSEP和 Rp2分别为0.246 9和0.951 2, 识别准确率达到100%, 在五种建模方法中性能最优。