一種基于運動想象的腦-機接口時空濾波器迭代算法

段 放 高小榕

(清華大學(xué)醫(yī)學(xué)院生物醫(yī)學(xué)工程系，北京 100084)

引言

腦-機接口(brain-computer interface，BCI)是一種不經(jīng)過外圍肌肉運動，通過直接提取大腦活動特征來指揮外部設(shè)備的通信方式［1］。在腦-機接口系統(tǒng)中所使用的監(jiān)測大腦活動的腦功能成像方法多為頭皮腦電圖(electroencephalogram，EEG)信號，其具有對使用者無損傷的優(yōu)勢。

運動想象是在肢體不實際運動的情況下對肢體某個部分的潛在想象［2］。基于運動想象的腦-機接口系統(tǒng)的生理依據(jù)是使用者在想象不同肢體運動時，自主控制其大腦運動感覺皮層與所想象肢體對應(yīng)區(qū)域μ節(jié)律或β節(jié)律的變化。如何提取μ節(jié)律或β節(jié)律的能量差異這種特征是運動想象腦-機接口系統(tǒng)的關(guān)鍵問題。

共空間模式(common spatial patterns，CSP)已被證實是運動想象腦-機接口系統(tǒng)特征提取的有效方法［3］，其生成的空域濾波器可有效提高信號的信噪比。但CSP只能在信號進行了比較合理的時域濾波條件下，才能有效地發(fā)揮其性能［3］。因此通常在使用CSP前對信號的時域帶通濾波變得十分關(guān)鍵。一般情況下，常采用性能不佳的較寬通帶濾波(如7～30 Hz)，或依據(jù)受試腦電信號特點對濾波器人為離線設(shè)計，這必然限制了系統(tǒng)的性能提升。

為了避免過多離線人工干預(yù)，很多機器學(xué)習(xí)方法被引入到CSP的改進中來。例如共空頻模式(common spatial spectral pattern，CSSP)、稀疏空頻模式(common sparse spectral spatial pattern，CSSSP)、濾波器庫空域模式(filter bank common spatial pattern，F(xiàn)BCSP)等［4－8］。這些方法都通過一定程度的遍歷搜索確定出腦電信號的有效頻率帶寬，但由于涉及遍歷搜索計算量均較大或優(yōu)化求解不確定性高等問題，很難在線實現(xiàn)。

本研究提出一種分別對空域濾波器和時域濾波器進行迭代優(yōu)化的特征提取方法，避免了遍歷搜索的時間不確定性，并給出了收斂的必然性證明。減少了特征提取中的人工干預(yù)，有利于在線應(yīng)用。

1 計算方法

1.1 共空間模式

對于兩類多導(dǎo)信號xi(t)∈RC，C為導(dǎo)聯(lián)數(shù)，i=1，2。信號的協(xié)方差矩陣為:

l為兩類信號的長度。CSP生成一組空間加權(quán)系數(shù)w∈RC，使經(jīng)過空域加權(quán)濾波后的信號yi(t)=wTxi(t)更具有辨識力。其信號能量為:

所謂辨識力體現(xiàn)在分別最大化兩者能量比E1/E2或E2/E1。不失一般性，設(shè)為:

這個優(yōu)化問題可通過求解廣義特征值問題Σ1w=λΣ2w獲得結(jié)果，其中λ即為兩信號能量比。限于篇幅，具體的推導(dǎo)過程可參考文獻［3］。

由于通常CSP只在經(jīng)過特定頻率濾波的情況下才能發(fā)揮其有效性，需分別引入一個給定長度的有限沖激響應(yīng)濾波器(Finite Impulse Response filter，F(xiàn)IR filter)H。確定一個有效的H是很關(guān)鍵的，Dornhege等人也提出了使用 FIR濾波器與 CSP結(jié)合最大化能量比的思路(CSSSP)。但他沒有給出FIR濾波器的求解方法，而采用線性梯度搜索的方法通過窮舉的方式確定FIR濾波器。即使對長度較短的 FIR 濾波器也無疑計算量是巨大的［5，7］，因此不適用于在線系統(tǒng)。而較短的濾波器也無法提供相對窄的濾波頻帶，CSP在較寬的濾波頻帶下性能是有限的，因此CSSSP并不能得到足夠好的效果。

