王健鋒， 黃曉霞

(上海海事大學(xué) 信息工程學(xué)院， 上海 201306)

0 引 言

腦機(jī)接口(Brain-Computer Interface，BCI)的目的是建立一種不依賴外周神經(jīng)和肌肉組織的新型通信系統(tǒng).該系統(tǒng)通過(guò)測(cè)量腦信號(hào)，如腦磁圖(magnetoencephalography, MEG)，分析人類的意圖，并將該意圖轉(zhuǎn)換為外部設(shè)備的控制信號(hào).[1-2]因此，BCI可以給外周神經(jīng)或肌肉組織受損的患者提供一種新的與外部世界交流的方式.BCI分為植入式BCI和非植入式BCI.非植入式BCI具有無(wú)創(chuàng)傷性、易于攜帶、價(jià)格低廉等優(yōu)點(diǎn)，為BCI的實(shí)際應(yīng)用和推廣提供可能，成為BCI研究領(lǐng)域的一個(gè)熱點(diǎn).在非植入式BCI研究中，腦電圖(electroencephalogram, EEG)和MEG是兩種常用的輸入信號(hào).對(duì)EEG信號(hào)的研究已經(jīng)取得很大成就，但是EEG信號(hào)的采集容易受到頭蓋骨的影響，導(dǎo)致其空間分辨率不高，誤差很大，使EEG的應(yīng)用范圍受到限制.MEG不受頭蓋骨的干擾，具有較高的時(shí)空分辨率和更低的失真率，能夠直接反映腦內(nèi)場(chǎng)源的活動(dòng)狀態(tài)，受到不少BCI研究者的青睞.[3-4]

在2008年第4屆國(guó)際BCI競(jìng)賽中，第3組數(shù)據(jù)使用的就是MEG數(shù)據(jù)，許多研究機(jī)構(gòu)參加這次競(jìng)賽，其中HAJIPOUR等獲得該組競(jìng)賽的第一名，他們用統(tǒng)計(jì)特征、自回歸模型和小波變換的方法進(jìn)行特征提取，用有監(jiān)督的遺傳算法進(jìn)行特征降維，用線性判別分析(Linear Discrimination Analysis, LDA)和支持向量機(jī)(Support Vector Machine, SVM)進(jìn)行分類，所得MEG的平均分類準(zhǔn)確率為46.9%.我國(guó)燕山大學(xué)的李靜等獲得第二名，他們用快速傅里葉變換和主成分分析(Principal Component Analysis, PCA)法進(jìn)行特征提取，用Fisher線性判別(Fisher Linear Discrimination, FLD)方法進(jìn)行特征降維，用LDA進(jìn)行分類，所得MEG的平均分類準(zhǔn)確率為25.1%.2011年SABRA等[5]對(duì)該MEG進(jìn)行研究，用統(tǒng)計(jì)特征、功率譜密度和離散小波變的方法進(jìn)行特征提取，分別用LDA和SVM進(jìn)行分類，獲得33.3%的平均分類準(zhǔn)確率.雖然這些特征處理方法能夠獲得不錯(cuò)的結(jié)果，但都沒(méi)有兼顧分類準(zhǔn)確率和處理方法簡(jiǎn)易性，即分類準(zhǔn)確率高時(shí)，處理方法比較復(fù)雜，處理方法簡(jiǎn)單時(shí)，分類準(zhǔn)確率卻不高.

在不明顯降低分類準(zhǔn)確率的情況下，兼顧MEG識(shí)別率和處理方法的簡(jiǎn)易性，本文提出一種基于經(jīng)驗(yàn)?zāi)Ｊ椒纸?Empirical Mode Decomposition, EMD)和Hilbert變換的MEG特征提取和分類方法.這是一種時(shí)頻數(shù)據(jù)的處理方法，非常適用于非線性和非穩(wěn)定性信號(hào)的處理，在軸承故障診斷、肌電信號(hào)分析和基于EEG的BCI研究方面都能取得不錯(cuò)的效果，但目前該方法還沒(méi)有被應(yīng)用到基于MEG的BCI研究中.為驗(yàn)證EMD和Hilbert變換對(duì)基于MEG的BCI的處理性能，使用EMD和Hilbert變換對(duì)MEG進(jìn)行特征提取，然后用PCA法對(duì)提取到的特征向量進(jìn)行特征降維，最后用SVM進(jìn)行分類.

