CN109283169A - 一种鲁棒的拉曼光谱峰识别方法 - Google Patents

Abstract

本发明提供了一种鲁棒的拉曼光谱峰识别方法，包括计算拉曼光谱数据的一阶导数；记录拉曼光谱数据的一阶导数的正负号发生改变的位置；根据拉曼光谱数据的一阶导数的正负号发生改变的位置，判断拉曼光谱数据的峰谷类型值的位置；所述峰谷类型值包括左侧谷值、峰值、右侧谷值及中间谷值；根据拉曼光谱数据的左侧谷值、峰值、右侧谷值及中间谷值的位置，利用高斯函数作为曲线拟合数学模型，采用最小中值二乘法对拉曼光谱数据进行拟合，得到特征峰的位置及强度。本发明的有益效果在于：该方法能够避免现有技术对拉曼光谱数据预处理产生的畸变，提高了拉曼光谱分析的精度，同时能够对拉曼光谱进行定量分析。

Description

一种鲁棒的拉曼光谱峰识别方法

技术领域

本发明涉及光谱峰识别领域，尤其是指一种鲁棒的拉曼光谱峰识别方法。

背景技术

拉曼光谱分析是一种常用的物质分析与检测手段，拉曼光谱分析可分为定性分析与定量分析；定性分析的基本原理为：不同位置的拉曼光谱发光峰位置不同，通过观察拉曼光谱在特定位置是否有发光峰能够判断是否存在某种物质。定量分析的基本原理为：拉曼光谱峰的强度(特征峰的面积)与相应物质的含量之间存在正比线性关系，所以，分析拉曼光谱峰的强度能够得到该种物质的含量。所以，对于拉曼光谱分析而言，关键在于确定光谱峰的位置以及光谱峰强度。然而，拉曼光谱常常会受到噪声干扰，主要原因是拉曼光谱信号强度较弱，仅约为原始激发信号的10^-8。由此，造成了拉曼光谱分析的困难。

目前，针对该问题的一种解决思路为首先对拉曼光谱进行降噪预处理，而后对光谱峰进行识别。例如：对拉曼光谱进行平滑处理。该方法能够很好的去除随机噪声，提高信噪比，但是，在平滑过程中降低了特征峰的峰值，使得难以识别一些峰值小的特征峰。利用小波变换进行降噪处理，该方法能够较好的保留特征峰，寻峰准确率较高，但是，计算量大，且小波基的选取对分析结果有显著影响。

针对该问题的另一种解决方案为：首先估计噪声的标准差，然后利用光谱的一阶导数获得可能的峰值与谷值，最后对比峰值与相邻的谷值。如果峰值与谷值之差大于r倍(通常取大于6)的噪声标准差时，则认为是拉曼光谱峰。该方法能够在不对拉曼光谱做降噪预处理的情况下，自动实现特征峰识别。但是，受噪声的影响，该方法得到的特征峰位置的误差较大，并且未能获得特征峰的强度，难以用于定量分析。

发明内容

本发明所要解决的技术问题是：克服目前拉曼光谱分析时难以获得准确的特征峰位置及获得的特征峰强度不足的问题。

为了解决上述技术问题，本发明采用的技术方案为：一种鲁棒的拉曼光谱峰识别方法，包括步骤，

S10、计算拉曼光谱数据的一阶导数；

S20、记录拉曼光谱数据的一阶导数的正负号发生改变的位置；

S30、根据拉曼光谱数据的一阶导数的正负号发生改变的位置，判断拉曼光谱数据的峰谷类型值的位置；所述峰谷类型值包括左侧谷值、峰值、右侧谷值及中间谷值，并将峰谷类型值作为元素组成一个向量z；

S40、根据拉曼光谱数据的左侧谷值、峰值、右侧谷值及中间谷值的位置，利用高斯函数作为曲线拟合数学模型，采用最小中值二乘法对拉曼光谱数据进行拟合，得到特征峰的位置及强度。

