CN117216475A - 一种脑电信号去噪方法及系统 - Google Patents

Abstract

本发明提出一种脑电信号去噪方法及系统，涉及脑电信号去噪的技术领域，方法包括：首先将获取的脑电信号代入变分模态分解模型中，对变分模态分解模型进行优化，得到优化后的变分模态分解模型；然后优化后的变分模态分解模型对脑电信号进行变分模态分解成模态分量；接着计算模态分量和原始脑电信号之间的相关系数，根据所述相关系数判断出有效模态分量，将有效模态分量进行重构，进一步对重构信号进行小波阈值去噪，得到经过阈值处理的小波系数并小波重构，获得去噪干净的脑电信号。本发明有效提高了信号分解质量和信号的去噪效果。

Description

一种脑电信号去噪方法及系统

技术领域

本发明涉及脑电信号去噪的技术领域，特别涉及一种脑电信号去噪方法及系统。

背景技术

脑电信号是一种记录大脑电活动的非侵入性生物电信号，具有重要的临床和研究价值，在生物医疗领域得到广泛应用，但由于幅值较小，容易受到伪迹的干扰，特别是眼电伪迹和肌电伪迹等噪声，它们很难直接消除，这些噪声的存在会对实验和研究的结果产生影响，因此研究脑电信号的去噪方法具有重要意义。

目前，研究脑电信号的去噪方法主要有基于经验模态分解、变分模态分解和小波阈值去噪三种方法，经验模态分解是一种适于非线性和非平稳信号信号分解方法，可以根据信号的特性进行分解和去噪，无需假设或模型拟合，但在信号分解过程中，不同频率成分的模态函数相互重叠，难以准确分离，导致模态混叠，同时由于缺乏固定的边界处理策略，信号的边界处可能会出现振荡或不稳定的模态函数，导致边界效应；变分模态分解在分解过程中采用优化算法，通过解决优化问题来确定每个模态函数的中心频率和带宽，能够自动确定分解的模态函数的个数，并且具有自适应的带宽，可以更好地捕捉信号的局部特征，但需要预设参数来控制信号分解的精细程度和噪声抑制能力，若错误地选择参数可能导致信号分解的质量下降，出现模态过分解或模态混叠现象，无法充分去除噪声或保留感兴趣的信号成分；小波阈值去噪利用小波变换将信号转换到时频域，根据设定的阈值准则选择保留或丢弃小波系数，以实现噪声抑制，但在在低信噪比的情况下，对脑电信号进行小波阈值去噪可能会出现信号的伪影或失真，影响后续脑电信号分析的可靠性。

现有专利文献公开了一种信号去噪方法，获取含噪电能质量信号；选取排列熵作为遗传算法的自适应度函数，通过遗传算法调用变分模态分解，对变分模态分解的惩罚因子α与分解模态数k进行迭代寻优，确定最优参数；使用变分模态分解将信号数据分解为k个模态分量，通过相关系数确定有效模态分量与噪声模态分量；对于改进小波阈值，提出参数可调的阈值函数，并且将小波能量熵的概念引入阈值函数中；使用改进小波阈值对噪声模态分量作去噪处理，并选择有效模态分量与去噪处理后的噪声模态分量进行重构，得到消噪后的电能质量扰动信号，但在面对复杂信号时，利用遗传算法对变分模态分解的参数进行迭代寻优，无法全面或精细化探索解空间中的多个局部最优解，全局搜索性能和局部搜索能力差，而且遗传算法的迭代寻优过程是无法自适应调整系数，导致信号分解质量和去噪效果差。

发明内容

为解决当前信号去噪方法存在信号分解质量和去噪效果差的问题，本发明提出一种脑电信号去噪方法及系统，采用了蜣螂优化算法和变分模态分解联合小波阈值去噪对脑电信号进行去噪处理，有效提高了信号分解质量和信号的去噪效果。

为了达到上述技术效果，本发明的技术方案如下：

一种脑电信号去噪方法，包括以下步骤：

S1.获取脑电信号，将脑电信号代入变分模态分解模型中；

S2.对变分模态分解模型进行优化，得到优化后的变分模态分解模型；

S3.基于优化后的变分模态分解模型，对脑电信号进行变分模态分解，得到K个模态分量；

S4.计算K个模态分量中的模态分量与原始脑电信号之间的相关系数，根据所述相关系数判断K个模态分量中的模态分量是否为噪声模态分量，若是，则去除噪声模态分量；否则，记该模态分量为有效模态分量，将有效模态分量进行重构，得到重构信号；

S5.对重构信号进行小波阈值去噪，得到经过阈值处理的小波系数并小波重构，获得去噪干净的脑电信号。

优选地，对变分模态分解模型进行优化，具体过程为：将最小排列熵作为蜣螂优化算法的适应度函数，利用蜣螂优化算法，优化变分模态分解模型中的调节参数，得到最佳参数组合；最小排列熵的计算过程如下：

最小排列熵的计算过程如下：

