基于BWO-SVM算法的广陈皮陈化年份高光谱鉴别模型

吕石磊, 王宏炜, 李 震*, 周 旭, 赵 静

1. 华南农业大学电子工程学院(人工智能学院), 广东 广州 510642

2. 人工智能与数字经济广东省实验室(广州), 广东 广州 510330

3. 国家柑橘产业技术体系机械化研究室, 广东 广州 510642

引 言

陈皮为芸香科柑橘属植物的干燥成熟果皮, 其中以广陈皮的质量为优。 作为我国岭南地区的道地药材, 广陈皮是陈化中药的典型代表, 被列为广东三宝之首[1]。 源于茶枝柑的广陈皮相比于普通陈皮, 具有疏肝利胆、 理气健脾等独特功效, 具有“越陈越香”的特征, 这也造成不同陈化年份的广陈皮在售价上有巨大的区别。 一些不法商家为了谋取巨大的经济利益进行年份造假的行为时有发生, 为有效避免消费者的损失, 如何应用简便、 准确、 非破坏性的手段鉴别出广陈皮的真实陈化年份具有重要意义。

国内外有很多学者都对陈皮陈化年份的鉴别进行了相关研究。 Choi等采用高效液相色谱与二极管阵列检测(HPLC-DAD)技术, 通过分析陈皮样本中的酚类化合物含量来推断陈化年份[2]; 柴利萍采用可见分光光度法测定不同年份陈皮样本的黄酮含量, 以此来构建陈化年份鉴别模型[3]。 在光谱学的应用领域, 杨方等在陈皮表面依附带有银纳米颗粒的金膜, 利用该特殊基底能够增强拉曼光谱信号的特性, 通过采集与分析陈皮样本表面的增强拉曼光谱数据来推断陈化年份[4]; 余梅等采集了不同陈化年份陈皮的近红外光谱, 并采用主成分分析(PCA)、 差异成分分析(SIMCA)等方法构建了不同年份陈皮的鉴别模型[5]。 由于陈皮成分复杂, 各光谱波段间的信息关联度较差, 因此运用光谱学知识对陈皮进行年份鉴别需要通过改进光谱数据处理方法提高鉴别准确率。 作为一种无损鉴别技术, 高光谱成像能够将样本的光谱信息与空间信息相结合, 进而在样本成分分析方面具有较好的技术优势, 已被广泛应用在地表遥感、 水质检测、 植被监测等领域[6]。

针对现阶段广陈皮陈化年份鉴别的不足和空白, 应用高光谱成像技术, 以不同陈化年份的广陈皮样本作为研究对象, 采用不同的光谱预处理算法、 特征波段选择方式和优化的鉴别模型对广陈皮的陈化年份进行鉴别, 并提出一种结合黑寡妇优化算法和支持向量机(BWO-SVM)的高光谱快速无损鉴别方法, 为广陈皮陈化年份鉴别提供了新的技术参考。

1 实验部分

1.1 材料与设备

采用的广陈皮出产于广东省江门市的新宝堂陈皮有限公司(网购), 样本均为茶枝柑的干燥成熟果皮, 样品陈化年份分别为5年、 10年、 15年、 20年, 不同年份的广陈皮在外观上很难进行人工区分。

采用的高光谱采集设备是北京卓立汉光仪器有限公司的GaiaField便携式高光谱成像系统, 如图1所示。 该系统主要由面阵CCD探测器, 高光谱镜头(分辨率696×256 Pix), 均匀光源(4个100 W卤素灯, 位于采样台四周), 光谱仪(波长范围385～1 014 nm, 共256个波段)等组成。

图1 高光谱成像系统结构

图2 5年广陈皮样本ROI取样

1.2 数据采集

由于广陈皮样本形状各异、 表面凹凸不平, 采集的高光谱图像会有较大的差异; 同时光谱采集过程中相机存在暗电流等噪声影响, 为了使采集效果达到最佳, 必须调整相机的镜头参数。 设置高光谱成像系统的相机曝光时间为27 ms, 镜头到样本的距离为42 cm, 平台移动速度为1.38 mm·s-1[7]。

高光谱需要在黑暗密闭的环境中进行采集, 使用100 W卤素灯模拟太阳光, 调整卤素灯方位角为45°且距离样本区垂直高度为40 cm。 采样前需对广陈皮样本擦拭去除表面的灰尘。 在相同条件下, 样本内侧朝向相机镜头放置于反射率接近于0的黑色绒布上, 光谱仪镜头尽量对准样本的中心位置。 共获得341个广陈皮样本的高光谱数据, 包括86个5年样本数据、 87个10年样本数据、 82个15年样本数据和86个20年样本数据。

