CN1903119A - 基于稀疏性约束的脑电源定位方法 - Google Patents

Abstract

基于稀疏性约束的脑电源定位方法，属于生物信息技术领域，本发明包括以下步骤：1)确定传递矩阵A；2)通过多道脑电信号记录系统获取实际记录的脑电信号Y，进行预处理，确定要进行源分析的时刻；3)初始化源向量x_k－1，k＝1，迭代终止误差ε和一个最大的迭代次数；4)更新权重对角矩阵：W_k＝(diag(x_k－1))；5)加入稀疏性约束，利用优化方法直接估计q_k：argmin‖b－AW_kq_k‖₂+λ‖q_k‖_p；6)更新源信息向量：x_k＝W_kq_k；7)迭代终止条件判断：比较更新前后的源分布变化，当‖x_k－x_k－1‖≤ε或迭代次数超出设定的最大迭代次数限制时，迭代终止，x_k即为源的最终定位估计结果；否则，k＝k+1，转步骤2)，继续迭代；其中，λ≠0，W_k≠I。本发明将l_p (p≤1)模约束融合到FOCUSS的迭代过程中，获得更加稳定的稀疏源定位结果。

Description

基于稀疏性约束的脑电源定位方法

技术领域：

本发明属于生物信息技术领域，涉及一种对脑电源进行定位的方法，主要应用于人脑功能及与人脑相关疾病的研究与诊断。

背景技术：

根据头表观测电位反演定位脑电源的空间活动位置是脑电研究中的一个重要问题，其本质上是一个非线性优化逆问题。为了简化计算的复杂性，在脑电源的反演定位中，常用一线性方法去逼近非线性问题，从而可把脑电逆问题归结为如下的一个线性系统(C.M.Michel，M.Murray，EEG source imaging，Clinical neurophysiology，115，2195-2222，1997)，

        b＝Ax                                (1)

其中b为头表电极记录到的电位，是一M×1的向量，M为头表记录电极数目；A为传递矩阵(Lead Field)，是一M×N的矩阵，N是源活动解空间的维数；x为一N×1的向量，是待进行空间定位的源信息向量。在脑电逆问题中，通常M＜＜N，因此，式(1)所示的线性系统是一个严重欠定的系统。对上面的方程，有无穷多组解x满足观测到的电压分布b，所以为了获得符合实际生理特性的解，人们往往根据不同的需要和研究目的对上面的逆问题施加合适的约束(尧德中.脑功能探测的电学理论与方法.北京：科学出版社，2003，195-243)。在早期的研究中，主要采用了最小模解(minimum norm solution，MNS)，在此基础上发展到了当前常用的加权模解(weightedminimum-norm solution，WMN)，其中最具代表性的是LORETA(Lowresolution electromanetic tomography)方法，它是一种采用Laplacian微分算子对解空间进行加权的方法。采用WMN方法，通常只能得到源活动的块状模糊区域，达不到神经科学研究需要的比较精细的源定位(R.D.Pascual-Marqui，Review of methods for solving the EEG inverseproblem，Int.J Bioelectromagnetism，1，75-86，1999；C.M.Michel，M.Murray，EEG source imaging，Clinical neurophysiology，115，2195-2222，1997)

由于人脑在执行某种任务或存在某种疾病时，其对应的主要神经活动区域是局部性和稀疏性的，因此，如果把这两个特性作为欠定系统求解时的约束，理论上可以获得更加符合脑电源活动生理特性的解(C.M.Michel，M.Murray，EEG source imaging，Clinicalneurophysiology，115，2195-2222，1997；尧德中.脑功能探测的电学理论与方法.北京：科学出版社，2003，195-243)。在当前的脑电逆问题技术中，已建立了两种方式来获得源的稀疏解：一是被称为FOCUSS(Focal underdetermined system solver)算法的一种迭代WMN方法；二是直接求解l_p(p≤1)模约束的解。FOCUSS通过逐步迭代，使解空间的大多数元素趋向零，从而使解的能量局部化(I.F.Gorodnitsky，B.D.Rao，Sparse signal reconstruction from limited datausing FOCUSS：A re-weighted minimum norm algorithm，IEEE Trans.S.P.，45，600-616，1997)。它通常是对LORETA等方法获得的模糊解进一步迭代，来获取能量局部化的稀疏解。FOCUSS的基本迭代过程可以叙述为：采用一线形变换x＝Wq，当前迭代的加权矩阵W_k取为由前一次迭代结果元素x_k-1构成的对角阵，记为W_k＝(diag(x_k-1))，q为一中间变量。FOCUSS可以由下面3步迭代来完成：

            1⁰ W_k＝(diag(x_k-1))

            2⁰ q_k＝(AW_k)⁺b

            3⁰ x_k＝W_kq_k

其中(AW_k)⁺表示对矩阵AW_k进行广义逆的求解。

从这里的迭代过程可以看出，在FOCUSS的迭代过程中有一个求解矩阵逆的过程，对矩阵逆的估计好坏，在很大程度上决定了FOCUSS算法的稳定性和对噪声的免疫能力。当前对FOCUSS的改进，主要是用各种正则化技术，包括奇异值截断和各种形式的Tikhonov正则化技术改善矩阵逆的特性。第2种采用l_p(p≤1)模约束的方法，就是把式(1)表示的逆问题转化为如下的问题(C.Silva，J.C.Maltez，E.Trindade，Evaluation of L1 and L2 minimum norm performances on EEGlocalizations，Clinical neurophysiology，115，1657-1668，2004)，

         arg min‖b-Ax‖₂+λ‖x‖_p              (2)

然后直接对(2)式采用优化方法进行求解，以获得一个稀疏的源分布x。在脑电逆问题中，通常p＝1。在这里，稀疏性约束是通过一个类似于Tikhonov正则化的项引入的，它既有稀疏性约束的作用，也有正则化的作用，从而使得整个计算过程是稳健的。但是，该技术没有迭代过程，使它得到的解的稀疏性不如FOCUSS方法。

上述的FOCUSS和l_p(p≤1)模算法都是求取源的稀疏解，但两种技术各有特点。本技术立足脑电源活动的稀疏性生理特性，将l_p(p≤1)模约束融合到FOCUSS的迭代过程之中，获得了更加稳健的稀疏源定位结果。

发明内容

本发明所要解决的技术问题是，提供一种基于l_p模约束的迭代WMN脑电源定位方法(简称LIPSS)。该方法将l_p模稀疏性约束融合到FOCUSS的迭代过程之中，以获得更加稳定而稀疏的源定位结果。

本发明解决所述技术问题采用的技术方案是，基于稀疏性约束的脑电源定位方法，包括以下步骤：

1)确定传递矩阵A；

2)通过多道脑电信号记录系统获取实际记录的脑电信号Y，进行预处理，确定要进行源分析的时刻；

3)初始化源向量x_k-1，k＝1，迭代终止误差ε和一个最大的迭代次数；

4)更新权重对角矩阵：W_k＝(diag(x_k-1))；

5)加入稀疏性约束，利用优化方法直接估计q_k：arg min‖b-AW_kq_k‖₂+λ‖q_k‖_p；

6)更新源信息向量：x_k＝W_kq_k；

7)迭代终止条件判断：比较更新前后的源分布变化，当‖x_k-x_k-1‖≤ε或迭代次数超出设定的最大迭代次数限制时，迭代终止，x_k即为源的最终定位估计结果；否则，k＝k+1，转步骤2)，继续迭代；

其中，λ≠0，W_k≠I(I为单位矩阵)。

所述步骤1中，包括以下分步骤：

a1)、对待测对象的头部进行MRI或CT扫描，获取头部解剖结构的影像信息；

a2)、提取步骤a1)所得的影像信息中的大脑部分，然后对大脑分割，再提取大脑部分的源功能区(主要包括灰质、海马、小脑等部位)；

a3)、以一定精度的网格将步骤a2)所得的脑源功能区进行剖分，确定解空间网格(包括解空间的维数和各个网格的空间位置序号)；

a4)、确定多道脑电信号记录系统各个电极的空间位置信息；

a5)、确定脑电源的模型；

a6)、利用步骤a3)到步骤a5)中确定的解空间网格、电极位置信息和脑电源模型，利用正演方法计算传递矩阵A，具体方法如下：在每个解空间位置上放置单位的源，利用数值计算方法计算该单位源在电极位置处产生的电位分布，该电位分布构成传递矩阵中的一列，以此类推，当把所有解空间遍历放置单位源后，就可以获得传递矩阵A。

特别的，p取值为1。

所述步骤2)中的预处理包括：滤波、基线校正、眼电去除。

本发明的有益效果是，相比较以前的方法而言，本发明采用LPISS技术，从脑电源活动的稀疏性生理事实出发，将l_p(p≤1)模约束融合到FOCUSS的迭代过程中，能获得更加稳定的稀疏源定位结果。

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

附图说明

图1：在NSR＝15％下对两个源的LORETA，FOCUSS和LPISS定位结果示意图。彩色的方形区域代表着估计的源位置，不同的颜色代表着不同的能量；标志有蓝色十字叉的彩色方形区域表示着估计源位置和模拟源位置间的重叠部分；在绿色圆圈内的蓝色十字叉表示模拟源和估计源间的不重叠部分。

其中a为LORETA定位结果，左(-23.00，-51.00，51.67)(mm)，右(-43.00，-41.00，61.67)(mm)。

b为FOCUSS定位结果，左(-23.00，-51.00，51.67)(mm)，右(-43.00，-41.00，61.67)(mm)。

c为LPISS定位结果，左(-23.00，-51.00，51.67)(mm)，右(-43.00，-41.00，61.67)(mm)。

图2：不同噪声条件下对两个源的定位性能指标图。其中a为源1，b为源2。

图3：深浅混合源的定位结果示意图。彩色的方形区域代表着估计的源位置，不同的颜色代表着不同的能量；标志有蓝色十字叉的彩色方形区域表示着估计源位置和模拟源位置间的重叠部分；在绿色圆圈内的蓝色十字叉表示着模拟源和估计源间的不重叠部分。其中a为LORETA定位结果，a1：(-83.00，19.00，31.67)(mm)；a2：(-23.00，29.00，-8.33)(mm)；a3：(-63.00，39.00，-48.33)(mm)。

b为FOCUSS定位结果，b1：(-83.00，19.00，31.67)(mm)；b2：(-23.00，29.00，-8.33)(mm)；b3：(-63.00，39.00，-48.33)(mm)。

c为LPISS定位结果。c1：(-83.00，19.00，31.67)(mm)；c2：(-23.00，29.00，-8.33)(mm)；c3：(-63.00，39.00，-48.33)(mm)。

图4：LPISS对IOR数据的源定位结果示意图。其中a为右额叶眼区，b为丘脑，c为右小脑，d为左枕叶。

具体实施方式

本发明采用LPISS技术，具体的说，包括以下步骤：

步骤1.确定传递矩阵A，包括以下分步骤：

a1)对待测对象的头部进行MRI或CT扫描，获取头部解剖结构的影像信息；

a2)提取步骤a1)所得的影像信息中的大脑部分，然后对大脑分割，再提取大脑部分的源功能区(主要包括灰质、海马、小脑等部位)；

a3)以一定精度的网格将步骤a2)所得的脑源功能区进行剖分，确定解空间网格(包括解空间的维数和各个网格的空间位置序号)；

a4)确定多道脑电信号记录系统各个电极的空间位置信息；

a5)确定脑电源的模型；

a6)利用步骤a3)到步骤a5)中确定的解空间网格、电极位置信息和脑电源模型，利用正演方法计算传递矩阵A，具体方法如下：在每个解空间位置上放置单位的源，利用数值计算方法计算该单位源在电极位置处产生的电位分布，该电位分布构成传递矩阵中的一列，以此类推，当把所有解空间遍历放置单位源后，就可以获得传递矩阵A；

步骤2.通过多道脑电信号记录系统获取实际记录的脑电信号Y，并进行必要的如：滤波、基线校正、眼电去除等基本的预处理，并确定要进行源分析的时刻；

步骤3.初始化源向量x_k-1，k＝1，迭代终止误差ε和一个最大的迭代次数；

步骤4.更新权重对角矩阵：W_k＝(diag(x_k-1))；

步骤5.加入稀疏性约束，利用优化方法直接估计q_k：

arg min‖b-AW_kq_k‖₂+λ‖q_k‖_p；

步骤6.更新源信息向量：x_k＝W_kq_k；

步骤7.迭代终止条件判断：比较更新前后的源分布变化，当‖x_k-x_k-1‖≤ε或迭代次数超出设定的最大迭代次数限制时，迭代终止，x_k即为源的最终定位估计结果；否则，k＝k+1，转步骤2，继续迭代。

优选的取值是，p＝1；在第3步中采用基于BFGS变尺度的优化方法对l_p(p≤1)模问题进行求解。

上述方案中，步骤1.的步骤a3)中所述的一定精度的网格，综合考虑计算精度和效率，一般取10mm/格；步骤1.的步骤a4)中所述的多道脑电信号记录系统可以是标准的32道、64道、128道及256道电极的脑电信号记录系统；步骤1.的步骤a5)中所述的脑电源模型通常为点电荷模型或偶极子模型；步骤1.的步骤a6)中所述的数值计算方法可以是边界元算法或有限元算法；步骤3中的源的初始解设为LORETA的解；步骤5中的p＝1，采用基于BFGS变尺度的优化方法对l_p(p≤1)模问题进行求解。步骤1实际上能根据现有技术实现。

作为应用效果，基于一真实的标准头模型，利用模拟源对算法进行了定位性能分析比较，并应用到对返回抑制(Inhibition of return，IOR)脑电实验数据的源定位分析。传递矩阵A是按这样的方式进行具体计算：通过扫描获得的MRI头模型，把偶极子源活动位置限定在大脑的灰质、海马及其他可能的源活动区域，通过10mm网格剖分离散成910个位置，在每一个网格位置上分别放置沿X，Y，Z方向的偶极子，采用标准128道电极系统，利用边界元(Boundary Elements Method，BEM)方法计算获得传递矩阵A。在实施时，LPISS的最大迭代次数为100，算法迭代终止误差ε为1E-5。

在评估算法性能时，主要采用了如下的三个性能指标：1)空间定位误差，E_localization，定义为模拟源和在包围模拟源的一个球状区域内的，具有最大能量的估计源的位置间的距离；2)源能量误差，E_energy，可以用来评估算法对源能量的恢复能力，其计算定义为： $E_{\mathrm{energy}}=|1-\frac{\left|\right|J_{\max }\left|\right|}{\left|\right|J_{\mathrm{simu}}\left|\right|}|,$ 其中J_simu是模拟源的偶极矩，J_max是在相应模拟源周围的一定半径的球形邻域内的具有最大能量的偶极子的偶极矩(C.M.Michel，M.Murray，EEG source imaging，Clinical neurophysiology，115，2195-2222，1997)；3)在感兴趣区域(region of interest，ROI)内的归一化模糊指数(normalized burring index，NBI)，NBI可以用来度量算法的空间分辨能力，NBI的定义为(D Yao，B.He，A Self-CoherenceEnhancement Algorithm and its Application to Enhancing 3D SourceEstimation from EEGs，Annals of Biomedical Engineering，29(11)，1019-1027，2001)，

${\mathrm{NBI}}_i=\sqrt{\underset{k}{\mathrm{\Σ}}{\left|\right|r_k-r_i\left|\right|}^2{J}^2\left(k\right)}/\sqrt{\underset{k}{\mathrm{\Σ}}{\left|\right|r_k-r_i\left|\right|}^2/\underset{k}{\mathrm{\Σ}}1}$

在本工作中，计算NBI时，采用的都是球形ROI。在NBI定义中的各符号意义如下：下标i代表着待计算NBI的球形ROI的中心点的网格序号，对模拟假设的源分布的NBI计算时，i就是模拟源的真实空间位置所对应的网格序号；对估计源分布的NBI进行计算时，i选为在模拟源的一定球形ROI内具有最大能量的估计偶极子位置所对应的网格序号。下标k是指在ROI内的所有有效网格节点序号。r_k，r_i分别是序号为k，i的网格点所对应的空间位置坐标。NBI可以表征在ROI区域内的源的分布情况。通常，源在ROI内分布越平缓，NBI越接近于1；如果在ROI内源呈现尖锐分布，则NBI接近于0。显然，离散分布的源，具有较小的NBI。一个优秀的源定位算法应该具有较小的E_localization和E_energy，同时估计的源分布的NBI应该和真实源分布NBI较相似。

1.定位性能的模拟实验分析比较：

方法：首先测试了算法在不同噪声条件下的源定位能力。将两个偶极矩分别为(1.90，0.60，0.40)和(0.40，2.30，0.00)的偶极子，放置在两个孤立的位置(-73.00，-51.00，31.67)(mm)和(-43.00，-41.00，61.67)(mm)。头表观测电位是通过BEM正演计算，并施加不同NSR的高斯噪声后获得。对不同的NSR，分别采用LORETA，FOCUSS和LPISS对源进行定位。在本工作中，NSR定义为噪声标准差和信号标准差的比值。计算E_energy，E_localization和NBI时采用的ROI的半径都是25mm，不同的是计算E_energy，E_localization时的球形ROI的中心在模拟源位置，计算NBI时的球形ROI的中心在估计源(最大能量的偶极子)位置。在NSR＝15％时的定位结果显示在附图1中，附图2显示了在所有模拟噪声情况下的定位性能指数。

其次测试了算法对深浅混合源的定位结果。在本模拟实验中，三个偶极子的偶极矩强度分别为(6.00，2.70，1.90)，(-8.00，3.00，1.00)和(5.80，1.00，2.00)，将它们分别放置在三个离散的位置(-83.00，19.00，31.67)(mm)，(-23.00，29.00，-8.33)(mm)和(-63.00，39.00，-48.33)(mm)，其中第一个和第三个源是两个浅源，第二个是深源。同样，也是用BEM正演方法获得头表电位分布，并给其施加NSR＝15％的高斯噪声。也是利用LORETA、FOCUSS和LPISS三种方法对其进行定位，按上面同样的方法计算各自相应的E_localization、E_energy和NBI指数。结果显示在附图3和表1中。

                                表1.在NSR＝15％下三个源定位结果的性能指标

	Source 1	Source 2	Source 3
				LORETA    FOCUSS   LPISS	LORETA    FOCUSS   LPISS	LORETA   FOCUSS   LPISS
	Eenergy	0.8980    0.4419   0.2502	0.9749    0.7566   0.1748				0.9152   0.2192   0.3721
Elocalization(mm)				10.000    0        0	22.3607   22.3607  0	17.3205  0        10
	NBI	0.5784    0.1045   0.1378	0.6311    0.0732   0.1189				0.4794   0.1178   0.2574

从上面的定位结果可以看出，LPISS的E_energy是最小的，E_localization的整体性能优于LORETA和FOCUSS，NBI性能与FOCUSS相当，优于LORETA。这些结果显示了LIPSS良好的稀疏源定位能力，且在与当前代表性的方法(LORETA，FOCUSS)的比较中，显示最好的综合性能。

2.IOR脑电数据的源定位结果：

方法：对Posner设计的经典IOR实验进行修改后(主要是刺激的修改)，采用修改后的实验刺激被试来获取和IOR有关的脑电数据。脑电数据用128道的EGI系统在250Hz的采样率下记录。采用vertex(Cz)电极作为参考。记录后，离线处理丢弃包含有过大的眼电伪迹的数据。在15个被试中，13个被试的数据符合记录要求。根据线索和靶呈现位置的不同组合(RR：right cue-location and righttarget-location；RL：right cue-location and left target-location；LL：Leftcue-location and left target-location；LR：Left cue-location and righttarget-location)，将被试的记录数据段(Epoch)分为4类：RR，RL，LL和LR。每一个数据段长1.2秒，起始于靶出现的前200ms，持续到靶出现后的1000ms。13个有效被试的数据按照线索和靶位置的组合进行总体平均，获得4类诱发电位：RR，RL，LL和LR。在本工作中，展示的是对RR诱发电位进行的源定位结果。

在IOR实验中，由靶刺激诱发的ERP在靶出现后200ms时有一个明显的负波，主要针对该时刻的成分进行源定位测试。在定位时，也是采用前面的标准真实头模型，将EGI的电极坐标配准到标准头模型上，采用BEM计算传递矩阵A，它的维数为128×2730，然后采用LPISS对200ms时刻的数据进行源定位。LPISS的定位结果如附图4所示。

IOR是一个复杂的脑功能区活动过程，当前已有的研究显示IOR是由很多脑功能区域共同作用的结果。在利用LPISS的定位结果中，在右额叶眼区(right frontal eye field，RFEF)、丘脑(thalamus)，右小脑(right cerebellum)区域和左枕叶(left occipital)部位都定位到了较强的活动，和当前IOR研究中利用fMRI等其他方法定位到的基本活动区域是大致相符合的。

Claims (4)

1、基于稀疏性约束的脑电源定位方法，包括以下步骤：

1)确定传递矩阵A；

2)通过多道脑电信号记录系统获取实际记录的脑电信号Y，进行预处理，确定要进行源分析的时刻；

3)初始化源向量x_k-1，k＝1，迭代终止误差ε和一个最大的迭代次数；

4)更新权重对角矩阵：W_k＝(diag(x_k-1))；

5)加入稀疏性约束，利用优化方法直接估计q_k：arg min‖b-AW_kq_k‖2+λ‖q_k‖_p；

6)更新源信息向量：x_k＝W_kq_k；

7)迭代终止条件判断：比较更新前后的源分布变化，当‖x_k-x_k-1‖≤ε或迭代次数超出设定的最大迭代次数限制时，迭代终止，x_k即为源的最终定位估计结果；否则，k＝k+1，转步骤2)，继续迭代；

其中，λ≠0，W_k≠I。

2、如权利要求1所述的基于稀疏性约束的脑电源定位方法，其特征在于，所述步骤1中，包括以下分步骤：

a1)、对待测对象的头部进行MRI或CT扫描，获取头部解剖结构的影像信息；

a2)、提取步骤a1)所得的影像信息中的大脑部分，然后对大脑分割，再提取大脑部分的源功能区(主要包括灰质、海马、小脑等部位)；

a3)、以一定精度的网格将步骤a2)所得的脑源功能区进行剖分，确定解空间网格(包括解空间的维数和各个网格的空间位置序号)；

a4)、确定多道脑电信号记录系统各个电极的空间位置信息；

a5)、确定脑电源的模型；

a6)、利用步骤a3)到步骤a5)中确定的解空间网格、电极位置信息和脑电源模型，利用正演方法计算传递矩阵A，具体方法如下：在每个解空间位置上放置单位的源，利用数值计算方法计算该单位源在电极位置处产生的电位分布，该电位分布构成传递矩阵中的一列，以此类推，当把所有解空间遍历放置单位源后，就可以获得传递矩阵A。

3、如权利要求1所述的基于稀疏性约束的脑电源定位方法，其特征在于，p＝1。

4、如权利要求1所述的基于稀疏性约束的脑电源定位方法，其特征在于，所述步骤2)中的预处理包括：滤波、基线校正、眼电去除。

Images