CN103961091B - 基于双树复小波样本熵的运动想象脑电信号特征提取方法 - Google Patents

Abstract

本发明提出了一种基于双树复小波样本熵的运动想象脑电信号特征提取方法。首先，利用双树复小波变换将脑电信号分解在不同的频段，依据运动想象脑电信号中ERD/ERS现象，抽取有用信号的频段进行重构；然后，利用样本熵提取出脑电信号特定频段的非线性特征。本方法可以作为脑电数据信号分析中一种有效的特征处理方法，具有一定的可行性，可以获得比较高的识别率，它为BCI的特征提取提供了新的思路。

Description

基于双树复小波样本熵的运动想象脑电信号特征提取方法

技术领域

本发明属于脑电信号处理领域，涉及一种脑电信号特征提取方法，特别涉及一种用于脑-机接口中运动想象脑电信号识别的特征提取方法。

背景技术

脑-机接口(Brain Computer Interface，BCI)是一种不依赖于外周神经系统及肌肉组织的参与，人脑与计算机或其它外部设备之间建立起直接交流和控制的通道，是一种全新的通讯和控制方式。脑电信号是由脑神经细胞自发性、节律性的生理电活动产生的，是大脑意识过程中神经细胞电生理活动的外在反映，具有较高的时间分辨率，因而脑电图（EEG）成为BCI最重要的信号获取手段。

BCI系统中，通过脑电信号控制辅助外设的方法有多种：按照操作方式的自动化程度可分为单步、半自动、全自动等；按照脑电信号产生机理可分为诱发脑电和自发脑电两大类。其中，自发脑电由受试者自主产生，具有灵活性和可控制性，相对诱发脑电来说是更为自然和实用的方式。运动想象脑电信号是自发脑电信号的一种，在脑-机接口研究和应用中受到广泛的关注。

基于运动想象脑电信号的多模式识别是目前BCI主要应用手段之一。以人类大脑为对象的头皮脑电信号研究表明，它主要由各种节律性电活动组成，与运动想象紧密联系的一种电生理现象是事件相关去同步（ERD）和事件相关同步（ERS）现象。在进行基于运动想象脑电信号的脑-机接口研究中，将运动想象脑电信号分解在不同频段上进行特征提取具有一定的必要性。将脑电信号分解在不同频段上常见方法有小波变换，小波包变换等。这些方法尽管也取得了不错的效果，但是离散小波分解后的信号在相邻尺度的过渡频带上存在着能量泄漏现象，同时小波变换在二抽取的过程中会产生较大的混叠现象，这些缺陷在特征提取时会造成一定程度的假象，影响后续模式分类器的识别率。Kingsbury等人在1988年提出双树复小波变换(Dual-Tree Complex Wavelet Transform，DTCWT)，DTCWT是离散小波变换的衍生，可以有效的克服混叠和能量泄漏，而且还具有时移不变性、多维方向选择性、完全重构性等诸多优点，特征提取的效果明显好于小波分析。

脑电信号是一种随机的非线性信号，非线性特点比较明显，随着非线性理论的发展，很多非线性方法已经被广泛用于脑电信号的特征提取，例如，Pincus于二十世纪九十年代提出量化时间序列复杂度的近似熵算法。近似熵为相似向量由m维增加至m+1维时继续保持其相似性的条件概率，也是当维数变化的时候时间序列中产生新模式的概率大小，所以从统计的角度来区别时间过程的复杂性。但是近似熵中存在统计量的不一致性，针对这一不足，一种近似熵的改进方法—样本熵被Richman和Moorman提出。样本熵不仅具备近似熵的所有优点，而且避免了统计量的不一致性。样本熵是时间序列复杂度的一种度量，在实际应用过程中，与Lyapunov指数、信息熵、关联维数、K熵等非线性动力学方法相比，样本熵因为只需较短的数据就能够得出稳健的估计值，同时还有较好的抗噪和抗干扰能力，又可用于随机成分和确定性成分组成的混合信号中，分析效果优于简单统计参数，不需要对原始信号进行粗粒化等特点，比较适合对生物信号进行分析。发明内容

本发明提出一种双树复小波样本熵的特征提取方法。该方法通过双树复小波变换，把采集到的运动想象脑电信号进行分解，抽取出对应于运动想象脑电信号中ERD和ERS现象的节律波信号，进行重构，然后对该信号利用样本熵方法进行特征提取。实验表明，双树复小波样本熵的特征提取方法的具有一定的可行性，可以获得比较高的识别率。

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

