赵 伟， 何 俊， 刘泉澄， 侯森林， 邓 琥

(1.西南科技大学 信息工程学院，四川 绵阳 621010； 2.西南科技大学 四川天府新区创新研究院,四川 成都 610299)

0 引 言

中药材传统经验鉴别法在鉴别产地时，对个人能力依赖较大，缺乏客观性。随着科技的进步，现代鉴别方法因能客观反映药材的内部质量而被广泛应用于中草药鉴别[1,2]。近年来，太赫兹技术也被用于中药材检测，众多学者还将机器学习与太赫兹光谱结合[3～9]实现中药材的鉴别，刘燕德等人[6]利用标准正态变量变换和主成分分析(principal component analysis,PCA)及支持向量机(support vector machine,SVM)3种算法结合太赫兹光谱，实现了对4种贝母的识别；Liang J等人[7]利用PCA-SVM粒子群优化(particle swarm optimization,PSO)算法对模型进行参数优化，实现了对黄芩的产地鉴别；Yang S等人[8]将全波段的太赫兹吸收光谱作为卷积神经网络(convolutional neural network,CNN)的输入量,实现了对3种产地咖啡豆的识别；Pan S等人[9]利用改进的杜鹃搜索(cuckoo search,CS)算法SPCS和SVM算法结合太赫兹光谱，实现了对3种年限人参的识别精度。

上述识别模型具有较高的识别率，但却缺乏普适性，即不能处理未知药材名称的产地识别问题。本文依据中医药的整体观思想选取了0.2～1.5 THz全波段的数据进行分析，提出了将太赫兹光谱结合机器学习的中药材产地识别方法。

1 材料与方法

1.1 太赫兹时域光谱系统

本文使用的太赫兹时域光谱(THz-time-domain spectrum,THz-TDS)系统如图1所示，系统由飞秒激光器(MaiTai，Spectra-Physics)和THz-TDS仪组成。激光器的中心波长为800 nm，带宽为10 nm，经立方体分束器分为泵浦光和探测光，当泵浦光照射光导天线时可产生太赫兹脉冲，其通过不同样品[17]后可产生不同的THz-TDS，THz-TDS仪为美国Zomega公司的Z-3型THz-TDS系统，该系统光谱范围是0.1～2.5 THz，Z-3系统设置扫描开始时间为50 ps，扫描长度为60 ps，扫描速度为0.122 mm/s，为了减少水分对光谱的干扰，将Z-3THz-TDS仪放置在密封箱中，在实验期间通入干燥空气，将相对湿度控制在4 %RH以下，并将重复4次扫描获得到的平均光谱作为一次有效的THz-TDS，经傅里叶变换及Dorney等人[10]提出的光学参数提取方法可以得到其太赫兹吸收光谱。

图1 透射式THz-TDS系统光路原理

1.2 样品来源与制备

实验所用中药材由绵阳市食品药品检验所提供，其中附片有江油、云南、汉中3种产地，其余6种药材为：大黄(甘肃和青海2种产地)、党参、陈皮、麦冬、天麻，实验前先将样品干燥3 h，然后粉碎并过200目筛网，之后将压片机的压力设置为5 T，压力保持时间为3 min，压制成直径约13 mm，厚度约1 mm的圆片。附片每种产地压制35片，共计105片，采集样品正反两面的太赫兹光谱，共得到210组附片的太赫兹吸收光谱数据；而其他品种药材各压制5片，共计30片，共得到60组太赫兹吸收光谱数据。

2 算法原理

2.1 分类算法

采用Python 3.7中的Sklearn学习库实现随机森林(random forest,RF)[11]、SVM和偏最小二乘判别分析(partial least square discriminant analysis,PLS-DA)等模型的建立，其操作较为便捷。其中，SVM采用线性核(linear kernel)函数，SVM的惩罚参数、RF的基分类器数量、PLS-DA的主成分数均通过绘制学习率曲线获得最佳值。另外，为保证训练的准确度，利用train_test_split()函数将不同产地附片的210组光谱按7︰3比例随机划分为训练集和测试集。

2.2 模糊判别分析算法

本文中的模糊判别分析(fuzzy discriminant analysis，FDA)算法是在给定标准库情况下，利用最大隶属度原则将待识别对象分类到对应标准库的方法。

1)样本库和细节点特征库[12]

将附片每种产地的50组原始吸收系数数据作为样本库。每种产地附片在每个频率采样点下的特征均值和标准差，用A(i)表示第i种产地药材的特征均值,C(i)表示第i种产地药材的标准差,即

A(i)=[mi1,…,mij,…,miN]

C(i)=[ci1,…,cij,…,ciN]

(1)

2)细节特征匹配

