为实现对站立、坐下和静止状态的识别，结合统计差异的网络特征参数和机器学习算法，本文提出了一种该类动作意图识别方法，其技术路线图如图1所示。

图  1  本文技术路线图

13名被试者完成此次实验，被试者均为在校学生，身体健康，没有任何的精神疾病病史，均为右利手，实验之前和被试者签署了实验知情协议书。实验过程中被试者坐在显示器前，显示器显示计时时间，要求被试者先静坐1分钟，然后被试者起立，由静坐转为站立，被试者站立持续1分钟，之后被试者坐下，由站立转为静坐，静止持续1分钟，依次反复交替完成上述动作，如图2所示^[3]。

图  2  站立−静坐实验过程示意图

在该过程中全程记录被试者的脑电信号，持续10 min为实验的一个Session，每名被试者总共完成6个Session的数据采集，各Session中间休息1～2 min。

整个动作实验过程中的数据采集使用SAGA 32+脑电采集系统(荷兰TMSi公司)，在线采样频率为1024 Hz，在线采样频段为0～200 Hz，数据记录中采用平均参考模式，额顶处为接地电极。该系统除了采集30通道的脑电信号之外，还有4组附加通道，通过附加通道接入了两组外接传感器，其中一组是惯性测量单元(inertial measurement unit， IMU)传感器，将其固定在被试者左上腿位置，另一组是4路足底压力传感器(FSR型压力传感器)，以足底压力鞋垫的形式固定在被试者左脚下方，通过IMU和足底压力传感器数据，来确定起立和坐下的起始时间点，实现和脑电数据的同步，传感器安装位置如图3所示。

图  3  足底压力鞋垫和传感器位置固定示意图

实验中记录的信号包括30通道脑电信号、4路足底压力传感器数据以及IMU数据，所有数据的预处理均在EEGLAB中完成。对所有被试者的处理逐个进行，首先将每名被试者所有Session数据整合到一起，然后对数据降采样到250 Hz，并完成0.1～48.0 Hz的带通滤波，使用Cleanline工具完成对信号的线性滤波，根据SAGA系统电极位置完成对所有脑电信号通道的电极定位，并剔除坏电极和无关通道。在完成数据读入的基本处理后，对IMU的各路数据进行低通滤波(Butterworth 滤波器，截止频率10 Hz)，并借助下肢上IMU和脚后跟的足底压力信号确定每个站立动作和坐下动作的起始时间点，作为EEG信号事件的标记点(将事件标记为起立、坐下和静态3个状态)。根据该标记点完成对所有数据的截断，为了讨论动作前后的脑电特征，每个Epoch的长度设置为动作起始点之前4.5 s到之后的2.5 s，最后得到356个Epoch数据段，所有的Epoch数据段都进行了基线校准。然后通过视觉检查剔除受伪迹影响较大的Epoch数据段，利用运行独立成分分析(run independent component analysis, RUNICA)操作完成对所有信号的ICA分解，并基于半自动的独立成分分析选择(semi-automated selection of independent component analysis, SASICA)工具包完成对伪迹成分的判断，并剔除其中的伪迹成分，主要包括眨眼、眼动、肌肉伪迹以及由于动作过程引起的其他伪迹成分。在ICA剔除伪迹之后，如有需要，对坏电极进行插值处理，最后对所有数据进行平均重参考，从而完成数据的预处理过程。

行走是人类一项复杂的认知活动，涉及视觉判断、运动控制、空间感知和规划等多个方面，需要整个大脑不同脑区间的协调和配合。在一个动作执行之前，通常大脑就表现出一定的动作意图，通过讨论动作执行前后相当时间内不同脑区间的相互依赖关系，有助于全面了解动作执行过程中大脑的活动状态，并有望发现针对动作意图检测的新特征。基于多层脑功能网络的复杂网络分析方法为此类动作过程对应脑区间依赖关系的研究提供了工具。

表征脑电信号之间相互依赖关系的度量方法有很多种^[14]，其中同步似然分析(sychronization likehood, SL)被广泛应用于非线性和非平稳信号间依赖关系的分析中，特别适合于脑电信号的分析，对EEG信号中的线性和非线性特征都具有较好的表征能力，且对容积传导等问题具有较好的鲁棒性^[15-16]。SL是神经生理学数据中估计广义同步性的常用指标，它的理论依据是系统在某一时刻的状态可通过一个插值向量来表示，当该状态再次出现时，可以用一个相似的插值向量来对其进行表示。因此，本文采用SL来构建不同通道脑电信号之间的相互依赖关系，从而构建多层脑功能网络。分为以下4个步骤。

