独立分量生物医学运用探析

1独立分量分析

独立分量分析(ICA)是统计信号处理近年来的一项发展。顾名思义,这是一种分解技术,其特点是把信号分解成若干相互独立的成分。主分量分析(PCA)和奇异值分解(SVD)是人们较熟悉的分解信号的线性代数方法,ICA与它们的主要不同之处表现在:

(1)后者只要求分解出来的各分量互相正交(不相关),但并不要求它们互相独立。用统计信号处理的语言来表达,即:后者只考虑二阶统计特性,而前者则要更全面考虑其概率密度函数的统计独立性。

(2)后者按能量大小排序来考虑被分解分量的重要性。这样的分解虽然在数据压缩和去除弱噪声方面有其优点,但分解结果往往缺乏明确的生理意义。前者虽然分解出的分量其能量大小存在不确定性,但当测量值确实是由若干独立信源混合而成时,分解结果往往具有更好的生理解释。由于测得的生理信号往往是若干独立成分的加权迭加(例如,诱发脑电总是被自发脑电所淹没,而且常伴随有心电、眼动、头皮肌电等干扰),此ICA是一项值得注意的分解方法。

此外,神经生理研究认为,人类对认知、感知信息的前期处理有“去冗余”的特点。ICA在这方面也表现出类似特性,因为互相独立的分量之间互信息是最少的。ICA是伴随着盲信号处理,特别是盲信源分离发展起来。其研究热潮方兴未艾,也正在引起生物医学工程界的注意,IEEETransBME正在组织出版以它为重点的专辑。就国际范围看,以下几个研究单位目前工作比较领先:(1)美国加州大学生物系计算神经生物学实验室,(2)日本Riken脑科学研究所脑信息研究室,(3)芬兰赫尔辛基工业大学计算机及信息科学实验室,目前发表有关文献较多的刊物有IEEETrans的SP和NN以及NeuralComputation等。本文目的是对ICA的原理、算法及应用作一简述,以引起国内同行对它的关注。将侧重于概念说明,而不追求数学上的严谨性。

2原理

2.1问题的提法,s-(n)是一组互相独立的信源,A是混合矩阵,x-(n)是观察记录,即x-(n)=As-(n)。问题的任务是:在A阵未知且对s-(n)除独立性外无其它先验知识的情况下,求解混矩阵B,使得处理结果y-(n)=Bx-(n)中各分量尽可能互相独立,且逼近s(n)。容易理解,解答不是唯一的,它至少受以下条件的限制:(1)比例不定性:s-(n)中某一分量大K倍时,只要使相应的A阵系数减小K倍,x-(n)便保持不变。

因此,求解时往往把s-(n)假设成具有单位协方差阵,即s-中各分量均值为零,方差为1,且互相独立。(2)排序不定性:y-与s-中各分量排序可以不同。因为只要对调B阵中任意两行,y-中相应元素的位置也便对调。(3)s-(n)中至多只能有一个高斯型信源:这是因为高斯信源的线性组合仍是高斯型的,因此混合后便无法再区别。(4)信源数目N只能小于或等于观测通道数M。N>M情况目前尚未解决。以下讨论设M=N。因此,y-(n)只是在上述条件下对s-(n)的逼近。换名话说,任务的实质是优化问题,它包括两个主要方面:优化判据(目标函数)和寻优算法。

2.2目标函数

这一领域的研究者已经从不同角度提出了多种判据。其中以互信息极小判据(MinimizationofMutualInformation,简记MMI)和信息或熵极大判据(Informax或MaximizationofEntropy,简记ME)应用最广。由于最基本的独立性判据应由概率密度函数(probabilitydensityfunction,简记pdf)引出,而工作时pdf一般是未知的,估计它又比较困难,因此通常采用一些途径绕过这一困难。

常用的方法有两类:①把pdf作级数展开,从而把对pdf的估计转化为对高阶统计量的估计;②在图1的输出端引入非线性环节来建立优化判据。后一作法实际上隐含地引入了高阶统计量。(1)互信息极小判据:统计独立性的最基本判据如下:令p(y-)是y-的联合概率密度函数,pi(yi)是y-中各分量的边际概率密度函数。当且仅当y-中各分量独立时有:p(y-)=∏Ni=1pi(yi)因此用p(y-)与∏i=1pi(yi)间的Kullback-Leibler散度作为独立程度的定量度量:I(y-)=KL[p(y-),∏Ni=1pi(yi)]=∫p(y-)log[p(y-)∏Ni=1pi(yi)]dy-(1)显然,I(y-)0,当且仅当各分量独立时I(y-)=0。因此,互信息极小判据的直接形式是:在y-=Bx-条件下寻找B,使(1)式的I(y-)极小为了使判据实际可用,需要把I(y-)中有关的pdf展成级数。

