湯 偉 耿逸飛
(陜西科技大學電氣與控制工程學院 陜西 西安 710021) (陜西科技大學工業(yè)自動化研究所 陜西 西安 710021)
腦機接口(Brain-Computer Interface,BCI)是利用大腦所產(chǎn)生的信號來建立人腦與外部設備的信息控制的交互技術[1],在諸多方向有著顯著的研究意義,如康復醫(yī)療[2]、游戲娛樂[3]等。運動想象是人對自身進行控制的意念活動。當人在想象單側肢體運動時,對側大腦皮層處于激活狀態(tài),相關的節(jié)律信號幅值會降低;而同側處于阻滯狀態(tài),相關的節(jié)律信號幅值會升高。上述現(xiàn)象分別稱為事件相關去同步(Event-Related Desynchronization,ERD)以及事件相關同步[4](Event-Related Synchronization,ERS)。因此,在ERS/ERD的基礎上,我們能夠?qū)崿F(xiàn)對左右手運動想象任務的分類。由于大腦系統(tǒng)十分復雜,所產(chǎn)生的信號隨機性強且不穩(wěn)定,非線性特征顯著,因此如何有效提取腦電特征成了一個難題。
目前,腦電信號常用的特征提取方法有:單一時域或頻域方法、時頻域結合方法和非線性分析方法等。其中前兩種方法適用于分析平穩(wěn)信號,應用于非平穩(wěn)的腦電信號時效果不佳[4]。由于腦電信號具有顯著的非線性特征,非線性分析方法十分適用于腦電信號,常用的方法有:香農(nóng)熵、自回歸(Auto Regressive,AR)模型參數(shù)、分形維數(shù)和Lempel-Ziv復雜度等[5]。孫會文等[6]提出的結合HHT與AR模型參數(shù)的腦電分類方法,分類識別率達81.08%,適應性良好。羅志增等[7]提出的基于多尺度Lempel-Ziv復雜度的特征提取方法,得到了較高的分類識別率,達87.87%。上述方法從非線性角度對腦電進行識別,都表現(xiàn)出了較高的識別率,但是特征選取考慮的都是單一角度,忽略了腦電信號在時頻域所包含的有效信息[8]。因此本文提出一種時頻域特征與非線性特征相融合的提取方法,從多個角度完備地描述腦電信號。
由于原始腦電中夾雜與運動想象無關的信號頻段,本文引入雙樹復小波(Dual-tree Complex Wavelet Transform,DTCWT)對C3、C4通道的腦電信號進行最優(yōu)頻段選取,這樣既去除了無關腦電頻段,又進行頻率分析,同時提高了分類準確率。在DTCWT分解基礎上,選取Hilbert瞬時能量譜、邊際能量譜和Lempel-Ziv復雜度作為腦電信號特征,最后采用線性判別分析對其實行分類。仿真實驗表明,本文提出的基于DTCWT的多特征融合方法,相較于使用相同數(shù)據(jù)集的文獻[9-10]方法,具有更大的最大互信息與更高的平均識別率。
本文使用的數(shù)據(jù)為被測試者在進行左右手運動想象任務時采集的腦電信號,源自第三屆國際腦-機接口(BCI Competition Ⅲ)的dataset Ⅲb。EEG信號按照國際10-20導聯(lián)系統(tǒng)的C3、C4、Cz通道采集獲得。該數(shù)據(jù)集包含3名受試者的實驗結果,分別為O3VR、S4和X11。O3VR為400個試次,S4和X11皆為1 080個試次,每次實驗時間為8 s,采樣頻率為125 Hz。單次實驗過程如圖1所示。
圖1 單次實驗過程
受試者安靜坐在顯示器前,根據(jù)屏幕提示進行左右手運動想象任務。0~2 s,屏幕空白,受試者進行放松;2~3 s,屏幕投放“+”,準備運動想象任務;3~9 s,屏幕顯示左右方向箭頭,受試者根據(jù)對應箭頭進行運動想象任務。t=9 s時,實驗結束。
1.2.1 雙樹復小波變換
雙樹復小波不僅具有離散小波變換的多尺度分析的能力,同時有效解決了復小波變換中計算復雜的問題,相對于傳統(tǒng)離散小波,可以較為精確地描述信號的細節(jié)信息,已成功應用于軸承故障分析和面部識別等領域[11]。
DTCWT采用實部樹和虛部樹對信號進行分解與重構,實部樹和虛部樹采用不同的實數(shù)濾波器組,從而得到實部系數(shù)與虛部系數(shù),使得實部樹與虛部樹的信息互補,從而保證了信號分析的完備性。另外,DTCWT在多尺度分解中運用了二分法來降低冗余度,有效提高了算法的計算速度。
設腦電信號c(n)的采樣頻率為fs,用DTCWT將腦電信號分為L層,則能得到CAL、CDL、CDL-1、…、CD1等L+1子帶的復系數(shù)。重構后得到的信號分量AL、DL、DL-1、…、D1所對應的頻帶范圍分別為[0,fs/2L+1]、[fs/2L+1,fs/2L]、[fs/2L,fs/2L-1]、…、[fs/22,fs/2]。因此,DTCWT進行L層分解可得L+1個重構分量。
1.2.2 時頻-非線性特征提取
在對腦電信號進行特征提取時,先進行DTCWT分解,得到一系列子頻帶,ERS/ERD現(xiàn)象主要發(fā)生在α(8~14 Hz)和β(14~30 Hz)節(jié)律,根據(jù)DTCWT的分解規(guī)則,分解尺度定為4層。子帶和節(jié)律信號的關系如表1所示。
表1 子帶信號和節(jié)律信號的對應關系
考慮到ERS/ERD現(xiàn)象,本實驗只對D3與D2分量進行特征值求取,這樣既能精確且完備地提取有效信息,又能降低特征向量維度。分別計算每個數(shù)據(jù)段的Hilbert瞬時能量譜(Instantaneous Energy Spectrum,IES)、邊際能量譜(Marginal Energy Spectrum,MES)和Lempel-Ziv復雜度三個特征參數(shù),從而構成時頻域-非線性特征向量。對于單次實驗數(shù)據(jù),對C3和C4兩個電極需要計算4個IES參數(shù)、4個MES參數(shù)、4個Lempel-Ziv復雜度參數(shù),共計12個參數(shù)值組成12維特征向量。
3個特征參數(shù)的定義如下:
(1) 對于一個給定信號C(t),進行Hilbert變換[12]:
(1)
可得到解析信號:
Z(t)=C(t)+jY(t)=A(t)ejφ(t)
(2)
式中:A(t)和φ(t)分別為瞬時幅值和瞬時相位。
由A(t)和φ(t)求取瞬時頻率ω(t):
(3)
則信號幅度在時-頻域的分布情況,即Hilbert譜為:
(4)
式中:Re為取實部。
根據(jù)式(4)可求取IES和MES:
(5)
式中:[ω1,ω2]為頻率范圍;[t1,t2]為時間范圍。
IES表現(xiàn)為時域能量特征,MES表現(xiàn)為頻域能量特征。
(2) Lempel-Ziv復雜度描述了序列與隨機序列的相似度,復雜度越大,序列新模式增加的速率越快,越接近隨機序列[13]。設序列S(s1,s2,…,sn)的Lempel-Ziv復雜度為z(n),當n→∞時,z(n)會趨于定值n/logln,l為粗粒化段數(shù)(傳統(tǒng)二值化時,l=2),則定義歸一化Lempel-Ziv復雜度為:
(6)
為了驗證上述三個特征能夠較好地反映ERS/ERD現(xiàn)象,從流程的第1 s開始,采用滑動窗來提取特征,每次滑動1個采樣點,直到最后1 s結束,從而獲得單次實驗任務的特征值時間序列。將同一運動想象任務的所有實驗所得的特征值疊加平均,得到平均特征值時間序列。圖2(a)為被試者O3VR左手運動想象的平均瞬時能量圖,可以看出,在4~5 s的時間段內(nèi),被試者進行左手想象任務時,C3通道較C4通道具有更高的能量分布;圖2(b)為被試者O3VR右手運動想象的平均瞬時能量圖,進行右手想象任務時,C4通道較C3通道具有更高的能量分布,符合ERS/ERD原理。圖3為左右手運動想象的平均邊際譜幅值圖,圖3(a)為左手,圖3(b)為右手。可以看出在8~15 Hz頻段,兩種不同動作的運動想象下,皆是C3通道具有較高能量,不符合ERS/ERD現(xiàn)象,但在15~30 Hz頻段,信號的特征符合ERS/ERD現(xiàn)象。若從時域能量角度統(tǒng)一分析信號特征,低能高頻信號所含有效信息會被高能低頻信號覆蓋,從而導致時域能量特征的敏感度下降、被試識別率差的結果。因此,本文利用DTCWT將信號拆分為成多個子頻帶分量,保證了有效信息的完整保留,提高了被試識別率。
(a)
(b)圖2 C3、C4的平均瞬時能量
(a)
(b)圖3 C3、C4的平均邊際譜幅值
圖4為左右手運動想象的平均Lempel-Ziv復雜度幅值圖。圖4(a)為想象左手運動,在3~5 s階段,同側C3電極復雜度減小,對側C4電極復雜度增加;圖4(b)為想象右手運動,亦出現(xiàn)了類似現(xiàn)象,故Lempel-Ziv復雜度符合ERS/ERD現(xiàn)象。
(a)
(b)圖圖4 C3、C4的平均Lempel-Ziv復雜度幅值
因此,綜合上述分析,三種不同觀測角度的特征都適合用來分析左右手想象運動。
腦電特征提取過程如圖5所示。
圖5 特征提取過程
1.2.3 分類方法及性能指標
線性判別分類器(Linear Discriminant Analysis,LDA)的基本思想是將訓練樣本投影后,使得類內(nèi)方差最小,類間方差最大,即使同種類別數(shù)據(jù)的投影點盡量靠近,不同類別數(shù)據(jù)的類別中心間距離盡可能的大。
LDA實現(xiàn)步驟[14]如下:
(1) Fisher規(guī)則函數(shù):
(7)
式中:向量Ω為投影方向;Sb是投影前兩類別的類間;SΩ是類內(nèi)離散度矩陣。
(2) 使用拉格朗日乘數(shù)法求解極值,約束條件為函數(shù)的分母為常數(shù)矩陣C:
L(Ω,λ)=ΩTSbΩ-λ(ΩTSΩΩ-C)
(8)
式中:λ為拉格朗日乘子。
(3) 將式(8)對Ω求偏導數(shù),得到最佳投影方向:
(9)
本文引入分類準確率和互信息兩個指標來評價本文方法的有效性[15]。
分類準確率定義如下:
(10)
互信息定義如下:
MI=0.5log2(RSN,t)+1
(11)
式中:RSN,t為信噪比。RSN,t表達式如下:
(12)
實驗環(huán)境為PC計算機:Windows 10 64位系統(tǒng),8 GB內(nèi)存,主頻2.60 GHz。首先,依據(jù)運動想象所對應的大腦區(qū)域,選取C3和C4兩個通道作為輸入信號。對所選通道的腦電信號進行DTCWT變換得到不同頻率范圍的分量信號(見表1),根據(jù)ERS/ERD現(xiàn)象,D2與D3重構分量接近α和β節(jié)律,以便于精確地對EEG進行特征提取。
然后,對D2與D3分量計算瞬時能量譜IMS、邊際能量譜MES和Lempel-Ziv復雜度三個特征參數(shù),從而構成12維的時頻域-非線性特征向量。
最后,采用分類器對腦電特征向量進行分類,LDA相較于支持向量機和神經(jīng)網(wǎng)絡,它不需要設置任何參數(shù),不依賴隱含層節(jié)點的選擇,在腦電識別方面具有良好的分類效果。因此選用LDA分類器來進行腦電想象運動的識別。本文采用4折交叉驗證方式進行50次隨機實驗,將50次實驗的平均識別率作為最終實驗結果。由于實驗時間的4~8 s為受試者的表現(xiàn)力最佳時期,因此選取該時間段進行實驗驗證。
為了方便,把瞬時能量譜IES記作特征F1,邊際能量譜MES記作特征F2,Lempel-Ziv復雜度記作特征F3。受試者O3VR記作S1,受試者S4記作S2,受試者X11記作S3。采用不同特征進行分類,包括采用單一特征和采用多特征進行腦電分類,所得結果如表2所示。
表2 不同被試者采用不同特征的腦電分類結果
可以看出,相同實驗任務下,不同受試者的識別率不盡相同,是每個受試者對實驗的敏感程度不同導致的。F1+F2為時頻域特征,F3為非線性特征,同一受試者,采用時域特征進行運動想象分類的識別率沒有采用非線性特征的識別率高,且三個受試者皆是如此,表明非線性特征能夠更好地表征腦電信號中所蘊含的有效信息。本文提出的時頻-非線性的組合特征(F1+F2+F3),相較于單一特征,識別率均得到了顯著提升,這是由于組合特征對腦電信息進行了多角度的描述,完備地表征了腦電信號。表2實驗數(shù)據(jù)表明,多種特征組合的分類效果明顯優(yōu)于單一特征。
表3給出了本文方法與近期使用相同數(shù)據(jù)集的其他方法所取得的分類結果對比。文獻[9]從非線性角度選取了相空間特征進行運動想象識別。文獻[10]選用了小波特征并使用粒子群對分類器進行了優(yōu)化。本文則是將時頻域特征以及非線性特征結合組成時頻域-非線性特征向量,從多個角度描述信號信息。由表3可知,本文方法的平均識別率均優(yōu)于其他兩種方法。
表3 不同方法所得到的最高識別率(%)
表4給出了不同方法下得到的最大互信息,其中包括BCI2003競賽前三名獲勝者以及文獻[9-10]的實驗結果。實驗結果表明,本文方法相較于其他實驗方法具有顯著優(yōu)勢。
表4 不同方法的最大互信息
腦電信號是一種隨機非線性信號,傳統(tǒng)方法只采用單一特征或只考慮了時頻域角度,而忽略了非線性這一特性,導致系統(tǒng)識別率較低。本文針對這一問題,提出一種基于DTCWT的運動想象腦電特征提取方法。首先,運用DTCWT近似平移不變性和有效的抗頻譜混疊的特性,對腦電信號進行分解與重構,篩選得到與ERS/ERD現(xiàn)象更加匹配的信號分量。然后對篩選后的分量提取Hilbert瞬時能量譜IMS、邊際能量譜MES和Lempel-Ziv復雜度三個特征,形成時頻域-非線性特征向量。最后用LDA對特征進行分類識別。實驗結果表明,引入非線性特征——Lempel-Ziv復雜度后,相對單一特征或時頻域特征,識別率都得到了一定程度的提升,最高分類準確率為89.84%,平均識別率為82.12%,最大互信息值為0.88,表明了非線性分析結合時頻域分析的優(yōu)勢。