1.2 共迭代時空模式

給出一種分別整合空域濾波器、FIR濾波器的機器學(xué)習(xí)方法來提高CSP的性能。原CSP問題為兩個形式相同的最大化能量比問題。

1.2.1 問題提出

圖8為母線電流響應(yīng)法的實驗波形。首先,對開關(guān)組合S1和S2,S5和S6分別進行開通和關(guān)斷,測得的電流響應(yīng)分別為i3和i2,如圖8(a)所示;然后,對開關(guān)組合S3和S4,S5和S6分別進行開通和關(guān)斷,測得的電流響應(yīng)分別是i1和i2,如圖8(b)所示。

對于給定長度的 FIR濾波器h(n)，n=0，…，L－ 1，向前補0l－ L項得到 H=h(n)，n=0，…，l－1，H ∈ Rl。則其濾波器循環(huán)行列式為 A，A ∈ Rl×l，其中:

經(jīng)過 FIR濾波器濾波后的信號 zi=xiA。原CSP問題為最大化能量比問題表述為:

將濾波前的信號帶入式(1)得到:

式(2)為ICSTP所想解決的優(yōu)化問題。

1.2.2 優(yōu)化求解

在機器學(xué)習(xí)中最優(yōu)解往往會引起過擬合的不穩(wěn)定問題，所以可以不追求全局最優(yōu)解。下面給出利用基于迭代的局部最優(yōu)解求解方法。

首先嘗試在某一個初始FIR濾波器H的基礎(chǔ)上優(yōu)化空域濾波器，式(2)轉(zhuǎn)化為

這是一個原始的CSP問題，可通過解廣義特征值問題求得空域濾波器 w。對于原始信號 xi(t)，經(jīng)過空域濾波后的信號為yi(t)=wxi(t)，代入式(2)得到:

得 yiA=HRi，其中 Ri∈ Rl×l為 yi的逆向循環(huán)矩陣。

對于長度為 L的濾波器，H中前 l－L項為0。現(xiàn)將 H長度還原為，只需考慮 Ri的后行，即 Qi∈RL×l，式(2)優(yōu)化轉(zhuǎn)化為

可以發(fā)現(xiàn)經(jīng)時空濾波器互換以上優(yōu)化問題與傳統(tǒng)CSP優(yōu)化問題形式相同，因此次優(yōu)化可使用優(yōu)化傳統(tǒng)CSP方法的求解廣義特征值問題的方法求得新的濾波器H。將新的H變換為A代入式(3)再求解空域濾波器w。

通過式(3)與式(4)交替迭代求得最終需要得到的FIR濾波器H與空域濾波器w。出于腦電信號的特點，我們可以使用任意7～30 Hz FIR濾波器作為初始濾波器。迭代次數(shù)的選擇可以通過以下兩點由使用者直接確定:即到達預(yù)先設(shè)置的迭代次數(shù)或相鄰兩次迭代的空域濾波器w的變化很小。對于運動想象問題有兩個能量比最大化問題E1/E2或E2/E1，分別對應(yīng)于增強左手或右手想象信號能量。因此我們最終根據(jù)兩個問題分別返回1組對應(yīng)的空域濾波器與 FIR 濾波器解 {wi，Hi}，i=1，2。整體算法流程見圖1。

1.2.3 算法的收斂性

作為迭代算法，其收斂性是非常重要的。下面給出這種迭代求解方法收斂性的證明。

每次對式(3)的優(yōu)化是在式(4)的優(yōu)化所得結(jié)果H基礎(chǔ)上進行的，而對式(4)的優(yōu)化則是在式(3)的優(yōu)化所得結(jié)果w基礎(chǔ)上進行的。每次優(yōu)化均極大化兩信號能量比，易知對于第i次迭代有

2 實際腦電信號評估

圖1 算法流程Fig.1 Diagram of algorithm

