CN107727634B - 一种激光拉曼光谱解谱处理方法 - Google Patents

Abstract

本发明涉及一种激光拉曼光谱解谱处理方法，包括以下步骤：1)对拉曼光谱原始信号进行基于自适应迭代重加权惩罚最小二乘法的基线校正处理和基于惩罚最小二乘法的平滑处理；2)对经所述步骤1)处理后的信号进行基于连续小波变换的峰检测。与现有技术相比，本发明具有速度快、峰检测精度高等优点。

Description

一种激光拉曼光谱解谱处理方法

技术领域

本发明涉及信号处理方法，尤其是涉及一种激光拉曼光谱解谱处理方法。

背景技术

拉曼光谱可以被看作是一项“指纹”技术，因为它可以提供非常丰富的结构信息。因此拉曼光谱可以被用作物质的定性识别。并且拉曼光谱具有制样简单，不破坏样品，在几乎所有的环境下都可以采集。

由于拉曼光谱具有上述的优点，故在化学品成分分析中被广泛应用。但是拉曼光谱激光源是可见光，所以有易产生噪声，荧光干扰严重的缺点。这些缺点会影响对样品的定性定量分析，然而现在的硬件技术无法避免这些缺点，所以这时就需要使用数学算法对拉曼光谱图进行后期的处理以达到过滤噪声和荧光的目的。

中国专利CN 103217409B公开了一种拉曼光谱的预处理方法，其使用基于小波变换的自适应阈值去噪声，采用非对称最小二乘的基线校正算法去除荧光背景，存在以下不足：首先，平滑参数需要优化以便得到最优结果；其次，非对称参数对于所有的基线数据点都是一成不变的，这样基线可能会出现负值部分。

发明内容

本发明的目的就是为了克服上述现有技术存在的缺陷而提供一种激光拉曼光谱解谱处理方法。

本发明的目的可以通过以下技术方案来实现：

一种激光拉曼光谱解谱处理方法，包括以下步骤：

1)对拉曼光谱原始信号进行基于自适应迭代重加权惩罚最小二乘法的基线校正处理和基于惩罚最小二乘法的平滑处理；

2)对经所述步骤1)处理后的信号进行基于连续小波变换的峰检测。

进一步地，所述基于自适应迭代重加权惩罚最小二乘法的基线校正处理中，所述自适应迭代重加权惩罚最小二乘法的表达式为：

式中，Q为原始基线与拟合基线保真度与粗糙度间的平衡，x＝[x₁,x₂,...,x_m]为原始信号向量，z＝[z₁,z₂,...,z_m]为校正向量，m为向量长度，λ为粗糙度系数，i、j为计数变量，t为迭代次数，w为权重向量，由自适应迭代方法获得，表达式为：

式中，d^t为在第t次迭代时原始信号向量与第t-1次迭代获得的校正向量的负值的系数，权重向量的初始值为1。

进一步地，所述基于自适应迭代重加权惩罚最小二乘法的迭代结束条件包括达到最大迭代次数或满足表达式|d^t|<0.001×|x|。

进一步地，所述基于惩罚最小二乘法的平滑处理中，引入一帽子矩阵或平滑矩阵获得最优平滑。

进一步地，所述基于连续小波变换的峰检测具体为：

201)通过连接在各放缩级别的连续小波变换系数的极大值识别脊线；

202)识别所述脊线对应的峰；

203)对所述峰的参数进行优化。

进一步地，所述步骤1)中，先进行所述基线校正处理，后进行所述平滑处理，或者，

先进行所述平滑处理，后进行所述基线校正处理。

进一步地，所述步骤2)还包括：估计拉曼光谱图中峰的宽度。

与现有技术相比，本发明具有以下有益效果：

(1)基于自适应迭代重加权惩罚最小二乘法的基线校正不需要任何的用户介入和预处理信息，全部自动化。

(2)基于惩罚最小二乘法的光谱平滑具有快速，可以连续控制平滑度并且可以进行交叉验证得到最客观的平滑值。

