CN104473635A - 混合小波和共同空间模式的左右手运动想象脑电信号特征提取方法 - Google Patents

Abstract

本发明涉及一种混合小波和共同空间模式的左右手运动想象脑电信号特征提取方法，属于脑电信号识别控制技术领域。该方法混合小波和共同空间模式算法对左右手运动想象脑电信号进行处理，提取出能够表征左右手运动想象脑电信号的时频空特征，具体包括以下步骤：1)使用小波变换对左右手运动想象脑电信号进行分解，得到能反映其时频特征的小波系数；2)使用CSP算法对能反映左右手运动想象脑电信号时频特征的小波系数D2进行分析，得到该小波系数的最优空间滤波器；3)使用该最优空间滤波器提取出能反映左右手运动想象脑电信号时频空特征的特征值。本方法能够解决现有技术中的脑电信号特征提取方法不能同时提取出脑电信号的时频特征和空域特征的问题。

Description

混合小波和共同空间模式的左右手运动想象脑电信号特征提取方法

技术领域

本发明属于脑电信号识别控制技术领域，涉及一种混合小波和共同空间模式的左右手运动想象脑电信号特征提取方法。

背景技术

脑-机接口(brain-computer interface,BCI)是不依赖于大脑外周神经与肌肉系统，在人脑和计算机或外部设备之间建立起来的一种通信系统。它能够为肢体残疾但思维意识正常的患者提供一种新型的对外信息交流手段，并在残疾人康复、正常人辅助控制、娱乐等领域有着广泛的应用前景。

特征提取是BCI研究的关键技术之一，常用的方法有FFT、AR、AAR、小波变换、共同空间模式(Common Spatial Pattern,CSP)等方法。FFT、AR和AAR通过把幅度随时间变化的脑电信号变换为脑电功率随频率变化的谱图，从而提取出脑电信号的频域特征。然而这些方法都只适合分析平稳信号，对于非平稳的脑电信号，具有很大的局限性。小波变换是一种典型的时频分析法，具有多分辨率特性，在时域和频域均具有良好的分辨率，适用于非平稳信号处理。因此，小波变换非常合适分析脑电信号。CSP算法是一种利用两类脑电信号的估计协方差矩阵来设计每类脑电信号的最优空间滤波器，进而提取脑电信号的空域特征的方法。

在上述脑电信号的特征提取方法中，虽然小波变换非常适合分析脑电信号，并能提取出脑电信号的时频特性，但却无法提取出脑电信号的空域特征，而CSP算法虽然能提取出脑电信号的空域特征，但又不能提取出脑电信号的时频特征。因此，在脑电信号的特征提取中，使用小波变换和CSP算法都不能充分地提取出脑电信号的特征。

发明内容

有鉴于此，针对现有的脑电信号特征提取方法不能同时提取出脑电信号的时频特征和空域特征的问题，本发明提供了一种混合小波和共同空间模式的左右手运动想象脑电信号特征提取方法。

为达到上述目的，本发明提供如下技术方案：

一种混合小波和共同空间模式的左右手运动想象脑电信号特征提取方法，该方法混合小波和共同空间模式算法对左右手运动想象脑电信号进行处理，提取出能够表征左右手运动想象脑电信号的时频空特征；具体包括以下步骤：

步骤一：使用小波变换对左右手运动想象脑电信号进行分解，得到能反映其时频特征的小波系数；

步骤二：使用CSP算法对能反映左右手运动想象脑电信号时频特征的小波系数D₂进行分析，得到该小波系数的最优空间滤波器；

步骤三：使用该最优空间滤波器提取出能反映左右手运动想象脑电信号时频空特征的特征值。

进一步，在步骤一中，小波变换对脑电信号的分解层数将视脑电信号的采样频率而定，根据脑电信号的采样频率，使用小波变换提取出左右手运动想象脑电信号的β节律所对应的小波系数。

进一步，所述共同空间模式算法是对左右手运动想象脑电信号的β节律所对应的小波系数进行分析，提取出的特征值能反映出左右手运动想象脑电信号的时频空特征。

本发明的有益效果在于：本发明提供了一种混合小波和共同空间模式的左右手运动想象脑电信号特征提取方法，该方法能够解决现有技术中的脑电信号特征提取方法不能同时提取出脑电信号的时频特征和空域特征的问题。

附图说明

为了使本发明的目的、技术方案和有益效果更加清楚，本发明提供如下附图进行说明：

图1为左右手运动想象脑电信号的特征提取方法框图；

图2为Emotiv脑电信号采集仪；

图3为Emotiv脑电采集仪的电极安放位置；

