CN102824172B - 一种脑电特征提取方法 - Google Patents

Abstract

本发明公开了一种脑电特征提取方法，包括以下步骤：去除背景脑电数据和待处理脑电数据中的伪迹，分别获得背景脑电数据和待处理脑电数据的有效频段，分别将背景脑电数据和待处理脑电数据的有效频段分为若干数据段；对每个数据段分别提取时频特征和形态特征，得到每个数据段相应的时频特征值和形态特征值；利用背景脑电数据每个数据段的时频特征值和形态特征值，计算得到特征值的频率分布函数；利用背景脑电数据的频率分布函数，得到待处理脑电数据每个数据段的时频特征值出现的概率和形态特征值出现的概率；利用特征值的概率计算得到IMF-VoE特征值。利用IMF-VoE特征值可以快速有效识别脑电信号的变化特征，监测大脑状态。

Description

一种脑电特征提取方法

技术领域

本发明涉及脑电数据分析领域，具体涉及一种脑电特征提取方法。

背景技术

脑电数据信号是携带有人脑特征或状态的信息载体，人脑是一个开放的、时变的和非线性的系统，其产生的信号也是时变的，非线性的，同时脑电数据信号在测量后会产生随机误差，并且脑电信号还会受到个体差异的硬性，因此，对于脑电数据信号的分析成为难题。

申请号为200910196746.3的发明公开了一种脑电波分析方法，应用了经典的时频域分析和主成份分析方法解决了脑电信号特征提取的难题，成功地提取了与人体紧张，疲劳与放松息息相关的时频域参数，并且将其映射到主成分空间中去，还运用了支持向量机在主成份空间中高效地分析非线性关系，提高了判读的精确性和有效性。但是，该发明计算方法复杂，能够处理的频率范围有限，不能满足更宽频段的需求。

申请号为200710163302.0的发明公开了一种脑电超慢波分析方法及其数据处理装置，它通过对0.5-50Hz范围内的脑电波进行系列分析，提取出其中超慢波，并分析超慢谱系功率，用以反映所对应的脑神经递质的功能状况，为建立一种无创伤脑神经递质功能检测方法提供了条件。但是，该发明对于脑电数据特征具有针对性，当特征值不一致时效果较差。

通过对脑电数据信号的监测和分析，可以得到大脑的状态，如深度睡眠，大脑疲劳，大脑警觉等状态，同时，也为人类与计算机的交互控制提供了新的手段。

对脑电数据信号进行特征提取，得到可识别的，可靠程度高的特征值是脑电数据分析的基础环节，具有重要意义。

发明内容

本发明提供了一种脑电特征提取方法，该方法通过提取脑电数据得到时频特征和形态特征，并将时频特征和形态特征相结合，进而反应大脑状态，计算复杂程度低，实时性好。

一种脑电特征提取方法，包括以下步骤：

(1)去除背景脑电数据和待处理脑电数据中的伪迹，分别获得背景脑电数据的有效频段和待处理脑电数据的有效频段，再分别将背景脑电数据的有效频段和待处理脑电数据的有效频段分为若干数据段；

(2)对所得的每个数据段分别提取时频特征和形态特征，得到每个数据段相应的时频特征值和形态特征值；

(3)利用背景脑电数据每个数据段的时频特征值和形态特征值，计算得到时频特征值的频率分布函数和形态特征值的频率分布函数；

(4)利用背景脑电数据的时频特征值的频率分布函数和形态特征值的频率分布函数，得到待处理脑电数据每个数据段的时频特征值出现的概率和形态特征值出现的概率；

(5)利用下式计算得到IMF-VoE特征值，

$\mathrm{IMF}-\mathrm{VoE}\left(t\right)={\log }_2{\mathrm{\Π}}_{c=1}^n{P}_c;$

其中，t为数据段的序数；

n为特征值的数量；

P_c为步骤(4)中待处理脑电数据每个数据段的时频特征值出现的概率或形态特征值出现的概率；

c为特征值的序数。

首先，对脑电数据进行预处理，将原始的脑电数据通过带通滤波去除伪迹，选取有效频段进行分析和处理；优选地，所述步骤(1)中带通滤波后获得的有效频段频率为1.6～70Hz。带通滤波可以采用2阶巴特沃兹滤波器(Butterworth)，有效频段的频率可以根据需要通过滤波参数进行选择。

