基于FTIR技术和稀疏线性判别分析的秦艽种类鉴别
李四海1, 余晓晖2, 赵磊2, 晋玲2,*
1. 甘肃中医药大学信息工程学院, 甘肃 兰州 730000
2. 甘肃省高校中(藏)药化学与质量研究省级重点实验室, 甘肃 兰州 730000
*通讯联系人 e-mail: zyxyjl@163.com

作者简介: 李四海, 1972年生, 甘肃中医药大学信息工程学院副教授 e-mail: lshroom@163.com

摘要

傅里叶变换红外光谱通常包含有大量的波长变量点, 对其进行定性分析需要建立稳健的、 可解释性的分类模型。 稀疏线性判别分析(SLDA)是一种较为新颖和有效的机器学习算法, 常用于高维度、 小样本数据的变量筛选和判别分析, SLDA通过在线性判别分析中引入正则项, 使分类器训练过程和变量选择过程同时完成, 不同判别方向上载荷系数的稀疏性则增强了模型的可解释性。 采集甘肃不同产地的秦艽样本94个, 其中麻花秦艽( Gentiana straminea Maxim)30个, 黄管秦艽( Gentiana officinalis)28个, 大叶秦艽( Gentiana macrophylla Pall)36个, 利用傅里叶变换红外光谱法获得所有样本的光谱图。 取其中70个样本构成训练集, 剩余24个为测试集。 使用训练集建立SLDA模型, 对2个判别方向上不为0的载荷系数个数进行网格化寻优, 得到了最优的参数空间。 利用建立的SLDA模型对测试样本进行预测, 其分类准确率达到100%, 实现了对三种秦艽的快速、 准确鉴别。 实验结果表明, 与PLS-DA方法相比, SLDA模型在分类准确率、 稀疏性及可解释性方面均具有一定优势, 是一种新颖、 有效的光谱定性分析方法。

关键词: 秦艽; 傅里叶变换红外光谱; 正则化; 稀疏线性判别分析; 变量选择
中图分类号:TP391.4 文献标识码:A
Identification of Gentiana Macrophylla by FTIR Technology and Sparse Linear Discriminant Analysis
LI Si-hai1, YU Xiao-hui2, ZHAO Lei2, JIN Ling2,*
1. School of Information Science & Engineering, Gansu University of Traditional Chinese Medicine, Lanzhou 730000, China;
2. Key Laboratory of Chemistry and Quality for Traditional Chinese Medicines of the College of Gansu Province, Lanzhou 730000, China
Abstract

Fourier transform infrared(FTIR) spectrum usually includes a large number of wavelength variables and the qualitative analysis of FTIR spectrum needs to establish a stable and interpretable classification model. Sparse linear discriminant analysis (SLDA), a relatively new and effective machine learning algorithm, is commonly used for variable selection and discriminant analysis of high-dimensional settings, in which the number of wavelength variable is very large and the number of observations is limited. By introducing regularization items into linear discriminant analysis, the classifier training and variable selection are performed simultaneously in SLDA, and the sparsity of load coefficients in different discriminant directions increases the interpretability of the model. A total of 94 samples of Gentiana macrophylla, including 30 Gentiana straminea Maxims, 28 Gentiana officinalis and 36 Gentiana macropylla Pall, were collected. FTIR spectrum of all samples was obtained by Fourier transform infrared spectroscopy method. 70 of the samples were selected as the training set, the remaining as the test set. Based on the training set, the SLDA model was established through the grid optimization of the number of non-zero loading coefficients in the two discriminant directions, and the optimal parameter space was obtained. According to the model parameters, the prediction accuracy of the test set was 100%, and thus the rapid and accurate identification of the three kinds of Gentiana macrophylla was realized. The experimental results showed that the SLDA model was superior to PLS-DA method in terms of classification accuracy, sparseness and interpretability. SLDA will be a novel and effective method for spectroscopy qualitative analysis.

Key words: Gentiana macrophylla; FTIR; Regularization; Sparse linear discriminant analysis; Variable selection
引言

秦艽(Radix Gentiana macrophylla)是龙胆科龙胆属秦艽组植物, 为我国常用大宗中药材, 具有较高的药用价值, 是治疗风湿关节痛、 结核病、 潮热、 黄疸等症的主药之一。 秦艽组植物全世界约20种, 我国至少有17种, 其中大叶秦艽分布最广, 甘肃、 陕西、 山西和四川是秦艽的道地产区[1]

