面向腦機接口個體差異的濾波器組設計方法及其模型研究

李子波，唐榮年，謝小峰

（海南大學機電工程學院，海南 ?？?570228）

腦機接口（Brain?Computer Interfaces，BCI）將腦信號轉換為可傳達用戶意圖的輸出［1］，為大腦控制外部設備提供了一種新穎的連接方式，例如輪椅，假肢和計算機［2］.此項研究能夠為患有嚴重運動障礙的人開發(fā)出一種新的交流和控制的方式，有助于提高運動障礙者的生活［3］.不同的BCI系統(tǒng)設計出幾種腦電信號（Electroencephalogram，EEG）模式，其中基于運動想象（Motor Imagery，MI）的BCI系統(tǒng)可以將想象身體運動轉換為控制命令［4］.在處理運動想象的EEG信號方法中，共空域模式（Common Spatial Pattern，CSP）是非常有效的特征提取方法之一［5-6］.CSP算法不僅可以從運動想象的信號中提取出最具有判別性的特征，而且可以提供EEG信號通道選擇的權重信息［7］，但CSP算法的效果受濾波器篩選的頻帶信息影響.值得注意的是，在運動想象任務中，在特定的頻帶范圍蘊含非常重要的EEG信號信息，例如alpha波段（8～13 Hz），beta波段（14～30 Hz）.因此，選擇合適的頻帶信息尤為重要.

針對CSP算法和CSP變體算法的研究受到廣泛的關注，濾波器組共空域模式（Filter Bank CSP，F(xiàn)BC?SP）使用零相位的切比雪夫Ⅱ型濾波器將原始信息濾波為多個子帶，并使用互信息系數(shù)自動選擇特征［7］；子帶共空域模式（Sub?Band CSP，SBCSP）與FBCSP相似，在信號的多個子帶上使用Gabor濾波器，然后基于子帶分數(shù)融合技術選擇特征［8］；自適應濾波器組共空域模式（Adaptive Filter Bank CSP，AFBCSP）可以基于Fisher比率來自適應選擇特定于對象的區(qū)分性頻帶［9］.雖然上述研究方法在濾波器組設計時都包含了al?pha和beta波段，但是濾波器組劃分出所有的子頻帶都是相等權重，即與運動想象相關的頻帶信息沒有被突出表達.由于個體的差異，導致不同實驗對象的重要頻帶范圍會有波動，如圖1所示，展示了2個實驗對象在想象右手運動的時頻圖存在差異.圖1紅色框選擇范圍表示重要頻帶，由于個體區(qū)別使得alpha波段和beta波段在每個個體上的表達存在差異，從而影響信息的表達.

圖1 不同對象在運動想象時出現(xiàn)的頻段不一致現(xiàn)象

綜上所述，如何突出表達與運動想象相關的頻帶信息以及消除由于個體差異導致重點頻帶范圍波動尤為重要.筆者針對上述問題設計出一種新穎的濾波器組設計方法，不僅能夠突出表達alpha和beta波段的信息，而且還能夠在一定程度上改善由于個體差異導致的重要頻帶范圍波動的影響，提升模型的穩(wěn)定性.

1 試驗方法

1.1數(shù)據(jù)來源實驗中使用的n道EEG數(shù)據(jù)為BCI競賽第IV期的運動圖像數(shù)據(jù).BCI競賽IV的數(shù)據(jù)集IIa來自9名受試者，執(zhí)行4種類型的運動想象任務（右手、左手、腳和舌頭想象運動）.記錄的信號由22個腦電圖通道組成，通道構型如圖2a所示.實驗方案如圖2b所示：

圖2 電極分布與時序方案圖

1）在初始時間（0～2 s），發(fā)出短的聲學警告音；

2）2 s后，屏幕上出現(xiàn)了一個指向左、右、下或上的箭頭提示，時間為2～3.25 s，激發(fā)受試者進行運動想象任務，直到6 s固定十字從屏幕上消失；

3）有1.5 s短暫休息.

處理數(shù)據(jù)的時間間隔限定在被試執(zhí)行運動想象任務的3.75～5.75 s之間.對于每個受試者和運動想象任務，有72項訓練和72項測試試驗.因此，每個受試者的訓練試驗和測試試驗總數(shù)都為288次，EEG信號采集，采樣頻率為250 Hz.

1.2分析方法設計一種新的濾波器組方法，如圖3預處理部分所示，將原始的EEG信號劃分為多個子頻帶.在粗特征提取部分，使用CSP算法分別對多個子頻帶提取特征.雖然該濾波器組豐富了數(shù)據(jù)表達，同時也會導致信息冗余.采取3種精特征篩選方法（組稀疏（Group Lasso）、基于極限學習機的自編碼器（ELM_AE）、互信息（Mutual Information））進一步處理，最后使用LDA建立模型.

