何 群，杜 碩，王煜文，陳曉玲，謝 平

(燕山大學(xué) 河北省測試計量技術(shù)與儀器重點實驗室, 河北 秦皇島 066004)

1 引 言

腦機接口(brain-computer interface, BCI)能夠在人腦和電腦或者是外部設(shè)備之間建立直接通路；借此，人無需身體上的動作就可以將命令和控制指令直接發(fā)送給外部設(shè)備[1]。運動想象(motor imagery, MI)分類研究是BCI技術(shù)的重要組成部分; 但由于MI腦電信號的非平穩(wěn)、非線性以及低信噪比的特點[2,3]，其特征提取與分類識別的效果均不穩(wěn)定，傳統(tǒng)人工確定最優(yōu)時段及最優(yōu)頻段的方法又難免造成信息遺漏進(jìn)而導(dǎo)致識別率的降低；因此基于腦電信號的MI分類研究成為難點問題。

迄今為止，研究者已經(jīng)基于受試者進(jìn)行MI時在mu頻帶以及beta頻帶所發(fā)生的事件相關(guān)去同步/事件相關(guān)同步(event-related desynchronization/event-related synchronization, ERD/ERS)[4]現(xiàn)象, 提出了利用時-頻特征、共空間模式(common spatial pattern, CSP)以及非線性動力學(xué)特征在內(nèi)的多種特征提取方法[5～7]用于MI腦電信號的特征提?。?同時也將支持向量機(support vector machine, SVM)、線性判別分析(linear discriminant analysis, LDA)和logistic分類器[8～10]等多種機器學(xué)習(xí)方法用于MI分類中。 Bashar S K等[11]利用多元經(jīng)驗?zāi)B(tài)分解(multivariate empirical decomposition, MEMD)將多通道腦電信號分解后，選取最優(yōu)模態(tài)進(jìn)行短時傅里葉變換(short-time Fourier transform, STFT)，將頻譜峰值作為MI特征輸入k-最近鄰(knearest neighbor, KNN)分類器對MI信號進(jìn)行分類。李明愛[12]等利用主成分分析(principal component analysis, PCA)進(jìn)行特征提取，以IGHSOM神經(jīng)網(wǎng)絡(luò)為分類器識別想象模式。此外，有學(xué)者將深度網(wǎng)絡(luò)運用于MI分類問題。Ren Y F[13]等利用CSP算法得到特征，并以此作為輸入來訓(xùn)練一個4層的深度信念網(wǎng)絡(luò)(deep belief network, DBN)，得到MI的分類結(jié)果。 Lu N等[14]獲得腦電信號的頻域表示, 再以此作為輸入來訓(xùn)練一個4層的頻率深度信念網(wǎng)絡(luò)(frequency deep belief network, FDBN)得到MI的分類結(jié)果。Tabar Y R等人[15]使用卷積神經(jīng)網(wǎng)絡(luò)(convolution neural network, CNN)提取特征后使用堆疊自編碼器(stacked autoencoders, SAE)進(jìn)行MI分類。針對MI分類問題，現(xiàn)有的方法已經(jīng)取得了一定的效果，但腦電信號的識別率依然是限制MI腦機接口發(fā)展的瓶頸；其根源在于現(xiàn)有特征提取及模式識別方法難以兼顧受試者個體差異性，同時手動提取最優(yōu)想象時段以及最優(yōu)頻段的方法難以避免信息遺漏導(dǎo)致的識別率降低問題。

變分模態(tài)分解(variational mode decomposition, VMD)[16]可將腦電信號分解為多個窄帶分量，各分量包含了信號不同時-頻尺度上的特征，同時解決了模態(tài)混疊的問題。深度信念網(wǎng)絡(luò)(deep brief network, DBN)具有自動學(xué)習(xí)有效特征的能力[17]，從而達(dá)到特征降維的目的，已有研究者基于這一特性將其應(yīng)用到睡眠分期、情緒識別和癲癇診斷等腦電信號處理領(lǐng)域[18]。