进一步的，所述步骤S30具体包括，

当一阶导数由零变为正数时，判断该位置为左侧谷值的位置；

当一阶导数由正数变为负数时，判断该位置为峰值的位置；

当一阶导数由负数变为零且在零处持续一段时间，判断该位置为右侧谷值的位置；

当一阶导数由负数变为正数时，判断该位置为中间谷值的位置。

进一步的，所述步骤S40具体包括，

对向量z中的元素从左往右进行判断；

当判断到向量z中的第一个元素为峰值时，将原光谱数据的最左侧位置作为左侧谷值；

当判断到向量z中的最后一个元素为峰值时，将原光谱数据的最右侧位置作为右侧谷值；

从左往右找到第一个峰值，获取峰值的数据片段，进行光谱拟合。

进一步的，当从左往右找到第一个峰值后，还包括对所述峰值相邻元素的元素类型进行判断的步骤，具体包括，

当第一个峰值的相邻两个元素的类型分别为左侧谷值及右侧谷值时，将左侧谷值与右侧谷值之间的数据作为待拟合的单特征峰数据片段；

采用最小中值二乘法对单特征峰数据片段进行拟合。

进一步的，当第一个峰值的左右两侧元素中存在中间峰时，该特征峰属于多峰叠加的一个子峰；在向量z中，由第一个峰值的位置向左直至寻找到左侧谷值，作为该多峰叠加的最左侧位置；由第一个峰值的位置向右直至寻找到右侧谷值，作为该多峰叠加的最右侧位置；

将左侧谷值与右侧谷值之间的数据作为待拟合的多峰叠加数据片段；

采用最小中值二乘法对多个特征峰的数据片段进行拟合。

进一步的，所述单个特征峰的数学模型为：

其中，I为拉曼光谱强度，I_base为基底强度，R_c为特征峰的中心，R为拉曼频移，w为半高宽，S为高斯峰面积。

进一步的，所述单个特征峰的最小中值二乘法的目标函数为：

其中，I_0i为原始拉曼光谱数据片段中的第i个元素，上式表示的意义为寻求

一组参数I_base、R_c、w、S，使得由该组参数拟合得到的高斯光谱与原始拉曼光谱数据片段之差的绝对值的中值最小。

进一步的，所述多个特征峰的数学模型为：

其中，j为多个峰的索引，n为特征峰的个数，I_0j、R_cj、w_j、S_j分别为第j个特征峰的基底强度、中心、半高宽及面积；特征峰的个数为该多峰叠加数据片段中中间谷值的个数加一。

进一步的，所述多个特征峰的最小中值二乘法的目标函数为：

其中，上式表示的意义为寻求一组参数I_basej、R_cj、w_j、S_j，j从1到n，使得由该组参数拟合得到的高斯光谱与原始拉曼光谱数据片段之差的绝对值的中值最小。

本发明的技术方案通过对拉曼光谱数据的一阶导数来判断拉曼光谱数据的峰谷类型及峰谷值的位置，根据峰谷值的位置，利用高斯函数作为曲线拟合数学模型，采用最小中值二乘法对拉曼光谱数据进行拟合，得到特征峰的位置及强度；该方法能够在不需要预先对拉曼光谱进行降噪预处理的情况下，获得特征峰位置及强度，避免了现有技术对拉曼光谱数据预处理产生的畸变，提高了拉曼光谱分析的精度，同时拉曼光谱能够用于定量分析。

附图说明

下面结合附图详述本发明的具体结构。

图1为本发明一具体实施例的数据处理总体框图；

图2为本发明一具体实施例的普峰识别流程图；

图3为本发明一具体实施例的一阶前向差分值与峰谷值对应关系图；

图4为本发明一具体实施例的退火算法计算基于最小中值二乘法的高斯线型拟合流程图；

图5为本发明一具体实施例的鲁棒的拉曼光谱峰识别方法流程图。

