CN109524112A

CN109524112A - 一种基于部分定向相干法的脑功能网络构建方法

Info

Publication number: CN109524112A
Application number: CN201811601572.XA
Authority: CN
Inventors: 张启忠; 胡月静; 石鹏; 袁长敏; 杨晨; 席旭刚; 罗志增
Original assignee: Hangzhou Dianzi University
Current assignee: Hangzhou Dianzi University
Priority date: 2018-12-26
Filing date: 2018-12-26
Publication date: 2019-03-26

Abstract

本发明提出了一种基于部分定向相干法的脑功能网络构建方法；本发明首先采集脑电信号，对脑电数据进行小波硬阈值去噪、用独立成分分析方法去除心电、眼电。其次，对脑电信号运用多轨道自回归模型进行建模，并通过拉普拉斯变换将其转化到频域上，得到线性部分定向相干值。然后，利用AIC准则确定模型的阶数，采用显著性水平的方法计算阈值，根据阈值确定二值化矩阵。最后，根据二值化矩阵建立脑功能网络连接图。功能网络连接图以采集脑电通道位置作为节点位置。数据显示使用该方法能够鉴别患者发病间期分别与人体视觉、体觉及精神功能相关脑网络结点的变化情况。

Description

一种基于部分定向相干法的脑功能网络构建方法

技术领域

本发明属于信号特征分析领域，涉及一种基于部分定向相干法的脑网络特征分析方法。

背景技术

癫痫是由于脑神经元超同步化异常放电引起的一种常见的脑部疾病，在神经系统疾病中排名第二，发病率仅次于脑卒中。癫痫发作是一种慢性长期反复性的疾病，具有突发性，在病发期间会造成脑功能暂时紊乱。在国外癫痫的患病率约为3‰～10‰，在我国的患病率有4‰～9‰，到目前为止，全球大约有5千万患者正在经历病痛的折磨，其中有超过40％的是儿童患者。一般情况下1岁以内的患病率最高，随着年龄的增长会逐渐降低。癫痫发作的临床表现可分为全面强直阵挛发作(大发作)、部分发作、失神发作(小发作)、癫痫持续状态等几种，常伴有不自主肌肉抽动，意识短暂性丧失症状，如果单次癫痫超过30分钟或频繁发作视为一种需要抢救的急症。这不仅给患者的身体、心理都带来了极大的痛苦，还给家庭带来了精神上的折磨及经济上的负担。而且在儿童癫痫的治疗方面，普遍存在误诊率高、复发率高和治疗无法降低复发率等问题。

头皮脑电信号(Electroencephalogram,EEG)是由大脑神经元突触后电位同步综合形成的生物电信号，反映了大脑运行状态和神经细胞活动情况。神经系统疾病往往会产生病变的EEG，如癫痫疾病会引起大脑部分区域出现异常性放电，使得EEG中掺杂了称为癫痫样波的异常信息。由于头皮脑电EEG具有无创性且蕴含了大量能反映患者生理、病理状态的生物电信息，因而，EEG成为了一种常用的疾病信号检测手段。从而使通过对EEG的分析研究各种神经疾病形成机理，对各器官功能的影响及寻求治愈方法成了各国专家学者的研究热点。如在癫痫疾病的EEG研究中，为探究神经元

异常放电的原因，HL Chan等提出了一种基于相空间特征相似性的尖峰分类方法，在每个尖峰时间点计算可以提供源活动流的PDC值，结果显示前端峰值具有高度集中的分布，但后面的峰值集中度有所下降。从而得出癫痫放电最初是由小面积的皮质神经元诱导出来的。但是遗憾的是此研究在只停留在理论实验上，还需要更多的临床研究来将映射模式与临床发现联系起来。王刚等基于PDC分析提出从方向和强度的角度研究大脑区域的信息流，分析癫痫发作前后因果关系的变化。方法适用于癫痫发作的检测与病灶定位。癫痫发病时间预测对减少患者的伤害有较大的意义。Oscrio等采用时频分析方法，将预测时间提高到15.5s，预测率达到92％。2017年最新研究成果，JesusMartinez-del-Rincon等人将非线性分类器应用于癫痫发作检测的脑电分析，该方法利用了脑电数据的潜在非线性特性，它的提出主要有两方面的贡献：一是通过内核技巧使用非线性分类器，二是提出 Bag-of-Words(BoW)模型来提取输入的非线性特征，改善的BoW后系统性能通过公开的数据库进行验证，实验显示平均F1-measure比现有的交叉数据集实验的方法改善了10％。

