CN103110418B - 一种脑电信号特征提取方法 - Google Patents

Abstract

本发明提供一种脑电信号特征提取方法，通过小波重构、加窗水平可视图复杂网络转换和复杂网络分析的方法计算网络平均路径长度和聚类系数的特征提取算法，计算脑电信号构成的复杂网络的平均路径长度和网络聚类系数，从而实现脑电信号及不同节律脑电信号的混沌时间序列信号特性分析。本发明的效果该方法是将一维混沌时间序列转化为复杂网络，通过网络特征参数的分析，揭示脑电信号的分形特性，从一个全新的角度对脑电信号这种复杂的非线性随机信号进行刻画。该脑电信号特征可以应用于精神疾病自动诊断和脑机接口***的特征识别模块。该方法能够有效区分癫痫发作期和非发作期的脑电信号，经Mann–Whitney检验p<0.01,可以应用于癫痫脑电自动识别。

Description

一种脑电信号特征提取方法

技术领域

本发明属于脑电信号处理方法，特别是一种脑电信号特征提取方法。

背景技术

脑电信号是脑神经细胞电生理活动在大脑皮层或头皮表面的总体反映，包含大量的生理与病理信息，是非线性的混沌时间序列。脑电信号不仅为脑疾病提供诊断和治疗的依据；而且成为人脑语言、记忆、学习和思维等任职功能的重要研究方法；在工程应用方面,通过脑电信号识别与控制的实现脑-计算机接口也成为研究的热点。脑电信号分析和特征提取是诊断***和控制***中为模式识别提供客观参数的重要环节。

脑电信号特征提取的方法包括时域方法、频域方法、时-频分析方法以及近年来发展的非线性分析方法。

专利（申请号：200810056838.7）公开了一种基于小波变换和BP神经网络的脑电特征提取方法，以想象动作思维引起的能量变换作为想象左右手运动思维特征，应用于脑机接口***。

专利（申请号：200910196746.3）公开了一种脑电分析方法，应用时频分析和主成分分析提取脑电信号时频参数，并将其映射到主成分空间，运用支持向量机分析非线性关系。

专利（申请号：201210235865.7）公开了一种基于定量脑电图的脑电特征提取方法，将采集的信号进行预处理，然后提取功率耦合系数和功率谱密度的不对称性系数，最后利用支持向量机对静态脑电信号进行分类。

专利（申请号：201210308790.0）公开了一种脑电信号特征提取方法。该方法通过数据预处理消除信号伪迹和选择有效频段后，将脑电数据分段并提取时频特征和形态特征，利用时频特征和形态特征的频率分布函数计算特征值IMF-VoE。

专利（申请号：201210245452.7）公开了一种针对抑郁症PSD患者脑电信号的特征提取方法，通过功率谱不对称性对脑电特征进行分类识别，从而对被测试者是否抑郁及抑郁程度进行判别。

专利（申请号：201210189995.1）公开了一种基于脑功能网络邻接矩阵分解的脑电特征提取方法，对多通道脑电进行相关性分析建立脑功能网络，然后进行基于邻接矩阵奇异值分解的特征提取，最后运用支持向量机进行分类。

前五种方法均为基于时频域分析的特征提取方法，注重于激活脑区的脑电信号功率的定量及定性分析，不能体现脑电信号混沌性特性。第六种方法虽然采用了网络分析的方法计算了不同脑区脑电信号的相关性，不能反映脑电这种混沌信号的在单一脑区表现的时间序列随机信号内在规律的分形特性。

本发明针对现有方法无法揭示脑电信号混沌特性的问题，提出一种小波变换和时间序列复杂网络分析相结合的方法，不仅能够分析脑电信号和不同节律脑电信号作为混沌时间序列信号的内在规律的分形特性。

发明内容

针对现有脑电信号视频域特征提取方法不能体现脑电信号混沌性特性且复杂网络方法建网复杂的现状，本发明的目的是提供一种脑电信号特征提取方法，通过将小波重构、加窗水平可视图复杂网络转换和复杂网络分析计算网络平均路径长度和聚类系数的特征提取算法，从而实现脑电信号及不同节律脑电信号的混沌信号网络特性分析。该方法可以作为神经疾病自动诊断***和脑机接口控制***的特征提取环节的一个组成部分。

为实现上述目的，本发明采用的技术方案是提供一种脑电信号特征提取方法，该方法包括有以下步骤：

①小波重构

