基于高光谱成像技术的宁夏枸杞产地溯源鉴别

袁伟东，姜洪喆,2，杨诗雨，张 聪，周 禹，周宏平,2,*

（1.南京林业大学机械电子工程学院，江苏 南京 210037；2.南京林业大学 林业资源高效加工利用协同创新中心，江苏 南京 210037）

宁夏枸杞（Lycium barbarumL.）是茄科枸杞属的多年生灌木，营养丰富且具有药用价值，被称为“药食同源”水果[1]。近代药理研究表明，枸杞子具有延缓衰老、降血糖、降血脂、降血压、免疫调节、益精明目、滋阴补肾等多方面功效[2-3]。我国枸杞种植面积广泛，主要分布在宁夏、内蒙古、新疆、甘肃等地，其中宁夏枸杞是唯一被收录于《中华人民共和国药典》的枸杞品种[4]。随着人们生活水平的不断改善以及保健意识的增强，宁夏枸杞因具有肉质饱满、甘甜、色艳、含糖量低等优点受到消费者的喜爱。但品质优良的宁夏枸杞产量有限，导致一部分不良商贩滥用其他产区的枸杞冒充宁夏优质枸杞以假乱真，严重扰乱宁夏枸杞交易市场。

传统的枸杞产地溯源检测方法包括高效液相色谱法[5]、碳同位素分析[6]和电子鼻[7]，而这些技术成本高昂、耗时且具有破坏性。因此，迫切需要一种经济、高效、快速且无损的检测方法[8]。近红外光谱作为一种快速、无损检测食品化学质量属性的技术已被广泛使用[9-10]。然而，近红外光谱是一种点测量技术，在某些情况下由于样品的异质结构，获得的光谱信息可能不足或不能代表整个样品。高光谱成像技术将光谱和成像技术相融合，同时获取光谱和空间信息的特殊优势备受食品行业研究人员的青睐[11-12]。目前，高光谱成像技术已广泛应用于农产品的产地鉴别[13]、掺假[14]、含油率[15]等多方面的检测。然而高光谱数据包含大量的共线性和冗余信息，导致模型将面临更多的计算负担和所开发模型的不稳定和复杂性。特征波长的选择对于便携式仪器的开发具有重要意义，特征波长选择有利于挑选出共线性最小、冗余最少且包含主要信息的变量。选择方法通常是从原始波长中选择少数关键波长，以便提高模型的运行效率和可解释性，去除非信息变量波段，产生更简单、更稳健的校正模型。

在本研究中，采用400～1 000 nm波段范围内的高光谱成像系统结合化学计量学方法对不同产地枸杞进行检测鉴别。应用多种特征波长选择方法选取特征波长建立判别模型并对比结果，为消除干扰信号和光散射的影响还采用了多种光谱预处理方法，旨在提供一种能够快速、无损鉴别枸杞产地的检测方法。

1 材料与方法

1.1 材料

优质枸杞样品（宁杞1号）分别购于宁夏银川、甘肃靖远、内蒙古巴彦淖尔、青海西宁和新疆乌鲁木齐5 个具有代表性的生产基地。

1.2 仪器与设备

用于测试的高光谱成像检测系统由暗箱、升降平台、成像光谱仪（SOC710VP）、4 盏75 W卤钨灯和一台计算机组成，如图1所示。光谱范围374.98～1 038.79 nm，分辨率为4.68 nm，在256 个波段上进行成像。

图1 高光谱成像系统示意图Fig.1 Schematic diagram of hyperspectral imaging system

1.3 方法

1.3.1 样品准备

挑选大小、颜色均匀的枸杞，并剔除表面有明显缺陷、有虫眼、干枯的劣质枸杞。为了减少温度对样品鉴别的影响，所有的枸杞样品在进行成像前均贮存于4 ℃的冰箱中，贮存时间为12 h。每种产地枸杞均采集450 粒样品，5 种产地共收集了2 250 粒枸杞样品高光谱图像。每种产地枸杞随机选取300 个样品分配进入训练集，其余150 个样品进入预测集。

