基于HCSP和模糊熵的腦電信號分類

于沐涵，陳 峰

(南通大學(xué) 電氣工程學(xué)院，江蘇 南通 226019)

0 引 言

腦-機(jī)接口(brain-computer interface，BCI)[1]通過分析由運動想象產(chǎn)生的腦電信號就能讀取人腦發(fā)出的動作指令，從而替代大腦神經(jīng)與人體肌肉直接與外部設(shè)備建立連接通道，幫助肢體活動受限的人操作和控制一些輔助設(shè)備。然而由于腦電信號自身固有的非線性和非平穩(wěn)性，這給腦電信號特征提取與識別帶來了巨大的挑戰(zhàn)，同時也是運動想象BCI技術(shù)一直處于實驗室階段，不能走入實際應(yīng)用的主要原因[2]。

目前，研究人員已提出多種運動想象腦電信號的特征提取[3-7]方法。其中小波變換和小波包變換都是基于傅里葉變換的特征提取算法，不能在時域和頻域同時具有較高的分辨率。共同空間模式提取的特征缺乏相關(guān)的頻域信息，且需要大量的電極，無法應(yīng)用到少通道信號中[8]。希爾伯特-黃變換在時域和頻域都具有很高的分辨率，但是在EMD分解后需要人為選取IMF分量，由此可能導(dǎo)致重構(gòu)信號混入噪聲或丟失對分類有用的信息[9]。針對上述特征提取方法中的缺陷，提出了希爾伯特-黃變換與共同空間模式相結(jié)合(HCSP)的方法，結(jié)合自購實驗裝置EPOC，對表征手部運動想象腦電特征的C3、C4、CZ三通道分別進(jìn)行經(jīng)驗?zāi)B(tài)分解，選取分類效果最明顯的前三階IMF及其組合構(gòu)成15維的信號，應(yīng)用希爾伯特變換求取瞬時幅值后，利用CSP特征提取，聯(lián)合計算空間濾波器濾波后信號的AR模型參數(shù)和模糊熵[10,11]構(gòu)造特征向量。針對BCI competition 2003的數(shù)據(jù)，分別采用希爾伯特-黃變換、CSP和本文提出的方法進(jìn)行特征提取，再用線性分類器[12]進(jìn)行分類。利用本文提出的方法其分類準(zhǔn)確率較單獨使用HHT提高了5個百分點，較單獨使用CSP方法提高了15個百分點，分類結(jié)果表明此方法避免了人為選取IMF帶來的不確定因素，在減少分類導(dǎo)聯(lián)數(shù)的同時還能增加相關(guān)的頻域信息，可以提取對分類更為有效的信息。

1 腦電特征提取方法

1.1 HHT簡介

希爾伯特-黃變換是近年來信號處理領(lǐng)域的一個重大突破，自90年代末以來，已經(jīng)廣泛應(yīng)用于地震數(shù)據(jù)、氣候數(shù)據(jù)、語音信號、圖像信號等數(shù)據(jù)的分析，該方法在時域和頻域都具有較高的分辨率，完全適應(yīng)于腦電信號這類具有非線性以及非平穩(wěn)性特點的信號處理，HHT主要分為兩個部分：①經(jīng)驗?zāi)B(tài)分解(EMD)；②Hilbert譜分析[6](HSA)。

經(jīng)驗?zāi)B(tài)分解是希爾伯特-黃變換開始的第一步，也是整個特征提取最重要一步。利用經(jīng)驗?zāi)B(tài)分解可將一復(fù)雜多分量信號分解為一些簡單分量之和[13]，這些簡單分量必須滿足任意時刻上下包絡(luò)線相對于時間軸對稱，且其極值點數(shù)和過零點數(shù)相差最多為1，滿足上述條件的分量稱作固有模態(tài)函數(shù)(IMF)。

EMD的分解過程描述如下：

(1)對任一實信號x(t)，確定x(t)的所有極大值點和極小值點。

(2)對所有的極值點用三次樣條函數(shù)擬合成一條光滑的曲線形成信號的上下包絡(luò)線，計算其均值曲線，計算x(t)和均值曲線的差值h(t)。

(3)考察差值h(t)是否滿足IMF的條件。如果滿足條件直接作為IMF分量，如果不滿足條件，重新執(zhí)行步驟(1)、步驟(2)，直到滿足條件為止，并作為第一個IMF分量c1(t)，求得x(t)與IMF分量的差值r(t)。

