基于小波?传递熵的多时频尺度间脑肌电耦合分析方法
【专利摘要】本发明公开了一种基于小波?传递熵的多时频尺度间脑肌电耦合分析方法，所述方法包括脑肌电信号同步采集部分和信号处理部分，脑肌电信号同步采集部分包括脑电信号采集和肌电信号采集；信号处理部分包括信号预处理和脑肌电的小波?传递熵分析方法。本发明具有可应用性、可采纳性，在康复医学领域具有重要的应用价值。
【专利说明】
基于小波-传递熵的多时频尺度间脑肌电耦合分析方法
技术领域
[0001] 本发明涉及神经康复工程及运动机制研究领域，具体涉及一种基于小波-传递熵 的多时频尺度间脑肌电耦合分析方法。
【背景技术】
[0002] 脑电（electroencephalogram, EEG)和肌电（electromyographic，EMG)信号分别包 含躯体运动控制信息和肌肉对大脑控制意图的功能响应信息，脑肌电信号间的多尺度耦合 信息反映多层次的皮层-肌肉功能親合(Functional corticomuscular coupling,FCMC)的 连接信息。目前，脑肌电同步特征研究主要基于相干分析，获取大脑运动意识驱动与肌肉运 动响应之间的功能联系特征，但传统的相干性分析并不能体现耦合方向特征。为更好理解 大脑皮层与相应肌肉间的功能交互及信息传递特性，格兰杰因果性分析被应用于脑肌电同 步研究，发现脑肌电间存在双向（下行EEG-EMG、上行EMG-EEG)耦合联系。但由于脑肌电之 间耦合模型未知且脑肌电信号间的功能耦合存在着非线性因果关系，基于既定模型的格兰 杰因果分析方法不能有效描述脑肌电非线性耦合特征。传递熵具有不依赖于既定模型且实 现非线性定量分析的特点，能够有效的估计皮层-肌肉之间的功能耦合强度及信息转移方 向。因此，传递熵模型对于估计皮层-肌肉之间的功能耦合强度及信息传递方向特征、揭示 运动过程中皮层和肌肉间的运动控制及响应机制具有可行性。2015年，本发明作者在《基于 多尺度传递熵的脑肌电信号耦合分析》中曾提出多尺度传递熵方法，并基于该方法研究不 同时间尺度上脑肌电信号耦合特征。但随着粗粒化尺度的增加，序列长度减少，可能使得 熵估计不准确。移动均值化的提出克服了这一弊端，使得各个尺度的时间序列长度保持相 同。但以上研究还存在一些不足:粗粒化与移动均值化方法仅是将脑电和肌电信号进行时 间尺度化，不能够刻画出脑电和肌电的时频域特性以及不同时频尺度间脑肌电信号的非线 性耦合及信息传递。
【发明内容】
[0003] 本发明的目的在于提供一种能够发现皮层肌肉间非线性相关性、深入探索大脑皮 层与肌肉之间的耦合及信息传递特征的基于小波-传递熵的多时频尺度间脑肌电耦合分析 方法。
[0004] 本发明所述方法的步骤如下：
[0005] 步骤1，采用64导Neuroscan设备同步采集脑电信号和肌电信号；
[0006] 步骤2,利用Neuroscan设备数据处理软件对采集到的脑电信号和肌电信号分别去 除基线漂移、溢出、眼动和工频干扰；
[0007] 步骤3,采用Daubechies类db4小波基函数对脑电和肌电信号进行频谱分解，分析 脑电信号和肌电信号不同时频尺度间的同步特性，定量描述脑肌间非线性耦合和信息传递 特征；
[0008] 步骤4,对脑肌间非线性耦合和信息传递特征进行运动功能分析。
[0009] 进一步的，在步骤1中，脑电电极采用国际10-20系统标准，以双耳的乳突作为参 考，从32导头皮脑电采集设备记录对应运动的脑电信号;使用Synamp2设备采集肌电信号， 电极顺着肌纤维方向，放置在肌腹位置处。
[0010] 进一步的，
[0011] 所述步骤3的具体方法如下：
[0012] 基于步骤2中预处理后的实测数据构建脑电信号x(t)和肌电信号y(t)两组时间序 列;对脑电信号x(t)进行小波变换，将小波变换应用于尺度划分之中；
[0013]首先构造小波函数，公式如下：
[0014] Wj，k(t)=2-」/2Ψ(2-jt_k) (I)
[0015] 式中，Ψ(?)为母小波;k为Ψ(?)纵坐标方向的平移量，j表示信号的层数，j，kez， Z为整数集;尺度参数为平移参数2?; t为时间指标；
[0016] 然后对脑电信号x(t)进行7层频谱分解，得到小波变换系数
[0017]
(2)
[0018] 小波系数C^k按照频率范围从高到低的顺序排列，提取第3、4、5、6、7层系数重构出 gamma(32~64Hz)、beta(16~32Hz)、alpha(8~16Hz)、theta(4~8Hz^Pdelta(l~4HzMj^ 段的信号：
[0019]
(3)
[0020] 式中，C =--"为滤波器的品质因数； '布'宽
[0021 ]重构信号Xi(t)、X2(t)、X3(t)、X4(t)及X5(t)分别对应脑电功能频带delta、theta、 alpha、beta 和 gamma 的信号；
[0022]对于肌电信号y(t)进行上述相同的小波变换过程，得到ym(t)(m=l，2,3,4,5)，* 别对应肌电delta、theta、alpha、beta和gamma频段的信号；
[0023]基于传递熵计算方法，构造 x(t)到y(t)的小波-传递熵WTE^，公式如下：
[0024]
(4)
[0025] 式中，u为预测时间；p( ·)为变量之间的联合概率X、分别表示脑电和肌电 delta、theta、alpha、beta、gamma分量的延迟向量;Mu为< 的预测序列；
[0026] WTEX-表示脑电信号EEG的Xi (t)分量到肌电信号t)分量间的传递熵 值；同理信号y (t)到X (t)的小波-传递熵WTEex的表达式为：
[0027]
(5)
[0028] 式中，<"为< 的预测序列;WTEpx表示Ei?^ym(t)到分量到EEG的Xi(t)分量间的传 递熵值;传递熵值越大，说明在此频段间皮层肌肉耦合越强;反之亦然。
[0029] 本发明与现有技术相比有益效果如下：本发明方法利用小波-传递熵分析脑肌电 信号信息传递特性，定量描述脑电信号与肌电信号频段间非线性的耦合及信息传递特征， 有助于探索大脑皮层与肌肉之间的功能联系，研究运动控制反馈机制及运动障碍病理机 制，建立基于脑肌电信号的康复状态评价指标，构建康复机器人运动状态及患者生理状态 评价机制，可以获得可观的社会效益和经济效益。
【附图说明】
[0030] 图1为本发明中Neuroscan设备的结构简图。
[0031] 图2为本发明方法的工作流程图。
[0032] 图3为受试者C4通道脑电信号的小波分解后时频结果图。
[0033] 图4为受试者指浅屈肌处肌电信号的小波分解后时频结果图。
[0034]图5为受试者的脑肌电信号小波-传递熵分析结果图。
[0035]附图标号：1_脑电电极、2-电极帽、3-脑肌电采集仪、4-肌电导联线、5-肌电电极。
【具体实施方式】
[0036]脑电信号和肌电信号十分微弱，具有非线性、非平稳性和频域特性突出等特点。在 运动过程中，神经系统与肌肉之间的交互控制机制可以通过脑肌电信号的同步耦合分析体 现。小波分解能够提取出脑肌电信号的特定时频数据段，而且传递熵可以刻画信号间非线 性的耦合及信息传递特征，本发明通过研究脑肌电间的小波-传递熵分析，获得不同运动状 态下大脑皮层与肌肉之间信息传递关系，进而研究运动功能障碍产生的生理机制。
[0037] 实施例1:
[0038]如图2所示，方法步骤如下：
[0039] 步骤1，采用64导Neuroscan设备同步采集脑电信号和肌电信号。
[0040] Neuroscan设备的结构如图1所示，由脑电电极、电极帽、脑肌电采集仪、肌电导联 线、肌电电极连接组成。
[0041 ]脑电信号采集:脑电电极采用国际标准10-20电极放置标准，通过电极帽2实现脑 电电极1与头皮接触。在手部静态握力输出运动下进行脑肌电信号同步采集实验。M1、M2导 联分别连接到左右耳后乳突作为参考电极，接地电极布置在头顶正中，从32导头皮脑电采 集设备中选择C3、C4以及CPZ区记录对应运动的脑电信号。
[0042] 肌电信号采集：使用Synamp 2设备采集指浅屈肌（f I exor d i g i t orum superficialis，FDS)处的肌电信号，首先用酒精擦拭被测部位的皮肤表面，去除皮肤表面 油脂和皮肩，然后顺着肌纤维方向将肌电电极5粘贴在肌腹位置处，并将肌电导联线4适当 固定尽量减少运动过程中导联线晃动的干扰。
[0043]步骤2,利用Neuroscan设备数据处理软件对采集到的脑电信号和肌电信号分别去 除基线漂移、溢出、眼动和工频干扰；
[0044] 步骤3,采用Daubechies类db4小波基函数对脑电和肌电信号进行频谱分解，基于 步骤2中预处理后的实测数据构建脑电信号x(t)和肌电信号y(t)两组时间序列。以脑电信 号x(t)为例，将小波变换应用于尺度划分之中。首先构造小波函数，公式如下：
[0045]
(1)
[0046] 式中，ψ(υ为母小波;k为ψ(υ纵坐标方向的平移量，j表示信号的层数，j，kez， Z为整数集;尺度参数为平移参数2?; t为时间指标。
[0047] 然后对脑电信号x(t)进行7层频谱分解，得到小波变换系数
[0048]
(2)
[0049]基十上式得到的小波糸数Cj,k按照频率范围从高到低的顺序排列，提取第3、4、5、 6、7层系数重构出8&臟&(32~64]^)、匕6七&(16~32]^)、&]^]1&(8~16]^)、1：]16七&(4~8]^)和 delta(l~4Hz)频段的信号：
[0050]
(3) Γ ? ^ i广滤波器中心频率a u m
[0051] 式中，C =-4一 -为滤波器的品质因数； 市苋
[0052] 重构信号Xi(t)、X2(t)、X3(t)、X4(t)及X5(t)分别对应脑电功能频带delta、theta、 alpha、beta 和gamma 的信号。
[0053] 对于肌电信号y(t)进行上述相同的小波变换过程，得到ym(t)(m=l，2,3,4,5)，* 别对应肌电delta、theta、alpha、beta和gamma频段的信号。
[0054] 基于传递熵计算方法，构造 x(t)到y(t)的小波-传递熵WTEq，公式如下：
[0055]
(4)
[0056] 式中，u为预测时间；p( ·)为变量之间的联合概率;<、<分别表示脑电和肌电 del ta、theta、alpha、beta、gamma分量的延迟向量;.V二为#的预测序列；
[0057] WTEX-表示脑电信号EEG的Xi (t)分量到肌电信号t)分量间的传递熵 值。同理信号y (t)到X (t)的小波-传递熵WTEpx的表达式为：
[0058]
(5)
[0059] 式中，为'的预测序列;WTEpx表示Ei?^ym(t)到分量到EEG的Xi(t)分量间的传 递熵值;传递熵值越大，说明在此频段间皮层肌肉耦合越强;反之亦然。
[0060] 基于上述指标，计算手部静态握力输出运动下，不同耦合方向上、不同时频尺度间 的WTE值，即能够定量描述脑电信号EEG和肌电信号EMG多时频尺度间的非线性同步耦合特 征。
[0061] 为验证本发明所述的脑肌电信号小波-传递熵分析方法的可行性和有效性，募集8 名健康的受试者进行手部静态握力输出实验，受试者相关信息如表1所示。按照本发明所述 的脑肌电采集与分析过程，同步采集受试者恒力输出运动下的脑肌电信号，进行分析并研 究受试者运动运动过程中皮层肌肉间的耦合与信息传递机制。
[0062] 本次实验采集左手指浅屈肌（flexor digitorum superficialis，FDS)处的肌电 信号与对侧C4通道脑电信号，并计算WTE值。
[0063] 图3和图4分别为受试者脑电、肌电信号经过小波分解后各时频尺度上的时频域结 果(左侧是时域图，右侧是频域图），从中可以看出脑电和肌电信号经过小波分解后可以得 至丨Jde Ita、theta、alpha、beta 和gamma 频段的信号。
[0064] 图5为8位受试者脑肌间耦合的小波-传递熵分析后的平均值。从中可以看出，在静 态握力输出过程中皮层肌肉间beta频段的耦合强度最为显著，而且不同耦合方向上各时频 尺度间的耦合强度也存在差异，为探究神经肌肉功能耦合机制提供了理论研究方法。
[0066]表1受试者相关信息。
【主权项】
1. 一种基于小波-传递赌的多时频尺度间脑肌电禪合分析方法，其特征在于，所述方法 的步骤如下： 步骤1，采用64导化uroscan设备同步采集脑电信号和肌电信号； 步骤2,利用化uroscan设备数据处理软件对采集到的脑电信号和肌电信号分别去除基 线漂移、溢出、眼动和工频干扰； 步骤3,采用DaubecMes类化4小波基函数对脑电和肌电信号进行频谱分解，分析脑电 信号和肌电信号不同时频尺度间的同步特性，定量描述脑肌间非线性禪合和信息传递特 征； 步骤4,对运动状态下的脑肌间非线性禪合和信息传递特征进行分析。2. 根据权利要求1所述的基于小波-传递赌的多时频尺度间脑肌电禪合分析方法，其特 征在于:步骤1中，脑电电极采用国际10-20系统标准，W双耳的乳突作为参考，从32导头皮 脑电采集设备记录对应运动的脑电信号;使用Synamp2设备采集肌电信号，电极顺着肌纤维 方向，放置在肌腹位置处。3. 根据权利要求1所述的基于小波-传递赌的多时频尺度间脑肌电禪合分析方法，其特 征在于，所述步骤3的具体方法如下： 基于步骤2中预处理后的实测数据构建脑电信号x(t)和肌电信号y(t)两组时间序列； 对脑电信号x(t)进行小波变换，将小波变换应用于尺度划分之中； 首先构造小波函数，公式如下： Wj'k(t) = 2-j/2w(2_jt-k) (1) 式中，ψ(υ为母小波;k为Ψ(υ纵坐标方向的平移量，j表示信号的层数，j，kez，z为 整数集;尺度参数为户，平移参数2?; t为时间指标； 然后对脑电信号x(t)进行7层频谱分解，得到小波变换系数口） 小波系数Cj,k按照频率范围从高到低的顺序排列，提取第3、4、5、6、7层系数重构出 gamma(32~64Hz)、beta(16~32Hz)、al地a(8~16Hz)、theta(4~細z)和delta(l~4Hz)频 段的信号：餅 式中为滤波器的品质因数； 重构信号Xi(t)、X2(t)、X3(t)、X4(t)及X日（t)分别对应脑电功能频带delta、theta、 alpha、beta 和 gamma 的信号； 对于肌电信号y(t)进行上述相同的小波变换过程，得到ym(t)(m=l，2,3,4,5)，分别对 应肌电delta、theta、al地a、be1:a和gamma频段的信号； 基于传递赌计算方法，构造 x(t)到y(t)的小波-传递赌WTEx·^，公式如下：(4) 式中，U为预测时间;p( ·)为变量之间的联合概率;1;、沁分别表示脑电和肌电delta、 theta、al地a、be化、gamma分量的延迟向量；.>寬,为.)'Γ的预测序列； WTE^^y则表示脑电信号EEG的xi(t)分量到肌电信号EMG的ym(t)分量间的传递赌值;同理 信号y (t)到X (t)的小波-传递赌WTE^-x的表达式为：(5) 式中，故。为别勺预测序列;WTE^-x表示EMG的ym(t)到分量至化EG的xi(t)分量间的传递赌 值;传递赌值越大，说明在此频段间皮层肌肉禪合越强;反之亦然。
【文档编号】A61B5/0476GK106073702SQ201610362111
【公开日】2016年11月9日
【申请日】2016年5月27日
【发明人】谢平, 杨芳梅, 张园园, 陈晓玲, 吴晓光, 张晋铭, 王霄
【申请人】燕山大学