CN105147288A - 脑磁源强度定位方法 - Google Patents

Abstract

本发明公开了一种脑磁源强度定位方法，包括：通过引入时域平滑正则算子构造双参数混合正则化代价函数，然后根据广义交叉验证原则选取双正则化参数，通过单正则项在源信号矩阵中的比重对原始代价函数进行求解，得到脑磁源强度和位置确定的解矩阵，从而完成定位过程。本发明提出的脑磁源强度定位方法具有：总体均方误差小，且噪声越大时优势越明显；各时刻均方误差基本稳定在同一个水平，且受噪声影响小，根据本发明的方法能够重建得到时空准确且时域平滑的脑内神经信号，实现了脑磁源的精确定位。

Description

脑磁源强度定位方法

技术领域

本发明涉及生物信息技术领域，特别是一种脑磁源强度定位方法。

背景技术

根据头脑表观的磁场强度反演定位磁源的空间活动位置时脑磁研究中的一个重要问题，其本质上是一个非线性优化逆问题，为了简化计算的复杂性，在脑磁源的反演定位中，常用一线性方法去逼近非线性问题。现有技术中具体是采用脑磁源成像技术去进行脑磁源定位的，传统最小范数估计法(MNE，minimumnormestimate)，是最具代表性的源成像方法，主要原理是，基于大脑在特定时刻只有局部神经元活动的前提，对欠定线性方程增加l₂范数约束，求解一副能量最小的电流密度分布图像，常用的方式为：

假设脑外有m个通道的MEG信号，脑内有n个均匀分布的源信号，那么在i时刻，脑内源信号与MEG信号的关系可以用以下离散化的线性模型表示：

b_i＝Ax_i+e_i

其中，b_i为第i时刻大小为m×1的MEG测量信号；x_i为第i时刻脑内源信号，大小为n×1；e_i是第i时刻和_bi同维度的噪声信号；A为引导场矩阵，代表脑内源信号与MEG测量信号的映射关系，大小为m×n。当矩阵A已知时，即可由b_i求出脑磁逆问题的解x_i。但矩阵A是病态的，其条件数，即最大特征值与最小特征值之比很大，因此直接求逆不合适，通常转化为求解最小二次泛函的问题。不难得知，上式最小二乘解对噪声e_i非常敏感，MEG测量信号中很小的噪声将对解产生很大的扰动，造成无用解。引入Tikhonov正则化技术来减小噪声对对解的影响。在i时刻，脑磁逆问题求解转化为求解下式的最小值问题：

$f=\arg \min \left\{\right||{\mathrm{Ax}}_i-b_i|{|}_2^2+{\mathrm{\λ}}^2\left|\right|{\mathrm{Rx}}_i\left|{|}_2^2\right\}$

等式右边第一项表示测量数据和估计数据的拟合，第二项为正则项，表示解的先验信息，其中R为约束解空间的正则算子，当m个信号通道噪声均匀一致时，R取单位矩阵I，λ为正则化参数，调节拟合项和正则项在两项之间达到平衡。上式对应的解的形式为：

${\hat{x}}_i=A^T{({\mathrm{AA}}^T+{\mathrm{\λ}}^2{R}^{T}R)}^{-1}{b}_i$

由此可见，代价函数估算出来的源信号各个时刻之间是相互独立的。

但这种脑磁源方法的不足之处在于：1、对连续的MEG测量信号，没有全局考量噪声的影响，导致估算出来的源信号位置和强度的准确性差；

2、相邻时刻之间的估算结果在时域上有跳变，即时域上不平滑，不符合神经元定向传导的性质。

发明内容

针对上述技术问题，本发明中提出了一种脑磁源强度定位方法，该方法不同于传统最小范数估计算法(minimumnormestimate,MNE)，通过引入时域平滑正则算子构造双参数混合正则化代价函数，然后根据广义交叉验证(generalizedcross-validationcriterion，GCV)原则选取双正则化参数，求解代价函数，从而得出最优解，实现了脑磁源的精确定位。

为了实现根据本发明的这些目的和其它优点，提供了一种脑磁源强度定位方法，包括以下步骤：

