張聰聰，常湛源，2，李傳江

(1. 上海師范大學(xué)信息與機(jī)電工程學(xué)院，上海 200234；2. 上海智能教育大數(shù)據(jù)工程技術(shù)研究中心，上海 200234)

1 引言

腦機(jī)接口(brain computer interface，BCI)技術(shù)不依賴于神經(jīng)和肌肉，僅通過(guò)大腦可以為運(yùn)動(dòng)障礙者提供一個(gè)重要的新的向外部世界發(fā)送信息和控制命令的選擇[1-3]，分為侵入式和非侵入式兩種，侵入式腦機(jī)接口通常直接植入到大腦的灰質(zhì)，因而所獲取的神經(jīng)信號(hào)的質(zhì)量比較高但其容易引發(fā)免疫反應(yīng)，非侵入式腦機(jī)接口因其裝置方便佩戴于人體而成為研究重點(diǎn)[4]，在BCI系統(tǒng)的腦電活動(dòng)中，運(yùn)動(dòng)想象不依賴于外部事件，以其反應(yīng)主動(dòng)運(yùn)動(dòng)意圖的優(yōu)勢(shì)成為研究的熱點(diǎn)。

BCI系統(tǒng)的關(guān)鍵技術(shù)是對(duì)想象運(yùn)動(dòng)腦思維活動(dòng)所產(chǎn)生的腦電信號(hào)進(jìn)行綜合分析，提取不同類別腦電信號(hào)的特征，對(duì)特征進(jìn)行分類識(shí)別運(yùn)動(dòng)意圖。用于表示腦電信號(hào)運(yùn)動(dòng)意圖的常見(jiàn)特征是時(shí)域、頻域、空域特征，其中時(shí)域特征有幅值和幅值能量等，頻域特征有功率譜密度(Power Spectral Density，PSD)、自回歸模型系數(shù)等，空域特征有無(wú)監(jiān)督的主成分分析(Principal Component Analysis，PCA)、獨(dú)立成分分析(Independent Component Analysis，ICA)和有監(jiān)督的共空間模式(Common Spatial Patterns，CSP)等方法獲取[5-6]。根據(jù)腦電信號(hào)的非線性、非平穩(wěn)的特性，時(shí)頻域分析的方法應(yīng)用越來(lái)越廣泛，例如希爾伯特-黃變換(Hilbert-Huang Transform，HHT)[7]、小波變換(Wavelet Transform，WT)[8]和小波包變換(Wavelet Packet Transform，WPT)[9]。

近年來(lái)，多特征融合的方法取得較多成果，與單一特征相比，多個(gè)特征結(jié)合的方法取得了更高的分類精度[10-11]。例如，在運(yùn)動(dòng)想象分類領(lǐng)域，文獻(xiàn)[12]將HHT提取的邊際譜和瞬時(shí)能量譜結(jié)合近似熵進(jìn)行多特征研究，所有被試者平均識(shí)別率高于80%；在情感分類領(lǐng)域，文獻(xiàn)[13]結(jié)合6種時(shí)域情感特征進(jìn)行分類研究，并將小波特征和信息熵組成復(fù)合特征，都達(dá)到較高識(shí)別率；在腦卒分類領(lǐng)域，文獻(xiàn)[14]將層次理論、小波包能量和模糊熵相結(jié)合，得到的分類精度遠(yuǎn)高于單一基準(zhǔn)特征。在HHT時(shí)-頻分析中，首先通過(guò)EMD將原始信號(hào)分解為IMF分量，EMD分解的IMF擺脫了傅里葉變換的局限性，與小波變換和小波包變換相比，EMD不需要基函數(shù)，克服了小波基函數(shù)無(wú)自適應(yīng)性的問(wèn)題，對(duì)于一段未知信號(hào)可直接分解。但EMD的缺點(diǎn)是模態(tài)混疊，EEMD更好地解決了模態(tài)混疊問(wèn)題[15]，但EEMD分解在多次試驗(yàn)集成平均后，重構(gòu)的各個(gè)模式分量中仍然含有一定幅值的殘留噪聲。因此，本文引入CEEMDAN方法來(lái)克服EMD分解和EEMD分解效率低和模態(tài)混疊的問(wèn)題，CEEMDAN方法的分解過(guò)程具有較好的完備性，并可以提高信號(hào)重構(gòu)的準(zhǔn)確性[16]。熵度量了一個(gè)系統(tǒng)或一段信息的不確定性，將樣本熵(Sample Entropy，SampEn)、近似熵(Approximate Entropy，ApEn)、模糊熵等非線性動(dòng)力學(xué)特征引入心電、肌電、腦電等信號(hào)的數(shù)據(jù)處理具有很大優(yōu)越性[17]。模糊熵(Fuzzy Entropy，F(xiàn)uzzyEn)通過(guò)模糊函數(shù)得到的更多細(xì)節(jié)也使得其比樣本熵和近似熵更精確地定義了熵，此外，模糊熵的相對(duì)一致性更強(qiáng)，對(duì)數(shù)據(jù)長(zhǎng)度的依賴性更小[18]。因此，本文選用模糊熵作為運(yùn)動(dòng)想象信號(hào)的非線性動(dòng)力學(xué)特征。

