基于中红外光谱数据对中药材的鉴别分析

Identification of Chinese Medicinal Materials by Using Mid-Infrared Spectroscopy Data

摘要: 文章首先选取中红外光谱数据差异较为明显的波数区域，并使用K-means聚类分析法对其进行分析，从而鉴别出药材种类；再次选取同种药材在不同产地的中红外光谱波峰处数值，对这些数据针对产地和需鉴定药材种类进行均值处理，并对相应波段的差值绝对值的方差进行估计，方差最小的则确定为目标药材的产地；最后，采用暴力枚举法对数据进行分析，并利用方差估计进一步进行判断，得出对应种类和产地。

Abstract: In this paper, the wavenumber regions with obvious differences in mid-infrared spectral data were firstly selected and analyzed by k-means clustering analysis, so as to identify the medicinal materials. The values at the mid-infrared spectrum peaks of the same medicinal materials in different producing areas were selected again, and these data were averaged according to the origin and the species of medicinal materials to be identified. The variance of the absolute difference values of the corresponding bands was estimated, and the one with the smallest variance was determined as the origin of the target medicinal materials. Finally, the violence enumeration method was used to analyze the data, and variance estimation was used to further judge the corresponding species and origin.

1. 问题背景

本文以2021年高教社杯全国大学生数学建模竞赛E题为依托进行分析。中药材的种类及产地对中药材的药性有很大影响，因此，对中药材的种类及产地的判定十分重要。不同中药材表现的光谱特征差异较大，即使来自不同产地的同一药材，因其无机元素的化学成分、有机物等存在的差异性，在近红外、中红外光谱的照射下也会表现出不同的光谱特征，因此可以利用这些特征来鉴别中药材的种类及产地。本文主要目标是通过中红外光谱数据进行数据分析，从而对药材的种类和产地进行鉴别分析 [1]。对所给药材的中红外光谱数据进行数据处理，筛出异常数据，选取波峰差异明显区域，利用K-means聚类分析法，得到药材分类结果；接下来，根据已知药材种类数据，并依据相应波段的差值绝对值方差进行估计分析，得出目标药材的种类判定结果 [2] [3] [4]；最后，依据给定种类和产地的中药材的近红外光谱数据，利用暴力枚举法判定未知药材的种类及产地，并根据方差进行进一步精细判定。

2. 药材种类的鉴别模型

2.1. 异常数据处理

题中附件1给出了所需判断种类的药材中红外数据，由于所给数据量较大，因此对数据在EXCEL进行排序进行处理后发现NO64、NO136、NO201这三个编号数据与整体数据差异较大，因此本文首先将该异常数据清洗掉。

2.2. 模型建立与分析

选择波数区域(652-662)、(1144-1154)、(3242-3252)数据差异较为明显三个波段的数据进行均值处理，建立欧式距离分析处理较为合适使用K-means聚类分析法对三个区域的数据进行分析。

建立欧式距离的表达式如下：

$d_{ij}=\sqrt{{\left(x_i-x_j\right)}^2+{\left(y_i-y_j\right)}^2}$ (1)

$\left(i,j=652,\cdots ,662;1144,\cdots ,1154;3242,\cdots ,3252\right)$

其中： $d_{ij}$ 表示第i个样本和第j个样本之间的距离。

按照吸光性由低到高分为三类，利用K-means聚类分析法，对上述三个区域的数据进行迭代处理，得到10组结果。因此，分析得出所给药材数据中有10种药材。

3. 同种药材的不同产地的鉴别模型

3.1. 数据分析

根据题中提供的数据显示，同种药材在不同产地的中红外光谱数据曲线图变化趋势相同，但变化量有所差异。根据题目中给出的数据按照产地构建中红外光谱数据曲线图，如下如图1。

注：该图为产地1的药材中红外光谱数据图，其他产地同理，故不列入正文。

Figure 1. Mid-infrared spectrum data of crude drugs from place 1

图1. 产地1的药材中红外光谱数据图

其次，对中红外光谱数据曲线图选取波峰处数据作为对比参照，观察图中波峰区间进行分组，选取6组波峰，分别为(551-650)、(1011-1110)、(1371-1470)、(1601-1700)、(2901-3000)、(3251-3350)。因此，本考虑对已知药材的中红外光谱图的数据与其不同产地药材的中红外光谱数据的方差数据进行分析，从而鉴别同种药材的不同产地。

3.2. 模型建立与分析

考虑对已知药材的中红外光谱图的数据与其不同产地药材的中红外光谱数据的方差数据进行分析，从而鉴别同种药材的不同产地。

分别对该区域内的数据按照产地和种类分别进行均值处理，首先按照产地分类进行均值处理(表1)：

$\overline{O{P}_{n\partial k}}=\frac{{\displaystyle \underset{m=1}{\overset{100}{\sum }}O{P}_{n\partial km}}}{100}$ (2)

