一种基于贝叶斯估计的多层次多尺度层析成像方法与流程

本发明属于脑科学与信息技术交叉学科领域，涉及一种脑磁图(magnetoencephalography，meg)、脑电场(electroencephalogram，eeg)成像的方法，具体涉及一种基于贝叶斯估计的多层次多尺度层析成像方法。

背景技术：

人脑活动的成像在认知神经科学的研究中起着重要的作用，其可以促进人类对自身复杂活动时神经元工作机制的理解，在临床医学应用中，脑成像同样扮演着重要角色，特别对脑疾病患者比如脑肿瘤和癫痫手术的诊断、指导以及病灶的切除。然而，当前对脑源位置定位和活动时间序列估计是一个极具挑战的问题，因为潜在脑活动的源数目可以为成千上万，而脑外传感器数目仅仅只有数百个。脑磁图(meg)和脑电场(eeg)是目前存在的两种无损的对脑电生理活动进行测量和成像方法，由于meg和eeg分别测量脑内一簇神经元活动的电流在脑外产生的磁场和电场，其具有互补性。因此，基于meg和eeg对脑活动的动态研究有着重要意义。本发明主要针对当前层析方法无法对多尺度源活动进行准确成像的问题，结合脑功能结构以及解剖结构对源空间的区域分割，提出了基于贝叶斯的多尺度源活动层析成像方法。近些年，众多学者通过meg和eeg对脑源活动进行分析，提出了许多脑电磁源成像的方法，但是针对诸如一簇多个相邻偶极子活动或者是一簇多个相邻偶极子与独立偶极子同时活动的问题，始终无法找到有效解决方法。

技术实现要素：

为了解决上述技术问题，本发明提出了一种全新的多层次多尺度meg、eeg脑源重构通用模型的方法，可以看作是源重构方法的多层次多尺度扩展，旨在实现对多尺度源活动的重构。

本发明所采用的技术方案是：一种基于贝叶斯估计的多层次多尺度层析成像方法，其特征在于，包括以下步骤：

步骤1：多层次多尺度脑源概率模型的生成；

步骤2：脑源活动分布及其区域活动分布的估计；

步骤3：贝叶斯超参数协方差分布估计方法；

步骤4：据meg/eeg多尺度源活动层析成像。

本发明首先将脑源空间分解为相互独立的体元，每个体元作为一个潜在的活动源。然后将所有的体元根据人脑解剖结构或者功能区分解为若干个区域，每个区域对应着不同的分布参数。假设每个体元的活动分布由两部分组成：本身固有活动分布和区域活动分布，同时假设位于同一个区域的体元有着相同的区域活动协方差成分，每个源的活动与传感器采集数据之间通过标准leadfield导向矩阵来联系，而leadfield是由大脑的结构和导电性来决定。基于贝叶斯及其凸函数理论的方法对脑外采样数据进行分析，估计所有体元本征固有活动分布以及区域活动分布的协方差成分。基于上述模型框架提出两种方法：一种是体元最终活动的协方差成分只通过区域活动决定称为tling-champagne，另一种是体元的活动是由体元本征固有的活动和区域活动分布相加来决定，称为tree-champagne；最后通过模拟数据和真实脑数据对其进行性能分析。

相对于现有技术，本发明的有益效果是：针对现有的技术始终无法有效解决如一簇多个相邻偶极子活动或者是一簇多个相邻偶极子与独立偶极子同时活动的问题，本发明提出一个全新的多层次多尺度meg、eeg脑源重构通用模型。该模型首先将脑源空间分解为相互独立的体元，每个体元作为一个潜在的活动源。然后将所有的体元根据人脑解剖结构或者功能区分解为若干个区域，每个区域对应着不同的分布参数，以实现对多尺度源活动的重构。

附图说明

图1为本发明实施例的流程图。

具体实施方式

为了便于本领域普通技术人员理解和实施本发明，下面结合附图及实施例对本发明作进一步的详细描述，应当理解，此处所描述的实施示例仅用于说明和解释本发明，并不用于限定本发明。

基于贝叶斯估计的多层次多尺度层析成像方法主要包括对采样数据meg/eeg模型的生成，脑源活动分布及其区域活动分布的估计，贝叶斯超参数协方差分布估计和多尺度源活动层析成像方法的求解。