通过脑电图机采集的脑电信号，对脑电信号进行四级小波重构，第一级分解将采集的脑电信号从频率中点分开为低频、高频两部分，第二级分解将第一级分解得到的低频成分从其频率中点分开为低频、高频两部分，第三级分解将第二级分解得到的低频成分从其频率中点分开为低频、高频两部分，第四级分解将第三级分解得到的低频成分从其频率中点分开为低频、高频两部分，再通过第四级小波重构即可得到脑电信号的五种生理节律；

②加窗水平可视图复杂网络转换

脑电信号的采样频率为f，令时间窗时间长度为t₁,时间窗滑动距离为t₂，则时间窗内采样点个数为N_w=t₁×f，时间窗滑动的采样点个数为N_s=t₂×f，令时间窗内的脑电信号时间序列为{x_i}i=1,2,...,N_w，令i和j是复杂网络中的两个节点，如果对于所有的x_nn∈(i,j)，x_i,x_j>x_n，则认为i和j有链接，否则没有连接，从而建立N_w×N_w维复杂网络邻接矩阵M，其中，当i和j有链接时，M(i,j)=1；否则M(i,j)=0；

对步骤①中的脑电信号数据和所有小波重构获得的各节律脑电数据，或者根据需要选取通过小波重构获得的各节律脑电信号中的任意节律的脑电信号，按照上述方法尽心网络转换，得到脑电数据和各节律脑电数据对应的复杂网络邻接矩阵M；

③复杂网络参数特征提取；

对步骤②得到的复杂网络邻接矩阵M提取复杂网络的网络特性路径长度L和聚类系数C；

网络的平均路径长度定义为：

$L=\frac{1}{N(N-1)}\underset{i=2}{\overset{N}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{i-1}{\mathrm{\&Sigma;}}}{d}_{\mathrm{ij}}---\left(3\right)$

即网络中任意两节点的距离的平均值，其中N为网络节点个数，其中，d_ij定义为连接两个节点i和j的最短路径上的边数，网络的聚类系数定义为所有节点聚类系数的平均值：

$C=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{C}_i---\left(4\right)$

其中N为网络节点个数，其中，C_i即为节点i的聚类系数，是当节点i的度为k_i时，k_i个节点间实际存在的边数R_i与最多存在的边数E_i=k_i(k_i-1)/2的比值；

将步骤②得到的脑电信号和各节律脑电信号对应的网络矩阵按照公式（1）和公式（2）计算，得到相应的网络特征参数：平均路径长度L和平均聚类系数C；

④通过步骤③，实时获得脑电信号及各节律脑电信号时间窗内数据的网络特征参数平均路径长度L和平均聚类系数C，平均路径长度L和平均聚类系数C作为脑电信号的特征值。

本发明的效果该方法是将一维混沌时间序列转化为复杂网络，通过网络特征参数的分析，揭示脑电信号的分形特性，从一个全新的角度对脑电信号这种复杂的非线性随机信号进行刻画。该脑电信号特征可以应用于精神疾病自动诊断和脑机接口***的特征识别模块。该方法能够有效区分癫痫发作期和非发作期的脑电信号，经Mann–Whitney检验p<0.01,可以应用于癫痫脑电自动识别。

附图说明

图1为本发明的特征提取算法原理框图；

图2为利用本发明方法的癫痫脑电信号；

图3为利用本发明方法的癫痫发作前期的脑电信号和小波重构的各节律脑电信号；

图4为利用本发明方法的癫痫脑电信号对应的网络平均路径长度；

图5为利用本发明方法的癫痫脑电信号对应的网络聚类系数。

具体实施方式

如图1所示，本发明的一种脑电信号特征提取方法，该方法包括有以下步骤：小波重构、加窗水平可视图复杂网络转换和复杂网络参数特征提取。

该方法处理的脑电数据可以是通任何脑电采集装置获得的脑电数据。该算法可以在计算机等通用数据处理平台进行，也可以集成于脑机接口***等专用的数据处理装置进行。对于利用本方法处理的脑电数据可以在任何平台上使用，可以集成于精神疾病自动诊断***，也可以集成于脑机接口***，也可以应用于任何需要进行脑电信号特征值提取的需求中。

下面对这四个组成部分进行逐一说明：

1.小波重构

通过脑电图机采集的脑电信号，对脑电信号进行四级小波重构，基本过程为：第一级分解将原信号从频率中点分开为低频、高频两部分，第二级分解将第一级分解得到的低频成分从其频率中点分开为低频、高频两部分，第三级分解将第二级分解得到的低频成分从其频率中点分开为低频、高频两部分，第四级分解将第三级分解得到的低频成分从其频率中点分开为低频、高频两部分。再通过第四级小波重构即可得到脑电信号的五种生理节律。