步骤(1).抽取出对应频段的运动想象脑电信号。将采集到的运动想象脑电信号通过双数复小波变换进行分解，抽取出对应于运动想象脑电信号中ERD和ERS现象的节律波信号，再对其进行重构。

步骤(2).对信号进行特征提取。对含有有用成分的各层重构信号利用样本熵方法进行特征提取。

本发明与已有的运动想象脑电特征提取方法相比，具有如下特点：

1、利用双树复小波变换方法，在信号分解与重构过程中有实虚两平行的包含小波变换的双正交滤波器的小波树构成，能够实现交替奇偶滤波，有效地弥补复小波金字塔算法的重构性差的缺点，同时兼有计算效率高，数据冗余少的优点。

2、非线性动力学方法—样本熵，能够分析脑电信号的非线性特征。同时，运动想象脑电信号中的ERD和ERS现象在一些频段较为显著，更容易获取运动想象脑电信号的特征向量。

本发明可以有效表征运动想象时的EEG特征变换。因此，基于双树复小波样本熵的特征提取方法可以作为脑电数据信号分析中一种有效的特征处理方法，它为BCI的特征提取提供了新的思路，在脑-机接口领域具有广阔的应用前景。

附图说明

图1脑电信号特征提取流程图

图2双树复小波变换的分解和重构过程

具体实施方式

下面结合附图描述本发明基于双树复小波样本熵的运动想象脑电信号特征提取方法。

图1为脑电信号特征提取流程图，其实施主要包括以下几个步骤：

（1）利用双树复小波变换将采集到的运动想象脑电信号分解在不同的频段，依据运动想象脑电信号中ERD和ERS现象，抽取有用频段的信号进行重构；

（2）采用样本熵方法对重构后的不同频段的脑电信号提取非线性特征。

下面逐一对各步骤进行详细说明。

步骤一，利用双树复小波变换提取出对应频段的运动想象脑电信号。

（一）双树复小波变换的基本原理

DTCWT采用了二叉树结构的两路小波变换，将复小波的实部和虚部分离开，由两组并行的实数滤波器组获取实部和虚部的小波变换系数。DTCWT的主要思想是：在第一层分解时，要确保虚部树中的二抽取刚好采样得到实部树在二抽取过程中丢掉的信息，所以要在虚部树前面增加一个延迟器使得实部滤波器与虚部滤波器之间的延迟刚好是一个采样间隔，在之后的各层分解过程中，要求实虚两树对应滤波器的幅频响应相等，相频响应有半个采样周期的群延迟。同时，实虚滤波器采用双正交变换保证相位为线性的，两树滤波器长度分别是奇数长度和偶数长度，且每树不同层次之间交替采用奇偶滤波器，保证两树呈好的对称性。

由于双树复小波变换是基于两个并行的小波变换，因而根据小波分析的相关理论，在图2双树复小波变换的分解和重构过程中，↓2表示下取样算子，↑2表示上取样算子，实部树的小波系数及尺度系数分别如下：

$d_J^{\mathrm{Re}}\left(n\right)=2^{j/2}{\∫}_{-\∞}^{+\∞}x\left(t\right){\mathrm{\ψ}}_h(2^{j}t-n)\mathrm{dt}---\left(1\right)$

$d_J^{\mathrm{Re}}\left(n\right)=2^{J/2}{\∫}_{-\∞}^{+\∞}x\left(t\right){\mathrm{\φ}}_h(2^{J}t-n)\mathrm{dt}---\left(2\right)$

式(1)及(2)中，n为采样次数，j为比例因子，J为最大尺度且j＝1,2,…,J,x(t)为采样信号。

同理，虚部树的小波系数及尺度系数分别如下：

$d_j^{\mathrm{Im}}\left(n\right)=2^{j/2}{\∫}_{-\∞}^{+\∞}x\left(t\right){\mathrm{\ψ}}_g(2^{j}t-n)\mathrm{dt}---\left(3\right)$

$c_J^{\mathrm{Im}}\left(n\right)=2^{J/2}{\∫}_{-\∞}^{+\∞}x\left(t\right){\mathrm{\φ}}_g(2^{J}t-n)\mathrm{dt}---\left(4\right)$