图4为小波变换对左右手运动想象脑电信号的分解框图。

具体实施方式

下面将结合附图，对本发明的优选实施例进行详细的描述。

在本方法中，首先使用小波变换对左右手运动想象脑电信号进行分解，得到能反映其时频特征的小波系数；然后使用CSP算法对能反映左右手运动想象脑电信号时频特征的小波系数进行分析，进而提取出能反映左右手运动想象脑电信号时频空特征的特征值。

具体步骤如下：

步骤一：对左右手运动想象脑电信号进行小波分解，得到能反映其时频特征的小波系数。

人在想象单侧手运动时，其对侧相应初级感觉运动皮层区的脑电μ节律(8～12Hz)和β节律(14～30Hz)节律幅值降低，这种现象称为事件相关去同步(event-related desynchronization,ERD)；而同侧脑电μ节律和β节律幅度升高，称为事件相关同步(event—relatedsynchronization，ERS)。根据这一特征，可使用μ节律和β节律来分析左右手运动想象脑电信号。而小波变换能把信号的整个频带划分为多个子频带，因此可使用小波变换来分析左右运动想象脑电信号。

设x(n)表示实验采集的脑电信号，则x(n)的离散小波变换定义为

$C_{j,k}=2^{-j/2}\underset{n=-\&infin;}{\overset{\&infin;}{\mathrm{\&Sigma;}}}x\left(n\right){\stackrel{\&OverBar;}{\mathrm{\&psi;}}}_{j,k}(2^{-j}n-k)=<x\left(n\right),{\mathrm{\&psi;}}_{j,k}>,j,k\&Element;Z,---\left(1\right)$

式中，ψ(n)为小波基函数，j、k分别代表频率分辨率和时间平移量。采用Mallat算法，对信号进行有限层分解，即

$x\left(n\right)=A_L+\underset{j=1}{\overset{L}{\mathrm{\&Sigma;}}}{D}_j.---\left(2\right)$

式中，L为分解层数，A_L为低通逼近分量，D_j为不同尺度下的细节分量。设信号x(n)的采样频率为f_s，则(1)式中的A_L、D_L、D_L-1、…、D₁各分量所对应的子频带依次为

$[0,\frac{f_s}{2^{L+1}}],[\frac{f_s}{2^{L+1}},\frac{f_s}{2^L}],[\frac{f_s}{2^L},\frac{f_s}{2^{L-1}}],\&CenterDot;\&CenterDot;\&CenterDot;,[\frac{f_s}{2^2},\frac{f_s}{2}].---\left(3\right)$

将信号进行小波分解时，分解的层数将视具体信号的有用成分和采样频率而定。若左右手运动想象脑电信号的采样频率为128Hz，则可对脑电信号进行2层分解，即x(n)＝A₂+D₂+D₁。

根据(3)式知D₂(16～32Hz)在脑电信号的β节律频带范围附近，而小波系数能表达信号在时域和频域的能量分布，因此，小波系数D₂能反映出左右手运动想象脑电信号的时频域特征。

步骤二：使用CSP算法对能反映左右手运动想象脑电信号时频特征的小波系数D₂进行分析，得到该小波系数的最优空间滤波器。

CSP算法统一模型的目标函数为：

$J\left(w_1\right)=\frac{w_1^T{X}_1^T{X}_1{w}_1}{w_1^T{X}_2^T{X}_2{w}_1}=\frac{w_1^T{C}_1{w}_1}{w_1^T{C}_2{w}_1},---\left(4\right)$

其中，X₁,X₂∈R^n×N表示两类脑电信号(如左手和右手运动想象脑电信号)，n是脑电信号采集时每个通道的采样点数，N是脑电信号的通道数，C₁,C₂∈R^N×N为这两类脑电信号经零均值化处理后的估计协方差矩阵(这里的估计协方差是所有训练样本的协方差矩阵的均值)，T表示矩阵的转置，w₁∈R^N×1为脑电信号X₁的空间滤波器。显然，使得J(w₁)最大的空间滤波器w₁，就是使得X₁和X₂两类脑电信号的方差比最大的最优空间滤波器。

观察式(4)，知

$J\left(w_1\right)=\frac{w_1^T{C}_1{w}_1}{w_1^T{C}_2{w}_1}=\frac{k^2{w}_1^T{C}_1{w}_1}{k^2{w}_1^T{C}_2{w}_1}=J\left({\mathrm{kw}}_1\right),---\left(5\right)$

其中k为任意非零实数。显然，存在一个k，可使得从而求使得J(w₁)最大的空间滤波器w₁，就可等价于在的条件下，求使得