无论是从发病原因，病灶定位还是从癫痫发作预测，科研工作者都做了很好的研究及总结，但对于癫痫患者不同期脑网络特征与脑部分区功能之间关系分析的文献相对较少，有待更进一步深入。如果能够在研究儿童患者脑网络的特征上有所发现，则对预测、治疗、防范儿童癫痫病都有重要的意义。基于此，本文用部分定向相干法建立不同时期的脑功能网络连接图，分析发病间期与发病期的特征变化，在复杂网络静态特征的基础上研究大脑各功能分区的相关性变化趋势，从而达到分析癫痫儿童在发病期脑功能网络特征变化的目的。

发明内容

本发明针对现有技术的不足，提出了一种基于部分定向相干法的脑功能网络构建方法。

提出了基于部分定向相干法的脑功能网络构建方法。首先采集脑电信号，对脑电数据进行小波硬阈值去噪、用独立成分分析方法去除心电、眼电。其次，对脑电信号运用多轨道自回归模型(Multivariate autoregressive model，MVAR)进行建模，并通过拉普拉斯变换将其转化到频域上，得到线性部分定向相干值(Partial directional coherence，PDC)。然后，利用AIC准则确定模型的阶数，采用显著性水平的方法计算阈值，根据阈值确定二值化矩阵。最后，根据二值化矩阵建立脑功能网络连接图。功能网络连接图以采集脑电通道位置作为节点位置。数据显示使用该方法能够鉴别患者发病间期分别与人体视觉、体觉及精神功能相关脑网络结点的变化情况。标明了儿童患者在癫痫发作时与视觉、体觉、思维等能力相关的脑区脑功能的变化趋势，与其它方法相比，此方法揭示了脑功能网络在节点之间的信息流向及强度，更符合人类大脑的信息传递。此方法有望在研究儿童癫痫脑网络功能特征时提供一个新的视角。

为了实现以上目的，本发明方法主要包括以下步骤：

1.基于部分定向相干法的脑功能网络构建方法，其特征在于该方法包括如下步骤：

步骤(1)，采集脑电信号X_i(n)，对脑电数据进行小波硬阈值去噪、用独立成分分析方法去除心电、眼电。

步骤(2)，将步骤(1)中的信号运用多轨道自回归模型(Multivariateautoregressive model，MVAR)进行建模。

具体方法如下：

假设一个MVAR为：

式中，X(n)＝(x_i(n),x₂(n),...,x_d(n))^T，是一组离散的数据采样点，它表示的是在时刻n时采集到的数集。i，d表示通道数，P表示建立模型的阶数，由AIC 准则来确定。a_id为第i通道对第d通道影响状态系数，。E(n)为偏差矩阵，存放的是噪声信号。令：

式中p_m是MVAR的最大预测模型参数，上式可以写成：

Y＝DA

在公式两边同时左乘D^T即可求出系数矩阵A，这样就可以求得MVAR模型的系数：

A_d＝(D^TD)^-1D^TY

A是一个的MVAR模型d×d预测系数矩阵，表示预测时间序列，而E(n)也是一列长度为d的列向量，表示预测误差序列。

步骤(3)通过拉普拉斯变换将其转化到频域上，得到线性部分定向相干值(Partial directional coherence，PDC)。

具体方法如下：

在建立了MVAR之后，通过最小二乘法拟合了MVAR的系数A(r)，根据 M维单位矩阵和系数矩阵的傅里叶变换之间的差异，通过拉普拉斯变换将其转化到频域上：

式中p是MVAR的最大预测阶数，f是频率，那么从第j道到第i道的PDC值可以表示为：

步骤(4)利用AIC准则确定模型的阶数。AIC准则的计算公式为：

AIC＝2k-2ln(L)

其中，k是估计的参数个数，L是似然函数。公式右边的第一项描述了拟合的优劣性，当后边的估计得参数个数增加时，这一项的值会减少。而模型的复杂度就由第二项来决定，参数个数越小，模型就越优。

步骤(5)，采用显著性水平的方法计算阈值T，当PDC值大于T的时候，可以认为两个通道之间是具有相关性，用1表示；当PDC值小于T，则认为这两通道间不具有相关性用0表示，这样就得到二值化的矩阵。

步骤(6)，根据二值化的矩阵建立脑功能网络连接图。功能网络连接

图以采集脑电通道位置作为节点位置。

本发明与已有的诸多脑电信号的方法相比，具有如下特点：

由于脑电信号是非线性、非平稳的随机信号，此方法是基于图论的方法来研究患者的全局特征，在建立脑功能网络的基础上，提取度和聚类系数两个特征。特别是对基于PDC方法建立的有向图来说，顶点的出边条数称为该顶点的出度，顶点的入边条数称为该项点的入度，此方法揭示了信息传递的方向性。聚类系数表示图形中相同类节点聚集的系数。这样就整体刻画了所建立的脑功能网络的特性。

