基于小波‑传递熵的多时频尺度间脑肌电耦合分析方法

摘要

本发明公开了一种基于小波‑传递熵的多时频尺度间脑肌电耦合分析方法，所述方法包括脑肌电信号同步采集部分和信号处理部分，脑肌电信号同步采集部分包括脑电信号采集和肌电信号采集；信号处理部分包括信号预处理和脑肌电的小波‑传递熵分析方法。本发明具有可应用性、可采纳性，在康复医学领域具有重要的应用价值。

说明书

技术领域

本发明涉及神经康复工程及运动机制研究领域，具体涉及一种基于小波-传递熵的多时频尺度间脑肌电耦合分析方法。

背景技术

脑电(electroencephalogram,EEG)和肌电(electromyographic,EMG)信号分别包含躯体运动控制信息和肌肉对大脑控制意图的功能响应信息，脑肌电信号间的多尺度耦合信息反映多层次的皮层-肌肉功能耦合(Functional corticomuscular coupling,FCMC)的连接信息。目前，脑肌电同步特征研究主要基于相干分析，获取大脑运动意识驱动与肌肉运动响应之间的功能联系特征，但传统的相干性分析并不能体现耦合方向特征。为更好理解大脑皮层与相应肌肉间的功能交互及信息传递特性，格兰杰因果性分析被应用于脑肌电同步研究，发现脑肌电间存在双向(下行EEG→EMG、上行EMG→EEG)耦合联系。但由于脑肌电之间耦合模型未知且脑肌电信号间的功能耦合存在着非线性因果关系，基于既定模型的格兰杰因果分析方法不能有效描述脑肌电非线性耦合特征。传递熵具有不依赖于既定模型且实现非线性定量分析的特点，能够有效的估计皮层-肌肉之间的功能耦合强度及信息转移方向。因此，传递熵模型对于估计皮层-肌肉之间的功能耦合强度及信息传递方向特征、揭示运动过程中皮层和肌肉间的运动控制及响应机制具有可行性。2015年，本发明作者在《基于多尺度传递熵的脑肌电信号耦合分析》中曾提出多尺度传递熵方法，并基于该方法研究不同时间尺度上 脑肌电信号耦合特征。但随着粗粒化尺度的增加，序列长度减少，可能使得熵估计不准确。移动均值化的提出克服了这一弊端，使得各个尺度的时间序列长度保持相同。但以上研究还存在一些不足：粗粒化与移动均值化方法仅是将脑电和肌电信号进行时间尺度化，不能够刻画出脑电和肌电的时频域特性以及不同时频尺度间脑肌电信号的非线性耦合及信息传递。

发明内容

本发明的目的在于提供一种能够发现皮层肌肉间非线性相关性、深入探索大脑皮层与肌肉之间的耦合及信息传递特征的基于小波-传递熵的多时频尺度间脑肌电耦合分析方法。

本发明所述方法的步骤如下：

步骤1，采用64导Neuroscan设备同步采集脑电信号和肌电信号；

步骤2，利用Neuroscan设备数据处理软件对采集到的脑电信号和肌电信号分别去除基线漂移、溢出、眼动和工频干扰；

步骤3，采用Daubechies类db4小波基函数对脑电和肌电信号进行频谱分解，分析脑电信号和肌电信号不同时频尺度间的同步特性，定量描述脑肌间非线性耦合和信息传递特征；

步骤4，对脑肌间非线性耦合和信息传递特征进行运动功能分析。

进一步的，在步骤1中，脑电电极采用国际10-20系统标准，以双耳的乳突作为参考，从32导头皮脑电采集设备记录对应运动的脑电信号；使用Synamp2设备采集肌电信号，电极顺着肌纤维方向，放置在肌腹位置处。

进一步的，

所述步骤3的具体方法如下：

基于步骤2中预处理后的实测数据构建脑电信号x(t)和肌电信号y(t)两组时间序列；对脑电信号x(t)进行小波变换，将小波变换应用于尺度划分之中；

首先构造小波函数，公式如下：

Ψ_j,k(t)＝2^-j/2Ψ(2^-jt-k)>

式中，Ψ(t)为母小波；k为Ψ(t)纵坐标方向的平移量，j表示信号的层数，j,k∈Z，Z为整数集；尺度参数为2^j，平移参数2^jk；t为时间指标；

然后对脑电信号x(t)进行7层频谱分解，得到小波变换系数

$C_{j,k}={\int }_{-\infty }^{\infty }x\left(t\right){\Psi }_{j,k}^{*}\left(t\right)dt---\left(2\right)$

小波系数C_j,k按照频率范围从高到低的顺序排列，提取第3、4、5、6、7层系数重构出gamma(32～64Hz)、beta(16～32Hz)、alpha(8～16Hz)、theta(4～8Hz)和delta(1～4Hz)频段的信号：