圖3 濾波器組設計方法結構圖

1.2.1濾波器組設計方法基于運動想象的多通道的腦電信號可以表示

L和N分別表示采樣點個數(shù)和通道數(shù)，

設計了一種全新的濾波器組方法，以4 Hz為間隔將4～42 Hz頻帶范圍信號濾波為150個子頻帶，如圖4a所示，即［4～8 Hz，4～12 Hz，…，4～40 Hz，3～7 Hz］；［6～10 Hz，6～14 Hz，…，6～42 Hz，35～39 Hz］；…；［38～42 Hz，3～7 Hz，…，3～35 Hz，3～39 Hz］.從圖4b可以看出，alpha和beta波段占據(jù)很大比例，使重要頻帶信息得以突出表達，并且2個核心頻帶周圍的信息也被強調，可以在一定程度減緩由個體差異產生的核心波段波動的影響.

圖4 濾波器組結構圖與頻帶占比圖

1.2.2CSP算法粗特征提取通過上述濾波器組將原始腦電信號分解為多個子頻帶，使用空域濾波器矩陣W分別從每個子頻帶中提取特征.對于二分類問題.CSP算法可以在最大化一類方差的同時最小化另外一類的方差，因此空域濾波器可以通過最大化最小化目標函數(shù)得到［7］

其中，Ci表示屬于第i類EEG數(shù)據(jù)的協(xié)方差矩陣的算術平均值.

F表示CPS特征

1.2.33種特征精篩選方法多個子頻帶的CSP粗特征提取會增加冗余信息，因此分別采取3種降維和特征選擇的方法進行特征精篩選.

自編碼器算法降維 極限學習機（ELM）是單隱藏層前饋網(wǎng)絡（SLFNS），盡管采用了隨機投影，但在大多數(shù)情況下，ELM的性能都優(yōu)于反向傳播神經(jīng)網(wǎng)絡［10］.與一般的多層感知機相比，ELM具有非梯度優(yōu)化的特點，并且計算成本通常要低得多［11］，得益于其優(yōu)秀性能，近年來針對ELM的研究被人們廣為關注［12-14］.由于多頻帶所提取出的特征信息的高維度以及存在冗余，采用基于ELM的自動編碼器網(wǎng)絡對多頻帶CPS

特征進行編碼壓縮，ELM_AE的表達形式為

其中，h(fi)∈RNh×1表示隱含層為Nh的輸出，a和b表示激活函數(shù)的參數(shù)，在［-1，1］范圍內隨機生成.AE網(wǎng)絡的解碼部分

其中，β=H?F解出.因此可以獲得自編碼器降維的特征S_AE

互信息熵特征選擇 給定初始數(shù)據(jù)集合F（n維特征），尋找一個子集S屬于F并且擁有從原始數(shù)據(jù)集合篩選出的K個最大互信息的特征［7］.2個隨機變量的互信息定義

對于d維隨機變量X，可以表示為X={x1，x2，…，xd}，其信息熵計算

對于隨機變量X和Y的條件熵計算

以互信息熵值作為篩選依據(jù)，選擇前K個互信息熵的最大的特征構成特征子集S_MI.

組稀疏特征選擇 通過在目標函數(shù)中加入L2或者L1范數(shù)正則化約束條件，可以更容易獲得稀疏解，是一種特別適用于高維問題的高效算法［15］.由于稀疏正則化的優(yōu)秀性能，其被廣泛應用于圖像的處理［16］以及神經(jīng)網(wǎng)絡的結構優(yōu)化以及對EEG信號的特征篩選.

其中，fil表示對第l類的第i個樣本，wil表示對應于fil的權重信息，yil是fil的響應，cl是第l類的截距，x表示待求解.通過對上式目標函數(shù)的求解，完成對多頻帶CSP特征的篩選得到特征集合S_GL.

1.2.4LDA分類模型LDA模型用于對低維數(shù)據(jù)的分類，LDA是一種有監(jiān)督的方式，通過線性變化最大化類間距離的同時最小化類內距離［17］，通過降維處理和特征選擇方法對特征進行篩選，因此用LDA建立分類模型.

使用LDA分別對上述3種精特征篩選方法提取到的特征集合S_AE，S_GL，S_MI分別進行建立辨識模型.

2 結果與分析

