基于可见_近红外高光谱图像的药品快速鉴别研究

基于可见_近红外高光谱图像的药品快速鉴别研究

                                            四川双利合谱科技有限公司

0 前言

化橘红又名化皮、化州橘红，为芸香科植物化州柚的未成熟果实的外层果皮。前者习称“毛橘红”，后者习称“光七爪”、“光五爪”。化橘红不仅具有治咳化啖、健胃**、醒酒功能，而还是人体美容的*佳原料，有广阔的市场前景。研究表明，挥发油，黄酮类化合物，多糖以及香豆素类化合物等是化橘红的主要有效成分。吴宋夏证实了化橘红中的黄酮苷具有镇痛、**、解热的作用；香豆素类化合物具有抗氧化、**抗病毒之功效。不同的品种有效成分的含量不一样，功效不一样，并且在价格上也相差较大，以**皮的效果*佳。因此市场上存在许多用化橘红的**果、伪品果、伪品皮冒充**皮，损害了消费者利益，也冲击了种植优良品种的农民们的利益。

目前对**皮常用的鉴别方法主要有性状鉴定、显微鉴定、高效液相色谱法。这些方法虽然各有优势，但是存在不同程度上的主观性强、需要预处理、实验过程复杂等缺点，不能满足市场快速、可靠检测的需要。本研究分别利用400-1000nm、1000-2500nm两款成像高光谱相机获取**皮、伪品皮、**果、伪品果四种化橘红成分的高光谱信息，利用光谱指数(DVI、NDVI)、偏*小二乘判别分析(PLS-DA)和极限学习法(ELM)三种方法分别构建四种不同成分的鉴别模型，并用独立样本数据对不同的模型进行验证。

基于可见_近红外高光谱图像的药品快速鉴别研究

1 材料与方法

1.1 材料

试验用的化橘红四种不同成分**皮、伪品皮、**果、伪品果由广东省食品药品职业技术学校提供，其中**皮样本32个、**果10个、伪品果11个，伪品皮7个。样本经粉碎均匀后，各取5g放置于培养皿上，标号，用于高光谱相机的光谱采集。

1.2 高光谱图像采集

高光谱图像数据采集采用四川双利合谱科技有限公司的 GaiaSorter高光谱分选仪系统（V10E、N25E-SWIR）。该系统主要由高光谱成像仪、面阵列相机、卤素灯光源、暗箱、计算机组成，如图1。实验仪器参数设置如表1。

图1 GaiaSorter 高光谱分选仪

表1   GaiaSorter 高光谱分选仪系统参数

在进行高光谱图像采集时，需要设置相机曝光时间，平台移动速度以及物镜之间的距离。这 3 个参数相互影响，图像调节的目的是使采集的图像大小合适，清晰，不变形失真。经过反复尝试，物镜高度设置为 31 cm，曝光时间设置为10ms，平台移动速度分别设置为 6.0 mm/s（400-1000 nm）、16mm/s（1000-2500 nm）。图像采集软件采用四川双利合谱科技有限公司提供的高光谱成像系统采集软件完成。图像处理采用 ENVI5.3 软件进行处理。在进行图像处理之前，先要对采集的光谱图像进行图像校正，图像校正公式如下：

式中，R_ref 是校正过的图像，DN_raw 是原始图像，DN_white为白板校正图像，DN_dark 是黑板校正图像。高光谱图像的光谱与图像之间有着对应的关系，据此，在**皮、伪品皮、**果、伪品果四种样本上选取 25×25 像素点的感兴趣区域，以感兴趣区域所有像素的光谱平均值作为该样本的平均光谱。

1.3  光谱噪声去除

试验得到光谱含有由仪器和试验条件等引起的噪声，对这些噪声的处理有助于减少噪声对光谱分析的影响，突出光谱的有效信息。SG 平滑算法可以有效消减光谱数据中的随机噪声，消噪效果受平滑点数的影响，本文中选择SG 二次多项式 7 点平滑对光谱数据进行处理（何勇，2013）。

1.4  特征波长选择

光谱信息之间存在大量的冗余和共线性信息特征，对光谱有效信息的提取产生了较大的干扰，且大量光谱数据造成模型复杂、计算量大的问题。本文采用连续投影算法（successive projections algorithm，SPA）进行特征波长的选择，以减少信息冗余和共线性的影响，简化模型，减少计算量。

SPA 是一种特征变量前向选择算法，在光谱特征波长中取得了广泛的应用。本文采用 SPA 算法对去噪处理后的光谱进行特征波长选择。进行SPA 计算时，以建模集的光谱数据和类别赋值为输入，设置选择特征波长数的范围为 5～30。

