[发明专利]一种基于脑源空间分割的多相关源扫描成像方法

权利要求书

1.一种基于脑源空间分割的多相关源扫描成像方法，其特征在于，包括以下步骤：

步骤1：多层区域联合源重构概率模型的生成；

步骤1的具体实现包括以下子步骤：

步骤1.1：假设传感器在t时刻采集脑外MEG或EEG的数据表示为y(t)＝[y₁(t)，y₂(t)，…，y_M(t)]^T，其中y_i(t)表示传感器在t时刻采集的脑外数据，i＝1，2，...，M，M为采集点的通道数目；通过对脑活动区域进行离散分割，同时假设每个体元都为潜在脑源活动，故在t时刻的脑源活动为N为分割的潜在脑源活动数目，位于第i个元素在t时刻的源活动s_i(t)表示为d_c维度的向量这里d_c表示源活动方向个数；

则在对位置i的脑源活动扫描求解时，假设整个脑活动是由第i个元素的活动以及其它虚拟源相加，此模型为：

y(t)＝l_is_i(t)+g^\iv^\i(t)+ε (1)

这里为位置i的leadfield矩阵，l_i的第k列代表位置i的源在k方向为单位源活动时，脑外传感器采集其活动产生信号的强度；g^\i＝[g₁，g₂，...，g_Q-1]是其余Q-1个虚拟源活动对应的leadfield矩阵，这些虚拟源包含了除i源以外所有脑活动以及背景干扰和噪声，Q为虚拟源和真实源i活动的数目总和，v^\i(t)＝[v₁(t)，v₂(t)，…，v_Q-1(t)]^T是虚拟源对应的在时刻t的活动强度，其第j个虚拟源活动参数ε为相互独立且其分布服从正态分布的高斯噪声信号；

步骤1.2：根据脑解剖或者功能结构，将所有虚拟源无重合分割为R个区域，假设区域j包含p_j个虚拟源；则公式(1)表示为：

g_j，m和v_j，m(t)分别对应于位于区域j中编号为m的虚拟源，表示属于区域j中所有虚拟源产生的脑外采数据求和；如果对所有位置变量赋予概率分布，此模型即为一个概率模型；

步骤1.3：假设公式(2)的噪声ε为已知，假设源活动参数和噪声活动都为高斯分布，他们之间相互独立且在任意时刻都有着同样的随机概率分布；

假设源活动i的先验协方差分布为d_c×d_c的矩阵φ_i，下标i标明此次扫描是对位置i的源活动s_i进行求解，故有：

虚拟源v^\i(t)活动的先验分布假设为：

这里ψ_j为d_c×d_c的矩阵，其表示位于区域j虚拟源有着相同分布的协方差矩阵ψ_j；

如果第i个虚拟源位于区域j中，定义第i个虚拟源协方差矩阵为γ_i＝ψ_j，故得：

这里：

如果将位置i的脑源活动看做为仅有一个脑源活动i的区域，对公式(2)化简得通用模型如下：

其中：

公式(7)的模型中，只有第一个区域中的活动源为本次扫描求解的真实源活动，故f₁＝l_i和x₁＝s_i即为本次计算的关键，而其余源活动为伪源；

步骤1.4：假设y(t_k)表示为y(k)，x(t_k)表示为x(k)，整个时间窗口的脑源活动表示为X＝[x(t₁)，x(t₂)，…，x(t_K)]，脑外传感器采集数据表示为Y＝[y(t₁)，y(t₂)，…，y(t_K)]，假设其先验分布为：

其中先验协方差矩阵Ω为d_cQ×d_cQ的块对角矩阵：

这里，Ω_i为d_c×d_c的矩阵，当i＝1时，Ω_i＝φ₁对应于位置i的源活动先验协方差矩阵，当i＝2，3，…，Q时，Ω_i＝γ_i-1对应于第i-1虚拟源活动的先验协方差矩阵；由于已经假设噪声ε的分布已知，为故条件概率分布p(Y|X)表示成：

为了估计源X的分布，需要求解模型的后验概率分布p(X|Y)，定义为：