请见图1，本发明提供的一种基于贝叶斯估计的多层次多尺度层析成像方法，包括以下步骤：

步骤1：多层次多尺度脑源概率模型的生成；

具体步骤包括：

步骤1.1：假设采样数据meg/eeg是由脑内所有体元(网格点)的活动和所有脑区的活动以及传感器背景干扰噪声(包括脑、生物磁信息、环境干扰以及传感器噪声)组成。其数据采集的生成模型描述为：

其中，是传感器在时刻t采集时间序列，dy为传感器采集点的数目，n是脑区分割产生的潜在脑源活动的数目，为第i个体元(潜在源活动的位置)leadfield矩阵，li的第k列表示假设第i个体元为单位强度活动偶极子时(固定方向为k)产生的外面传感器所能采集到的数据，ε表示噪声和干扰信息，为了简单计算，这里假设噪声与干扰为不相关的高斯分布si(t)表示第i个体元在时刻t的活动。

步骤1.2：预先根据脑解剖结构或者功能分成r个区域(片)，第j个脑区包含pj个体元活动，并且分割脑区之间相互独立互不重叠，且每一个体元只属于一个脑区。假设第j个区域的在时刻t的平均活动为可以抽象为一个dc维度的神经元电流偶极子，对于第j个脑区的投影：

si(t)＝gjzj(t)+vi(t)(2)

其中gi为j脑区分布和第i个体元在时刻t的活动si(t)之间的加权矩阵，本发明假设vi(t)为体元i在t时刻的本征固有活动强度，与其它的体元活动之间具有相互独立性；zj(t)表示第j个区域的在时刻t的平均活动。

将公式(2)带入公式(1)中可得：

其中，表示为位于j区域的所有体元对应的leadfield矩阵求和，有ωj表示第j个区域的所有体元。

定义扩展(体元扩展)leadfield矩阵为h，有：

这里当i＝1,…,n时，有hi＝li；当i＝n+1,…,n+r时，有定义时刻t的扩展体元活动序列为：

x(t)＝[v1(t),…,vn(t),z1(t),…,zr(t)]^t

＝[x1(t),x2(t),…,xr+n(t)]^t(5)

当i＝1,…,n时，有xi(t)＝vi(t)；当i＝n+1,…,n+r时，有xi(t)＝zi(t)，将公式(4)和公式(5)带入公式(3)可得：

y(t)＝hx(t)+ε(6)

公式(6)即为发明所提出的通用新模型(多层次多尺度脑源概率模型)。

对公式的简化表达，在tk采样数据y(tk)表示为yk，同样x(tk)表示为xk。方法对源活动的重构在时域和空间域进行，故有时域活动序列为x＝[x1,x2,…,xk]，采样数据为y＝[y1,y2,…,yk]。

定义γi为xi的dc×dc先验协方差矩阵，且定义γ为dc(n+r)×dc(n+r)的块对角矩阵为：

公式(6)中模型的先验分布假设为：

将噪声和干扰活动分布带入可得条件概率分布为：

这里∑ε为固定的已知噪声和干扰活动协方差矩阵。

步骤2：脑源活动分布及其区域活动分布的估计；

为了求解源活动x，首先定义后验概率分布p(x|y)为：

利用贝叶斯估计，后验概率分布的方差和均值表示为：

公式(12)的后验均值也可以写成：

上式中：

为模型数据的协方差矩阵。故公式(13)可以表示为：

故公式(13)的空间域滤波更新规则表示为：

公式(16)为源活动时间序列更新规则。

步骤3：贝叶斯超参数协方差分布估计方法；

具体步骤包括：

步骤3.1：根据上述公式(12)和(16)对脑源活动xk进行贝叶斯估计，然而只预先估计超参数γ的值，才能利用公式(12)或公式(16)进行估计。超参数γ可以通过最大化边缘分布p(y|γ)的估计值来获取。此边缘分布p(y|γ)为：

尽管超参数γ的值可以通过最大化边缘分布p(y|γ)(公式17)来获得，但事实上，很难对公式(17)的最大化求解，因为log|∑y|的存在。由于log|∑y|为超参数γ的凸函数，故使用dc×dc的辅助矩阵参数λj,(j＝1,…,n+r)得关系如下：

公式(18)始终满足，然后构造辅助代价函数为：

故始终有：

