脑机接口(brain computer interface, BCI)系统可以不依赖大脑神经通路，通过采集人的脑电信号发送至计算机，借用机器学习算法进行意图识别，控制外部设备，实现人与外部环境的交互^[1-2]. 脑机接口技术在运动康复、改善、替代方面具有重大的应用前景，给运动障碍患者带来康复的希望. 近年来，国内外关于脑机接口的研究蓬勃发展，通过采集脑电图(electroencephalogram, EEG)信号不同成分进行研究已经实现对外骨骼^[3]、机械臂^[4]、智能轮椅^[5]等外部设备的控制. 目前，比较成熟的关于EEG信号实现BCI控制的研究主要有稳态视觉诱发电位(steady-state visual evoked potentials, SSVEP)^[6]、P300^[7]、运动想象(motor imagery, MI)脑电信号^[8]等.

现有的BCI系统虽然能有效控制外部设备实现人机交互，但也存在一些问题. 最主要的问题是系统实时性低，响应时间长，无法实现自然控制. 由于被试需要通过眼睛接收外部的视觉刺激，在使用过程中会产生视觉疲劳，严重的甚至会诱发眼部疾病. 在实际运动发生时刻的−1 500 ms前后，大脑已经被激活，会产生8~10 μV负向偏移的运动相关皮质电位 (movement related cortical potentials , MRCP )^[9]. 运动前的这段脑电信号中包含有能够预测接下来要发生运动的相关信息，而且不依赖外部刺激^[10]，这可以解决目前BCI系统存在的问题，实现对外部设备更自然、更流畅的控制，是未来BCI系统的发展趋势. 基于MRCP的意图识别也面临一些问题，如信号较微弱、抗干扰性差. 传统特征提取方法（如频带能量、小波变换）建立在EEG信号是线性的假设基础之上，所得结果却不尽如人意^[11-12]. 大脑神经元是高度非线性的组织，脑电信号是神经元放电的时空累积. 当被试产生运动意图时，大脑的神经元集群放电模式会发生变化，貌似不可重复、不可预测的脑电信号中，可能隐含着简单的规律，这就是大脑的混沌现象^[13].

针对传统BCI系统分类准确率相对较低、实时性较差的问题，在传统线性特征提取方法的基础上，基于大脑的混沌现象，将非线性混沌特征样本熵引入MRCP的研究中. 通过计算大脑状态复杂度的变化，表征不同大脑活动的差异，实现对手部运动意图预测. 基于手部运动意图的提前预测，提供提高BCI系统实时性的控制策略.

如图1所示为运动相关皮质电位采集实验过程及设备. 理论上需要大脑参与的运动都有可能会产生相关的MRCP信号，尤其是非节律性的动作（如上肢的运动），实验采集的MRCP可能为多种运动意图产生的MRCP叠加. 为了保证所提取的MRCP是由手部抓握运动意图产生的，在实验过程中，需要在不依赖视听觉刺激的情况下采集被试完全自主运动时的EEG信号. 通过肌电(electromyogram，EMG)信号检测自主运动的发生时刻. 被试保持平静，思想放松，要避免过度紧张以免脑电信号产生较大的波动，尽量减少吞咽、眨眼运动. 被试自主决定进行抓握运动开始时间，抓握频率不能太高，2次运动间隔时间至少为4 000 ms. 每组实验结束后，被试休息一定时间，进入下组实验.

本系统的脑电电极设定遵从标准的10−20国际电极布局系统. 在实验中，所有电极的阻抗都保持在10 kΩ以下，肌电信号用于检测手部抓握运动的发生. EEG信号采集设备(奥地利NeuroScan，主要包括64导Quik-Cap脑电帽、脑电信号放大器型号SynAmps² 、控制盒、电源等设备)和EMG信号采集设备(美国Noraxon TeleMyo System-1400A)的采样频率均为1 000 Hz.

实验共采集5名志愿者(年龄24±2岁)的脑肌电信号，每名被试实验前均未接受过手部抓握动作的训练. 针对手部抓握运动进行运动意图的预测，实际发生的动作为手掌的抓握、舒展.

