自动化学报  2017, Vol. 43 Issue (2): 302-312   PDF    
无监督的猕猴运动皮层锋电位信号CKF解码
薛明龙, 吴海锋, 曾玉     
云南民族大学电气信息工程学院 昆明 650500
摘要: 如何通过猕猴运动皮层的神经元锋电位信号估计其手指移动位置是一神经解码问题,现存方法解决该问题大多采用有监督训练,需要通过训练数据得到神经元锋电位信号与手指移动位置的关系,因此其估计性能依赖于训练数据.本文提出了一种无监督解码方法,该方法基于状态空间模型(State space model,SSM),利用神经网络得到神经元锋电位数与手指移动位置的关系权值,再用逐次状态估计方法去估计手指移动的位置.为减少训练的复杂度和提高估计准确度,采用一种非线性的积分卡尔曼滤波(Cubature Kalman filtering,CKF)来完成神经网络的训练和手指位置的逐次状态估计.与传统方法相比,该方法的最大特点是无监督,可以由神经元锋电位簇向量直接估计手指移动位置,而无需有监督训练.实验结果显示,当采用较少的有监督数据,现存方法与本文方法相比有较大的估计误差;当采用较多的有监督数据,现存方法才具有与本文方法相近似的估计误差.
关键词: 神经解码     状态空间模型     无监督训练     积分卡尔曼滤波    
Unsupervised CKF Decoding for Macaque Motor Cortical Spikes
XUE Ming-Long, WU Hai-Feng, ZENG Yu     
School of Electrical and Information Technology, Yunnan Minzu University, Kunming 650500
Received: 2016-01-21, Accepted: 2016-05-31.
Foundation Item: Supported by National Natural Science Foundation of China (61262091), the 17th Batches of Young and Middle-aged Leaders in Academic and Technical Reserved Talents Project of Yunnan Province (2014HB019), the Project of Scientific Research Foundation of Yunnan Provincial Department of Education (2014Z093), and the Project of Postgraduate Innovation Foundation of Yunnan Minzu University (2016YJCXS03)
Author brief: XUE Ming-Long Master student at the School of Electrical and Information Technology, Yunnan Minzu University. His research interest covers machine learning and neural network;
ZENG Yu Lecturer at Yunnan Minzu University. Her main research interest is radio frequency identification (RFID) technology
Corresponding author. WU Hai-Feng Professor at Yunnan Minzu University. His main research interest is radio frequency identification (RFID) technology. Corresponding author of this paper
Recommended by Associate Editor TIAN Jie
Abstract: How to estimate macaque's moving finger position through neuron spikes in his mortor cortex is a problem about neural decoding. For the problem, most of existing methods use a supervised training algorithm and require supervised data to obtain the relationship between the spikes and the finger's moving position. Therefore, the performance of the existing methods depends on the training data. This paper proposes an unsupervised decoding method, which, based on a state space model (SSM), adopts neural networks to obtain the weights between the neuron spikes and the finger's moving position, and then estimates the finger's position through sequential state estimators. To reduce computational complexity and enhance estimation accuracy, a nonlinear cubature Kalman filter (CKF) is used to train the neural network and estimate the sequential moving positions. Compared with the existing methods, the proposed method's advantage is to be unsupervised. It could estimate the finger's position only through the spike vector instead of the supervised training data. Experiment results show that the existing methods have more estimation errors than the proposed method when a small amount of supervised data is adopted, and that the existing ones have similar estimation errors only when more supervised data adopted.
Key words: Neural decoding     state space model (SSM)     unsupervised training     cubature Kalman filter (CKF)    

神经编码是将外部世界映射到脑活动的一个过程[1], 外部世界经由视觉、听觉或嗅觉等感知信息传递到大脑, 大脑在得到这些信息后作出一系列反应, 这些反应体现在身体的骨骼或肌肉的活动中[2].可以用神经生理记录等分析仪器采集相应的脑活动信号, 然后将骨骼或肌肉的运动与脑活动建立一一对应关系, 从而完成神经编码[3].例如, 可以记录肢体运动时的神经元锋电位(Neuron spike), 然后找到肢体运动方向、速度与神经元锋电位的对应关系[4-8].神经解码是神经编码的逆过程, 它将脑活动映射到外部世界.例如, 由运动皮层中神经元锋电位的变化去估计手指运动的方向、速度[9-12], 通过辅助眼区的神经元锋电位估计眼睛移动的方向[13], 通过脑磁图(Magnetoencephalography, MEG) 信号估计手腕的移动[14].神经编码和解码的原理可用到生物医学工程等领域, 在有听力障碍的患者耳中植入人工耳蜗, 将声音信号编码为可刺激听觉神经的数字信号, 从而使患者具有感知声音的能力[15]; 用神经解码的方法让伤残人士通过大脑直接去控制设备的移动, 例如鼠标[10]或机械手臂等[16].

本文主要研究通过猕猴运动皮层中神经元锋电位数去估计其手指移动的位置, 这是一种典型的神经编码和解码问题.较早的编码和解码方法采用线性方法[17-18], 即神经脉元冲数与运动状态的关系为正比.线性方法最大的优点是计算简单, 但其估计的准确性较低.目前, 较为流行的是采用状态空间模型(State space model, SSM) 方法解决该问题. SSM模型已被广泛地运用到神经科学领域中, Czanner等建立了一个广义线性SSM模型, 这种模型比传统模型更能准确描述神经元的活动特性[19]; Malik等利用SSM模型, 采用EM (Expectation maximum) 算法准确估计神经元锋电位激发率(Fire rate) 函数[20]; Brown等通过SSM把老鼠海马区的神经锋电位序列(Spike train) 解码为位置信息[21]; Shanechi等建立了一个SSM模型, 并运用最优反馈控制模型来解码猴子运动的状态[22].

本文提出了一种猕猴运动皮层信号的无监督积分卡尔曼滤波解码方法(Unsupervised cabuture-Kalman-filtering decoding, UCKD). UCKD也是基于SSM模型, 对多个神经元锋电位数所构成的簇向量(Population vector)[23-24]进行解码.该解码方法采用神经网络得到神经元锋电位数与手指移动位置的关系权值, 同时采用逐次状态估计方法[25]估计手指移动的位置.为减少训练的复杂度和提高估计准确度, 采用一种非线性的积分卡尔曼滤波(Cubature Kalman filtering, CKF)[26]完成神经网络的训练和手指位置的逐次状态估计.与传统方法相比, 该方法的最大特点是无监督, 可以由记录的神经元锋电位簇向量直接估计手指移动的位置, 不需要先用有监督的训练数据对手指移动位置进行神经编码.当有监督的数据缺失或者不完备时, 只能依靠无监督神经解码.实验结果验证, 当采用较少的有监督数据, 现存的有监督方法与本文方法相比有较大的估计误差.

1 相关工作

由神经元锋电位簇向量估计猕猴手指移动的方法中, 较早提出的方法是线性拟合[17-18]可表示为

$\begin{align} {{\boldsymbol{Y}}_k} = {\boldsymbol{L}}({{\boldsymbol{s}}_k}) \end{align} $ (1)

其中, ${{\boldsymbol{Y}}_k} = {[{{y_{1k}}, {y_{2k}}, \cdots, {y_{Mk}}}]^{\rm T}}$为第${k \times \Delta t}$时刻猕猴的手指移动信息矢量, ${k = 1, 2, \cdots, K}$, ${\Delta t}$为采样时间, ${y_{mk}}$为猕猴第${m}$个手指移动信息, $m = 1, 2, \cdots, M$, 例如横坐标、纵坐标等信息, ${{\boldsymbol{s}}_k} =[{s_{1k}}$, ${s_{2k}}$, $\cdots $, ${s_{Nk}}]^{\rm T}$为神经元锋电位簇向量, ${s_{Nk}}$为第${N}$个电极所得到的神经元锋电位数, ${\boldsymbol{L}}( \cdot )$为一线性函数.

