一种基于多通道皮层脑电记录的大尺度皮层网络信息流个体化三维动态可视化方法

技术领域

本发明涉及神经工程领域，尤其是涉及一种基于多通道皮层脑电记录的大尺度皮层网 络信息流个体化三维动态可视化方法。

背景技术

人脑是一个开放的复杂巨系统，生理学和解剖学的研究表明，大脑皮层内神经元的连 接具有局部性和全局性。这种大脑皮层连接的结构，既保证了局部区域内的信息的高效交 互，又确保了脑区间的信息传输。大脑通过位于不同空间位置脑区的信息整合来完成认知 功能，即使对于单一的认知任务，也同时需要激活多个空间分离的脑功能区域。先前的很 多实验和理论研究一致表明：大脑在完成认知功能加工过程中既要实现自下而上的信息加 工，又要完成自上而下的调节与控制功能。脑区的神经元群在不同的认知任务下动态的相 互作用并且形成了不同的神经环路。先前的研究主要集中在定位认知功能特异脑区，强调 用大脑区域来定位认知功能，孤立的关注脑区与特定认知功能的联系；近年来随着观测技 术和研究手段的不断发展，越来越多的研究者逐渐着眼于大脑不同脑区在完成某一认知功 能时的相互作用，从脑网络的层次，从信息流传导的层次，将大脑作为一个整体进行功能 整合来对大脑进行深入、细致、全面的研究。

虽然大脑不同脑区对应于不同的认知功能，但要完成一项认知功能，往往需要大脑多 个不同的脑区相互协调、交互作用，共同动态的构成一个复杂的网络来发挥其功能，即大 脑的功能执行总是依赖于多个脑区之间广泛地交互。因此，从脑网络的角度来研究大脑不 同脑区的神经生理动力学特性是极为必要的。神经生理学实验在大脑大尺度上得到了大量 皮层脑电（ECoG）数据，这种ECoG数据具有极高的时空分辨率，为我们研究脑功能提 供了新的窗口，但是对于大脑网络不同模块间信息交互的度量和量化还是很浅显的，要求 发展新的计算和可视化方法来分析和研究大脑网络不同模块间的信息交互。

发明内容

本发明的目的在于提供一种基于多通道皮层脑电记录的大尺度皮层网络信息流个体化 三维动态可视化方法，以便从大尺度脑网络信息传导的层次研究不同脑区间的信息交互， 为大脑网络不同模块间信息交互的度量和量化提供新的计算和可视化方法。

为实现上述目的，本发明提出一种基于多通道皮层脑电记录的大尺度皮层网络信息流 个体化三维动态可视化方法，包括以下步骤：

步骤S1，对采集到的皮层脑电（ECoG）数据进行预处理，以选取感兴趣通道；

步骤S2，依据所选取的感兴趣ECoG数据，计算在不同时间和频率的大脑皮层网络信 息流矩阵；

步骤S3，对经过显著性检验的信息输入、输出流在三维立体脑结构上进行三维动态可 视化。

所述步骤S2具体包括以下步骤：

步骤2.1：利用基于多变量自回归模型的自适应有向传递函数ADTF计算对所述感兴 趣ECoG数据在不同时频点的传递矩阵；

多变量自回归模型的描述如下：

$X\left(t\right)=\underset{i=1}{\overset{p}{\Sigma }}A(\tau ,t)X(t-\tau )+E\left(t\right)---\left(1\right)$

其中，X(t)是信号矩阵，A(τ，t)是延迟时间τ的时变系数矩阵，E(t)是和信号X(t)独立的 白噪声，p是可以使用施瓦茨贝叶斯准则（SBC）来确定的模型的阶；为了在频域研究信号 之间的信息交互，使用傅立叶变换将（1）式方程变换到频域：

A(t，f)X(f)=E(f)   （2）

其中，

$A(t,f)=-\underset{\tau =0}{\overset{p}{\Sigma }}{A}_{\tau }\left(t\right)e^{-j2\mathrm{\pi f\tau }}---\left(3\right)$

