一种神经肌肉信息交互模型构建及参数辨识优化方法

本发明涉及神经信息处理及神经动力学领域，特别涉及一种基于神经肌肉信息交互模型及参数辨识优化方法。

背景技术：

大脑运动皮层的脑电信号(electroencephalogram,eeg)和相应肢体肌肉组织的肌电信号(electromyography,emg)分别反映了运动控制和肌肉功能响应信息，所以脑电和肌电信号的同步特征分析能够体现出大脑皮层与肌肉组织之间的功能联系特征，进而体现神经肌肉控制系统的功能状态。eeg和emg信号间的神经肌肉功能耦合(functionalcorticomuscularcoupling,fcmc)现象能够反映运动控制系统中大脑与肌肉间不同层面的信息交互及传递，有助于从系统层面揭示运动控制系统神经网络的协同工作方式，进一步理解大脑皮层对于信息处理和选择运动产生模式的能力，为运动功能评价提供全新视角。研究受神经肌肉功能耦合机制约束下的模型构建及参数辨识方法，有助于从微观层面揭示上位运动神经元与肢体动作电位间的信息交互关系。

目前，大多是从功能耦合的角度探究大脑皮层和躯体肌肉之间的信息交互，但基于脑电和肌电信号的神经肌肉功能耦合分析只是通过数学方法客观计算推断与之相应的运动控制系统的内在生理机制，未能直接表征运动控制系统内部结构和功能的真实变化。而系统建模为解决这一问题提供了很好的途径，通过对控制系统内在结构建模，以可变参量的形式模拟仿真控制系统内在的变化因子，实现对系统内部特性的深度解析。

目前，基于仿真模型探究人体生理内在变化机制的研究中，神经元群模型或运动神经元模型已被广泛应用，但都局限于探究大脑神经元的实际作用机制关系。在此基础上，本发明基于神经元群模型及下肢动作神经元模型，建立大脑与肌肉间的运动控制和感觉反馈的闭环控制环路模型，基于eeg和emg信号获取神经肌肉耦合强度、耦合方向及延迟时间等多层次特征指标，构建基于多层次耦合特征指标的目标函数，并基于无迹卡尔曼滤波(unscentedkalmanfilter,ukf)方法和混合粒子群算法进行模型参数有效辨识及优化选取，实现神经肌肉信息交互模型的构建，客观、定量地描述神经肌肉系统的运动功能状态。

技术实现要素：

针对现有技术存在的问题，本发明提供一种神经肌肉信息交互模型构建及参数辨识优化方法，神经肌肉信息交互模型构建及参数辨识优化方法具体实施步骤如下：

s1、分别构建运动控制通路模型部分和感觉反馈通路模型部分：

s11、构建运动控制通路模型部分：

s111、给定引起运动控制模块中电位u变化的外界信号，获得与运动控制模块中电位u相关的局部兴奋性平均突触增益参数he和兴奋性细胞膜平均时间常数和树突平均时间常数的和εe，通过相关参数得到运动控制模块的模型函数表达式如下：

式中，u(t)为模型输入信号，表示运动控制模块中动作电位的变化，x(t)代表运动控制模块中局部场电位，代表由运动控制模块传至感觉反馈模块的兴奋信号状态量；

s112、根据步骤s111中x0(t)的变化，通过控制感觉反馈模块中元群募集α与γ，引起感觉反馈模块中元群的变化，输出感觉反馈模块中单元动作电位y(t)，其函数关系如下：

y(t)＝[α(αs,αd)+γ(γs,γd)]g[x0(t+τ)](2)

式中，g[·]为x(t)和y(t)间的具有一定延时的信息传递函数，x0(t)代表由运动控制模块传至感觉反馈模块的兴奋信号，τ为延时因子，α(·)和γ(·)分别为α与γ感觉反馈模块元群的增益函数，αs,αd,γs,γd分别为静态和动态感觉反馈模块中的元群；

s12、构建感觉反馈通路模型部分：

s121、输入感觉反馈模块中的外界感觉信号y0(t)，通过中介体ia和ib传回运动控制模块，产生感觉反馈模块至运动控制模块的反馈传入信号z(t)，其具体的表达式为：