$x_i\left(t\right)=C\underset{-\infty }{\overset{\infty }{\Sigma }}\underset{-\infty }{\overset{\infty }{\Sigma }}{C}_{j_i,k}{\Psi }_{j,k}\left(t\right),(i=1,2,3,4,5)---\left(3\right)$

式中，为滤波器的品质因数；

重构信号x₁(t)、x₂(t)、x₃(t)、x₄(t)及x₅(t)分别对应脑电功能频带delta、theta、alpha、beta和gamma的信号；

对于肌电信号y(t)进行上述相同的小波变换过程，得到y_m(t)(m＝1,2,3,4,5)，分别对应肌电delta、theta、alpha、beta和gamma频段的信号；

基于传递熵计算方法，构造x(t)到y(t)的小波-传递熵WTE_x→y，公式如下：

${\mathrm{WTE}}_{x\to y}=\underset{y_{t+u}^m,y_t^m,x_t^i}{\Sigma }p(y_{t+u}^m,y_t^m,x_t^i)log\frac{p(y_{t+u}^m,y_t^m,x_t^i)p\left(y_t^m\right)}{p(y_{t+u}^m,y_t^m)p(y_t^m,x_t^i)}---\left(4\right)$

式中，u为预测时间；p(·)为变量之间的联合概率；分别表示脑电和肌电delta、theta、alpha、beta、gamma分量的延迟向量；为的预测序列；

WTE_x→y则表示脑电信号EEG的x_i(t)分量到肌电信号EMG的y_m(t)分量间的传递熵值；同理信号y(t)到x(t)的小波-传递熵WTE_y→x的表达式为：

${\mathrm{WTE}}_{y\to x}=\underset{x_{t+u}^i,x_t^i,y_t^m}{\Sigma }p(x_{t+u}^i,x_t^i,y_t^m)log\frac{p(x_{t+u}^i,x_t^i,y_t^m)p\left(x_t^i\right)}{p(x_{t+u}^i,x_t^i)p(x_t^i,y_t^m)}---\left(5\right)$

式中，为的预测序列；WTE_y→x表示EMG的y_m(t)到分量到EEG的x_i(t)分量间的传递熵值；传递熵值越大，说明在此频段间皮层肌肉耦合越强；反之亦然。

本发明与现有技术相比有益效果如下：本发明方法利用小波-传递熵分析脑肌电信号信息传递特性，定量描述脑电信号与肌电信号频段间非线性的耦合及信息传递特征，有助于探索大脑皮层与肌肉之间的功能联系，研究运动控制反馈机制及运动障碍病理机制，建立基于脑肌电信号的康复状态评价指标，构建康复机器人运动状态及患者生理状态评价机制，可以获得可观的社会效益和经济效益。

附图说明

图1为本发明中Neuroscan设备的结构简图。

图2为本发明方法的工作流程图。

图3为受试者C4通道脑电信号的小波分解后时频结果图。

图4为受试者指浅屈肌处肌电信号的小波分解后时频结果图。

图5为受试者的脑肌电信号小波-传递熵分析结果图。

附图标号：1-脑电电极、2-电极帽、3-脑肌电采集仪、4-肌电导联线、5-肌电电极。

具体实施方式

脑电信号和肌电信号十分微弱，具有非线性、非平稳性和频域特性突出等特点。在运动过程中，神经系统与肌肉之间的交互控制机制可以通过脑肌电信号的同步耦合分析体现。小波分解能够提取出脑肌电信号的特定时频数据段，而且传递熵可以刻画信号间非线性的耦合及信息传递特征，本发明通过研究脑肌电间的小波-传递熵分析，获得不同运动状态下大脑皮层与肌肉之间信息传递关系，进而研究运动功能障碍产生的生理机制。

实施例1：

如图2所示，方法步骤如下：

步骤1，采用64导Neuroscan设备同步采集脑电信号和肌电信号。

Neuroscan设备的结构如图1所示，由脑电电极、电极帽、脑肌电采集仪、肌电导联线、肌电电极连接组成。

脑电信号采集：脑电电极采用国际标准10-20电极放置标准，通过电极帽2实现脑电电极1与头皮接触。在手部静态握力输出运动下进行脑肌电信号同步采集实验。M1、M2导联分别连接到左右耳后乳突作为参考电极，接地电极布置在头顶正中，从32导头皮脑电采集设备中选择C3、C4以及CPZ区记录对应运动的脑电信号。

肌电信号采集：使用Synamp2设备采集指浅屈肌(flexor digitorum superficialis,FDS)处的肌电信号，首先用酒精擦拭被测部位的皮肤表面，去除皮肤表面油脂和皮屑，然后顺着肌纤维方向将肌电电极5粘贴在肌腹位置处，并将肌电导联线4适当固定尽量减少运动过程中导联线晃动的干扰。