令A₀(t)=-I，方程（2）改写为：

X(f)=A^-1(t，f)E(f)=H(t，f)E(f)   （4）

其中，初始值A₀(t)为负的单位矩阵，H(t，f)是随时间变化的转移矩阵，H_ij(t，f)表示在 时间t和频率f从信号x_j到信号x_i的信息流大小；则归一化的ADTF值可由下式计算：

${\mathrm{ADTF}}_{\mathrm{ij}}(t,f)=\frac{{\left|H_{\mathrm{ij}}(t,f)\right|}^2}{{\Sigma }_{k=1}^n{\left|H_{\mathrm{ik}}(t,f)\right|}^2}---\left(5\right)$

其中，n是脑电通道的总数，一般脑电系统为64通道，i和j表示通道编号，则ADTF_ij(t，f) 就表示了在时间t和频率f时由通道j流向通道i的信息流大小；

步骤2.2：利用偏相关函数计算对所述感兴趣ECoG数据在不同时频点的偏相关系数；

首先计算所述感兴趣ECoG数据的交叉谱C(t，f)，计算公式如下：

C(t，f)=H(t，f)∑H^*-1(t，f)   （6）

其中，∑表示白噪声矩阵E(t)的方差，上标*表示对所述转移矩阵H(t，f)取共轭，上标-1 表示对所述转移矩阵H(t，f)进行转置；

则偏相关系数的计算公式为：

${\Gamma }_{\mathrm{ij}}(t,f)=\frac{M_{\mathrm{ij}}^2(t,f)}{M_{\mathrm{ii}}(t,f)M_{\mathrm{jj}}(t,f)}---\left(7\right)$

其中，Γ_ij(t，f)表示在时间t和频率f时信号x_j和信号x_i的偏相关系数，M_ij(t，f)表示所述 交叉谱C(t，f)的代数余子式；

步骤2.3：依据自适应直接传递函数和偏相关系数计算自适应有向直接传递函数 AdDTF，得到在不同时间、频率点感兴趣脑区间的直接信息流强度和方向；

所述自适应有向直接传递函数AdDTF的计算公式如下：

AdDTF(t，f)=ADTF(t，f)Γ_ij(t，f)   （8）

其中，AdDTF_ij(t，f)值表示了在时间t和频率f时信号x_j和信号x_i间的直接信息流强度和 方向，反映了大尺度皮层网络信息的交互与整合；

因此，流入或流出某一通道的归一化直接信息流计算方法为：

${\mathrm{flow}}_{\mathrm{ij}}=\frac{{\Sigma }_{t=t_1}^{t_2}{\Sigma }_{f=f_1}^{f_2}{\mathrm{AdDTF}}_{\mathrm{ij}}(t,f)}{(t_2--t_1)(f_2-f_1)}---\left(9\right)$

其中，t₁和t₂分别对应所选取数据的起始时刻和终止时刻，f₁和f₂分别对应所选取数据 的起始频率和终止频率；

流出某一通道的直接信息输出流为：

$\mathrm{out}{\mathrm{flow}}_j=\frac{{\Sigma }_{i=1}^n{\mathrm{flow}}_{\mathrm{ij}}}{n-1},\mathrm{fori}\ne j---\left(10\right)$

流入某一通道的直接信息输入流为：

${\mathrm{inflow}}_i=\frac{{\Sigma }_{j=1}^n{\mathrm{flow}}_{\mathrm{ij}}}{n-1},\mathrm{forj}\ne i---\left(11\right)$

其中，公式（10）和（11）中n为通道总数，脑电系统为64通道，i和j表示通道编号；

步骤2.4：信息流统计显著性检验。

进一步，所述步骤S1具体包括以下步骤：

步骤1.1：对所述ECoG数据进行0.3-150Hz滤波和50Hz陷波处理，以消除低频飘移 和工频干扰；