(3)改善了基于非对称最小二乘法的传统基线校正方法的缺陷，传统方法对基线区域的点的校正参数是不变的，本发明方法可以根据前一个拟合点与其原始点的差值来更改校正参数以便更好地拟合。

(4)基于连续小波变换的峰检测算法可以自动地并且同时考虑峰形和峰高对峰进行检测，最大的降低了峰检测假阳性的概率。

附图说明

图1为本发明的流程示意图；

图2为自适应迭代重加权惩罚最小二乘法基线校正的运行结构图；

图3为自适应迭代重加权惩罚最小二乘法平滑的运行结构图；

图4为墨西哥帽函数；

图5为原始拉曼光谱图；

图6为单独使用自适应迭代重加权惩罚最小二乘法基线校正的拉曼光谱图；

图7为单独使用惩罚最小二乘法平滑的拉曼光谱图；

图8为同时使用自适应迭代重加权惩罚最小二乘法基线校正和惩罚最小二乘法平滑的拉曼光谱图。

具体实施方式

下面结合附图和具体实施例对本发明进行详细说明。本实施例以本发明技术方案为前提进行实施，给出了详细的实施方式和具体的操作过程，但本发明的保护范围不限于下述的实施例。

如图1所示，本实施例提供一种激光拉曼光谱解谱处理方法，具体步骤如下：

1)对拉曼光谱原始信号进行基于自适应迭代重加权惩罚最小二乘法的基线校正处理和基于惩罚最小二乘法的平滑处理，基线校正处理和平滑处理的先后顺序可任意调整；

2)对经所述步骤1)处理后的信号进行基于连续小波变换的峰检测。

如图2所示，步骤1)中所述对拉曼光谱原始信号进行基于自适应迭代重加权惩罚最小二乘法的算法的具体步骤如下：

2-1)自适应迭代重加权惩罚最小二乘法的表达式为：

式中，Q为原始基线与拟合基线保真度与粗糙度间的平衡，x＝[x₁,x₂,...,x_m]为原始信号向量，z＝[z₁,z₂,...,z_m]为校正向量，m为向量长度，λ为粗糙度系数，i、j为计数变量，t为迭代次数，w为权重向量，由自适应迭代方法获得。

2-2)在迭代开始，给定w一个初始值即：w⁰＝1。迭代开始之后，在每一个迭代步骤t，w均可由下列表达式得到：

向量d^t是在t迭代步骤，原始信号向量与第t-1次迭代获得的校正向量的负值的系数。如果第i个数据点比在(t-1)步迭代时得到的z^t-1大时，这个数据点可以被看作是峰上的一点，所以此时的权重可以设置为零以便使其不进入下一步迭代。在本发明中，这种方法可以在权重向量w中自动地逐步排除峰上的点并保留基线上的点。

2-3)迭代会在达到最大迭代次数或满足下面的条件时结束：

|d^t|<0.001×|x|。

如图3所示，步骤2)中所述对进行完1)的拉曼光谱信号进行基于惩罚最小二乘法的平滑的算法根据原始曲线分量，依次经稀疏矩阵、微分矩阵和平滑矩阵获得粗平滑向量，再经优化获得最终平滑向量。具体步骤如下：

3-1)将2-1)中的加权系数去除，即得到峰平滑的数学表达式：

(I+αD'D)z＝y

式中I为单位矩阵；D为微分矩阵；z为平滑后光谱的向量；Δz＝Dz；y为原始曲线向量，α为平滑度。

3-2)依次将y中的一个数据取出，然后对剩下的数据进行平滑，此时对于取出的y_i得到一个预测值

以此类推，可以计算交叉验证的标准误差：

3-3)对于一个最优的平滑，也就是一个最优的λ值，s_cv应为其最小值。为了加快计算速度，引入H，帽子矩阵或平滑矩阵，其对角线上数值的平均数为

(I+αD‘D)H＝I

H矩阵的列可以通过平滑其所对应的单位矩阵找到。

基于惩罚最小二乘法的平滑的算法在Matlab中进行编码，如下：

m＝length(y)；

E＝speye(m)；

D＝diff(E,d)；

