CN107126193B - 基于滞后阶数自适应选择的多变量因果关系分析方法 - Google Patents

Abstract

本发明涉及一种基于滞后阶数自适应选择的多变量因果关系分析方法，旨在为了准确地预测基于运动想象脑电信号的各个脑区域之间的因果影响。目前基于自回归模型的传统Granger因果关系方法缺乏考虑多变量时间序列中存在的时滞依赖结构以及模型系数对因果性的影响。本项目首先获取多通道运动想象脑电信号，其次采用改进的后向时间选择算法估计回归模型中每个变量的最优滞后阶数，建立多通道脑电信号的动态回归模型，然后利用模型的残差和系数来定义多变量之间的条件因果测度，可以有效地提高真实因果关系的估计性能。该方法在因效性脑功能网络、皮层肌肉耦合分析领域具有广阔的应用前景。

Description

基于滞后阶数自适应选择的多变量因果关系分析方法

技术领域

本发明属于脑功能网络分析领域，涉及动态回归模型中每个变量滞后阶数的自适应选择和多变量因果关系分析方法。

背景技术

人脑作为世界上中最为复杂的动力学系统之一，其皮层由150-330亿个神经细胞组成，这些神经元间的相互连接构成了一个拥有相当复杂结构和功能的脑网络，在神经学和生理学方面，大脑的主要功能是控制和支配身体的各个器官，因此，这个复杂而庞大的脑网络使它具有高级的信息处理和认知表达功能，如：语言、情感、记忆、认知等，并对来自人体内部和周围环境的信息进行存储、处理、加工和整合，探究大脑的结构和功能，加速脑科学领域的研究不仅可以提高对大脑疾病的预防、诊断和治疗，更可以推动人工智能的发展。

近年来，已有越来越多的学者注意到研究脑网络结构对分析大脑各个区域之间功能连接的重要性。以Stam为代表的一些研究人员最先将“小世界”网络理论应用于分析疾病产生原因，这种“小世界”属性可以使大脑完成某种功能时消耗尽可能少的资源，并且使得大脑在完成复杂功能时需要的网络连接较少，实现大脑的最优连接。脑网络在疾病诊断和治疗方面具有很好的前景，例如在阿尔茨海默氏症(Alzheimer’s Disease)抑郁症、癫痫、精神分裂等疾病的早期检测和诊断中有一定的研究进展，这些研究基于功能磁共振成像、脑电图(EEG)、脑磁图、结构磁共振成像等多种成像方法。Alexander等人系统地探讨了脑区之间的功能连通性、复杂的网络拓扑和结构上(欧式)的距离连接，发现精神分裂症患者功能网络组织的空间和拓扑紊乱可能产生于过度“修剪”的短距离功能连接。Supekar采用区域脑网络构建方法分别对阿尔茨海默氏症患者和健康受试者的功能磁共振成像数据进行研究，发现在脑网络聚类系数上阿尔茨海默氏症患者明显低于健康受试者，证明了该患者脑网络的局部效率较低。这些研究表明，基于复杂网络理论的脑功能网络在探索大脑信息处理和传递机制方面具有巨大的优势和潜力。

脑网络可分为结构性网络(Structure Network)、功能性网络(FunctionalNetwork)、效用性网络(Effective Network)。结构性网络可通过神经元之间的解剖性连接或成像技术来确定，反映了大脑的生理结构。功能性网络是描述网络节点之间存在的统计性连接关系，可采用互相关、互信息量等量化方法构建，属于无向网络。效应性网络相较于反映节点间关联性的功能性连接，不仅反映统计意义上的连通性，而且能够判断信息在节点之间的传播方向。已有许多学者进行了有效的连通性建模方法的研究，例如应用于协方差结构方程模型的神经成像数据、非线性系统识别技术以及确定性状态空间模型的贝叶斯估计。但是，这些有效连接方法都需要预先选择相关区域，假设关联性影响方向。这些预先指定的模型可用于假设在进行认知任务时，不同脑区之间的关联性。然而，该类方法隐含的问题是不适合的模型可能导致错误的结论。近年来人们将推测变量之间因果信息流向的格兰杰因果关系(Granger causality,GC)方法应用于研究脑网络的效应性连接，该方法使得时间序列优先从数据中的信息识别关联性方向，提供了在行为和认知任务过程中，大脑区域间的方向性相互作用与信息传递机制。但是建立在线性回归模型上的传统GC分析方法中，线性回归模型中所有变量具有相同的滞后阶数，没有考虑到多变量时间序列中通常存在的时滞依赖结构，限制了模型的估计性能。