脑电信号极其微弱，掺杂着大量的噪声、眼电、工频干扰等伪迹信号. 在离线数据处理前需要预处理采集的脑电信号. 主要预处理流程包括：重参考、伪迹去除、基线校正、50 Hz工频滤波、0.1~30 Hz带通滤波. 预处理采用Matlab中的eeglab工具箱完成，在采集EEG信号实验中以GND接地电极作为参考，采集头皮上信号和大地的差值，同时记录身体其他部位产生的电位信号. 为了提高信号质量，将参考电极转换为左右乳突M1、M2电极进行重参考.

实验需要采集运动前的脑电信号，用于提取运动准备电位，进行被试的手部抓握运动意图的分析预测. 实际抓握运动发生时刻前的1 500 ms到抓握运动的整个时间段的EEG信号作为运动准备时期，此时被试大脑已经处于准备状态. 实际运动发生时刻前的3 000 ms到前1 500 ms为平静时期，此时被试大脑处于静息状态，无任何运动意图的产生. 须确定运动的起始时刻t作为0 ms时刻. 检测指浅屈肌的EMG信号，当肌电信号超过手部抓握最大肌电信号幅值A阈值10%时，便认为肌肉已经被激活，实际运动已经发生. 其中EEG信号与EMG信号的采集处理流程图如图2所示.

现有研究尚未明确大脑运动意图的产生机理，相关大脑的激活区域也未明确. 采集的脑电信号所使用的通道数量，直接影响采集数据大小.若选取过多的通道进行离线数据分析，不仅包含大量冗余信息，而且会导致算法运行速度过慢，不符合BCI系统要求，因此须确定激活脑区. 采集被试全通道EEG信号进行平均叠加，分析时序关系和各个时刻激活的相关区域，观察被试者从无任何运动意图平静时期到运动准备状态，再到运动执行时期整个过程内，不同时刻大脑激活程度的变化.

如图3所示为手部抓握运动脑地形图及通道位置.从−3 000 ms开始绘制脑地形图，每图间隔200 ms，观察大脑的激活脑区，在−3 000 ~−2 200 ms时，EEG信号的幅值未出现显著变化，此时大脑处于静息状态；从−2 000 ms开始，可观察到EEG信号的幅值出现负向偏移，大脑脑区开始激活. 从脑地形图可以看出FC5、FC3、 FC1、FC2、 FC4、 FC6、 C5、 C3、C1、 C2、 C4、C6、 CP5、 CP3、 CP1、CP2、 CP4、 CP6这18个通道在自主运动时，有较大的负向偏移振幅，这些活动反映运动意图的预测过程.

这18个通道所在脑区(后顶叶和额叶皮质区)与运动规划和运动意图关联较大. 同时，内侧前额叶皮质区是重要的功能区，对目标任务的执行起着重要作用，反映运动决策的过程.

对于长度为N的时间序列 $\boldsymbol{x}=\left[x\right(1), x(2),\cdots , $ $ x(N\left)\right]$，计算线性特征幅值A的均值：

(1) $ \stackrel-{A}=\frac{1}{N}\sum _{i=1}^{N}{x}{\left(i\right)}. $

式中： $ x\left(i\right) $为时间序列。

计算大脑信号系统的样本熵复杂度之前需要对大脑系统进行相空间重构. 相空间重构是通过时间延迟嵌入法，获得与脑电信号具有相同拓扑性质的几何对象，揭示可能蕴含在原序列的信息. 对于N个数据组成的一维脑电信号时间序列 $ \boldsymbol{x} $，通过相空间重构为 $ M\times m $的矩阵空间，其中M为m维嵌入后的点数， $ M=N-(m-1)t $. 重构为M个m维的时间序列，其轨迹矩阵X可表示为

(2) $ \boldsymbol{X}=\left[\begin{array}{c}{\boldsymbol{X}}_{1}\\ \begin{array}{c}{\boldsymbol{X}}_{2}\\ \begin{array}{c}\vdots \\ {\boldsymbol{X}}_{\mathit{M}}\end{array}\end{array}\end{array}\right]\left[ {\begin{array}{*{20}{c}} {{{x}}_1}&{x_{1 + \tau }}&\cdots &{x_{1 + \left( {m - 1} \right)\tau }}\\ {{{x}}_2}&{x_{2 + \tau }}&\cdots &{x_{2 + \left( {m - 1} \right)\tau }}\\ \vdots &\vdots &\; &\vdots\\ {x}_{M}&{x}_{M+\tau }&\cdots &{x}_M+\left(m-1\right)\tau \end{array}} \right]. $

