基于DTCWT的運(yùn)動想象腦電特征提取

湯 偉 耿逸飛

(陜西科技大學(xué)電氣與控制工程學(xué)院 陜西 西安 710021) (陜西科技大學(xué)工業(yè)自動化研究所 陜西 西安 710021)

0 引 言

腦機(jī)接口(Brain-Computer Interface,BCI)是利用大腦所產(chǎn)生的信號來建立人腦與外部設(shè)備的信息控制的交互技術(shù)[1],在諸多方向有著顯著的研究意義,如康復(fù)醫(yī)療[2]、游戲娛樂[3]等。運(yùn)動想象是人對自身進(jìn)行控制的意念活動。當(dāng)人在想象單側(cè)肢體運(yùn)動時(shí),對側(cè)大腦皮層處于激活狀態(tài),相關(guān)的節(jié)律信號幅值會降低;而同側(cè)處于阻滯狀態(tài),相關(guān)的節(jié)律信號幅值會升高。上述現(xiàn)象分別稱為事件相關(guān)去同步(Event-Related Desynchronization,ERD)以及事件相關(guān)同步[4](Event-Related Synchronization,ERS)。因此,在ERS/ERD的基礎(chǔ)上,我們能夠?qū)崿F(xiàn)對左右手運(yùn)動想象任務(wù)的分類。由于大腦系統(tǒng)十分復(fù)雜,所產(chǎn)生的信號隨機(jī)性強(qiáng)且不穩(wěn)定,非線性特征顯著,因此如何有效提取腦電特征成了一個難題。

目前,腦電信號常用的特征提取方法有:單一時(shí)域或頻域方法、時(shí)頻域結(jié)合方法和非線性分析方法等。其中前兩種方法適用于分析平穩(wěn)信號,應(yīng)用于非平穩(wěn)的腦電信號時(shí)效果不佳[4]。由于腦電信號具有顯著的非線性特征,非線性分析方法十分適用于腦電信號,常用的方法有:香農(nóng)熵、自回歸(Auto Regressive,AR)模型參數(shù)、分形維數(shù)和Lempel-Ziv復(fù)雜度等[5]。孫會文等[6]提出的結(jié)合HHT與AR模型參數(shù)的腦電分類方法,分類識別率達(dá)81.08%,適應(yīng)性良好。羅志增等[7]提出的基于多尺度Lempel-Ziv復(fù)雜度的特征提取方法,得到了較高的分類識別率,達(dá)87.87%。上述方法從非線性角度對腦電進(jìn)行識別,都表現(xiàn)出了較高的識別率,但是特征選取考慮的都是單一角度,忽略了腦電信號在時(shí)頻域所包含的有效信息[8]。因此本文提出一種時(shí)頻域特征與非線性特征相融合的提取方法,從多個角度完備地描述腦電信號。

由于原始腦電中夾雜與運(yùn)動想象無關(guān)的信號頻段,本文引入雙樹復(fù)小波(Dual-tree Complex Wavelet Transform,DTCWT)對C3、C4通道的腦電信號進(jìn)行最優(yōu)頻段選取,這樣既去除了無關(guān)腦電頻段,又進(jìn)行頻率分析,同時(shí)提高了分類準(zhǔn)確率。在DTCWT分解基礎(chǔ)上,選取Hilbert瞬時(shí)能量譜、邊際能量譜和Lempel-Ziv復(fù)雜度作為腦電信號特征,最后采用線性判別分析對其實(shí)行分類。仿真實(shí)驗(yàn)表明,本文提出的基于DTCWT的多特征融合方法,相較于使用相同數(shù)據(jù)集的文獻(xiàn)[9-10]方法,具有更大的最大互信息與更高的平均識別率。

1 實(shí)驗(yàn)數(shù)據(jù)與方法

1.1 實(shí)驗(yàn)數(shù)據(jù)

本文使用的數(shù)據(jù)為被測試者在進(jìn)行左右手運(yùn)動想象任務(wù)時(shí)采集的腦電信號,源自第三屆國際腦-機(jī)接口(BCI Competition Ⅲ)的dataset Ⅲb。EEG信號按照國際10-20導(dǎo)聯(lián)系統(tǒng)的C3、C4、Cz通道采集獲得。該數(shù)據(jù)集包含3名受試者的實(shí)驗(yàn)結(jié)果,分別為O3VR、S4和X11。O3VR為400個試次,S4和X11皆為1 080個試次,每次實(shí)驗(yàn)時(shí)間為8 s,采樣頻率為125 Hz。單次實(shí)驗(yàn)過程如圖1所示。

圖1 單次實(shí)驗(yàn)過程