具体实施方式

为详细说明本发明的技术内容、构造特征、所实现目的及效果，以下结合实施方式并配合附图详予说明。

实施例一

参考图5，一种鲁棒的拉曼光谱峰识别方法，包括步骤，

S10、计算拉曼光谱数据的一阶导数；

S20、记录拉曼光谱数据的一阶导数的正负号发生改变的位置；

S30、根据拉曼光谱数据的一阶导数的正负号发生改变的位置，判断拉曼光谱数据的峰谷类型值的位置；所述峰谷类型值包括左侧谷值、峰值、右侧谷值及中间谷值，并将峰谷类型值作为元素组成一个向量z；

S40、根据拉曼光谱数据的左侧谷值、峰值、右侧谷值及中间谷值的位置，利用高斯函数作为曲线拟合数学模型，采用最小中值二乘法对拉曼光谱数据进行拟合，得到特征峰的位置及强度。

本技术方案中，能够在不需要预先对拉曼光谱进行降噪预处理的情况下，获得特征峰位置及强度，避免了现有技术对拉曼光谱数据预处理产生的畸变，提高了拉曼光谱分析的精度。

进一步的，所述步骤S30具体包括，

当一阶导数由零变为正数时，判断该位置为左侧谷值的位置；

当一阶导数由正数变为负数时，判断该位置为峰值的位置；

当一阶导数由负数变为零且在零处持续一段时间，判断该位置为右侧谷值的位置；

当一阶导数由负数变为正数时，判断该位置为中间谷值的位置。

进一步的，所述步骤S40具体包括，

对向量z中的元素从左往右进行判断；

当判断到向量z中的第一个元素为峰值时，将原光谱数据的最左侧位置作为左侧谷值；

当判断到向量z中的最后一个元素为峰值时，将原光谱数据的最右侧位置作为右侧谷值；

从左往右找到第一个峰值，获取峰值的数据片段，进行光谱拟合。

进一步的，当从左往右找到第一个峰值后，还包括对所述峰值相邻元素的元素类型进行判断的步骤，具体包括，

当第一个峰值的相邻两个元素的类型分别为左侧谷值及右侧谷值时，将左侧谷值与右侧谷值之间的数据作为待拟合的单特征峰数据片段；

采用最小中值二乘法对单特征峰数据片段进行拟合。

进一步的，当第一个峰值的左右两侧元素中存在中间峰时，该特征峰属于多峰叠加的一个子峰；在向量z中，由第一个峰值的位置向左直至寻找到左侧谷值，作为该多峰叠加的最左侧位置；由第一个峰值的位置向右直至寻找到右侧谷值，作为该多峰叠加的最右侧位置；

将左侧谷值与右侧谷值之间的数据作为待拟合的多峰叠加数据片段；

采用最小中值二乘法对多个特征峰的数据片段进行拟合。

进一步的，所述单个特征峰的数学模型为：

其中，I为拉曼光谱强度，I_base为基底强度，R_c为特征峰的中心，R为拉曼频移，w为半高宽，S为高斯峰面积。

进一步的，所述单个特征峰的最小中值二乘法的目标函数为：

其中，I_0i为原始拉曼光谱数据片段中的第i个元素，上式表示的意义为寻求一组参数I_base、R_c、w、S，使得由该组参数拟合得到的高斯光谱与原始拉曼光谱数据片段之差的绝对值的中值最小。

进一步的，所述多个特征峰的数学模型为：

其中，j为多个峰的索引，n为特征峰的个数，I_0j、R_cj、w_j、S_j分别为第j个特征峰的基底强度、中心、半高宽及面积；特征峰的个数为该多峰叠加数据片段中中间谷值的个数加一。

进一步的，所述多个特征峰的最小中值二乘法的目标函数为：

其中，上式表示的意义为寻求一组参数I_basej、R_cj、w_j、S_j，j从1到n，使得由该组参数拟合得到的高斯光谱与原始拉曼光谱数据片段之差的绝对值的中值最小。

