CN108577835A - 一种基于微状态的脑功能网络构建方法 - Google Patents

Abstract

本发明涉及脑功能网络研究技术领域，更具体而言，涉及一种基于微状态的脑功能网络构建方法。本发明从构建EEG脑功能网络出发，根据EEG微状态相关方面的理论，即EEG微状态是反映大脑活动的最基本状态，提出了EEG微状态也能够应用在脑功能网络的构建上。本发明解决了传统构建脑功能网络的时间开销过大问题，通过使用微状态重构后的时间序列，减少了原始EEG时间序列中大量的冗余信息。通过将本发明应用在医院采集的精神疾病工作记忆的EEG数据集，构建了精神分裂症病人和正常人的EEG功能脑网络，证明了用微状态重构时间序列构建的脑功能网络与使用原始EEG信号构建的脑功能相比，能够更加准确的反映人脑网络的特性。

Description

一种基于微状态的脑功能网络构建方法

技术领域

本发明涉及脑功能网络研究技术领域，更具体而言，涉及一种基于微状态的脑功能网络构建方法。

背景技术

复杂网络作为近年来一种新兴的数据分析方法，被应用于各个方面。由于大脑是一个十分复杂的系统，不同神经元之间相互联系构成了一个复杂的脑网络，且这个网络具有非线性、混沌等方面的特点，因此人们逐渐使用复杂网络理论对大脑进行研究。通常脑网络可以分为结构性网络以及功能性网络。结构性网络主要是以大脑的脑区为节点，以它们之间的连接关系作为边，常见于功能性磁共振成像(functional magnetic resonanceimaging，fMRI)网络。功能性脑网络同样是以脑区为节点，以它们之间的相关性作为边，通常脑电图技术(Electroencephalography，EEG)网络则是功能网络。

人脑电活动可以通过EEG获得，具有高时间分辨率以及相当好的空间分辨率。微状态作为一种越来越受欢迎的分析方法，通常被用来分析大脑活动的变化。1987年，Lehmann根据EEG脑地形图的拓扑变化，首次提出了微状态的概念。认为微状态能够反映出全局脑功能状态变化是不连续以及非线性的。2009年，Lehmann根据实验确定一个微状态持续的时间为80-120ms，且不能平滑地转换为另一个状态，是突然变化的。微状态揭示了人脑活动90％的变化可以由四个微状态之间的转换来解释。这四类微状态之间的转换可以代表任何个体中大脑活动的变化。

目前使用EEG构建脑功能网络通常是以电极作为网络的节点，以各个电极在一段时间内的相互联系作为边。然而，由于EEG信号的时间分辨率是毫秒级，且EEG信号是一段连续的脑电信号，所以必然会存在大量的冗余信息。这就会导致计算通道之间联系的时间开销会特别大。微状态能够反映一段时间内大脑整体活动的状态，揭示大脑活动的变化。目前对微状态的研究主要是基于微状态之间的转换以及每个微状态出现的次数，对微状态其他方面的研究很少见。

发明内容

为了克服现有技术中所存在的不足，本发明提供一种基于微状态的脑功能网络构建方法，从复杂脑功能网络构建的角度出发，降低了构建脑功能网络时间复杂度，研究不同微状态时间序列构建的脑网络与微状态本身之间联系。

为了解决上述技术问题，本发明所采用的技术方案为：

一种基于微状态的脑功能网络构建方法，包括以下步骤：

S1、对EEG信号进行预处理，减少电信号的噪声与伪迹，提高信噪比，为后续的特征提取和模式识别结果提供高质量的数据，包括：

1)将EEG原始信号当中的四个活动电极去除掉；

2)将EEG数据的参考电极转换为全脑平均参考电极，脑记录系统高密集分布，电极密度高，故采用“平均参考”，平均参考的优点，在整个球面上求和时，根据欧姆定律，外侧正电流与内侧副电流之和为0，使得其他电极的信息更加明显，便于后续处理和分析；

3)根据采集EEG信号时所依据的实验范式对EEG信号进行分段，分为保持、编码和检索三个阶段；

4)除去步骤3)分段EEG信号中的眼电伪迹；

5)将步骤4)中去眼点后数据进行滤波，分为四个频段；