1.3.2 高光谱图像采集与校准

高光谱图像采集软件为Hyper Scanner_2.0.127，该系统由北京安洲科技有限公司提供。经过多次调试，样品表面到高光谱镜头的距离为50 cm，积分时间设定为30 ms。为了捕捉到清晰的枸杞形状，图像分辨率设置为696×520像素，采集结束后得到一张696×520×256像素的三维图像。

为消除光照不均匀和暗电流的影响，对原始图像Ir进行反射率校正，根据式（1）计算Ic：

式中：Ic为校正后的高光谱图像；Ir为原始高光谱图像；Istd为校准后标准板的反射率；Is为标准板的反射率（约30%的反射率）。

1.3.3 反射光谱数据的提取和预处理

感兴趣区域（region of interest，ROI）光谱信息提取的主要步骤如图2所示。为了选择单个枸杞样品作为ROI，首先对原始高光谱图像进行分水岭分割、腐蚀、膨胀处理（图2a），二值化生成掩膜并分割得到单个枸杞的原始光谱图像（图2b）。然后选取图像中反射率最高和最低的两个波段（885.89 nm和417.49 nm）进行波段运算，从而得到枸杞样品与背景特征具有高对比度的灰度图像（图2c）。在生成的图像中，通过将反射阈值设置为0.15得到二值化的掩膜图像（图2d），应用掩膜图像将背景去除得到仅有枸杞的高光谱图像（图2e）。提取ROI的平均光谱作为对应样品的反射光谱（图2f）。由于336.2～400 nm和1 000～1 038.79 nm波段上的光谱强度较低，其数据波形抖动较大，包含很多噪声信号，因此在后续分析中只保留400～1 000 nm之间的光谱数据，共计231 个波段。

图2 高光谱信息提取过程Fig.2 Information extraction process from hyperspectral image

光谱数据中除了包含样品的组分信息外，还夹杂着环境噪声、杂散光、样品背景等不相关信息。因此为了减小或消除这些无用信息的影响，提高检测模型的精度和稳定性，建模前需要对原始光谱进行标准正态变量（standard normal variate，SNV）、SNV结合去趋势、归一化反射光谱（normalized reflectance spectrum，NR）、一阶导数和二阶导数等预处理。其中，一阶导数和二阶导数是基于二阶多项式拟合和五点移动窗口的Savitzky-Golay平滑实现。

1.3.4 特征波长选择方法

高光谱数据中包含大量的冗余信息，需要提取出最小冗余新变量消除数据间的共线性和重叠问题，提高模型的运行效率和模型的可解释性。

连续投影算法（successive projections algorithm，SPA）是一种前向变量循环选择方法。SPA在向量空间中使用简单的投影操作，以选择具有最小相关性的变量子集[16]。从一个波长开始，每次循环都计算它在未选中波长上的投影，最终将投影向量最大的波长引入循环后的波长组合。在本研究中，变量数的范围设置为10～50，最小均方根误差（root mean square error，RMSE）所对应的数量即为最终的选择结果。

竞争性自适应重加权算法（competitive adaptive reweighted sampling，CARS）是一种将蒙特卡罗采样与偏最小二乘回归系数相结合的特征变量选择方法，采用了“适者生存”的理论。其原理是通过自适应重加权采样方法和指数衰减函数去除偏最小二乘回归系数中绝对值最小的波长点，得到多个变量子集后，选择交叉验证均方根误差最小的变量子集作为特征变量[17]。

粒子群优化算法（particle swarm optimization，PSO）是一种基于种群的全局智能优化算法。在PSO中，每个优化问题的解都被抽象作为多维空间一个没有质量和体积的粒子。所有粒子都有一个评价函数决定的适应值（以F＝R2作为适应度函数，其中R2表示决定系数），粒子们通过相互作用发现复杂搜索空间中的最优区域[15]。本研究中，粒子种群数量大小设置为20，迭代次数设置为1 000。

