付榮榮， 楊 陽， 于 寶， 劉 沖， 張 馳

(1.燕山大學(xué) 電氣工程學(xué)院，河北 秦皇島 066004；2.東北大學(xué) 機(jī)械工程與自動(dòng)化學(xué)院，遼寧 沈陽 110819；3.大連理工大學(xué) 生物醫(yī)學(xué)工程學(xué)院，遼寧 大連 116024)

1 引 言

腦電信號(hào)(electroencephalography，EEG)是由大腦神經(jīng)元活動(dòng)產(chǎn)生的信號(hào)，其中包含了豐富的大腦狀態(tài)信息，為實(shí)現(xiàn)腦-機(jī)接(brain-computer interface，BCI)需要對(duì)腦電信號(hào)進(jìn)行有效解碼，解碼過程包括對(duì)EEG信號(hào)進(jìn)行特征提取和模式分類。傳統(tǒng)的BCI多通道EEG信號(hào)特征提取方法有AR模型、時(shí)頻濾波方法、共空間模式、主成分分析和2D主成分分析(2DPCA)等[1～3]，這些從時(shí)頻域或時(shí)空域提取特征方法都已取得不錯(cuò)的效果。但是，這些二階特征提取方法必然會(huì)忽略頻域或空間域的特征，限制了所提取特征被準(zhǔn)確識(shí)別的上限。近年來，Lu H等[4]提出了多線性主成分分析(multi-linear principal component analysis，MPCA)等高階特征提取方法，這種高階特征提取方法不需要將原始數(shù)據(jù)展開成向量形式，保留了數(shù)據(jù)的結(jié)構(gòu)特性。隨后，一些研究小組開始研究時(shí)空頻高階特征提取方法，文獻(xiàn)[5,6]提出了一種基于小波和獨(dú)立分量分析的時(shí)間頻率空間EEG特征的提取方法，分別用離散小波變換和獨(dú)立分量分析提取時(shí)頻域特征和空域特征，有效地克服了傳統(tǒng)的時(shí)頻特征提取方法空間信息描述不足等問題；王金甲等[7]首先利用復(fù)Morlet小波分析方法生成張量數(shù)據(jù)，然后進(jìn)行張量降維并提取特征，克服了主成分分析維度小的局限性?？梢姼唠A特征提取方法在單次運(yùn)動(dòng)想象腦電信號(hào)的特征提取中是極為有效的，但小波分解方法構(gòu)建高階張量時(shí)存在小波參數(shù)難以確定的缺陷。

針對(duì)上述問題，本文提出了一種新的特征提取方法，這種方法通過集合經(jīng)驗(yàn)?zāi)B(tài)分解[8～11](Ensemble Empirical Mode Decomposition，EEMD)提取EEG頻域信息，構(gòu)造包含頻域、時(shí)域、空間域的三維EEG信號(hào)；并結(jié)合MPCA對(duì)EEG信號(hào)進(jìn)行降維。這種新的高階特征提取方法無需將原始腦電信號(hào)張量展開為向量，能夠分別從時(shí)域、頻域和空間域提取腦電信號(hào)的多模態(tài)特征。通過與使用小波分解構(gòu)建EEG高階張量并使用EEMD進(jìn)行降維的特征提取方法和直接使用2DPCA[12]方法進(jìn)行降維的特征提取方法進(jìn)行對(duì)比，驗(yàn)證了本文提出方法的有效性，為單次運(yùn)動(dòng)想象腦電信號(hào)的準(zhǔn)確識(shí)別奠定了基礎(chǔ)。

2 實(shí)驗(yàn)方法

2.1 腦電信號(hào)的采集

實(shí)驗(yàn)共7名受試者參與腦電信號(hào)的采集，年齡在22～27歲之間，實(shí)驗(yàn)前已被告知實(shí)驗(yàn)過程并簽署知情同意文件。實(shí)驗(yàn)采用虛擬碗球?qū)嶒?yàn)，這是一種動(dòng)態(tài)復(fù)雜約束特性的腦電實(shí)驗(yàn)，實(shí)驗(yàn)環(huán)境包括演示刺激系統(tǒng)和腦電信號(hào)采集系統(tǒng)兩部分。演示刺激系統(tǒng)為一臺(tái)在PsychToolbox環(huán)境下運(yùn)行虛擬碗球系統(tǒng)的計(jì)算機(jī)，腦電采集系統(tǒng)為Emotiv的EPOC+。