步骤1)通过核磁共振仪对待测试人头部进行扫描，生成脑MR结构图像；通过脑磁仪在所述待检测人头部设置m个信号采集点，在时段a内连续采集k个时刻的信号，生成脑磁图MEG信号b，其中，MEG信号b以m×k维矩阵[b₁...b_i...b_k]表示，b_i为某一时刻的m×1维列向量，i为正整数，且1≤i≤k-1，将MR结构图像与MEG信号b的空间位置进行配准，得到真实几何头模型；

步骤2)将所述头模型上的大脑皮层设定为球模型，其上均匀分布有n个位置确定的等效磁偶极子，即n个脑磁源信号，进而确定引导场矩阵A；用源信号矩阵x表示在时段a内n个脑内源信号的强度随时间的变化，通过MEG信号b与源信号矩阵x如下关系式计算出源信号矩阵x，

b＝Ax+e1)

其中，x为n×k维矩阵[x₁...x_i...x_k]，x_i为某一时刻的n×1维源信号列向量，e为m×k维噪声信号矩阵，n远大于m；

步骤3)构造双参数正则化代价函数，使得求解所述代价函数得出的解矩阵x_opt在整个时段a内是所有解中全局能量最小，且在解矩阵x_opt中相邻两个源信号列向量的强度平滑变化，代价函数为：

$x_{opt}=\arg \min \left\{\underset{i=1}{\overset{k}{\Σ}}\right||{\mathrm{Ax}}_i-b_i|{|}_2^2+{{\mathrm{\λ}}_1}^2\underset{i=1}{\overset{k}{\Σ}}\left|\right|x_i|{|}_2^2+{{\mathrm{\λ}}_2}^2\underset{i=1}{\overset{k-1}{\Σ}}||x_{i+1}-x_i|{|}_2^2\}---2)$

其中，x_opt为n×k维矩阵，等式右边第一项表示测量数据和估计数据的拟合项，第二项为时域平滑约束项，采用广义交叉验证方法自动选取正则化参数λ₁和λ₂；

步骤4)将求解式1)的脑磁逆问题转化为求解式2)的最小值问题，求解式2)得到x_opt，x_opt中每一个元素表示某一时刻对应该位置处的源信号强度，将x_opt中每一列向量中n个元素匹配到大脑皮层n个精确位置上，即完成了任意时刻每个位置上脑磁源信号强度的定位。

优选的，所述步骤1)中，脑磁仪采集到的数据经过去眼电、滤波以及基线校准后得到所述MEG信号b。

优选的，采用边界元或有限元方法结合所述头模型求解正问题获取m×n维的所述引导场矩阵A，其反映MEG信号b与源信号矩阵x映射关系。

优选的，所述脑磁仪设置有148个信号采集通道。

优选的，所述步骤3)中，根据在a时段内第_i+1时刻的源信号解向量x_i+1满足如下条件：x_i+1＝x_i+Δη_i，Δη_i→0，来构造所述时域平滑约束项，联合正则化参数λ₂将所述时域平滑约束项构造成 ${{\mathrm{\λ}}_2}^2\underset{i=1}{\overset{k-1}{\Σ}}\left|\right|x_{i+1}-x_i|{|}_2^2.$

优选的，所述步骤3)中，具体地，构造以下方程：

$GCV=\frac{km\left|\right|\stackrel{\‾}{A}\widehat{\stackrel{\‾}{x}}-\stackrel{\‾}{b}|{|}_2^2}{{\[trace(I_{km}-f_{{\mathrm{\λ}}_1{\mathrm{\λ}}_2})\]}^2}---3)$

其中，k与m都为正整数，I_km为单位矩阵， $f_{{\mathrm{\λ}}_1{\mathrm{\λ}}_2}=\stackrel{\‾}{A}{\stackrel{\‾}{A}}^T{(\stackrel{\‾}{A}{\stackrel{\‾}{A}}^T+{{\mathrm{\λ}}_1}^{2}I+{{\mathrm{\λ}}_2}^2{L}^{T}L)}^{-1}{\stackrel{\‾}{A}}^T,$ 通过遗传算法求式3)的最小值来确定λ₁和λ₂，I_k和I为单位矩阵。

优选的，还包括步骤5)，通过确定对应的单正则项的最优解和对应的单正则项的最优解在源信号矩阵中所占的比重求解式2)，从而得到x_opt。

优选的，所述步骤5)中，先引用Kronecker积将式2)转化成如下形式：

