CN117056671A - 一种基于emd的拉曼光谱降噪方法 - Google Patents

Abstract

本发明公开了一种基于EMD的拉曼光谱降噪方法，包括以下步骤：步骤S1：获取含噪声的拉曼光谱数据；步骤S2：对含噪声的拉曼光谱数据进行经验模态分解；步骤S3：对经EMD分解的IMF噪声含量利用IMF的排列熵进行判别；步骤S4：检测n维本征模态分量拉曼峰信号；步骤S5：软阈值滤波去噪与分解峰边界处理；步骤S6：Savi tzky‑Golay滤波处理；步骤S7：信号重组得到降噪后拉曼光谱数据。本发明提供的基于EMD的拉曼光谱降噪方法，具有操作简单，拉曼特征峰信息高保真度等优点。

Description

一种基于EMD的拉曼光谱降噪方法

技术领域

本发明涉及拉曼光谱降噪技术领域，具体是一种基于EMD(Empirical ModeDecomposition，经验模态分解)的拉曼光谱降噪方法。

背景技术

现有的拉曼光谱降噪方法有两种降噪途径。

第一种是信号拟合降噪。Savitzky-Golay滤波和Whittaker滤波是在拉曼光谱降噪上应用较为广泛的信号拟合类的降噪方法。

Savitzky-Golay滤波通过拟合多项式来平滑光谱数据，通过这种平滑性，可以去除光谱中的高频噪声和突变，使数据更加平稳。

Whittaker滤波是基于最小二乘曲线拟合的思想，通过在滤波过程中平衡光滑度和数据拟合度来实现平滑效果。

第二种是信号分解降噪。信号分解降噪类方法试图将光谱中的有效信号与噪声信号分离。傅里叶变换(Fourier Transformation)阈值滤波、小波变换(Wavelet Transform)滤波和经验模态分解(Empirical Mode Decomposition)均属于信号分解降噪。

傅里叶变换阈值滤波是从频域的角度对信号进行分解。通过噪声的频域特点，构造特定的滤波器(通常为低通或带通滤波器)对频域信号进行滤波操作，达到选择性的增强或抑制信号的作用，最后对滤波后的频域信号进行逆傅里叶变换得到降噪后的拉曼光谱数据。

小波变换阈值滤波是使用选定的小波基函数，将信号分解为不同尺度和频率的小波基函数的线性组合，得到不同尺度下的小波系数。低频小波系数表示信号的平滑成分，高频小波系数则表示信号的细节成分。依据选择的小波基函数与分解后的小波系数选择合适的阈值进行滤波处理，最后将滤波后的小波系数重构得到降噪后的拉曼光谱数据。

经验模态分解是从信号的极值出发构建满足特定条件的本征模态分量(Intrinsic Mode Functions)分解信号。将含噪光谱数据进行EMD分解，得到n个本征模态分量(IMF)分量和一个残余(residual)分量，构建并应用一定的筛选规则，寻找并筛除高噪声占比的IMF，利用筛选后的所有分量重构信号，得到降噪后的拉曼光谱数据。

现有技术的缺点/不足：

尽管Savi tzky-Golay滤波在拉曼光谱降噪中是一种常用的方法，但是SG滤波仍存在一些不足。平滑参数的选择问题：对于不同的拉曼光谱信号，想要得到较为理想的降噪数据，通常需要依据信号特点选择不同的降噪参数，如滤波窗口和拟合阶数，不合适的参数选择会导致“过平滑”与“欠平滑”。“过平滑”是平滑过度导致的，虽然大幅度降低了环境及暗电流等随机噪声，但是导致拉曼峰形变形，从而影响到信号中的峰值和细节信息，相当于额外引入了另一种形式的噪声；“欠平滑”则是平滑效果太差导致的，滤波后的信号保留了过多的噪声，影响拉曼信号后续的使用；边界处理问题：由于SG滤波是基于滑动窗口的，故在两端边界小于选定的窗口的范围内，无法做到与其它范围内相同的处理。

Whittaker滤波同SG滤波一样属于信号拟合类降噪，同样存在SG滤波类似的问题。平滑参数选择问题：平滑度是一个影响降噪效果的重要参数，它控制着平滑度和数据拟合度之间的权衡。过大或过小的平滑度参数都可能导致滤波效果不理想。较大的平滑度参数可能会过度平滑信号，丧失细节信息，而较小的平滑度参数可能无法有效降噪；拟合阶数的选择：通常选择越高的拟合阶数，保留的信号细节就会越多，但保留的细节信息可能是有效信号，也有可能是噪声信号。拟合误差：Whi ttaker滤波通过最小化数据拟合项和平滑度约束项来实现滤波效果。然而，这种拟合过程可能会引入一定的拟合误差，尤其是当原始信号中存在复杂的噪声或非线性成分时。异常值处理问题：Whi ttaker滤波对异常值敏感。当存在极端值或离群点时，滤波结果可能会受到异常值的干扰，导致滤波效果不理想。

