1.1.1 数学模型

TPA是一种双群体搜索算法^[16], 模仿团队两个小组(精英组和普通组)的学习和探索过程, 并设计合理的成员更新规则, 逐步提升其评价值以达到全局最优。 实验表明, 该算法能够在较少的计算量前提之下快速寻找全局最优。

算法模型多变量无约束最大化问题可表示为

Y=min{f(x)}x=[x1, x2, …, xn]Txi∈[ai, bi], i=1, 2, …, n(1)

式(1)中, 向量x代表一个成员, 即为包含多个波长的波段; x_i表示该成员的第i个能力因素, 即该波段中的第i个波长点; a_i和b_i分别表示x_i的上下边界值; 而函数f(x)代表该成员x的评价值。 通过更新成员以逐步提升或降低评价值来寻求最优成员。

1.1.2 分组规则

将整条光谱波段均分为P个波段同时确定其评价值, 评价值为该波段与测得的含量理化值进行PLS建模得到的校正集均方差(RMSEC)和相关系数(r)为变量的函数值, 本工作设定的评价值为

f(x)=r1+RMSEC(2)

每一个波段就相当于向量x。 将P个波段按评价值从大到小排列分成N+M个波段(N和M均为整数), N个评价值较高的波段组成精英组{x_e1, x_e2, …, x_en}, M个评价值较低的波段组成普通组{x_p1, x_p2, …, x_pm}。 通过团队进步过程更新组内成员。

1.1.3 行为定义

需要产生一条新生波段x_r, 新生波段可从任意一组产生, 其波长点从当前组随机波段同一波长点中继承。 若新生波段出身自精英组, 且该新生波段的第i个波长点如果是在精英组第n个波段中产生, 那么该波长点需继承精英组第n个波段中第i个波长点; 若新生波段出身自普通组, 且该新生波段的第i个波长点如果是在普通组第m个波段中产生, 那么该波长点需继承普通组第m个波段中第i个波长点。 继承下来的新生波段通过设定概率选择一次学习或者探索行为以更新自身的波长点, 才能成为候选波段x_c。

新生成员x_r如果选择进行学习行为, 则需要向参照目标方向调节。 参照方向产生于精英组和普通组, 分别称为精英组样板e_e和普通组样板e_p, 且样板值取所在组波段波长的平均值。 若新生波段产生于普通组, 则其需向精英组样板调节; 若新生波段产生于精英组, 则向普通组样板的反方向调节, 如式(3)所示。

xc=(1-γ)xr+γeexc=(1+γ)xr-γep(3)

式(3)中, γ 为区间[0, 1]内随机数。 若x_c某个波长点范围越界, 则改用其边界值。

新生成员如果进行探索行为, 则其各波长点x_ri(i=1, 2, …, n)将做随机改变。 并且探索强度逐步减小。 两组新生波段x_r经过探索行为生成x_c的表达式如式(4)和式(5)所示。

xc=[xc1, xc2, …, xcn]T(4)

xci=xri+γite, p(bi-xri), mi=0xri-γite, p(xri-ai), mi=1(5)

式(5)中, γ _i是区间(0, 1)的随机数, m_i为0和1二值随机整数。 收缩系数t_{e, p}为

te, p=1-kKαe, p(6)

式(6)中, K为算法最大迭代次数, k表示当前累计的迭代次数。 收缩指数α _{e, p}表示当新生波段继承自精英组时选取α _e, 当新生波段继承自普通组时选取α _p。

1.1.4 更新规则

若候选波段x_c的评价值高于精英组末位x_ewst的评价值, 则x_c进入精英组, 同时精英组再次进行评价值排序, 精英组末位x_ewst不进入普通组直接淘汰。 这是因为TPA设置了学习行为, 加强了算法的定向搜索和局部搜索能力, 遭到淘汰的精英组末位x_ewst如果进入普通组的话容易导致算法陷入局部最优。 若候选波段x_c的评价值劣于x_ewst但优于普通组末位x_pwst, 还需检查x_c是否由探索行为得到, 若是, 则x_c进入普通组, 淘汰x_pwst。 若不是, 直接丢弃x_c。 这是因为学习行为产生高评价值候选波段的可能性比较大, 且趋同性强, 容易导致普通组波段同化, 降低全局寻优能力。

在TPA算法中, 继承于精英组的新生波段x_r在选择学习行为时, 需要向普通组样板e_p的反方向调节。 但在迭代过程中, 精英组和普通组的成员在不断进行更新优化, 从而在迭代后期, 普通组的样板值e_p也随着精英组样板值增大, 在选择学习行为时, 并不能为继承于精英组的新生波段x_r提供一个良好的调节方向, 容易陷入局部最优。 因此对TPA算法进行如下改进:

1.2.1 改进分组规则

在对P个成员进行分组时, 由原来的N+M模式改为N+M+L模式, 评价值依次降序排列。 N为评价值高的波段组成的精英组{x_e1, x_e2, …, x_en}, M为评价值适中的波段组成的普通组{x_p1, x_p2, …, x_pm}, 新增添的L组为评价值最低的波段组成的垃圾回收组{x_g1, x_g2, …, x_gl}。 再通过团队进步过程更新成员。

1.2.2 学习行为重定义

新增添的L组不参与继承、 学习以及探索行为。 新生波段进行学习行为时, 仍需要向参照目标方向调节, 参照方向更改为分别产生于精英组和垃圾回收组, 称为精英组样板e_e和垃圾回收组样板e_g, 样板值依旧取所在组波段波长的平均值。 即若新生波段产生于普通组, 则其需向精英组样板调节; 若新生波段产生于精英组, 就向垃圾回收组样板的反方向调节, 如式(7)所示。

xc=(1-γ)xr+γeexc=(1+γ)xr-γeg(7)

式(7)中, γ 为区间[0, 1]内随机数。 若x_c某个波长点范围越界, 则改用其边界值。

1.2.3 更新规则的修改

在进行波段更新的时候需要将低评价值的波段回收进L组, 使L组的波段评价值极低。 新的波段更新规则在精英组和普通组不做修改, 当候选波段x_c的评价值低于垃圾回收组首位(垃圾回收组中评价值最高)x_gbst的评价值, 则x_c进入垃圾回收组, 同时垃圾回收组中的x_gbst波段遭到淘汰, 促使垃圾回收组评价值一直保持低的状态, 为继承于精英组的新生波段进行学习行为时提供一个正确的更新方向, 有效地避免算法陷入局部最优。

Step1: 初始化基本参数。 包含总波段成员个数, 精英组、 普通组以及垃圾回收组的波段个数, 新生波段选择学习或探索行为的概率, 精英组和普通组的收缩系数, 迭代次数。

Step2: 分组。 将各波段按照与其相应的理化值PLS建模得到的评价值函数大小降序排列, 评价值函数如式(2)所示, 依次放入精英组、 普通组和垃圾回收组, 使得精英组中评价值最高, 普通组次之, 垃圾回收组最低。

Step3: 产生新生波段。 新生波段随机选择从精英组或普通组中产生, 其波长点从当前组随机波段同一波长点中继承。

Step4: 新生波段选择学习或探索行为生成候选波段。 生成一个随机数, 若随机数满足选择学习行为的概率, 则需进行学习行为。 在进行学习行为之前, 需判断该新生波段是从精英组还是普通组产生, 以便选择不同的学习行为, 如式(7)所示。 若随机数满足选择探索行为的概率, 则需进行探索行为。 在进行探索之前, 仍需判断该新生波段是从精英组还是普通组产生, 从而为探索行为选择不同的收缩系数, 如式(5)和式(6)所示。 若行为过程中某个波长点范围越界, 则改用其边界值。

Step5: 波段更新。 候选波段的评价值需要跟三个小组中的波段评价值进行对比更新。 更新规则参阅本文1.1.4和1.2.3。

Step6: 循环Step3—Step5进行迭代更新。 迭代结束, 选出精英组中评价值最高的波段为最终筛选波段。

2.1.1 数据来源

为考察新提出的变量筛选算法对建模预测的效果, 将其应用在一组标准玉米近红外光谱数据集。 该光谱数据集引用自eigenvector网站上开源的玉米样本光谱数据集, 网址https://eigenvector.com/resources/data-sets/。 该数据集为80个玉米样品用mp5spec仪器扫描得到的光谱数据, 并用化学方法测定了其淀粉和蛋白质含量。 以间隔为2 nm在波长范围为1 100~2 498 nm上收集(700个波长点)。

2.1.2 剔除异常数据及样本划分

考虑到仪器测量光谱数据时因误差而得到异常光谱, 会影响模型性能, 故先用马氏距离剔除光谱中异常数据。 图1为各个数据样本点到数据中心点的马氏距离分布图。 马氏距离最远的两个样本(75和77号样本)作为异常点剔除, 图2为剔除异常样本之后的原始光谱图。

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 样本马氏距离分布图Fig.1 Mahalanobis distance distribution of samples

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 剔除异常样本的光谱图Fig.2 Spectra of abnormal samples removed

