莫云 李智 趙紫寧 張紹榮

1 桂林電子科技大學(xué) 電子工程與自動(dòng)化學(xué)院,廣西 桂林 541004; 2 桂林航天工業(yè)學(xué)院 電子信息與自動(dòng)化學(xué)院,廣西 桂林 541004

運(yùn)動(dòng)想象腦電信號(hào)是大腦活動(dòng)自發(fā)產(chǎn)生的腦電信號(hào)，不需要外部的刺激，在腦機(jī)接口(brain computer interface，BCI)中得到了廣泛應(yīng)用。實(shí)現(xiàn)穩(wěn)定可靠的BCI系統(tǒng)關(guān)鍵在于腦電信號(hào)的特征提取。但是運(yùn)動(dòng)想象腦電信號(hào)隨機(jī)性強(qiáng)，信噪比低，容易受生理和非生理噪聲干擾[1]，從運(yùn)動(dòng)想象腦電信號(hào)中提取有效的特征非常困難。因此，特征提取一直是BCI系統(tǒng)研究的重點(diǎn)和難點(diǎn)。

共空域模式(common spatial pattern，CSP)可以有效地捕獲在運(yùn)動(dòng)想象過程中產(chǎn)生的事件相關(guān)同步(event related synchronization，ERS)和事件相關(guān)去同步(event related desynchronization，ERD)現(xiàn)象，從而提取具有判別性的特征，因此CSP在運(yùn)動(dòng)想象腦電(electroencephalogram，EEG)解碼中得到了廣泛的應(yīng)用[2]。但是CSP的性能依賴于最優(yōu)頻帶和時(shí)間窗的選擇[3]。針對(duì)CSP這兩方面的缺陷，已有大量的研究工作對(duì)CSP進(jìn)行了改進(jìn)。Ang等人[4]對(duì)CSP進(jìn)行了改進(jìn)，提出了濾波器組共空域模式(filter bank common spatial patterns，F(xiàn)BCSP)，此方法將原始的腦電信號(hào)濾波成多個(gè)不重疊的子帶成分，再對(duì)每個(gè)子帶成分用CSP算法得到能量特征，然后使用互信息的方法選取最具有判別性的特征用于分類。后面還有一系列的文章對(duì)FBCSP進(jìn)行改進(jìn)和優(yōu)化。近期，Kumar等人[5]對(duì)這一系列的文章做了較為詳細(xì)的論述。文獻(xiàn)[6]使用滑動(dòng)時(shí)間窗技術(shù)選擇最佳的時(shí)間窗。文獻(xiàn)[7]同樣使用滑動(dòng)時(shí)間窗技術(shù)，但是時(shí)間窗寬度可變，通過移動(dòng)時(shí)間窗的起點(diǎn)和終點(diǎn)選擇最優(yōu)的時(shí)間窗，此種方法在訓(xùn)練階段非常耗時(shí)。文獻(xiàn)[3]和文獻(xiàn)[8]聯(lián)合考慮了時(shí)間窗和頻帶的選擇。值得指出的是，以上方法需要把腦電信號(hào)濾波成多個(gè)頻帶和時(shí)間窗，工作量大，計(jì)算耗時(shí)。

本文從另外一個(gè)角度出發(fā)，對(duì)CSP特征提取方法進(jìn)行改進(jìn)。傳統(tǒng)的CSP特征提取方法在進(jìn)行CSP變換之后提取信號(hào)的對(duì)數(shù)方差作為特征。本文假設(shè)腦電信號(hào)在CSP變換之后，有比對(duì)數(shù)方差更具判別性的特征?；诖讼敕?，我們提出了基于CSP變換的新特征提取方法。首先，對(duì)原始的腦電信號(hào)進(jìn)行8～30 Hz帶通濾波，濾除與運(yùn)動(dòng)想象無(wú)關(guān)的腦電成分。其次，對(duì)濾波的信號(hào)進(jìn)行CSP變換。然后，對(duì)CSP變換后的信號(hào)提取對(duì)數(shù)方差、自回歸(auto regressive，AR)系數(shù)、帶通功率和小波能量以及融合特征作為特征。最后，使用Fisher線性判別分析(Fisher Linear Discriminant Analysis，F(xiàn)LDA)、最小絕對(duì)值收縮和選擇算子(least absolute shrinkage and selection operator，LASSO)、融合LASSO(fused LASSO)和彈性網(wǎng)絡(luò)(elastic net，EN)四種分類器進(jìn)行分類。實(shí)驗(yàn)結(jié)果表明基于CSP變換的新特征提取方法提高了分類準(zhǔn)確率，優(yōu)于傳統(tǒng)的CSP特征提取方法。

1 方法