式中：m为嵌入维数，τ为时间延迟.

m、τ值的确定，在本研究中，采用计算量小、计算准确的C-C算法^[14-15]. 算法流程如下。1)定义相关积分函数 $ C\left(m,N,r,t\right) $，其中N为EEG信号时间序列长度，t是时间尺度，r为半径. 2)将时间序列分解为t个不相交子序列，当 $ N\to \infty $，计算统计量 $ S\left(m,r,t\right) $. 3)定义△ $ S\left(m,r,t\right) $为最大 $ S\left(m,r,t\right) $和最小 $ S\left(m,r,t\right) $的差量，S_cor (t) 为△S (t) 与 $ \stackrel{-}{S}\left(t\right) $绝对值之和. 4)基于统计结论，计算 $ \stackrel{-}{S}\left(t\right) $，△ $ \stackrel{-}{S}\left(t\right) $， $ {S}_{{\rm{cor}}}\left(t\right) $，确定τ，从而确定m.

样本熵作为时间序列规律性度量^[16]，衡量时间序列中产生新模式概率大小. 产生新模式概率越大，序列越复杂^[17]. 样本熵计算流程如下.

1)相空间重构m维空间向量：

(3) $ \left.\begin{array}{c} {\boldsymbol{X}}_{\mathit{m}}\left(\mathit{i}\right)=\left[x\left(i\right),x\left(i+1\right),\cdots ,x\left(i+m-1\right)\right]; \\ 1\leqslant i\leqslant N-m+1. \end{array} \right\}$

式中： $ {\mathit{X}}_{\mathit{m}}\left(\mathit{i}\right) $表示第i个点相空间位置.

2)定义向量 $ \boldsymbol{X}\left(\mathit{i}\right) $与 $ \boldsymbol{X}\left(\mathit{j}\right) $的距离 $d\left[{\boldsymbol{X}}_{\mathit{m}}\left({i}\right),{\boldsymbol{X}}_{\mathit{m}}\left({j}\right)\right]$为两者对应元素的最大差值绝对值. 即

(4) $ \left.\begin{array}{c} d\left[\left.{\boldsymbol{X}}_{\mathit{m}}\left({i}\right),{\boldsymbol{X}}_{\mathit{m}}({j}\right)\right]=\mathrm{max}\;\left|x\left(i+k\right)-x\left(j+k\right)\right|;\\ k=\mathrm{0,1},\cdots ,m-1,\;i\ne j. \end{array} \right\} $

3)对于给定的 $ {\boldsymbol{X}}_{\mathit{m}}\left(\mathit{i}\right),{统}{计} $d[X_m(i), X_m(j)]≤r的 $ j(1\leqslant j\leqslant N-m,j\ne i) $的数目，记做B_i.对于1≤i≤N−m，定义：

(5) $ {B}_{i}^{m}\left(r\right)={B}_{i}/{(N-m)}. $

式中： $ {B}_{i}^{m}\left(r\right) $为第i个点与其他点的距离小于阈值r的个数占所有点数的比例.

4)定义所有点数比例的平均值为

(6) $ {B}^{m}={(N-m+1)}^{-1}{\sum }_{i=1}^{N-m}{B}_{i}^{m}. $

5)增加维数到m+1，当N有限值时，计算样本熵：

(7) $ S=-\mathrm{ln}\;\left[{{B}^{m+1}\left(r\right)}/{{B}^{m}\left(r\right)}\right]. $

从式 (7) 可以看到，样本熵的值依赖于r的取值，其中 $r=(0.10～0.25)D$， $ D $为 $ \left\{x\right(i\left)\right\} $的标准差^[18]，样本熵算法流程图如图4所示.

单一特征并不能有效揭示大脑在产生运动意图时潜在的动力学行为，为此分别从幅值及复杂度提取原始EEG信号特征，实现特征互补. 根据脑地形图结果选取18个通道EEG信号数据提取特征，将2个18维特征矩阵线性拼接成36维的特征矩阵进行特征融合，输入分类器中. 其中特征向量表达式为

