段 放 高小榕

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

引言

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

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

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

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

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

1 計(jì)算方法

1.1 共空間模式

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

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

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

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

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

1.2 共迭代時(shí)空模式

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

1.2.1 問題提出

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

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

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

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

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

1.2.2 優(yōu)化求解

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

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

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

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

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

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

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

1.2.3 算法的收斂性

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

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

2 實(shí)際腦電信號(hào)評(píng)估

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

為檢驗(yàn)ICSTP的運(yùn)動(dòng)想象特征提取能力，使用人工選擇濾波頻段的CSP算法進(jìn)行特征提取作為比較對(duì)象，評(píng)價(jià)ICSTP算法的性能。

2.1 人工進(jìn)行頻段選擇

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

2.2 特征分類器

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

2.3 腦電實(shí)驗(yàn)數(shù)據(jù)

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

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

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

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

2.5 數(shù)據(jù)實(shí)驗(yàn)結(jié)果

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

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

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

3 討論

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

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

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

除正確率外，算法的計(jì)算時(shí)間也非常重要。ISSPL算法每次進(jìn)行頻域?yàn)V波器計(jì)算需求解一個(gè)二次規(guī)劃問題。此二次規(guī)劃問題的維度為2×頻率窗長 +樣本個(gè)數(shù)。而本算法只需解決一個(gè)廣義特征值問題，矩陣的階次為FIR濾波器長度。消耗時(shí)間較少。同時(shí)ISSPL算法中的正則化項(xiàng)的數(shù)值也需要通過交叉驗(yàn)證來確認(rèn)。對(duì)于迭代算法每次迭代過程的正則化項(xiàng)是否應(yīng)該使用相同的值也仍然是一個(gè)開放的問題。使用一臺(tái)配置為Inter Core2 2.8 GHz CPU、4 GB內(nèi)存的普通個(gè)人電腦考察算法計(jì)算時(shí)間。文獻(xiàn)［5］中所描述CSSSP算法中單次優(yōu)化求解需要約1 min，由于CSSSP需要交叉驗(yàn)證，最終約需要20 min。對(duì)于ISSPL算法，忽略算法的正則化系數(shù)的選取時(shí)間，僅在確定一個(gè)系數(shù)的條件下進(jìn)行ISSPL算法運(yùn)算。5次迭代使用時(shí)間約為80 s。本算法經(jīng)過5次迭代選取特征后在設(shè)計(jì)分類器，總耗時(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分布圖。可以發(fā)現(xiàn)其腦電在13 Hz與14.5 Hz的r2值出現(xiàn)兩個(gè)峰值。在經(jīng)過ICSTP迭代學(xué)習(xí)后知道了出現(xiàn)這種現(xiàn)象的原因。即這個(gè)受試的在進(jìn)行左側(cè)想象和右側(cè)想象時(shí)對(duì)應(yīng)的μ節(jié)律是不同的。比較圖3中(b)與(c)，可以發(fā)現(xiàn)ICSTP方法計(jì)算出的FIR濾波器很好地吻合了r2分布圖的結(jié)果。強(qiáng)化左側(cè)想象的FIR濾波器與強(qiáng)化右側(cè)想象的FIR濾波器分別與r2分布圖中兩個(gè)辨識(shí)度最高的頻段相吻合。

4 結(jié)論

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

(致謝 感謝本實(shí)驗(yàn)室吳畏與閆錚提供的重要建議。感謝本實(shí)驗(yàn)室王毅軍提供了研究數(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.