傅里叶变换阈值滤波和小波变换滤波都是尝试使用基函数取分解信号，不同点在于傅里叶变换使用正弦函数，而小波变换可以自由选择要使用的基函数，如Daubechies小波、Haar小波、Morlet小波等，并且可以使基函数以不同尺度和平移去分解，具有更高的自由性，与此同时也带来了更高的使用门槛与特定信号的特殊处理方法，模型的泛用性较差；傅里叶变换阈值滤波假设信号是周期性的。当信号不满足这个假设时，傅里叶变换可能导致频谱泄漏，即频谱成分在频域中被错误地分布到其他频率上。频谱泄漏可能会导致滤波结果中的误差和伪影。

现有的经验模态分解降噪对拉曼光谱的有效拉曼峰信号影响大，导致半峰宽、峰值强度、峰中心波长等产生畸变，严重降低了拉曼光谱中的有效信息含量，影响拉曼光谱的后续数据分析工作。现有的EMD降噪，会将EMD分解出的前k(k<n)个含噪声占比高的IMF分量直接抛弃，这虽然能滤除大部分高频噪声，但对于多个拉曼峰密集的情况，前k个IMF分量也会包含很多拉曼峰信号，这会一并将其滤除，从而导致信号缺失。由于EMD分解过程需要保证每个IMF满足零均值条件，是产生模态混叠的原因之一，而这一条件会导致一个拉曼峰信号被分解到多个IMF中，这也包含前k个IMF分量，直接滤除也会影响拉曼峰的峰值强度。同时在多个IMF滤除的情况下，可能会导致半峰宽畸变与峰中心波长偏移。

发明内容

本发明的目的在于提供一种基于EMD的拉曼光谱降噪方法，以解决现有技术中的问题。

为实现上述目的，本发明提供如下技术方案：一种基于EMD的拉曼光谱降噪方法，包括如下步骤：步骤S1：获取含噪声的拉曼光谱数据；步骤S2：对含噪声的拉曼光谱数据进行经验模态分解；步骤S3：对经EMD分解的IMF噪声含量利用IMF的排列熵进行判别；步骤S4：检测n维本征模态分量拉曼峰信号；步骤S5：软阈值滤波去噪与分解峰边界处理；步骤S6：Savitzky-Golay滤波处理；步骤S7：信号重组得到降噪后拉曼光谱数据。

进一步地，所述步骤S2对含噪声的拉曼光谱数据进行EMD分解，获得对应的n维本征模态分量和一维分解残余分量，将每个一维本征模态分量记为IMF_i,其中i＝1,2,3…n，将一维分解残余分量记为IMF_n+1。

进一步地，所述步骤S3利用式1计算每个一维本征模态分量和一维分解残余分量的排列熵P_i，结合特定的排列熵计算窗口为P_i设定一个合适的T_P，筛选计算得出每个P_i，若P_i＞T_P，转入步骤S4，否则IMF_i”＝IMF_i，继续处理下一个IMF分量；

进一步地，所述步骤S4包括：寻找IMF_i中所有峰值信息P_i,j,其中j＜L，结合当前分量的排列熵P_i与当前分量的最大峰值M_i＝max(P_i,j)，依据式2计算得出当前分量的有效拉曼峰峰值判定阈值R_i；

进一步地，所述步骤S5包括：依次取相邻的峰值信息P_i,j与P_i,j+1，与R_i进行对比判断，并且记录P_i,j对应的信号下标I_j与P_i,j+1；并将处理后的IMF_i重组为IMF_i'。

进一步地，所述步骤S6包括：对IMF_i'进行拟合窗口为5、拟合阶数为2的特定平滑参数的Savitzky-Golay滤波处理，得到IMF_i”。

进一步地，所述步骤S7将所有IMF_i”按照式3进行重组，得到降噪后的拉曼光谱数据S'；

进一步地，所述步骤S5对应的信号下标I_j+1分为三种情况处理：

a.|P_i,j|＞|P_i,j+1|＞R_i或|P_i,j+1|＞|P_i,j|＞R_i：此时IMF_i从I_j到I_j+1为有效拉曼峰的分解信息，直接保留该范围信号；

