CN101498661B - 高精度分辨中药材品种、产地及生长方式的红外光谱特征提取方法 - Google Patents

Abstract

本发明披露了一种高精度分辨中药材品种、产地及生长方式的红外光谱特征提取方法。此方法取得用于训练的每个品种的中药材样本的红外光谱后再将其预处理，根据产地及生长方式将样本分类，并计算出各品种的类间类内差异比谱函数，借以取得各样本的特征向量，再将这些特征向量降维及进行分布优化，从而取得类间距离大而类内距离小的训练样本特征向量集。在该特征向量集形成过程中所获得的类间类内差异比谱函数以及降维矩阵和分布优化变换矩阵即可被用于送检样本的红外光谱特征提取，从而使该送检样本得到自动识别。

Description

高精度分辨中药材品种、产地及生长方式的红外光谱特征提取方法

技术领域

本发明属于对中药材品种、产地及生长方式的自动识别方法的技术领域。

背景技术

中药是中华民族在悠久的历史进程中与疾病斗争的科学结晶和向人类贡献的璀璨瑰宝，疗效突出。然而由于当代全球医疗科技的迅猛发展，中药的传统优势地位正在受到冲击。目前我国中药市场上存在着的一些质量问题，已对中医药的信誉、疗效和人们的生命健康带来了不良的后果。

中药材的质量是中药质量的基础，中药材质量鉴别在中药疗效、安全及发展中扮演着极其重要的作用。目前，中药材的鉴别主要有二大方法：(1)基于经验的性状鉴别，及(2)理化鉴别方法。前者主要直接透过观察药材的形状、颜色、气味、尺寸及剖面而鉴定其质量，因此对观察者要求极其严格，必须经过长期的经验积累。但即使这样，此方法人为误差依然较大。相比之下，理化鉴别以分离或以微区分析为基础，如色谱指纹图谱的应用对于解决质量的整体控制起到了一定的作用。但这种技术相对复杂、费时、成本高，而且属于灵敏、精细的微观方法，重复性较差。因此，这种方法不适合从宏观角度，定性的快速鉴定中药材质量。

中药材的质量决定于所含化学成分和各成分的比例，因而与它的产地及生长方式有着密不可分的关系。通过找出中药材的真实来源地及生长方式，从宏观上为中药材质量鉴别提供了一条快速而简易的可行方法。中药材种类、产地及生长方式的不同，将会造成中药材红外光谱的差异。我们知道，用红外光谱鉴定化合物，其优点是简便、迅速和可靠；同时样品用量少、可回收；对样品也无特殊要求，无论气体、液体和固体均可进行检测。而自从20世纪70年代后期傅里叶变换技术引入红外光谱领域以来，化合物红外光谱的获得已变得相当容易、快速和精确，这就为通过红外光谱鉴定中药材质量提供了一个崭新的平台。将红外光谱法与计算机和模式识别技术相结合，将使中药材质量鉴定手段产生巨大的飞跃。根据模式识别原理，通过红外频谱识别中药材的品种、产地及生长方式的过程可分为训练及判别二个基本阶段。在训练阶段中，先利用红外光谱仪获取相当数量的已知真实产地和生长方式的中药材样本的红外光谱图，再利用信号处理技术对各光谱图进行预处理、特征提取和选择。经过训练阶段，这些中药材样本的红外光谱原始测量数据映像为高维空间中的特征向量，不同种类、产地和生长方式的中药材样本的特征向量应形成不同的集群。在判别阶段中，送检的中药材红外光谱进行一系列与训练阶段一致的处理，被映像到训练样本集的特征空间中，再用分类算法判别其种类、产地和生长方式。上述方法的核心关键是如何从中药材红外光谱图提取有效特征，更具体地说，是从红外光谱的哪些频率(波数)点处对光谱幅度进行抽样用以提取有效特征。目前已有一些科学工作者在这方面作了有益的探索，例如凭药物理论知识和经验在某些特征峰所在处提取特征，或者通过贪婪搜索法搜索特征谱段提取特征。这些方法有下列局限：需要有关人员具有较高的药物知识水平和丰富的经验积累，效率不高而缺少普遍适用性，不便于推广到品种、产地及生长方式较多情况下的中药材鉴别。

为了使基于红外频谱的对多种中药材的品种、产地及生长方式的自动识别技术尽早推广应用，异常迫切需要寻找出一个既具有普遍适用性，又能服务于高精度鉴别中药材品种、产地及生长方式的红外光谱特征提取方法。

