CN103610461B - 基于双密度小波邻域相关阈值处理的脑电信号消噪方法 - Google Patents

Abstract

本发明涉及一种基于双密度小波邻域相关阈值处理的脑电信号消噪方法。当前大都采用经典离散小波变换结合传统阈值法对脑电信号进行消噪，而已有的基于经典小波变换结合传统阈值法的消噪方法都存在不足之处。本发明首先从大脑皮层上采集脑电信号，然后用双密度小波正变换对脑电信号进行分解，得到多层的信号高频系数，根据小波系数的局部统计依赖性，运用邻域相关阈值处理算法进行收缩，最后将收缩后的小波系数进行重构得到消噪后的信号。本发明根据脑电信号特性和干扰噪声特点，以信噪比为目标函数，采用网格寻优法对邻域相关阈值处理算法中的三个可调参数进行寻优，进而有效地平滑噪声并保留EEG的细节特征。

Description

基于双密度小波邻域相关阈值处理的脑电信号消噪方法

技术领域

本发明属于生物医学信号消噪领域，涉及一种基于双密度小波邻域相关阈值处理的脑电信号消噪方法。

背景技术

脑电信号(Electroencephalogram,EEG)是中枢神经系统产生的生物电活动而衍射到脑皮的募集信号，人在主动思维或受到不同的感觉刺激时，EEG的特征具有明显差异。对EEG进行分析，可以获得大量的生理、心理及病理信息。EEG是一种非线性非平稳信号,而且非常微弱,其幅值仅为μV级，极易受到噪声干扰，如脉冲干扰、工频干扰、呼吸干扰、眼动干扰、心电干扰、肌电干扰、头皮电极的抖动等，这些干扰产生的噪声混杂在EEG中，使得EEG不能正确反映大脑的真实生物电活动特征，给医生诊断带来误诊，也同样给研究人员对EEG的进一步深入研究，例如基于EEG的脑机接口(BCI)中运动意识任务的模式识别带来困难。因此，消除EEG中的噪声就显得非常必要。

小波变换因其具有良好的时频局部化特性、多分辨率特性、去相关性等特点，在生物医学信号消噪领域得到了广泛应用，也是传统的脑电信号消噪方法。经典离散小波变换虽功能强大，但它采用的是临界采样，导致经典离散小波变换具有数据敏感性，输入数据在时域上微小的变化会对离散小波变换系数产生不可预测的结果，即使输入有很小的移位都可能导致在各个尺度上的能量分布有明显的变化,这将降低BCI系统的EEG控制信号鲁棒性。小波系数阈值处理策略的选取是影响基于小波变换阈值消噪方法降噪效果的主要因素之一。信号经小波分解后，每个子空间内的小波系数具有一定的统计依赖性，利用这种局部依赖性，可以有效地提高信号的降噪效果。但传统软、硬阈值算法只是单独对每一个小波系数进行阈值处理，并未考虑小波系数的局部依赖性，造成系数估计值偏差过大，进而使得传统小波阈值消噪方法对EEG的消噪性能降低。

综上所述，已有的基于经典小波变换和传统软阈值法的脑电信号消噪方法存在不足之处。

发明内容

本发明的目的就是针对现有技术的不足，提出一种基于双密度小波邻域相关阈值处理的脑电信号消噪新方法。双密度小波变换（Double-DensityDiscreteWaveletTransform，DD-DWT）是一种新的小波分析方法，它的滤波器组是由一个尺度函数和两个彼此之间偏移0.5个单位的小波函数构成，同时使用过采样代替临界采样，因而它很好的客服了经典离散小波变换的局限性，具有运算速度快、近似的平移不变性、完美的重构性和有限的冗余性等特点，且与经典小波变换相比能够更准确地描述信号的细节特征，进而提高信号特征的鲁棒性；邻域相关阈值处理算法在小波系数的估计计算中融入了邻域内小波系数的信息，可以减小系数估计值的偏差，进而有效地平滑噪声并保留EEG的细节特征。

实现本发明方法的主要思路是：利用双密度小波对采集到的EEG分解，得到多层的信号高频系数；根据小波系数的局部统计依赖性，运用邻域相关阈值处理算法进行收缩，将收缩后的小波系数进行重构得到消噪后的信号。

为了实现以上目的，本发明方法主要包括以下步骤：

步骤一.获取人大脑运动想象脑电信号样本数据。

步骤二.将步骤一获取的运动想象脑电信号进行双密度小波分解，得到低频小波系数和各层高频小波系数，分别记为C_J(l)和w_j,h1(l)、w_j,h2(l)。其中J表示最大分解尺度；j表示分解尺度；h₁、h₂表示两个高通滤波器；l表示小波系数序列的标号，其中w_j,h1(l)、w_j,h2(l)分别表示子空间(j,h₁)和子空间(j,h₂)内的小波系数序列。