(4)將差值r(t)作為原始信號重復(fù)上述的篩選，直至剩余分量為一個單調(diào)信號或者只有一個極值點為止，本文所用的停止條件：當(dāng)兩次篩選間的標(biāo)準(zhǔn)差SD小于0.2時即認(rèn)為篩選已經(jīng)結(jié)束。表達(dá)式為

(1)

在利用EMD分解得到各個IMF分量后，對每一個IMF分量做Hilbert變換即

(2)

構(gòu)造解析信號Zi(t)

Zi(t)=Ci(t)+jH[Ci(t)]=ai(t)ejθi(t)

(3)

進(jìn)而求得解析信號的瞬時相位和瞬時幅值

(4)

(5)

在傳統(tǒng)的HHT的特征提取的過程中通常是選取前幾階IMF分量的疊加來進(jìn)行分析，本文采用共同空間模式算法來選取對特征提取最有效的分量。

1.2 CSP簡介

CSP算法是近年來應(yīng)用于多通道分類的一種主流的分析方法，最初是應(yīng)用于BCI系統(tǒng)的腦電異常檢測，到后來逐步應(yīng)用于運動想象腦電信號的分類中。其核心思想請參考文獻(xiàn)[14]。設(shè)一次實驗的腦電信號表示為一個N×T維的矩陣X

(6)

式中：N是采樣通道數(shù)，T表示每個通道的采樣點數(shù)。Xl和Xr分別是想象左右手的腦電數(shù)據(jù)。CSP運算過程如下：

(1)計算想象左右手動作的空間協(xié)方差矩陣Cl、Cr

(7)

(8)

(2)計算出所有實驗的平均協(xié)方差矩陣組成混合空間協(xié)方差矩陣

(9)

對C進(jìn)行特征值分解，式中Σ為特征值對角矩陣，其中UO為特征值所對應(yīng)的特征向量。

(10)

(11)

(12)

(4)根據(jù)矩陣同時對角化原理：Sl、Sr應(yīng)該具有相同的特征向量，且特征值矩陣的和為單位陣

Sl=UλUT

(13)

Sr=U(I-λ)UT

(14)

(5)Ul和Ur分別是特征值矩陣λ和Ι-λ中最大特征值所對應(yīng)的特征向量，利用白化矩陣P構(gòu)造空間濾波器Wl、Wr

(15)

(16)

(6)則原始的腦電信號經(jīng)過濾波器濾波，求得特征向量Zl、Zr

Zl=Wl×X

(17)

Zr=Wr×X

(18)

(7)進(jìn)一步提取出2維特征向量

(19)

Var(X)是計算X的方差。

利用矩陣對角化原理，CSP算法可實現(xiàn)空間濾波器尋優(yōu)，利用該最優(yōu)濾波器獲取更明顯的特征向量，但CSP算法缺乏相關(guān)的頻域信息而且只有在輸入導(dǎo)聯(lián)數(shù)達(dá)到一定數(shù)量時，才能達(dá)到較佳的分類效果，這也限制了共同空間模式算法的應(yīng)用，制約了BCI系統(tǒng)的發(fā)展。

1.3 HHT和CSP相結(jié)合的特征提取算法

本文提出一種HHT與CSP相結(jié)合的特征提取算法，首先采用EMD方法將3個導(dǎo)聯(lián)的原始腦電信號分解為眾多簡單分量的疊加，提取每個通道對運動想象分類有效的IMF分量及它們的組合組成5維重構(gòu)信號，可假定15個通道的腦電信號來做CSP特征提取，算法流程如圖1所示。

圖1 特征提取算法流程

研究表明，運動想象腦電信號主要發(fā)生在8 Hz～30 Hz，本文在特征提取之前首先對原始腦電信號進(jìn)行8 Hz～30 Hz的帶通濾波，采用butterworth濾波器，設(shè)置通帶截止頻率8 Hz～30 Hz，阻帶截止頻率5 Hz～35 Hz，通帶衰減0.5 db，阻帶衰減50 db。然后對濾波后的信號進(jìn)行經(jīng)驗?zāi)B(tài)分解，分解后的各階IMF如圖2所示。