S2、计算脑电信号的微状态，包括：

1)计算脑电数据的总体场功率(Global field power，简称GFP)：

式中，GFP为总体场功率，x_n为第n个通道的电压观测值，N为通道总数，为所有通道的电压平均值；

2)根据步骤1)所得的总体场功率时间序列，求取它的极大值所对应的时间点，将这些时间点在原始EEG时间序列中所对应的电位信息单独拿出来放在一起，得到总体场功率极大值时间点处所对应的时间序列，因为此时的脑电信号具有最大的信噪比，而位于两个极大值之间的信号，它们的地形图有着相似的结构，GFP极大值点处的地形图能够代替周围时间点；

3)根据K-均值算法(K-means)，对每个被试的总体场功率极大值时间序列进行计算，每个被试的总体场功率极大值时间序列得到四种微状态；

4)将每个被试计算得到的四种微状态重新放入步骤3)当中计算，得到所有被试的四个基本微状态；

5)以步骤4)得到的四个基本微状态为标准，对每个被试的GFP极大值时间序列进行标记，将它们标记为四个微状态中的一个，依据为如下公式：

式中，d_kt表示时间点t处所对应的电位值与第k种微状态之间的距离，V_t表示时间点t处各个通道的电位值，Γ_k代表第k种微状态；

6)据步骤5)得到每个被试标记好的微状态时间序列，对该序列做进一步处理，得到最终的微状态时间序列；

S3、构建EEG脑功能网络，包括：

1)以头皮电极作为网络的节点；

2)以头皮电极之间的相位锁(PLV)作为节点之间的边，边的定义如下：

其中，PLV为相位锁，表示任意两个信号之间平均相位差的绝对值，表示信号x与信号y在t时刻的瞬时相位差值，N为通道数；表示两个信号在各个时间点瞬时相位平均值；

3)计算任意两个头皮电极之间的PLV值，这些PLV值放在一起，可以得到一个N*N的关联矩阵，即脑电信号的脑功能网络，其中N为通道数；

4)将步骤3)中得到的关联矩阵选取合适的阈值构建脑网络，将其转化为二值矩阵，构建无权无向网络。

进一步地，所述EEG信号预处理采用矩阵实验室(Matrix Laboratory，MATLAB)平台上的EEGLAB工具箱。

进一步地，所述S1中步骤4)去除眼电伪迹使用EEGLAB包中的基于独立成分分析(Independent Component Analysis，ICA)算法进行处理。

进一步地，所述S2中步骤3)具体计算算法如下：

1)初始化N_μ个微状态：从GFP极大值时间序列中随机选取N_μ个时间点，将它们的脑地形图作为初始微状态；

2)给定Γ_k，令k＝1…N_μ，满足||Γ_k||＝1，且

其中，Γ_k为第k种微状态，Γ_k为N_μ个时间点所对应通道电位值的归一化；

3)令t＝1…N_t，

其中，N_t为GFP极大值个数，L_t表示每个GFP极大值所对应的微状态，V_t表示每个时刻各通道的电位值，V_t'为向量V_t的转置；

4)令k＝1…N_μ，在||X||＝1的约束条件下，计算计算

其中，S_k为协方差矩阵，X为协方差矩阵S_k的特征向量，为S_k的最大特征值所对应的特征向量；

5)计算

式中，σ₀ ²为噪声方差，N_s为电极数；

6)如果Γ_k即为最终的微状态，若果算法回到步骤3)。

进一步地，所述S2中步骤6)将标记好的微状态时间序列具体处理方法为：判断相同状态所持续时间是否为80ms-120ms，若不满足，则去掉，若满足，则保留，从而得到最终的微状态时间序列。当每一类微状态所持续时间为80ms-120ms时，可以断定该状态是一个稳定的微状态，否则，该状态仅是半稳定状态。

进一步地，所述S3中步骤2)中PLV值的范围为[0，1]。PLV值为1时，表示两个信号相位同步；PLV值为0时，表示两个信号相位不同步。

进一步地，所述频段频率4-7Hz为theta波，频率7-14Hz为alpha波，频率14-20Hz为beta1波，频率20-30Hz为beta2波。