$J\left(w_1\right)=k^2{w}_1^T{C}_1{w}_1---\left(6\right)$

最大的w₁。

使用拉格朗日乘数法求J(w₁)的最大值。设拉格朗日函数为

$L({\mathrm{\&lambda;}}_1,w_1)=k^2{w}_1^T{C}_1{w}_1-{\mathrm{\&lambda;}}_1(k^2{w}_1^T{C}_2{w}_1-1),---\left(7\right)$

其中λ₁为拉格朗日乘子。

对式(7)求关于w₁和λ₁的一阶偏导数，并使之为零，得到

$\frac{\&PartialD;L}{{\&PartialD;w}_1}={2k}^2{w}_1^T{C}_1-{2\mathrm{\&lambda;}}_1{k}^2{w}_1^T{C}_2=0,---\left(8\right)$

$\frac{\&PartialD;L}{{\&PartialD;\mathrm{\&lambda;}}_1}=1-k^2{w}_1^T{C}_2{w}_1=0.---\left(9\right)$

由式(8)得

$C_2^{-1}{C}_1{w}_1={\mathrm{\&lambda;}}_1{w}_1.---\left(10\right)$

从式(10)知，λ₁和w₁是矩阵的特征值和特征向量，因此矩阵的所有特征值可记为λ₁₁,λ₁₂,…,λ_1N(λ₁₁≥λ₁₂≥…≥λ_1N)，其对应的特征向量分别为w₁₁,w₁₂,…,w_1N。

由式(8)和式(9)得

$J\left(w_1\right)=\frac{k^2{w}_1^T{C}_1{w}_1}{k^2{w}_1^T{C}_2{w}_1}=k^2{w}_1^T{C}_1{w}_1={\mathrm{\&lambda;}}_1{k}^2{w}_1^T{C}_2{w}_1={\mathrm{\&lambda;}}_1.---\left(11\right)$

由式(11)知，J(w₁)的值等于λ₁，即等于矩阵的特征值。因此脑电信号X₁的空间滤波器就为矩阵的特征值对应的特征向量，而最优空间滤波器为矩阵的最大特征值λ₁₁对应的特征向量w₁₁。

同理，若目标函数为

$J\left(w_2\right)=\frac{w_2^T{X}_2^T{X}_2{w}_2}{w_2^T{X}_1^T{X}_1{w}_2}=\frac{w_2^T{C}_2{w}_2}{w_2^T{C}_1{w}_2},---\left(12\right)$

则，可求得J(w₂)的值等于矩阵的特征值λ₂，脑电信号X₂的空间滤波器就为矩阵的特征值对应的特征向量w₂₁,w₂₂,…,w_2N，而最优空间滤波器为矩阵的最大特征值对应的特征向量w₂₁。

步骤三：使用该最优空间滤波器提取出能反映左右手运动想象脑电信号时频空特征的特征值。

设X_L,X_R∈R^n×N分别表示左右手运动想象脑电信号β节律所对应的小波系数D₂，C_L,C_R∈R^N×N分别为X_L和X_R经零均值化处理后的估计协方差矩阵。由步骤2)知，左手运动想象脑电信号的空间滤波器为矩阵C_R ^-1C_L的特征向量；右手运动想象脑电信号的空间滤波器为矩阵C_L ^-1C_R的特征向量。

设矩阵C_R ^-1C_L和C_L ^-1C_R的特征向量排列方式是按对应特征值的降序方式排列，则可取矩阵C_R ^-1C_L和C_L ^-1C_R的第一列特征向量就是左右手运动想象脑电信号的最优的空间滤波器，分别记为w_L,w_R∈R^N×1。记一次实验任务的脑电信号β节律所对应的小波系数D₂为X∈R^n×N，则X通过两个滤波器后就变为Z_L＝Xw_L和Z_R＝Xw_R。

显然，当X为X_L时，有

(Z_L)^TZ_L＞(Z_R)^TZ_R，        (13)

当X为X_R时，有

(Z_L)^TZ_L＜(Z_R)^TZ_R.         (14)

根据以上分析，可得到左右手运动想象脑电信号的时频空特征为

$f=\lg \frac{{\left(Z_L\right)}^T{Z}_L}{{\left(Z_R\right)}^T{Z}_R}.---\left(15\right).$

图1为左右手运动想象脑电信号的特征提取框图。在该框图中，首先使用小波变换对左右手运动想象脑电信号进行分解，得到能反映其时频特征的小波系数；然后使用CSP算法对能反映左右手运动想象脑电信号时频特征的小波系数进行分析，进而提取出能反映左右手运动想象脑电信号时频空特征的特征值。

