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

于沐涵，陳 峰

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

0 引 言

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

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

1 腦電特征提取方法

1.1 HHT簡介

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

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

EMD的分解過程描述如下：

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

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

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

(4)將差值r(t)作為原始信號重復上述的篩選，直至剩余分量為一個單調信號或者只有一個極值點為止，本文所用的停止條件：當兩次篩選間的標準差SD小于0.2時即認為篩選已經結束。表達式為

(1)

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

(2)

構造解析信號Zi(t)

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

(3)

進而求得解析信號的瞬時相位和瞬時幅值

(4)

(5)

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

1.2 CSP簡介

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

(6)

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

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

(7)

(8)

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

(9)

對C進行特征值分解，式中Σ為特征值對角矩陣，其中UO為特征值所對應的特征向量。

(10)

(11)

(12)

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

Sl=UλUT

(13)

Sr=U(I-λ)UT

(14)

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

(15)

(16)

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

Zl=Wl×X

(17)

Zr=Wr×X

(18)

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

(19)

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

利用矩陣對角化原理，CSP算法可實現空間濾波器尋優，利用該最優濾波器獲取更明顯的特征向量，但CSP算法缺乏相關的頻域信息而且只有在輸入導聯數達到一定數量時，才能達到較佳的分類效果，這也限制了共同空間模式算法的應用，制約了BCI系統的發展。

1.3 HHT和CSP相結合的特征提取算法

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

圖1 特征提取算法流程

研究表明，運動想象腦電信號主要發生在8 Hz～30 Hz，本文在特征提取之前首先對原始腦電信號進行8 Hz～30 Hz的帶通濾波，采用butterworth濾波器，設置通帶截止頻率8 Hz～30 Hz，阻帶截止頻率5 Hz～35 Hz，通帶衰減0.5 db，阻帶衰減50 db。然后對濾波后的信號進行經驗模態分解，分解后的各階IMF如圖2所示。

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

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

圖2 腦電信號EMD后各階IMF

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

1.4 AR模型特征提取

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

圖4 各階AR模型分類準確率

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

1.5 模糊熵特征提取

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

(1)設N點采樣序列為

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

(20)

(2)對上述時間序列按序號連續順序進行相空間重構，組成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)定義函數

(27)

(28)

(6)定義模糊熵為

(29)

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

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

(30)

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

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

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

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

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

2 實驗結果

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

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

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

3 結束語

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

[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].太原:太原理工大學,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).[鄭戍華,閆琛,王向周.一種重復二分CSP4類運動想象腦電信號特征提取算法[J].北京理工大學學報,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).[孫會文,伏云發,熊馨,等.基于HHT運動想象腦電模式識別研究[J].自動化學報,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].杭州:浙江大學,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).[楊幫華,陸文宇,何美燕,等.腦機接口中基于WPD和CSP的特征提取[J].儀器儀表學報,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).[李聰改.熵在癲癇信號檢測中的對比研究與應用[D].太原:太原理工大學,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].華中科技大學學報(自然科學版),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].廣州:華南理工大學,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).[林海波,龔璐,張毅,等.基于改進HHT和樣本熵的腦電信號特征提取[J].計算機工程與設計,2015,36(6):1608-1613.]