圖1給出了本文方法的整體流程圖。首先，對(duì)輸入的原始EEG信號(hào)進(jìn)行預(yù)處理，具體預(yù)處理過程將在下文詳細(xì)介紹。其次，對(duì)預(yù)處理的數(shù)據(jù)進(jìn)行CSP變換。然后進(jìn)行特征提取，所提取特征包括對(duì)數(shù)方差、AR系數(shù)、帶通功率和小波能量，并計(jì)算4種特征的融合特征。最后，采用FLDA、LASSO、fused LASSO和EN分類器進(jìn)行分類。

圖1 整體流程圖

1.1 預(yù)處理流程

對(duì)輸入的EEG原始信號(hào)進(jìn)行預(yù)處理，詳細(xì)處理方法如下：

(1)使用共同平均參考(common average reference，CAR)對(duì)原始信號(hào)進(jìn)行空間濾波，降低容積傳導(dǎo)效應(yīng)的影響。

(2)使用6階的Butterworth濾波器對(duì)每個(gè)通道的腦電信號(hào)進(jìn)行8～30 Hz的帶通濾波，濾除與運(yùn)動(dòng)想象無(wú)關(guān)的腦電成分。

(3)提取單試次數(shù)據(jù)。數(shù)據(jù)集選取的時(shí)間窗為0.5～3.5 s，其中0 s表示運(yùn)動(dòng)想象任務(wù)開始的時(shí)間。

1.2 CSP變換

CSP尋找一組空間濾波器，使投影濾波后的腦電信號(hào)在一類運(yùn)動(dòng)想象任務(wù)的方差最大，而在另一類任務(wù)的方差最小。由于在給定頻帶內(nèi)的腦電信號(hào)方差對(duì)應(yīng)該頻帶內(nèi)的信號(hào)功率，因此CSP也可以認(rèn)為是基于頻帶功率特征實(shí)現(xiàn)對(duì)腦電信號(hào)的識(shí)別。

空間濾波器w的計(jì)算是通過對(duì)兩類樣本協(xié)方差矩陣同時(shí)進(jìn)行對(duì)角化，具體如下：

(1)

(2)

式(2)中，trace(·)表示矩陣的跡的求解；Nk表示第k類任務(wù)的樣本數(shù)，即第k類任務(wù)的單試次數(shù)據(jù)個(gè)數(shù)；D(k,n)∈C×T代表第k類任務(wù)的第n個(gè)試次的數(shù)據(jù)，其中C為腦電信號(hào)通道總數(shù)，T為各通道樣點(diǎn)數(shù)。

公式(1)可轉(zhuǎn)化為式(3)廣義特征值的求解問題[9]。

(3)

對(duì)新的單試次數(shù)據(jù)，空間投影信號(hào)為：

Z=WTD

(4)

傳統(tǒng)的CSP特征提取方法，提取空間濾波信號(hào)的對(duì)數(shù)方差作為特征，具體如下：

(5)

式(5)中，var(·)表示信號(hào)的方差。最后通過公式(5)可以計(jì)算得到單試次數(shù)據(jù)的特征向量，表示為x=[f1,f2,…,f2m]。

1.3 新特征提取方法

在1.2節(jié)中，通過式(4)得到空間濾波信號(hào)Z后，分別提取信號(hào)Z的AR系數(shù)、帶通功率和小波能量作為特征，并計(jì)算4種特征(包含對(duì)數(shù)方差)的融合特征。

(1)AR系數(shù)

自回歸(Autoregressive Model，AR)模型是利用自身做回歸變量的過程，即用前期若干時(shí)刻隨機(jī)變量的線性組合來(lái)描述后期某時(shí)刻隨機(jī)變量的線性回歸模型，是時(shí)間序列中一種常見的形式。AR模型通常表示為

(6)

式(6)中，y(t)是隨機(jī)時(shí)間序列。φi是模型參數(shù)，εt是白噪聲序列，p為AR模型階數(shù)。

(2)帶通功率

帶通功率使用文獻(xiàn)[11]的方法計(jì)算，功率譜估計(jì)使用Welch方法，窗口函數(shù)使用海明窗，濾波頻帶為：8～12 Hz、12～16 Hz、16～20 Hz、20～24 Hz和24～30 Hz。信號(hào)序列的帶通功率計(jì)算如下：

(7)

式(7)中，pi,j[t]表示第i個(gè)通道，第j個(gè)頻帶的信號(hào)序列的帶通功率。計(jì)算各個(gè)頻帶的平均功率作為特征，記為pi,j。由于本文選擇了5個(gè)濾波頻帶，所以每個(gè)通道的帶通功率特征維數(shù)為5。

(3)小波能量