1) 对来自通道C的EEG信号$ {x_{C,i}} $(i =1,2,···N，N为样本点数)进行重建，得到m维空间嵌入之后的信号：

式中，l为时间延迟。同理，对通道D的EEG信号进行重构得到$ {{\boldsymbol{X}}_{D,i}} $。

2) 对通道C上的第i个样本，设定各向量之间彼此的距离小于s的概率为：

式中，$ |\cdot| $表示两个相空间矢量的欧氏距离；$ \theta $为Heaviside阶梯函数，满足当$ x \leqslant 0 $时，$ \theta (x) = 0 $，当$ x > 0 $时，$ \theta (x) = 1 $；对时间窗的控制通过$ {\omega _1} $和$ {\omega _2} $来控制，且满足$ {\omega _1} \ll {\omega _2} \ll N $。

3) 定义通道C上在采样点i处的临界距离为：$ {s_{C,i}}:P_{C,i}^{{s_{C,i}}} = {p_{{\rm{ref}}}} $，其中有$ {p_{{\rm{ref}}}} \leqslant 1 $；同样可确定通道D中信号的邻接距离：$ {s_{D,i}}:P_{D,i}^{{s_{D,i}}} = {p_{{\rm{ref}}}} $。

4) 在完成上述步骤之后，就可以实现对通道C和通道D之间EEG信号同步似然指数SL的计算，其表示的是通道C信号和D信号同时重现其状态向量的可能性，即当$ {{\boldsymbol{X}}_{C,j}} $和给定$ {{\boldsymbol{X}}_{C,i}} $之间的距离小于邻接距离$ {s_{C,i}} $，同时$ {{\boldsymbol{X}}_{D,j}} $和$ {{\boldsymbol{X}}_{D,i}} $之间的距离也小于$ {s_{{D} ,i}} $的可能性。本文在计算SL时，对应参数依次设置为：l=10、m=10、$ {\omega _1} = 100 $、$ {\omega _2} = 400 $、$ {p_{{\rm{ref}}}} = 0.01 $。

对某段时间内脑功能网络拓扑结构动态关系的一种表征就是对该时间段内信号按一定时间窗进行截断，依次构建功能网络并形成超邻接矩阵，然后借助多层网络分析方法进行特征分析。为了通过脑功能网络来表征站立、坐下以及静止状态的脑信号特征，首先提取了和运动意图相关的delta慢波信号，如图1所示，然后根据站立和坐下以及静止状态的时间点，将3类状态对应的delta波段信号分割为[−4.5 s, 2.5 s]的Epoch数据段(长度为7 s)，并以3.5 s的时长作为时间窗，0.5 s的步长对每个Epoch依次进行截断，从而得到8个时间段的子窗口(依次为w1,w2,···,w8)，这样对一个Epoch样本进行操作，就可以得到30×875×8的数据样本，然后对每个窗口内以及每对窗口间使用SL来计算两个通道间的相互依赖关系。图4为层间和层内各通道SL计算示意图，通过层间和层内各通道之间SL计算最终可以得到多层网络结构并以超邻接矩阵的形式表示。

图  4  多层脑功能网络构建

图5为3种状态下所有被试、所有Epoch数据段对应超邻接矩阵的总平均显示，该矩阵为240×240的SL值矩阵表示，对角线位置为8个子窗口(w1～w8)的窗内连接关系(层内连接)，而上三角或下三角矩阵为不同子窗口之间的连接关系(层间连接)，可以看到层内连接强度要远高于层间连接。由于图示超邻接矩阵中SL为总平均值，3种状态之间的差异并不明显。为此，图6计算了该超邻接矩阵中SL的总平均值，可以看到3种状态之间，SL均值在动作发生期间要高于静止状态，且坐下和起立动作对应均值接近，但是统计分析并没有发现3种状态之间存在显著差异。

图  5  3种状态对应EEG信号SL超邻接矩阵

图  6  3种状态SL超邻接矩阵总平均值

为了就整个动作过程对应的特征量进行量化分析，本文重点讨论4类多层网络特征参数，首先是层−层相关系数(layer-layer correlation, LLC)，将其定义为两个网络层之间的皮尔逊相关系数^[17-18]，用于表示各网络层之间，特别是不同子窗口对应网络层之间的相关性程度。为了就多层网络中不同层之间整合和分离程度进行度量，文献[17]引入了一个度量指标，该指标量化了第m层中的链路关系仍然存在于第n层中的条件概率(conditional probability, CP)，将其定义为：

式中，$ {{\boldsymbol{A}}_m} $(m=1,2,···,8)为层内邻接矩阵。如果P值较高，表示各层之间整合关系较强，而如果各层之间关系相对比较独立，则P值较低。此外，不同层之间的分离特征也可以用多层模块值来进行量化，文献[19-20]将其定义为：