$\left(\partial =1,2,\cdots ,6;k=550,1000,1370,1600,2900,3250;n=1,2,\cdots ,11\right)$

$\partial$ ：波峰区间分组；

$\overline{O{P}_{n\partial k}}$ ：针对第n个产地的第 $\partial$ 组的均值处理结果；

$O{P}_{n\partial km}$ ：针对第n个产地的第 $\partial$ 组的第 $k+m$ 个数据。

其次，按照产地分类进行均值处理：

$\overline{N{O}_{\partial \omega k}}=\frac{{\displaystyle \underset{m=1}{\overset{100}{\sum }}N{O}_{\partial \omega km}}}{100}$ (3)

$\left(\omega =3,14,38,48,58,71,79,86,89,110,134,152,227,331,618\right)$

$\overline{N{O}_{\partial \omega k}}$ ：针对第 $\omega$ 种药材的第 $\partial$ 组的均值处理结果；

$N{O}_{\partial \omega km}$ ：针对第 $\omega$ 种药材的的第 $\partial$ 组的每个数据。

随后，为了估计第 $\omega$ 个未知产地药材与第n个产地的逼近效果，对相应波段的差值绝对值的方差进行估计：

相应波段的差值的绝对值：

$D_{\partial \omega n}=\left|\overline{N{O}_{\partial \omega k}}-\overline{O{P}_{n\partial k}}\right|$ (4)

相应波段的差值的平均值：

$\overline{D_{\omega n}}=\frac{{\displaystyle \underset{\partial =1}{\overset{6}{\sum }}{D}_{\partial \omega n}}}{6}$ (5)

相应波段的差值的方差：

$S{D}_{\omega n}^2=\frac{{\displaystyle \underset{\partial =1}{\overset{6}{\sum }}\left(D_{\partial \omega n}-\overline{D_{\omega n}}\right)}}{6}$ (6)

Table 1. Data variance processing results for 11 producing areas

表1. 针对11个产地的数据方差处理结果

最终判断下列序号对应的未知产地药材的产地结果如下表2：

Table 2. Origin of 15 kinds of medicinal materials

表2. 15种药材产地

4. 未知药材的不同产地的鉴别模型

4.1. 模型建立

以给定的A、B、C三种种类的药材的近红外光谱数据作为参考依据，判定未知药材的种类和产地。利用暴力枚举法在C++中建立循环遍历模型，首先筛选出同一光数处的A、B、C三种种类的药材数据，选取药材种类所在列作为主要关键字，进行循环遍历，得到药材分类，其次选择药材产地所在列作为次要关键字再次进行循环遍历，最终得到目标药材的种类和产地。算法程序如下图2：

Figure 2. Flow chart of violence enumeration method

图2. 暴力枚举法流程图

根据问题要求，根据给出的几种药材的近红外光谱数据判断出其种类和产地。首先以药材种类所在列为主要关键字，药材产地所在列为次要关键字，对数据进行排列，并利用暴力枚举法对数据进行分析。

4.2. 药材种类分析及结果

首先以药材种类作为分类，在C++中进行循环遍历，如果目标在某一个光数处不满足数据所给出的范围，那么返回false，如果满足所给范围，则循环至下一个光数处，直至遍历整个数组。若一直能满足数组范围，说明目标在每一处光数处都能满足此种类的吸光度，所以返回true，证明此目标满此子种类。得到以下布尔型值(表3)：

Table 3. Values of Boolean types

表3. 布尔型的值

即，A种类药材为NO94、109、140；B种类药材有NO278、308、330；C种类药材为NO347。

4.3. 药材产地分析及结果

在同种药材类别下对产地进行分类，按照附件所给数据，A种类共分为5个产地，为A_a $\left(a=1,\cdots ,5\right)$ ；B种类共分为16个产地，为B_b $\left(b=1,\cdots ,16\right)$ ；C种类共分为4个产地，为C_c $\left(c=1,\cdots ,4\right)$。将目标药材循环遍历数组，根据程序运行返回false或者true，得到产地初步分类(表4)：

Table 4. Preliminary classification of origin

表4. 产地初步分类

由于NO94、109、140、347出现多个产地的判定结果，因此为了精确判断出产地，求出该目标在每一处光数处与数组中相应产地的中位数的方差，最小方差则为对应产地，得到药材最后的产地分类(表5)：

Table 5. Final classification of origin

表5. 产地最终分类

5. 总结

利用K-means聚类分析法对数据进行迭代处理，可以按照特征对数据进行分类，同时利用方差公式对种类进行吻合分析，最后利用暴力枚举法以及方差估计进一步得到种类和产地的最优匹配结果。在对未知药材的种类进行判定时，对相关数据进行了均值简化处理，虽不失精确性，但在接下来的改进中我们可以对给出的全体样本数据直接进行方差处理，这样能够更加精确地判断出产地和种类。

基金项目

2021年度安徽高校自然科学基金(KJ2021A1447)。