MEG源定位方法及系统

技术领域

本发明涉及一种MEG源定位方法及系统。

背景技术

目前，商用的脑磁图仪(magnetoencephalography，MEG)主要基于超导量子干涉仪(SQUID)原理，SQUID传感器具有极强的弱磁探测能力，其探测灵敏度可达1fT√Hz。但存在几个缺陷：一是需液氦冷却，由于传感器阵列置于液氦杜瓦内部，导致多通道MEG体积较为庞大笨重；二是杜瓦的真空隔离层(一般大于20mm)使得传感器探测单元无法密切贴合所有患者的头皮，致使SQUID传感器的探测能力降低(磁场强度与信号源、传感器的距离平方成反比)。

基于原子磁力计的新型脑磁图扫描仪(MEG)，探测器单元更贴近头皮，探测信号具有更高的信噪比；此外，结合柔性头盔设计，实现了被试者在自然状态下的脑磁信号探测。但是，受限于高昂的原子磁力计成本和有限的头表皮空间，如何利用极少数量探测器单元进行脑磁源信号定位，对于基于原子磁力计的MEG系统设计极其重要。

现有的源定位处理流程主要基于传统的商用脑磁图仪(MEG)，源定位方法主要包含：LCMV算法、迭代类定位算法。传统的LCMV算法实时性较好，但无法满足高精度定位需求；传统的迭代类定位算法定位精度高，但无法满足MEG系统的实时性定位需求。

发明内容

有鉴于此，有必要提供一种MEG源定位方法及系统，其能够解决极少数探测器条件下脑磁源精确定位问题，提高MEG系统的空间分辨率，满足MEG系统的实时性定位需求。

本发明提供一种MEG源定位方法，该方法包括如下步骤：a.利用柔性头盔和由原子磁力计组成的探测器单元，按任务态对被试者进行单任务多区域的数据采集，并将采集的多通道数据按刺激器的初始位置对齐；b.将采集的上述多通道数据，按采集区域分组进行归一化处理；c.利用参考探测器通道数据，通过数值模拟方法计算实时动态的噪声系数，并对归一化处理后的数据进行降噪处理；d.根据降噪处理后的数据获取多通道数据的频谱信息，并从中提取单频段时域信号；e.利用被试者的T1磁共振图像计算得到其个性化T1-MRI模板，将MEG原子磁力计与T1-MRI模板坐标系配对，并计算得到被试者正向传递函数；f.利用步骤d提取的多个所述单频段时域信号和步骤e计算得到的所述被试者正向传递函数，估算得到多个单频段时域数据的源空间信号分布；g.将步骤f估算得到的上述源空间信号分布分解成信号部分和噪声部分，并计算得到分解矩阵和混合矩阵；h.利用步骤g计算得到的所述分解矩阵和混合矩阵，对步骤g中的信号部分进行分解，得到其空域主成分和时域主成分，并结合步骤f得到的单频段时域数据的源空间信号分布，准确定位被试者的脑磁源信号分布。

进一步地，所述的步骤c具体包括：

所述计算实时动态的噪声系数为通过合成梯度法获取，公式如下：

s＝σ-ξ·b

其中，s为降噪后结果输出，σ为测量探测器的输出，b

静息态时，测量探测器与参考探测器输出结果应该一致；利用最小二乘法，得到使s最小的情况下的实时动态噪声系数为：