1.5 光谱指数

光谱指数的产生来源于植被指数，植被指数是指利用卫星不同波段探测数据组合而成的，能反映植物生长状况的指数。植物叶面在可见光红光波段有很强的吸收特性，在近红外波段有很强的反射特性，这是植被遥感监测利用卫星不同波段探测数据组合而成的，能反映植物生长状况的指数。植物叶面在可见光红光波段有很强的吸收特性，在近红外波段有很强的反射特性，这是植被遥感监测的物理基础，通过这两个波段测值的不同组合可得到不同的植被指数。光谱指数是通过任意两波段组合或三波段组合成各种光谱指数，如NDVI、DVI等，探寻*佳的波段组合用于各个领域的模型构建等。

归一化植被指数(Normalized difference vegetable index)

差值植被指数(Difference vegetable index)

其中，λ1和λ2代表任意波长的反射率，波段范围为400-1000 nm与1000-2500 nm。

1.6 判别分析方法

偏*小二乘法判别分析(PLS-DA ,Partial least squares discrimination analysis)是一种用于判别分析的多变量统计分析方法。判别分析是一种根据观察或测量到的若干变量值，来判断研究对象如何分类的常用统计分析方法。其原理是对不同处理样本（如观测样本、对照样本）的特性分别进行训练，产生训练集，并检验训练集的可信度（Luna et al., 2013）。本文分别基于全光谱、特别波长光谱建立 PLS-DA 判别分析模型，通过建立光谱数据与类别特征之间的回归模型，进行判别分析。

1.7 极限学习机

极限学习机（extreme learning machine）ELM是一种简单易用、有效的单隐层前馈神经网络SLFNs学习算法。2004年由南洋理工大学黄广斌副教授提出。传统的神经网络学习算法（如BP算法）需要人为设置大量的网络训练参数，并且很容易产生局部*优解。极限学习机只需要设置网络的隐层节点个数，在算法执行过程中不需要调整网络的输入权值以及隐元的偏置，并且产生**的*优解，因此具有学习速度快且泛化性能好的优点。本文中隐含层神经元个数从 1 到 520（288）以步长 1 进行寻优，以*小训练误差下的神经元个数为 ELM 模型隐含层神经元个数。

1.8  评价指标

    回归模型得到的样本的预测值不是整数，需要设置阈值以判断样本的归

属。本文中阈值设置为 0.5 ，预测值小数点大于或等于0.5则加1归整，小于0.5则减1归整。总体识别精度是指正确识别的个数除以总数，**皮识别精度是指**皮正确识别的个数除以**皮的总数，**皮识别错误率指数被错误分为**皮的个数除以**皮的总数。

基于可见_近红外高光谱图像的药品快速鉴别研究

2  结果与分析

2.1  化橘红不同成分的原始光谱曲线

本试验采用的V10E 相机获取的是400-1000 nm波长范围共520个波段的可见/近红外光谱数据，N25E-SWIR相机获取的是1000-2500 nm波长范围共288个波段的近红外光谱数据，全部样本的原始光谱图如图 2所示，**皮、伪品皮、**果、伪品果的光谱比较图如图3所示。

从图2和图3可以看到，总体而言，无论是400-1000 nm或1000-2500 nm波长范围内，**皮的光谱反射率值低于其他三种成分的光谱曲线，从曲线变化趋势来看四种不同成分并没有十分明显的差异。本研究按照Kennard-Stone 算法将样本分成建模集和预测集，其中建模集 38 个样本，预测集32个样本。**皮、伪品皮、**果、伪品果分别赋值为 1、2、3、4（表 2），不同化橘红成分建模集和预测集样本的划分如表 2所示。

图2 全部化橘红样本的原始反射光谱图

图3 化橘红不同成分反射光谱曲线图

表2  化橘红不同成分类别赋值与建模集合检验集样本划分

2.2  主成分分析定性分析

对化橘红的四种不同成分的光谱数据进行主成分分析（principal component analysis，PCA） 可知，400-1000 nm范围内，**主成分（principle component 1, PC1）的贡献率为88.36%，PC2 的贡献率为7.24%，PC1 和 PC2 累计贡献率为95.6%，PC1 和 PC2 能够解释绝大部分的变量；在1000-2500 nm范围内，**主成分（principle component 1, PC1）的贡献率为93.27%，PC2 的贡献率为3.72%，PC1 和 PC2 累计贡献率为97.0%，PC1 和 PC2 也能够解释绝大部分的变量；图4分别为400-1000 nm与1000-2500 nm范围内，**主成分与**主成分的散点分布图。