记脑电信号的时间序列经变分模态分解得到的其中一个模态分量为{X(i),i＝1,2,…,N}，N为正整数，对其进行相空间重构,得到重构矩阵Y如下：

其中，d表示嵌入维数，τ表示延时因子，M表示重构矩阵Y分量的个数；将重构矩阵Y中的第j个分量(x(j),x(j+τ),…,x(j+(d-1)τ)按照数值的大小进行升序排列，得到符号序列S(q)＝(j₁,j₂,…,j_d)；其中q＝1,2,…,r，r表示正整数，r≤d！，j₁,j₂,…,j_d表示各元素在原重构矩阵分量中的索引号；重构矩阵中的每个分量是d维空间，映射到d维符号序列S(q)，共有d！种排列方式；

计算每一种d维符号序列出现的概率{P1,P2,…,Pd}，基于计算出的概率，计算所述时间序列的排列熵H_p(d)如下：

其中，P_j表示j维符号序列出现的概率；排列熵的值越小，说明时间序列越简单、越规则，反之，排列熵的值越复杂，越不规则；将最小值的排列熵作为最小排列熵。

优选地，利用蜣螂优化算法，优化变分模态分解模型中的调节参数，具体过程如下：

S21.初始化蜣螂种群大小和变分模态分解的调节参数范围；

S22.利用排列熵作为适应度函数，计算蜣螂种群的适应度值；

S23.通过模拟蜣螂行为，进行蜣螂的位置迭代更新；

S24.判断更新后蜣螂的位置是否位于原蜣螂的位置，若是，利用更新后蜣螂的位置代替原蜣螂的位置，执行步骤S25；否则，保持原蜣螂的位置不变，执行步骤S25；

S25.判断当前迭代次数是否小于最大迭代次数，若是，则返回步骤S22，否则，将当前蜣螂的位置作为最优位置，当前蜣螂种群的适应度值为最优适应度值，输出最优位置和最优适应度值对应的调节参数为最佳参数组合。

优选地，所述调节参数为分解个数K和惩罚因子α。

优选地，所述蜣螂行为包括滚球行为、跳舞行为、卵球繁殖行为、觅食行为和偷窃行为。

优选地，所述通过模拟蜣螂行为，进行蜣螂的位置迭代更新，具体更新过程为：

利用滚球行为进行蜣螂的位置迭代更新，负责滚球的蜣螂的位置更新如下：

x_i(t+1)＝x_i(t)+α×k×x_i(t-1)+b×Δx

△x＝|x_i(t)-X^w|

其中，x_i表示第i个蜣螂的位置信息，t表示当前迭代次数，α∈(0,1)表示一个随机数，k代表一个常数，表示偏转系数，b表示自然系数，b取值为-1或1，X^w表示全局最差位置，△x表示光强的变化；

当遇到障碍物时，蜣螂利用跳舞行为重新定位方向，此时蜣螂的位置更新如下：

x_i(t+1)＝x_i(t)+tan(θ)|x_i(t)-x_i(t-1)|

tan(θ)表示正切函数，θ表示偏转系角度，当θ等于0或π或2π时，不更改蜣螂的位置。

当蜣螂执行卵球繁殖行为时，利用边界选择策略来模拟雌性蜣螂产卵的区域，边界选择策略定义为：

Lb^*＝max(X^*×(1-R),Lb)

Ub^*＝min(X^*×(1+R),Ub)

其中，Lb^*表示蜣螂产卵时区域的下界，Ub^*表示蜣螂产卵时区域的下界，表示X^*当前最佳位置，Lb表示优化问题的下界，Ub表示优化问题的上界，R表示惯性权值，其计算表达式为：

R＝1-t/T_max

其中，T_max表示最大迭代次数；

在迭代过程中，蜣螂产卵的卵球的位置随着产卵的区域动态变化：

B_i(t+1)＝X^*+b₁×(B_i(t)-Lb^*)+b₂×(B_i(t)-Ub^*)

其中，B_i表示第i个卵球的位置，b₁和b₂表示1×D的随机向量，D表示优化问题的维数；

当小蜣螂进行觅食行为时，建立最优觅食区域来引导小蜣螂进行觅食，其中最优觅食区域定义为：

Lb^b＝max(X^b×(1-R),Lb)

Ub^b＝min(X^b×(1+R),Ub)

其中，X^b表示全局最佳位置，Lb^b表示最佳觅食区域的下限，Ub^b表示最佳觅食区域的上限；小蜣螂的位置更新为：

x_i(t+1)＝x_i(t)+C₁×(x_i(t)-Lb^b)+C₂×(x_i(t)-Ub^b)

其中，C₁表示服从正态分布的随机数；C₂∈(0,1)表示一个随机向量；

当蜣螂具有偷盗行为时，具有偷盗行为的蜣螂的位置更新如下：

x_i(t+1)＝X^b+S×g×(|x_i(t)-X^*|+|x_i(t)-X^b|)

其中，X^b表示争夺食物的最佳地点，S表示一个常数，g表示服从均值为0、方差为1的正态分布的随机向量。

优选地，基于优化后的变分模态分解模型，对脑电信号f进行变分模态分解，得到K个具有中心频率ω_k的有限带宽模态分量，具体过程包括：

S31.对脑电信号f进行希尔伯特变换，得到脑电信号f对应的解析信号，对解析信号进行单边谱处理，得到每个模态分量u_k(t)的单边谱：

其中，t表示时间，δ(t)表示关于时间t的狄拉克分布，j表示复数，k表示模态分量个数；

S32.将S31得到的单边谱与预估中心频率混合，调制出各模态信号带宽，各模态信号带宽受约束的变分表达式如下：

S33.基于步骤S32的变分表达式，构建目标约束函数：

其中，u_k表示本征模态分量，w_k表示各个模态分量的中心频率，s.t表示约束于，f为原始脑电信号，表示对模态分量个数k求偏导；

S34.引入二次惩罚因子α和拉格朗日乘子λ，构造增广拉格朗日函数L，将目标约束函数的约束变分问题转变成非约束变分问题，L的具体表达式如下：

其中，{u_k}表示本征模态分量u_k的集合，{w_k}表示各个模态分量的中心频率的集合，λ表示拉格朗日乘子算子，f(t)表示关于时间t的脑电信号，λ(t)表示关于时间t的拉格朗日乘子算子，<λ(t),f(t)-∑_ku_k(t)>表示拉格朗日乘子项；

利用交替方向乘子算法，得到更新后的各个模态分量、中心频率和拉格朗日乘子，更新后的模态分量的计算表达式为：

其中，表示更新后的本征模态分量，/>表示傅里叶变换的结果，/>表示本征模式分量在傅里叶域中的迭代值，/>表示更新后的拉格朗日乘子算子，w表示频率参数，k表示模态分量的个数；

更新后的中心频率的计算表达式为：

其中，为更新后的中心频率；

更新后的拉格朗日乘子的计算表达式为：

其中，表示更新后的拉格朗日乘子，/>表示拉格朗日乘子算子值在傅里叶域中的迭代值，τ表示更新参数。

S35.判断分解得到的模态分量是否满足下式：

其中，ε为判断阈值；若是，则循环分解模态分量完成；否则，返回步骤S34。

优选地，所述计算K个模态分量中的模态分量与原始脑电信号之间的相关系数，具体包括：

S41.获取原始脑电信号的自相关函数R_ss和各模态分量的互相关函数R_J；

S42.基于R_ss和R_J，计算模态分量与原始脑电信号之间的相关系数：

其中，r_J表示模态分量与原始脑电信号之间的相关系数，J表示各模态分量的序号，r_J的值越大，模态分量与原始脑电信号之间的相关性越强，相反，r_J的值越小，模态分量与原始脑电信号之间的相关性越弱。

优选地，从计算得出的模态分量与原始脑电信号之间的相关系数中获取数值最大的相关系数，基于数值最大的相关系数，计算临界值ρ：

其中，表示数值最大的相关系数；

将计算出的相关系数大于或等于ρ的模态分量记为有效模态分量，将计算出的相关系数小于ρ的模态分量记为噪声模态分量。

本发明还提出了一种脑电信号去噪系统，所述系统包括：

信号获取模块，用于获取脑电信号，将脑电信号代入变分模态分解模型中；

优化模块，用于对变分模态分解模型进行优化，得到优化后的变分模态分解模型；

分解模块，用于根据优化后的变分模态分解模型，对脑电信号进行变分模态分解，得到K个模态分量；

计算模块，用于计算K个模态分量中的模态分量与原始脑电信号之间的相关系数，根据所述相关系数判断K个模态分量中的模态分量是否为噪声模态分量，若是，则去除噪声模态分量；否则，记该模态分量为有效模态分量，将有效模态分量进行重构，得到重构信号；

去噪模块，用于对重构信号进行小波阈值去噪，得到经过阈值处理的小波系数并小波重构，获得去噪干净的脑电信号。

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

本发明提出一种脑电信号去噪方法及系统，首先将获取的脑电信号代入变分模态分解模型中，其次对变分模态分解模型进行优化，得到优化后的变分模态分解模型，使得优化后的变分模态分解模型能够实现自适应调整，再结合优化后的变分模态分解模型对脑电信号进行变分模态分解成模态分量，提升了脑电信号的信噪比，减弱脑电信号中因存在大量的噪声对小波阈值算法准确性和可靠性的影响；接着计算模态分量和原始脑电信号之间的相关系数，根据所述相关系数判断出有效模态分量，将有效模态分量进行重构，进一步对重构信号进行小波阈值去噪，得到经过阈值处理的小波系数并小波重构，获得去噪干净的脑电信号，减少了计算复杂度，有效提高了信号分解质量和信号的去噪效果。

附图说明

图1表示本发明实施例中提出的一种脑电信号去噪方法的流程示意图；

图2表示本发明实施例中提出的蜣螂优化算法的流程示意图；

图3表示本发明实施例中提出的一种脑电信号去噪的结构图。

具体实施方式

附图仅用于示例性说明，不能理解为对本专利的限制；

为了更好地说明本实施例，附图某些部位会有省略、放大或缩小，并不代表实际尺寸，“上”“下”等部位方向的描述非对本专利的限制；

对于本领域技术人员来说，附图中某些公知内容说明可能省略是可以理解的；

附图中描述位置关系的用于仅用于示例性说明，不能理解为对本专利的限制；

下面结合附图和实施例对本发明的技术方案做进一步的说明。

实施例1

如图1所示，本实施例提出了一种脑电信号去噪方法，包括以下步骤：

S1.获取脑电信号，将脑电信号代入变分模态分解中；

S2.对变分模态分解模型进行优化，得到优化后的变分模态分解模型；

在步骤S2中，对变分模态分解模型进行优化，具体过程为：将最小排列熵作为蜣螂优化算法的适应度函数，利用蜣螂优化算法，优化变分模态分解模型中的调节参数，得到最佳参数组合；最小排列熵的计算过程如下：

记脑电信号的时间序列经变分模态分解得到的其中一个模态分量为{X(i),i＝1,2,…,N}，N为正整数，对其进行相空间重构,得到重构矩阵Y如下：

其中，d表示嵌入维数，τ表示延时因子，M表示重构矩阵Y分量的个数；将重构矩阵Y中的第j个分量(x(j),x(j+τ),…,x(j+(d-1)τ)按照数值的大小进行升序排列，得到符号序列S(q)＝(j₁,j₂,…,j_d)；其中q＝1,2,…,r，r表示正整数，r≤d！，j₁,j₂,…,j_d表示各元素在原重构矩阵分量中的索引号；重构矩阵中的每个分量是d维空间，映射到d维符号序列S(q)，共有d！种排列方式；

计算每一种d维符号序列出现的概率{P1,P2,…,Pd}，基于计算出的概率，计算所述时间序列的排列熵H_p(d)如下：

其中，P_j表示j维符号序列出现的概率；排列熵的值越小，说明时间序列越简单、越规则，反之，排列熵的值越复杂，越不规则；将最小值的排列熵作为最小排列熵；利用蜣螂优化算法，优化变分模态分解模型中的调节参数，具体过程如下：

S21.初始化蜣螂种群大小和变分模态分解的调节参数范围；

在步骤S21中，设置蜣螂优化算法的初始化参数，蜣螂种群大小设置为10，最大迭代次数为20，初始化分解个数K和惩罚因子α的范围如下：

K＝[5,10],α＝[500,2000]

S22.利用排列熵作为适应度函数，计算蜣螂种群的适应度值；

S23.通过模拟蜣螂行为，进行蜣螂的位置迭代更新；

S24.判断更新后蜣螂的位置是否位于原蜣螂的位置，若是，利用更新后蜣螂的位置代替原蜣螂的位置，执行步骤S25；否则，保持原蜣螂的位置不变，执行步骤S25；

S25.判断当前迭代次数是否小于最大迭代次数，若是，则返回步骤S22，否则，将当前蜣螂的位置作为最优位置，当前蜣螂种群的适应度值为最优适应度值，输出最优位置和最优适应度值对应的调节参数为最佳参数组合；

在步骤S25中，对获取的脑电信号的时间序列进行变分模态分解后，得到K个模态分量，分别计算K个模态分量的排列熵值，并求和，得到K个模态分量总的排列熵值，将其作为优化算法的适应度值；其中最优适应度值的获取过程为：分别计算每次位置更新后相应的适应度值之和，比较所有适应度值的大小，将最小的适应度值作为最优适应度值，最优适应度值对应的分解个数K和惩罚因子α作为最优参数组合；

S3.基于优化后的变分模态分解模型，对脑电信号进行变分模态分解，得到K个模态分量；

S4.计算K个模态分量中的模态分量与原始脑电信号之间的相关系数，根据所述相关系数判断K个模态分量中的模态分量是否为噪声模态分量，若是，则去除噪声模态分量；否则，记该模态分量为有效模态分量，将有效模态分量进行重构，得到重构信号；

S5.对重构信号进行小波阈值去噪，得到经过阈值处理的小波系数并小波重构，获得去噪干净的脑电信号。

在本实施例中，首先将获取的脑电信号代入变分模态分解中，其次将最小排列熵作为蜣螂优化算法的适应度函数，利用蜣螂优化算法，优化变分模态分解的调节参数，得到最佳参数组合，利用蜣螂优化算法与变分模态分解结合不仅可以抑制模态混叠，还可以自适应地选取调节参数；再结合最佳参数组合对脑电信号进行变分模态分解成模态分量，提升了脑电信号的信噪比，减弱脑电信号中因存在大量的噪声对小波阈值算法准确性和可靠性的影响；接着计算模态分量和原始脑电信号之间的相关系数，根据所述相关系数判断出有效模态分量，将有效模态分量进行重构，进一步对重构信号进行小波阈值去噪，得到经过阈值处理的小波系数并小波重构，获得去噪干净的脑电信号，减少了计算复杂度，有效提高了信号分解质量和信号的去噪效果。

实施例2

本实施例进一步对蜣螂的位置迭代更新进行说明，在蜣螂的位置迭代更新过程中，所述蜣螂行为包括滚球行为、跳舞行为、卵球繁殖行为、觅食行为和偷窃行为；

所述通过模拟蜣螂行为，进行蜣螂的位置迭代更新，具体更新过程为：

利用滚球行为进行蜣螂的位置迭代更新，负责滚球的蜣螂的位置更新如下：

x_i(t+1)＝x_i(t)+α×k×x_i(t-1)+b×Δx

△x＝|x_i(t)-X^w|

其中，x_i表示第i个蜣螂的位置信息，t表示当前迭代次数，α∈(0,1)表示一个随机数，k代表一个常数，表示偏转系数，b表示自然系数，b取值为-1或1，X^w表示全局最差位置，△x表示光强的变化；

当遇到障碍物时，蜣螂利用跳舞行为重新定位方向，此时蜣螂的位置更新如下：

x_i(t+1)＝x_i(t)+tan(θ)|x_i(t)-x_i(t-1)|

tan(θ)表示正切函数，θ表示偏转系角度，当θ等于0或π或2π时，不更改蜣螂的位置。

当蜣螂执行卵球繁殖行为时，利用边界选择策略来模拟雌性蜣螂产卵的区域，边界选择策略定义为：

Lb^*＝max(X^*×(1-R),Lb)

Ub^*＝min(X^*×(1+R),Ub)

其中，Lb^*表示蜣螂产卵时区域的下界，Ub^*表示蜣螂产卵时区域的下界，表示X^*当前最佳位置，Lb表示优化问题的下界，Ub表示优化问题的上界，R表示惯性权值，其计算表达式为：

R＝1-t/T_max

其中，T_max表示最大迭代次数；

在迭代过程中，蜣螂产卵的卵球的位置随着产卵的区域动态变化：

B_i(t+1)＝X^*+b₁×(B_i(t)-Lb^*)+b₂×(B_i(t)-Ub^*)

其中，B_i表示第i个卵球的位置，b₁和b₂表示1×D的随机向量，D表示优化问题的维数；

当小蜣螂进行觅食行为时，建立最优觅食区域来引导小蜣螂进行觅食，其中最优觅食区域定义为：

Lb^b＝max(X^b×(1-R),Lb)

Ub^b＝min(X^b×(1+R),Ub)

其中，X^b表示全局最佳位置，Lb^b表示最佳觅食区域的下限，Ub^b表示最佳觅食区域的上限；小蜣螂的位置更新为：

x_i(t+1)＝x_i(t)+C₁×(c_i(t)-Lb^b)+C₂×(x_i(t)-Ub^b)

其中，C₁表示服从正态分布的随机数；C₂∈(0,1)表示一个随机向量；

当蜣螂具有偷盗行为时，具有偷盗行为的蜣螂的位置更新如下：

x_i(t+1)＝X^b+S×g×(|x_i(t)-X^*|+|x_i(t)-X^b|)

其中，X^b表示争夺食物的最佳地点，S表示一个常数，g表示服从均值为0、方差为1的正态分布的随机向量。

实施例3

本实施例对脑电信号进行变分模态分解的过程做进一步说明，基于优化后的变分模态分解模型，对脑电信号f进行变分模态分解，得到K个具有中心频率ω_k的有限带宽模态分量，具体过程包括：

S31.对脑电信号f进行希尔伯特变换，得到脑电信号f对应的解析信号，对解析信号进行单边谱处理，得到每个模态分量u_k(t)的单边谱：

其中，t表示时间，δ(t)表示关于时间t的狄拉克分布，j表示复数，k表示模态分量个数；

S32.将S31得到的单边谱与预估中心频率混合，调制出各模态信号带宽，各模态信号带宽受约束的变分表达式如下：

S33.基于步骤S32的变分表达式，构建目标约束函数：

其中，u_k表示本征模态分量，w_k表示各个模态分量的中心频率，s.t表示约束于，f为原始脑电信号，表示对模态分量个数k求偏导；

S34.引入二次惩罚因子α和拉格朗日乘子λ，构造增广拉格朗日函数L，将目标约束函数的约束变分问题转变成非约束变分问题，L的具体表达式如下：

其中，{u_k}表示本征模态分量u_k的集合，{w_k}表示各个模态分量的中心频率的集合，λ表示拉格朗日乘子算子，f(t)表示关于时间t的脑电信号，λ(t)表示关于时间t的拉格朗日乘子算子，<λ(t),f(t)-∑_ku_k(t)>表示拉格朗日乘子项；

利用交替方向乘子算法，得到更新后的各个模态分量、中心频率和拉格朗日乘子，更新后的模态分量的计算表达式为：

其中，表示更新后的本征模态分量，/>表示傅里叶变换的结果，/>表示本征模式分量在傅里叶域中的迭代值，/>表示更新后的拉格朗日乘子算子，w表示频率参数，k表示模态分量的个数；

更新后的中心频率的计算表达式为：

其中，为更新后的中心频率；

更新后的拉格朗日乘子的计算表达式为：

其中，表示更新后的拉格朗日乘子，/>表示拉格朗日乘子算子值在傅里叶域中的迭代值，τ表示更新参数。

S35.判断分解得到的模态分量是否满足下式：

其中，ε为判断阈值；若是，则循环分解模态分量完成；否则，返回步骤S34。

本实施例还进一步计算K个模态分量中的模态分量与原始脑电信号之间的相关系数，具体包括：

S41.获取原始脑电信号的自相关函数R_ss和各模态分量的互相关函数R_J；

在步骤S41中，R_ss和R_J是利用下式计算得出，具体计算公式如下：

S42.基于R_ss和R_J，计算模态分量与原始脑电信号之间的相关系数：

其中，r_J表示模态分量与原始脑电信号之间的相关系数，J表示各模态分量的序号，r_j的值越大，模态分量与原始脑电信号之间的相关性越强，相反，r_J的值越小，模态分量与原始脑电信号之间的相关性越弱。

从计算得出的模态分量与原始脑电信号之间的相关系数中获取数值最大的相关系数，基于数值最大的相关系数，计算临界值ρ：

其中，表示数值最大的相关系数；

将计算出的相关系数大于或等于ρ的模态分量记为有效模态分量，将计算出的相关系数小于ρ的模态分量记为噪声模态分量。

将有效模态分量进行重构，得到重构信号y(t)；选择合适的小波基，确定分解层数，对重构信号进行小波分解，得到分解后的各尺度小波分解系数；选择合适的阈值，对小波系数进行阈值处理；根据阈值处理后的小波系数，进行小波重构，得到去噪后的信号；在小波阈值去噪的过程中：

选择db3作为小波基函数，分解层数为5层，对重构信号y(t)进行小波变换。

选取sqtwolog阈值函数，阈值为对每层分解层的小波系数进行小波阈值去噪。

其中，w_j,k表示小波分解系数；W_j,k表示阈值处理后的小波系数，α和β均表示小波阈值参数，均取1，对处理后的N层系数进行小波重构得到去噪后的信号。

实施例4

参见图3，本实施例提出了一种脑电信号去噪系统300，所述系统包括：

信号获取模块301，用于获取脑电信号，将脑电信号代入变分模态分解模型中；

优化模块302，用于对变分模态分解模型进行优化，得到优化后的变分模态分解模型；

分解模块303，用于根据优化后的变分模态分解模型，对脑电信号进行变分模态分解，得到K个模态分量；

计算模块304，用于计算K个模态分量中的模态分量与原始脑电信号之间的相关系数，根据所述相关系数判断K个模态分量中的模态分量是否为噪声模态分量，若是，则去除噪声模态分量；否则，记该模态分量为有效模态分量，将有效模态分量进行重构，得到重构信号；

去噪模块305，用于对重构信号进行小波阈值去噪，得到经过阈值处理的小波系数并小波重构，获得去噪干净的脑电信号。

在本实施例中，首先将获取的脑电信号代入变分模态分解模型中，其次对变分模态分解模型进行优化，得到优化后的变分模态分解模型，使得优化后的变分模态分解模型能够实现自适应调整，再结合优化后的变分模态分解模型对脑电信号进行变分模态分解成模态分量，提升了脑电信号的信噪比，减弱脑电信号中因存在大量的噪声对小波阈值算法准确性和可靠性的影响；接着计算模态分量和原始脑电信号之间的相关系数，根据所述相关系数判断出有效模态分量，将有效模态分量进行重构，进一步对重构信号进行小波阈值去噪，得到经过阈值处理的小波系数并小波重构，获得去噪干净的脑电信号，减少了计算复杂度，有效提高了信号分解质量和信号的去噪效果。

显然，本发明的上述实施例仅是为清楚地说明本发明所作的举例，而并非是对本发明的实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动。这里无需也无法对所有的实施方式予以穷举。凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明权利要求的保护范围之内。

Claims (10)

1.一种脑电信号去噪方法，其特征在于，包括以下步骤：

S1.获取脑电信号，将脑电信号代入变分模态分解模型中；

S2.对变分模态分解模型进行优化，得到优化后的变分模态分解模型；

S3.基于优化后的变分模态分解模型，对脑电信号进行变分模态分解，得到K个模态分量；

S4.计算K个模态分量中的模态分量与原始脑电信号之间的相关系数，根据所述相关系数判断K个模态分量中的模态分量是否为噪声模态分量，若是，则去除噪声模态分量；否则，记该模态分量为有效模态分量，将有效模态分量进行重构，得到重构信号；

S5.对重构信号进行小波阈值去噪，得到经过阈值处理的小波系数并小波重构，获得去噪干净的脑电信号。

2.根据权利要求1所述的脑电信号去噪方法，其特征在于，对变分模态分解模型进行优化，具体过程为：将最小排列熵作为蜣螂优化算法的适应度函数，利用蜣螂优化算法，优化变分模态分解模型中的调节参数，得到最佳参数组合；最小排列熵的计算过程如下：

记脑电信号的时间序列经变分模态分解得到的其中一个模态分量为{X(i),i＝1,2,…,N}，N为正整数，对其进行相空间重构,得到重构矩阵Y如下：

其中，d表示嵌入维数，τ表示延时因子，M表示重构矩阵Y分量的个数；将重构矩阵Y中的第j个分量(x(j),x(j+τ),…,x(j+(d-1)τ)按照数值的大小进行升序排列，得到符号序列S(q)＝(j₁,j₂,…,j_d)；其中q＝1,2,…,r，r表示正整数，r≤d！，j₁,j₂,…,j_d表示各元素在原重构矩阵分量中的索引号；重构矩阵中的每个分量是d维空间，映射到d维符号序列S(q)，共有d！种排列方式；

计算每一种d维符号序列出现的概率{P1,P2,…,Pd}，基于计算出的概率，计算所述时间序列的排列熵H_p(d)如下：

其中，P_j表示j维符号序列出现的概率；排列熵的值越小，说明时间序列越简单、越规则，反之，排列熵的值越复杂，越不规则；将最小值的排列熵作为最小排列熵。

3.根据权利要求2所述的脑电信号去噪方法，其特征在于，利用蜣螂优化算法，优化变分模态分解模型中的调节参数，具体过程如下：

S21.初始化蜣螂种群大小和变分模态分解的调节参数范围；

S22.利用排列熵作为适应度函数，计算蜣螂种群的适应度值；

S23.通过模拟蜣螂行为，进行蜣螂的位置迭代更新；

S24.判断更新后蜣螂的位置是否位于原蜣螂的位置，若是，利用更新后蜣螂的位置代替原蜣螂的位置，执行步骤S25；否则，保持原蜣螂的位置不变，执行步骤S25；

S25.判断当前迭代次数是否小于最大迭代次数，若是，则返回步骤S22，否则，将当前蜣螂的位置作为最优位置，当前蜣螂种群的适应度值为最优适应度值，输出最优位置和最优适应度值对应的调节参数为最佳参数组合。

4.根据权利要求3所述的脑电信号去噪方法，其特征在于，所述调节参数为分解个数K和惩罚因子α。

5.根据权利要求4所述的脑电信号去噪方法，其特征在于，所述蜣螂行为包括滚球行为、跳舞行为、卵球繁殖行为、觅食行为和偷窃行为。

6.根据权利要求5所述的脑电信号去噪方法，其特征在于，所述通过模拟蜣螂行为，进行蜣螂的位置迭代更新，具体更新过程为：

利用滚球行为进行蜣螂的位置迭代更新，负责滚球的蜣螂的位置更新如下：

x_i(t+1)＝x_i(t)+α×k×x_i(t-1)+b×Δx

△x＝|x_i(t)-X^w|

其中，x_i表示第i个蜣螂的位置信息，t表示当前迭代次数，α∈(0,1)表示一个随机数，k代表一个常数，表示偏转系数，b表示自然系数，b取值为-1或1，X^w表示全局最差位置，△x表示光强的变化；

当遇到障碍物时，蜣螂利用跳舞行为重新定位方向，此时蜣螂的位置更新如下：

x_i(t+1)＝x_i(t)+tan(θ)|x_i(t)-x_i(t-1)|

tan(θ)表示正切函数，θ表示偏转系角度，当θ等于0或π或2π时，不更改蜣螂的位置。

当蜣螂执行卵球繁殖行为时，利用边界选择策略来模拟雌性蜣螂产卵的区域，边界选择策略定义为：

Lb^*＝max(X^*×(1-R),Lb)

Ub^*＝min(X^*×(1+R),Ub)

其中，Lb^*表示蜣螂产卵时区域的下界，Ub^*表示蜣螂产卵时区域的下界，表示X^*当前最佳位置，Lb表示优化问题的下界，Ub表示优化问题的上界，R表示惯性权值，其计算表达式为：

R＝1-t/T_max

其中，T_max表示最大迭代次数；

在迭代过程中，蜣螂产卵的卵球的位置随着产卵的区域动态变化：

B_i(t+1)＝X^*+b₁×(B_i(t)-Lb^*)+b₂×(B_i(t)-Ub^*)

其中，B_i表示第i个卵球的位置，b₁和b₂表示1×D的随机向量，D表示优化问题的维数；

当小蜣螂进行觅食行为时，建立最优觅食区域来引导小蜣螂进行觅食，其中最优觅食区域定义为：

Lb^b＝max(X^b×(1-R),Lb)

Ub^b＝min(X^b×(1+R),Ub)

其中，X^b表示全局最佳位置，Lb^b表示最佳觅食区域的下限，Ub^b表示最佳觅食区域的上限；小蜣螂的位置更新为：

x_i(t+1)＝x_i(t)+C₁×(x_i(t)-Lb^b)+C₂×(x_i(t)-Ub^b)

其中，C₁表示服从正态分布的随机数；C₂∈(0,1)表示一个随机向量；

当蜣螂具有偷盗行为时，具有偷盗行为的蜣螂的位置更新如下：

x_i(t+1)＝X^b+S×g×(|x_i(t)-X^*|+|x_i(t)-X^b|)

其中，X^b表示争夺食物的最佳地点，S表示一个常数，g表示服从均值为0、方差为1的正态分布的随机向量。

7.根据权利要求6所述的脑电信号去噪方法，其特征在于，基于优化后的变分模态分解模型，对脑电信号f进行变分模态分解，得到K个具有中心频率ω_k的有限带宽模态分量，具体过程包括：

S31.对脑电信号f进行希尔伯特变换，得到脑电信号f对应的解析信号，对解析信号进行单边谱处理，得到每个模态分量u_k(t)的单边谱：

其中，t表示时间，δ(t)表示关于时间t的狄拉克分布，j表示复数，k表示模态分量个数；

S32.将S31得到的单边谱与预估中心频率混合，调制出各模态信号带宽，各模态信号带宽受约束的变分表达式如下：

S33.基于步骤S32的变分表达式，构建目标约束函数：

其中，u_k表示本征模态分量，w_k表示各个模态分量的中心频率，s.t表示约束于，f为原始脑电信号，表示对模态分量个数k求偏导；

S34.引入二次惩罚因子α和拉格朗日乘子λ，构造增广拉格朗日函数L，将目标约束函数的约束变分问题转变成非约束变分问题，L的具体表达式如下：

其中，{u_k}表示本征模态分量u_k的集合，{w_k}表示各个模态分量的中心频率的集合，λ表示拉格朗日乘子算子，f(t)表示关于时间t的脑电信号，λ(t)表示关于时间t的拉格朗日乘子算子，<λ(t),f(t)-∑_ku_k(t)>表示拉格朗日乘子项；

利用交替方向乘子算法，得到更新后的各个模态分量、中心频率和拉格朗日乘子，更新后的模态分量的计算表达式为：

其中，表示更新后的本征模态分量，/>表示傅里叶变换的结果，/>表示本征模式分量在傅里叶域中的迭代值，/>表示更新后的拉格朗日乘子算子，w表示频率参数，k表示模态分量的个数；

更新后的中心频率的计算表达式为：

其中，为更新后的中心频率；

更新后的拉格朗日乘子的计算表达式为：

其中，表示更新后的拉格朗日乘子，/>表示拉格朗日乘子算子值在傅里叶域中的迭代值，τ表示更新参数。

S35.判断分解得到的模态分量是否满足下式：

其中，ε为判断阈值；若是，则循环分解模态分量完成；否则，返回步骤S34。

8.根据权利要求7所述的脑电信号去噪方法，其特征在于，所述计算K个模态分量中的模态分量与原始脑电信号之间的相关系数，具体包括：

S41.获取原始脑电信号的自相关函数R_ss和各模态分量的互相关函数R_J；

S42.基于R_ss和R_J，计算模态分量与原始脑电信号之间的相关系数：

其中，r_J表示模态分量与原始脑电信号之间的相关系数，J表示各模态分量的序号，r_J的值越大，模态分量与原始脑电信号之间的相关性越强，相反，r_J的值越小，模态分量与原始脑电信号之间的相关性越弱。

9.根据权利要求8所述的脑电信号去噪方法，其特征在于，从计算得出的模态分量与原始脑电信号之间的相关系数中获取数值最大的相关系数，基于数值最大的相关系数，计算临界值ρ：

其中，表示数值最大的相关系数；

将计算出的相关系数大于或等于ρ的模态分量记为有效模态分量，将计算出的相关系数小于ρ的模态分量记为噪声模态分量。

10.一种脑电信号去噪系统，其特征在于，所述系统包括：

信号获取模块，用于获取脑电信号，将脑电信号代入变分模态分解模型中；

优化模块，用于对变分模态分解模型进行优化，得到优化后的变分模态分解模型；

分解模块，用于根据优化后的变分模态分解模型，对脑电信号进行变分模态分解，得到K个模态分量；

计算模块，用于计算K个模态分量中的模态分量与原始脑电信号之间的相关系数，根据所述相关系数判断K个模态分量中的模态分量是否为噪声模态分量，若是，则去除噪声模态分量；否则，记该模态分量为有效模态分量，将有效模态分量进行重构，得到重构信号；

去噪模块，用于对重构信号进行小波阈值去噪，得到经过阈值处理的小波系数并小波重构，获得去噪干净的脑电信号。