DBN能夠?qū)Ω呔S特征降維得到MI的最優(yōu)特征，從而避免手動確定最優(yōu)想象時段及最優(yōu)頻段導(dǎo)致的識別率降低問題。因此，本文提出使用VMD對腦電信號分解；利用希爾伯特變換提取邊際譜(marginal spectrum, MS)、特征頻帶下的瞬時能譜(instantaneous energy spectrum, IES)以及時-頻聯(lián)合特征(joint time-frequency features, JT-F)來涵蓋完整MI時程下的時頻域特征，并融合3種特征；引入DBN對融合后的高維特征降維得到MI的最優(yōu)特征并實現(xiàn)MI腦電模式的識別，避免了手動確定最優(yōu)時段及最優(yōu)頻段導(dǎo)致的識別率降低問題。

2 基于VMD的高維時頻特征提取

由于腦電信號具有非線性、非平穩(wěn)和低信噪比的特點，經(jīng)驗?zāi)B(tài)分解(empirical mode decomposition, EMD)及局部均值分解(local mean decomposition, LMD)等自適應(yīng)分解方法被廣泛應(yīng)用于腦電信號特征的提取[19, 20]。

但現(xiàn)有的自適應(yīng)分解框架存在多種缺陷[16]：(1)缺乏數(shù)學(xué)理論支撐；(2)多數(shù)分解方法基于遞歸篩選，誤差因多次遞歸被放大，進(jìn)而出現(xiàn)模態(tài)混疊現(xiàn)象；(3)缺乏噪聲魯棒性。

VMD是一種新的自適應(yīng)分解方法，與EMD及LMD等遞歸式分解方式不同；VMD假設(shè)每個模態(tài)為帶寬有限，中心頻率不同的窄帶分量，將模態(tài)的估計問題轉(zhuǎn)化為求解變分問題，改善了EMD及LMD等方法的模態(tài)混疊問題。

VMD算法中，定義本征模態(tài)分解函數(shù)(intrinsic mode function，IMF)為調(diào)幅-調(diào)頻信號，記為：

uk(t)=Ak(t)cos(φk(t))

(1)

(2)

式中：{u1,…,uK}及{ω1,…,ωK}分別為分解得到的K個模態(tài)分量及每個模態(tài)對應(yīng)的中心頻率。

式(2)引入增廣拉格朗日函數(shù)得到式(3)：

(3)

式中：α為懲罰參數(shù)；λ為拉格朗日乘子。

(4)

采用Parseval/Plancherel傅里葉等距變換將式(4)轉(zhuǎn)化到頻域中，得到各模態(tài)分量及中心頻率的的頻域更新分別為：

(5)

(6)

同時利用式(7)更新：

(7)

本文首先對MI腦電信號進(jìn)行8～30 Hz帶通濾波，然后對濾波后的腦電信號進(jìn)行VMD得到K個IMF分量。

為提取運動想象腦電信號的時頻特征，對分解得到的各IMF分量進(jìn)行希爾伯特變換：

(8)

Zk(t)=uk(t)+jYk(t)=ak(t)ejθk(t)

(9)

式中： PV是柯西主值；t′為時域中的某一時刻；ak(t)是瞬時幅值;θk(t)是瞬時相位。則瞬時頻率定義為：

(10)

則原始腦電信號可以表示為各IMF分量的和值形式：

(11)

定義希爾伯特譜為：

(12)

本文按照式(13)～式(15)分別求取希爾伯特邊際譜(MS)，特征頻帶下的瞬時能譜(IES)及時-頻聯(lián)合特征(JT-F)。

(13)

(14)

(15)