进一步地，所述S3中步骤4)阈值取稀疏度12％～40％，步长为2％。

进一步地，所述无权无向网络模型的规模为60*60。

与现有技术相比，本发明所具有的有益效果为：

本发明涉及一种基于微状态的脑功能网络构建方法，利用微状态能够反映一段时间内大脑状态的特性，将微状态与脑功能网络相结合，使用微状态时间序列构建脑功能网络，在构建EEG脑功能网络时，由于大脑的活动状态可以由微状态之间相互的转换来表示，因此用微状态时间序列代替EEG原始时间序列，在保留EEG信号中主要信息的同时，对数据中冗余信息进行了剔除，降低了构建EEG脑网络的时间成本。

附图说明

图1为四种微状态；

图2为本发明与传统方法的脑网络拓扑图；

图3为本发明与传统分类准确率比较；

图4为保持阶段α频段下正常人与病人聚类系数的脑地形分布图；

图5为EEG信号分段实验范式。

其中图2中，a为传统PLV>0.8的脑网络(病人1)，b本发明PLV>0.8的脑网络(病人1)，c为传统PLV>0.8的脑网络(正常人)，d为本发明PLV>为0.8的脑网络(正常人)，e为传统稀疏度10所形成的网络(病人1)，f为本发明稀疏度10所形成的网络(病人1)，g为传统稀疏度10的脑网络(正常人)，h为微状态稀疏度10的脑网络(正常人)；

图4中A1为微状态-正常人，A2为传统-正常人，A为微状态-传统；B1为微状态-病人，B2为传统-病人，B为微状态-传统；C1为微状态，C2为传统。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

一种基于微状态的脑功能网络构建方法，包括以下步骤：

S1、对EEG信号进行预处理包括：

1)1.去除无用电极：matlab平台上的EEGLAB包将原始信号当中的四个活动电极去除掉；

2)参考电极的转换：使用EEGLAB包对EEG原始数据的参考电极转换为全脑平均参考电极；

3)分段：使用EEGLAB包，根据任务实验范式图对EEG信号进行分段，分为保持，编码，检索三个阶段，实验范式如图5，其中编码阶段为5秒，保持阶段为3秒，检索阶段为2.5秒；

4)去眼电：使用EEGLAB包中的ICA进行处理，用ICA算法分离伪迹与脑电信号，去除眼电的成分；

5)分频段：将步骤4)中去眼点后的数据进行滤波，分为theta(4-7Hz)，alpha(7-14Hz)，beta1(14-20Hz)，beta2(20-30Hz)四个频段；

S2、计算脑电信号的微状态

微状态是大脑活动的头皮表面电位变化地形，被认为是大脑活动的基本状态，通常为A、B、C、D四类。它的计算方法是由如下步骤实现的：

1)计算脑电数据的总体场功率(Global field power，简称GFP)：

式中，GFP为总体场功率，x_n为第n个通道的电压观测值，N为通道总数，为所有通道的电压平均值；

2)根据步骤1)所得到的GFP时间序列，求取它的极大值所对应的时间点，得到GFP极大值时间点处所对应的时间序列，这是因为此时的脑电信号具有最大的信噪比，而位于两个极大值之间的信号，它们的地形图有着相似的结构，GFP极大值点处的地形图能够代替周围时间点；

3)根据改进的K-means算法，对每个被试的GFP极大值时间序列进行计算，得到四个基本的微状态；

4)将每个被试计算得到的四种微状态重新放入步骤3)当中计算，得到本次实验中所有被试的四个基本微状态，具体算法如下：

①初始化N_μ个微状态：从GFP极大值时间序列中随机选取N_μ个时间点，将它们的脑地形图作为初始微状态；

②给定Γ_k，令k＝1…N_μ，满足||Γ_k||＝1，且

其中，Γ_k为第k种微状态，Γ_k为N_μ个时间点所对应通道电位值的归一化；

③令t＝1…N_t，

其中，N_t为GFP极大值个数，L_t表示每个GFP极大值所对应的微状态，V_t表示每个时刻各通道的电位值，V_t'为向量V_t的转置；

④令k＝1…N_μ，在||X||＝1的约束条件下，计算计算

