一种基于电生理与血液微循环信号的脑功能耦合分析方法

技术领域

本发明涉及信号处理技术领域，更具体地，涉及一种基于电生理与血液微循环信号的脑功能耦合分析方法。

背景技术

脑卒中是导致运动功能障碍的病因之一。脑卒中会导致高级中枢神经元受损，下运动神经元失去控制，反射活动活跃等。而神经系统参与调节人体几乎所有的功能活动，中枢神经系统包括脑和脊髓，外周神经系统包括肢体运动神经和支配心血管及心肺等内脏器官活动的自主神经。大脑通过脊髓传导通路接收感觉信息并发出运动指令控制肢体运动，通过自主神经调节心脏搏动、呼吸运动等；外周神经通过发放动作电位支配骨骼肌、心肌、呼吸肌等肌肉的运动。另一方面，血液循环的正常运行为全身每个组织器官的活动提供氧气和营养物质并排出二氧化碳和代谢物。循环系统功能损伤往往会导致神经系统功能损伤，并且神经的代谢可通过神经的血流和血氧状况表现出来。

人体各系统、组织和器官之间存在普遍的耦合(Coupling)关系，各系统、组织和器官的协同工作和相互影响决定了人体行为与功能的实现。由于运动过程中人体的神经肌肉耦合规律尚不明晰，传统的运动训练康复疗法在运动功能康复治疗中效果欠佳。为了科学制定运动功能障碍患者的运动功能康复方案，理解运动过程中的耦合规律，从信息学的角度来进行研究是一种十分科学的方法。神经血管单元(Neurovascular Unit，NVU)是观测神经-血液循环关系的重要研究单元，而电生理信号与血液微循环信号又是检测神经与血管活动的两种重要信号。同步采集电生理与血液微循环信号，对其进行预处理和尺度化分析后，再采用能够量化非线性耦合特征的参数进行分析，探究耦合机制，这对发展运动功能综合干预方法，重建患者神经肌肉功能，开发新型康复机器人系统实现应用，提升偏瘫卧床患者的运动功能康复效果具有显著意义。

在现有技术中，分析多模信号耦合关系的方法有很多种，按分析的耦合关系类型区分，可分为线性和非线性分析方法；按分析的功能区分，可分为无向和定向分析方法。

例如，相干性是一种线性的、无向的分析方法，用于度量两个时间序列之间的线性相关程度，定义为互谱密度对自谱密度的归一化。由于其算法原理简单和稳定性强，相干性被广泛用来描述生理信号的耦合关系。但是，相干性需要对两个信号进行频域分析，因此需要使用频域算法，而频域算法的计算较为复杂，需要占用较大的计算资源和时间。相干性无法区分信号的振幅和相位，且应用前提是假设分析的信号是平稳的。传统相干性分析不足以准确描述信号之间的非线性特征，而且可能会造成耦合强度的高估。

互信息是信息论中的一种信息度量工具，用于表示一个随机变量中蕴含关于另外一个随机变量的信息量。它是一种非线性、无向的分析方法，用于计算两个时间序列之间的统计独立性。但是，互信息计算需要较长的数据，并且估计互信息时需要依赖联合概率密度估计。

格兰杰因果是一种线性的、定向的分析方法，能够描述双向信息传递，计算相对简单，但是它只能测得时间序列之间的线性同步因果耦合关系以及信息流向，并且需要时间序列平稳，对于非线性、高阶等复杂的因果耦合关系难以定量分析。并且格兰杰因果方法无法处理噪声干扰等问题，其模型需要大量的数据，并且需要对数据进行适当的预处理，才能获得较好的估计效果。

经验模态分解是一种自适应的信号分解方法，不论是非线性信号还是非平稳信号都能够进行有效的分解，可以更好的描述信号的瞬态行为和频谱特征。但经验模态分解的结果可能存在模态混叠现象，不同模态之间会相互干扰或重叠，影响信号的分析。

