运动控制的关键在于运动皮层和运动所涉及的肌肉之间的相互作用^[1]. 在肌肉持续收缩中，皮层肌肉功能耦合(functional cortical muscular coupling, FCMC)突出表现在alpha频段(8~13 Hz)上^[2]；在肌肉维持恒力输出时，FCMC主要表现在beta频段(13~35 Hz)上^[3]；在肌肉持续输出较强力量时，gamma频段(35~60 Hz)上的FCMC耦合强度最大^[4]. 这些数据表明，不同运动系统的FCMC在不同频段上起着不同的作用.

FCMC分析主要采用相干性方法. 这种方法仅能描述线性耦合关系，不能描述脑肌电信号间存在的非线性耦合关系. 为了更准确地理解FCMC，Schreiber^[5]提出传递熵(transfer entropy, TE)分析方法. TE不需要假设任何模型，但在实际情况中没有考虑瞬时效应，这会影响因果关系的度量^[6]. 为了解决上述问题，Faes等^[7]基于有效条件熵估计器，并遵循用于多元时间序列的非均匀嵌入的顺序过程，提出补偿传递熵(compensated TE, cTE)方法. 其中嵌入向量满足应用于系统当前状态的熵的最小化准则，随着嵌入维数的增加来补偿熵估计偏差，保证了结束时的最小熵率. 谢平等^[8]将TE与变分模态分解(variational mode decomposition, VMD)结合应用于脑肌间耦合分析，马鹏刚等^[9]将TE与多元经验模态分解(multivariate empirical mode decomposition, MEMD)结合构建MEMD-TE分析模型. VMD、MEMD是经验模态分解(empirical mode decomposition, EMD)的衍生. EMD适用于非线性和非平稳序列，但经过该方法分解后，信号存在边界效应及模态混叠现象. MEMD可以较好地消除模态混叠现象^[10]，但在捕获多变量数据动态方面有更好的方法. Yang等^[11]引入自适应投影多元经验模态分解方法(adaptive-projection intrinsically transformed MEMD, APITMEMD)，可以解决多通道脑肌电数据中的相关性和功率不平衡问题. APITMEMD具有MEMD的缓解模式混合的优点，还能产生更少的本征模态函数(intrinsic mode function, IMF)分量.

本研究提出新的多尺度补偿传递熵(multiscale cTE, McTE)方法, 用于脑肌电不同尺度上的耦合分析. 通过研究恒力输出上肢运动过程中不同尺度间的神经肌肉耦合，定量描述脑肌电信号间的耦合特征和信息传递方向，深入理解运动控制机制，为探索康复运动功能评价提供理论依据.

TE可以度量2个生物系统间的信息流向和转移量. 考虑2个表示物理系统 $ x、y $随时间变化的随机过程 $ {\boldsymbol{X}}、{\boldsymbol{Y}} $，令 ${x_t}$、 ${y_t}$为时间 $t$时访问 ${\boldsymbol{X}}$和 $ \boldsymbol{Y} $的随机变量， ${{\boldsymbol{x}}_{n:t}}$为从时间 $n$到时间 $t$的所有 $x$样本组成的向量，则定义 ${\boldsymbol{X}}$到 $ \boldsymbol{Y} $的TE为 $ {T_{{\boldsymbol{X}} \to {\boldsymbol{Y}}}} $，其中“ $ \to $”代表传递方向.

(1) $ T_{{\boldsymbol{X}} \rightarrow {\boldsymbol{Y}}} = \sum p\left(y_{t}, \boldsymbol{y}_{1: t-1}, \boldsymbol{x}_{1: t-1}\right) \ln \frac{p\left(y_{t}, \boldsymbol{y}_{1: t-1}, \boldsymbol{x}_{1: t-1}\right) p\left(\boldsymbol{y}_{1: t-1}\right)}{p\left(y_{t}, \boldsymbol{y}_{1: t-1}\right) p\left(\boldsymbol{y}_{1: t-1}, \boldsymbol{x}_{1: t-1}\right)}. $

式中： $p( \cdot )$为联合概率。根据条件熵的定义，可以得出传递熵与条件熵(conditional entropy，CE) $H( \cdot ) $的等价关系：

(2) $ T_{{\boldsymbol{X }}\rightarrow {\boldsymbol{Y}}}=H\left(y_{t} \mid \boldsymbol{y}_{1: t-1}\right)-H\left(y_{t} \mid \boldsymbol{x}_{1: t-1}, \boldsymbol{y}_{1: t-1}\right). $