$\begin{array}{c}{x}_{opt}=\arg \min \left\{\right||\stackrel{\‾}{A}\stackrel{\‾}{x}-\stackrel{\‾}{b}|{|}_2^2+{{\mathrm{\λ}}_1}^2\left|\right|\stackrel{\‾}{x}\left|\right|+{{\mathrm{\λ}}_2}^2\left|\right|L\stackrel{\‾}{x}\left|{|}_2^2\right\}\\ =\arg \min \left\{\frac{1}{2}\underset{j=1}{\overset{2}{\Σ}}\left(\right||\stackrel{\‾}{A}\stackrel{\‾}{x}-\stackrel{\‾}{b}|{|}_2^2+2{{\mathrm{\λ}}_j}^2\left|\right|Q_j\stackrel{\‾}{x}\left|{|}_2^2\right)\right\}\end{array}---4)$

其中：

$\stackrel{\‾}{A}=\left[\begin{array}{cccc}A& 0& \mathrm{...}& 0\\ 0& A& \mathrm{...}& 0\\ .& .& & .\\ .& .& & .\\ .& .& & .\\ 0& 0& \mathrm{...}& A\end{array}\right]$

$\stackrel{\‾}{x}={\[x^T\left(1\right),x^T\left(2\right),\mathrm{...},x^T\left(k\right)\]}^T$

$\stackrel{\‾}{b}={\[b^T\left(1\right),b^T\left(2\right),\mathrm{...},b^T\left(k\right)\]}^T$

Q₁＝I，Q₂＝L

求解λ_j对应的单正则项下的解： ${\stackrel{\‾}{x}}_{\mathrm{\λ}j}={\[{x_{\mathrm{\λ}j}}^T\left(1\right),{x_{\mathrm{\λ}j}}^T\left(2\right),\mathrm{...},{x_{\mathrm{\λ}j}}^T\left(k\right)\]}^T,$ j＝1,2

最后由以下方程得出所述解矩阵x_opt：

$x_{opt}=a{\stackrel{\‾}{x}}_{\mathrm{\λ}1}+(1-a){\stackrel{\‾}{x}}_{\mathrm{\λ}2}$

其中 $a=\frac{-{c_2}^T(c_1-c_2)}{\left|\right|c_1-c_2|{|}_2^2},c_j=(2{{\mathrm{\λ}}_j}^2{Q}_j-\underset{r=1}{\overset{2}{\Σ}}{{\mathrm{\λ}}_r}^2{Q}_r){\stackrel{\‾}{x}}_{\mathrm{\λ}j},$ j＝1,2，r＝1,2。

本发明至少包括以下有益效果：

1、实现了时序脑磁信号源估算时正则参数的自动定位，脑磁源定位更快更精确；

2、减弱了噪声对源估算的影响；

3、减弱了各个时刻脑磁源之间的跳变，保证了脑磁源在整个时域上的平滑性；

4、重建的脑磁源强度随时间的变化过程更加逼近真实的神经元传导性质。

本发明的其它优点、目标和特征将部分通过下面的说明体现，部分还将通过对本发明的研究和实践而为本领域的技术人员所理解。

附图说明

图1为本发明的脑磁源强度定位方法的流程示意图；

图2(a)为一实施例中两个仿真源信号的信号发生位置示意图；

图2(b)为一实施例中仿真MEG测量信号的示意图(信噪比为6dB时)；

图3(a)为6ms时无噪声仿真信号在脑皮层上脑磁源信号强度的成像图；

图3(b)为19ms时无噪声仿真信号在脑皮层上脑磁源信号强度的成像图；

图4(a)传统MNE方法计算出各时刻原始仿真信号与估算信号的均方误差示意图；

图4(b)时域平滑约束条件下计算出各时刻原始仿真信号与估算信号的均方误差示意图；

图4(c)传统MNE方法与本发明的定位方法的总体均方误差比较示意图；

图5(a)左脑仿真信号与采用传统MNE方法得出的估算信号之间的吻合度对比示意图；

图5(b)右脑仿真信号与采用传统MNE方法得出的估算信号之间的吻合度对比示意图；

图6(a)左脑仿真信号与采用本发明的时域平滑约束下得出的估算信号之间的吻合度对比示意图；