附图说明

图1为本发明的实施流程图；

图2为本发明实施例的EEG采集图；

图3为本发明实施例发病间期各频率段的因果关系图；

图4为本发明实施例在发病期各频率段的因果关系图；

图5为本发明实施例23通道在发病间期与发病期平均入度值的变化图；

图6为本发明实施例23通道在发病间期与发病期平均出度值的变化图。

图7为本发明实施例T8-P8区在发病间期与发病期各个频率段的出度值比较图。

具体实施方式

下面结合附图对本发明的实施例作详细说明：本实施例在以本发明技术方案为前提下进行实施，给出了详细的实施方式和具体的操作过程，但本发明的保护范围不限于下述的实施例。

如图1所示，本实施例包括如下步骤：

步骤(1)，采集脑电数据及预处理，本实例采用麻省理工学院提供的儿童癫痫数据(CHB-MIT)库和申请者与医院合作采集的癫痫信号进行对照实验，所有的信号均由国际标准的10-20电极分布系统采样得到。数据预处理包括小波去噪、心电、眼电剔除等。图2为本发明实施例的EEG采集图；

具体方法如下：

假设一个MVAR为：

式中p_m是MVAR的最大预测模型参数，上式可以写成：

Y＝DA

只要在公式两边同时左乘D^T即可求出系数矩阵A，由于D^TD是一个对称方,假如|D^TD|≠0，那么D^TD为可逆矩阵，这样就可以求得MVAR模型的系数：

A＝(D^TD)^-1D^TY

A是一个的MVAR模型d×d预测系数矩阵，表示预测时间序列，而E(n) 也是一列长度为d的列向量，表示预测误差序列。

具体方法如下：

图3为本发明实施例发病间期各频率段的因果关系图；图4为本发明实施例在发病期各频率段的因果关系图。

步骤(4)利用AIC准则确定模型的阶数。AIC准则的计算公式为：

AIC＝2k-2ln(L)

确定了模型阶数，接下来可以求取各通道间的相干系数。为描述方便，将各通道进行编号，同时列出了电极所放置位置，如表1所示。

表1通道间序号及对应脑电图极的位置

步骤(5)，采用显著性水平的方法计算阈值T，仍然以儿童A为例，其阈值为T＝0.1273，当PDC值大于0.1273的时候，可以认为两个通道之间是具有相关性，用1表示；当PDC值小于0.1273，则认为这两通道间不具有相关性用0表示，这样就得到二值化的矩阵。

步骤(6)，根据二值化的矩阵建立脑功能网络连接图。功能网络连接图以采集脑电通道位置作为节点位置。图3和图4分别为儿童A在发病间期及发病期各频率段的脑功能网络连接图。

分析显示，所建立的脑功能网络连接图符合复杂网络的定义：小世界性和无标度性，因而可用复杂网络的特征分析方法研究其特征。首先从度的角度来分析患者的脑功能网络连接特性的变化，度是指在因果关系图中节点υi连接的边数。特别是对基于PDC方法建立的有向图来说，顶点的出边条数称为该顶点的出度，顶点的入边条数称为该项点的入度，此方法揭示了信息传递的方向性。针对139次患病实例，分别求取发病间期和发病期复杂网络图每个节点的入度、出度平均值，将两个时期的结果求差，所得结果为正数表明在发病期此通道的出度值上升，若为负数说明在发病期此通道的出度值下降，相干性减弱。图5、图6给出的是各通道在发病间期与发病期平均入度变化情况和平均出度值的变化情况。

由数据可计算得入度变化值的方差为2.3004，出度变化值的方差为29.2648，相较于出度变化，入度值的变化不大，因而本文重点关注网络的出度。患者发病间期与发病期的出度的变化表现为：第2通道，第15通道，第22通道和第23通道的值增加，其它通道的出度值变化不明显或有明显减小趋势。第15通道和第23 通道由表1可知其测试电极都置于同一脑区，即T8-P8处。图7给出的是T8-P8 处23通道(节点)所有病例在发病间期与发病期各个频率段的出度值变化比较图。

在所记录的23个患病儿童，共计139次发病中，有101次癫痫发作时呈现(a) 图趋势，在癫痫病发中期大脑的T8-P8区的出度值明显上升，与其它通道的信息交互频率提高，相关性呈现上升状态，有32次的癫痫发作时呈现(c)图趋势，在发病期T8-P8区出度值无明显变化，但在各个频率段内出度值也有升有降，还有其余6 次病发状态呈现(b)图趋势，其出度值明显减小，与其它各通道的信息交互能力下降，相关性减小，但相对处于其它脑功能分区的节点其下降幅度较小。