本文使用小波包分解(wavelet packet decomposition，WPD)對(duì)信號(hào)Z的每個(gè)通道進(jìn)行小波分解，WPD的詳細(xì)推導(dǎo)公式可以參考文獻(xiàn)[12]。針對(duì)所選數(shù)據(jù)集，選用‘db4’小波基，進(jìn)行3層小波分解。對(duì)各個(gè)子帶的小波系數(shù)計(jì)算小波能量。子帶j的小波能量計(jì)算公式為：

(8)

其中，N為相應(yīng)子帶小波系數(shù)的個(gè)數(shù)，di,j為第j個(gè)子帶的第i個(gè)小波系數(shù)。

(4)融合特征計(jì)算

把腦電信號(hào)每個(gè)通道的特征(包括對(duì)數(shù)方差、AR系數(shù)、帶通功率和小波能量)逐行進(jìn)行連接，最終得到一個(gè)特征向量，該特征向量的維數(shù)等于每通道的特征維數(shù)乘以通道總數(shù)。

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

采用一個(gè)公開的BCI競(jìng)賽數(shù)據(jù)集驗(yàn)證所提出方法的性能，并選用FLDA、LASSO、fused LASSO和EN分類器進(jìn)行實(shí)驗(yàn)。其中，LASSO、fused LASSO和EN分類器具有特征選擇功能。LASSO、fused LASSO和EN分類器模型參數(shù)通過10重交叉驗(yàn)證進(jìn)行選擇，LASSO和EN分類器超參數(shù)備選集合均為：λ∈0.1×{1,2,…,30}，fused LASSO分類器的超參數(shù)備選集合為：λ1,λ2∈0.1×{1,2,…,30}。FLDA參考文獻(xiàn)[13]實(shí)現(xiàn)，LASSO、fused LASSO和EN使用SLEP工具箱[14]實(shí)現(xiàn)。

2.1 數(shù)據(jù)說(shuō)明

下面將簡(jiǎn)要介紹所用到的數(shù)據(jù)集：

BCI競(jìng)賽III(2005)數(shù)據(jù)集IVa[15]。該數(shù)據(jù)集共包括118個(gè)電極通道，采樣率為100 Hz。5個(gè)健康被試(aa、al、av、aw、ay)分別執(zhí)行左、右手和右腳3類運(yùn)動(dòng)想象任務(wù)。由于競(jìng)賽只提供了右手和右腳的運(yùn)動(dòng)想象腦電數(shù)據(jù)，因此只對(duì)這兩類任務(wù)數(shù)據(jù)進(jìn)行分類。每個(gè)被試包含280個(gè)單試次的腦電數(shù)據(jù)，其中aa、al、av、aw、ay各個(gè)被試的訓(xùn)練集樣本個(gè)數(shù)分別為168、224、84、56、28，其余為測(cè)試集。

2.2 實(shí)驗(yàn)結(jié)果

在使用AR系數(shù)作為特征時(shí)，AR系數(shù)的階數(shù)選擇對(duì)最終的分類結(jié)果影響比較大。因此，我們通過實(shí)驗(yàn)選擇AR系數(shù)特征的最優(yōu)階數(shù)。圖2給出了使用AR系數(shù)作為特征時(shí)，采用FLDA等4種分類方法得到的整個(gè)數(shù)據(jù)集平均分類準(zhǔn)確率隨AR階數(shù)(1至20階)的變化曲線。從圖2中可以明顯地看出，對(duì)于4種分類器，平均分類準(zhǔn)確率均在階數(shù)為4左右時(shí)，達(dá)到最高值。其中，F(xiàn)LDA、LASSO和fused LASSO分類器在4階時(shí)，分類準(zhǔn)確率最高；EN在3階時(shí)，分類準(zhǔn)確率最高?？偲骄诸悳?zhǔn)確率先隨著階數(shù)的增加而增加，在階數(shù)為4時(shí)達(dá)到峰值，之后隨著階數(shù)的增加準(zhǔn)確率連續(xù)下降。通過圖2我們可以看出，當(dāng)使用AR系數(shù)作為特征時(shí)，AR階數(shù)為4時(shí)，效果最佳。因此，在本文的后續(xù)實(shí)驗(yàn)中，AR階數(shù)選擇為4。

圖2 平均準(zhǔn)確率隨AR階數(shù)的變化曲線

為了準(zhǔn)確地比較各種特征提取方法的分類效果，表1-4分別展示了在不同分類方法下每個(gè)被試對(duì)應(yīng)不同特征提取方法得到的分類準(zhǔn)確率。表中準(zhǔn)確率最高的數(shù)據(jù)用黑體加粗標(biāo)注。

表1 (FLDA)分類準(zhǔn)確率

表2 (LASSO)分類準(zhǔn)確率

表3 (fused LASSO)分類準(zhǔn)確率

表4 (EN)分類準(zhǔn)確率

