基于小波-雙譜分析的視覺誘發(fā)腦電特征提取?

嚴(yán) 娜,喬曉艷,李 鵬

(山西大學(xué)電子信息技術(shù)系,山西太原 030006)

0 引 言

視覺誘發(fā)腦電是指對視覺系統(tǒng)施加適當(dāng)?shù)拇碳r所引起的大腦電位變化.其幅值約為0.5 μ V ～50 μ V,頻率范圍0.5 Hz～30 Hz,信噪比為0～10 dB,屬于非平穩(wěn)和非高斯的微弱信號.誘發(fā)腦電反映了大腦在決策和判定過程中的認(rèn)知功能,成為臨床和實驗室腦功能測試的通常選擇,也被大量應(yīng)用于腦機(jī)接口(BCI)系統(tǒng)中作為反映大腦活動的特征信號.近年來,基于P300的BCI系統(tǒng)仍在不斷研究中[1,2],2008年瑞士聯(lián)邦洛桑理工大學(xué)BCI小組使用圖片輪換的方案,要求受試者默數(shù)目標(biāo)刺激圖片的數(shù)目獲取視覺誘發(fā)特征信號,并成功構(gòu)建腦機(jī)接口系統(tǒng)[1].同時,視覺誘發(fā)特征信號在測謊技術(shù)[3]、犯罪意識檢測、病態(tài)人格測試中均得到廣泛應(yīng)用[4].視覺誘發(fā)腦電的特征提取方法主要包括相干平均法、功率譜分析、AAR模型法、獨立分量分析、支持向量機(jī)、高階譜分析等.疊加平均技術(shù)能有效提取誘發(fā)腦電特征,但需要多次刺激疊加,容易引起受試者神經(jīng)系統(tǒng)疲勞,而且該方法忽略了每次刺激之間誘發(fā)電位的變異,不能反應(yīng)誘發(fā)腦電的逐次變化,從而影響特征提取準(zhǔn)確性;功率譜分析法利用了信號頻域的能量信息,但損失了時域信息;AAR模型法更適合分析平穩(wěn)隨機(jī)信號,而EEG信號是高度非平穩(wěn)的信號.此外,AAR模型法對偽跡很敏感;獨立分量分析能從單次(或少次)刺激中提取出誘發(fā)腦電,但各個獨立分量所蘊含的物理意義有待進(jìn)一步研究;支持向量機(jī)只能提取特征,不能提取信號,因而丟失了部分信息;小波變換具有多尺度、相對帶寬恒定特點,適當(dāng)選擇基本小波可在時域和頻域都具有表征信號局部特性的能力[5],但是小波分析不能有效獲取信號的相位信息,也不能抑制加性高斯噪聲.高階累積量定義的高階譜分析可以從更高階的概率結(jié)構(gòu)來表征信號,完全抑制了信號中的高斯噪聲,彌補二階統(tǒng)計量的缺陷[6].但高階譜分析的數(shù)據(jù)來自信號時域信息,丟失了有用的頻域信息,無法進(jìn)行時頻域分析.為了更有效、更精確并能從少次刺激中獲得視覺誘發(fā)腦電特征,本文采用基于小波-雙譜分析相結(jié)合的方法提取誘發(fā)腦電.

1 算法基本理論

1.1 相干平均

相干平均是腦電信號研究中提取腦誘發(fā)電位的傳統(tǒng)方法,其利用對齊疊加平均的原理,可獲得倍信噪比增益[7,8].

設(shè)各次測量信號xi(t)、誘發(fā)電位si(t)和自發(fā)腦電(可視為噪聲)ni(t)的關(guān)系表示如下

相干平均之后信噪比為

式中:σ2表示各次測量信號方差,σ2/N表示N次疊加平均誘發(fā)響應(yīng)方差.可見,信噪比提高了倍.

1.2 小波變換

小波變換是一種時頻分析方法,適用于非平穩(wěn)信號處理,其變焦距特性,容易將特征間的差異突出表現(xiàn).小波分析描述非平穩(wěn)信號是在時間-尺度平面上,將信號分解成一系列小波函數(shù)的疊加,小波窗口大小隨頻率改變.在低頻段,時間分辨率較低而頻率分辨率較高;在高頻段,時間分辨率較高而頻率分辨率較低.由于小波變換的多分辨率特點,適合提取非平穩(wěn)的腦電信號特征.

在多分辨分析中,尺度函數(shù)為φ(t),對應(yīng)的小波函數(shù)為 ψ(t),它們滿足二尺度差分方程

式中:h0(k)對應(yīng)正交低通濾波器系數(shù);h1(k)對應(yīng)正交高通濾波器系數(shù).

由 φj,k和 ψj,k各自的正交性,h0(k),h1(k)可由下式求得