1 數(shù)據(jù)描述

本文選用第4屆國(guó)際BCI競(jìng)賽的第3組數(shù)據(jù).該數(shù)據(jù)記錄兩名健康的右利手實(shí)驗(yàn)者S1和S2在執(zhí)行手腕運(yùn)動(dòng)期間的腦部活動(dòng)信息.采集數(shù)據(jù)時(shí)，兩名實(shí)驗(yàn)者坐在帶有MEG采集設(shè)備的椅子上，保持放松，把肘部放在枕墊上以防止上臂和肩運(yùn)動(dòng)，頭部用小枕頭固定.要求實(shí)驗(yàn)者僅僅使用右利手和手腕從中心位置向與其夾角為90°的4個(gè)方向移動(dòng)一個(gè)操縱桿，運(yùn)動(dòng)幅度是4.5 cm.在每次實(shí)驗(yàn)中，方向由實(shí)驗(yàn)者自己選擇，目標(biāo)位置在菱形的4個(gè)角上，分別表明左、右、遠(yuǎn)離實(shí)驗(yàn)者、朝向?qū)嶒?yàn)者.

該MEG數(shù)據(jù)記錄頭部運(yùn)動(dòng)區(qū)域上的10個(gè)信道(LC21，LC22，LC23，LC31，LC32，LC41，LC42，RC41，ZC01和ZC02)的信號(hào)(包含運(yùn)動(dòng)開(kāi)始前0.4 s至開(kāi)始后0.6 s的數(shù)據(jù))，并通過(guò)頻帶0.5～100 Hz的帶通濾波器以采樣率400 Hz進(jìn)行重新采樣.本實(shí)驗(yàn)所使用的MEG數(shù)據(jù)包括：實(shí)驗(yàn)者S1的160次訓(xùn)練數(shù)據(jù)(每個(gè)目標(biāo)方向分別進(jìn)行40次實(shí)驗(yàn))和74次測(cè)試數(shù)據(jù)；實(shí)驗(yàn)者S2的160次訓(xùn)練數(shù)據(jù)(每個(gè)目標(biāo)方向分別進(jìn)行40次實(shí)驗(yàn))和73次測(cè)試數(shù)據(jù).

2 方法描述

2.1 基于EMD和Hilbert變換的特征提取

2.1.1 EMD

EMD主要思想是從復(fù)雜信號(hào)里分離出有限個(gè)本征模函數(shù)(Intrinsic Mode Function，IMF)，這些IMF具有不同的頻率，并且滿足：①其極值點(diǎn)個(gè)數(shù)與過(guò)零點(diǎn)個(gè)數(shù)相同或相差一個(gè)；②其上下包絡(luò)關(guān)于時(shí)間軸局部對(duì)稱.[6-8]EMD具體分解過(guò)程如下.

(1)對(duì)任一信號(hào)x(t)，首先確定x(t)上所有的極大值點(diǎn)和極小值點(diǎn)，然后用三次樣條函數(shù)將所有極大值點(diǎn)擬合成原數(shù)據(jù)序列的上包絡(luò)線，將所有極小值點(diǎn)擬合成原數(shù)據(jù)序列的下包絡(luò)線.

(2)計(jì)算上下包絡(luò)線的均值曲線m1(t)，將原數(shù)據(jù)序列x(t)減去m1(t)可得到一個(gè)新的數(shù)據(jù)序列h1(t)，即

(1)

(3)判斷h1(t)是否滿足IMF定義.若滿足，信號(hào)x(t)第1個(gè)IMFc1(t)=h1(t)；若不滿足，將h1(t)作為原信號(hào)，重復(fù)上述步驟，執(zhí)行k次，直到h1k(t)滿足IMF的定義為止，即得信號(hào)x(t)第1個(gè)IMFc1(t)=h1k(t).

(4)將第1個(gè)IMFc1(t)從x(t)中分離出來(lái)，得到一個(gè)去掉高頻分量的剩余信號(hào)r1(t)，即

(2)