可以看到Q值受到参数$ \gamma $和$ \omega $的控制，它们分别决定了网络结构中的时间连通性和空间连通性。式中，$ \mu $是多层网络中的整体链路权重；$ {\gamma _l} $对应于第l层的拓扑分辨参数；$ {\omega _{jml}} $是节点 j 在第 l 层和第 m层之间的时间耦合参数；$ k_i^l = \displaystyle\sum\limits_j {w_{ij}^l} $是节点 i 在第l 层上的连接强度。节点i和节点j在同一个模块时，$ \delta ({{\boldsymbol{A}}_{il}},{{\boldsymbol{A}}_{jl}}) $和$ \delta ({{\boldsymbol{A}}_{il}},{{\boldsymbol{A}}_{jm}}) $的取值为1，否则取值为0。

最后，多层网络还可以使用多层参与系数(multilayer participant coefficient, MPC)来量化其特征。该系数与熵有关，主要用来量化给定节点i的链路是均匀分布在各个层中，还是主要集中在一个或几个层中^[17-18]。节点i的MPC值被定义为：

式中，$ {o_i} = \displaystyle\sum\limits_{l = 1}^M {k_i^l} $是节点i跨越M个层的整体连接强度。当节点i的所有连接边都位于同一个层时，MPC_i=0，仅当节点i在M个层中的每个层上具有完全相同的边数时，MPC_i=1。将整个多层网络的MPC值定义为各个节点该值的总平均为：

上述SL值构成的超邻接矩阵表征了层内和层间的相互依赖关系，是对同一时间段或不同时间段各通道之间相互依赖关系的一种表征，在这种连接关系中也存在着大量的弱连接边，可能是由于噪声或其他非特征性相互作用造成的，并不能有效表征不同动作状态下脑区间的真实依赖关系。因此，在一般性问题的讨论过程中，需要设定一个合适的阈值对上述超邻接矩阵中的连接边权值进行筛选，保留在不同状态间具有一定区分度的重要连接边。为此，首先将3种状态对应SL超邻接矩阵权值展平，并组合在一起构成一个一维向量，然后通过选择该向量取值分布图上的上分位点作为最终的阈值点，图7为该向量的频率直方图分布以及确定的阈值(阈值T = 0.1)。

对上述3种状态对应所有超邻接矩阵使用本阈值进行重点权值边的筛选，发现由于层间连接都较弱，保留下的重点连接边都表征了层内各通道之间的连接关系。因此，对各层内连接以单层网络的方法进行对比分析，重点考虑了各层网络的度、聚集系数、全局效率、同配系数、传递性系数和节点的介数中心性^[14,21]。网络中最基本的一个参数就是节点度(degree, D)，定义为与该节点相连的边的权值之和，并将网络中所有节点度值的平均作为网络的度(average degree, AD)。同时，主要对网络的聚集系数(clustering coefficient, C)、全局效率(global efficiency, GE)、网络的同配系数(assortativity coefficient, ASS)、网络的传递性系数(transitivity, Tra)以及节点的介数中心性(betweenness centrality, BC)等网络参数展开分析。

图  7  3种状态对应SL超邻接矩阵值统计分布

在得到如图5所示各状态对应的超邻接矩阵之后，本文分别计算了各Epoch片段对应的超邻接矩阵的4种多层网络特征参数，图8为对应的Q值和MPC值，可以看出Q值保持在0.8左右，而MPC值在0.9以上，两组参数值都较高，表明不管是静止状态，还是起立和坐下动作，各网络层之间的连接关系较为密切，整个动作过程中大脑活动表现出很强的连续性。对3种状态之间的统计分析也发现，Q值在3种状态间有较为明显的差异(t-test, p<0.05)，而MPC值并没有表现出显著性差异。

图9分别对应3种状态下各网络层之间条件概率CP矩阵和层−层相关系数LLC矩阵。可以看到各层之间的CP值都普遍较高，均在0.6以上，表明各子窗口对应网络层之间存在较强的整合关系，不管是相邻时间窗的网络层(如w1和w2)，还是间隔较远的时间窗的网络层(如w1和w8)，各层之间均存在较为紧密的联系，但是3种状态之间并无显著差异。相比之下，图9中两层之间相关性值情况有所不同，可以看到两层之间的LLC值随着两层间距离的增加而降低，表现出相邻层有较高的LLC值，而间隔较远的层之间有较低的LLC值。为了能够找到可应用于不同运动状态分类的特征量，对LLC值同样进行了统计分析，图9中标记为两种状态之间存在统计性显著差异的情况(t-test，p<0.05，其中 “*” 为起立 vs. 坐下， “★” 为起立 vs. 静态，“▲” 为坐下 vs. 静态)。可以看到，对于总平均后的LLC矩阵值并没有明显不同，但是统计分析发现在w4和w6,w7,w8之间，以及w6和w7,w8之间的相关性存在统计性显著差异，这可能和站立和坐下动作开始的时间点有一定关系。特别是从w4子窗口开始，出现了起立和坐下的动作起始点，因此和后续子窗口的网络在层相关性上表现出显著性差异。