基于小波變換的分析技術(shù)主要包括小波分解和小波重構(gòu)兩部分.設(shè)原始信號為f(n),如果aj(k),dj(k)是多分辨率分析中的離散逼近系數(shù)和離散細(xì)節(jié)系數(shù),h0(k),h1(k)是滿足尺度方程的兩個正交濾波器,則小波分解算法可表示為

令j由0逐漸增大,可得到多分辨率的逐級實現(xiàn)[9].

小波重構(gòu)算法是其分解算法的逆過程,也可由濾波器組實現(xiàn),即

由于小波函數(shù)具有不唯一性,選擇不同的基本小波會產(chǎn)生不同的分析結(jié)果.Daubechies系列小波是工程上應(yīng)用最廣泛、最成熟的緊支集正交實小波函數(shù)族,簡稱dbN小波系(其中N為小波序號這一系列),小波共有49個(db1～db49);其特點是:支集長度L=2N,消失矩階數(shù)p=N;隨著序號N的增大,時間局部性變差,但頻域局部性變好.其它類似的緊支集正交小波,例如coif小波,Sym小波特性均不如db小波,因此,通過比較分析,選擇db7小波作為腦電信號分析的小波函數(shù).

1.3 雙譜分析

雙譜即三階累積量譜,是具有幅值和相位的復(fù)值譜.在高階譜分析中,雙譜具有高階譜的所有特性,并且它的階次最低,計算量也是最小的,雙譜可以較好地反映信號的特征信息.

設(shè)x(n)為零均值、三階實平穩(wěn)隨機(jī)序列,其三階相關(guān)函數(shù)為[10]

其雙譜表達(dá)式為

采用非參數(shù)直接雙譜估計,先計算觀測序列的FFT,再求頻域相關(guān),算法具體步驟描述如下[7]:

設(shè)x(0),x(1),…,x(N-1)為一組測試數(shù)據(jù),fs為采樣頻率.在雙譜域內(nèi),若 ω1和 ω2軸的頻率采樣點數(shù)為N0,則頻率抽樣間隔為Δ0=fs/N0.

1)將所給數(shù)據(jù)分成K段,每段M個觀測樣本,且每段數(shù)據(jù)之間允許重疊,對每段數(shù)據(jù)減去該段的均值,使每段成為零均值序列.

2)計算每段的DFT系數(shù)

3)根據(jù)DFT系數(shù),分別求出每段數(shù)據(jù)的雙譜估計,即

式中:λ1,λ2分別對應(yīng) ω1,ω2軸的 DFT變換后的點,L1表示平滑點數(shù)且N0和L1應(yīng)選擇為滿足的值.

4)對各段雙譜估計的結(jié)果進(jìn)行統(tǒng)計平均

2 實驗范式與數(shù)據(jù)描述

采用Donchin提供的視覺誘發(fā)Oddball腦-機(jī)接口實驗范式,該數(shù)據(jù)被BCI2005競賽選為標(biāo)準(zhǔn)數(shù)據(jù).實驗范式構(gòu)成如圖1所示,實驗中,受試者需注視此6×6字符矩陣中某個字符,目標(biāo)是通過視覺刺激誘發(fā)腦電特征,識別受試者正在注視的字符,從而完成拼寫單詞的任務(wù).

對需要拼寫的單詞中每個字母而言,字符矩陣以 2.5 s為周期顯示,其中,前400 ms字符矩陣沒有被點亮,即無顯示,接下來字符矩陣中的行和列以5.7 Hz的頻率隨機(jī)加亮,即每次加亮100 ms,加亮完成后,矩陣變黑75 ms.6×6字符矩陣共需進(jìn)行12次行、列加亮,其中會有兩次包含目標(biāo)字符(字符所在的特定行和特定列),這些包含目標(biāo)字符的刺激引起的誘發(fā)腦電響應(yīng)與不包含目標(biāo)刺激的誘發(fā)腦電響應(yīng)是不同的.每一個字符重復(fù)15次試驗,以保證獲得足夠強(qiáng)度的誘發(fā)電位.數(shù)據(jù)通過0.1 Hz～60 Hz的帶通濾波和 240 Hz頻率采樣后,以Matlab數(shù)據(jù)格式*.mat提供.有兩位受試者A和B參與了實驗,各有一個包含85個字符的訓(xùn)練文件和包含100個字符的測試文件.所有數(shù)據(jù)文件以單精度存儲,每個文件內(nèi)都包含64導(dǎo)聯(lián)的誘發(fā)腦電數(shù)據(jù).

圖1 視覺誘發(fā)腦控拼寫器界面Fig.1 The interface of visual evoked brain control speller

3 視覺誘發(fā)腦電特征提取