步骤三.计算各分解子空间的Donoho阈值，记为λj,hk：

${\mathrm{\λ}}_{j,h_k}=\mathrm{median}\left(\right|w_{j,h_k}\left|\right)*\sqrt{2\mathrm{InM}}/0.6745$

其中，M为子空间(j,h_k)内的小波系数数目。

步骤四.定义加权阈值缩放因子λ_j，计算式为：

λ_j＝β/(1+In(j))

其中，j为分解尺度。

步骤五.对Donoho阈值进行缩放得到各分解子空间的阈值，记为T_j,hk：

T_j,hk＝λ_j*λ_j,hk

根据分解尺度对频率相对高子空间的阈值进行放大、频率相对低子空间的阈值进行缩小，进而增强对EEG中高频噪声的抑制，同时保留低频有用信息。

步骤六.计算以w_j,hk(l),k＝1,2为中心，大小为2m+l的操作邻域窗口内小波系数的均值,记为：

${\stackrel{\‾}{w}}_{j,h_k}\left(l\right)=\frac{1}{2m+1}{\mathrm{\Σ}}_{m_1=-m}^m{w}_{j,h_k}(l+m_1)---\left(1\right)$

步骤七.使用邻域相关阈值处理函数对高频系数w_j,hk(l)进行收缩，其函数构造为:

式(2)在小波系数的估计计算中考虑了相邻小波系数的影响。当一个较大的小波系数周围领域的小波系数都相对较小，那么这个较大的小波系数受噪声污染严重的可能性极大,由式(1)、(2)知，邻域相关处理函数能够很好地平滑孤立的受噪声污染的小波系数，能对混杂在EEG中的脉冲干扰噪声、心电干扰噪声等起到有效的抑制作用。此外，当时，，从而很好地克服了软阈值中与w_j,hk(l)之间具有恒定偏差的缺点，有效地保留了信号的一些重要局部信息；当K→∞时，它与硬阈值函数相当，但是避免了硬阈值函数在阈值±T处不连续的缺点。因此，只要选取适合的参数K，就可以既保证重构信号的整体连续性、抑制振荡，又能提高重构信号的精度。

本发明以信噪比为目标函数，使用网格搜索法对m、β、K三个可调参数进行了寻优，当m＝2、β＝1.15、K＝4时，算法的消噪效果最为理想。参数寻优的具体步骤如下：

(1)根据真实EEG波形特点，采用Matlab生成采样率Fs(与脑电信号采样频率一致)、时间t=0:1/Fs:4、频率范围为2～30Hz的标准信号 $s\left(t\right)={\mathrm{\Σ}}_{m=2}^{13}{2}^{-\mathrm{rand}\left(1\right)}\sin (2\mathrm{\πnt}+2*\mathrm{rand}\left(1\right)*\mathrm{\π})+{\mathrm{\Σ}}_{n=14}^{30}\frac{\mathrm{rand}\left(1\right)}{5}\cos (2\mathrm{\πnt}+2*\mathrm{rand}\left(1\right)*\mathrm{\π}),$ 与EEG的主要节律（δ、θ、α、β节律）对应。

(2)对标准信号加噪声x(t)＝awgn(s(t),n,'measured')，即标准信号加噪后的信噪比为n（dB）。

(3)将x(t)进行双密度小波分解，得到各层高频小波系数，记为w_j,h1(l)、w_j,h2(l)；计算各分解子空间的Donoho阈值，记为λ_j,hk。

(4)设定m∈[0，(M-1)/2]，M为最大分解子空间内小波系数数目，步长Δm=1；设定β∈[0，20]，步长Δβ=0.05；设定K∈[1，60]，步长ΔK=0.05。然后，根据步骤四到步骤七计算每种参数组合(m，β，K)的消噪信号x以及信噪比SNR_（m,β,K）。信噪比定义为：

式中，N是信号长度，s(i)是原始标准信号，为去噪后的信号。通常信噪比越大信号的整体消噪效果越好，因此，只要找出SNR_（m,β,K）的最大值就可以寻得最优的(m，β，K)参数组合。

步骤八.将低频系数和经收缩后的高频系数进行双密度小波逆变换，得到消噪后的EEG。

本发明与已有的诸多脑电信号消噪方法相比，具有如下特点：

由于双密度小波变换运算速度快，且具有出色的信号特征表达能力，特别适合用于对实时性和稳健性要求高的在线BCI系统中。邻域相关阈值处理算法在小波系数的估计计算中融入了邻域内小波系数的信息，可以减小系数估计值的偏差，进而有效地平滑噪声并保留EEG的细节特征。

附图说明

图1为本发明的实施流程图；

