认知神经科学领域脑电复杂性测度方法的新进展

1、认知神经科学领域脑电复杂性测度方法的新进展     1引言从Berger(1929)发现脑电(electroencephalogram,EEG)开始1，脑电信号中有效信息的提取一直是困扰研究者的难题。传统方法主要有脑电地形图(EEG  mapping)和谱分析(spectral  analysis)两类。脑电地形图只能粗略地描述人在认知加工过程中各脑区的激活程度。在脑电频域和时域特征(frequency  and  time  domain &#

2、160;features)分析中，数字信号的线性处理方法已得到广泛应用，如事件相关电位(event-related  potential,ERP)。然而实际记录的脑波很难满足线性分析方法的要求（如低信噪化、脑电信号平稳等）2，且认知神经科学通常采用的平均叠加法会导致有用信息的大量损失，因此线性分析方法在很大程度上限制了认知电位时空模式研究的发展。大量研究表明人脑是一个结构和功能高度复杂的系统，而脑电信号是神经细胞生物电活动在时间和空间上的非线性耦合3。从80年代中期开始，许多研究者用非线性混沌动力学理论发展了一些脑电信号复杂性测度的算法4，如分型维数(fractal

3、0; dimension)和Lyapunov指数(L-exponential)等5。由于这些方法无需作锁时(time-locked)和锁相(phaselocked)处理，在早期的研究中得到了广泛的应用。然而这些方法要求的数据量较大、对取样信号的平稳度要求较高5，再者混沌动力学中讨论的对象是混沌吸引子，并且不同的研究者在相似的实验条件下所得到的结果变异较大，脑电信号是否具有低维混沌特性从而受到了质疑6，因此上述方法可能并不适合于人脑这种各向异性的空间扩展系统。随着非线性理论的发展，脑电复杂性测度分析方法进一步得到完善。目前常用的脑电复杂性测度算法主要有K,c复杂度（包括K,c复杂度及其

4、各种改进算法和信息传输矩阵(Information  Transmission  Matrix,ITM)和近似熵(Approximate  Entropy,ApEn)。它们对脑电信号的取样量及其平稳度的要求较低，且无需考虑其是否具有低维混沌特性，从而成为刻画脑电信号非线性变化特征的有效手段2。本文就上述方法、特点及其应用作一简要介绍。2基于K,c复杂度的分析方法Kolmogorov(1965)提出用产生给定0、1序列最少的计算机程序的比特数作为序列的复杂性度量，这种刻画序列复杂性的方法称为算法复杂性(Algorithm &#

5、160;complexity)2。Lempel和Ziv以复制和添加两个简单操作为核心，对序列的复杂性作了进一步描述。他们定义的复杂性是一个时间序列随其长度的增长出现新模式的速率，表现了序列接近随机的程度，能反映一个动力学系统的动态特征7。在此基础上Kaspar和Schuster发展了随机序列复杂性测度的算法8，Wu等人(1991)则首先将这种算法引入脑电信号的分析中，作为反映大脑信息加工活动的有序程度的指标9。2.1K,c复杂度k,c复杂度的计算步骤如下：(1)粗粒化预处理(coarse  graining  preprocessing)。对于一给定序列

6、X=(X,1,X,2,X,n)，首先求得这个序列的平均值，再重构该序列。令大于平均值的X,i为1，小于平均值的X,i为0。将序列(X,1,X,2,X,n)转化为一个字符串形式的0、1序列(s,1,s,2,s,n)。(2)在S=(s,1,s,2,s,m)后加一个或一串字符Q（Q=s,m+1或Q=s,m+1,s,m+2,s,m+k），得到字符串SQ=(s,1,s,2,s,m,s,m+1)或SQ=(s,1,s,2,s,m,s,m+1,s,m+2,s,m+k)，令SQv为SQ减去最后一个字符所得到的字符串。如果Q属于SQv中的“字句”（即两点间的字符串），那么把Q加在S后称之“复制”；反之则称为“插入

7、”，即用一个"."把Q与S前后分开。再把"."前面的所有的字符看成S，重复如上步骤。(3)如上所述，得到用"."分成段的字符串，分成的段的数目就定义为“复杂度”C(n)；(4)根据Lempel和Ziv的研究，对几乎所有的X属于0,1的C(n)都会趋向一个定值b(n)（见公式）。附图以b(n)来对C(n)进行归一化后得到一个相对复杂度c(n)=C(n)/b(n)，称之为Kolmogorov复杂度(K,c)。K,c复杂度反映了时间序列的随机程度，如果时间序列是周期性的，那么K,c就会随时间序列的增加而趋向于0；如果时间序列是随机的，则K,