本发明针对传统GC方法存在的该问题，提出了基于动态回归模型滞后阶数自适应选择的多变量因果关系方法，用来分析多导联运动想象EEG信号中以C3、C4和Cz为典型代表的不同脑区之间的因果影响。

发明内容

本发明的目的就是针对传统GC方法没有考虑多变量时间序列中存在的时滞依赖结构以及模型系数对因果关联性的影响，提出了一种基于动态回归模型滞后阶数自适应选择的多变量因果关系方法。

脑功能网络是个复杂且稀疏的抽象网络，其构建首先要定义网络节点。对于多通道EEG信号，往往把每个EEG导联(通道)对应的电极所测量的区域定义为一个节点，其电活动为若干时间序列，然后计算这些时间序列之间的相关系数，各节点之间相关系数的大小反映出对应脑区之间的功能连接强度，采用改进的后向时间选择(modified backward-in-time-selection,mBTS)算法动态选取回归模型中每个变量的滞后阶数，然后利用模型的残差和系数重新定义变量之间的因果关系。

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

步骤(1)获取多通道运动想象脑电信号，具体是：采用多导联电极帽采集不同运动想象实验范式下的脑电信号；

步骤(2)采用改进的后向时间选择算法建立多通道脑电信号的动态回归模型；

具体是：给定P个通道、长度为N的脑电信号时间序列{X_j,t:j＝1,2,...,P；t＝1,2,…,N}，首先采用改进的后向时间选择(mBTS)算法估计出每个变量X_j的最优滞后阶数，然后考虑所有变量的历史信息对X_j的影响，X_j的动态回归模型可表示为

其中，u_j,t为模型的回归估计残差，a_jp,l是模型系数，p＝1,2,...,P，l＝1,2,...,k_p；k_p是由mBTS算法估计出的变量X_p的最优滞后阶数；

步骤(3)以步骤(2)得到的动态回归模型为基础，采用模型的残差和系数构建条件因果测度，描述多变量之间因果关系的强度；

具体是：建立X_i→X_j的条件因果测度为

其中，k_p为变量X_p的最优滞后阶数，m＝max(k_p|p＝1,2,…,P)，p＝1,2,...,P，i∈{1,2,...,P}且i≠j；同理，X_j→X_i的条件因果测度表示为

步骤(4)采用步骤(3)的条件因果测度分析运动想象任务时不同区域脑电信号之间的因果关系。

本发明与已有的多种GC分析方法相比，具有如下特点：

1、采用mBTS算法估计回归模型中每个变量的最优滞后阶数

传统的GC分析方法建立在回归(Autoregressive,AR)模型之上，该模型中所有变量具有相同的阶数，这种向量形式极大地限制了模型的估计，而且在模型的预测中矢量时间序列的最佳拟合对于分量时间序列而言不是最优的。此外，由于所有参数同时被估计，当它们的数量相对于可用数据的数量较大时，估计变得数值不稳定。本发明将AR模型的矢量形式分解成标量方程，构建动态回归(Dynamic regression,DR)模型，采用mBTS算法估计DR模型中每个变量的最优滞后阶数，避免了每个分量具有相同阶数的限制。

2、采用模型系数和残差项构建新的条件因果关系测度，定量描述多变量之间的相关性

在基于AR模型的GC方法中，定义的因果关系测度只包括残差项。然而，实际上影响变量间因果关系的因素包括模型的系数和残差项。本发明采用新的因果关系测度定义，不但包括DR模型的残差和系数对预测因果关系准确性的影响，而且在考虑两两时间序列之间的因果关系时加入了第三时间序列的影响，在一定程度上抑制了间接因果关系对直接因果关系的影响，有效提高了对真实因果关系的估计性能。

附图说明

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

具体实施方式

下面结合附图详细描述本发明基于DR模型滞后阶数自适应选择的多变量因果关系方法，图1为实施流程图。

如图1，本发明方法的实施主要包括六个步骤：(1)获取多通道运动想象脑电信号样本数据；(2)建立多通道脑电信号的DR模型，采用改进的后向时间选择算法估计回归模型中每个变量的最优滞后阶数；(3)以步骤(2)得到的DR模型为基础，采用模型的残差和系数构建条件因果测度，描述多变量之间因果关系的强度；(4)采用步骤(3)得到的条件因果测度分析运动想象任务时不同区域脑电信号之间的因果关系。