在正式實(shí)驗(yàn)之前，參加實(shí)驗(yàn)的7名受試者需要先進(jìn)行虛擬碗球?qū)嶒?yàn)的模擬練習(xí)。在整個(gè)實(shí)驗(yàn)中，受試者正坐在椅子上，目視前方的演示刺激系統(tǒng)電腦屏幕。

演示刺激系統(tǒng)電腦屏幕上顯示內(nèi)容為虛擬碗球系統(tǒng)，受試者需要控制碗在水平面上從起點(diǎn)到達(dá)終點(diǎn)，碗里的球此時(shí)會(huì)由于慣性在碗中晃動(dòng)，受試者同時(shí)還必須將球約束在碗中，如果球從碗中掉出，實(shí)驗(yàn)會(huì)被立即終止。按壓電腦鍵盤的左右鍵就能控制碗左右移動(dòng)，在按壓左右鍵的同時(shí)，受試者還需要想象對(duì)應(yīng)的左右運(yùn)動(dòng)。在實(shí)驗(yàn)過程中，通過Emotiv的EPOC+對(duì)腦電信號(hào)進(jìn)行采集。信號(hào)采集示意圖如圖1所示。

圖1 EEG信號(hào)采集示意圖Fig.1 EEG signal acquisition diagram

2.2 腦電信號(hào)的預(yù)處理

通過腦電采集設(shè)備采集到的腦電信號(hào)非常微弱而且含有噪聲，需要對(duì)信號(hào)進(jìn)行預(yù)處理[13]。由于運(yùn)動(dòng)想象誘發(fā)的ERD/ERS[14]現(xiàn)象主要集中在Mu節(jié)律和Beta節(jié)律，即要分析的信號(hào)頻率范圍為8～30 Hz，故采用8～30 Hz的6階巴特沃斯帶通濾波器對(duì)采集到的原始數(shù)據(jù)進(jìn)行濾波。經(jīng)過Emotiv的EPOC+能夠采集到14導(dǎo)聯(lián)的腦電數(shù)據(jù)，采樣頻率為128 Hz。經(jīng)過預(yù)處理后將14通道的腦電數(shù)據(jù)[15]劃分為2個(gè)類別，將數(shù)據(jù)逐次截取每一類別的前1 s數(shù)據(jù)，整理為多個(gè)四維張量結(jié)構(gòu)(樣本點(diǎn)數(shù)目×導(dǎo)聯(lián)數(shù)目×實(shí)驗(yàn)次數(shù)×類別)。

2.3 對(duì)照實(shí)驗(yàn)

為了驗(yàn)證新提出的特征提取方法的有效性，設(shè)計(jì)了3種不同特征提取方法的對(duì)照實(shí)驗(yàn)。實(shí)驗(yàn)1為本文提出的特征提取方法，采用EEMD構(gòu)建腦電信號(hào)高階張量，采用EEMD進(jìn)行降維，流程圖如圖2所示；實(shí)驗(yàn)2采用復(fù)Morlet小波分解構(gòu)建腦電信號(hào)高階張量，采用MPCA進(jìn)行降維；實(shí)驗(yàn)3直接對(duì)采集到的EEG信號(hào)進(jìn)行2DPCA進(jìn)行降維。特征提取后都采用Fisher線性判別分析分類方法取得分類準(zhǔn)確率。由于EEG數(shù)據(jù)集采集難度大，數(shù)據(jù)集較小，為獲得更佳的模型，本實(shí)驗(yàn)采用10折交叉驗(yàn)證將數(shù)據(jù)隨機(jī)分為90%的訓(xùn)練集和10%的測(cè)試集進(jìn)行分類。

圖2 本文提出的特征提取方法流程圖Fig.2 The flow chart of feature extraction method proposed