其中，S_k为协方差矩阵，X为协方差矩阵S_k的特征向量，为S_k的最大特征值所对应的特征向量；

⑤计算

式中，σ₀ ²为噪声方差，N_s为电极数；

⑥如果Γ_k即为最终的微状态，如果算法回到步骤3)；

在本实施例中，计算所得最终微状态见图1；

5)以步骤4)得到的四个基本微状态为标准，对每个被试的GFP极大值时间序列进行标记，将它们标记为四个微状态中的一个，依据为如下公式：

式中，d_kt表示时间点t处所对应的电位值与第k种微状态之间的距离，V_t表示时间点t处各个通道的电位值，V_t'为向量V_t的转置，Γ_k代表第k种微状态；

6)根据步骤5)得到每个被试标记好的微状态时间序列，对该序列做进一步处理，判断相同状态所持续时间是否为80ms-120ms，若不满足，则去掉，满足，则保留，从而得到最终的微状态时间序列。这是因为，当每一类微状态所持续时间为80ms-120ms时，可以断定该状态是一个稳定的微状态，否则，该状态仅是半稳定状态；

S3、构建EEG脑功能网络

通常EEG脑功能网络是依据不同电极之间的相互关系来构建，通过量化各个通道之间的联系，来反映整个脑网络，具体步骤如下：

1)节点的选择：以头皮电极作为网络的节点；

2)边的定义：EEG脑功能网络通常是以头皮电极之间的相互关系作为节点之间的边，常用的指标有互信息，相干性，相位同步性，同步似然性等。本发明所使用的是PLV(相位锁)，具体定义如下：

其中，PLV为相位锁，表示任意两个信号之间平均相位差的绝对值，表示信号x与信号y在t时刻的瞬时相位差值，N为通道数；

PLV值介于0到1之间。值为1时，表示两个信号相位同步；值为0时，表示两个信号相位不同步；

3)计算两两头皮电极之间的PLV，可以得到一个60*60的关联矩阵，该矩阵即为脑电信号的脑功能网络；

4)将步骤3)中得到的关联矩阵选取合适的阈值构建脑网络，将其进行二值化处理，得到二值矩阵，也就得到了节点之间连边的关系，构建无权无向网络，网络模型的规模为60*60。

在本实施例中，所述S3中步骤4)阈值取稀疏度12％～40％，步长为2％。

在医院采集的精神疾病工作记忆的EEG数据集，构建了34例精神分裂症病人和34例正常人的EEG功能脑网络，结果分析说明如下：

(1)时间开销

表1.本发明与传统的时间开销之比

由表1可知，本发明在时间开销上远远优于传统构建脑网络的方法，在理论上，使用本发明的时间开销为传统方法的其中n为数据规模。而在实际实验当中，使用本发明的时间开销为20min，传统的时间开销为120min，时间开销为1：6。这是因为在计算过程中要不断的打开与保存文件，会有额外的时间开销，所以在实际实验当中，本发明的时间开销仅为传统方法的1/6。

(2)网络结构

如图2所示，传统方法构建的脑网络拓扑结构与本发明构建的脑网络相比，在相同的阈值下，本发明所构建的网络结构更加清晰，网络的边数也更多。在图a-d中我们可以看到，在阈值为PLV>0.8时，无论是传统网络还是本发明构建的网络，病人脑网络的边主要分布在右前额以及枕区，而正常人脑网络的边则主要分布在右前额，顶区以及枕区。但是相比于传统网络，本发明构建网络的节点之间联系更加紧密，区域更加集中，在接下来对脑网络的分析也更加具有优势。而在图e-h中，也即相同稀疏度下的脑网络当中，可以观察到本发明构建的脑网络，边的分布更加集中，而传统脑网络的边则更加分散。其中，病人脑网络的边主要分布在右前额以及枕区，正常人脑网络的边主要分布在前额，顶区以及枕区。但是两者网络的整体分布趋势还是大致相同的。所以，本发明构建的脑网络能够在保留原始EEG时间序列主要信息的同时，去除大量不必要的冗余信息。且使得原本联系紧密的头皮电极变得更加紧密，排除了数据中冗余信息的干扰。

(3)属性值分析