b.|P_i,j|＞R_i＞|P_i,j+1|或|P_i,j+1|＞R_i＞|P_i,j|：此时IMF_i从I_j到I_j+1的范围中存在部分有效拉曼峰的分解信息，结合IMF构造特性，在该范围内寻找零点下标I_z；|P_i,j|＞R_i＞|P_i,j+1|时，将IMF_i中I_j到I_z范围的信号置0；P_i,j+1|＞R_i＞|P_i,j|时，将IMF_i中I_z到I_j+1范围的信号置0；

c.R_i＞|P_i,j|＞|P_i,j+1|或R_i＞|P_i,j+1|＞|P_i,j|：此时IMF_i从I_j到I_j+1的范围内不存在有效拉曼峰的分解信息，全部做置0处理。

与现有技术相比，本发明的有益效果是：

1、本申请使用简单，无需对不同的拉曼光谱进行不同的降噪处理；

2、避免了模态混叠而导致的有效拉曼峰信息丢失问题；

3、在降低噪声的同时，很好的保持了拉曼峰细节信息，包括峰值、半峰宽、峰中心波长等。

附图说明

附图用来提供对本发明的进一步理解，并且构成说明书的一部分，与本发明的实施例一起用于解释本发明，并不构成对本发明的限制。在附图中：

图1是本发明的拉曼光谱降噪流程图；

图2是现有EMD分解筛除降噪、Savitzky-Golay滤波降噪与申请方法降噪结果对比图；

图3是本发明获取含噪声的对乙酰氨基酚拉曼光谱数据S图；

图4是本发明将RES记为IMF_n+1，分解后信号图；

图5是本发明IMF₁与IMF₁'的对比图；

图6是本发明IMF₁与IMF₁”对比图；

图7是本发明降噪后的拉曼光谱数据S'图。

具体实施方式

为使本发明实施方式的目的、技术方案和优点更加清楚，下面将结合本发明实施方式中的附图，对本发明实施方式中的技术方案进行清楚、完整地描述，显然，所描述的实施方式是本发明一部分实施方式，而不是全部的实施方式。基于本发明中的实施方式，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施方式，都属于本发明保护的范围。因此，以下对在附图中提供的本发明的实施方式的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施方式。

请参阅图1-7，本发明实施例中，提供一种基于EMD的拉曼光谱降噪方法，包括：S1：获取含噪声的拉曼光谱数据；S2：经验模态分解；S3：噪声含量判别；S4：n维本征模态分量(IMF)拉曼峰信号检测；S5：IMF软阈值滤波去噪与分解峰边界处理；S6：IMF Savitzky-Golay滤波处理；S7：信号重组得到降噪后拉曼光谱数据。

各步骤具体包括实现如下：

步骤S1：获取含噪声的对乙酰氨基酚拉曼光谱数据S，其信号长度为L；

步骤S2：对含噪声的拉曼光谱数据进行EMD分解，获得对应的n维本征模态分量(IMF)和一维分解残余分量(RES)，将每一维IMF记为IMF_i,其中i＝1,2,3…n，将RES记为IMF_n+1；

步骤S3：利用式1计算每个IMF及RES分量的排列熵P_i，结合特定的排列熵计算窗口为P_i设定一个合适的T_P，筛选计算得出的每个P_i，计算排列熵的时候，存在一个计算窗口的参数，不同的计算窗口会得到不同的排列熵取值范围，具体是PE∈[0,log₂ D！]，其中D是窗口大小。由于排列熵范围的会随排列熵的计算窗口的改变而发生变化，故对含噪声分量的判定阈值也应该随之改变，在该算法中，这个计算窗口和排列熵阈值是固定的。可选地，窗口数值为5，排列熵阈值取值为1.3。若P_i＞T_P，转入步骤S4，否则IMF_i”＝IMF_i，继续处理下一个IMF分量；

步骤S4：寻找IMF_i中所有峰值信息P_i,j,其中j＜L，结合当前分量的排列熵P_i与当前分量的最大峰值M_i＝max(P_i,j)，依据式2计算得出当前分量的有效拉曼峰峰值判定阈值R_i；