由于在协方差相等的概率分布中高斯分布的熵值最大,因此展开时常用同协方差的高斯分布作为参考标准。例如,采用Gram-Charlier展开时有:P(yi)PG(yi)=1+13!k2yih3(y-i)+14!k4yih4(yi)+…式中PG(yi)是与P(yi)具有同样方差(σ2=1)和均值(μ=0)的高斯分布。k3yi、k4yi是yi的三、四阶累计量(cumulant),hn(yi)是n阶Hermit多项式。此外还有许多其他展开办法,如Edgeworth展开,利用负熵(Negentropy)等。不论采用何种展开方式,经推导后总可把式(1)近似改成k3、k4的函数:I(y)=F(k3y-,k4y-,B)(1)’F(·)的具体形式多种多样,视推导时的假设而异。

这样就得到互信息判据的实用近似形式:在y-=Bx-条件下寻找B,使式(1)的I(y-)极小(2)Infomax判据:这一判据的特点是在输出端逐分量地引入一个合适的非线性环节把yi转成ri(如图2)。可以证明,如果gi(·)取为对应信源的累积分布函数cdf(它也就是概率密度函数的积分),则使r-=(r1…rN)T的熵极大等效于使I(y-)极小,因此也可达使y-中各分量独立的要求。从而得到Infomax判据:在选定适当gi(·)后,寻找B使熵H(r-)极大需要指出的是,虽然理论上gi(·)应取为各信源的cdf,但实践证明此要求并不很严格,有些取值在0～1之间的单调升函数也可以被采用,如sigmoid函数、tanh(·)等。估计H(r-)固然也涉及pdf,但由于其作用已通过gi(·)引入,所以可以不必再作级数展开而直接用自适应选代寻优步骤求解。文献中还提出了一些其他判据,如极大似然、非线性PCA等,但它们本质上都可统一在信息论的框架下,所以不再一一列举[1]。

3处理算法优化算法

可大致分为两类,即批处理与自适应处理。

3.1批处理批处理比较成熟的方法有两类。较早提出的是成对旋转法[2],其特点是把优化过程分解成两步。先把x-(n)经W阵加以“球化”得z-(n),使z-(n)T=IN,即:各分量不相关且方差为1,然后再寻找合适的正交归一阵U达到使y-各分量独立的目的。前一步类似于PCA,后一步则可利用Givens旋转,根据目标函数,将z-中各分量两两成对反复旋转直到收敛。这种方法计算量较大。1999年,Gadoso提出几种方法对它作了进一步改进[3],其中包括:Maxkurt法、JADE法、SHIBBS法等,限于篇幅,本文不再叙述。近年来,提出的另一类方法是所谓“固定点”法(FixedPointMethod)[4,5],其思路虽来源于自适应处理,但最终算法属于批处理。

简单地说,通过随机梯度法调节B阵来达到优化目标时,有:B(k+1)=B(k)+ΔB(k)ΔB(k)=-μεkB(k)式中k是选代序号,εk是瞬时目标函数。当到达稳态时必有[E是总集均值算子]:E[ΔB(k)]=0(2)如果ΔB(k)与B(k)有关,就可由(2)式解出B的稳态值。不过由于(2)式总是非线性方程,因此求解时仍需要采用数值方法(如牛顿法、共轭梯度法等)迭代求解。实践证明,不论是收敛速度还是计算量,此法均优于前一种方法,而且它还可以根据需要逐次提取最关心的yi,因此是一类值得注意的方法。

3.2结合神经网络的自适应处理结合神经网络的自适应处理算法的框图。1994年Cichocki提出的调节算法是:B(k+1)=B(k)+ΔB(k)ΔB(k)=μk[I-Ψ(y-k)ΦT(y-k)]B(k)式中Ψ、Φ都是N维矢量,其各元素都是单调升的非线性函数:Ψ(yk)=sgnyk·y2k,ΦTy-k=3tanh(10yk)所得结果虽令人鼓舞,但是方法是经验性的。其后学者们从理论上沿着这一方向作了更深入的讨论,并发展出多种算法。概括地说,主要发展有以下几点:

(1)引入自然梯度(或相对梯度)。按照最陡下降的随机梯度法推导出的系数调节公式往往具有如下一般形式:ΔB(k)=μk[B-T(k)-Ψ(y-k)x-Tk]式中的Ψ(y-k)视具体算法而异。Infomax法中Ψ(·)由所选用的g(·)决定;MMI法中则与yk的三、四阶矩有关。B-T(k)是矩阵求逆再转置,它的计算量很大。Amari[7]在1998年提出将最陡下降梯度改为“自然梯度”,两者间关系是:[自然梯度]=[最陡下降梯度]·BT(k)B(k)于是有:ΔB(k)=μk[B-T(k)-Ψ(y-k)x-Tk]BT(k)B(k)=μk[I-Ψ(y-k)y-Tk]B(k)由于此式避免了矩阵求逆,因此计算量明显降低且收敛加快。目前,这一作法已被普遍接受。

(2)引入自然梯度后,采用不同的优化判据得出的调节公式虽各有千秋,但大致都可表示为如下的“串行更新”形式:B(k+1)=B(k)+ΔB(k)=[I+H(y-k)]B(k)只是H(y-k)的具体形式各不相同。串行矩阵更新的算法还具有一些理论上值得注意的性质,如均匀特性(uniformproperty)和等变性(equivariant)等[8,9]。

