本发明涉及一种脑区脉冲神经信号的预测方法，包括如下步骤：1)、脑区脉冲神经信号通道预处理；2)、后继脑区脉冲神经信号生成概率模型建模；3)、模型准确性度量；4)、利用数值梯度下降进行模型优化；5)迭代计算完成后继脑区脉冲神经信号预测。本预测方法将脉冲神经信号的点过程的天然特性纳入预测模型的优化目标中，从而提高模型对神经脉冲序列的预测能力。
1.一种脑区脉冲神经信号的预测方法，其特征在于，包括如下步骤：
1)、脑区脉冲神经信号通道预处理：首先将原始脑区脉冲神经信号采用时间槽划分完成离散化；筛除所有脑区脉冲神经信号通道中脉冲发放率过高或过低的信号通道；所述筛除方法是指：首先将每一路脉冲神经信号通道按照时间轴划分成50％训练集、20％验证集和30％测试集，检测此信号通道的三个集合中脉冲发放率是否均在2Hz至40Hz之间，满足此条件便通过筛选；
2)、后继脑区脉冲神经信号生成概率模型建模：利用互信息获取与后继脑区的每一路信号通道互信息值最高的若干前继脑区信号通道作为模型自变量，然后针对每一个待预测的后继脑区信号通道中的每一个时间槽设置滑动窗口来选取模型的输入自变量；得到模型的输入自变量之后，以非齐次泊松点过程通用线性模型进行后继脑区脉冲神经信号生成概率模型建模；
3)、模型准确性度量：通过步骤2)中的模型得到后继脑区待预测脉冲神经信号的脉冲发放概率序列；利用离散时间重放缩柯尔莫诺夫-斯米尔诺夫统计量作为度量模型效果的量纲，计算模型的效果度量值；
4)、利用数值梯度下降进行模型优化：将步骤3)中得到的模型的效果度量值作为目标函数，然后对目标函数求导得到数值梯度，采用拟牛顿法进行最优参数优化计算；
5)迭代计算：在求解最优参数计算的过程中，迭代执行步骤3)和步骤4)，直至优化后的参数对应的模型准确性达到局部极值点或达到迭代次数上限，将最终确定的参数代入模型中，并利用此模型完成后继脑区脉冲神经信号预测。
2.根据权利要求1所述的脑区脉冲神经信号的预测方法，其特征在于，所述步骤1)中离散化是指：采取10毫秒作为时间槽宽度对原始脑区脉冲神经信号进行划分，将时间槽内存在脉冲神经信号的时间槽记为1，反之记为0，以此完成离散化。
3.根据权利要求1所述的脑区脉冲神经信号的预测方法，其特征在于，所述步骤2)中获取模型自变量的方法：计算前继脑区与后继脑区各路脉冲神经信号两两互信息值之后，针对后继脑区每一路脉冲神经信号，选取与之互信息值最高的前8路前继脑区信号通道作为此后继脑区信号通道的模型自变量；互信息值为：其中，xi为前继脑区第i路信号通路，yj为后继脑区第j路信号通路，p(xt,yt)是事件xt与yt同时发生的联合概率，p(xt)和p(yt)分别表示事件xt与yt各自发生的概率，由于原始脑区脉冲神经信号已经被离散化，因此xt和yt的取值范围均在{0，1}之间。
4.根据权利要求3所述的脑区脉冲神经信号的预测方法，其特征在于，所述步骤2)中选取模型的输入自变量的方法：针对每一个待预测的后继脑区信号通道中的每一个时间槽设置一个时间跨度为90毫秒的滑动窗口，用来截取与之互信息值最高的前8路前继脑区脉冲信号通道，作为模型的输入自变量，此滑动窗口的末端与当前时间槽保持齐平。
5.根据权利要求4所述的脑区脉冲神经信号的预测方法，其特征在于，所述步骤2)中以非齐次泊松点过程通用线性模型进行后继脑区脉冲神经信号生成概率模型建模，其具体形式为：其中，n(t)表示第t个时间槽内发放脉冲的个数，由于原始脑区脉冲神经信号采取离散化处理，因此n(t)取值范围为{0,1}；P(n(t))表示在第t个时间槽内发放n(t)个脉冲的概率；λ(t|x(t))表示第t个时间槽内发放脉冲的条件概率密度函数；x(t)是由滑动窗口采集得到的前继脑区脉冲神经信号，Δ是时间槽的长度；
随后将λ(t)与自变量x(t)联系在一起，采用以指数函数为连接函数的通用线性模型，具体形式为：其中，λ(t)＝λ(t|x(t))，θ是待确定的模型参数向量。
6.根据权利要求5所述的脑区脉冲神经信号的预测方法，其特征在于，所述步骤3)中利用离散时间重放缩柯尔莫诺夫-斯米尔诺夫统计量作为度量模型效果的量纲，计算模型的效果度量值，其具体形式是：I、首先为了消除离散时间采样带来的影响，在λ(t)的基础之上生成一个新序列q(t)，其具体形式为：q(t)＝-log(1-λ(t)Δ) (4)
II、随后将序列q(t)按照待预测神经元产生的脉冲序列进行切割，然后对相邻的两个脉冲信号之间的q(t)值进行积分，并对离散时间处理产生的偏差进行修正，具体形式为：其中，ti指待预测神经元的第i个脉冲信号发生的时间，ti+1指待预测神经元的第i+1个脉冲信号发生的时间，δ(i)是一个按照以下形式采集得到的随机变量：式(6)中的r(i)是一个在[0,1]之间均匀分布的随机变量；
根据柯尔莫诺夫-斯米尔诺夫检验理论，式(7)产生的集合{z(i)}的分布应遵循在[0,
IV、将集合{z(i)}按照数值由小到大进行重排列，并与标准[0,1]之间的均匀分布同时绘图于坐标系中进行比较，横轴为集合{z(i)}，纵轴为标准均匀分布，取纵轴上二者最大数值差作为效果度量值f(y|θ)。
7.根据权利要求6所述的脑区脉冲神经信号的预测方法，其特征在于，所述步骤3)中计算模型的效果度量值时，针对每一个待预测神经元信号通路，计算过程均执行10遍并取平均值。
8.根据权利要求6所述的脑区脉冲神经信号的预测方法，其特征在于，所述步骤4)中模型优化是指：通过对目标函数f(y|θ)进行近似求导直接计算目标函数f(y|θ)的数值梯度，求解数值梯度的具体形式为：其中，y是后继脑区待预测脉冲神经信号通道，θ是待确定的模型参数向量，f(y|θ)是目标函数，θa与θb在第i维上分别表示为 以及 其余维度的值与θ保持一致，此处ε的取值为1e-5。
一种脑区脉冲神经信号的预测方法
技术领域
[0001] 本发明涉及脉冲神经信号预测领域，具体涉及一种脑区脉冲神经信号的预测方法。
背景技术
[0002] 随着微电极阵列技术的不断发展和现代医学对脑功能研究的不断深入，大脑神经元脉冲信号能够以极高的时空分辨率被采集和分析，这为科研人员们探索大脑信号传递机制和各个脑区的功能联系提供了便利。
[0003] 大脑可以依照不同的功能分区粗略划分为若干脑区。利用信号前继脑区的脑信号来预测信号后继脑区的脑信号，是了解两个在物理分布上相连的脑区之间功能联系的一大途径。由于大脑神经元以脉冲信号为基本单位进行神经元间通信的工作原理，单个神经元所产生的脑信号本身可视为一个时间序列点过程。然而，现有的多数模型大多未将脑信号的点过程特性纳入模型的考量中，这对准确预测目标脑区的脉冲信号、判断相邻脑区之间的功能联系带来一定的挑战。
发明内容
[0004] 本发明的目的在于针对现有技术的不足，提供一种脑区脉冲神经信号的预测方法，以通用线性模型为基础、以离散时间重放缩柯尔莫诺夫-斯米尔诺夫统计量为优化目标、以数值梯度为优化手段的脉冲神经信号预测方法，主要解决的问题是将脉冲神经信号的点过程的天然特性纳入预测模型的优化目标中，从而提高模型对神经脉冲序列的预测能力。
[0005] 本发明解决上述技术问题所提供的技术方案为：
[0006] 一种脑区脉冲神经信号的预测方法，包括如下步骤：
[0007] 1)、脑区脉冲神经信号通道预处理：首先将原始脑区脉冲神经信号采用时间槽划分完成离散化；筛除所有脑区脉冲神经信号通道中脉冲发放率过高或过低的信号通道；
[0008] 2)、后继脑区脉冲神经信号生成概率模型建模：利用互信息获取与后继脑区的每一路信号通道互信息值最高的若干前继脑区信号通道作为模型自变量，然后针对每一个待预测的后继脑区信号通道中的每一个时间槽设置滑动窗口来选取模型的输入自变量；得到模型的输入自变量之后，以非齐次泊松点过程通用线性模型进行后继脑区脉冲神经信号生成概率模型建模；
[0009] 3)、模型准确性度量：通过步骤2)中的模型得到后继脑区待预测脉冲神经信号的脉冲发放概率序列；利用离散时间重放缩柯尔莫诺夫-斯米尔诺夫统计量作为度量模型效果的量纲，计算模型的效果度量值；
[0010] 4)、利用数值梯度下降进行模型优化：将步骤3)中得到的模型的效果度量值作为目标函数，然后对目标函数求导得到数值梯度，采用拟牛顿法进行最优参数优化计算；
[0011] 5)迭代计算：在求解最优参数计算的过程中，迭代执行步骤3)和步骤4)，直至优化后的参数对应的模型准确性达到局部极值点或达到迭代次数上限，将最终确定的参数代入模型中，并利用此模型完成后继脑区脉冲神经信号预测。
[0012] 所述步骤1)中离散化是指：采取10毫秒作为时间槽宽度对原始脑区脉冲神经信号进行划分，将时间槽内存在脉冲神经信号的时间槽记为1，反之记为0，以此完成离散化。
[0013] 所述步骤1)中筛除方法是指：首先将每一路脉冲神经信号通道按照时间轴划分成50％训练集、20％验证集和30％测试集，检测此信号通道的三个集合中脉冲发放率是否均在2Hz至40Hz之间，满足此条件便通过筛选。
[0014] 所述步骤2)中获取模型自变量的方法：计算前继脑区与后继脑区各路脉冲神经信号两两互信息值之后，针对后继脑区每一路脉冲神经信号，选取与之互信息值最高的前8路前继脑区信号通道作为此后继脑区信号通道的模型自变量。互信息值为：
[0015]
[0016] 其中，xi为前继脑区第i路信号通路，yj为后继脑区第j路信号通路，p(xt,yt)是事件xt与yt同时发生的联合概率，p(xt)和p(yt)分别表示事件xt与yt各自发生的概率，由于原始脑区脉冲神经信号已经被离散化，因此xt和yt的取值范围在{0，1}之间。
[0017] 所述步骤2)中选取模型的输入自变量的方法：针对每一个待预测的后继脑区信号通道中的每一个时间槽设置一个时间跨度为90毫秒的滑动窗口，用来截取与之互信息值最高的前8路前继脑区脉冲信号通道，作为模型的输入自变量，此滑动窗口的末端与当前时间槽保持齐平。
[0018] 所述步骤2)中以非齐次泊松点过程通用线性模型进行后继脑区脉冲神经信号生成概率模型建模，其具体形式为：
[0019]
[0020] 其中，n(t)表示第t个时间槽内发放脉冲的个数，由于原始脑区脉冲神经信号采取离散化处理，因此n(t)取值范围为{0,1}；P(n(t))表示在第t个时间槽内发放n(t)个脉冲的概率；λ(t|x(t))表示第t个时间槽内发放脉冲的条件概率密度函数；x(t)是由滑动窗口采集得到的前继脑区脉冲神经信号，Δ是时间槽的长度；
[0021] 随后将λ(t)与自变量x(t)联系在一起，采用以指数函数为连接函数的通用线性模型，具体形式为：
[0022]
[0023] 其中，λ(t)＝λ(t|x(t))。
[0024] 所述步骤3)中利用离散时间重放缩柯尔莫诺夫-斯米尔诺夫统计量作为度量模型效果的量纲，计算模型的效果度量值，其具体形式是：
[0025] I、首先为了消除离散时间采样带来的影响，在λ(t)的基础之上生成一个新序列q(t)，其具体形式为：
[0026] q(t)＝-log(1-λ(t)Δ) (4)
[0027] II、随后将序列q(t)按照待预测神经元产生的脉冲序列进行切割，然后对相邻的两个脉冲信号之间的q(t)值进行积分，并对离散时间处理产生的偏差进行修正，具体形式为：
[0028]
[0029] 其中，ti指待预测神经元的第i个脉冲信号发生的时间，ti+1指待预测神经元的第i+1个脉冲信号发生的时间，δ(i)是一个按照以下形式采集得到的随机变量：
[0030]
[0031] 式(6)中的r(i)是一个在[0,1]之间均匀分布的随机变量；
[0032] III、随后对ζ(i)进行尺度放缩，得到：
[0033] z(i)＝1-e-ζ(i) (7)
[0034] 根据柯尔莫诺夫-斯米尔诺夫检验理论，式(7)产生的集合{z(i)}的分布应遵循在[0,1]之间的均匀分布；
[0035] IV、将集合{z(i)}按照数值由小到大进行重排列，并与标准[0,1]之间的均匀分布同时绘图于坐标系中进行比较，横轴为集合{z(i)}，纵轴为标准均匀分布，取纵轴上二者最大数值差作为效果度量值f(y|θ)。
[0036] 所述步骤3)中计算模型的效果度量值时，针对每一个待预测神经元信号通路，计算过程均执行10遍并取平均值。
[0037] 所述步骤4)中模型优化是指：通过对目标函数f(y|θ)进行近似求导直接计算目标函数f(y|θ)的数值梯度，求解数值梯度的具体形式为：
[0038]
[0039] 其中，y是后继脑区待预测脉冲神经信号通道，θ是待确定的模型参数向量，f(y|θ)是目标函数，θa与θb在第i维上分别表示为 以及 其余维度的值与θ保持一致，此处ε的取值为1e-5。
[0040] 同现有技术相比，本发明的有益效果体现在：将脉冲神经信号的点过程的天然特性纳入预测模型的优化目标中，从而提高模型对神经脉冲序列的预测能力。
附图说明
[0041] 图1是实施例中后继脑区脉冲神经信号预测的整体流程图；
[0042] 图2是实施例中后继脑区脉冲神经信号预测的各个变量之间的依赖图。
具体实施方式
[0043] 下面结合具体的实施例以及附图对本发明作进一步的说明。
[0044] 如图1和2所示，脑区脉冲神经信号的预测方法，具体步骤如下：
[0045] (1)脑区脉冲神经信号通道预处理
[0046] 首先使用微电极阵列采集原始脑区脉冲神经信号，然后将原始脑区脉冲神经信号采用时间槽划分完成离散化，具体方法：采取10毫秒作为时间槽宽度对原始脑区脉冲神经信号进行等距划分，将时间槽内存在脉冲神经信号的时间槽记为1，反之记为0，以此完成离散化。
[0047] 随后筛除所有脑区脉冲神经信号通道中脉冲发放率过高或过低的信号通道；具体方法：首先将每一路脉冲神经信号通道按照时间轴划分成50％训练集、20％验证集和30％测试集，检测此信号通道的三个集合中脉冲发放率是否均在2Hz至40Hz之间，满足此条件便通过筛选。
[0048] (2)后继脑区脉冲神经信号生成概率模型建模
[0049] 计算前继脑区与后继脑区各路脉冲神经信号两两互信息值之后，针对后继脑区每一路脉冲神经信号，选取与之互信息值最高的前8路前继脑区信号通道作为此后继脑区信号通道的模型自变量。互信息值为：
[0050]
[0051] 其中，xi为前继脑区第i路信号通路，yj为后继脑区第j路信号通路，p(xt,yt)是事件xt与yt同时发生的联合概率，p(xt)和p(yt)分别表示事件xt与yt各自发生的概率，由于原始脑区脉冲神经信号已经被离散化，因此xt和yt的取值范围在{0，1}之间。
[0052] 然后针对每一个待预测的后继脑区信号通道中的每一个时间槽设置一个时间跨度为90毫秒的滑动窗口，用来截取与之互信息值最高的前8路前继脑区脉冲信号通道，作为模型的输入自变量，此滑动窗口的末端与当前时间槽保持齐平。
[0053] 得到模型的输入自变量之后，以非齐次泊松点过程通用线性模型进行后继脑区脉冲神经信号生成概率模型建模，其具体形式为：
[0054]
[0055] 其中，n(t)表示第t个时间槽内发放脉冲的个数，由于原始脑区脉冲神经信号采取离散化处理，因此n(t)取值范围为{0,1}；P(n(t))表示在第t个时间槽内发放n(t)个脉冲的概率；λ(t|x(t))表示第t个时间槽内发放脉冲的条件概率密度函数(简写为λ(t))；x(t)是由滑动窗口采集得到的前继脑区脉冲神经信号，Δ是时间槽的长度；
[0056] 随后将λ(t)与自变量x(t)联系在一起，采用以指数函数为连接函数的通用线性模型，具体形式为：
[0057]
[0058] 其中，λ(t)＝λ(t|x(t))。
[0059] (3)模型准确性度量
[0060] 通过步骤2)中的模型得到后继脑区待预测脉冲神经信号的脉冲发放概率序列；利用离散时间重放缩柯尔莫诺夫-斯米尔诺夫统计量作为度量模型效果的量纲，计算模型的效果度量值，其具体形式是：
[0061] I、首先为了消除离散时间采样带来的影响，在λ(t)的基础之上生成一个新序列q(t)，其具体形式为：
[0062] q(t)＝-log(1-λ(t)Δ) (4)
[0063] II、随后将序列q(t)按照待预测神经元产生的脉冲序列进行切割，然后对相邻的两个脉冲信号之间的q(t)值进行积分，并对离散时间处理产生的偏差进行修正，具体形式为：
[0064]
[0065] 其中，ti指待预测神经元的第i个脉冲信号发生的时间，ti+1指待预测神经元的第i+1个脉冲信号发生的时间，δ(i)是一个按照以下形式采集得到的随机变量：
[0066]
[0067] 式(6)中的r(i)是一个在[0,1]之间均匀分布的随机变量；
[0068] III、随后对ζ(i)进行尺度放缩，得到：
[0069] z(i)＝1-e-ζ(i) (7)
[0070] 根据柯尔莫诺夫-斯米尔诺夫检验理论，式(7)产生的集合{z(i)}的分布应遵循在[0,1]之间的均匀分布；
[0071] IV、将集合{z(i)}按照数值由小到大进行重排列，并与标准[0,1]之间的均匀分布同时绘图于坐标系中进行比较，横轴为集合{z(i)}，纵轴为标准均匀分布，取纵轴上二者最大数值差作为效果度量值f(y|θ)。效果度量值f(y|θ)随着模型的参数θ的变化而变化，因此可看作是关于θ的一个函数，也即模型的目标函数(此处y是待预测后继脑区脉冲神经信号通道)。由于在计算过程中引入了随机变量，因此对此效果度量值的计算过程需要进行若干遍并取平均值来消除影响。针对每一个待预测神经元信号通路，计算过程均执行10遍并取平均值。
[0072] (4)利用数值梯度下降进行模型优化
[0073] 将步骤3)中得到的模型的效果度量值作f(y|θ)为目标函数，然后对目标函数求导得到数值梯度，采用拟牛顿法进行最优参数优化计算；
[0074] 通过对目标函数f(y|θ)进行近似求导直接计算目标函数f(y|θ)的数值梯度，求解数值梯度的具体形式为：
[0075]
[0076] 其中，y是后继脑区待预测脉冲神经信号通道，θ是待确定的模型参数向量，f(y|θ)是目标函数，θa与θb在第i维上分别表示为 以及 其余维度的值与θ保持一致，此处ε的取值为1e-5。
[0077] 利用数值梯度近似求得目标函数的导数值之后，采用拟牛顿法进行最优化计算。由于数值梯度计算过程较为耗时，利用拟牛顿法进行优化相比于最速下降法而言，能够更快收敛至局部极值，从而尽量减少迭代次数，以达成减少计算量的目的。
[0078] (5)迭代计算
[0079] 在求解最优参数计算的过程中，迭代执行步骤3)和步骤4)，直至优化后的参数对应的模型准确性达到局部极值点或达到迭代次数上限。达到局部极值点的标准是模型的目标函数值不再随着新一轮迭代而降低，预设的迭代次数上限是400次，若已达到迭代次数上限，则直接退出参数优化过程，不再进行新一轮参数迭代。将最终确定的参数代入模型中，并利用此模型完成后继脑区脉冲神经信号预测。
[0080] 效果比较
[0081] 本发明的试验结果在浙江大学求是高等研究院所采集的猕猴(编号为B04)的PMd和M1脑区的神经元脉冲信号上进行了实验，该信号的背景是猕猴执行四方向center-out任务时采集而得，信号有效时长为636.4秒，包括127路PMd脑区和102路M1脑区的神经元脉冲信号通路，经过筛选后得到96路有效PMd信号通路和57路有效M1信号通路。
[0082] 本发明的实验效果同时与其他几种脉冲信号模型得到的效果进行比较，包括线性回归模型(LR)、Spike-Triggered Average(STA)、Spike-Triggered Convaiance(STC)、传统通用线性模型(GLM)等。表1是本方法与其余几类方法对于脉冲信号预测效果的横向比较，其中DTR-KS是离散时间重放缩柯尔莫诺夫-斯米尔诺夫检验(Discrete Time Rescaling Kolmogorov Smirnov test)的简称，其值在[0,1]之间，数值越接近0表示该模型的预测效果越好。
[0083] 表1本方法与其他几类预测方法的效果比较
[0084]
