SPA算法^[18]是一种基于变量信息的变量降维技术, 它利用向量的投影分析来寻找含有最低冗余信息的变量组合, 能够有效地消除光谱波长共线性、 奇异性和不稳定性的影响, 使向量间的共线性达到最小, 减少建模所用变量的个数, 降低模型复杂度。 对于光谱矩阵X_{n× m}(n为样本数, m为光谱变量数), 首先设定待选特征变量个数H, 然后执行以下步骤:

(1)初始迭代t=1时, 在光谱矩阵中任选一列向量x_j, 记为x_k(0), k(0)为所选变量的最初位置(j=k(0), 1≤ j≤ m), 则其他剩余变量位置的集合定义为s

s={j, 1≤j≤m, j∉{k(0), ..., k(H-1)}}(1)

(2)计算剩余列向量x_j(j∈ s)在所选向量x_k(t-1)构成的正交向量空间中的投影

P=I-xk(t-1)(xk(t-1))T(xk(t-1))Txk(t-1)(2)

xj=Pxj(3)

其中, I为单位矩阵; P为投影算子。

(3)提取最大投影值变量arg[max(‖ Px_j‖ )], j∈ s, 加入所选变量集。

(4) t=t+1, 如果t< H, 则返回步骤(2)循环计算。

当循环终止时, 得到的变量集合{x_k(0), x_k(1), ..., x_k(H)}就是选取的特征波长集合。 因为迭代的第一个变量x_k(0)是随机选取的, 因此令光谱中的每个波长都作一次初始变量, 进行上述迭代, 每次迭代选取H个变量, 即可得到n× H维矩阵X={X₁, X₂, ..., X_n}^T, 此矩阵为基于n个初始变量的迭代所选取的n个候选变量集。 然后对每个变量集进行PLS交叉验证, 得到交叉验证均方根误差RMSECV_j(1≤ j≤ n), 取最小的RMSECV所对应的k(0)和所选出的变量组合, 即为最终筛选出的最优组合。

SPA方法在特征波长选择方面具有一定的优势, 但存在两方面的不足: 一是候选变量个数的确定没有标准, H过小会造成所选变量不能包含光谱中的大部分关键信息, 而由于变量之间的共线性, H又不能设置过大以至于超过独立变量的个数, 因此需要寻找取得最优H值的方法; 二是将每一个光谱变量都作为初始变量进行迭代, 得到备选变量组合, 之后要对每组变量进行PLS交叉验证, 当光谱变量较多时, 算法的效率不高, 需要寻找合适的降维方法, 在SPA之前对光谱数据进行预降维。 针对H确定的问题, 采取的改进方法如下:

(1)令H在某个合理的范围[M, N]内变化, 得到N-M组最小RMSECV_i(1≤ i≤ N-M)和对应的变量组合。

(2)从中选取最小RMSECV_i所对应的H以及变量组合。

这样获取的H值具有区域范围的最优性, 称之为自适应(Self Adaption)SPA, 简称SPA_sa。

针对SPA预降维的问题, 提出一种改进的平均影响值(MIV)方法。 MIV方法^[19]最初由Dombi提出, 用于表示神经网络中权重矩阵的变化, 且可用于评估输入变量对神经网络模型性能的影响, 通常应用于质谱分析、 生物医学中^{[20, 21]}。 本研究中的最终建模方法采用BP神经网络, 使用MIV方法可以精确给出波长变量对种子活力等级影响程度的排序, 然后使用SPA方法从排序靠前的波长中进一步进行优选, 达到提高算法效率的目的。 MIV的具体步骤如下, 对于光谱矩阵X_{n× m}(n为样本数, m为光谱变量数)有:

(1)首先采用全谱BP建模, 训练出预测准确率超过90%的模型。

(2)将训练集T第i个样本的第j个波长的吸光度分别增加和减少10%, 得到新的样本集 Tj+i和 Tj-i。

(3)将 Tj+i和 Tj-i作为全谱模型的输入进行预测, 得到的平均预测值代入式(4), 计算出样本变动对输出的影响变化值I Vji(Impact Value)。