发明内容

技术问题：

本发明的目的是给出一个既具有普遍适用性，又能服务于高精度鉴别多种中药材品种、产地及生长方式的红外光谱特征提取方法，建立一个操作方便、显示内容丰富、有利于开发新算法的红外光谱识别系统。为了使初始特征提取有利于红外光谱分类，本发明提出了一个能综合反映红外光谱中各频率(波数)点处的类间样本差异与类内样本差异的“类间类内差异比”谱函数(也称为VoM/MoV函数)。为了有利于在特征空间中采用简单而高效的算法判别送检样本的类别，本发明提出了一个实现训练样本特征向量在多维空间中分布最优化的迭代算法。

技术方案：

本发明以下列条件作基础：对于要求鉴别的每一种中药材(例如丹参、人参等)，已经获取了用傅里叶变换红外光谱仪测得的、已知产地和生长方式的相当数量样本的红外光谱。

对于要求鉴别的每一种中药材，根据获取的红外光谱的特点进行数据的预处理(包括归一化，去噪，求导等)，然后求取类间类内差异比谱函数(VoM/MoV函数)；根据VoM/MoV函数指示的位置，从红外光谱的指纹区抽取经过预处理的波形幅度值作为初始特征；再通过K-L变换进行维数压缩并通过一种迭代算法实现训练样本特征向量在多维空间中的分布最优化。对送检的中药材红外光谱，进行一系列与每一种中药材训练阶段一致的处理，使之映像到训练样本集的特征空间中，再用分类算法判别该中药材的种类、产地和生长方式。为了操作方便、显示内容丰富、有利于开发新算法，本发明专门研制一套由Visual C++语言编写的软件系统。该软件系统不但使红外光谱数据文件的调用、转换和处理非常容易，而且具有科学计算可视化的功能，可以将每一步骤的结果以彩色图形的方式及时显示。

本发明提供一个既具有普遍适用性，又能高精度分辨中药材品种、产地及生长方式的红外光谱特征提取方法。可以高精度、快速且同时检测多品种中药材的产地及生长方式，从而可从宏观上定性检测药材的质量及真假。

基于所发明的红外光谱特征提取方法，可将在训练步骤中的分类结果存于远程服务器数据库，并对受测药材光谱特征进行网上查询，可开发出基于Web的快速及有效的中药材品种、产地及生长方式在线检测系统。

给定相当数量的已知质量级数的中药材红外光谱图，本发明可快速测定受检未知药材的真假及质量级数。

本发明可用于对食品及饮品类的真、劣、假进行快速准确的检测。

附图说明

图1是丹参的训练集样本的产地及生长方式显示。

图2是丹参的测试集样本的产地及生长方式显示。

图3是丹参训练集样本在30维特征空间中x₁-x₂平面上的投影显示。

图4是对丹参某个测试样本的识别结果。

图5是图4结果在30维特征空间中x₇-x₈平面上的投影显示。

图6是系统对丹参错误识别的整体显示。

图7是对人参的某个样本的识别结果。

图8是对不在识别范围内(其品种未对系统进行过训练)的某个样本的测试结果。

图9是图8的测试样本与人参训练集样本在30维特征空间中x₁-x₂平面上的投影显示。与其它对系统进行过训练的中药材相比，在特征空间中人参与该样本最接近。

图10是丹参训练集样本红外光谱经归一化和平滑之后的一阶导数谱的″类间类内差异比谱函数″的一段。此段包含了指纹区。其余部分已被切除。

图11是实现训练样本特征向量在多维空间中分布最优化的算法流图。

具体实施方式

由于所获取的红外光谱原始数据文件有可能是以不同的格式存储的(例如：*.asc或*.sp)，并且存储的光谱幅度值可能是“百分透过率(percent transmission)”，也可能是“吸光度(absorbance)”。为了方便后续处理，本方案第一步是将用于训练的每个数据文件中的光谱幅度值读取出来按统一的格式保存到单独的文件目录中(保存前，如果某幅度值是吸光度，则转换为百分透过率)，并将中药材的品种、产地和生长方式信息以统一规定的代码添加到文件名中。

为了突出谱线的变化部分，消除各谱线基线不同带来的影响，我们吸取了别人的经验，对所有红外光谱求取一阶导数，然后从一阶导数谱中提取特征。由于各个原谱在纵轴方向的幅度范围也往往有明显差异，而且往往迭加上了很多并不代表样本特征的很稠密的小幅度波纹。因此在求取一阶导数前我们对原谱进行某种归一化以及平滑化操作。

