CN112932505A - 一种基于时频能量的符号传递熵及脑网络特征计算方法 - Google Patents

Abstract

本发明公开了一种基于时频能量的符号传递熵及脑网络特征计算方法，首先，基于共平均参考对采集的运动想象脑电信号(MI‑EEG)进行预处理；然后，对各导联MI‑EEG进行连续小波变换，求得其时‑频‑能量矩阵，并将与运动想象密切相关的频带内各频率所对应的时间‑能量序列依次拼接，得到该导联的一维时频能量序列；进而，计算任意两个导联时频能量序列之间的符号传递熵，构建大脑连通性矩阵，并使用皮尔逊特征选择算法优化矩阵元素；最后，计算脑功能网络的度和中间中心性，构成特征向量，用于MI‑EEG的分类。结果表明，本发明可以有效地提取MI‑EEG的频域特征和非线性特征，相比于传统的基于脑功能网络的特征提取方法具有明显的优势。

Description

一种基于时频能量的符号传递熵及脑网络特征计算方法

技术领域

本发明属于脑电信号处理领域，涉及一种基于时频能量的符号传递熵及脑网络特征计算方法，应用于脑-机接口系统中运动想象脑电信号(MI-EEG)的特征提取。具体涉及：基于连续小波变换和符号传递熵构建动态脑功能网络，并结合皮尔逊特征选择算法优选网络特征，用于不同运动想象脑电信号的识别。

背景技术

人类的大脑是一个复杂而密集的网络，由数十亿相互连接的神经元组成。近年来，基于图论的复杂网络分析方法被应用于神经科学中，利用复杂网络基本原理可以进行大脑属性分析，以及发现脑网络节点间潜在的信息传递关系。越来越多的研究表明，利用脑功能网络中的度量来区分运动想象脑电信号具有一定的可行性。在脑功能网络中，节点通常是EEG电极、而边定义为节点之间的连通性。计算连通性的指标可分为两类：功能连接和有效连接。功能连接度量(如互相关(CC)、互信息(MI)等)用于捕获不同电极之间的相互依赖性，但无法提供具有流动方向的信息；而有效连接度量(如格兰杰因果关系(GC)、定向传递函数(DTF)等)，可描述节点之间的因果相互作用和信息传输方向。因此，使用有效连接度量更容易发现运动想象过程中节点之间隐藏和被忽视的连接性。然而，GC和DTF依赖于线性回归模型，这使得它们难以准确识别出人脑这样的非线性动力系统中存在的因果关系。

传递熵(TE)是一种基于非参数统计的信息理论量，它不需要预先假设交互作用的模型(无论是线性的还是非线性的)，可用于衡量两个非线性向量之间信息交互，是研究大脑信息流的理想方法。研究人员已将其应用于各种领域，如识别听觉皮层数据之间的信息传递，定位癫痫患者的焦点，以及分析不同驾驶状态下的信息流模式等。然而，传统的传递熵需要大量的数据来估计概率分布函数，计算速度较慢且易受噪声影响，这限制了其在实际脑功能网络建模中的应用。此外，现有技术直接通过计算原始脑电信号之间的连通性来构建脑功能网络，这难以获取脑电信号的频域特征，使得基于脑功能网络特征提取方法的最终分类准确率一般较低。

发明内容

针对上述现有技术的不足，本发明提出了一种基于时频能量的符号传递熵及脑网络特征计算方法。首先，基于共平均参考(CAR)对采集的MI-EEG进行预处理，去除空间噪声；然后，对各通道MI-EEG进行连续小波变换，求得其时-频-能量矩阵，并将与运动想象密切相关的频带内各频率所对应的时间-能量序列依次拼接，得到该通道的一维时频能量序列；进而，计算任意两个通道时频能量序列之间的符号传递熵，构建大脑连通性矩阵，并使用皮尔逊特征选择算法优化矩阵元素；最后，计算脑功能网络的节点度和节点中间中心性，构成特征向量，并采用支持向量机(SVM)对特征向量进行特征评估，根据分类准确率确定最优时间段、频带以及符号的阶次，提取出最优参数下的脑电特征。

本发明的具体实现包括以下步骤：