数据分析处理软件包括: SpecVIEW、 ENVI5.3、 The Unscrambler X10.4、 Matlab 2016a、 Python3.10等。

1.3 数据处理

使用SpecVIEW软件对广陈皮的高光谱样本进行镜头校准和反射率校准, 利用自带的镜头校准修复失真和透视; 利用标准反射板的反射率校正数据, 得到更准确的样本高光谱反射图像[8]。 高光谱图像反射率校正如式(1)所示。

R=(R0-RB)/(RW-RB)

(1)

式(1)中,R为校正后高光谱图像,R0为原始图像,RB为遮挡镜头得到的黑板图像,RW为白板反射图像。

使用ENVI5.3软件对采集到的广陈皮高光谱数据进行处理。 在每个样本的四角位置和中间位置分别取点, 总计选取5个感兴趣区域(region of interest, ROI), 每个ROI选择的范围是20×20 Pix, 因此每个样本可得到2 000(400×5)条光谱数据, 将数据导入到Matlab 2016a中进行均值化处理。

1.4 高光谱数据预处理算法

高光谱成像数据受到来自广陈皮自身因素干扰(样本表面平整度、 水分含量等)和其他外界因素干扰(环境、 暗电流等), 因此采集的样本数据中通常包含着无关信息和噪声[9]。 为了降低以上干扰对数据采集的影响, 需要进行光谱预处理。

1.4.1 多项式平滑算法

多项式平滑算法(savitzky-golay smoothing, SG)是基于最小二乘原理的平滑算法, 其原理是利用中心点以及周边若干点进行最小二乘拟合, 尽可能减少噪音对有用信息的影响, 同时提高光谱数据的平滑性。

1.4.2 多元散射校正算法

多元散射校正算法(multiple scattering correction, MSC)通过理想光谱修正数据的基线偏移现象, 可以有效的消除不同散射水平的光谱差异, 从而增强光谱与数据之间的相关性, 分离散射光谱中的物理光散射信息以及化学光吸收信息, 校正样品中的光谱信息。

1.4.3 去趋势算法

去趋势算法(detrend)将光谱图象分解成静态稳定分量和非线性趋势项, 并消减非线性趋势项, 从而得到静态稳定分量, 其优势是减少相关性分析中数据畸变的概率, 防止低频成分上翘掩盖主频成分, 影响模型鉴别精度。

基于各类预处理算法的广陈皮高光谱数据如图3所示, 其中基于SG_detrend算法处理后的光谱曲线有相对明显的波峰和波谷, 因此更适合于处理广陈皮高光谱数据。

图3 基于各类预处理算法的广陈皮高光谱数据

1.5 降维算法

高光谱具有波长多、 谱带重叠等特征, 在处理高光谱信息时可能出现运算资源的浪费和运算时间过长等问题。 同时无关的波段对建立模型还会产生负面影响, 所以需要通过光谱降维去除冗余波段, 这样既提高模型的鉴别准确率, 同时也降低模型过拟合风险和存储运算压力。

1.5.1 连续投影算法

连续投影算法(successive projections algorithm, SPA)是一种使矢量空间共线性最小化的前向变量选择算法, 在降低共线信息和获取有效信息等方面有较好的效果[10]。 在全波段的光谱数据中选择出提供主要贡献的波段, 减少参与鉴别的光谱维度, 同时降低噪声和其他冗余信息的影响。

将基于SG_detrend算法处理后的广陈皮高光谱数据作为输入样本, 经SPA算法筛选特征波段后, 得到18个特征波段, 如图4所示。

图4 SPA算法筛选特征变量过程

1.5.2 竞争性自适应重加权算法

竞争性自适应重加权算法(competitive adaptive reweighted sampling, CARS)以波段变量回归系数的绝对值当参考指标, 选取数值大的变量, 筛除数值小的变量, 其中变量个数由衰减指数法(exponentially decreasing function, EDF)确定; 在未被选取的波段变量中, 采用自适应加权算法(adaptive reweighted sampling, ARS)再次选取回归系数绝对值相对较大的波段[11]。

将基于SG_detrend预处理算法的广陈皮高光谱数据作为输入样本, 设置蒙特卡洛迭代次数为50; 结果表明, 当迭代次数为8时, 筛选效果最好, 经CARS算法筛选后得到50个特征波段, 如图5所示。

图5 CARS筛选特征变量过程

1.5.3 逐步回归算法

逐步回归算法(stepwise regression, SR)通过逐个引入贡献相对大的显著自变量来增加回归方程的准确性, 同时排除非显著变量, 减少冗余信息干扰, 最终得到最优的回归模型。 由于CARS算法筛选出的特征波段数量较大, 需要进一步降维, 采用SR算法对CARS算法筛选后的广陈皮高光谱特征波段进一步筛选, 最终得到22个特征波段, 如图6所示。