该方法先通过一组有监督的数据拟合出${\boldsymbol{L}}( \cdot )$, 然后将第${k}$个时刻的簇向量${\boldsymbol{s}_{k}}$代入${\boldsymbol{L}}( \cdot )$后估计出手指第${k}$个时刻的移动信息.线性拟合方法简单, 线性函数容易得到, 然而却存在以下问题.通过${\boldsymbol{L}}( \cdot )$求得反函数${{\boldsymbol{s}}_k} = {{\boldsymbol{L}}^{{ -1}}}({{\boldsymbol{Y}}_k})$时, 发现${\pmb {s}_{k}}$可能为负数, 但${\boldsymbol{s}_{k}}$只可能为$\geq 0$的数.

最大后验概率(Maximum a posteriori, MAP) 估计能避免该问题, 表示为[1]

$\begin{align} {{{ \hat{\pmb Y}}}_k} = \mathop {\max }\limits_{{{\boldsymbol{Y}}_k} \in {\mathbb{Y}^M}} \{ P({{\boldsymbol{s}}_k}|{{\boldsymbol{Y}}_k})P({{\boldsymbol{Y}}_k})\} \end{align}$ (2)

其中, $P({{\boldsymbol{s}}_k}|{{\boldsymbol{Y}}_k})$是在条件${{\boldsymbol{Y}}_k}$${\boldsymbol{s}_{k}}$的观测值的概率, $P({{\boldsymbol{Y}}_k})$${{\boldsymbol{Y}}_k}$的概率, ${{\mathbb{Y}^M}}$${{\boldsymbol{Y}}_k}$的搜索集合. MAP方法利用广义线性函数(General linear model, GLM)[1], 可以避免估计出来的${\boldsymbol{s}_{k}}$为负数的情况.然而同线形拟合一样, 最大后验概率估计将${{\boldsymbol{Y}}_k}$, ${k = 1, 2, \cdots, K}$看作为相互独立的数据.但实际上, 手指移动是一个具有一定速度和方向的连续过程, 特别当${\Delta t}$较小, 例如${10^{ - 2}}$ s时, ${{\boldsymbol{Y}}_k}$$ {{\boldsymbol{Y}}_{k - 1}}$是相关的.因此, 该方法的估计准确性较低.

SSM方法采用逐次状态估计的方法表示${{\boldsymbol{Y}}_k}$$ {{\boldsymbol{Y}}_{k -1}}$的相关性, 其状态方程和观测方程为[1, 9-12]

$\begin{align} {{\boldsymbol{Y}}_k} = {\boldsymbol{h}}({{\boldsymbol{Y}}_{k - 1}}) + {{\boldsymbol{\omega }}_k} \end{align}$ (3a)
$\begin{align} {{\boldsymbol{s}}_k} = {\boldsymbol{f}}({{\boldsymbol{Y}}_k}) + {{\boldsymbol{v}}_k} \end{align}$ (3b)

其中, ${\boldsymbol{h}}( \cdot )$表示自变量为${{\boldsymbol{Y}}_{k - 1}}$的函数, ${\boldsymbol{f}}( \cdot )$表示自变量为${{\boldsymbol{Y}}_k}$的函数, ${{\boldsymbol{\omega }}_k}$是过程噪声矢量, ${{\boldsymbol{v}}_k}$是观测噪声矢量.

通过式(3a) 和(3b) 估计${{\boldsymbol{Y}}_k}$, 需要知道函数${\boldsymbol{h}}( \cdot )$${\boldsymbol{f}}( \cdot )$, 若均为线性函数, 可以利用卡尔曼滤波方法估计${{\boldsymbol{Y}}_k}$ [11]; 若为非线性函数, 可以利用扩展卡尔曼滤波或粒子滤波方法[9].但是, 无论采用何方法估计${{\boldsymbol{Y}}_k}$, 必须知道这两个函数.传统的SSM方法需要通过有监督的训练数据来拟合得到${\boldsymbol{h}}( \cdot )$${\boldsymbol{f}}( \cdot )$, 因此其估计性能必然会受到拟合的影响.当拟合函数较精确, 其估计效果就会较好, 反之则会较差.本文提出一种不依赖有监督数据的无监督方法, 通过神经网络训练出神经元锋电位${\boldsymbol{s}_{k}}$和手指移动信息$ {{\boldsymbol{Y}}_k}$间的关系, 然后利用SSM估计出手指移动信息.

2 系统模型

猕猴的手指移动估计问题描述如下:首先对一只猕猴进行训练, 使其能够做"弹珠台"运动; 一旦目标物出现, 就用手指指向目标物. 图 1是猕猴的手指在一个长为L、宽为W的平面内的运动轨迹, 猕猴的手指初始位置在A点, 当目标物出现在B点, 猕猴的手指就移动到B点, 随后B点目标物消失, C点出现目标物, 猕猴的手指又移向C点, 再至D点, 反复向目标物出现的位置运动.同时给猕猴的脑部连接一组电极, 在${k \times \Delta t}$时刻, ${k = 1, 2, \cdots, K}$记录猕猴运动皮层的神经元锋电位数${\boldsymbol{s}_{k}}$, 然后通过该神经元锋电位数估计猕猴手指移动信息${{\boldsymbol{Y}}_k}$.采用SSM模型估计手指移动首先需要确定状态方程和观测方程.

图 1 猕猴手指移动轨迹图 Figure 1 An example for the trajectory of a macaque's finger
2.1 状态方程

以运动信息${{\boldsymbol{Y}}_k}$中手指移动纵坐标${y_k}$为例确定状态方程(3a).用${y_k}$代替${{\boldsymbol{Y}}_k}$后, 状态方程变为${y_k} = h({y_{k - 1}})+{\omega _k}$, 其中${h({y_{k - 1}})}$${\omega_k}$需要确定.由于纵坐标${y_k} \in [0,L]$, 则两个相邻时刻纵坐标间差值${{y_k} -{y_{k - 1}}}=\Delta y \in [{-L, L}$].将$[-L, L]$等分成${2L}$个区间${{\cal L}_1}$, ${{\cal L}_2}, \cdots, {{\cal L}_{2L}}$, 并让${\Delta t}$足够小, 例如${10^{ - 2}}$ s, 统计出${\Delta y}$落在各个区间的概率p $(\Delta y \in {{\cal L}_i})$, ${i}= 1, 2, \cdots, {2L}$. 图 2给出了其概率分布图, 数据来源于公开数据(详见第3.1节系统设置), 其中${\Delta t}=70\, {\rm ms}$, ${L}=15\, {\rm cm}$.由图 2可以看出, ${\Delta y_k}$小于0.5 cm的概率约为60 %, 小于2.5 cm的概率约为97 %, 其概率分布曲线近似于零均值的高斯分布, 由此可得到相应的状态方程为

$\begin{align} {y_k} = {y_{k - 1}} + {\omega _k} \end{align}$ (4a)

其中, ${\omega _k}$是方差为${Q_\omega }$、均值为0的高斯分布噪声.

图 2 相邻两时刻纵坐标间差值${\Delta y}$的分布概率图 Figure 2 An example for the distribution of ${\Delta y}$
2.2 测量方程

同样, ${y_k}$代替${{\boldsymbol{Y}}_k}$后, 观测方程(3b) 变为

$\begin{align} {{\boldsymbol{s}}_k}={\boldsymbol{f}}({y_k})+{{\boldsymbol{v}}_k}\end{align}$ (4b)

其中${\boldsymbol{f}}({y_k})$${{\boldsymbol{v}}_k}$需要确定.由文献[1, 17-18]可知, ${\boldsymbol{f}}( \cdot )$近似为一种线性函数.

将式(4a) 代入${{\boldsymbol{s}}_k} = {\boldsymbol{f}}({y_k}) + {{\boldsymbol{v}}_k}$中, 有

$\begin{align} {{\boldsymbol{s}}_k} - {{\boldsymbol{s}}_{k - 1}} = {\boldsymbol{\gamma }}{\omega _k} + {{\boldsymbol{v}}_k} - {{\boldsymbol{v}}_{k - 1}} \end{align}$ (5)

其中, ${\boldsymbol{\gamma}}$为一常数向量.令第${n}$个电极上两相邻时刻神经元锋电位数差值$\Delta \boldsymbol{s}_n=s_{nk} - s_{n(k - 1)}$, 则$\Delta {\boldsymbol{s}_n} \in$ $[-S, S]$, 其中S为电极所能采集到的最大锋电位数.将$[-S, S]$等分成$[2S]$个区间${\cal S}_1, {\cal S}_2, \cdots, {\cal S}_{2S}$, 统计出${\Delta \boldsymbol{s}}$落在各个区间的概率p $(\Delta {\boldsymbol{s}_n}\in {{\cal S}_i})$, ${i}= 0, 1, 2, \cdots, {2S}$. 图 3给出了其概率分布图, 数据来源于公开数据(详见第3.1节系统设置), 其中${\Delta t} = 70$ ms.从图 3可以看出, $\Delta{\boldsymbol{s}_n}$为0的概率约为30 %, $\Delta{\boldsymbol{s}_n}$为1或-1的概率约为80 %, 其概率分布曲线也近似于零均值的高斯分布.由于${\omega_k}$$\Delta{\boldsymbol{s}_n}$均近似为零均值高斯分布, 则由式(5) 可知, ${{\boldsymbol{v}}_k}$也应是零均值的高斯分布的噪音矢量.而且, 由于各电极上各神经元锋电位数相互独立, 因此各电极上的观测噪声也相互独立.最终, 得到观测方程${\boldsymbol{f}}({y_k}) + {{\boldsymbol{v}}_k}$中, 有$ {{\boldsymbol{s}}_k} = {\boldsymbol{f}}({y_k}) + {{\boldsymbol{v}}_k} $.其中, ${\boldsymbol{f}}( \cdot )$为一线性函数, $[{v_{1k}}, {v_{2k}}, \cdots, {v_{Nk}}]^{\rm{T}}$${\boldsymbol{0}}$均值、方差矩阵为${{ {R}}_v}$的高斯分布噪声矢量, ${v_{nk}}$为零均值、方差为$\sigma _{nk}^2$的高斯分布噪声, ${{\boldsymbol{R}}_v} = {\rm diag}\{[\sigma _{1k}^2$, $\sigma _{2k}^2$, $\cdots, \sigma _{Nk}^2]^{\rm{T}}\}$.

图 3 相邻两时刻纵坐标间差值${\Delta s}$的分布概率图 Figure 3 An example for the distribution of ${\Delta s}$

通过状态方程(4a) 和观测方程(4b) 估计猕猴手指移动的位置${y_k}$, 需要知道函数${\boldsymbol{f}}( \cdot)$的具体表达形式.传统的方法通过有监督的训练方法, 即知道神经元锋电位数目${{\boldsymbol{s}}_k}$和对应的猕猴手指移动位置${y_k}$, 然后通过该数据拟合出${\boldsymbol{f}}( \cdot)$.本文提出一种UCKD方法, 不需要监督数据就能估计出${y_k}$.

3 无监督的猕猴运动皮层信号的CKF解码 3.1 权值训练与位置估计

由上一节所述, 估计手指移动的纵坐标$y_k $需要先得到函数${{\pmb f}}(\cdot )$, 为此, 将式(4b) 变为

$\begin{align} {{\pmb s}}_k ={{\pmb f}}({{\pmb y}}_k, {{\pmb w}}_k )+{{\pmb v}}_k \end{align}$ (6)

其中, ${ {\pmb a}}_k =[a_{1k}, a_{2k}, \cdots, a_{Nk}]^{\rm T}$, ${ {\pmb b}}_k=[b_{1k}, b_{2k}, \cdots$, $b_{Nk}]^{\rm T}$, ${ {\pmb y}}_k =[y_k\ 1]^{\rm T}$, ${ {\pmb f}}(\cdot )$表示${ {\pmb s}}_k $$y_k $成线性关系, 为

$\begin{align} \label{eq2} {{\pmb f}}({{\pmb y}}_k, {{\pmb w}}_k )={{\pmb a}}_k \times y_k +{{\pmb b}}_k \end{align}$ (7)

将式(7) 视为单层感知器的神经网络, 如图 4所示, 其中, 激活函数$\varphi ({{\pmb s}}_k )={{\pmb s}}_k $.利用SSM方法训练权值向量${{\pmb w}}_k $, 于是给出一组关于${{\pmb w}}_k $的状态空间方程

$\begin{align} \label{eq3}&{ { {\hat{\pmb w}}}}_k ={ { {\hat{\pmb w}}}}_{k-1} +{{\pmb \theta }}_{k} \end{align}$ (8a)
$\begin{align} \label{eq4}& {{\pmb s}}_k ={{\pmb f}}( {\hat{y}}_k, { { {\hat{\pmb w}}}}_k )+{{\pmb v}}_k \end{align}$ (8b)

其中, ${{\pmb \theta }}_{k}$为权值噪声, 假定其均值为${{\pmb 0}}$、方差矩阵为${ { R}}_w $, $ {\hat{y}}_k $${\rm { {\hat{\pmb w}}}}_k $$y_k $${{\pmb w}}_k $的估计值.得到状态空间方程(8) 后, 采用逐次状态估计方法来训练权值${{\pmb w}}_k$, 从而得到${{\pmb f}}(\cdot )$.

图 4 权值训练神经网络信号流 Figure 4 Weight training neural network signals flow

值得注意的是, 通过式(8) 训练${\rm { {\hat{\pmb w}}}}_k $需要知道$ {\hat{y}}_k $的信息, 我们仍然通过一组状态空间方程求解$ {\hat{y}}_k $.由式(4a) 和式(6) 构成状态方程和观测方程, 并用$ {\hat{y}}_k $代替$y_k $, 有以下状态空间方程

$\begin{align} \label{eq5}&{\hat{y}}_k = {\hat{y}}_{k-1} +\omega _k \end{align}$ (9a)
$\begin{align} \label{eq6} &{{\pmb s}}_k ={{\pmb f}}( {\hat{y}}_k, {\hat {\pmb w}}_k )+{{\pmb v}}_k \end{align}$ (9b)

同样, 由式(9) 估计$y_k $也需要${\rm { {\hat{\pmb w}}}}_k $. UCKD采用的卡尔曼类的逐次状态方法(详见第4.2节), 其实是在高斯分布的噪声下, 使得条件概率期望值${\rm E}[p(y_k$, ${{\pmb w}}_k \vert \omega _k, {{\pmb v}}_k )]$最大.其过程是, 先固定状态$ {y}_k $, 求解权值${{\pmb {w}}}_k $, 然后再固定权值${\rm { {\hat{\pmb w}}}}_k $, 求解状态$ {y}_k $, 反复迭代直到收敛, 如图 5所示.该过程与EM算法类似[27], 由Jensen不等式可知, 若$p(y_k, {{\pmb w}}_k \vert\omega _k, {{\pmb v}}_k )$为凸函数, 那么有${\rm E}[p(y_k, {{\pmb w}}_k \vert\omega _k, {{\pmb v}}_k )]\ge p[{\rm E}(y_k, {{\pmb w}}_k \vert \omega _k $, ${{\pmb v}}_k )]$, 经过反复训练, ${\rm E}[p(y_k, {{\pmb w}}_k \vert\omega _k, {{\pmb v}}_k )]$将接近其下限值, 致使算法达到收敛.

图 5 位置估计与权值训练 Figure 5 Position estimation and weight training
3.2 积分卡尔曼滤波