對于式(13)，為避免人工選擇造成的時、頻域信息遺漏，將式中的T定為運動想象結(jié)束時刻，由于C3與C4通道分別分布在大腦的左側(cè)和右側(cè)，CZ通道則分布在大腦的正中央，所以受試者在進(jìn)行左右手MI實驗時，僅C3與C4通道腦電信號能夠體現(xiàn)出ERS/ERD現(xiàn)象。因此，本文在進(jìn)行特征提取時，不對CZ通道腦電信號進(jìn)行處理, 可以得到完整MI時程下C3、C4通道的邊際譜特征F1：

F1={MSC3,MSC4}

(16)

式(14)中的ω1及ω2分別對應(yīng)MI的mu及beta特征頻帶，分別取兩組對應(yīng)值，其中mu頻帶下的ωmu1與ωmu2分別為8 Hz和13 Hz，beta頻帶下的ωbeta1與ωbeta2分別為16 Hz和28 Hz，則得到兩個特征頻帶下的瞬時能量特征IESmu及IESbeta，至此得到特征頻帶下C3、C4通道的瞬時能量特征F2為：

F2={IESmuC3,IESmuC4,IESbetaC3,IESbetaC4}

(17)

計算JT-F特征，則首先對時間及頻率分別加窗，定義twindow1及twindow2分別為時間窗的起止時間，ωwindow1及ωwindow2分別為頻率窗的起止頻率。本文定義0.25 s時間窗將完整MI時程分割為n個時段，以2 Hz頻率窗將8～30 Hz的特征頻帶分割為11個子頻帶，然后以式(15)計算第i(≤n)個時段第j(j≤11)個頻段下的JT-Fi,j值，則JT-F特征可以表示為：

J={JT-F1,1,JT-F2,1,…JT-Fn,1,…JT-F1,11,

JT-F2,11,…JT-Fn,11}

(18)

則C3、C4通道的時頻聯(lián)合特征F3為：

F3={JC3,JC4}

(19)

融合3種特征形成高維特征向量F={F1,F2,F3}作為MI特征。

3 基于DBN的高維特征分類

DBN由底層的若干個限制玻爾茲曼機(restricted boltzmann machines，RBMS)和頂層的分類器組成，相鄰2層即為一個RBM,各層內(nèi)部相互獨立，但數(shù)據(jù)可通過特定的激活函數(shù)按照RBM的學(xué)習(xí)規(guī)則在層間相互轉(zhuǎn)換，圖1為5層DBN網(wǎng)絡(luò)結(jié)構(gòu)圖。

圖1 DBN結(jié)構(gòu)圖Fig.1 The structure of DBN

圖1中，每個RBM由可視層(v)和隱含層(h)組成，層與層之間通過權(quán)重w彼此互聯(lián)，但層內(nèi)無連接，假設(shè)v層包括n個可見單元，h層包括m個隱含單元。則1組確定狀態(tài)的RBM系統(tǒng)的能量為：

(20)

式中：vi為第i個可見單元狀態(tài)；hj為第j個隱含單元狀態(tài)；對?i,j,vi∈{0,1},hj∈{0,1}，θ={wij,aj,bj}為RBM模型的參數(shù)；wij為可見單元i與隱含單元j間的連接權(quán)重；ai為可見單元i的偏置；bj為隱含單元j的偏置。定義可視層和隱含層的聯(lián)合概率密度為：

P(v,h;θ)=e-E(v,h;θ)/Z(θ)

(21)

(22)

由式(22)，隱含節(jié)點與可視節(jié)點被激活的概率為：

(23)

式中:σ為激活函數(shù)。

為使RBM良好擬合訓(xùn)練樣本，將訓(xùn)練目標(biāo)定為最大化似然函數(shù)：

(24)

Hinton利用對比散度(contrastive divergence，CD)方法，只需少數(shù)次狀態(tài)轉(zhuǎn)移就可以得到對可視層的較好估計。該算法首先計算所有隱含層單元條件概率，并使用Gibbs抽樣確定隱含單元狀態(tài)；再計算可視層狀態(tài)，完成可視層的重構(gòu)，并求解式(24)。

最后得出權(quán)重、偏置的更新過程為：

(25)

