基于高光谱成像技术鉴别大米品种

王靖会, 程娇娇, 刘洋, 常佳乐, 王朝辉

(1.吉林农业大学信息技术学院， 长春 130118； 2.吉林农业大学食品科学与工程学院， 长春 130118)

大米是生活中必不可少的主食之一，不同品种大米除外观存在一定差异，其内部的淀粉、蛋白质、脂肪酸等含量也不同，故而品质及价格差别很大，导致市场上利用形似优质米的普通米冒充或掺杂优质米现象频发，造成了巨大的经济损失，同时给大米的国际贸易、地理标志产品保护以及消费者权益保护带来阻碍。2019年出台的《关于深入实施优质粮食工程的意见》指导文件中提出“要促进农民增收、企业增效、消费者得实惠”[1]，更确定了大米品种鉴别的重要性。

目前，常用大米品种鉴别方法有电泳法、分子标记法、拉曼光谱技术、人工嗅觉技术、近红外光谱技术等，其中人工嗅觉技术误差较大，电泳法、分子标记法、拉曼光谱技术和近红外光谱技术存在样品准备工作繁琐、耗时长、不能无损检测等缺陷。高光谱成像技术具有“图谱合一”的特点，能够同时获得样本的图像和光谱信息。图像信息可以反映样本的大小、形状等外部特征，光谱信息可以充分表达样本内部的物理结构和化学组分，其将光谱信息和图像信息融合使用，更能表达品种的特征，这决定了高光谱成像技术在农产品检测领域中的优势地位[2-8]。但高光谱图像属于高维非线性数据，存在数据量大、噪声、冗余及多重共线性等问题，另外，大米品种鉴别的样本来源大都是地域分散、外观差异大的品种，而现在市场上通常为外观相似的品种掺伪，增加了鉴别难度[9-12]。因此，光谱数据处理和建模方法的选择以及计算效率的提升成为利用高光谱图像技术鉴别品种的关键。常用的建模方法有神经网络、决策树、支持向量机等，其中支持向量机(support vector machine, SVM)对解决非线性及高维数据问题效果较好，但其分类精度受参数影响较大。乌鸦搜索算法(crow search algorithm, CSA)是一种智能优化算法[13]，具有收敛速度快、全局搜索能力强、输入变量少等优点，在优化约束工程设计的参数选择上也有比较广泛的应用。本研究采集外观相近的6个品种共600粒大米的高光谱图像，利用支持向量机算法分别建立基于全波段、特征波段、图像纹理特征以及光谱-纹理特征融合的大米品种鉴别模型，利用CSA算法进行优化，通过对比分析探索实现大米品种快速无损检测的新的方法。

1 材料与方法

1.1 数据来源

1.1.1样本采集 在吉林省通化地区采集龙稻20、中科发5号、中科804、稻花香、龙洋11、长粒香共6个外观近似的大米品种作为研究对象(标号为R1～R6)，经晾晒、除杂、砻谷、精米等预处理后，挑选色泽均匀、颗粒完整的大米，每个品种100粒，共600粒样本，利用Kennard-Stone算法[14-15]把样本按照近似2∶1的比例划分成训练集(66)和测试集(34)。

1.1.2高光谱数据采集及校正 利用高光谱图像采集系统(图1)采集大米高光谱数据，为确保光源能稳定照射，在采集前预热光谱仪。该系统由成像光谱仪(Imspector V10E-QE，Spectral Imaging Ltd.，Oulu，Finland)、C8484-05G相机(Hamamatsu Photonics，Japan)、V23-f/2.403 060 3镜头(Specim，Finland)、P/N 913 0线光源、2 900 ER控制器(Illuminati on Technologies，USA)、GZ02DS20可升降样品台(光正公司)及PSA200-11-X电控位移台(卓立汉光公司)组成。将大米样本按照5×5网格分布，以单粒随机摆放至载物台黑板上，经多次实验确定最佳参数：物距为13.5 cm，曝光时间为7 ms，位移台移动速度为3.4 mm·s-1，波长范围为400～1 000 nm，波段数为477个，采集图像时对每个样本都进行三次扫描。

对采集到的大米高光谱图像按照式(1)进行校准，以降低外界环境如光源照射不均匀、传感器响应不同等造成的样品差异，高光谱图像采集和黑白校正均由高光谱采集软件Spectracube 2.75b完成。

(1)