ξ＝(b

任务态脑磁信号去噪后的结果则为：

out＝σ

其中，σ

进一步地，所述的步骤e具体包括：加载被试者T1磁共振扫描图像；

计算个性化被试者T1-MRI模板图像：利用的被试者的T1磁共振图像计算得到其个性化模板T1-MRI模板；

MEG原子磁力计与T1-MRI模板坐标系配对：将MEG原子磁力计的坐标信息与T1-MRI模板的基准点配对，并计算得到其在个性化模板中的坐标信息；

计算被试者正向传递函数：利用个性化的T1-MRI模板和MEG原子磁力计的坐标信息，计算得到被试者正向传递函数。

进一步地，所述的步骤f具体包括：

利用步骤d提取的多个所述单频段时域信号和步骤e计算得到的所述被试者正向传递函数，通过向量化的Beamformer方法估算得到多个单频段时域数据的源空间信号分布；其中，给定任意空间位置r，估算得到的源空间信号分布为：

s(t,r)＝W

其中，s(t,r)为t时刻、空间位置r估算得到的源空间信号分布；b(t)为t时刻向量化的多通道数据；W(r)为向量化空间滤波系数，其向量化的表达式如下：

W(r)＝(L

上式中，L为步骤e中计算得到的被试者正向传递函数，C为向量化的多通道数据的协方差矩阵，T为矩阵的转置操作符。

进一步地，所述的步骤g具体包括：

通过矩阵的奇异值分解，将步骤f估算得到的源空间信号分布分解成信号部分和噪声部分，并通过独立成分分析计算分解矩阵和混合矩阵；

其中，所述的奇异值分解公式定义如下：

对于时域子空间主成分

本发明提供一种MEG源定位系统，该系统包括：数据采集模块、归一化处理模块、降噪处理模块、信号提取模块、传递函数计算模块、源空间信号分布估算模块、矩阵计算模块、定位模块，其中：所述数据采集模块用于利用柔性头盔和由原子磁力计组成的探测器单元，按任务态对被试者进行单任务多区域的数据采集，并将采集的多通道数据按刺激器的初始位置对齐；所述归一化处理模块用于将采集的上述多通道数据，按采集区域分组进行归一化处理；所述降噪处理模块用于利用参考探测器通道数据，通过数值模拟方法计算实时动态的噪声系数，并对归一化处理后的数据进行降噪处理；所述信号提取模块用于根据降噪处理后的数据获取多通道数据的频谱信息，并从中提取单频段时域信号；所述传递函数计算模块用于利用被试者的T1磁共振图像计算得到其个性化T1-MRI模板，将MEG原子磁力计与T1-MRI模板坐标系配对，并计算得到被试者正向传递函数；所述源空间信号分布估算模块用于利用所述信号提取模块提取的多个所述单频段时域信号和所述传递函数计算模块计算得到的所述被试者正向传递函数，估算得到多个单频段时域数据的源空间信号分布；所述矩阵计算模块用于将所述源空间信号分布估算模块估算得到的上述源空间信号分布分解成信号部分和噪声部分，并计算得到分解矩阵和混合矩阵；所述定位模块用于利用所述矩阵计算模块计算得到的所述分解矩阵和混合矩阵，对所述矩阵计算模块中的信号部分进行分解，得到其空域主成分和时域主成分，并结合所述源空间信号分布估算模块得到的单频段时域数据的源空间信号分布，准确定位被试者的脑磁源信号分布。

进一步地，所述的降噪处理模块具体用于：

所述计算实时动态的噪声系数为通过合成梯度法获取，公式如下：

s＝σ-ξ·b

其中，s为降噪后结果输出，σ为测量探测器的输出，b

静息态时，测量探测器与参考探测器输出结果应该一致；利用最小二乘法，得到使s最小的情况下的实时动态噪声系数为：

ξ＝(b

任务态脑磁信号去噪后的结果则为：

out＝σ

其中，σ

进一步地，所述的传递函数计算模块具体用于：

加载被试者T1磁共振扫描图像；

计算个性化被试者T1-MRI模板图像：利用的被试者的T1磁共振图像计算得到其个性化模板T1-MRI模板；

MEG原子磁力计与T1-MRI模板坐标系配对：将MEG原子磁力计的坐标信息与T1-MRI模板的基准点配对，并计算得到其在个性化模板中的坐标信息；

计算被试者正向传递函数：利用个性化的T1-MRI模板和MEG原子磁力计的坐标信息，计算得到被试者正向传递函数。

进一步地，所述的源空间信号分布估算模块具体用于：

利用所述信号提取模块提取的多个所述单频段时域信号和所述传递函数计算模块计算得到的所述被试者正向传递函数，通过向量化的Beamformer方法估算得到多个单频段时域数据的源空间信号分布；其中，给定任意空间位置r，估算得到的源空间信号分布为：

s(t,r)＝W

其中，s(t,r)为t时刻、空间位置r估算得到的源空间信号分布；b(t)为t时刻向量化的多通道数据；W(r)为向量化空间滤波系数，其向量化的表达式如下：

W(r)＝(L

上式中，L为所述传递函数计算模块中计算得到的被试者正向传递函数，C为向量化的多通道数据的协方差矩阵，T为矩阵的转置操作符。

进一步地，所述的矩阵计算模块具体用于：

通过矩阵的奇异值分解，将所述源空间信号分布估算模块估算得到的源空间信号分布分解成信号部分和噪声部分，并通过独立成分分析计算分解矩阵和混合矩阵；

其中，所述的奇异值分解公式定义如下：

对于时域子空间主成分

本发明一种MEG源定位方法及系统，其能够解决极少数探测器条件下脑磁源精确定位问题，提高MEG系统的空间分辨率，满足MEG系统的实时性定位需求。本申请利用少量的探测器单元(一般5～30个)，基于柔性头盔，进行多通道数据采集，显著减低设备硬件成本，且被试者可在自然移动条件下进行数据采集；且同时结合其空域主成分和时域主成分以及单频段时域数据的源空间分布进行联合分析，准确定位被试者的脑磁源信号分布；本申请具有较好的经济效益。

附图说明

图1为本发明MEG源定位方法的流程图；

图2为本发明MEG源定位系统的硬件架构图；

图3为本发明一实施例被试者光阻断实验源定位结果示意图。

具体实施方式

下面结合附图及具体实施例对本发明作进一步详细的说明。

参阅图1所示，是本发明MEG源定位方法较佳实施例的作业流程图。

步骤S1，利用柔性头盔和由原子磁力计组成的探测器单元，按任务态对被试者进行单任务多区域的数据采集，并将采集的多通道数据按刺激器的初始位置对齐。具体而言：

在本实施例中，将柔性头盔固定于被试者头部，记录所有探测器单元卡槽的坐标位置；将一定数量的原子磁力计(5-30个)插入柔性头盔部分区域的卡槽内，记录卡槽编号位置，同时按任务态对被试者进行单任务多通道的数据采集；然后，通过改变原子磁力计插入卡槽位置，按任务态对被试者进行3～5次单任务多通道数据的重复采集；并将采集的多通道数据按刺激器的初始位置对齐。

步骤S2，将采集的上述多通道数据，按采集区域分组进行归一化处理。具体而言：

设定采集的上述多通道数据的数据范围分别为(T

如果假设第i次重复试验的探测器测量信号为T

步骤S3，利用参考探测器通道数据，通过数值模拟方法计算实时动态的噪声系数，并对归一化处理后的数据进行降噪处理。具体而言：

其中，所述计算实时动态的噪声系数为通过合成梯度法获取，公式如下：

s＝σ-ξ·b

其中，s为降噪后结果输出，σ为测量探测器的输出，b

静息态时，测量探测器与参考探测器输出结果应该一致；利用最小二乘法，得到使s最小的情况下的实时动态噪声系数为：

ξ＝(b

任务态脑磁信号去噪后的结果则为：

out＝σ

其中，σ

步骤S4，根据降噪处理后的数据获取多通道数据的频谱信息，并从中提取单频段时域信号。具体而言：

通过对降噪处理后的数据进行快速傅里叶变换获取多通道数据的频谱信息，并在频域空间内，对获取的所述多通道数据的频谱信息分割得到感兴趣的频段，对各个感兴趣的频段分别进行逆傅里叶变换获取对应频段的时域数据。

步骤S5，加载被试者T1磁共振扫描图像，利用被试者的T1磁共振图像计算得到其个性化T1-MRI模板，将MEG原子磁力计与T1-MRI模板坐标系配对，并计算得到被试者正向传递函数。

具体包括：

加载被试者T1磁共振扫描图像；

计算个性化被试者T1-MRI模板图像：利用的被试者的T1磁共振图像(T1-MRI)计算得到其个性化模板T1-MRI模板；

MEG原子磁力计与T1-MRI模板坐标系配对：将MEG原子磁力计的坐标信息与T1-MRI模板的基准点配对，并计算得到其在个性化模板中的坐标信息；

计算被试者正向传递函数：利用个性化的T1-MRI模板和MEG原子磁力计的坐标信息，计算得到被试者正向传递函数。

步骤S6，利用步骤S4提取的多个所述单频段时域信号和步骤S5计算得到的所述被试者正向传递函数，估算得到多个单频段时域数据的源空间信号分布。具体而言：

利用步骤S4提取的多个所述单频段时域信号和步骤S5计算得到的所述被试者正向传递函数，通过向量化的Beamformer方法估算得到多个单频段时域数据的源空间信号分布。其中，给定任意空间位置r，估算得到的源空间信号分布为：

s(t,r)＝W

其中，s(t,r)为t时刻、空间位置r估算得到的源空间信号分布；b(t)为t时刻向量化的多通道数据；W(r)为向量化空间滤波系数，其向量化的表达式如下：

W(r)＝(L

上式中，L为步骤S5中计算得到的被试者正向传递函数，C为向量化的多通道数据的协方差矩阵，T为矩阵的转置操作符。

步骤S7，将步骤S6估算得到的上述源空间信号分布分解成信号部分和噪声部分，并计算得到分解矩阵和混合矩阵。具体而言：

通过矩阵的奇异值分解，将步骤S6估算得到的源空间信号分布分解成信号部分和噪声部分，并通过独立成分分析计算分解矩阵和混合矩阵；

其中，所述的奇异值分解公式定义如下：

对于时域子空间主成分

步骤S8，利用步骤S7计算得到的所述分解矩阵和混合矩阵，对步骤S7中的信号部分进行分解，得到其空域主成分和时域主成分，并结合步骤S6得到的单频段时域数据的源空间信号分布，准确定位被试者的脑磁源信号分布。具体而言：

利用步骤S7计算得到的分解矩阵和混合矩阵，对步骤S7中的信号部分

在本实施例中，对步骤S7中的时域子空间主成分进行独立成分分析，得到分解矩阵，对分解矩阵求逆得到混合矩阵；利用混合矩阵对源空间信号进行分解，计算得到其独立时间成分的Map图，利用独立时间成分的Map图对步骤S6单频段时域数据的源空间信号分布进行验证，将源空间信号分布一致性好的作为备选的源空间信号估计。

参阅图2所示，是本发明MEG源定位系统10的硬件架构图。该系统包括：数据采集模块101、归一化处理模块102、降噪处理模块103、信号提取模块104、传递函数计算模块105、源空间信号分布估算模块106、矩阵计算模块107、定位模块108。

所述数据采集模块101用于利用柔性头盔和由原子磁力计组成的探测器单元，按任务态对被试者进行单任务多区域的数据采集，并将采集的多通道数据按刺激器的初始位置对齐。具体而言：

在本实施例中，将柔性头盔固定于被试者头部，记录所有探测器单元卡槽的坐标位置；将一定数量的原子磁力计(5-30个)插入柔性头盔部分区域的卡槽内，记录卡槽编号位置，同时按任务态对被试者进行单任务多通道的数据采集；然后，通过改变原子磁力计插入卡槽位置，按任务态对被试者进行3～5次单任务多通道数据的重复采集；并将采集的多通道数据按刺激器的初始位置对齐。

所述归一化处理模块102用于将采集的上述多通道数据，按采集区域分组进行归一化处理。具体而言：

所述归一化处理模块102设定采集的上述多通道数据的数据范围分别为(T

如果假设第i次重复试验的探测器测量信号为T

所述降噪处理模块103用于利用参考探测器通道数据，通过数值模拟方法计算实时动态的噪声系数，并对归一化处理后的数据进行降噪处理。具体而言：

其中，所述计算实时动态的噪声系数为通过合成梯度法获取，公式如下：

s＝σ-ξ·b

其中，s为降噪后结果输出，σ为测量探测器的输出，b

静息态时，测量探测器与参考探测器输出结果应该一致；利用最小二乘法，得到使s最小的情况下的实时动态噪声系数为：

ξ＝(b

任务态脑磁信号去噪后的结果则为：

out＝σ

其中，σ

所述信号提取模块104用于根据降噪处理后的数据获取多通道数据的频谱信息，并从中提取单频段时域信号。具体而言：

所述信号提取模块104通过对降噪处理后的数据进行快速傅里叶变换获取多通道数据的频谱信息，并在频域空间内，对获取的所述多通道数据的频谱信息分割得到感兴趣的频段，对各个感兴趣的频段分别进行逆傅里叶变换获取对应频段的时域数据。

所述传递函数计算模块105用于加载被试者T1磁共振扫描图像，利用被试者的T1磁共振图像计算得到其个性化T1-MRI模板，将MEG原子磁力计与T1-MRI模板坐标系配对，并计算得到被试者正向传递函数。

具体包括：

加载被试者T1磁共振扫描图像；

计算个性化被试者T1-MRI模板图像：利用的被试者的T1磁共振图像(T1-MRI)计算得到其个性化模板T1-MRI模板；

MEG原子磁力计与T1-MRI模板坐标系配对：将MEG原子磁力计的坐标信息与T1-MRI模板的基准点配对，并计算得到其在个性化模板中的坐标信息；

计算被试者正向传递函数：利用个性化的T1-MRI模板和MEG原子磁力计的坐标信息，计算得到被试者正向传递函数。

所述源空间信号分布估算模块106用于利用信号提取模块104提取的多个所述单频段时域信号和传递函数计算模块105计算得到的所述被试者正向传递函数，估算得到多个单频段时域数据的源空间信号分布。具体而言：

利用信号提取模块104提取的多个所述单频段时域信号和传递函数计算模块105计算得到的所述被试者正向传递函数，通过向量化的Beamformer方法估算得到多个单频段时域数据的源空间信号分布。其中，给定任意空间位置r，估算得到的源空间信号分布为：

s(t,r)＝W

其中，s(t,r)为t时刻、空间位置r估算得到的源空间信号分布；b(t)为t时刻向量化的多通道数据；W(r)为向量化空间滤波系数，其向量化的表达式如下：

W(r)＝(L

上式中，L为传递函数计算模块105计算得到的被试者正向传递函数，C为向量化的多通道数据的协方差矩阵，T为矩阵的转置操作符。

所述矩阵计算模块107用于将源空间信号分布估算模块106估算得到的上述源空间信号分布分解成信号部分和噪声部分，并计算得到分解矩阵和混合矩阵。具体而言：

通过矩阵的奇异值分解，将源空间信号分布估算模块106估算得到的源空间信号分布分解成信号部分和噪声部分，并通过独立成分分析计算分解矩阵和混合矩阵；

其中，所述的奇异值分解公式定义如下：

对于时域子空间主成分

所述定位模块108用于利用矩阵计算模块107计算得到的所述分解矩阵和混合矩阵，对矩阵计算模块107中的信号部分进行分解，得到其空域主成分和时域主成分，并结合源空间信号分布估算模块106得到的单频段时域数据的源空间信号分布，准确定位被试者的脑磁源信号分布。具体而言：

利用矩阵计算模块107计算得到的分解矩阵和混合矩阵，对矩阵计算模块107中的信号部分

在本实施例中，对矩阵计算模块107中的时域子空间主成分进行独立成分分析，得到分解矩阵，对分解矩阵求逆得到混合矩阵；利用混合矩阵对源空间信号进行分解，计算得到其独立时间成分的Map图，利用独立时间成分的Map图对源空间信号分布估算模块106单频段时域数据的源空间信号分布进行验证，将源空间信号分布一致性好的作为备选的源空间信号估计。

值得说明的是，在本申请实施例中：

所述的刺激器：一般由外部装置接入，包括：视觉刺激器、味觉刺激器、听觉刺激器以及触觉刺激器等；

所述的测量探测器和所述的参考探测器均为由原子磁力计组成的探测器单元，置于柔性头盔上测量人大脑磁场的为测量探测器，置于头部外围测量环境磁场用于补偿的为参考探测器；进一步地，在本申请实施例将测量探测器和参考探测器按照不同的信号标记区分为任务态探测器、任务态参考探测器等。

本发明另一实施例以5个原子磁力计通道，重复3次光阻断实验采集的多通道数据，验证本发明达到的技术效果，如图3所示。本发明另一实施例可以准确定位到光阻断实验中，脑磁源信号在后枕叶部位活动较强。

本申请利用少量的探测器单元(一般5～30个)，基于柔性头盔，进行多通道数据采集，显著减低设备硬件成本，且被试者可在自然移动条件下进行数据采集；且同时结合其空域主成分和时域主成分以及单频段时域数据的源空间分布进行联合分析，准确定位被试者的脑磁源信号分布；本申请具有较好的经济效益。

虽然本发明参照当前的较佳实施方式进行了描述，但本领域的技术人员应能理解，上述较佳实施方式仅用来说明本发明，并非用来限定本发明的保护范围，任何在本发明的精神和原则范围之内，所做的任何修饰、等效替换、改进等，均应包含在本发明的权利保护范围之内。