图4  400-1000nm（左）、1000-2500nm（右）**主成分与**主成分的得分散点分布图

图5 主成分组合识别**皮（左：400-1000 nm，右：1000-2500 nm）

由图 4 可知，在400-1000 nm与1000-2500 nm范围内，各不同成分之间较难直接区分出来，得分图中均有重合的地方。从图5亦可以看出无论是400-1000 nm或1000-2500 nm，主成分的成分组合均未能完全从**皮、伪品皮、**果、伪品果四种样本中识别出**皮，从图5可知，识别出的**皮亦有假的**皮。因此需要对光谱数据进行进一步的分析和处理，以鉴别化橘红不同成分。

2.3  光谱指数

本研究使用的成像高光谱仪V10E范围是400-1000nm，N25E-SWIR范围是1000-2500 nm，通过不同波段的重新组合，形成不同形式的光谱指数(归一化光谱指数，NDVI；差值光谱指数，DSI)。

运用Matlab软件编程，将两两组合的所有波段构建的光谱指数与各类别赋值，计算相应的决定系数(R²)，绘制决定系数图。图中颜色从蓝色到红色变化，图像颜色越红，表示决定系数越大，图像颜色越蓝，表示决定系数越小。图6为化橘红各成分分类赋值与NDVI、DVI决定系数（R2）二维图。从图6可知，在400-1000 nm范围内， DVI模型决定系数*高的两波段在绿光范围内，分别是538.91 nm和543.75nm；NDVI模型**系数*高的两波段989.76 nm和670.39 nm。在1000-2500 nm范围内，NDVI和DVI两个光谱指数模型决定系数*高的两波段组合均为是1820.72 nm和1787.12 nm。综合400-1000 nm与1000-2500 nm

图 6化橘红各成分分类赋值与NDVI、DVI决定系数（R²）二维图

两组数据构建NDVI和DVI模型时，NDVI和DVI两个光谱指数模型决定系数*高的两波段组合与1000-2500 nm范围内相同。比较400-1000 nm和1000-2500 nm的决定系数图可知，在1000-2500 nm范围内，NDVI与DVI构建的模型决定系数较高。

根据图6可知，在400-1000 nm范围内，DVI_{(538.91/543.75)}与NDVI_{(989.76/670.39)}与类别赋值构建的模型决定系数*高，分别为0.643和0.640。在1000-2500 nm范围内，DVI_{(1820.72 /1787.12)}与NDVI_{(1820.72 /1787.12)}与类别赋值构建的模型决定系数*高，分别为0.861和0.834。图7分别为DVI_{(538.91/543.75)}、NDVI_{(989.76/670.39)} 、DVI_{(1820.72 /1787.12)}、NDVI_{(1820.72 /1787.12)}与类别赋值的线性拟合散点分布图，从图中可知，在1000-2500 nm范围内构建的NDVI、DVI光谱指数与类别赋值拟合度*高，且变化曲线显著。 

                图7  DVI、NDVI光谱指数与类别赋值的散点分布图

运用独立的数据，分别对图7中DVI_{(538.91/543.75)}、NDVI_{(989.76/670.39)} 、DVI_{(1820.72 /1787.12)}、NDVI_{(1820.72 /1787.12)}与类别赋值构建的模型进行检验，检验结果如图8和表3所示。根据查表可知，1-15与28-32为**皮，15-25号为16-20**果、21-24为伪品果，25-27为伪品皮。从图8可知，**皮、**果、伪品果、伪品皮之间均存在不同程度的错误识别。表3基于DVI、NDVI光谱指数模型检验化橘红样本的精度评价表。表中分别统计了总体识别精度、**皮识别精度和**皮识别错误率。从表3中可知，总体识别精度、**皮识别精度*高的均为DVI_{(1820.72 /1787.12)}构建的模型，分别是66%和75%，**皮识别错误率*低的则为DVI_{(1820.72 /1787.12)}、NDVI_{(1820.72 /1787.12)} 构建的模型，均为5%。

图 8 基于NDVI、DVI检验集预测结果对比图

表2 基于光谱指数模型检验化橘红样本的精度评价

2.4  特征波长选择

本文以建模集样本的光谱数据和类别赋值为输入，利用连续投影算法SPA选择特征波长。选出的特征波长的个数如表 3所示。从表 3 可知，400-1000 nm范围内所选择的特征波段为15个，1000-2500 nm范围内所选择的特征波段较少，为5个。