通常采用逐次状态估计的方法求解权值和位置信息[25], 例如线性卡尔曼滤波(Kalman filter, KF)、扩展卡尔曼滤波(Extended Kalman filter, EKF) 和粒子滤波(Particle filter, PF) 等.然而, KF要求方程涉及的函数均为线性函数, 而${{\pmb f}}(\cdot )$仅是关于$y_k $的线性函数而非${{\pmb w}}_k $; EKF可处理非线性函数, 但它是一阶泰勒展开式近似, 且需求导数, 复杂度较高; PF也能处理非线性函数, 但需要大量粒子, 耗时较长.

本文采用CKF[26]求解式(8) 的权值训练和式(9) 的位置估计, 同上述方法相比, CKF不仅可以处理非线性过程, 而且保留了函数的二阶信息, 同时无需求导, 因此具有较高的近似度和较快的计算时间. CKF基于三阶球面径向积分法则, 近似地求解高斯非线性函数的积分, 例如, 对于一个n维高斯分布向量${{\pmb x}}\in { {\bf R}}^n$, 其近似的积分规则可表示为[26]

$\int_{{{\bf R}}^n} {{{\pmb h}}({{\pmb x}})} {\cal N}\left( {{{\pmb x}}; {{\pmb \mu }}, {{\pmb \Sigma }}} \right){\rm d {\pmb x}}\approx \frac{1}{2n}\sum\limits_{i=1}^{2n} {{{\pmb h}}({{\pmb \mu }}+\sqrt {\rm { \Sigma }} {{\pmb \xi }}_i )}$

其中, ${{\pmb h}}({{\pmb x}})$表示关于向量${{\pmb x}}$的函数, ${\cal N}\left( {{\boldsymbol{x}}; {\boldsymbol{\mu }}, {\boldsymbol{\Sigma }}} \right)$表示为均值为${{\pmb \mu }}$、方差矩阵为${ { \Sigma }}$n维高斯分布密度函数, $2n$个数值积分点为

${\xi _i} = {\begin{cases} \sqrt n \, {{\boldsymbol{e}}_i}, &i = 1, 2, \cdots, n \\[1mm] - \sqrt n\, {{\boldsymbol{e}}_{i - n}}, &i = n + 1, n + 2, \cdots, 2n \end{cases}}$

${{\pmb e}}_i \in {\bf{R}}^{n}$表示为n维单位向量中的第i个列向量, 即除第i元素为1, 其余为0.

算法1和算法2分别给出了用CKF求解权值${{\pmb w}}_k $和纵坐标$y_k $的算法过程, 均包含了时间更新和观测更新两个步骤.在算法中引入了符号$k\vert k-1$$k\vert k$. $ {y}_{k\vert k-1} $表示给定观测值${{\pmb s}}_1 $, ${{\pmb s}}_2, {\cdots}, {{\pmb s}}_{k-1} $下对$y_k $的估计值.由此, $y_k $的最终估计值表示为$ {y}_k = {y}_{k\vert k} $.算法中以下参数需要确定:

1) 遗忘因子$\lambda $, 取值范围10$^{-2}$~10$^{-5}$;

2) ${{\pmb {w}}}_{k\vert k} $的初值${{\pmb {w}}}_0 $, 可在-0.5~0.5之间随机产生;

3) $ {y}_{k\vert k} $的初始值$ {y}_0 $, 为0~1之间的随机数;

4) 状态噪声$\omega _k $的方差$Q_\omega $, 取值范围为10$^{-1} \sim$ 10$^{0}$;

5) 状态噪声${{\pmb \theta }}_k $的方差矩阵${ { R}}_w $, 取值为$\delta {{ I}}_{2N\times 2N} $, 其中$\delta $约为10$^{-1}$ $\sim $ 10$^{-2}$;

6) 观测噪声矢量${{\pmb v}}_k $的方差矩阵${ { R}}_v $, 由式(5) 可近似为${{\pmb s}}_k -{{\pmb s}}_{k-1} $的方差, 有

$\begin{align} { R}_v \approx{\rm diag}\left\{ \sum\limits_{k=1}^{K-1} \dfrac{({{\pmb s}}_{k+1} -{{\pmb s}}_k )^2}{K-1}\right\} \end{align}$ (10)

算法 1.  训练${\boldsymbol{w}_k}$的CKF算法

初始条件.  ${{\boldsymbol{P}}_{w, 0|0}} =\delta{ I}$, ${ I}$为单位矩阵, $\delta$约为$10^{ - 1}\;\sim 10^{ - 2}$, ${{{ \hat{\pmb w}}}_{0|0}} = {{{ \hat{\pmb w}}}_0}$

步骤 1.  时间更新

步骤 1.1.  分解

$\begin{align}{{\boldsymbol{P}}_{w, k|k}} = {Z_{w, k|k}}Z_{w, k|k}^{\rm T} \end{align}$

步骤 1.2.  计算数值积分点$(i=1, 2, \cdots, 2n)$

${{\boldsymbol{W}}_{i, k{\rm{|}}k}} = {Z_{w, k|k}}{\xi _i} \cdot {\lambda ^{ - \frac{1}{2}}} + {{ \hat{{\pmb w}}}_{k|k}}$

步骤 1.3.  计算传播积分点$(i=1, 2, \cdots, 2n)$

${\boldsymbol{W}}_{i, k + 1{\rm{|}}k}^ * = {{\boldsymbol{W}}_{i, k{\rm{|}}k}} $

步骤 1.4.  计算预测状态

${{{ \hat{\pmb w}}}_{k + 1|k}} = \dfrac{1}{2}{{n}}\sum\limits_{i = 1}^{2n} {{\boldsymbol{W}}_{i, k + 1{\rm{|}}k}^ * }$

步骤 1.5.  计算预测错误协方差

$\begin{align} {{\boldsymbol{P}}_{w, k + 1|k}} = &\ \frac{{1}}{2}{{n}}\sum\limits_{i = 1}^{2n} {{\boldsymbol{W}}_{i, k + 1{\rm{|}}k}^ * } {\boldsymbol{W}}_{i, k + 1{\rm{|}}k}^{ * {\rm T}}\, -\\ &\ {{{ \hat{\pmb w}}}_{k + 1|k}}{{ \hat{\pmb w}}}^{\rm T}_{k + 1|k} + {{ {R}}_w} \end{align}$

步骤 2.  测量更新

步骤 2.1.  分解

${{\boldsymbol{P}}_{w, k + 1|k}} = {Z_{w, k + 1|k}}Z_{w, k + 1|k}^{\rm T}$

步骤 2.2.  计算数值积分点$(i=1, 2, \cdots, 2n)$

${{\boldsymbol{W}}_{i, k + 1{\rm{|}}k}} = {Z_{w, k + 1|k}}{\xi _i} + {{{\hat{\pmb w}}}_{k + 1|k}}$

步骤 2.3.  计算传播积分点$(i=1, 2, \cdots, 2n)$

${{\boldsymbol{S}}_{i, k + 1{\rm{|}}k}} = {\mathop{\boldsymbol{f}}\nolimits} ({{\boldsymbol{Y}}_{k + 1{\rm{|}}k + 1}}, {{\boldsymbol{W}}_{i, k + 1{\rm{|}}k}})$

步骤 2.4.  计算预测观测值

${{\hat{ \boldsymbol{s}}}_{k + 1|k}} = \frac{{1}}{2}{{n}}\sum\limits_{i = 1}^{2n} {{{\boldsymbol{S}}_{i, k + 1{\rm{|}}k}}} $

步骤 2.5.  计算新息协方差矩阵

$\begin{align*} {{\boldsymbol{P}}_{SS, k + 1|k}} = &\ \frac{{1}}{2}{{n}}\sum\limits_{i = 1}^{2n} {{\boldsymbol{S}}_{i, k + 1|k}^{} {\boldsymbol{S}}_{i, k + 1|k}^{\rm T}} \, -\\ &\ {{{ \hat{\pmb s}}}_{k + 1|k}}{{{ \hat{\pmb s}}}^{\rm T}}_{k + 1|k} + {{ {R}}_v} \end{align*}$