8、c趋向于1。2.2C,1和C,2复杂度D'Alessandro和Politi认为K,c复杂度只反映了时间序列的随机化程度，并不能完全反映大脑认知功能复杂性的实质10。X,u发展了复杂度C,1和C,2算法11。附图在时间序列中有长度为n-1的子序列但没有长度为n的子序列(S,1S,2S,3S,n)，则称(S,1S,2S,3S,n-1S,n)为长度为n的禁止字。记N,f(n)为时间序列中的禁止字数目，那么C,2的计算见公式。附图2.3C,0复杂度K,c、C,1、C,2算法中过粗粒化(over-coarse)的预处理可能会导致原始信号中信息的大量丢失，不恰当的粗粒化甚至会改变原始时间序列的动

9、力学特性，例如，有可能将随机时间序列改变成周期时间序列。为了消除这种潜在的危险，Chen等人定义了一种新的复杂度算法C,012。C,0复杂度假设任何复杂运动的时间序列都是由规则运动时间序列和随机运动时间序列组成。因此C,0复杂度的定义就为时间序列随机运动时序和时间轴所围区域的面积与整个复杂运动时间序列和时间轴所围面积之比，具体的计算步骤如下：(1)利用快速傅立叶变换(Fast  Fourier  Transform,FFT)计算原始时间序列的功率谱和平均值；(2)只有那些振幅比平均值大的波谱成分才被保留，其余的均被置为0；(3)然后对这个新的波谱进行FF

10、T反转，从而得到一个新的时间序列；将此序列作为原始时间序列的规则成分(regular  component)，而原始时间序列与规则成分之差称为无规则成分(disorder  component)；(4)无规则成分的面积与原始时间序列面积的比值记为复杂度C,0。可见周期信号的C,0值为0，白噪声(white  noise)的C,0为1。2.4信息传输矩阵ITMXu等人根据互信息论(the  mutual  information  theory)提出每一个电极的EEG序列都可以

11、重建一个m维的相空间2。在第i个电极处，取一段从t,0开始、时间窗长(time  window)为1024ms的脑电数据x,it,0,x,i(t,0+1),x,i(t,0+2),x,i(t,0+1023)。据此可以计算向量x,i(t),x,i(t+1),x,i(t,0+2)及其头落在相空间三维子空间中的概率，并且可以计算其熵(entropy)HX,i(t,0)、HX,j(t,0+k,)（电极j的t,0+k,的熵）以及联合熵HX,i(t,0),X,j(t,0+k,)，其中为1ms。因此从第i个电极到第j个电极之间的延迟为k,的信息传输可以由公式决定：  I

12、T,ij(t,0,k,)=HX,i(t,0+k,)-HX,i(t,0+k,) 确定t,0，且k的取值范围是0到511，得到信息传输的时间序列。用复杂度计算这个时间序列，得到从第i个电极到第j个电极在区间t,0,t,0+511之间的信息传输活动程度的特征指标。由此可以构建一个由n×n=n2个值的矩阵，其第i行第j列的值为C,i，j(t,0)。Xu称该矩阵为信息传输矩阵13。信息传输矩阵是一种直观表示不同脑皮层之间信息传递量的指标：第i行表示从第i个电极向其他电极位置的信息传递量（包括第i个电极本身），第j列表示在第j个电极处接收到的信息量。以为步长，逐渐增加t,0的值，重复上面的步骤，

13、就可以得到一系列信息传输矩阵，以此表征大脑信息传输复杂度的动力学过程。3近似熵分析方法近似熵是一种不需要进行粗粒化的脑电复杂性测度分析方法，该方法于1991年由Pincus提出14，并在脑电分析领域中得到广泛应用。ApEn的定义和算法如下：(1)对于一给定的时间序列(1),(2),(N)，按顺序将其组成一个m维的向量集X(i)，即X(i)=(i),(i+1),(i+m-1)(i=1,2,3,N-m+1)；(2)计算向量X(i)与其余向量X(j)之间的距离dX(i),(X(j)并将最大值定义为最大反应成分距离，见公式(5)：  DX(i),(X(j)=maxx(i+k)-x(j+k)(k=0,m-1)(3)定义一个阈值r(r0)，对于每一个i值，记录满足条件dX(i),X(j)r的个数。把这个值与N-m的比值定义为Cm,i(r)，见公式：  Cm,i(r)=dX(i),X(j)r的数目(N-m+1)(4)对每