受試者安靜坐在顯示器前,根據(jù)屏幕提示進(jìn)行左右手運(yùn)動想象任務(wù)。0～2 s,屏幕空白,受試者進(jìn)行放松;2～3 s,屏幕投放“+”,準(zhǔn)備運(yùn)動想象任務(wù);3～9 s,屏幕顯示左右方向箭頭,受試者根據(jù)對應(yīng)箭頭進(jìn)行運(yùn)動想象任務(wù)。t=9 s時(shí),實(shí)驗(yàn)結(jié)束。

1.2 實(shí)驗(yàn)原理

1.2.1 雙樹復(fù)小波變換

雙樹復(fù)小波不僅具有離散小波變換的多尺度分析的能力,同時(shí)有效解決了復(fù)小波變換中計(jì)算復(fù)雜的問題,相對于傳統(tǒng)離散小波,可以較為精確地描述信號的細(xì)節(jié)信息,已成功應(yīng)用于軸承故障分析和面部識別等領(lǐng)域[11]。

DTCWT采用實(shí)部樹和虛部樹對信號進(jìn)行分解與重構(gòu),實(shí)部樹和虛部樹采用不同的實(shí)數(shù)濾波器組,從而得到實(shí)部系數(shù)與虛部系數(shù),使得實(shí)部樹與虛部樹的信息互補(bǔ),從而保證了信號分析的完備性。另外,DTCWT在多尺度分解中運(yùn)用了二分法來降低冗余度,有效提高了算法的計(jì)算速度。

設(shè)腦電信號c(n)的采樣頻率為fs,用DTCWT將腦電信號分為L層,則能得到CAL、CDL、CDL-1、…、CD1等L+1子帶的復(fù)系數(shù)。重構(gòu)后得到的信號分量AL、DL、DL-1、…、D1所對應(yīng)的頻帶范圍分別為[0,fs/2L+1]、[fs/2L+1,fs/2L]、[fs/2L,fs/2L-1]、…、[fs/22,fs/2]。因此,DTCWT進(jìn)行L層分解可得L+1個重構(gòu)分量。

1.2.2 時(shí)頻-非線性特征提取

在對腦電信號進(jìn)行特征提取時(shí),先進(jìn)行DTCWT分解,得到一系列子頻帶,ERS/ERD現(xiàn)象主要發(fā)生在α(8～14 Hz)和β(14～30 Hz)節(jié)律,根據(jù)DTCWT的分解規(guī)則,分解尺度定為4層。子帶和節(jié)律信號的關(guān)系如表1所示。

表1 子帶信號和節(jié)律信號的對應(yīng)關(guān)系

考慮到ERS/ERD現(xiàn)象,本實(shí)驗(yàn)只對D3與D2分量進(jìn)行特征值求取,這樣既能精確且完備地提取有效信息,又能降低特征向量維度。分別計(jì)算每個數(shù)據(jù)段的Hilbert瞬時(shí)能量譜(Instantaneous Energy Spectrum,IES)、邊際能量譜(Marginal Energy Spectrum,MES)和Lempel-Ziv復(fù)雜度三個特征參數(shù),從而構(gòu)成時(shí)頻域-非線性特征向量。對于單次實(shí)驗(yàn)數(shù)據(jù),對C3和C4兩個電極需要計(jì)算4個IES參數(shù)、4個MES參數(shù)、4個Lempel-Ziv復(fù)雜度參數(shù),共計(jì)12個參數(shù)值組成12維特征向量。

3個特征參數(shù)的定義如下:

(1) 對于一個給定信號C(t),進(jìn)行Hilbert變換[12]:

(1)

可得到解析信號:

Z(t)=C(t)+jY(t)=A(t)ejφ(t)

(2)

式中:A(t)和φ(t)分別為瞬時(shí)幅值和瞬時(shí)相位。

由A(t)和φ(t)求取瞬時(shí)頻率ω(t):

(3)

則信號幅度在時(shí)-頻域的分布情況,即Hilbert譜為:

(4)

式中:Re為取實(shí)部。

根據(jù)式(4)可求取IES和MES:

(5)

式中:[ω1,ω2]為頻率范圍;[t1,t2]為時(shí)間范圍。

IES表現(xiàn)為時(shí)域能量特征,MES表現(xiàn)為頻域能量特征。

(2) Lempel-Ziv復(fù)雜度描述了序列與隨機(jī)序列的相似度,復(fù)雜度越大,序列新模式增加的速率越快,越接近隨機(jī)序列[13]。設(shè)序列S(s1,s2,…,sn)的Lempel-Ziv復(fù)雜度為z(n),當(dāng)n→∞時(shí),z(n)會趨于定值n/logln,l為粗粒化段數(shù)(傳統(tǒng)二值化時(shí),l=2),則定義歸一化Lempel-Ziv復(fù)雜度為:

(6)