3 算法原理

3.1 EEMD方法

對(duì)單次運(yùn)動(dòng)想象的腦電信號(hào)，其數(shù)據(jù)形式為樣本點(diǎn)×導(dǎo)聯(lián)數(shù)目，分別表示腦電信號(hào)的時(shí)域信息和空間域信息。為了獲得更加全面的腦電信息，需要對(duì)采集到的時(shí)域信息進(jìn)行時(shí)頻變換。本方法采用EEMD進(jìn)行時(shí)頻變換，EEMD方法依據(jù)數(shù)據(jù)自身的時(shí)間尺度特征進(jìn)行信號(hào)分解，無需設(shè)定預(yù)先設(shè)定的任何基函數(shù)，對(duì)于處理非線性非平穩(wěn)時(shí)間序列非常有效。經(jīng)過EEMD處理后，腦電數(shù)據(jù)同時(shí)包含了腦電的時(shí)域信息，頻域信息和空間域信息。

EEMD的原理是先將頻譜均勻分布的白噪聲加入原信號(hào)中作為新的待分解信號(hào)，這樣不同時(shí)間尺度的信號(hào)就自動(dòng)分布到了合適的參考尺度。再將信號(hào)進(jìn)行經(jīng)驗(yàn)?zāi)B(tài)分解(empirical mode decompose-tion，EMD)[16～18]處理，由于白噪聲零均值的特性，最后取均值得到逼近真實(shí)的模態(tài)。具體步驟如下：

(1) 往信號(hào)x(t)中輸入白噪聲信號(hào)k·σx·n(t)，其中：n(t)為歸一化白噪聲，σx為信號(hào)標(biāo)準(zhǔn)差，k為比例系數(shù)，得到信噪混合函數(shù)，

y(t)=x(t)+k·σx·n(t)

(1)

(2) 確定信號(hào)y(t)的所有局部極值點(diǎn)。

(3) 利用3次樣條插值將所有局部極值點(diǎn)連接起來分別形成上、下包絡(luò)線包絡(luò)所有的數(shù)據(jù)點(diǎn)。

(4) 上、下包絡(luò)線的平均值記為m1(t)，

h1(t)=y(t)-m1(t)

(2)

(5) 如果h1(t)滿足在整個(gè)數(shù)據(jù)段內(nèi)，極值點(diǎn)的個(gè)數(shù)和過零點(diǎn)的個(gè)數(shù)必須相等或相差最多不超過一個(gè)且上、下包絡(luò)線相對(duì)于時(shí)間軸對(duì)稱，則h1(t)就是y(t)的第一個(gè)本征模態(tài)函數(shù)(IMF)分量，否則y(t)=h1(t)，重復(fù)(2)至(4)。

(6) 將h1(t)從y(t)中分離出來得到，

r1(t)=y(t)-h1(t)

(3)

重復(fù)(2)～(5)。

(7) 最終將y(t)分解成IMF的組合，

(4)

式中rm(t)為分解m次后的剩余分量。

3.2 MPCA方法

得到了包含腦電的時(shí)域信息，頻域信息和空間域信息的腦電信號(hào)后，需要對(duì)腦電信號(hào)進(jìn)行降維，從而剔除對(duì)腦電信號(hào)分類影響很小甚至誤導(dǎo)分類的部分特征(噪聲)。

本方法對(duì)腦電數(shù)據(jù)使用MPCA進(jìn)行降維。MPCA算法的目的是找到多個(gè)模態(tài)的投影矩陣1U(1)T，21U(2)T，31U(3)T，…，n1U(n)T來對(duì)原始張量數(shù)據(jù)進(jìn)行降維，并盡量保持原張量組的離散度。由于本研究腦電信號(hào)提取了頻域、時(shí)域、空間域三維信息，以三階張量為例，具體步驟如下：

(1) 對(duì)所述張量樣本X∈1Uc×f×t×s進(jìn)行中心化處理，其中，c為空間域信號(hào)，f為頻域信號(hào)，t為時(shí)域信號(hào)，s為實(shí)驗(yàn)次數(shù)。