图  8  不同状态超邻接矩阵对应Q值和MPC值

图  9  不同状态对应网络层−层相关系数矩阵

为了突出上述网络连接中的重点连接边，将阈值设定为0.1，完成对上述3种状态所有Epoch数据网络权值的筛选，大于该阈值的连接边保留，小于该阈值的连接边权值设置为0。权值筛选后的网络最终为各层内连接网络，层间连接权值相对较小，因此均被忽略。图10为3种状态各子窗口对应阈值筛选后总平均网络连接拓扑图。可以看出，3种状态在8个时间窗上拓扑结构变化比较平稳，对站立和坐下状态，均从w4开始，网络连接较之前的w1、w2、w3开始变得稀疏，并且一直持续到最终的w8，这和图10中网络层之间相关系数显著性差异结果一致。本文认为这主要是由于在整个时间窗的移动过程中，站立和坐下动作的起始点最开始均出现在w4窗口，而在站立和坐下过程中的动作电位相应地影响了不同脑区间的这种依赖关系，且在动作完成之后的一段时间内，这种依赖关系得到了一定程度上的保持。从图中也可看出，在静止状态下，网络拓扑结构保持平稳，不同脑区间的拓扑关系几乎没有发生变化，整个网络拓扑结构的变化同起立及坐下状态有一定区别，但无法通过该动态拓扑图来区分起立和坐下状态。

图  10  3种状态各子窗口对应脑功能网络拓扑连接

为了便于进行量化分析，确定可用于不同状态识别的新特征量，本文重点分析了3种状态下各窗口对应的网络特征参数。如表1所示，分别计算了聚集系数C、全局效率E、平均度AD、同配系数Ass、介数中心性Bet和网络传递性系数Tra这6类基本的网络特征参数。其中，“*”表示t-test，p <0.05。由表可知，动作状态对应网络特征参数在8个递进子窗上整体呈现出变小的趋势，但是变化差异很小，而在静止状态网络参数变化没有明显的趋势性。对上述各网络参数在3种状态之间的统计分析表明，只有w7窗口对应的网络传递系数Tra表现出显著的统计性差异。

通过上述对各网络特征参数的对比和统计分析，本文选择了多层网络特征参数Q、有显著差异的层间相关系数LLC以及w7子窗口对应网络传递性系数Tra构成特征向量：

使用支持向量机(support vector machine, SVM)、逻辑回归(logistic regression, LR)、线性判别分析(linear discriminat analysis, LDA)和朴素贝叶斯(navive bayes, NB)4类机器学习算法对该特征向量进行测试，完成对两类状态的识别，在机器学习分类过程中将70%的数据作为训练集，剩余30%作为测试集。图11为5折交叉验证得到的平均分类准确率。可以看到，针对两类状态的分类准确率均在71%以上，其中SVM分类器的效果最好，针对站立和坐下、站立和静态以及坐下和静态的分类准确率分别达到74.4%，77.0%和76.9%，证明了该方法的有效性。

图  11  不同分类器对应分类结果

本文通过设计静坐到站立的动作转换实验，完成对被试者在起立、坐下和静止状态下对应不同脑区脑电信号的同步似然分析，通过构建动作执行前后8个时间窗口上的脑功能网络，以超邻接矩阵表示的多层网络形式展开对不同动作状态下网络拓扑结构的特征分析。基于对多层网络特征和各子窗口网络参数的统计分析，确定对不同动作状态敏感的网络特征参数，并使用SVM、LR、LDA以及NB等分类算法实现对起立、坐下和静止状态的分类识别，分类准确率最高达到77%，证明了基于动作前后动态网络特征进行下肢运动意图检测和识别的可行性。针对多层动态脑功能网络的分析结果表明，站立和坐下动作的发生会弱化脑区间的相互依赖关系，从而导致网络拓扑连接结构变得逐渐稀疏，该结果也为基于脑电的动作意图检测提供了新思路，对开展基于脑机接口的下肢助行和康复机器人的开发具有一定的参考价值。

致谢：本文得到了兰州交通大学天佑青年托举人才计划的基金资助，同时也感谢博睿康开放课题(BRKOT-LZJTU-20220329C)对文章的支持。