集合经验模态分解能够克服经验模态分解带来的随机性问题，在一定程度上消除了模态混叠问题，由于它的自适应性和灵活性，能够适应不同类型的信号和噪声，提高了信号分解的适用范围和准确性。但该种分解方法需要进行多次噪声扰动和分解，噪声水平过高或分解次数过少都可能会影响分解结果的质量。

经分析，尽管已有多种方法能够实现信号的时频分析，然而对于非周期性信号，传统的傅里叶变换不能给出有效的表示和分析，并且存在精度问题，在信号的高频率部分，由于正弦和余弦函数的振幅和相位不确定，可能会导致信号处理的误差和失真，并且针对大规模的信号处理时，计算复杂度较高，成本较大。经验模态分解方法虽然能自适应地分解信号，能够更好的描述信号的频谱特征，但是它分解的结果可能存在模态混叠问题，影响信号的分析。

综上，现有的分析生理信号的耦合关系方法有相干性、互信息、格兰杰因果等。传统的相干性方法是一种线性的、无向的分析方法，相干性无法区分信号的振幅和相位，且应用前提是假设分析的信号是平稳的。传统相干性分析不足以准确描述信号之间的非线性特征，而且可能会造成耦合强度的高估。互信息是一种非线性、无向的分析方法，用于计算两个时间序列之间的统计独立性。但是，互信息只能衡量两个变量之间的相关性，并不能准确地反映两个变量之间的因果关系，而且互信息计算需要较长的数据，估计互信息时需要依赖联合概率密度估计。格兰杰因果是一种线性的、定向的分析方法，但是它只能测得时间序列之间的线性同步因果耦合关系以及信息流向，并且需要时间序列平稳，对于非线性、高阶等复杂的因果耦合关系难以定量分析。并且格兰杰因果方法无法处理噪声干扰等问题，其模型需要大量的数据，并且需要对数据进行适当的预处理，才能获得较好的估计效果。

发明内容

本发明的目的是克服上述现有技术的缺陷，提供一种基于电生理与血液微循环信号的脑功能耦合分析方法。该方法包括以下步骤：

在不同运动状态下，同步采集受试者的脑电信号和fNIRS信号并进行预处理，得到预处理的脑电信号和脑氧信号；

对于所述脑电信号进行完全自适应噪声集合经验模态分解，获得多个不同子频带分量；

对所述脑电信号的多个不同子频带分量和所述脑氧信号的时间序列进行符号化，获得不同子频带分量的符号化序列以及脑氧信号的符号化序列；

计算所述不同子频带分量的符号化序列与所述脑氧信号的符号化序列两两之间的传递熵；

基于所述传递熵分析所述脑电信号与所述脑氧信号的耦合程度。

与现有技术相比，本发明的优点在于，针对目前神经系统中电生理信号与血液微循环信息耦合机制模糊不清，难以对脑卒中偏瘫卧床患者提出多角度协同康复干预手段的问题，本发明提出一种基于电生理与血液微循环信号的耦合关系分析方法，从不同尺度上来探究电生理与血液微循环信号的非线性耦合特征及信息传递，为进一步探索运动功能障碍的产生机理和提出协同康复干预手段提供依据。

通过以下参照附图对本发明的示例性实施例的详细描述，本发明的其它特征及其优点将会变得清楚。

附图说明

被结合在说明书中并构成说明书的一部分的附图示出了本发明的实施例，并且连同其说明一起用于解释本发明的原理。

图1是根据本发明一个实施例的基于电生理与血液微循环信号的脑功能耦合分析方法的流程图；

图2是根据本发明一个实施例的基于电生理与血液微循环信号的脑功能耦合分析方法的实现过程示意图。

图3是根据本发明一个实施例的同步采集信号通道示意图；

图4是根据本发明一个实施例的原始脑电信号示意图；

图5是根据本发明一个实施例的预处理后的脑电信号示意图；

图6是根据本发明一个实施例的原始fNIRS信号示意图；

图7是根据本发明一个实施例的预处理后的fNIRS信号示意图；

图8是根据本发明一个实施例的经过完全自适应噪声集合经验模态分解后的脑电模态以及对应频谱图；