步骤 2.6.  计算互协方差

$\begin{align} {{\boldsymbol{P}}_{ws, k + 1|k}} = \frac{1}{2}{{n}}\sum\limits_{i = 1}^{2n} {{\boldsymbol{W}}_{i, k + 1|k}^{}{\boldsymbol{S}}_{i, k + 1|k}^{\rm T}} - {{{ \hat{\pmb w}}}_{k + 1|k}}{{{\hat{\pmb s}}}^{\rm T}}_{k + 1|k}\end{align}$

步骤 2.7.  计算卡尔曼增益

${\boldsymbol{G}_{w,k + 1}} = {\boldsymbol{P}_{ws,k + 1|k}}\boldsymbol{P}_{ss,k + 1|k}^{ - 1}$

步骤 2.8.  计算滤波状态值

${{{\hat{\pmb w}}}_{k + 1|k + 1}} = {{{ \hat{\pmb w}}}_{k + 1|k}} + {{\boldsymbol{G}}_{w, k + 1}}({{\boldsymbol{s}}_{k + 1}} - {{\boldsymbol{ s}}_{k + 1|k}})$

步骤 2.9.  计算新息协方差矩阵

${{\boldsymbol{P}}_{w, k + 1|k + 1}} = {{\boldsymbol{P}}_{w, k + 1|k}} - {{\boldsymbol{G}}_{w, k + 1}}{{\boldsymbol{P}}_{ss, k + 1|k}}{\boldsymbol{G}}_{w, k + 1}^{\rm T}$

算法 2.  估计${\pmb y_k}$的CKF算法

初始条件.${P_{y, 0|0}}=1$, ${ y_{0|0}} = { y_0}$

步骤 1.  时间更新

步骤 1.1.  分解

${P_{y, k|k}} = {Z_{y, k|k}}Z_{y, k|k}^{\rm T}$

步骤 1.2.  计算数值积分点$(i=1, 2, \cdots, 2n)$

${{\boldsymbol{Y}}_{i, k{\rm{|}}k}} = {Z_{y, k|k}}{\xi _i} + {\hat{ y}_{k|k}}$

步骤 1.3.  计算传播积分点$(i=1, 2, \cdots, 2n)$

${\boldsymbol{Y}}_{i, k + 1{\rm{|}}k}^ * = {{\boldsymbol{Y}}_{i, k{\rm{|}}k}}$

步骤 1.4.  计算预测状态

${ \hat{y}_{k + 1|k}} = \frac{{1}}{2}{{n}}\sum\limits_{i = 1}^{2n} {{\boldsymbol{Y}}_{i, k + 1{\rm{|}}k}^ * }$

步骤 1.5.  计算预测错误协方差

$\begin{align*} {P_{y, k + 1|k}} =&\ \frac{{1}}{2}{{n}}\sum\limits_{i = 1}^{2n} {{\boldsymbol{Y}}_{i, k + 1{\rm{|}}k}^ * } {\boldsymbol{Y}}_{i, k + 1{\rm{|}}k}^{ *{\rm T}}\, -\\ &\ { \hat{y}_{k + 1|k}} \hat{y}^{\rm T}_{k + 1|k} + {Q_\omega } \end{align*}$

步骤 2.  测量更新

步骤 2.1.  分解

${P_{y, k + 1|k}} = {Z_{y, k + 1|k}}Z_{y, k + 1|k}^{\rm T}$

步骤 2.2.  计算数值积分点$(i=1, 2, \cdots, 2n)$

${{\boldsymbol{Y}}_{i, k + 1{\rm{|}}k}} = {Z_{y, k + 1|k}}{\xi _i} + { \hat{y}_{k + 1|k}}$

步骤 2.3.  计算传播积分点$(i=1, 2, \cdots, 2n)$

${{\boldsymbol{S}}_{i, k + 1{\rm{|}}k}} = {\boldsymbol{f}}({{\boldsymbol{Y}}_{i, k + 1{\rm{|}}k}}, {{\boldsymbol{w}}_{k + 1{\rm{|}}k}})$

步骤 2.4.  计算预测观测值

${{\boldsymbol{ s}}_{k + 1|k}} = \frac{{1}}{2}{{n}}\sum\limits_{i = 1}^{2n} {{{\boldsymbol{S}}_{i, k + 1{\rm{|}}k}}}$

步骤 2.5.  计算新息协方差矩阵

$\begin{align*}{{\boldsymbol{P}}_{SS, k + 1|k}} =&\ \frac{{1}}{2}{{n}}\sum\limits_{i = 1}^{2n} {{\boldsymbol{S}}_{i, k + 1|k}^{} {\boldsymbol{S}}_{i, k + 1|k}^{\rm T}} \, -\\ &\ {{{ \hat{\pmb s}}}_{k + 1|k}}{{{ \hat{\pmb s}}}^{\rm T}}_{k + 1|k} + {{ {R}}_v} \end{align*}$

步骤 2.6.  计算互协方差

$\begin{align*}{{\boldsymbol{P}}_{ys, k + 1|k}} =&\ \frac{{1}}{2}{{n}}\sum\limits_{i = 1}^{2n} {{\boldsymbol{Y}}_{i, k + 1|k}^{} {\boldsymbol{S}}_{i, k + 1|k}^{\rm T}}\, -\\ &\ { \hat{y}_{k + 1|k}}{{{ \hat{\pmb s}}}^{\rm T}}_{k + 1|k}\end{align*}$

步骤 2.7.  计算卡尔曼增益

${{\boldsymbol{G}}_{y, k + 1}} = {{\boldsymbol{P}}_{ys, k + 1|k}}{\boldsymbol{P}}_{ss, k + 1|k}^{ - 1}$

步骤 2.8.  计算滤波状态值

${ \hat{y}_{k + 1|k + 1}} = { \hat{y}_{k + 1|k}} + {{\boldsymbol{G}}_{y, k + 1}}({{\boldsymbol{s}}_{k + 1}} - {{{\hat{\pmb s}}}_{k + 1|k}})$

步骤 2.9.  计算新息协方差矩阵

${P_{y, k + 1|k + 1}} = {P_{y, k + 1|k}} - {{\boldsymbol{G}}_{y, k + 1}}{{\boldsymbol{P}}_{ss, k + 1|k}}{\boldsymbol{G}}_{y, k + 1}^{\rm T} $
3.3 正则化位置估计

通过状态空间方程(8a), (8b) 和(9a), (9b) 求解纵坐标$y_k $并不具有唯一解, 假定${ \hat {y}}'_k $也为$y_k $的一个估计值, 且与期望值$y_k $成线性关系, 令

$\begin{align} \label{eq11} \hat {y}_k =c{ \hat {y}}'_k +d \end{align}$ (11)

其中, cd为任意常数.将式(11) 代入式(8a), 式(8b) 和式(9a), 式(9b) 中, 有

$\begin{align} \label{eq12} &{ { { \hat{\pmb w}}'}}_k ={\rm { { \hat{\pmb w}}'}}_{k-1} +{{\pmb \theta }}_{k} \end{align}$ (12a)
$\begin{align} \label{eq13} &{{\pmb s}}_k ={{\pmb f}}({ \hat{y}}'_k, {\rm { { \hat{\pmb w}}'}}_k )+{{\pmb v}}_k \end{align}$ (12b)
$\begin{align} \label{eq14} &{ \hat {y}}'_k ={ \hat {y}}'_{k-1} +\omega _k \end{align}$ (13a)
$\begin{align} \label{eq15} &{{\pmb s}}_k ={{\pmb f}}({ \hat{y}}'_k, {\rm { { \hat{\pmb w}}'}}_k )+{{\pmb v}}_k \end{align}$ (13b)