图6(b)右脑仿真信号与采用本发明的时域平滑约束下得出的估算信号之间的吻合度对比示意图。

具体实施方式

下面结合附图对本发明做进一步的详细说明，以令本领域技术人员参照说明书文字能够据以实施。

应当理解，本发明所使用的诸如“具有”、“包含”以及“包括”术语并不配出一个或多个其它元件或其组合的存在或添加。

本发明的脑磁源信号强弱的定位方法不用在疾病诊断过程中。

如图1所示，本发明提供了一种脑磁源强度定位方法，具体包括以下步骤：

步骤1)通过核磁共振仪对待测试人头部进行扫描，生成脑MR结构图像；通过脑磁仪在所述待检测人头部设置m个信号采集点，在时段a内连续采集k个时刻的信号，包括第1时刻到第k时刻，生成脑磁图MEG信号b，其中，MEG信号b以m×k维矩阵[b₁...b_i...b_k]表示，b_i为某一时刻的m×1维列向量，i为正整数，且1≤i≤k-1，将MR结构图像与MEG信号b的空间位置进行配准，得到真实几何头模型；

步骤2)将所述头模型上的大脑皮层设定为球模型，其上均匀分布有n个位置确定的等效磁偶极子，即n个脑磁源信号，每个脑磁源信号大小即磁偶极子强度，脑磁源方向设定为与脑皮层表面垂直的方向，结合所述头模型，确定引导场矩阵A；用源信号矩阵x表示在时段a内n个脑内源信号的强度随时间的变化，通过MEG信号b与源信号矩阵x如下关系式计算出源信号矩阵x，

b＝Ax+e1)

其中，x为n×k维矩阵[x₁...x_i...x_k]，其中，x_i为某一时刻的n×1维源信号列向量，即某一时刻式1)的解，e为m×k维噪声信号矩阵，为已知量；当矩阵A已知时，即可由b求出脑磁逆问题的解x，也就是说求解式1)的逆问题即可得出所述源信号矩阵x，一般n远大于m，由此根据式1)求解得出的源信号矩阵x不唯一，且矩阵A是病态的，其条件数，即最大特征值与最小特征值之比很大，因此直接求逆不合适，通常转化为求解最小二次泛函的问题，因此进入步骤3)；

步骤3)构造双参数正则化代价函数，使得求解所述代价函数得出的解矩阵x_opt在整个时段a内是所有解中能量最小，从而确定唯一解，且在解矩阵x_opt中相邻两个源信号列向量的强度平滑变化，消除相邻两个解向量之间的过渡依然具有跳变性，避免解向量在时域上不平滑和出现不规则震荡，使得源信号矩阵信号强度随时间的变化符合神经元定向传导的性质，代价函数为：

$x_{opt}=\mathrm{argmin}\left\{\underset{i=1}{\overset{k}{\Σ}}\right||{\mathrm{Ax}}_i-b_i|{|}_2^2+{\mathrm{\λ}}^2\underset{i=1}{\overset{k}{\Σ}}\left|\right|x_i|{|}_2^2+{\mathrm{\μ}}^2\underset{i=1}{\overset{k-1}{\Σ}}||x_{i+1}-x_i|{|}_2^2\}---2)$

其中，x_opt为n×k维矩阵，等式右边第一项表示测量数据和估计数据的拟合项，第二项为时域平滑约束项，采用广义交叉验证方法自动选取正则化参数λ₁和λ₂；

步骤4)将求解式1)的脑磁逆问题转化为求解式2)的最小值问题，求解式2)得到x_opt，x_opt中每一个元素表示某一时刻对应该位置处的源信号强度，将x_opt中每一列向量中n个元素匹配到大脑皮层n个精确位置上，从而可以的在时段a上对n个位置确定的脑磁源的信号强弱进行准确分析，n个精确位置上的脑磁源信号强度随时间的变化就可以重建获知，也就是完成了脑磁源信号强弱在脑模型上的快速准确定位，任意时刻每个位置上脑磁源信号的强度即可获知，从而可以对脑磁源信号进行重建并成像，得到源图像，重建的脑磁源强度随时间的变化过程减弱了各个时刻脑磁源之间的跳变，保证了脑磁源在整个时域上的平滑性，从而更加逼近真实的神经元传导性质。