前期研究表明，只有前三階分量的頻率范圍屬于運動想象的節(jié)律段，選取前三階IMF分量及它們的組合，進(jìn)行希爾伯特變換提取瞬時幅值，由于腦電特征的瞬時幅值波動較大，本文采用1S的滑動時間窗，滑動步長為一個采樣點，提取平均瞬時能量并進(jìn)行歸一化處理來表征分類特征信息。圖3是歸一化后的能量隨時間變化的曲線。

選擇與手部運動想象最相關(guān)的3個導(dǎo)聯(lián)(C3、CZ、C4)，進(jìn)行EMD分解，每個導(dǎo)聯(lián)提取前三階固有模態(tài)分量及前三階分量的組合，組成一15維向量E15×T作為原始腦電信號。在運動想象信號的研究中，ERD/ERS現(xiàn)象并非存在于整個想象過程中。由圖3可知，在4 s至7 s期間，C3、C4通道的能量差異最大，可以認(rèn)為在此時間段的ERD/ERS現(xiàn)象最明顯。本文選取此時間段內(nèi)200個采樣點數(shù)據(jù)(即601-800點)，并對E15×200進(jìn)行CSP特征提取。本文選取特征值矩陣中最大和最小特征值對應(yīng)的特征向量構(gòu)成空間濾波器，對E15×200進(jìn)行濾波處理，濾波后的信號為Z2×200，進(jìn)一步提取出2維特征向量fp。

圖2 腦電信號EMD后各階IMF

圖3 左右手運動想象腦電能量曲線

1.4 AR模型特征提取

進(jìn)而對經(jīng)空間濾波器濾波后的信號Z2×200提取自回歸(auto regressive，AR)模型參數(shù)，本文分別使用3、5、6、8、10、12、14、16階AR模型提取AR模型參數(shù)，分別求取最后的分類準(zhǔn)確率[6]，如圖4所示。

圖4 各階AR模型分類準(zhǔn)確率

鑒于在階次為6時對應(yīng)了最大的分類準(zhǔn)確率，本文基于Burg算法提取6階AR模型系數(shù)組成12維的特征向量：{C3AR1，…，C3AR6，C4AR1，…，C4AR6}。

1.5 模糊熵特征提取

腦電信號是一個非線性信號，它還具有一定的非線性特征[15]。熵是一種非線性動力學(xué)參數(shù)，可以用于衡量時間序列中產(chǎn)生新模式概率的大小，得出腦電信號的復(fù)雜度。在對腦電信號進(jìn)行非線性檢測時，我們通常提取信號的近似熵、樣本熵、模糊熵作為信號的非線性特征[10]。近似熵(approximate entropy，ApEn)是Pincus在20年代末提取的一種非線性指標(biāo)，可以對離散時間序列進(jìn)行估計，計量信號的復(fù)雜性及不規(guī)則性質(zhì)。在近似熵的基礎(chǔ)上Richman和Moorman提出了一種時間序列復(fù)雜性測量方法——樣本熵(sample entropy，SampEn)，與近似熵相比，樣本熵在計算過程中對數(shù)據(jù)長度的敏感性降低，擴(kuò)大了非線性檢測的適用范圍，而且計算過程略去了自匹配和自適應(yīng)模板算法，消除了了計算過程帶來的誤差。本文使用模糊熵提取腦電信號的非線性特征，模糊熵(fuzzy entropy，F(xiàn)uzzyEn)是Chen等針對樣本熵和近似熵提取過程使用二值函數(shù)進(jìn)行相似性度量導(dǎo)致熵值不連續(xù)而提出的一種改進(jìn)算法[11]。模糊熵在計算向量的相似度時用模糊隸屬度函數(shù)取代Heaviside函數(shù)，克服了二值函數(shù)方法所導(dǎo)致的缺乏連續(xù)性問題。該算法不僅繼承了樣本熵算法的優(yōu)點，而且對時間序列的長度的依賴性更小，對含噪信號的魯棒性更加優(yōu)良。模糊熵算法的具體步驟如下：

(1)設(shè)N點采樣序列為

{u(i):i=1～N}

(20)

(2)對上述時間序列按序號連續(xù)順序進(jìn)行相空間重構(gòu)，組成m維和m+1維的矢量即

(21)

(22)

