程昭立，王曉琳，李紅權(quán)，李鴻云，伏云發(fā)，3

1(昆明理工大學(xué) 信息工程與自動化學(xué)院，云南 昆明 650500)

2(昆明理工大學(xué) 腦認知與腦機智能融合創(chuàng)新團隊，云南 昆明 650500)

3(云南省計算機技術(shù)應(yīng)用重點實驗室，云南 昆明 650500)

1 引 言

腦機接口(Brain Computer Interface，BCI)是一種變革性的人機交互，屬于人機交互的子領(lǐng)域，旨在大腦與外設(shè)之間建立一條通路，而不受大腦外周神經(jīng)和肌肉的控制[1，2].該技術(shù)可為運動障礙患者在特定情況下提供可選的、新的通信或控制方式，目前，基于運動想象的BCI是一種新興的可用于運動障礙康復(fù)的技術(shù)[3].

腦卒中是一種腦組織損傷的疾病，造成該疾病的原因是腦部血管突然破裂或因血管阻塞導(dǎo)致血液不能流入大腦.近些年來我國每年新發(fā)腦卒中病例250萬，其中50%～80%的腦卒中患者會遺留有不同程度的肢體運動功能障礙[4]，從而嚴重影響患者的生活自理能力和社會參與能力，給社會和家庭帶來沉重的負擔(dān).已有研究表明，運動想象和實際肢體運動能激活相同的大腦皮層運動區(qū)[5]，由于人的大腦具有極強的可塑性，運動想象將有助于腦卒中偏癱患者大腦受損神經(jīng)的康復(fù)[6].

綜上所述，本研究針對運動協(xié)調(diào)想象，設(shè)計了一種較新的實驗范式，包括步行想象和空閑狀態(tài)，前者為被試想象正常的步行動作，后者為被試保持清醒放松的狀態(tài).HHT是一種適用于非線性非平穩(wěn)信號的有效時頻分析方法[7-9]，考慮到腦電信號的非線性非平穩(wěn)性，已有研究把HHT應(yīng)用到MI-BCI中，且具有較好的效果，但是該方法尚沒有在復(fù)合肢體參與的MI中驗證其有效性.AR模型已用于MI-BCI中，該模型是一種時間序列分析方法，模型參數(shù)凝聚了重要信息，同時研究表明，AR模型的自回歸參數(shù)對狀態(tài)變化反映很敏感[10-12].因此，AR模型的自回歸參數(shù)作為特征向量來表征被試狀態(tài)相關(guān)的腦電變化具有較好的效果[13，14].本研究擬采用HHT提取平均瞬時能量特征，并與AR模型系數(shù)結(jié)合構(gòu)造特征向量，最后采用支持向量機(Support Vector Machine，SVM)對步行想象和空閑狀態(tài)兩類任務(wù)進行分類，以驗證其有效性.

2 材料及方法

2.1 被試

本研究招募了18名被試參與本實驗，均為男性，本科以上學(xué)歷，右利手，所有被試身體健康且無任何心理疾病，年齡在22～28歲之間，此前均無參與過相關(guān)實驗.由于不同被試進行運動想象時激活腦區(qū)的強度有較大差異，故在實驗前，依據(jù)運動覺-視覺想象問卷(The Kinesthetic and Visual Imagery Questionnaire，KVIQ)[15]對每名被試的運動想象能力進行評估.本實驗獲得昆明理工大學(xué)醫(yī)學(xué)院醫(yī)學(xué)倫理委員會的批準并遵照《赫爾辛基宣言》，每名被試在實驗前都自愿簽署了知情同意書.

2.2 實驗范式

運動協(xié)調(diào)想象至少包含兩種肢體參與協(xié)作，是復(fù)合肢體運動想象，其運動感覺的心理排演和傳統(tǒng)單側(cè)肢體簡單運動想象的心理排演是不一樣的.傳統(tǒng)單側(cè)肢體簡單運動想象可能激活有限數(shù)目的腦區(qū)，而運動協(xié)調(diào)想象可能需要更多的腦區(qū)參與，協(xié)同相關(guān)的感知覺資源.其中，上肢和手的協(xié)調(diào)運動可以訓(xùn)練動作的精確性、反應(yīng)速度以及動作的節(jié)奏性，而下肢的協(xié)調(diào)運動可以練習(xí)正確的步態(tài)和上下肢動作的配合、協(xié)調(diào)等.為探索腦卒中偏癱患者運動障礙康復(fù)的方法，針對運動協(xié)調(diào)想象，本實驗設(shè)計了一種新的步行想象范式，被試想象正常的步行動作，步行動作涉及人體四肢的協(xié)調(diào)配合，同時又設(shè)計了任務(wù)以識別被試的空閑狀態(tài).