目前, 秦艽主要依靠野生资源满足用药需求, 但长期的滥采滥挖导致了秦艽资源的严重匮乏, 一些伪品秦艽也充斥于药材市场, 导致市场用药品种混乱。 因此, 开展秦艽种类、 产地及真伪鉴别对于保证秦艽质量及临床疗效具有重要意义。 秦艽的鉴定方法主要包括传统生药学鉴定和现代鉴定方法两大类[2]。 传统生药学鉴定是最常用的方法, 但该类方法存在主观性强、 周期长, 成本高、 无法对中药材质量进行整体评价等问题。 在现代鉴定方法中, 运用傅里叶变换红外光谱技术结合化学计量学方法对秦艽种类进行定性分析, 具有简单、 快速、 无破坏性、 可综合评价等优点。

在FTIR光谱的定性分析中, 应用最为广泛的是主成分分析(PCA)及偏最小二乘判别分析(PLS-DA)。 De Luca M[3]等利用FTIR技术结合PLS-DA对摩洛哥四种产地的橄榄油进行了聚类分析。 吴喆[4]等对不同采收期滇重楼的FTIR光谱进行平滑及小波去噪预处理, 然后结合化学计量学中的偏最小二乘判别分析方法对不同年限的滇重楼进行了聚类分析。 偏最小二乘判别分析的主要缺点在于每个主成分都是原始波长变量的线性组合, 模型中包含的变量过多, 一方面会使模型过于复杂, 泛化能力下降; 另一方面, 由于绝大多数波长变量的载荷系数不为0, 导致无法有效判断波长变量的重要性, 模型的可解释性不强。

稀疏判别分析是在Fisher’ s线性判别分析的基础上, 通过在目标函数中引入L1和L2正则项, 使模型既具有稀疏性又有良好的稳健性。 该方法的优点在于建立的模型简单, 能够自动筛选出重要变量, 模型判别准确率高。 稀疏判别分析已被成功应用于基因微阵列数据分析[5, 6], 但在FTIR光谱判别分析方面还未见报道。

本文利用傅里叶变换红外光谱仪获得甘肃产三种秦艽的FTIR光谱, 建立稀疏线性判别分析模型, 筛选出重要的波长变量点, 实现三种秦艽的快速、 准确鉴别, 为说明本方法的有效性, 同时建立了三种秦艽的偏最小二乘判别分析模型, 对两种模型的分类准确率及模型性能进行了对比分析。

1 稀疏线性判别分析
1.1 原理

假设X为原始光谱矩阵, XRn× p, pn。 共n个样本, p个波长变量; yn× 1为类标签向量, y∈ {1, 2, …, K}, K为类别数。 Ck为第k类样本的集合, k=1, 2, …, K。 对X进行中心化使其每一列均满足: 均值为0、 标准差为1。 将类别向量y编码为n× K稀疏矩阵形式

Y={Yik=1|iCk, i=1, 2, , n; k=1, 2, , K}

Fisher’ s线性判别分析的基本原理是找到一系列相互正交的判别方向, 使得原始数据在各判别方向上的类间方差最大、 类内方差最小。

引入最优得分矩阵θ k, 将Y投影到相互正交的各个判别方向将分类问题转化为回归问题, 对第k个判别方向, 问题转化为求解如下的优化问题[7]

minβk, θk1nYθk-Xβk2s.t 1nθTkYTYθk=1; θTkYTYθl=0 l< k(1)

式(1)中, θ kK维列向量, β kp维列向量。

为获得模型的稀疏解, 对式(1)添加L1正则项

minβk, θk1nYθk-Xβk2+λβk1s.t 1nθTkYTYθk=1; θTkYTYθl=0 l< k(2)

式(2)为lasso线性回归问题, ‖ β k1为所有回归系数的绝对值之和, λ 为L1正则项的系数。 通过施加L1正则可将模型的解限定在一个L1-ball中[8], 使大多数变量在判别方向上的载荷系数为0, 有利于获得稀疏解, 使模型具有更好的可解释性。 但L1正则也存在如下问题: (1)选择波数变量的上界为min(n, p)。 (2)不能有效解决变量选择中的组效应问题: 如果一个光谱区间的多个波数变量具有同等重要程度, lasso也只从中选择一个变量, 即lasso不能有效解决变量选择的组效应问题。