式中，l[·]为z(t)和y0(t)间的信息传递函数，l是牵张反射增益，s是交互兴奋增益，r是交互抑制增益，b是ib的增益，g是运动控制模块的中间抑制性元群的增益，v'和分别为来自于中介体ia和ib传入运动控制模块的放电频率，g'(t)是随机的已知向量；

s122、基于步骤s11中运动控制模块的函数关系，得到x2(t)与z(t)关系如下：

式中，hi为抑制性平均突触增益参数，εi为抑制性细胞膜平均时间常数和树突平均时间常数的和，z(t)为感觉反馈模块至运动控制模块的反馈传入信号，为运动控制模块中部分反馈信号状态量；

s13、综合步骤s11和步骤s12中运动控制模块与感觉反馈模块的作用关系，得到运动控制模块信号与感觉反馈模块信号的表达关系分别为：

式中，n1＝[αγ00]，n2＝[(αs,αd)(γs,γd)00]^t分别为未知参数向量，he为运动控制模块中局部兴奋性平均突触增益参数，εe为兴奋性细胞膜平均时间常数和树突平均时间常数的和，hi为运动控制模块中局部抑制性平均突触增益参数，εi为抑制性细胞膜平均时间常数和树突平均时间常数的和，s是交互兴奋增益，r是交互抑制增益，b是ib的增益，v'和分别为来自于中介体ia和ib传入运动控制模块的放电频率，α(·)和γ(·)分别为感觉反馈模块元群α与γ的增益函数；v(t)＝g[x0(t+τ)]分别为随机的已知向量；

s2、基于多层次耦合特征的参数辨识及优化：

s21、基于多层次耦合特征指标构建参数辨识的目标函数：

s211、对采集的数据进行预处理，并用非线性多尺度脑肌电同步耦合分析方法建立耦合强度vc、耦合方向tse和时间延迟td多层次特征指标；

s2111、基于时频一致性方法，获取耦合强度vc；

s2112、基于传递谱熵，获取耦合方向tse；

s2113、基于时延信息传递，获取时间延迟td；

s212、利用方程误差参数辨识方法构建基于多层次耦合特征指标的目标函数，目标函数的具体表达式为：

其中，ε(k)代表k时刻模型输出与系统输出的偏差，是模型参数估计值，是模型参数估计值下所有时刻的偏差和；

s22、基于无迹卡尔曼滤波(ukf)以及步骤s212构建的目标函数对神经肌肉信息交互模型进行参数辨识：

s221、令脑电和肌电信号表达式的输入信号u(t)服从高斯分布；

s222、基于步骤s211建立的多层次特征指标，利用增广状态向量构建状态方程，将模型中的未知参数向量定义为系统中新增状态向量，具体表示为：

x'(k)＝[x(k)h(k)ξ(k)](7)

其中，x为特征指标，h＝[m1m2]和ξ＝n1·n2分别表示式(5)中的运动控制模块参数和感觉反馈模块参数；

新的状态方程如下：

式中，x'(k+1)为n维随机状态向量序列，z(k+1)为n维系统可观测输出变量，包含x(t)和y(t)，f(·)和h(·)分别为n维和m维函数，u(k)为系统输入变量；

s223、利用步骤212中的目标函数并基于无迹卡尔曼滤波中采样策略逼近非线性分布的估计原则，对步骤s222中的未知参数h和ξ进行辨识；

s23、基于混合粒子群算法对参量h＝[hvchtsehtd]及ξ＝[ξvcξtseξtd]进行最优选取。

可优选的是，运动控制模块中的元群信号x(t)和感觉反馈模块中的单元动作电位y(t)的变化分别与外界信号的刺激有关。

可优选的是，在步骤s223中，根据无迹卡尔曼滤波中采样策略逼近非线性分布的估计原则，当状态值x(k+1)满足时，辨识过程结束。

可优选的是，多层次耦合特征指标下的目标函数参数集合，是根据目标函数的输出有效逼近基于实测信号计算的u(k)值。

可优选的是，步骤s23中的最优选取，是基于对实测控制系统神经肌肉功能耦合(fcmc)分析指标的逼近优化原则，以实测运动控制模块和感觉反馈模块中信号能量谱为目标。

