CN102743166A - 一种事件相关电位源定位方法 - Google Patents

Abstract

本发明公开了一种事件相关电位源定位方法，包括以下步骤：事件相关电位提取，按刺激或反应对单试次的脑电进行叠加平均，得到待源定位的事件相关电位数据；假设θ为待求解的源的电位强度，每一项对应皮层上的一个位置，根据每项对应位置的国际自动解剖标记号将θ分为k类，分别对应k个磁共振静息态标准网络，按照静息态标准网络，构建θ的先验分布；计算每个网络的强度参数γ_i，头表记录噪声的协方差强度α；根据以上得到的α，得出皮层的源分布：θ＝α(αL^TL+C^-1)^-1L^TY。本发明利用了功能磁共振的静息态脑网络的空间分布作为先验信息，成像结果有效整合了脑电的高时间分辨率和功能磁共振的高空间分辨率。

Description

一种事件相关电位源定位方法

技术领域

本发明涉及生物信息技术领域，特别涉及脑电信号源电位技术。

背景技术

脑电图因其具有无创、高时间分辨率、能反映大尺度的神经电活动等特点，在脑功能研究与诊断中具有十分重要的地位。脑电是由脑内神经细胞群的电生理活动所产生的电势，经容积导体(由皮层、颅骨、脑膜及头皮等组织构成)传导后，在头皮表面的综合表现。通过对脑电数据的反演研究，可以定量地提供脑内神经活动源的位置、强度及分布情况，这就是脑电逆问题。脑电逆问题在神经科学基础研究和临床应用中具有重要意义。

脑电逆问题是一个不适定问题，其本质上是一个非线性优化问题。但由于非线性问题的复杂性，人们为了简化计算的复杂性，在脑电源的反演定位中，常用一线性方法去逼近非线性问题。因此从采用的方法上，可以把脑电源定位方法分为线性反演和非线性反演方法两类。

非线性反演方法不包括模型中源位置的先验假设，把大脑的一个相对集中小区域的电活动模拟为一个等效偶极子源。这个偶极子代表了该区域神经电活动的位置和方向信息。假设1到5个偶极子，计算产生的头表电活动与记录的电活动差异，可以通过非线性优化迭代得到偶极子的位置的方向参数。

线性反演又叫源成像方法，目前有最小模解(Minimum NormModel，MNM)，加权最小模解(weighted Minimum Norm，WMN)，低分辨层析成像(Low-resolution electromagnetic tomography，LORETA)，多重稀疏先验(Multiple Sparse Prior model，MSP)以及动态统计参数成像(dynamic Statistical Parametric Mapping，dSPM)等。LORETA和MSP采用了解剖上的空间邻接关系作为先验，这两种方法中临近的区域被认为有相似的神经活动。dSPM则采用功能磁共振的激活信息作为先验。专利申请人于2011年前后国际上首次提出采用功能磁共振的连通信息作为先验的做法：网络源成像(Network-based source imaging，NESOI)[Lei X，Xu P，Luo C，Zhao J，Zhou D，Yao D(2011)：fMRIFunctional Networks for EEG Source Imaging.Human Brain Mapping32：1141-1160.]。该方法根据脑区的功能磁共振是否有时间相关性，将脑区划分为大的模块，作为脑电源定位的先验。这些功能网络与空间邻接信息(LORETA采用)或功能激活信息(dSPM采用)不同，它们往往覆盖了空间上相隔很远的多个脑区。本方法是在网络源成像的基础上，将功能连接信息具体化和规范化，使用数个静息态脑网络作为先验对事件相关电位进行源成像。

发明内容

本发明的目的是提供一种引入静息态脑网络信息，更准确的定位事件相关电位的产生源的事件相关电位源定位方法。

本发明解决所述技术问题采用的技术方案是，通过贝叶斯理论，使用数个已知的静息态脑网络的空间分布信息和头表电位分布信息，得到每个网络对头表电位的贡献强度，求出事件相关电位的源分布。包括以下步骤：