表 3在400-1000 nm与1000-2500 nm范围内SPA 算法选择的特征波长个数

2.5  偏*小二乘判别分析

分别将建模集合检验集的全波段光谱及特征波段光谱作为偏*小二乘法判别分析输入变量，从而获取预测样本类别的赋值。图9为检验集与预测值的类别赋值图，并针对图9进行总体识别精度、**皮识别精度、 **皮识别错误率的统计，如表4所示。从表4中可知，总体识别精度、**皮识别精度*高的均为在1000-2500 nm范围内的PLS-DA构建的模型，分别是78%和90%，**皮识别错误率*低的则为1000-2500 nm范围内的PLS-DA、PLS-DA 构建的模型，均为5%。从表4可以看出，无论是400-1000 nm或1000-2500 nm，基于全波段的PLS-DA模型总体识别率和**皮识别率均高于基于特征波段的PLS-DA模型，而**皮的错误识别率，无论是400-1000 nm或1000-2500 nm范围内，基于全波段的PLS-DA模型与基于特征波长的PLS-DA模型的错误识别率相同，400- 1000 nm范围错误识别率均为10%，1000-2500范围错误识别率则为5%。

图 9 基于PLS-DA检验集预测结果对比图

表4基于PLS-DA模型检验化橘红样本的精度评价

基于可见_近红外高光谱图像的药品快速鉴别研究

2.6  极限学习机

分别将建模集合检验集的全波段光谱及特征波长光谱作为极限学习机的输入变量，从而获取预测样本类别的赋值。图10为检验集与预测值的类别赋值图，并针对图10进行总体识别精度、**皮识别精度、 **皮识别错误率的统计，如表5所示。从表5中可知，总体识别精度、**皮识别精度*高的均为在1000-2500 nm范围内的ELM与ELM-SPA构建的模型，分别是84%和95%，**皮识别错误率*低的则为1000-2500 nm范围内的ELM与ELM-SPA 构建的模型，均为5%。从表5可知，在400-1000 nm范围内，基于特征波段光谱的ELM模型总体识别率与**皮识别率均高于基于全波段的ELM模型，对于**皮的识别错误率，基于全波段与基于特征波段的错误识别率相同；在1000-2500 nm范围内，无论是基于全波段光谱或基于特征波段的ELM模型，其总体识别率、**皮识别率、**皮识别错误率均相同，分别为84%、95%和5%。 

图 10 基于ELM检验集预测结果对比图

表4 基于ELM模型检验化橘红样本的精度评价

2.7  光谱指数模型，PLS-DA 模型，和 ELM 模型的比较

综合对比光谱指数模型，PLS-DA 模型，和 ELM 模型的识别效果可知，无论是光谱指数模型，PLS-DA 模型或ELM 模型，基于1000-2500 nm范围内构建的模型，其预测值的总体识别率、**皮识别率均高于400-1000 nm范围内的模型，且**皮的识别错误率也低于400-1000nm范围内的模型。在光谱指数模型、PLS-DA 模型和 ELM 模型的模型中，ELM 模型的识别准确性*高，其次是PLS-DA模型，*后是光谱指数模型。基于特征波段光谱的PLS-DA模型其识别准确性低于基于全波段光谱的PLS-DA的模型，但是基于特征波段光谱的ELM模型在400-1000 nm范围内，其识别准确性高于基于全波段光谱的ELM模型，在1000-25000 nm范围内，其识别准确性与基于全波段光谱的ELM模型相同。

基于可见_近红外高光谱图像的药品快速鉴别研究

3  结论与讨论

基于V10E与N25E-SWIR两款成像高光谱相机，分别获取**皮、**果、伪品皮、伪品果四种化橘红成分400-1000 nm与1000-2500 nm范围的光谱反射率，采用 SG 平滑算法对提取出的光谱数据进行去噪处理，同时采用 SPA 算法对去噪后的光谱提取特征波长，并分别基于全波段光谱、特征波段光谱建立 PLS-DA 判别模型和 ELM 模型，同时采用全波段循环，探寻*佳的NDVI、DVI两个光谱指数构建判别模型，用于鉴别**皮、**果、伪品皮、伪品果，取得了比较好的识别效果。基于特征波段光谱与全波段光谱建立的 ELM 模型取得了*佳效果，总体识别精度、**皮识别精度、**皮识别错误率分别为84%、95%和5%。在实际运用中，考虑到识别时间与成分，基于SPA算法提取的特征波段构建的ELM模型效果*佳。本论文研究结果为高光谱成像技术在药品真伪等检测中的应用提供了可行性。