上述技术方案中，所述步骤1)中所述脑磁仪设置有m个采集通道，采集到的数据在脑磁仪中经过去眼电、滤波和基线校准等预处理步骤后，即可得到所需的MEG信号。

上述技术方案中，所述步骤2)中，采用边界元或有限元等方法结合所述头模型求解正问题获取m×n维的所述引导场矩阵A，其反映MEG信号b与源信号矩阵x映射关系，本实施例中采用的边界元方法。

上述技术方案中，所述步骤3)中，根据在a时段内第_i+1时刻的源信号解向量x_i+1满足如下条件：x_i+1＝x_i+Δη_i，Δη_i→0，来构造所述时域平滑约束项，联合正则化参数λ₂将所述时域平滑约束项构造成引入了时域平滑约束项，使得减弱了噪声对源信号矩阵x_opt的影响，从而减弱了各个时刻脑磁源之间的跳变，保证了脑磁源在整个时域上的平滑性。

上述技术方案中，所述步骤3)中，采用广义交叉验证GCV方法确定式2)中的两个正则化参数λ₁和λ₂，具体地，构造以下方程：

$GCV=\frac{km\left|\right|\stackrel{\‾}{A}\widehat{\stackrel{\‾}{x}}-\stackrel{\‾}{b}|{|}_2^2}{{\[trace(I_{km}-f_{{\mathrm{\λ}}_1{\mathrm{\λ}}_2})\]}^2}---3)$

其中，k与m都为正整数，I_km为单位矩阵， $f_{{\mathrm{\λ}}_1{\mathrm{\λ}}_2}=\stackrel{\‾}{A}{\stackrel{\‾}{A}}^T{(\stackrel{\‾}{A}{\stackrel{\‾}{A}}^T+{{\mathrm{\λ}}_1}^{2}I+{{\mathrm{\λ}}_2}^2{L}^{T}L)}^{-1}{\stackrel{\‾}{A}}^T,$ 通过遗传算法求式3)的最小值来确定λ₁和λ₂，I_k和I为单位矩阵，这种自动方法实现了正则参数一次性的自动定位，使得脑磁源定位更快更精确。

上述技术方案中，还包括步骤5)，通过确定对应的单正则项的最优解和对应的单正则项的最优解在源信号矩阵中所占的比重求解式2)，从而得到x_opt。

具体求解方法为：

先引用Kronecker积将式2)转化成如下形式：

$\begin{array}{c}{x}_{opt}=\arg \min \left\{\right||\stackrel{\‾}{A}\stackrel{\‾}{x}-\stackrel{\‾}{b}|{|}_2^2+{{\mathrm{\λ}}_1}^2\left|\right|\stackrel{\‾}{x}\left|\right|+{{\mathrm{\λ}}_2}^2\left|\right|L\stackrel{\‾}{x}\left|{|}_2^2\right\}\\ =\arg \min \left\{\frac{1}{2}\underset{j=1}{\overset{2}{\Σ}}\left(\right||\stackrel{\‾}{A}\stackrel{\‾}{x}-\stackrel{\‾}{b}|{|}_2^2+2{{\mathrm{\λ}}_j}^2\left|\right|Q_j\stackrel{\‾}{x}\left|{|}_2^2\right)\right\}\end{array}---4)$

其中：

$\stackrel{\‾}{A}=\left[\begin{array}{cccc}A& 0& \mathrm{...}& 0\\ 0& A& \mathrm{...}& 0\\ .& .& & .\\ .& .& & .\\ .& .& & .\\ 0& 0& \mathrm{...}& A\end{array}\right]$

$\stackrel{\‾}{x}={\[x^T\left(1\right),x^T\left(2\right),\mathrm{...},x^T\left(k\right)\]}^T$

$\stackrel{\‾}{b}={\[b^T\left(1\right),b^T\left(2\right),\mathrm{...},b^T\left(k\right)\]}^T$

Q₁＝I，Q₂＝L

求解λ_j对应的单正则项下的解： ${\stackrel{\‾}{x}}_{\mathrm{\λ}j}={\[{x_{\mathrm{\λ}j}}^T\left(1\right),{x_{\mathrm{\λ}j}}^T\left(2\right),\mathrm{...},{x_{\mathrm{\λ}j}}^T\left(k\right)\]}^T,$ j＝1,2