考虑到FTIR光谱中变量数p远远大于样本数n, 为有效解决高维光谱数据的变量选择问题, 在式(2)中继续添加L2正则项

minβk, θk1nYθk-Xβk2+γβTkΩβk+λβk1s.t 1nθTkYTYθk=1, θTkYTYθl=0 l< k(3)

式(3)中, γ 为非负惩罚系数, Ω 为正定矩阵, 当Ω =I时, 式(3)为弹性网回归问题[9]

式(3)可通过交替迭代算法求解: 首先固定θ k, 求解弹性网回归问题得到β k; 然后固定β k, 更新θ k

对于固定的β k, 令Dπ =YTY/n , 则式(3)变为式(4)

minθk{Yθk-Xβk2}s.t θTkDπθk=1,  θTkDπθl=0l< k(4)

θ k的无约束估计为 Dπ-1YTk, 考虑到与前k-1个判别方向的正交性, 故θ kI-Qk QTkDπ 投影算子上的投影即为需要计算的θ k

1.2 稀疏线性判别分析算法

稀疏线性判别分析算法寻找q=K-1个相互正交的判别方向并返回判别方向上的系数矩阵β p× q

Step1: 将类别向量y编码为n× K矩阵形式

Yn×K={Yik=1|iCk, i=1, 2, , n; k=1, 2, , K}

Step2: Dπ =YTY/n

Step3: 初始化判别方向上的类别得分矩阵θ k及回归系数矩阵β p× q, θ k的列向量为相互正交的单位向量。

Step4: for k=1:q, 使用交替迭代算法求解(θ k, β k)。

(1) 固定θ k, 使用最小角回归算法[10, 11]求解如下的弹性网回归问题, 计算β k

minβk1nYθk-Xβk2+γβTkΩβk+λβk1

(2) 固定β k, 更新θ k

Qk=θ k(:, 1: q-1), 即Qk为前q-1个判别方向上的类别得分矩阵, 则

θ˙k=(I-QkQTkDπ)Dπ-1YTXβk

θk=θ˙k/θ˙TkDπθ˙k

(3)重复以上两步直至β k收敛或达到最大迭代次数。

Step5: 将原始光谱数据X投影到q个判别方向: Xn× q=(1, 2, …, q), 进行LDA线性判别分析。

2 实验部分
2.1 仪器与试药

ALPHA-T傅里叶变换红外光谱仪(配有DTGS检测器, 德国Bruker公司), 溴化钾(美国PIKE公司)等试剂均为分析纯。 94个秦艽样本分别采自甘肃古浪、 天祝、 卓尼等地, 药材经甘肃中医药大学中药资源教研室晋玲教授鉴定。

2.2 样品MIR指纹图谱采集

将各秦艽样品干燥粉碎, 过120目筛, 压片前放入干燥器中放置48 h。 取样品粉末约10 mg, 置于玛瑙研钵中, 与约100 mg的KBr粉末混合, 研磨均匀后压片, 制成厚度为1 mm透明晶片作为供试品。 将供试品放入ALPHA-T型傅里叶变换红外光谱仪, 以DTGS检测器在4 000~400 cm-1光谱范围内扫描, 分辨率为4 cm-1, 扫描信号累加次数为16次, 求取平均光谱。 每批样品重复检测3次。 扫描实时扣除H2O和CO2的干扰, 获得一维FTIR谱图。 傅里叶变换红外光谱如图1所示。

图1 样本的傅里叶变换红外光谱Fig.1 FTIR spectrum of all samples

2.3 数据预处理

94个秦艽样本中麻花秦艽30个, 黄管秦艽28个, 大叶秦艽36个, 原始光谱矩阵大小为94× 2 558。

对原始光谱数据进行如下预处理: (1) 为消除基线漂移及背景噪声的影响, 首先进行基线校正及平滑预处理, 采用Savitsky-Golay平滑滤波, 窗口大小为13, 多项式阶数为3。 (2) 将每个光谱波长变量进行中心化处理, 使其均值为0, 方差为1。

2.4 模型的建立