W＝spdiags(w,0,m,m)；

C＝chol(W+lambda*D’*D)；

z＝C\(C’\(w*y))；

步骤2)中，进行峰检测的条件有很多，比如信噪比、峰强度阈值、峰形、脊线、极大值、峰宽等等。本发明是使用信噪比和脊线来作为峰检测条件，用连续小波变换作为算法。连续小波变换是对信号函数与经过缩放与平移的小波母函数的乘积在整个时间域的积分，其公式如下：

S(t)是信号，a是缩放系数，b是位移系数。Ψ(t)是小波母函数，Ψ_a,b(t)即为经过缩放和平移的小波函数。结果C(a,b)是一个小波系数的二维矩阵(2D)。

由于小波系数反映了信号s和Ψ_a,b(t)间的相似程度，所以小波母函数需要具有拉曼光谱峰最基本的特点。本发明中选择了“墨西哥帽”函数作为小波母函数，如图4所示，其数学表达为：

当用此方法进行峰探测时，连续小波变换系数在每一次缩放后都在峰中心周围有一个极大值。极大值的值随着小波尺寸的增大而增大，直到小波尺寸与峰宽为最佳匹配时，之后极大值的值又慢慢变小。当我们把对连续小波变换系数的放大倍数作为第三维度放到连续小波变换二维系数图中时，在峰位置就会出现一条清晰的脊线。所以峰检测算法此时就包含三个步骤：

1.通过连接极大值来识别脊线。

脊线可以通过连接在各放缩级别的连续小波变换系数的极大值来实现。首先，找到在各放缩级别的极大值。然后使用一个滑动窗口，其大小和在此缩放级别的小波变换支持区域的大小成比例。接下来将这些极大值练成线，这些线就是识别出来的脊线。

2.识别出基于脊线对应的峰。

通过下面这三条法则来识别主峰：

a)在脊线中对应最大强度的缩放值，也就是那个和峰宽成比例的缩放值，应该被给定一定范围。

b)信噪比要大于某一阈值。信号强度定义为在一定缩放范围内，在脊线上最大的连续小波变换系数。噪声定义为当缩放系数等于1时，连续小波变换系数95％的分位数。

c)脊线的强度应该大于某一阈值。

3.优化峰参数。

为了精确的识别出峰，在进行完整体的峰检测后，可以对感兴趣的拉曼光谱区域以及使用一些特定的连续小波变换缩放系数来再次对拉曼光谱进行峰检测。此时的峰检测是通过附加计算完成的。

采用本发明所述处理方法进行处理后的拉曼谱图如图5、图6、图7、图8所示。

在某些实施例中，步骤2)还可估算出拉曼光谱图中峰的宽度，使用的算法是基于哈尔小波函数的微分运算，具体地：

5-1)根据哈尔小波函数的特点，一组信号的n次导数可以通过应用n次连续小波变换来实现。哈尔小波函数的数学表达式为：

5-2)峰宽估计步骤如下：

5-2-1)使用在峰检测中同样的缩放系数对此哈尔小波进行连续小波变换。二维连续小波变换系数以M*N的矩阵表示。

5-2-2)然后对此矩阵中所有值取绝对值。

5-2-3)对于在峰检测阶段检测到的每一个峰都有两个参数：峰指数和峰尺寸。二维连续小波变换矩阵中对应峰尺寸的行被用来从峰指数中寻找每个区域的极小值。

5-2-4)如果极小值不存在，那么峰的起点或终点就是三倍于其峰尺寸的最小值或下一个峰指数。如果极小值存在，那么峰起点或终点就是最近的那个极小值。

5-2-5)重复步骤5-2-3)～5-2-4)，直到得到所有的峰宽。

以上详细描述了本发明的较佳具体实施例。应当理解，本领域的普通技术人员无需创造性劳动就可以根据本发明的构思作出诸多修改和变化。因此，凡本技术领域中技术人员依本发明的构思在现有技术的基础上通过逻辑分析、推理或者有限的实验可以得到的技术方案，皆应在由权利要求书所确定的保护范围内。

Claims (5)