其中，i=1，2，…N-m，其中um(i)和um+1(i)分別為m維和m+1維矢量的均值。

(23)

(24)

其中，i，j=1～N-m,i≠j。

(25)

(26)

(5)定義函數(shù)

(27)

(28)

(6)定義模糊熵為

(29)

當(dāng)序列長度為有限值時，可以得出模糊熵的估計值為

F(m,n,r,N)=lnΦm(n,r)-lnΦm+1(n,r)

(30)

上述函數(shù)中，m是相空間重構(gòu)的維數(shù)，r是相似容限度，相似容限r(nóng)選擇不當(dāng)會對模糊熵的估計造成影響，r的取值范圍一般為0.1～0.25倍采樣序列標(biāo)準(zhǔn)差，本文為了獲取更大足夠的序列信息，m和n的取值均為2，r為標(biāo)準(zhǔn)差的1/5。如圖5所示是想象左右手運動C3、C4通道的平均模糊熵歸一化值。

圖5 想象左右手運動EEG平均模糊熵值

提取二維特征向量{FuzzyC3,FuzzyC4}，表1為分別使用近似熵、樣本熵、模糊熵的分類結(jié)果。

表1 不同非線性特征提取方式下的分類準(zhǔn)確率/%

從上表可知，相對于近似熵和樣本熵，模糊熵更適合用于衡量信號的復(fù)雜性，提取的非線性特征的分類正確率也更高。從理論上和實踐上證明了模糊熵優(yōu)于近似熵和樣本熵。

2 實驗結(jié)果

為驗證本文所述方法的有效性，對BCI competition 2003的左右手運動想象腦電數(shù)據(jù)進(jìn)行了特征提取與模式分類。將對上述CSP算法提取的二維特征向量、AR模型提取的12維特征向量和模糊熵算法提取的2維特征向量組成最終的16維特征向量{fp，C3AR1，…，C3AR6，C4AR1，…，C4AR6，F(xiàn)uzzyC3，F(xiàn)uzzyC4}，形成腦電信號的多特征值提取。對提取的特征向量采用線性分類器(LDA)進(jìn)行分類。線性分類器執(zhí)行速度快，推廣能力強(qiáng)，是腦機(jī)接口模式識別中最常用分類算法，在對實時性要求極高的在線系統(tǒng)中應(yīng)用尤其廣泛。LDA算法即將一個d維的樣本投影到一條直線上，使得在這條直線上的樣本投影能最大尺度的被分開，降低了原始信號的維度。LDA分類器的核心就是找到這一最佳投影線，能將兩類樣本最好地分開。表2給出了分別使用傳統(tǒng)HHT，傳統(tǒng)CSP和改進(jìn)HCSP的算法下得到的分類準(zhǔn)確率。

表2 不同特征提取方法下的正確識別率/%

可見當(dāng)輸入信號導(dǎo)聯(lián)數(shù)較少，單一的CSP特征提取算法的分類效果最差，無法提取到有效特征向量。使用單一的HHT算法的分類準(zhǔn)確率有了明顯提高，表明HHT算法在分析非線性非平穩(wěn)信號時，在時域和頻域都可以具有較高的分類正確率，可以有效提取運動想象腦電的特征模式。本文在分析了HHT和CSP算法的優(yōu)缺點的基礎(chǔ)上提出了將兩者結(jié)合的特征提取算法，從表2可知本文提出的HCSP特征提取方法是有效的，可以提取更為有效的分類特征。

3 結(jié)束語

本文利用HHT與CSP混合特征提取方法，獲取EEG時頻空域信息，既解決了人為選取IMF分量時可能混入噪聲抑或丟失有用信息的問題，同時減少了輸入信號的導(dǎo)聯(lián)數(shù)。HHT算法側(cè)重于提取腦電信號中具有可分類信息的瞬時幅值，CSP算法能夠構(gòu)造適用于分類的空間濾波器，同時提取AR模型參數(shù)和模糊熵。利用BCI competition 2003中提供的數(shù)據(jù)使用本文提出的特征提取方法，獲得了良好的效果。

[1]Long Jinyi,Li Yuanqi,Wang Hongtao,et al.A hybrid brain computer interface to control the direction and speed of a simulated or real wheelchair[J].Neural Systems and Rehabilitation Engineering,2012,20(5):720-729.

