一种神经元群模型的参数和状态估计方法与流程

本发明涉及神经信号处理技术领域，具体涉及一种神经元群模型的参数和状态估计方法。

背景技术：

大脑是一种由许多神经元群组成的复杂动力学系统，它具有非线性、不平稳性的特征。大脑皮层中包含大量细胞体、树突和轴突组成的神经元，这些神经细胞相互作用产生突触后电位的总和就是脑电信号。脑电图(electroencephalography,EEG)对大脑功能状态的变化有很高的灵敏度，能反映出大脑皮质中不同神经细胞的电生理活动。数学模型已广泛的应用于神经科学领域，并为该领域的探索与研究奠定了基础。神经元群模型的建立为大脑生理学机制的研究提供了便利条件。但是当用神经元群模型进行脑电信号仿真时，其中的生理学参数值是无法测量的，因此要想更深刻的理解脑电信号的机理，需要对脑电信号的生理学参数进行研究。

在构造每种神经元群模型时，其中的生理学参数是需要确定的，前人发现神经元群模型中的生理学参数是根据解刨学信息进行估计，通常采用经验值进行脑电信号仿真。之后有人利用优化算法中的遗传算法可以估计出模型生理学参数值，但是遗传算法编程过程比较复杂，其中的交叉概率和变异概率需要根据经验值确定，并且搜索速度较慢，要想获得较精确的最优解需要较长的运行时间。

技术实现要素：

(一)要解决的技术问题

本发明要解决的技术问题是：如何实现对神经元群模型的参数和状态的高效估计。

(二)技术方案

为了解决上述技术问题，本发明提供了一种神经元群模型的参数和状态估计方法，包括以下步骤：

将单通道脑电信号输入UKF中，对Jansen模型进行转换，从而估计神经元模型中的生理学参数：抑制性突触增益B、抑制性时间常数b、模型的输入均值μ，其他模型参数都被视为常数；被估计的三个参数用θ＝[B bμ]^T参数集表示；

UKF算法将Jansen模型定义为如下离散时间状态变量：

xt+1＝f(xt,θt)+εt (1)

yt＝Cxt+et (2)

其中，xt为Jansen模型中表示突触后膜电位的6个状态变量，f是非线性Jansen模型状态方程函数，yt代表脑信号的观测值，C为观测矩阵，θt为模型的参数集；公式(1)中的附加项εt表示模型输入高斯白噪声的标准差、et表示脑电信号观测噪声，即测量噪声，这两个参数都被定义为独立分布高斯过程，把这两个公式看作为常数项附加一个误差值，即一阶自回归过程，反映出参数变化比状态变化慢的特点；最终想要得到Jansen模型状态的变化和参数的变化，在UKF估计之前，将观测噪声添加到Jansen模型输出中，得到测量的脑电信号yt，然后再对模型进行参数估计；

用UKF算法估计状态变量，在每次更新数据时对模型预测进行校正，目标是估计模型状态xt和模型的参数集θt，因此将变量的增广矩阵定义为：

将需要估计的变量增广矩阵表示为：

其中，需要估计的参数集用表示，模型的六个状态变量的估计用表示，T表示向量的转置，UKF算法是一个时间更新过程，从t时刻推出t+1时刻状态，UKF算法通过确定后验增广状态估计找到时变参数估计，其中后验增广状态估计是t时刻神经元群模型输出yt的期望值，即：

UKF算法中的状态变量包括Jansen模型输出的六个状态变量和三个待估计的参数变量：抑制性突触增益B、抑制性时间常数b和模型的输入均值μ，因此，原来六维的非线性状态空间转换成九维的非线性增广状态空间，这里，UKF算法直接应用于九维增广状态空间模型，可同时估计模型的状态和参数，当对增广状态进行估计时，模型中其他参数的初始值都设为标准值。

优选地，六维的非线性状态空间转换成九维的非线性增广状态空间之后，状态变量在不同时刻t的UKF算法执行步骤如下：

(1)求采样点：获得2n个Sigma采样点集：