1.一种激光拉曼光谱解谱处理方法，其特征在于，包括以下步骤：

1)对拉曼光谱原始信号进行基于自适应迭代重加权惩罚最小二乘法的基线校正处理和基于惩罚最小二乘法的平滑处理，所述基于惩罚最小二乘法的平滑处理中，引入一帽子矩阵或平滑矩阵获得最优平滑；

2)对经所述步骤1)处理后的信号进行基于连续小波变换的峰检测，并估计拉曼光谱图中峰的宽度；

步骤1)中，对拉曼光谱信号进行基于惩罚最小二乘法的平滑处理是根据原始曲线分量，依次经稀疏矩阵、微分矩阵和平滑矩阵获得粗平滑向量，再经优化获得最终平滑向量，具体步骤如下：

3-1)将自适应迭代重加权惩罚最小二乘法中的加权系数去除，得到峰平滑的数学表达式：

(I+αD'D)z＝y

式中I为单位矩阵；D为微分矩阵；z为平滑后光谱的向量；y为原始曲线向量，α为平滑度，D'为D的导数；

3-2)依次将y中的一个数据取出，然后对剩下的数据进行平滑，此时对于取出的y_i得到一个预测值

以此类推，计算交叉验证的标准误差：

其中，m为向量长度；

3-3)对于一个最优的平滑，也就是一个最优的λ值，s_cv应为其最小值，引入帽子矩阵或平滑矩阵H，其对角线上数值的平均数为

(I+αD'D)H＝I

H矩阵的列通过平滑其所对应的单位矩阵找到；

所述拉曼光谱图中峰的宽度通过基于哈尔小波函数的微分运算估计获得，具体地：

5-1)根据哈尔小波函数的特点，一组信号的n次导数通过应用n次连续小波变换来实现，哈尔小波函数的数学表达式为：

5-2)峰宽估计：

5-2-1)使用在峰检测中同样的缩放系数对哈尔小波进行连续小波变换，二维连续小波变换系数以M*N的矩阵表示；

5-2-2)然后对所述矩阵中所有值取绝对值；

5-2-3)基于所述矩阵中对应峰尺寸的行从峰指数中寻找每个区域的极小值；

5-2-4)如果极小值不存在，那么峰的起点或终点就是三倍于其峰尺寸的最小值或下一个峰指数，如果极小值存在，那么峰起点或终点就是最近的那个极小值；

5-2-5)重复步骤5-2-3)～5-2-4)，直到得到所有的峰宽。

2.根据权利要求1所述的激光拉曼光谱解谱处理方法，其特征在于，所述基于自适应迭代重加权惩罚最小二乘法的基线校正处理中，所述自适应迭代重加权惩罚最小二乘法的表达式为：

式中，Q为原始基线与拟合基线保真度与粗糙度间的平衡，x＝[x₁,x₂,...,x_m]为原始信号向量，x_i为其中一个元素，z＝[z₁,z₂,...,z_m]为校正向量，z_i为其中一个元素，m为向量长度，λ为粗糙度系数，i、j为计数变量，t为迭代次数，w为权重向量，由自适应迭代方法获得，表达式为：

式中，d^t为在第t次迭代时原始信号向量与第t-1次迭代获得的校正向量的负值的系数，权重向量的初始值为1。

3.根据权利要求2所述的激光拉曼光谱解谱处理方法，其特征在于，所述基于自适应迭代重加权惩罚最小二乘法的迭代结束条件包括达到最大迭代次数或满足表达式|d^t|<0.001×|x|。

4.根据权利要求1所述的激光拉曼光谱解谱处理方法，其特征在于，所述基于连续小波变换的峰检测具体为：

201)通过连接在各放缩级别的连续小波变换系数的极大值识别脊线；

202)识别所述脊线对应的峰；

203)对所述峰的参数进行优化。

5.根据权利要求1所述的激光拉曼光谱解谱处理方法，其特征在于，所述步骤1)中，先进行所述基线校正处理，后进行所述平滑处理，或者，

先进行所述平滑处理，后进行所述基线校正处理。

Images