作为优选，采用滑动时间窗的方法将有效频段分为若干数据段，滑动时间窗的长度为1s，滑动步长为0.2s。

数据分段的目的在于，确定数据处理的最小单元。滑动时间窗的长度以及滑动步长可以根据需要选取，时间窗的长度越短，滑动步长越短，划分的数据段越多，则提取得到的时频特征值和形态特征值越精确，但是，相应的计算量也大。

作为优选，所述步骤(2)中对每个数据段采用经验模态分解法提取时频特征，取前三个固有模态函数，利用下式计算得到每个数据段的时频特征值VoIMF：

${\mathrm{VoIMF}}_N\left(t\right)=\frac{1}{n-1}{\mathrm{\Σ}}_{i=1}^n{(x_i-\stackrel{\‾}{x})}^2,N=\mathrm{1,2,3};$

其中，x_i为每个数据段中每个数据点的固有模态函数值；

为每个数据段中所有数据点的固有模态函数平均值；

i为每个数据段的数据点的序数；

n为每个数据段的数据点的个数；

N为固有模态函数x的序数；

t为数据段的序数。

脑电数据经过经验模态分解(Empirical Mode Decomposition，简称EMD)，被转化为包含了原脑电数据信号的不同时频尺度特征的固有模态函数(Intrinsic Mode Function，简称IMF)，能够反映脑电数据在不同时频尺度上的特征。

一般来说，序数较小的固有模态函数代表频率较高的信号分量，N可以根据需要选择所需的固有模态函数的数目。当N为3时，即选取前三个固有模态函数，分别得到VoIMF₁、VoIMF₂和VoIMF₃。

对于每一个数据段分别计算VoIMF₁、VoIMF₂和VoIMF₃。例如计算VoIMF₁时，选取第一个固有模态函数，计算该数据段中每个数据点的第一个固有模态函数值x_i，计算该数据段中所有数据点的第一个固有模态函数值x_i的平均值得到

利用式

计算得到VoIMF₁，式中，i为每个数据段的数据点的序数；n为每个数据段的数据点的个数；t为数据段的序数。

作为优选，所述步骤(2)中的对每个数据段提取形态特征的步骤如下：

a、对每个数据段进行均值滤波，得到平滑后的数据；

b、求取平滑后的数据段的极大值，连接极大值得到上包络线E_upper；求取平滑后的数据段的极小值，得到下包络线E_lower；利用下式求取上包络线和下包络线的之间的包络范围：

Envelope_Range(t)＝E_upper(t)-E_lower(t)

其中，t为数据段序数；

Envelope_Range(t)为第t个数据段的包络范围；

利用下式计算每个数据段的形态特征值VoE，

$\mathrm{VoE}\left(t\right)=\frac{1}{n-1}{\mathrm{\Σ}}_{i=1}^n{(\mathrm{Envelope}\_\mathrm{Range}{\left(t\right)}_i-\stackrel{\‾}{\mathrm{Envelope}\_\mathrm{Range}\left(t\right)})}^2;$

其中，t为数据段的序数；

n为每个数据段数据点的个数；

i为每个数据段中数据点的序数；

Envelope_Range(t)_i为第t个数据段第i个数据点的包络范围值；

为第t个数据段所有数据点包络范围的平均值。

利用均值滤波(即线性滤波)，对每个数据段进行平滑，去除频率较高的极值，得到平滑过的数据段，均值滤波可以采用滑动时间窗平均法，时间窗的长度和滑动步长可以根据需要进行选择。

作为优选，所述步骤(3)中采用频率分布直方图归一化的方法得到时频特征值的频率分布函数和形态特征值的频率分布函数。

利用背景脑电数据每个数据段的时频特征值(时频特征值的数目依据选择的固有模态函数的数目而定，例如固有模态函数为3个，则相应有3个时频特征值VoIMF₁、VoIMF₂和VoIMF₃)和形态特征值，画出每个特征值(包括3个时频特征值和1个形态特征值)的频率分布直方图，归一化后，得到时频特征值的频率分布函数和形态特征值的频率分布函数，利用得到的频率分布函数，得到待处理脑电数据的时频特征值的概率和形态特征值的概率，利用下式计算得到IMF-VoE特征值，