图9是根据本发明一个实施例的脑氧信号至各模态脑电信号的符号化传递熵变化曲线示意图；

图10是根据本发明一个实施例的各模态脑电信号至脑氧信号的符号化传递熵变化曲线示意图。

具体实施方式

现在将参照附图来详细描述本发明的各种示例性实施例。应注意到：除非另外具体说明，否则在这些实施例中阐述的部件和步骤的相对布置、数字表达式和数值不限制本发明的范围。

以下对至少一个示例性实施例的描述实际上仅仅是说明性的，决不作为对本发明及其应用或使用的任何限制。

对于相关领域普通技术人员已知的技术、方法和设备可能不作详细讨论，但在适当情况下，所述技术、方法和设备应当被视为说明书的一部分。

在这里示出和讨论的所有例子中，任何具体值应被解释为仅仅是示例性的，而不是作为限制。因此，示例性实施例的其它例子可以具有不同的值。

应注意到：相似的标号和字母在下面的附图中表示类似项，因此，一旦某一项在一个附图中被定义，则在随后的附图中不需要对其进行进一步讨论。

本发明公开了一种基于电生理与血液微循环信号的脑功能耦合分析方法，采用完全自适应噪声集合经验模态分解-符号化传递熵定量分析信息传递以及耦合程度。在下文中，电生理信号以脑电信号为例，血液微循环信号以脑氧信号为例进行描述。脑电信号是脑神经组织的电生理活动在大脑皮层表面的总体反映。脑氧信号可以反映血液中氧气的含量。

结合图1和图2所示，所提供的基于电生理与血液微循环信号的脑功能耦合分析方法包括以下步骤：

步骤S110，在不同运动状态下，同步采集受试者的脑电信号和fNIRS信号。

首先，操作员连接脑电设备、fNIRS(功能性近红外光谱)设备和同步盒等硬件设备。受试者按照操作员指令完成不同实验任务，在实验过程中同步采集脑电信号和fNIRS信号，从而获得在不同运动状态下受试者的脑电信号和fNIRS信号。图3是同步采集信号通道示意图。

例如，招募了3名健康受试者，均无运动功能障碍或神经系统疾病。实验采用的脑电设备为NeuroScan，采样频率为1000Hz。脑氧采集设备为NirScan，采样频率为17Hz。受试者按照操作员指令分别完成上肢手摇车运动抗阻任务、下肢脚踏车抗阻任务和上下肢联动抗阻任务，每个任务分别设置0、5、10三个档位的阻力，每个任务持续时间为60s。为避免受试者疲劳，每组任务结束后休息2min后再开展下一组实验。在实验过程中同步采集受试者脑电信号和fNIRS信号，并以数字信号存储在电脑中。

步骤S120，对采集到的脑电信号和fNIRS信号进行预处理，得到预处理的脑电信号和脑氧信号。

具体地，对采集到的脑电信号进行预处理，包括去除基线漂移、降采样、利用FIR数字滤波器进行0.1-50Hz的带通滤波、去除50Hz的工频信号、使用独立成分分析(ICA)算法去除脑电信号中的眼电和肌电干扰成分，最终获得有效的脑电信号，并选择C3，C4通道的脑电数据用作后续的耦合计算。图4是原始脑电信号示意图。图5是预处理后的脑电信号示意图。

对采集到的fNIRS信号进行预处理，包括将光密度信号通过比尔朗伯定律转化为三种血流动力参数的时间序列：含氧血红蛋白浓度，脱氧血红蛋白浓度和总血红蛋白浓度，上述三种均为相对浓度变化。使用0.01-0.2Hz的带通滤波器滤波，将滤波后的信号进行分段。使用含氧血红蛋白相对浓度表征大脑皮层激活程度，用以后续和脑电子频带分量耦合计算的脑氧信号，并将其按照左侧运动皮层和右侧运动皮层进行划分。fNIRS信号转化为大脑中含氧血红蛋白浓度变化，以此代表脑氧信号。图6原始fNIRS信号示意图。图7是预处理后的fNIRS信号示意图。