式中，I0为校正前光谱图像；W为反射率是100%的全白标定图像；B为盖上镜头盖并关闭电源之后获取的反射率为0的全黑标定图像。

1.2 感兴趣区域选择及光谱预处理

采用ENVI5.0选取6个品种大米样本的感兴趣区域(region of interest，ROI)，一个单粒米为一个ROI，所有像素的平均光谱表示该粒米的样本信息。对同一粒样本的三幅图像分别提取ROI然后取平均，每个单粒大米可以得到一条代表样本的平均光谱曲线，600粒大米获得一个600×477的光谱矩阵。利用软件MATLAB2017a中多元散射校正(multiplicative scatter correction, MSC)、标准正态分布(standard normal variate, SNV)和二阶导数(2ndderivative,2ND)三种算法对样本光谱数据进行预处理[16-17]，并建立模型。根据模型的判别效果确定最佳预处理方法，以消除基线漂移，背景噪声及其他随机噪声的影响。预处理算法及模型建立均在软件MATLAB2017a中完成。

1.3 特征波长提取

原始高光谱图像包含477个波段，波段较多且相邻波段间的相关性较高，导致光谱信息大量冗余，因此需要降低光谱维度以提高模型的计算效率。本文使用连续投影和主成分分析方法实现特征波长提取，算法及模型建立均在软件MATLAB2017a中完成。

1.3.1连续投影算法 连续投影算法(successive projections algorithm, SPA)[18]首先设定特征波长数目的范围，以一个波长开始，在迭代中合并波长，达到波长数N，用N个特征波长代替全波段，其中每迭代一次都要对波长组合进行多元线性回归分析，得到验证集的预测标准偏差(root mean square error of calibration,RMSEC)，最小的RMSEC值为最优值，对应的波段即为最优特征波段组合。

1.3.2主成分分析 利用主成分分析(principal component analysis，PCA)对高光谱数据降维[19-20]，根据权值系数的绝对值与其对应的波长贡献率成正比[21]，从前几个累计贡献率大的主成分权值系数曲线中提取特征波长，即权值系数曲线中波峰波谷的位置。

1.4 图像纹理特征提取

高光谱图像技术相较于其他传统光谱技术的优势是可以提供样品的形状、大小、颜色、纹理等丰富的图像信息，纹理作为图像的重要特征之一，是分类的重要依据之一[22]。灰度共生矩阵(gray level co-occurrence matrix, GLCM)是基于共生矩阵的二阶统计量的纹理分析技术，它测量特定灰度级别的某两个像素在指定方向上和一定距离的联合分布概率。将像素间的距离参数d设为1，计算0 °、45 °、90 °和135 °四个方向的灰度共生矩阵，选取熵(entropy)、对比度(contrast)、相关性(correlation)、协同性(homogeneity)、二阶矩(second moment)及协方差(variance)作为纹理信息表征量。灰度共生矩阵算法利用软件MATLAB2017a完成。

1.5 支持向量机

支持向量机(SVM)旨在找到解决分类问题的最优超平面，使得间隔最大以实现不同类别数据的划分，在高光谱图像技术分类中使用较广泛[23-26],与其他机器学习算法相比泛化能力强。核函数的选择是训练样本前比较重要的步骤，RBF核函数适用范围较广，无论是小样本、大样本还是低维、高维等情况均可，已在许多研究中得到成功应用。在使用RBF-SVM分类时需要设置合适的参数C和g，C是惩罚因子，g是核函数中的gamma函数设置，两个数值过大或过小都会影响模型的泛化能力，利用MATLAB2017a建立模型。

1.6 乌鸦搜索算法

乌鸦搜索算法(CSA)是近年来兴起的一种智能优化算法，通过对乌鸦生活习性的模拟来找寻目标最优解，即合适的适应度函数选择，本文将分类准确率作为适应度函数。利用CSA算法对SVM模型的参数C和g进行优化，寻找最优组合使得模型分类效果最佳。算法参数设置：飞行长度为 2.5，意识概率为0.1，种群数量的大小为15，最大迭代次数为100。惩罚因子C 的取值范围是[0.01，20 000]，核参数g 的取值范围是[0.01，100]。利用软件MATLAB2017a实现优化算法，流程如图3所示。

1.7 混淆矩阵

混淆矩阵(confusion matrix)也称误差矩阵，是一种衡量分类器分类准确程度的方法。在计算多标签混淆矩阵时，将多分类问题转化二分类，即one-V-Srest，某一类为正样本，其余为负样本。模型准确率由混淆矩阵计算得出，可以衡量模型对样本的预测能力。