群體智能優(yōu)化算法因具有結(jié)構(gòu)簡(jiǎn)單、易于實(shí)現(xiàn)等特點(diǎn)，被廣泛應(yīng)用于復(fù)雜問(wèn)題的求解中?；依撬惴?Grey Wolf Optimizer，GWO)自提出以來(lái)，就因具有良好的性能而引起了眾多學(xué)者的廣泛關(guān)注。在函數(shù)優(yōu)化方面，已經(jīng)證明GWO的收斂速度和求解精度優(yōu)于粒子群優(yōu)化算法(Particle Swarm Optimization，PSO)、引力搜索算法(Gravitational Search Algorithm，GSA)、差分進(jìn)化算法(Differential Evolution，DE)、進(jìn)化規(guī)劃(evolutionary programming，EP)和進(jìn)化策略(Evolution Strategy，ES)[19]。在各種分類器中，支持向量機(jī)(Support Vector Machines，SVM)因其對(duì)非平穩(wěn)數(shù)據(jù)的通用性和魯棒性而被廣泛應(yīng)用于腦機(jī)接口中[20]。

為提高運(yùn)動(dòng)想象分類正確率，為腦電信號(hào)模式識(shí)別研究提供一種新的識(shí)別方法和研究思路，本文引入CEEMDAN方法應(yīng)用于運(yùn)動(dòng)想象腦電信號(hào)特征提取，對(duì)信號(hào)進(jìn)行加時(shí)間窗滑動(dòng)改進(jìn)分解，獲取更多的處理數(shù)據(jù)，增加腦電信號(hào)的利用率，并結(jié)合時(shí)-頻-空域-非線性動(dòng)力學(xué)的多特征組合的方法進(jìn)行處理分析，分類器選用GWO-SVM。結(jié)果表明，多特征融合互補(bǔ)，融合識(shí)別率高于單一特征，加時(shí)間窗的CEEMDAN分解方法可以提高分類正確率。

2 實(shí)驗(yàn)數(shù)據(jù)

