獨立分量生物醫學運用探析
前言：本站為你精心整理了獨立分量生物醫學運用探析范文，希望能為你的創作提供參考價值，我們的客服老師可以幫助你提供個性化的參考范文，歡迎咨詢。
獨立分量分析(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)采用二階以上的累計量為判據時,對高斯噪聲是透明的;但對非高斯噪聲情況如何,尚待研究。應當強調任何方法都不是萬能的。合理的途徑是把它結合到一個多方法、分層次的信號處理框架中發揮自己的優勢。