式中:〈·〉data為原始數(shù)據(jù)的期望；〈·〉recon為重構(gòu)數(shù)據(jù)的期望；ε為學(xué)習(xí)率。針對參數(shù)ε，Hinton提出學(xué)習(xí)動量項參數(shù)Momentum，使收斂更為精準(zhǔn)，Momentum定義為：

(26)

式中:ρ為動量學(xué)習(xí)率。

則權(quán)重偏置的更新法則變?yōu)椋?/p>

(27)

在圖1所示的DBN中，上一個RBM的隱含層即為下一個RBM的可視層，上一個RBM的輸出即為下一個RBM的輸入，輸入數(shù)據(jù)經(jīng)過低層RBM學(xué)習(xí)后的輸出結(jié)果作為高一層RBM的輸入，直至最后一層；因此高層RBM能夠?qū)W習(xí)到更抽象，更稀疏的特征表示。

以上過程為前向的無監(jiān)督學(xué)習(xí)過程，利用數(shù)據(jù)的標(biāo)簽，采用反向傳播方法從網(wǎng)絡(luò)的softmax分類器開始對網(wǎng)絡(luò)參數(shù)進(jìn)行反向微調(diào)，可以進(jìn)一步優(yōu)化DBN網(wǎng)絡(luò)參數(shù)，避免網(wǎng)絡(luò)陷入局部最優(yōu)。

針對運動想象高維細(xì)節(jié)特征F，本文構(gòu)建4層DBN模型對其進(jìn)行特征降維得到運動想象的最優(yōu)特征，進(jìn)而實現(xiàn)對受試者左右手MI意圖的識別。至此，本文完成了如圖2所示的基于VMD與DBN的腦電信號特征提取與識別框架。

圖2中，利用DBN中的4層網(wǎng)絡(luò)結(jié)構(gòu)對高維特征F降維得到其稀疏表達(dá)，降低MI高維特征中的噪聲、非線性和非相關(guān)信息，達(dá)到特征降維的目的，從而達(dá)到更好的分類效果。

圖2 基于VMD與DBN的特征提取與識別框架Fig.2 The framework of feature extraction and recognition based on VMD and DBN

4 實驗研究與結(jié)果分析

4.1 受試對象與實驗流程

為驗證本文提出的方法的普適性及有效性，實驗采用BCI 2003競賽Dataset Ⅲ數(shù)據(jù)集[21](S1)、BCI 2005競賽Dataset Ⅲb數(shù)據(jù)集[22](S2～S4)以及BCI 2008競賽Dataset IIb數(shù)據(jù)集[23](S5～S7)。

各數(shù)據(jù)集信息如表1所示。

表1 數(shù)據(jù)集信息Tab.1 The information of datasets

信號采集過程中，由于采集設(shè)備原因，個別時刻信號未被成功采集，這些數(shù)據(jù)點在信號中表現(xiàn)為“NAN”點，為提高信號信噪比，實驗中對數(shù)據(jù)進(jìn)行預(yù)處理工作，對數(shù)據(jù)中的“NAN”數(shù)據(jù)點置零，同時為減小數(shù)據(jù)處理的工作量，對受試者S5～S7對應(yīng)數(shù)據(jù)進(jìn)行了降采樣處理，將采樣頻率降為125 Hz。由于運動想象的ERD/ERS現(xiàn)象主要出現(xiàn)在mu節(jié)律(8～13 Hz)與beta節(jié)律(16～28 Hz)，因此對原始信號進(jìn)行8～30 Hz帶通濾波，所用濾波器為6階巴特沃斯濾波器，設(shè)置阻帶截止頻率分別為6 Hz和32 Hz。

實驗中將VMD分解模態(tài)數(shù)K定為7，選擇sigmoid函數(shù)作為DBN網(wǎng)絡(luò)的激活函數(shù)，隱藏層神經(jīng)元個數(shù)定為[400，200]，動量學(xué)習(xí)率定為0.1，批處理容量為全體訓(xùn)練集，各受試者學(xué)習(xí)率、預(yù)訓(xùn)練次數(shù)及微調(diào)次數(shù)等信息如表2所示。