對于步行想象，一次完整試驗的時序如圖1所示，在進行試驗前，先告知被試試驗任務(wù)為進行步行想象以及保持空閑狀態(tài)，步行想象為在大腦中想象正常的行走，空閑狀態(tài)為全身心清醒放松的狀態(tài).為了提高被試對兩種任務(wù)的適應(yīng)度，在第一個trial開始前，被試先觀看正常步行動作視頻，并認真體會肢體感覺，以強化步行動作在大腦中的預(yù)演，以及練習(xí)保持空閑狀態(tài)，每名被試僅在第一個trial開始前進行兩種任務(wù)的預(yù)演.在t=0s時，屏幕上呈現(xiàn)“固定十字”，提示被試準備開始實驗，持續(xù)時間為3s；在t=3s時，屏幕上的固定十字消失，且呈現(xiàn)試驗任務(wù)的提示，提示被試準備進行對應(yīng)試驗任務(wù)，持續(xù)時間為1s；在t=4s時，提示消失，屏幕為空屏，要求被試執(zhí)行試驗任務(wù)，持續(xù)時間為4s；t=8s時，任務(wù)結(jié)束，屏幕呈現(xiàn)星形提示，提示被試進行休息，持續(xù)時間為2s.休息結(jié)束，開始下一個trial.每名被試對兩種任務(wù)分別進行80個trials，兩種任務(wù)則為160個trials.

圖1 實驗范式

2.3 實驗設(shè)置

本研究使用的實驗設(shè)備是北京中科新拓儀器有限責(zé)任公司生產(chǎn)的NT9200系列EEG采集設(shè)備，電極帽為32通道，采樣率設(shè)為1000Hz，帶通濾波設(shè)為0.1-100Hz之間，同時采用50Hz陷波濾波器避免電源線干擾.電極帽的電極位置按照國際標準10-20導(dǎo)聯(lián)系統(tǒng)設(shè)置[16]，參考電極為A1和A2，所有電極的阻抗都調(diào)至10kΩ以下.在本實驗中，采集了FC3、FCz、FC4、C3、Cz、C4、CP3、CPz、CP4、O1、Oz、O2共12個通道的腦電數(shù)據(jù)進行研究.電極布局示意圖如圖2所示.

圖2 電極布局示意圖

2.4 數(shù)據(jù)分析

2.4.1 數(shù)據(jù)預(yù)處理

現(xiàn)有研究表明，運動想象能夠激活相關(guān)大腦感覺皮層運動區(qū)的腦電活動，并且主要在μ(8-13Hz)和β(14-30Hz)節(jié)律段表現(xiàn)最為顯著[17，18].腦電信號是一種很微弱的信號，幅值大概在10-100uV之間，很容易受到其他電極或者環(huán)境造成的干擾，因此需要對腦電信號進行預(yù)處理，從而提高信噪比[19，20].本文首先采用橢圓濾波器對原始腦電數(shù)據(jù)進行8-30Hz的數(shù)字帶通濾波，通帶衰減為0.5dB，阻帶衰減為50dB.然后對濾波后的腦電信號進行基線漂移校正，以消除腦電信號相對于基線的偏離.最后采用獨立成分分析(Independent Component Analysis，ICA)去除眼電、肌電等偽跡.

2.4.2 希爾伯特-黃變換

特征提取是腦機接口系統(tǒng)中的一個關(guān)鍵部分，通過特征提取可以將表征被試進行運動想象任務(wù)的大腦活動特征量提取出來.對于大腦活動中差異顯著的特征量，可以降低模式識別分類器的錯誤分類率，從而有效提高特征分類正確率.

希爾伯特-黃變換(Hilbert-Huang Transform，HHT)由美國工程院院士黃鍔于1998年提出[21]，HHT是一種適用于處理非線性、非平穩(wěn)信號的自適應(yīng)信號處理方法，它能夠很好地描述信號中存在的局部振蕩頻率分量[22].HHT主要由兩個部分組成：1)經(jīng)驗?zāi)B(tài)分解(Empirical Mode Decomposition，EMD)；2)Hilbert譜分析(Hilbert Spectral Analysis，HSA).