将样本划分为训练集和测试集。 训练集中麻花秦艽23个, 黄管秦艽21个, 大叶秦艽26个, 共70个样本; 测试集中麻花秦艽7个, 黄管秦艽7个, 大叶秦艽10个, 共24个样本。

稀疏线性判别模型中需要设置的参数包括: L2正则项系数delta、 判别方向个数及不同判别方向上载荷系数值不为0的波长变量个数、 收敛误差、 迭代次数。 其中delta用于权衡模型的稀疏性和分类能力之间的相对重要程度, 对模型性能影响较大。 delta值越大, 模型越稀疏, 波长变量载荷系数也越大, 但模型的分类能力会有所下降, 一般设置delta为10-5时会在二者之间取得较好的平衡。

实验过程中设置delta为10-5, 收敛误差为10-5, 迭代100次, 判别方向个数为2, 各判别方向上载荷系数值不为0的波长变量个数则通过网格搜索法进行寻优确定。

稀疏线性判别分析算法通过交替迭代, 求解出2个判别方向上所有波长变量的载荷系数矩阵β , 大小为2 558× 2, 然后将经过预处理后的光谱投影到2个判别方向, 将光谱的维度由2 558维降为2维, 有效实现了高维光谱数据的低维表达, 最后对该低维空间数据建立线性判别分析模型实现对三种秦艽的鉴别。 算法实现平台为MATLAB R2016b。 为说明本方法的有效性, 同时建立了偏最小二乘判别分析模型, 对两种模型的性能进行对比分析, PLS-DA算法在SIMCA-P软件中实现。

3 结果与讨论
3.1 参数网格化寻优

为评估稀疏性对模型分类性能的影响, 对两个判别方向上非0载荷系数的个数进行网格化寻优, 实验时设置的搜索范围分别为[20, 50], [20, 50], 步长均为1, 生成网格化参数数据分别进行稀疏线性判别分析, 模型在训练集上的最优参数空间如图2所示。

图2 非0载荷系数个数与分类准确率的关系Fig.2 Relationship between classification accuracy and the number of non-zero loading coefficients

从图2可得, 模型在训练集上的分类准确率较高, 最低为92.81%。 随着第1判别方向上非0载荷系数个数增加, 分类准确率呈阶梯状稳步上升, 当非0载荷系数增加到45个以上时, 分类准确率能稳定保持在100%。 图2同时表明, 模型分类准确率受第1判别方向上非0载荷系数个数的影响较大, 主要原因在于稀疏线性判别分析算法将光谱数据投影到不同的判别方向, 而第1判别方向是方差最大的方向, 所携带的原始光谱信息最多, 第1判别方向上非0系数过少, 则不能很好地提取光谱的重要特征。 仿真实验表明, 第1判别方向上非0载荷系数的个数应取较大值。

3.2 分类结果

根据训练集上参数的寻优结果, 最终确定2个判别方向上非0载荷系数个数分别取50和30, 其他参数设置如2.4, 建立稀疏线性判别分析模型。

70个训练样本在2个判别方向的投影散点图如图3所示。

图3 训练样本在2个判别方向上的投影Fig.3 Projection of the training samples onto the two discriminant directions

从图3可以看出, 三种秦艽样本之间的线性边界非常明显, 表明算法对光谱特征的提取非常有效, 其中麻花秦艽与其他两类秦艽的差异较大。

进一步考察模型在测试集上的预测能力。 根据模型在训练集上所得到的2个判别方向的载荷系数矩阵, 对24个测试样本进行投影, 得到测试集的投影散点图, 如图4所示。 可以看出, 三种秦艽在平面上的投影差异明显, 相互之间的线性边界很清晰, 非常易于分类。 选择最简单的线性判别分类器, 对投影得到的24× 2数据矩阵进行判别分析, 最终得到模型在测试集上的分类准确率为100%。

图4 测试样本在2个判别方向上的投影Fig.4 Projection of the test samples onto the two discriminant directions

为进一步说明模型的稀疏效果, 图5显示了原始光谱中2 558个波长变量点在2个判别方向上的载荷系数。

图5 2个判别方向上的载荷系数Fig.5 Loading coefficients onto the two discriminant directions