表1的分類方法為FLDA。在表1中可以看出，AR系數(shù)特征的平均分類準(zhǔn)確率最高，并且優(yōu)于傳統(tǒng)的CSP特征提取方法，即對(duì)數(shù)方差特征。其他特征提取方法低于對(duì)數(shù)方差特征，其中帶通功率特征的平均分類準(zhǔn)確率最低。

表2的分類方法為L(zhǎng)ASSO。在表2中可以看出，AR系數(shù)特征的平均分類準(zhǔn)確率最高，AR系數(shù)特征和融合特征都優(yōu)于對(duì)數(shù)方差特征。但是，融合特征的平均分類準(zhǔn)確率幾乎和對(duì)數(shù)方差特征的一致。

表3的分類方法為fused LASSO。在表3中可以看出，AR系數(shù)特征的平均分類準(zhǔn)確率最高，并且優(yōu)于對(duì)數(shù)方差特征。其他特征提取方法低于對(duì)數(shù)方差特征。

表4的分類方法為EN。在表4中可以看出，AR系數(shù)特征的平均分類準(zhǔn)確率最高，其次是融合特征，這兩種特征提取方法的分類結(jié)果都優(yōu)于對(duì)數(shù)方差特征。其他特征提取方法低于對(duì)數(shù)方差特征。

對(duì)比表1-4的結(jié)果，無(wú)論使用哪種分類器，帶通功率特征的分類效果都是最差的，而AR系數(shù)特征則是最佳的。當(dāng)分類方法具備特征選擇功能時(shí)，融合特征優(yōu)于對(duì)數(shù)方差特征，比如LASSO和EN。但是使用fused LASSO分類器時(shí)，融合特征的分類準(zhǔn)確率低于對(duì)數(shù)方差特征。其主要原因是：fused LASSO分類器可以選擇具有局部光滑特性的特征，但是我們所提取融合特征沒有這樣的特性。

為了更加直觀地查看不同特征提取方法使用不同分類器所取得的分類效果，圖3給出了不同特征提取方法使用不同分類器的情況下取得的平均類準(zhǔn)確率。圖3中可以清晰地看到AR系數(shù)特征效果最佳，其次是對(duì)數(shù)方差和融合特征。帶通功率特征提取方法的分類效果最差，且在不同分類方法的分類結(jié)果中差異較大。

圖3 不同特征提取方法的平均準(zhǔn)確率

3 討論

當(dāng)使用AR系數(shù)作為特征時(shí)，AR階數(shù)對(duì)準(zhǔn)確率的影響比較大，所以選擇合適的AR階數(shù)比較重要。在使用不同數(shù)據(jù)集進(jìn)行實(shí)驗(yàn)時(shí)，最佳的AR階數(shù)有可能不同，需要實(shí)驗(yàn)進(jìn)行選擇。從圖3來(lái)看，AR系數(shù)在分類性能上相較于對(duì)數(shù)方差、帶通功率、小波能量和融合特征，對(duì)各分類器有更好的適應(yīng)性。

另外，AR系數(shù)特征使用4種分類方法得到的平均準(zhǔn)確率都是最高的，且都維持在82%～83%之間。驗(yàn)證了使用新特征來(lái)代替?zhèn)鹘y(tǒng)的對(duì)數(shù)方差特征的可行性和有效性。

表4中使用EN作為分類器，AR系數(shù)雖然選擇了4階，沒有選擇準(zhǔn)確率最高的3階，但仍在所有特征中保持最高的分類準(zhǔn)確率，且在aa、al、av三個(gè)被試中均取得最高的分類準(zhǔn)確率，進(jìn)一步驗(yàn)證了提取AR系數(shù)作為特征的優(yōu)勢(shì)。

4 總結(jié)

基于CSP變換，本文使用AR系數(shù)、帶通功率、小波能量和融合特征代替?zhèn)鹘y(tǒng)CSP提取方差，并選用FLDA、LASSO、fused LASSO和EN四種分類方法進(jìn)行分類。實(shí)驗(yàn)結(jié)果表明AR系數(shù)作為新的特征可以代替對(duì)數(shù)方差特征取得更好的分類效果，其他特征提取方法則效果不佳。AR系數(shù)特征顯著地提高了腦電信號(hào)的分類準(zhǔn)確率，驗(yàn)證了新特征提取方法的可行性。相比其他CSP改進(jìn)方法，本文提出的方法計(jì)算量少，實(shí)時(shí)性較好，可以有效地提高BCI系統(tǒng)的性能。

在未來(lái)的工作中，我們將進(jìn)一步研究新的特征提取方法。另外，本文只選取了一個(gè)數(shù)據(jù)集作為實(shí)驗(yàn)，可能得到的分類結(jié)果具有一定的片面性。在后面的工作中，我們將使用更多的數(shù)據(jù)集進(jìn)行實(shí)驗(yàn)，進(jìn)一步驗(yàn)證方法的有效性。