本实施例中，能够在避免对原始数据预处理的情况下对拉曼光谱进行分析，因为降噪预处理会使得特征峰产生畸变，所以本发明能够提高拉曼光谱分析精度；同时能够用于对拉曼光谱进行定量分析。

实施例二

如图1所示，本发明提供的一种鲁棒的鲁棒的拉曼光谱峰识别方法主要包含两步：首先初步找到拉曼光谱的峰值和谷值，然后利用高斯数学模型作为特征峰模型，采用最小中值二乘算法对拉曼特征峰进行拟合，计算特征峰的位置以及强度。

如图2所示，详细的拉曼光谱特征峰识别具体包含以下步骤：

步骤1，确定可能的拉曼光谱特征峰。

步骤1.1，计算拉曼光谱数据的一阶导数(即一阶前向差分)，得到可能得峰值和谷值，将峰值和谷值组成一个新的向量z。如图3所示，对于两个高斯峰组合的光谱(图3(a))，计算得到的一阶前向差分如图3(b)。一阶前向差分中反应了所有峰值与谷值的位置；定义图中用圆所示位置为左侧谷值处，该位置处一阶前向差分由零变为正值；定义图中三角形所示位置为峰值处，该位置处一阶前向差分由正数变为负数；定义菱形位置处为中间谷值处，该位置处一阶前向差分由负数变为正数；定义图中五边形位置处为右侧谷值处，该位置处一阶前向差分由负数变为零。根据这些性质，判断峰谷类型。

步骤1.2，记录所有一阶导数的正负号发生改变的位置，并记录每一个z中元素对应位置的类型，分为左侧谷值、峰值、右侧谷值及中间谷值。

步骤2，利用高斯函数作为曲线拟合数学模型，结合最小中值二乘法拟合，计算特征峰的位置及强度。

步骤2.1，当z中，第一个元素为峰值时，将原光谱数据的最左侧位置作为左侧谷值。

当z中，最后一个元素为峰值时，将原光谱数据的最右侧位置作为右侧谷值。

步骤2.2，在z中，从左往右找到第一个类型为峰值的元素，获取当前特征峰的数据片段，进行光谱拟合。

定义，当z中相邻三个元素对应位置的类型分别为左侧谷值、峰值和右侧谷值时，该段特征峰为单峰；当z中相邻的左侧谷值和右侧谷值中间存在中间谷值时，该段特征峰为多峰叠加，峰的数量为中间谷值数量加一。

当该峰值的相邻两个元素的类型分别为左侧谷值及右侧谷值时，该特征峰为单峰，并将左侧谷值与右侧谷值之间的数据作为待拟合的单峰数据片段。

当该峰值的左右两侧元素中存在中间谷值时，该特征峰属于多峰叠加的一个子峰。在z中，由该峰位置处向左直至寻找左侧谷值，即为该多峰叠加的最左侧位置；由该峰位置处向右直至寻找到右侧谷值，即为该多峰叠加的最右侧位置，将左侧谷值与右侧谷值之间的数据作为待拟合的多峰叠加数据片段。

步骤2.3，采用最小中值二乘法对单特征峰数据片段进行拟合。

高斯函数作为单个特征峰的数学模型为：

其中，I为拉曼光谱强度(因变量)，I_base为基底强度，R_c为特征峰的中心，R为拉曼频移(自变量)，w为半高宽，S为高斯峰面积。

最小中值二乘法的目标函数表示为：

其中，I_0i为原始拉曼光谱数据片段中的第i个元素。上式表示的意义为寻求一组参数I_base、R_c、w、S，使得由该组参数拟合得到的高斯光谱与原始拉曼光谱数据片段之差的绝对值的中值最小。

多个特征峰的数学模型为：