进一步，在所建立的因果关系图中，从复杂网路的聚类系数特性探讨T患者发病间期与发病期的脑网络特征。聚类系数表示图形中相同类节点聚集的系数，在图形理论中，网络中的一个节点v_i有k_i条边将和它的其他节点相连，这k_i个节点就是节点v_i的邻居。显然在这k_i个节点之间最多只有k_i(k_i-1)/2条边，那么节点上实际存在的边数E_i与总的可能的边数的比值，就是这个节点的聚类系数C_i，也即为：

C_i＝2E_i/(k_i(k_i-1))

整个网络的聚类系数C就是所有节点的聚类系数C_i的平均值。很明显， 0≤C≤1，C＝0当且仅当所有的节点均为孤立节点，即没有任何连接边；C＝1当且仅当网络是全局耦合的，即网络中任意两个节点都直接相连。对于脑功能网络节点而言，节点的聚类系数值越大，说明该节点与其相邻节点越类似，该节点与其邻居节点协同完成某任务目标的能力越强。

表2各脑功能分区节点的聚类系数均值对照表

建立用于求取聚类系数阈值的患者脑功能网络连接图，阈值的选择依然参照上文建立邻连接矩阵的方法。表2为10个患病儿童发病间期与发病期求得的全部节点、各脑功能分区节点的聚类系数均值对照表。数据显示，额叶、枕叶及顶叶区的节点，其平均聚类系数原来处于较大值的仍保持较大值，其它原来较小的则明显增大。但处于大脑左、右他分界处的FZ-CZ、 CZ-PZ不同，其聚类系数是下降的。

Claims

1.一种基于部分定向相干法的脑功能网络构建方法，其特征在于，该方法包括如下步骤：

步骤(1)，采集脑电信号X_i(n)，对脑电数据进行小波硬阈值去噪、用独立成分分析方法去除心电、眼电；

步骤(2)，将步骤(1)中的脑电信号运用多轨道自回归模型MVAR进行建模；

具体算法如下：

假设一个MVAR为：

式中，X(n)＝(x_i(n),x₂(n),...,x_d(n))^T，是一组离散的数据采样点，它表示的是在时刻n时采集到的数集；i，d表示通道数，P表示建立模型的阶数，由AIC准则来确定；a_id为第i通道对第d通道影响状态系数；E(n)为偏差矩阵，E(n)^T＝(e_p+1 ^T,e_p+2 ^T,...,e_p+d ^T)是噪声信号；令：

式中p_m是MVAR的最大预测模型参数，表示最大模型参数下的系数矩阵，上式写成：

Y＝DA

在公式两边同时左乘D^T即可求出系数矩阵A，求得MVAR模型的系数：

A_d＝(D^TD)^-1D^TY

A是一个的MVAR模型d×d预测系数矩阵，表示预测时间序列，而E(n)是列长度为d的偏差矩阵，表示预测误差序列，X_i(n-i)表示时刻n时第i通道的系数；

步骤(3)通过拉普拉斯变换将其转化到频域上，得到线性部分定向相干值PDC；

具体算法如下：

在建立了MVAR之后，通过最小二乘法拟合了MVAR的系数A(r)，根据M维单位矩阵和系数矩阵的傅里叶变换之间的差异，通过拉普拉斯变换将其转化到频域上：

式中p是MVAR的最大预测阶数，f是频率，那么从第j道到第i道的PDC值π_ij(f)表示为：

其中，H表示共轭转置，A_ij(f)是傅里叶变换后的系数矩阵A(r)对应的第i行第j行的系数；a_j(f)是j列的系数

步骤(4)利用AIC准则确定模型的阶数；AIC准则的计算公式为：

AIC＝2k-2ln(L)

其中，k是估计的参数个数，L是似然函数；公式右边的第一项描述了拟合的优劣性，当估计的参数个数增加时，这一项的值会减少；而模型的复杂度就由第二项来决定，参数个数越小，模型就越优；

步骤(5)，采用显著性水平的方法计算阈值T，当线性部分定向相干值PDC值大于T的时候，认为两个通道之间是具有相关性，用1表示；当PDC值小于T，则认为这两通道间不具有相关性用0表示，这样就得到二值化的矩阵；

步骤(6)，根据二值化矩阵建立脑功能网络连接图；功能网络连接图以采集脑电通道位置作为节点位置。