其中P表示变量的方差，n＝9；

(2)将点集代入方程：把2n个Sigma点集带入状态方程中，并分别计算这些点集的状态值Xt和预测值Yt；

Xt＝f(Xt-1,θt-1)+εt-1 (7)

Yt＝CXt+et (8)

(3)求估计值：求出Sigma点的状态估计和测量估计

(4)求协方差：根据状态估计和测量估计求出系统的协方差：

(5)求增益矩阵：根据系统的协方差求卡尔曼增益矩阵K：

(6)更新状态：根据增益矩阵更新系统的状态信息：

为更新后的9维状态空间变量，其中包括要估计Jansen模型的状态变量和参数变量，公式(11)中增广状态的协方差矩阵Q被定义为：

其中，α1＝10^-3，α2＝10^-7，分别表示三个需要估计参数的初始值，协方差矩阵Q是九阶对角矩阵，对角线上的元素依次是三个参数的方差和Jansen模型六个状态变量的方差，Δt(Aaσ)²是过程噪声，表示输入到Jansen模型中的高斯白噪声的方差，即公式(1)中的εt的平方；公式(12)中R表示观测噪声，恒等于(0.2σsim)²，即公式(2)中的et的平方，其中σsim是Jansen模型输出信号的标准差，从而UKF算法不断把协方差递归，通过减小模型产生的EEG信号yt和估计的模型输出信号之间的数据预测误差估计得到模型的参数。

优选地，表示需要估计参数的初始值：22mv。

优选地，表示需要估计参数的初始值：50s^-1。

优选地，表示需要估计参数的初始值：220。

优选地，σ为5.74。

本发明又提供了一种所述的方法在Wending模型中的应用，将所述方法应用于估计Wending模型的状态和生理学参数。

本发明还提供了一种所述的方法在Wending模型中的应用，将所述方法应用于估计丘脑模型中的状态和生理学参数。

(三)有益效果

本发明采用UKF算法对神经元群模型的参数和状态进行估计，其中，采用高效率的递归算法对神经元群模型的参数进行估计，利用滤波算法追踪模型中生理学参数的变化，为分析神经科学和临床应用中脑电信号潜在的生理机制提供可靠的依据。从仿真结果可以看出，在含有噪声的情况下，UKF算法可以估计出随时间变化的状态和生理学参数。UKF算法克服了遗传算法运行效率低、参数难以确定等缺陷，同时能估计出随时间变化的生理学参数。

具体实施方式

为使本发明的目的、内容、和优点更加清楚，下面结合实施例，对本发明的具体实施方式作进一步详细描述。

无迹卡尔曼滤波器算法(Unscented Kalman Filter，UKF)是一种非线性滤波方法，利用一定的滤波准则和统计量对系统进行估计，它是用无迹变换UT处理均值和协方差的非线性传递问题，摒弃了对非线性函数进行线性化的传统做法，是对非线性函数的概率密度分布进行近似，用一系列确定样本来逼近状态的后验概率密度。本发明就是采用UKF算法对神经元群模型的参数和状态进行估计。

下面以Jansen模型为例，利用UKF算法估计Jansen模型中的参数和状态。

将单通道脑电信号输入UKF中，对Jansen模型进行转换，从而估计神经元模型中的生理学参数：抑制性突触增益B、抑制性时间常数b、模型的输入均值μ，其他模型参数都被视为常数。被估计的三个参数用θ＝[B bμ]^T参数集表示。

UKF算法将Jansen模型定义为如下离散时间状态变量：

xt+1＝f(xt,θt)+εt (1)

yt＝Cxt+et (2)

其中xt为Jansen模型中表示突触后膜电位的6个状态变量，f是非线性Jansen模型状态方程函数，yt代表脑信号的观测值，C为观测矩阵，θt为模型的参数集。由于非线性系统动力学上的不确定性，公式(1)中的附加项εt表示模型输入高斯白噪声的标准差、et表示脑电信号观测噪声(测量噪声)，这两个参数都被定义为独立分布高斯过程，可以把这两个公式看作为常数项附加一个误差值，即一阶自回归过程，反映出参数变化比状态变化慢的特点。最终想要得到Jansen模型状态的变化和参数的变化。在UKF估计之前，将观测噪声添加到Jansen模型输出中，得到测量的脑电信号yt，然后再对模型进行参数估计。