步骤1.2：对滤波处理后的ECoG数据进行坏通道检验处理，去除坏通道，避免不可 靠数据的影响；

步骤1.3：对去除坏通道后的ECoG数据进行转换参考处理，采用平均参考；

步骤1.4：对转换参考后的ECoG数据进行基线校正处理，去除基线漂移，以使ECoG 数据趋于平稳；

步骤1.5：对基线校正后的ECoG数据进行伪迹检测处理，去除头、面部肌肉抖动干 扰的数据片段；

步骤1.6：从预处理后的ECoG数据中选取记录电极落在感兴趣脑区的ECoG数据。

本发明所述的基于多通道皮层脑电记录的大尺度皮层网络信息流个体化三维动态可视 化方法，通过ECoG数据预处理、选取感兴趣通道、利用自适应传递函数及偏相关系数计 算大尺度皮层网络信息流矩阵，再用替代数据的方法检验皮层网络信息流交互矩阵的统计 显著性，消除了有随机性引起的无意义信息传导，提高了信息流计算的精确性和鲁棒性。 最后利用箭头编码信息流向，利用箭头的粗细和颜色编码信息流强度，在个体化的三维立 体脑结构上实现动态信息流可视化。所述方法从大尺度脑网络信息传导的层次研究不同脑 区间的信息交互，为大脑皮层网络不同模块间信息交互与整合的度量和量化提供了新的计 算和可视化方法，是皮层脑电数据分析方法的有益补充，对于深入探索脑功能机制有着重 要启示。

附图说明

图1为本发明多通道皮层脑电记录的大尺度皮层网络信息流个体化三维动态可视化方 法流程图；

图2为根据本发明一实施例在特定时间和频率的大脑皮层网络信息流矩阵示意图；

图3为一实施例使用本方法得到的大尺度皮层网络信息流个体化3D动态可视化结果。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附 图，对本发明进一步详细说明。

参照图1，本发明所述的一种基于多通道皮层脑电记录的大尺度皮层网络信息流个体 化三维动态可视化方法，具体实施步骤如下：

步骤S1：对采集到的皮层脑电（ECoG）数据进行预处理，以选取感兴趣通道。

所述步骤S1具体包括以下步骤：

步骤1.1：对所述ECoG数据进行0.3-150Hz滤波和50Hz陷波处理，以消除低频飘移 和工频干扰。

步骤1.2：对滤波处理后的ECoG数据进行坏通道检验处理，去除坏通道，避免不可靠 数据的影响。

步骤1.3：对去除坏通道后的ECoG数据进行转换参考处理，采用平均参考。

步骤1.4：对转换参考后的ECoG数据进行基线校正处理，去除基线漂移，以使ECoG 数据趋于平稳。

步骤1.5：对基线校正后的ECoG数据进行伪迹检测处理，去除头、面部肌肉抖动干扰 的数据片段。

步骤1.6：从预处理后的ECoG数据中选取记录电极落在感兴趣脑区的ECoG数据。所 述感兴趣脑区的选取视不同的认知任务加工不同而不同，比如研究语言认知加工过程的 ECoG数据，则感兴趣脑区一般选取Broca区、Wernicke区及运动区；如果研究学习和记 忆，则对应的感兴趣脑区一般为大脑皮层联络区、海马及其邻近结构、杏仁核、丘脑和脑 干网状结构。

步骤S2:依据所选取的感兴趣ECoG数据，计算在不同时间和频率的大脑皮层网络信 息流矩阵。

所述步骤S2具体包括以下步骤：

步骤2.1：利用基于多变量自回归模型的自适应有向传递函数（ADTF）计算对所述感 兴趣ECoG数据在不同时频点上的传递矩阵。

自适应有向传递函数是一个随时间变化的多变量自回归方法，它对于大尺度皮层网络 的非稳态信号间的信息交互的估计非常有用。多变量自回归模型的描述如下：

$X\left(t\right)=\underset{i=1}{\overset{p}{\Sigma }}A(\tau ,t)X(t-\tau )+E\left(t\right)---\left(1\right)$