步骤1：MI-EEG信号预处理。

在预处理阶段，首先对原始MI-EEG信号进行CAR滤波去除空间噪声，并截取运动想象期数据；采用滑动时间窗法将数据划分为多个时间段子信号，分时段对各子信号进行后续分析，窗长为p，步长为q。

步骤2：获取EEG信号的时频能量序列

步骤2.1设

为第i导预处理后的脑电信号，L代表导联数，N_t为采样点数，R代表实数集。对X_i进行Morlet连续小波变换，得到其小波系数矩阵W_i：

式中，

表示第i导脑电信号在m频率下n时刻的小波系数值，N_f为小波分解得到的总频率点数，f₁和f₂为小波分解后感兴趣频带的左、右边界频率，C为虚数集。

步骤2.2截取小波系数矩阵W_i中的第f₁至f₂行，并由公式(2)计算各个元素

的小波能量值

式中，

为

的模长。得到X_i在感兴趣频带内的时频能量矩阵E_i：

其中，N′_f＝f₂-f₁+1为小波分解后感兴趣频带内频率的个数。

步骤2.3将时频能量矩阵E_i各行依次拼接，求得第i导脑电信号的时频能量序列

式中，N_c＝N′_f×N_t为序列

的长度，将(4)式重写为：

步骤3：基于符号传递熵的脑功能网络构建。

步骤3.1对各导联时频能量序列

进行符号化，获得符号序列。

设

为

中D个元素的集合，D≥2，λ∈N^*为延迟因子，t＝1，2，…，N_c--(D-1)λ，N^*为正整数集。当λ＝1时，根据各元素的大小排列情况定义符号

用符号

对

进行标注，得到

的符号序列，记作：

式中，N_m＝N_c-D+1为符号序列M_i的长度。

步骤3.2以EEG电极作为网络的节点，电极时频能量序列之间的符号传递熵作为边，计算有效连通性矩阵，构建脑功能网络。

传递熵用于估计两个序列之间的因果关系或信息流。符号传递熵(STE)是TE的一种变体，计算速度快，抗噪声能力强。符号序列M_i和M_j之间的符号传递熵由下式计算：

式中，

表示符号

和

同时出现的概率；

表示在

和

出现的前提下，

的条件概率；

为仅给定

的情况下，

的条件概率，τ∈N^*为时间滞后。计算任意两个导联时频能量序列之间的符号传递熵，得到脑功能网络的有效连通性矩阵A：

式中，a_ij＝STE(M_i→M_j)表示导联i到导联j的符号传递熵，i，j＝1，2，…，L。

步骤3.3使用皮尔逊特征选择算法计算A中元素a_ij与类别标签之间的相关性，将相关性小于阈值的元素置零，得到新的连通性矩阵A^*∈R^L×L。

步骤4：特征计算与融合。

根据A^*计算脑功能网络特征参数：度、中间中心性。

导联i的度d_i的计算式如下：

d_i＝∑_ja_ij (10)

导联的度越大表明该导联越重要。计算所有导联的度，构成特征向量F₁：

F₁＝[d₁，d₂，…，d_L]∈R^1×L (11)

导联i的中间中心性b_i可由公式(12)计算:

式中，

为导联j与导联k之间的特征路径长度，l_jk(i)表示导联j和导联k经过导联i的特征路径长度。b_i值越高，表明该导联承载的信息量越大，对脑功能网络的影响也越大。计算每个导联的中间中心性，构成特征向量F₂：

F₂＝[b₁，b₂，…，b_L]∈R^1×L (13)

将特征向量F₁和F₂串行融合，获得特征向量F：

F＝[F₁ F₂]∈R^1×2L (14)

与现有技术相比，本发明具有以下优点：

(1)本发明针对传统的脑功能网络构建方法易丢失脑电信号在频域上的信息这一不足，首先，对脑电信号进行连续小波变换，得到其时-频-能量图；之后，将与运动想象密切相关的频带内各频率对应的时间-能量序列依次拼接，得到各导联的一维时频能量序列，并将其作为节点的输入信号来构建脑功能网络，可以更好地体现不同运动想象脑电信号在频域上的差异。