由于誘發(fā)腦電信號的微弱性,在腦電采集過程中常常會存在各種干擾和噪聲,同時在誘發(fā)腦電測量過程中,不可避免地會記錄到自發(fā)腦電、眼電、肌電等多種電生理信號的偽跡以及由于基線漂移引起的記錄誤差,通常在提取特征之前,需要對采集的腦電信號進(jìn)行預(yù)處理,以便抑制和消除測量信號中的干擾和噪聲,而又不損失腦電信息.由于BCI CompetionⅢ提供的數(shù)據(jù)是未經(jīng)處理的原始數(shù)據(jù),本文對腦電信號的預(yù)處理包括了導(dǎo)聯(lián)選擇,共平均參考,帶通濾波處理,偽跡去除,疊加平均等幾個部分.首先,對數(shù)據(jù)進(jìn)行共平均參考(CAR),獲得一個相對理想的參考電極,選擇10個最優(yōu)導(dǎo)聯(lián),通過 0.5 Hz～30 Hz的帶通濾波器,使信號更加平滑,按照每次刺激起始點后1s進(jìn)行數(shù)據(jù)分割,并將數(shù)據(jù)重排和降采樣,在以上預(yù)處理的基礎(chǔ)上,對數(shù)據(jù)進(jìn)行15次疊加平均;然后,按照2.2中描述的多分辨率分析原理,選用db7母函數(shù)進(jìn)行6層小波分解,并對第6層的細(xì)節(jié)系數(shù)進(jìn)行重構(gòu)濾波并白化處理,去掉相關(guān)性.最后,采用2.3中描述的雙譜分析方法分別對靶刺激和非靶刺激下重構(gòu)的誘發(fā)腦電細(xì)節(jié)信號進(jìn)行提取特征.

4 結(jié)果與分析

實驗數(shù)據(jù)采用BCI2000腦機(jī)接口系統(tǒng)和64導(dǎo)聯(lián)的腦電放大器記錄得到,本文僅使用Fz,Fc1,Fc2,C3,Cz,C4,Pz,PO7,PO8,FP1共10個導(dǎo)聯(lián),因為在這些導(dǎo)聯(lián)上P300響應(yīng)較為明顯.

首先,對原始標(biāo)準(zhǔn)實驗的腦電數(shù)據(jù)按照上述預(yù)處理方法進(jìn)行處理.因訓(xùn)練集中共有85個字符,本文以Cz導(dǎo)聯(lián)上的第一個字符為例,按照每次刺激起始點后1s內(nèi)的數(shù)據(jù)進(jìn)行處理,圖2是共平均參考后的結(jié)果,圖3顯示的是15次疊加平均的結(jié)果.經(jīng)過以上預(yù)處理之后,信噪比可提高30%,該腦電數(shù)據(jù)用于后續(xù)的特征提取中.

圖2 Cz導(dǎo)聯(lián)上使用共平均參考進(jìn)行重參考處理Fig.2 Re-reference handling used CAR on Cz channel

圖3 Cz導(dǎo)聯(lián)上原始腦電信號15次疊加平均Fig.3 15 times average?stack of?the original EEG on Cz channel

其次,為了便于比較,本文對記錄的視覺誘發(fā)腦電信號進(jìn)行雙譜分析,同樣以Cz導(dǎo)聯(lián)為例,圖4和圖5分別顯士出了對第一個字符的靶刺激與非靶刺激進(jìn)行雙譜分析的結(jié)果.由圖中可知不論腦功能處于何種狀態(tài),腦電信號各分量均發(fā)生平方相位耦合現(xiàn)象,說明腦電是一種非常復(fù)雜的非高斯隨機(jī)過程,因此采用高階譜技術(shù)處理腦電信號,能獲得更多有關(guān)腦功能狀態(tài)的有用信息.

圖4 雙譜處理三維圖(靶刺激)Fig.4 Three-dimensional map of bispectrum processing(target stimulus)

圖5 雙譜處理三維圖(非靶刺激)Fig.5 Three-dimensional map of bispectrum processing(non-target stimulus)

最后,采用本文提出的小波-雙譜分析方法,即以db7小波母函數(shù)進(jìn)行6層小波分解和重構(gòu),并對細(xì)節(jié)系數(shù)進(jìn)行白化處理,最后再進(jìn)行雙譜估計,本文處理的數(shù)據(jù)長度為60,將其數(shù)據(jù)分成6段,重疊率為50%,圖6和圖7分別顯示了對第一個字符的靶刺激和非靶刺激特征提取的結(jié)果.

圖6 小波-雙譜處理三維圖(靶刺激)Fig.6 Three-dimensional map of Wavelet-bispectrum processing(target stimulus)

圖7 小波-雙譜處理三維圖(非靶刺激)Fig.7 Three-dimensional map of wavelet-bispectrum processing(non-target stimuli)