步骤S130，对预处理的脑电信号进行完全自适应噪声集合经验模态分解，将其分解为多个不同子频带分量。

例如，使用完全自适应噪声集合经验模态分解算法(Complete EnsembleEmpirical Mode Decomposition with Adaptive Noise，CEEMDAN)将脑电信号分解为不同频段的子分量，求第一阶IMF(本征模函数)，将正负成对的高斯白噪声加入原始信号y(t)中，得到新信号为：

y(t)+(-1)

其中，t表示时刻，q＝1和2，ε为幅值，n

对新信号进行EMD(经验模态分解)分解，得到第一阶IMF分量u

其中，E表示经验模态分解，

对产生的N个模态分量进行总体平均计算，表示为：

其中，

则余量信号为：

对该新信号进行EMD分解，得到对应的第一阶IMF分量u

则余量信号表示为：

其中，r

重复上述操作，直到获得的余量信号为单调信号或者极值点不超过两个的时候，余量信号不能再继续分解，算法结束，此时得到的本征模态分量为K，则原始信号可表达为：

完全自适应噪声集合经验模态分解可以有效去除信号中的噪声，提高了信号的质量和准确性。CEEMDAN分解后的输出是不同频率和幅度的子频带分量，代表了原始信号在不同时间尺度和频率上的本质特征。CEEMDAN分解后的输出是不同频率和幅度的子频带分量，代表了原始信号在不同时间尺度和频率上的本质特征。图8是经过完全自适应噪声集合经验模态分解后的脑电模态以及对应频谱图。

步骤S140，将脑电信号的各个子频带分量以及脑氧信号时间序列进行符号化，将其转变为符号化观测值。

具体地，对于分解得到的脑电信号子带分量，将脑电信号的各个子频带分量与脑氧信号时间序列进行符号化，表示为：

其中，S(i)为符号化后的序列，x(i)为原始时间序列，delta是递增量，取值为(max(x)-min(x))/nums，min(x)表示x(i)的最小值，max(x)表示x(i)的最大值，nums表示为该时间序列的符号化观测值。

在该步骤中，对脑电信号的各个子频带分量进行符号化，得到各个子频带分量的符号化序列，并且将脑氧信号进行符号化，得到脑氧信号的符号化序列。时间序列在符号化时能够捕捉大尺度的特征，并且符号化序列保存了更重要的动力学特征。

步骤S150，针对脑电信号子带分量以及脑氧信号的符号化序列，利用传递熵计算各脑电信号子带分量与左右脑区脑氧信号的耦合关系。

具体地，基于上述分解得到的脑电信号子带分量以及脑氧信号的符号化序列，利用传递熵计算各脑电信号子带分量与左右脑区脑氧信号的耦合关系，表示如下：

X和Y为两个不同的时间序列，信号Y到信号X传输信息的传递熵定义为：

其中，x

在该步骤中，选取不同的子频带分量与脑氧信号的符号化序列两两之间进行传递熵计算，观察其耦合方向及强度。传递熵是一种非线性、定向的分析方法，本身包含着因果方向和动态强度，能够很好地量化生理信号的非线性特征，可用作研究两个时间序列间的非线性信息传递。传递熵因具有无模型、非线性、定量分析等优势，因而能够用于探索神经生理信号间的非线性耦合特征。将时间序列符号化之后再进行传递熵的计算，可以避免传递熵空间提取的计算复杂度。

综上，本发明实施例采用完全自适应噪声集合经验模态分解结合符号化传递熵的优势，有利于处理电生理与血液微循环信号，定量描述在不同运动状态下脑电信号和脑氧信号的特征频段的耦合特性，为进一步探究电生理与血液微循环之间的耦合机制提供新的研究方向。

步骤S160，绘制符号化传递熵变化曲线。