為檢驗ICSTP的運動想象特征提取能力，使用人工選擇濾波頻段的CSP算法進行特征提取作為比較對象，評價ICSTP算法的性能。

2.1 人工進行頻段選擇

CSP算法的目標(biāo)是使數(shù)據(jù)更有辨識力，因此要選擇有足夠辨識力的頻段作為CSP方法的預(yù)先濾波頻段。利用r2值來評價兩類數(shù)據(jù)在頻域上是否有辨識力［1］，使用r2較高的頻段作為 CSP的預(yù)先濾波頻段。

2.2 特征分類器

模式識別中特征提取與特征分類是兩個連貫的過程。本算法作為一個特征提取的方法，仍需要在特征提取后應(yīng)用相對應(yīng)的特征分類方法。對運動想象的腦電數(shù)據(jù)分類中線性分類器具有足夠好的性能［10］，因此使用線性軟間隔支持向量機(Support Vector Machine，SVM)。SVM 分類器源于統(tǒng)計學(xué)習(xí)原理，由于其基于結(jié)構(gòu)風(fēng)險最小化的特點使其有相對強的抗過擬合能力［11］，這樣較有利于小樣本學(xué)習(xí)問題。之所以使用線性核的SVM，是為了防止樣本學(xué)習(xí)在高維非線性空間產(chǎn)生的過擬合現(xiàn)象。

2.3 腦電實驗數(shù)據(jù)

圖2 數(shù)據(jù)處理流程Fig.2 Diagram of data processing

用以檢驗算法的數(shù)據(jù)是本實驗室采集的運動想象腦電數(shù)據(jù)。7組腦電數(shù)據(jù)分別從年齡在21～24歲之間的7位志愿者處采集，數(shù)據(jù)是在有視覺反饋提示的情況下獲得的。每次試驗時間長度為8 s。在2 s的無提示準(zhǔn)備時間后屏幕上出現(xiàn)方向提示，提示時間1 s。之后受試通過左手或右手的運動想象控制屏幕上光標(biāo)根據(jù)提示進行垂直運動。運動想象時間持續(xù)5 s。5 s想象后根據(jù)光標(biāo)，最終位置給出正確與否的提示。每個受試分兩天共記錄240次試驗數(shù)據(jù)，其中左手想象右手想象各120次用于離線分析。整個實驗采用 BioSemi ActiveTwo腦電放大器系統(tǒng)記錄數(shù)據(jù)。根據(jù)國際10/20系統(tǒng)，在第一運動區(qū)與輔助運動區(qū)設(shè)置了32個電極。參考電極為左耳垂，采樣頻率256 Hz。用以離線分析的數(shù)據(jù)的時間窗選定為試驗的3.5～7.5 s。為CSP算法與ICSTP算法的公平性使用全部32導(dǎo)數(shù)據(jù)用以檢驗算法而不再各自單獨挑選導(dǎo)聯(lián)。更具體實驗時序圖和腦電電極安放的實際位置圖可參考文獻［7］。

2.4 ICSTP設(shè)置與迭代條件

進行一次時空優(yōu)化作為一次迭代。在整個數(shù)據(jù)實驗過程設(shè)定10次迭代為最終迭代次數(shù)。同時為檢驗收斂速度，對于實驗所用32導(dǎo)數(shù)據(jù)，設(shè)定基于空域濾波器的迭代停止條件運動想象人為設(shè)計濾波器帶寬通常小于5 Hz［12］。為獲得小于5 Hz的通帶寬度，設(shè)定FIR濾波器約為1/5采樣頻率，長度為50點。

2.5 數(shù)據(jù)實驗結(jié)果

對數(shù)據(jù)進行了10×3隨機分組劃分實驗。即將兩類數(shù)據(jù)平均分為相等的3份，每份中兩類數(shù)據(jù)數(shù)量相等。每次測試將其中兩份作為訓(xùn)練集，最后一份作為測試集。同時給出寬帶濾波(7～35 Hz)和人工挑選頻段的后CSP結(jié)果，說明頻段選擇的重要性。作為一種相似的迭代算法，給出文獻［7］中對相同數(shù)據(jù)使用ISSPL算法的實驗結(jié)果作比較。分類正確率與迭代次數(shù)見表1，數(shù)據(jù)處理流程見圖2。