(5)

其中，

(6)

(2) 計(jì)算初始協(xié)方差矩陣，

(7)

式中：Xm(n)為Xmn模展開后的矩陣，并對(duì)所述初始協(xié)方差矩陣進(jìn)行特征分解，取最大的d′個(gè)特征值對(duì)應(yīng)的特征向量組成投影矩陣1U(n)(n=1，2，3)得到初始化的降維投影矩陣。

潑水節(jié)是傣族人民的傳統(tǒng)節(jié)日，傣族人民喜水，不僅與氣候環(huán)境文化傳統(tǒng)有關(guān)，更是一種風(fēng)俗習(xí)慣，傣族是一個(gè)熱愛水的民族。在傣族社會(huì)中,水不僅僅是一種自然的物質(zhì),同時(shí)也被賦予了豐富的文化內(nèi)涵,從而形成了自身的水文化。水文化作為一種重要的社會(huì)傳統(tǒng),在保持傣族居住區(qū)人與自然環(huán)境的和諧起到了積極的作用。歷史上傣族和水就有一種特殊的關(guān)系,并且傣族擁有生動(dòng)的水文化,因此大多數(shù)傣族被其他民族生動(dòng)地稱為“水傣”。根據(jù)傣族人類起源的傳說,天空和大地都起源于水,而且認(rèn)為人類生命的一半也是由水創(chuàng)造。因此水是無比偉大和圣潔的。

(3) 對(duì)步驟(2)得到的所述初始化降維投影矩陣進(jìn)行局部?jī)?yōu)化，包括如下具體步驟：

① 進(jìn)行投影，

(8)

式中左下標(biāo)1、2、3表示1模、2模、3模乘積。

② 計(jì)算總散度，

(9)

③ 對(duì)n=1、2、3，計(jì)算投影后張量n模展開后協(xié)方差矩陣的特征值，取d′個(gè)所述特征值對(duì)應(yīng)的特征向量組成新的投影矩陣1U(n)(n=1，2，3)，用所述新的投影矩陣更新Sm并計(jì)算新的ψSk；

④ 判斷ψSk-ψSk-1<η是否成立，式中k為優(yōu)化迭代次數(shù)，η為自定義的閾值，如果ψSk-ψSk-1<η，得到最終的投影矩陣11U(1)T，21U(2)T，31U(3)T，計(jì)算

(10)

得到的投影矩陣21U(2)T，21U(2)T，31U(3)T即為MPCA要求的投影矩陣。對(duì)原始數(shù)據(jù)進(jìn)行投影就得到了降維后的腦電數(shù)據(jù)。

3.3 Fisher線性判別分析

降維后的腦電數(shù)據(jù)是三階多維腦電信號(hào)，需要先將其展開為一階多維信號(hào)，然后使用分類器進(jìn)行分類。本實(shí)驗(yàn)采用Fisher線性判別分析[13]方法對(duì)數(shù)據(jù)進(jìn)行分類。線性判別分析的思想為：使投影后類內(nèi)間距最小，類間方差間距最大。具體步驟如下：

(1) 計(jì)算出兩類數(shù)據(jù)特征向量的均值向量

(11)

式中Nj為第j類特征向量對(duì)應(yīng)的數(shù)據(jù)點(diǎn)數(shù)。

(2) 將二者方差加起來求出類內(nèi)散度矩陣

(12)

(3) 計(jì)算類間散度矩陣Sb

Sb=(m0-m1)(m0-m1)T

(13)

(5) 計(jì)算投影后的數(shù)據(jù)點(diǎn)Y。

Y=WTX

(14)

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

4.1 分類結(jié)果評(píng)價(jià)

表1為實(shí)驗(yàn)分類結(jié)果，ACC(Accuracy)為準(zhǔn)確度，AUC(Area Under Curve)為受試者工作特性ROC(receiver operating characteristic)曲線下的面積，ACC值和AUC值越接近1代表分類結(jié)果越好。

表1 實(shí)驗(yàn)分類結(jié)果Tab.1 Classification results of experiment