基于步骤S150的计算结果，可进一步绘制符号化传递熵变化曲线，以观察不同运动状态、不同尺度下符号化传递熵的变化，并统计分析观察符号化传递熵值的显著变化。图9脑氧信号至各模态脑电信号的符号化传递熵变化曲线示意图。图10是各模态脑电信号至脑氧信号的符号化传递熵变化曲线示意图。

综上，本发明提出的完全自适应噪声集合经验模态分解-符号化传递熵方法，对于分析电生理信号与血液微循环信号的耦合关系有很大优势。完全自适应噪声集合经验模态分解能够解决模态混叠问题，引入不同尺度的噪声扰动，在分解过程中将干扰分配到不同的本征模态函数中。这种方法支持对非线性、非平稳和非高斯信号的分解和分析，具有较强的适应性和普适性，能够提高分解结果的准确性和重复性。

综上所述，本发明基于电生理与血液微循环信号，从不同尺度上来探究电生理与血液微循环信号的非线性耦合特征及信息传递。在不同运动态下同步采集人体的电生理与血液微循环信号，使用完全自适应噪声集合经验模态分解方法对划分脑电信号的频域尺度，并采用能够量化非线性特征的符号化传递熵参数分析信号间的耦合关系，探究耦合机制。这种方法支持对非线性、非平稳和非高斯信号的分解和分析，具有较强的适应性和普适性，能够提高分解结果的准确性和重复性。相对于现有技术，本发明具有以下优势：

1)本发明采用完全自适应噪声集合经验模态分解方法分解预处理后的脑电信号，加入经EMD分解后含辅助噪声的本征模态分量，而不是将高斯白噪声信号直接添加在原始信号中，提升了信号处理的普适性。

2)本发明利用完全自适应噪声集合经验模态分解得到脑电信号的子带信号，将脑电信号进行频带划分，并将时间序列进行符号化。并且在得到第一阶本征模态分量时就对产生的模态分量进行总体平均，得到最终的第一阶本征模态分量，然后对余量信号重复分解步骤。这种方式能够捕捉不同尺度的特征，并且符号化序列保存了重要的动力学特征。

3)本发明采用符号化传递熵量化脑电信号与脑氧信号的信息传递，较传统的相干性方法而言，符号化传递熵是一种非线性、定向的分析方法，包含着因果方向和动态强度，可用作研究两个时间序列间的非线性信息传递，探索信号的非线性耦合特征，对于探索神经生理信号的耦合机制有很大助益。

4)经过论证和预实验，本发明能够为进一步探索运动功能障碍的产生机理和提出协同康复干预手段提供依据。

本发明可以是系统、方法和/或计算机程序产品。计算机程序产品可以包括计算机可读存储介质，其上载有用于使处理器实现本发明的各个方面的计算机可读程序指令。

计算机可读存储介质可以是可以保持和存储由指令执行设备使用的指令的有形设备。计算机可读存储介质例如可以是但不限于电存储设备、磁存储设备、光存储设备、电磁存储设备、半导体存储设备或者上述的任意合适的组合。计算机可读存储介质的更具体的例子(非穷举的列表)包括：便携式计算机盘、硬盘、随机存取存储器(RAM)、只读存储器(ROM)、可擦式可编程只读存储器(EPROM或闪存)、静态随机存取存储器(SRAM)、便携式压缩盘只读存储器(CD-ROM)、数字多功能盘(DVD)、记忆棒、软盘、机械编码设备、例如其上存储有指令的打孔卡或凹槽内凸起结构、以及上述的任意合适的组合。这里所使用的计算机可读存储介质不被解释为瞬时信号本身，诸如无线电波或者其他自由传播的电磁波、通过波导或其他传输媒介传播的电磁波(例如，通过光纤电缆的光脉冲)、或者通过电线传输的电信号。

这里所描述的计算机可读程序指令可以从计算机可读存储介质下载到各个计算/处理设备，或者通过网络、例如因特网、局域网、广域网和/或无线网下载到外部计算机或外部存储设备。网络可以包括铜传输电缆、光纤传输、无线传输、路由器、防火墙、交换机、网关计算机和/或边缘服务器。每个计算/处理设备中的网络适配卡或者网络接口从网络接收计算机可读程序指令，并转发该计算机可读程序指令，以供存储在各个计算/处理设备中的计算机可读存储介质中。