其中, ${\pmb {w}}'_k ={[{{\pmb {a}}'_k}^{\rm T} \ \ {\pmb {b}'_k}^{\rm T}]}^{\rm T}$, ${\pmb {a}}'_k =c\, {\pmb a}_k$, $\pmb {b}'_k =d\, {\pmb {a}}'_k +{{\pmb b}}_k$, ${{\pmb f}}(\cdot )$表示为

$\begin{align} {{\pmb f}}( \hat {\pmb y}'_k, \hat{\pmb w}'_k )=\pmb {a}'_k \times { {y}}'_k +\pmb {b}'_k \end{align}$ (14)

将新的状态空间方程(12a), (12b) 和(13a), (13b) 与原状态空间方程(8a), (8b) 和(9a), (9b) 相比, 发现两者具有相同的形式.由此可知, 只要与期望值$y_k $成线性关系的均是其估计值, 它并不具有唯一解.注意到式(11)~(14) 给出了必要性的推导, 没有给出充分性, 即如果是式(8a), 式(8b) 和式(9a), 式(9b) 的解, 那么一定与期望值成线性关系.但是, 经多次实验结果显示, 通过UKCD得到的解均与期望值成线性关系.为限定估计范围, 将其规定在$-1$~1之间, 这个带有约束条件的估计问题可表示为

$\begin{align} \label{eq7} \hat{y}_k =u({{\pmb s}}_k ) \end{align}$ (15)
$\begin{align} \, {\rm s.t. }\ \ -1\le {y}_k \le 1 \end{align}$ (16)

其中, $u(\cdot )$函数表示为通过式(8a), 式(8b) 和式(9a), 式(9b) 以及算法1和算法2估计$y_k $.为满足式(15), 在每一次估计结果$y_k $上加一个双曲正切的激活函数, 使其不会超出区间[$-1$, 1], 即

$\begin{align} \hat{y}_k =\tanh ( {y}_k ) \end{align}$ (17)

由上所述, 估计值$ \hat{y}_k $与期望值$y_k $成线性关系, 对$ \hat {y}_k $做一线性变化可得其实际状态$\tilde {y}_k $, 表示为

$\begin{align} \tilde {y}_k =\frac{\left( {y_{\max } -y_{\min } } \right)\left( { \hat{y}_k - \hat{y}_{\min } } \right)}{ \hat{y}_{\max } - \hat{y}_{\min } }+y_{\min } \end{align}$ (18)

其中, $y_{\max } $$y_{\min } $为猕猴手指移动位置的最大值和最小值, 由于猕猴手指的移动范围已经被限定, 因此通常为已知, 例如$y_{\max } =L$, $y_{\min } =0$; $ \hat{y}_{\max } $$ \hat{y}_{\min } $为估计值$ \hat{y}_k $, $k=1, 2, \cdots, K$中的最大值和最小值, 若K足够大, 则$ \hat{y}_k $遍历, 那么$ \hat{y}_{\max } =1$, $ \hat{y}_{\min } = -1$.

3.4 平滑

在估计的结果中, 某些区间的数据中包含较多的高频分量.为提高估计性能, 采用五点三次移动平均的平滑方法[28]消除高频分量对估计结果的影响.

3.5 UCKD算法

算法 3.   UCKD算法

给定观测值: ${ {{{\pmb s}}_1, {{\pmb s}}_2, \cdots, {{\pmb s}}_K}}$

步骤 1.  初始化${\rm { {\hat{\pmb w}}}}_{0\vert 0} ={\rm {{\hat{\pmb w}}}}_0 $$ \hat {y}_{0\vert 0} = \hat {y}_0$, 并确定各噪声方差$Q_\omega $, ${ { R}}_w $${ { R}}_v$.

步骤 2.  令$k=0, 1, \cdots, K$

步骤 2.1.  由${ { {\hat{\pmb w}}}}_{k\vert k} $$ \hat{y}_{k\vert k} $, 通过式(9) 和算法2计算$ \hat{y}_{k+1\vert k+1} $;

步骤 2.2.  由$ \hat{y}_{k+1\vert k+1} $${\rm { {\hat{\pmb w}}}}_{k\vert k} $, 通过式(8) 和算法1计算${ { {\hat{\pmb w}}}}_{k+1\vert k+1} $;

步骤 2.3.  使$ \hat{y}_{k+1\vert k+1}=\tanh (\hat{y}_{k+1\vert k+1})$, 并令$ \hat{y}_{k+1} = \hat{y}_{k+1\vert k+1}$, ${ { {\hat{\pmb w}}}}_{k+1}={ { {\hat{\pmb w}}}}_{k+1\vert k+1}$;

步骤 2.4. 由式(17) 对$ \hat{y}_{k+1} $线性变换至$\tilde {y}_{k+1} $;

步骤 2.5. 对$\tilde {y}_{k+1} $做五点三次平滑;

步骤 2.6. 若$k=K-1$, 结束, 否则, $k=k+1$, 回到步骤2.1.

4 实验结果 4.1 系统设置

实验对记录的猕猴运动皮层神经元锋电位数及其对应手指移动位置的数据进行处理, 数据来源于文献[1]的公开数据1.具体过程为:猕猴手指在长为L, 宽为W的平面内运动, 当目标物出现, 猕猴的手指就指向目标物; 猕猴的脑部连接N个电极, 每隔$\Delta t$的采样时间就记录猕猴运动皮层的峰值电位数, 记录时间共为$K\Delta t$; 与此同时, 在每$\Delta t$时刻记录猕猴手指在移动平面内的坐标值以供验证.实验参数L=25 cm, $W=18$ cm; 采样周期$\Delta t=70$ ms; 数据长度$K=3\, 000$; 每一组神经元簇向量由$N=42$个电极记录不同的神经元的峰值电位数.

1http://booksite.elsevier.com/0780123838360

我们评价猕猴手指移动位置的估计性能采用以下方式:从数据中随机抽取一段进行有监督的训练或拟合, 再随机抽取一段进行验证, 评价的指标采用均方根误差(Rooted mean square error, RMSE) $e_r $, 定义为

$\begin{align} \label{eq1} e_r =\sqrt {\frac{1}{K}\sum\limits_{k=1}^K {(y_k - \hat{y}_k )^2} } \end{align}$ (19)

将本文的UCKD方法与传统的线性拟合方法[17-18]和SSM方法进行对比, 其中SSM方法的逐次状态估计器采用KF滤波[11-12], 涉及的参数设置:

1) 线性拟合方法

a) 训练:将数据分为两部分, 一半用于进行有监督的拟合, 另一半用于测试;

b) 拟合:采用最小二乘拟合得到式(1) 的线性函数, 并用该拟合函数估计位置.

2)  KF滤波方法

a) 训练:与线性拟合相同, 一半数据用于有监督的拟合, 另一半用于测试;

b) 拟合:采用最小二乘拟合得到${{\pmb f}}(y_k )$, 然后通过式(4) 和卡尔曼滤波方法估计$y_k$;

c) 初始值$ {y}_0 =0$;

d) 状态噪声$\omega _k $的方差

$Q_\omega =\frac{1}{K-1}\sum _{k=1}^K (y_{k+1} -y_k)^2$

其中, K为训练数据长度;

e) 观测噪声矢量${{\pmb v}}_k $的方差矩阵

${{\pmb R}}_v ={\rm diag}\left\{\sqrt {\frac{1}{K-1}({{\pmb s}}-{\rm { {\hat{\pmb s}}}})({{\pmb s}}-{\rm { {\hat{\pmb s}}}})^{\rm T}} \right\}$

其中, ${{\pmb S}}=[{{\pmb s}}_1, {{\pmb s}}_2, \cdots, {{\pmb s}}_K]$, ${ { {\hat{\pmb s}}}}_k ={{\pmb f}}(y_k )$.

3) 本文方法