迭代保留信息变量算法（iteratively retaining informative variables，IRIV）是一种基于二进制矩阵重排过滤器提出的特征波长选择算法，将所有变量分为强信息变量、弱信息变量、无信息变量和干扰变量4 类[18]。通过多次迭代的方式，对于出现比例较多的变量赋予较高的权重，并去除无信息变量和干扰变量灯对模型无关的变量，最后保留有用的信息变量，经过反向消除获得特征变量。随着变量选择方法的不断发展，变量方法联用组合使用也逐渐受到重视。梁琨等[19]通过CARS-IRIV筛选高光谱数据特征变量建立库尔勒香梨可溶性固形物含量预测模型，在提高预测精度的同时简化了模型的运算。CARS组合IRIV充分利用不同算法之间的互补性，构建更少、更有效的变量模型，对设备的开发和基础建设具有重要的指导意义。

1.3.5 模型的建立与评价

偏最小二乘判别分析（partial leasts quares discriminant analysis，PLS-DA）是光谱分析中常用的定性建模方法，基于最小二乘算法开发的监督线性算法，适用于解决变量较多且存在多重线性的问题[20]。PLS-DA将光谱变量视为独立输入（X变量），并将不同的类别信息分配为1、2、3、4、5等（即Y变量）。为了校准最优PLS-DA模型，本研究采用10折交叉验证法并依据最小均方根误差确定PLS-DA模型的最佳潜在变量（latent variables，LVs），在全光谱模型和简化模型中将LVs最大数量分别设置为40和30。

为了验证上述建立模型的性能，本研究中采用训练集、交叉验证集和预测集的分类准确率（correct classification rate，CCR）评估分类模型的性能，CCR计算如式（2）所示：

式中：CCR为分类准确率（正确分类的样本占全部样本的比率）/%；TP为真阳性；TN为真阴性；FN为假阴性；FP为假阳性。

为了进一步评估所选模型的结果，基于混淆矩阵的灵敏度、特异性和Kappa系数评价模型的分类性能，计算分别如式（3）～（6）所示：

式中：Pe为预期的随机准确率；P0为CCR；M为类别总数；i为类别数。

具体来说，灵敏度（又称真正率或召回率）量化了模型识别正样本的能力，特异性（又称真负率）则是模型在正确识别伪样本方面的效能。Kappa系数用于衡量分类模型的性能与随机分类的性能之间的差异。当Kappa系数位于0.8～1之间，表示模型的预测性能与真实值之间几乎完全一致，表明模型预测性能具有极高的可靠性和稳健性。

2 结果与分析

2.1 原始光谱

图3显示了从校准高光谱图像中提取的5 种不同产地枸杞样品的光谱曲线以及平均光谱曲线。5 种不同产地枸杞样品之间的光谱曲线表现出相似的轮廓，主要源于枸杞内部组织共性，而反射率的强度差异主要受到其内部化学成分含量影响。多糖、黄酮、总糖和多酚是枸杞主要的活性成分，它们也是衡量其内部质量的主要特征指标，这些化学成分在不同产地的枸杞样品中含量不同[21]。从图3b可以看出，在可见光区域内的400～550 nm波段，不同产地枸杞的平均光谱反射强度几乎相同，表明它们在蓝色和绿色的颜色分量上几乎没有差异。然而，在550～700 nm的波段区域，各产地枸杞的平均光谱反射强度开始有所不同，且呈现出明显的上升趋势，表明红色波段内存在的颜色分量差异较为明显。在波长450 nm附近的反射峰主要与枸杞中的酚类物质阿魏酸相关[22]。550 nm波长处的波谷是总糖的有效波长[21]，而在波长560 nm附近的吸收峰与枸杞表面叶绿素和类胡萝卜素的吸收带相关[23]。在近红外区域，黄酮的有效波段为800～900 nm，而在860 nm附近吸收峰为枸杞黄酮的有效波长[18]。在910～960 nm波段可归因于水或碳水化合物的O—H拉伸模式的第二泛音[24]。由于5 个不同产地枸杞光谱之间存在很多交叉和重叠，因此将光谱与化学计量学方法相结合进行深度分析，并作出准确判别。