[2]ZHAO Limin.Research on EEG feature extraction and classi-fication based on movement imagination[D].Taiyuan:Taiyuan University of Technology,2016(in Chinese).[趙利民.基于運動想象的腦電信號特征提取與分類方法研究[D].太原:太原理工大學(xué),2016.]

[3]Atyabi A,Shic F,Naples A.Mixture of autoregressive mode-ling orders and its implication on single trial EEG classification[J].Expert Systems with Applications,2016,65:164-180.

[4]Zhou Z,Wan B.Wavelet packet-based independent component analysis for feature extraction from motor imagery EEG of complex movements[J].Clinical Neurophysiology Official Journal of the International Federation of Clinical Neurophysiology,2012,123(9):1779-1788.

[5]ZHENG Shuhua,YAN Chen,WANG Xiangzhou.A repeated bisection CSP feature extraction algorithm of four-class motor imagery EEG[J].Transactions of Beijing Institute of Techno-logy,2016(8):844-850(in Chinese).[鄭戍華,閆琛,王向周.一種重復(fù)二分CSP4類運動想象腦電信號特征提取算法[J].北京理工大學(xué)學(xué)報,2016(8):844-850.]

[6]SUN Huiwen,FU Yunfa,XIONG XIN,et al.Identification of EEG induced by motor imagery based on Hilbert-Huang transform[J].Acta Automatica Sinica,2015,41(9):1686-1692(in Chinese).[孫會文,伏云發(fā),熊馨,等.基于HHT運動想象腦電模式識別研究[J].自動化學(xué)報,2015,41(9):1686-1692.]

[7]JI Yu.Research of P300 processing algorithm based on independent component analysis[D].Hangzhou:Zhejiang University,2013(in Chinese).[計瑜.基于獨立分量分析的P300腦電信號處理算法研究[D].杭州:浙江大學(xué),2013.]

[8]YANG Banghua,LU Wenyu,HE Meiyan,et al.Novel feature extraction method for BCI based on WPD and CSP[J].Chinese Journal of Scientific Instrument,2012,33(11):2560-2565(in Chinese).[楊幫華,陸文宇,何美燕,等.腦機(jī)接口中基于WPD和CSP的特征提取[J].儀器儀表學(xué)報,2012,33(11):2560-2565.]

[9]Park C,Looney D,Rehman NU,et al.Classification of motor imagery BCI using multivariate empirical mode decomposition[J].IEEE Transactions on Neural Systems & Rehabilitation Engineering,2013,21(1):10-22.

[10]LI Conggai.The comparation and application of entropy in the detection of epilepsy[D].Taiyuan:Taiyuan University of Technology,2015(in Chinese).[李聰改.熵在癲癇信號檢測中的對比研究與應(yīng)用[D].太原:太原理工大學(xué),2015.]

[11]TIAN Jing,LUO Zhizeng.Motor imagery EEG feature extraction based on fuzzy entropy[J].Journal of Huazhong University of Science and Technology(Nature Science Edition),2013,41(s1):92-94(in Chinese).[田京,羅志增.基于模糊熵的運動想像腦電信號特征提取[J].華中科技大學(xué)學(xué)報(自然科學(xué)版),2013,41(s1):92-94.]

[12]LI Lijun.Feature extraction and classification of imaginary movements in EEG[D].Guangzhou:South China University of Technology,2012(in Chinese).[李麗君.基于運動想象的腦電信號特征提取及分類算法研究[D].廣州:華南理工大學(xué),2012.]

[13]Bajaj V,Pachori RB.Classification of seizure and non-seizure EEG signals using empirical mode decomposition[J].IEEE Transactions on Information Technology in Biomedicine,2012,16(6):1135-1142.

[14]Song X,Yoon SC,Perera V.Adaptive common spatial pattern for single-trial EEG classification in multisubject BCI[C]//International IEEE/EMBS Conference on Neural Engineering.IEEE,2013:411-414.

[15]LIN Haibo,GONG Lu,ZHANG Yi,et al.Feature extraction of EEG signal based on improved HHT and sample entropy[J].Computer Engineering and Design,2015,36(6):1608-1613(in Chinese).[林海波,龔璐,張毅,等.基于改進(jìn)HHT和樣本熵的腦電信號特征提取[J].計算機(jī)工程與設(shè)計,2015,36(6):1608-1613.]