步骤3.2：通过优化超参数γ的值来获取辅助代价函数值的同时，此超参数υ也增加边缘分布logp(y|γ)的值。因此，超参数υ可以通过辅助函数进行估计：

公式(19)对超参数γ求导可得：

设置公式(22)右边等于0，有

由于求解辅助变量矩阵为半正定矩阵故有：

公式(24)即为超参数υ的更新规则。

由于超平面为凸函数log|∑y|的切平面，故辅助参数矩阵λj通过函数log|∑y|斜率来求取，故有辅助参数矩阵λj的更新规则为

对超参数γj的求取是通过循环迭代公式(16)公式(24)和公式(25)来实现，理论上来说，每次迭代过程都保证降低(或保持不变)辅助代价函数的值。每次迭代求取的计算为n+r的线性复杂度，故其此计算过程的计算代价是可以接受的。但是由此辅助函数求取的迭代过程收敛速度比基于em方法的速度快很多。

步骤4：基于多尺度源活动层析成像方法(即tree_champagne和tiling_champagne两种不同的源重构计算方法)，通过计算上述的迭代更新过程直至收敛，可以同时得到体元本征协方差和脑区协方差的信息；

具体步骤包括：

步骤4.1：由于要得到体元本征协方差脑区协方差信息，就要超参数之间的关系。这里用γ^v表示体元的本征协方差矩阵，γ^r表示分割区域区的协方差矩阵，其与求取超参数γ之间的关系：

上公式中

基于上述分析，利用体元本征分布协方差矩阵以及区域平均协方差矩阵提出两种不同的源重构计算方法。

步骤4.2：tree_champagne重构方法。假设第i个体元(待求解源活动位置)的协方差是由i个源的本征固有协方差与其所属于的脑区的协方差之和，故根据上述理论，tree_champagne求取的源的活动时间序列为：

式中，第i个体元的协方差分布为：

其中体元i被分割在第j个脑区中。

步骤4.3：tiling_champagne重构方法。假设第i个体元的分布的协方差矩阵为其所在脑区分布的平均协方差，故位于同一个脑区中的所有的体元有着同样的协方差矩阵，体元i的时间序列为：

这里第i个体元的协方差为

这里体元i被分割在第j个脑区中。

本发明提供的一种基于全新的多层次多尺度meg、eeg脑源重构通用模型，基于此模型的方法可以看作是champagne多层次多尺度扩展，其实现了对多尺度脑源活动的重构。面对更复杂的不同尺度脑源活动(如一簇多个相邻偶极子活动或者是一簇多个相邻偶极子与独立偶极子同时活动的情况)问题时，本发明首先将脑源空间分解为相互独立的体元，每个体元作为一个潜在的活动源。然后将所有的体元根据人脑解剖结构或者功能区分解为若干个区域，每个区域对应着不同的分布参数。另外，采用基于贝叶斯及其凸函数理论的方法对脑外采样数据进行分析，估计所有体元本征固有活动分布以及区域活动分布的协方差成分。

本发明是基于贝叶斯及其凸函数理论的方法对脑外采样数据进行分析，估计所有体元本征固有活动分布以及区域活动分布的协方差成分。根据上述模型框架提出两种方法tree_champagne和tiling_champagne是基于模型公式(1)的两个源重构方法，其计算过程如下：首先通过利用变分贝叶斯方法(vbfa)对刺激前的背景数据进行噪声协方差∑ε估计，设置超参数γ的初始值，利用公式(11)和公式(12)计算后验分布的均值和方差。然后利用公式(24)和公式(25)计算超参数γ以及辅助矩阵λ的值；接着利用公式(16)来对模型源活动时间序列进行更新，重复以上迭代过程直至收敛，得到贝叶斯后验分布超参数γ的值公式(24)；最后通过公式(28)以及公式(29)来计算tree_champagne的估计时间序列值以及协方差，公式(30)和(31)计算tiling_champagne的估计值。

具体实施时，本领域技术人员可采用计算机软件技术实现以上流程的自动运行。

应当理解的是，本说明书未详细阐述的部分均属于现有技术。

应当理解的是，上述针对较佳实施例的描述较为详细，并不能因此而认为是对本发明专利保护范围的限制，本领域的普通技术人员在本发明的启示下，在不脱离本发明权利要求所保护的范围情况下，还可以做出替换或变形，均落入本发明的保护范围之内，本发明的请求保护范围应以所附权利要求为准。