4.2 實驗結(jié)果與分析

4.2.1 特征組合及對比

為比較不同特征提取方法對識別率的影響，實驗設(shè)計3個參照組與所提方法進(jìn)行對比，分別為：

(1)采用4層DBN網(wǎng)絡(luò)對邊際譜特征F1降維并識別MI腦電信號；(2)采用4層DBN網(wǎng)絡(luò)對瞬時能量譜特征F2降維并識別MI腦電信號；(3)采用4層DBN網(wǎng)絡(luò)對時-頻聯(lián)合特征F3降維并識別MI腦電信號。表3呈現(xiàn)了各受試者在3個參照組與本文所提方法下的最大識別率。

表2 各受試者參數(shù)信息Tab.2 The parameters information of subjects

表3 不同特征下的最大識別率Tab.3 The maximum recognition rates of different features (%)

由表3中發(fā)現(xiàn)，采用邊際譜特征所得平均最大識別率為60.51%，采用瞬時能量譜特征所得平均最大識別率為73.03%，采用時-頻聯(lián)合特征所得平均最大識別率為77.59%，采用融合特征所得平均最大識別率為79.00%。

其中，邊際譜特征取得的分類效果相對較差，這是因為實驗范式的設(shè)定導(dǎo)致MI信號強弱隨時間變化，完整MI時程下的邊際譜特征的時域分辨率較低，不能準(zhǔn)確表征MI，但邊際譜特征的頻域分辨率高，對最后的識別結(jié)果有輔助作用；與邊際譜特征相反，特征頻帶下的瞬時能譜則具有優(yōu)秀的時域分辨率和較差的頻域分辨率，取得了更加良好的識別效果；由于兼顧了時、頻分辨率，3種單一特征中，時-頻聯(lián)合特征取得了最優(yōu)的識別效果；而將3種特征融合后，由于時、頻分辨率均較高，MI的識別率得到提升，同時不同受試者識別率間的標(biāo)準(zhǔn)差降低，分類框架的穩(wěn)定性得到了進(jìn)一步提升。

4.2.2 不同分解方法對比

實驗對預(yù)處理后的腦電信號進(jìn)行VMD分解，S1的C3通道腦電信號VMD分解結(jié)果及其分量頻譜圖如圖3。

實驗中作為對比，將腦電信號進(jìn)行EMD分解，S1的C3通道腦電信號EMD分解結(jié)果及分量頻譜圖如圖4。

在圖3及圖4中，由于對腦電信號預(yù)先進(jìn)行了濾波處理，自適應(yīng)分解的IMF分量頻率集中在8～30 Hz范圍內(nèi)。由圖3及圖4可以看出，EMD分解得到的分量頻帶混疊現(xiàn)象嚴(yán)重，而VMD的分解結(jié)果均為窄帶分量，不會出現(xiàn)相同成分信息出現(xiàn)在不同分量的情況。實驗中發(fā)現(xiàn)EMD分解所得結(jié)果，特征頻帶主要分布在前3階分量中。

圖3 VMD分解結(jié)果及各分量頻譜圖Fig.3 The decomposition results of VMDand spectral diagrams of components

圖4 EMD分解結(jié)果及各分量頻譜圖Fig.4 The decomposition results of EMD and spectral diagrams of components

本文選取EMD分解的前3階分量，提取相同高維特征，利用DBN對2種方法所得特征進(jìn)行特征降維與分類，分別對2種分解方法所得識別率進(jìn)行比較，識別結(jié)果如表4所示。表4中可以看出，由于更優(yōu)的分解效果，相比于用EMD提取的特征，VMD所提取的特征具有更佳的識別率。

表4 EMD與VMD所得識別率比較Tab.4 Comparison of recognition rates between EMD and VMD (%)

4.2.3 不同方法對比