本研究所采用的數(shù)據(jù)集為BCI Competition II的Data set III競(jìng)賽數(shù)據(jù)，由格拉茨理工大學(xué)(TU Graz)生物醫(yī)學(xué)工程研究所醫(yī)學(xué)信息學(xué)系提供(http：∥bbci.de/competition/ii/download/)。受試者為一位25歲女性，受試者以放松狀態(tài)坐在椅子上，試驗(yàn)任務(wù)是通過(guò)顯示器隨機(jī)的左右提示來(lái)想象左手和右手運(yùn)動(dòng)獲得反饋數(shù)據(jù)。實(shí)驗(yàn)共有7組，每組40次，所有實(shí)驗(yàn)都在同一天進(jìn)行，中間休息幾分鐘。數(shù)據(jù)給出了280個(gè)9s長(zhǎng)的試驗(yàn)。其中t=2時(shí)，聲刺激提示試驗(yàn)開(kāi)始，顯示器“+”顯示1s；然后t=3時(shí)，箭頭(左或右)顯示為左右運(yùn)動(dòng)提示，同時(shí)要求受試者開(kāi)始想象運(yùn)動(dòng)，即運(yùn)動(dòng)想象時(shí)間為3s～9s。圖1為單次試驗(yàn)方案設(shè)計(jì)。單次試驗(yàn)中考慮到受試者從看到提示到開(kāi)始想象運(yùn)動(dòng)之間會(huì)有反應(yīng)延時(shí)，同時(shí)持續(xù)運(yùn)動(dòng)想象時(shí)間長(zhǎng)會(huì)導(dǎo)致疲勞，容易引入較多噪聲，因此，本文中截取中間4s～8s的信號(hào)作為研究對(duì)象。

試驗(yàn)采用Ag/AgCl電極，通過(guò)C3、Cz、C4三通道獲得反饋數(shù)據(jù)，其中C3和C4電極位于大腦的初級(jí)感覺(jué)皮層運(yùn)動(dòng)功能區(qū)，反應(yīng)受試者在想象左手和右手運(yùn)動(dòng)時(shí)的狀態(tài)變化，因此，本文實(shí)驗(yàn)分析選取C3和C4通道。試驗(yàn)電極位置分布如圖2所示。信號(hào)的采樣頻率為128 Hz，原始信號(hào)經(jīng)過(guò)0.5～30 Hz的帶通濾波器。圖2為試驗(yàn)電極位置。在本文中，為增加試驗(yàn)數(shù)據(jù)處理次數(shù)得到更加合理的結(jié)果，采用k折交叉驗(yàn)證的思想，將280次試驗(yàn)分為4組，每組70次，選取其中的3組作為訓(xùn)練集剩下一組為測(cè)試集，選取4次，第五組試驗(yàn)數(shù)據(jù)為原始140次訓(xùn)練和140次測(cè)試數(shù)據(jù)，最終試驗(yàn)結(jié)果為5組試驗(yàn)結(jié)果的平均值。

圖2 試驗(yàn)電極位置

3 方法

3.1 CEEMDAN方法

EMD方法依據(jù)數(shù)據(jù)自身的時(shí)間尺度特征來(lái)進(jìn)行信號(hào)分解，無(wú)須預(yù)先設(shè)定基函數(shù)。EEMD方法把高斯白噪聲引入要分析的信號(hào)中。由于零均值噪聲的特性，噪音經(jīng)過(guò)多次的平均計(jì)算后會(huì)相互抵消，這樣集成均值的計(jì)算結(jié)果就可以直接視作最終結(jié)果，顯著提高EMD算法的穩(wěn)定性。CEEMDAN在EEMD的基礎(chǔ)上，在分解的每一個(gè)階段自適應(yīng)加入均值為0、方差為1的高斯白噪聲序列，計(jì)算唯一的余量信號(hào)來(lái)得到各個(gè)IMF分量。具體實(shí)現(xiàn)過(guò)程如下：

1) 將原始信號(hào)f(t)與白噪聲ni(t)相結(jié)合，ε為信噪比控制系數(shù)，構(gòu)造噪聲信號(hào)fi(t)，得到N個(gè)噪聲信號(hào)fi(t)=f(t)+εni(t)，i=1，2，…，N；利用EMD分解各fi(t)，得到IMF1的ci1(t)及其殘余項(xiàng)ri(t)，計(jì)算ci1(t)的平均值，得到CEEMDAN的第一個(gè)IMF分量c1(t)：

(1)

2) 從f(t)中減去c1(t)得到剩余項(xiàng)R1(t)：

R1(t)=f(t)-c1(t)

(2)

4) 重復(fù)3)得到

fnj-1(t)=Rj-1(t)+εj-1Ej-1(ni(t))

(3)

(4)

Rj(t)=Rj-1(t)-cj(t)

(5)