其中，X(t)是信号矩阵，A(τ，t)是延迟时间τ的时变系数矩阵，E(t)是和信号X(t)独立的 白噪声，p是可以使用施瓦茨贝叶斯准则（SBC）来确定的模型的阶。为了在频域研究信号 之间的信息交互，使用傅立叶变换将（1）式方程变换到频域：

A(t，f)X(f)=E(f)   （2）

其中，

$A(t,f)=-\underset{\tau =0}{\overset{p}{\Sigma }}{A}_{\tau }\left(t\right)e^{-j2\mathrm{\pi f\tau }}---\left(3\right)$

令A₀(t)=-I，方程（2）可改写为：

X(f)=A^-1(t，f)E(f)=H(t，f)E(f)   （4）

其中，初始值A₀(t)为负的单位矩阵，H(t，f)是随时间变化的转移矩阵，H_ij(t，f)表示在 时间t和频率f从信号x_j到信号x_i的信息流大小。则归一化的ADTF值可由下式计算：

${\mathrm{ADTF}}_{\mathrm{ij}}(t,f)=\frac{{\left|H_{\mathrm{ij}}(t,f)\right|}^2}{{\Sigma }_{k=1}^n{\left|H_{\mathrm{ik}}(t,f)\right|}^2}---\left(5\right)$

其中，n是脑电通道的总数，一般脑电系统为64通道，i和j表示通道编号，则ADTF_ij(t，f) 就表示了在时间t和频率f时由通道j流向通道i的信息流大小。

步骤2.2：利用偏相关函数计算对所述感兴趣ECoG数据在不同时频点的偏相关系数。

首先计算所述感兴趣ECoG数据的交叉谱C(t，f)，计算公式如下：

C(t，f)=H(t，f)∑H^*-1(t，f)   （6）

其中，∑表示白噪声矩阵E(t)的方差，上标*表示对所述转移矩阵H(t，f)取共轭，上标-1 表示对所述转移矩阵H(t，f)进行转置。

则偏相关系数的计算公式为：

${\Gamma }_{\mathrm{ij}}(t,f)=\frac{M_{\mathrm{ij}}^2(t,f)}{M_{\mathrm{ii}}(t,f)M_{\mathrm{jj}}(t,f)}---\left(7\right)$

其中，Γ_ij(t，f)表示在时间t和频率f时信号x_j和信号x_i的偏相关系数，M_ij(t，f)表示所述 交叉谱C(t，f)的代数余子式。

步骤2.3：依据自适应直接传递函数和偏相关系数计算自适应有向直接传递函数 （AdDTF），得到在不同时间、频率点感兴趣脑区间的直接信息流强度和方向。

所述自适应有向直接传递函数（AdDTF）的计算公式如下：

AdDTF(t，f)=ADTF(t，f)Γ_ij(t，f)   （8）

其中，AdDTF_ij(t，f)值表示了在时间t和频率f时信号x_j和信号x_i间的直接信息流强度和 方向，反映了大尺度皮层网络信息的交互与整合。

因此，流入或流出某一通道的归一化直接信息流计算方法为：

${\mathrm{flow}}_{\mathrm{ij}}=\frac{{\Sigma }_{t=t_1}^{t_2}{\Sigma }_{f=f_1}^{f_2}{\mathrm{AdDTF}}_{\mathrm{ij}}(t,f)}{(t_2--t_1)(f_2-f_1)}---\left(9\right)$

其中，t₁和t₂分别对应所选取数据的起始时刻和终止时刻，f₁和f₂分别对应所选取数据 的起始频率和终止频率。

流出某一通道的直接信息输出流为：

$\mathrm{out}{\mathrm{flow}}_j=\frac{{\Sigma }_{i=1}^n{\mathrm{flow}}_{\mathrm{ij}}}{n-1},\mathrm{fori}\ne j---\left(10\right)$

流入某一通道的直接信息输入流为：

