基于功能磁共振成像的时变约束脑电图或脑磁图溯源方法与流程

本发明涉及脑电图、脑磁图及磁共振成像领域，尤其涉及一种基于功能磁共振成像的时变约束脑电图或脑磁图溯源方法。

背景技术：

脑电图(electroencephalography，eeg)或脑磁图(magnetoencephalography，meg)有很高的时间分辨率，可以达到毫秒量级，并且他们的信号主要来源于突触后电位，是对大脑活动的直接测量。但是由于脑电图电极数或脑磁图传感器数目远远小于大脑内电流源的数目，因此想要唯一地确定脑神经电流源的位置需要加入外加的约束，而这些外加的约束往往是缺乏生理学依据的。采用脑电图或脑磁图信号研究脑神经活动，确定神经电流源的逆问题算法中具有代表性的如最小二次模算法(minimumnormestimation，mne)、深度加权最小范数法(depth-weightedminimumnorm，dwmn)等，但这些方法始终不能摆脱空间分辨率低的问题。

基于血氧水平依赖(blood-oxygen-level-dependence，bold)的功能磁共振成像(functionalmagneticresonanceimaging，fmri)方法自1992年ogawa等小组在人体上成功后得到了空前的发展。bold功能磁共振成像现在已经成为认知神经科学中重要的研究手段。由于血氧动力学效应，功能磁共振成像信号的时间分辨率低于脑电图或脑磁图信号，但是功能磁共振成像的空间分辨率要优于脑电图或者脑磁图，可以达到毫米的量级。

采用功能磁共振成像约束脑电图或脑磁图溯源方法确定大脑内电流源位置，不但有效克服了脑电图或者脑磁图数据空间分辨率不足的问题,而且在分析中同时引入了神经活动电磁信息和代谢信息，提高了分析结果的可信度。

1998年，liu等第一次将功能磁共振成像的激活结果用于约束脑电图或脑磁图的源定位，并提出了功能磁共振成像约束最小二次模算法(fmri-weightedminimumnormestimation，fmne)。但是由于脑电图或脑磁图信号反映的是大脑神经活动的电磁特征，功能磁共振成像信号反映的是大脑神经活动的血氧代谢信息，脑电图或脑磁图电流源和功能磁共振成像激活区域间关系并不可知。用功能磁共振成像约束脑电图或者脑磁图溯源的这种方法依赖于两种假设：1、功能磁共振成像激活区位置存在可被脑电或者脑磁记录到的神经活动；2、神经活动可以引起对应的位置产生可探测到的血氧信号。但是并不是所有情况都满足这两种假设，举例说明：如果一个区域内有血氧代谢信号，但是不能记录到对应的电磁信号，这种情况可能是因为神经活动产生的是封闭场，故而违反了第一个假设，会造成功能磁共振成像额外源的存在；如果能在头皮记录到某个神经活动产生的电磁信号但是却测量不到血氧代谢信号，这种情况可能因为这些神经元放电时间太短而造成很少的血氧消耗，这种情况违反了第二个假设，可能会造成功能磁共振成像丢失源的存在。

因此，我们需要一种更好的方法来解决两者位置不匹配所产生的溯源结果不准确的问题。

技术实现要素：

为了解决上述不足，本发明提出一种基于功能磁共振成像的时变约束脑电图或脑磁图溯源方法，将源的协方差矩阵认为是随时间变化，以此为先验条件来约束脑电图或脑磁图溯源，本方法适用于所有的功能磁共振成像约束脑电图或脑磁图溯源分析。

为了解决上述技术问题，本发明采用以下技术方案：

一种基于功能磁共振成像的时变约束脑电图或脑磁图溯源方法，步骤包括：

记录脑电图或脑磁图信号并进行溯源，得到初始的脑电图或脑磁图源强度值，根据所述初始的脑电图或脑磁图源强度值得到脑电图或脑磁图源强度矩阵；

记录功能磁共振成像数据并据以得到功能磁共振成像激活图，根据所述功能磁共振成像激活图得到功能磁共振成像权重矩阵；

根据所述功能磁共振成像权重矩阵和所述脑电图或脑磁图源强度矩阵获得源的协方差矩阵，将该源的协方差矩阵设置为随时间变化；

根据所述初始的脑电图或脑磁图源强度值计算交互通讯效应(crosstalkeffect)和归一化部分曲线下面积(normalizedpartialareaundercurvematrix)，并根据该交互通讯效应和归一化部分曲线下面积修正所述源的协方差矩阵；

根据所述修正的源的协方差矩阵进行溯源，得到最终的脑电图或脑磁图源强度值，完成溯源。

进一步地，根据所述脑电图或脑磁图信号、所述修正的源的协方差矩阵进行溯源所采用的线性方程为：