基于相同的腦電數(shù)據(jù)，將本文所述方法與其它方法進(jìn)行比對，表5及表6給出了比較結(jié)果。

表5中，空白處表示文獻(xiàn)[13]未采用此位受試者的MI腦電數(shù)據(jù)驗證所提方法。表5中，文獻(xiàn)[24]選用功率譜密度、Hjorth參數(shù)、AR模型系數(shù)和小波系數(shù)等多種特征對運動想象分類。文獻(xiàn)[25]利用總體經(jīng)驗?zāi)B(tài)分解(ensemble empirical mode decomposition, EEMD)提取邊際譜、瞬時能譜及近似熵作為特征輸入線性判別分類器識別運動想象腦電信號。

表6中，文獻(xiàn)[26]采用濾波器組共空間模式(filter bank common spatial pattem, FBCSP)提取MI特征，利用多分類器降低訓(xùn)練集與測試集之間的分布差異，提升了MI識別率；文獻(xiàn)[27]則利用新型相空間重構(gòu)算法計算非線性動力學(xué)特征，選用集成分類器識別想象類別。此類人工確定MI最優(yōu)時段及最優(yōu)頻段的方法難免遺漏細(xì)節(jié)信息造成識別率的降低，本文利用DBN網(wǎng)絡(luò)對高維時、頻特征降維得到最優(yōu)特征并識別MI模式，明顯提高了分類準(zhǔn)確率，且S2、S4及S7的識別率提升明顯。

文獻(xiàn)[13]同樣將深度框架應(yīng)用于MI問題中，利用CSP算法提取特征，并構(gòu)建2層DBN網(wǎng)絡(luò)對特征降維并識別腦電模式，從表5的比較可以看出本文所構(gòu)建的方法在特征提取與識別框架方面優(yōu)于文獻(xiàn)[13]所提的方法。

表5 S1～S4所得識別率比較Tab.5 The comparison of recognition rates of S1～S4 (%)

表6 S5～S7所得識別率比較Tab.6 The comparison of recognition rates of S5～S7 (%)

BCI 2003競賽Dataset Ⅲ數(shù)據(jù)集中，組織者將互信息與識別率同時定為判斷指標(biāo)；BCI 2005競賽Dataset Ⅲb數(shù)據(jù)集中，組織者將互信息作為判斷指標(biāo)；BCI 2008競賽Dataset IIb數(shù)據(jù)集中，組織者將Kappa值與識別率作為判斷指標(biāo)；實驗中將本文方法所得分類結(jié)果與BCI競賽冠軍指標(biāo)進(jìn)行對應(yīng)比較，比較結(jié)果見表7。表7中空白處表示此屆比賽組織者未采取此項指標(biāo)作為評判標(biāo)準(zhǔn)。

表7 BCI競賽冠軍結(jié)果與所提方法結(jié)果比較Tab.7 Comparison of results between the BCI competition’s winner and the proposed method

從表7中可以看出：對于S2、S4、S5～S7，本文所述方法所得指標(biāo)均優(yōu)于BCI冠軍所得結(jié)果；對于S1，本文方法所得互信息低于BCI冠軍結(jié)果，但所得識別率優(yōu)于BCI冠軍結(jié)果。

5 結(jié) 論

本文應(yīng)用VMD方法對腦電信號分解，利用希爾伯特變換提取邊際譜、特征頻帶下的瞬時能譜和時-頻聯(lián)合特征表征MI，引入4層DBN網(wǎng)絡(luò)對融合后的高維特征降維得到最優(yōu)特征并分類，避免了人工確定最優(yōu)時段及最優(yōu)頻段導(dǎo)致的信息遺漏，提升了MI腦電信號的識別率。運用所述方法采用腦機接口競賽數(shù)據(jù)集進(jìn)行分類實驗，平均準(zhǔn)確率達(dá)到79%。本文方法有效提升了MI腦電信號的識別率，為MI腦機接口的應(yīng)用奠定了基礎(chǔ)。