步骤S5：IMF软阈值滤波去噪是为了消除高频低幅值的噪声，这里是指有非常大的把握确定筛选出来的信号就是无意义的噪声信号，同时因为是直接消除这部分信号，故这里的阈值计算是相对保守的，防止将有效拉曼信号误判成了噪声，后面的Savitzky-Golay滤波，是对这一步保守策略的弥补，抑制可能被放行的相对较高幅值的噪声信号。分解峰边界处理具体包括依次取相邻的峰值信息P_i,j与P_i,j+1，与R_i进行对比判断，并且记录P_i,j对应的信号下标I_j与P_i,j+1；对应的信号下标I_j+1分为三种情况处理：

a.|P_i,j|＞|P_i,j+1|＞R_i或|P_i,j+1|＞|P_i,j|＞R_i：此时IMF_i从I_j到I_j+1为有效拉曼峰的分解信息，直接保留该范围信号；

b.|P_i,j|＞R_i＞|P_i,j+1|或|P_i,j+1|＞R_i＞|P_i,j|：此时IMF_i从I_j到I_j+1的范围中存在部分有效拉曼峰的分解信息，结合IMF构造特性，在该范围内寻找零点下标I_z。|P_i,j|＞R_i＞|P_i,j+1|时，将IMF_i中I_j到I_z范围的信号置0；P_i,j+1|＞R_i＞|P_i,j|时，将IMF_i中I_z到I_j+1范围的信号置0；

c.R_i＞|P_i,j|＞|P_i,j+1|或R_i＞|P_i,j+1|＞|P_i,j|：此时IMF_i从I_j到I_j+1的范围内不存在有效拉曼峰的分解信息，全部做置0处理；

将经过处理后的IMF_i重组为IMF_i',以IMF₁为例，IMF₁与IMF₁'的对比如图5，在去除IMF₁分量中高频且低于阈值的噪声信号的同时，保留特征峰的边界信号；

步骤S6：对IMF_i'进行拟合窗口为5、拟合阶数为2的特定平滑参数的Savitzky-Golay滤波处理，得到IMF_i”，以IMF₁为例，IMF₁与IMF₁”的对比如图6，显示单条IMF分量处理前后的数据变化，引入Savitzky-Golay滤波拟合分解后的有效拉曼特征峰，能有效抑制步骤S5放行的高频但幅值较高的噪声信号，并且提高特征峰与非有效信号过渡区域的信号平滑度；

步骤S7：将所有IMF_i”按照式3进行重组，得到降噪后的拉曼光谱数据S'；

所述经EMD分解的IMF噪声含量判别是利用IMF的排列熵进行划分。

对于含噪声IMF的有效拉曼峰信号判定是综合IMF的排列熵和IMF的峰值进行计算判定。

对于IMF中非有效拉曼特征峰信号的剔除，综合信号分解过滤噪声与有效拉曼特征峰信号的拟合保留。

通过对分解后的信号进行有效拉曼特征峰信号判定与对不同信号采取对应处理方法，以实现对拉曼光谱的无参、高特征峰信息降噪；无参是指使用该算法时不需要额外设定任何参数；高特征峰信息是指在降噪的同时，对于有效拉曼特征峰影响小，能有效的保留特征峰信息。

经验模态分解(EMD)：是一种信号处理技术，用于将非线性和非平稳信号分解成一系列称为本征模态分量(Intrins ic Mode Funct ions，简称IMF)的振动模式。EMD方法不需要先验假设或模型，可以适应不同类型和非线性、非平稳信号的降噪需求。EMD将信号分解成局部振动模式，有助于对局部特征进行处理和分析并且分解后的分量能够较好地保留信号的时域和频域特征，避免了传统方法中可能引入的伪迹和失真。但EMD在某些情况下可能会出现模态混叠现象，即一个特定信号被分解到两个或多个IMF之中，或一个IMF中表现出多种不同的时域和频域特征，难以清晰地区分和处理。

排列熵(Permutat ion Entropy)：是一种用于分析时间序列的非线性动力学特征的方法。通过比较顺序序列中D各连续值的所有D！种排列来计算，可以用来描述顺序序列中的复杂度，特别是在动力学或观测噪声的情况下。

Savi tzky-Golay滤波：主要思想是利用最小二乘多项式拟合来实现平滑滤波。它的原理是在滑动窗口内进行多项式拟合，通过在每个窗口内计算多项式拟合系数来实现平滑效果。这些系数是通过最小二乘法来确定的，使得拟合多项式在窗口内与原始数据的平方误差最小。然后，使用这些系数对窗口内的数据进行加权平均，得到滤波后的数值。