${\mathrm{inflow}}_i=\frac{{\Sigma }_{j=1}^n{\mathrm{flow}}_{\mathrm{ij}}}{n-1},\mathrm{forj}\ne i---\left(11\right)$

其中，公式（10）和（11）中n为通道总数，一般脑电系统为64通道，i和j表示通道 编号。

步骤2.4：信息流统计显著性检验。

为了检验采用自适应有向直接传递函数（AdDTF）计算得到的信息流的统计显著性， 使用替代数据方法进行统计检验。采用一种改进的AAFT迭代生成算法，即IAAFT算法 生成替代数据。IAAFT算法生成的替代数据以其与原始数据的概率密度分布和傅里叶频谱 极好的匹配度，在数据的非线性统计检验中被广泛采用。依据IAAFT算法重复2000次产 生一经验分布，在设定的显著性水平P=0.05下进行显著性检验，以去除随机性引起的无意 义信息流。

本实施例方法共选取11个感兴趣通道，利用所述步骤S2计算得到大脑皮层网络信息 流矩阵如图2所示。

步骤S3，对经过显著性检验的信息输入、输出流在三维立体脑结构上进行三维动态 可视化。

所述步骤S3具体包括以下步骤：

步骤3.1：综合利用多模态医学影像数据，对颅内电极进行三维立体定位，以得到电极 三维坐标。

首先通过X线平片建立颅内电极与脑组织结构的联系，然后对扫描脑部获取的MRI 影像进行分割与重建得到单一被试个体的真实脑模型，再把提取得到的电极坐标标定到重 建得到的真实脑模型上，完成颅内电极的个体化三维立体定位。

步骤3.2：根据重建得到的颅内电极三维坐标，进行大尺度皮层网络信息流交互模式的 三维动态可视化。

利用从一个通道指向另一个通道的箭头来编码信息流向，利用箭头的颜色和大小来编 码信息流的强度。取经过显著性检验的信息流交互值，即可在个体化的三维立体脑结构上 描绘出信息流的交互模式图，进而可以呈现大尺度皮层脑电数据ECoG在不同时间、频率 时的信息流交互模式动态图，如图3a所示。

步骤3.3：根据重建得到的颅内电极三维坐标，进行大尺度皮层网络信息输出、输入流 模式的三维动态可视化。

首先，为了获得流经一通道的信息输出、输入流，分别采用所述步骤2.3中公式（10） 和（11）来计算其输出、输入流。然后利用球体的颜色和大小来编码信息输出、输入流的 大小，取经过显著性检验的信息输出、输入流值，即可在个体化的三维立体脑结构上描绘 出信息输出、输入流的模式图，进而可以呈现大尺度皮层脑电数据ECoG在不同时间、频 率时的信息输出、输入流模式动态图，如图3b和图3c所示。

从以上基于真实皮层脑电数据ECoG的大尺度皮层网络信息流个体化三维动态可视 化结果可以看出，本发明所述的基于多通道皮层脑电记录的大尺度皮层网络信息流个体化 三维动态可视化方法，通过ECoG数据预处理、选取感兴趣通道、利用自适应传递函数及 偏相关系数计算大尺度皮层网络信息流矩阵，再用替代数据的方法检验皮层网络信息传导 矩阵的统计显著性，最后利用可视化技术，实现了大尺度皮层脑网络信息流个体化的三维 动态可视化。所述方法从大尺度脑网络信息传导的层次研究不同脑区间的信息交互，为大 脑皮层网络不同模块间信息交互与整合的度量和量化提供了新的计算和可视化方法，是皮 层脑电数据分析方法的有益补充，对于深入探索脑功能机制有着重要启示。

从事本领域研究工作的工程技术人员应当意识到，以上所述的实施例是为了帮助读者 理解本发明的原理，应被理解为发明的保护范围并不局限于这样的特别陈述和实施例。凡 是根据上述描述做出各种可能得等同替换或改变，均被认为属于本发明的权利要求的保护 范围。