本次实验分别计算了原始脑电数据构建脑网络的属性值与本发明构建脑网络的属性值。由于需要找出正常人与病人之间的差异性指标，因此需要对计算出来的网络属性进行检验。本实验所采用的检验方法为KS检验，这是因为KS检验是比较一个频率分布f(x)与理论分布g(x)或者两个观测值分布的检验方法。其原假设H0:两个数据分布一致或者数据符合理论分布。D＝max|f(x)-g(x)|，当实际观测值D≥D(n,α)则拒绝H0，否则则接受H0假设。KS检验不需要知道数据的分布情况，是一种非参数检验方法。在样本量比较小的时候，KS检验作为非参数检验在分析两组数据之间分布是否不同时相当常用。由于本试验中样本数量相对较小，所以适合用KS检验来找出差异性比较显著的特征。根据KS检验，找到了每个网络中各个属性值的特异性节点。把他们作为特征放入SVM分类器当中，分类结果如图3，图中柱形从左到右依次为：微状态编码-α-聚类系数，传统编码-α-聚类系数；微状态编码-α-局部效率，传统编码-α-局部效率；微状态编码-θ-聚类系数，传统编码-θ-聚类系数；微状态编码-θ-局部效率，传统编码-θ-局部效率；微状态编码-α-聚类系数，传统编码-α-聚类系数；微状态编码-α-局部效率，传统编码-α-局部效率；微状态编码-θ-聚类系数，传统编码-θ-聚类系数；微状态编码-θ-局部效率，传统编码-θ-局部效率；微状态编码-α-聚类系数，传统编码-α-聚类系数；微状态编码-α-局部效率，传统编码-α-局部效率；微状态编码-θ-聚类系数，传统编码-θ-聚类系数；微状态编码-θ-局部效率，传统编码-θ-局部效率。相比于传统构建的脑网络，本发明的分类准确率有了大大提升。这说明本发明，在有效减少时间开销的同时，能够更加准确的反映大脑属性。进一步说明了本发明的有效性，微状态可以应用于脑网络的构建，为EEG数据分析提供了一个新的思路。

图4为保持阶段α频段下正常人与病人聚类系数的脑地形分布图。我们可以看到，无论是传统方法还是本发明，正常人的聚类系数在右前额以及枕区的值更高，而病人的聚类系数则只在右前额较高。这说明正常人右前额以及枕区的神经元更加活跃，而病人枕区的神经元不太活跃。而对于同一个被试来说，在传统脑功能网络以及本发明构建的脑功能网络当中，它们的聚类系数分布整体上是一致的。且相比较传统方法，本发明构建的脑网络聚类系数地形图的分布层次更加分明。

实验结果表明，本发明使用微状态时间序列构建脑网络的方法可以在减少时间开销的同时保留EEG原始时间序列中的主要信息，且构建的脑功能网络更加准确，能够更好地反映大脑整体的活动状态。用微状态重构时间序列构建的脑功能网络与使用原始EEG信号构建的脑功能相比，能够更加准确的反映人脑网络的特性。使用微状态时间序列构建的脑功能网络在大大减少构建网络时间复杂度的同时，网络的拓扑属性并没有发生很大的变化。并且，正常人与病人之间的差异更加明显，差异性节点更加集中。对研究精神疾病发病机制具有重要意义。

上面仅对本发明的较佳实施例作了详细说明，但是本发明并不限于上述实施例，在本领域普通技术人员所具备的知识范围内，还可以在不脱离本发明宗旨的前提下作出各种变化，各种变化均应包含在本发明的保护范围之内。

Claims (9)

1.一种基于微状态的脑功能网络构建方法，其特征在于，包括以下步骤：

S1、对EEG信号进行预处理，包括：

1)将EEG信号当中的四个活动电极去除掉；

2)将EEG信号的参考电极转换为全脑平均参考电极；

3)对EEG信号进行分段，分为保持、编码和检索三个阶段；

4)除去步骤3)分段EEG信号中的眼电伪迹；

5)将步骤4)中去眼点后数据进行滤波，分为四个频段；

S2、计算脑电信号的微状态，包括：

1)计算脑电数据的总体场功率：

式中，GFP为总体场功率，x_n为第n个通道的电压观测值，N为通道总数，为所有通道的电压平均值；

