王永軒 邱天爽 劉 蓉
(大連理工大學(xué) 電子信息與電氣工程學(xué)部,大連 116024)
誘發(fā)電位(evoked potential,EP)反映了中樞神經(jīng)系統(tǒng)的感覺通路及皮層區(qū)域的神經(jīng)電活動情況,是臨床醫(yī)學(xué)診斷神經(jīng)系統(tǒng)損傷及病變的重要手段之一,也是實現(xiàn)腦-機(jī)接口(brain-computer interface,BCI)的重要腦電信號之一[1-2]。由頭皮表面測得的EP信號是和自發(fā)腦電(EEG)等背景噪聲同時記錄到的,信噪頻譜重疊,而且信噪比很低,所以提取難度較大。目前,臨床醫(yī)學(xué)通常使用的方法是疊加平均法[3-4],通過很多的刺激次數(shù)來得到足夠的信噪比。然而,疊加平均法中太多的刺激次數(shù)會增大每次EP信號之間的差異,而且反復(fù)刺激容易引起神經(jīng)疲勞甚至產(chǎn)生暴發(fā)偽跡,因此少次提取問題成為一個重要的研究熱點和難點。此外,EP信號是由特定位置的大腦皮層細(xì)胞產(chǎn)生的反應(yīng),傳導(dǎo)到頭皮不同位置時相當(dāng)于經(jīng)過了不同的系統(tǒng),很難用線性模型描述,從而 EP信號的單導(dǎo)提取也顯得尤為重要。
獨立分量分析(independent component analysis,ICA)是為解決盲源分離問題而產(chǎn)生的一種有效算法[5-8]。它從多路觀測信號出發(fā),對已知信息量很少的源信號進(jìn)行估計,最后得到相互獨立的源信號的近似值。ICA算法處理的對象是相互獨立的源信號經(jīng)未知線性組合而產(chǎn)生的一組混合信號,在實際使用時,一般假設(shè)觀測信號路數(shù)大于或等于獨立源的個數(shù)。顯然在EP信號單導(dǎo)提取任務(wù)中無法利用多導(dǎo)信號提供的空間信息,因此不能直接使用ICA算法。為了實現(xiàn)這一目的,文獻(xiàn)[9-10]均引入了一路噪聲(EEG)作為虛擬通道,將單路觀測信號擴(kuò)展為兩路,以滿足ICA算法的應(yīng)用條件。這個模型建立的前提必須是虛擬噪聲通道與觀測信號中的噪聲具有很高的相關(guān)性,但實驗表明EEG作為一種隨機(jī)信號,不同時間段的EEG之間的相關(guān)性很弱,即虛擬噪聲通道與觀測信號中的噪聲匹配度很低,因此這個模型很難消除EP信號中的噪聲。文獻(xiàn)[11,12]假設(shè)單路觀測信號是由 EP信號和多個噪聲分量疊加而成,引入的多路虛擬噪聲通道分別代表了這些噪聲分量。雖然模型中假設(shè)多路虛擬噪聲通道是獨立源信號,但是ICA算法要求輸出信號最大程度符合獨立性準(zhǔn)則,所以引入的虛擬通道之間的實際相關(guān)程度會束縛ICA算法的效果。針對以上問題,本文在EP信號的單導(dǎo)少次提取中,將少次信噪混合信號作為虛擬觀測通道,構(gòu)建了新的虛擬通道模型,增加了ICA算法中解混矩陣的自由度,從而提升了ICA的應(yīng)用性能。
ICA算法處理的信號是多路觀測信號,將多路觀測信號用向量表示為
式中,M表示觀測信號路數(shù),x的各個分量xm表示的是第m路觀測信號,即
式中,N為采樣點數(shù)。在ICA算法中把觀測信號以及源信號都看作是隨機(jī)變量,每個采樣點是隨機(jī)變量的不同取值,在算法中的各種統(tǒng)計量正是利用這些采樣點來進(jìn)行計算。多路觀測信號的形成是多路未知源信號的未知線性混合,數(shù)學(xué)模型為
式中,s= [s1,s2,…,sM]T為未知源信號向量,A 為M×M維未知混合矩陣。ICA算法的應(yīng)用前提是s的各個分量sm相互獨立,且A是滿秩的。在這個假設(shè)條件下根據(jù)觀測信號x得到源信號的估計^s,因此^s是在獨立性準(zhǔn)則條件下的最佳估計。通過尋找一個解混矩陣BM×M,得到多路輸出信號
如果y的各個分量ym在最大程度符上合獨立性準(zhǔn)則,那么就可以認(rèn)為y是 s的一個最佳估計。顯然如果解混矩陣B是混合矩陣A的逆矩陣,即
那么ICA算法將準(zhǔn)確還原多路源信號。由于ICA算法假設(shè)源信號是統(tǒng)計獨立的,因此需要給出獨立性準(zhǔn)則以及反映獨立性準(zhǔn)則的目標(biāo)函數(shù),當(dāng)目標(biāo)函數(shù)達(dá)到極值時,就可以認(rèn)為輸出信號是獨立的,也就是源信號的估計。比較經(jīng)典的算法包括FastICA算法和 Infomax算法等[13]。
單次EP信號觀測模型可表示為
式中,x1= [x1(1),x1(2),…,x1(N)]為單次含噪信號,s1為純凈EP信號,v1表示 EEG等噪聲。為了應(yīng)用ICA算法實現(xiàn)單導(dǎo)誘發(fā)電位提取,需要引入虛擬通道,將一路觀測信號擴(kuò)展為多路。文獻(xiàn)[11-12]的方法是將刺激前記錄的EEG按刺激間隔分段截取作為虛擬噪聲通道,然后與單次含噪信號一起構(gòu)成多路觀測信號,從而建立如下模型
式中,A 為未知混合矩陣,s= [s1,v2,…,vM]T為獨立源向量。雖然模型中假設(shè)多路虛擬噪聲通道是獨立源信號,但是ICA算法要求輸出信號最大程度符合獨立性準(zhǔn)則,所以引入的虛擬通道之間的實際相關(guān)程度會束縛ICA算法的效果。事實上由于未知混合矩陣的元素給定值很多,所以解混矩陣的自由度會很低,這樣ICA輸出分量的獨立性也會受到限制。
在本算法模型中,引入的是虛擬觀測通道,也就是將刺激后記錄的單導(dǎo)少次含噪信號形成多路觀測信號,同時假設(shè)多路觀測信號中的噪聲均由有限個獨立噪聲模式組成,即
式中,uj,(j=2,3,…,M) 為獨立噪聲模式,aij為第j種噪聲模式在第i路噪聲中的權(quán)重。另外根據(jù)EP信號的鎖時性[1],可假設(shè)單導(dǎo)相鄰少次記錄中 EP信號基本不變,這樣虛擬通道模型可用矩陣表示為
式中,A 為未知混合矩陣,s= [s1,u2,u3,…,uM]T為獨立源向量。在這個模型中未知混合矩陣假定元素很少,提高了解混矩陣的自由度,從而 ICA輸出分量的獨立程度得到了提高。因為模型中的每路觀測信號之間(即每次含噪信號之間)的相關(guān)性不強(qiáng),所以觀測信號矩陣可以保證行滿秩,也就是混合矩陣一定是滿秩的。
仿真實驗選取4路獨立源信號,分別代表 EP信號(通過疊加平均法得到的貓的體感誘發(fā)電位)、工頻干擾(50 Hz正弦波)、自發(fā)腦電(高斯分布隨機(jī)信號)和亞高斯噪聲(周期方波)。4路源信號經(jīng)過隨機(jī)矩陣 A4×4混合得到4路觀測信號,最后通過fastICA算法提取得到4路輸出信號,源信號、觀測信號和提取信號如圖1所示。其中觀測信號的混合矩陣 A4×4由 Matlab隨機(jī)生成:
從圖1中可以看到,源信號淹沒在觀測信號中,而提取后得到的信號與源信號非常相似,同時可以發(fā)現(xiàn)提取后得到的信號的順序與源信號相比發(fā)生了變化,這是由 ICA的不確定性造成的??傮w說來,對于符合ICA假設(shè)條件的源信號與混合方式得到的結(jié)果是比較理想的,即使在較低信噪比情況下提取結(jié)果與原始信號仍有很高的相關(guān)系數(shù)。圖1中四路信號的相關(guān)系數(shù)分別為0.9994(EP)、-0.9999(正弦)、-0.9899(高斯)和0.9999(方波)。
在虛擬通道ICA仿真實驗中分別對兩種模型進(jìn)行了仿真。第1種是文獻(xiàn)[11,12]的模型,其中第1路觀測信號是EP信號與噪聲的混合信號,信噪比為0 dB,另外3路為其他時間段記錄的噪聲,觀測信號與提取信號如圖2所示。第2種為所給出的模型,4路信號均為含噪信號,其中的EP信號是相同的,但噪聲為不同時間段的記錄,信噪比均為0 dB,觀測信號與提取信號如圖3所示。
從圖2中可以看出4路輸出信號中第一路為EP信號,但信噪比不但沒有上升反而略有下降,這是因為4路觀測信號中的噪聲之間的相關(guān)程度很弱,彼此無法互相削弱抵消。其他3路輸出信號與觀測信號完全相同,這是由模型中混合矩陣的約束造成的,解混矩陣與混合矩陣有著完全相同的型式。而從圖3中可以看到EP信號出現(xiàn)在輸出信號中的第1路,信噪比得到了明顯的提高。
圖1 4路源信號、觀測信號和提取信號。(a)EP信號;(b)正弦波;(c)高斯噪聲;(d)周期方波;(e)~(h)相應(yīng)的4路觀測信號;(i)~(l)相應(yīng)的4路提取信號Fig.1 Four channels of source signals,observed signals and extracted signals.(a)EP signal;(b)Sine wave;(c)Gaussian noise;(d)Cycle square wave;(e)~(h)Four channels of observed signals;(i)~(l)Four channels of extracted signals
圖2 第1種虛擬通道模型的觀測信號與提取信號。(a)~(d)4路觀測信號;(e)~(h)相應(yīng)的4路提取信號Fig.2 The observed and extracted signals with the first virtual-channel ICA model.(a) ~ (d)Four channels of observed signals;(e)~(h)Four channels of extracted signals
實驗中使用的視覺誘發(fā)電位(visual evoked potential,VEP)的原始數(shù)據(jù)由美國 Neuroscan公司32導(dǎo)腦電設(shè)備Scan4.3采集得到。對12位受試者進(jìn)行模式翻轉(zhuǎn)視覺誘發(fā)電位(pattern reversal visual evoked potential,PRVEP)測試。在測試過程中進(jìn)行多次刺激,經(jīng)過機(jī)器自動識別與人工篩選去除了有明顯偽跡的記錄。使用傳統(tǒng)的疊加平均法得到每導(dǎo)記錄中的 VEP,圖4為受試者 S1的 Cz,O1,Oz,O2四路導(dǎo)聯(lián)的VEP波形。
從圖4中可以看到 O1,Oz,O2等3個導(dǎo)聯(lián)均在刺激后的100 ms左右出現(xiàn)了P100峰值,且每導(dǎo)記錄中的VEP波形非常接近,而Cz導(dǎo)聯(lián)沒有出現(xiàn)明顯的P100峰值,這反映了VEP主要出現(xiàn)在頭皮枕骨處這一現(xiàn)象[1]。
選擇受試者 S1的 Cz,O1,Oz,O2等4個導(dǎo)聯(lián)單次記錄作為觀測數(shù)據(jù)使用ICA算法進(jìn)行處理,觀測信號與提取信號如圖5所示,由圖5可見在4路提取信號中任何一路都沒有明顯出現(xiàn)期望的P100峰值。
圖3 第2種虛擬通道模型的觀測信號與提取信號。(a)~(d)4路觀測信號;(e)~(h)相應(yīng)的4路提取信號Fig.3 The observed and extracted signals with the second virtual-channel ICA model.(a) ~ (d)Four channels of observed signals;(e)~(h)Four channels of extracted signals
出現(xiàn)這種結(jié)果的原因在于,雖然多導(dǎo)腦電信號在形式上符合ICA的模型,但實際上并不真正符合。以兩路觀測信號為例進(jìn)行分析,假設(shè)兩路信號分別為xi=si+vi,i=1,2。對于相鄰導(dǎo)聯(lián),其中的EP信號非常接近(見圖4),因此觀測信號可以簡化為xi=s+vi(先不考慮EEG等噪聲的相近程度)。
ICA算法的任何一路輸出實際上就是多路輸入的線性組合,即y=k1x1+k2x2,通過ICA算法可以得到在獨立性準(zhǔn)則條件下的加權(quán)系數(shù)k1和k2。然而無論k1和k2是通過什么準(zhǔn)則確定的,y和s的最大相關(guān)系數(shù)是有限的。計算 y和s的相關(guān)系數(shù) ry,然后通過對k1和k2求偏導(dǎo)可得最大相關(guān)系數(shù)為
式中,rv為v1和v2的相關(guān)系數(shù),λ為信噪比。顯然rv越小越好,而實驗表明相鄰導(dǎo)聯(lián)中的EEG等噪聲具有較大的正相關(guān)性(圖5中 O1,Oz,O2這3個導(dǎo)聯(lián)的波形比較接近),所以無論加權(quán)系數(shù)k1和k2如何選擇,也就是ICA中的解混矩陣如何選擇,都無法得到很好的結(jié)果。
圖4 疊加平均法得到的VEP波形。(a)Cz導(dǎo)聯(lián);(b)O1導(dǎo)聯(lián);(c)Oz導(dǎo)聯(lián);(d)O2導(dǎo)聯(lián)Fig.4 The VEP waveforms with ensemble average method.(a)The Cz channel;(b)The O1 channel;(c)The Oz channel;(d)The O2 channel
使用受試者S1的Oz導(dǎo)聯(lián)的連續(xù)4次記錄構(gòu)建虛擬通道ICA模型,觀測信號與提取信號如圖6所示。圖6中第3路提取信號中出現(xiàn)了P100峰值的輪廓,雖然與圖4所示的疊加平均結(jié)果相比仍然有很大噪聲,但原來完全淹沒在噪聲里的EP信號得到了提取,背景噪聲得到了一定程度的去除。如果以圖4所示的疊加平均結(jié)果作為純凈EP信號,經(jīng)計算可知提取信號與觀測信號相比信噪比提高了約為12 dB。此外實驗也證明了隨著虛擬通道數(shù)量(即刺激次數(shù))的增加,EP信號提取效果也會更好。圖7(b)為使用10路虛擬通道(連續(xù)10次記錄)經(jīng)ICA算法得到的結(jié)果,此時信噪比提高約為20 dB。
圖5 使用Cz,O1,Oz,O2四導(dǎo)數(shù)據(jù)進(jìn)行多導(dǎo)聯(lián) ICA的VEP提取。(a)Cz導(dǎo)聯(lián);(b)O1導(dǎo)聯(lián);(c)Oz導(dǎo)聯(lián);(d)O2導(dǎo)聯(lián);(e)~(h)相應(yīng)的4路提取信號Fig.5 VEP extraction with the Cz,O1,Oz and O2 channels based on ICA method.(a)The Cz channel;(b)The O1 channel;(c)The Oz channel;(d)The O2 channel;(e)~(h)Four channels of extracted signals
圖6 使用Oz導(dǎo)聯(lián)連續(xù)4次數(shù)據(jù)進(jìn)行虛擬通道ICA的VEP提取。(a)~(d)觀測信號;(e)~(h)相應(yīng)的提取信號Fig.6 VEP extraction with four times Oz channel observed signals based on virtual-channel ICA method.(a)-(d)The observed signals;(e)-(h)The extracted signals
圖8為12位受試者分別在4路和10路虛擬通道ICA方法下得到的 Cz,O1,Oz,O2等 4個導(dǎo)聯(lián)VEP相關(guān)系數(shù)的統(tǒng)計結(jié)果(均值與方差)。首先使用疊加平均法分別得到每位受試者 Cz,O1,Oz,O2各個導(dǎo)聯(lián)的VEP作為標(biāo)準(zhǔn)信號。然后使用4路和10路虛擬通道ICA方法提取VEP,并計算與相應(yīng)的標(biāo)準(zhǔn)信號的相關(guān)系數(shù)以及每類相關(guān)系數(shù)的均值與方差。通過統(tǒng)計分析可知在顯著性水平α=0.01條件下4路和10路虛擬通道方法的 O1,Oz,O2導(dǎo)的VEP提取差別是顯著的。此外,在實驗中還觀察到Cz導(dǎo)聯(lián)中的VEP很微弱且成倒相現(xiàn)象,反映在統(tǒng)計結(jié)果中Cz導(dǎo)聯(lián)的相關(guān)系數(shù)明顯較小,這說明VEP并不主要出現(xiàn)在頂葉處。
圖7 Oz導(dǎo)聯(lián)VEP波形比較。(a)4路虛擬通道ICA的提取結(jié)果;(b)10路虛擬通道ICA的提取結(jié)果;(c)疊加平均法結(jié)果Fig.7 The extracted VEPs of Cz channel with different methods.(a)Four virtual channels ICA method;(b)Ten virtual channels ICA method;(c)The estimated result with ensemble average method.
圖8 4路和10路虛擬通道 ICA方法下得到的 Cz,O1,Oz,O2四導(dǎo)聯(lián) VEP相關(guān)系數(shù)的統(tǒng)計結(jié)果。顯著性差異由*表示(**P<0.01)Fig.8 The statistical results of estimated VEP correlation coefficients of the Cz,O1,Oz and O2 channels based on four and ten virtual channels ICA methods.Statistically significant differences denoted by asterisks(**P <0.01).
ICA模型的解混矩陣相當(dāng)于一個空間濾波器組,任何一路輸出都是多路觀測信號的加權(quán)組合。若多導(dǎo)觀測信號中的EEG等噪聲具有較強(qiáng)的正相關(guān)性,則任何一路輸出都不可能減弱或消除噪聲,所以很難得到期望的EP信號。如果能將觀測信號在形式上組成ICA算法所需要的觀測矩陣,而噪聲之間具有較弱的相關(guān)性甚至負(fù)相關(guān)性,那么ICA是有可能達(dá)到目的的。這樣空間信息可以來源于廣義的信號空間而并不局限于物理空間(多傳感器得到的多路觀測信號)。因此,本研究針對單導(dǎo)EP信號提取問題,將單導(dǎo)多次觀測信號以虛擬通道的觀點構(gòu)成多路觀測信號,實現(xiàn)了EP信號的有效提取。仿真實驗表明當(dāng)源信號與混合方式符合ICA假設(shè)條件時可以得到很好的結(jié)果,即使在較低信噪比情況下提取結(jié)果與原始信號的相關(guān)系數(shù)仍然很高。但是在多導(dǎo)VEP提取中,并沒得到預(yù)期的結(jié)果。這說明應(yīng)用ICA對多導(dǎo)信號進(jìn)行處理,其前提是多導(dǎo)觀測信號必須符合ICA的假設(shè)條件。在多導(dǎo)腦電信號處理中,因為每導(dǎo)腦電信號都是EP與EEG的疊加,而ICA算法卻要求它們都是來自某些獨立源的某種線性組合,所以直接使用ICA方法處理多導(dǎo)腦電信號以實現(xiàn)信號分離是不適合的。而將單導(dǎo)觀測信號以虛擬通道的觀點構(gòu)成多路觀測信號以后,近似符合了ICA的應(yīng)用條件,可以實現(xiàn)EP信號的有效提取。從空間濾波器角度分析可得出結(jié)論,基于虛擬通道的ICA算法實際上實現(xiàn)了獨立性準(zhǔn)則下的加權(quán)疊加平均效果。
[1]潘映輔.臨床誘發(fā)電位學(xué)(第二版)[M].北京:人民衛(wèi)生出版社,2000.
[2]Wang Yijun,Gao Xiaorong,Hong Bo,et al.Brain-computer interfaces based on visual evoked potentials [J]. IEEE Engineering in Medicine and Biology Magazine,2008,27(5):64-71.
[3]Davila CE,Mobin MS.Weighted averaging of evoked potentials[J].IEEE Trans Biomed Eng,1992,39:338-345.
[4]Bezerianos A,Laskaris N,F(xiàn)otopoulos S,et al.Data dependent weighted averages for recording of evoked potential signals[J].Electroencephalography and Clinical Neurophysiology,1995,96:468-471.
[5]Aapo H,Erkki O.Independent component analysis:algorithms and applications[J].Neural Networks,2000,13(4-5):411-430.
[6]Zou Ling,Zhang Yingchun.Single-trial evoked potentials study by combining wavelet denoising and principal component analysis methods[J].Clin Neurophysiol,2010,27:17-24.
[7]Hu L,Zhang GZ,Hung YS,et al.Single-trial detection of somatosensory evoked potentials by probabilistic independent component analysis and wavelet filtering[J].Clin Neurophysiol,2011,122:1429-1439.
[8]Villa NC, James CJ. Independent component analysis for auditory evoked potentials and cochlear implant artifact estimation[J].IEEE Trans Biomed Eng,2011,58:348-354.
[9]師黎,鐘麗輝,王端.基于虛擬通道ICA-WT大鼠視覺誘發(fā)電位少次提取[J].中國生物醫(yī)學(xué)工程學(xué)報,2010,29(3):379-383.
[10]李洪,孫云蓮.基于EMD虛擬通道的ICA算法在信號消噪中的應(yīng)用[J].中國生物醫(yī)學(xué)工程學(xué)報,2007,30(5):33-36.
[11]季忠,金濤,楊炯明,等.虛擬噪聲通道在基于ICA消噪過程中的應(yīng)用 [J]. 中國機(jī)械工程,2005,16(4):350-364.
[12]劉欽團(tuán),邱飛岳,李浩君.基于虛擬通道的 ICA的 P-VEP提取方法的研究[J].計算機(jī)工程與科學(xué),2010,32(7):137-139.
[13]周宗潭,董國華,徐昕,等.獨立分量分析[M].北京:電子工業(yè)出版社,2007.