(5)將r1(t)作為原信號(hào)，重復(fù)上述步驟，可得到其他IMF，直到剩余信號(hào)是一個(gè)單調(diào)函數(shù)或只有一個(gè)極值時(shí)篩選過(guò)程結(jié)束.最終原始信號(hào)x(t)可表示為

(3)

式中：ci(t)代表每個(gè)IMF；rn(t)代表剩余信號(hào).

2.1.2 Hilbert變換

經(jīng)過(guò)EMD后，對(duì)每個(gè)IMF進(jìn)行Hilbert變換，可得到原始信號(hào)的Hilbert譜.[9-10]具體步驟如下.

對(duì)每個(gè)IMF作Hilbert變換得

(4)

(5)

這里省略剩余信號(hào)rn(t)，Re表示取實(shí)部，其中每個(gè)IMF的幅值和相位隨時(shí)間改變.式(5)等號(hào)右部稱為信號(hào)x(t)的Hilbert振幅譜，簡(jiǎn)稱Hilbert譜，記為

(6)

該公式將信號(hào)的時(shí)間-頻率-幅值的三維關(guān)系表達(dá)出來(lái)，其含義是瞬時(shí)振幅在頻率-時(shí)間上的分布.

對(duì)Hilbert譜取平方可得Hilbert能量譜H2(w,t).在此基礎(chǔ)上，定義瞬時(shí)能量密度

(7)

很顯然,瞬時(shí)能量密度EI是一個(gè)關(guān)于時(shí)間的函數(shù)，反映能量在時(shí)間上的波動(dòng).

2.2 基于PCA法的特征降維

PCA法是揭示大樣本、多變量數(shù)據(jù)或樣本間內(nèi)在關(guān)系的一種方法，旨在利用降維的思想，把多指標(biāo)轉(zhuǎn)化為少數(shù)幾個(gè)綜合指標(biāo)，降低觀測(cè)空間的維數(shù)，以獲取最主要的信息.通常把轉(zhuǎn)化生成的綜合指標(biāo)稱為主成分，每個(gè)主成分都是原始特征的線性組合，并且各個(gè)主成分之間互不相關(guān)，這就使得主成分比原始特征具有某些更好的性能.[11-12]

在實(shí)際的降維過(guò)程中，根據(jù)前m個(gè)主成分的累積方差貢獻(xiàn)率決定特征向量的維數(shù).當(dāng)前m個(gè)主成分的累積方差貢獻(xiàn)率大于給定的上限值時(shí)，可取前m個(gè)主成分作為特征向量，達(dá)到降維目的.

2.3 基于SVM的分類

SVM是一種基于統(tǒng)計(jì)學(xué)習(xí)理論的機(jī)器學(xué)習(xí)方法，它首先選擇某種非線性映射算法將輸入向量映射到一個(gè)高維特征空間，在該高維特征空間中利用結(jié)構(gòu)風(fēng)險(xiǎn)最小化原則和間隔的概念，構(gòu)造最優(yōu)線性決策函數(shù)，求解最優(yōu)分類超平面，達(dá)到分類的目的.[13-14]

3 實(shí)驗(yàn)過(guò)程

本文所提出的算法包括預(yù)處理、特征提取、特征降維和分類等4個(gè)步驟，其實(shí)驗(yàn)流程見(jiàn)圖1.

圖1 算法實(shí)驗(yàn)流程

3.1 預(yù)處理過(guò)程

本實(shí)驗(yàn)所采用的MEG數(shù)據(jù)在頻帶0.5～7 Hz范圍內(nèi)包含的運(yùn)動(dòng)方向信息最為明顯[15]，故將這10個(gè)信道的MEG數(shù)據(jù)通過(guò)頻帶為7 Hz的低通濾波器進(jìn)行濾波處理.

3.2 特征提取過(guò)程

使用EMD對(duì)預(yù)處理后的MEG數(shù)據(jù)進(jìn)行分解，得到各信道中各樣本的IMF.其中某信道單個(gè)樣本經(jīng)EMD后得到的結(jié)果見(jiàn)圖2.從該圖可以看出，經(jīng)EMD所得的剩余信號(hào)逐漸趨于平滑，故棄之不用，只選取經(jīng)EMD得到的IMF.

圖2 某信道單個(gè)樣本經(jīng)EMD后的結(jié)果

