基于粒子群算法思想和模拟退火策略的光谱特征波长选择方法

摘要

本发明公开了一种基于粒子群算法思想和模拟退火策略的光谱特征波长选择方法，以光谱数据波长变量个数为码长进行粒子群初始化和温度初始化，使用粒子中值为“1”二进制位对应的光谱波长变量建立偏最小二乘回归模型，以交叉验证均方根误差为适应度函数，选出每个粒子的个体历史最优和整个种群的群体历史最优，计算每个粒子中值为“1”二进制位的移动速度并更新位置，对更新位置的粒子生成扰动解，最终采用交叉验证均方根误差最小的粒子作为群体历史最优位置，该粒子值为“1”二进制位的编号即为优选的光谱特征波长变量。本发明实现了粒子值为“1”二进制位的指导性寻优，融合了模拟退火策略，具有高效的特征波长选择能力。

说明书

技术领域

本发明涉及光谱分析领域，尤其涉及一种基于粒子群算法思想和模拟退火策略的光谱特征波长选择方法。

背景技术

光谱分析技术具有无损、快速、低成本等优势，通过快速获取大量、准确的光谱数据结合已经建立的定量校正模型即可实现样品理化指标的快速检测。但是在采集光谱数据时，由于采集的光谱数据量很大，光谱数据中含有的样品背景、高频噪声等无关信息很难使用预处理方法完全消除，以采集的光谱全部波长变量建模时，计算量大，波长冗余严重，不仅增加了模型的复杂程度，还严重影响了模型的预测精度。因此，有必要通过特征波长优选，有效消除光谱中不相干和共线性的波长变量对建模精度和效率的影响。

目前国内外常用的特征波长选择方法有区间偏最小二乘法、反向区间偏最小二乘法、连续投影算法、竞争性自适应重加权采样算法等，同时遗传算法、模拟退火算法、随机蛙跳算法等智能优化算法在光谱特征波长选择方面得到了有效应用。相比于其它智能优化算法，粒子群算法具有实现简单、收敛迅速、参数调整少、稳定性高等优点，在特征波长选择方面应用越来越广泛。然而用于求解波长选择问题的二进制粒子群算法虽然具有很强的搜索能力，但随着算法迭代搜索次数的增加粒子将迅速失去多样性，导致算法丧失全局搜索能力，不能收敛于全局最优解。因此，亟待对二进制粒子群算法进行改进，并研究改进后的算法在光谱特征波长选择方面的应用效果具有重要意义。

发明内容

为了克服现有二进制粒子群算法在光谱波长选择过程中容易出现早熟收敛的问题，本发明提出了一种基于粒子群算法思想和模拟退火策略的光谱特征波长选择方法，基于粒子群算法思想实现了粒子值为“1”二进制位的指导性进化，融合了模拟退火策略进行扰动解的Metropolis选择复制，充分利用算法强大搜索能力的同时有效避免早熟收敛，进而实现光谱特征波长的高效优选。

本发明采用的技术方案包括以下步骤：

(1)输入光谱数据集和化学浓度数据，以Savitzky-Golay平滑、标准正态变量变换、多元散射校正、一阶导数及其组合对光谱数据进行预处理，以校正集样本光谱数据建立全谱下的偏最小二乘回归模型，通过对比交叉验证均方根误差最小确定光谱数据的最佳预处理方法；

(2)以预处理后的光谱数据为波长选择对象，依据光谱数据特点进行粒子种群的初始化和温度初始化。随机产生N个M长的二进制序列，N为种群规模，M为预处理后的光谱波长变量个数。一个二进制序列为一个粒子，序列内所有值为“1”二进制位的组合作为该粒子所选中的特征波长。初始化每个粒子的初始整体惯性移动速度、个体历史最优适应度函数值以及群体历史最优适应度函数值；

(3)对每一代种群中的每个粒子，将其二进制序列中值为“1”二进制位对应的波长选中参与建模运算，使用K折交叉验证计算其偏最小二乘回归模型的交叉验证均方根误差，并将其作为本算法的适应度函数值；

(4)每个粒子各自将该次计算得到的适应度函数值和个体历史最优值比较，将数值较小值更新为该粒子的个体历史最优适应度函数值，将该值所对应的粒子位置记录为个体历史最优位置。将每个粒子的个体历史最优适应度函数值和群体历史最优适应度函数值进行比较，将数值最小值更新为群体历史最优适应度函数值，将该值所对应的粒子位置记录为群体历史最优位置；

(5)每个粒子根据个体历史最优位置和群体历史最优位置进行指导性进化，变换各粒子中所有值为“1”二进制位的当前位置，使粒子二进制位逐渐向具有较好适应度函数值的方向移动；

(6)对移动二进制位“1”后的每一个粒子随机选取W位进行二进制“位变异”生成扰动解，并采用Metropolis判别准则接受扰动解。若接受扰动解，需要对扰动解中的“1”变“0”位速度置零、“0”变“1”位速度重新初始化。最后，计算种群中每一个粒子的整体惯性移动速度，并执行降温操作和运行次数递增操作；