图6 SR算法筛选特征变量

1.6 分类鉴别模型

分类鉴别模型是模式识别与机器学习等领域的重要组成部分, 其表征了模型输入与输出之间的特定映射关系, 能够通过对已有样本数据的学习与训练来实现对广陈皮高光谱数据的精准分类鉴别。

1.6.1 偏最小二乘判别分析模型

偏最小二乘判别分析(partial least-regression, PLS)模型通过将因变量和自变量投影到新的多维空间中, 形成基于最大化自变量相干性的线性回归模型。 PLS将原始光谱中包含的信息放置在排序在前面的若干主成分中, 所选择的主成分之间互不干扰, 这样既解决了高光谱数据中的共线性难题, 同时减少了噪声影响。

1.6.2 支持向量机模型

支持向量机(support vector machine, SVM)模型采用线性方程组得到支持向量, 基于最小化经验误差和最大化几何边缘区的方式来降低结构分类风险和过拟合风险, 进而增加分类模型的稳定性。 但是SVM模型的错误惩罚系数c和核系数g大多基于经验数据获取, 通用性较差, 因此本研究结合样本数据特征进一步优化SVM模型参数。

1.6.3 粒子群算法优化SVM模型

粒子群算法(particle swarm optimization, PSO)的灵感来自于生物体的社会行为, 利用一群粒子通过学习不断调整位置和速度来寻找最佳的个体和全局最优解。 设置PSO算法的种群规模为20, 迭代次数为100, 独立运行次数为20; 将SVM模型系数作为PSO算法的可行解域, 基于最优的适应度值更新粒子位置和速度, 采用当前的全局最优解训练SVM模型, 直至达到最大迭代次数。

1.6.4 蝗虫算法优化SVM模型

蝗虫算法(grasshopper optimization algorithm, GOA) 是一种元启发式仿生优化算法[12], 模拟蝗虫在自然界中的种群迁移和觅食行为, 具有较高的搜索效率和较快的收敛速度, 并且算法本身特殊的自适应机制能够很好地平衡全局和局部搜索过程, 具有较好的效果。 设置GOA算法的种群规模为20, 迭代次数为100, 独立运行次数为20; 基于GOA算法优化SVM模型的过程与1.6.3节一致。

1.6.5 黑寡妇优化算法优化SVM模型

黑寡妇蜘蛛对能够表现个体状态的性信息素较为敏感, 高性信息素的配偶能够提供足够的生育能力, 低性信息素的个体因倾向于吞噬同类而被种群集体排斥和遗弃。 黑寡妇优化算法(black widow optimization algorithm, BWO)是基于黑寡妇蜘蛛独特行为启发的一种新型元启发式算法[13], 具有收敛速度快、 寻优效率高等优势。 在BWO算法中, 黑寡妇蜘蛛能够根据上一次迭代过程中最优解所对应的位置和本次迭代过程中其他黑寡妇蜘蛛所在的位置来调整自身的下一个寻优位置。 基于BWO算法优化SVM模型的计算流程如图7所示, 具体步骤如下所述。

图7 BWO优化SVM流程

(1)BWO算法参数: 种群规模为20, 迭代次数为100, 独立运行次数为20, 将作为黑寡妇蜘蛛的性信息素警戒值设为0.3, 个体螺旋式移动决定阈值设为0.7; SVM模型参数: 错误惩罚系数c和核系数g的取值范围设为[0.000 1, 1 000]。

(2)采用均方根误差(root mean squared error, RMSE)作为适应度函数fitness, 如式(2)所示。

fitness=argmin(RMSESVMpredict)

(2)

对所有的黑寡妇蜘蛛对应的适应度值进行排序, 得到当前最优的SVM模型参数。

(3)计算并更新当前迭代过程中所有黑寡妇蜘蛛对应的性信息素值, 如式(3)所示。 其中fitness(i)是当前黑寡妇蜘纂的适应度, pheromone(i)是当前黑寡妇蜘蛛的性信息素值。

(3)

(4)

(5)

(6)计算并更新黑寡妇群体适应度值, 若达到算法停止条件则输出最优解, 否则重复步骤2—步骤6。

2 结果与讨论

2.1 基于PLS模型的高光谱数据预处理算法性能分析

采用PLS模型检验不同的高光谱数据预处理算法性能, 指标包括决定系数(determination coefficients,R2)、 RMSE和PLS模型鉴别准确率[14]。 其中,R2值越大表征回归效果越好, RMSE越趋近于0则表示鉴别结果越靠近实际结果。 预处理算法性能结果如表1所示。 从表中可知, 各类高光谱数据预处理算法均对PLS模型的鉴别结果有一定程度的提升作用, 其中SG_detrend预处理算法的性能最佳。

