基于非负块稀疏贝叶斯学习的瞬态脑电源定位方法及系统与流程

1.本发明涉及脑电图(eeg)溯源技术领域，特别涉及一种基于非负块稀疏贝叶斯学习(nnbsbl)的瞬态脑电源定位方法及系统。


背景技术：

2.脑电图(eeg)是伴随大脑活动产生的可被记录的电位信号，对于研究大脑活动、评估大脑功能、临床以及心理学研究、脑机接口等具有十分重要的意义。eeg包含瞬态eeg与稳态eeg两种形式，其中瞬态eeg体现为持续时间较短的具有一定特征的eeg波形，而稳态eeg则在较长时间段内与特定频率相关联。瞬态eeg可由诱发产生，比如各种形式的瞬态诱发电位，也可以自发产生，比如癫痫发作时的棘(尖)慢复合波等。
3.eeg的源定位问题，也被称为eeg溯源问题，是通过利用仅在头皮上放置的多通道电极所采集的eeg信号，设计算法逆推出产生这些可被记录的eeg信号的电流偶极子源在大脑皮层上的位置。eeg源定位方法，目前大体可分为基于偶极子匹配的方法和基于分布式源模型的方法两类。其中基于偶极子匹配的方法需要预设偶极子源的个数，然而实际中源的个数无法被预知或者准确地估计出来，故目前eeg源定位的主流研究都是在分布式源模型下开展的。
4.分布式源模型利用一个预设的离散网格将大脑皮层完全覆盖，网格的每个顶点就代表了eeg源在皮层上的可能的位置，并且分布式源模型利用有限元法或者边界元法建立了对应的导程场矩阵以表征由皮层网格各顶点到头皮上的各记录电极上的映射关系。在分布式源模型下，eeg源定位是一个严重的病态问题，主要解决思路为空间滤波器方法和断层成像方法。空间滤波器法基于抑制不在扫描范围内的所有位置上的信号的原理，设计空间滤波器扫描可能存在偶极子源的每个顶点位置并估计该位置上的强度，但该方法不能很好地考虑偶极子源之间存在相关性的情况，且定位结果的分辨率较低。相对地，断层成像方法正成为近些年来eeg源定位中最受关注的一类方法，这类方法通过设计正则化约束来试图解决该严重的欠定问题，以标准化低分辨率脑电磁成像(sloreta)为典型方法。近些年来，贝叶斯推断思想被引入该领域以解决eeg源定位问题，而各种断层成像方法被证明可视为使用了不同的源分布先验信息(对应于正则化约束)的贝叶斯推断方法。目前绝大部分eeg源定位贝叶斯推断方法都是在假设偶极子源在皮层上稀疏分布这一特性来设计的，即导出了稀疏贝叶斯学习(sbl)方法，其中具有代表性的是champagne算法(david p.wipf,et al.robust bayesian estimation of the location,orientation,and time course ofmultiple correlated neural sources using meg,neuroimage,2010)。这些方法大多数假设产生eeg的电流偶极子源只存在于皮层网格的少数顶点上，即满足元稀疏(atom-wise sparsity)，然而在实际中，eeg信号应源于神经功能核团，每个核团包含网格上相邻的多个顶点，且核团内部每个顶点上的电流偶极子源信号之间应具有相关性。在这个考虑下，电流偶极子源在皮层网格上的元稀疏特征应被块稀疏(block-wise sparsity)所替代。目前，在champagne算法的基础上，结合块稀疏先验，已发展出了tree-champagne算法
(chang cai,et al.hierarchical multiscale bayesian algorithm for robust meg/eeg source reconstruction,neuroimage,2018)和smooth-champagne算法(chang cai,et al.robust empirical bayesian reconstruction ofdistributed sources for electromagnetic brain imaging,ieee tmi,2020)等。这些结合块稀疏先验的champagne算法面临着需要预设或者预知噪声功率的问题，为此研究者们又提出了先估计噪声功率后估计信号协方差矩阵的两步块稀疏贝叶斯学习(bsbl-2s)算法(alejandro ojeda,fast and robust block-sparse bayesian learning for eeg source imaging,neuroimage,2018)。在eeg源定位问题建模过程中，噪声往往不仅包括信号获取过程中线路的电噪声，还包括与感兴趣的eeg成分无关的自发eeg成分，故很难保证信号采集过程中每个电极通道上的噪声功率都相等，即噪声很可能是空间非均匀分布的；此外，在一个迭代算法中同时估计噪声功率与信号的协方差矩阵，两者在估计过程中会相互之间引入误差。因此，设计一种结合块稀疏先验的瞬态eeg源定位sbl算法，在不需要预知或预设各电极通道上的噪声功率的条件下，只在算法中考虑对瞬态eeg信号的偶极子源位置进行估计，并且在算法中通过某种方式免除对噪声功率的估计环节，预期可提升瞬态eeg源定位的算法性能。
5.目前在本领域有一些已公开的eeg源定位方法的专利申请，如：
6.发明专利申请(基于稀疏性约束的脑电源定位方法，申请号200610021584.6)将lp范数正则化约束融合到focuss算法中，从而给出了一种eeg源定位方法，该方法未考虑电流偶极子源在皮层上的块稀疏特征，且lp-focuss算法的分辨率和估计精度较低。
7.发明专利申请(一种基于格兰杰因果性的脑电源定位方法，申请号201410636555.5)通过分析每个电极通道与其他电极通道间的格兰杰因果性进行源定位，以期解决脑电逆问题解的不稳定性和非唯一性，然而该方法所确定的源的位置在头皮电极处的下方皮层，未引入导程场矩阵的映射作用，因此源定位精度将受到影响。
8.发明专利申请(一种基于lstm的脑电信号源定位方法，申请号201910178711.0)通过训练出一个lstm神经网络，将多通道eeg数据输入该神经网络来直接求得源的位置，然而该方法无论在神经网络的训练还是测试环节都未考虑电流偶极子源在皮层上的块稀疏分布，且输出二值结果非常容易导致eeg源定位结果的高假阳性或假阴性。
9.发明专利申请(基于稀疏贝叶斯学习的稳态诱发响应脑源定位方法，申请号202010001431.5)给出了一种针对稳态eeg的sbl源定位方法，然而该方法针对的是稳态eeg而非瞬态eeg，且仍未考虑电流偶极子源在皮层上的块稀疏分布。
10.发明专利申请(基于种群多样性控制的脑电源定位方法、系统及设备，申请号202011096444.1)，提出一种基于种群多样性控制的量子粒子群算法，通过控制算法的种群多样性，提高算法全局遍历寻优能力，该方法仍属于偶极子匹配方法的范畴，需要预设电流偶极子源的个数，限制了其实用性。
11.综合当前研究现状，断层成像方法是近年来eeg源定位中最受关注的一类方法，其中的sbl方法应用于瞬态eeg源定位问题具有很好的优势，但结合实际中电流偶极子源在皮层上的块稀疏分布先验特征，在源定位过程中无需预知、预设或者估计噪声信息，进而采取措施提升源定位估计精度，尚需提出切实可行的方法。
12.综上所述，传统技术大多存在以下技术问题：
13.1.未充分考虑电流偶极子源在皮层上的块稀疏分布先验特征。
14.2.需要预设噪声功率。
15.3.需要预设电流偶极子源的个数。
16.4.不能解决空间非均匀噪声问题，且在算法中同时估计噪声功率会影响电流偶极子源的定位精度。
17.5.定位目标往往是各顶点上的时变信号本身而非信号的功率，未利用信号功率的非负性来提升eeg源定位的效果。


