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

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

1 引言

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

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

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

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

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

2 實驗數(shù)據(jù)

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

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

圖2 試驗電極位置

3 方法

3.1 CEEMDAN方法

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

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

(1)

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

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

(2)

4) 重復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)不超過兩個極值點時，算法結束，Rj(t)即為最終的殘余分量R(t)，最終f(t)表示為

(6)

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

3.2 Hilbert變換

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

(7)

得到解析信號zi(t)

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

(8)

(9)

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

(10)

(11)

3.3 共空間模式特征提取

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

1) 選取一次試驗腦電數(shù)據(jù)構成N×T的矩陣EL或ER，其中N為通道數(shù)，T為一次試驗采樣點個數(shù)，EL表示左手運動的腦電數(shù)據(jù)，ER表示右手運動的腦電數(shù)據(jù)，由此得到協(xié)方差矩陣

(12)

(13)

C=BΛBT

(14)

(15)

其中Λ是特征值矩陣；B是對應的特征向量矩陣。

(16)

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

SL=UΛLUTSR=UΛRUT

(17)

其中ΛL和ΛR是特征值矩陣，U是對應的特征向量?？勺CSL和SR的特征向量相等，ΛL與ΛR之和是單位矩陣，SL的最大特征值對應SR的特征值最小。

5) 由白化矩陣構造空間濾波器W

W=UTP

(18)

6) 對腦電信號濾波得到Z，進而得到特征向量fp

ZN*T=WN*NEN*T

(19)

(20)

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

3.4 模糊熵

算法實現(xiàn)過程如下：

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)

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

(26)

(27)

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

(28)

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

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

(29)

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

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

4 實驗

4.1 CEEMDAN分解

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

(30)

圖3 C3和C4通道腦電信號CEEMDAN分解及其頻譜分析

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

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

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

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

4.4 模糊熵(FuzzyEn)特征提取

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

圖4 想象左手和右手運動C3和C4通道時域振幅對比

圖5 想象左手和右手運動C3和C4通道邊際譜對比

圖6 想象左手和右手運動C3和C4通道模糊熵值對比

4.5 實驗結果與分析

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

表1 分類結果對比

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

5 結論

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