①小波变换

设x(n)表示实验采集的脑电信号，则x(n)的离散小波变换定义为

$C_{j,k}=2^{-j/2}\underset{n=-\&infin;}{\overset{\&infin;}{\mathrm{\&Sigma;}}}x\left(n\right){\stackrel{\&OverBar;}{\mathrm{\&psi;}}}_{j,k}(2^{-j}n-k)=<x\left(n\right),{\mathrm{\&psi;}}_{j,k}>,j,k\&Element;Z,---\left(16\right)$

式中，ψ(n)为小波基函数，j、k分别代表频率分辨率和时间平移量。采用Mallat算法，对信号进行有限层分解，即

$x\left(n\right)=A_L+\underset{j=1}{\overset{L}{\mathrm{\&Sigma;}}}{D}_j.---\left(17\right)$

式中，L为分解层数，A_L为低通逼近分量，D_j为不同尺度下的细节分量。设信号x(n)的采样频率为f_s，则(16)式中的A_L、D_L、D_L-1、…、D₁各分量所对应的子频带依次为

$[0,\frac{f_s}{2^{L+1}}],[\frac{f_s}{2^{L+1}},\frac{f_s}{2^L}],[\frac{f_s}{2^L},\frac{f_s}{2^{L-1}}],\&CenterDot;\&CenterDot;\&CenterDot;,[\frac{f_s}{2^2},\frac{f_s}{2}].---\left(18\right)$

将信号进行小波分解时，分解的层数将视具体信号的有用成分和采样频率而定。图2中的Emotiv脑电信号采集仪的采样频率为128Hz，图3为Emotiv脑电采集仪的电极安放位置，因此可对脑电信号进行2层分解，即x(n)＝A₂+D₂+D₁，其分解框图如图4所示。

根据图4知D₂(16～32Hz)在脑电信号的β节律频带范围附近，而小波系数能表达信号在时域和频域的能量分布，因此，小波系数D₂能反映出左右手运动想象脑电信号的时频域特征。

②CSP算法

CSP算法统一模型的目标函数为：

$J\left(w_1\right)=\frac{w_1^T{X}_1^T{X}_1{w}_1}{w_1^T{X}_2^T{X}_2{w}_1}=\frac{w_1^T{C}_1{w}_1}{w_1^T{C}_2{w}_1},---\left(19\right)$

其中，X₁,X₂∈R^n×N表示两类脑电信号(如左手和右手运动想象脑电信号)，n是脑电信号采集时每个通道的采样点数，N是脑电信号的通道数，C₁,C₂∈R^N×N为这两类脑电信号经零均值化处理后的估计协方差矩阵(这里的估计协方差是所有训练样本的协方差矩阵的均值)，T表示矩阵的转置，w₁∈R^N×1为脑电信号X₁的空间滤波器。显然，使得J(w₁)最大的空间滤波器w₁，就是使得X₁和X₂两类脑电信号的方差比最大的最优空间滤波器。

观察式(19)，知

$J\left(w_1\right)=\frac{w_1^T{C}_1{w}_1}{w_1^T{C}_2{w}_1}=\frac{k^2{w}_1^T{C}_1{w}_1}{k^2{w}_1^T{C}_2{w}_1}=J\left({\mathrm{kw}}_1\right),---\left(20\right)$

其中k为任意非零实数。显然，存在一个k，可使得从而求使得J(w₁)最大的空间滤波器w₁，就可等价于在的条件下，求使得

$J\left(w_1\right)=k^2{w}_1^T{C}_1{w}_1---\left(21\right)$

最大的w₁。

使用拉格朗日乘数法求J(w₁)的最大值。设拉格朗日函数为

$L({\mathrm{\&lambda;}}_1,w_1)=k^2{w}_1^T{C}_1{w}_1-{\mathrm{\&lambda;}}_1(k^2{w}_1^T{C}_2{w}_1-1),---\left(22\right)$

其中λ₁为拉格朗日乘子。

对式(22)求关于w₁和λ₁的一阶偏导数，并使之为零，得到

$\frac{\&PartialD;L}{{\&PartialD;w}_1}={2k}^2{w}_1^T{C}_1-{2\mathrm{\&lambda;}}_1{k}^2{w}_1^T{C}_2=0,---\left(23\right)$

$\frac{\&PartialD;L}{{\&PartialD;\mathrm{\&lambda;}}_1}=1-k^2{w}_1^T{C}_2{w}_1=0.---\left(24\right)$

由式(23)得

$C_2^{-1}{C}_1{w}_1={\mathrm{\&lambda;}}_1{w}_1.---\left(25\right)$