图2为本发明的实施例的实验范式示意图；

图3为本发明实施例的原始含噪声的脑电信号图；

图4为本发明实施例的消噪后的脑电信号图；

图5为本发明实施例的原始含噪声的脑电信号频谱图；

图6为本发明实施例的消噪后的脑电信号频谱图；

具体实施方式

下面结合附图对本发明的实施例作详细说明：本实施例在以本发明技术方案为前提下进行实施，给出了详细的实施方式和具体的操作过程。

如图1所示，本实施例包括如下步骤：

步骤一.获取人大脑运动想象脑电信号样本数据，具体是：本发明采用美国Neuroscan公司的Scan4.3系统采集样本数据，采样频率为250Hz，精度为32bit。总共10名受试者，均为身体健康的24±1.6岁的在校大学生，在头脑清醒的情况下接受测试。脑电电极按照国际标准导联10-20系统放置，以左侧乳突为参考电极，右侧乳突为接地电极。数据采集实验示意图如图2所示。采集过程由Presentation软件控制完成,受试者坐在扶椅上双目注视屏幕,每一次实验持续9s。在0～4s,受试者保持休息状态；在第4s时，显示器上出现一个十字光标，持续1s,同时伴随一个提示声音提示受试者做好准备；在第5s时,十字光标由一个指示左、右、上、下方向的箭头或圆提示符代替,受试者根提示符想象左手、右手、左脚、右脚或舌头运动，直到第9s实验完成。总共采集50组EEG数据,即每名受试者每类运动意识任务各执行一次。本发明选取其中一段典型样本C3通道EEG作为实施例进行详细说明，图3是选取的EEG，原始脑电信号的频谱图如图4所示。。

步骤二.将步骤一获取的运动想象脑电信号进行双密度小波分解，得到低频小波系数和各层高频小波系数，分别记为C_J(l)和w_j,h1(l)、w_j,h2(l)。其中J表示最大分解尺度，本实施实例的J取2；j表示分解尺度；h₁、h₂表示两个高通滤波器；l表示小波系数序列的标号，其中w_j,h1(l)、w_j,h2(l)分别表示子空间(j,h₁)和子空间(j,h₂)内的小波系数序列。

步骤三.计算各分解子空间的Donoho阈值，记为λ_j,hk：

${\mathrm{\λ}}_{j,h_k}=\mathrm{median}\left(\right|w_{j,h_k}\left|\right)*\sqrt{2\mathrm{InM}}/0.6745$

其中，M为子空间(j,h_k)内的小波系数数目。

步骤四.定义加权阈值缩放因子λ_j，计算式为：

λ_j＝β/(1+In(j))

其中，j为分解尺度。

步骤五.对Donoho阈值进行缩放得到各分解子空间的阈值，记为T_j,hk：

T_j,hk＝λ_j*λ_j,hk

根据分解尺度对频率相对高子空间的阈值进行放大、频率相对低子空间的阈值进行缩小，进而增强对EEG中高频噪声的抑制，同时保留低频有用信息。

步骤六.计算以w_j,hk(l),k＝1,2为中心，大小为2m+l的操作邻域窗口内小波系数的均值,记为：

${\stackrel{\‾}{w}}_{j,h_k}\left(l\right)=\frac{1}{2m+1}{\mathrm{\Σ}}_{m_1=-m}^m{w}_{j,h_k}(l+m_1)$

步骤七.使用邻域相关阈值处理函数对高频系数w_j,hk(l)进行收缩，其函数构造为:

对高频系数收缩前需确定m、β、K三个可调参数，本发明以信噪比为目标函数，使用网格搜索法对m、β、K三个可调参数进行了寻优，当m＝2、β＝1.15、K＝4时，算法的消噪效果最为理想。参数寻优的具体步骤如下：

(1)根据真实EEG波形特点，采用Matlab生成采样率250Hz(与脑电信号采样频率一致)、时间t=0:1/250:4、频率范围为2～30Hz的标准信号 $s\left(t\right)={\mathrm{\Σ}}_{m=2}^{13}{2}^{-\mathrm{rand}\left(1\right)}\sin (2\mathrm{\πnt}+2*\mathrm{rand}\left(1\right)*\mathrm{\π})+{\mathrm{\Σ}}_{n=14}^{30}\frac{\mathrm{rand}\left(1\right)}{5}\cos (2\mathrm{\πnt}+2*\mathrm{rand}\left(1\right)*\mathrm{\π}),$ 与EEG的主要节律（δ、θ、α、β节律）对应。

(2)对标准信号加噪声x(t)＝awgn(s(t),n,'measured')，即标准信号加噪后的信噪比为n（dB）。

(3)将x(t)进行2尺度双密度小波分解，得到各层高频小波系数，记为w_j,h1(l)、w_j,h2(l)；计算各分解子空间的Donoho阈值，记为λ_j,hk。