假定某一品种中药材按产地和生长方式可分为K类，第k类的样本数为N_k，每个红外光谱图的波数数目为n。(如果这种中药材红外光谱图的波数范围不完全一致，例如有的谱图的波数范围为4000cm^-1到400cm^-1，另外的谱图的波数范围为4000cm^-1到650cm^-1，则取其公有范围内的光谱图数据作为后续处理的信号。在上述情况下，n＝4000-650+1＝3351)。我们用d(i，j，k)表示红外光谱图经过预处理后第k类中的第j个谱图的第i个幅度值，其中，i＝1，2，…，n；j＝1，2，…，N_k；k＝1，2，…，K。总的训练样本数为

$N=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{N}_k---\left(1\right)$

为了导出“类间类内差异比”谱函数的定义，我们先给出几个术语：

·第k类均值谱

$M(i,k)=\frac{1}{N_k}\underset{j=1}{\overset{N_k}{\mathrm{\Σ}}}d(i,j,k)---\left(2\right)$

·第k类方差谱：

$V(i,k)=\frac{1}{N_k}\underset{j=1}{\overset{N_k}{\mathrm{\Σ}}}{[d(i,j,k)-M(i,k)]}^2---\left(3\right)$

·类内方差均值谱(MoV，即Mean of Variances)：

$M_V\left(i\right)=\frac{1}{K}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}V(i,k)---\left(4\right)$

·类内均值方差谱(VoM，即Variance of Means)：

$V_M\left(i\right)=\frac{1}{K}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{[M(i,k)-M_M\left(i\right)]}^2---\left(5\right)$

式中，

$M_M\left(i\right)=\frac{1}{K}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}M(i,k).---\left(6\right)$

由式(4)和(5)，我们可以针对该中药材定义一种特殊的谱图——“类间类内差异比”谱函数(或者称为VoM/MoV函数)：

$R\left(i\right)=\frac{V_M\left(i\right)}{M_V\left(i\right)}---\left(7\right)$

在统计的意义上，V_M(i)和M_V(i)分别描述了这种中药材的红外光谱在点i处的类间差异和类内差异的大小。从分类的目的出发，V_M(i)越大越好，M_V(i)越小越好。在VoM/MoV函数值较大的i处取红外光谱幅度值来提取特征，将可以使特征空间中不同类别间的特征向量平均距离大而相同类别间的特征向量平均距离小，有利于分类。我们目前的方案是在求出VoM/MoV函数之后，用一种算法寻找出它在红外光谱指纹区的若干个局部极大点取这些点处的光谱幅度值，将每个光谱图转化为n₁维特征向量。一般取n₁为60至90。

实际上我们能够获得的训练集样本数不多。在此情况下，为了解决模式识别技术中所谓的维数灾难问题，本方案采用K-L变换将n₁维特征向量降维为n₂维特征向量(取n₂=30左右)。在此过程中，我们得到一个n₂×n₁的降维矩阵。

使训练样本特征向量在多维空间中分布最优化，以便能在特征空间中采用简单而高效的算法对未知类别的样本进行分类，是本发明除“类间类内差异比”谱函数之外的另一创新点。

假定我们用X(j，k)表示经过降维后第k类中的第j个特征向量，其中，j＝1，2，…，N_k，k＝1，2，…，K。各类均值向量和协方差矩阵分别为

$M_k=\frac{1}{N_k}\underset{j=1}{\overset{N_k}{\mathrm{\Σ}}}X(j,k)---\left(8\right)$

$C_k=\frac{1}{N_k}\underset{j=1}{\overset{N_k}{\mathrm{\Σ}}}[X(j,k)-M_k]{[X-(j,k)-M_k]}^T---\left(9\right)$

令

$\stackrel{\‾}{C}=\frac{1}{N}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{N}_k{C}_k---\left(10\right)$

是各类协方差矩阵的以(N_k/N)加权的平均的协方差矩阵。它的各特征值往往彼此相差很大，各类训练样本在特征空间中呈现很扁平的分布，使得样本间的欧氏距离难以反映样本间的相似性。我们通过一种迭代算法对特征向量进行变换。在有限的几步迭代后就可以使新

的各特征值几乎完全相等。这时，从平均意义上就可以认为各类样本在特征空间中的分布是优化的。样本数越多的类，其样本越呈“球形分布”。这样就可以在对未知类别的样本进行分类操作时用简单的分类算法。迭代过程将产生一个n₂×n₂的分布优化变换矩阵。