本发明的工作原理是：对于不同的拉曼光谱需要设定不同的降噪参数；降噪引起的拉曼峰发生畸变，导致有效信息丢失或错误；EMD分解降噪在拉曼光谱上由于模态混叠而导致的拉曼峰分解信息丢失。本发明使用简单，无需对不同的拉曼光谱进行不同的降噪处理；避免了模态混叠而导致的有效拉曼峰信息丢失问题，图2中对比现有EMD分解筛除降噪与本发明申请方法降噪结果体现出该优点；在降低噪声的同时，很好的保持了拉曼峰细节信息，包括峰值、半峰宽、峰中心波长等，图2中对比Savi tzky-Golay滤波降噪与本发明申请方法降噪结果体现出该优点。

最后应说明的是：以上所述仅为本发明的优选实施例而已，并不用于限制本发明，尽管参照前述实施例对本发明进行了详细的说明，对于本领域的技术人员来说，其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (8)

1.一种基于EMD的拉曼光谱降噪方法，其特征在于，包括以下步骤：

步骤S1：获取含噪声的拉曼光谱数据；

步骤S2：对含噪声的拉曼光谱数据进行经验模态分解；

步骤S3：对经EMD分解的IMF噪声含量利用IMF的排列熵进行判别；

步骤S4：检测n维本征模态分量拉曼峰信号；

步骤S5：软阈值滤波去噪与分解峰边界处理；

步骤S6：Savitzky-Golay滤波处理；

步骤S7：信号重组得到降噪后拉曼光谱数据。

2.根据权利要求1所述的一种基于EMD的拉曼光谱降噪方法，其特征在于，所述步骤S2对含噪声的拉曼光谱数据进行EMD分解，获得对应的n维本征模态分量和一维分解残余分量，将每个一维本征模态分量记为IMF_i,其中i＝1,2,3…n，将一维分解残余分量记为IMF_n+1。

3.根据权利要求2所述的一种基于EMD的拉曼光谱降噪方法，其特征在于，所述步骤S3利用式1计算每个一维本征模态分量和一维分解残余分量的排列熵P_i，结合特定的排列熵计算窗口为P_i设定一个合适的T_P，筛选计算得出每个P_i，若P_i＞T_P，转入步骤S4，否则IMF_i”＝IMF_i，继续处理下一个IMF分量；

4.根据权利要求2所述的一种基于EMD的拉曼光谱降噪方法，其特征在于，所述步骤S4包括：寻找IMF_i中所有峰值信息P_i,j,其中j＜L，结合当前分量的排列熵P_i与当前分量的最大峰值M_i＝max(P_i,j)，依据式2计算得出当前分量的有效拉曼峰峰值判定阈值R_i；

5.根据权利要求2所述的一种基于EMD的拉曼光谱降噪方法，其特征在于，所述步骤S5包括：依次取相邻的峰值信息P_i,j与P_i,j+1，与R_i进行对比判断，并且记录P_i,j对应的信号下标I_j与P_i,j+1；并将处理后的IMF_i重组为IMF_i'。

6.根据权利要求2所述的一种基于EMD的拉曼光谱降噪方法，其特征在于，所述步骤S6包括：对IMF_i'进行拟合窗口为5、拟合阶数为2的特定平滑参数的Savitzky-Golay滤波处理，得到IMF_i”。

7.根据权利要求3所述的一种基于EMD的拉曼光谱降噪方法，其特征在于，所述步骤S7将所有IMF_i”按照式3进行重组，得到降噪后的拉曼光谱数据S'；

8.根据权利要求5所述的一种基于EMD的拉曼光谱降噪方法，其特征在于：所述步骤S5对应的信号下标I_j+1分为三种情况处理：

a.|P_i,j|＞|P_i,j+1|＞R_i或|P_i,j+1|＞|P_i,j|＞R_i：此时IMF_i从I_j到I_j+1为有效拉曼峰的分解信息，直接保留该范围信号；

b.|P_i,j|＞R_i＞|P_i,j+1|或|P_i,j+1|＞R_i＞|P_i,j|：此时IMF_i从I_j到I_j+1的范围中存在部分有效拉曼峰的分解信息，结合IMF构造特性，在该范围内寻找零点下标I_z；|P_i,j|＞R_i＞|P_i,j+1|时，将IMF_i中I_j到I_z范围的信号置0；|P_i,j+1|＞R_i＞|P_i,j|时，将IMF_i中I_z到I_j+1范围的信号置0；

c.R_i＞|P_i,j|＞|P_i,j+1|或R_i＞|P_i,j+1|＞|P_i,j|：此时IMF_i从I_j到I_j+1的范围内不存在有效拉曼峰的分解信息，全部做置0处理。