张灵枝, 黄 艳, 于英杰, 林 刚, 孙威江
(1. 福建农林大学园艺学院,福建 福州 350002;2. 福建农林大学安溪茶学院,福建 泉州 362400;3. 中国茶叶流通协会,北京 100801;4. 福建融韵通生态科技有限公司,福建 福州 350025;5. 福建农林大学福建省茶产业工程技术研究中心,福建 福州 350002;6. 福建农林大学海峡两岸特色作物安全生产省部共建协同创新中心,福建 福州 350002)
茶叶富含茶多酚、氨基酸、生物碱等多种对人体健康有益的成分,是世界上最受欢迎的三大无酒精饮料之一。 根据加工工艺及发酵程度不同,茶叶可分为绿茶(未发酵)、白茶(微发酵)、黄茶(轻发酵)、乌龙茶(半发酵)、红茶(全发酵)以及黑茶(后发酵)[1-3]六大类。 六大茶类识别主要通过专业人士对茶叶外形、风味进行感官审评,结果易受审评人员的身体、精神状况等影响,审评主观性和经验性强。 国际层面,茶叶评价术语尚未实现规范化、统一化和标准化, 化学分析的国际标准制定严重滞后,使得国外非专业人士对六大茶类的分辨更加困难[4-6],严重制约着茶叶国际贸易和流通。 因此,不少学者在化学分析技术、计算机视觉技术、光谱技术等领域开展了多茶类判别研究。 Ning 等以儿茶素、咖啡碱等生化成分含量结合Fisher 判别分析对六大茶类进行逐步判别, 其最佳模型的识别正确率达到88.30%以上[7];Wu 等以茶多酚、咖啡碱等生化成分含量结合逐步线性回归判别分析,建立绿茶、白茶、乌龙茶以及红茶的判别函数,识别正确率达97.80%[8];Jiang 等基于氨基酸、儿茶素等生化成分含量结合线性判别分析, 对六大茶类可实现100.00%准确鉴别[9];Zhang 等基于非靶向代谢组学识别绿茶、 白茶以及黄茶,找出了关键分类指标,分类识别正确率可达100.00%[10]。上述研究虽然建立了优异的判别函数和模型,但存在特征成分检测耗时长、成本高、分析复杂等问题,不易进行快速无损鉴别,难以在国际上实现产业化应用。 Wu 等通过采集红茶、绿茶、乌龙茶的外形图片, 利用机器学习和计算机视觉技术,结合K 近邻(K-nearest neighbor,KNN)分类器构建判别模型,识别正确率可达94.7%[11]。Ning 等将高光谱成像技术结合化学计量学建立的Lib-SVM 模型对绿茶、黄茶、白茶、黑茶、乌龙茶的识别正确率可达98.39%[12]。 由于茶叶的年份差异导致的波动性较大、高光谱仪普及率低等问题,上述两项技术并未实现产业化应用。
近红外光谱作为绿色分析技术[13],具有高效、便捷、准确性高等优点,已在食品[14-16]、炼油[17-18]、药物[19-20]等领域广泛应用。本研究中以NIRS 结合机器学习算法,探索六大茶类鉴定的可行性,包括以下4个步骤:1)获取六大茶类在3 600~12 500 cm-1波段的近红外光谱;2) 应用最小最大归一化(minmax scaler,Minmax)、连续小波变换(continuous wavelet transform,CWT)、 标准正态变换 (standard normal variate,SNV)及多元散射校正(multiplicative scatter correction,MSC)4 种预处理算法, 建模分析他们对OS 的去噪和散射校正性能;3) 比较主成分分析(principal component analysis,PCA)、线性判别分析(linear discriminant analysis,LDA) 和连续投影算法(successive projections algorithm,SPA)3 种方法提取茶叶光谱特征的能力;4)建立基于RF 和SVM 分类器的茶类鉴别模型,实现对六大茶类的快速、无损识别,为近红外光谱技术在茶类识别的产业应用上奠定理论基础和科学依据。
共收集370 份六大茶类样品,包括122 种乌龙茶、110 种绿茶、55 种红茶、34 种黑茶、29 种白茶、20 种黄茶,来自中国福建、广东、台湾等地区,以及日本、斯里兰卡等国家。 所有茶叶样品均具有该茶类正常的商品外形及特有的色、香、味,无异味、无劣变、无污染、无非茶类夹杂物、无任何添加剂,满足实验材料要求,具体信息详见表1。
表1 各茶类样品信息Table 1 Sample information of each tea category
MPA 型傅里叶变换近红外光谱仪:德国布鲁克光谱仪器公司;高速粉碎机:上海鼎广机械设备有限公司;CFJ-II 茶叶筛分机: 杭州大吉光电仪器有限公司。
1.2.1 样品制备茶叶粉碎后,置于茶叶筛分机中过80 目筛,取筛下茶粉5 g,密封编号,放于4 ℃冰箱中备用。
1.2.2 光谱采集使用MPA 型傅里叶变换近红外光谱仪采集样品光谱信息。 仪器工作时温度控制在25 ℃,相对湿度<70%。 光谱采集工作流参数:波数为3 600~12 500 cm-1,光谱扫描次数为64 次,分辨率为8.0 cm-1。 为确保近红外光谱检测数据的可靠性,每个样品扫描3 次,取平均光谱作为原始光谱数据进行后续分析。
1.2.3 光谱预处理采用Minmax 算法增强数据;选用CWT 算法校正基线漂移并消除高频噪声;使用SNV 及MSC 算法校正散射, 消除因茶粉粒径的不均匀、光程不恒定等因素所带来的影响。
1.2.4 特征提取为提升模型性能、 运算效率,选用PCA、LDA 及SPA 方法进行光谱数据的特征提取,降低数据维数。
1.2.5 模型构建与评价数据挖掘分类器广泛应用于NIRS 数据的分析与利用, 不存在始终保持最优效果的分类器,因此使用多种分类器建模更利于优质模型的构建[21]。 本文中基于RF、SVM 两种分类器,结合不同预处理、特征提取方法,优化模型参数,探究六大茶类最佳识别模型构建流程。
为确保所建模型的适用性, 将数据按照3∶1 的比例划分为训练集和验证集两个子集,样本数分别为277 和93 个,其中训练集用于模型训练,验证集用于测试模型稳健性。使用RA、AUC 以及混淆矩阵作为模型精度及性能的评价指标。
1.2.6 数据处理软件数据处理软件包括MATLAB 2016a、Origin 2019b、Excel 及Python。
370 份茶叶样品在3 600~12 500 cm-1波段的近红外光谱如图1(a)所示,各波段吸光度变化趋势趋于一致。 随着波数的增加,吸光度总体呈现下降趋势,变化范围处于0.249~2.196。
图1 370 份茶叶样品近红外光谱和六大茶类平均光谱图Fig. 1 Near-infrared spectra of 370 tea samples and average spectra of six major tea types
六大茶类的平均光谱在大多波段范围内趋势一致(见图1(b)),吸收光谱趋于平行,各茶类对应谱图于3 750 cm-1处初步分开,于3 750~9 000 cm-1处波动最为明显。受C—H+C—H 组合频伸缩、变形振动影响,4 000 cm-1处出现明显波峰;参照相关文献[22-25],4 500~4 950 cm-1的谱带变化可归因于N—H+O—H 的组合频;5 000~5 300 cm-1的谱带变化主要受O—H+O—H 的组合频振动影响;5 600~6 300 cm-1的谱带变化受C—H 和S—H 的一级倍频振动影响;6 700~7 400 cm-1的谱带变化主要受N —H 的一级倍频振动影响;8000~9000cm-1中出现的波峰可能与CHCH及C—H的二级倍频振动有关。
3 750~9 000 cm-1中各平均光谱能基本分开,说明六大茶类样品吸光度在该波段增减性不同,即光谱信息与茶类间具有相关性。 首段谱带中各平均光谱间交叉重叠现象频发,说明该波段体现的与产地相关的有效信息较少,信噪比低;末端谱带波动同原始光谱一致,趋于平缓,无明显波峰波谷,特征信息不显,因此结合Meng 等结论[23],本文中的模型构建使用3 750~9 000 cm-1波段的光谱数据。
为校正光谱采集过程中因环境、光程不恒定及样品粒径差异等因素所带来的误差, 从数据增强、基线校正、 散射校正等角度, 使用Minmax、CWT、SNV 及MSC 算法对OS 进行预处理。 CWT 中小波参数(小波基、分解尺度)的选择至关重要,直接决定了后续模型的优劣,经比对分析后,本研究中选择应用最为广泛的db(daubechies)族中的db4 小波基,分解尺度定为100[26]。
从光谱变化情况(见图2)可知,4 种处理方式都使光谱形态发生了较大改变。Minmax(见图2(a))将光谱吸光度凝练到-1.0~1.0,增强了数据,消除了数据量纲及取值范围的影响,后续可使所建模型收敛速度加快,提高模型性能。 应用CWT 进行光谱预处理如图2(b)所示,其形态变换程度为4 种预处理方式中最大,基线漂移、背景干扰、噪音现象等得到明显消除,谱峰更清晰,差异信息段更明显。 受茶粉颗粒大小不均、产生的散射影响,采用SNV 与MSC进行预处理(见图2(c)和图2(d)),处理后光谱中的散射干扰得到明显消除,特征信息更加突出。 相较于OS,预处理可有效消除光谱中因光散射、基线漂移等造成的信号干扰,但处理后的光谱图仍无法直观分辨茶类间的差异,这可能是由于不同茶类在内含物的组成与含量上具有较多的相似性。 为了进一步评估各处理对建模结果的影响,将各处理所得的光谱吸光度分别作为模型的输入变量,后续基于模型评价指标判断预处理效果。
图2 光谱预处理效果Fig. 2 Effect of spectral pretreatment
近红外光谱的连续波数中存在大量的冗余信息,其与特征信息之间存在很强的相关性[27]。 通过选取特征向量或波数降低数据维数,可保留原始数据中的主要特征信息, 减少后续处理的计算任务。本研究中采用PCA、LDA 和SPA 方法进行数据降维。
2.3.1 PCAPCA 是一种常用于降低大数据集维数的无监督特征提取方式,能从大量数据中提取出特征,转换为仍包含绝大部分有效信息却拥有较小维数的数据集, 最大程度保留原始数据信息, 因而PCA 是一种最优、最为常用的方法[28]。 使用PCA 方法对3 750~9 000 cm-1中OS 及经Minmax、CWT、SNV 和MSC 预处理后的4 种光谱进行降维, 截取前15 个主成分特征值与累积贡献度。 结果如表2所示, 以特征值大于1、 累积贡献度大于80%为原则, 筛选模型输入主成分个数。 OS 以及Minmax、CWT、SNV 和MSC 预处理后光谱分别筛出6、11、13、12、12 个主成分,累积贡献度分别达到99.89%、99.71%、99.67%、99.82%、99.82%,符合原则。基于筛选的主成分构建模型。
表2 PCA 特征值及累积贡献度Table 2 PCA feature values and cumulative contribution
2.3.2 LDALDA 是一种有监督的特征提取方法[29],在茶叶领域常作为分类器使用。而利用LDA 进行光谱特征提取、降低数据维数,并结合分类器建立茶类识别模型的研究,目前尚未见相关报道。 LDA 最多可使数据矩阵降至类别数减1 的维数,降低维数的同时不过多丢失原始信息,LDA 降维后所得维数将被用于后续建模。
2.3.3 SPASPA 是一种前向循环特征提取方法,其可以通过连续投影的方式从原始光谱矩阵中提取有效预测响应变量的信息,最大限度减少光谱变量之间的共线性效应,达到所选响应变量预测能力的最大化。 其主要通过将波数投影于其他波数,比较投影向量大小, 波数间投影向量大者为待选波数,最终投影向量最大且与特征集内波数共线性最小的波数选入特征集合[30]。 特征波数的数目由校准集内部完全交叉验证的均方根误差 (root mean square error,RMSE)确定,与最小RMSE 对应的特征波数数目和特征波数为最佳值[31-32]。
由SPA 方法提取的特征向量如图3 所示。 SPA中4 种预处理后的光谱信息的RMSE 迭代下降曲线分别如图3(a)、图3(c)、图3(e)、图3(g)和图3(i)所示。 从图中可以看出,当选择特定数量的波数时,RMSE 达到最小值; 而后RMSE 虽仍波动下降,但降幅很小且导致所选波数增加,没有必要为了追求微小的RMSE 而增加维数。因此,最终从OS 及经Minmax、CWT、SNV 和MSC 预处理光谱中获得的特征波数数目分别为15、16、19、15、7 个 (如图3 所示,特征波数具体信息如表3)。
图3 连续投影算法(SPA)提取特征波数Fig. 3 Extraction of feature wavenumbers by successive projections algorithm (SPA)
表3 通过SPA 筛选出的特征波数Table 3 Feature wavenumbers selected by SPA
SVM 是近年来茶叶中应用最广、效果最好的机器学习方法之一。 它是一种利用核函数将n 维输入向量映射到K 维特征空间(K>n),从而通过高维特征空间进行分类的算法[33]。为提高模型质量,本文中所有SVM 模型皆基于高斯 (radial basis function,RBF)核函数,该核函数可降低训练过程的计算复杂度,在一般平滑假设下具有良好性能;与此同时,惩罚参数及gamma 参数最优值的确定也至关重要,SVM 模型精度取决于这两个参数的组合。根据初步试算结果,将惩罚参数取值定为1×103、1×104、1×105、1×106,gamma 参数取1×100、1×10-1、1×10-2、1×10-3。将3 750~9 000 cm-1的OS 数据作为输入量,结合模型的识别正确率进行参数优化。 结果如表4 所示,当惩罚参数为1×106、gamma 参数为1×10-2时模型具最佳识别正确率,后续模型皆基于此参数构建。
表4 SVM 模型参数优化Table 4 Parameter optimization of SVM model
RF 是一种有监督的集成分类算法, 主要是为解决单一决策树可能出现的很大误差和过拟合的问题,在分类问题中表现优异,具有成为各情况下效果最优分类器的巨大潜力[34]。 该模型由许多的决策树组成,但每一棵决策树之间没有关联,得到森林之后对新样本进行判断或预测时,将由森林中的每一棵决策树分别进行判断, 分辨该样本属于哪类,比对出选择数最多的类别,从而对样本类别做出判断, 因此该模型中树木数量的选择极为重要。OS 数据试算后, 选定树木数量为1~100 进行参数优化,结果如图4 所示,树木数量为70 时,RF 模型具最佳识别正确率,后续模型皆基于此参数构建。
图4 RF 模型参数优化Fig. 4 Parameter optimization of RF model
RA 常被用于模型预测能力的评估,AUC 则常被用于模型泛化能力的评价, 取值区间为0.5~1.0,值的大小与模型质量呈正相关。 因此, 采用RA、AUC 联合评估模型性能。 此外,为直观呈现所建模型对各茶类识别性能的优劣,引入混淆矩阵对模型预测结果进行评价。
由表5 可知, 各茶类NIRS 数据结合不同预处理、 特征提取方法及数据挖掘分类器最终获得40个茶类识别模型, 识别正确率介于59.14%~100.00%,AUC 处于0.70~1.00, 大多数模型识别正确率高于70%,模型性能良好。 OS 结合RF、SVM 所得模型识别正确率分别为69.89%及92.47%,光谱经预处理、特征提取后,建立的RF 模型绝大多数识别正确率、AUC 显著提升,模型精度、泛化能力均得到改善, 最佳模型OS-LDA-RF 的识别正确率可达94.62%,AUC 可达0.96;SVM 结合预处理后光谱建模效果欠佳,多出现识别准度下降的问题,不同特征提取方法中效果最佳的为LDA,结合不同预处理后的光谱数据皆优化了模型性能, 基于SVM 建立的茶类识别模型中最佳的是OS-LDA-SVM,识别正确率为100.00%,AUC 为1.00。 总体而言,不同分类器结合不同预处理、特征提取方法所取得的效果也不尽相同,在茶类识别模型构建中,RF 适合与不同化学计量学方法结合,多数预处理、特征提取方法对RF 模型性能优化效果显著;SVM 适合基于去除头尾信息匮乏波段的原始光谱结合特征提取方法进行建模,所得模型特征数可得到简化,提升运算速度及模型质量。
观察最佳模型混淆矩阵可知 (见图5),OSLDA-RF 对乌龙茶、 黑茶、 白茶的识别正确率为100.00%,误判发生于绿茶、红茶及黄茶的识别中,三者识别正确率分别为89.28%、92.86%以及80.00%;OS-LDA-SVM 对各茶类的识别正确率皆为100.00%,混淆矩阵中数值皆处于对角线,表明识别效果优异,模型质量好。
图5 最佳模型混淆矩阵Fig. 5 Optimal model confusion matrix
对性能最优模型OS-LDA-SVM 进行三维空间可视化(见图6)。 可发现,识别正确率达100.00%的茶类识别模型中,绿茶与黄茶的光谱特征于三维空间上的分布极为接近,可能与其加工工艺的高度相似有关,黄茶仅比绿茶多了闷黄工艺;还可能受黄茶样本数量较少导致特征信息不显的影响,黄茶在六大茶类中占比最小,仅在我国四川、湖南等少数省份小规模生产,后期将通过逐年增加黄茶样本数量的方式,强化黄茶光谱特征,提高模型性能。 其他茶类的光谱特征在三维空间分布差异大,可能与加工工艺、茶树品种不同有关。
图6 OS-LDA-SVM 模型三维空间效果Fig. 6 Three-dimensional space effect of OS-LDA-SVM model
基于茶类的370 个近红外光谱数据,剔除认为不含茶类相关信息的波段后, 使用Minmax、CWT、SNV 及MSC 进行预处理,PCA、LDA 及SPA 进行特征提取,最后基于RF、SVM 构建茶类识别模型。 主要结论如下:
1)茶类识别模型构建中,不同预处理对不同分类器的效果不尽相同。对RF 模型而言,预处理对模型性能提升效果显著, 而相同处理于SVM 模型中效果欠佳,模型质量多不如原始光谱模型。
2)通过特征提取优化模型,效果显著。 PCA、LDA、SPA 方法皆大幅度降低了数据维数,提高了模型运算效率。 其中LDA 效果最佳,与不同预处理方法、分类器结合所得模型皆质量优异。
3)RF、SVM 皆适用于茶类识别模型构建, 相较RF 模型,SVM 模型总体效果略胜一筹。 RF 模型中最佳模型为OS-LDA-RF,RA 为94.62%, 对乌龙茶、 黑茶及白茶的RA 可达100.00%,AUC 为0.96,模型性能优异、稳定;OS-LDA-SVM 为SVM 模型中的最优模型,不同茶类的RA 皆达100.00%,AUC 为1.00,模型质量高、泛化能力好。
不同茶类的近红外光谱数据, 经适当预处理、特征提取方法的选择后,结合RF、SVM 分类器可挖掘出近红外光谱中茶类识别相关的关键信息,构建出高识别正确率的茶类识别模型。 基于此,针对光谱特征接近的茶类,可逐年扩大样本量,优化与验证茶类识别模型性能,提高模型适用性。 除此之外,本研究后续将利用互联网技术构建六大茶类的近红外光谱数据库,搭建在线茶类识别平台,开展远程茶类识别,以期早日在国际市场上实现产业化应用,促进我国茶产业的高质量发展。