对送检的中药材红外光谱，进行与每一种中药材训练阶段一致的处理，使之映像到每个训练样本集的特征空间中。然后，用合适的分类算法判别该中药材的种类、产地和生长方式。本方案目前采用的分类算法是改进了的近邻法。

本发明在Visual C++6.0编程的软件系统中具体实施。该系统具有科学计算可视化功能，有友好用户界面。

我们先后共取得丹参、人参、天麻和淫羊霍4种中药材的中红外光谱图，该系统目前可全部鉴别它们。其中丹参的样本数最多，红外光谱图也较规范，因此下面主要以丹参为例来说明具体实施方式。我们获得了940个丹参样本的中红外光谱。这些样本来自27个产地，生长条件分为采挖野生(有些又进一步细分为阴坡和阳坡)，野生茎叶，采挖栽培(有些又进一步细分为一年和二年)，栽培茎叶，购买野生，购买栽培等。我们的任务是：既要对丹参的产地进行区分，又要对同一产地的丹参生长条件进行区分。类别总数有49类。训练样本总数与类别总数相比是很少的，而且各类样本数目严重不均。因此，必须在信号的预处理和特征提取上狠下功夫，尽可能使同一类别的特征向量相距较近而不同类别的特征向量相距较远，然后针对小样本情况考虑分类方法。

我们将样本划分为870个训练样本和70个测试样本。图1及图2显示了丹参的训练集样本和测试集样本的某些信息。图1在地图上显示了不同训练集样本的产地，该地图的下方有丹参样本的其它数据，如生长条件等。图2显示了用于测试的丹参样本的产地、生长条件及红外光谱的文件名称。

对于870个训练样本，我们进行下列操作：

预处理：对各样本光谱的幅度进行归一化，用长度为13的三角形窗进行平滑，再求一阶导数谱。

求一阶导数谱的“类间类内差异比谱函数”，将此函数的幅度归一化后，切去对分类无意义的部分，保留包含了红外光谱图指纹区的部分，如图10所示。

求这一段中的局部极大点(此例中有64个)，将其保存。在这些局部极大点处提取一阶导数谱的幅度，得到870个64维的特征向量。

求出这些特征向量总体的协方差矩阵及其特征向量和特征值，按K-L变换降维的原则构造出一个30x64的变换矩阵，将各特征向量降到30维。(之所以是30，一是因为变换后的新的特征向量总体的协方差矩阵的主对角线元素之和下降不多，仅比变换前下降6％左右；二是可以保证对70个测试样本有较高的正确识别率)。

对于降维后的训练样本特征向量，用图11所示的迭代算法实现其在多维空间中的分布最优化。每次迭代循环中包含下列步骤：(1)由训练集所有样本的特征向量计算出每一类的均值向量及协方差矩阵；(2)将每一类的协方差矩阵以样本数的比率加权后求和，得到平均协方差矩阵；(3)计算出该平均协方差矩阵的各个本征值及对应的本征向量；(4)将最大本征值和最小本征值的差与一个既定的微小阈值(此例中是0.001)比较；(5)如果最大本征值和最小本征值的差小于该既定阈值，则转向步骤(6)；否则，将各本征值的倒数组成的对角矩阵左乘本征向量矩阵的转置矩阵，构成一个变换矩阵；(6)利用此变换矩阵对各特征向量进行线性变换，所得结果作为新的特征向量，转向步骤(1)；(7)结束迭代，将各中间过程的变换矩阵相乘得到总的分布优化变换矩阵。

图3显示了经过空间分布优化的30维的丹参训练集样本特征向量在一个平面(x₁-x₂平面)上的投影。图中用连线标出了某个类别的五个样本。可以看出，同类别的样本在特征空间里的距离是较近的。

训练集样本特征向量经过空间分布优化后，即可利用一定的分类算法对未知类别的样本进行识别了。如果系统对样本的正确识别率高，则将处理过程和对分类有用的数据保存备用。否则，要重新考虑从预处理到分类算法的各个环节。

本发明采用的分类算法是改进的近邻法。图4显示了对丹参某个测试样本的识别结果。图的右上方可以看到该样本的红外光谱。图的下方显示了最接近该测试样本的二十个训练集样本。图5显示了这个测试样本的特征向量及其二十个最接近的训练集样本的特征向量在一个平面上(x₇-x₈平面)的投影。

图6显示了系统对中药丹参的错误识别率。对于未参与训练的70个不同产地、不同生长条件的丹参测试样本，只有4个被错误识别。这说明本发明提供的特征提取和分类方法是有效的。