其中，j为多个峰的索引，n为特征峰的个数，I_0j、R_cj、w_j、S_j分别为第j个特征峰的基底强度、中心、半高宽及面积，特征峰的个数为该多峰叠加数据片段中中间谷值的个数加一。

最小中值二乘法的目标函数为：

上式表示的意义为寻求一组参数I_basej、R_cj、w_j、S_j，(j从1到n)，使得由该组参数拟合得到的高斯光谱与原始拉曼光谱数据片段之差绝对值的中值最小。

单峰可作为多峰的一种特殊类型对待，在多峰的数学模型中n取1即为单峰模型。

步骤2.4，求解该优化问题，确定特征峰的参数，给出特征峰的位置以及强度。求解方法可采用模拟退火算法、遗传算法、蚁群算法。

如图4所示，本实施例中，以模拟退火算法进行多峰高斯模型的最小二乘拟合为例进行说明。

步骤2.4.1，设定模拟退火算法中所需参数，设定初始温度T，降温系数a，总迭代次数N，迭代次数n＝1以及收敛阈值thres。

步骤2.4.2，随机生成待拟合的高斯模型中的参数初始值x。这里x为向量，向量中包含所有高斯模型中待拟合的参数。

计算高斯光谱I(x)，计算高斯光谱与原始拉曼光谱数据之间的差值，并计算差值的绝对值的中值I_m0。

步骤2.4.3，随机扰动产生新的解x_new，产生方法为在x的基础上增加一个随机数；对于x为矢量而言，则是矢量中的每个变量增加一个随机数，计算高斯光谱I(x_new)。计算高斯光谱与原始拉曼光谱数据之间的差值。并计算差值的绝对值的中值I_{m_new}。

步骤2.4.4，比较I_m0和I_{m_new}，计算：

dI＝I_{m_new}-I_m0

如果dI<0，则接受新解，令x＝x_new，以及I_m0＝I_{m_new}，进一步判断是否达到迭代次数上限，并将迭代步数加1；

否则，以一定的概率接受该新解，该概率通过公式，

p＝exp(-dI/T)

计算得到，并将计算结果与一个0～1之间的随机数比较。如果p大于该随机数，则接受该解；否则，判断是否达到迭代次数上限，并将迭代步数加1。

如果达到迭代次数上限，则判断是否满足收敛条件，收敛条件为I_m0<thres。

如果未达到迭代次数上限，则跳转到步骤2.4.3。

如果满足，则输出此时的解x。

否则，降低温度，

T＝aT

并将迭代步数重置为1，跳转到步骤2.4.3。

本实施例的技术方案通过对拉曼光谱数据的一阶导数来判断拉曼光谱数据的峰谷类型及峰谷值的位置，根据峰谷值的位置，利用高斯函数作为曲线拟合数学模型，采用最小中值二乘法对拉曼光谱数据进行拟合，得到特征峰的位置及强度；该方法能够在不需要预先对拉曼光谱进行降噪预处理的情况下，获得特征峰位置及强度，避免了现有技术对拉曼光谱数据预处理产生的畸变，提高了拉曼光谱分析的精度，同时能够对拉曼光谱进行定量分析。

此处第一、第二……只代表其名称的区分，不代表它们的重要程度和位置有什么不同。

此处，上、下、左、右、前、后只代表其相对位置而不表示其绝对位置。以上所述仅为本发明的实施例，并非因此限制本发明的专利范围，凡是利用本发明说明书及附图内容所作的等效结构或等效流程变换，或直接或间接运用在其他相关的技术领域，均同理包括在本发明的专利保护范围内。

Claims (9)

1.一种鲁棒的拉曼光谱峰识别方法，其特征在于：包括步骤，

S10、计算拉曼光谱数据的一阶导数；

S20、记录拉曼光谱数据的一阶导数的正负号发生改变的位置；

