认知神经科学领域脑电复杂性测度方法的新进展，本文主要内容关键词为：复杂性论文,认知论文,新进展论文,神经论文,领域论文，此文献不代表本站观点，内容供学术参考，文章仅供参考阅读下载。
中图分类号：B84
文献标识码：A
文章编号：1006-6020(2002)-03-0051-05
1 引言
从Berger(1929)发现脑电(electroencephalogram,EEG)开始[1]，脑电信号中有效信息的提取一直是困扰研究者的难题。传统方法主要有脑电地形图(EEG mapping)和谱分析(spectral analysis)两类。脑电地形图只能粗略地描述人在认知加工过程中各脑区的激活程度。在脑电频域和时域特征(frequency and time domain features)分析中，数字信号的线性处理方法已得到广泛应用，如事件相关电位(event-related potential,ERP)。然而实际记录的脑波很难满足线性分析方法的要求（如低信噪化、脑电信号平稳等）[2]，且认知神经科学通常采用的平均叠加法会导致有用信息的大量损失，因此线性分析方法在很大程度上限制了认知电位时空模式研究的发展。
大量研究表明人脑是一个结构和功能高度复杂的系统，而脑电信号是神经细胞生物电活动在时间和空间上的非线性耦合[3]。从80年代中期开始，许多研究者用非线性混沌动力学理论发展了一些脑电信号复杂性测度的算法[4]，如分型维数(fractal dimension)和Lyapunov指数(L-exponential)等[5]。由于这些方法无需作锁时(time-locked)和锁相(phaselocked)处理，在早期的研究中得到了广泛的应用。然而这些方法要求的数据量较大、对取样信号的平稳度要求较高[5]，再者混沌动力学中讨论的对象是混沌吸引子，并且不同的研究者在相似的实验条件下所得到的结果变异较大，脑电信号是否具有低维混沌特性从而受到了质疑[6]，因此上述方法可能并不适合于人脑这种各向异性的空间扩展系统。
随着非线性理论的发展，脑电复杂性测度分析方法进一步得到完善。目前常用的脑电复杂性测度算法主要有K[,c]复杂度（包括K[,c]复杂度及其各种改进算法和信息传输矩阵(Information Transmission Matrix,ITM)和近似熵(Approximate Entropy,ApEn)。它们对脑电信号的取样量及其平稳度的要求较低，且无需考虑其是否具有低维混沌特性，从而成为刻画脑电信号非线性变化特征的有效手段[2]。本文就上述方法、特点及其应用作一简要介绍。
2 基于K[,c]复杂度的分析方法
Kolmogorov(1965)提出用产生给定0、1序列最少的计算机程序的比特数作为序列的复杂性度量，这种刻画序列复杂性的方法称为算法复杂性(Algorithm complexity)[2]。Lempel和Ziv以复制和添加两个简单操作为核心，对序列的复杂性作了进一步描述。他们定义的复杂性是一个时间序列随其长度的增长出现新模式的速率，表现了序列接近随机的程度，能反映一个动力学系统的动态特征[7]。在此基础上Kaspar和Schuster发展了随机序列复杂性测度的算法[8]，Wu等人(1991)则首先将这种算法引入脑电信号的分析中，作为反映大脑信息加工活动的有序程度的指标[9]。
2.1 K[,c]复杂度
k[,c]复杂度的计算步骤如下：
(1)粗粒化预处理(coarse graining preprocessing)。对于一给定序列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后称之“复制”；反之则称为“插入”，即用一个"."把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[,c]趋向于1。
2.2 C[,1]和C[,2]复杂度
D'Alessandro和Politi认为K[,c]复杂度只反映了时间序列的随机化程度，并不能完全反映大脑认知功能复杂性的实质[10]。X[,u]发展了复杂度C[,1]和C[,2]算法[11]。
在时间序列中有长度为n-1的子序列但没有长度为n的子序列(S[,1]S[,2]S[,3]…S[,n])，则称(S[,1]S[,2]S[,3]…S[,n-1]S[,n])为长度为n的禁止字。记N[,f](n)为时间序列中的禁止字数目，那么C[,2]的计算见公式③。
2.3 C[,0]复杂度
K[,c]、C[,1]、C[,2]算法中过粗粒化(over-coarse)的预处理可能会导致原始信号中信息的大量丢失，不恰当的粗粒化甚至会改变原始时间序列的动力学特性，例如，有可能将随机时间序列改变成周期时间序列。为了消除这种潜在的危险，Chen等人定义了一种新的复杂度算法C[,0][12]。
C[,0]复杂度假设任何复杂运动的时间序列都是由规则运动时间序列和随机运动时间序列组成。因此C[,0]复杂度的定义就为时间序列随机运动时序和时间轴所围区域的面积与整个复杂运动时间序列和时间轴所围面积之比，具体的计算步骤如下：
(1)利用快速傅立叶变换(Fast Fourier Transform,FFT)计算原始时间序列的功率谱和平均值；
(2)只有那些振幅比平均值大的波谱成分才被保留，其余的均被置为0；
(3)然后对这个新的波谱进行FFT反转，从而得到一个新的时间序列；将此序列作为原始时间序列的规则成分(regular component)，而原始时间序列与规则成分之差称为无规则成分(disorder component)；
(4)无规则成分的面积与原始时间序列面积的比值记为复杂度C[,0]。
可见周期信号的C[,0]值为0，白噪声(white noise)的C[,0]为1。
2.4 信息传输矩阵ITM
Xu等人根据互信息论(the mutual information theory)提出每一个电极的EEG序列都可以重建一个m维的相空间[2]。在第i个电极处，取一段从t[,0]开始、时间窗长(time window)为1024ms的脑电数据[x[,i]t[,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)H[X[,i](t[,0])]、H[X[,j](t[,0]+k[,τ])]（电极j的t[,0]+k[,τ]的熵）以及联合熵H[X[,i](t[,0]),X[,j](t[,0]+k[,τ])]，其中τ为1ms。因此从第i个电极到第j个电极之间的延迟为k[,τ]的信息传输可以由公式④决定：
IT[,ij](t[,0],k[,τ])=H[X[,i](t[,0]+k[,τ])]
-H[X[,i](t[,0]+k[,τ])]
④
确定t[,0]，且k的取值范围是0到511，得到信息传输的时间序列。用复杂度计算这个时间序列，得到从第i个电极到第j个电极在区间[t[,0],t[,0]+511]之间的信息传输活动程度的特征指标。
由此可以构建一个由n×n=n[2]个值的矩阵，其第i行第j列的值为C[,i]，j(t[,0])。Xu称该矩阵为信息传输矩阵[13]。信息传输矩阵是一种直观表示不同脑皮层之间信息传递量的指标：第i行表示从第i个电极向其他电极位置的信息传递量（包括第i个电极本身），第j列表示在第j个电极处接收到的信息量。以为步长，逐渐增加t[,0]的值，重复上面的步骤，就可以得到一系列信息传输矩阵，以此表征大脑信息传输复杂度的动力学过程。
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)之间的距离d[[X(i),(X(j)]]并将最大值定义为最大反应成分距离，见公式(5)：
D[X(i),(X(j)]=max[｜x(i+k)-x(j+k)｜](k=[0,m-1])
⑤
(3)定义一个阈值r(r＞0)，对于每一个i值，记录满足条件d[X(i),X(j)]＜r的个数。把这个值与N-m的比值定义为C[m,i](r)，见公式⑥：
C[m,i](r)=｛d[X(i),X(j)]＜r的数目｝／(N-m+1)
⑥
(4)对每一个可能的i值，计算C[m,i](r)的对数，求这些对数的平均值，定义为Φ[m](r)，见公式⑦：
可以证明该极限存在且极值为1。因此，ApEn可以表示向量集随着m增大产生新模式的概率，产生新模式的概率越大ApEn值就越大，即时间序列的复杂度越大。实际上，N不可能取无穷大，所以通常只能在N足够大的时候对ApEn进行估计（见公式⑧），而且ApEn的值还依赖于m和r。
ApEn(m,r,N)=Φ[m](r)-Φ[m+1](r)
⑧
根据经验，Pincus建议m取2，r取0.1～0.2倍原始数据的标准差。从而只需用很短时间序列（约1000个数据点）就足以估算出可靠的ApEn值。这种方法特别适用于生物电这类极其不稳定信号的分析。
上述几种复杂度分析方法中，由于K[,c]计算简单，易于理解而应用最为广泛。基于K[,c]发展起来的C[,1]和C[,2]改进了脑电信号出现新模式的检测方法，能够更为深刻地描述脑电动力学系统的复杂性本质。然而由于它们都需要对脑电时间序列进行粗粒化预处理，从而可能丢失脑电信号中有意义的信息。C[,0]强调脑电信号由规则运动部分和随机运动部分组成，其算法避免了粗粒化处理，有助于减少脑电信号中有效信息的损失。基于复杂度的信息传输矩阵提供了描述脑电信息传输量更为直观有效的指标。近似熵适合于研究短数据、抗干扰能力强，且不需对原始数据进行粗粒化处理而成为生物电信号分析的重要方法。然而近似熵需要事前设定m和r两个参数，且相对运算量较大。可见上述方法各具特点，应根据具体的研究目的和实验条件进行合理选择。
4 应用及展望
利用非线性动力学复杂性测度研究脑电信号，其实质是把测定的时间序列的复杂度作为衡量该时间序列所含信息量的指标，分析人在不同状态下脑电信号的时空模式，以揭示脑的认知功能[3]，因此复杂度在认知科学、临床等领域得到了越来越广泛的应用。如大脑成熟度评估、情绪的变化、各种思维方式的对比等都采用了脑电复杂性测度方法[15]。有研究者发现，正常人在不同的状态下脑电复杂度的变化表现出一定的规律性。如在睁眼状态下复杂度高于闭眼状态下的复杂度，在执行任务时额叶大脑活动区域复杂度降低[16]，而对帕金森、精神分裂症病人的研究表明其脑电复杂度的变化趋势与正常人相反[17]。另有研究模拟了高空飞行不同程度的缺氧条件，发现脑电复杂度对脑缺氧十分敏感，可以作为一项临床诊断的指标[18]。此外脑电复杂度也被广泛地应用于麻醉深度测定、中风病人的脑电活动特征分析、encephalopathy、Creutzfeldt-Jakob病症监测和癫痫发作预测等[19,20]。
就目前的研究而言，大多局限于离线(offline)数据的分析，尚未突破实时(on-line)分析的技术难点，如何有效地利用和发展这些方法以获取反映认知加工过程的动态指标将成为方法学研究的重点所在。