(7)重复步骤(3)～(6)，若此时已达到了算法运行的最大运行次数，则结束循环，输出此时的群体历史最优适应度函数值和群体历史最优位置，该群体历史最优位置对应粒子中值为“1”二进制位的编号即为优选的光谱特征波长变量；

作为本发明的进一步改进，所述步骤(2)中每个粒子的初始整体惯性移动速度为V

V

其中，round()函数为四舍五入取整函数，sign()函数为归一函数，即正数归为1，负数归为-1；r

作为本发明的进一步改进，所述步骤(5)中粒子的移动遵循以下规律：

①粒子的移动即为粒子中的每一个值为“1”二进制位的移动；

②第i代进化时种群中某个粒子的第j个值为“1”二进制位的移动速度为

其中，

③关于V

V

V

其中，V

作为本发明的进一步改进，所述步骤(6)中对扰动解“0”变“1”位速度重新初始化，并计算粒子的整体惯性移动速度。设第i代进化时扰动解第j个“0”变“1”位的初始化速度为

其中，round()函数为四舍五入取整函数，sign()函数为归一函数，即正数归为1，负数归为-1；r

与现有技术相比，本发明的有益效果在于：

相比于传统二进制粒子群算法在更新粒子的速度位置时需要对粒子内所有二进制位进行操作(判定该位是否为“1”)，本发明将每个粒子中的值为“1”二进制位看作一个独立个体，仅对值为“1”二进制位进行指导性移位操作，基于粒子群算法思想使其朝着适应度更好的区段移动，算法运行更简便、高效。通过融合模拟退火策略对移位后的粒子二进制位进行扰动，有效解决了二进制粒子群算法迭代搜索后期因种群多样性丧失易导致早熟收敛的问题。本发明通过构建粒子群算法思想与模拟退火策略有机结合的光谱特征波长选择方法，充分发挥了粒子群算法的高效搜索能力，解决了二进制粒子群算法早熟收敛和模拟退火算法搜索效率低的问题，能够有效实现光谱数据的特征波长优选。

附图说明

图1为本发明的流程示意图；

图2为实施例的土壤样本原始光谱图；

图3为个体历史最优位置示意图；

图4为个体历史最优位置左侧无“1”情况示意图；

图5为个体历史最优位置“1”重叠情况示意图；

图6为群体历史最优位置示意图；

图7为粒子二进制位“1”的位移示意图；

图8为不同波长选择方法进化过程对比图。

具体实施方式

下面结合附图和具体实施例对本发明作进一步详细说明。

本发明实施例取自网站Quality&Technology土壤开源数据集，由可见近红外光谱数据和化学指标数据构成，共计108个样本。样本光谱波长范围为400～2500nm，采样间隔为2nm，波长变量个数为1050个，样品光谱数据如图2所示。本发明采用土壤有机质含量为因变量进行特征波长选择和光谱建模分析以证明本发明方法的有效性。具体实施步骤如下：

S1：输入光谱数据集和化学浓度数据，采用随机选择法按3:1的比例划分校正集和验证集，得到81个校正集样本和27个验证集样本。以Savitzky-Golay平滑(SG)、标准正态变量变换(SNV)、多元散射校正(MSC)、一阶导数(FD)其组合对光谱数据进行预处理，以校正集样本光谱数据建立全谱下的偏最小二乘回归模型，通过对比10折交叉验证均方根误差最小确定光谱数据的最佳预处理方法为Savitzky-Golay平滑。本实施例光谱数据不同预处理方法对比结果如表1所示。

表1不同预处理方法对比结果

S2：以预处理后的光谱数据为对象，进行光谱波长变量优选，初始化粒子种群和各项参数。随机产生N个M长的二进制序列作为初始种群，其中N为粒子种群规模，取待优化波长变量个数的五分之一上取整数，M为预处理后的光谱波长变量个数。本实施例设定N为210，M为1050。依据公式(1)将每个粒子的初始整体惯性移动速度初始化为[-10,-1]∪[1,10]区间内的随机整数：

V

其中，round()函数为四舍五入取整函数，sign()函数为归一函数，即正数归为1，负数归为-1；r

以初始种群每个粒子所有值为“1”二进制位的组合作为选中的特征波长，使用10折交叉验证计算偏最小二乘回归模型的交叉验证均方根误差作为初始种群每个粒子的适应度函数，按t

S3：对每一代种群中的每个粒子，将其二进制序列中值为“1”二进制位对应的波长选中参与建模运算，使用10折交叉验证计算其偏最小二乘回归模型的交叉验证均方根误差，并将其作为当前粒子的适应度函数值。

S4：将当前粒子的适应度函数值与其个体历史最优值和群体历史最优值比较，取较小者为新的最优值，并记录此时该粒子的个体历史最优位置和种群的群体历史最优位置。

S5：每个粒子根据粒子的当前位置与速度和个体历史最优位置与群体历史最优位置的差距，通过粒子内每个值为“1”二进制位的移动来计算并更新各粒子的速度和位置。具体步骤如下：

第一步，计算粒子趋向于个体历史最优的速度时，对该粒子内每个“1”，计算其距离个体历史最优位置中该“1”的位编号左右最近的“1”的距离R