a) 神经网络设置:输入层1, 输出层84, 激活函数$\varphi ({{\pmb s}}_k )={{\pmb s}}_k $;

b) 初始值

$\begin{align*} & \hat{y}_0 =0\\ &{{ {\hat{\pmb w}}}}_{0\vert 0} =[{\rm rand}(1, 84)-0.5]\left(\frac{84}{3}\right)^{\frac{1}{2}} \end{align*}$

其中, rand ($m, n)$为0~1之间均匀分布的$m\times n$矩阵;

c) 各方差矩阵初始值

$\begin{align*} &P_{y, 0\vert 0} =1\\ & {{ P}}_{w, 0\vert 0} =0.05\times{{ I}}_{84\times 84} \end{align*}$

d) 状态噪声$\omega _k $的方差$Q_\omega =0.8$;

e) 观测噪声矢量${{\pmb v}}_k $的方差矩阵${ { R}}_v $由式(10) 得到;

f) 遗忘因子$\lambda =0.005$;

g) 权值噪声${{\pmb \theta }}_k $的方差矩阵${ { R}}_w =10^{-2} ~\times~ { I}_{84\times 84}$.

需要注意的是, 本文方法中初始值${{\pmb {w}}}_{0\vert 0} $, 状态噪声方差, 权值噪声方差矩阵以及遗忘因子等参数的选择并不唯一, 其取值范围参见第4.2节, 在该取值范围内取值对算法的性能没有本质影响.

4.2 结果分析

图 6是线性拟合方法对猕猴手指移动的纵坐标估计曲线, 随机抽取350个数据训练进行线性拟合(训练数据长度采取350原因在于, 即使长度增加也不能显著改善算法性能), 然后对剩余数据进行测试.图 6给出了其中100个测试数据的结果图.图中Linear表示线性拟合, Linerfitsmooth表示线性拟合后再进行五点三次平滑, Actual表示手指实际运动的纵坐标.从图 6可以看出, 线性拟合方法能够将猕猴手指移动位置的变化趋势估计出来, 例如, 当约第10~20数据点间实际轨迹向下移动, 估计曲线也向下移动; 当约第30~40点间实际轨迹向上移动, 估计曲线也向上移动.但是该曲线即使经过平滑后在相邻时刻间的波动也比较大; 当约第50~60点间估计曲线有上下波动, 主要原因还是线性拟合把相邻时刻猕猴手指移动的过程看作独立所致.

图 6 线性拟合位置估计与手指移动实际位置曲线对比 Figure 6 Linear fitting position estimation compared with fingers moving curve of actual position

图 7给出了KF方法对猕猴手指移动纵坐标的估计, 其中随机抽取350个数据进行拟合得到线性函数, 然后再对剩余数据进行了估计, 图中仍然显示了其中100个点的数据.由于该方法采用了SSM模型, 将手指的移动当作非独立的过程, 因此不仅能估计手指移动的趋势, 而且曲线也较线性拟合光滑.例如, 约第50~60点间KF估计曲线的上下波动程度就比线性拟合的缓慢.

图 7 KF位置估计与手指移动实际位置曲线对比 Figure 7 KF algorithm position estimation compared with fingers moving curve of actual position

图 8给出了UCKD算法对猕猴手指移动纵坐标的估计, 由于该方法不需要有监督的数据, 因此只是随机抽取了一段长350的数据进行训练, 在估计过程中对该数据训练了20次, 然后再对该数据进行估计.图 8显示了数据中第101~200个数据的估计结果.从图 8可以看出, 虽然UCKD并没有采用有监督的数据, 但是也估计出了手指移动的趋势, 且UCKD方法还采用了五点三次平滑方法, 其曲线也比较平滑.另外, 注意到在约第50~60点间的估计曲线有较大误差.原因在于, 无论本文的无监督方法还是有监督的KF方法都是一种逐次状态估计, 即需要通过前一次的估计状态和当前的观测值去修正当前的估计状态, 它带来的效应就是当状态的趋势出现反转(由下降突变成上升或上升变下降) 时, 估计的曲线就会出现一定的误差或震荡, 20~30点正是趋势出现逐步改变时.图 7中有监督的KF方法在这些点上也有较大的震荡, 而且第10点和第80点附近的误差都很大, 而这些点恰恰是趋势发生改变的时候.

图 8 UCKD位置估计与手指移动实际位置曲线对比 Figure 8 UCKD algorithm position estimation compared with fingers moving curve of actual position

图 9给出了UCKD算法随训练次数增加RMSE收敛的过程图.图中用长度为3 000的数据进行了训练, 训练结束后再求这3 000个估计点的RMSE值.从图 9可以看出, 在第2次训练之后RMSE就稳定在2.1左右. 图 10是在一次训练下, UCKD随训练数据长度从1 200变化到3 000时, 对整个3 000点取得估计的RMSE曲线.从图 10可以看出, RMSE由1 200时的4.5 cm降到了3 000的2.2 cm左右, 说明随着数据长度增加, 估计误差逐渐减小.从图 9图 10可以发现, UCKD算法性能与训练次数和数据的长度有关系, 训练次数或数据长度的增加均会使估计值越来越精确. 表 1给出了5组交叉验证的RMSE, 首先随机抽取长度为500的数据进行UCKD训练, 然后再随机抽取另外500个数据进行测试.由表 1可以看出, RMSE随着训练次数的增加而减小, 并在第6次训练时趋向稳定.

图 9 UCKD训练次数RMSE变化图 Figure 9 UCKD's RMSE correlated with the number of training
图 10 UCKD算法随数据长度变化的位置估计RMSE曲线图 Figure 10 UCKD algorithm with data length change of position estimation RMSE curve
表 1 5次训练的RMSE (cm) Table 1 Five times the RMSE of training (cm)

为进一步评判上述三种方法的性能, 图 11给出了它们的RMSE曲线, 其中线性拟合方法和KF方法分别抽取了长度为50~350的有监督数据进行函数拟合, 然后再通过测试数据进行估计, 反复进行15次实验, 取平均值后得到最后RMSE曲线. UCKD不需要有监督数据, 对随机抽取的长度为350的数据分别进行10次和20次训练得到权值后, 再得到该数据估计的平均RMSE指标, 实验次数也是15次.从图 11可以看出, 线性拟合方法和KF方法的RMSE和训练数据的长度有关, 随着训练数据的长度由50增加至150, 线性拟合和KF方法分别从7 cm减小至2.2 cm和从5 cm减小至2 cm左右.这说明两种方法的估计性能依赖于有监督训练数据的长度, 当数据越多时, 拟合的函数越精确, 估计也越精确.对于UCKD来说, 训练次数为10次和20次的RMSE分别约为2.2 cm和1.9 cm左右, RMSE水平约等于线性拟合和KF方法采用有监督训练数据长度为150个时取得的水平.说明UCKD虽然没经过有监督数据的训练, 但是一样能取得线性拟合和KF方法取得的估计性能.

图 9~11可以看出, KF算法采用350个有监督数据进行拟合, 其RMSE曲线就收敛到约2 cm; 而UCKD需要对350个无监督数据训练10次, 或者对3 000个数据训练一次后才得到相似的RMSE性能.从时间上看, UCKD算法收敛时间约为KF或线性拟合算法的$3\, 000/350 =8.5$倍.因此, UCKD虽然无需有监督的数据训练, 但却以更多的无监督训练数据或训练次数为代价.

图 11 三种方法的RMSE变化图 Figure 11 Three algorithm of RMSE variation
5 讨论