IVji=BPfull(Tj+i)-BPfull(Tj-i)(4)

(4)对各样本所得的I Vji求平均, 即为各波长的平均影响值MIV_j, 如式(5)所示

MIVj=1N∑I=1n|IVji|(5)

(5)对得到的平均影响值进行排序。

MIV方法可以得到波长变量的平均影响值排序, 但如何确定序列前多少个波长作为有效波长不同文献的方法不一。 李大虎等^[22]采用相对贡献率来表征某个因素的MIV值对于全部因素MIV总和的百分比。 在其研究中, 全部因素只有9个, 以相对贡献率超过10%作为特征筛选的依据可以有效地实现BP神经网络输入特征的筛选。 但对于数目众多的近红外全谱数据来说, 每个波长的相对贡献率很小且相互之间的差异不大, 以相对贡献率作为筛选指标不够妥当。 因此, 提出相对距离比δ _j作为选择标准

δj=2×|MIVj-MVI¯|MIVmax-MIVmin(6)

当δ _j大于某个数值时的变量即选入备选变量集。 其中, MIV¯为全部波长的平均MIV, MIV_max为最大MIV, MIV_min为最小MIV。

将优化的MIV算法MIV_opt对全谱BP模型进行平均影响值排序及选择, 得到的光谱变量作为SPA_sa算法的输入, 由此降低SPA算法的循环次数, 提高波长选择的效率。 耦合MIV_opt-SPA_sa算法是一种基于变量信息的波长选择方法, 在MIV_opt阶段, 将非线性的BP模型预测用于平均影响值排序, 剔除了非线性模型无关的波长变量; 在SPA_sa阶段, 使用相对全谱变量较少的波长变量, 并采用线性的PLS作为交叉验证, 剔除了与线性模型无关的波长变量, 并进行自适应的筛选, 得到最优的波长变量组合。

实验用种子样本购自大庆市萨尔图萨中种子公司, 为黑龙江省农垦科学院作物所玉米育种研究室杂交培育的垦粘一号玉米品种。 实验用仪器如图1所示。 (a)是近红外光谱采集仪器, 为德国Bruker公司Tango近红外光谱仪, 采用积分球漫反射测量方式, 分辨率为8 cm^-1, 样品和背景的扫描时间均为32 s, 谱区范围11 550~3 950 cm^-1, 每条光谱采集的数据点数为1 845个; (b)是科文KW-TH型种子老化(恒温恒湿)实验箱; (c)是上海菁华公司JA2003N高精度电子天平, 精确到1 mg; (d)是近红外光谱仪配套的IN312-SHD0型量杯。 光谱分析及建模采用的软件主要使用挪威CAMO公司的UnscrambX10.3和美国MathWorks公司的Matlab R2020。

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 近红外光谱仪(a), 老化实验箱(b), 电子秤(c)和IN312-SHD0量杯(d)Fig.1 Near infrared spectrometer (a), Aging test box (b), Electronic scale (c) and IN312-SHD0 measuring cup (d)

种子在自然条件下的贮藏时间越长, 种子的活力和生活力下降的越快。 研究表明, 人工加速老化与自然老化对种子内部物质含量及结构的影响差异不大, 且发芽情况相近。 将种子置于干燥(湿度低于10%)、 低温(温度10~20℃)的环境中保存备用, 实验前对种子进行筛选, 清除干瘪、 瘦小、 损伤以及腐坏的种子, 选出健康、 饱满的种子总计5 000 g。 Tango采集颗粒状样本时, 要求样本的容量要达到量杯容量的2/3以上, 以此确定每个样本的种子质量为(37.0± 0.3) g, 用高精度电子秤量出。 将种子共分为5组(D0, D2, D4, D6, D8), D0组样本15个, 不进行老化处理; 其余各组每组样本13个, 进行不同程度的老化处理, 将样本装入尼龙袋中并编号。 根据《国际种子检验规程》中对玉米种子人工加速老化测定的规定, 将样本放入高温高湿老化箱中, 薄层平铺于老化箱的网架上进行老化, 温度设为41 ℃, 相对湿度设为99%, 5组样本的老化时间分别为0, 48, 96, 144和192 h。