5) 直到獲取的Rj(t)不超過(guò)兩個(gè)極值點(diǎn)時(shí)，算法結(jié)束，Rj(t)即為最終的殘余分量R(t)，最終f(t)表示為

(6)

其中M為分解的IMF的數(shù)量，R(t)為殘余分量。

3.2 Hilbert變換

對(duì)IMF分量信號(hào)進(jìn)行Hilbert譜分析是腦電信號(hào)時(shí)-頻域分析的有效方法，對(duì)CEEMDAN分解得到的IMF分量代入式(7)進(jìn)行希爾伯特變換

(7)

得到解析信號(hào)zi(t)

zi(t)=ci(t)+jyi(t)=Ai(t)ejθi(t)

(8)

(9)

進(jìn)一步求得希爾伯特瞬時(shí)能量譜(Instantaneous Energy Spectrum) IES(t)和邊際能量譜(Marginal Energy Spectrum)MES(ω)：

(10)

(11)

3.3 共空間模式特征提取

共空間模式算法的基本原理是利用矩陣的對(duì)角化，找到一組最優(yōu)空間濾波器進(jìn)行投影，使得兩類信號(hào)的方差差異最大化，從而得到具有較高區(qū)分度的特征向量。算法具體流程如下：

1) 選取一次試驗(yàn)?zāi)X電數(shù)據(jù)構(gòu)成N×T的矩陣EL或ER，其中N為通道數(shù)，T為一次試驗(yàn)采樣點(diǎn)個(gè)數(shù)，EL表示左手運(yùn)動(dòng)的腦電數(shù)據(jù)，ER表示右手運(yùn)動(dòng)的腦電數(shù)據(jù)，由此得到協(xié)方差矩陣

(12)

(13)

C=BΛBT

(14)

(15)

其中Λ是特征值矩陣；B是對(duì)應(yīng)的特征向量矩陣。

(16)

4) 對(duì)白化后的矩陣分別做特征值分解

SL=UΛLUTSR=UΛRUT

(17)

其中ΛL和ΛR是特征值矩陣，U是對(duì)應(yīng)的特征向量。可證SL和SR的特征向量相等，ΛL與ΛR之和是單位矩陣，SL的最大特征值對(duì)應(yīng)SR的特征值最小。

5) 由白化矩陣構(gòu)造空間濾波器W

W=UTP

(18)

6) 對(duì)腦電信號(hào)濾波得到Z，進(jìn)而得到特征向量fp

ZN*T=WN*NEN*T

(19)

(20)

其中var(X)是計(jì)算樣本X的方差。

3.4 模糊熵

算法實(shí)現(xiàn)過(guò)程如下：

i=1，2，…，N-m+1

(21)

(22)

j=1，2，…，N-m+1，且j≠i

(23)

|u(j+p-1)-u0(j)|)

(24)

(25)

進(jìn)而定義φm(r)和φm+1(r)

(26)

(27)

4) 因此可以定義時(shí)間序列的模糊熵FuzzyEn(m，r)為：

(28)

對(duì)于有限長(zhǎng)的數(shù)據(jù)集，模糊熵估計(jì)為

FuzzyEn(m，r，N)=lnφm(r)-lnφm+1(r)

(29)

3.5 灰狼算法優(yōu)化的支持向量機(jī)

支持向量機(jī)是在人機(jī)交互中識(shí)別生理模式的常用分類器，其具有多功能性、魯棒性和免費(fèi)專用工具箱的大量可用性等優(yōu)點(diǎn)。SVM的目的是尋找一個(gè)超平面來(lái)對(duì)樣本進(jìn)行分割，分割的原則是間隔最大化，最終轉(zhuǎn)化為一個(gè)凸二次規(guī)劃問(wèn)題來(lái)求解。c為懲罰參數(shù)，是調(diào)節(jié)分類間隔和誤分類點(diǎn)個(gè)數(shù)的系數(shù)，c值增大時(shí)對(duì)誤分類的懲罰增大，c值減小時(shí)對(duì)誤分類的懲罰減小，最佳原則是使分類間隔盡量大，同時(shí)使誤分類點(diǎn)的個(gè)數(shù)盡量小。核函數(shù)g的作用是與線性分類方法相結(jié)合來(lái)解決非線性問(wèn)題，核函數(shù)的選擇直接影響性能。本文中將兩者相結(jié)合，利用灰狼算法選擇SVM的最佳懲罰參數(shù)c和最佳核函數(shù)g，各等級(jí)狼的初始位置設(shè)置和狼群數(shù)量以及最大迭代次數(shù)分別設(shè)置為2，20，100。