步骤2，利用Neuroscan设备数据处理软件对采集到的脑电信号和肌电信号分别去除基线漂移、溢出、眼动和工频干扰；

步骤3，采用Daubechies类db4小波基函数对脑电和肌电信号进行频谱分解，基于步骤2中预处理后的实测数据构建脑电信号x(t)和肌电信号y(t)两组时间序列。以脑电信号x(t)为例，将小波变换应用于尺度划分之中。首先构造小波函数，公式如下：

Ψ_j,k(t)＝2^-j/2Ψ(2^-jt-k)>

式中，Ψ(t)为母小波；k为Ψ(t)纵坐标方向的平移量，j表示信号的层数，j,k∈Z，Z为整数集；尺度参数为2^j，平移参数2^jk；t为时间指标。

然后对脑电信号x(t)进行7层频谱分解，得到小波变换系数

$C_{j,k}={\int }_{-\infty }^{\infty }x\left(t\right){\Psi }_{j,k}^{*}\left(t\right)dt---\left(2\right)$

基于上式得到的小波系数C_j,k按照频率范围从高到低的顺序排列，提取第3、4、5、6、7层系数重构出gamma(32～64Hz)、beta(16～32Hz)、alpha(8～16Hz)、theta(4～8Hz)和delta(1～4Hz)频段的信号：

$x_i\left(t\right)=C\underset{-\infty }{\overset{\infty }{\Sigma }}\underset{-\infty }{\overset{\infty }{\Sigma }}{C}_{j_i,k}{\Psi }_{j,k}\left(t\right),(i=1,2,3,4,5)---\left(3\right)$

式中，为滤波器的品质因数；

重构信号x₁(t)、x₂(t)、x₃(t)、x₄(t)及x₅(t)分别对应脑电功能频带delta、theta、>

对于肌电信号y(t)进行上述相同的小波变换过程，得到y_m(t)(m＝1,2,3,4,5)，分别对应肌电delta、theta、alpha、beta和gamma频段的信号。

基于传递熵计算方法，构造x(t)到y(t)的小波-传递熵WTE_x→y，公式如下：

${\mathrm{WTE}}_{x\to y}=\underset{y_{t+u}^m,y_t^m,x_t^i}{\Sigma }p(y_{t+u}^m,y_t^m,x_t^i)log\frac{p(y_{t+u}^m,y_t^m,x_t^i)p\left(y_t^m\right)}{p(y_{t+u}^m,y_t^m)p(y_t^m,x_t^i)}---\left(4\right)$

式中，u为预测时间；p(·)为变量之间的联合概率；分别表示脑电和肌电delta、theta、alpha、beta、gamma分量的延迟向量；为的预测序列；

WTE_x→y则表示脑电信号EEG的x_i(t)分量到肌电信号EMG的y_m(t)分量间的传递熵值。同理信号y(t)到x(t)的小波-传递熵WTE_y→x的表达式为：

${\mathrm{WTE}}_{y\to x}=\underset{x_{t+u}^i,x_t^i,y_t^m}{\Sigma }p(x_{t+u}^i,x_t^i,y_t^m)log\frac{p(x_{t+u}^i,x_t^i,y_t^m)p\left(x_t^i\right)}{p(x_{t+u}^i,x_t^i)p(x_t^i,y_t^m)}---\left(5\right)$

式中，为的预测序列；WTE_y→x表示EMG的y_m(t)到分量到EEG的x_i(t)分量间的传递熵值；传递熵值越大，说明在此频段间皮层肌肉耦合越强；反之亦然。

基于上述指标，计算手部静态握力输出运动下，不同耦合方向上、不同时频尺度间的WTE值，即能够定量描述脑电信号EEG和肌电信号EMG多时频尺度间的非线性同步耦合特征。

为验证本发明所述的脑肌电信号小波-传递熵分析方法的可行性和有效性，募集8名健康的受试者进行手部静态握力输出实验，受试者相关信息如表1所示。按照本发明所述的脑肌电采集与分析过程，同步采集受试者恒力输出运动下的脑肌电信号，进行分析并研究受试者运动运动过程中皮层肌肉间的耦合与信息传递机制。

本次实验采集左手指浅屈肌(flexor digitorum superficialis,FDS)处的肌电信 号与对侧C4通道脑电信号，并计算WTE值。

图3和图4分别为受试者脑电、肌电信号经过小波分解后各时频尺度上的时频域结果(左侧是时域图，右侧是频域图)，从中可以看出脑电和肌电信号经过小波分解后可以得到delta、theta、alpha、beta和gamma频段的信号。

图5为8位受试者脑肌间耦合的小波-传递熵分析后的平均值。从中可以看出，在静态握力输出过程中皮层肌肉间beta频段的耦合强度最为显著，而且不同耦合方向上各时频尺度间的耦合强度也存在差异，为探究神经肌肉功能耦合机制提供了理论研究方法。

表1 受试者相关信息。