下面逐一对各步骤进行详细说明。

步骤一：获取多通道运动想象脑电信号样本数据

采用美国NeuroScan公司的40导脑电帽采集运动想象脑电信号。受试者按要求佩戴好脑电帽后坐在轮椅上，保持安静、自然，注视实验环境中设定的情景提示。实验范式如下：想象左手操控轮椅控制杆向左、想象右手操控轮椅控制杆向右，分别对应轮椅左转、右转的控制运动形式，在实施过程中还可根据实验的具体情况对实验模式的设计做适当修正。

步骤二：建立多通道运动想象脑电信号的动态回归(DR)模型

给定P个通道、长度为N的脑电信号时间序列{X_j,t:j＝1,2,...,P；t＝1,2,…,N}，首先采用改进的后向时间选择(mBTS)算法估计出每个变量X_j的最优滞后阶数，然后考虑所有变量的历史信息对X_j的影响，X_j的DR模型可表示为：

其中，u_j,t为模型的回归估计残差，a_jp,l是模型系数，p＝1,2,...,P，l＝1,2,...,k_p；k_p是由mBTS算法估计出的变量X_p的最优滞后阶数。而传统的AR模型表示如下：

其中，m是AR模型的阶数。可以看出，AR模型中所有变量具有相同的滞后阶数，而基于滞后阶数自适应选择的动态回归模型中，不同变量具有各自的最优滞后阶数。

步骤三：以基于mBTS的动态回归模型为基础，采用模型的残差和系数构建条件因果关系测度

在步骤二得到的DR模型中，哪一项在所有项中占的比例越大，则该项对X_j有较大的因果关联性，从而建立X_i→X_j的条件因果测度为

其中，k_p为变量X_p的最优滞后阶数，m＝max(k_p|p＝1,2,…,P)，p＝1,2,...,P，i∈{1,2,...,P}且i≠j。而在传统的格兰杰因果关系方法中，除了公式(3)所示的全模型之外，还需要根据不考虑变量X_i的历史数据对X_j的影响来建立X_j的简模型，表示如下

进而以全模型和简模型为基础，X_i→X_j的条件格兰杰因果测度被定义为

比较公式(3)和公式(5)，可以看出:传统的条件格兰杰因果测度只涉及模型的残差项，而本发明建立的条件因果测度同时考虑了模型的残差和系数对预测因果关系准确性的影响。

根据公式(3)的定义，同理可得X_j→X_i的条件因果测度，可表示为

步骤四：采用步骤三得到的条件因果测度分析运动想象任务时不同区域脑电信号之间的因果作用及影响程度。

本实施方案采用Graz科技大学提供的脑-机接口国际竞赛数据集Datasets2b。该数据集包含9名受试者，共采集了C3、C4、Cz三个通道的EEG数据，包括左手和右手两类运动想象数据。每个受试者采集5组，前三组为训练数据集，后两组是测试数据集。其中，前两组训练数据集是无视觉反馈的数据，各包含120次运动想象数据，而第三组训练数据集是有视觉反馈的数据，包含160次运动想象数据。考虑到测试数据的分类结果未知，因此本文选取前三组训练数据集进行分析，共有27组实验数据。对于每个受试者的每组实验，将同类任务中C3、Cz、C4三个通道的运动想象数据分别进行叠加和求平均运算，分别得到每个受试者在两类任务中各个通道的平均数据，用于因果关系分析。

采用文献“Hu S,Wang H,Zhang J,et al.Comparison analysis:grangercausality and new causality and their applications to motor imagery.IEEETransactions on Neural Networks and Learning Systems,2016,27(7):1429-1444”提出的指标来评价不同因果分析方法的性能，具体的实验步骤如下：首先，采用AIC算法确定回归模型阶次最合适的取值，得出最佳阶次m＝12。然后，分别运用基于AR模型的条件格兰杰因果方法(CGC)、新型条件因果方法(NCC)、基于DR模型的格兰杰因果方法(CGCI)、本发明方法(NCCI)四种因果分析方法，计算每次实验左右手想象运动中C3→C4、C4→C3、C3→Cz、Cz→C3、C4→Cz、C4→Cz的因果值，其次分别统计出Cz→C3的因果值大于的因果值的次数、Cz→C4的因果值大于C4→Cz的因果值的次数、C4→C3的因果值大于C3→C4的因果值的次数，然后统计出这些次数在27次实验中所占的比值(分别记为Cz→C3＞C3→Cz、Cz→C4＞C4→Cz、C4→C3＞C3→C4)，作为因果分析方法预测变量间因果性能的评判标准。该比值越大，表示该方法的因果分析性能越好。实验结果如表1和表2所示。