表1 预处理算法性能分析

2.2 基于降维算法的广陈皮高光谱特征波段筛选

针对基于SG_detrend预处理算法的广陈皮高光谱数据, 分别采用SPA和CARS_SR两类降维算法来筛选特征波段, 如表2所示, 其中SPA算法筛选得到18个特征波段, CARS_SR筛选得到22个特征波段。

表2 广陈皮高光谱特征波段筛选分析

2.3 广陈皮分类鉴别模型性能分析

使用蒙特卡洛算法将各年份广陈皮高光谱数据按比例3∶1随机分为训练集和测试集, 共得到训练集258个, 测试集83个。 错误惩罚系数c和核系数g是决定SVM模型性能的重要参数, 本研究分别采用了PSO算法、 GOA算法和BWO算法来优化SVM模型参数。 其中, 参数c和g的取值范围为[0.000 1, 1 000]; 各优化算法的种群规模为20, 迭代次数为100, 独立运行次数为20。 经各算法优化得到的SVM模型最优参数分别为: (PSO-SVM)c=59.931 6,g=0.010 3; (GOA-SVM)c=93.975 9,g=0.506 6; (BWO-SVM)c=80.270 1,g=0.159 5。

为综合分析不同分类鉴别模型性能, 本研究设计了包含PLS、 PSO-SVM、 GOA-SVM和BWO-SVM四种算法的广陈皮陈化年份鉴别实验。 其中, 输入样本数据均采用了SG_detrend预处理算法, 特征波段分别采用了未经降维算法处理的全波段(FULL)、 SPA特征波段集合和CARS_SR特征波段集合。 不同分类鉴别模型的实验结果如表3所示。 首先, 从降维算法的角度看, 与全波段相比, SPA算法和CARS_SR算法均能够在一定程度上剔除冗余信息, 进而提升鉴别模型的各项性能; 其次, 从分类鉴别模型的角度看, 与PLS模型相比, SVM模型能够得到更好的鉴别结果, 尤其是基于全波段的FULL-BWO-SVM模型对广陈皮年份鉴别的准确率(95.18%)高于基于CARS_SR降维算法处理后的SPA-PLS模型(91.57%); 最后, 从优化算法的角度看, 与PSO算法和GOA算法相比, BWO算法优化SVM模型的实验结果具有明显优势, 也验证了BWO算法寻优效率高、 不易于陷入局部最优等特点。 从表中可知, SG_detrend-CARS_SR-BWO-SVM模型性能最好, 对广陈皮陈化年份的鉴别准确率达到97.59%, RMSE为0.060 2,R2为0.952 9, 该误差主要来自高陈化年份的广陈皮样本数据, 验证结果如图8所示; 该模型比次优模型SG_detrend-SPA-BWO-SVM提升了1.20个百分点, 比对照模型CARS_SR-PSO-SVM提升了4.82个百分点。

表3 不同分类鉴别模型的实验结果

图8 SG_detrend-CARS_SR-BWO-SVM模型验证结果

3 结 论

采集四类陈化年份的广陈皮高光谱信息作为数据集, 结合预处理算法、 降维算法, 通过构建PLS、 PSO-SVM、 GOA-SVM、 BWO-SVM四种分类模型对广陈皮的陈化年份进行鉴别, 得出如下结论:

(1)采用SG_MSC算法和SG_detrend算法对广陈皮高光谱数据进行预处理, 其中SG_detrend算法较SG_MSC算法更能有效降低环境因素干扰, 消除光谱差异, 进而突出广陈皮自身的特征信息。

(2)采用SPA算法和CARS_SR算法对广陈皮高光谱数据进行特征提取, 其中CARS_SR算法较SPA算法更能有效增强高光谱数据和陈化年份之间的关联性, 进而降低分类鉴别模型输入变量的复杂度。

(3)采用PSO算法、 GOA算法和BWO算法优化SVM模型, 以此来得到分类鉴别模型的最优控制参数, 并与PLS模型进行对照实验; 结果表明, SVM模型较PLS模型具有更好的鉴别结果, BWO算法较其他优化算法具有更好的寻优性能。

(4)SG_detrend-CARS_SR-BWO-SVM模型对广陈皮陈化年份的鉴别准确率达到97.59%, RMSE=0.060 2,R2=0.952 9; 下一步工作将侧重研究广陈皮的不同品种、 产地等因素对陈化年份与高光谱数据之间的影响, 并且加强对高年份广陈皮的鉴别能力, 提升分类鉴别模型的泛化能力。