图3 5 种不同产地的枸杞样品光谱曲线Fig.3 Spectral curves of L.barbarum samples from five different geographical origins

2.2 基于全光谱的判别模型

对原始全光谱应用不同的预处理方法，并采用原始和预处理全光谱建立PLS-DA枸杞产地判别模型。为了确定最佳的预处理方法，本研究基于5 个产地的枸杞数据建立模型，结果如表1所示。结果表明，无论采用原始光谱还是预处理光谱PLS-DA模型，分类准确率均大于90%，表明所建立的模型可以轻松区分不同产地的枸杞样本。在PLS-DA模型中，LVs的选择会严重影响最终结果。当LVs过多时，会造成过拟合。相反，如果LVs的数量太少，会丢失一些有用的信息，降低模型的准确性。因此，通过交叉验证计算最小预测残差平方和以确定最优LVs[25]。图4显示了不同LVs对原始光谱PLS-DA判别模型性能的影响。还可以观察到，NR和SNV预处理均提升了模型的性能，而SNV＋去趋势、一阶导数和二阶导数略微降低了性能。其中，基于NR预处理全光谱的PLSDA模型表现最佳，最佳LV为30，训练集分类准确率为95.5%，交叉验证集分类准确率为91.9%，预测集分类准确率为93.1%。NR预处理可有效抑制部分光照差异的影响和消除无关信息，这一结论与徐新刚等[26]的研究结果一致。因此，在随后的分析中最终选择NR预处理光谱用于枸杞产地溯源鉴别。

表1 基于不同预处理方法全光谱PLS-DA模型性能Table 1 Performance of full-spectrum PLS-DA models based on different pre-processing methods

图4 基于不同LVs原始光谱PLS-DA判别模型性能Fig.4 Performance of the raw spectral PLS-DA discriminant model based on different LVs

本研究对枸杞产地多元化鉴别需求进行深入分析，采用原始光谱和NR预处理光谱构建PLS-DA产地溯源判别模型，研究结果如表2所示。可以发现，随着枸杞样本产地的增加，无论采用原始光谱还是预处理光谱PLS-DA模型分类准确率总体呈下降趋势。当模型输入枸杞产地数量从2增加到5时，经预处理后模型的训练集分类准确率从99.2%下降到95.5%，交叉验证集分类准确率从94.8%下降到91.9%，测试集分类准确率从98.3%下降到93.1%。

表2 基于PLS-DA构建枸杞产地多元化判别模型Table 2 Multivariate identification model based on PLS-DA for geographical origins of L.barbarum

2.3 特征波长选择

全光谱数据包含大量的冗余信息，特征波长的选择有利于降低数据维度，这对于开发实时多光谱检测系统非常需要[27]。在本研究中，SPA、CARS、PSO、IRIV和CARS＋IRIV对4 组光谱集筛选出的特征波长分布如图5所示（为了便于观察，同时将5 种特征选择分布情况置于一张图中，其中纵坐标值按比例增加）。当枸杞产地数量为2时，相较于其他产地数量所筛选的特征波段较少，其中5 种方法选择的特征波长都相对分散不连续。当模型输入为5 个产地枸杞样本数据时所选的特征波长相对较多，但仅占全光谱的14.3%～42.4%。可以发现不同的特征变量选择方法将选取不同数量的特征波长，因此确定最优的变量选择方法对于构建高质量判别模型至关重要。枸杞产地溯源的进一步分类将使用基于选择的特征波长进行PLS-DA建模。

图5 不同方法选择的特征波长分布Fig.5 Distribution of characteristic wavelengths selected by different methods

2.4 简化模型

基于不同特征波长选择方法选取的特征波长构建多元化枸杞产地溯源PLS-DA判别模型，结果如表3所示。与全光谱PLS-DA模型性能相比，结果有所下降。由于特征波长极大简化了分类模型，可能造成这种性能下降。整体来看，基于SPA筛选特征波长建立的模型性能表现较差，主要原因是所筛选的波长主要集中在900～1 000 nm，无法提供全部的有效信息。虽然使用CARS和IRIV选择的波长建模取得了优异的结果，由于选择的波长较多不利于光谱检测系统的开发，对于PSO-PLS-DA和CARS＋IRIV-PLS-DA模型表现出令人满意的结果，所选的特征波长仅占全波长的15.6%～27.7%，在二元分类模型中预测集分类准确率分别为96.0%和97.7%，在三元分类模型中预测集分类准确率分别为90.0%和90.9%，在四元分类模型中预测集分类准确率分别为86.7%和89.2%，在五元分类模型中仍取得84.1%和87.1%的分类结果。

表3 基于不同波长选择方法的最优PLS-DA建模结果Table 3 Results of optimal PLS-DA models based on different wavelength selection methods

为了进一步探索CARS＋IRIV-PLS-DA模型鉴别枸杞产地溯源的能力，图6给出了该简化模型预测集的混淆矩阵、灵敏度和特异性以及Kappa系数计算结果。在混淆矩阵中，纵坐标表示实际类，横坐标表示预测类。主对角线内的值表示正确分类的样本，主对角线外的值表示错误分类的样本。结果表明，内蒙古和新疆枸杞样本发生错误分类数量最少（灵敏度均大于94%），可能由于内蒙古和新疆枸杞相较于其他产地枸杞含糖量高，具有较好的区分性。还可以观察到宁夏枸杞易被错误分类成内蒙古枸杞，在二元分类模型中由于模型简单宁夏枸杞的识别率高达96.7%。随着输入产地数量的增加，宁夏枸杞识别率和Kappa系数整体呈下降趋势，在五元分类模型中宁夏枸杞仍取得了82.7%的识别率。4 组简化模型的Kappa系数均超过0.83，说明分类模型具有较强的稳定性和鲁棒性。这些结果表明，CARS＋IRIV-PLS-DA模型在没有任何化学或物理信息的情况下识别宁夏枸杞产地溯源具有巨大的潜力。

3 结论

鉴别枸杞产地的传统方法耗时且具有破坏性和主观性，因此，本研究旨在开发一种基于高光谱成像技术，从而快速、无损识别宁夏枸杞产地的检测方法。通过高光谱成像采集宁夏、甘肃、内蒙古、青海和新疆共2 250 个枸杞样本，并从ROI提取光谱数据。研究发现相比于SNV、SNV结合去趋势、一阶导数和二阶导数，NR可以更好地降低光谱噪声和散射效应。为了降低数据维度并进一步减少建模时间，采用SPA、CARS、PSO、IRIV和CARS＋IRIV选择特征波长，并基于特征波长建立PLS-DA判别模型。随着枸杞产地数量的增加，模型性能呈下降趋势。当仅输入两个枸杞产地数量时，全光谱模型分类准确率高达98.3%，鉴于实用性最佳简化模型CARS＋IRIV-PLS-DA分类准确率高达97.7%。当输入为5 个枸杞产地数量，简化模型CARS＋IRIV-PLS-DA仍能获得87.1%的分类准确率和0.839的Kappa系数。综合分析表明，高光谱成像技术（400～1 000 nm）结合化学计量学方法可以作为一种快速、无损的检测方法鉴别宁夏枸杞的真伪性。