2)根据步骤1)所得的总体场功率的时间序列，求取它的极大值所对应的时间点，将这些时间点在原始EEG时间序列中所对应的电位信息单独拿出来放在一起，得到总体场功率极大值时间点处所对应的时间序列；

3)根据K-均值算法，对每个被试的总体场功率极大值时间序列进行计算，每个被试的总体场功率极大值时间序列得到四种微状态；

4)将每个被试计算得到的四种微状态重新放入步骤3)当中计算，得到所有被试的四个基本微状态；

5)以步骤4)得到的四个基本微状态为标准，对每个被试的GFP极大值时间序列进行标记，将它们标记为四个微状态中的一个，依据为如下公式：

式中，d_kt表示时间点t处所对应的电位值与第k种微状态之间的距离，V_t表示时间点t处各个通道的电位值，V_t'为向量V_t的转置，Γ_k代表第k种微状态；

6)据步骤5)得到每个被试标记好的微状态时间序列，对该序列做进一步处理，得到最终的微状态时间序列；

S3、构建EEG脑功能网络，包括：

1)以头皮电极作为网络的节点；

2)以头皮电极之间的相位锁作为节点之间的边，边的定义如下：

其中，PLV为相位锁，表示任意两个信号之间平均相位差的绝对值，表示信号x与信号y在t时刻的瞬时相位差值，N为通道数；表示两个信号在各个时间点瞬时相位平均值；

3)计算任意两个头皮电极之间的PLV值，这些PLV值放在一起，可以得到一个N*N的关联矩阵，N为通道数；

4)将步骤3)中得到的关联矩阵选取合适的阈值构建脑网络，将其转化为二值矩阵，构建无权无向网络。

2.根据权利要求1所述的一种基于微状态的脑功能网络构建方法，其特征在于：所述EEG信号预处理采用矩阵实验室平台上的EEGLAB工具箱。

3.根据权利要求1所述的一种基于微状态的脑功能网络构建方法，其特征在于：所述S1中步骤4)去除眼电伪迹使用EEGLAB包中的基于独立成分分析算法进行处理。

4.根据权利要求1所述的一种基于微状态的脑功能网络构建方法，其特征在于：所述频段频率4-7Hz为theta波，频率7-14Hz为alpha波，频率14-20Hz为beta1波，频率20-30Hz为beta2波。

5.根据权利要求1所述的一种基于微状态的脑功能网络构建方法，其特征在于：所述S2中步骤3)具体计算算法如下：

1)初始化N_μ个微状态：从GFP极大值时间序列中随机选取N_μ个时间点，将它们的脑地形图作为初始微状态；

2)给定Γ_k，令k＝1…N_μ，满足||Γ_k||＝1，且

其中，Γ_k为第k种微状态，Γ_k为N_μ个时间点所对应通道电位值的归一化；

3)令t＝1…N_t，

其中，N_t为GFP极大值个数，L_t表示每个GFP极大值所对应的微状态，V_t表示每个时刻各通道的电位值，V_t'为向量V_t的转置；

4)令k＝1…N_μ，在||X||＝1的约束条件下，计算计算

其中，S_k为协方差矩阵，X为协方差矩阵S_k的特征向量，为S_k的最大特征值所对应的特征向量；

5)计算

式中，σ₀ ²为噪声方差，N_s为电极数；

6)如果则Γ_k代表最终的微状态，如果算法回到步骤3)。

6.根据权利要求1所述的一种基于微状态的脑功能网络构建方法，其特征在于，所述S2中步骤6)将标记好的微状态时间序列具体处理方法为：判断相同状态所持续时间是否为80ms-120ms，若不满足，则去掉；若满足，则保留，从而得到最终的微状态时间序列。

7.根据权利要求1所述的一种基于微状态的脑功能网络构建方法，其特征在于：所述S3中步骤2)中PLV值的范围为[0，1]。

8.根据权利要求1所述的一种基于微状态的脑功能网络构建方法，其特征在于：所述S3中步骤4)阈值取稀疏度12％～40％，步长为2％。

9.根据权利要求1所述的一种基于微状态的脑功能网络构建方法，其特征在于：所述无权无向网络模型的规模为60*60。