用UKF算法估计状态变量，在每次更新数据时对模型预测进行校正，目标是估计模型状态xt和模型的参数集θt，因此这里将变量的增广矩阵定义为：

需要估计的变量增广矩阵表示为：

其中需要估计的参数集用表示，模型的六个状态变量的估计用表示(最后输出的参数估计值为这六个状态变量和这三个参数)，T表示向量的转置。UKF算法是一个时间更新过程，从t时刻推出t+1时刻状态。UKF算法通过确定后验增广状态估计找到时变参数估计，其中后验增广状态估计是t时刻神经元群模型输出yt的期望值，即：

UKF算法中的状态变量包括Jansen模型输出的六个状态变量和三个待估计的参数变量：抑制性突触增益B、抑制性时间常数b和模型的输入均值μ。因此，原来六维的非线性状态空间转换成九维的非线性增广状态空间。这里，UKF直接应用于九维增广状态空间模型，而不是原始的Jansen模型，并且相比于其他估计方法(例如最大似然估计)，UKF可以同时估计模型的状态和参数。当对增广状态进行估计时，模型中其他参数的初始值都设为标准值。状态变量在不同时刻t的UKF算法基本步骤如下：

(1)求采样点获得2n个Sigma采样点集：

其中P表示变量的方差；

(2)点集代入方程把2n个Sigma点集带入状态方程中，并分别计算这些点集的状态值Xt和预测值Yt。

Xt＝f(Xt-1,θt-1)+εt-1 (7)

Yt＝CXt+et (8)

(3)求估计值求出Sigma点的状态估计和测量估计

(4)求协方差根据状态估计和测量估计求出系统的协方差：

(5)求增益矩阵根据系统的协方差求卡尔曼增益矩阵K：

(6)更新状态根据增益矩阵更新系统的状态信息：

为更新后的9维状态空间变量，其中包括要估计Jansen模型的状态变量和参数变量，公式(11)中增广状态的协方差矩阵Q被定义为：

其中α1＝10^-3，α2＝10^-7，分别表示三个需要估计参数的初始值：22mv、50s^-1、220。协方差矩阵Q是九阶对角矩阵，对角线上的元素依次是三个参数的方差和Jansen模型六个状态变量的方差。Δt(Aaσ)²是过程噪声，表示输入到Jansen模型中的高斯白噪声的方差，即公式(1)中的εt的平方，这里σ等于5.74；公式(12)中R表示测量噪声，恒等于(0.2σsim)²，即公式(2)中的et的平方，其中σsim是Jansen模型输出信号的标准差，过程噪声和测量噪声一般假设不随系统状态变化而变化。简而言之，UKF算法是不断的把协方差递归，通过减小模型产生的EEG信号yt和估计的模型输出信号之间的数据预测误差从而估计模型的参数。

Jansen神经元群模型原理框图如图1所示，可以产生不同频率的窄带信号，体现单一细胞群的动力学特性。观测的脑电信号yt是由Jansen模型的输出加观测噪声构成。

对Jansen模型进行状态估计的结果如图2所示，图2a)表示Jansen模型输出信号的真实值和估计值变化曲线。截取49s～50s的信号，观察估计值的变化趋势，红色曲线表示Jansen模型的输出真实值，黑色曲线表示脑电信号估计值的分布情况即yt。图中2b)表示0s～100s附加噪声的Jansen模型的输出信号yt。由仿真结果可以看出，模型输出信号状态的估计值主要分布在真实值附近处，可以追踪随着时间变化的状态值。