用于执行本发明操作的计算机程序指令可以是汇编指令、指令集架构(ISA)指令、机器指令、机器相关指令、微代码、固件指令、状态设置数据、或者以一种或多种编程语言的任意组合编写的源代码或目标代码，所述编程语言包括面向对象的编程语言—诸如Smalltalk、C++、Python等，以及常规的过程式编程语言—诸如“C”语言或类似的编程语言。计算机可读程序指令可以完全地在用户计算机上执行、部分地在用户计算机上执行、作为一个独立的软件包执行、部分在用户计算机上部分在远程计算机上执行、或者完全在远程计算机或服务器上执行。在涉及远程计算机的情形中，远程计算机可以通过任意种类的网络—包括局域网(LAN)或广域网(WAN)—连接到用户计算机，或者，可以连接到外部计算机(例如利用因特网服务提供商来通过因特网连接)。在一些实施例中，通过利用计算机可读程序指令的状态信息来个性化定制电子电路，例如可编程逻辑电路、现场可编程门阵列(FPGA)或可编程逻辑阵列(PLA)，该电子电路可以执行计算机可读程序指令，从而实现本发明的各个方面。

这里参照根据本发明实施例的方法、装置(系统)和计算机程序产品的流程图和/或框图描述了本发明的各个方面。应当理解，流程图和/或框图的每个方框以及流程图和/或框图中各方框的组合，都可以由计算机可读程序指令实现。

这些计算机可读程序指令可以提供给通用计算机、专用计算机或其它可编程数据处理装置的处理器，从而生产出一种机器，使得这些指令在通过计算机或其它可编程数据处理装置的处理器执行时，产生了实现流程图和/或框图中的一个或多个方框中规定的功能/动作的装置。也可以把这些计算机可读程序指令存储在计算机可读存储介质中，这些指令使得计算机、可编程数据处理装置和/或其他设备以特定方式工作，从而，存储有指令的计算机可读介质则包括一个制造品，其包括实现流程图和/或框图中的一个或多个方框中规定的功能/动作的各个方面的指令。

也可以把计算机可读程序指令加载到计算机、其它可编程数据处理装置、或其它设备上，使得在计算机、其它可编程数据处理装置或其它设备上执行一系列操作步骤，以产生计算机实现的过程，从而使得在计算机、其它可编程数据处理装置、或其它设备上执行的指令实现流程图和/或框图中的一个或多个方框中规定的功能/动作。

附图中的流程图和框图显示了根据本发明的多个实施例的系统、方法和计算机程序产品的可能实现的体系架构、功能和操作。在这点上，流程图或框图中的每个方框可以代表一个模块、程序段或指令的一部分，所述模块、程序段或指令的一部分包含一个或多个用于实现规定的逻辑功能的可执行指令。在有些作为替换的实现中，方框中所标注的功能也可以以不同于附图中所标注的顺序发生。例如，两个连续的方框实际上可以基本并行地执行，它们有时也可以按相反的顺序执行，这依所涉及的功能而定。也要注意的是，框图和/或流程图中的每个方框、以及框图和/或流程图中的方框的组合，可以用执行规定的功能或动作的专用的基于硬件的系统来实现，或者可以用专用硬件与计算机指令的组合来实现。对于本领域技术人员来说公知的是，通过硬件方式实现、通过软件方式实现以及通过软件和硬件结合的方式实现都是等价的。

以上已经描述了本发明的各实施例，上述说明是示例性的，并非穷尽性的，并且也不限于所披露的各实施例。在不偏离所说明的各实施例的范围和精神的情况下，对于本技术领域的普通技术人员来说许多修改和变更都是显而易见的。本文中所用术语的选择，旨在最好地解释各实施例的原理、实际应用或对市场中的技术改进，或者使本技术领域的其它普通技术人员能理解本文披露的各实施例。本发明的范围由所附权利要求来限定。