假设细节点特征库中包含B种产地的附片，现对待测样本Q的模糊判别过程做进一步说明，假设待识别药材吸收光谱记为FQ，则FQ=[fQ1,…,fQj,…,fQN]。若待测药材被分类器预测为第i种药材，则将待测光谱与细节点特征库中的第i种产地药材匹配，但由于样品厚度、系统噪声等因素的影响，FQ不可能恒等于第i种药材的均值矢量A(i),而是以高斯分布的形式分布在A(i)附近。计算待测样本属于中药材r的隶属度drj，则

(2)

其中，隶属度drj在0～1之间。求得待测药材与细节点特征库中某种产地隶属度后，根据最大隶属度原则以及设定的阈值即可验证分类结果是否正确，若隶属度大于设定阈值则分类结果正确并将其输出，否则分类错误输出结果为：未知药材。其识别流程如图2所示。

图2 3种模型的识别流程

2.3 分类模型阈值的选择与模型评价

本文将等错误率(equal error rate,EER)处的值，即误识率(FAR)和拒误率(FRR)相等处的值设定为识别模型的阈值。FAR=(NFA/NIRA)×100 %，其中，NIRA为类间测试的总次数，而类间测试集由不同品种的中药材样本组成，NFA为错误接受的样本数。FRR=(NFR/NGRA)×100 %，NGRA为类内测试的总次数，而类内测试集由不同产地的附片样本组成，NFR为错误拒绝的样本数。

在进行类内测试实验时，随着阈值的增大，错误拒绝的样本数量越来越多，即待测光谱必须与目标光谱高度相似才能被正确识别；而对于类间测试实验，随着阈值的增大，错误接受的样本数量越来越少，其他药材的光谱本就与附片光谱差异较大，更容易被识别为非附片药材。因此，将EER作为识别模型的阈值，此时识别模型在对类内样本的识别和对类间的样本识别，其错误率均最低且相等，如图3。

图3 FAR/FRR—阈值曲线示意

识别准确率A(Accuracy)是识别正确样本与总样本的比例，A=((TP+TN)/(TP+FP+TN+FN))×100 %。其中，TP为真阳性，TN为真阴性，FP为假阳性，FN为假阴性。

3 结果与分析

3.1 太赫兹光谱分析

由于1.5 THz后信噪比变差，所以仅研究了0.2～1.5 THz波段的吸收光谱。图4(a)～(c)为来自于汉中、江油、云南3种产地附片各70组太赫兹吸收光谱；在与附片相同的检测系统及条件下检测其余6种药材，得到的60组太赫兹原始吸收光谱如图4(d)所示，将其与部分附片光谱数据作为类间测试集的样本；从图中可以看出，仅靠吸收系数的大小难以区分不同产地的附片药材，而使用机器学习可能是实现中药材产地鉴别的关键。

图4 样品的太赫兹吸收光谱

3.2 建立识别模型

3.2.1 识别模型的建立与对比分析

本文将所有附片70 %的太赫兹吸收光谱数据作为输入量输入到SVM识别模型、RF识别模型、PLS-DA中进行训练。各分类器最佳参数值如表1所示。

表1 3种分类器的最佳参数

根据表2可绘制出如图5所示的3种模型在不同阈值下的识别准确率，从其中可以看出,PLS-DA结合FDA(PLS-DA-FDA)模型的准确率在任何阈值下均是最高的。

表2 3种识别模型的类内测试实验结果

图5 3种模型在不同阈值下的识别准确率

将类内测试集(即未知产地的附片)的63个测试样本输入识别模型，同时将从附片中随机抽取10个样本(3汉中+3江油+4云南)和其余6种药材各10组光谱组成类间测试集输入3种识别模型。图6中类内测试集对应FAR曲线，类间测试集对应FRR曲线，同一识别模型的2条曲线的交点称作EER点，将其所对应的值设定为识别模型在确认分类结果时的阈值。从图中可以看出，PLS-DA-FDA模型的EER最低约2 %，优于其他模型，此时阈值为0.16。

图6 FAR/FRR—阈值曲线

3.2.2 PLS-DA-FDA模型的实验结果

将分类模型的阈值设置为0.16。从表3中可以看出,该模型对未知产地的附片平均准确率为95.23 %,训练集准确率为95.89 %。

表3 PLS-DA-FDA模型的类内测试实验结果

PLS-DA-FDA模型进行类间测试，从附片的每种产地中随机抽取共计10个样本，与其余6种药材各10个测试样本，共计70个样本组成类间测试集，并将其输入模型进行验证，其结果如图7所示，其中“☆”表示真实标签，“●”表示预测标签，从图中可看出，仅有2个其他品种的药材被误识为汉中制附片，未知药材名称时的识别准确率为97.14 %。

图7 PLS-DA-FDA模型在类间测试集上的识别结果

4 结 论

用6种非附片的其他药材和附片的太赫兹光谱数据来检验模型对类间样本的识别，用附片的数据验证模型对3种产地的鉴别能力。结果表明：太赫兹光谱结合PLS-DA和模糊判别方法是实现附片产地快速识别的有效方法。