(2)本发明在计算脑功能网络连通性矩阵上使用了符号传递熵，相比于传统的连通性度量，符号传递熵更适合反映出不同脑电极之间的因果交互作用，且这种交互作用是有向的、非线性的，符合人脑的非线性动力学特性。因此，使用符号传递熵构建脑网络能发现运动想象过程中各导联之间隐藏和被忽视的连接性，从而为脑电信号的识别提供了更多的有用特征。

(3)本发明相比于传统的脑功能网络特征提取方法，能更有效地获取运动想象脑电信号的频域特征和非线性特征，同时，结合皮尔逊特征选择算法优选网络参数，进一步提高了脑电信号的分类正确率。在BCI 2000公开数据集的测试中，109名受试者取得了92.28％的平均分类准确率，单个受试者的最高分类准确率为100.00％，相比于其他传统的基于脑功能网络的特征提取方法具有明显的优势。

附图说明

图1为本专利提出方法的流程图。

图2为电极位置分布图。

图3为D＝3时序列符号化的示意图。

图4为本发明实施例有效连通性矩阵图。

图5为本发明实施例中特征参数度的脑地形图。

图6为本发明实施例中特征参数中间中心性的脑地形图。

具体实施方式

本发明具体实验的软件环境为：Windows 10(64位)、Matlab R2017a。

本发明实施例MI-EEG数据来自于BCI 2000公开数据集，使用标准10-20系统分布下的64个电极采集脑电数据，电极分布位置如图2所示。脑电信号采样频率为160Hz，经过1-50Hz滤波及50Hz陷波滤波。该数据集共包含109个受试者，想象任务为左手或右手运动，每名受试者共进行约45次实验。每次实验持续约8秒,其中0至4秒为运动想象期。

基于上述MI-EEG数据集，本发明具体的实施步骤如下所述：

步骤1：信号预处理。

对原始MI-EEG信号进行CAR滤波去除空间噪声，并截取0至4秒的运动想象期数据，采用滑动时间窗法(窗长为2s，步长为1s)将数据划分为三个时间段：T₁(0s-2s)、T₂(1s-3s)和T₃(2s-4s)，每个时间段数据包含320个采样点。

步骤2：获取EEG信号的时频能量序列

步骤2.1分别对各个时间段内的脑电信号进行Morlet连续小波变换，求得其小波系数矩阵：W_i∈C^69×320。

步骤2.2选取α(8-13Hz)、β(13-30Hz)和α+β(8-30Hz)三个频带进行分析，由公式(2)计算小波系数的小波能量值，得到各通道感兴趣频带内的时-频-能量矩阵：E_i∈R^12×320。

步骤2.3，将E_i各行依次进行拼接，得到各个导联的一维时频能量序列

步骤3：基于符号传递熵的脑功能网络构建。

步骤3.1依据公式(6)将各导联的时频能量序列按幅值大小排列情况进行符号化，得到

的符号序列，记为：

步骤3.2以EEG电极为网络的节点，并由公式(8)(式中的τ设为1)计算任意两个导联时频能量序列之间的符号传递熵作为边，得到脑功能网络的有效连通性矩阵：A∈R^64×64。

步骤3.3使用皮尔逊特征选择算法优选矩阵中的元素，将与类别标签之间的相关性小于阈值(阈值设为0.24)的矩阵元素置零，得到新的连通性矩阵：A^*∈R^64×64。

步骤4：特征计算与融合。

由公式(10)和公式(12)分别计算A^*的网络特征参数：度和中间中心性，构成特征向量F₁＝[d₁，d₂，…，d₆₄]∈R^1×64，F₂＝[b₁，b₂，…，b₆₄]∈R^1×64。特征串行融合，得到最终特征向量：F＝[F₁ F₂]∈R^1×128。

为验证本发明提出的基于时频能量的符号传递熵及脑网络特征计算方法的有效性和实用性，使用SVM分类器对本发明计算出的特征向量进行了模式分类并进行了参数寻优。由于每个受试者实验次数较少，使用了留一法对样本进行测试，结果如表1所示。最终得到的最优时间段为T₂，最优频带为β(13-30Hz)，最优的符号阶次为D＝4，此参数下109名受试者