(a)、事件相关电位提取，按刺激或反应对单试次的脑电进行叠加平均，得到待源定位的事件相关电位数据，用Y表示，Y为m行1列的向量，m表示电极个数；

(b)、假设θ为待求解的源的电位强度，为n行1列的向量，n为解空间的大小，即大脑皮层分布的偶极子个数。对θ的每一项θ_i，都有一个对应的范围在1到116之间的国际自动解剖标记号(Automatedanatomical labeling，AAL[Tzourio-Mazoyer N，Landeau B，Papathanassiou D，Crivello F，Etard O，Delcroix N，Mazoyer B，Joliot M(2002)：Automated anatomical labeling of activations in SPM using amacroscopic anatomical parcellation of the MNI MRI single-subject brain.Neuroimage 15：273-289.])，解剖标记号的划定是Tzourio-Mazoyer等人通过对一个具有解剖代表性的健康受试者进行多次MRI扫描T1加权像进行空间平均化之后做出的。他们将脑沟描绘出，据此做为分界标志在脑的两个半球各划分45个区域，由此给90个区域分别人为指定了一个标记号。标记号相同的脑区构成了一个功能模块，在认知活动中发挥相似的功能。

通过对39例受试者的静息态功能磁共振数据的独立成分分解，可得到相似时间序列波形的K个皮层网络。K个网络实际覆盖了源空间的绝大部分区域，根据AAL和K个皮层网络重合的大小，AAL相同的脑区被唯一的划分到K个网络之一。建议的K取值为9，该值是根据smith等(2009)研究发现的10个静息态脑网络除去1个小脑网络(不在源定位的解空间上)后的结果，即：视觉网络中部、视觉网络侧部、视觉网络腹侧、默认模式网络、运动感觉皮层网络、听觉网络、决策网络、偏左侧的额顶叶网络、偏右侧的额顶叶网络。我们根据39例静息态功能磁共振数据构造了9个网络对应的AAL标记表，即所谓的静息态标准网络，AAL和静息态标准网络对应关系如下：

表1AAL和静息态标准网络对应关系

不在以上9个网络中的其余的偶极子被对应的划分在第十个网络中。

标准化后的静息网络空间分布见图2。需要强调的是，在实际源定位中由于大脑皮层分布的每个偶极子对应的AAL标记号已知，可以直接通过上表将其划分到不同的静息态标准网络，而不再需要静息态功能磁共振扫描。

(c)、按照静息态标准网络，构建θ的先验分布，即θ满足均值为0(n行1列的全零向量)，协方差为C(n行n列)的正态分布。C是协方差矩阵V_i的加权和：

其中γ_i为协方差矩阵的强度参数，V_i为n行n列的协方差矩阵，主对角线上的每一项对应一个激活位置。

V_i(1≤i≤K)根据第i个静息态标准网络确定，即V_i主对角线上的第j项对应的位置如果在第i个静息态标准网络中，则为1，否则为0。大脑皮层中不被K个网络覆盖的区域构成第K+1个协方差矩阵，即V_K+1其主对角线上的第j项对应的位置如果不在任何一个静息态标准网络中，则为1，否则为0。可以看出，V_i只提供空间位置信息，而每个网络的激活强度用参数γ_i控制。

(d)、对事件相关电位Y源定位的过程就是解线性方程组

Y＝Lθ+ε

其中ε满足正态分布的向量，该分布的均值为0(m行1列的全零向量)，协方差为1/α×I_m，I_m为项数为m行m列的单位矩阵，α为权重系数；θ为待求解的皮层源的电位强度，L为从电流从皮层传递到头表的传递矩阵，Y为事件相关电位，均为已知。总的来说，未知的项包括：每个网络的强度参数γ_i，头表记录噪声的协方差强度α，皮层源的电位分布：θ。