表1数据集Datasets 2b中多次左手运动想象所得Cz→C3＞C3→Cz、Cz→C4＞C4→Cz和C4→C3＞C3→C4的比列

表2数据集Datasets 2b中多次右手运动想象所得Cz→C3＞C3→Cz、Cz→C4＞C4→Cz和C3→C4＞C4→C3的比列(％)

在表1和2中，对CGC、NCC、CGCI和NCCI方法进行了比较和分析，可以看出：

1)想象左、右手运动时，四种方法都表现出同样的趋势，即Cz对C3的因果关系强于C3对Cz的因果关系，Cz对C4的因果关系强于C4对Cz的因果关系。想象左手运动时，C3对C4的因果关系弱于C4对C3的因果关系，而想象右手运动时，C3对C4的因果关系强于C3对C4的因果关系，此结果与文献“Hu S,Wang H,Zhang J,et al.Comparison analysis:grangercausality and new causality and their applications to motor imagery.IEEETransactions on Neural Networks and Learning Systems,2016,27(7):1429-1444”一致。

2)想象左手运动时，NCCI得出的Cz→C3＞C3→Cz比例(93％)比NCC高4％，比CGCI高12％，比CGC高23％；NCCI得出的Cz→C4＞C4→Cz比例(89％)分别比CGCI高8％，比NCC高15％，比CGC高26％；NCCI得出的C4→C3＞C3→C4的比例(74％)分别比CGCI、NCC、CGC高出11％、15％、22％。可见，在想象左手运动时，本发明提出的NCCI方法更能清楚地确定因果关系影响。

3)想象右手运动时，NCCI得出的Cz→C3＞C3→Cz比例(96％)远远高于其他三种方法，Cz→C4＞C4→Cz比例(74％)略低于NCC，但在四种方法中准确率位居第二，C3→C4＞C4→C3比例(67％)分别比CGCI、NCC、CGC高15％、8％、11％。综合而言，在想象右手运动时，本发明提出的NCCI方法在一定程度上提高了CGCI和NCC的预测准确率。

4)对比左、右手想象运动的结果，可以发现：想象右手运动时，NCC、CGC和NCCI得出的Cz→C3＞C3→Cz比例分别为96％、81％和70％，明显高于想象左手运动时Cz→C4＞C4→Cz的比例(分别为89％、74％和63％)。考虑到竞赛数据Datasets 2b的所有受试者惯用右手，这一结果表明受试者想象左、右手运动时，大脑有效连接具有不对称性，与文献“Gao Q,Duan X,Chen H.Evaluation of effective connectivity of motor areas duringmotor imagery and execution using conditional Granger causality.NeuroImage,2011,54(2):1280-1288”所揭示的结论一致。

Claims (1)

1.基于滞后阶数自适应选择的多变量因果关系分析方法，其特征在于该方法包括如下步骤：

步骤(1)获取多通道运动想象脑电信号，具体是：采用多导联电极帽采集不同运动想象实验范式下的脑电信号；

步骤(2)采用改进的后向时间选择算法建立多通道脑电信号的动态回归模型；

具体是：给定P个通道、长度为N的脑电信号时间序列{X_j,t:j＝1,2,...,P；t＝1,2,…,N}，首先采用改进的后向时间选择(mBTS)算法估计出每个变量X_j的最优滞后阶数，然后考虑所有变量的历史信息对X_j的影响，X_j的动态回归模型可表示为

其中，u_j,t为模型的回归估计残差；a_jp,l是模型系数，p＝1,2,...,P，l＝1,2,...,k_p；k_p是由mBTS算法估计出的变量X_p的最优滞后阶数；

步骤(3)以步骤(2)得到的动态回归模型为基础，采用模型的残差和系数构建条件因果测度，描述多变量之间因果关系的强度；

具体是：建立X_i→X_j的条件因果测度为

其中，通道号i∈{1,2,...,P}且i≠j；k_i为变量X_i的最优滞后阶数；m＝max(k_p|p＝1,2,…,P)，p＝1,2,...,P；同理，X_j→X_i的条件因果测度表示为

步骤(4)采用步骤(3)的条件因果测度分析运动想象任务时不同区域脑电信号之间的因果关系。

Images