$\mathrm{IMF}-\mathrm{VoE}\left(t\right)={\log }_2{\mathrm{\Π}}_{c=1}^n{P}_c;$

其中，t为数据段的序数；

n为特征值的数量；

P_c为待处理脑电数据每个数据段的时频特征值出现的概率或形态特征值出现的概率；

c为特征值的序数。

所述特征值的数量包括时频特征值和形态特征值的总数量；所述特征值的序数是对时频特征值和形态特征值依次编号。

本发明一种脑电特征提取方法，将时频特征和形态特征相结合，计算复杂度低，实时性好，可以用于快速识别脑电数据信号的特征变化，从而检测大脑的状态。

附图说明

图1为本发明一种脑电特征提取方法包络范围计算示意图；

图2为本发明一种脑电特征提取方法IMF-VoE特征值检测结果。

具体实施方式

下面结合附图，对本发明一种脑电特征提取方法做详细描述：

一种脑电特征提取方法，包括以下步骤：

(1)利用带通滤波去除脑电数据伪迹，获得有效频段，利用滑动时间窗将有效频段均匀分为若干数据段；

首先，对脑电数据进行预处理，将原始的脑电数据通过带通滤波去除伪迹，选取有效频段进行分析和处理，带通滤波采用2阶巴特沃兹滤波器(Butterworth)，滤波参数为1.6～70Hz，即得到的有效频段频率为1.6～70Hz。

将得到的有效频段均匀分为若干数据段，确定数据处理的最小单元，采用滑动时间窗的方法对得到的有效频段进行分解，滑动时间窗的长度为1s，滑动步长为0.2s。

例如，采样频率为1000Hz，滑动时间窗的长度为1s，则每个时间窗内包含1000个数据点，滑动步长为0.2s，即下一个时间窗的起点为当前时间窗起点后的0.2s，即当前时间窗的第201个数据点为下一个时间窗的第1个数据点。

(2)对步骤(1)所得的每个数据段分别提取时频特征和形态特征，得到每个数据段的时频特征值和形态特征值；

提取时频特征采用现有的经验模态分解方法，首先，将数据段分解成为具有不同时频特征的固有模态函数，再通过计算固有模态函数的方差得到时频特征值。

经验模态分解步骤如下：

首先，求取第t个数据段中所有极大值点，并利用三次样条插值函数拟合形成上包络线；求取第t个数据段中所有的极小值点，并利用三次样条插值函数拟合形成下包络线；计算上包络线和下包络线的平均值，得到第t个数据段的平均数据序列；

其次，第t个数据段的原始数据序列减去平均数据序列，得到第t个数据段的最终数据序列；

最后，若得到的第t个数据段的最终数据序列仍然存在负的极大值和正的极小值，说明这还不是一个本征模函数，重复前述步骤，直至最终数据序列满足固有模态函数要求。

一般来说，序数较小的固有模态函数代表频率较高的信号分量，根据数据处理经验，选取前三个固有模态函数作为有效时频特征，利用下式计算得到时频特征值；

${\mathrm{VoIMF}}_N\left(t\right)=\frac{1}{n-1}{\mathrm{\Σ}}_{i=1}^n{(x_i-\stackrel{\‾}{x})}^2,N=\mathrm{1,2,3};$

其中，x_i为每个数据段中每个数据点的固有模态函数值；

为每个数据段中所有数据点的固有模态函数平均值；

i为每个数据段的数据点的序数；

n为每个数据段的数据点的个数；

N为固有模态函数x的序数；

t为数据段的序数。

N为固有模态函数x的序数，这里N取3，即可以得到3个时频特征值，分别为VoIMF₁、VoIMF₂、VoIMF₃。

对每个数据段均做时频特征的提取，得到相应的时频特征值VoIMF₁、VoIMF₂、VoIMF₃。

对每个数据段提取形态特征的步骤如下：

1)采用时间窗平均法，时间窗长度为0.02s，步长为0.02s；将数据窗中点的数据值赋值为时间窗均值，再通过三次样条插值方法得到平滑过的数据段。