用基于Hilbert變換的方法計(jì)算每個(gè)樣本的IMF的瞬時(shí)能量密度，并將其組合為MEG的特征向量，用于區(qū)分MEG的運(yùn)動(dòng)方向.將圖2所對(duì)應(yīng)樣本的IMF進(jìn)行Hilbert變換后，計(jì)算其瞬時(shí)能量譜，能量隨時(shí)間的變化情況見(jiàn)圖3.

圖3 樣本IMF的瞬時(shí)能量譜

3.3 特征降維過(guò)程

提取MEG數(shù)據(jù)的特征后，得到單個(gè)MEG樣本單個(gè)通道的瞬時(shí)能量譜為398維，每個(gè)樣本共有10個(gè)通道，若將單個(gè)樣本的10個(gè)通道的瞬時(shí)能量譜組合為該樣本的特征向量，特征向量的維度將高達(dá)3 980.為避免復(fù)雜的計(jì)算且提高分類準(zhǔn)確率，在分類之前對(duì)特征向量進(jìn)行特征降維.在本實(shí)驗(yàn)中，使用PCA法將特征向量的維度降為159，然后根據(jù)主成分的累積方差貢獻(xiàn)率的原則取其前26個(gè)主成分，并將歸一化處理過(guò)的26個(gè)主成分作為MEG的特征向量.

3.4 分類過(guò)程

本實(shí)驗(yàn)使用基于線性核函數(shù)的SVM對(duì)獲得的特征向量進(jìn)行分類.為獲得較好的分類效果，對(duì)SVM的懲罰因子C進(jìn)行尋優(yōu).圖4是S1和S2的分類準(zhǔn)確率隨C的變化情況.從圖4可知，當(dāng)C達(dá)到22時(shí)，S1和S2的平均分類準(zhǔn)確率達(dá)到最高值，為44.3%。

圖4 S1和S2分類準(zhǔn)確率隨C的變化

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

通過(guò)本文提出的算法獲得的分類結(jié)果：S1和S2的分類準(zhǔn)確率以及二者的平均分類準(zhǔn)確率分別為36.5%，52.1%和44.3%.而第4屆國(guó)際BCI競(jìng)賽各個(gè)競(jìng)賽小組所獲得的分類結(jié)果見(jiàn)表1.其中競(jìng)賽第一名分別提取MEG的時(shí)域特征、頻域特征和時(shí)頻特征作為的聯(lián)合特征，競(jìng)賽第二名分別提取MEG信號(hào)的時(shí)域特征和頻域特征作為信號(hào)的聯(lián)合特征.

表1 第4屆國(guó)際BCI競(jìng)賽結(jié)果 %

采用本實(shí)驗(yàn)所使用的算法獲得的測(cè)試集的平均分類準(zhǔn)確率比競(jìng)賽第二名所獲得的平均分類準(zhǔn)確率高19.2%，雖沒(méi)能超過(guò)競(jìng)賽第一名所獲得的分類準(zhǔn)確率，但只與競(jìng)賽第一名所獲得的平均分類準(zhǔn)確率相差2.6個(gè)百分點(diǎn).本實(shí)驗(yàn)特征提取算法使用的是EMD和Hilbert變換，只需提取信號(hào)的時(shí)頻特征，比競(jìng)賽第一名所使用的特征提取算法簡(jiǎn)單且易于實(shí)現(xiàn)；本實(shí)驗(yàn)使用PCA法進(jìn)行特征降維，與競(jìng)賽第一名所使用的遺傳算法相比，時(shí)間復(fù)雜度降低很多；本實(shí)驗(yàn)分類器只選用線性SVM，使算法復(fù)雜度和時(shí)間消耗降低.本算法在奔騰處理器、主頻2.1 GHz和1 G內(nèi)存的Window XP下處理整個(gè)MEG數(shù)據(jù)的運(yùn)行時(shí)間是281 s，而競(jìng)賽第一名所使用的算法在奔騰處理器、主頻3 GHz和1 G內(nèi)存的Windows XP下的運(yùn)行時(shí)間是1 043 s，在運(yùn)行環(huán)境基本相同且計(jì)算機(jī)性能并不占優(yōu)勢(shì)的情況下，本算法的運(yùn)行時(shí)間顯著減少.總體來(lái)說(shuō)，在MEG分類準(zhǔn)確率不明顯降低的情況下，本算法的實(shí)現(xiàn)比較簡(jiǎn)單，且時(shí)間性能優(yōu)于競(jìng)賽第一名所使用的MEG數(shù)據(jù)處理算法.