采用Kennard-Stone(KS)方法^[18]将剩余的78个样本点分成校正集和预测集。 划分结果为校正集样本50个, 预测集样本为28个。 校正集和预测集淀粉和蛋白质含量值统计如表1。 由表1可知, 校正集样本与预测集样本的平均值和标准差相差不大, 通过KS方法划分数据集保证了校正集样本均匀分布。

表1

Table 1

表1(Table 1)

表1 校正集和预测集中淀粉和蛋白质含量值g/100 g统计 Table 1 Statistics of starch and protein contents in correction set and prediction set 最大值 最小值 平均值 标准差 淀粉校正集 66.472 0 62.826 0 64.809 9 0.799 7 淀粉预测集 65.841 0 63.099 0 64.432 3 0.813 3 蛋白质校正集 9.595 0 7.654 0 8.569 2 0.485 2 蛋白质预测集 9.711 0 8.112 0 8.848 8 0.498 0

表1 校正集和预测集中淀粉和蛋白质含量值g/100 g统计 Table 1 Statistics of starch and protein contents in correction set and prediction set

实验使用的是一台戴尔计算机, 处理器是Intel(R) Core i5-9400, CPU主频为2.90 GHz, 操作系统为Windows10, 所有计算均在MATLAB 2016a中进行。 为了验证iTPA算法的有效性和优越性, 分别用全谱波长、 GA、 TPA以及主成分分析(principal component analysis, PCA)算法对玉米的淀粉和蛋白质含量的建模效果进行了对比。 进行GA算法时, 由于原始光谱波长点数众多, 如对波长点进行优选组合, 运算效率将会非常低, 因此将原光谱均分若干个子区间, 用GA算法进化迭代获取最大适应度值所对应的优选子区间组合。 根据基因选出的波段建立PLS模型, 计算出模型的相关系数(r)和校正集预测均方根误差(RMSEC), 适应度函数F跟iTPA算法评价值函数f(x)保持一致为

F=r1+RMSEC(8)

将原光谱在1 100~2 498 nm范围之间共700个波长点数划分35个等距区间, 即遗传编码长度为35, 每一个基因包含20个波长点数。 设定群体个数为50个, 交叉概率为0.85, 变异概率为0.1, 最大迭代数为200代。 在种群进化过程中寻找最大迭代次数内进化过程中最优适应度个体。

同样的, 采用TPA算法时将原光谱等分35个波段区间, 每个波段内含20个波长点数, 即每个成员对应20个能力因素。 按照算法经验, 精英组由15个波段组成, 普通组由20个波段组成, 该分组模式算法性能最优。 设定选择学习概率为0.35, 循环尝试不同的收缩指数, 当精英组收缩指数为20时, 算法预测效果最优, 普通组收缩指数设定为精英组的0.5倍。 迭代次数为1 000。 迭代完毕之后, 取精英组评价值最高的成员为优选成员, 优选成员中包含的能力因素就是所要筛选的波长点。 如出现重复波长点则去除即可。 因此, 采用TPA算法筛选出来的波长点数最多不超过20个。

对于iTPA算法, 同样分为35个成员, 按照算法经验, 其中精英组10个成员, 普通组10个成员, 垃圾回收组15个成员, 该分组模式算法效果最优, 其余设定条件与原TPA算法保持一致。 由于GA算法、 TPA算法、 iTPA算法都具有的随机性, 因此将以上三种算法分别运行50次求平均。

改进算法iTPA优化了TPA算法的更新方向, 促使算法很大程度上避免陷入局部最优。 图3为用TPA和iTPA筛选之后的变量分别与淀粉和蛋白质含量测试50次进行建模分析的预测均方根误差值。

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 TPA和iTPA算法运行50次(a)淀粉预测RMSEP值和 (b)蛋白质预测RMSEP值Fig.3 Prediction RMSEP values of starch (a) and protein (b) after running TPA and iTPA for 50 times

由图3可知, iTPA算法由于避免了陷入局部最优, 使得整体预测效果得到了明显的提升。 表2和表3分别为各算法测试淀粉含量和蛋白质含量的各项性能数据。

表2

Table 2

表2(Table 2)

表2 不同变量筛选方法对玉米淀粉含量的预测结果 Table 2 Prediction results of corn starch content by different variable screening methods 波长数 Rc RMSEC Rp RMSEP T/s F-PLS 700 0.905 7 0.335 7 0.871 7 0.391 4 0.028 GA-PLS 350 0.958 2 0.226 4 0.926 2 0.300 6 21.980 TPA-PLS 17.64 0.945 3 0.258 1 0.901 3 0.345 3 0.680 iTPA-PLS 17.55 0.944 1 0.260 9 0.907 9 0.334 4 0.676 PCA 20 0.899 6 0.377 0 0.792 7 0.438 3 0.467