(2)

式中，TP为真正例，即被判定为正样本且实际上也是正样本；TN为真反例，即被判定为负样本且实际上也是负样本；FP为假正例，即被判定为正样本实际上是负样本；FN为假反例,即被判定为负样本实际上是正样本。

图2 乌鸦搜索算法流程Fig.2 Crow search algorithm flow

2 结果与分析

2.1 大米预处理后的光谱分析

为消除设备以外的噪声影响，本研究采用MSC、SNV、2ND三种方法对大米原始高光谱进行预处理。将三种预处理后的光谱数据作为SVM分类模型的输入，其分类准确率如表1所示。可以看出，在同一模型下，未经预处理的原始光谱数据的模型整体分类准确率最低，测试集分类准确率仅为75.98%，而经过预处理后的模型整体分类准确率均有所提高，其中基于MSC预处理后的模型分类准确率最高，对大米品种分类效果最好，训练集和测试集的模型分类效果均达90%以上，表明利用MSC可以很好的消除光谱数据的散射影响，修正大米样本光谱间的相对基线偏移和平移现象。

表1 全波段分类准确率对比Table 1 Full-band classification accuracy

对降噪后的大米样本反射率求平均值，得到6个品种大米平均光谱曲线，图3为MSC预处理后的平均光谱反射率曲线，虽然大米品种不同，但光谱曲线走势、波峰、波谷位置基本相同，反射率不同，说明其内部组成成分相同，但成分含量不同。在849 nm处的吸收峰主要是由大米内部水分吸收引起的，是水分子中的O-H 基团在该波段的倍频的振动和延伸。

图3 六个品种大米ROI平均光谱曲线Fig.3 ROI average spectral curve of six varieties of rice

2.2 基于特征波长的模型分析

2.2.1特征波长的确定 由图3可以看出，品种R1、R4、R5和品种R2、R6均有重合部分，基于全波段光谱数据建立的分类模型会受到影响，同时全波段下的光谱数据信息冗余量大，在模型计算量及计算效率方面很不友好，因此需要对全波段光谱数据降维，以特征波长代替全波段。采用SPA进行光谱数据特征波长选取，将波长数的范围设定在5～30，从RMSEV值随波长数量变化的曲线(图4)可以看出，波长数在16之前，RMSEV值处于下降状态，而后趋于平稳，经多元线性回归分析计算可得RMSEV最小值为0.072 842，对应的特征波长个数为16个。选定的16个最优波长分布如图5所示，分别为415.63、426.84、429.25、461.91、470.59、520.32、580.36、861.89、881.17、892.75、910.76、927.48、939.06、969.65、975.46、978.56 nm。

图4 SPA的RMSEC最小值Fig.4 Minimum RMSEC of SPA

图5 SPA选取的特征波长Fig.5 Featured wavelength selected by SPA

2.2.2主成分分析 PCA主成分个数的选取一般有两种标准：一是累计贡献率≥85%以上，二是特征值≥1。经过计算，前三个主成分的累计贡献率已达97.11%(表2)，很好地表达了原变量信息。在前三个主成分权值系数曲线中，其波峰波谷的位置即为特征波长，选取的7个波长分别为429.84、712.26、737.7、761.93、791.34、867.03、914.62 nm。

表2 前三个主成分累计贡献率Table 2 Cumulative contribution rate of the top three principal components

2.2.3模型分类结果分析 将两种降维方法提取的特征波长分别输入SVM分类模型，分类结果如表3所示，虽然基于SPA和PCA提取特征波长的模型整体分类准确率相较于全波段的降低了，但是输入变量维度减少了96%以上，模型训练时间有了大幅度降低。而利用PCA降维的模型训练时间稍低于SPA降维的模型训练时间，但测试集分类准确率相较于SPA却降低了2.94%，表明利用SPA算法选取特征波长效果较好，在保证模型分类准确率的同时大幅度提升模型运算效率。

表3 模型分类结果Table 3 Classification result based on model

2.3 基于纹理特征建立模型