由于任意信號都可以看作是由一系列的本征振蕩模態(tài)組成，則經(jīng)驗?zāi)B(tài)分解能夠根據(jù)信號特征自適應(yīng)地將任意復(fù)雜信號分解成一系列本征模函數(shù)(Intrinsic Mode Function，IMF)，且EMD具有自適應(yīng)性、正交性、完備性以及對IMF分量的調(diào)制特性[23].經(jīng)驗?zāi)B(tài)分解需要滿足兩個條件[24]：1)初始信號中極值點的數(shù)量和過零點數(shù)量之差不超過1；2)初始信號上下包絡(luò)關(guān)于時間軸局部對稱.

EMD分解過程如下[25，26]：

1)求出原始數(shù)據(jù)x(t)中所有的極大值點和極小值點，并用Spline插值擬合形成原數(shù)據(jù)的上包絡(luò)線，以及將所有的極小值點通過Spline插值擬合形成數(shù)據(jù)的下包絡(luò)線，上包絡(luò)線和下包絡(luò)線的均值記為h1，將原始數(shù)據(jù)x(t)減去該平均包絡(luò)h1，得到一個新的數(shù)據(jù)序列k1：

x(t)-h1=k1

(1)

如果k1滿足IMF的兩個條件，則k1為x(t)的第1個分量.

2)若k1不滿足IMF的兩個條件，則將k1作為原始數(shù)據(jù)，重復(fù)步驟1)，直到循環(huán)第i次，得到k1(i-1)-h1i=k1i，使得k1i滿足IMF條件.記d1=k1i，則d1為原始數(shù)據(jù)x(t)的第1個IMF分量.

3)從x(t)中分離出d1，得到：

r1=x(t)-d1

(2)

將r1作為原數(shù)據(jù)重復(fù)步驟1)-步驟2)，求得x(t)的第2個IMF分量d2，重復(fù)循環(huán)n次，得到n個IMF分量.

(3)

若rn成為單調(diào)函數(shù)且不能再從中提取IMF分量時，則循環(huán)結(jié)束.由式(2)和式(3)得到：

(4)

因此，對于任何信號x(t)，都可以將其分解為n個IMF分量和一個殘量rn之和.其中，n個IMF分量d1，d2，…，dn分別包含了數(shù)據(jù)信號從高到低不同頻段的成分，且其都是平穩(wěn)的.

經(jīng)EMD分解后，再進行Hilbert譜分析，下式對每個IMF分量進行Hilbert譜變換：

(5)

解析信號為：

zj(t)=dj(t)+iyj(t)=aj(t)eiθj(t)

(6)

進而，可以分別得到瞬時幅值和瞬時相位為：

(7)

(8)

本研究采用HHT變換提取平均瞬時能量再結(jié)合AR模型系數(shù)構(gòu)造特征向量.首先，對每一種任務(wù)下的每個腦電信號進行EMD分解，然后，對每一個IMF分量用Hilbert變換進行譜分析，計算出經(jīng)HHT變換后幅值的平均瞬時能量，再結(jié)合AR模型系數(shù)構(gòu)造特征向量，最后采用SVM對測試集分類.

2.4.3 AR模型

對于一般的隨機信號，其AR模型可以描述如下：

其中，y(i)是信號的第i個采樣值，φi是AR模型的第j個系數(shù)，n(i)是白噪聲殘差，p是AR模型的階數(shù).

本研究首先通過AIC準則對AR模型進行定階，得到的AR模型最佳階次為6.然后采用Burg算法提取6階的AR模型系數(shù)AR1，…，AR6，由C3、Cz和C4這3個通道的AR模型系數(shù)組成18維特征向量：{C3AR1，…，C3AR6，CzAR1，…，CzAR6，C4AR1，…，C4AR6}.

最后把C3、Cz和C4這3個通道由AR模型系數(shù)構(gòu)成的18維特征向量與alpha和beta頻帶能量構(gòu)成的2維特征組合成20維的特征向量，這20維的特征向量由2.4.4節(jié)的支持向量機方法進行分類.整個特征提取過程如圖3所示.

圖3 特征提取過程

2.4.4 支持向量機

支持向量機(Support Vector Machine，SVM)[27]通過尋求結(jié)構(gòu)風(fēng)險最小化來提升算法的泛化能力，使得置信范圍和經(jīng)驗風(fēng)險達到最小，從而能夠在樣本量較少的情況下，也能獲得較高的分類準確率.