技术实现要素：

18.本发明要解决的技术问题是提供一种在瞬态脑电源定位算法的建模阶段充分考虑电流偶极子源在皮层上的块稀疏分布先验特征；在无需预设或预知噪声功率的前提下进行瞬态脑电源定位；在无需预设电流偶极子源的个数的前提下进行瞬态脑电源定位；在瞬态脑电源定位算法中绕开对噪声功率的估计，从而提升源定位效果，并且同时适用于空间非均匀噪声场景；以及在瞬态脑电源定位算法中充分利用信号功率的非负性来提升源定位效果的基于非负块稀疏贝叶斯学习的瞬态脑电源定位方法。
19.为了解决上述问题，本发明提供了一种基于非负块稀疏贝叶斯学习的瞬态脑电源定位方法，其包括以下步骤：
20.s1、利用m个电极采集对应于瞬态响应的多通道eeg，采样率为fs，共记录t个采样点，并对其预处理获得m
×
t维预处理后的多通道eeg数据y；同时，将皮层划分成n
roi
个脑区，根据这些脑区将导程场矩l阵划分成各脑区上的子导程场矩阵ln,n＝1,2,...,n
roi
，并计算各脑区上的空间平滑核矩阵kn,n＝1,2,...,n
roi
；
21.s2、利用预处理后的多通道eeg数据y计算样本协方差矩阵利用样本协方差矩阵各脑区上的子导程场矩阵ln和各脑区上的空间平滑核矩阵kn,n＝1,2,...,n
roi
，构造加权样本数据矩阵z以及平滑脑区基矩阵m＝1,2,...,m；
22.s3、将非负脑区功率稀疏支撑向量α初始化为α
init
，并设置迭代参数；
23.s4、利用加权样本数据矩阵z、平滑脑区基矩阵m＝1,2,...,m以及前一非负脑区功率稀疏支撑向量α计算非负脑区功率系数γ的后验均值向量μ
γ
和后验协方差矩阵σ
γ
；
24.s5、根据非负脑区功率系数γ的后验均值向量μ
γ
和后验协方差矩阵σ
γ
计算非负脑区功率稀疏支撑向量α的更新值α
new
；
25.s6、判断是否满足迭代停止条件：迭代次数n≥nm或||α
new-α||2/||α||2≤ε，其中||
·
||2表示l2范数，若不满足则令α＝α
new
；若是，则回到步骤s4继续进行迭代运算；否则结束迭代，利用加权样本数据矩阵z、平滑脑区基矩阵m＝1,2,...,m以及当前非负脑区功率稀疏支撑向量α的更新值α
new
计算非负脑区功率系数γ的一阶矩《γ》，并联合各脑区上的空间平滑核矩阵kn,n＝1,2,...,n
roi
计算得到皮层神经元电活动的空间协方差矩阵σ
x
(《γ》)。
26.作为本发明的进一步改进，所述计算各脑区上的空间平滑核矩阵kn,n＝1,2,...,n
roi
包括：
27.第n个脑区上的空间平滑核矩阵的计算公式为n＝1,2,...,n
roi
，其中，hn表示从全脑离散拉普拉斯平滑矩阵h中抽取的对应于第n个脑区的子阵，(
·
)
t
表示矩阵转置，(
·
)-1
表示矩阵求逆。
28.作为本发明的进一步改进，所述全脑离散拉普拉斯平滑矩阵h的计算如下：
29.对于导程场矩阵l的列所对应的nv个大脑皮层网格顶点，利用这些顶点的位置和顶点两两之间的连接情况，所得nv×nv
维全脑离散拉普拉斯平滑矩阵h的计算方法为：其第(i,i)个元素当i≠j且第i个顶点和第j个顶点间有连接时h的第(i,j)个元素当i≠j且第i个顶点和第j个顶点间无连接时h的第(i,j)个元素h[i,j]＝0，其中表示与第i个顶点有连接的所有邻居顶点的距离的平均值，表示第i个顶点与跟它有连接的所有邻居顶点的距离的倒数的平均值，ni表示与第i个顶点有连接的所有邻居顶点的个数，d
ij
表示在第i个顶点和第j个顶点间有连接的情况下这两个顶点间的距离。
[0030]
作为本发明的进一步改进，步骤s2中，所述加权样本数据矩阵z＝[z1,z2,...,zm]，其第m列的计算公式为：
[0031][0032]
其中，第m列样本协方差向量的协方差矩阵其中(
·
)
t
表示矩阵转置，表示的第(m,m)个元素，表示的第m列，jm为第m个权重矩阵，其计算公式为j1＝[0
m-1
,i
m-1
]，jm＝[i
m-1
[:,1:m-1],0
m-1
,i
m-1
[:,m:m-1]],2≤m≤m-1，jm＝[i
m-1
,0
m-1
]，0
m-1
表示一个(m-1)
×
1维全1向量，i
m-1
表示一个(m-1)
×
(m-1)维单位矩阵；
[0033]
平滑脑区基矩阵m＝1,2,...,m的计算公式为：
[0034][0035]
其中，第n个空间平滑成分子阵n＝1,2,...,n
roi
，an[:,m]表示矩阵an的第m列。
[0036]
作为本发明的进一步改进，步骤s3中，非负脑区功率稀疏支撑向量α初始化为：
[0037][0038]
其中(
·
)
t
表示矩阵转置，||
·
||f表示frobenius范数，表示一个n
roi
×
1维全1向量。
[0039]
作为本发明的进一步改进，步骤s3中的迭代参数包括：误差比率阈值ε和最高迭代次数nm；误差比率阈值ε设置为大于0且小于1；最高迭代次数nm设置为一个大于1的整数。
[0040]
作为本发明的进一步改进，步骤s4中：
[0041]
所述非负脑区功率系数γ的后验均值向量μ
γ
如下：
[0042][0043]
其中，(
·
)
t
表示矩阵转置；
[0044]
所述非负脑区功率系数γ的后验协方差矩阵σ
γ
如下：
[0045][0046]
其中，diag(
·
)表示对角化矩阵，(
·
)-1
表示矩阵求逆。
[0047]
作为本发明的进一步改进，步骤s5中：
[0048]
非负脑区功率稀疏支撑向量α的更新值α
new
的第n个元素计算公式为：
[0049][0050]
其中，μn＝μ
γ
[n]是μ
γ
的第n个元素，是σ
γ
第(n,n)个元素的平方根，增补误差函数
[0051]
作为本发明的进一步改进，步骤s6中：
[0052]
所述非负脑区功率系数γ的一阶矩计算公式为：
[0053][0054]
其中，更新的非负脑区功率系数γ的后验协方差矩阵其中diag(
·
)表示对角化矩阵，(
·
)-1
表示矩阵求逆；
[0055]
所述皮层神经元电活动的空间协方差矩阵σ
x
(《γ》)的计算公式为：
[0056][0057]
其中，《γn》,n＝1,2,...,n
roi
表示《γ》的第n个元素。
[0058]
本发明还提供了一种基于非负块稀疏贝叶斯学习的瞬态脑电源定位系统，其包括以下模块：
[0059]
多通道eeg数据获取模块，用于利用m个电极采集对应于瞬态响应的多通道eeg，采样率为fs，共记录t个采样点，并对其预处理获得m
×
t维预处理后的多通道eeg数据y；
[0060]
空间平滑核矩阵计算模块，用于将皮层划分成n
roi
个脑区，根据这些脑区将导程场
矩l阵划分成各脑区上的子导程场矩阵ln,n＝1,2,...,n
roi
，并计算各脑区上的空间平滑核矩阵kn,n＝1,2,...,n
roi
；
[0061]
加权样本数据矩阵及平滑脑区基矩阵构造模块，用于利用预处理后的多通道eeg数据y计算样本协方差矩阵利用样本协方差矩阵各脑区上的子导程场矩阵ln和各脑区上的空间平滑核矩阵kn,n＝1,2,...,n
roi
，构造加权样本数据矩阵z以及平滑脑区基矩阵m＝1,2,...,m；
[0062]
初始化和迭代参数设置模块，用于将非负脑区功率稀疏支撑向量α初始化为α
init
，并设置迭代参数；
[0063]
非负脑区功率系数的后验均值向量和后验协方差矩阵计算模块，用于利用加权样本数据矩阵z、平滑脑区基矩阵m＝1,2,...,m以及前一非负脑区功率稀疏支撑向量α计算非负脑区功率系数γ的后验均值向量μ
γ
和后验协方差矩阵σ
γ
；
[0064]
非负脑区功率稀疏支撑向量的更新值计算模块，用于根据非负脑区功率系数γ的后验均值向量μ
γ
和后验协方差矩阵σ
γ
计算非负脑区功率稀疏支撑向量α的更新值α
new
；
[0065]
判断模块，用于判断是否满足迭代停止条件：迭代次数n≥nm或||α
new-α||2/||α||2≤ε，其中||
·
||2表示l2范数，若不满足则令α＝α
new
；若是，则回到步骤s4继续进行迭代运算；否则结束迭代；
[0066]
皮层神经元电活动的空间协方差矩阵计算模块，用于在迭代结束后，利用加权样本数据矩阵z、平滑脑区基矩阵m＝1,2,...,m以及当前非负脑区功率稀疏支撑向量α的更新值α
new
计算非负脑区功率系数γ的一阶矩《γ》，并联合各脑区上的空间平滑核矩阵kn,n＝1,2,...,n
roi
计算得到皮层神经元电活动的空间协方差矩阵σ
x
(《γ》)。
[0067]
本发明的有益效果：
[0068]
(1)本发明在瞬态脑电源定位算法的建模阶段即充分考虑电流偶极子源在皮层上的块稀疏分布先验特征，从而使瞬态脑电源定位结果更符合实际；
[0069]
(2)本发明无需预设或预知噪声功率即可进行瞬态脑电源定位；
[0070]
(3)本发明无需预设电流偶极子源的个数即可进行瞬态脑电源定位；
[0071]
(4)本发明直接绕开对噪声功率的估计，从而提升源定位效果，并且同时适用于空间非均匀噪声场景；
[0072]
(5)本发明充分利用了信号功率的非负性来提升瞬态脑电源定位效果。
[0073]
上述说明仅是本发明技术方案的概述，为了能够更清楚了解本发明的技术手段，而可依照说明书的内容予以实施，并且为了让本发明的上述和其他目的、特征和优点能够更明显易懂，以下特举较佳实施例，并配合附图，详细说明如下。
附图说明
[0074]
图1是本发明优选实施例中基于非负块稀疏贝叶斯学习的瞬态脑电源定位方法的流程图。
[0075]
图2是本发明优选实施例中基于非负块稀疏贝叶斯学习的瞬态脑电源定位方法的a
′
度量指标随着snr变化的曲线。
[0076]
图3是本发明优选实施例中基于非负块稀疏贝叶斯学习的瞬态脑电源定位方法的a
′
度量指标随着采样向量数变化的曲线。
[0077]
图4是本发明优选实施例中基于非负块稀疏贝叶斯学习的瞬态脑电源定位方法的a
′
度量指标随着激活脑区数变化的曲线。
[0078]
图5是本发明优选实施例中基于非负块稀疏贝叶斯学习的瞬态脑电源定位方法在同方差噪声条件下源定位rmse随snr变化的曲线。
[0079]
图6是本发明优选实施例中基于非负块稀疏贝叶斯学习的瞬态脑电源定位方法在同方差噪声条件下源定位rmse随采样向量数变化的曲线。
[0080]
图7是本发明优选实施例中基于非负块稀疏贝叶斯学习的瞬态脑电源定位方法在异方差噪声条件下源定位rmse随snr变化的曲线。
[0081]
图8是本发明优选实施例中基于非负块稀疏贝叶斯学习的瞬态脑电源定位方法在异方差噪声条件下源定位rmse随采样向量数变化的曲线。
[0082]
图9是本发明用于算法性能验证的经预处理和平均后的多通道p300数据在cz电极通道上的波形示意图。
[0083]
图10是本发明基于非负块稀疏贝叶斯学习的瞬态脑电源定位方法用于实际p300数据源定位结果在皮层上的分布图。
具体实施方式
[0084]
下面结合附图和具体实施例对本发明作进一步说明，以使本领域的技术人员可以更好地理解本发明并能予以实施，但所举实施例不作为对本发明的限定。
[0085]
本发明旨在基于非负块稀疏贝叶斯学习，提出了一种瞬态脑电源定位方法。本发明首先采集瞬态多通道eeg，计算样本协方差矩阵，对矩阵每列加权，可表示为laplace平滑后各脑区上非负块稀疏表示。然后设置迭代停止条件以及非负脑区功率稀疏支撑向量的初始值。而后基于非负高斯分布迭代更新非负脑区功率向量的后验均值与协方差并由此更新非负脑区功率稀疏支撑向量。最后，利用最新的非负脑区功率稀疏支撑向量给出源定位结果。本发明利用协方差向量的脑区非负块稀疏表示，结合非负高斯后验分布的期望与方差，给出瞬态eeg源定位nnbsbl方法。
[0086]
如图1所示，为本发明优选实施例中基于非负块稀疏贝叶斯学习(nnbsbl)的瞬态脑电源定位方法，其包括以下步骤：
[0087]
s1、利用m个电极采集对应于瞬态响应的多通道eeg，采样率为fs，共记录t个采样点，并对其预处理获得m
×
t维预处理后的多通道eeg数据y；同时，将皮层划分成n
roi
个脑区，根据这些脑区将导程场矩l阵划分成各脑区上的子导程场矩阵ln,n＝1,2,...,n
roi
，并计算各脑区上的空间平滑核矩阵kn,n＝1,2,...,n
roi
；
[0088]
具体地，所述计算各脑区上的空间平滑核矩阵kn,n＝1,2,...,n
roi
包括：
[0089]
第n个脑区上的空间平滑核矩阵的计算公式为n＝1,2,...,n
roi
，其中，hn表示从全脑离散拉普拉斯平滑矩阵h中抽取的对应于第n个脑区的子阵，(
·
)
t
表示矩阵转置，(
·
)-1
表示矩阵求逆。
[0090]
其中，所述全脑离散拉普拉斯平滑矩阵h的计算如下：
[0091]
对于导程场矩阵l的列所对应的nv个大脑皮层网格顶点，利用这些顶点的位置和顶点两两之间的连接情况，所得nv×nv
维全脑离散拉普拉斯平滑矩阵h的计算方法为：其第(i,i)个元素当i≠j且第i个顶点和第j个顶点间有连接时h的第(i,j)个元素当i≠j且第i个顶点和第j个顶点间无连接时h的第(i,j)个元素h[i,j]＝0，其中表示与第i个顶点有连接的所有邻居顶点的距离的平均值，表示第i个顶点与跟它有连接的所有邻居顶点的距离的倒数的平均值，ni表示与第i个顶点有连接的所有邻居顶点的个数，d
ij
表示在第i个顶点和第j个顶点间有连接的情况下这两个顶点间的距离。
[0092]
可选的，所述预处理包括基线校正、带通滤波以及伪迹剔除等。
[0093]
s2、利用预处理后的多通道eeg数据y计算样本协方差矩阵利用样本协方差矩阵各脑区上的子导程场矩阵ln和各脑区上的空间平滑核矩阵kn,n＝1,2,...,n
roi
，构造加权样本数据矩阵z以及平滑脑区基矩阵m＝1,2,...,m；
[0094]
具体地，所述加权样本数据矩阵z＝[z1,z2,...,zm]，其第m列的计算公式为：
[0095][0096]
其中，第m列样本协方差向量的协方差矩阵其中(
·
)
t
表示矩阵转置，表示的第(m,m)个元素，表示的第m列，jm为第m个权重矩阵，其计算公式为j1＝[0
m-1
,i
m-1
]，jm＝[i
m-1
[:,1:m-1],0
m-1
,i
m-1
[:,m:m-1]],2≤m≤m-1，jm＝[i
m-1
,0
m-1
]，0
m-1
表示一个(m-1)
×
1维全1向量，i
m-1
表示一个(m-1)
×
(m-1)维单位矩阵；
[0097]
平滑脑区基矩阵m＝1,2,...,m的计算公式为：
[0098][0099]
其中，第n个空间平滑成分子阵n＝1,2,...,n
roi
，an[:,m]表示矩阵an的第m列。
[0100]
s3、将非负脑区功率稀疏支撑向量α初始化为α
init
，并设置迭代参数；
[0101]
具体地，非负脑区功率稀疏支撑向量α初始化为：
[0102][0103]
其中(
·
)
t
表示矩阵转置，||
·
||f表示frobenius范数，表示一个n
roi
×
1维全1向量。
[0104]
具体地，迭代参数包括：误差比率阈值ε和最高迭代次数nm；误差比率阈值ε设置为大于0且小于1，且数值尽可能取小；最高迭代次数nm设置为一个大于1的整数。
[0105]
s4、利用加权样本数据矩阵z、平滑脑区基矩阵m＝1,2,...,m以及前一非负脑
区功率稀疏支撑向量α计算非负脑区功率系数γ的后验均值向量μ
γ
和后验协方差矩阵σ
γ
；
[0106]
具体地，所述非负脑区功率系数γ的后验均值向量μ
γ
如下：
[0107][0108]
其中，(
·
)
t
表示矩阵转置；
[0109]
所述非负脑区功率系数γ的后验协方差矩阵σ
γ
如下：
[0110][0111]
其中，diag(
·
)表示对角化矩阵，(
·
)-1
表示矩阵求逆。
[0112]
s5、根据非负脑区功率系数γ的后验均值向量μ
γ
和后验协方差矩阵σ
γ
计算非负脑区功率稀疏支撑向量α的更新值α
new
；
[0113]
具体地，非负脑区功率稀疏支撑向量α的更新值α
new
的第n个元素计算公式为：
[0114][0115]
其中，μn＝μ
γ
[n]是μ
γ
的第n个元素，是σ
γ
第(n,n)个元素的平方根，增补误差函数
[0116]
s6、判断是否满足迭代停止条件：迭代次数n≥nm或||α
new-α||2/||α||2≤ε，其中||
·
||2表示l2范数，若不满足则令α＝α
new
；若是，则回到步骤s4继续进行迭代运算；否则结束迭代，利用加权样本数据矩阵z、平滑脑区基矩阵m＝1,2,...,m以及当前非负脑区功率稀疏支撑向量α的更新值α
new
计算非负脑区功率系数γ的一阶矩〈γ》，并联合各脑区上的空间平滑核矩阵kn,n＝1,2,...,n
roi
计算得到皮层神经元电活动的空间协方差矩阵σ
x
(〈γ》)。即所要求的瞬态脑电源定位结果。
[0117]
具体地，所述非负脑区功率系数γ的一阶矩计算公式为：
[0118][0119]
其中，更新的非负脑区功率系数γ的后验协方差矩阵其中diag(
·
)表示对角化矩阵，(
·
)-1
表示矩阵求逆；
[0120]
所述皮层神经元电活动的空间协方差矩阵σ
x
(《γ》)的计算公式为：
[0121]
[0122]
其中，《γn》,n＝1,2,...,n
roi
表示《γ》的第n个元素。
[0123]
本发明在瞬态脑电源定位算法的建模阶段即充分考虑电流偶极子源在皮层上的块稀疏分布先验特征，从而使瞬态脑电源定位结果更符合实际；本发明无需预设或预知噪声功率即可进行瞬态脑电源定位；本发明无需预设电流偶极子源的个数即可进行瞬态脑电源定位；本发明直接绕开对噪声功率的估计，从而提升源定位效果，并且同时适用于空间非均匀噪声场景；本发明充分利用了信号功率的非负性来提升瞬态脑电源定位效果。
[0124]
本发明基于非负块稀疏贝叶斯学习的瞬态脑电源定位方法涉及的理论推导过程如下：
[0125]
(1)瞬态eeg源定位问题的概率生成模型
[0126]
假设某种自发的或者由某种范式的刺激诱发的瞬态eeg信号，由头皮上分布的多导联电极采集，记录电极的通道数为m，记录时长为t
rec
，采样频率为fs。数字化的eeg采样数据随后被通过一个包括基线校正、工频干扰抑制以及伪迹剔除在内的预处理模块。预处理后的数据可表示为一个数据矩阵y＝[y1,y2,...,y
t
]，其中y
t
为第t个采样向量，该向量的第m个元素由第m个电极通道采集得到，m＝1,2,...,m，并且采样向量的个数为t＝t
recfs
。预处理后的采样向量y
t
的构成可由下式表示：
[0127]yt
＝lx
t
n
t
,t＝1,2,...,t,
ꢀꢀꢀ
(1)
[0128]
其中：l为导程场矩阵，其各元素表示由大脑皮层网格各顶点体素上的电流偶极子传导到头皮上各电极处的权重，即电流偶极子由皮层传导到头皮上的映射关系；x
t
表示与所采集的eeg相对应的皮层上的电流偶极子源活动，其每个元素为导程场矩阵对应列所指定的皮层网格顶点体素上的电流偶极子源强度；n
t
是噪声成分，有传感器上的电噪声和非定位目标的自发脑电成分构成。我们在分布式源模型下研究瞬态eeg源定位问题，使用一个包含nv个顶点的网格完全覆盖大脑皮层表面，每一个顶点位置(体素)表示电流偶极子源可能存在的位置。在分布式源模型下，由于锥体细胞的排列特点，各体素上的电流偶极子的方向是固定的(垂直于皮层表面)，因此融合了偶极子方向信息后式(1)中的导程场矩阵l包含nv列，nv维向量x
t
表示这nv个体素上的电流偶极子强度，其中不存在源的体素上的电流偶极子强度为0。当噪声向量假设服从均值为零向量、协方差矩阵为λ的高斯分布时，可给出y
t
的条件分布表示为：
[0129]
p(y
t
|x
t
,λ)＝n(y
t
|lx
t
,λ),t＝1,2,...,t,
ꢀꢀꢀ
(2)
[0130]
其中n(
·
|μ,σ)表示均值向量为μ、协方差矩阵为σ的高斯分布的概率密度函数。噪声随机向量n
t
经常被简单地建模为同方差的，在本研究中为了使问题更贴近于实际，n
t
被假设为异方差的，即其协方差矩阵λ是对角矩阵且对角元素可以各不相等。
[0131]
由式(2)，eeg逆问题的求解目标在于，在未知噪声协方差矩阵λ的条件下，由预处理后的eeg数据矩阵y，推算出各采样时刻上的电流偶极子强度向量所排列而成的矩阵x＝[x1,x2,...,x
t
]。很显然，由于nv＞＞m，该问题是一个严重的病态问题：基矩阵列数远大于行数的欠定问题。由于电流偶极子源在皮层上的稀疏分布特性，各种稀疏先验对应的正则化措施被引入到了分布式源模型下的eeg源定位求解方法中，特别是在稀疏贝叶斯学习(sbl)这一类方法中的应用尤为突出。不同于大多数sbl方法假设产生eeg的电流偶极子源只存在于皮层网格的少数顶点上，即满足元稀疏(atom-wise sparsity)，近年来电流偶极子源在皮层网格上的元稀疏特征渐渐被块稀疏(block-wise sparsity)所替代。块稀疏假
设的出发点在于，eeg信号应源于神经功能核团，每个核团包含网格上相邻的多个顶点，这些神经功能核团内的皮层顶点位置是由某种功能脑区分割方式所确定的。块稀疏包含两层含义：块内源之间的相关性和块层面上的稀疏性。其具体的实现方式为：参考某种大脑图谱(atlas)对皮层进行功能区域分割，获得不同功能的块，它们包含的顶点数可以各不相同，块内源之间的相关性可通过某种空间平滑核来定义。以destrieux大脑图谱为例，可将皮层分割为n
roi
＝147个感兴趣区域(roi)并对皮层全覆盖，不同的roi之间满足块稀疏性。在每一个roi中，定义对应的空间平滑核矩阵kn,n＝1,2,...,n
roi
，即块内各顶点上源之间的协方差矩阵，用以规定块内源之间的相关性。由此，可得零均值高斯分布随机向量x
t
的协方差矩阵表示为：
[0132][0133]
其中，各roi上的功率系数γn,n＝1,2,...,n
roi
为非负乘性系数，而γ
nkn
则为第n个roi块内源的协方差矩阵，γn为零或者很小意味着该roi块内不存在电流偶极子源。类似于已有算法bsbl-2s的设计过程，在本研究中，空间平滑核矩阵定义为：n＝1,2,...,n
roi
，其中hn表示从全脑离散拉普拉斯平滑矩阵h中抽取的对应于第n个脑区的子阵，(
·
)
t
表示矩阵转置，(
·
)-1
表示矩阵求逆。全脑离散拉普拉斯平滑矩阵h的计算过程为：对于导程场矩阵l的列所对应的nv个大脑皮层网格顶点，利用这些顶点的位置和顶点两两之间的连接情况，所得nv×nv
维全脑离散拉普拉斯平滑矩阵h的计算方法为，其第(i,i)个元素当i≠j且第i个顶点和第j个顶点间有连接时h的第(i,j)个元素当i≠j且第i个顶点和第j个顶点间无连接时h的第(i,j)个元素h[i,j]＝0，其中表示第i个顶点与跟它有连接的所有邻居顶点的距离的平均值，表示第i个顶点与跟它有连接的所有邻居顶点的距离的倒数的平均值，ni表示与第i个顶点有连接的所有邻居顶点的个数，d
ij
表示在第i个顶点和第j个顶点间有连接的情况下这两个顶点间的距离。由式(3)可见，t个采样向量共用一个块稀疏支撑集，这要求在整个eeg数据获取时间范围内皮层上的电流偶极子源需保持不变。已有研究表明，eeg信号可划分成多个似平稳的微状态，每个似平稳微状态的持续时长为80-120ms。因此，在本研究中，t个采样向量的总采样持续时长需≤120ms。
[0134]
目前已有的基于bsbl的eeg源定位方法都是在式(2)与式(3)所构成的概率生成模型基础上而推导得出的，这其中需要知晓噪声协方差矩阵λ的信息，包括三种方式：噪声协方差矩阵λ本身已知、噪声协方差矩阵λ由基线数据估算或者噪声协方差矩阵λ与同时进行估计。在本研究中，我们另辟蹊径绕开对噪声协方差矩阵λ的需求：考虑eeg数据的样本协方差矩阵的分布而非如式(2)所示的直接eeg数据本身的分布，该考虑的分布模型不再直接以噪声协方差矩阵λ为条件。y
t
的协方差矩阵定义为：
[0135][0136]
其中，ln是导程场矩阵l的子阵，由对应于第n个roi块的所有列构成，并且an定义为若对角的噪声协方差矩阵表示为其对角元素可相等或不相等，可以对ry的第m列左乘一个对应的加权矩阵jm,m＝1,2,...,m，从而获得无噪向量即：
[0137][0138]
其中，加权矩阵的构造方法为：j1＝[0
m-1
,i
m-1
],jm＝[i
m-1
,0
m-1
],以及jm＝[i
m-1
[:,1:m-1],0
m-1
,i
m-1
[:,m:m-1]],2≤m≤m-1；0
m-1
表示(m-1)
×
1维全零向量；i
m-1
表示(m-1)
×
(m-1)维单位矩阵。式(5)中的无噪向量可进一步转化为线性表示方式，即：
[0139][0140]
其中，构造的基矩阵通过式(6)中的线性表示形式，复杂的块稀疏转化为简单的元稀疏形式，其中m个向量共用同一个稀疏支撑集，且共用的稀疏向量γ是非负的。
[0141]
在实际中，只能通过有限个样本向量来估计多通道eeg信号的协方差矩阵ry，即样本协方差矩阵其估计方法为基于已有文献，可得到实信号的样本协方差矩阵的概率分布为：
[0142][0143]
其中，vec(
·
)表示对矩阵向量化，即将矩阵的所有列依次堆叠构成一个长向量。由样本协方差矩阵进一步定义可得到：
[0144][0145]
其中的协方差矩阵定义为若定义以及由于当m≠n时有zm|γ条件独立于zn|γ，可得z的条件分布表示为：
[0146][0147]
其中，
[0148][0149]
非负稀疏向量γ控制着每个roi块内电流偶极子源的功率，其中γ中数值较大的
元素意味着对应roi块内存在偶极子源。为了进一步在所研究算法中利用到γ的非负稀疏特性，使用一个基于非负参数向量α的非负高斯先验对其进行建模，可得：
[0150][0151]
其中，n