式（1）、（2）、（3）及（4）中，ψ_h、φ_h、ψ_g、φ_g表示小波变换函数，满足下列关系： ${\mathrm{\φ}}_h\left(t\right)=\sqrt{2}\underset{n}{\mathrm{\Σ}}{h}_0\left(n\right)\mathrm{\φ}(2t-n),$ ${\mathrm{\ψ}}_h\left(t\right)=\sqrt{2}\underset{n}{\mathrm{\Σ}}{h}_1\left(n\right)\mathrm{\φ}(2t-n),$ ${\mathrm{\φ}}_g\left(t\right)=\sqrt{2}\underset{n}{\mathrm{\Σ}}{g}_0\left(n\right)\mathrm{\φ}(2t-n),$ ${\mathrm{\ψ}}_g\left(t\right)=\sqrt{2}\underset{n}{\mathrm{\Σ}}{g}_1\left(n\right)\mathrm{\φ}(2t-n).$ 其中，h₀(n)和h₁(n)分别表示实数部低通滤波器和高通滤波器，而g₀(n)和g₁(n)分别表示虚数部低通滤波器和高通滤波器。

综合实虚部相关系数，可得双树复小波变换分解过程的小波系数和尺度系数分别如下：

$d_j^{\left(C\right)}\left(n\right)=d_j^{\mathrm{Re}}\left(n\right)+{\mathrm{id}}_j^{\mathrm{Im}}\left(n\right)---\left(5\right)$

$c_J^{\left(C\right)}\left(n\right)=c_J^{\mathrm{Re}}\left(n\right)+{\mathrm{ic}}_J^{\mathrm{Im}}\left(n\right)---\left(6\right)$

（二）选取有用频段的信号，并对其进行重构

在利用双树复小波变换对采集到的脑电信号进行分解时，分解的层数将视具体信号的有用成分和采样频率而定。设采集脑电信号的采样频率为f_s，采用双树复小波变换对信号进行分解的层次为L，对低频子带复系数cA_L和高频子带复系数cD_L，cD_L-1，……，cD₁进行系数重构，根据分解原理可以得到L+1个重构信号，频段范围由低到高依次为[0,f_s/2^L+1]，[f_s/2^L+1,f_s/2^L]，[f_s/2^L,f_s/2^L-1]，……，[f_s/2²,f_s/2]。其中f(D_l)∈[f_s/2^l+1,f_s/2^l]Hz，l∈1,2,...,L,f(A_L)∈[0,f_s/2^L+1]。然后选取有用信号的alpha节律和beta节律波对应的频率范围作为重构信号。

双树复小波变换重构过程的小波系数d_j(t)和尺度系数c_J(t)如下：

$d_j\left(t\right)=2^{(j-1)/2}[\underset{n=-\∞}{\overset{\∞}{\mathrm{\Σ}}}{d}_j^{\mathrm{Re}}\left(n\right){\mathrm{\ψ}}_h(2^{j}t-n)+\underset{k=-\∞}{\overset{\∞}{\mathrm{\Σ}}}{d}_j^{\mathrm{Im}}\left(n\right){\mathrm{\ψ}}_g(2^{j}t-k)]---\left(7\right)$

$c_J\left(t\right)=2^{(J-1)/2}[\underset{n=-\∞}{\overset{\∞}{\mathrm{\Σ}}}{c}_J^{\mathrm{Re}}\left(n\right){\mathrm{\φ}}_h(2^{J}t-n)+\underset{k=-\∞}{\overset{\mathrm{Im}}{\mathrm{\Σ}}}{c}_J^{\mathrm{Im}}\left(n\right){\mathrm{\φ}}_g(2^{J}t-k)]---\left(8\right)$

由式(7)和(8)可得双树复小波变换多尺度分解下重构后的信号：

x(t)＝d_j(t)+c_J(t) (9)

步骤二，利用样本熵提取出重构后的脑电信号的非线性特征。特征选取是模式识别的核心问题，它不仅影响的到分类器的设计，还关系到分类的有效性。本发明运用样本熵方法对重构的不同频段的脑电信号提取非线性特征。

采用样本熵的快速算法，在保持样本熵优点的情况下，计算速度更快。

（一）快速样本熵的计算步骤