最后由以下方程得出所述解矩阵x_opt：

$x_{opt}=a{\stackrel{\‾}{x}}_{\mathrm{\λ}1}+(1-a){\stackrel{\‾}{x}}_{\mathrm{\λ}2}$

其中 $a=\frac{-{c_2}^T(c_1-c_2)}{\left|\right|c_1-c_2\left|\right|},c_j=(2{{\mathrm{\λ}}_j}^2{Q}_j-\underset{r=1}{\overset{2}{\Σ}}{{\mathrm{\λ}}_r}^2{Q}_r){\stackrel{\‾}{x}}_{\mathrm{\λ}j},$ j＝1,2，r＝1,2。从而得出所述源信号解矩阵x_opt，求解过程更加简化方便，且源信号解矩阵x_opt更加逼近真实的神经元传导性质，便于认知神经学、脑神经疾病的研究。

以下结合几个对比例来进一步说明本发明的脑磁源强度定位方法优于基于传统MNE算法的脑磁源强度定位方法。

为了利用分析和说明，采用电脑模拟的方式来更加清楚、准确地说明本发明相对于传统定位方法的有益效果，具体对比方式如下：

对比例1

模拟148通道的脑磁图仪器，设定真实几何头模型内部共有7850个均匀分布的格点，代表7850个源的位置，由正弦指数函数生成模拟两个仿真源信号，分辨率为1000Hz，时长为40ms，且分别在第6ms和19ms处达到能量峰值，如图2(a)所示，为模拟真实脑磁信号，生成MEG仿真数据时，根据信噪比定义添加方差为σ²高斯白噪声，分别得到4dB～12dB的MEG仿真数据，图2(b)所示为仿真的148通道6dB测量信号波形图。

在脑皮层选取两个活化位置，坐标分别为(-39.4982，-36.6656，56.8917)和(36.0071，-18.8000，58.9000)，两个位置分别对应左脑和右脑感觉区。6ms时达到能量峰值的源信号被放置在(-39.4982，-36.6656，56.8917)处，19ms时达到能量峰值的源信号则被放置在(36.0071，-18.8000，58.9000)处，如3(a)和(b)所示。

本实施例中，在用于脑电/脑磁信号分析的开源软件eConnectome平台上完成了上述仿真数据的设计，并在此基础上对发明提出的方法进行了实验验证。具体地，我们将MEG仿真数据导入eConnectome后执行数据预处理(Preprocessing)，包括BaselineCorrection(以1～4ms为基准线)和Filtering(50Hz陷波滤波器)，采用真实几何头模型和边界元方法求解正问题获取转换矩阵A，然后分别用传统MNE方法和式2)基于时域平滑约束方法对经过预处理的数据进行脑磁源重建，对比两者的实验结果。

实验结果考察两个方面：一是考察数据精确度参数均方误差，二是考察两个活化位置的估算信号与原始模拟信号的吻合情况。

采用均方误差MSE来评价重建方法的精确度：

$MSE=\left|\right|x-\hat{x}|{|}_2^2/m$

其中m是MEG信号通道数，本发明m取值148。x和分别为原始仿真信号和求逆得到的估算信号。

分别计算了信噪比为4dB、6dB、8dB、10dB和12dB时对应的各时刻均方误差和总体均方误差，如图4所示。由图4(a)可以看出，采用传统MNE方法的各时刻均方误差波动较大，而且噪声越大时，均方误差波动越大；图4(b)显示本发明提出的基于时域平滑约束算法各时刻均方误差基本稳定在一个水平，而且受噪声影响小；图4(c)显示基于时域平滑约束算法的总体均方误差远小于传统MNE方法，且噪声越大，优势越明显。

对比例2

图5(a)和图5(b)分别展示了脑皮层上两个活化位置(-39.4982，-36.6656，56.8917)和(36.0071，-18.8000，58.9000)采用传统MNE方法的估算信号和仿真信号之间的吻合情况。发现各个时刻间独立求逆使得解在时域上不规则振荡，且某些时刻与真实值相去甚远。图6(a)和图6(b)显示引入双参数正则化增加时域平滑约束项后，估算信号基本复原了仿真信号变化趋势，而且分别在6ms和19ms处具有能量峰值，解决了信号在相邻时刻的跳变问题。需要注意的是，估算信号的幅度小于真实信号，是因为式4)中第二项是能量约束项，也就是说所求的估算信号是所有解中能量最小的解，这是重建算法本身决定的，MNE算法也存在同样的现象。