2.1CSP粗特征提取結果以實驗對象一的左右手運動想象為例，挑選出的18個子頻帶和4～42 Hz總頻帶的CSP特征的大腦拓撲圖，如圖5所示.左手/右手運動想象主要分布在C3和C4通道所在的腦域，4～42 Hz和18個頻帶的空域濾波器所提取出的特征有效的區(qū)分出C3和C42個部分，說明了CSP算法在這些子頻帶粗提取出的特征是有效的為下述的信號處理提供了保障.

圖5 2008年BCI競賽數(shù)據(jù)集中的18個子頻帶以及4～42頻帶對應CSP腦部拓撲圖（以左右想象為例）

2.23種精特征篩選結果分析

2.2.1自編碼器篩選結果圖6b～e為二維CSP特征圖.圖6a，f～k為自編碼器壓縮提取后的二維特征圖.圖6f～k所示4個類別兩兩之間的精篩選特征信息呈現(xiàn)出較為明顯的可分性能，并與圖6b～eCSP粗提取出的特征相比，圖6a所示自編碼器處理后的特征更加可分.表明經(jīng)過自編碼器處理后能夠在CSP粗提取特征的基礎上精提取出更加有效的特征信息.

圖6 CSP和AE二維特征圖

2.2.2組稀疏互信息熵篩選結果圖7展示了組稀疏和互信息熵2種特征選擇方法在CSP粗特征提取后精篩選特征的位置分布.2種篩選方法在0～150維所選擇的特征分布相似，但是互信息熵在150～300維未選擇特征，可能丟失部分重要頻段的信息，而組稀疏則選擇了更多重要的特征.

圖7 組稀疏及互信息熵特征篩選圖

2.3預測結果為了驗證本文所設計濾波器組的效果，將其與常規(guī)的濾波組方法進行比較，常規(guī)濾波器組是將4～40 Hz的EEG信號等權重劃分為［4～8 Hz，8～12 Hz，…，36～40 Hz］.表1中A表示本文所提出的4～42 Hz濾波器組方法，B表示4～40 Hz均等濾波器組方法.

表1 分類結果（Kappa值）對比

表1為3種精特征篩選方法分別在2種濾波器組結構下的分類結果.結果表明在相同的精特征篩選方法下，本文所提出的濾波器組方法的結果都高于4～40 Hz均等濾波器組的結果，說明本文所提出的濾波器組設計方法的有效性.秩和檢驗的結果顯示2種方法無統(tǒng)計顯著性，這可能是由于實驗樣本量少以及部分對象結果極值導致的.

圖8分析了3種精特征篩選方法在本文所提的濾波器組的效果.除了在實驗對象8以外，組稀疏的識別效果均是最優(yōu)，這可能是由于本文提出的濾波器組方法提升了數(shù)據(jù)的維度，而組稀疏在處理高維數(shù)據(jù)有著優(yōu)秀的性能.因此，本文以結合組稀疏所建立的模型作為最終的識別方法.

圖8 3種精篩選方法建模結果圖

為了進一步說明本文所提方法性能，將其與競賽前三名以及近年基于該數(shù)據(jù)集的結果進行比較.

表2為本文方法在2008年BCI競賽數(shù)據(jù)集上的分類結果.由表2可以看出，本文基于濾波器組設計方法的模型的分類精度達到0.58，高于競賽前三名的預測Kappa平均值0.57（P=0.97），0.52（P=0.571）和0.31（P=0.005），但分類精度略低于SSDT_KNN（KNN結合特定對象的決策樹模型）的0.61（P=0.85）［18］和ESVL（集成支持向量學習）［19］的0.60（P=0.596）.SSDT_KNN將數(shù)據(jù)投影到目標空間，并訓練4個KNN模型用于構建決策樹模型；ESVL需要計算先驗概率和多個特征，訓練3個SVM分類器組合用以預測.相較于SSDT_KNN和ESVL需要訓練多個分類器，本文設計的方法則不需要.結果表明了本文提出的濾波器組設計方法能夠突出強調與運動想象相關的頻帶信息表達，并且可以在一定程度消除個體差異導致的頻帶波動問題，提供了一種穩(wěn)定的識別方法.

3 小結

面向腦機接口的運動想象問題，筆者分析研究了運動想象的頻帶選擇問題，設計出一種新的濾波器組方法.針對該濾波器組的特性，按照粗特征提取和精特征篩選的方法測試了該濾波器組的效果，最終預測Kappa平均值為0.58，具有優(yōu)秀的分類性能.實驗結果表明，所提出的方法能夠突出強調alpha和beta波段的信息，并且能夠在一定程度上消除由于個體差異導致關鍵頻帶范圍波動的問題.