式(2)中定义的TE没有考虑瞬时效应。在过程 $ \boldsymbol{X}$具有内部存储器结构(即 $x_{t}$可由 $\boldsymbol{x}_{1: t-1} $解释)并且与过程 $ \boldsymbol{Y}$ 瞬时相关(即 $y_{t} $ 可由 $x_{t} $ 解释)的情况下，即使 $ \boldsymbol{x}_{1: t-1}$ 不能解释 $y_{t} $ ， 式(2)中得到的TE也有意义. 这表示即使所测的2个系统不耦合，也能够由式(2)得到 ${\boldsymbol{X}}$、 $ \boldsymbol{Y} $的信息传输. 针对这个问题，在TE计算的2项CE中考虑瞬时效应，定义 $ {\boldsymbol{X}} $到 $ {\boldsymbol{Y}} $的补偿TE为 $ {F_{{\boldsymbol{X}} \to{\boldsymbol{ Y}}}} $.

(3) $ \begin{gathered} {F_{{\boldsymbol{X}} \to {\boldsymbol{Y}}}} = H({y_t}|{x_t},{{\boldsymbol{y}}_{1:t - 1}}) - {\text{ }}H({y_t}|{{\boldsymbol{x}}_{1:t}},{{\boldsymbol{y}}_{1:t - 1}}). \hfill \\ \end{gathered} $

构建时间序列EEG信号 ${\boldsymbol{X}} = [{x_1},{x_2}, \cdots {x_i}, \cdots ,{x_M}]$和EMG信号 ${\boldsymbol{Y}} = [{y_1},{y_2}, \cdots {y_i}, \cdots ,{y_M}]$. 假设多元信号 $ {\boldsymbol{x}}(t) $的协方差矩阵为 ${\boldsymbol{C}} = E({{\boldsymbol{x}}^{\rm T}}(t){\boldsymbol{x}}(t))$，信号的第一主成分的方向由协方差矩阵的特征分解确定，公式为

(4) $ {\boldsymbol{C}} = {\boldsymbol{V\varLambda }}{{\boldsymbol{V}}^{\rm T}}. $

式中： ${\boldsymbol{\varLambda }} $为特征值矩阵， ${\boldsymbol{\varLambda }} = {\text{diag}}\left[ {{\lambda _1},{\lambda _2}, \cdot \cdot \cdot ,{\lambda _n}} \right]$； ${\boldsymbol{V}} $为特征向量矩阵， ${\boldsymbol{V}} = \left[ {{{\boldsymbol{v}}_1},{{\boldsymbol{v}}_2}, \cdots, {{\boldsymbol{v}}_n}} \right]$；最大特征值 ${\lambda _1}$所对应的特征向量 ${{\boldsymbol{v}}_1}$规定为第一主成分方向，规定向量 ${{\boldsymbol{v}}_{{\rm{o}}1}}$为 ${{\boldsymbol{v}}_1}$反方向的向量，即 ${{\boldsymbol{v}}_{{\rm{o}}1}} = - {{\boldsymbol{v}}_1}$；用 ${{\boldsymbol{v}}_1}$、 ${{\boldsymbol{v}}_{{\rm{o}}1}}$重新定位先前由均匀投影方案生成的所有方向向量，公式为

(5) $ {\boldsymbol{\hat x}}_{{{\boldsymbol{v}}_1}}^k = \frac{{{\boldsymbol{x}}_{{{\boldsymbol{v}}_1}}^k + \alpha {{\boldsymbol{v}}_1}}}{{|{\boldsymbol{x}}_{{{\boldsymbol{v}}_1}}^k + \alpha {{\boldsymbol{v}}_1}|}}, $

(6) $ {\boldsymbol{\hat x}}_{{{\boldsymbol{v}}_{{\rm{o}}1}}}^k = \frac{{{\boldsymbol{x}}_{{{\boldsymbol{v}}_{{\rm{o}}1}}}^k + \alpha {{\boldsymbol{v}}_{{\rm{o}}1}}}}{{|{\boldsymbol{x}}_{{{\boldsymbol{v}}_{{\rm{o}}1}}}^k + \alpha {{\boldsymbol{v}}_{{\rm{o}}1}}|}}. $

式中： ${{\boldsymbol{\hat x}}^k}$为第 $k$( $k = 1,2, \cdots ,K$)个重置的方向向量； ${{\boldsymbol{x}}^k}$为均匀投影向量； $ \alpha \in [0,1] $由多元信号的功率不平衡程度决定， $\alpha = 1$表示各个通道之间的高功率不平衡. APITMEMD可以很好地处理信号在空间的重要性采样. 在迭代过程中，沿着自适应投影向量迭代 $ {\boldsymbol{x}}(t) $，并通过MEMD获得局部均值，具体APITMEMD算法步骤参见文献[12].