對(duì)比3種實(shí)驗(yàn)的AUC和ACC結(jié)果如圖3所示。

由表1可以看出，實(shí)驗(yàn)1(本文提出的EEMD+MPCA方法)中8個(gè)受試者EEG信號(hào)的識(shí)別ACC均值為0.933 4，AUC均值為0.950 2都高于實(shí)驗(yàn)2(小波+MPCA)和實(shí)驗(yàn)3(2DPCA)的結(jié)果，由圖3能夠直觀地看出實(shí)驗(yàn)1獲得的AUC和ACC在每個(gè)實(shí)驗(yàn)對(duì)象上都是最高的，即實(shí)驗(yàn)1的方法比實(shí)驗(yàn)2和實(shí)驗(yàn)3能取得更好的分類效果。這個(gè)結(jié)果說明，采用EEMD構(gòu)建腦電信號(hào)高階張量，采用MPCA進(jìn)行降維的特征提取方法相對(duì)于采用復(fù)Morlet小波分解構(gòu)建腦電信號(hào)高階張量，采用MPCA進(jìn)行降維和直接對(duì)采集到的EEG信號(hào)進(jìn)行2DPCA進(jìn)行降維的特征提取方法而言能夠取得更好的特征，在單次運(yùn)動(dòng)想象腦電信號(hào)識(shí)別中能夠取得更準(zhǔn)確的識(shí)別結(jié)果。

圖3 3種實(shí)驗(yàn)的AUC和ACC結(jié)果對(duì)比圖Fig.3 Comparison of AUC and ACC results of three experiments

4.2 適用性評(píng)價(jià)

由于腦電信號(hào)的個(gè)體差異，同一方法對(duì)不同受試者腦電信號(hào)的分類可能取得不同的效果，因此，分別計(jì)算了3種實(shí)驗(yàn)AUC和ACC的方差如表2所示。

表2 3種實(shí)驗(yàn)AUC和ACC的方差Tab.2 Variance of AUC and ACC in three experiments

由表2可以看出：實(shí)驗(yàn)1的AUC方差為 0.009 151，ACC為0.008 744都遠(yuǎn)遠(yuǎn)小于實(shí)驗(yàn)2和實(shí)驗(yàn)3的方差；由圖3也可以看出實(shí)驗(yàn)1每個(gè)受試者的AUC和ACC結(jié)果差別很小，而實(shí)驗(yàn)2和實(shí)驗(yàn)3每個(gè)受試者的AUC和ACC結(jié)果參差不齊，說明實(shí)驗(yàn)1的方法對(duì)于不同的受試者有更好的適用性。該結(jié)果說明：采用EEMD構(gòu)建腦電信號(hào)高階張量和采用MPCA進(jìn)行降維的特征提取方法相對(duì)于采用復(fù)Morlet小波分解構(gòu)建腦電信號(hào)高階張量和采用MPCA進(jìn)行降維，直接對(duì)采集到的EEG信號(hào)進(jìn)行2DPCA進(jìn)行降維的特征提取方法具有更高的適用性，在對(duì)不同個(gè)體的單次運(yùn)動(dòng)想象腦電信號(hào)識(shí)別中能夠取得穩(wěn)定的識(shí)別效果。

5 討 論

本文提出的采用EEMD構(gòu)建腦電信號(hào)高階張量，MPCA進(jìn)行降維的特征提取方法可以提取出對(duì)EEG信號(hào)識(shí)別最有效的特征。通過100次10折交叉驗(yàn)證取平均值得到的結(jié)果避免隨機(jī)性，保證了分類結(jié)果的準(zhǔn)確性。最終得到本特征提取方法結(jié)合Fisher線性判別分析能達(dá)到93.34%的準(zhǔn)確率，比采用Morlet小波分解構(gòu)建腦電信號(hào)高階張量，MPCA進(jìn)行降維，結(jié)合Fisher線性判別分析分類的方法準(zhǔn)確率高4.75%；比直接對(duì)采集到的EEG信號(hào)進(jìn)行2DPCA進(jìn)行降維，結(jié)合Fisher線性判別分析分類的方法準(zhǔn)確率高2.6%。并且，對(duì)8名受試者的分類準(zhǔn)確率的方差也最小僅為0.008 744，比采用Morlet小波分解構(gòu)建腦電信號(hào)高階張量，MPCA進(jìn)行降維，結(jié)合Fisher線性判別分析分類的方法準(zhǔn)確率方差減小72.69%，比直接對(duì)采集到的EEG信號(hào)進(jìn)行2DPCA進(jìn)行降維，結(jié)合Fisher線性判別分析分類的方法準(zhǔn)確率方差減小23.86%。這說明該方法對(duì)不同受試者的單次運(yùn)動(dòng)想象腦電信號(hào)都能取得很好的識(shí)別效果。