SVM算法的關(guān)鍵在于[28]：1)當線性可分時，利用所求的支持向量來尋求一個最大間隔的超平面；2)當線性不可分時，通過核函數(shù)將在低維空間中線性不可分的特征向量映射到高維空間中使得線性可分.本實驗采集的腦電信號正是樣本量較少且在低維空間中線性不可分的，因此將SVM算法應(yīng)用于本實驗.

在本研究中，對C3、Cz、C4這3個通道的腦電數(shù)據(jù)進行分析，實驗數(shù)據(jù)為80個訓(xùn)練樣本和80個測試樣本，以及訓(xùn)練樣本和測試樣本的標簽.首先用訓(xùn)練樣本和標簽訓(xùn)練模型，然后用訓(xùn)練好的模型對測試樣本進行標簽預(yù)測，最后將預(yù)測標簽與真實標簽進行對比，從而得到分類正確率.

3 結(jié) 果

表1為12個電極通道每個頻段的t檢驗值和步行想象的方差分析p值，方差分析因子有3個：第1個是試驗任務(wù)，包括步行想象及空閑狀態(tài)；第2個是腦電信號的頻段范圍(delta，theta，alpha，beta，gamma)；第3個是腦電電極通道(12個通道).顯著性水平設(shè)為p<0.05.兩種視覺想象任務(wù)的t檢驗結(jié)果如表1所示，表中還包括了方差分析得出的每個電極位置的p值和df值.從表中可以看出，最顯著的差異出現(xiàn)在C3、Cz和C4的α頻帶(p<0.001).除此之外，還進行了t檢驗，以確定在哪個通道有顯著性的差異，結(jié)果表明，位置C3、Cz、C4的p值小于0.01，這表明步行想象和空閑狀態(tài)兩種任務(wù)之間存在99%以上的顯著差異.

表1 12個電極通道每個頻段的t檢驗值和方差分析p值

本文計算各個IMF的方差貢獻率與相關(guān)系數(shù)以進行IMF篩選，如表2所示.從表2可以看出，前4階IMF的方差貢獻率之和達到了96.47%，相關(guān)系數(shù)也很高.因此，本文對前4階IMF采用Hibert變換并求取平均能量值作為能量特征.

表2 所有IMF分量的方差貢獻率與相關(guān)系數(shù)

表3呈現(xiàn)了采用HHT和AR模型組合方法提取特征(試驗任務(wù)開始后4s期間)并利用SVM分類器取的平均、最高和最低分類率.

表3 采用HHT和AR模型組合方法提取特征取得的分類正確率

表4呈現(xiàn)了分別采用CSP、AR、HHT和AR模型組合方法提取特征(試驗任務(wù)開始后4s期間：4～8s)并利用SVM取得的平均、最高和最低分類率.

表4 分別采用CSP、AR、HHT和AR模型組合方法提取特征并利用SVM取得的平均、最高和最低分類率

圖4為步行想象腦電信號經(jīng)過EMD后得到的IMF，由圖3可以看出，腦電信號分解成了8個IMF分量.

圖5呈現(xiàn)了分別采用CSP提取特征、AR模型提取特征、HHT和AR模型結(jié)合提取特征所得的分類正確率隨時間變化的曲線.圖4(a)結(jié)果表明，測試集在0～1.5s間的分類曲線的趨勢是向上的，且包含一個極大值點，在1.5～2.5s之間波動幅度比較大，在2.5～4s趨于平穩(wěn)，并維持在80%附近.圖4(b)結(jié)果表明，測試集在整個0～4s期間的分類曲線波動幅度較大，在0～1.5s振蕩向上，包含一個極大值點，1.5～3s之間有較大幅度下降，包含一個極小值點，3～4s之間又有較大幅度上升，并在4s達到最大值.圖4(c)結(jié)果表明，測試集在0～4間的分類曲線整體趨于平穩(wěn)，0～1.2s之間趨于上升，之后便趨于平穩(wěn)，并維持在83%附近.

圖4 步行想象腦電信號EMD后各階IMF

圖5 分類正確率隨時間變化的曲線

4 討 論

迄今為止，大多數(shù)運動想象針對的是單側(cè)肢體的研究，而傳統(tǒng)的單側(cè)肢體簡單的運動想象不易被患者理解并且康復(fù)效果有限.與傳統(tǒng)單側(cè)肢體簡單運動想象相比，運動協(xié)調(diào)想象至少包含兩種肢體參與協(xié)作，是復(fù)合肢體運動想象，其運動感覺的心理排演和傳統(tǒng)單側(cè)肢體簡單運動想象的心理排演是不一樣的.因此，有別于傳統(tǒng)單側(cè)肢體運動想象研究，本研究針對運動協(xié)調(diào)想象，設(shè)計了一種新的實驗范式，包括想象正常的步行以及保持清醒放松的狀態(tài)，易于被試理解.