通过APITMEMD方法分别分解EEG信号和EMG信号得到 $K$个IMF分量. 由式(3)准确描述EEG、EMG信号各个尺度上的耦合强度和信息传递特征.

构造EEG信号第 ${k_1}$个IMF分量 ${{\boldsymbol{x}}^{{k_1}}}$到EMG信号第 ${k_2}$个IMF分量 ${{\boldsymbol{Y}}^{{k_2}}}$的多尺度补偿传递熵 $Z{}_{{\rm E}{\rm EG} \to {\rm E}{\rm MG}}$，公式为

(7) $ \begin{split} Z{}_{{\rm E}{\rm EG} \to {\rm E}{\rm MG}} = H(y_t^{{k_2}}|x_t^{{k_1}},{\boldsymbol{y}}_{1:t - 1}^{{k_2}}) - H(y_t^{{k_2}}|{\boldsymbol{x}}_{1:t}^{{k_1}},{\boldsymbol{y}}_{1:t - 1}^{{k_2}}). \hfill \\ \end{split} $

$Z{}_{{\rm E}{\rm EG} \to {\rm E}{\rm MG}}$的生理意义为脑皮层与肌肉组织在某尺度上从脑电到肌电方向的信息传递量和耦合强度. 同理，肌电信号第 ${k_2}$个IMF分量 ${{\boldsymbol{Y}}^{{k_2}}}$到脑电信号第 ${k_1}$个IMF分量 ${{\boldsymbol{X}}^{{k_1}}}$的多尺度补偿传递熵 $Z{}_{{\rm E}{\rm MG} \to {\rm E}{\rm EG}}$公式为

(8) $ \begin{split} Z{}_{{\rm E}{\rm MG} \to {\rm E}{\rm EG}} = H(x_t^{{k_1}}|y_t^{{k_2}},{\boldsymbol{x}}_{1:t - 1}^{{k_1}}) - H(x_t^{{k_1}}|{\boldsymbol{y}}_{1:t}^{{k_2}},{\boldsymbol{x}}_{1:t - 1}^{{k_1}}). \hfill \\ \end{split} $

$Z{}_{{\rm E}{\rm MG} \to {\rm E}{\rm EG}}$表示EMG的 ${k_2}$分量到EEG的 ${k_1}$分量的多尺度补偿传递熵值. $Z$的值越大，表明在此尺度上EEG信号和EMG信号间的耦合强度越强；反之亦然.

本实验数据包含6名成年人，其中5名受试者均无上肢运动功能障碍历史，被依次标记为S1~S5，1名受试者四肢有轻微运动功能障碍，被标记为S6. 所有受试者均为右利手且精神状态良好，均被告知实验详情，由本人签署知情同意书，自愿参加此次数据采集实验. 为了防止运动疲劳对实验采集数据的影响，受试者前1 d未进行剧烈运动. 实验数据采集装置为BrainAmp DC系统，同步采集32路脑电信号、左右手肱二头肌(biceps brachii，BB)和尺侧腕屈肌(flexor carpi ulnaris，FCU)的肌电信号，上肢抓握实验环境如图1所示. 在采集实验前，用酒精清洁受试者被测部位的皮肤，以减少阻抗. 采样时间为5 s，采样频率为1 000 Hz，每次抓握动作后休息15 s，S1~S5完成5、10、20 kg的左右手抓握各5次，S6由于身体原因，只能完成左右手5、10 kg的抓握.

使用APITMEMD和MEMD算法对受试者抓握运动时采集的信号进行分解并对比. 如图2所示为受试者S6的3个脑肌电通道信号的2种算法分解结果. 图中，V为信号振幅，N为点数，C4为脑中央通道. 可见， MEMD分解产生的IMF分量个数相同；APITMEMD分解产生的IMF分量个数也相同，但比MEMD分解产生的少.

如图3所示为APITMEMD方法对受试者S1肌电信号BB的时频分析结果. 图中， $f$为频率，A为幅度. 可见，该方法分解出的IMF频率分量是从高到低排列的，有效缓解了模式混叠.