構(gòu)建EEG信號(hào)高階張量時(shí)，Morlet小波分解方法需要調(diào)節(jié)不同的小波參數(shù)來適應(yīng)不同人的EEG信號(hào)才能取得更準(zhǔn)確的分類結(jié)果，因此在不同人身上的識(shí)別準(zhǔn)確率出現(xiàn)很大的差異性，甚至同一個(gè)人的不同次實(shí)驗(yàn)的識(shí)別準(zhǔn)確率也會(huì)有較大的差異。

EMD分解方法由于具有自適應(yīng)性，能夠根據(jù)數(shù)據(jù)本身來分解數(shù)據(jù)，但是會(huì)出現(xiàn)模態(tài)混疊，同一個(gè)人不同次實(shí)驗(yàn)的分解長(zhǎng)度不能對(duì)齊，不能用來構(gòu)建EEG信號(hào)高階張量。而EEMD分解方法在EMD分解信號(hào)中加入白噪聲，繼承了EMD優(yōu)點(diǎn)的同時(shí)，解決了模態(tài)混疊現(xiàn)象，能夠?qū)EG根據(jù)頻率自適應(yīng)分解。因此，在構(gòu)建EEG信號(hào)高階張量時(shí)，EEMD方法要優(yōu)于Morlet小波分解方法。在整個(gè)EEG信號(hào)的分類識(shí)別中，基于EEMD方法構(gòu)建腦電信號(hào)高階張量，結(jié)合MPCA進(jìn)行降維的特征提取方法提取了EEG信號(hào)時(shí)域，頻域和空間域的多模態(tài)特征[19～20]，而直接采用2DPCA的方法忽略了EEG信號(hào)的頻域信息，未能取得更好的識(shí)別效果。

本文將EEMD方法引入EEG信號(hào)高階張量的構(gòu)建，結(jié)合MPCA進(jìn)行特征提取，使得在特征提取的過程中，不需要將原始高階張量展開為向量，保留了腦電信號(hào)的時(shí)域，頻域和空間域信息。并結(jié)合Fisher線性判別分析構(gòu)成的算法組合用于單次運(yùn)動(dòng)想象腦電信號(hào)識(shí)別。采集腦電數(shù)據(jù)數(shù)據(jù)，設(shè)計(jì)對(duì)照實(shí)驗(yàn)，并通過100次10折交叉驗(yàn)證的方式避免了實(shí)驗(yàn)結(jié)果的隨機(jī)性，驗(yàn)證了新提出方法的有效性。

6 結(jié) 論

本研究提出基于EEMD方法構(gòu)建腦電信號(hào)高階張量，結(jié)合MPCA進(jìn)行降維的特征提取方法。通過與基于Morlet小波構(gòu)建腦電信號(hào)高階張量結(jié)合MPCA降維的特征提取方法和直接對(duì)腦電信號(hào)進(jìn)行2DPCA降維的特征提取方法進(jìn)行對(duì)比，結(jié)果表明：

(1) 本研究提出的特征提取方法能夠提取更明顯的特征，結(jié)合Fisher線性判別分析能夠取得更高的分類準(zhǔn)確率；

(2) 本研究提出的特征提取方法相對(duì)于上述其他兩種特征提取方法對(duì)不同的受試者具有更好的適用性，能夠更好的應(yīng)用于不同個(gè)體。