由上述几个对比例的实验结果可知，本发明提出的基于时域平滑约束的双参数MEG时序信号逆问题求解方法而进行的脑磁源强度定位方法要优于传统MNE方法。

可以理解的是，本发明的脑磁源强度定位方法也可用于根据脑电信号来进行源定位。

尽管本发明的实施方案已公开如上，但其并不仅仅限于说明书和实施方式中所列运用，它完全可以被适用于各种适合本发明的领域，对于熟悉本领域的人员而言，可容易地实现另外的修改，因此在不背离权利要求及等同范围所限定的一般概念下，本发明并不限于特定的细节和这里示出与描述的图例。

Claims (8)

1.一种脑磁源强度定位方法，其特征在于，包括以下步骤：

步骤1)通过核磁共振仪对待测试人头部进行扫描，生成脑MR结构图像；通过脑磁仪在所述待检测人头部设置m个信号采集点，在时段a内连续采集k个时刻的信号，生成脑磁图MEG信号b，其中，MEG信号b以m×k维矩阵[b₁...b_i...b_k]表示，b_i为某一时刻的m×1维列向量，i为正整数，且1≤i≤k-1，将MR结构图像与MEG信号b的空间位置进行配准，得到真实几何头模型；

步骤2)将所述头模型上的大脑皮层设定为球模型，其上均匀分布有n个位置确定的等效磁偶极子，即n个脑磁源信号，进而确定引导场矩阵A；用源信号矩阵x表示在时段a内n个脑内源信号的强度随时间的变化，通过MEG信号b与源信号矩阵x如下关系式计算出源信号矩阵x，

b＝Ax+e1)

其中，x为n×k维矩阵[x₁...x_i...x_k]，x_i为某一时刻的n×1维源信号列向量，e为m×k维噪声信号矩阵，n远大于m；

步骤3)构造双参数正则化代价函数，使得求解所述代价函数得出的解矩阵x_opt在整个时段a内是所有解中全局能量最小，且在解矩阵x_opt中相邻两个源信号列向量的强度平滑变化，代价函数为：

$x_{opt}=\arg \min \left\{\underset{i=1}{\overset{k}{\Σ}}\left|\right|{\mathrm{Ax}}_i-b_i|{|}_2^2+{{\mathrm{\λ}}_1}^2\underset{i=1}{\overset{k}{\Σ}}\left|\right|x_i|{|}_2^2+{{\mathrm{\λ}}_2}^2\underset{i=1}{\overset{k-1}{\Σ}}\left|\right|x_{i+1}-x_i|{|}_2^2\right\}---2)$

其中，x_opt为n×k维矩阵，等式右边第一项表示测量数据和估计数据的拟合项，第二项为时域平滑约束项，采用广义交叉验证方法自动选取正则化参数λ₁和λ₂；

步骤4)将求解式1)的脑磁逆问题转化为求解式2)的最小值问题，求解式2)得到解矩阵x_opt，x_opt中每一个元素表示某一时刻对应该位置处的源信号强度，将x_opt中每一列向量中n个元素匹配到大脑皮层n个精确位置上，即完成了任意时刻每个位置上脑磁源信号强度的定位。

2.如权利要求1所述的脑磁源强度定位方法，其特征在于，所述步骤1)中，脑磁仪采集到的数据经过去眼电、滤波以及基线校准后得到所述MEG信号b。

3.如权利要求2所述的脑磁源强度定位方法，其特征在于，所述步骤2)中，采用边界元或有限元方法结合所述头模型求解正问题获取m×n维的所述引导场矩阵A，其反映MEG信号b与源信号矩阵x映射关系。

4.如权利要求1所述的脑磁源强度定位方法，其特征在于，所述脑磁仪设置有148个信号采集通道。

5.如权利要求1所述的脑磁源强度定位方法，其特征在于，所述步骤3)中，根据在a时段内第i+1时刻的源信号解向量x_i+1满足如下条件：x_i+1＝x_i+Δη_i，Δη_i→0，来构造所述时域平滑约束项，联合正则化参数λ₂将所述时域平滑约束项构造成 ${{\mathrm{\λ}}_2}^2\underset{i=1}{\overset{k-1}{\Σ}}\left|\right|x_{i+1}-x_i|{|}_2^2.$

6.如权利要求5所述的脑磁源强度定位方法，其特征在于，所述步骤3)中，具体地，构造以下方程：

$GCV=\frac{km\left|\right|\stackrel{\‾}{A}\widehat{\stackrel{\‾}{x}}-\stackrel{\‾}{b}|{|}_2^2}{{\[trace(I_{km}-f_{{\mathrm{\λ}}_1{\mathrm{\λ}}_2})\]}^2}---3)$

其中，k与m都为正整数，I_km为单位矩阵， $f_{{\mathrm{\λ}}_1{\mathrm{\λ}}_2}=\stackrel{\‾}{A}{\stackrel{\‾}{A}}^T{(\stackrel{\‾}{A}{\stackrel{\‾}{A}}^T+{{\mathrm{\λ}}_1}^{2}I+{{\mathrm{\λ}}_2}^2{L}^{T}L)}^{-1}{\stackrel{\‾}{A}}^T,$ 通过遗传算法求式3)的最小值来确定λ₁和λ₂，I_k和I为单位矩阵。

7.如权利要求6所述的脑磁源强度定位方法，其特征在于，还包括步骤5)，通过确定对应的单正则项的最优解和对应的单正则项的最优解在源信号矩阵中所占的比重求解式2)，从而得到x_opt。

8.如权利要求7所述的基于时域平滑约束的脑磁源定位方法，其特征在于，所述步骤5)中，先引用Kronecker积将式2)转化成如下形式：

$\begin{array}{c}{x}_{opt}=\arg \min \left\{\left|\right|\stackrel{\‾}{A}\stackrel{\‾}{x}-\stackrel{\‾}{b}|{|}_2^2+{{\mathrm{\λ}}_1}^2|\left|\stackrel{\‾}{x}\right|{|}_2^2+{{\mathrm{\λ}}_2}^2\left|\right|L\stackrel{\‾}{x}|{|}_2^2\right\}\\ =\arg \min \left\{\frac{1}{2}\underset{j=1}{\overset{2}{\Σ}}\left(\left|\right|\stackrel{\‾}{A}\stackrel{\‾}{x}-\stackrel{\‾}{b}|{|}_2^2+2{{\mathrm{\λ}}_j}^2|\left|Q_j\stackrel{\‾}{x}\right|{|}_2^2\right)\right\}\end{array}---4)$

其中：

$\stackrel{\‾}{A}=\left[\begin{array}{cccc}A& 0& ...& 0\\ 0& A& ...& 0\\ \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}.\\ .\\ .\end{array}& \begin{array}{c}\\ \\ \end{array}& \begin{array}{c}.\\ .\\ .\end{array}\\ 0& 0& ...& A\end{array}\right]$

$\stackrel{\‾}{x}={\[x^T\left(1\right),x^T\left(2\right),...,x^T\left(k\right)\]}^T$

$\stackrel{\‾}{b}={\[b^T\left(1\right),b^T\left(2\right),...,b^T\left(k\right)\]}^T$

Q₁＝I，Q₂＝L

求解λ_j对应的单正则项下的解： ${\stackrel{\‾}{x}}_{\mathrm{\λ}j}={\[{x_{\mathrm{\λ}j}}^T\left(1\right),{x_{\mathrm{\λ}j}}^T\left(2\right),...,{x_{\mathrm{\λ}j}}^T\left(k\right)\]}^T,$ j＝1,2，最后由以下方程得出所述解矩阵x_opt：

$x_{opt}=a{\stackrel{\‾}{x}}_{\mathrm{\λ}1}+(1-a){\stackrel{\‾}{x}}_{\mathrm{\λ}2}$

其中 $a=\frac{-{c_2}^T(c_1-c_2)}{\left|\right|c_1-c_2|{|}_2^2},c_j=(2{{\mathrm{\λ}}_j}^2{Q}_j-\underset{r=1}{\overset{2}{\Σ}}{{\mathrm{\λ}}_r}^2{Q}_r){\stackrel{\‾}{x}}_{\mathrm{\λ}j},j=1,2,r=1,2.$