表1中黑體字部分為正確率最高的情況?？梢钥吹絻H有兩名受試在人工選擇濾波頻段后進行CSP的結(jié)果好于ICSTP。而CSP優(yōu)于ICSTP的情況是在受試運動想象能力較強的情況下，且兩者結(jié)果差距很小。而對于兩名想象運動能力相對較弱的受試，ICSTP的結(jié)果明顯優(yōu)于 CSP較多。同時本算法與ISSPL相比仍保持向?qū)Ω叩恼_率。

表1 分類正確率±方差與迭代次數(shù)±方差Tab.1 Classification accuracies±standard deviation and iteration times±standard deviation

3 討論

可以發(fā)現(xiàn)平均迭代次數(shù)均不超過4次。所有的數(shù)據(jù)實驗過程均未達到預(yù)先設(shè)定的迭代10次的迭代次數(shù)最大值。在很少的迭代次數(shù)之后空域濾波器已無明顯變化，算法收斂速度快。

3.1 與ISSPL相比較的優(yōu)勢

同樣作為迭代算法，本算法在正確率上也普遍高于ISSPL。這可能是由于本算法在收斂性上能夠得到保證。ISSPL中的將分類器設(shè)計與頻率窗的學(xué)習(xí)相結(jié)合，采用支持向量機計算頻率窗函數(shù)同時確定決策函數(shù)分類面，其空域濾波器的設(shè)計采用CSP最大能量比的思想，頻域濾波器期望取得盡量大的間隔。迭代算法的兩個迭代步驟沒有統(tǒng)一的優(yōu)化目標(biāo)，這樣算法的收斂性難以證明。本算法在兩個步驟均最大化能量比，已證明了算法收斂。本算法與ISSPL還有一個不同點是本算法首先進行特征提取，將特征維度壓縮到2維后再進行分類器學(xué)習(xí)。而ISSPL在學(xué)習(xí)頻率窗函數(shù)的同時進行分類函數(shù)學(xué)習(xí)。分類決策函數(shù)所需要確定的系數(shù)與頻率窗長度相同，維數(shù)高，相比本算法可能容易出現(xiàn)過擬合。

除正確率外，算法的計算時間也非常重要。ISSPL算法每次進行頻域濾波器計算需求解一個二次規(guī)劃問題。此二次規(guī)劃問題的維度為2×頻率窗長 +樣本個數(shù)。而本算法只需解決一個廣義特征值問題，矩陣的階次為FIR濾波器長度。消耗時間較少。同時ISSPL算法中的正則化項的數(shù)值也需要通過交叉驗證來確認(rèn)。對于迭代算法每次迭代過程的正則化項是否應(yīng)該使用相同的值也仍然是一個開放的問題。使用一臺配置為Inter Core2 2.8 GHz CPU、4 GB內(nèi)存的普通個人電腦考察算法計算時間。文獻［5］中所描述CSSSP算法中單次優(yōu)化求解需要約1 min，由于CSSSP需要交叉驗證，最終約需要20 min。對于ISSPL算法，忽略算法的正則化系數(shù)的選取時間，僅在確定一個系數(shù)的條件下進行ISSPL算法運算。5次迭代使用時間約為80 s。本算法經(jīng)過5次迭代選取特征后在設(shè)計分類器，總耗時約15 s?？梢姳舅惴ǜ欣谠诰€系統(tǒng)。

3.2 特征提取效果

圖3 頻域特征。(a)各導(dǎo)聯(lián)頻域 分布圖;(b)左手想象特征濾波器的幅頻響應(yīng);(c)右手想象特征濾波器幅頻響應(yīng)Fig.3 Feature in frequency domain.(a)r2-value of all channels in frequency domain;(b)magnitude response of left-hand movement imagery filter;(c)magnitude response of right-hand movement imagery filter