4 結(jié)束語(yǔ)

實(shí)驗(yàn)結(jié)果表明，本文提出的基于EMD和Hilbert變換的特征提取和分類算法不僅能夠用于MEG數(shù)據(jù)的分類，而且能夠獲得較高的分類準(zhǔn)確率和較優(yōu)的時(shí)間性能，綜合性能較好，具有更高的實(shí)際應(yīng)用價(jià)值.由于本研究還在進(jìn)行中，該算法還有許多地方可以改進(jìn)，下一步的目標(biāo)是進(jìn)一步提高該算法的分類準(zhǔn)確率.

參考文獻(xiàn)：

[1] SARDOUIE S H, SHAMSOLLAHI M B. Selection of efficient features for discrimination of hand movements from MEG using a BCI competition IV data set[J]. Frontiersin Neuroscience, 2012(6): 42-47.

[2] 高上凱. 基于節(jié)律性腦電信號(hào)的腦-機(jī)接口[J]. 生命科學(xué), 2008, 20(5): 722-724.

[3] 吳婷, 顏國(guó)正, 楊幫華. 基于小波包分解的腦電信號(hào)特征提取[J]. 儀器儀表學(xué)報(bào), 2008, 28(12): 2230-2234.

[4] KRISHNA S, VINAY K C, RAJA K B. Efficient MEG signal decoding of direction in wrist movement using curve fitting (EMDC)[C]// 2011 Int Conf Image Inform Proc (ICIIP). IEEE, 2011: 1-6.

[5] SABRA N I, WAHED M A. The use of MEG-based brain computer interface for classification of wrist movements in four different directions[C]// 2011 28th Natl Radio Sci Conf (NRSC). IEEE, 2011: 1-7.

[6] LEI Y, LIN J, HE Z,etal. A review on empirical mode decomposition in fault diagnosis of rotating machinery[J]. Mech Systems & Signal Processing, 2013, 35(1): 108-126.

[7] LI Lin, JI Hongbing. Signal feature extraction based on an improved EMD method[J]. Measurement, 2009, 42(5): 796-803.

[8] MENDEZ O, CORTHOUT J, Van HUFFEL S,etal. Automatic screening of obstructive sleep apnea from the ECG based on empirical mode decomposition and wavelet analysis[J]. Physiol Measurement, 2010, 31(3): 273-289.

[9] HUANG N E, WU M L, QU W,etal. Applications of Hilbert-Huang transform to non-stationary financial time series analysis[J]. Appl Stochastic Models in Business & Industry, 2003, 19(3): 245-268.

[10] PANOULAS K I, HADJILEONTIADIS L J, PANAS S M. Hilbert-Huang spectrum as a new field for the identification of EEG event related de-synchronization for BCI applications[C]// 30th Annu Int Conf Eng Med & Biol Soc (EMBS). IEEE, 2008: 3832-3835.

[11] 張穎. 基于PCA的模糊BP網(wǎng)絡(luò)建模方法在藻類生長(zhǎng)狀態(tài)預(yù)測(cè)中的應(yīng)用[J]. 上海海事大學(xué)學(xué)報(bào), 2011, 32(4): 76-79.

[12] JOHNSTONE I M, LU A Y. On consistency and sparsity for principal components analysis in high dimensions[J]. J Am Stat Assoc, 2009, 104(486): 682-693.

[13] 安博文, 李丹, 龐然. 基于 SVM 分類器的集裝箱箱號(hào)識(shí)別法[J]. 上海海事大學(xué)學(xué)報(bào), 2011, 32(1): 25-29.

[14] 馮超, 賀俊吉, 史立. 基于支持向量機(jī)的轎車車型識(shí)別[J]. 上海海事大學(xué)學(xué)報(bào), 2011, 32(3): 85-89.

[15] WALDERT S, PREISSL H, DEMANDT E,etal. Hand movement direction decoded from MEG and EEG[J]. J Neuroscience, 2008, 28(4): 1000-1008.