脑机接口(BCI)是大脑与外界信息交互的桥梁，目前主要通过头皮或者皮层内传感器采集脑电信号(EEG)，对信号进行解码，进而完成意图检，其中运动想象(MI)识别在BCI中占有非常重要的地位．EEG中包含了丰富的运动信息，可以从中提取相关运动想象特征．文献[1]设计了情绪刺激实验，分别构建了μ节律和β节律的脑网络，以区分左右手想象任务．合理选择EEG信号频段对运动想象解码有重要影响，但目前大部分运动想象研究基于脑电信号的频域特征，忽略了人脑处理运动想象是一个连续时变过程的事实，脑网络随时间发展产生的动态变化仍有很大的探索空间．目前，有关网络的时间动态分析也已经在意识障碍等脑部神经疾病研究中得到了应用．文献[2]在对意识障碍与健康人分类研究中指出时间的动态连接变化对于相关分类有着重要意义．文献[3]基于时变特性进行多层脑网络拓扑，使用统计指标的多视角多层脑网络拓扑属性计算方法进行脑疾病分类．综上，时间测度脑电的动态结构变化与频率测度一样，提供了运动想象识别的重要特征．同时，关于功能脑网络的研究也出现了新的动态，文献[4-5]介绍了使用新型的多层网络模型对人的大脑进行空间描述所具备的潜在优势．文献[6]分析了多层网络中网络模块的时间衍化、动态重构及功能网络的灵活性，并且整理了相关数学定义及常用的多层网络评估手段．文献[7]依据时间序列构建多层网络对人体运动姿态改变进行了网络拓扑研究，发现了不同姿态的运动过程中，不同时刻的功能脑网络产生了显著拓扑差异．文献[8]使用神经成像技术和EEG建立了功能脑网络，比较了脑创伤患者于对照组的平衡神经相关性，揭示了不同运动状态下脑区功能连接具有不同的聚类组合．结合上述内容，为解决运动想象识别问题中特征维度单一、缺少运动想象时间维度动态过程分析的问题，本研究提出一种基于多层时变功能脑网络的运动想象特征提取和识别方法．运用多层网络模型，将运动想象脑电数据按特定时间窗口与步长定义为一系列时间窗口的串联，构建多层时变功能脑网络．通过量化分析网络特征与网络拓扑，确定多层网络中的核心层，选择其核心单层网络特征描述人脑在不同运动想象状态下的空间功能结构连接，选择多层网络特征描述随时间产生的动态变化与分离整合特征．最终联合两类特征进行运动想象识别任务，在公开数据集上取得了较好的识别成果．1 基本理论1.1　方法路线图为实现不同MI的分类任务，本研究构筑多层时变网络，通过网络拓扑与异层相关性分析确定核心差异网络层，结合具有统计差异的网络特征参数进行机器学习，最终达到分类目标．其理论方法路线如图1所示．10.13245/j.hust.240450.F001图1基于多层时变网络的MI识别方法路线图首先依据滑动窗口方法，确定窗口长度与步长，保证每个相邻窗口有部分重叠且窗口长度大于步长，将单次运动想象采集划分为多个连续的时间窗口；然后将不同时间窗口，按照相似性计算不同窗口和通道间的功能连接，构成超邻接矩阵，通过阈值筛选，保留连接强度较大，信息交流紧密部分，去除弱连接与干扰信号，组成多层时变功能脑网络．通过对不同网络层的拓扑分析与网络评价指标分析，确定差异最明显的核心网络层．最终选取核心网络层的网络特征，结合多层网络特征组合成为分类特征，完成机器学习及识别任务．1.2　网络模型1.2.1　多层模型定义多层网络在张量形式上可定义为M=(G,C);G={Gα|α∈N};Gα=(Vα,Eα),式中：M为多层网络；G为一组图；C为连接不同图节点一组边；N为不同网络层的集合；Gα为α层上的图；Vα为Gα的节点集合；Eα为连接边的集合．在不同层α和β上，连接边集合Call可以表示为Call={Eα∈VαVβ|α≠β}．同时，多层网络可以被降维展开成对应的超邻接矩阵，即A={Zijαβ}，式中：A为一个矩阵；Zijαβ为第α层的i节点与第β层的j节点之间的连接．对于层数为q的多层网络，A是一个主对角线有q块，负责层内连接；非对角线上有q(q-1)块描述不同层间连接的对称矩阵，可得A=E11E12…E1qE21E22…E2q︙Eq1Eq2…Eqq.1.2.2　单层网络特征a．聚类系数网络中用聚类系数(CC)来表示局部信息的处理，平均聚类系数反映了网络中两个节点之间存在连接边的可能性，可以表达为Li=∑i≠j,m≠iaijajmamiki(ki-1)，式中：Li为节点i的聚类系数；ki为节点i的单层网络节点度；aij为节点i与节点j之间的连接边．对网络所有节点聚类系数求平均即得到网络的聚类系数L，可表示为L=1N∑i=1NLi．b．节点度节点度是衡量一个节点与其他节点共享的实际链路数，网络的平均节点度K 就是所有网络节点度的平均值，是描述网络网络规模大小的重要指标，可表示为K=1N∑i=1Nki．c．平均路径长度平均路径长度(CPL)即任意两个节点之间距离的平均值，它表征了单层网络中某一层的传输信息的能力，可表示为D=1N(N-1)∑i≠jdij，式中：D为网络所有节点的平均路径长度；dij为节点i与j的距离．1.2.3　多层网络特征a．多层聚类系数在多层网络中，同一个节点在不同层上有不一样的聚类系数，单层聚类系数难以表征整个网络聚类系数，因此在多层网络中节点i的多层聚类系数(MCC)可以表示为Li,1=∑μ∑k≠μ ∑i≠j,m≠iaijajmamiq-1∑μkiαkiα-1，式中：Li,l 为第l层中节点i的多层聚类系数；aijα为节点i，j在第α层上的连接；kiα为节点i在第α层上的节点度．该MCC表达式反映了跨越不同频段的局部信息处理，若节点i为中心的边可以在层间形成三角，则节点i的多路复用聚类系数为1，反之为0．对多层网络所有节点求平均可得网络的多层聚类系数U，表示为U=1N∑iCi,1．b．多层参与系数多层网络中的等效度量是重叠度，即每个层该节点的节点度按权相加，这在多层网络中是一个定值，但是不同层中有不同的度分布，因此引入多层参与系数(MPC)来量化节点对不同层信息交流的参与程度，可表示为pi=qq-11-∑L=1M(kiα/oi)2，式中：pi为节点i的多层参与系数；oi为节点i的重叠度，即节点i在所有层上的节点度之和；kiα/oi衡量了节点度在每层所占比例．pi表征了节点i的连接在各层是否均匀，与网络多层聚类系数一样，对多层网络所有节点求平均可得网络的多层参与系数P，可表示为P=1N∑ipi．c．异层相关系数异层相关系数(CP)，用于表示不同子窗口对应网络层之间的相关性程度．为了对多层网络中不同层之间整合和分离程度进行度量，引入一个度量指标H，该指标量化了第α层中的连接关系仍然存在于第β层中的条件概率，可表示为H=∑ijaijαaijβ/∑ijaijβ，式中a α为α层内邻接矩阵．若CP较高，则表示各层之间整合关系较强；若CP值较低，则各层之间关系相对比较独立．d．多层模块值对于多层时变网络，引入一个时间相关的参数将层与层之间联系起来，多层模块值能够量化不同层之间的分离特征，它表征了一段时间中网络聚类集群分类变化的程度，其模块化计算公式为Q=12θaijα-γnkiαkjm2zαδ+εωjαβδ，式中：θ为多层网络中整体链路的权重；zα为第α层的所有节点度之和；ωjαβ为节点j的α层与β层时间耦合系数，它与结构分辨率参数γ通常被设置为1，当节点i与节点j在一个模块时，δ和ε取值为1，否则为0．2 网络建立及特征分析2.1　实验数据本研究的MI数据来源于BCI Competition Ⅳ Data SetsⅠ[9]，该数据集共有4名被试者，分别被标号为A，E，F，G，分别测量了左手和脚运动想象，本研究选取了该4名被试者数据作为实验数据．整个实验中，每名被试者进行左手、脚各100次运动想象，单次运动想象持续8 s，连续进行10次运动想象后休息30 s．共有59个采集电极，符合国际10-20导联标准．所采信号采样频率均为100 Hz，并且经过50 Hz的陷波滤波．整个MI采样持续8 s，前2 s保持静息，2～6 s为MI阶段，被试者根据屏幕显示箭头提示做出相应MI，最后6～8 s保持黑屏，被试者处于休息状态，结束该次采集．2.2　数据处理本研究使用Matlab的EEGLAB工具进行数据预处理．首先对数据集所提供的采集电极空间信息，根据10-20标准进行通道定位，得到正确的空间分布，便于后期拓扑处理；然后对59个通道数据进行电极共平均参考，单个电极减去所有电极的平均值消除噪声，该方法可使参考电极不受位置影响，消除电极距离不相等、电压可比性与双侧半球不对等的问题．在采集电极中，为使左右功能脑网络在拓扑结构中更加明显，选用其中41个导联，重新排列为1～41号，最终电极分布如图2所示．10.13245/j.hust.240450.F002图2电极分布图之后使用FIR滤波器，对信号进行8～14 Hz带通滤波，得到μ节律信号，同时依据需要，在每段运动想象任务的200与600 ms位置重新标记，截取作为运动想象任务数据段．利用独立成分分析法完成对所有信号的ICA分解，并且剔除其中的伪迹成分．最终可以得400组2～6 s的运动想象脑电数据，本实验选择窗口长度为1 s，步长为0.5 s的滑动窗口，并将每段运动想象脑电信号划分为7个时间子窗．2.3　多层时变功能脑网络建立锁相值(PLV)被广泛地应用于非线性和非平稳信号间依赖关系的分析中，对EEG信号中的线性和非线性特征都有很强的表征能力．本研究采用PLV来构建不同通道EEG信号之间的相互依赖关系，其计算公式为R=eiφxy(t)t=cosφxy(t)t2+sinφxy(t)t2;φxy(t)=φx(t)-φy(t),式中：R为PLV的值；∙表示刻度为时间的平均值；φxy(t)为x(t)和y(t)瞬时相位差．PLV取值范围为0~1，其中：1表示具有稳定的相位差即相位同步；0则表示无同步．每个运动想象任务所包含的4 s信号均已划分成为了7个滑动窗口(W1，W2，…，W7)，将脑电信号采集电极作为节点，节点的PLV作为边，构建每个时间窗口下的层内连接与每个时间窗口之间的层间连接．通过锁相值的连接可以获得多个通道信号相关矩阵，矩阵中的值为任意两个通道信号之间的锁相值．因此，每位被试者的每次运动想象均会得到7个描述自身时间窗口层内连接的矩阵，构成超邻接矩阵的主对角线，同时会得到42个描述不同滑动时间窗口之间的层间连接的矩阵，组成超邻接矩阵的其他部分．将数据平均后作为分析对象，会掩盖实验过程中单体的噪声，故此本次选择被试者F第1次左手与第6次脚运动想象作为拓扑样本说明，该被试者选择是基于最接近实验群体平均值的个体做出分析，使其具有代表性同时具有个体的不确定性．图3所示为被试者F第1次左手与第6次脚运动想象W4相关矩阵及超邻接矩阵．10.13245/j.hust.240450.F003图3被试F左手与脚运动想象W4窗口相关矩阵与超邻接矩阵观察图3(a)和图3(b)可知：在W4窗口下，左手和脚的运动想象相关矩阵中有着显著的颜色分布差异，每个窗口中左手和脚的空间功能连接差异十分明显，主要集中于通道1～10与28～41之间．对比两者可见，脚运动想象通道1～10内部间连接与通道28～41之间有着更深的颜色分布即有着更强的连接，W4窗口中运动想象脚动作有着更强的连接关系．对比图3(c)和图3(d)两张超邻接矩阵可以看出：无论是相对应的时间窗口，还是不同窗口之间，脚运动想象均有着更深的颜色分布，说明脚运动想象总体上有着更紧密的连接，更显著的动态变化．但无论是左手还是脚，不同窗口之间即不同层之间同样包含大量功能连接，这是单层网络所无法反映的．本实验中基于PLV构成的超邻接矩阵表征了不同时间窗口的层内与层间相互依赖关系，在这种连接关系中也存在着大量的弱连接，它通常是由噪声或者其他非特征性相互作用造成的，并不能表示该动作脑功能网络的结构关系与信息传递．因此，设置一个合适的阈值对上述超邻接矩阵中连接边的权值进行筛选，能更好表征功能脑网络．对两种运动想象任务对应的超邻接矩阵权值展平，并组合成为一个一维向量，通过选择该向量取值分布图上的上分点位0.7作为最终的阈值点，即阈值T=0.7．2.4　网络拓扑与网络特征参数计算根据所选阈值，完成重点连接边筛选，大于阈值的保留，小于该阈值的连接边设置为0，层间连接权值相对较小，因此阈值筛选后的网络最终为各滑动窗口的层连接网络．图4所示为对应两种MI多层网络各层之间的条件概率异层相关系数的值．10.13245/j.hust.240450.F004图4被试者F左手与脚MI中各网络层之间的异层相关系数可以看出：无论是左手还是脚的运动想象，各个层之间的异层相关系数值均普遍较高，大部分都在0.5以上，前后连接窗口可以达到0.9以上，表明不管是相邻时间窗的网络层(如W2和W3)，还是间隔较远的时间窗的网络层(如W1和W7)，均存在较为紧密的联系，有很强的整合关系．比较图4(a)和图4(b)可知：总体上脚的运动想象有着更高的异层相关系数，代表着脚运动想象功能脑网络中有着更为明显的功能连接组合特征，符合先验的运动想象事实．单独观察左手和脚运动想象的异层相关系数矩阵可见，左手运动想象的W4窗口和脚的W4窗口均有着最高的异层相关系数总和，代表该两者的W4网络层与其他网络层有着更紧密的联系与变化．为进一步分析时变网络动态变化及功能结构连接变化，对每个时间子窗建立网络拓扑并按时间排序．图5为其连接筛选后总平均网络连接拓扑图．10.13245/j.hust.240450.F005图5被试者F第1次左手及第6次脚运动想象脑网络拓扑图可以看出：左手MI总体的连接相比脚MI更为稀疏，两种状态在7个时间窗口上的动态变化较平稳，整个过程中产生了两次明显聚类改变，将整个过程分为三段聚类．左手MI，W3与W5时间窗口产生了两次较明显的聚类分布改变，将整个过程分为W1～W2(右斜稀疏)、W3～W4(左斜密集)和W5～W7(右斜密集)；脚MI，在W4窗口前后产生了明显的聚类改变，也将整个过程分成了W1～W3(团状密集)、W4(右斜稀疏)和W5～W7(左斜密集)三类，W4窗口均处于两类范式下聚类组合变化的关键时刻．对比两者可以看出：两组功能脑网络的动态变化与聚类组合有着明显的差别，左手网络聚类分布呈现疏-密-密变化，而在同一时间段中脚网络呈现密-疏-密的不同变化趋势，聚类分布中左手网络呈现右斜-左斜-右斜的变化，脚网络则是中心-右斜-左斜的变化过程．在两种状态下，W4窗口在动态连接与聚类分布的改变中都有着最为显著的网络拓扑差异，结合上面所提到的W4有着最高的异层相关系数，有着更显著的功能连接组合特征，因此选择W4为本次拓扑的网络核心层．为进一步量化分析，确定提取核心网络特征与层间网络特征是否具备显著特征，本研究求出不同被试不同MI的平均节点度K、平均聚类系数L、多层参与系数P、多层聚类系数U及多层模块度Q，并最终求其平均值，结果如表1所示．可以看出：无论是MCC还是MPC，两者均较高，表明无论是左手还是脚的运动想象，在整个过程中大脑的活动表现较强的连续性，这与图4所表示的各时间窗对应的网络层之间均有较强的整合关系一致，而较高的Q值表明了在整个运动想象过程中功能脑网络有着显著的聚类集群分类变化，以及明显的动态连接变化与聚类重组，进一步支撑了图5的变化分析．10.13245/j.hust.240450.T001表 1四名被试者不同运动想象任务W4窗口下单层与多层网络特征参数被试任务窗口KLPUQA手W414.9770.6750.6870.6160.732脚W414.8810.6730.8100.6320.826E手W413.9270.6810.7280.6330.774脚W413.5710.6020.7590.6200.699F手W411.0990.6040.8310.5440.759脚W410.5150.5920.8620.5390.835G手W411.7450.6190.8540.5480.682脚W410.9280.5990.8610.5370.729选取被试者F运动想象的核心变化层W4，其相关网络特征中的平均节点度、平均聚类系数用以描述核心层的空间连接特征，同时选择构筑方式相似的MPC与MCC来描述多层网络的动态变化与分离组合特征．挑选出来的网络特征分布如图6所示．10.13245/j.hust.240450.F006图6被试者F单层与多层网络特征分布各个特征分布具有一定的区分度，不同维度下特征长度不同，但同类运动想象特征分布趋势大致相同，可见不同的运动想象在本研究建立的多层时变网络上具有一定的区分度．3 基于多层时变功能脑网络特征的运动想象分类任务3.1　运动想象实验分类结果本研究构建的多层时变功能脑网络，通过上述的网络拓扑分析与特征参数的统计分析，选择多层网络参数P和U描述网络的动态变化，选择核心层W4的K和L描述网络层的空间功能结构，构建特征向量Fall．为了实验对比，将单次运动想象任务2~6 s时间段的脑电信号求平均构建脑功能网络，并且运用其K，L和D组建成特征向量Fs．构筑两种特征向量分别为：Fs=[L,K,D];Fall=[L4,K4,U,P].选用支持向量机分类器(SVM)，采用5折交叉验证法进行验证．从A，E，F，G两名被试者的样本数据中，随机选取160组左手和右手运动想象样本，合计640组作为训练集，剩下160组作为测试集．重复进行20次分类，取20次实验的平均正确识别率．再重新选择640组数据作为新训练集，剩下160组作为测试集，重新实验共重复3次，识别结果如表2所示．10.13245/j.hust.240450.T002表2不同特征向量的识别率特征向量识别率分类一分类二分类三平均Fs82.50±1.0383.75±1.1781.33±1.4282.53±1.82Fall87.50±0.8690.53±0.5389.38±0.8989.14±1.08%由表2知：使用Fs特征向量的平均识别率为82.53%，方差为1.82%，最高识别率为83.75%；当使用Fall特征向量进行分类时平均识别率为89.14%，高出Fs特征向量6.61%，方差为1.08%，其中最高识别率达到了90.53%．比较两者结果，多层特征拥有更高的识别率，与更小的方差，表明Fall在获得了高识别率的同时具备更好的稳定性，验证了运动想象过程中不同时间窗口中的层间连接拥有单层网络特征所不能表征的动态信息和功能结构特征，联合使用单层网络特征与层间网络特征能实现更好的分类效果．3.2　与相关文献方法对比为了验证本研究的有效性，将本研究与现有常用方法在同一数据集上进行对比，选取同样的SVM分类器，进行10折交叉验证，结果如表3所示．10.13245/j.hust.240450.T003表3不同特征提取方法的平均识别率对比特征提取方法平均识别率本研究89.14滤波器组共空间[10]84.79希尔伯特-黄变换[11]86.36双数复小波变换[12]88.27%3.3　其他运动想象范式分类效果为了验证本研究对于不同范式下运动想象分类的适用性，改变被试范式为左手-右手运动想象，运用相同方法进行数据处理、分析、特征选择，最终进行基于SVM的运动想象分类．数据来源于同一数据集的被试者B．经处理后，其相关网络特征如表4所示．10.13245/j.hust.240450.T004表4被试者B被试左右手运动想象任务不同窗口下单层与多层网络特征参数任务窗口KLPUQ左手W217.9750.7140.9100.5760.687右手W216.7250.6960.9270.6920.753依据拓扑分析，选用W3窗口中的K与L，联合P与U组合成为分类所用的特征向量，选用支持向量机分类器，采用5折交叉验证法进行验证．最终所获得的平均识别率为86.67%由表3可知：滤波器组共空间模式，希尔伯特-黄变换，双数复小波变换方法的正确率分别为84.79%，86.36%和88.27%，均低于本研究所提出的多层时变网络特征提取方法，证明了本研究提出的时间测度提供了重要的运动想象特征，可以有效提高左手和脚运动想象的识别准确率．由表4可知：本方法不仅针对于左手和脚运动想象有着较好的分类效果，而且在左右手分类上也具备一定的准确率，证明了本方法在不同范式的运动想象任务分类中具有较好的适用性．4 结语本研究提出一种基于脑电信号的多层时变功能脑网络的运动想象特征提取方法，并使用左手和脚的运动想象分类进行验证．首先确定合适的窗口长度与步长，将信号分为多个部分重叠且连续的时间窗口，依据锁相值计算节点间功能连接，构建多层时变功能脑网络；然后运用网络拓扑与统计参数量化分析，选定特征最明显的核心差异层，计算其K，L，P和U，组成多层特征向量，分别描述网络的空间功能连接、动态变化与分离整合特性．分类结果表明：本研究所提取特征向量相较单层网络特征有更好的识别率，能说明时间测度下脑电的动态变化提供了重要的运动想象特征，充分反映运动想象中人脑是一个动态改变、分离整合过程这一事实．通过与其他功能脑网络特征提取方法对比，本研究具有更高的识别率，证明了本研究的有效性．