脑电信号的五种生理节律分别为delta节律（约0-4Hz）、theta节律（约4-8Hz）、alpha节律（约8-12Hz）、beta节律（约13-30Hz）以及gamma节律（约30-60Hz），所以将已经通过0-60Hz的低通滤波器的脑电信号经过四级小波重构即可得到脑电信号的五种生理节律。

2.加窗的水平可视图复杂网络转换

水平可视图方法是一种将一维时间序列转化为复杂网络的建网方法。该方法将时间序列的每个数据点看作复杂网络的一个节点，如果两节点之间满足可视性准则两节点间存在连接，即如果两节点能相互水平看到对方即认为两节点间存在连接。可视性准则描述为：如果(t_a,x_a)和(t_b,x_b)两节点间的任一点(t_c,x_c)均满足x_a,x_b>x_c，则这两个节点存在连接。

因为对于时间序列两个数据点间的时间间隔越长其相关性越小，为了消除网络的时间不相关性冗余及体现网络拓扑随时间的变化关系，本发明的方法采用了加窗的水平可视图方法。在矩形时间窗内利用水平可视性准则判断两节点间是否存在连接，时间窗每前进一步重复上面的计算，得到相应的网络邻接矩阵。

基本过程为：脑电的采样频率为f，令时间窗时间长度为t₁,时间窗滑动距离为t₂，则时间窗内采样点个数为N_w=t₁×f，时间窗滑动的采样点个数为N_s=t₂×f。令时间窗内的脑电信号时间序列为{x_i}，i=1,2,...,N_w。令i和j是复杂网络中的两个节点。如果对于所有的x_nn∈(i,j)，x_i,x_j>x_n，则认为i和j有链接，否则没有连接。从而建立N_w×N_w维复杂网络邻接矩阵M。其中，当i和j有链接时，M(i,j)=1；否则M(i,j)=0。

3.复杂网络参数特征提取

提取复杂网络M的网络特性路径长度L和聚类系数C。

网络的平均路径长度定义为

$L=\frac{2}{N(N-1)}\underset{i=2}{\overset{N}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{i-1}{\mathrm{\&Sigma;}}}{d}_{\mathrm{ij}}---\left(5\right)$

即网络中任意两节点的距离的平均值，其中N为网络节点个数。其中，d_ij定义为连接两个节点i和j的最短路径上的边数。

网络的聚类系数定义为所有节点聚类系数的平均值

$C=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{C}_i---\left(6\right)$

其中N为网络节点个数。其中，C_i即为节点i的聚类系数，是当节点i的度为k_i时，k_i个节点间实际存在的边数R_i与最多存在的边数E_i=k_i(k_i-1)/2的比值。

4.特征值

通过步骤3得到的网络特性路径长度L和聚类系数C，即为对应时间窗内脑电信号的特征值。

下面以该提取方法在癫痫脑电信号分析中的应用为例对本发明的方法进行说明。

1.实验数据

实验数据采用医院提供的20例包含癫痫强直性阵挛发作前后C4导的三分钟数据脑电数据，数据采样频率256Hz，已通过数据预处理去除了伪迹和0-60Hz低通滤波，如图2所示。其中0-80秒为癫痫发作前期脑电信号，80-160秒为癫痫发作期脑电信号，160-180秒为癫痫发作后期脑电信号。

2.算法实现

2.1小波重构

采用长度为8的Daubechies最小非对称相移滤波器（Daubechies leastasymmetric filter of length8）对脑电信号进行四级分解得到EEG信号的5种生理节律成分。实验使用MATLAB7.10.0版本的Wavelet工具箱中的wavedec函数，得到各个生理频段成分的小波系数；然后使用waverec函数对各个生理节律成分进行重构，将小波系数通过离散小波逆变换转换为各生理节律的时间序列。小波重构效果如图3所示，纵坐标的单位为微伏，横坐标中的数字是截取的数据点的序号。

2.2加窗水平可视图复杂网络转换

分别对癫痫脑电信号及五种节律脑电信号（delta节律、theta节律、alpha节律、beta节律和gamma节律）进行加窗水平可视图复杂网络转换。具体方法如下。

令时间窗长度约为6秒,时间窗点数N_w=1536点，；令每个时间窗的间隔时间约为2秒,时间窗滑动距离为512点。因此180秒脑电数据生成87个网络邻接矩阵。时间窗起始时间t_i和时间窗编号N_i的对应关系为t_i=2N_i。

2.3网络特征参数提取

