于沐涵，陳 峰

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

0 引 言

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

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

1 腦電特征提取方法

1.1 HHT簡(jiǎn)介

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

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

EMD的分解過(guò)程描述如下：

(1)對(duì)任一實(shí)信號(hào)x(t)，確定x(t)的所有極大值點(diǎn)和極小值點(diǎn)。

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

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

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

(1)

在利用EMD分解得到各個(gè)IMF分量后，對(duì)每一個(gè)IMF分量做Hilbert變換即

(2)

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

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

(3)

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

(4)

(5)

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

1.2 CSP簡(jiǎn)介

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

(6)

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

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

(7)

(8)

(2)計(jì)算出所有實(shí)驗(yàn)的平均協(xié)方差矩陣組成混合空間協(xié)方差矩陣

(9)

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

(10)

(11)

(12)

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

Sl=UλUT

(13)

Sr=U(I-λ)UT

(14)

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

(15)

(16)

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

Zl=Wl×X

(17)

Zr=Wr×X

(18)

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

(19)

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

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

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

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

圖1 特征提取算法流程

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

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

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

圖2 腦電信號(hào)EMD后各階IMF

圖3 左右手運(yùn)動(dòng)想象腦電能量曲線

1.4 AR模型特征提取

進(jìn)而對(duì)經(jīng)空間濾波器濾波后的信號(hào)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時(shí)對(duì)應(yīng)了最大的分類準(zhǔn)確率，本文基于Burg算法提取6階AR模型系數(shù)組成12維的特征向量：{C3AR1，…，C3AR6，C4AR1，…，C4AR6}。

1.5 模糊熵特征提取

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

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

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

(20)

(2)對(duì)上述時(shí)間序列按序號(hào)連續(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)序列長(zhǎng)度為有限值時(shí)，可以得出模糊熵的估計(jì)值為

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)會(huì)對(duì)模糊熵的估計(jì)造成影響，r的取值范圍一般為0.1～0.25倍采樣序列標(biāo)準(zhǔn)差，本文為了獲取更大足夠的序列信息，m和n的取值均為2，r為標(biāo)準(zhǔn)差的1/5。如圖5所示是想象左右手運(yùn)動(dòng)C3、C4通道的平均模糊熵歸一化值。

圖5 想象左右手運(yùn)動(dòng)EEG平均模糊熵值

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

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

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

2 實(shí)驗(yàn)結(jié)果

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

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

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

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

本文利用HHT與CSP混合特征提取方法，獲取EEG時(shí)頻空域信息，既解決了人為選取IMF分量時(shí)可能混入噪聲抑或丟失有用信息的問(wèn)題，同時(shí)減少了輸入信號(hào)的導(dǎo)聯(lián)數(shù)。HHT算法側(cè)重于提取腦電信號(hào)中具有可分類信息的瞬時(shí)幅值，CSP算法能夠構(gòu)造適用于分類的空間濾波器，同時(shí)提取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).[趙利民.基于運(yùn)動(dòng)想象的腦電信號(hào)特征提取與分類方法研究[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類運(yùn)動(dòng)想象腦電信號(hào)特征提取算法[J].北京理工大學(xué)學(xué)報(bào),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).[孫會(huì)文,伏云發(fā),熊馨,等.基于HHT運(yùn)動(dòng)想象腦電模式識(shí)別研究[J].自動(dòng)化學(xué)報(bào),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).[計(jì)瑜.基于獨(dú)立分量分析的P300腦電信號(hào)處理算法研究[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é)報(bào),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).[李聰改.熵在癲癇信號(hào)檢測(cè)中的對(duì)比研究與應(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).[田京,羅志增.基于模糊熵的運(yùn)動(dòng)想像腦電信號(hào)特征提取[J].華中科技大學(xué)學(xué)報(bào)(自然科學(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).[李麗君.基于運(yùn)動(dòng)想象的腦電信號(hào)特征提取及分類算法研究[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和樣本熵的腦電信號(hào)特征提取[J].計(jì)算機(jī)工程與設(shè)計(jì),2015,36(6):1608-1613.]