表1最优参数及分类准确率

的平均分类准确率为92.28％，单个受试者的最高分类准确率为100.00％。

Claims (3)

1.一种基于时频能量的符号传递熵及脑网络特征计算方法，其特征在于：该方法包括以下几个步骤，

步骤1：运动想象脑电信号MI-EEG预处理；

使用共平均参考CAR对原始MI-EEG进行空间滤波，并截取运动想象期数据；采用滑动时间窗法将数据划分为多个时间段的子信号，分时段对各子信号进行后续分析，窗长为p，步长为q；

步骤2：获取EEG信号的时频能量序列

步骤3：基于符号传递熵的脑功能网络构建；

步骤4：特征计算与融合；

根据有效连通性矩阵A^*计算脑功能网络特征参数：度、中间中心性；

导联i的度d_i的计算式如下：

d_i＝∑_ja_ij (10)

计算所有导联的度，构成特征向量F₁：

F₁＝[d₁，d₂，…，d_L]∈R^1×L (11)

导联i的中间中心性b_i可由公式(12)计算：

式中，

为导联j与导联k之间的特征路径长度，l_jk(i)表示导联j和导联k经过导联i的特征路径长度；计算每个导联的中间中心性，构成特征向量F₂：

F₂＝[b₁，b₂，…，b_L]∈R^1×L (13)

将特征向量F₁和F₂串行融合，获得特征向量F：

F＝[F₁ F₂]∈R^1×2L (14)。

2.根据权利要求1所述的一种基于时频能量的符号传递熵及脑网络特征计算方法，其特征在于：步骤2中，步骤2.1设

为第i导预处理后的脑电信号，L代表导联数，N_t为采样点数，R为实数集；对X_i进行Morlet连续小波变换，得到其小波系数矩阵W_i：

式中，

表示第i导脑电信号在m频率下n时刻的小波系数值，N_f为小波分解得到的总频率点数，f₁和f₂为小波分解后感兴趣频带的左、右边界频率，C为虚数集；

步骤2.2截取小波系数矩阵W_i中的第f₁至f₂行，并由公式(2)计算各个元素

的小波能量值

式中，

为

的模长；得到X_i在感兴趣频带内的时频能量矩阵E_i：

其中，N′_f＝f₂-f₁+1为小波分解后感兴趣频带内频率的个数；

步骤2.3将时频能量矩阵E_i各行依次拼接，求得第i导脑电信号的时频能量序列

式中，N_c＝N′_f×N_t为序列

的长度，将(4)式重写为：

3.根据权利要求1所述的一种基于时频能量的符号传递熵及脑网络特征计算方法，其特征在于：步骤3中，步骤3.1对各导联的时频能量序列

进行符号化，获得符号序列；

设

为

中D个元素的集合，D≥2，λ∈N^*为延迟因子，t＝1，2，…，N_c-(D-1)λ，N^*为正整数集；当λ＝1时，根据各元素的大小排列情况定义符号

用符号

对

进行标注，得到

的符号序列，记作：

式中，N_m＝N_c-D+1为符号序列M_i的长度；

步骤3.2以EEG电极作为网络的节点，电极时频能量序列之间的符号传递熵作为边，计算有效连通性矩阵，构建脑功能网络；

传递熵TE用于估计两个序列之间的因果关系或信息流；符号传递熵(STE)是TE的一种变体，计算速度快，抗噪声能力强；符号序列M_i和M_j之间的符号传递熵由下式计算：

式中

表示符号

和

同时出现的概率；

表示在

和

出现的前提下，

的条件概率；

为仅给定

的情况下，

的条件概率，τ∈N^*为时间滞后；计算任意两个导联时频能量序列之间的符号传递熵，得到脑功能网络的有效连通性矩阵A：

式中，a_ij＝STE(M_i→M_j)表示导联i到导联j的符号传递熵，i，j＝1，2，…，L；

步骤3.3使用皮尔逊特征选择算法计算A中元素a_ij与类别标签之间的相关性，将相关性小于阈值的元素置零，得到新的连通性矩阵A^*∈R^L×L。

Images