2)求取平滑后的数据段的极大值，连接极大值得到上包络线E_upper；求取平滑后的数据段的极小值，得到下包络线E_lower；利用下式求取上包络线和下包络线的之间的包络范围：

Envelope_Range(t)＝E_upper(t)-E_lower(t)

其中，t为数据段序数；

Envelope_Range(t)为第t个数据段的包络范围；

利用下式计算每个数据段的形态特征值VoE，

$\mathrm{VoE}\left(t\right)=\frac{1}{n-1}{\mathrm{\Σ}}_{i=1}^n{(\mathrm{Envelope}\_\mathrm{Range}{\left(t\right)}_i-\stackrel{\‾}{\mathrm{Envelope}\_\mathrm{Range}\left(t\right)})}^2;$

其中，t为数据段的序数；

n为每个数据段数据点的个数；

i为每个数据段中数据点的序数；

Envelope-Range(t)_i为第t个数据段第i个数据点的包络范围值；

为第t个数据段所有数据点包络范围的平均值。

如图1所示，图1中上面的图为原始脑电数据，下面的图为计算包络范围的示意图，点划线和虚线之间的范围即为包络范围。

(3)计算步骤(2)中得到的每个数据段的时频特征值的概率和形态特征值的概率；

将时频特征和形态特征相结合，即通过将时频特征值VoIMF和形态特征值VoE相结合，得到IMF-VoE特征值，利用IMF-VoE特征值得到大脑的状态，如深度睡眠、疲劳等。

IMF-VoE特征值的计算步骤如下：

A、选取静息状态下的脑电数据(即背景脑电数据)，通过步骤(1)和步骤(2)得到背景脑电数据的时频特征值VoIMF和形态特征值VoE，利用频率分布直方图计算时频特征值VoIMF和形态特征值VoE的频率分布函数，具体方法如下：

A-1、选取相对平稳的背景信号，最少时间长度为10min；

A-2、利用滑动时间窗法将背景信号平均分为若干数据段，计算每个数据段的时频特征值VoIMF₁、VoIMF₂、VoIMF₃和形态特征值VoE特征值；其中，VoIMF₁即利用经验模态分析得到的第一固有模态函数，VoIMF₂即利用经验模态分析得到的第二固有模态函数，VoIMF₃即利用经验模态分析得到的第三固有模态函数；

A-3、对于每一个特征值(包括3个时频特征值和1个形态特征值)，分别确定其分布区间[a，b]，设定组数m，则组距为h＝(a-b)/m，分别做特征值的频率分布直方图；若某个特征值的频率为0，则将其置为一个小于1的值以方便计算，一般取0.01-0.1之间。

A-4、对频率分布直方图进行归一化，得到背景脑电数据每一个特征值的频率分布函数；

B、得到背景脑电数据的特征值的分布函数后，对于待处理脑电数据，IMF-VoE特征值的计算方法如下：

B-1、计算待处理脑电数据的时频特征值VoIMF₁、VoIMF₂、VoIMF₃和形态特征值VoE；

B-2、利用背景信号特征值的分布函数计算每个特征值的出现概率P_c；

B-3、利用下述公式计算得到IMF-VoE特征值：

$\mathrm{IMF}-\mathrm{VoE}\left(t\right)={\log }_2{\mathrm{\Π}}_{c=1}^n{P}_c;$

其中，t为数据段的序数；

n为特征值的数量，特征值分别为时频特征值VoIMF₁、VoIMF₂、VoIMF₃和形态特征值VoE，n为4；

P_c为每个特征值出现的概率；包括时频特征值VoIMF₁的概率、VoIMF₂的概率、VoIMF₃的概率和形态特征值VoE的概率；

c为特征值的序数。

通过IMF-VoE特征值对头皮脑电信号进行分析，结果表明IMF-VoE特征值对脑电信号的时频特征变化和形态特征变化敏感。

当IMF-VoE特征值小于一定阈值，表明待处理脑电数据段属于背景脑电数据为小概率事件，即待处理脑电数据有异于背景脑电数据。

慢波是睡眠监测等大脑状态分析过程中重要的参考依据，慢波发生时，脑电信号的时频特征和形态特征均发生变化，因此通过IMF-VoE特征值，可以监测大脑状态。