表2 不同变量筛选方法对玉米淀粉含量的预测结果 Table 2 Prediction results of corn starch content by different variable screening methods

表3

Table 3

表3(Table 3)

表3 不同变量筛选方法对玉米蛋白质含量的预测结果 Table 3 Prediction results of corn protein content by different variable screening methods 波长数 Rc RMSEC Rp RMSEP T/s F-PLS 700 0.951 7 0.147 4 0.930 7 0.178 9 0.028 GA-PLS 346 0.979 1 0.097 8 0.973 5 0.111 4 22.120 TPA-PLS 19.64 0.975 5 0.105 6 0.969 6 0.119 3 0.678 iTPA-PLS 19.60 0.975 9 0.101 9 0.970 4 0.117 7 0.666 PCA 20 0.956 9 0.141 7 0.891 6 0.236 2 0.469

表3 不同变量筛选方法对玉米蛋白质含量的预测结果 Table 3 Prediction results of corn protein content by different variable screening methods

从表2可得, 在玉米淀粉含量的预测上, iTPA相比于全谱PLS(F-PLS), 变量个数从700个减少到均值为17.55个(50次运算求平均), 模型的校正均方根误差RMSEC从0.335 7降到0.260 9左右, 校正集预测精度提升22.3%。 模型的预测均方根误差RMSEP从0.391 4下降到0.334 4左右, 预测集预测精度提升14.6%。 相比原TPA算法有了提升, TPA算法的RMSEP为0.345 3。 GA算法总体预测效果是最佳, 但其筛选出来的平均波长点数高达350个和算法运算时间21.98 s, 远远大于iTPA算法的17.55个和0.676 s。 通过对原光谱进行PCA降维, 筛选出20个波长点(贡献率已达0.999 98), 以此与筛选出近乎相同波长数的iTPA算法作比较。 经PCA算法筛选波长之后得到的RMSEP为0.792 7, 预测效果远不如iTPA算法。 因此, iTPA算法应用在淀粉含量光谱数据集中能在保持预测能力的前提下大幅度削减波长点数, 有效地减小建模的计算量, 同时算法速度更快。

从表3可得, 对于玉米蛋白质的预测, iTPA算法预测均方根误差RMSEP为0.117 7, 预测相关系数R_p为0.970 4, 校正均方根误差RMSEC为0.101 9, 校正相关系数R_c为0.975 9, 预测效果比全谱PLS和原TPA算法都有提升, 全谱PLS的RMSEP为0.178 9, TPA算法的RMSEP为0.119 3。 GA算法总体预测效果最佳, 但其筛选出来的平均波长点数高达346个, 以及算法运算时间22.12 s, 远远大于iTPA的19.60个和0.666 s。 通过对原光谱进行PCA降维, 筛选出20个波长点(贡献率已达0.999 997)。 经PCA算法筛选波长之后得到的RMSEP为0.236 2, 预测效果远不及iTPA算法。 因此, iTPA算法应用在蛋白质含量光谱数据集中, 同样能在保持预测能力的前提下大幅度削减波长点数, 有效地减小建模的计算量, 同时算法速度更快。

图4表示iTPA算法分别在玉米淀粉和蛋白质数据集上运行50次后变量被选取的频率。

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 iTPA算法运行50次后(a)淀粉光谱变量被选取的频率和(b)蛋白质光谱变量被选取的频率Fig.4 Frequencies of spectral variable selections for starch (a) and protein (b) after running iTPA for 50 times

由图4(a)可知, 淀粉近红外光谱被iTPA算法筛选后的信息变量区域主要分布在1 540~1 546, 1 576~1 588, 1 724~1 730和1 766~1 772 nm等区域, 而这些区域与淀粉中O—H的伸缩振动一级倍频以及C—H的伸缩振动一级倍频的频率一致, 这与本次研究中淀粉的化学性质相一致。

由图4(b)可知, 蛋白质近红外光谱被iTPA算法筛选后的信息变量区域主要分布在1 920, 1 958~1 962, 2 050, 2 106~2 110, 2 180~2 182和2 242~2 244 nm等区域, 而这些区域与蛋白质中C＝O的伸缩振动一级倍频、 N—H键对称、 不对称、 伸缩振动以及AmideⅠ (Ⅱ , Ⅲ )(分别为酰胺分子中羰基与胺基不同的耦合方式)的频率一致。 同时筛选的变量也有分布在1 700 nm附近, 这些区域主要对应着C—H伸缩振动一级倍频, 这是因为蛋白质中含有C元素, 这与本次研究中蛋白质的化学性质相一致。

因此iTPA算法能够有效地消除光谱信息中的一些干扰波长点, 达到筛选出有效波长点目的。