然后用UKF算法估计随时间变化的参数值。首先将Jansen模型中的生理学参数抑制性突触增益B、抑制性时间常数b和模型的输入均值μ设置成随时间变化的函数，初始值设置为标准值B＝22mv、b＝50s^-1、μ＝220。UKF算法对Jansen模型的参数估计结果如图3所示，图3a)是将观测信号通过1阶500Hz的巴特沃斯滤波器，经过滤波得到的脑电信号，图3b)、3c)、3d)分别为抑制性突触增益B的估计、抑制性时间常数b的估计和模型的输入均值μ的估计，红线表示参数的真实值随时间变化的曲线，黑线表示根据UKF算法对参数的估计值曲线。由仿真结果可以看出，参数的估计值一直分布在真实值附近处，UKF算法可以估计出随时间变化的参数值。

本发明将UKF算法应用在三种不同的神经元群模型中，对Jansen模型中的状态和生理学参数进行估计，用同样的原理和步骤还可以估计以下两种神经元群模型，即Wending模型和丘脑模型中的状态和生理学参数。

Wending模型的原理框图如图4所示，利用UKF算法对Wending模型进行仿真。首先将Wending模型中的生理学参数：兴奋性突触增益A、慢速抑制性突触增益B和快速抑制性突触增益G设置成随时间变化的函数，初始值设置为标准值A＝3.25mv、B＝22mv、G＝10mv。利用UKF算法对Wending模型的状态估计如图5所示，图5a)表示Wending模型输出信号的真实值和估计值变化曲线。截取49s～50s的信号，观察估计值的变化趋势，红色曲线表示模型的输出信号真实值，黑色点线表示脑电信号估计值的分布情况即yt。由仿真结果可以看出，模型输出信号状态的估计值主要分布在真实值附近处。图中5b)表示0s～100s附加噪声的Wending模型的输出信号，即观测信号yt。

利用UKF算法对Wending模型的参数估计结果如图6所示，图6a)是将观测信号通过1阶500Hz的巴特沃斯滤波器，得到滤波后的脑电信号，图6b)、6c)、6d)分别为兴奋性突触增益A的估计、慢速抑制突触增益B的估计和快速抑制突触增益G的估计，红线表示真实的参数随时间的变化，黑线表示根据UKF算法对参数的估计值。由仿真结果可以看出，无迹卡尔曼滤波算法可以同时估计出随时间变化的三个生理学参数，但是生理学参数快速抑制突触增益G的估计值比真实值偏低，存在一定的误差。

丘脑模型的原理框图如图7所示。利用UKF算法对丘脑模型的状态估计结果如图8所示，图8a)表示丘脑模型输出信号的真实值和估计值变化曲线。截取49s～50s的仿真信号，观察估计值的变化趋势，红色曲线表示模型的输出信号真实值，黑色点线表示脑电信号估计值的分布情况即yt。由仿真结果可以看出，模型输出信号状态的估计值主要分布在真实值附近处。图中8b)表示是0s～100s附加噪声的丘脑模型的输出信号，即观测信号yt。

然后运用UKF算法对丘脑模型的参数进行估计。同样，首先将三个参数设置为随着时间变化的函数。对丘脑模型的参数估计结果如图9所示，图9a)是将观测信号通过1阶500Hz的巴特沃斯滤波器，得到滤波后的脑电信号，图9b)、9c)、9d)分别为兴奋性突触增益He的估计、抑制突触增益Hi的估计和抑制性时间常数τi的估计，红线表示真实的参数随时间的变化，黑线表示根据UKF算法对参数的估计值。由仿真结果可以看出，随着参数真实值的上升或下降，估计值也随之变化，UKF算法可以同时估计出随时间变化的三个生理学参数的分布，图中生理学参数抑制突触增益Hi的估计值比真实值偏低，存在一定的误差。

从仿真结果可以看出，在含有噪声的情况下，UKF算法可以估计出随时间变化的状态和生理学参数。UKF算法克服了遗传算法运行效率低、参数难以确定等缺陷，同时能估计出随时间变化的生理学参数，因此能够更深刻的理解脑电信号中参数的生理学意义。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明技术原理的前提下，还可以做出若干改进和变形，这些改进和变形也应视为本发明的保护范围。