通过猕猴运动皮层神经元脉冲信号来估计其手指移动位置是一种典型的神经解码问题, 本文针对该问题提出了一种无监督积分卡尔曼滤波解码方法, 称为UCKD方法.该解码方法采用神经网络得到神经元锋电位数与手指移动位置的关系权值, 然后再利用权值关系估计手指移动的位置.为此, 构建了两组SSM方程, 利用逐次状态估计方法来得到权值和移动位置.由于该SSM方程并不是线性, 为减少训练的复杂度和提高估计准确度, 逐次状态估计方法采用的是CKF滤波.同时, 由于方程并不具有唯一解, 往往与期望值成线性关系, 还需要将该结果线性变换至手指实际移动的变化区间内, 最后, 为了提高估计精度, 采用了五点三次平滑来减少RMSE值.该方法的最大特点是无监督训练, 在只知道神经元锋电位信号的情况下依旧能够估计猕猴手指移动的位置.相反, 现存的线性拟合以及SSM方法需要知道神经元锋电位信号以及与其对应的手指移动位置信息, 通过这些有监督的训练数据才能估计猕猴手指移动位置.从实验结果看, 当采用较少的有监督数据, 现存方法与本文方法相比有较大的估计误差; 当采用较多的有监督数据, 现存方法才具有与本文的方法相近似的估计误差.另外, 本文的UCKD方法是一种可同时进行训练和位置估计的方法, 随着神经元锋电位信号的长度增大, 其估计准确度会逐渐提高, 因此无需对同样一段数据进行反复训练而使计算时间增加.另外, 本算法需要一些先验知识.例如, 神经元锋电位信号与猕猴手指移动之间成何关系.我们假设其为线性关系[2], 然后针对该假设设计单层的神经网络进行训练.另外, 由于本文算法采用无监督学习方式, 相比于有监督算法, 其需要较多的训练次数或数据长度才能使算法收敛, 达到较低的估计误差.

参考文献
1 Wallisch P, Lusignan M E, Benayoun M D, Baker T I, Dickey A S, Hatsopoulos N G. MATLAB For Neuroscientists (2nd edition). London: Elsevier, 2014.
2 Gabbiani F, Cox S J. Mathematics for Neuroscientists. London: Elsevier, 2010.
3 Kass R E, Eden U T, Brown E N. Analysis of Neural Data. New York: Springer, 2014.
4 Hatsopoulos N G, Ojakangas C L, Paninski L, Donoghue J P. Information about movement direction obtained from synchronous activity of motor cortical neurons. Proceedings of the National Academy of Sciences of the United States of America, 1998, 95 (26): 15706–15711. DOI:10.1073/pnas.95.26.15706
5 Hatsopoulos N G, Xu Q Q, Amit Y. Encoding of movement fragments in the motor cortex. The Journal of Neuroscience, 2007, 27 (19): 5105–5114. DOI:10.1523/JNEUROSCI.3570-06.2007
6 Moran D W, Schwartz A B. Motor cortical representation of speed and direction during reaching. Journal of Neurophysiology, 1999, 82 (5): 2676–2692.
7 Reynaud-Bouret P, Rivoirard V, Grammont F, Tuleau-Malot C. Goodness-of-fit tests and nonparametric adaptive estimation for spike train analysis. Journal of Mathematical Neuroscience, 2014, 4 (1): 1–41. DOI:10.1186/2190-8567-4-1
8 Georgopoulos A P, Kettner R E, Schwartz A B. Primate motor cortex and free arm movements to visual targets in three-dimensional space. Ⅱ. Coding of the direction of movement by a neuronal population. The Journal of Neuroscience, 1988, 8 (8): 2928–2937.
9 Brockwell A E, Rojas A L, Kass R E. Recursive Bayesian decoding of motor cortical signals by particle filtering. Journal of Neurophysiology, 2004, 91 (4): 1899–1907. DOI:10.1152/jn.00438.2003
10 Brockwell A E, Kass R E, Schwartz A B. Statistical signal processing and the motor cortex. Proceedings of the IEEE, 2007, 95 (5): 881–898. DOI:10.1109/JPROC.2007.894703
11 Wu W, Shaikhouni A, Donoghue J P, Black M J. Closed-loop neural control of cursor motion using a Kalman filter. In:Proceedings of the 26th Annual International Conference of the IEEE Engineering in Medicine and Biology Society. San Francisco, CA:IEEE, 2004. 4126-4129
12 Wu W, Gao Y, Bienenstock E, Donoghue, J P, Black M J. Bayesian population decoding of motor cortical activity using a Kalman filter. Neural Computation, 2006, 18 (1): 80–118. DOI:10.1162/089976606774841585
13 Ebitz R B, Platt M L. Neuronal activity in primate dorsal anterior cingulate cortex signals task conflict and predicts adjustments in pupil-linked arousal. Neuron, 2015, 85 (3): 628–640. DOI:10.1016/j.neuron.2014.12.053
14 Wang W, Sudre G P, Xu Y, Kass R E, Collinger J L, Degenhart A D, Bagic A I, Weber D J. Decoding and cortical source localization for intended movement direction with MEG. Journal of Neurophysiology, 2010, 104 (5): 2451–2461. DOI:10.1152/jn.00239.2010
15 Campbell J, Sharma A. Visual cross-modal re-organization in children with cochlear implants. PLoS One, 2016, 11 (1): e0147793. DOI:10.1371/journal.pone.0147793
16 Hochberg L R, Serruya M D, Friehs G M, Mukand J A, Saleh M, Caplan A H, Branner A, Chen D, Penn R D, Donoghue J P. Neuronal ensemble control of prosthetic devices by a human with tetraplegia. Nature, 2006, 442 (7099): 164–171. DOI:10.1038/nature04970
17 Serruya M D, Hatsopoulos N G, Paninski L, Fellows M R, Donoghue J P. Instant neural control of a movement signal. Nature, 2002, 416 (6877): 141–142. DOI:10.1038/416141a
18 Warland D K, Reinagel P, Meister M. Decoding visual information from a population of retinal ganglion cells. Journal of Neurophysiology, 1997, 78 (5): 2336–2350.
19 Czanner G, Eden U T, Wirth S, Yanike M, Suzuki W A, Brown E N. Analysis of between-trial and within-trial neural spiking dynamics. Journal of Neurophysiology, 2008, 99 (5): 2672–2693. DOI:10.1152/jn.00343.2007
20 Malik W Q, Hochberg L R, Donoghue J P, Brown E N. Modulation depth estimation and variable selection in state-space models for neural interfaces. IEEE Transactions on Biomedical Engineering, 2015, 62 (2): 570–581. DOI:10.1109/TBME.2014.2360393
21 Brown E N, Frank L M, Tang D D, Quirk M C, Wilson M A. A statistical paradigm for neural spike train decoding applied to position prediction from ensemble firing patterns of rat hippocampal place cells. The Journal of Neuroscience, 1998, 18 (18): 7411–7425.
22 Shanechi M M, Wornell G W, Williams Z M, Brown E N. Feedback-controlled parallel point process filter for estimation of goal-directed movements from neural signals. IEEE Transactions on Neural Systems and Rehabilitation Engineering, 2013, 21 (1): 129–140. DOI:10.1109/TNSRE.2012.2221743
23 Sanger T D. Neural population codes. Current Opinion in Neurobiology, 2003, 13 (2): 238–249. DOI:10.1016/S0959-4388(03)00034-5
24 Zhang K C, Sejnowski T J. Accuracy and learning in neuronal populations. Progress in Brain Research, 2001, 130 : 333–342. DOI:10.1016/S0079-6123(01)30022-5
25 Haykin S O. Neural Networks and Learning Machines (3rd edition). New Jersey:Pearson Education, Inc., 2009.
26 Arasaratnam I, Haykin S. Cubature Kalman filters. IEEE Transactions on Automatic Control, 2009, 54 (6): 1254–1269. DOI:10.1109/TAC.2009.2019800
27 Friedman N. The Bayesian structural EM algorithm. In:Proceedings of the 14th Conference on Uncertainty in Artificial Intelligence. San Francisco, CA, USA:Morgan Kaufmann Publishers Inc., 1998. 129-138
28 Chou Y L. Statistical Analysis. Holt International, 1975.