图7显示了系统对人参的某个样本的识别结果。结果是正确的。

图8和图9显示了系统对一个在丹参、人参、天麻和淫羊霍4种中药材范围之外的某个样本的识别结果。由图8中的数值可以看出，该样本与所有品种的训练样本在特征空间中的距离都超出了一定范围，因此本系统将该样本视作不可识别的。图9是该样本与人参训练集样本特征向量的2维显示。与其它三种对系统进行过训练的中药材相比，在特征空间中人参与该样本最接近，但距离仍明显偏大。

一个以红外光谱特征识别样本类别的系统可以通过服务器与互联网连接，也可直接与傅立叶变换红外光谱仪连接。远程用户可将需要识别的红外光谱数据通过互联网传送至系统。本地用户可用光谱仪直接将待识别的实物样品的红外光谱数据输入系统。系统中的专用软件将根据训练阶段建立的数据库和探索出的算法对输入的红外光谱进行一系列处理，最后完成识别任务，并将结果以文字和图形形式显示给用户。

Claims (4)

1.一种提取一个对象品种的红外光谱特征的方法，包括以下步骤：

(1)根据最少一个分类条件，将该品种的对象分成多个类别，并获取属于该品种的多个已知类别的样本；

(2)将该多个样本以随机方式分成一个训练样本集及一个测试样本集，该训练样本集包括多个训练样本，该测试样本集包括多个测试样本；

(3)测定该多个样本的红外光谱；

(4)对该多个训练样本的红外光谱进行一系列预处理，获得该多个训练样本的一阶导数谱；

(5)从该多个训练样本的一阶导数谱，计算该品种的类间类内差异比谱函数；类间类内差异比谱函数的定义为类内均值方差谱除以类内方差均值谱；类内均值方差谱及类内方差均值谱从该品种各类训练样本的一阶导数谱的幅度值计算出；

(6)将该类间类内差异比谱函数归一化；

(7)求该类间类内差异比谱函数的多个局部极大点；

(8)在该多个局部极大点提取该多个训练样本的一阶导数谱的幅度值，作为该多个训练样本的原始特征向量；

(9)对该多个训练样本的原始特征向量进行降维处理，得到一个降维矩阵；

(10)对该多个降维后的特征向量进行多维特征空间中的分布优化，从而得到一个分布优化变换矩阵，包括以下步骤：

A.由该多个训练样本的特征向量计算出该多个类别的均值向量及协方差矩阵；

B.将该每一个类别的协方差矩阵以样本数的比率加权后求和，得到平均协方差矩阵；

C.计算出该平均协方差矩阵的多个本征值及对应的本征向量；

D.将最大本征值和最小本征值的差与一个既定的微小阈值进行比较；

E.如果最大本征值和最小本征值的差小于该既定的微小阈值，则转向步骤H；否则，将各本征向量合并成为一个本征向量矩阵，并用相应的各本征值的倒数组成的对角矩阵左乘该本征向量矩阵的转置矩阵，其结果作为中间过程的变换矩阵；

F.利用该中间过程的变换矩阵对该多个特征向量进行线性变换，所得结果作为新的特征向量；

G.重复步骤A至步骤F，直至最大本征值和最小本征值的差比该既定的微小阈值小；以及

H.将该各中间过程的变换矩阵依次相乘，以最后一个在最左及第一个在最右，得到总的分布优化变换矩阵；

在步骤A至步骤G的迭代中，每一次迭代均会产生一个中间过程变换矩阵，而且均会改善训练样本特征向量在该多维特征空间中的分布；这种迭代一直进行到该平均协方差矩阵的最大本征值和最小本征值的差小于该既定阈值为止；

(11)确定距离阈值，设计分类算法，并根据该距离阈值及该分类算法对该多个测试样本进行试分类；以及

(12)调节该距离阈值及该分类算法，重复步骤(4)到步骤(11)，直至该试分类所测得的正确识别率高过一个既定的比率，然后将该类间类内差异比谱函数的多个局部极大点、降维矩阵、分布优化变换矩阵、距离阈值及分类算法储存。

2.如权利要求1所述的方法，其特征在于所述对象为中药材，并且所述的分类条件至少包括产地及生长条件。

3.如权利要求1所述的方法，其特征在于所述降维处理为K-L变换法。

4.如权利要求1所述的方法，其特征在于所述预处理包括归一化、平滑运算及求导运算。

Images