$ {\boldsymbol{F}}={[{\boldsymbol{F}}}_{{\rm{AC1}}},\cdots,{{\boldsymbol{F}}}_{{\rm{AC6}},}{{\boldsymbol{F}}}_{{\rm{ACP1}}},\cdots,{{\boldsymbol{F}}}_{{\rm{ACP6}},}{{\boldsymbol{F}}}_{{\rm{AFC1}}},\cdots,{{\boldsymbol{F}}}_{{\rm{AFC6}}}, $

(8) $ {{\boldsymbol{F}}}_{{\rm{SC1}}},\cdots,{{\boldsymbol{F}}}_{{\rm{SC6}},}{{\boldsymbol{F}}}_{{\rm{SCP1}}},\cdots,{{\boldsymbol{F}}}_{{\rm{SCP6}},}{{\boldsymbol{F}}}_{{\rm{SFC1}}},\cdots,{{\boldsymbol{F}}}_{{\rm{SFC6}}}] . $

式中：F为融合后特征，各元素为单个通道特征向量，A、S分别为幅值特征、复杂度特征，C、CP、FC分别为3种通道。

本研究对手部抓握运动意图的整个过程进行时频分析，分析其不同时刻t下不同频段f的能量密度，在时频域上研究平静时期与运动准备时期之间存在的功率谱差异. 计算整个过程的事件相关谱扰动(event-related spectral perturbation, ERSP)，测量所采集不同通道、不同信号成分的平均功率谱变化.

通过事件相关频谱扰动图，在较宽频率范围内实现对随时间变化的平均事件相关频谱功率变化的可视化，以C1通道为例，结果如图5所示，其他通道均表现出相似的变化模式. 分析表明当产生不同运动意图时，脑电信号在时域及时频域上存在可区分的特征，可通过时域频域特征识别运动意图. 在初始平静状态时，大脑未产生任何运动意图大脑平均功率谱稳定. 当0 ms时刻，产生被试实际的手部抓握运动时，大脑对侧脑区EEG会波幅变大，伴随着频带能量的上升，出现事件相关同步化(event related synchronization，ERS). 相反，大脑同侧脑区EEG波幅减小，被称为时间相关去同步(event-related desynchronization，ERD)^[19]. 在−1 500 ms时刻左右，被试产生运动意图，处于运动准备状态，大部分大脑区域被激活，激活程度有差异，主要集中脑区为额叶皮质区. 与−3 000~−1 500 ms阶段相比，此时EEG信号出现负向偏移，功率谱密度降低，出现ERD现象，差异主要集中在0~13 Hz低频信号内.

利用C-C算法计算5个被试在不同通道(C3、C4)下，运动准备时期的EEG信号相空间重构的所有实验次数T嵌入维数及时间延迟，以被试Ⅱ的分析结果为例，结果如图6所示. 同时根据其他被试的结果可知，嵌入维数的取值为[2, 5]，时间延迟的取值为[4, 10]. 选择的相空间m=2，τ=5.

不同大脑状态下脑电信号相空间重构轨迹如图7所示。由图可知，在相空间中存在相对稳定、收敛的极限环. 重构向量包含与原来系统具有相同拓扑性质的等价吸引子，大脑系统的所有可能状态都对应空间中的独立点. 具体表现为局部不稳定，但整体上的规律性^[20]. 不同大脑运动意图在各个维度上的延伸程度以及总体的聚集程度各不相同，利用样本熵计算运动轨迹的复杂度以获取不同运动意图的差异信息.

EEG信号不同节律包含不同大脑活动的动力学信息，与神经元的活动模式密切相关^[21]. 因此，通过非线性方法分析得到整个运动产生过程的不同频段脑电信号样本熵S在不同时刻t变化, 如图8所示. 由图可知，当大脑产生运动意图时Delta波段(0.1~4.0 Hz)、Theta波段(4.0~8.0 Hz)、Alpha波段(8.0~13.0 Hz)大脑信号复杂度出现明显的下降，与频带能量变化出现相反的变化，Beta (13.0~30.0 Hz)信号未显示出明显的变化.