通过贝叶斯估计，结合公式(9)，公式(11)以及公式(12)得后验概率的精度和均值为：

后验均值表示为：

其中：

∑_y＝∑_ε+fΩf^T (16)

为模型数据的协方差矩阵，公式(15)的结果表示为：

公式(15)的j源活动的更新公式为：

公式(18)即为多层区域联合源重构概率模型的时间序列更新规则；

步骤2：多层区域联合源超参数分布的估计；

步骤2的具体实现包括以下子步骤：

步骤2.1：在贝叶斯估计中，超参数Ω的估计通过最大化边缘概率分布p(Y|Ω)来进行，边缘概率分布为：

步骤2.2：对公式(19)的代价函数进行修改，使用与Champagne方法相同的辅助函数来进行加速；由于log|∑_y|是超参数Ω的凸函数，故引入d_c×d_c辅助参数矩阵Λ_j，j＝1，2，…，Q，满足：

这里Λ₀为标量值，然后定义辅助代价函数为：

其始终满足：

故如果超参数Ω的值能使辅助代价函数的值增长，其超参数Ω也使边缘概率logp(Y|Ω)值增加，故问题转化为对辅助代价函数求最大值的优化，有：

超参数的更新规则通过辅助代价函数对其求导得：

设置公式(25)为0，得：

由于Ω的需要满足公式(23)，故得结果：

公式(27)即为超参数Ω的更新规则；

步骤2.3：求取辅助参数Λ_j的更新规则；

Λ_j的更新规则通过对log|∑_y|求超平面斜率获得：

这里Ω_j为超参数Ω的块对角矩阵元素，下面需要求解每个分割区域的平均协方差矩阵，因为假设每个区域内的所有元素有着相同的协方差矩阵，故对于区域j的协方差矩阵有：

然后将本次求解的公式(18)，公式(27)，公式(28)以及公式(29)的值作为下次迭代循环的初始值，继续迭代运算，知道辅助函数或边缘分布收敛；

步骤3：多层区域联合源虚拟源的确定；

步骤3的具体实现过程是：

在对位置i的源活动进行求取时，假设虚拟源及其leadfield矩阵与非位于i的源活动及其leadfield矩阵相同，即g^\i＝l^\i，这里l^\i＝[l₁，…，l_i-1，l_i+1，…，l_N]，N是源空间分割成元素的数目，同时有Q＝N；在此假设下，如果设置所有虚拟源都有着相同的分布以及相等的协方差矩阵，即对其不进行区域分割R＝1，所有的虚拟源都隶属于同一个区域，基于多层区域联合源的模型称为SOCHU_MN；假如将虚拟源分割区域的数量R＞1，隶属于不同区域的虚拟源有着不同的分布及协方差矩阵，同一区域的虚拟源有着同样的分布，基于此设定的多层区域联合源称为SOCHU_tiles；

步骤4：多层区域联合源余最小方差自适应波束形成方法联系；

步骤4的具体实现包括以下子步骤：

步骤4.1：最小方差自适应波束形成MVAB方法通过对位置i的偶极子源在时间为t时的活动进行重建，这里W_i为元素i的空间域滤波器定义为：

由公式(1)模型y(t)＝l_is_i(t)+g^\iv^\i(t)+ε得出，对于位置i的元素进行源活动估计时，其干扰与噪声之和为z_n＝g^\iv^\i(t)+ε，采样信号y(t)的协方差矩阵可以表示为∑_y＝∑_ε+fΩf^T；假设已经估计模型参数为∑_ε和φ_i，则多层区域联合源对位置i元素的源活动时间序列为其中：

这里φ_i为位置i源活动的先验分布值，公式(32)与公式(31)在经过变换后是相同；

步骤4.2：假设在无限数据采样点的情况下，采样数据的协方差矩阵满足使用矩阵求逆引理，R_yy的逆可以表示为因此有：

下面使用Θ的估计值为其假设位置i的源活动先验分布方差约等于0；故有：

最后一步使用公式(33)。