从式(25)知，λ₁和w₁是矩阵的特征值和特征向量，因此矩阵的所有特征值可记为λ₁₁,λ₁₂,…,λ_1N(λ₁₁≥λ₁₂≥…≥λ_1N)，其对应的特征向量分别为w₁₁,w₁₂,…,w_1N。

由式(23)和式(24)得

$J\left(w_1\right)=\frac{k^2{w}_1^T{C}_1{w}_1}{k^2{w}_1^T{C}_2{w}_1}=k^2{w}_1^T{C}_1{w}_1={\mathrm{\&lambda;}}_1{k}^2{w}_1^T{C}_2{w}_1={\mathrm{\&lambda;}}_1.---\left(26\right)$

由式(26)知，J(w₁)的值等于λ₁，即等于矩阵的特征值。因此脑电信号X₁的空间滤波器就为矩阵的特征值对应的特征向量，而最优空间滤波器为矩阵的最大特征值λ₁₁对应的特征向量w₁₁。

同理，若目标函数为

$J\left(w_2\right)=\frac{w_2^T{X}_2^T{X}_2{w}_2}{w_2^T{X}_1^T{X}_1{w}_2}=\frac{w_2^T{C}_2{w}_2}{w_2^T{C}_1{w}_2},---\left(27\right)$

则，可求得J(w₂)的值等于矩阵的特征值λ₂，脑电信号X₂的空间滤波器就为矩阵的特征值对应的特征向量w₂₁,w₂₂,…,w_2N，而最优空间滤波器为矩阵的最大特征值对应的特征向量w₂₁。

设X_L,X_R∈R^n×N分别表示左右手运动想象脑电信号β节律所对应的小波系数D₂，C_L,C_R∈R^N×N分别为X_L和X_R经零均值化处理后的估计协方差矩阵。由步骤二知，左手运动想象脑电信号的空间滤波器为矩阵C_R ^-1C_L的特征向量；右手运动想象脑电信号的空间滤波器为矩阵C_L ^-1C_R的特征向量。

设矩阵C_R ^-1C_L和C_L ^-1C_R的特征向量排列方式是按对应特征值的降序方式排列，则可取矩阵C_R ^-1C_L和C_L ^-1C_R的第一列特征向量就是左右手运动想象脑电信号的最优的空间滤波器，分别记为w_L,w_R∈R^N×1。记一次实验任务的脑电信号β节律所对应的小波系数D₂为X∈R^n×N，则X通过两个滤波器后就变为Z_L＝Xw_L和Z_R＝Xw_R。

显然，当X为X_L时，有

(Z_L)^TZ_L＞(Z_R)^TZ_R，          (28)

当X为X_R时，有

(Z_L)^TZ_L＜(Z_R)^TZ_R.          (29)

根据以上分析，可得到左右手运动想象脑电信号的时频空特征为

$f=\lg \frac{{\left(Z_L\right)}^T{Z}_L}{{\left(Z_R\right)}^T{Z}_R}.---\left(30\right).$

最后说明的是，以上优选实施例仅用以说明本发明的技术方案而非限制，尽管通过上述优选实施例已经对本发明进行了详细的描述，但本领域技术人员应当理解，可以在形式上和细节上对其作出各种各样的改变，而不偏离本发明权利要求书所限定的范围。

Claims (3)

1.一种混合小波和共同空间模式的左右手运动想象脑电信号特征提取方法，其特征在于：混合小波和共同空间模式算法对左右手运动想象脑电信号进行处理，提取出能够表征左右手运动想象脑电信号的时频空特征；具体包括以下步骤：

步骤一：使用小波变换对左右手运动想象脑电信号进行分解，得到能反映其时频特征的小波系数；

步骤二：使用CSP算法对能反映左右手运动想象脑电信号时频特征的小波系数D2进行分析，得到该小波系数的最优空间滤波器；

步骤三：使用该最优空间滤波器提取出能反映左右手运动想象脑电信号时频空特征的特征值。

2.根据权利要求1所述的一种混合小波和共同空间模式的左右手运动想象脑电信号特征提取方法，其特征在于：在步骤一中，小波变换对脑电信号的分解层数将视脑电信号的采样频率而定，根据脑电信号的采样频率，使用小波变换提取出左右手运动想象脑电信号的β节律所对应的小波系数。

3.根据权利要求1所述的一种混合小波和共同空间模式的左右手运动想象脑电信号特征提取方法，其特征在于：所述共同空间模式算法是对左右手运动想象脑电信号的β节律所对应的小波系数进行分析，提取出的特征值能反映出左右手运动想象脑电信号的时频空特征。