CN112914578B - Meg源定位方法及系统 - Google Patents
Meg源定位方法及系统 Download PDFInfo
- Publication number
- CN112914578B CN112914578B CN202110073793.XA CN202110073793A CN112914578B CN 112914578 B CN112914578 B CN 112914578B CN 202110073793 A CN202110073793 A CN 202110073793A CN 112914578 B CN112914578 B CN 112914578B
- Authority
- CN
- China
- Prior art keywords
- calculating
- matrix
- source
- tested person
- meg
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 47
- 238000012546 transfer Methods 0.000 claims abstract description 54
- 238000012545 processing Methods 0.000 claims abstract description 32
- 210000004556 brain Anatomy 0.000 claims abstract description 25
- 238000001228 spectrum Methods 0.000 claims abstract description 17
- 238000010606 normalization Methods 0.000 claims abstract description 14
- 239000011159 matrix material Substances 0.000 claims description 96
- 238000000354 decomposition reaction Methods 0.000 claims description 43
- 238000004364 calculation method Methods 0.000 claims description 40
- 238000002156 mixing Methods 0.000 claims description 18
- 238000005259 measurement Methods 0.000 claims description 17
- 238000000605 extraction Methods 0.000 claims description 13
- 238000012880 independent component analysis Methods 0.000 claims description 13
- FGUUSXIOTUKUDN-IBGZPJMESA-N C1(=CC=CC=C1)N1C2=C(NC([C@H](C1)NC=1OC(=NN=1)C1=CC=CC=C1)=O)C=CC=C2 Chemical compound C1(=CC=CC=C1)N1C2=C(NC([C@H](C1)NC=1OC(=NN=1)C1=CC=CC=C1)=O)C=CC=C2 FGUUSXIOTUKUDN-IBGZPJMESA-N 0.000 claims description 8
- 238000011068 loading method Methods 0.000 claims description 8
- GNFTZDOKVXKIBK-UHFFFAOYSA-N 3-(2-methoxyethoxy)benzohydrazide Chemical compound COCCOC1=CC=CC(C(=O)NN)=C1 GNFTZDOKVXKIBK-UHFFFAOYSA-N 0.000 claims description 6
- 230000000284 resting effect Effects 0.000 claims description 6
- 238000004088 simulation Methods 0.000 claims description 6
- 230000004807 localization Effects 0.000 claims description 5
- 230000001629 suppression Effects 0.000 abstract 1
- 238000002582 magnetoencephalography Methods 0.000 description 41
- 238000012360 testing method Methods 0.000 description 7
- 238000001514 detection method Methods 0.000 description 5
- 238000010586 diagram Methods 0.000 description 5
- 238000004422 calculation algorithm Methods 0.000 description 4
- 210000003128 head Anatomy 0.000 description 4
- 241000238366 Cephalopoda Species 0.000 description 3
- 230000000903 blocking effect Effects 0.000 description 3
- 241000712899 Lymphocytic choriomeningitis mammarenavirus Species 0.000 description 2
- 238000004458 analytical method Methods 0.000 description 2
- 238000013461 design Methods 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 238000002474 experimental method Methods 0.000 description 2
- 229910052734 helium Inorganic materials 0.000 description 2
- 239000001307 helium Substances 0.000 description 2
- SWQJXJOGLNCZEY-UHFFFAOYSA-N helium atom Chemical compound [He] SWQJXJOGLNCZEY-UHFFFAOYSA-N 0.000 description 2
- 239000007788 liquid Substances 0.000 description 2
- 210000004761 scalp Anatomy 0.000 description 2
- 238000001816 cooling Methods 0.000 description 1
- 230000001339 gustatory effect Effects 0.000 description 1
- 238000003384 imaging method Methods 0.000 description 1
- 238000005305 interferometry Methods 0.000 description 1
- 238000002955 isolation Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 210000000869 occipital lobe Anatomy 0.000 description 1
- 230000035945 sensitivity Effects 0.000 description 1
- 230000000007 visual effect Effects 0.000 description 1
Classifications
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/0033—Features or image-related aspects of imaging apparatus classified in A61B5/00, e.g. for MRI, optical tomography or impedance tomography apparatus; arrangements of imaging apparatus in a room
- A61B5/0035—Features or image-related aspects of imaging apparatus classified in A61B5/00, e.g. for MRI, optical tomography or impedance tomography apparatus; arrangements of imaging apparatus in a room adapted for acquisition of images from more than one imaging mode, e.g. combining MRI and optical tomography
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/0033—Features or image-related aspects of imaging apparatus classified in A61B5/00, e.g. for MRI, optical tomography or impedance tomography apparatus; arrangements of imaging apparatus in a room
- A61B5/004—Features or image-related aspects of imaging apparatus classified in A61B5/00, e.g. for MRI, optical tomography or impedance tomography apparatus; arrangements of imaging apparatus in a room adapted for image acquisition of a particular organ or body part
- A61B5/0042—Features or image-related aspects of imaging apparatus classified in A61B5/00, e.g. for MRI, optical tomography or impedance tomography apparatus; arrangements of imaging apparatus in a room adapted for image acquisition of a particular organ or body part for the brain
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/05—Detecting, measuring or recording for diagnosis by means of electric currents or magnetic fields; Measuring using microwaves or radio waves
- A61B5/055—Detecting, measuring or recording for diagnosis by means of electric currents or magnetic fields; Measuring using microwaves or radio waves involving electronic [EMR] or nuclear [NMR] magnetic resonance, e.g. magnetic resonance imaging
Landscapes
- Health & Medical Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
- Radiology & Medical Imaging (AREA)
- Physics & Mathematics (AREA)
- Medical Informatics (AREA)
- Surgery (AREA)
- Engineering & Computer Science (AREA)
- Biomedical Technology (AREA)
- Heart & Thoracic Surgery (AREA)
- Biophysics (AREA)
- Molecular Biology (AREA)
- Pathology (AREA)
- Animal Behavior & Ethology (AREA)
- General Health & Medical Sciences (AREA)
- Public Health (AREA)
- Veterinary Medicine (AREA)
- High Energy & Nuclear Physics (AREA)
- Neurology (AREA)
- Measurement And Recording Of Electrical Phenomena And Electrical Characteristics Of The Living Body (AREA)
Abstract
本发明涉及一种MEG源定位方法，包括：结合柔性头盔装置和原子磁力计对被试者进行单任务多区域的数据采集，对多区域多通道数据进行归一化处理；计算实时动态的噪声系数，对多通道数据进行噪声抑制；针对降噪后的数据获取多通道数据的频谱信息，并在频域空间内分割感兴趣的频段，对各个感兴趣的频段获取对应频段的时域数据；利用的被试者的T1磁共振模板、MEG原子磁力计坐标信息计算被试者正向传递函数；估计多频段时域数据的源空间分布，准确定位被试者脑磁源信号分布。本发明还涉及一种MEG源定位系统。本发明能够解决极少数探测器条件下脑磁源精确定位问题，提高MEG系统的空间分辨率，满足MEG系统的实时性定位需求。
Description
技术领域
本发明涉及一种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＝σ-ξ·bref
其中，s为降噪后结果输出，σ为测量探测器的输出，bref为参考探测器的输出，ξ为实时动态噪声系数；
静息态时，测量探测器与参考探测器输出结果应该一致；利用最小二乘法，得到使s最小的情况下的实时动态噪声系数为：
ξ＝(bTb)-1bTσ
任务态脑磁信号去噪后的结果则为：
out＝σact-ξ·bref·act
其中，σact为任务态探测器的输出，bref·act为任务态参考探测器的输出。
进一步地，所述的步骤e具体包括：加载被试者T1磁共振扫描图像；
计算个性化被试者T1-MRI模板图像：利用的被试者的T1磁共振图像计算得到其个性化模板T1-MRI模板；
MEG原子磁力计与T1-MRI模板坐标系配对：将MEG原子磁力计的坐标信息与T1-MRI模板的基准点配对，并计算得到其在个性化模板中的坐标信息；
计算被试者正向传递函数：利用个性化的T1-MRI模板和MEG原子磁力计的坐标信息，计算得到被试者正向传递函数。
进一步地，所述的步骤f具体包括：
利用步骤d提取的多个所述单频段时域信号和步骤e计算得到的所述被试者正向传递函数，通过向量化的Beamformer方法估算得到多个单频段时域数据的源空间信号分布；其中，给定任意空间位置r，估算得到的源空间信号分布为：
s(t,r)＝WT(r)·b(t)
其中，s(t,r)为t时刻、空间位置r估算得到的源空间信号分布；b(t)为t时刻向量化的多通道数据；W(r)为向量化空间滤波系数，其向量化的表达式如下：
W(r)＝(LT(r)C-1L(r))-1·LT(r)C-1
上式中，L为步骤e中计算得到的被试者正向传递函数，C为向量化的多通道数据的协方差矩阵，T为矩阵的转置操作符。
进一步地，所述的步骤g具体包括：
通过矩阵的奇异值分解，将步骤f估算得到的源空间信号分布分解成信号部分和噪声部分，并通过独立成分分析计算分解矩阵和混合矩阵；
其中，所述的奇异值分解公式定义如下：
对于时域子空间主成分通过独立成分分析寻找得到分解矩阵H与混合矩阵H-1对源空间信号进行分解，计算得到其未知的源信号分量。
本发明提供一种MEG源定位系统，该系统包括：数据采集模块、归一化处理模块、降噪处理模块、信号提取模块、传递函数计算模块、源空间信号分布估算模块、矩阵计算模块、定位模块，其中：所述数据采集模块用于利用柔性头盔和由原子磁力计组成的探测器单元，按任务态对被试者进行单任务多区域的数据采集，并将采集的多通道数据按刺激器的初始位置对齐；所述归一化处理模块用于将采集的上述多通道数据，按采集区域分组进行归一化处理；所述降噪处理模块用于利用参考探测器通道数据，通过数值模拟方法计算实时动态的噪声系数，并对归一化处理后的数据进行降噪处理；所述信号提取模块用于根据降噪处理后的数据获取多通道数据的频谱信息，并从中提取单频段时域信号；所述传递函数计算模块用于利用被试者的T1磁共振图像计算得到其个性化T1-MRI模板，将MEG原子磁力计与T1-MRI模板坐标系配对，并计算得到被试者正向传递函数；所述源空间信号分布估算模块用于利用所述信号提取模块提取的多个所述单频段时域信号和所述传递函数计算模块计算得到的所述被试者正向传递函数，估算得到多个单频段时域数据的源空间信号分布；所述矩阵计算模块用于将所述源空间信号分布估算模块估算得到的上述源空间信号分布分解成信号部分和噪声部分，并计算得到分解矩阵和混合矩阵；所述定位模块用于利用所述矩阵计算模块计算得到的所述分解矩阵和混合矩阵，对所述矩阵计算模块中的信号部分进行分解，得到其空域主成分和时域主成分，并结合所述源空间信号分布估算模块得到的单频段时域数据的源空间信号分布，准确定位被试者的脑磁源信号分布。
进一步地，所述的降噪处理模块具体用于：
所述计算实时动态的噪声系数为通过合成梯度法获取，公式如下：
s＝σ-ξ·bref
其中，s为降噪后结果输出，σ为测量探测器的输出，bref为参考探测器的输出，ξ为实时动态噪声系数；
静息态时，测量探测器与参考探测器输出结果应该一致；利用最小二乘法，得到使s最小的情况下的实时动态噪声系数为：
ξ＝(bTb)-1bTσ
任务态脑磁信号去噪后的结果则为：
out＝σact-ξ·bref·act
其中，σact为任务态探测器的输出，bref·act为任务态参考探测器的输出。
进一步地，所述的传递函数计算模块具体用于：
加载被试者T1磁共振扫描图像；
计算个性化被试者T1-MRI模板图像：利用的被试者的T1磁共振图像计算得到其个性化模板T1-MRI模板；
MEG原子磁力计与T1-MRI模板坐标系配对：将MEG原子磁力计的坐标信息与T1-MRI模板的基准点配对，并计算得到其在个性化模板中的坐标信息；
计算被试者正向传递函数：利用个性化的T1-MRI模板和MEG原子磁力计的坐标信息，计算得到被试者正向传递函数。
进一步地，所述的源空间信号分布估算模块具体用于：
利用所述信号提取模块提取的多个所述单频段时域信号和所述传递函数计算模块计算得到的所述被试者正向传递函数，通过向量化的Beamformer方法估算得到多个单频段时域数据的源空间信号分布；其中，给定任意空间位置r，估算得到的源空间信号分布为：
s(t,r)＝WT(r)·b(t)
其中，s(t,r)为t时刻、空间位置r估算得到的源空间信号分布；b(t)为t时刻向量化的多通道数据；W(r)为向量化空间滤波系数，其向量化的表达式如下：
W(r)＝(LT(r)C-1L(r))-1·LT(r)C-1
上式中，L为所述传递函数计算模块中计算得到的被试者正向传递函数，C为向量化的多通道数据的协方差矩阵，T为矩阵的转置操作符。
进一步地，所述的矩阵计算模块具体用于：
通过矩阵的奇异值分解，将所述源空间信号分布估算模块估算得到的源空间信号分布分解成信号部分和噪声部分，并通过独立成分分析计算分解矩阵和混合矩阵；
其中，所述的奇异值分解公式定义如下：
对于时域子空间主成分通过独立成分分析寻找得到分解矩阵H与混合矩阵H-1对源空间信号进行分解，计算得到其未知的源信号分量。
本发明一种MEG源定位方法及系统，其能够解决极少数探测器条件下脑磁源精确定位问题，提高MEG系统的空间分辨率，满足MEG系统的实时性定位需求。本申请利用少量的探测器单元(一般5～30个)，基于柔性头盔，进行多通道数据采集，显著减低设备硬件成本，且被试者可在自然移动条件下进行数据采集；且同时结合其空域主成分和时域主成分以及单频段时域数据的源空间分布进行联合分析，准确定位被试者的脑磁源信号分布；本申请具有较好的经济效益。
附图说明
图1为本发明MEG源定位方法的流程图；
图2为本发明MEG源定位系统的硬件架构图；
图3为本发明一实施例被试者光阻断实验源定位结果示意图。
具体实施方式
下面结合附图及具体实施例对本发明作进一步详细的说明。
参阅图1所示，是本发明MEG源定位方法较佳实施例的作业流程图。
步骤S1，利用柔性头盔和由原子磁力计组成的探测器单元，按任务态对被试者进行单任务多区域的数据采集，并将采集的多通道数据按刺激器的初始位置对齐。具体而言：
在本实施例中，将柔性头盔固定于被试者头部，记录所有探测器单元卡槽的坐标位置；将一定数量的原子磁力计(5-30个)插入柔性头盔部分区域的卡槽内，记录卡槽编号位置，同时按任务态对被试者进行单任务多通道的数据采集；然后，通过改变原子磁力计插入卡槽位置，按任务态对被试者进行3～5次单任务多通道数据的重复采集；并将采集的多通道数据按刺激器的初始位置对齐。
步骤S2，将采集的上述多通道数据，按采集区域分组进行归一化处理。具体而言：
设定采集的上述多通道数据的数据范围分别为(T1 min,T1 max)、(T2 min,T2 max)、…(Ti min,Ti max)，其中，i为多通道数据采集编号，Ti min为第i次重复试验多通道数据的最小阈值，Ti max为第i次重复试验多通道数据的最大阈值；利用多通道数据的阈值范围，对多次采集的多通道数据做归一化处理。
如果假设第i次重复试验的探测器测量信号为Ti，则归一化的测量数值Ti norm：
步骤S3，利用参考探测器通道数据，通过数值模拟方法计算实时动态的噪声系数，并对归一化处理后的数据进行降噪处理。具体而言：
其中，所述计算实时动态的噪声系数为通过合成梯度法获取，公式如下：
s＝σ-ξ·bref
其中，s为降噪后结果输出，σ为测量探测器的输出，bref为参考探测器的输出，ξ为实时动态噪声系数。
静息态时，测量探测器与参考探测器输出结果应该一致；利用最小二乘法，得到使s最小的情况下的实时动态噪声系数为：
ξ＝(bTb)-1bTσ
任务态脑磁信号去噪后的结果则为：
out＝σact-ξ·bref·act
其中，σact为任务态探测器的输出，bref·act为任务态参考探测器的输出。
步骤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)＝WT(r)·b(t)
其中，s(t,r)为t时刻、空间位置r估算得到的源空间信号分布；b(t)为t时刻向量化的多通道数据；W(r)为向量化空间滤波系数，其向量化的表达式如下：
W(r)＝(LT(r)C-1L(r))-1·LT(r)C-1
上式中，L为步骤S5中计算得到的被试者正向传递函数，C为向量化的多通道数据的协方差矩阵，T为矩阵的转置操作符。
步骤S7，将步骤S6估算得到的上述源空间信号分布分解成信号部分和噪声部分，并计算得到分解矩阵和混合矩阵。具体而言：
通过矩阵的奇异值分解，将步骤S6估算得到的源空间信号分布分解成信号部分和噪声部分，并通过独立成分分析计算分解矩阵和混合矩阵；
其中，所述的奇异值分解公式定义如下：
对于时域子空间主成分通过独立成分分析寻找得到分解矩阵H与混合矩阵H-1对源空间信号进行分解，计算得到其未知的源信号分量。
步骤S8，利用步骤S7计算得到的所述分解矩阵和混合矩阵，对步骤S7中的信号部分进行分解，得到其空域主成分和时域主成分，并结合步骤S6得到的单频段时域数据的源空间信号分布，准确定位被试者的脑磁源信号分布。具体而言：
利用步骤S7计算得到的分解矩阵和混合矩阵，对步骤S7中的信号部分进行分解，得到其空域主成分和时域主成分，并结合步骤S6估算得到的单频段时域数据的源空间分布进行分析，准确定位被试者的脑磁源信号分布。
在本实施例中，对步骤S7中的时域子空间主成分进行独立成分分析，得到分解矩阵，对分解矩阵求逆得到混合矩阵；利用混合矩阵对源空间信号进行分解，计算得到其独立时间成分的Map图，利用独立时间成分的Map图对步骤S6单频段时域数据的源空间信号分布进行验证，将源空间信号分布一致性好的作为备选的源空间信号估计。
参阅图2所示，是本发明MEG源定位系统10的硬件架构图。该系统包括：数据采集模块101、归一化处理模块102、降噪处理模块103、信号提取模块104、传递函数计算模块105、源空间信号分布估算模块106、矩阵计算模块107、定位模块108。
所述数据采集模块101用于利用柔性头盔和由原子磁力计组成的探测器单元，按任务态对被试者进行单任务多区域的数据采集，并将采集的多通道数据按刺激器的初始位置对齐。具体而言：
在本实施例中，将柔性头盔固定于被试者头部，记录所有探测器单元卡槽的坐标位置；将一定数量的原子磁力计(5-30个)插入柔性头盔部分区域的卡槽内，记录卡槽编号位置，同时按任务态对被试者进行单任务多通道的数据采集；然后，通过改变原子磁力计插入卡槽位置，按任务态对被试者进行3～5次单任务多通道数据的重复采集；并将采集的多通道数据按刺激器的初始位置对齐。
所述归一化处理模块102用于将采集的上述多通道数据，按采集区域分组进行归一化处理。具体而言：
所述归一化处理模块102设定采集的上述多通道数据的数据范围分别为(T1 min,T1 max)、(T2 min,T2 max)、…(Ti min,Ti max)，其中，i为多通道数据采集编号，Ti min为第i次重复试验多通道数据的最小阈值，Ti max为第i次重复试验多通道数据的最大阈值；利用多通道数据的阈值范围，对多次采集的多通道数据做归一化处理。
如果假设第i次重复试验的探测器测量信号为Ti，则归一化的测量数值Ti norm：
所述降噪处理模块103用于利用参考探测器通道数据，通过数值模拟方法计算实时动态的噪声系数，并对归一化处理后的数据进行降噪处理。具体而言：
其中，所述计算实时动态的噪声系数为通过合成梯度法获取，公式如下：
s＝σ-ξ·bref
其中，s为降噪后结果输出，σ为测量探测器的输出，bref为参考探测器的输出，ξ为实时动态噪声系数。
静息态时，测量探测器与参考探测器输出结果应该一致；利用最小二乘法，得到使s最小的情况下的实时动态噪声系数为：
ξ＝(bTb)-1bTσ
任务态脑磁信号去噪后的结果则为：
out＝σact-ξ·bref·act
其中，σact为任务态探测器的输出，bref·act为任务态参考探测器的输出。
所述信号提取模块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)＝WT(r)·b(t)
其中，s(t,r)为t时刻、空间位置r估算得到的源空间信号分布；b(t)为t时刻向量化的多通道数据；W(r)为向量化空间滤波系数，其向量化的表达式如下：
W(r)＝(LT(r)C-1L(r))-1·LT(r)C-1
上式中，L为传递函数计算模块105计算得到的被试者正向传递函数，C为向量化的多通道数据的协方差矩阵，T为矩阵的转置操作符。
所述矩阵计算模块107用于将源空间信号分布估算模块106估算得到的上述源空间信号分布分解成信号部分和噪声部分，并计算得到分解矩阵和混合矩阵。具体而言：
通过矩阵的奇异值分解，将源空间信号分布估算模块106估算得到的源空间信号分布分解成信号部分和噪声部分，并通过独立成分分析计算分解矩阵和混合矩阵；
其中，所述的奇异值分解公式定义如下：
对于时域子空间主成分通过独立成分分析寻找得到分解矩阵H与混合矩阵H-1对源空间信号进行分解，计算得到其未知的源信号分量。
所述定位模块108用于利用矩阵计算模块107计算得到的所述分解矩阵和混合矩阵，对矩阵计算模块107中的信号部分进行分解，得到其空域主成分和时域主成分，并结合源空间信号分布估算模块106得到的单频段时域数据的源空间信号分布，准确定位被试者的脑磁源信号分布。具体而言：
利用矩阵计算模块107计算得到的分解矩阵和混合矩阵，对矩阵计算模块107中的信号部分进行分解，得到其空域主成分和时域主成分，并结合源空间信号分布估算模块106估算得到的单频段时域数据的源空间分布进行分析，准确定位被试者的脑磁源信号分布。
在本实施例中，对矩阵计算模块107中的时域子空间主成分进行独立成分分析，得到分解矩阵，对分解矩阵求逆得到混合矩阵；利用混合矩阵对源空间信号进行分解，计算得到其独立时间成分的Map图，利用独立时间成分的Map图对源空间信号分布估算模块106单频段时域数据的源空间信号分布进行验证，将源空间信号分布一致性好的作为备选的源空间信号估计。
值得说明的是，在本申请实施例中：
所述的刺激器：一般由外部装置接入，包括：视觉刺激器、味觉刺激器、听觉刺激器以及触觉刺激器等；
所述的测量探测器和所述的参考探测器均为由原子磁力计组成的探测器单元，置于柔性头盔上测量人大脑磁场的为测量探测器，置于头部外围测量环境磁场用于补偿的为参考探测器；进一步地，在本申请实施例将测量探测器和参考探测器按照不同的信号标记区分为任务态探测器、任务态参考探测器等。
本发明另一实施例以5个原子磁力计通道，重复3次光阻断实验采集的多通道数据，验证本发明达到的技术效果，如图3所示。本发明另一实施例可以准确定位到光阻断实验中，脑磁源信号在后枕叶部位活动较强。
本申请利用少量的探测器单元(一般5～30个)，基于柔性头盔，进行多通道数据采集，显著减低设备硬件成本，且被试者可在自然移动条件下进行数据采集；且同时结合其空域主成分和时域主成分以及单频段时域数据的源空间分布进行联合分析，准确定位被试者的脑磁源信号分布；本申请具有较好的经济效益。
虽然本发明参照当前的较佳实施方式进行了描述，但本领域的技术人员应能理解，上述较佳实施方式仅用来说明本发明，并非用来限定本发明的保护范围，任何在本发明的精神和原则范围之内，所做的任何修饰、等效替换、改进等，均应包含在本发明的权利保护范围之内。
Claims (6)
1.一种MEG源定位方法，其特征在于，该方法包括如下步骤：
a.利用柔性头盔和由原子磁力计组成的探测器单元，按任务态对被试者进行单任务多区域的数据采集，并将采集的多通道数据按刺激器的初始位置对齐；
b.将采集的上述多通道数据，按采集区域分组进行归一化处理；
c.利用参考探测器通道数据，通过数值模拟方法计算实时动态的噪声系数，并对归一化处理后的数据进行降噪处理；
d.根据降噪处理后的数据获取多通道数据的频谱信息，并从中提取单频段时域信号；
e.利用被试者的T1磁共振图像计算得到其个性化T1-MRI模板，将MEG原子磁力计与T1-MRI模板坐标系配对，并计算得到被试者正向传递函数；
f.利用步骤d提取的多个所述单频段时域信号和步骤e计算得到的所述被试者正向传递函数，估算得到多个单频段时域数据的源空间信号分布；
g.将步骤f估算得到的上述源空间信号分布分解成信号部分和噪声部分，并计算得到分解矩阵和混合矩阵；
h.利用步骤g计算得到的所述分解矩阵和混合矩阵，对步骤g中的信号部分进行分解，得到其空域主成分和时域主成分，并结合步骤f得到的单频段时域数据的源空间信号分布，准确定位被试者的脑磁源信号分布；
其中，所述的步骤c具体包括：
所述计算实时动态的噪声系数为通过合成梯度法获取，公式如下：
s＝σ-ξ·bref
其中，s为降噪后结果输出，σ为测量探测器的输出，bref为参考探测器的输出，ξ为实时动态噪声系数；
静息态时，测量探测器与参考探测器输出结果应该一致；利用最小二乘法，得到使s最小的情况下的实时动态噪声系数为：
ξ＝(bTb)-1bTσ
任务态脑磁信号去噪后的结果则为：
out＝σact-ξ·bref·act
其中，σact为任务态探测器的输出，bref·act为任务态参考探测器的输出；
所述的步骤f具体包括：
利用步骤d提取的多个所述单频段时域信号和步骤e计算得到的所述被试者正向传递函数，通过向量化的Beamformer方法估算得到多个单频段时域数据的源空间信号分布；其中，给定任意空间位置r，估算得到的源空间信号分布为：
s(t,r)＝WT(r)·b(t)
其中，s(t,r)为t时刻、空间位置r估算得到的源空间信号分布；b(t)为t时刻向量化的多通道数据；W(r)为向量化空间滤波系数，其向量化的表达式如下：
W(r)＝(LT(r)C-1L(r))-1·LT(r)C-1
上式中，L为步骤e中计算得到的被试者正向传递函数，C为向量化的多通道数据的协方差矩阵，T为矩阵的转置操作符。
2.如权利要求1所述的方法，其特征在于，所述的步骤e具体包括：
加载被试者T1磁共振扫描图像；
计算个性化被试者T1-MRI模板图像：利用的被试者的T1磁共振图像计算得到其个性化模板T1-MRI模板；
MEG原子磁力计与T1-MRI模板坐标系配对：将MEG原子磁力计的坐标信息与T1-MRI模板的基准点配对，并计算得到其在个性化模板中的坐标信息；
计算被试者正向传递函数：利用个性化的T1-MRI模板和MEG原子磁力计的坐标信息，计算得到被试者正向传递函数。
3.如权利要求2所述的方法，其特征在于，所述的步骤g具体包括：
通过矩阵的奇异值分解，将步骤f估算得到的源空间信号分布分解成信号部分和噪声部分，并通过独立成分分析计算分解矩阵和混合矩阵；
其中，所述的奇异值分解公式定义如下：
对于时域子空间主成分通过独立成分分析寻找得到分解矩阵H与混合矩阵H-1对源空间信号进行分解，计算得到其未知的源信号分量。
4.一种MEG源定位系统，其特征在于，该系统包括：数据采集模块、归一化处理模块、降噪处理模块、信号提取模块、传递函数计算模块、源空间信号分布估算模块、矩阵计算模块、定位模块，其中：
所述数据采集模块用于利用柔性头盔和由原子磁力计组成的探测器单元，按任务态对被试者进行单任务多区域的数据采集，并将采集的多通道数据按刺激器的初始位置对齐；
所述归一化处理模块用于将采集的上述多通道数据，按采集区域分组进行归一化处理；
所述降噪处理模块用于利用参考探测器通道数据，通过数值模拟方法计算实时动态的噪声系数，并对归一化处理后的数据进行降噪处理；
所述信号提取模块用于根据降噪处理后的数据获取多通道数据的频谱信息，并从中提取单频段时域信号；
所述传递函数计算模块用于利用被试者的T1磁共振图像计算得到其个性化T1-MRI模板，将MEG原子磁力计与T1-MRI模板坐标系配对，并计算得到被试者正向传递函数；
所述源空间信号分布估算模块用于利用所述信号提取模块提取的多个所述单频段时域信号和所述传递函数计算模块计算得到的所述被试者正向传递函数，估算得到多个单频段时域数据的源空间信号分布；
所述矩阵计算模块用于将所述源空间信号分布估算模块估算得到的上述源空间信号分布分解成信号部分和噪声部分，并计算得到分解矩阵和混合矩阵；
所述定位模块用于利用所述矩阵计算模块计算得到的所述分解矩阵和混合矩阵，对所述矩阵计算模块中的信号部分进行分解，得到其空域主成分和时域主成分，并结合所述源空间信号分布估算模块得到的单频段时域数据的源空间信号分布，准确定位被试者的脑磁源信号分布；
其中，所述的降噪处理模块具体用于：
所述计算实时动态的噪声系数为通过合成梯度法获取，公式如下：
s＝σ-ξ·bref
其中，s为降噪后结果输出，σ为测量探测器的输出，bref为参考探测器的输出，ξ为实时动态噪声系数；
静息态时，测量探测器与参考探测器输出结果应该一致；利用最小二乘法，得到使s最小的情况下的实时动态噪声系数为：
ξ＝(bTb)-1bTσ
任务态脑磁信号去噪后的结果则为：
out＝σact-ξ·bref·act
其中，σact为任务态探测器的输出，bref·act为任务态参考探测器的输出；
所述的源空间信号分布估算模块具体用于：
利用所述信号提取模块提取的多个所述单频段时域信号和所述传递函数计算模块计算得到的所述被试者正向传递函数，通过向量化的Beamformer方法估算得到多个单频段时域数据的源空间信号分布；其中，给定任意空间位置r，估算得到的源空间信号分布为：
s(t,r)＝WT(r)·b(t)
其中，s(t,r)为t时刻、空间位置r估算得到的源空间信号分布；b(t)为t时刻向量化的多通道数据；W(r)为向量化空间滤波系数，其向量化的表达式如下：
W(r)＝(LT(r)C-1L(r))-1·LT(r)C-1
上式中，L为所述传递函数计算模块中计算得到的被试者正向传递函数，C为向量化的多通道数据的协方差矩阵，T为矩阵的转置操作符。
5.如权利要求4所述的系统，其特征在于，所述的传递函数计算模块具体用于：
加载被试者T1磁共振扫描图像；
计算个性化被试者T1-MRI模板图像：利用的被试者的T1磁共振图像计算得到其个性化模板T1-MRI模板；
MEG原子磁力计与T1-MRI模板坐标系配对：将MEG原子磁力计的坐标信息与T1-MRI模板的基准点配对，并计算得到其在个性化模板中的坐标信息；
计算被试者正向传递函数：利用个性化的T1-MRI模板和MEG原子磁力计的坐标信息，计算得到被试者正向传递函数。
6.如权利要求5所述的系统，其特征在于，所述的矩阵计算模块具体用于：
通过矩阵的奇异值分解，将所述源空间信号分布估算模块估算得到的源空间信号分布分解成信号部分和噪声部分，并通过独立成分分析计算分解矩阵和混合矩阵；
其中，所述的奇异值分解公式定义如下：
对于时域子空间主成分通过独立成分分析寻找得到分解矩阵H与混合矩阵H-1对源空间信号进行分解，计算得到其未知的源信号分量。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110073793.XA CN112914578B (zh) | 2021-01-20 | 2021-01-20 | Meg源定位方法及系统 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110073793.XA CN112914578B (zh) | 2021-01-20 | 2021-01-20 | Meg源定位方法及系统 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN112914578A CN112914578A (zh) | 2021-06-08 |
| CN112914578B true CN112914578B (zh) | 2024-02-09 |
Family
ID=76164535
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202110073793.XA Active CN112914578B (zh) | 2021-01-20 | 2021-01-20 | Meg源定位方法及系统 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN112914578B (zh) |
Families Citing this family (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113951887B (zh) * | 2021-12-22 | 2022-03-25 | 北京航空航天大学杭州创新研究院 | 一种频谱匹配独立成分分析方法及系统 |
| CN114052668B (zh) * | 2022-01-17 | 2022-06-17 | 北京航空航天大学杭州创新研究院 | 一种基于脑磁图数据的脑功能分析方法 |
| CN114065825B (zh) * | 2022-01-17 | 2022-04-19 | 北京航空航天大学杭州创新研究院 | 一种基于结构相似性的脑磁meg源定位方法 |
| CN115018018B (zh) * | 2022-08-05 | 2022-11-11 | 北京航空航天大学杭州创新研究院 | 一种抑制背景噪声的双重空间滤波方法 |
Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2000010454A1 (en) * | 1998-08-24 | 2000-03-02 | Ctf Systems Inc. | Functional brain imaging from magnetoencephalographic data |
| US6697660B1 (en) * | 1998-01-23 | 2004-02-24 | Ctf Systems, Inc. | Method for functional brain imaging from magnetoencephalographic data by estimation of source signal-to-noise ratio |
| CN105147289A (zh) * | 2015-08-18 | 2015-12-16 | 高家红 | 基于原子磁力计的meg系统及方法 |
| CN105147288A (zh) * | 2015-07-23 | 2015-12-16 | 中国科学院苏州生物医学工程技术研究所 | 脑磁源强度定位方法 |
| CN105212895A (zh) * | 2015-08-24 | 2016-01-06 | 中国科学院苏州生物医学工程技术研究所 | 动态脑磁源定位方法 |
| CN106798557A (zh) * | 2017-02-27 | 2017-06-06 | 上海理工大学 | 基于原子磁力传感的脑磁信息检测分析方法 |
| CN106859599A (zh) * | 2017-01-23 | 2017-06-20 | 上海理工大学 | 基于全光原子磁力检测的脑磁图系统及获取方法 |
Family Cites Families (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2015021070A1 (en) * | 2013-08-05 | 2015-02-12 | The Regents Of The University Of California | Magnetoencephalography source imaging for neurological functionality characterizations |
| FR3025037B1 (fr) * | 2014-08-22 | 2016-09-30 | Commissariat Energie Atomique | Procede de localisation d'une activite cerebrale associee a une tache |
| CN110728704B (zh) * | 2019-11-13 | 2022-12-06 | 北京航空航天大学 | 一种基于mri和opm的脑磁多模态影像配准系统及方法 |
-
2021
- 2021-01-20 CN CN202110073793.XA patent/CN112914578B/zh active Active
Patent Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6697660B1 (en) * | 1998-01-23 | 2004-02-24 | Ctf Systems, Inc. | Method for functional brain imaging from magnetoencephalographic data by estimation of source signal-to-noise ratio |
| WO2000010454A1 (en) * | 1998-08-24 | 2000-03-02 | Ctf Systems Inc. | Functional brain imaging from magnetoencephalographic data |
| CN105147288A (zh) * | 2015-07-23 | 2015-12-16 | 中国科学院苏州生物医学工程技术研究所 | 脑磁源强度定位方法 |
| CN105147289A (zh) * | 2015-08-18 | 2015-12-16 | 高家红 | 基于原子磁力计的meg系统及方法 |
| CN105212895A (zh) * | 2015-08-24 | 2016-01-06 | 中国科学院苏州生物医学工程技术研究所 | 动态脑磁源定位方法 |
| CN106859599A (zh) * | 2017-01-23 | 2017-06-20 | 上海理工大学 | 基于全光原子磁力检测的脑磁图系统及获取方法 |
| CN106798557A (zh) * | 2017-02-27 | 2017-06-06 | 上海理工大学 | 基于原子磁力传感的脑磁信息检测分析方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN112914578A (zh) | 2021-06-08 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN112914578B (zh) | Meg源定位方法及系统 | |
| Taulu et al. | Applications of the signal space separation method | |
| US6697660B1 (en) | Method for functional brain imaging from magnetoencephalographic data by estimation of source signal-to-noise ratio | |
| Zetter et al. | Requirements for coregistration accuracy in on-scalp MEG | |
| Gross et al. | Properties of MEG tomographic maps obtained with spatial filtering | |
| US5526811A (en) | Apparatus and process for determining the sources of biomagnetic activity | |
| US5263488A (en) | Method and apparatus for localization of intracerebral sources of electrical activity | |
| Wikswo Jr et al. | The future of the EEG and MEG | |
| US9560986B2 (en) | Magnetometer for medical use | |
| JP5361131B2 (ja) | 直交仮想チャネルを使用したマルチチャネル測定データの分析 | |
| US9370309B2 (en) | Magnetoencephalography system and method for 3D localization and tracking of electrical activity in brain | |
| Vrba et al. | Fetal MEG redistribution by projection operators | |
| Nolte et al. | Current multipole expansion to estimate lateral extent of neuronal activity: a theoretical analysis | |
| Xie et al. | Benchmarking for on-scalp MEG sensors | |
| Sekihara et al. | Maximum-likelihood estimation of current-dipole parameters for data obtained using multichannel magnetometer | |
| WO2022069896A1 (en) | Magnetoencephalography method and system | |
| EP0477434B2 (en) | Analysis of biological signals using data from arrays of sensors | |
| Singh et al. | Neuromagnetic localization using magnetic resonance images | |
| CN115132349B (zh) | 一种事件相关脑活动时空频动态定位系统及方法 | |
| Primin et al. | A method and an algorithm to reconstruct the spatial structure of current density vectors in magnetocardiography | |
| CN113951886A (zh) | 一种脑磁纹生成系统及测谎决策系统 | |
| WO2000010454A1 (en) | Functional brain imaging from magnetoencephalographic data | |
| CN113951887A (zh) | 一种频谱匹配独立成分分析方法及系统 | |
| Mohseni et al. | A new beamforming-based MEG dipole source localization method | |
| Mino et al. | Neuromagnetic source analysis with a 64-channel SQUID system and MR imaging |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |