基于四阶累积量矩阵子空间分解的脑电逆问题算法
作者：尧德中
单位：尧德中(电子科技大学 自动化系，成都 610054)
关键词：脑电图；逆问题；四阶累积量；特征分解算法
生物医学工程学杂志000215 摘要 根据头表观测电位反演脑电源的空间位置是脑电研究中的一个重要问题。基于四阶累积量矩阵的子空间分解，提出了一种新的脑电逆问题算法。与现有的基于二阶统计量矩阵的子空间分解算法相比，二阶方法存在对空间相关噪音敏感的问题，而四阶方法可以抑制在时间方向上呈高斯分布，在空间方向相关的噪音，从而显示了四阶方法的较好性能。
Electroencephalography Inverse Problem by Subspace Decomposition of the Fourth-order Cumulant Matrix
, 百拇医药
Yao Dezhong
(Department of Automation University of Electronic Science and Technology, Chengdu 610054)
Abstract It is an important topic in electroencephalography (EEG) research to localize the EEG activity sources from the scalp recordings. In this paper, based on the fourth-order cumulant matrix, a new sub-space decomposition algorithm is proposed for the EEG inverse problem. As the second-order moments (cumulants) has the drawback of being sensitive to the noise covariance. Using the fourth-order cumulants we need not know the noise covariances, as long as the noise is Gaussian. Computer simulation study on a three-layer concentric sphere head model shows its better performance than the two-order cumulate method in depressing the spatial coherent Gaussian noise.
, 百拇医药
Key words Electroencephalography (EEG) Inverse problem Fourth-order cumulant Subspace decomposition
1 引 言
脑功能的无损测量在当前引起了人们广泛的兴趣。这种测量通常采用不同的感觉刺激，以研究人脑某些特定的高级认知功能，或脑功能的损伤。脑功能的无损测量手段包括脑电(EEG)，脑磁(MEG)和功能磁共振(fMRI)等，其中EEG的费用远低于其它方法，且具有毫秒级的时间分辨率，因而一直受到相关学者的高度重视。
脑电逆问题的研究方法大致可以分为两类[1]：一类是成像方法[2，3]，具体有基于等效分布源、基于边界元、有限元和基于球谐分析皮层电位成像方法，以及假设源呈三维体分布的低分辨层析方法；另一类是偶积源定位方法[4～6]，具体有非线性最小二乘法[4]，全局优化算法[5]和子空间分解算法或称多信号分类算法[6，7]。在这些方法中，大都是以某一时刻的电位观测值为已知信息，唯有子空间分解算法是直接建立在一段观测记录上的，从而较好的同时利用了观测记录中的时间和空间信息[6，7]，受到人们较广泛的重视。
, 百拇医药
Mosher等[6]首先将子空间分解算法用于脑磁逆问题的研究，Sekihara等[7]将其进一步用于同时进行了“任务实验”和“控制实验”的MEG逆问题研究，Koles等[8]将其与共同空间图像分解相结合，研究了自发脑电中癫痫波的定位。在这些工作中，子空间分解是建立在观测记录的二阶统计量-互相关阵的基础上的，其效果受空间相关噪音影响较为明显[6～8]。鉴于高阶统计量对任何形式的高斯过程的不敏感性，本作者发展了基于四阶累积量矩阵的子空间分解脑电逆问题方法，并在数值实验中显示了一定的效果。
2 算 法
2.1 四阶累积量的基本概念
脑电记录是一个典型的多传感信号。假设M个空间实信号构成向量X
XT＝[X1(t),X2(t),…，XM(t)](1)
，空间一、二、三和四阶矩为[9，10]
, 百拇医药
μ1(k1)=E{Xk1(t)}
μ2(k1,k2)=E{Xk1(t)Xk2(t)}
μ3(k1,k2,k3)=E{Xk1(t)Xk2(t)Xk3(t)}
μ4(k1,k2,k3,k4)=E{Xk1(t)Xk2(t)Xk3(t)Xk4(t)} (2)
, 百拇医药
显然，空间二阶矩即空间互相关函数，空间一阶矩即数学期望。
空间互四阶累积量为[10，11]
k4(k1,k2,k3,k4)=μ4(k1,k2,k3,k4)-(μ3(k2,k3,k4)μ1(k1)+μ3(k1,k3,k4)μ1(k2)+μ3(k1,k2,k4)μ1(k3)+μ3(k1,k2,k3)μ1(k4))-
, 百拇医药
μ2(k1,k4)μ2(k2,k3)+μ2(k1,k3)μ2(k2,k4)+μ2(k1,k2)μ2(k3,k4))+2(μ2(k3,k4)μ1(k1)μ1(k2)+
μ2(k2,k4)μ1(k1)μ1(k3)+μ2(k2,k3)μ1(k1)μ1(k4)+μ2(k1,k4)μ1(k2)μ1(k3)+μ2(k1,k3)μ1(k2)μ1(k4)+
, http://www.100md.com
μ2(k1,k2)μ1(k3)μ1(k4))-6μ1(k1)μ1(k2)μ1(k3)μ1(k4) (3)
对于零均值信号，上式简化为
k4(k1,k2,k3,k4)=μ4(k1,k2,k3,k4)-μ2(k1,k4)μ2(k2,k3)-
, 百拇医药
μ2(k1,k3)μ2(k2,k4)-μ2(k1,k2)μ2(k3,k4) (4)
(4)式是空间谱估计中常用的四阶累积量表达式[9，11]。对于脑电信号，零均值假设不一定符合实际，尤其对于诱发脑电信号，为此，在本工作中将以上面(3)式的一般形式作出发点。
累积量有许多优良性质可以利用，下面列出有关的三项性质[10]：
(1)若{αi}ni=1为常量，{xi}ni=1为随机变量，则累积量满足以下关系
(5)
, http://www.100md.com
(2)若随机变量{yi}ni=1与随机变量{xi}ni=1之间统计独立，则
cum(x1+y1，…，xn+yn)=cum(x1,…,xn)+
cum(y1，…，yn) (6)
式中：cum(cumulant)表示累积量。
(3)四阶累积量完全抑制高斯噪音，若{xi}ni=1为高斯随机变量，则
, http://www.100md.com
cum(xk1,xk2,xk3，xk4)=0.0 (7)
2.2 脑电逆问题模型
脑电问题可归纳为以下的模型：
(8)
式中：y(t)∈RM×1为头表M个电极的观测记录；x(t)∈RL×1为头内L个独立源的强度函数；e(t)∈RM×1为M个电极处的噪音。A(α1,…，αj，…，αL)为源到电极的传递函数矩阵，它与头模型有关，本工作中，我们采用的是同心的三层球头模型[12]。
一般地，可假设脑电模型中的噪音为高斯噪音，以前基于二阶统计量的子空间分解算法要求噪音是道间独立的，在空间相关噪音(空间色噪)情况下，该类算法的性能将大大降低[8，9，11]。上述四阶累积量的性质表明((7)式)，基于四阶的算法能完全抑制高斯噪音。即只要噪音的时间过程是高斯的，无论它在道间是独立的还是相关的，都不会对基于四阶的算法产生影响。
, http://www.100md.com
根据四阶累积量的性质和上述脑电模型，我们有
(9)
此式表明，四阶累积量可完全抑制高斯噪音。
假设k4(k1,k2,k3,k4)位于一个M2×M2矩阵C的(k1-1)M+k2行(k3-1)M+k4列， 则(9)式可用矩阵表示为
C＝BSBT (10)
式中：
(11)
, 百拇医药
式中：为Kronecker积；T为转置；Diag{}为以括号内元素为对角元素的对角阵。
对(8)式表示的模型，我们列出以下二个限制条件：
条件1： M>L即观测电极数多于未知独立源的数目，电极位置互不重叠，因而在传递矩阵A(α1，…，αj，…，αL)中，当i不等于j时，αj与αi相互独立，从而有rank(A)=L，rank(B)=L。
条件2： 源的时间过程是非高斯的，因而rank(S)=L,rank(BSBT)=L。
对C作奇异值分解得
(12)
, 百拇医药
其中：
σ2i为矩阵C的本征值，且有λi>0，i=1，…，L；λi≈0，i=L+1，…，M2。这里假设与噪音对应的本征值为零，对于实际情况，它可能是一个小量。正交矩阵U和V由相应的本征矢量构成，根据本征值的分布情况，可将正交矩阵U或V作信号与噪音子空间分解，如U＝[S｜G]＝[u1，…，uL｜uL+1，…，uM2]，其中S＝[u1，…，uL]和G＝[uL+1，…，uM2]分别张成信号子空间和噪音子空间。显然
CG＝BSBTG＝0 (13)
, 百拇医药
因而
GTBSBTG＝(BTG)TS(BTG)＝0 (14)
又由于S是非奇异的，从而有[9]
BTG=0；GTB=0；BTGGTB=0 (15)
(15)式表明
(16)
根据(16)式，子空间分解算法的具体实现为，扫描解空间的每一点i，计算以下目标函数或其倒数：
(17)
, 百拇医药
在源实际存在的地方，该目标函数将趋于零，在其他地方则不趋于零。因此，可以将该目标函数据的三维分布图上的极值点作为对真实源位置的定位。
对于脑电逆问题，位于某个位置的偶积子在整个观测时段内的偶积矩取向可能是固定的，也可能是旋转的。参照文献[6]中二阶子空间分解方法，对于偶积矩可旋转的情况，可把目标函数或称定位算子定义为：
(18)
式中：F为矩阵的Frobenious模；Ai为位于体内位置i处，沿三个基矢量方向的单位偶积子总的传递函数矩阵。对于偶积矩取向是固定的情况，设
为该固定取向偶积子的单位取向矢量，此时的定位算子可定义为：
(19)
, 百拇医药
(19)式最后一个等式的分母表示取大括号内矩阵束的最小 广义本征值。在本仿真研究中，我们只考虑了偶积矩取向是固定的情况，对于取向可为旋转的情况，可用(18)式处理。
综上所述，脑电逆问题的四阶累积量阵子空间分解算法由以下几个步骤构成：
(1)求观测记录的四阶累积量C；
(2)对C作奇异值分解(SVD)得C＝U∑VT；
(3)根据奇异值∑的分布情况和一些生理知识完成信、噪空间的分解，即对矩阵进行分割。在本仿真研究中，我们认为L是已知的，对于实际应用，如何有效地确定L还有待深入的研究。令U＝[S｜G]，其中S由对应∑前面L个较大奇异值的本征向量构成，G由后面M2-L个本征向量构成。对每一个可能的源位置计算J(i)，根据J(i)分布即得到源的位置信息，在实际上就是把J(i)分布的前L个较大值解释为真实源位置。
, 百拇医药
3 数值计算结果
3.1 仿真方法描述
本研究采用的三层同心球头模型的参数为[12]，头皮半径1，颅骨半径0.92，脑半径0.87。有关的正演计算是根据相应的球谐级数解来进行的。脑电活动的时间过程用正弦函数模拟，通过调节正弦函数的频率和初相位，即得到相似但不相同的时间过程，从而可以方便地构造出非高斯且相互独立的源组合。正弦函数的形式取为
f(t)=sin(2πfmt+v)=sin(2π/Nsi+j) i=1，…，N (20)
如果dt=4ms，则Ns=13、14对应的频率为19.2 Hz和17.8 Hz，在本工作中我们取j=0，π＝3.14159。
由于四阶累积量方法的矩阵大小为M2×M2，M太大将导致计算上的困难，在本工作中，我们取M＝11，它们是近似均匀地分布于上半头表的。从头内任意一点到球头表面的电位正演计算公式见文献[12]。研究中所用的噪音为独立同分布的高斯白噪音。对应M道记录的白噪音可以是分别运行M次白噪程序得到的，此时它们在道间是统计独立的，即下面所说的空间不相关噪音。如果是把同一段白噪数据用于所有的M道，即构成空间相关噪音。
, 百拇医药
本研究中的信噪比定义为所有M道噪音数据的标准差与所有M道号的标准差之比，因此，具体到各道的信噪比是不同的，信号强的道，信噪比高，这与实际情况是类似的。
以下计算结果是以二、四阶方法相对照的形式给出的，二阶方法的详细情况请参见文献[6-9]。
3.2 计算结果
3.2.1 空间不相关噪音的影响 假设有单位强度的偶积源1和2位于球头模型中(-0.7，0.4，0.3)和(-0.7，-0.4，0.3)，其径向坐标为0.86，因而它们是近似处于皮层上的，偶积矩呈径向。信噪比为40%，M＝11，NS＝13、14，N＝500。图1和图2分别给出了四阶和二阶方法的结果。它们表明，在空间(道间)不相关噪音的情况下，四阶和二阶方法具有相似的效果。
, 百拇医药
图1 空间不相干噪音情况下四阶累积量法重建结果
图上由小到大的圆圈分别为球头模型的上半球的不同横断面。假设坐标原点位于球心，X轴为从左向右，Z轴向上，则这些横端面所对应的位置为：从上面的左面向右向下分别为Z＝0.8、0.7、0.6、0.5、0.4、0.3、0.2、0.1、0.0
Fig 1 The reconstructed result based on fourth-order cumulant when spatial incoherent noise exists
The circles in the graphs represent respectively different cross sections of the upper part of the three concentric sphere head model. The Cartesian coordinate origin is at the center of the model, the positive x-axis direction is from left to right, and z-axis from bottom to top of the model,that is, z-axis was vertical to the circles plane. The positions of the croos sections respectively is at z=0.8，0.7，0.6、0.5、0.4、0.3、0.2、0.1、0.0, whose order is from top to bottom and left to right in each row
, 百拇医药
3.2.2 空间相关噪音的影响 在其它条件与上面完全相同的情况下，只是加在每一道上的噪音为同一噪音数据，因而为空间相关噪音。图3和图4分别给出了四阶和二阶方法的结果。它们表明，在空间(道间)相关噪音的情况下，二阶方法导致了一个源的丢失，而四阶方法仍较好地实现了对二个源位置的成像。从而用数值算例证明了理论预言的两者在性能上的差异。
图2 空间不相干噪音情况下二阶方法重建结果
Fig 2 The reconstructed result based on second-order cumulant when spatial incoherent noise exists
图3 空间相干噪音情况下四阶累积量法重建结果
, 百拇医药
Fig 3 The reconstructed result based on fourth-order cumulant when spatial coherent noise exists
4 结束语
脑电逆问题是脑认知科学研究和脑功能病变损伤定位研究的一个重要方面，目前虽然发表了许多研究成果，但由于脑电逆问题本身的非线性、奇异性和复杂性，有关的实际应用研究表明，这一问题还远远没有解决。脑电逆问题的子空间分解算法，以正在迅速发展的现代谱估计技术为依托，具有很大的发展空间。本研究表明，子空间分解算法利用11个电极记录就可能完成对两个源的成像，因而具有相当的实际意义。与现有的二阶方法相比，作者提出的基于四阶累积量的方法在克服空间相关噪音方面取得了一定的进展。但是，如何将这一方法推向实际还有大量的工作要做，例如，在本仿真研究中，我们假设信号为正弦信号，因而是非高斯信号，而噪音的时间过程为高斯分布，从而发挥了四阶方法的优势，但要把本方法用于实际问题，就需要明确实际脑电信号的非高斯程度等，这方面的工作需要与具体的应用问题结合进行研究，也是我们下一步努力的方向。
, 百拇医药
图4 空间相干噪音情况下二阶方法重建结果
Fig 4 The reconstructed result based on second-order cumulant when spatial coherent noise exists
国家自然科学基金资助项目(39770215，3998009)
参考文献
1，尧德中，孟繁盛，胡 晓.脑电逆问题研究进展.电子科技导报，1998；11∶31
2，Yao D, Luo B. Theory of the EEG coritcal imaging techniques. Chinese J Biomed Eng (English edition), 1996;5(3)∶128
, http://www.100md.com
3，Yao D. Imaging of EEG by spherical harmonic analysis. Journal of Electronics, 1996;13(1)∶82(“电子科学学刊”中文版，1995；17(5)∶456)
4，张 彤，杨福生，唐庆玉.利用脑电空间高频分量分离深浅源偶积子.电子学报，1996；24(7)∶83
5，Uutela K, Hamalainen M, Salmelin R. Global optimization in the localization of neuromagnetic sources. IEEE Trans BME, 1998;45(6)∶716
6，Mosher JC, Lewis PS, Leahy RM. Multiple dipole modeling and localization from spatio-temporal MEG data. IEEE Trans BME,1992;39∶541
, http://www.100md.com
7，Sekilhara K, Poeppel D, Marantz A et al.MEG covariance difference analysis:a method to extract target source activities by using task and control measurements. IEEE Trans BME,1998;45(1)∶87
8，Koles ZJ. Trends in EEG source localization. Electroenceph Clin Neurophysiol 1998;106∶127
9，张贤达.现代信号处理.北京：清华大学出版社，1995∶70-72,125-128
10，Mendel JM. Tutorial on higher-order statistics (Spectra) in signal processing and system theory: theoretical results and some applications. Proceedings of the IEEE, 1991;79(3)∶278
11，Porat B, Friedlander B. Direction finding algorithm based on high-order statistics. IEEE Trans SP, 1991;39(9)∶2016
12，Rush S, Driscoll DA. EEG electrode sensitivity-An application of reciprocity. IEEE Trans BME, 1969;16∶15
(收稿：1999-04-05), 百拇医药
单位：尧德中(电子科技大学 自动化系，成都 610054)
关键词：脑电图；逆问题；四阶累积量；特征分解算法
生物医学工程学杂志000215 摘要 根据头表观测电位反演脑电源的空间位置是脑电研究中的一个重要问题。基于四阶累积量矩阵的子空间分解，提出了一种新的脑电逆问题算法。与现有的基于二阶统计量矩阵的子空间分解算法相比，二阶方法存在对空间相关噪音敏感的问题，而四阶方法可以抑制在时间方向上呈高斯分布，在空间方向相关的噪音，从而显示了四阶方法的较好性能。
Electroencephalography Inverse Problem by Subspace Decomposition of the Fourth-order Cumulant Matrix
, 百拇医药
Yao Dezhong
(Department of Automation University of Electronic Science and Technology, Chengdu 610054)
Abstract It is an important topic in electroencephalography (EEG) research to localize the EEG activity sources from the scalp recordings. In this paper, based on the fourth-order cumulant matrix, a new sub-space decomposition algorithm is proposed for the EEG inverse problem. As the second-order moments (cumulants) has the drawback of being sensitive to the noise covariance. Using the fourth-order cumulants we need not know the noise covariances, as long as the noise is Gaussian. Computer simulation study on a three-layer concentric sphere head model shows its better performance than the two-order cumulate method in depressing the spatial coherent Gaussian noise.
, 百拇医药
Key words Electroencephalography (EEG) Inverse problem Fourth-order cumulant Subspace decomposition
1 引 言
脑功能的无损测量在当前引起了人们广泛的兴趣。这种测量通常采用不同的感觉刺激，以研究人脑某些特定的高级认知功能，或脑功能的损伤。脑功能的无损测量手段包括脑电(EEG)，脑磁(MEG)和功能磁共振(fMRI)等，其中EEG的费用远低于其它方法，且具有毫秒级的时间分辨率，因而一直受到相关学者的高度重视。
脑电逆问题的研究方法大致可以分为两类[1]：一类是成像方法[2，3]，具体有基于等效分布源、基于边界元、有限元和基于球谐分析皮层电位成像方法，以及假设源呈三维体分布的低分辨层析方法；另一类是偶积源定位方法[4～6]，具体有非线性最小二乘法[4]，全局优化算法[5]和子空间分解算法或称多信号分类算法[6，7]。在这些方法中，大都是以某一时刻的电位观测值为已知信息，唯有子空间分解算法是直接建立在一段观测记录上的，从而较好的同时利用了观测记录中的时间和空间信息[6，7]，受到人们较广泛的重视。
, 百拇医药
Mosher等[6]首先将子空间分解算法用于脑磁逆问题的研究，Sekihara等[7]将其进一步用于同时进行了“任务实验”和“控制实验”的MEG逆问题研究，Koles等[8]将其与共同空间图像分解相结合，研究了自发脑电中癫痫波的定位。在这些工作中，子空间分解是建立在观测记录的二阶统计量-互相关阵的基础上的，其效果受空间相关噪音影响较为明显[6～8]。鉴于高阶统计量对任何形式的高斯过程的不敏感性，本作者发展了基于四阶累积量矩阵的子空间分解脑电逆问题方法，并在数值实验中显示了一定的效果。
2 算 法
2.1 四阶累积量的基本概念
脑电记录是一个典型的多传感信号。假设M个空间实信号构成向量X
XT＝[X1(t),X2(t),…，XM(t)](1)
，空间一、二、三和四阶矩为[9，10], 百拇医药
μ1(k1)=E{Xk1(t)}
μ2(k1,k2)=E{Xk1(t)Xk2(t)}
μ3(k1,k2,k3)=E{Xk1(t)Xk2(t)Xk3(t)}
μ4(k1,k2,k3,k4)=E{Xk1(t)Xk2(t)Xk3(t)Xk4(t)} (2)
, 百拇医药
显然，空间二阶矩即空间互相关函数，空间一阶矩即数学期望。
空间互四阶累积量为[10，11]
k4(k1,k2,k3,k4)=μ4(k1,k2,k3,k4)-(μ3(k2,k3,k4)μ1(k1)+μ3(k1,k3,k4)μ1(k2)+μ3(k1,k2,k4)μ1(k3)+μ3(k1,k2,k3)μ1(k4))-
, 百拇医药
μ2(k1,k4)μ2(k2,k3)+μ2(k1,k3)μ2(k2,k4)+μ2(k1,k2)μ2(k3,k4))+2(μ2(k3,k4)μ1(k1)μ1(k2)+
μ2(k2,k4)μ1(k1)μ1(k3)+μ2(k2,k3)μ1(k1)μ1(k4)+μ2(k1,k4)μ1(k2)μ1(k3)+μ2(k1,k3)μ1(k2)μ1(k4)+
, http://www.100md.com
μ2(k1,k2)μ1(k3)μ1(k4))-6μ1(k1)μ1(k2)μ1(k3)μ1(k4) (3)
对于零均值信号，上式简化为
k4(k1,k2,k3,k4)=μ4(k1,k2,k3,k4)-μ2(k1,k4)μ2(k2,k3)-
, 百拇医药
μ2(k1,k3)μ2(k2,k4)-μ2(k1,k2)μ2(k3,k4) (4)
(4)式是空间谱估计中常用的四阶累积量表达式[9，11]。对于脑电信号，零均值假设不一定符合实际，尤其对于诱发脑电信号，为此，在本工作中将以上面(3)式的一般形式作出发点。
累积量有许多优良性质可以利用，下面列出有关的三项性质[10]：
(1)若{αi}ni=1为常量，{xi}ni=1为随机变量，则累积量满足以下关系
(5), http://www.100md.com
(2)若随机变量{yi}ni=1与随机变量{xi}ni=1之间统计独立，则
cum(x1+y1，…，xn+yn)=cum(x1,…,xn)+
cum(y1，…，yn) (6)
式中：cum(cumulant)表示累积量。
(3)四阶累积量完全抑制高斯噪音，若{xi}ni=1为高斯随机变量，则
, http://www.100md.com
cum(xk1,xk2,xk3，xk4)=0.0 (7)
2.2 脑电逆问题模型
脑电问题可归纳为以下的模型：
(8)式中：y(t)∈RM×1为头表M个电极的观测记录；x(t)∈RL×1为头内L个独立源的强度函数；e(t)∈RM×1为M个电极处的噪音。A(α1,…，αj，…，αL)为源到电极的传递函数矩阵，它与头模型有关，本工作中，我们采用的是同心的三层球头模型[12]。
一般地，可假设脑电模型中的噪音为高斯噪音，以前基于二阶统计量的子空间分解算法要求噪音是道间独立的，在空间相关噪音(空间色噪)情况下，该类算法的性能将大大降低[8，9，11]。上述四阶累积量的性质表明((7)式)，基于四阶的算法能完全抑制高斯噪音。即只要噪音的时间过程是高斯的，无论它在道间是独立的还是相关的，都不会对基于四阶的算法产生影响。
, http://www.100md.com
根据四阶累积量的性质和上述脑电模型，我们有
(9)此式表明，四阶累积量可完全抑制高斯噪音。
假设k4(k1,k2,k3,k4)位于一个M2×M2矩阵C的(k1-1)M+k2行(k3-1)M+k4列， 则(9)式可用矩阵表示为
C＝BSBT (10)
式中：
(11), 百拇医药
式中：为Kronecker积；T为转置；Diag{}为以括号内元素为对角元素的对角阵。
对(8)式表示的模型，我们列出以下二个限制条件：
条件1： M>L即观测电极数多于未知独立源的数目，电极位置互不重叠，因而在传递矩阵A(α1，…，αj，…，αL)中，当i不等于j时，αj与αi相互独立，从而有rank(A)=L，rank(B)=L。
条件2： 源的时间过程是非高斯的，因而rank(S)=L,rank(BSBT)=L。
对C作奇异值分解得
(12), 百拇医药
其中：
σ2i为矩阵C的本征值，且有λi>0，i=1，…，L；λi≈0，i=L+1，…，M2。这里假设与噪音对应的本征值为零，对于实际情况，它可能是一个小量。正交矩阵U和V由相应的本征矢量构成，根据本征值的分布情况，可将正交矩阵U或V作信号与噪音子空间分解，如U＝[S｜G]＝[u1，…，uL｜uL+1，…，uM2]，其中S＝[u1，…，uL]和G＝[uL+1，…，uM2]分别张成信号子空间和噪音子空间。显然CG＝BSBTG＝0 (13)
, 百拇医药
因而
GTBSBTG＝(BTG)TS(BTG)＝0 (14)
又由于S是非奇异的，从而有[9]
BTG=0；GTB=0；BTGGTB=0 (15)
(15)式表明
(16)根据(16)式，子空间分解算法的具体实现为，扫描解空间的每一点i，计算以下目标函数或其倒数：
(17), 百拇医药
在源实际存在的地方，该目标函数将趋于零，在其他地方则不趋于零。因此，可以将该目标函数据的三维分布图上的极值点作为对真实源位置的定位。
对于脑电逆问题，位于某个位置的偶积子在整个观测时段内的偶积矩取向可能是固定的，也可能是旋转的。参照文献[6]中二阶子空间分解方法，对于偶积矩可旋转的情况，可把目标函数或称定位算子定义为：
(18)式中：F为矩阵的Frobenious模；Ai为位于体内位置i处，沿三个基矢量方向的单位偶积子总的传递函数矩阵。对于偶积矩取向是固定的情况，设
为该固定取向偶积子的单位取向矢量，此时的定位算子可定义为： (19), 百拇医药
(19)式最后一个等式的分母表示取大括号内矩阵束的最小 广义本征值。在本仿真研究中，我们只考虑了偶积矩取向是固定的情况，对于取向可为旋转的情况，可用(18)式处理。
综上所述，脑电逆问题的四阶累积量阵子空间分解算法由以下几个步骤构成：
(1)求观测记录的四阶累积量C；
(2)对C作奇异值分解(SVD)得C＝U∑VT；
(3)根据奇异值∑的分布情况和一些生理知识完成信、噪空间的分解，即对矩阵进行分割。在本仿真研究中，我们认为L是已知的，对于实际应用，如何有效地确定L还有待深入的研究。令U＝[S｜G]，其中S由对应∑前面L个较大奇异值的本征向量构成，G由后面M2-L个本征向量构成。对每一个可能的源位置计算J(i)，根据J(i)分布即得到源的位置信息，在实际上就是把J(i)分布的前L个较大值解释为真实源位置。
, 百拇医药
3 数值计算结果
3.1 仿真方法描述
本研究采用的三层同心球头模型的参数为[12]，头皮半径1，颅骨半径0.92，脑半径0.87。有关的正演计算是根据相应的球谐级数解来进行的。脑电活动的时间过程用正弦函数模拟，通过调节正弦函数的频率和初相位，即得到相似但不相同的时间过程，从而可以方便地构造出非高斯且相互独立的源组合。正弦函数的形式取为
f(t)=sin(2πfmt+v)=sin(2π/Nsi+j) i=1，…，N (20)
如果dt=4ms，则Ns=13、14对应的频率为19.2 Hz和17.8 Hz，在本工作中我们取j=0，π＝3.14159。
由于四阶累积量方法的矩阵大小为M2×M2，M太大将导致计算上的困难，在本工作中，我们取M＝11，它们是近似均匀地分布于上半头表的。从头内任意一点到球头表面的电位正演计算公式见文献[12]。研究中所用的噪音为独立同分布的高斯白噪音。对应M道记录的白噪音可以是分别运行M次白噪程序得到的，此时它们在道间是统计独立的，即下面所说的空间不相关噪音。如果是把同一段白噪数据用于所有的M道，即构成空间相关噪音。
, 百拇医药
本研究中的信噪比定义为所有M道噪音数据的标准差与所有M道号的标准差之比，因此，具体到各道的信噪比是不同的，信号强的道，信噪比高，这与实际情况是类似的。
以下计算结果是以二、四阶方法相对照的形式给出的，二阶方法的详细情况请参见文献[6-9]。
3.2 计算结果
3.2.1 空间不相关噪音的影响 假设有单位强度的偶积源1和2位于球头模型中(-0.7，0.4，0.3)和(-0.7，-0.4，0.3)，其径向坐标为0.86，因而它们是近似处于皮层上的，偶积矩呈径向。信噪比为40%，M＝11，NS＝13、14，N＝500。图1和图2分别给出了四阶和二阶方法的结果。它们表明，在空间(道间)不相关噪音的情况下，四阶和二阶方法具有相似的效果。
, 百拇医药
图1 空间不相干噪音情况下四阶累积量法重建结果
图上由小到大的圆圈分别为球头模型的上半球的不同横断面。假设坐标原点位于球心，X轴为从左向右，Z轴向上，则这些横端面所对应的位置为：从上面的左面向右向下分别为Z＝0.8、0.7、0.6、0.5、0.4、0.3、0.2、0.1、0.0
Fig 1 The reconstructed result based on fourth-order cumulant when spatial incoherent noise exists
The circles in the graphs represent respectively different cross sections of the upper part of the three concentric sphere head model. The Cartesian coordinate origin is at the center of the model, the positive x-axis direction is from left to right, and z-axis from bottom to top of the model,that is, z-axis was vertical to the circles plane. The positions of the croos sections respectively is at z=0.8，0.7，0.6、0.5、0.4、0.3、0.2、0.1、0.0, whose order is from top to bottom and left to right in each row
, 百拇医药
3.2.2 空间相关噪音的影响 在其它条件与上面完全相同的情况下，只是加在每一道上的噪音为同一噪音数据，因而为空间相关噪音。图3和图4分别给出了四阶和二阶方法的结果。它们表明，在空间(道间)相关噪音的情况下，二阶方法导致了一个源的丢失，而四阶方法仍较好地实现了对二个源位置的成像。从而用数值算例证明了理论预言的两者在性能上的差异。
图2 空间不相干噪音情况下二阶方法重建结果
Fig 2 The reconstructed result based on second-order cumulant when spatial incoherent noise exists
图3 空间相干噪音情况下四阶累积量法重建结果
, 百拇医药
Fig 3 The reconstructed result based on fourth-order cumulant when spatial coherent noise exists
4 结束语
脑电逆问题是脑认知科学研究和脑功能病变损伤定位研究的一个重要方面，目前虽然发表了许多研究成果，但由于脑电逆问题本身的非线性、奇异性和复杂性，有关的实际应用研究表明，这一问题还远远没有解决。脑电逆问题的子空间分解算法，以正在迅速发展的现代谱估计技术为依托，具有很大的发展空间。本研究表明，子空间分解算法利用11个电极记录就可能完成对两个源的成像，因而具有相当的实际意义。与现有的二阶方法相比，作者提出的基于四阶累积量的方法在克服空间相关噪音方面取得了一定的进展。但是，如何将这一方法推向实际还有大量的工作要做，例如，在本仿真研究中，我们假设信号为正弦信号，因而是非高斯信号，而噪音的时间过程为高斯分布，从而发挥了四阶方法的优势，但要把本方法用于实际问题，就需要明确实际脑电信号的非高斯程度等，这方面的工作需要与具体的应用问题结合进行研究，也是我们下一步努力的方向。
, 百拇医药
图4 空间相干噪音情况下二阶方法重建结果
Fig 4 The reconstructed result based on second-order cumulant when spatial coherent noise exists
国家自然科学基金资助项目(39770215，3998009)
参考文献
1，尧德中，孟繁盛，胡 晓.脑电逆问题研究进展.电子科技导报，1998；11∶31
2，Yao D, Luo B. Theory of the EEG coritcal imaging techniques. Chinese J Biomed Eng (English edition), 1996;5(3)∶128
, http://www.100md.com
3，Yao D. Imaging of EEG by spherical harmonic analysis. Journal of Electronics, 1996;13(1)∶82(“电子科学学刊”中文版，1995；17(5)∶456)
4，张 彤，杨福生，唐庆玉.利用脑电空间高频分量分离深浅源偶积子.电子学报，1996；24(7)∶83
5，Uutela K, Hamalainen M, Salmelin R. Global optimization in the localization of neuromagnetic sources. IEEE Trans BME, 1998;45(6)∶716
6，Mosher JC, Lewis PS, Leahy RM. Multiple dipole modeling and localization from spatio-temporal MEG data. IEEE Trans BME,1992;39∶541
, http://www.100md.com
7，Sekilhara K, Poeppel D, Marantz A et al.MEG covariance difference analysis:a method to extract target source activities by using task and control measurements. IEEE Trans BME,1998;45(1)∶87
8，Koles ZJ. Trends in EEG source localization. Electroenceph Clin Neurophysiol 1998;106∶127
9，张贤达.现代信号处理.北京：清华大学出版社，1995∶70-72,125-128
10，Mendel JM. Tutorial on higher-order statistics (Spectra) in signal processing and system theory: theoretical results and some applications. Proceedings of the IEEE, 1991;79(3)∶278
11，Porat B, Friedlander B. Direction finding algorithm based on high-order statistics. IEEE Trans SP, 1991;39(9)∶2016
12，Rush S, Driscoll DA. EEG electrode sensitivity-An application of reciprocity. IEEE Trans BME, 1969;16∶15
(收稿：1999-04-05), 百拇医药