4 實(shí)驗(yàn)

4.1 CEEMDAN分解

本文中將原始信號(hào)各通道直接CEEMDAN分解改進(jìn)為設(shè)置時(shí)間窗滑動(dòng)分解，時(shí)間窗設(shè)置為1s，信噪比控制系數(shù)ε為0.1，加入噪聲的次數(shù)N為100次，原始4s的腦電信號(hào)由時(shí)間窗分割為4段分別進(jìn)行CEEMDAN分解，在分解具有較高頻率分辨率的基礎(chǔ)上增加各時(shí)間窗上腦電數(shù)據(jù)的利用率。如圖3所示，為某次試驗(yàn)C3通道和C4通道信號(hào)CEEMDAN分解的前5階IMF分量及其頻譜。可以看出各IMF分量的頻率從高到底排列，后幾階IMF為低頻，而運(yùn)動(dòng)執(zhí)行相關(guān)電位的頻率范圍為8～30HZ，因此，本文試驗(yàn)中直接舍棄低頻偽跡，可以看做是自適應(yīng)高通濾波，使信號(hào)處理更加高效。為得到更準(zhǔn)確的分析信號(hào)，對(duì)保留的IMF分量根據(jù)式(30)計(jì)算分量與原始分解信號(hào)的相關(guān)系數(shù)r(其中：c為IMF分量，x為原始信號(hào)，n=128)，選取相關(guān)系數(shù)最大的分量做后續(xù)處理。

(30)

圖3 C3和C4通道腦電信號(hào)CEEMDAN分解及其頻譜分析

4.2 基于Hilbert變換的瞬時(shí)能量差和邊際能量差特征提取

將選取的IMF分量進(jìn)行Hilbert變換求取瞬時(shí)能量譜和邊際譜。如圖4所示，(a)為某次試驗(yàn)想象左手運(yùn)動(dòng)信號(hào)幅值變化，(a)為某次試驗(yàn)想象右手運(yùn)動(dòng)信號(hào)幅值變化，與想象左手或右手運(yùn)動(dòng)產(chǎn)生的事件相關(guān)去同步(event-related desynchronization，ERD)和事件相關(guān)同步(event related synchronization，ERS)現(xiàn)象相一致。圖5為運(yùn)動(dòng)想象的邊際譜，由圖5可以看出左手和右手運(yùn)動(dòng)邊際能量差異明顯。因此選取C3和C4通道的瞬時(shí)能量差和邊際能量差能更好的表示左手和右手運(yùn)動(dòng)想象特征。結(jié)合瞬時(shí)能量差和邊際能量差得到時(shí)-頻組合特征向量F1。

4.3 共空間模式(CSP)的特征提取

將CEEMDAN按時(shí)間窗滑動(dòng)分解的IMF分量篩選和計(jì)算相關(guān)系數(shù)之后，保留每個(gè)時(shí)間窗分解后相關(guān)系數(shù)最大的IMF組合成新的多通道信號(hào)，記為X(IMFC31，IMFC32，IMFC33，IMFC34，IMFC41，IMFC42，IMFC43，IMFC44)。將組合而成的信號(hào)X經(jīng)空間濾波器(spatial filters)W映射處理，得到空域特征記為F2。

4.4 模糊熵(FuzzyEn)特征提取