本发明与现有技术相比，具有如下优点：

在神经肌肉系统信息交互仿真建模层面，将神经元群模型和运动神经元模型有机整合，构建神经肌肉信息交互模型，构建目标函数并基于无迹卡尔曼滤波方法对仿真模型参数进行有效辨识，再以脑电和肌电信号能量谱为目标，基于混合粒子群算法对所辨识参量进行优化选取，从微观上为运动功能评价提供研究新思路。

附图说明

图1为本发明一种神经肌肉信息交互模型构建及参数辨识优化方法的总体研究方案；

图2为本发明一种神经肌肉信息交互模型构建及参数辨识优化方法的参数辨识及优化策略。

具体实施方式

为详尽本发明之技术内容、结构特征、所达成目的及功效，以下将结合说明书附图进行详细说明。

一种神经肌肉信息交互模型构建及参数辨识优化方法，如图1所示。以运动控制过程中大脑信息与肌肉信息交互为基础，建立大脑与肌肉之间的信息交互仿真模型，构建目标函数并基于无迹卡尔曼滤波方法对仿真模型参数进行有效辨识，再以脑电和肌电信号能量谱为目标，基于混合粒子群算法对所辨识参量进行优化选取，实现神经肌肉信息交互模型的构建。

本发明是这样实现的：

s1、分别构建运动控制通路模型部分和感觉反馈通路模型部分。

s2、基于多层次耦合特征的参数辨识及优化。

下面将对本发明进行详细说明。

s1、分别构建运动控制通路模型部分和感觉反馈通路模型部分。

神经肌肉信息交互模型主要分为运动控制通路建模和感觉反馈通路建模两部分。

s11、构建运动控制通路模型部分：

s111、给定引起运动控制模块中电位u变化的外界信号，获得与运动控制模块中电位u相关的局部兴奋性平均突触增益参数he、兴奋性细胞膜平均时间常数和树突平均时间常数的和εe，通过相关参数得到运动控制模块的模型函数表达式如下：

式中，u(t)为模型输入信号，表示运动控制模块中动作电位的变化，x(t)代表运动控制模块中局部场电位，代表由运动控制模块传至感觉反馈模块的兴奋信号状态量；

s112、根据步骤s111中x0(t)的变化，通过控制感觉反馈模块中元群募集α与γ，引起感觉反馈模块中元群的变化，输出感觉反馈模块中单元动作电位y(t)，其函数关系如下：

y(t)＝[α(αs,αd)+γ(γs,γd)]g[x0(t+τ)](2)

式中，g[·]为x(t)和y(t)间的具有一定延时的信息传递函数，x0(t)代表由运动控制模块传至感觉反馈模块的兴奋信号，τ为延时因子，α(·)和γ(·)分别为α与γ感觉反馈模块元群的增益函数，αs,αd,γs,γd分别为静态和动态感觉反馈模块中的元群。

s12、构建感觉反馈通路模型部分：

s121、输入感觉反馈模块中的外界感觉信号y0(t)，通过中介体ia和ib传回运动控制模块，产生感觉反馈模块至运动控制模块的反馈传入信号z(t)，其具体的表达式为：

式中，l[·]为z(t)和y0(t)间的信息传递函数，l是牵张反射增益，s是交互兴奋增益，r是交互抑制增益，b是ib的增益，g是运动控制模块的中间抑制性元群的增益，v'和分别为来自于中介体ia和ib传入运动控制模块的放电频率，g'(t)是随机的已知向量；

s122、基于步骤s11中基于运动控制模块的函数关系，得到x2(t)与z(t)关系如下：

式中，hi为抑制性平均突触增益参数，εi为抑制性细胞膜平均时间常数和树突平均时间常数的和，z(t)为感觉反馈模块至运动控制模块的反馈传入信号，为运动控制模块中部分反馈信号状态量。

s13、综合步骤s11和步骤s12中运动控制模块与感觉反馈模块的作用关系，得到运动控制模块信号与感觉反馈模块信号的表达关系分别为：