為了驗(yàn)證上述三個特征能夠較好地反映ERS/ERD現(xiàn)象,從流程的第1 s開始,采用滑動窗來提取特征,每次滑動1個采樣點(diǎn),直到最后1 s結(jié)束,從而獲得單次實(shí)驗(yàn)任務(wù)的特征值時(shí)間序列。將同一運(yùn)動想象任務(wù)的所有實(shí)驗(yàn)所得的特征值疊加平均,得到平均特征值時(shí)間序列。圖2(a)為被試者O3VR左手運(yùn)動想象的平均瞬時(shí)能量圖,可以看出,在4～5 s的時(shí)間段內(nèi),被試者進(jìn)行左手想象任務(wù)時(shí),C3通道較C4通道具有更高的能量分布;圖2(b)為被試者O3VR右手運(yùn)動想象的平均瞬時(shí)能量圖,進(jìn)行右手想象任務(wù)時(shí),C4通道較C3通道具有更高的能量分布,符合ERS/ERD原理。圖3為左右手運(yùn)動想象的平均邊際譜幅值圖,圖3(a)為左手,圖3(b)為右手。可以看出在8～15 Hz頻段,兩種不同動作的運(yùn)動想象下,皆是C3通道具有較高能量,不符合ERS/ERD現(xiàn)象,但在15～30 Hz頻段,信號的特征符合ERS/ERD現(xiàn)象。若從時(shí)域能量角度統(tǒng)一分析信號特征,低能高頻信號所含有效信息會被高能低頻信號覆蓋,從而導(dǎo)致時(shí)域能量特征的敏感度下降、被試識別率差的結(jié)果。因此,本文利用DTCWT將信號拆分為成多個子頻帶分量,保證了有效信息的完整保留,提高了被試識別率。

(a)

(b)圖2 C3、C4的平均瞬時(shí)能量

(a)

(b)圖3 C3、C4的平均邊際譜幅值

圖4為左右手運(yùn)動想象的平均Lempel-Ziv復(fù)雜度幅值圖。圖4(a)為想象左手運(yùn)動,在3～5 s階段,同側(cè)C3電極復(fù)雜度減小,對側(cè)C4電極復(fù)雜度增加;圖4(b)為想象右手運(yùn)動,亦出現(xiàn)了類似現(xiàn)象,故Lempel-Ziv復(fù)雜度符合ERS/ERD現(xiàn)象。

(a)

(b)圖圖4 C3、C4的平均Lempel-Ziv復(fù)雜度幅值

因此,綜合上述分析,三種不同觀測角度的特征都適合用來分析左右手想象運(yùn)動。

腦電特征提取過程如圖5所示。

圖5 特征提取過程

1.2.3 分類方法及性能指標(biāo)

線性判別分類器(Linear Discriminant Analysis,LDA)的基本思想是將訓(xùn)練樣本投影后,使得類內(nèi)方差最小,類間方差最大,即使同種類別數(shù)據(jù)的投影點(diǎn)盡量靠近,不同類別數(shù)據(jù)的類別中心間距離盡可能的大。

LDA實(shí)現(xiàn)步驟[14]如下:

(1) Fisher規(guī)則函數(shù):

(7)

式中:向量Ω為投影方向;Sb是投影前兩類別的類間;SΩ是類內(nèi)離散度矩陣。

(2) 使用拉格朗日乘數(shù)法求解極值,約束條件為函數(shù)的分母為常數(shù)矩陣C:

L(Ω,λ)=ΩTSbΩ-λ(ΩTSΩΩ-C)

(8)

式中:λ為拉格朗日乘子。

(3) 將式(8)對Ω求偏導(dǎo)數(shù),得到最佳投影方向:

(9)

本文引入分類準(zhǔn)確率和互信息兩個指標(biāo)來評價(jià)本文方法的有效性[15]。

分類準(zhǔn)確率定義如下:

(10)

互信息定義如下:

MI=0.5log2(RSN,t)+1

(11)

式中:RSN,t為信噪比。RSN,t表達(dá)式如下:

(12)

2 實(shí)驗(yàn)與結(jié)果分析

2.1 實(shí)驗(yàn)過程

實(shí)驗(yàn)環(huán)境為PC計(jì)算機(jī):Windows 10 64位系統(tǒng),8 GB內(nèi)存,主頻2.60 GHz。首先,依據(jù)運(yùn)動想象所對應(yīng)的大腦區(qū)域,選取C3和C4兩個通道作為輸入信號。對所選通道的腦電信號進(jìn)行DTCWT變換得到不同頻率范圍的分量信號(見表1),根據(jù)ERS/ERD現(xiàn)象,D2與D3重構(gòu)分量接近α和β節(jié)律,以便于精確地對EEG進(jìn)行特征提取。