S30、根据拉曼光谱数据的一阶导数的正负号发生改变的位置，判断拉曼光谱数据的峰谷类型值的位置；所述峰谷类型值包括左侧谷值、峰值、右侧谷值及中间谷值，并将峰谷类型值作为元素组成一个向量z；

S40、根据拉曼光谱数据的左侧谷值、峰值、右侧谷值及中间谷值的位置，利用高斯函数作为曲线拟合数学模型，采用最小中值二乘法对拉曼光谱数据进行拟合，得到特征峰的位置及强度。

2.如权利要求1所述的鲁棒的拉曼光谱峰识别方法，其特征在于：所述步骤S30具体包括，

当一阶导数由零变为正数时，判断该位置为左侧谷值的位置；

当一阶导数由正数变为负数时，判断该位置为峰值的位置；

当一阶导数由负数变为零且在零处持续一段时间，判断该位置为右侧谷值的位置；

当一阶导数由负数变为正数时，判断该位置为中间谷值的位置。

3.如权利要求1所述的鲁棒的拉曼光谱峰识别方法，其特征在于：所述步骤S40具体包括，

对向量z中的元素从左往右进行判断；

当判断到向量z中的第一个元素为峰值时，将原光谱数据的最左侧位置作为左侧谷值；

当判断到向量z中的最后一个元素为峰值时，将原光谱数据的最右侧位置作为右侧谷值；

从左往右找到第一个峰值，获取峰值的数据片段，进行光谱拟合。

4.如权利要求3所述的鲁棒的拉曼光谱峰识别方法，其特征在于：

当从左往右找到第一个峰值后，还包括对所述峰值相邻元素的元素类型进行判断的步骤，具体包括，

当第一个峰值的相邻两个元素的类型分别为左侧谷值及右侧谷值时，将左侧谷值与右侧谷值之间的数据作为待拟合的单特征峰数据片段；

采用最小中值二乘法对单特征峰数据片段进行拟合。

5.如权利要求4所述的鲁棒的拉曼光谱峰识别方法，其特征在于：

当第一个峰值的左右两侧元素中存在中间峰时，该特征峰属于多峰叠加的一个子峰；在向量z中，由第一个峰值的位置向左直至寻找到左侧谷值，作为该多峰叠加的最左侧位置；由第一个峰值的位置向右直至寻找到右侧谷值，作为该多峰叠加的最右侧位置；

将左侧谷值与右侧谷值之间的数据作为待拟合的多峰叠加数据片段；

采用最小中值二乘法对多个特征峰的数据片段进行拟合。

6.如权利要求5所述的鲁棒的拉曼光谱峰识别方法，其特征在于：所述单个特征峰的数学模型为：

其中，I为拉曼光谱强度，I_base为基底强度，R_c为特征峰的中心，R为拉曼频移，w为半高宽，S为高斯峰面积。

7.如权利要求6所述的鲁棒的拉曼光谱峰识别方法，其特征在于：所述单个特征峰的最小中值二乘法的目标函数为：

其中，I_0i为原始拉曼光谱数据片段中的第i个元素，上式表示的意义为寻求一组参数I_base、R_c、w、S，使得由该组参数拟合得到的高斯光谱与原始拉曼光谱数据片段之差的绝对值的中值最小。

8.如权利要求5所述的鲁棒的拉曼光谱峰识别方法，其特征在于：所述多个特征峰的数学模型为：

其中，j为多个峰的索引，n为特征峰的个数，I_0j、R_cj、w_j、S_j分别为第j个特征峰的基底强度、中心、半高宽及面积；特征峰的个数为该多峰叠加数据片段中中间谷值的个数加一。

9.如权利要求8所述的鲁棒的拉曼光谱峰识别方法，其特征在于：所述多个特征峰的最小中值二乘法的目标函数为：

其中，上式表示的意义为寻求一组参数I_basej、R_cj、w_j、S_j，j从1到n，使得由该组参数拟合得到的高斯光谱与原始拉曼光谱数据片段之差的绝对值的中值最小。