如图2所示，选取长度为10min的静息状态下的脑电做背景信号(图2中最上面的图为原始脑电数据)，当慢波脑电信号出现时，时频特征值VoIMF₁、VoIMF₂、VoIMF₃和形态特征值VoE迅速升高，这些特征值属于背景信号的概率降低，因此，IMF-VoE特征值显著减小。

因此，通过IMF-VoE特征值，可以对深度睡眠等不同于静息状态的脑电数据分析提供参考。

Claims (3)

1.一种脑电特征提取方法，包括获取背景脑电数据和待处理脑电数据，其特征在于，还包括以下步骤：

（1）去除背景脑电数据和待处理脑电数据中的伪迹，分别获得背景脑电数据的有效频段和待处理脑电数据的有效频段，再分别将背景脑电数据的有效频段和待处理脑电数据的有效频段分为若干数据段；

采用带通滤波去除伪迹，去除伪迹后获得的有效频段频率为1.6～70Hz；

采用滑动时间窗的方法将有效频段分为若干数据段，滑动时间窗的长度为1s，滑动步长为0.2s；

（2）对所得的每个数据段分别提取时频特征和形态特征，得到每个数据段相应的时频特征值和形态特征值；

（3）利用背景脑电数据每个数据段的时频特征值和形态特征值，计算得到时频特征值的频率分布函数和形态特征值的频率分布函数；

采用频率分布直方图归一化的方法得到时频特征值的频率分布函数和形态特征值的频率分布函数；

（4）利用背景脑电数据的时频特征值的频率分布函数和形态特征值的频率分布函数，得到待处理脑电数据每个数据段的时频特征值出现的概率和形态特征值出现的概率；

（5）利用下式计算得到IMF-VoE特征值，

$\mathrm{IMF}-\mathrm{VoE}\left(t\right)={\log }_2{\mathrm{\Π}}_{c=1}^n{P}_c;$

其中，t为数据段的序数；

n为特征值的数量；

P_c为步骤（4）中待处理脑电数据每个数据段的时频特征值出现的概率或形态特征值出现的概率；

c为特征值的序数。

2.如权利要求1所述脑电特征提取方法，其特征在于，所述步骤（2）中对每个数据段采用经验模态分解法提取时频特征，取前三个固有模态函数，利用下式计算得到每个数据段的时频特征值VoIMF：

$\mathrm{VoI}{\mathrm{MF}}_N\left(t\right)=\frac{1}{n-1}{\mathrm{\Σ}}_{i=1}^n{(x_i-\stackrel{\‾}{x})}^{2}N=\mathrm{1,2,3};$

其中，x_i为每个数据段中每个数据点的固有模态函数值；

为每个数据段中所有数据点的固有模态函数平均值；

i为每个数据段的数据点的序数；

n为每个数据段的数据点的个数；

N为固有模态函数x的序数；

t为数据段的序数。

3.如权利要求2所述脑电特征提取方法，其特征在于，所述步骤（2）中的对每个数据段提取形态特征的步骤如下：

a、对每个数据段进行均值滤波，得到平滑后的数据；

b、求取平滑后的数据段的极大值，连接极大值得到上包络线E_upper；求取平滑后的数据段的极小值，得到下包络线E_lower；利用下式求取上包络线和下包络线的之间的包络范围：

Envelope_Range(t)＝E_upper(t)-E_lower(t)

其中，t为数据段序数；

Envelpe_Range(t)为第t个数据段的包络范围；

利用下式计算每个数据段的形态特征值VoE，

$\mathrm{VoE}\left(t\right)=-\frac{1}{n-1}{\mathrm{\Σ}}_{i=1}^n{(\mathrm{Envelope}\_\mathrm{Rang}{\left(t\right)}_i-\stackrel{\‾}{\mathrm{Envelope}\_\mathrm{Range}\left(t\right)})}^2;$

其中，t为数据段的序数；

n为每个数据段数据点的个数；

i为每个数据段中数据点的序数；

Envelope_Range(t)_i为第t个数据段第i个数据点的包络范围值；

为第t个数据段所有数据点包络范围的平均值。

Images