1.本发明涉及脑功能网络分析技术领域，尤其涉及一种基于脑磁图数据的脑功能分析方法。


背景技术：

2.人脑作为世界上中最为复杂的动力学系统之一，其皮层由150-330亿个神经细胞组成，这些神经元间的相互连接构成了一个拥有相当复杂结构和功能的脑网络，在神经学和生理学方面，大脑的主要功能是控制和支配身体的各个器官，因此，这个复杂而庞大的脑网络使它具有高级的信息处理和认知表达功能，如：语言、情感、记忆、认知等，并对来自人体内部和周围环境的信息进行存储、处理、加工和整合，探究大脑的结构和功能，加速脑科学领域的研究不仅可以提高对大脑疾病的预防、诊断和治疗，更可以推动人工智能的发展。
3.脑磁图能够以较高的时空分辨率对大脑神经活动进行直接成像，为大脑连通分析提供详细的时空以及特定节律活动的信息，定位精细度高，能有效地捕捉深部源放电过程。为此，本发明提供一种基于脑磁图数据的脑功能分析方法。


技术实现要素：

4.（一）要解决的技术问题鉴于上述技术中存在的问题，本发明至少从一定程度上进行解决。为此，本发明的目的在于提出了一种基于脑磁图数据的脑功能分析方法，能够构建出更精确的脑功能网络。
5.（二）技术方案为了达到上述目的，本发明提供一种基于脑磁图数据的脑功能分析方法，包括以下步骤：步骤s1、获取被试者的各个数据通道的脑磁图数据；脑磁图数据基于用于被试者的serf磁强计采集；步骤s2、将脑磁图数据输入预先构建的溯源定位模型，获得脑磁图数据的信号源信息；步骤s3、将脑磁图数据输入基于贝叶斯的多元自回归模型中，根据期望最大化算法进行模型求解，获得模型系数，基于模型系数进行显著性检验连接数据通道节点，构建传感器级脑网络；和/或，将信号源信息和脑磁图数据输入基于贝叶斯的多元自回归模型中，根据期望最大化算法进行模型求解，获得模型系数，基于模型系数进行显著性检验连接信号源节点，构建信号源级脑网络。
6.可选地，基于贝叶斯的多元自回归模型可表示为：
其中，为模型系数；为矩阵的转置；的转置；，为n维脑磁图数据的联合平稳时间序列，可以包含有脑磁图数据的信号源信息；q为通道数。
7.可选地，步骤s1中，脑磁图数据为被试者处于磁屏蔽室中的脑磁图数据，步骤s1还包括：获取第一参考磁图数据和第二参考磁图数据，第一参考磁图数据基于用于磁屏蔽室的serf磁强计采集，第二参考磁图数据是在无被试者时基于用于被试者的serf磁强计采集；在步骤s1和步骤s2之间，还包括：步骤s12、对脑磁图数据进行预处理；预处理依次包括以下步骤：步骤a1、去除损坏的数据通道的脑磁图数据；步骤a2、令脑磁图数据先通过截止频率为40hz的低通滤波器，再通过截止频率为0.1hz的高通滤波器；步骤a3、根据第一参考数据对脑磁图数据进行回归处理，去除背景噪声；步骤a4、基于信号空间投影法，将第二参考数据在脑磁图数据上进行投影，去除背景噪声和生理伪迹；步骤a5、通过自适应阈值算法得到生理伪迹的时间窗，对时间窗内数据进行主成分分析，从得到的主成分中选择生理伪迹相关的成分，在脑磁图数据中将生理伪迹相关的成分去除；步骤a6、对脑磁图数据进行分段处理，去除坏段；步骤a7、以被试者静息态下的各个数据通道的脑磁图数据为基准，校正脑磁图数据；步骤a8、对脑磁图数据进行叠加平均。
8.可选地，步骤a2中，在脑磁图数据通过截止频率为0.1hz的高通滤波器后，还通过截止频率为50hz的陷波滤波器。
9.可选地，步骤s2中，溯源模型根据正向引导场矩阵和配准构建，将脑磁图数据输入预先构建的溯源定位模型中，根据单偶极子模型或分布式源逆向求解，获得脑磁图数据的信号源信息；正向引导场矩阵是对球模型或单壳模型进行正向问题建模，并利用有限元对大脑及边界进行模拟仿真后得到的；配准是传感器坐标系、大脑坐标系和mri坐标系之间的配准关系。
10.可选地，还包括：步骤s4、根据每一种刺激形式下的脑磁图数据，令偏差刺激的波
形减去标准刺激的波形，获得每一种刺激形式的差异波，根据每一种刺激形式的差异波的参数信息绘制响应的地形图。
11.可选地，步骤s4、根据静息态下的脑磁图数据，获得静息态下响应的地形图；步骤s5、对静息态下响应的地形图进行聚类，确定每一聚类出现的概率。
12.可选地，还包括：步骤s4、根据每一种刺激形式下的脑磁图数据的频谱，频谱密度，功率谱和功率谱密度，绘制每一种刺激形式下的频域图。
13.可选地，还包括：步骤s4、根据每一种刺激形式下的脑磁图数据，进行参数检验和非参数检验，获得数据均值和方差；根据数据均值、方差和预设的评价标准，获得评价结果。
14.（三）有益效果本发明的有益效果是：本发明提供的基于脑磁图数据的脑功能分析方法，在滤波处理脑磁图数据的基础上，依次结合利用第一参考数据对脑磁图数据进行回归处理，利用第二参考数据对脑磁图数据进行投影处理，以及主成分分析，能够有效去除生理伪迹和背景噪声，为后续的数据处理提供更为干净的脑磁图数据；并且能够更全面地对脑功能进行分析，包括时域分析、频域分析、微状态分析、网络分析和统计分析，其中在网络分析中基于贝叶斯的多元自回归模型能够构建出更精确的网络结构。
附图说明
15.本发明借助于以下附图进行描述：图1为根据本发明一个实施例的脑磁分析装置的结构示意图；图2为根据本发明一个实施例的基于脑磁图数据的脑功能分析方法的流程示意图；图3为构建的传感器级脑网络示意图；图4为构建的信号源级脑网络示意图。
16.【附图标记说明】1：屏蔽房；2：头戴式脑磁图阵列传感器系统；3：感官刺激系统；4：数据采集系统；5：数据处理系统。
具体实施方式
17.为了更好的解释本发明，以便于理解，下面结合附图，通过具体实施方式，对本发明作详细描述。
18.存在如图1所示的脑磁分析装置，该装置包括磁屏蔽系统、头戴式脑磁图阵列传感器系统2、感官刺激系统3、数据采集系统4和数据处理系统5。
19.磁屏蔽系统包括用于屏蔽地磁环境的被动屏蔽子系统以及用于抵消屏蔽房1内静态剩磁和梯度磁场的主动屏蔽子系统；被动屏蔽子系统采用多层高磁导率坡莫合金板构筑屏蔽房1墙壁，主动屏蔽子系统包括布置在屏蔽房1内被试者左右两侧的赫姆霍兹线圈，与
线圈相连的线圈驱动以及用于磁屏蔽室的serf磁强计构成的参考磁传感阵列，根据参考磁传感阵列提供的背景磁场信息通过线圈驱动控制补偿线圈内电流大小以响应背景磁场。
20.头戴式脑磁图阵列传感器系统2包括无磁脑电采集帽、用于被试者的serf磁强计构成的磁传感阵列、用来安放磁传感阵列的脑磁采集柔性单元以及磁传感阵列固定底座，无磁脑电采集帽使用单极电极记录刺激呈现时被试者的大脑神经元电活动并使用双极电极记录生物电伪影，无磁脑电采集帽部署在被试头部，头戴式脑磁图阵列传感器系统2与下一级数据采集系统4相连。
21.感官刺激系统3包括包括生成感官刺激的计算机以及感官刺激的配套呈现设备；生成感官刺激的计算机作为感官刺激的激励源能够在线执行刺激范式程序，生成高保真诊断实验所需刺激；感官刺激的配套呈现设备作为刺激传递中介将刺激有效地呈现给被试者；感官刺激系统3与下一级数据采集系统4相连。
22.数据采集系统4用于采集被试者的脑磁图数据和感官刺激信号，将感官刺激信号标记在脑磁图数据上。
23.数据处理系统5用于对数据采集系统4的数据进行处理，包括时域分析和网络分析。
24.基于上述脑磁分析装置，本发明提出一种基于脑磁图数据的脑功能分析方法，能够根据脑磁图数据更完善、更精确地进行脑功能分析。
25.如图2所示，本发明提出的基于脑磁图数据的脑功能分析方法，包括以下步骤：步骤s1、获取被试者的各个数据通道的脑磁图数据，以及获取第一参考磁图数据和第二参考磁图数据。
26.其中，脑磁图数据基于用于被试者的serf磁强计采集；第一参考磁图数据基于用于磁屏蔽室的serf磁强计采集；第二参考磁图数据是在无被试者时基于用于被试者的serf磁强计采集。
27.步骤s2、对脑磁图数据进行预处理。
28.其中预处理的步骤依次包括：步骤a1、去除损坏的数据通道的脑磁图数据。
29.作为一个示例，去除损坏的数据通道的脑磁图数据包括：去除没有采集到数据的数据通道，和去除超出预设幅值的数据通道的脑磁图数据。
30.步骤a2、令脑磁图数据先通过截止频率为40hz的低通滤波器，再通过截止频率为0.1hz的高通滤波器。
31.优选地，在脑磁图数据通过截止频率为0.1hz的高通滤波器后，还通过截止频率为50hz的陷波滤波器。去除掉50hz的工频噪声。
32.步骤a3、根据第一参考数据对脑磁图数据进行回归处理，去除背景噪声。
33.步骤a4、基于信号空间投影法，将第二参考数据在脑磁图数据上进行投影，去除背景噪声和生理伪迹。
34.步骤a5、通过自适应阈值算法得到生理伪迹的时间窗，对时间窗内数据进行主成分分析，从得到的主成分中选择生理伪迹相关的成分，在脑磁图数据中将生理伪迹相关的成分去除。
35.步骤a6、对脑磁图数据进行分段处理，去除坏段。
36.步骤a7、以被试者静息态下的各个数据通道的脑磁图数据为基准，校正脑磁图数据。其中，静息态是指在没有任何刺激的情形下被试者产生的脑磁图数据。
37.步骤a8、对脑磁图数据进行叠加平均。如此，能够去除随机伪迹。
38.脑磁图数据的噪声来源包括生理伪迹（如眼电，眼动和肌电产生的信号）和背景噪声，本发明在滤波的基础上，依次结合利用第一参考数据对脑磁图数据进行回归处理，利用第二参考数据对脑磁图数据进行投影处理，以及主成分分析，能够有效去除生理伪迹和背景噪声，为后续的数据处理提供更为干净的脑磁图数据。
39.步骤s3、将预处理的脑磁图数据输入预先构建的溯源定位模型，获得脑磁图数据的信号源信息。
40.具体地，步骤s3中，溯源模型根据正向引导场矩阵和配准构建，将脑磁图数据输入预先构建的溯源定位模型中，根据单偶极子模型或分布式源逆向求解，获得脑磁图数据的信号源信息；正向引导场矩阵是对球模型或单壳模型进行正向问题建模，并利用有限元对大脑及边界进行模拟仿真后得到的；配准是传感器坐标系、大脑坐标系和mri坐标系之间的配准关系。其中，信号源信息包括信号源的位置信息、方向信息和强度信息等。
41.由于本发明提出的溯源模型的解不是唯一的，而且求得的解是数值解而非解析解，因此可以得到相应的误差函数，从而得到更精确的解，即能够得到更精确的脑磁图数据的信号源信息。利用脑磁图数据的信号源信息后续可以结合功能核磁共振成像fmri进行相应位置的显示，和进行相应的基于源的统计分析和脑网络分析，进而得到更多信息，从而有利于脑功能分析。
42.步骤s4、对脑磁图数据进行时域分析：根据每一种刺激形式下的脑磁图数据，令偏差刺激的波形减去标准刺激的波形，获得每一种刺激形式的差异波，根据每一种刺激形式的差异波的参数信息绘制响应的地形图。
43.其中，每一种刺激形式包括预先定义的依时间顺序排列的标准刺激和偏差刺激；差异波的参数信息包括峰值幅度能量和潜伏期等。
44.步骤s5、对脑磁图数据进行微状态分析：从步骤s4中获得静息态下响应的地形图，对静息态下响应的地形图进行聚类，确定每一聚类出现的概率。
45.脑磁图数据的微状态是多数据通道脑磁图数据中地形图拓扑结构的准稳定时期；静息态脑磁图数据由少数交替的微状态所主导；各种神经精神类疾病选择性地影响脑磁图数据微状态；脑磁图数据微状态可反映特定fmri探测的静息态网络（rsns）。脑磁图数据微状态是一种有前景的神经生理学工具，可用于了解和评估健康、患病群体在毫秒时间尺度上的脑网络动力学。
46.本发明提出的方法能够基于脑磁图数据进行微状态分析，更全面地进行脑功能分析。
47.步骤s6、对脑磁图数据进行频域分析：根据每一种刺激形式下的脑磁图数据的频谱，频谱密度，功率谱和功率谱密度，绘制每一种刺激形式下的频域图。
48.如此，能够找出比较活跃的频率段，以在频域上进一步进行统计分析和脑网络分析。
49.步骤s7、对脑磁图数据进行网络分析：将脑磁图数据输入基于贝叶斯的多元自回归模型中，根据期望最大化算法进行模型求解，获得模型系数，基于模型系数进行显著性检
验连接数据通道节点，构建传感器级脑网络，如图3所示；和/或，将信号源信息和脑磁图数据输入基于贝叶斯的多元自回归模型中，根据期望最大化算法进行模型求解，获得模型系数，基于模型系数进行显著性检验连接信号源节点，构建信号源级脑网络，如图4所示。
50.具体地，基于贝叶斯的多元自回归模型可表示为：其中，为模型系数；为矩阵的转置；的转置；，为n维脑磁图数据的联合平稳时间序列，可以包含有脑磁图数据的信号源信息；q为通道数。
51.如此，能够更精确地对脑网络进行分析。
52.基于贝叶斯的多元自回归模型ba-mvar的构建过程，参见以下：假设向量是n维脑磁数据的联合平稳时间序列，假设阶数为q (数据通道数为q)，则多元自回归模型能够表示为以下形式：（1）其中，q代表模型延迟观测的最大值，表示延迟为r时的系数矩阵，它的形式如下：（2）其中表示对的线性影响。非零系数值可以看做时间序列j对时间序列i的影响。表示零均值、协方差矩阵为o2的多元高斯白噪声。
53.通常，多元自回归模型在估计系数时使用了最小二乘法。基于这个方法，可以把公式（1）写成：（3）
其中,公式右侧表示l2模，是的向量；是待估计的系数，定义是待估计的系数，定义；是延迟矩阵，其形式为：（4）对（3）式进行化简，当，我们有以下的结果（5）则模型系数能够表示为（6）估计出多元自回归模型mvar 系数后，网络构建将基于式（6）的计算结果，因此mvar系数的精确估计对于得到可靠的因果连接网络具有决定性作用。然而，在诸如功能磁共振成像这样的实际应用中，实验中采集的信号数据不可避免地会包含异常值的出现。基于最小二乘法的多元自回归分析（ls-mvar）在这样的实际应用中有着巨大的缺陷：噪声的影响会被放大，这一点从式（3）中的二次方项容易看出。这会导致估计出的mvar系数不准确从而进一步影响后续网络构建，产生更多伪因故连接，扭曲网络拓扑结构。而且，基于最小二乘法的mvar需要比较长的数据才能捕捉时间序列的动态变化。为了改善mvar系数的精确程度，本发明采用贝叶斯方法。具体如下：假设未知系数服从零均值的独立高斯分布，则有（7）其中n是的长度由于高斯白噪声，其中，则，有（8）贝叶斯公式用来描述两个概率之间的关系，形式如下：（9）则后验概率与成正比，由此可以得到最大后验概率表达式（10）对该函数取对数后得到目标函数
（11）为求得目标函数最大值，我们取，从而有（12）（13）（14）在本发明的基于贝叶斯的多元自回归模型中，未知变量是，可以通过期望最大化em算法迭代求解。未知变量中他们能够彼此估计出对方，具体地说：当已知，就能得到，当已知， 也就能得到。em算法会先给出的初值，然后有此得到；然后由这个值再估计出。重复这个过程，直到式（11）的目标函数达到收敛。
54.em算法的步骤如下：（1）e步：计算期望。根据未知量的目前已经估计的值，计算式（12）中的最大似然值；（2）m步：最大化。最大化e步中得到的最大似然值，得到的值，换句话说，就是计算式（13）和式（14）（3）重复步骤（1）（2），直到式子（11）收敛，则我们就得到未知量 的最终估计值。
55.基于ba-mvar估计出的模型系数，可通过显著性检验去构造因果网络。首先，估计出每个时间延迟中的系数矩阵的方差和均值，将的值转化为学生t分布的；第二，计算p=0.05的学生t分布式累积概率密度函数（tcdf）的值，就将其置为0，最后，对所有事件延迟上的求和。如果总的值大于0，这表明这两个节点之前的连接是显著的。
56.步骤s8、对脑磁图数据进行统计分析：根据每一种刺激形式下的脑磁图数据，进行参数检验和非参数检验，获得数据均值和方差；根据数据均值、方差和预设的评价标准，获得评价结果。
57.其中，参数检验可以是t检验和方差分析，非参数检验可以是置换检验。
58.综上，本发明提供的基于脑磁图数据的脑功能分析方法，在滤波处理脑磁图数据的基础上，依次结合利用第一参考数据对脑磁图数据进行回归处理，利用第二参考数据对脑磁图数据进行投影处理，以及主成分分析，能够有效去除生理伪迹和背景噪声，为后续的数据处理提供更为干净的脑磁图数据；并且能够更全面地对脑功能进行分析，包括时域分析、频域分析、微状态分析、网络分析和统计分析，其中在网络分析中基于贝叶斯的多元自回归模型能够构建出更精确的网络结构。
59.需要理解的是，以上对本发明的具体实施例进行的描述只是为了说明本发明的技
术路线和特点，其目的在于让本领域内的技术人员能够了解本发明的内容并据以实施，但本发明并不限于上述特定实施方式。凡是在本发明权利要求的范围内做出的各种变化或修饰，都应涵盖在本发明的保护范围内。