根据公式（1）和（2）对癫痫患者脑电数据及五种节律脑电数据生成的网络邻矩阵计算每个时间窗内网络的聚类系数和平均路径长度。网络聚类系数和网络平均路径长度即为该脑电信号的特征值。

3.结果与分析

图4和图5是一例癫痫脑电信号时间序列提出的网络平均路径长度和网络聚类系数。图4纵坐标是网络平均路径长度L，横坐标是时间窗编号N_i和对应的时间窗起始时间t_i，t_i单位为秒，脑电信号在癫痫发作前40秒开始的平均路径长度增加，在癫痫大发作结束后平均路径长度降低。图5纵坐标是网络聚类系数C，横坐标是时间窗编号N_i和对应的时间窗起始时间t_i，t_i单位为秒，脑电信号在癫痫发作时网络聚类系数增加，在癫痫大发作结束后聚类系数降低。

对20例癫痫脑电信号进行发作前期和发作期，发作后期和发作期的脑电及各节律脑电信号的Mann–Whitney检验，p均小于0.01，差异显著。

因为在癫痫发作前，发作中和发作后网络平均路径长度和网络平均聚类系数差异显著，所以该特征值可以用于癫痫脑电自动识别。

另外，在其他精神疾病及其他意识行为中，脑电信号的时间序列混沌特性都会发生相应变化，所以特征提取方法也可以集成于其他精神疾病自动诊断***和脑机接口***作为特征参数提取环节的一个组成部分。

Claims (1)

1.一种脑电信号特征提取方法，该方法包括有以下步骤：

①小波重构

通过脑电图机采集的脑电信号，采用长度为8的Daubechies最小非对称相移滤波器对脑电信号进行四级分解得到EEG信号的5种生理节律成分，即第一级分解将采集的脑电信号从频率中点分开为低频、高频两部分，第二级分解将第一级分解得到的低频成分从其频率中点分开为低频、高频两部分，第三级分解将第二级分解得到的低频成分从其频率中点分开为低频、高频两部分，第四级分解将第三级分解得到的低频成分从其频率中点分开为低频、高频两部分，然后使用MATLAB7.10.0版本的Wavelet工具箱中的wavedec函数，计算各个生理频段成分的小波系数；使用waverec函数对各个生理节律成分进行重构，将小波系数通过离散小波逆变换转换为各生理节律的时间序列；

②加窗水平可视图复杂网络转换

脑电信号的采样频率为f，令时间窗时间长度为t₁，时间窗起始时间t_i和时间窗编号N_i的对应关系为t_i＝2N_i，时间窗滑动距离为t₂，则时间窗内采样点个数为N_w＝t₁×f，时间窗滑动的采样点个数为N_s＝t₂×f，令时间窗内的脑电信号时间序列为{x_i}i＝1,2,...,N_w，令i和j是复杂网络中的两个节点，如果对于所有的x_nn∈(i,j)_，x_i,x_j>x_n，则认为i和j有链接，否则没有连接，从而建立N_w×N_w维复杂网络邻接矩阵M，其中，当i和j有链接时，M(i,j)＝1；否则M(i,j)＝0；

对步骤①中的脑电信号数据和所有小波重构获得的各节律脑电数据，或者根据需要选取通过小波重构获得的各节律脑电信号中的任意节律的脑电信号，按照上述方法尽心网络转换，得到脑电数据和各节律脑电数据对应的复杂网络邻接矩阵M；

③复杂网络参数特征提取；

对步骤②得到的复杂网络邻接矩阵M提取复杂网络的网络特性路径长度L和聚类系数C；

网络的平均路径长度定义为：

$L=\frac{2}{N(N-1)}\underset{i=2}{\overset{N}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{i-1}{\mathrm{\&Sigma;}}}{d}_{\mathrm{ij}}---\left(1\right)$

即网络中任意两节点的距离的平均值，其中N为网络节点个数，其中，d_ij定义为连接两个节点i和j的最短路径上的边数，网络的聚类系数定义为所有节点聚类系数的平均值：

$C=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{C}_i---\left(2\right)$

其中N为网络节点个数，其中，C_i即为节点i的聚类系数，是当节点i的度为k_i时，k_i个节点间实际存在的边数R_i与最多存在的边数E_i＝k_i(k_i-1)/2的比值；

将步骤②得到的脑电信号和各节律脑电信号对应的网络矩阵按照公式(1)和公式(2)计算，得到相应的网络特征参数：平均路径长度L和平均聚类系数C；

④通过步骤③，实时获得脑电信号及各节律脑电信号时间窗内数据的网络特征参数平均路径长度L和平均聚类系数C，平均路径长度L和平均聚类系数C作为脑电信号的特征值。