高光谱图像包含连续波长的灰度图像，如果从每个灰度图像中提取纹理特征，就会导致数据冗余难以处理，同时为了提高特征提取的稳健性[27]，本研究选取基于SPA、PCA提取的特征波长对应的大米高光谱灰度图像为研究对象，分别计算每个品种在四个方向上的灰度共生矩阵，对矩阵做归一化处理，得到每个品种四个方向的六个特征量值，然后取其平均作为样本的最终纹理信息,从最终每个品种得到6个纹理特征，如图6所示。从纹理特征平均值(表4)可以看出，不同品种的二阶矩值相近且较小，说明图像灰度分布不均。以纹理特征作为模型输入，大米品种分类模型的准确率如表5所示，可以看出，基于SPA特征波长下提取的纹理特征分类模型准确率明显高于PCA特征波长下的纹理特征分类模型准确率，两种组合方法下的品种分类模型整体鉴别率低于基于光谱数据建立的模型，表明单独利用纹理特征分类方法误差较大。

图6 R1在580.36 nm下的六个纹理特征Fig.6 Six texture features of R1 at 580.36 nm

表4 六个品种平均纹理特征Table 4 Average texture characteristics of six varieties

表5 基于纹理特征分类结果Table 5 Classification result based on texture characteristic

2.4 基于融合光谱和纹理数据特征建立模型

光谱特征波段与图像纹理特征在数值上有较大差异，直接将二者融合，大数值区间的属性会过分支配小数值区间的属性进而影响建模分类效果，因此在进行特征融合前需对数据做归一化处理。 分别将利用SPA、PCA提取的特征波长与6个纹理特征数据归一化后进行融合，将融合后的数据作为SVM分类器的输入，分类结果如表6所示，经过特征融合后的分类模型整体鉴别率优于单独使用光谱数据或纹理数据的分类模型鉴别率，融合特征后的模型结合了大米的外部和内部属性，可以更充分地说明大米品种的质量变化，表明经过特征融合的数据可以更充分的表达品种信息，补足了由单一特征分类的缺陷，提高分类模型的精度。

表6 基于融合数据分类结果Table 6 Classification result based on fusion data

2.5 优化模型分类结果分析

随着迭代次数的增加，每代最佳适应度(准确率)发生变化(图7)，从48代开始最佳适应度趋于稳定达到0.97。将CSA算法得到的最优(C，g)组合输出带入SVM分类模型计算分类准确率，结果如表7所示，优化后的两种数据融合后模型整体分类准确率均有所提高，能很好地区分6个大米品种，其中基于SPA数据融合下的模型优化前后的整体分类准确率高于基于PCA数据融合下的模型优化前后的整体分类准确率，表明在此实验中SPA的特征提取效果优于PCA，同时将乌鸦搜索算法与支持向量机结合具有合理性，为相关研究提供了思路。

图7 适应度曲线Fig.7 Fitness curve

表7 优化后分类结果Table 7 Classification result after optimization

3 讨论

传统的大米品种鉴别技术很难无损高效地进行大米品种鉴别，而高光谱成像技术可以从样本的外部特征和内里的化学成分着手进行大米品种的快速无损鉴别，但其获得的数据量大，冗余度高，导致实际在线检测时出现计算量大、数据处理不方便等问题，因而如何降低数据维度，提高模型分类准确率等是目前值得研究的问题。

本研究利用高光谱成像技术鉴别大米品种的，结果表明，利用SPA和PCA在进行特征提取是可行的，不同的方法原理在特征选择时存在差异，这与Yang等[18]利用SPA提取不同品种玉米的特征波段和樊阳阳等[21]利用PCA及SPA提取不同品种干枣特征波段的结论一致，同时本研究表明，利用特征波段代替全波段建立模型后的分类准确率降低，这是由于降维后丢失了一些关于品种的特征信息，但结合模型效率考虑，降维后的代价最小。单独使用光谱数据或图像数据进行大米品种鉴别的模型准确率低于图谱数据融合使用后的模型分类准确率且单独使用图像纹理特征进行大米品种鉴别准确率最低，这与黄蒂云等[5]利用高光谱成像技术进行脱绒棉种品种分类的结论一致，表明单纯依靠图像纹理鉴别外观近似的大米品种存在难度，而高光谱成像优势可以弥补单一特征的缺陷，提升鉴别准确率。为了进一步提高模型的分类性能，本研究利用CSA对SVM参数进行优化，优化后的模型分类准确率有所提高，表明该算法融合能够有效提高模型分类准确率。综上所述，以特征光谱融合纹理数据，利用SVM和CSA的融合算法，可以实现大米品种的快速无损鉴别，能够为其他研究提供思路和参考，接下来的研究目标是收集更多的粒形相似的大米品种，建立更加完善的大米品种鉴别模型。