由圖6,圖7可知:腦功能狀態(tài)的變化在雙譜中直接反映在腦電分量發(fā)生平方相位耦合的差異上,在大腦處于注意力高度集中時,雙譜中的三階能量幾乎集中在α波段的9 Hz的頻率分量上,且這一頻率分量的有序性迅速增強(qiáng),非靶刺激的腦電誘發(fā)信號的雙譜在 α波段和θ波段均出現(xiàn)明顯的非線性相位耦合現(xiàn)象,有幾個比較明顯的譜峰,但在 9 Hz分量(μ節(jié)律)上幅值遠(yuǎn)遠(yuǎn)低于靶刺激在該分量的幅值.基于高階累積量的雙譜分析方法,能將腦電中反映大腦狀態(tài)的信息提取出來,較好地保留了注意相關(guān)和非注意相關(guān)時腦功能狀態(tài)的差異性.

表1給出了靶刺激和非靶刺激在 μ節(jié)律處的幅值,從表中可以明顯看出,小波-雙譜分析方法與雙譜分析方法相比,提取的腦電特征更加明顯,靶刺激相對非靶刺激在μ節(jié)律上腦電信號幅值相差近16倍.

表1 靶刺激與非靶刺激在 9 Hz頻率上腦電幅值相位耦合現(xiàn)象,Tab.1 The EEG amplitude of target stimulus and non-target stimulus in the frequency 9 Hz

5 結(jié) 論

采用現(xiàn)代信號處理技術(shù)分析非平穩(wěn)和非線性的腦電信號并快速有效地提取特征進(jìn)而對腦電信號進(jìn)行自動識別與分類是腦-機(jī)接口技術(shù)研究的一個重要環(huán)節(jié).由于通常的腦電雙譜分析僅僅利用了腦電的時域數(shù)據(jù),丟失了頻域中大量有用的信息.因此,本文采用時頻分析中的小波變換方法對視覺誘發(fā)腦電數(shù)據(jù)進(jìn)行小波分解和重構(gòu),再對重構(gòu)的細(xì)節(jié)系數(shù)進(jìn)行白化和雙譜分析,而且應(yīng)用小波分解多分辨率特性,可以更有效地提取非平穩(wěn)的腦電特征,并且仿真試驗驗證了該算法具有很好的泛化能力,可以用于腦機(jī)接口系統(tǒng)中.

[1]Piccione F,Giorgi F,Tonin P,et al.P300-based brain computer interface:Reliability and performance in healthy and paralysed participants[J].Clinical Neurophysiology,2006,117(3):531-537.

[2]Lebedev M A,Nicolelis M A.Brain-machine interfaces:Past,present and future[J].Trends in Neurosciences,2006,29(9):536-546.

[3]Vahid Abootalebi.A new approach for EEG feature extraction in P300-based lie detection[J].Computer Methods and Programs in Biomedicine,2008,10-14.

[4]Lebedev MA,Nicolelis MA.Brain-machine interfaces:past,present and future[J].Trends in Neurosciences,2006,29(9):536-546.

[5]徐寧壽,張建華,曹正才.小波變換在視覺誘發(fā)腦電信號提取中的應(yīng)用[J].北京工業(yè)大學(xué)學(xué)報,2000,26(4):9-14.

Xu Ningchun,Zhang Jianhua,Cao Zhengcai.Application of wavelet transform technology in VEP signal extract[J].Journal of Beijing Polytechnic University,2000,26(4):9-14.(in Chinese)

[6]王群,樂建威,金頌揚,等.腦電信號高階譜分析技術(shù)的研究[J].中國醫(yī)療器械雜志,2009,33(2):79-82.

Wang Qun,Le Jianwei,Jin Songyang,et al.The study of EEG higher oder spectral analysis technology[J].Chinese Journal of Medical Instrumentation,2009,33(2):79-82.(in Chinese)

[7]黃日輝,李霆,阜艷,等.誘發(fā)腦電提取方法的研究進(jìn)展[J].現(xiàn)代電子技術(shù),2008,22(285):139-144.

Huang Rihui,Li Ting,Fu Yan,et al.Review of methods for extracting evoked potential[J].Modern Electronics Technique,2008,22(285):139-144.(in Chinese)

[8]王洪濤.視覺誘發(fā)電位腦機(jī)接口關(guān)鍵技術(shù)研究[J].重慶文理學(xué)院學(xué)報(自然科學(xué)版),2010,29(1):69-74.

Wang Hongtao.Researches of key techniques on brain-computer interface of visual evoked potential[J].Journal of Chongqing University of Arts and Science(Natural Science Edition),2010,29(1):69-74.(in Chinese)

[9]Eduardo Bayro-corrochano.The theory and use of the quaternion wavelet transform[J].Journal of Mathematical Imaging and Vision,2006,24:19-35.

[10]Huang Liyu,Wang Weirong,Singare S.Bispectrum quantification analysis of EEG and artificial neural network may classify ischemic states[J].Lecture Notes in Computer Science,2006,4233:533-542.