觀察圖3(a)為受試WW r2分布圖?？梢园l(fā)現(xiàn)其腦電在13 Hz與14.5 Hz的r2值出現(xiàn)兩個峰值。在經(jīng)過ICSTP迭代學(xué)習(xí)后知道了出現(xiàn)這種現(xiàn)象的原因。即這個受試的在進行左側(cè)想象和右側(cè)想象時對應(yīng)的μ節(jié)律是不同的。比較圖3中(b)與(c)，可以發(fā)現(xiàn)ICSTP方法計算出的FIR濾波器很好地吻合了r2分布圖的結(jié)果。強化左側(cè)想象的FIR濾波器與強化右側(cè)想象的FIR濾波器分別與r2分布圖中兩個辨識度最高的頻段相吻合。

4 結(jié)論

本研究提出了一種基于迭代優(yōu)化的運動想象特征提取方法。并且數(shù)據(jù)實驗證實了算法的有效性。這種算法能在無人工干預(yù)的情況下，快速計算出適用于不同受試的運動想象時域濾波器與空域濾波器。作為一種迭代算法，并給出了迭代收斂性的證明。并在數(shù)據(jù)實驗中證實了其收斂的快速性。與傳統(tǒng)的運動想象空頻濾波器特征提取方法相比有著計算速度快的優(yōu)勢。這無疑有利于在線系統(tǒng)的實現(xiàn)，特別是應(yīng)用于在線系統(tǒng)的在線更新問題上。

(致謝 感謝本實驗室吳畏與閆錚提供的重要建議。感謝本實驗室王毅軍提供了研究數(shù)據(jù))

［1］Wolpaw JR， BirbaumerN， McFarlandDJ， etal.Braincomputer interfaces for communication and control［J］.Clinical Neurophysiology，2002，113(6):767－791.

［2］De Vries S，Mulder T.Motor imagery and stroke rehabilitation:a critical discussion ［J］.Journal of Rehabilitation Medicine，2007，39(1):5 －13.

［3］Blankertz B，Tomioka R，Lemm R，et al.Optimizing spatial filters for robust EEG single-trial analysis［J］.IEEE Signal Proc Magzine，2008，25(1):41－56.

［4］Lemm S，Blankertz B，Curio G，et al.Spatio-spectral filters for improved classification of single trial EEG ［J］.IEEE Trans Biomed Eng，2005，52(9):1541－1548.

［5］Dornhege G，Blankertz B，Krauledat M，et al.Combined optimization of spatial and temporal filters for improving brain–computer interfacing［J］.IEEE Trans Biomed Eng，2006，53(11):2274－2281.

［6］Novi Q，Guan C，Dat TH，et al.Sub-band Common Spatial Pattern (SBCSP)forBrain-ComputerInterface［C］//Proceedings of 3rd International IEEE/EMBS Conference on Neural Engineering.Hawaii:IEEE，2007:296－300.

［7］Wu Wei，Gao Xiaorong，Hong Bo，et al.Classifying Single-Trial EEG During Motor Imagery by Iterative Spatio-Spectral Patterns Learning(ISSPL)［J］.IEEE Trans Biomed Eng，2008，55(6):1733－1743.

［8］Ang KK，Chin ZY，and Guan C.Filter Bank Common Spatial Pattern(FBCSP)in Brain-Computer Interface［C］//Proceedings of IJCNN 2008［C］.HongKong:IEEE，2008:2391－2398.

［9］Fukunaga F and Koontz W.Applications of the Karhunen-Loève expansion to feature selection and ordering［J］.IEEE Trans Computers，1970，19(5):311－318.

［10］Dornhege G，Blankertz B，Curio G，et al.Boosting bit rates in noninvasive EEG single-trial classifications by feature combination and multiclass paradigms［J］.IEEE Trans Biomed Eng，2004，51(6):993－1002.

［11］Vapnik VN.Statistical Learning Theory［M］.New York:Wiley，1998.

［12］Pfurtscheller G，Lopes Da Silva FH.Event-related EEG/MEG synchronization and desynchronization:basic principles［J］.Clin Neurophysiol，1999，110(10):1842－1857.