由图5可得: (1) 光谱的2 558个波长变量中, 在2个判别方向上具有非0载荷值的变量数分别为50和30, 这表明所建立的稀疏判别模型实现了光谱重要变量的自动选择。 (2) 稀疏模型有利于实现高维数据的低维可视化, 具有更好

的直观性。 (3) 2个判别方向上具有较大载荷系数值的波数分别为: 558.43, 619.38和955.29 cm-1, 这说明三种秦艽在上述光谱波数处具有明显差异。

3.3 与PLS-DA模型的对比分析

为进一步验证本方法的有效性, 同时建立三种秦艽的偏最小二乘判别分析模型, 对两种模型的性能进行对比分析。 使用相同的训练集和测试集, 对原始光谱进行一阶导数预处理, 建立PLS-DA判别模型, 取前10个主成分。 两种模型对24个测试样本的预测结果如表1所示。

表1 两种模型的性能对比 Table 1 Performance comparison of SLDA and PLS-DA

表1可以看出, PLS-DA模型提取前10个主成分时, 分类准确率为87.5%, 明显低于SLDA模型。 在载荷系数的稀疏性方面, SLDA模型通过设置波长变量点在不同判别方向上不为0的系数个数从而获得模型的稀疏解, PLS-DA模型波长变量点在不同主成分上的载荷系数不具有稀疏性, 这说明SLDA模型的可解释性要优于PLS-DA模型。

4 结 论

利用FTIR技术结合稀疏线性判别分析, 实现了三种秦艽的准确鉴别。 稀疏线性判别分析算法中的关键参数包括: L2正则项系数及不同判别方向上不为0的载荷系数个数。 通过网格化寻优, 得到了最优参数空间, 实验表明, 第1判别方向上非0载荷系数个数对模型性能影响较大。 同时建立了三种秦艽的PLS-DA判别分析模型, 仿真结果表明, SLDA模型在分类准确率、 低维可视化及稀疏性方面均优于PLS-DA模型。

光谱定性分析需要运用化学计量学及机器学习中的模式识别方法建立相应的分类模型, 模型的稳健性是评价模型性能的一个重要指标。 研究表明, 进行光谱波长变量筛选是提高模型稳健性的一个有效方法。 近年来, 在光谱的定性及定量分析方面, 出现了一些基于蒙特卡洛随机搜索策略的波长选择方法, 如竞争自适应重采样、 随机蛙跳算法等, 其目的也是通过选择重要变量, 建立更简单、 更稳健、 具有更好可解释性的分类模型。 本文提出的稀疏线性判别模型为光谱的定性分析提供了新的思路和方法。

The authors have declared that no competing interests exist.

参考文献
[1] HUANG Lu-lin, YANG Xiao, FENG Xian-hong, et al(黄璐琳, 杨晓, 丰先红, ). Modern Chinese Medicine(中国现代中药), 2011, 13(5): 40. [本文引用:1]
[2] XIONG Bo, GUO Shu-peng, HU Lin(熊波, 郭树鹏, 胡林). Chinese Journal of Experimental Traditional Medical Formulae(中国实验方剂学杂志), 2015, 21(17): 230. [本文引用:1]
[3] De Luca M, Terouzi W, Ioele G, et al. Food Chemistry, 2011, 124(3): 1113. [本文引用:1]
[4] WU Zhe, ZHANG Ji, JIN Hang, et al(吴喆, 张霁, 金航, ). Spectroscopy and Spectral Analysis(光谱学与光谱分析), 2017, 37(6): 1754. [本文引用:1]
[5] Shao J, Wang Y, Deng X, et al. The Annals of Statistics, 2011, 39(2): 1241. [本文引用:1]
[6] Kong H, Lai Z, Wang X, et al. Neurocomputing, 2016, 177: 198. [本文引用:1]
[7] Clemmensen L, Hastie T, Witten D, et al. Technometrics, 2011, 53(4): 406. [本文引用:1]
[8] Rasmussen M A, Bro R. Chemometrics and Intelligent Laboratory Systems, 2012, 119: 21. [本文引用:1]
[9] Zou H, Hastie T. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 2005, 67(2): 301. [本文引用:1]
[10] Efron B, Hastie T, Johnstone I, et al. The Annals of Statistics, 2004, 32(2): 407. [本文引用:1]
[11] Liu C, Yang S X, Deng L. Expert Systems with Applications, 2015, 42(22): 8497. [本文引用:1]