其中，是脑电图或脑磁图源强度；r是源的协方差矩阵；l是引导场矩阵；c是噪声协方差矩阵，根据刺激前基线数据来确定；x是记录的脑电图或脑磁图信号。其中r和c是未知量，x和l是已知量。

进一步地，根据所述脑电图或脑磁图信号进行溯源时，使用最小二次模算法、深度加权最小二次模算法得到所述初始的脑电图或脑磁图源强度值。

进一步地，根据所述初始的脑电图或脑磁图源强度值得到所述脑电图或脑磁图源强度矩阵所采用的公式为：

其中，是初始的脑电图或脑磁图源强度值；re是脑电图或脑磁图源强度矩阵。

进一步地，根据所述功能磁共振成像权重矩阵和所述脑电图或脑磁图源强度矩阵获得所述源的协方差矩阵所采用的公式为：

r(t)＝rfre(t)(3)

其中，r是源的协方差矩阵；rf是功能磁共振成像权重矩阵；re是脑电图或者脑磁图源强度矩阵。

进一步地，所述交互通讯效应的计算公式为：

其中，ξ是交互通讯效应；g是线性逆算子(inverseoperator)；l是引导场矩阵。

所述归一化部分曲线下面积的计算公式为：

其中，ρ是归一化部分曲线下面积；是第i个位置得到的初始的脑电图或脑磁图源强度值；q是通过所述初始的脑电图或脑磁图源强度值的t检验来确定的向量。

进一步地，所述向量q通过对任务状态和基线状态下的初始的脑电图或脑磁图源强度值做t检验，选用统计检验参数p值和多重比较校正方法。

进一步地，所述统计检验参数p<0.05，所述多重比较校正方法采用邦费罗尼校正(bonferronicorrected)。

进一步地，根据所述计算交互通讯效应和所述归一化部分曲线下面积修正所述源的协方差矩阵的方法是，根据所述计算交互通讯效应和所述归一化部分曲线下面积修正所述功能磁共振成像权重矩阵，进而修正所述源的协方差矩阵，采用的公式为：

r(t)＝rξ,ρre(t)(6)

其中，r是源的协方差矩阵；rξ,ρ是修正的功能磁共振成像权重矩阵，由交互通讯效应ξ和归一化部分曲线下面积ρ两个参数来确定。

进一步地，针对所述修正的功能磁共振成像权重矩阵rξ,ρ，首先根据所述交互通讯效应和所述归一化部分曲线下面积的计算公式计算在功能磁共振成像激活区域内的平均值ρave和ξave，然后在功能磁共振成像激活区域外满足条件：ρi＞ρave并且ξi＜ξave，满足此条件的区域定义为延展的激活区域。

进一步地，根据所述功能磁共振成像激活图得到所述功能磁共振成像权重矩阵的方法是，所述功能磁共振成像权重矩阵的对角元素在实际功能磁共振成像激活区域和延展的激活区域内则置为1，其余位置则置为0.1，非对角元素设置为0。

本方法通过功能磁共振成像约束脑电图或脑磁图信号溯源，由功能磁共振成像权重矩阵和脑电图或脑磁图源强度矩阵两部分共同决定源的协方差矩阵，能够解决功能磁共振成像额外源的影响，并且在一定程度上降低了功能磁共振成像丢失源的影响；通过交互通讯效应和归一化部分曲线下面积修正功能磁共振成像权重矩阵，进而修正源的协方差矩阵，能够更好地解决功能磁共振成像丢失源的影响。在上述的约束条件下，通过消除两方面的影响，达到溯源结果的准确定位。本方法适用于不论神经电活动位置和功能磁共振成像激活区域位置的匹配或者不匹配等所有情况。

附图说明

图1是本发明的溯源方法的流程图。

图2是单侧视觉刺激下得到的事件相关电位曲线图。

图3是单侧视觉刺激下得到的功能磁共振激活图。

图4是mne、fmne和fitc对图2中脑电信号溯源的溯源结果示意图。

具体实施方式

为使本发明的上述特征和优点能更明显易懂，下文特举实施例，并配合所附图作详细说明如下。

应指出的是，本发明公开的仅是具有代表性的一种较佳实施例。显然，本发明并不局限于所描述的具体结构、功能、器件和方法，也可以具有其他实施方式，或者是其他实施方式的组合。本发明中所描述的元素数目也可以设置为多个，除非明确限制为单数。此外，为避免其他例与本发明发生混淆，对于本领域中众所周知的一些技术特征和细节未进行描述。