然后,對D2與D3分量計(jì)算瞬時(shí)能量譜IMS、邊際能量譜MES和Lempel-Ziv復(fù)雜度三個特征參數(shù),從而構(gòu)成12維的時(shí)頻域-非線性特征向量。

最后,采用分類器對腦電特征向量進(jìn)行分類,LDA相較于支持向量機(jī)和神經(jīng)網(wǎng)絡(luò),它不需要設(shè)置任何參數(shù),不依賴隱含層節(jié)點(diǎn)的選擇,在腦電識別方面具有良好的分類效果。因此選用LDA分類器來進(jìn)行腦電想象運(yùn)動的識別。本文采用4折交叉驗(yàn)證方式進(jìn)行50次隨機(jī)實(shí)驗(yàn),將50次實(shí)驗(yàn)的平均識別率作為最終實(shí)驗(yàn)結(jié)果。由于實(shí)驗(yàn)時(shí)間的4～8 s為受試者的表現(xiàn)力最佳時(shí)期,因此選取該時(shí)間段進(jìn)行實(shí)驗(yàn)驗(yàn)證。

2.2 多特征和單一特征的對比

為了方便,把瞬時(shí)能量譜IES記作特征F1,邊際能量譜MES記作特征F2,Lempel-Ziv復(fù)雜度記作特征F3。受試者O3VR記作S1,受試者S4記作S2,受試者X11記作S3。采用不同特征進(jìn)行分類,包括采用單一特征和采用多特征進(jìn)行腦電分類,所得結(jié)果如表2所示。

表2 不同被試者采用不同特征的腦電分類結(jié)果

可以看出,相同實(shí)驗(yàn)任務(wù)下,不同受試者的識別率不盡相同,是每個受試者對實(shí)驗(yàn)的敏感程度不同導(dǎo)致的。F1+F2為時(shí)頻域特征,F3為非線性特征,同一受試者,采用時(shí)域特征進(jìn)行運(yùn)動想象分類的識別率沒有采用非線性特征的識別率高,且三個受試者皆是如此,表明非線性特征能夠更好地表征腦電信號中所蘊(yùn)含的有效信息。本文提出的時(shí)頻-非線性的組合特征(F1+F2+F3),相較于單一特征,識別率均得到了顯著提升,這是由于組合特征對腦電信息進(jìn)行了多角度的描述,完備地表征了腦電信號。表2實(shí)驗(yàn)數(shù)據(jù)表明,多種特征組合的分類效果明顯優(yōu)于單一特征。

2.3 多種識別方法結(jié)果比較

表3給出了本文方法與近期使用相同數(shù)據(jù)集的其他方法所取得的分類結(jié)果對比。文獻(xiàn)[9]從非線性角度選取了相空間特征進(jìn)行運(yùn)動想象識別。文獻(xiàn)[10]選用了小波特征并使用粒子群對分類器進(jìn)行了優(yōu)化。本文則是將時(shí)頻域特征以及非線性特征結(jié)合組成時(shí)頻域-非線性特征向量,從多個角度描述信號信息。由表3可知,本文方法的平均識別率均優(yōu)于其他兩種方法。

表3 不同方法所得到的最高識別率(%)

表4給出了不同方法下得到的最大互信息,其中包括BCI2003競賽前三名獲勝者以及文獻(xiàn)[9-10]的實(shí)驗(yàn)結(jié)果。實(shí)驗(yàn)結(jié)果表明,本文方法相較于其他實(shí)驗(yàn)方法具有顯著優(yōu)勢。

表4 不同方法的最大互信息

3 結(jié) 語

腦電信號是一種隨機(jī)非線性信號,傳統(tǒng)方法只采用單一特征或只考慮了時(shí)頻域角度,而忽略了非線性這一特性,導(dǎo)致系統(tǒng)識別率較低。本文針對這一問題,提出一種基于DTCWT的運(yùn)動想象腦電特征提取方法。首先,運(yùn)用DTCWT近似平移不變性和有效的抗頻譜混疊的特性,對腦電信號進(jìn)行分解與重構(gòu),篩選得到與ERS/ERD現(xiàn)象更加匹配的信號分量。然后對篩選后的分量提取Hilbert瞬時(shí)能量譜IMS、邊際能量譜MES和Lempel-Ziv復(fù)雜度三個特征,形成時(shí)頻域-非線性特征向量。最后用LDA對特征進(jìn)行分類識別。實(shí)驗(yàn)結(jié)果表明,引入非線性特征——Lempel-Ziv復(fù)雜度后,相對單一特征或時(shí)頻域特征,識別率都得到了一定程度的提升,最高分類準(zhǔn)確率為89.84%,平均識別率為82.12%,最大互信息值為0.88,表明了非線性分析結(jié)合時(shí)頻域分析的優(yōu)勢。