在搜寻最接近的“1”的时候，若其中某一个方向直至编码尽头也没有“1”，则第二步在计算该“1”向该方向的吸引速度时，设定该方向的吸引速度为0，如图4所示。若该粒子的“1”处于个体历史最优位置编码中的位数正好有“1”，则在第三步中计算该“1”的速度时，设定该速度为0，如图5所示；

第二步，根据R

第三步，根据已计算出的V

V

第四步，计算粒子趋向于群体历史最优的速度时，对该粒子内每个“1”，计算其距离群体历史最优位置中该“1”的位编号左右最近的“1”的距离R

在搜寻最接近的“1”的时候，若其中某一个方向直至编码尽头也没有“1”或者若该粒子的“1”处于群体历史最优位置编码中的位数正好有“1”，则处理办法同理第一步，分别为在第五步设定该“1”对该方向的吸引速度为0和在第六步设定该“1”的速度为0；

第五步，根据R

第六步，根据已计算出的V

V

第七步，根据公式(8)计算出该“1”的临时移动速度

第八步，对该二进制位“1”的移动速度

第九步，根据计算出的

第十步，对该粒子内的每个值为“1”二进制位重复第一至第九步，直至所有“1”移动完毕，得到该粒子的新位置编码，并将新位置编码作为该粒子更新位置后的位置编码。

S6：对移动二进制位“1”后的每一个粒子，构建邻域扰动解，基于Metropolis判别准则接受扰动解。在完成扰动解二进制位速度更新后，计算种群中每一个粒子的整体惯性移动速度，并执行降温操作和运行次数递增操作。具体步骤如下：

第一步，对二进制位“1”移位后的粒子X随机选取W位进行二进制“位变异”，即对随机选取的W位二进制位执行原“0”变“1”、原“1”变“0”操作，进而生成扰动解X′，其中W取码长的二十分之一上取整，本实施例W取53；

第二步，对X的扰动解X′，令Δf＝f(X)-f(X′)，若Δf≥0,则接受邻域解X′到下一代种群；若Δf＜0，则生成随机数r∈[0，1]，当r＜exp(Δf/t

第三步，若接受扰动解X′，需要更新扰动解中变异位的速度，将扰动解中的“1”变“0”位速度置零，根据公式(10)计算公式将“0”变“1”位速度重新初始化，其他值为“1”的二进制位移动速度保持不变；

第四步，更新完所有粒子的位置并执行Metropolis选择复制后，按根据公式(11)计算种群中每一个粒子的整体惯性移动速度，并按公式t

S7：重复S3～S6进行历史最优位置的迭代更新。若达到了算法运行的最大运行次数，则结束循环，输出此时的群体历史最优适应度函数值以及群体历史最优位置作为该算法优选的光谱特征波长。

为了评测本发明提出方法的寻优性能，将其与标准二进制粒子群算法、模拟退火算法、标准二进制粒子群算法联合模拟退火算法的波长选择性能进行对比。为了消除智能优化算法光谱波长选择的随机性对评测结果的影响，每种算法均运行5次，取5次优选结果的适应度函数平均值作为评价参数对算法的寻优性能进行评测。图8为本发明提出的基于粒子群算法思想和模拟退火策略的光谱特征波长选择方法(本方法)与标准二进制粒子群算法(BPSO)、模拟退火算法(SA)、标准二进制粒子群算法联合模拟退火算法(BPSO-SA)的波长选择进化过程对比图。图8中横坐标为波长选择算法运行次数，纵坐标为交叉验证均方根误差，实线为本发明提出的基于粒子群算法思想和模拟退火策略的光谱特征波长选择方法进化过程中得到的每一代群体历史最优适应度函数的平均值，虚线为标准二进制粒子群算法波长选择过程中得到的每一代群体历史最优适应度函数值的平均值，点划线为模拟退火算法波长选择过程中得到的每一代群体历史最优适应度函数的平均值，点虚线为标准二进制粒子群算法联合模拟退火算法波长选择过程中得到的每一代群体历史最优适应度函数的平均值。从图8可以看出，本发明提出的基于粒子群算法思想和模拟退火策略的光谱特征波长选择方法要优于标准二进制粒子群算法、模拟退火算法和标准二进制粒子群算法联合模拟退火算法特征波长选择的寻优性能。

为了评测本发明提出方法的建模精度和效率，将其与标准二进制粒子群算法、模拟退火算法、标准二进制粒子群算法联合模拟退火算法优选特征波长的建模性能进行对比，以5次运行的波长优选结果分别建立偏最小二乘回归模型，使用校正决定系数

表2不同方法建立的偏最小二乘回归模型性能指标

从表2可以看出，本方法得到的校正模型预测精度优于全谱和标准二进制粒子群算法、模拟退火算法、标准二进制粒子群算法联合模拟退火算法优选波长的建模精度，其验证集的

本实施例公布的是较佳的实施例，但并不局限于此，本领域的普通技术人员，极易根据上述实施例，领会本发明的精神，并做出不同的引申和变化，但只要不脱离本发明的精神，都在本发明的保护范围内。