(3)四阶累计量k4>0的超高斯信号和k4<0的欠高斯信号,其处理过程应当予以区别。采用同一算法效果往往不好。目前的办法多是在调节公式中引入一个开关。根据估计得k4的符号来切换不同算法,如扩展的Infomax法就是一例[10]。此法的系数调节公式是:ΔB(k)=μk[I-Ktanh(y-k)·y-Tk-y-ky-Tk]B(k)其中K是对角阵,其对角元素之值为+1或-1,视该信号分量k4>0或<0而定。为了实时应用,估计K4也可采用递归算法。总之,自适应算法是目前采用较广的方法。

4应用举例

4.1仿真计算为检验经ICA算法分解信源的能力,左图是一组源信号,它们对系统来说是未知的。这一组信号经混合后的观察信号作为(中图所示)ICA算法的输入,分解后的结果如右图所示。可以看到,除了波形的次序、极性和波幅发生变化之外,源信号的波形被很好地分解出来。一般情况下,临床脑电信号中既有超高斯成分(如诱发电位),也有亚高斯成分(如肌电和工频干扰)。为了检验扩展Infomax算法处理这类情况的能力,我们又用此法进行了如图6所示仿真实验。左图第一行是一段自发脑电信号,第二行是仿真的视觉诱发电位,第三行是肌电干扰。混合后的信号(图中第二列所示)经ICA分解得到如右图所示的结果。这一结果表明扩展ICA算法在同时存在超高斯和亚高斯信号的情况下,仍然能够很好地实现盲分解。但应指出:这一仿真结果并不说明通过ICA分解就能直接得到视觉诱发电位,因为还没有涉及头皮上的多导数据。

4.2实验VEP分析(1)多导脑电观察中VEP的增强:需要强调,把多导脑电作ICA分解后直接取出其中与VEP有关的成分,得到的并不是头皮电极处的VEP分量,因为它们只是分解出来的信源,而这些信源的位置并不在头皮上,为了得到电极处测量值中的VEP成分,需按下述步骤处理:用训练得的W阵直接对头皮上取得的多导脑电数据进行ICA分解,得到各独立分量组成的矩耻y=Bx(见图7a);再根据各分量的波形特征及产生时段,选择与VEP有关的一部分分量(例如在前300ms中具有较大幅度的分量),并将其余分量置0,得到新的独立分量矩阵y’;再反变换回头皮各电极处得x’=B-1-y’。这样才能得到去除噪声和干扰后各电极处的VEP。

采用这样的方法可显著地减少提取VEP所需要的累加次数。左图是经3次累加所得VEP,中图是经50次累加所得结果,右图则是用左图经图7中ICA处理后提取的VEP。比较中、右两图,两者波形趋势基本相同,但后者比前者其主要峰、谷显然更清楚,而累加次数由50减到3。(2)ICA分量的空间模式:把某一个ICA分量的瞬时值经B-1逆推回头皮各电极处得x-’后,就可以按断层图的插补方法得到该时该分量在头皮上的空间分布模式。这个空间分布模式也可以用更简单办法得到:只要把逆矩阵B-1中相应于某ICA分量的列中各元素的值赋与头皮各电极处,再作断层图插值,就可以表现该ICA分量在任意时刻的空间分布模式。也就是:x’i(t)=b’ijy’j(t),i=1～N式中b’ij是B-1的第i行第j列元素。

可见ICA分量y’j(t)在头皮各电极处的对应值等于用逆阵B-1第j列各元素来对y’j(t)加权。因此,列矢量b’j=[b’1,…,b’Nj]可以用来统一地表现任意时刻y’j的空间模式。

5总结与展望

本文粗略介绍了ICA的原理、算法和应用,可以看到ICA确是一个值得注意的研究方向,但其理论体系尚未完整,实际采用的处理方法多少还带有经验性。例如为什么对非线性特性gi的要求不甚严格就没有明确解释;又如算法的稳定性、收敛性在实践中是经常遇到的问题。从应用方面看也还有许多待开发的领域,例如如何应用于生理信号的模式识别与系统建模等。从生物医学信号分析的角度看,还有一些亟待深入的问题。例如:

(1)在以上分析中混合阵A被假设为恒定。这对静态的图像分析或固定信源是合理的;但在生理实际中,等效信源一般在空间并不固定,因而混合阵A应视为时变的,而且传导过程中还会引入容积导体的卷积及迟作用。这可能是实际生理信号分解结果不够理想的原因之一。

(2)一般公认,生理信号的非平稳性较强,而以上分析并没有考虑信号的非平稳性。

(3)生理研究表明脑内电活动往往是分区聚集的,活动区的部位随感觉、认知具体任务而定;而各区间的电活动经皮层下的神经网络存在着同步联系。因此各活动区的等效源的电活动不是完全独立的。采用ICA技术如何反映这种情况也值得研究。

(4)采用二阶以上的累计量为判据时,对高斯噪声是透明的;但对非高斯噪声情况如何,尚待研究。应当强调任何方法都不是万能的。合理的途径是把它结合到一个多方法、分层次的信号处理框架中发挥自己的优势。