使用积分球漫反射测量方式采集光谱数据, 为了扩大样本数目, 将每个样本重复装样3次(每次装样都要将样本翻动摇匀)测3条光谱取平均。 所有样本光谱采集的环境条件相同: 温度22 ℃, 相对湿度30%。 采样点数为1 845个, 开始波数为11 542.16 cm^-1, 结束波数为3 926.249 cm^-1, 采样间距为4.119 cm^-1。

采集共得到402条光谱数据, 如图2所示。 可以看出, 不同老化时间的样本光谱的整体趋势、 波峰位置高度相似, 属于高相似度样本分类问题。

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 402个玉米种子样本光谱图Fig.2 Spectra of 402 corn seed samples

测量的样品光谱中除了包含样品的真实信息还包括与仪器响应、 测试条件和光的散射等有关的背景信息^[23], 这些信息导致了光谱噪声和基线漂移。 因此, 在建立种子活力检测模型之前, 进行光谱预处理以削弱各种背景信息对真实光谱的影响、 降低模型的复杂度并提高模型的稳健性是十分必要的。 在进行预处理方法选择时, 首先使用高斯滤波(guassian filter, GS)、 卷积平滑(Savitzky-Golay, SG)平滑、 多元散射校正(multiplicative scatter correction, MSC)、 标准正态变量变换(standard normal variate, SNV)方法及其组合对原始光谱进行预处理, 然后使用Kennard-Stone(KS)法将预处理后的光谱按3:1划分为训练集和预测集, 最后建立全谱BP预测模型并根据模型性能确定种子活力分级所采用的预处理方法。

由表1可知, 在光谱平滑处理时, SG平滑消除随机噪声的效果较GS好, SNV消除样品颗粒大小和表面散射光的影响优于MSC, 组合预处理方法SG-SNV的模型表现最优, 训练集的准确率达到92.53%, 预测集的准确率达到90.18%。 组合预处理之后的光谱如图3所示, 与原始光谱相比减弱了噪声干扰和光谱散射问题, 吸收峰的位置更加清晰。 在波长8 210, 6 846, 5 182, 4 737和4 366 cm^-1处有5个显著的吸收峰。 吸收峰在光谱的低频部分更为频繁, 吸光度随着波数的减少而增加。

表1

Table 1

表1(Table 1)

表1 不同预处理方法的BP全谱模型建模结果 Table 1 Results of BP full spectrum prediction model with different pretreatment methods 预处理方法 训练集 预测集 准确率/% 交叉熵 准确率/% 交叉熵 GS 82.14 0.026 3 76.25 0.031 4 SG 89.22 0.018 7 83.34 0.024 8 MSC 83.52 0.025 1 79.56 0.029 1 SNV 87.69 0.018 9 83.47 0.023 8 SG-SNV 92.53 0.013 4 90.18 0.016 2

表1 不同预处理方法的BP全谱模型建模结果 Table 1 Results of BP full spectrum prediction model with different pretreatment methods

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 SG-SNV预处理之后的光谱曲线Fig.3 Spectra after SG-SNV pretreatment

首先利用经典SPA对全谱数据进行特征选择, 设置选择变量数20个, 采用PLS交叉验证择优, 其中校正集样本282个, 预测集样本120个。

图4显示了SPA根据交叉验证均方根误差(RMSECV)从1 845个变量中选出的20个波长变量的位置。 从图中可以看出, 选中的波长集中在光谱第一个波峰附近, 只有几个波长分布在其他几个吸收峰附近, 这些特征波长所携带的信息量明显缺失, 此时的交叉验证误差RMSECV为0.714 7, SPA波长选择的时间为12.125 3 s。 当设置选择变量数为26时的交叉验证误差RMSECV为0.685 0, SPA波长选择的时间为14.417 5 s; 当设置选择变量数为35时的交叉验证误差RMSECV为0.669 1, SPA波长选择的时间为16.824 3 s。 由此可见, 设置不同的变量数得到的最小RMSECV是变化的, 总的趋势是随着选择变量数的增加而减小, 但什么时候RMSECV达到最小, 很难通过不系统的抽样式的SPA选择来确定。

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 设置20个选择变量的SPA所选择的变量Fig.4 Selected variables of SPA with assigning value 20 to the number of variables selected