式中，n1＝[αγ00]，n2＝[(αs,αd)(γs,γd)00]^t分别为未知参数向量，he为运动控制模块中局部兴奋性平均突触增益参数，εe为兴奋性细胞膜平均时间常数和树突平均时间常数的和，hi为运动控制模块中局部抑制性平均突触增益参数，εi为抑制性细胞膜平均时间常数和树突平均时间常数的和，s是交互兴奋增益，r是交互抑制增益，b是ib的增益，v'和分别为来自于中介体ia和ib传入运动控制模块的放电频率，α(·)和γ(·)分别为感觉反馈模块元群α与γ的增益函数；v(t)＝g[x0(t+τ)]分别为随机的已知向量。

s2、如图1所示，基于神经肌肉多层次耦合特征的参数辨识及优化。首先基于脑电和肌电信号建立多层次特征指标，并构建基于多层次耦合特征指标的目标函数；然后，基于无迹卡尔曼滤波(ukf)方法对神经肌肉信息交互模型进行参数辨识，获取体现运动控制系统的关键模型参量指标；最后，以脑电和肌电信号能量谱为目标，基于混合粒子群算法实现对所辨识参量进行优化选取，实现神经肌肉信息交互模型构建。

s21、基于多层次耦合特征指标构建参数辨识的目标函数：

s211、对采集的脑电和肌电信号进行预处理，并用非线性多尺度脑肌电同步耦合分析方法建立耦合强度vc、耦合方向tse和时间延迟td多层次特征指标；

s2111、基于时频一致性方法，获取耦合强度vc；

s2112、基于传递谱熵，获取耦合方向tse；

s2113、基于时延信息传递，获取时间延迟td；

s212、利用方程误差参数辨识方法构建基于多层次耦合特征指标的目标函数，目标函数的具体表达式为：

其中，ε(k)代表k时刻模型输出与系统输出的偏差，是模型参数估计值，是模型参数估计值下所有时刻的偏差和。

s22、基于无迹卡尔曼滤波(ukf)以及步骤s212构建的目标函数对神经肌肉信息交互模型进行参数辨识：

s221、令脑电和肌电信号表达式的输入信号u(t)服从高斯分布；

s222、基于步骤s211建立的多层次特征指标，利用增广状态向量构建状态方程，将模型中的未知参数向量定义为系统中新增状态向量，具体表示为：

x'(k)＝[x(k)h(k)ξ(k)](7)

其中，x为特征指标，h＝[m1m2]和ξ＝n1·n2分别表示式(5)中的运动控制模块参数和感觉反馈模块参数；

新的状态方程如下：

式中，x'(k+1)为n维随机状态向量序列，z(k+1)为n维系统可观测输出变量，包含x(t)和y(t)，f(·)和h(·)分别为n维和m维函数，u(k)为系统输入变量；

s223、基于无迹卡尔曼滤波方法以及步骤s222中建立的新的状态方程，并利用步骤212中的目标函数对未知参数h和ξ进行辨识；即基于ukf中采样策略逼近非线性分布的估计原则，当状态值x(k+1)满足时，辨识过程结束，此时模型的输出有效逼近基于输入的脑电和肌电信号计算的u(k)值，得到多层次耦合特征指标下的模型参数集合。由s1可知神经群参数h和肌纤维参数ξ与模型各个未知参量存在函数关系。

s23、基于混合粒子群算法对参量h＝[hvchtsehtd]及ξ＝[ξvcξtseξtd]进行最优选取，基于对输入的脑电和肌电信号控制系统神经肌肉功能耦合(fcmc)分析指标的逼近优化原则，以输入的脑电和肌电信号运动控制模块和感觉反馈模块中信号能量谱为目标。

在本发明的一个优选实施例中，运动控制模块中的元群信号x(t)和感觉反馈模块中的单元动作电位y(t)的变化分别与外界刺激信号有关。

以下结合实施例对本发明一种神经肌肉信息交互模型构建及参数辨识优化方法做进一步描述：