传递矩阵L的计算是根据电场传播的麦克斯韦方程(Maxwell，1865)得到的，在这里对其计算方法进行简单介绍。L的维数和解空间离散的偶极子个数有关系，L中的一列通常是代表着在某一网格中心处沿垂直网格表面的方向，单位强度的偶极子在头表产生的电位分布。其计算方法是这样的：将源解空间按一定的精度进行三角片网格剖分后，在每一个节点(偶极子)上放置偶极子沿垂直网格表面方向的单位偶极子，然后采用一定的正向计算方法分别计算头表电位分布，这个分布构成L中的1列，对每个节点都重复类似的计算，就可以获得传递矩阵L。可见在采用偶极子模型时，一般情况下L的列数是剖分偶极子的个数，而行数为电极个数。

(e)、利用约束最大化似然(restricted maximum likelihood，ReML)算法求解参数γ_i和α。

(f)、将γ_i代入

得到源空间偶极子激活强度的协方差C，根据以上得到的α，得出皮层的源分布：

θ＝α(αL^TL+C^-1)^-1L^TY

从线性方程到解出θ的推导过程为：根据解线性方程组Y＝Lθ+ε可知θ的似然概率密度函数为：P(Y |θ，α)＝N(Lθ，α^-1I_m)，I_m为m行m列的单位矩阵，先验概率密度函数为：P(θ|γ)＝N(0，C)。根据贝叶斯公式，得到θ的后验概率密度函数为：P(θ|Y，γ，α)＝P(Y|θ，α)P(θ|γ)/∫P(Y|θ，α)P(θ|γ)。另一方面，由于似然概率密度和先验概率密度均为高斯分布，θ的后验概率密度函数也应为高斯分布：P(θ|Y，γ，α)＝N(μ，∑)，∑为θ后验概率协方差；对比两个后验概率密度函数得：μ＝α(αL^TL+C^-1)^-1L^TY，此即为θ的最大后验估计，μ为θ的均值。相应的，可以得出

$\begin{array}{c}\mathrm{\&gamma;}=\underset{\mathrm{\&gamma;}}{\max }P\left(Y\right|\mathrm{\&gamma;})=\underset{\mathrm{\&gamma;}}{\max }<\ln P\left(Y\right|\mathrm{\&theta;},\mathrm{\&alpha;})+P\left(\mathrm{\&theta;}\right|\mathrm{\&gamma;})>\\ \mathrm{\&alpha;}=\underset{a}{\max }<\ln P\left(Y\right|\mathrm{\&theta;},\mathrm{\&alpha;})>\end{array}$

将γ_i代入得到源空间偶极子激活强度的协方差C，根据以上得到的α，得出皮层的源分布：

θ＝α(αL^TL+C^-1)^-1L^TY

L为从电流从皮层传递到头表的传递矩阵。