受试者S4脑肌电信号经过APITMEMD分解得到的各IMF分量的频率范围如表1所示. 表中脑电信号C4与肌电信号BB的IMF降序排列(IMF1~IMF11所对应的带宽范围由高频到低频)，且IMF数量相同. 依据不同特征频段的频率范围，将分解得到的IMF依次划分，为分析脑肌电信号不同特征频段的耦合特性做准备.

为了对比受试者在 ${{\rm{EEG}}} \to {\rm{EMG}}$、 ${{\rm{EMG}}} \to {\rm{EEG}}$方向上各尺度的皮层肌肉耦合的差异，分别计算受试者左、右手不同力度下的McTE和MEMD-cTE的耦合强度并进行对比. 如图4 、5所示分别为受试者S3、S6运动时脑肌电耦合2种方法的对比图 (其他受试者结果与此类似). 图中，k为不同尺度频带；Z为耦合强度；标题中R代表右手，L代表左手；5、10、20代表受试者使用的抓握力度分别为5、10、20 kg. 可知，受试者的皮层肌肉功能耦合分为 ${{\rm{EEG}}} \to {\rm{EMG}}$(下行)方向和 ${{\rm{EMG}}} \to {\rm{EEG}}$(上行)方向，并且皮层肌肉间的耦合情况具有相似性. 同一受试者在不同的握力情况下，其EEG和EMG信号经过APITMEMD分解后的IMF分量个数有所差异. 第1个尺度(图中的大方框处)代表的是高gamma频段(50~75 Hz). 在这个尺度上无论是左手还是右手，无论握力是大是小，始终存在 ${\text{EEG}} \to {\text{EMG}}$的耦合强度大于 ${\text{EMG}} \to {\text{EEG}}$的耦合强度. 这体现了在运动过程中，信息传递有不同的传递通道，运动皮层需要传递更多的信息至肌肉进行控制. 在高gamma频段该方法结果曲线和MEMD-cTE曲线相似，但是在beta频段(图中的小方框处) 2种分解方法的分解结果具有一定差异性.

为了进一步研究受试者在 ${{\rm{EEG}}} \to {\rm{EMG}}$、 ${{\rm{EMG}}} \to {\rm{EEG}}$方向上各个尺度间的耦合特征，受试者S2、S6在不同运动情况下的各尺度间脑肌耦合McTE分别如图6、7所示. 图中，p为脑电信号IMF分量，q为肌电信号IMF分量. 可见，无论是左手还是右手，无论握力是大是小，beta频段(对应于图中脑肌电信号第6~8个IMF分量)的耦合强度最大，这与图4、5中第6~8个尺度耦合强度最大的结果相对应. 该结果与已有的静态握力输出实验模式下EEG-EMG耦合特征研究相符合^[13]，共同证实了在维持恒力输出下，脑肌电耦合主要表现在beta频段，运动皮层与肌肉间beta频段的振荡起主要作用，体现了皮层与肌肉间的信息交互.

研究静态握力下的EEG和EMG信号，深入理解中枢神经系统的控制机制. 皮层肌肉间的2种不同的分解算法：MEMD、 APITMEMD，都可以将脑肌电信号分解成不同的IMF分量，缓解模式混合. APITMEMD方法能产生比MEMD少的IMF分量，这与Yang等^[11]的研究结果一致. 针对皮层肌肉间信息流向的研究表明，脑电和肌电间的耦合是双向的(上行和下行方向)，体现了运动过程中皮层与相应肌肉在互相传递信息. 与此同时，本研究的实验结果进一步表明，下行方向beta频段的耦合强度在皮层肌肉耦合中占据主导地位，且在高gamma频段下行方向的耦合强度始终大于上行方向的耦合强度，这与Chen等^[14]得出的结论相一致. 受试者在不同尺度上的实验结果表明，在左、右手不同握力的情况下，beta频段的耦合特征最为明显，再次验证了受试者在不同方向上各尺度皮层肌肉耦合差异实验的结果.

针对传统耦合方法不能准确描述脑肌电信号在不同尺度上的耦合特征的问题，提出多尺度补偿传递熵分析方法. 通过实验表明，基于APITMEMD的生理信号分解比MEMD分解产生的IMF分量少；采用本研究所提方法，能够比MEMD-cTE更准确地描述脑肌电信号在beta频段上的耦合特征. 对于人体运动这一复杂过程，本研究选取的脑肌电信号通道数量、受试者数量均有限，且仅探究了脑肌电信号间的耦合特征，不具有全面性. 下一步计划招募更多的受试者，采集更多的脑肌电通道信号，从脑电耦合、肌间耦合与脑肌电耦合多层次深入研究人体运动的控制机制.