(
·
|μ,σ)表示均值向量为μ、协方差矩阵为σ的非负高斯分布的概率密度函数。
[0152]
由式(9)和式(10)，可得该概率生成模型的联合分布的概率密度函数，表示为
[0153]
p(z,γ；α)＝p(z|γ)p(γ；α).
ꢀꢀꢀ
(11)
[0154]
基于上述概率生成模型，利用我们所研究的nnbsbl算法，可首先推算出γ的最大后验估计，然后计算出σ
x
(γ)的最大后验估计，其对角线元素即对应着eeg源定位估计结果。
[0155]
(2)基于em的nnbsbl算法推导
[0156]
本研究所得的nnbsbl算法是基于em思想推导得来的，其中非负稀疏向量γ被视为隐变量向量。联合利用式(9)与式(10)，可推导出隐变量向量γ的后验分布仍符合非负高斯分布，其表示为：
[0157]
p(γ|z；α)＝n

(γ|μ
γ
,σ
γ
),
ꢀꢀꢀ
(12)
[0158]
其中，
[0159]
em中的e-step是用来计算已知后验概率分布p(γ|z；α)的条件下log p(z,γ；α)的期望，从而得到以α为自变量的代价函数q(α)，表示为：
[0160][0161]
其中，与分别是γ在已知后验概率分布条件下的一阶矩和二阶矩。《γ》与《γγ
t
》的计算公式以及证明过程已由之前所发表论文给出，表示为：
[0162][0163]
以及，
[0164][0165]
其中，其中，累积概率补余误差函数特别的，在计算过程中若erfc(
·
)函数值接近于0，可得《γn》
→
0，以及在m-step中，通过令代价函数q(α)的导数为0，可得到非负参数向量α的更新公式为：
[0166][0167]
在本发明所设计的nnbsbl算法中，当时，可直接令
[0168]
由式(8)与式(9)可知，在构造cm与zm的过程中，ry是未知的，因为计算ry仍需要知道噪声协方差矩阵λ。在前人的研究中，可知在构造cm与zm的过程中利用替代ry仍可获得可观的效果，这个思想也在本研究的nnbsbl算法中得到采纳，因此算法无需知道或者估计噪声协方差矩阵λ。在nnbsbl算法中唯一需要迭代更新的是非负参数向量α，其迭代的初始值由下式给出：
[0169][0170]
其中，1
nroi
表示n
roi
×
1维全1向量。迭代停止条件可由两种方式给出：基于代价函数q(α)或者基于非负参数向量α本身在两次相邻迭代过程更新得到的值。考虑到q(α)中的《γmγn》
m≠n
计算过程可能导致较大的计算量，推荐采用后者：当||α
(t)-α
(t-1)
||2/||α
(t-1)
||2＜ε或者预设的最大迭代次数达到时迭代停止，其中α
(t)
是在第t次迭代时α的更新值，ε
→
0是一个预设的数值很小的比例超参数。当迭代停止时，参考式(14)，可得γ的最大后验估计《γ》，进而可以计算出x
t
的最大后验协方差矩阵σ
x
(《γ》)，其对角线元素即对应值当前eeg信号的电流偶极子源在皮层上的功率分布，故可得出瞬态eeg的皮层源定位结果。
[0171]
为了验证本发明的有效性，下面进行基于非负块稀疏贝叶斯学习的瞬态脑电源定
位方法性能验证仿真实验：
[0172]
a.计算头模型和导程场矩阵
[0173]
本发明使用mni/icbm152mri模板分割头模型并计算与之对应的导程场矩阵。借助brainstorm软件，本发明使用边界元法计算由头皮、头骨和皮质组成的3层头模型，构造了包含1500个顶点的三角网格全覆盖皮质表面，这些顶点表示电流偶极子可能存在的皮层表面所有的空间位置。假设电流偶极子的方向垂直于大脑皮层表面，结合德国easycap公司生产的64导联电极帽“64ch-braincap-feba”(10-20国际标准导联系统，包含参考电极在内的eeg信号采集电极数为m＝62)所对应的电极分布信息，假设采用全脑平均作为参考，通过openmeeg软件计算出基于真实头模型的导程场矩阵l，其维度为62
×
1500。
[0174]
b.瞬态eeg信号仿真产生
[0175]
基于destrieux大脑图谱，将皮层分割为n
roi
＝147个感兴趣区域，并计算对应的空间平滑核矩阵kn,n＝1,2,...,n
roi
。选取其中k个脑区作为存在偶极子源的区域，在这些脑区上产生均值为零向量、协方差矩阵为γ
nkn
的高斯分布信号，采样点数为t，而未被选择的脑区上信号值统一置为0。不失一般性和便利性，在此假设存在偶极子源的区域上功率系数γn都相等且等于1。将导程场矩阵l与产生的信号矩阵s相乘，并加上噪声数据矩阵n，即获得瞬态eeg信号仿真数据y。噪声数据矩阵维度为62
×
t，服从零均值高斯分布，其协方差矩阵λ为对角阵，对角线元素可全部相等也可各不相同。信噪比定义为其中||
·
||f表示frobenius范数。为了使算法应用在仿真信号与实采信号上无差别，在应用算法前统一对瞬态eeg信号数据y作归一化处理。在本实施例中迭代过程中的参数设置为：误差比率阈值ε＝10-4
，最高迭代次数nm＝1000。
[0176]
c.瞬态eeg仿真信号下皮层源定位方法性能的a
′
度量指标对比
[0177]
首先验证本研究所提出的nnbsbl算法在皮层源定位的成功概率方面的性能，并与已有的典型算法及最新算法作对比。源定位成功概率由自由响应受试者工作特征(froc)曲线体现，froc曲线的纵轴为定位到正确源位置的脑区数占所有检测到存在源的脑区的比例，横轴为将不存在源的脑区检测为存在源的脑区的概率。一般通过froc曲线下的面积来衡量算法在皮层源定位的成功概率方面的优劣，曲线下面积越大意味着算法越容易成功定位到存在源的脑区且误定位到不存在源的脑区的概率小。一般利用a
′
度量来估计froc曲线下的面积，其定义为：a
′
＝(h
r-fr)/2 1/2，其中hr为击中比率，即被定位算法检测到的真实存在电流偶极子源的脑区占所有存在电流偶极子源的脑区的比率；fr为误检比率，即被定位算法检测到但实际并不存在电流偶极子源的脑区占所有不存在电流偶极子源的脑区的比率。定义一个脑区被源定位算法“检测”到存在偶极子源为：该脑区上的功率系数的最大后验估计值《γn》位于向量《γ》所有元素从大到小排列的前10％，并且《γn》的值不低于《γ》最大元素值的1％。
[0178]
下面开展实验计算所提出的nnbsbl算法在不同条件变化下的a
′
度量指标，在相同的实验条件下，还同时执行了三种对比算法并计算a
′
度量指标：(1)sloreta；(2)bsbl-2s(3)包含em噪声鲁棒估计的smooth-champagne算法(表示为champagneem)。图2展示了nnbsbl以及对比算法的a
′
度量指标随着snr变化的曲线，其中采样向量数t＝80，噪声为同
方差噪声，每个snr条件下经过了1000次蒙特卡洛实验结果平均，每个试次随机选择5个脑区。图3展示了nnbsbl以及对比算法的a
′
度量指标随着采样向量数变化的曲线，其中snr＝5db，噪声为同方差噪声，每个采样向量数条件下经过了1000次蒙特卡洛实验结果平均，每个试次随机选择5个脑区作为存在电流偶极子源的脑区。图4展示了nnbsbl以及对比算法的a
′
度量指标随着存在电流偶极子源的脑区(激活脑区)数变化的曲线，其中snr＝5db，采样向量数t＝80，噪声为同方差噪声，每个激活脑区数条件下经过了1000次蒙特卡洛实验结果平均。由图2、图3和图4可见，在这些设置的变化条件下，相对于其他三种算法，本发明的nnbsbl算法始终获得最高的a
′
度量值，意味着本发明在瞬态eeg源定位问题上以更高的概率定位到存在偶极子源的脑区，并且其误定位到不存在偶极子源的脑区的概率更小。
[0179]
d.瞬态eeg仿真信号下皮层源定位均方根误差(rmse)指标对比
[0180]
接下来，统计本研究所提出的nnbsbl算法在源定位误差方面的性能，并与bsbl-2s和champagneem算法作对比。皮层源定位均方根误差(rmse)定义为其中k为存在电流偶极子源的脑区个数，n
mc
为蒙特卡洛实验次数，rk为存在电流偶极子源的第k个脑区的中心点坐标，为第n次蒙特卡洛实验由源定位算法所定位到的与存在电流偶极子源的第k个脑区相匹配的脑区的中心点坐标，这里“匹配”定义为被源定位算法“检测”到的脑区集合中与真实存在电流偶极子源的脑区距离最近的脑区，当被源定位算法“检测”到的脑区数少于k时，该集合更改为《γn》的值最大的k个脑区。
[0181]
图5展示的是同方差噪声条件下源定位rmse随snr变化的曲线，其中采样向量数t＝80，蒙特卡洛实验试次1000次，每个蒙特卡洛试次随机选择k＝3个脑区作为存在电流偶极子源的脑区。图6展示的是同方差噪声条件下源定位rmse随采样向量数变化的曲线，其中snr＝5db，蒙特卡洛实验试次1000次，每个蒙特卡洛试次随机选择k＝3个脑区作为存在电流偶极子源的脑区。接下来研究算法在异方差噪声条件下的性能。异方差噪声的参数方法为：在每个蒙特卡洛试次中，第m个电极通道上的噪声方差为其中为噪声系数，ξm为在[0.5,1.5]区间内均匀分布的随机数。图7展示的是异方差噪声条件下源定位rmse随snr变化的曲线，其中采样向量数t＝80，蒙特卡洛实验试次1000次，每个蒙特卡洛试次随机选择k＝3个脑区作为存在电流偶极子源的脑区。图8展示的是异方差噪声条件下源定位rmse随采样向量数变化的曲线，其中snr＝5db，蒙特卡洛实验试次1000次，每个蒙特卡洛试次随机选择k＝3个脑区作为存在电流偶极子源的脑区。由图5、图6、图7与图8可见，无论是在同方差噪声还是异方差噪声条件下，本发明的nnbsbl算法在中高信噪比条件下相比对比算法都取得了更小的定位误差，且随着采样向量数的增加误差逐渐下降。
[0182]
综合以上仿真结果及分析，本发明的nnbsbl算法在性能上区别于其他典型算法或者最新算法的主要特征是：首先，nnbsbl算法对瞬态eeg的电流偶极子源的检测性能更佳，能够以更高的概率定位到存在偶极子源的脑区，并且其误定位到不存在偶极子源的脑区的概率更小；然后，nnbsbl算法在中、高信噪比条件下可取得了更小的定位误差，且定位误差随着采样向量数的增加还将逐渐下降。
[0183]
下面，将本发明基于非负块稀疏贝叶斯学习的瞬态脑电源定位方法应用于p300源定位。
[0184]
作为一种典型的瞬态诱发响应，p300在认知心理学、临床诊断和研究、脑机接口开发等领域应用广泛。本发明使用实际采集的被试的p300数据来验证本发明所提出算法的瞬态eeg源定位效果。
[0185]
实验过程描述如下。被试端坐在屏幕前保持专注，阅读完实验须知后点击鼠标左键开始实验。实验采用的图片呈现方式为出现/消失模式，图片为6
×
6字母表，每次仅出现字母表的随机行或者随机列，一个图片刺激的周期为800ms，其中出现125ms，消失675ms。当6行6列的字母表全部呈现一次视为一轮(共12次)。每次实验包含20轮，当字母“s”所在的行或者列出现的时候被试进行心理计数(每次实验正确计数应为40次)，共连续进行4次这样的实验，每次实验之间提供3分钟的休息间隙(字母“s”出现的总次数应为160次)。整个过程中同时记录的数据包括标准64导联分布电极采集到的电生理信号、用于去除眼电的垂直/水平眼电信号以及屏幕出现靶目标的时间标记。
[0186]
实验使用德国easycap公司生产的64导联电极帽“64ch-braincap-feba”，以所有电极的平均值为参考，使用德国brain products公司的brainamp(s/namp 13122084mrplus)放大器和brainvision recorder软件记录eeg信号。在eeg数据预处理过程中，通过analyzer2.0软件借助采集到的eog信号进行eog伪迹消除等预处理后，参考实验过程中所打的标记，对数据叠加平均得到最终的p300波形数据。图9所示为本实施例所采用的某被试最终平均得到的p300数据在cz电极通道上的波形示意图。对该多通道p300数据300ms处左右80ms长度数据，执行本发明的基于非负块稀疏贝叶斯学习的瞬态脑电源定位方法，得到的定位结果在皮层上的分布如图10所示。由图10可见，本发明所提出的瞬态eeg源定位方法，应用于真实的p300信号的皮层源定位任务时，其源定位结果与p300的综述论文“john polich,updating p300:an integrative theory ofp3a and p3b,clinical neurophysiology,2007”中与功能磁共振(fmri)同步开展所得出的p300皮层上源分布的结论非常吻合，表明了本发明所提出的瞬态eeg源定位方法不仅在仿真数据上表现优异，还在实际应用中证明了其有效性。
[0187]
本发明优选实施例还公开了一种基于非负块稀疏贝叶斯学习的瞬态脑电源定位系统，其包括以下模块：
[0188]
多通道eeg数据获取模块，用于利用m个电极采集对应于瞬态响应的多通道eeg，采样率为fs，共记录t个采样点，并对其预处理获得m
×
t维预处理后的多通道eeg数据y；
[0189]
空间平滑核矩阵计算模块，用于将皮层划分成n
roi
个脑区，根据这些脑区将导程场矩l阵划分成各脑区上的子导程场矩阵ln,n＝1,2,...,n
roi
，并计算各脑区上的空间平滑核矩阵kn,n＝1,2,...,n
roi
；
[0190]
加权样本数据矩阵及平滑脑区基矩阵构造模块，用于利用预处理后的多通道eeg数据y计算样本协方差矩阵利用样本协方差矩阵各脑区上的子导程场矩阵ln和各脑区上的空间平滑核矩阵kn,n＝1,2,...,n
roi
，构造加权样本数据矩阵z以及平滑脑区基矩阵m＝1,2,...,m；
[0191]
初始化和迭代参数设置模块，用于将非负脑区功率稀疏支撑向量α初始化为α
init
，
并设置迭代参数；
[0192]
非负脑区功率系数的后验均值向量和后验协方差矩阵计算模块，用于利用加权样本数据矩阵z、平滑脑区基矩阵m＝1,2,...,m以及前一非负脑区功率稀疏支撑向量α计算非负脑区功率系数γ的后验均值向量μ
γ
和后验协方差矩阵σ
γ
；
[0193]
非负脑区功率稀疏支撑向量的更新值计算模块，用于根据非负脑区功率系数γ的后验均值向量μ
γ
和后验协方差矩阵σ
γ
计算非负脑区功率稀疏支撑向量α的更新值α
new
；
[0194]
判断模块，用于判断是否满足迭代停止条件：迭代次数n≥nm或||α
new-α||2/||α||2≤ε，其中||
·
||2表示l2范数，若不满足则令α＝α
new
；若是，则回到步骤s4继续进行迭代运算；否则结束迭代；
[0195]
皮层神经元电活动的空间协方差矩阵计算模块，用于在迭代结束后，利用加权样本数据矩阵z、平滑脑区基矩阵m＝1,2,...,m以及当前非负脑区功率稀疏支撑向量α的更新值α
new
计算非负脑区功率系数γ的一阶矩《γ》，并联合各脑区上的空间平滑核矩阵kn,n＝1,2,...,n
roi
计算得到皮层神经元电活动的空间协方差矩阵σ
x
(《γ》)。
[0196]
本实施例中基于非负块稀疏贝叶斯学习的瞬态脑电源定位系统用于实现前述的基于非负块稀疏贝叶斯学习的瞬态脑电源定位方法，因此该系统的具体实施方式可见前文中的基于非负块稀疏贝叶斯学习的瞬态脑电源定位方法的实施例部分，所以，其具体实施方式可以参照相应的各个部分实施例的描述，在此不再展开介绍。
[0197]
另外，由于本实施例的基于非负块稀疏贝叶斯学习的瞬态脑电源定位系统用于实现前述的基于非负块稀疏贝叶斯学习的瞬态脑电源定位方法，因此其作用与上述方法的作用相对应，这里不再赘述。
[0198]
以上实施例仅是为充分说明本发明而所举的较佳的实施例，本发明的保护范围不限于此。本技术领域的技术人员在本发明基础上所作的等同替代或变换，均在本发明的保护范围之内。本发明的保护范围以权利要求书为准。