在想象左手和右手運(yùn)動(dòng)想象時(shí)，模糊熵衡量了C3和C4通道的信號(hào)序列復(fù)雜度。如圖6所示，為100個(gè)樣本中想象左手運(yùn)動(dòng)(a)和想象右手運(yùn)動(dòng)(b)的模糊熵值。本文為得到模糊熵特征，對(duì)每個(gè)通道信號(hào)計(jì)算模糊熵，嵌入維數(shù)m設(shè)置為2，將信號(hào)標(biāo)準(zhǔn)差歸一化后的相似容限值r設(shè)置為0.15，計(jì)算結(jié)果組合成特征集F3。

圖4 想象左手和右手運(yùn)動(dòng)C3和C4通道時(shí)域振幅對(duì)比

圖5 想象左手和右手運(yùn)動(dòng)C3和C4通道邊際譜對(duì)比

圖6 想象左手和右手運(yùn)動(dòng)C3和C4通道模糊熵值對(duì)比

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

將所有特征集組合成18維特征F={F1，F(xiàn)2，F(xiàn)3}，輸入到GWO-SVM分類器進(jìn)行分類，各單一特征和組合特征的分類結(jié)果見(jiàn)表1。

表1 分類結(jié)果對(duì)比

由平均分類準(zhǔn)確率可以看出組合特征結(jié)果明顯高于單一特征，加時(shí)間窗改進(jìn)CEEMDAN分解后的正確率也高于未改進(jìn)。融合特征的準(zhǔn)確率更高，驗(yàn)證了多種方法的融合特征可以對(duì)單一特征彌補(bǔ)不足，是提高分類正確率的重要思路。多種特征融合的分類準(zhǔn)確率中，加時(shí)間窗改進(jìn)處理比未改進(jìn)提高了2.142%；在所有與CEEMDAN分解相關(guān)的特征分類中，平均分類準(zhǔn)確率平均提高了1.673%；其中，特征F1和特征F2分類正確率提高最高，分別為3.002%和2.426%，說(shuō)明加時(shí)間窗分解的方法在一定情況下相當(dāng)于將兩通道數(shù)據(jù)增加為多通道，可以進(jìn)一步提高分類正確率。為描述多次試驗(yàn)結(jié)果的波動(dòng)，選用和數(shù)據(jù)單位一致的標(biāo)準(zhǔn)差來(lái)衡量。對(duì)于改進(jìn)后的單一特征和兩兩特征融合，雖然提高了分類正確率，但是結(jié)果波動(dòng)相比更大，而改進(jìn)后的多特征融合的方法在提高正確率的基礎(chǔ)上標(biāo)準(zhǔn)差幾乎沒(méi)有變化甚至略微減小，結(jié)果具有較強(qiáng)的穩(wěn)定性。因此，加時(shí)間窗改進(jìn)后的多特征融合方法對(duì)運(yùn)動(dòng)想象腦電信號(hào)分類正確率的提高具有重要意義。

5 結(jié)論

本文引入CEEMDAN方法應(yīng)用于運(yùn)動(dòng)想象腦電信號(hào)特征提取，對(duì)BCI Competition II數(shù)據(jù)集進(jìn)行處理分析，在對(duì)信號(hào)CEEMDAN分解時(shí)加時(shí)間窗滑動(dòng)改進(jìn)處理，并結(jié)合時(shí)-頻-空域-非線性動(dòng)力學(xué)模糊熵的多特征融合的方法進(jìn)行組合分析，融合特征的分類正確率達(dá)89.290%，改進(jìn)后的方法對(duì)正確率提高了2.142%，結(jié)果表明改進(jìn)的CEEMDAN分解方法可以提高分類正確率而且具有更高的穩(wěn)定性。因此，多特征融合的方法將各單個(gè)方法的缺點(diǎn)互補(bǔ)，進(jìn)一步提高分類正確率，為腦電信號(hào)模式識(shí)別研究提供了一種新的識(shí)別方法和研究思路。下一步的研究將引入其它新的機(jī)器學(xué)習(xí)方法應(yīng)用于腦電信號(hào)處理，將繼續(xù)研究多特征融合來(lái)進(jìn)一步提高腦電信號(hào)分類正確率。