(1)设重构后的信号时间序列{x_i}含有N个数据，分别为x(1),x(2)…x(N)。

(2)将序列{x_i}按顺序组成一组m维矢量，

X(i)＝[x(i),x(i+1),…,x(i+m-1)] (10)式中，i=1～N-m+1；

矢量X(i)与X(j)之间的距离dis[X(i),X(j)]用分量间的最大欧氏距离表示，即

$\mathrm{dis}[X\left(i\right),X\left(j\right)]=\underset{k=0,...,m-1}{\max }|x_{i+k}-x_{j+k}|---\left(11\right)$

式中，k=1～m-1,i,j=1～N-m+1;

(3)定义N×N的二值距离矩阵为D，D的第i行第j列为d_ij，给定阈值r（r>0），则：

$d_{\mathrm{ij}}=\left\{\begin{array}{c}0,\mathrm{dis}[X\left(i\right),X\left(j\right)]\&GreaterEqual;r\\ 1,\mathrm{dis}[X\left(i\right),X\left(j\right)]<r\end{array},,(i,j=1~N)\right.---\left(12\right)$

(4)对每一个i值统计dis[X(i),X(j)]＜r的数目，记为同理将维数加1，计算

$B_i^m\left(r\right)=\underset{j=1}{\overset{N-1}{\mathrm{\&Sigma;}}}{d}_{\mathrm{ij}}\&cap;d_{(i+1)(j+1)}...d_{(i+m-1)(j+m-1)}---\left(13\right)$

$B_i^{m+1}\left(r\right)=\underset{j=1}{\overset{N-2}{\mathrm{\&Sigma;}}}{d}_{\mathrm{ij}}\&cap;d_{(i+1)(j+1)}...d_{(i+m)(j+m)}---\left(14\right)$

(5)求所有的的平均值记为B^m(r)及的平均值B^m+1(r)，

$B^m\left(r\right)=\frac{1}{N-m+1}\underset{i=1}{\overset{N-m+1}{\mathrm{\&Sigma;}}}{B}_i^m\left(r\right)---\left(15\right)$

$B^{m+1}\left(r\right)=\frac{1}{N-m}\underset{i=1}{\overset{N-m}{\mathrm{\&Sigma;}}}{B}_i^{m+1}\left(r\right)---\left(16\right)$

(6)脑电信号时间序列{x_i}的样本熵可表示为：

SampEn(N,m,r)＝-ln[B^m+1(r)/B^m(r)] (17)（二）样本熵参数的选择

样本熵SampEn(N,m,r)的值与嵌入维数m，相似容限r，数据长度N都有关系，而参数m、r的选择是样本熵估计的关键。大多数情况下嵌入维数取m=2，因为随着嵌入维数m列的联合概的增大，序列进行动态重构时，详细信息就会越多，于此同时嵌入维数m需要长度为N=₁₀m～₂₀m的数据长度，在计算的过程中所要的时间也会越来越长。相似容限一般r取0.1～0.25SD(SD为原始数据的标准差)，相似容限r如果选取的过大，很多时间序列的细节信息会被丢失掉，如果过小，噪声对统计结果的影响就很显著。鉴于此，为了抑制噪声对信号的干扰，就需要使信号中噪声的幅值小于相似容限r。

Claims (1)

1.基于双树复小波样本熵的运动想象脑电信号特征提取方法，其特征在于：

步骤(1).抽取出对应频段的运动想象脑电信号：将采集到的运动想象脑电信号通过双数复小波变换进行分解，抽取出对应于运动想象脑电信号中ERD和ERS现象的节律波信号，再对其进行重构，具体是：

设采集脑电信号的采样频率为f_s，采用双树复小波变换对信号进行分解的层次为L，对低频子带复系数cA_L和高频子带复系数cD_L，cD_L-1，……，cD₁进行系数重构，根据分解原理得到L+1个重构信号，频段范围由低到高依次为[0,f_s/2^L+1]，[f_s/2^L+1,f_s/2^L]，[f_s/2^L,f_s/2^{L -1}]，……，[f_s/2²,f_s/2]；其中f(D_l)∈[f_s/2^l+1,f_s/2^l]Hz，l∈1,2,...,L,f(A_L)∈[0,f_s/2^L+1]；然后选取有用信号的alpha节律和beta节律波对应的频率范围作为重构信号；

双树复小波变换重构过程的小波系数d_j(t)和尺度系数c_J(t)如下：

$d_j\left(t\right)=2^{(j-1)/2}[\underset{n=-\&infin;}{\overset{\&infin;}{\mathrm{\&Sigma;}}}{d}_j^{\mathrm{Re}}\left(n\right){\mathrm{\&psi;}}_h(2^{j}t-n)+\underset{k=-\&infin;}{\overset{\&infin;}{\mathrm{\&Sigma;}}}{d}_j^{\mathrm{Im}}\left(n\right){\mathrm{\&psi;}}_g(2^{j}t-k)]$

$c_J\left(t\right)=2^{(J-1)/2}[\underset{n=-\&infin;}{\overset{\&infin;}{\mathrm{\&Sigma;}}}{c}_J^{\mathrm{Re}}\left(n\right){\mathrm{\&phi;}}_h(2^{J}t-n)+\underset{k=-\&infin;}{\overset{\mathrm{Im}}{\mathrm{\&Sigma;}}}{c}_J^{\mathrm{Im}}\left(n\right){\mathrm{\&phi;}}_g(2^{J}t-k)]$

由以上两式得双树复小波变换多尺度分解下重构后的信号：

x(t)＝d_j(t)+c_J(t)

其中n为采样次数，j为比例因子，J为最大尺度且j＝1,2,…,J,ψ_h、φ_h、ψ_g、φ_g表示小波变换函数，分别为实部树的小波系数及尺度系数，分别为虚部树的小波系数及尺度系数；

步骤(2).对信号进行特征提取：对含有有用成分的各层重构信号利用样本熵方法进行特征提取，具体是：

（一）快速样本熵的计算：

(1)设重构后的信号时间序列{x_i}含有N个数据，分别为x(1),x(2)…x(N)；

(2)将序列{x_i}按顺序组成一组m维矢量，

X(i)＝[x(i),x(i+1),…,x(i+m-1)]

式中，i=1～N-m+1；

矢量X(i)与X(j)之间的距离dis[X(i),X(j)]用分量间的最大欧氏距离表示，即

$\mathrm{dis}[X\left(i\right),X\left(j\right)]=\underset{k=0,...,m-1}{\max }|x_{i+k}-x_{j+k}|$

式中，k=1～m-1,i,j=1～N-m+1;

(3)定义N×N的二值距离矩阵为D，D的第i行第j列为d_ij，给定阈值r，r>0，则：

$d_{\mathrm{ij}}=\left\{\begin{array}{c}0,\mathrm{dis}[X\left(i\right),X\left(j\right)]\&GreaterEqual;r\\ 1,\mathrm{dis}[X\left(i\right),X\left(j\right)]<r\end{array},,(i,j=1~N)\right.$

(4)对每一个i值统计dis[X(i),X(j)]＜r的数目，记为同理将维数加1，计算

$B_i^m\left(r\right)=\underset{j=1}{\overset{N-1}{\mathrm{\&Sigma;}}}{d}_{\mathrm{ij}}\&cap;d_{(i+1)(j+1)}...d_{(i+m-1)(j+m-1)}$

$B_i^{m+1}\left(r\right)=\underset{j=1}{\overset{N-2}{\mathrm{\&Sigma;}}}{d}_{\mathrm{ij}}\&cap;d_{(i+1)(j+1)}...d_{(i+m)(j+m)}$

(5)求所有的的平均值记为B^m(r)及的平均值B^m+1(r)，

$B^m\left(r\right)=\frac{1}{N-m+1}\underset{i=1}{\overset{N-m+1}{\mathrm{\&Sigma;}}}{B}_i^m\left(r\right)$

$B^{m+1}\left(r\right)=\frac{1}{N-m}\underset{i=1}{\overset{N-m}{\mathrm{\&Sigma;}}}{B}_i^{m+1}\left(r\right)$

(6)脑电信号时间序列{x_i}的样本熵可表示为：

SampEn(N,m,r)＝-ln[B^m+1(r)/B^m(r)]

（二）样本熵参数的选择

样本熵SampEn(N,m,r)的值与嵌入维数m，相似容限r，数据长度N都有关系，选定嵌入维数取m=2，于此同时嵌入维数m需要长度为N=₁₀m～₂₀m的数据长度，相似容限r取0.1～0.25SD，SD为原始数据的标准差。