CN114052668B - 一种基于脑磁图数据的脑功能分析方法 - Google Patents

一种基于脑磁图数据的脑功能分析方法 Download PDF

Info

Publication number
CN114052668B
CN114052668B CN202210046874.5A CN202210046874A CN114052668B CN 114052668 B CN114052668 B CN 114052668B CN 202210046874 A CN202210046874 A CN 202210046874A CN 114052668 B CN114052668 B CN 114052668B
Authority
CN
China
Prior art keywords
data
magnetoencephalogram
magnetoencephalogram data
model
signal source
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
Application number
CN202210046874.5A
Other languages
English (en)
Other versions
CN114052668A (zh
Inventor
高阳
路浩
宁晓琳
房建成
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Hangzhou Innovation Research Institute of Beihang University
Original Assignee
Hangzhou Innovation Research Institute of Beihang University
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Hangzhou Innovation Research Institute of Beihang University filed Critical Hangzhou Innovation Research Institute of Beihang University
Priority to CN202210046874.5A priority Critical patent/CN114052668B/zh
Publication of CN114052668A publication Critical patent/CN114052668A/zh
Application granted granted Critical
Publication of CN114052668B publication Critical patent/CN114052668B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
    • A61B5/7235—Details of waveform analysis
    • A61B5/7264—Classification of physiological signals or data, e.g. using neural networks, statistical classifiers, expert systems or fuzzy systems
    • 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
    • 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/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
    • A61B5/7203—Signal processing specially adapted for physiological signals or for diagnostic purposes for noise prevention, reduction or removal
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
    • A61B5/7235—Details of waveform analysis
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
    • A61B5/7235—Details of waveform analysis
    • A61B5/725—Details of waveform analysis using specific filters therefor, e.g. Kalman or adaptive filters
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis

Landscapes

  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Surgery (AREA)
  • Molecular Biology (AREA)
  • Veterinary Medicine (AREA)
  • Public Health (AREA)
  • General Health & Medical Sciences (AREA)
  • Biophysics (AREA)
  • Pathology (AREA)
  • Biomedical Technology (AREA)
  • Heart & Thoracic Surgery (AREA)
  • Medical Informatics (AREA)
  • Animal Behavior & Ethology (AREA)
  • Artificial Intelligence (AREA)
  • Signal Processing (AREA)
  • Psychiatry (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Physiology (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Radiology & Medical Imaging (AREA)
  • Data Mining & Analysis (AREA)
  • Theoretical Computer Science (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Evolutionary Biology (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Operations Research (AREA)
  • Probability & Statistics with Applications (AREA)
  • Fuzzy Systems (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)

Abstract

本发明提供一种基于脑磁图数据的脑功能分析方法，包括：S1、获取被试者的各个数据通道的脑磁图数据；脑磁图数据基于用于被试者的SERF磁强计采集；S2、将脑磁图数据输入预先构建的溯源定位模型，获得脑磁图数据的信号源信息；S3、将脑磁图数据输入基于贝叶斯的多元自回归模型中，根据期望最大化算法进行模型求解，获得模型系数，基于模型系数进行显著性检验连接数据通道节点，构建传感器级脑网络；和/或，将信号源信息和脑磁图数据输入基于贝叶斯的多元自回归模型中，根据期望最大化算法进行模型求解，获得模型系数，基于模型系数进行显著性检验连接信号源节点，构建信号源级脑网络。能够构建出更精确的脑功能网络。

Description

一种基于脑磁图数据的脑功能分析方法
技术领域
本发明涉及脑功能网络分析技术领域，尤其涉及一种基于脑磁图数据的脑功能分析方法。
背景技术
人脑作为世界上中最为复杂的动力学***之一，其皮层由150-330亿个神经细胞组成，这些神经元间的相互连接构成了一个拥有相当复杂结构和功能的脑网络，在神经学和生理学方面，大脑的主要功能是控制和支配身体的各个器官，因此，这个复杂而庞大的脑网络使它具有高级的信息处理和认知表达功能，如：语言、情感、记忆、认知等，并对来自人体内部和周围环境的信息进行存储、处理、加工和整合，探究大脑的结构和功能，加速脑科学领域的研究不仅可以提高对大脑疾病的预防、诊断和治疗，更可以推动人工智能的发展。
脑磁图能够以较高的时空分辨率对大脑神经活动进行直接成像，为大脑连通分析提供详细的时空以及特定节律活动的信息，定位精细度高，能有效地捕捉深部源放电过程。为此，本发明提供一种基于脑磁图数据的脑功能分析方法。
发明内容
（一）要解决的技术问题
鉴于上述技术中存在的问题，本发明至少从一定程度上进行解决。为此，本发明的目的在于提出了一种基于脑磁图数据的脑功能分析方法，能够构建出更精确的脑功能网络。
（二）技术方案
为了达到上述目的，本发明提供一种基于脑磁图数据的脑功能分析方法，包括以下步骤：
步骤S1、获取被试者的各个数据通道的脑磁图数据；脑磁图数据基于用于被试者的SERF磁强计采集；
步骤S2、将脑磁图数据输入预先构建的溯源定位模型，获得脑磁图数据的信号源信息；
步骤S3、将脑磁图数据输入基于贝叶斯的多元自回归模型中，根据期望最大化算法进行模型求解，获得模型系数，基于模型系数进行显著性检验连接数据通道节点，构建传感器级脑网络；和/或，
将信号源信息和脑磁图数据输入基于贝叶斯的多元自回归模型中，根据期望最大化算法进行模型求解，获得模型系数，基于模型系数进行显著性检验连接信号源节点，构建信号源级脑网络。
可选地，基于贝叶斯的多元自回归模型可表示为：
Figure 441748DEST_PATH_IMAGE001
Figure 96720DEST_PATH_IMAGE002
Figure 708967DEST_PATH_IMAGE003
Figure 685275DEST_PATH_IMAGE004
Figure 916537DEST_PATH_IMAGE005
其中，
Figure 426015DEST_PATH_IMAGE006
为模型系数；
Figure 209164DEST_PATH_IMAGE007
为矩阵
Figure 436883DEST_PATH_IMAGE008
的转置；
Figure 206255DEST_PATH_IMAGE009
Figure 337285DEST_PATH_IMAGE010
，为N维脑磁图数据的联合平稳时间序列，可以包含有脑磁图数据的信号源信息；q为 通道数。
可选地，步骤S1中，脑磁图数据为被试者处于磁屏蔽室中的脑磁图数据，步骤S1还包括：获取第一参考磁图数据和第二参考磁图数据，第一参考磁图数据基于用于磁屏蔽室的SERF磁强计采集，第二参考磁图数据是在无被试者时基于用于被试者的SERF磁强计采集；
在步骤S1和步骤S2之间，还包括：步骤S12、对脑磁图数据进行预处理；预处理依次包括以下步骤：
步骤A1、去除损坏的数据通道的脑磁图数据；
步骤A2、令脑磁图数据先通过截止频率为40Hz的低通滤波器，再通过截止频率为0.1Hz的高通滤波器；
步骤A3、根据第一参考数据对脑磁图数据进行回归处理，去除背景噪声；
步骤A4、基于信号空间投影法，将第二参考数据在脑磁图数据上进行投影，去除背景噪声和生理伪迹；
步骤A5、通过自适应阈值算法得到生理伪迹的时间窗，对时间窗内数据进行主成分分析，从得到的主成分中选择生理伪迹相关的成分，在脑磁图数据中将生理伪迹相关的成分去除；
步骤A6、对脑磁图数据进行分段处理，去除坏段；
步骤A7、以被试者静息态下的各个数据通道的脑磁图数据为基准，校正脑磁图数据；
步骤A8、对脑磁图数据进行叠加平均。
可选地，步骤A2中，在脑磁图数据通过截止频率为0.1Hz的高通滤波器后，还通过截止频率为50Hz的陷波滤波器。
可选地，步骤S2中，溯源模型根据正向引导场矩阵和配准构建，将脑磁图数据输入预先构建的溯源定位模型中，根据单偶极子模型或分布式源逆向求解，获得脑磁图数据的信号源信息；正向引导场矩阵是对球模型或单壳模型进行正向问题建模，并利用有限元对大脑及边界进行模拟仿真后得到的；配准是传感器坐标系、大脑坐标系和MRI坐标系之间的配准关系。
可选地，还包括：步骤S4、根据每一种刺激形式下的脑磁图数据，令偏差刺激的波形减去标准刺激的波形，获得每一种刺激形式的差异波，根据每一种刺激形式的差异波的参数信息绘制响应的地形图。
可选地，步骤S4、根据静息态下的脑磁图数据，获得静息态下响应的地形图；
步骤S5、对静息态下响应的地形图进行聚类，确定每一聚类出现的概率。
可选地，还包括：步骤S4、根据每一种刺激形式下的脑磁图数据的频谱，频谱密度，功率谱和功率谱密度，绘制每一种刺激形式下的频域图。
可选地，还包括：步骤S4、根据每一种刺激形式下的脑磁图数据，进行参数检验和非参数检验，获得数据均值和方差；根据数据均值、方差和预设的评价标准，获得评价结果。
（三）有益效果
本发明的有益效果是：
本发明提供的基于脑磁图数据的脑功能分析方法，在滤波处理脑磁图数据的基础上，依次结合利用第一参考数据对脑磁图数据进行回归处理，利用第二参考数据对脑磁图数据进行投影处理，以及主成分分析，能够有效去除生理伪迹和背景噪声，为后续的数据处理提供更为干净的脑磁图数据；并且能够更全面地对脑功能进行分析，包括时域分析、频域分析、微状态分析、网络分析和统计分析，其中在网络分析中基于贝叶斯的多元自回归模型能够构建出更精确的网络结构。
附图说明
本发明借助于以下附图进行描述：
图1为根据本发明一个实施例的脑磁分析装置的结构示意图；
图2为根据本发明一个实施例的基于脑磁图数据的脑功能分析方法的流程示意图；
图3为构建的传感器级脑网络示意图；
图4为构建的信号源级脑网络示意图。
【附图标记说明】
1：屏蔽房；
2：头戴式脑磁图阵列传感器***；
3：感官刺激***；
4：数据采集***；
5：数据处理***。
具体实施方式
为了更好的解释本发明，以便于理解，下面结合附图，通过具体实施方式，对本发明作详细描述。
存在如图1所示的脑磁分析装置，该装置包括磁屏蔽***、头戴式脑磁图阵列传感器***2、感官刺激***3、数据采集***4和数据处理***5。
磁屏蔽***包括用于屏蔽地磁环境的被动屏蔽子***以及用于抵消屏蔽房1内静态剩磁和梯度磁场的主动屏蔽子***；被动屏蔽子***采用多层高磁导率坡莫合金板构筑屏蔽房1墙壁，主动屏蔽子***包括布置在屏蔽房1内被试者左右两侧的赫姆霍兹线圈，与线圈相连的线圈驱动以及用于磁屏蔽室的SERF磁强计构成的参考磁传感阵列，根据参考磁传感阵列提供的背景磁场信息通过线圈驱动控制补偿线圈内电流大小以响应背景磁场。
头戴式脑磁图阵列传感器***2包括无磁脑电采集帽、用于被试者的SERF磁强计构成的磁传感阵列、用来安放磁传感阵列的脑磁采集柔性单元以及磁传感阵列固定底座，无磁脑电采集帽使用单极电极记录刺激呈现时被试者的大脑神经元电活动并使用双极电极记录生物电伪影，无磁脑电采集帽部署在被试头部，头戴式脑磁图阵列传感器***2与下一级数据采集***4相连。
感官刺激***3包括包括生成感官刺激的计算机以及感官刺激的配套呈现设备；生成感官刺激的计算机作为感官刺激的激励源能够在线执行刺激范式程序，生成高保真诊断实验所需刺激；感官刺激的配套呈现设备作为刺激传递中介将刺激有效地呈现给被试者；感官刺激***3与下一级数据采集***4相连。
数据采集***4用于采集被试者的脑磁图数据和感官刺激信号，将感官刺激信号标记在脑磁图数据上。
数据处理***5用于对数据采集***4的数据进行处理，包括时域分析和网络分析。
基于上述脑磁分析装置，本发明提出一种基于脑磁图数据的脑功能分析方法，能够根据脑磁图数据更完善、更精确地进行脑功能分析。
如图2所示，本发明提出的基于脑磁图数据的脑功能分析方法，包括以下步骤：
步骤S1、获取被试者的各个数据通道的脑磁图数据，以及获取第一参考磁图数据和第二参考磁图数据。
其中，脑磁图数据基于用于被试者的SERF磁强计采集；第一参考磁图数据基于用于磁屏蔽室的SERF磁强计采集；第二参考磁图数据是在无被试者时基于用于被试者的SERF磁强计采集。
步骤S2、对脑磁图数据进行预处理。
其中预处理的步骤依次包括：
步骤A1、去除损坏的数据通道的脑磁图数据。
作为一个示例，去除损坏的数据通道的脑磁图数据包括：去除没有采集到数据的数据通道，和去除超出预设幅值的数据通道的脑磁图数据。
步骤A2、令脑磁图数据先通过截止频率为40Hz的低通滤波器，再通过截止频率为0.1Hz的高通滤波器。
优选地，在脑磁图数据通过截止频率为0.1Hz的高通滤波器后，还通过截止频率为50Hz的陷波滤波器。去除掉50Hz的工频噪声。
步骤A3、根据第一参考数据对脑磁图数据进行回归处理，去除背景噪声。
步骤A4、基于信号空间投影法，将第二参考数据在脑磁图数据上进行投影，去除背景噪声和生理伪迹。
步骤A5、通过自适应阈值算法得到生理伪迹的时间窗，对时间窗内数据进行主成分分析，从得到的主成分中选择生理伪迹相关的成分，在脑磁图数据中将生理伪迹相关的成分去除。
步骤A6、对脑磁图数据进行分段处理，去除坏段。
步骤A7、以被试者静息态下的各个数据通道的脑磁图数据为基准，校正脑磁图数据。其中，静息态是指在没有任何刺激的情形下被试者产生的脑磁图数据。
步骤A8、对脑磁图数据进行叠加平均。如此，能够去除随机伪迹。
脑磁图数据的噪声来源包括生理伪迹（如眼电，眼动和肌电产生的信号）和背景噪声，本发明在滤波的基础上，依次结合利用第一参考数据对脑磁图数据进行回归处理，利用第二参考数据对脑磁图数据进行投影处理，以及主成分分析，能够有效去除生理伪迹和背景噪声，为后续的数据处理提供更为干净的脑磁图数据。
步骤S3、将预处理的脑磁图数据输入预先构建的溯源定位模型，获得脑磁图数据的信号源信息。
具体地，步骤S3中，溯源模型根据正向引导场矩阵和配准构建，将脑磁图数据输入预先构建的溯源定位模型中，根据单偶极子模型或分布式源逆向求解，获得脑磁图数据的信号源信息；正向引导场矩阵是对球模型或单壳模型进行正向问题建模，并利用有限元对大脑及边界进行模拟仿真后得到的；配准是传感器坐标系、大脑坐标系和MRI坐标系之间的配准关系。其中，信号源信息包括信号源的位置信息、方向信息和强度信息等。
由于本发明提出的溯源模型的解不是唯一的，而且求得的解是数值解而非解析解，因此可以得到相应的误差函数，从而得到更精确的解，即能够得到更精确的脑磁图数据的信号源信息。利用脑磁图数据的信号源信息后续可以结合功能核磁共振成像fMRI进行相应位置的显示，和进行相应的基于源的统计分析和脑网络分析，进而得到更多信息，从而有利于脑功能分析。
步骤S4、对脑磁图数据进行时域分析：根据每一种刺激形式下的脑磁图数据，令偏差刺激的波形减去标准刺激的波形，获得每一种刺激形式的差异波，根据每一种刺激形式的差异波的参数信息绘制响应的地形图。
其中，每一种刺激形式包括预先定义的依时间顺序排列的标准刺激和偏差刺激；差异波的参数信息包括峰值幅度能量和潜伏期等。
步骤S5、对脑磁图数据进行微状态分析：从步骤S4中获得静息态下响应的地形图，对静息态下响应的地形图进行聚类，确定每一聚类出现的概率。
脑磁图数据的微状态是多数据通道脑磁图数据中地形图拓扑结构的准稳定时期；静息态脑磁图数据由少数交替的微状态所主导；各种神经精神类疾病选择性地影响脑磁图数据微状态；脑磁图数据微状态可反映特定fMRI探测的静息态网络（RSNs）。脑磁图数据微状态是一种有前景的神经生理学工具，可用于了解和评估健康、患病群体在毫秒时间尺度上的脑网络动力学。
本发明提出的方法能够基于脑磁图数据进行微状态分析，更全面地进行脑功能分析。
步骤S6、对脑磁图数据进行频域分析：根据每一种刺激形式下的脑磁图数据的频谱，频谱密度，功率谱和功率谱密度，绘制每一种刺激形式下的频域图。
如此，能够找出比较活跃的频率段，以在频域上进一步进行统计分析和脑网络分析。
步骤S7、对脑磁图数据进行网络分析：将脑磁图数据输入基于贝叶斯的多元自回归模型中，根据期望最大化算法进行模型求解，获得模型系数，基于模型系数进行显著性检验连接数据通道节点，构建传感器级脑网络，如图3所示；和/或，
将信号源信息和脑磁图数据输入基于贝叶斯的多元自回归模型中，根据期望最大化算法进行模型求解，获得模型系数，基于模型系数进行显著性检验连接信号源节点，构建信号源级脑网络，如图4所示。
具体地，基于贝叶斯的多元自回归模型可表示为：
Figure 25755DEST_PATH_IMAGE001
Figure 740770DEST_PATH_IMAGE002
Figure 48255DEST_PATH_IMAGE003
Figure 532326DEST_PATH_IMAGE004
Figure 424321DEST_PATH_IMAGE005
其中，
Figure 236419DEST_PATH_IMAGE006
为模型系数；
Figure 737807DEST_PATH_IMAGE007
为矩阵
Figure 76385DEST_PATH_IMAGE008
的转置；
Figure 372237DEST_PATH_IMAGE009
Figure 671631DEST_PATH_IMAGE010
，为N维脑磁图数据的联合平稳时间序列，可以包含有脑磁图数据的信号源信息；q为 通道数。
如此，能够更精确地对脑网络进行分析。
基于贝叶斯的多元自回归模型BA-MVAR的构建过程，参见以下：
假设向量
Figure 947017DEST_PATH_IMAGE011
是N维脑磁数据的联合平稳时间 序列，假设阶数为q (数据通道数为q)，则多元自回归模型能够表示为以下形式：
Figure 405680DEST_PATH_IMAGE012
（1）
其中，q代表模型延迟观测的最大值，
Figure 872434DEST_PATH_IMAGE013
表示延迟为r时的系数矩阵，它 的形式如下：
Figure 659124DEST_PATH_IMAGE014
（2）
其中
Figure 236736DEST_PATH_IMAGE015
表示
Figure 51370DEST_PATH_IMAGE016
对
Figure 689025DEST_PATH_IMAGE017
的线性影响。非零系数值可 以看做时间序列j对时间序列i的影响。
Figure 963012DEST_PATH_IMAGE018
表示零均值、协方差矩 阵为o2的多元高斯白噪声。
通常，多元自回归模型在估计系数时使用了最小二乘法。基于这个方法，可以把公式（1）写成：
Figure 609894DEST_PATH_IMAGE019
（3）
其中,公式右侧表示L2模，
Figure 777570DEST_PATH_IMAGE020
是
Figure 461492DEST_PATH_IMAGE021
的向量；
Figure 583294DEST_PATH_IMAGE022
是待估计的系数，定义
Figure 502708DEST_PATH_IMAGE023
Figure 790470DEST_PATH_IMAGE024
；
Figure 645294DEST_PATH_IMAGE021
是延迟矩阵，其形式为：
Figure 18506DEST_PATH_IMAGE025
（4）
对（3）式进行化简，当
Figure 243077DEST_PATH_IMAGE026
，我们有以下的结果
Figure 119766DEST_PATH_IMAGE027
（5）
则模型系数能够表示为
Figure 676649DEST_PATH_IMAGE028
（6）
估计出多元自回归模型MVAR 系数后，网络构建将基于式（6）的计算结果，因此MVAR系数的精确估计对于得到可靠的因果连接网络具有决定性作用。然而，在诸如功能磁共振成像这样的实际应用中，实验中采集的信号数据不可避免地会包含异常值的出现。基于最小二乘法的多元自回归分析（LS-MVAR）在这样的实际应用中有着巨大的缺陷：噪声的影响会被放大，这一点从式（3）中的二次方项容易看出。这会导致估计出的MVAR系数不准确从而进一步影响后续网络构建，产生更多伪因故连接，扭曲网络拓扑结构。而且，基于最小二乘法的MVAR需要比较长的数据才能捕捉时间序列的动态变化。为了改善MVAR系数的精确程度，本发明采用贝叶斯方法。具体如下：
假设未知系数
Figure 271578DEST_PATH_IMAGE006
服从零均值的独立高斯分布
Figure 267216DEST_PATH_IMAGE029
，则有
Figure 765456DEST_PATH_IMAGE030
（7）
其中n是
Figure 227661DEST_PATH_IMAGE006
的长度
由于高斯白噪声
Figure 309887DEST_PATH_IMAGE031
，其中
Figure 640374DEST_PATH_IMAGE032
，则
Figure 367022DEST_PATH_IMAGE033
，有
Figure 859183DEST_PATH_IMAGE034
（8）
贝叶斯公式用来描述两个概率之间的关系，形式如下：
Figure 930169DEST_PATH_IMAGE035
（9）
则后验概率
Figure 798768DEST_PATH_IMAGE036
与
Figure 379922DEST_PATH_IMAGE037
成正比，由此可以得到最大后验 概率表达式
Figure 308564DEST_PATH_IMAGE038
（10）
对该函数取对数后得到目标函数
Figure 99802DEST_PATH_IMAGE039
（11）
为求得目标函数最大值，我们取
Figure 913038DEST_PATH_IMAGE040
，从而有
Figure 240376DEST_PATH_IMAGE041
（12）
Figure 808761DEST_PATH_IMAGE042
（13）
Figure 618454DEST_PATH_IMAGE043
（14）
在本发明的基于贝叶斯的多元自回归模型中，未知变量是
Figure 704221DEST_PATH_IMAGE044
，可以通过期 望最大化EM算法迭代求解。未知变量中他们能够彼此估计出对方，具体地说：当
Figure 384601DEST_PATH_IMAGE006
已知，就 能得到
Figure 890931DEST_PATH_IMAGE045
，当
Figure 656762DEST_PATH_IMAGE046
已知， 也就能得到
Figure 811800DEST_PATH_IMAGE006
。EM算法会先给出
Figure 346686DEST_PATH_IMAGE045
的初值，然后有此得 到
Figure 522453DEST_PATH_IMAGE006
；然后由这个
Figure 916525DEST_PATH_IMAGE006
值再估计出
Figure 970194DEST_PATH_IMAGE046
。重复这个过程，直到式（11）的目标函数达到收敛。
EM算法的步骤如下：
（1）E步：计算期望。根据未知量的目前已经估计的值，计算式（12）中的最大似然值；
（2）M步：最大化。最大化E步中得到的最大似然值，得到
Figure 625166DEST_PATH_IMAGE046
的值，换句话说，就是 计算式（13）和式（14）
（3）重复步骤（1）（2），直到式子（11）收敛，则我们就得到未知量
Figure 237413DEST_PATH_IMAGE044
的最终 估计值。
基于BA-MVAR估计出的模型系数，可通过显著性检验去构造因果网络。首先，估计 出每个时间延迟中的系数矩阵的方差和均值，将
Figure 587623DEST_PATH_IMAGE047
的值转化为学生T分布的
Figure 943518DEST_PATH_IMAGE048
；第二，计 算P=0.05的学生T分布式累积概率密度函数（TCDF）的值，就将其置为0，最后，对所有事件延 迟上的
Figure 931024DEST_PATH_IMAGE048
求和。如果总的
Figure 714172DEST_PATH_IMAGE048
值大于0，这表明这两个节点之前的连接是显著的。
步骤S8、对脑磁图数据进行统计分析：根据每一种刺激形式下的脑磁图数据，进行参数检验和非参数检验，获得数据均值和方差；根据数据均值、方差和预设的评价标准，获得评价结果。
其中，参数检验可以是T检验和方差分析，非参数检验可以是置换检验。
综上，本发明提供的基于脑磁图数据的脑功能分析方法，在滤波处理脑磁图数据的基础上，依次结合利用第一参考数据对脑磁图数据进行回归处理，利用第二参考数据对脑磁图数据进行投影处理，以及主成分分析，能够有效去除生理伪迹和背景噪声，为后续的数据处理提供更为干净的脑磁图数据；并且能够更全面地对脑功能进行分析，包括时域分析、频域分析、微状态分析、网络分析和统计分析，其中在网络分析中基于贝叶斯的多元自回归模型能够构建出更精确的网络结构。
需要理解的是，以上对本发明的具体实施例进行的描述只是为了说明本发明的技术路线和特点，其目的在于让本领域内的技术人员能够了解本发明的内容并据以实施，但本发明并不限于上述特定实施方式。凡是在本发明权利要求的范围内做出的各种变化或修饰，都应涵盖在本发明的保护范围内。

Claims (8)

1.一种基于脑磁图数据的脑功能分析方法，其特征在于，包括以下步骤：
步骤S1、获取被试者的各个数据通道的脑磁图数据；所述脑磁图数据基于用于被试者的SERF磁强计采集；
步骤S2、将所述脑磁图数据输入预先构建的溯源定位模型，获得脑磁图数据的信号源信息；
步骤S3、将所述脑磁图数据输入基于贝叶斯的多元自回归模型中，根据期望最大化算法进行模型求解，获得模型系数，基于模型系数进行显著性检验连接数据通道节点，构建传感器级脑网络；和，
将所述信号源信息和所述脑磁图数据输入基于贝叶斯的多元自回归模型中，根据期望最大化算法进行模型求解，获得模型系数，基于模型系数进行显著性检验连接信号源节点，构建信号源级脑网络；
所述基于贝叶斯的多元自回归模型可表示为：
Figure 637538DEST_PATH_IMAGE001
Figure 57018DEST_PATH_IMAGE002
Figure 267419DEST_PATH_IMAGE003
Figure 947799DEST_PATH_IMAGE004
Figure 719709DEST_PATH_IMAGE005
其中，
Figure 360905DEST_PATH_IMAGE006
为模型系数；
Figure 109419DEST_PATH_IMAGE007
为矩阵
Figure 909884DEST_PATH_IMAGE008
的转置；
Figure 85651DEST_PATH_IMAGE009
Figure 214144DEST_PATH_IMAGE010
，为N 维脑磁图数据的联合平稳时间序列，当构建传感器级脑网络时，x(k)包括脑磁图数据，当构 建信号源级脑网络时，x(k)包括脑磁图数据和脑磁图数据的信号源信息；q为通道数。
2.根据权利要求1所述的基于脑磁图数据的脑功能分析方法，其特征在于，步骤S1中，所述脑磁图数据为被试者处于磁屏蔽室中的脑磁图数据，步骤S1还包括：获取第一参考磁图数据和第二参考磁图数据，所述第一参考磁图数据基于用于磁屏蔽室的SERF磁强计采集，所述第二参考磁图数据是在无被试者时基于用于被试者的SERF磁强计采集；
在步骤S1和步骤S2之间，还包括：步骤S12、对所述脑磁图数据进行预处理；所述预处理依次包括以下步骤：
步骤A1、去除损坏的数据通道的脑磁图数据；
步骤A2、令脑磁图数据先通过截止频率为40Hz的低通滤波器，再通过截止频率为0.1Hz的高通滤波器；
步骤A3、根据第一参考数据对脑磁图数据进行回归处理，去除背景噪声；
步骤A4、基于信号空间投影法，将第二参考数据在脑磁图数据上进行投影，去除背景噪声和生理伪迹；
步骤A5、通过自适应阈值算法得到生理伪迹的时间窗，对时间窗内数据进行主成分分析，从得到的主成分中选择生理伪迹相关的成分，在脑磁图数据中将生理伪迹相关的成分去除；
步骤A6、对脑磁图数据进行分段处理，去除坏段；
步骤A7、以被试者静息态下的各个数据通道的脑磁图数据为基准，校正脑磁图数据；
步骤A8、对脑磁图数据进行叠加平均。
3.根据权利要求2所述的基于脑磁图数据的脑功能分析方法，其特征在于，步骤A2中，
在脑磁图数据通过截止频率为0.1Hz的高通滤波器后，还通过截止频率为50Hz的陷波滤波器。
4.根据权利要求1所述的基于脑磁图数据的脑功能分析方法，其特征在于，步骤S2中，
所述溯源模型根据正向引导场矩阵和配准构建，将所述脑磁图数据输入预先构建的溯源定位模型中，根据单偶极子模型或分布式源逆向求解，获得脑磁图数据的信号源信息；
所述正向引导场矩阵是对球模型或单壳模型进行正向问题建模，并利用有限元对大脑及边界进行模拟仿真后得到的；所述配准是传感器坐标系、大脑坐标系和MRI坐标系之间的配准关系。
5.根据权利要求1所述的基于脑磁图数据的脑功能分析方法，其特征在于，还包括：
步骤S4、根据每一种刺激形式下的脑磁图数据，令偏差刺激的波形减去标准刺激的波形，获得每一种刺激形式的差异波，根据每一种刺激形式的差异波的参数信息绘制响应的地形图。
6.根据权利要求5所述的基于脑磁图数据的脑功能分析方法，其特征在于，还包括：
步骤S5、从步骤S4中获得静息态下响应的地形图，对静息态下响应的地形图进行聚类，确定每一聚类出现的概率。
7.根据权利要求1所述的基于脑磁图数据的脑功能分析方法，其特征在于，还包括：
步骤S4、根据每一种刺激形式下的脑磁图数据的频谱，频谱密度，功率谱和功率谱密度，绘制每一种刺激形式下的频域图。
8.根据权利要求1所述的基于脑磁图数据的脑功能分析方法，其特征在于，还包括：
步骤S4、根据每一种刺激形式下的脑磁图数据，进行参数检验和非参数检验，获得数据均值和方差；根据数据均值、方差和预设的评价标准，获得评价结果。
CN202210046874.5A 2022-01-17 2022-01-17 一种基于脑磁图数据的脑功能分析方法 Active CN114052668B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210046874.5A CN114052668B (zh) 2022-01-17 2022-01-17 一种基于脑磁图数据的脑功能分析方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210046874.5A CN114052668B (zh) 2022-01-17 2022-01-17 一种基于脑磁图数据的脑功能分析方法

Publications (2)

Publication Number Publication Date
CN114052668A CN114052668A (zh) 2022-02-18
CN114052668B true CN114052668B (zh) 2022-06-17

Family

ID=80231173

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210046874.5A Active CN114052668B (zh) 2022-01-17 2022-01-17 一种基于脑磁图数据的脑功能分析方法

Country Status (1)

Country Link
CN (1) CN114052668B (zh)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2024023268A1 (en) * 2022-07-29 2024-02-01 F. Hoffmann-La Roche Ag Characterisation of neurological dysfunction
CN115018018B (zh) * 2022-08-05 2022-11-11 北京航空航天大学杭州创新研究院 一种抑制背景噪声的双重空间滤波方法
CN115359236B (zh) * 2022-10-19 2023-03-28 之江实验室 一种基于机器学习的工作记忆任务脑磁图分类***
CN116863025B (zh) * 2023-09-05 2023-12-05 之江实验室 脑磁图数据的溯源重建方法、装置、电子装置和介质

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113827246A (zh) * 2021-11-25 2021-12-24 北京航空航天大学杭州创新研究院 脑磁数据采集分析方法及***

Family Cites Families (19)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CA2649906A1 (en) * 2006-04-21 2007-11-01 Regents Of The University Of Minnesota Activation and causal interaction of bioelectrical activity
CN102360501A (zh) * 2011-07-19 2012-02-22 中国科学院自动化研究所 一种基于核磁共振成像的脑区有效连接度构建方法
CN102509282B (zh) * 2011-09-26 2014-06-25 东南大学 一种融合结构连接的各脑区间的效能连接分析方法
US20180279939A1 (en) * 2015-06-09 2018-10-04 The Children's Medical Center Corporation Method and system for locating seizure focus from interictal data
CN105147288B (zh) * 2015-07-23 2018-06-08 中国科学院苏州生物医学工程技术研究所 脑磁源强度定位方法
CN105249964B (zh) * 2015-11-10 2017-12-22 东南大学 基于脑磁图和弥散张量成像的多模态脑功能重建评估方法
WO2018066715A1 (ja) * 2016-10-07 2018-04-12 国立大学法人筑波大学 脳波検出装置およびプログラム
JP6812022B2 (ja) * 2018-10-11 2021-01-13 株式会社国際電気通信基礎技術研究所 脳機能結合相関値の調整方法、脳機能結合相関値の調整システム、脳活動分類器のハーモナイズ方法、脳活動分類器のハーモナイズシステム、および脳活動バイオマーカシステム
WO2020075737A1 (ja) * 2018-10-11 2020-04-16 株式会社国際電気通信基礎技術研究所 脳機能結合相関値の調整方法、脳機能結合相関値の調整システム、脳活動分類器のハーモナイズ方法、脳活動分類器のハーモナイズシステム、及び脳活動バイオマーカシステム
CN110090017B (zh) * 2019-03-11 2021-09-14 北京工业大学 一种基于lstm的脑电信号源定位方法
CN111096745B (zh) * 2020-01-02 2021-03-30 苏州大学 基于稀疏贝叶斯学习的稳态诱发响应脑源定位方法
CN111543949B (zh) * 2020-05-13 2021-09-10 北京航空航天大学 一种基于脑磁图与脑电图的儿童asd诊断装置
CN111340142B (zh) * 2020-05-14 2020-08-14 南京慧脑云计算有限公司 一种癫痫脑磁图棘波自动检测方法与溯源定位***
CN112294339B (zh) * 2020-10-14 2022-12-09 中国科学院苏州生物医学工程技术研究所 基于种群多样性控制的脑电源定位方法、***及设备
CN112914578B (zh) * 2021-01-20 2024-02-09 季华实验室 Meg源定位方法及***
CN113143247A (zh) * 2021-04-29 2021-07-23 常州大学 一种大脑功能超网络的构建方法
CN113274037B (zh) * 2021-06-30 2022-08-26 中国科学院苏州生物医学工程技术研究所 一种动态脑功能网络的生成方法、***及设备
CN113576491A (zh) * 2021-07-26 2021-11-02 深圳市人民医院 基于静息态eeg频域特征及脑网络自动分析方法和***
CN113768464A (zh) * 2021-08-20 2021-12-10 西安交通大学 一种基于经验模态分解的癫痫发作时间段检测***及方法

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113827246A (zh) * 2021-11-25 2021-12-24 北京航空航天大学杭州创新研究院 脑磁数据采集分析方法及***

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
用于构建脑磁图网络的信号提取方法;杨春兰等;《北京工业大学学报》;20200710(第07期);全文 *
脑磁图脑功能连接网络癫痫棘波识别方法研究;张航宇等;《计算机工程与应用》(第08期);全文 *

Also Published As

Publication number Publication date
CN114052668A (zh) 2022-02-18

Similar Documents

Publication Publication Date Title
CN114052668B (zh) 一种基于脑磁图数据的脑功能分析方法
Wu et al. Bayesian Machine Learning: EEG\/MEG signal processing measurements
De Munck et al. Estimating stationary dipoles from MEG/EEG data contaminated with spatially and temporally correlated background noise
JP6884344B2 (ja) 脳内ネットワークの活動推定システム、脳内ネットワークの活動推定方法、脳内ネットワークの活動推定プログラム、および、学習済み脳活動推定モデル
US20160051162A1 (en) Method for locating a brain activity associated with a task
US20120232376A1 (en) Methods and systems for channel selection
Barton et al. Evaluating the performance of kalman-filter-based eeg source localization
US10945627B2 (en) Systems and methods for super-resolving electromagnetic localization and causality
Hansen et al. Unmixing oscillatory brain activity by EEG source localization and empirical mode decomposition
Jatoi et al. Brain source localization using reduced EEG sensors
Hecker et al. Long-short term memory networks for electric source imaging with distributed dipole models
Chang et al. Assessing recurrent interactions in cortical networks: Modeling EEG response to transcranial magnetic stimulation
Aram et al. Model-based estimation of intra-cortical connectivity using electrophysiological data
US20210186400A1 (en) Apparatuses, systems, and methods for suppression of artifacts in non-invasive electromagnetic recordings
Sun et al. SIFNet: Electromagnetic source imaging framework using deep neural networks
Dinh et al. Contextual minimum-norm estimates (CMNE): a deep learning method for source estimation in neuronal networks
Baillet Electromagnetic brain mapping using MEG and EEG
JP2021087761A (ja) 信号分離装置、プログラムおよび信号分離方法
Pflieger et al. Nonlinear analysis of multimodal dynamic brain imaging data
Treder et al. Source reconstruction of broadband EEG/MEG data using the frequency-adaptive broadband (FAB) beamformer
Wang et al. OPM-MEG bad channel identification method based on the improved box-isolation forest algorithm
Pascarella et al. Particle filtering, beamforming and multiple signal classification for the analysis of magnetoencephalography time series: a comparison of algorithms
Zaitcev EEG source imaging for improved control BCI performance
López et al. Spatio-temporal EEG brain imaging based on reduced Kalman filtering
Pankajbhai et al. Medical Signal Analysis Using Artificial Intelligence

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