下面采用SPA_sa进行波长变量选择, 设定H的变化范围为[1, 80], 校正集和预测集样本的设定保持不变, 变量的选择结果如图5(a)所示, 此时选中波长在各个吸收峰的附近都有分布, 提取的信息比较均衡。 在H变化过程中, 以选择的变量数为横坐标, 以RMSECV为纵坐标, 绘制RMSECV随变量数变化的趋势图, 如图5(b)所示。 当选择变量数增加时, RMSECV最小值是逐渐减小的, 当变量数达到47时, RMSECV值达到最小0.6217; 当变量数继续增加时, RMSECV趋向于稳定。 当变量数增加到接近80, 继续增加会引入与预测值无关的波长变量或具有较大噪声的变量, 此时RMSECV会急剧增加。

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 SPAsa特征波长选择结果(a)和SPAsa波长选择中RMSECV随设定变量的个数的变化(b)Fig.5 Selected variables of SPAsa (a) and the change of RMSECV with the setting number of selected variables in SPAsa (b)

由于SPA_sa需要在H的一定范围内反复进行SPA操作, 因此其运算时间相当于多次SPA的时间累加。 当光谱数据量较大时, 算法的运行时间较长, 因此对其进行预降维十分必要。 下面首先对11 542.16~3 926.249 cm^-1范围的光谱数据进行MIV平均影响值计算, MIV影响值随波长的分布如图6所示。

图6

Fig.6

	Figure Option ViewDownloadNew Window
	图6 全谱数据的MIV值分布Fig.6 Distribution of MIV values of full spectrum data

由图6可见, 不同波长变量的MIV数值差异较大, 为了去除与种子活力信息无关或相关性较小的波长, 根据式(6)计算各个波长对应的相对距离比, 选取的波长变量作为SPA_sa的备选光谱数据。 为了找到最佳的D, 建立选取的波长数、 BP模型预测的准确率随D变化的关系, 如图7所示。

图7

Fig.7

	Figure Option ViewDownloadNew Window
	图7 MIVopt预降维 (a): 选择的波长数随相对距离比的变化; (b): 预测准确率随相对距离比的变化Fig.7 MIVopt pre-dimensionality reduction (a): Number of selected wavelengths varies with MIV relative distance ratio; (b): Prediction accuracy varies with MIV relative distance ratio

从图7可以看出, D值的范围从0.05开始逐渐增加, 随着D值增加, 选择的波长数逐渐减少, BP模型的预测准确率开始时逐渐增大, 而后逐步降低。 当D值在0.35附近时, 预测准确率达到94.31%, 当D值取0.40时, 预测准确率变为94.28%。 让D值在[0.35, 0.40]范围内以步长0.01变化, 求得对应的预测准确率。 当D为0.37时达到预测准确率最大值94.79%, 此时选择的变量数目为644个。 以这644个变量作为SPA_sa的备选光谱数据, 设定SPA_sa的优选变量个数小于75个, 校正集和预测集样本数目保持282和120个不变, 经过SPA_sa共筛选出变量37个, 此时PLS模型的RMSECV为0.504 9, 运算时间为14.357 s。

如图8所示, 筛选出的特征波长主要集中在7个波峰附近。 在4 000~4 500 cm^-1波段, 以4 235为中心, 左右分布着4 165和4 358等几个特征波长, 此波段为玉米脂肪C— H基团的吸收峰。 在4 600~5 500 cm^-1波段, 以5 000为中心, 两侧各有一个波峰, 分布着4 787, 4 844, 4 976, 5 182和5 285等特征波长, 此波段为玉米蛋白质N— H基团及淀粉O— H基团的合频吸收峰区。 在5 500~7 500 cm^-1波段, 有两个波峰, 一个以5 620为中心, 两边分布着5 615和5 903等特征波长; 一个以6 719为中心, 两侧分布着6 397, 6 546, 6 854和6 941等特征波长, 此吸收峰是水份的倍频吸收区。 在8 000~9 000 cm^-1波段, 以8 333为中心, 分布着7 971, 8 045和8 910等特征波长, 此波段为玉米淀粉甲基C— H基团二级倍频的吸收谱带。 在9 000~11 000 cm^-1波段有一个波峰, 以10 030为中心, 两侧分布着9 437, 9 787, 10 200, 10 360和10 510等特征波长, 此波段是淀粉甲基C— H基团三级倍频及组合频的吸收谱带。

图8

Fig.8

	Figure Option ViewDownloadNew Window
	图8 MIVopt-SPAsa筛选出的特征波长分布Fig.8 The distribution of wavelengths selected by MIVopt-SPAsa

由上述分析可知, MIV_opt-SPA_sa优选得到的特征波长分布与玉米种子生化物质构成有着高度的一致性, 具有明显的物理意义, 可以体现玉米老化过程中种子内部物质组成的变化; 实现光谱数据的大幅度降维, 是一种有效的基于变量信息的特征提取方法。

BP神经网络能学习和存贮大量的输入-输出模式映射关系, 而无需揭示描述这种映射关系的数学方程。 其学习规则为梯度下降法, 通过反向传播不断调整网络的权值和阈值, 使网络的误差平方和最小^[24]。 本研究采用三层的BP网络, 隐层采用sigmoid激活函数, 输出层采用softmax损失函数, 训练算法采用比例共轭梯度反向传播算法。 样本集划分方法采用KS法, 训练集和预测集的比例为3:1, 将402个样本数据划分为训练集(282)、 预测集(120), 最大迭代次数设为1 000。

为了对比全谱、 MIV、 SPA_sa、 耦合MIV_opt-SPA_sa以及目前为许多学者青睐的CARS波长提取方法对模型性能的影响, 建立了5个BP模型, Full-BP、 MIV-BP、 SPA_sa-BP和MIV_opt-SPA_sa-BP以及CARS-BP。 待输入的光谱数据均经过SG-SNV预处理, 模型的最佳隐层节点数根据经验公式和数据实测综合确定。 所依据的经验公式为

m=n+l+α(7)

式(7)中, m为隐层节点数; n为输入层节点数; l为输出层节点数, α 为1~10之间的常数。

表征模型效率的评价指标是运算时间; 表征模型精度的指标有准确率和交叉熵(cross-entropy, CE), 每个模型都经过50次训练, 各个指标取平均来表征模型最终性能, 5种模型的性能对比如表2所示。

表2

Table 2

表2(Table 2)

表2 5种模型的性能对比 Table 2 Performance comparison of 5 models 模型 输入节点数 隐层节点数 输出层节点数 运算时间/s 训练准确率/% 预测准确率/% 最佳交叉熵 Full-BP 1 845 49 5 0.253 1 90.3 89.4 0.081 241 MIV-BP 644 30 5 24.523 96.7 95.3 0.052 381 SPAsa-BP 46 14 5 101.224 94.4 93.6 0.094 531 MIVopt-SPAsa-BP 37 11 5 14.382 99.2 99.1 0.007 892 CARS-BP 173 18 5 97.226 97.5 97.3 0.013 425

表2 5种模型的性能对比 Table 2 Performance comparison of 5 models

从表2可以看出, 全谱模型的总的运算时间最少, 其次是MIV_opt-SPA_sa-BP, 效率最低的是SPA_sa-BP。 这是因为具有相同输入变量的SPA_sa算法与BP的运算时间相比要大得多, 对于具有1 845个输入的BP来说, 其运算时间一般是几十毫秒, 而SPA_sa-BP一般要达到100 s左右, 而具有1 845个输入的MIV_opt-SPA_sa-BP的运算时间一般是十几秒。 从模型的准确性和稳健性来看, MIV_opt-SPA_sa-BP模型的准确率可达99%以上, 最佳交叉熵为0.007 892远远小于另外4个模型。