本論文在實驗范式上選擇步行想象和空閑狀態(tài)作為試驗任務(wù)對，科學(xué)假設(shè)是該兩種任務(wù)誘發(fā)的腦電特征具有可分性.由表1的結(jié)果可以看出，兩種任務(wù)誘發(fā)的腦電信號之間最顯著的差異出現(xiàn)在C3、Cz和C4的α頻帶；由表2的結(jié)果可以看出，前四階IMF分量的方差貢獻率達到了96.47%；由表3的結(jié)果可以看出，采用HHT和AR模型相結(jié)合的方法提取特征能夠較好地區(qū)分(平均分類正確率為80.51±4.06%)步行想象和空閑狀態(tài)兩種任務(wù).因此，本研究設(shè)計的實驗任務(wù)具有可分性.可以預(yù)見，如果能夠找到更好的特征提取和分類方法，還有可能提高分類正確率，這是我們未來要開展的一個工作.

在本文的研究中，共空間模式(CSP)特征提取與自回歸(AR)模型特征提取相比，最高分類正確率提高了5.32個百分點，平均分類正確率提高了10.29個百分點，這表明前者優(yōu)于后者；希爾伯特-黃變換(HHT)與AR模型相結(jié)合的特征提取和CSP特征提取相比，最高分類正確率提高了4.12個百分點，平均分類正確率提高了6.93個百分點，這表明該組合方法優(yōu)于CSP特征提取方法，更優(yōu)于AR模型.比較圖5的結(jié)果，HHT和AR模型相結(jié)合的方法比CSP具有更佳的分類正確率隨時間變化的曲線(其分類曲線整體上更趨于穩(wěn)定且呈上升狀態(tài)).

表5呈現(xiàn)了運動想象研究中實驗范式、特征提取方法、分類方法及分類正確率.本研究對步行想象和空閑狀態(tài)的平均分類正確率比Tariq M等人[29]對左右膝關(guān)節(jié)伸展運動想象的平均分類正確率高10.51個百分點，我們推測可能是由于想象左右膝關(guān)節(jié)伸展具有一定的難度，導(dǎo)致誘發(fā)的腦電信號可分性較差；本研究的平均分類正確率比Hsu W C等人[30]對左右腿踩踏運動想象的平均分類正確率低5.74個百分點，比唐智川等人[31]對左手和腳運動想象的平均分類正確率低7.47個百分點，比任士鑫等人[32]對腳運動想象和空想的平均分類正確率低1.49個百分點，我們推測可能是由于實驗范式有較大的差異，且特征提取和分類方法各不相同.

表5 運動想象研究中實驗范式、特征提取方法、分類方法及分類正確率

本研究進行的是離線實驗任務(wù)以及在健康被試中進行的實驗，我們未來的研究工作是：1)在線驗證和完善所提方法；2)進一步完善實驗范式，可以考慮將步行想象分解為數(shù)個動作；3)將此實驗應(yīng)用到腦卒中患者，以期步行想象能應(yīng)用到腦卒中患者運動功能障礙康復(fù)中.

5 結(jié) 論

有別于傳統(tǒng)的單側(cè)肢體簡單運動想象，運動協(xié)調(diào)想象至少包含兩種肢體參與協(xié)作，是復(fù)合肢體運動想象，其激活的腦區(qū)數(shù)目更多，且更易于被試理解并在大腦中預(yù)演.針對運動協(xié)調(diào)想象，同時為提高試驗任務(wù)對誘發(fā)的腦電特征之間的可分性，本研究提出一種新的且易于被試執(zhí)行的步行想象和空閑狀態(tài)實驗范式.本研究表明該試驗任務(wù)對誘發(fā)的腦電信號之間最顯著的差異出現(xiàn)在C3、Cz和C4這3個通道的α頻帶.通過不同的特征提取方法，并采用SVM對步行想象和空閑狀態(tài)兩種任務(wù)進行分類，發(fā)現(xiàn)HHT和AR模型的組合方法優(yōu)于CSP和AR模型，且平均分類正確率達到80.51±4.06%，具有良好的可分性.本研究可望為基于步行想象的BCI用于腦卒中偏癱患者運動障礙康復(fù)提供思路，同時也可望為構(gòu)建新型的在線下肢運動協(xié)調(diào)想象BCI提供策略.