大脑运动区域可以建模为由状态层次组成的混沌系统. 根据大脑的活动状况（例如，放松，睁眼，闭眼等），系统会跳至特定状态^[22]. 在产生运动意图期间，信号诱发振荡性爆发，从低能量混沌状态到高能量状态的转变. 对于每个运动意图，引起的振荡的幅度是不同的. 当被试处于静息时，脑电信号会相对无序. 当产生手部抓握运动意图时，被试集中精神，神经元放电会变得有序，导致样本熵的降低. 当产生实际运动时，大脑状态产生明显的振荡，随之而来是大脑信号复杂度的上升. 因此样本熵混沌特征变化可以呈现不同运动意图的特性. 基于样本熵的非线性分析可以研究大脑自主运动时的时空演化过程及神经机制^[23].

通过提取EEG信号幅值与神经复杂度特征实现对于平静时期、准备时期的大脑状态识别，预测手部运动意图. 2种大脑状态EEG信号采集时被试均无任何实际运动产生.

提取实际运动发生时刻−3 000 ~−1 500 ms时间段的EEG信号作为休息时期的运动意图，运动发生时刻−1 500 ~0 ms作为运动准备时期的运动意图，分别对2个时间段EEG信号提取特征进行离线分类识别. 输入2种不同的分类器1）线性判别分析分类器(linear discriminant analysis, LDA) 2）支持向量机(SVM)，得到不同被试的运动意图识别结果. LDA原理是通过寻找最优投影，将高维特征向量映射到一维，实现数据集分类的分离. SVM的原理是使分离的超平面之间的距离最大. SVM具有良好的泛化能力，并能选择最佳超平面作为决策边界. 经PCA降维后2种分类器预测结果如图9所示. 其中，被试Ⅱ意图预测的准确率p最高为88.46%，整体取得良好的分类准确率，基于LDA的分类识别p=80.75%±3.28%, 基于SVM分类识别p=77.37%±7.80%.

与其他研究团队的实验数据的离线分析结果对比如表1所示.基于特征融合的手部运动意图预测单通道数据运行时间为130 ms，总共18通道即一次的算法运行提取特征所需时间为2.34 s. 每分钟的假阳性(false positive, FP)检测为0.22±0.11 FP_S/min.

由于实验被试、动作、过程等因素的影响导致实验范式各不相同，同时关于算法运行时间延迟的结果较少，算法执行时间无法进行准确比较。通过比较不同研究中的假阳性率，希望借此说明本文算法的可行性.

为了验证实时控制设备的可行性，使用时间窗模拟实时发送数据流进行伪在线实验，以用于后续实时在线控制设备的实现. 基于实验采集的被试脑电信号数据进行伪在线实验分析，为了防止误判，采用投票机制， 以100 ms不重叠的滑动窗口从实际运动发生时刻的−2 000 ms开始进行运动意图的判别，只有当10个连续检测的窗口中有5个及5个以上的窗口判别存在意图时，才认为检测到被试发生了运动意图.

通过不同分类器提前预测不同被试手部运动意图，预测时间如表2所示. 由表可知，基于LDA可在实际运动发生时刻（−1455±312） ms检测到被试的手部抓握运动意图平均预测准确率为82.45%，基于SVM可在实际运动发生时刻的（−1608±270） ms检测到被试的手部抓握运动意图，平均预测准确率为79.4%.

不同被试预测准确率和提前预测到运动意图的时间存在差异，本研究认为这是不同被试存在个体差异导致的. 被试Ⅱ分类准确率最高，可能是实验中该被试大脑参与度最高，MRCP的幅值和神经复杂度的变化相对明显. 预测时间提前量是5个被试中最小的，即从运动意图产生到手部实际抓握运动发生的时间间隔最短，该被试的MRCP在幅值与神经复杂度的变化相对明显且快速.

针对现有BCI系统实时性低，输出存在较大时间滞后的问题， 基于大脑存在混沌现象，利用非线性方法表征在不同大脑状态下，神经元的活跃程度差异. 运动准备阶段的EEG信号神经复杂度明显低于大脑进行运动想象阶段的复杂度，基于此对手部运动意图的实现进行提前预测，提供解决BCI实时性低问题的可行性策略.限于EEG信号过于微弱、易受干扰，本研究提出实验范式的控制任务数量有限，不能满足实际BCI应用中多分类的场景. 在未来的研究中，可以研究引入如速度、方向的更多运动参数对EEG信号复杂度的影响，以扩展BCI的输出控制命令种类.