(4)设定m∈[0，(M-1)/2]，M为最大分解子空间内小波系数数目，步长Δm=1；设定β∈[0，20]，步长Δβ=0.05；设定K∈[1，60]，步长ΔK=0.05。然后，根据步骤四到步骤八计算每种参数组合(m，β，K)的消噪信号以及信噪比SNR_（mβK）。信噪比定义为：

式中，N是信号长度，s(i)是原始标准信号，为去噪后的信号。通常信噪比越大信号的整体消噪效果越好，因此，只要找出SNR_（m,β,K）的最大值就可以寻得最优的(m，β，K)参数组合。

步骤八.将低频系数和经收缩后的高频系数进行双密度小波逆变换，得到消噪后的EEG，如图5所示，消噪后的脑电信号的频谱图如图6所示。

Claims (1)

1.基于双密度小波邻域相关阈值处理的脑电信号消噪方法，该方法首先获取大脑运动想象脑电信号样本数据，其特征在于该方法对运动想象脑电信号样本数据进行如下步骤处理：

步骤一、将获取的运动想象脑电信号进行双密度小波分解，得到低频小波系数和各层高频小波系数，分别记为C_J(l)和其中J表示最大分解尺度；j表示分解尺度；h₁、h₂表示两个高通滤波器；l表示小波系数序列的标号，其中分别表示子空间(j,h₁)和子空间(j,h₂)内的小波系数序列；

步骤二、计算各分解子空间的Donoho阈值，记为

${\mathrm{\λ}}_{j,h_k}=median\left(\right|w_{j,h_k}\left|\right)*\sqrt{2InM}/0.6745$

其中，M为子空间(j,h_k)内的小波系数数目；

步骤三、定义加权阈值缩放因子λ_j，计算式为：

λ_j＝β/(1+In(j))

其中，β≥0、j为分解尺度；

步骤四、对Donoho阈值进行缩放得到各分解子空间的阈值，记为

$T_{j,h_k}={\mathrm{\λ}}_j*{\mathrm{\λ}}_{j,h_k}$

根据分解尺度对频率相对高子空间的阈值进行放大、频率相对低子空间的阈值进行缩小，进而增强对脑电信号中高频噪声的抑制，同时保留低频有用信息；

步骤五、计算以k＝1,2为中心，大小为2m+l的操作邻域窗口内小波系数的均值,记为

${\stackrel{\‾}{w}}_{j,h_k}\left(l\right)=\frac{1}{2m+1}{\mathrm{\Σ}}_{m_1=-m}^m{w}_{j,h_k}(l+m_1)$

步骤六、使用邻域相关阈值处理函数对高频系数进行收缩，其函数构造为:

以信噪比为目标函数，使用网格搜索法对m、β、K三个可调参数进行寻优，参数寻优的具体步骤如下：

(1)根据真实脑电信号波形特点，生成采样率Fs、时间t＝0:1/Fs:4、频率范围为2～30Hz的标准信号

$s\left(t\right)={\mathrm{\Σ}}_{m=2}^{13}{2}^{-rand\left(1\right)}\sin \left(\begin{array}{cc}2\mathrm{\π}nt+2*& rand\left(1\right)*\mathrm{\π}\end{array}\right)+{\mathrm{\Σ}}_{n=14}^{30}\frac{rand\left(1\right)}{5}\cos \left(2\mathrm{\π}nt+2*rand\left(1\right)*\mathrm{\π}\right),$ 与脑电信号的主要节律对应；

(2)对标准信号加噪声x(t)＝awgn(s(t),n,'measured')，即标准信号加噪后的信噪比为n(dB)；

(3)将x(t)进行双密度小波分解，得到各层高频小波系数，记为计算各分解子空间的Donoho阈值，记为

(4)设定m∈[0，(M-1)/2]，M为最大分解子空间内小波系数数目，步长Δm＝1；设定β∈[0，20]，步长Δβ＝0.05；设定K∈[1，60]，步长ΔK＝0.05；然后，根据步骤三到步骤六计算每种参数组合(m，β，K)的消噪信号以及信噪比SNR_(m,β,K)；信噪比定义为：

式中，N是信号长度，s(i)是原始标准信号，为去噪后的信号；通常信噪比越大信号的整体消噪效果越好，因此，只要找出SNR_(m,β,K)的最大值就可以寻得最优的(m，β，K)参数组合；

经计算，当m＝2、β＝1.15、K＝4时，消噪效果最好，因此确定m＝2、β＝1.15、K＝4作为最优的(m，β，K)参数组合；

步骤七、将低频系数和经收缩后的高频系数进行双密度小波逆变换，得到消噪后的EEG。