本发明中我们将常用的电极(共343个)的传递矩阵模板都在软件Fieldtrip中计算得到(http://fieldtrip.fcdonders.nl/download.php)，并存储下来。当需要任何电极时，只需查找该次实验的电极名，就可从传递矩阵模板中快速提取需要的传递矩阵。

传递矩阵其计算方法在电磁场研究中已成定论，参见(尧德中：脑功能探测的电学理论与方法2003年科学出版社)。这里简要描述如下：

对于我们构造的皮层三角片网格采用三层同心球模型进行包裹。此模型将人头部考虑为头皮层、头盖骨、脑脊液的3个同心球。在这个模型中，头盖骨的内外半径为8cm和8.5cm，而头表的半径为9.2cm。脑脊液和头皮层的电阻率设为2.22Ωm，而头盖骨的电阻率为80×2.22Ωm＝177Ωm，这些都是用来代表典型的生物生理数值。在确定好这些参数后，对于脑皮层任意位置有一个强度为v的偶极子，设其球坐标为

根据电场传播的麦克斯韦方程，该位置传导到头表电极的电位为：

其中δ为大脑的电导率，R为头皮层半径，

为Legendre函数，n为自然数1，2，3...∞。上式为一个位置的结果，矩阵L的一列代表着在某一网格中心处沿垂直网格表面的方向，单位强度的偶极子在头表产生的电位分布。重复上面的步骤，在每一个网格节点(偶极子)上放置偶极矩沿垂直网格表面方向的单位偶极子，就可以获得传递矩阵L。可见L的列数是剖分偶极子的个数。

本发明利用了功能磁共振的静息态脑网络的空间分布作为先验信息，成像结果有效整合了脑电的高时间分辨率和功能磁共振的高空间分辨率。采用了标准化的功能磁共振的静息态脑网络，不需进行功能磁共振扫描，就可达到高时空分辨率。与MNM，LORETA，MSP等其他线性源成像方法对比，定位误差更小，时间相似度更强，数据拟合度更高。

附图说明

以下结合附图和具体实施方式对本发明做进一步的说明。

图1为事件相关电位源成像的算法流程图；

图2为通过静息态标准网络构造协方差矩阵的示意图；

图3为人脑部偶极子的示意图；

图4为人脸和混拼人脸刺激产生的事件相关电位差异，上半部分为差异波的波形，下半部分为差异波的头表分布图，分别对应刺激呈现后100ms、170ms、240ms。

图5为本发明在人脸识别的源强度分布在皮层的投影(视角为从枕部往前看)：

图6为本方法采用4个网络的源定位结果；

图7为本方法采用2个网络的源定位结果。

具体实施方式

如图1所示，图1对本发明做的大体流程进行了描述。

一种事件相关电位源定位方法，包括以下步骤：

(a)、事件相关电位提取，按刺激或反应对单试次的脑电进行叠加平均，得到待源定位的事件相关电位数据，用Y表示，含m行1列，m为电极个数；

(b)、假设θ为待求解的源的电位强度，每一项对应皮层上的一个位置，有n行1列，根据每项对应位置的国际自动解剖标记号(Automated anatomical labeling)将θ分为k类，分别对应k个磁共振静息态标准网络，同一类的源的激活强度相近，不同类源的激活强度不同，其中k为自然数；

(c)、按照静息态标准网络，构建θ的先验分布，即θ满足均值为0(n行1列的全零向量)，协方差为C(n行n列)的正态分布；C是协方差矩阵V_i的加权和：

其中γ_i为协方差矩阵的强度参数，V_i(1≤i≤k)根据第i个静息态标准网络确定，大脑皮层中不被前面k个静息态标准网络覆盖的区域构成第k+1个协方差矩阵；

(d)、利用约束最大化似然算法(restricted maximum likelihood，ReML)计算每个网络的强度参数γ_i，头表记录噪声的协方差强度α；从线性方程到解出θ的推导过程为：根据解线性方程组Y＝Lθ+ε可知θ的似然概率密度函数为：P(Y|θ，α)＝N(Lθ，α^-1I_m)，I_m为m行m列的单位矩阵；先验概率密度函数为：P(θ|γ)＝N(0，C)。根据贝叶斯公式，得到θ的后验概率密度函数为：P(θ|Y，γ，α)＝P(Y|θ，α)P(θ|γ)/∫P(Y|θ，α)P(θ|γ)；另一方面，由于似然概率密度和先验概率密度均为高斯分布，θ的后验概率密度函数也应为高斯分布：P(θ|Y，γ，α)＝N(μ，∑)，∑为θ后验概率协方差；对比两个后验概率密度函数得：μ＝α(αL^TL+C^-1)^-1L^TY，此即为θ的最大后验估计，μ为θ的均值，相应的，可以得出

$\begin{array}{c}\mathrm{\&gamma;}=\underset{\mathrm{\&gamma;}}{\max }P\left(Y\right|\mathrm{\&gamma;})=\underset{\mathrm{\&gamma;}}{\max }<\ln P\left(Y\right|\mathrm{\&theta;},\mathrm{\&alpha;})+P\left(\mathrm{\&theta;}\right|\mathrm{\&gamma;})>\\ \mathrm{\&alpha;}=\underset{a}{\max }<\ln P\left(Y\right|\mathrm{\&theta;},\mathrm{\&alpha;})>\end{array}$

(e)、将γ_i代入

得到源空间偶极子激活强度的协方差C，根据以上得到的α，得出皮层的源分布：

θ＝α(αL^TL+C^-1)^-1L^TY

L为从电流从皮层传递到头表的传递矩阵。

其中k的值为9。

为了说明本发明所提及的一种脑电源定位方法的过程及效果，将采用一段真实事件相关电位数据Y，通过128道Active Two系统采集，另有四道位于双耳位置和眼部周围(眼电电极：VEOG，HEOG)，采样率1024Hz。

EEG的源空间由分布在大脑皮层的8196个偶极子构成，通过3个邻近偶极子的连接，构成如图3所示的三角片网格，三角片网络是在一个解剖学典型的人脑上得到的。采用文献[Henson RN，Mattout J，Phillips C，Friston KJ(2009)：Selecting forrd models for MEGsource-reconstruction using model-evidence.NeuroImage 46：168-176.]中的方法，在兼顾均匀性和局部解剖特异性的基础上全自动进行的。该网格从一个具有神经解剖典型性的磁共振结构图中提取。根据各个偶极子的坐标，可以找出偶极子所在的AAL标记号，再由表1得到对应的静息态标准网络号。将得到10个8196行8196列的协方差矩阵V_i，这些矩阵只有主对角线上有值，其他位置均为0。V_i根据第i个静息态标准网络确定，如图3中，通过第1号偶极子的坐标，我们知道该位置处于AAL的44号区，查询表1得知，它位于视觉网络中部，则V₁的1号位填1，其他V_i的该位置都填0；再如，第5100号偶极子的坐标处于AAL的17号区，查询表1得知，它位于听觉网络，则V₆的5100号位填1，其他V_i的该位置都填0。由于8196个偶极子与静息态标准网络的对应关系是固定的，在实际源定位中V_i被固定保存下来，而不再需要每次进行定位时都构造一遍V_i。

事件相关电位的预处理包括平均参考和伪迹去除。绝对值大于120μV的坏段将被剔除，344次试验中共去除80次实验。人脸和混拼人脸图片相关的事件相关电位分别进行叠加，并按刺激的-200ms到0ms进行基线校正。两种任务的差异ERP及相应的头表分布如图4所示。可以看到，在刺激呈现后的170ms附近有一个峰值，对应人脸识别特异的N170成分。以上都是事件相关电位提取的标准流程，图5中175ms的头表电位就是我们获得的需要定位的脑电数据：Y＝[7.66 7.35 6.61 5.82 4.38 4.09 2.94 -0.800 -2.67 -4.87 -4.55 -3.86 -4.05-3.48 -2.96 -1.81 1.21 2.58 4.38 2.91 1.09 -1.49 -2.68 -2.20 -2.29 -3.26-4.26 -4.72 -2.83 -0.440 2.37 3.28 7.24 7.13 4.78 2.79 0.210 -3.19 -5.15-5.01 -4.59 -3.47 -1.87 0.280 1.59 5.07 1.66 3.20 4.24 5.71 5.07 6.01 6.195.49 5.21 3.39 3.42 1.89 2.71 4.29 4.38 5.07 5.17 5.97 7.31 6.92 5.624.68 4.57 3.59 3.21 0.0600 4.39 5.09 6.15 6.14 5.59 6.42 2.51 4.39 5.605.06 5.60 5.66 6.35 6.85 7.01 6.20 5.81 5.34 4.81 4.81 5.45 4.60 3.597.69 7.58 6.89 5.57 4.59 3.28 3.51 2.97 3.20 3.37 3.97 4.63 4.81 5.746.62 7.18 5.82 5.68 6.09 4.91 4.06 3.67 2.77 2.71 0.390 1.52 2.70 2.864.57 0.770 -0.900 -2.20 -3.77]^T，即128行1列的电位分布数据，该数据是源定位的输入。

由偶极子的位置和电极位置，可以根据麦克斯韦方程计算得到128行8196列的传递矩阵L。将Y作为输入，借助约束最大化似然算法得到γ₁到γ₁₀的取值为[0.01 0 3.12 0 0 0 0 0 0 0]，即第3号网络视觉网络腹侧对该事件相关电位的产生具有主要贡献，另外得到α的值为0.1。将γ_i代入

得到源空间偶极子激活强度的协方差C，代入α，根据θ＝α(αL^TL+C^-1)^-1L^TY得出皮层的源分布：θ＝[-0.20-0.49 -0.63 -0.41 -0.54 -0.26 -0.16 -0.040 -0.18 0.010 0.12 0.060 -0.14-0.25 -0.30 -0.17 -0.59 -0.41 -0.37 -0.43 -0.44 -0.31 -0.060 0.030 -0.51-0.64 -0.46 -0.24 -0.83 -1.04 -0.87 -0.79 -0.86 -0.82 -0.53 0.020 0.46-0.71 -0.55 -0.85 -0.81 -0.78 -0.87 -0.79 -0.35 0.33 0.85 -0.43 -0.59 -0.500.46 -0.94 -0.14 -0.040 -0.050 -0.12 3.66 4.66 5.46 0.24 0.28 -0.060 0.110.13 0.29 0.72 -0.73 -0.75 -0.54 -0.27 -0.16 -0.14 -0.070 0.11 0.38 0.940.17 0.30 0.45 0.27 0.040 -0.070 -0.57 -0.48 -0.34 -0.26 -0.29 -0.34-0.24...]^T，是8196行1列的皮层电位分布数据。该数据是本方法的输出，将θ的数值取绝对值后在皮层上按灰度投影，会得到如图5所示的结果(从0开始随着数据的增加，颜色变化为白-黄-红-黑；既值越小越接近白色，值越大越接近黑色)。

为了与此进行对比，对同样的事件相关电位数据Y另外采用了将全脑划分为4个网络和划分成2个网络进行源定位的结果，其结果为图6和图7所示。

通过采用不同网络数量进行源定位以及与其它源定位方法的对比可以看出，本发明的方法成功的将双侧的基底节定位出来了，且划分的网络数量增多时，结果趋于更精确。

Claims (2)

1.一种事件相关电位源定位方法，其特征在于：所述方法包括以下步骤：

(a)、事件相关电位提取，按刺激或反应对单试次的脑电进行叠加平均，得到待源定位的事件相关电位数据，用Y表示，含m行1列，m为电极个数；

(b)、假设θ为待求解的源的电位强度，每一项对应皮层上的一个位置，有n行1列，根据每项对应位置的国际自动解剖标记号(Automated anatomical labeling)将θ分为k类，分别对应k个磁共振静息态标准网络，同一类的源的激活强度相近，不同类源的激活强度不同，其中k为自然数；

(c)、按照静息态标准网络，构建θ的先验分布，即θ满足均值为0(n行1列的全零向量)，协方差为C(n行n列)的正态分布；C是协方差矩阵V_i的加权和：

其中γ_i为协方差矩阵的强度参数，V_i(1≤i≤k)根据第i个静息态标准网络确定，大脑皮层中不被前面k个静息态标准网络覆盖的区域构成第k+1个协方差矩阵；

(d)、利用约束最大化似然算法(restricted maximum likelihood，ReML)计算每个网络的强度参数γ_i，头表记录噪声的协方差强度α；

$\begin{array}{c}\mathrm{\&gamma;}=\underset{\mathrm{\&gamma;}}{\max }P\left(Y\right|\mathrm{\&gamma;})=\underset{\mathrm{\&gamma;}}{\max }<\ln P\left(Y\right|\mathrm{\&theta;},\mathrm{\&alpha;})+P\left(\mathrm{\&theta;}\right|\mathrm{\&gamma;})>\\ \mathrm{\&alpha;}=\underset{a}{\max }<\ln P\left(Y\right|\mathrm{\&theta;},\mathrm{\&alpha;})>\end{array}$

(e)、将γ_i代入

得到源空间偶极子激活强度的协方差C，根据以上得到的α，得出皮层的源分布：

θ＝α(αL^TL+C^-1)^-1L^TY

L为从电流从皮层传递到头表的传递矩阵。

2.如权利要求1所述的事件相关电位源定位方法，其特征在于：所述k的值为9。

Images