本实施例提出了一种基于功能磁共振成像的动态约束脑电图或脑磁图溯源方法，亦称之为基于功能磁共振成像的时变约束法(fmriinformedtime-variantconstrains，fitc)。本例以单侧棋盘格视觉刺激下得到的脑电图或脑磁图信号和功能磁共振成像数据为例，详细说明溯源的具体方法，如图1所示。值得强调的是，该方法并不局限于视觉刺激，其他类型刺激任务可以按照相似的步骤处理。

脑电图或脑磁图信号和大脑内电流源的关系可以用线性方程(1)来进行溯源，利用最小二次模算法或其他溯源算法得到初始的脑电图或脑磁图电流源强度值根据功能磁共振成像数据利用现有技术得到功能磁共振成像激活图，再进而根据功能磁共振成像激活图得到功能磁共振成像权重矩阵rf。由于功能磁共振成像激活区域位置与脑电图或脑磁图电流源位置可能会出现不匹配的情况，为了解决这样由于这种不匹配情况产生的溯源问题，我们加入了一些约束条件。

首先解决功能磁共振成像额外源的影响，因为电流源的活动是随着时间发生变化的，它在某一时间的活动位置可能只是功能磁共振成像激活区域的一个子集，因此在这一时间点其他功能磁共振成像激活区域则表现为功能磁共振成像额外的源。为了解决这一问题，本方法的源的协方差矩阵r由功能磁共振成像权重矩阵rf和脑电图或脑磁图源强度矩阵re两部分决定，并且该源的协方差矩阵r是随时间变化的。从公式(3)就可以看出，本方法主要强调了那些在功能磁共振成像激活区内并且神经电流源强度很高的位置，所以可以解决功能磁共振成像额外源的影响，并且在一定程度上降低了功能磁共振成像丢失源的影响。

为了更好地解决功能磁共振成像丢失源的影响，功能磁共振成像权重矩阵rf通过公式(4)和公式(5)中的参数ξ和ρ来决定。首先计算在功能磁共振成像激活区域内的两个参数的平均值ρave和ξave，然后在功能磁共振成像激活区域外满足：ρi大于ρave并且ξi小于ξave，满足这两个条件的位置定义为延展的激活区域。修正的功能磁共振成像权重矩阵rξ,ρ的对角线元素在功能磁共振成像激活区域和延展的激活区域内设置为1，其余设置为0.1，非对角元素设置为0。

在计算归一化部分曲线下面积ρ时，通过对任务状态和基线状态下的初始的脑电图或脑磁图源强度值做t检验，p<0.05，邦费罗尼校正。此处需要说明的是，选用统计检验参数p值和多重比较校正方法不一定非要选择0.05和邦费罗尼校正，只要是计算向量q时采用统计检验的方式，都在本发明的保护范围内。

将经过修正后的源的协方差矩阵rξ,ρ再重新带入公式(1)中进行溯源，得到最终的脑电图或脑磁图源强度值完成溯源。

本发明所述的脑电图或脑磁图信号具体是指脑电图事件相关电位(eventrelatedpotential，erp)或脑磁图事件相关场(eventrelatedfield，erf)，图2是受试者在2赫兹棋盘格的刺激下得到的脑电事件相关电位曲线图，该视觉刺激只出现在视野的右下侧。从图中可以清楚看到与视觉刺激相关的n75、p100和n135等成分，这些成分是脑电波的峰，大概发生在69毫秒、107毫秒和128毫秒。

图3是受试者在相同视觉刺激下得到的功能磁共振激活图，从图中也可以看到包括初级视觉皮层v1以及v2、v3区域的激活。这些区域是人类视觉系统最常见的区域，现有技术中的相关研究指出在观看棋盘格刺激时会导致这三个区域的激活，本图结果与相关研究结论一致，由此可见本方法是可行的。

图4是对图2中脑电信号用三种方法在三个时间点下的溯源得到的源强度值的示意图，第一行是用最小二次模算法(mne)得到的溯源结果，虽然可见视觉区的激活，但是激活区域很宽泛，定位准确性差；第二行是用功能磁共振成像约束最小二次模算法(fmne)溯源得到的溯源结果，相对于最小二次模算法来说，功能磁共振成像约束最小二次模算法得到的激活区更加集中，但是在128ms时出现了错误的激活区。第三行是本发明阐述的方法(fitc)得到的溯源结果，相对于最小二次模算法和功能磁共振成像约束最小二次模算法，利用本发明的方法可以得到更好的定位准确性。

综上所述，利用本发明的方法处理功能磁共振成像约束脑电图溯源结果能够得到更好的溯源结果。此外，本方法适用于不论神经电活动位置和功能磁共振成像激活区域位置的匹配或者不匹配等所有情况，根据同步记录脑电图或者脑磁图和功能磁共振成像数据的广泛应用范围，本发明还可运用到认知神经科学中(比如语言、记忆等)以及精准医疗或者与大脑活动相关的个体差异研究中。