s1、构建神经肌肉信息交互模型。基于生物物理学理论，对信号的fcmc指标分析神经肌肉闭环通路的特性进行建模，将神经元群和肢体运动神经元模型有机整合，从神经元层面建立大脑与肌肉间的运动控制环路，即神经肌肉信息交互模型。神经肌肉信息交互模型主要分为运动控制通路建模和感觉反馈通路建模两部分。

s11、构建运动控制通路模型部分：

s111、根据输入的引起相应脑区兴奋神经元动作电位u变化的外界信号，获得与脑区兴奋神经元动作电位u相关的局部脑区兴奋性平均突出增益参数he、细胞膜平均时间常数和树突平均时间常数和εe的变化，通过上述相关参数得到基于神经元群的模型函数表达式如下：

式中，u(t)为模型输入信号，表示运动控制模块中动作电位的变化，x(t)代表运动控制模块中局部场电位，代表由运动控制模块传至感觉反馈模块的兴奋信号状态量；

s112、根据步骤s111知，大量神经元的电势变化将导致局部场电位x(t)(x1-x2的值)的改变，其中x1为兴奋性变化量，x2为抑制性反馈变化量，兴奋信号x0可以控制α与γ运动神经元募集，同时输出运动单元动作电位y(t)，其函数关系如下：

y(t)＝[α(αs,αd)+γ(γs,γd)]g[x0(t+τ)](2)

式中，g[·]为x(t)和y(t)间的具有一定延时的信息传递函数，x0(t)代表由运动控制模块传至感觉反馈模块的兴奋信号，τ为延时因子，α(·)和γ(·)分别为α与γ感觉反馈模块元群的增益函数，αs,αd,γs,γd分别为静态和动态感觉反馈模块中的元群。

外界刺激信号可以导致大脑神经元群信号x(t)及肌纤维运动单元动作电位y(t)的变化。

s12、构建感觉反馈通路模型部分：

s121、输入感觉反馈模块中的外界感觉信号y0(t)，通过纺锤体ia和高尔基腱器ib传入神经元传回大脑对动作调整，从而产生感觉反馈传入信号z(t)，其具体的表达式为：

式中，l[·]为z(t)和y0(t)间的信息传递函数，l是牵张反射增益，s是交互兴奋增益，r是交互抑制增益，b是ib的增益，g是运动控制模块的中间抑制性元群的增益，v'和分别为来自于中介体ia和ib传入运动控制模块的放电频率，g'(t)是随机的已知向量；

s122、基于步骤s11中基于运动控制模块的函数关系，得到x2(t)与z(t)关系如下：

式中，hi为抑制性平均突触增益参数，εi为抑制性细胞膜平均时间常数和树突平均时间常数的和，z(t)为感觉反馈模块至运动控制模块的反馈传入信号，为运动控制模块中部分反馈信号状态量。

s13、综合s11和s12中运动控制模块与感觉反馈模块的作用关系，得到运动控制模块信号与感觉反馈模块信号的表达关系分别为：

式中，n1＝[αγ00]，n2＝[(αs,αd)(γs,γd)00]^t分别为未知参数向量，he为运动控制模块中局部兴奋性平均突触增益参数，εe为兴奋性细胞膜平均时间常数和树突平均时间常数的和，hi为运动控制模块中局部抑制性平均突触增益参数，εi为抑制性细胞膜平均时间常数和树突平均时间常数的和，s是交互兴奋增益，r是交互抑制增益，b是ib的增益，v'和分别为来自于中介体ia和ib传入运动控制模块的放电频率，α(·)和γ(·)分别为感觉反馈模块元群α与γ的增益函数；v(t)＝g[x0(t+τ)]分别为随机的已知向量。

由上述建模过程可知，脑电仿真信号的输出既受外部输入的影响也受神经肌肉系统的自反馈信息的约束；肌电信号的输出受大脑外部输入的影响，也受神经肌肉系统各级单元的调制，并与大脑输入信息之间有一定的时延关系，可以通过对(5)式中的模型参数进行辨识估计，以实现外部刺激信号下神经肌肉闭环控制模型的有效构建。

s2、基于神经肌肉多层次耦合特征指标的参数辨识及优化。首先同步采集脑电和肌电信号，并进行预处理，再基于脑电和肌电信号建立多层次特征指标，构建基于多层次耦合特征的目标函数；然后，基于无迹卡尔曼滤波(ukf)方法对神经肌肉信息交互模型进行参数辨识，获取体现运动控制系统内在生理变化的关键模型参量指标；最后，以脑电和肌电信号能量谱为目标，基于混合粒子群算法实现对所辨识参量进行优化选取，具体辨识及优化策略如图2所示，进而实现神经肌肉信息交互模型的构建。

s21、基于多层次耦合特征指标构建参数辨识的目标函数：

s211、基于64导博睿康(neuracle)设备和10通道delsys肌电设备搭建脑电和肌电信号同步采集实验平台，实现多通道脑电和肌电信号同步采集。选取与运动及运动感觉相关额区、枕区、顶区以及颞区的64个脑电极位置。并对采集到的脑电和肌电数据进行预处理，基于自适应高通滤波进行基线校正，通过自适应陷波去除工频干扰，基于独立分量分析去眼电干扰。将预处理的数据用非线性多尺度脑肌电同步耦合分析方法建立耦合强度vc、耦合方向tse和时间延迟td多层次特征指标；

s2111、基于时频一致性方法，获取耦合强度vc；

s2112、基于传递谱熵，获取耦合方向tse；

s2113、基于时延信息传递，获取时间延迟td；

s212、利用方程误差参数辨识方法构建基于多层次耦合特征指标的目标函数，其基本思想是通过极小化模型输出与系统输出的误差来实现，目标函数的具体表达式为：

其中，ε(k)代表k时刻模型输出与系统输出的偏差，是模型参数估计值，是模型参数估计值下所有时刻的偏差和。

s22、基于无迹卡尔曼滤波(ukf)以及步骤s212构建的目标函数对神经肌肉信息交互模型进行参数辨识：

s221、为使神经信息交互模型及作用关系符合实际神经肌肉系统工作特性，令脑电和肌电信号表达式的输入信号u(t)服从高斯分布(e[u(t)]＝0,σ[u(t)]＝q)，建立外界刺激下神经肌肉控制系统仿真模型；

s222、基于步骤s211建立的多层次特征指标，利用增广状态向量构建状态方程，将模型中的未知参数向量定义为系统中新增状态向量，具体表示为：

x'(k)＝[x(k)h(k)ξ(k)](7)

其中，x为特征指标，h＝[m1m2]和ξ＝n1·n2分别表示式(5)中的神经群参数和肌纤维参数；

新的状态方程如下：

式中，x'(k+1)为n维随机状态向量序列，z(k+1)为n维系统可观测输出变量，包含x(t)和y(t)，f(·)和h(·)分别为n维和m维函数，u(k)为系统输入变量；

s223、基于无迹卡尔曼滤波方法以及步骤s222中建立的新的状态方程，并利用步骤212中的目标函数对未知参数h和ξ(与式(5)中m1和m2，n1和n2有关)进行辨识；

由步骤s1可知神经群参数h和肌纤维参数ξ与模型各个未知参量存在函数关系，因此基于ukf中采样策略逼近非线性分布的估计原则，对系统模型的未知参数h和ξ进行辨识，当状态值x(k+1)满足时，辨识过程结束，此时模型的输出有效逼近基于脑电和肌电信号计算的u(k)值，得到多层次耦合特征指标下的模型参数集合。

s23、基于混合粒子群算法对参量h＝[hvchtsehtd]及ξ＝[ξvcξtseξtd]进行最优选取，基于对神经肌肉控制系统fcmc分析指标的逼近优化原则，以脑电和肌电信号能量谱为目标，基于混合粒子群算法实现对神经群参数h＝[hvchtsehtd]及肌纤维参数ξ＝[ξvcξtseξtd]的最优选取，进而实现对正常神经肌肉信息交互模型的构建，并能获取体现大脑与肌肉信息交互内在微观表现的最优化特征指标h和ξ。

以上所述的实施例仅是对本发明的优选实施方式进行描述，并非对本发明的范围进行限定，在不脱离本发明设计精神的前提下，本领域普通技术人员对本发明的技术方案做出的各种变形和改进，均应落入本发明权利要求书确定的保护范围内。