基于VMD與CSP的腦電特征提取方法

劉 帥，烏日￣開西·艾依提

(新疆大學機械工程學院，新疆 烏魯木齊 830047)

1 引言

腦機接口(Brain Computer Interface， BCI)為人類與周圍環境進行交互提供了一種不依賴于外周神經和肌肉的新的輸出通道[1，2]。腦電圖(Electroencephalographic，EEG)具有無創、高時間分辨率、低成本等優點，被廣泛應用于腦機接口中。研究表明，當人在想象左手或右手運動時，大腦的感覺運動皮層Mu節律和Beta節律能量會產生變化，這一現象被稱為事件相關去同步(Event-related Desynchronization， ERD)和事件相關同步(Event-related Synchronization， ERS)[3]。近年來，運動想象腦機接口技術被用來訓練神經肌肉損傷的患者，幫助其恢復肢體的運動感知功能[4，5]。由于運動想象腦電信號是由大腦活動產生的，信噪比低，容易受到外界噪聲的影響，難以獲得較高的識別精度。如何對腦電信號進行有效的特征提取和分類成為一個難點問題。

目前，國內外研究者提出了若干運動想象腦電特征提取和分類算法，常用的特征提取算法有自回歸模型(Auto Regressive，AR)、小波包變換(Wavelet Package Transform，WPD)、功率譜密度(Power Spectral Density，PSD)、共空間模式(Common Spatial Pattern，CSP)等[6-8]。常用的分類算法有支持向量機(Support Vector Machine，SVM)、線性判別分析(Linear Discriminant Analysis，LDA)等[9]。楊默涵等[10]提出了一種基于總體經驗模態分解(Ensemble Empirical Mode Decomposition，EEMD)的運動想象腦電處理方法，將分解后的分量利用希爾伯特變換提取邊際譜及瞬時能譜，選取較優頻帶波幅特征作為特征，取得了較好的分類效果。Shalu Chaudhary等[11]采用兩種時頻分析方法將腦電信號轉化為圖像，輸入到基于AlexNet模型的深度卷積神經網絡(Deep Convolutional Neural Networks，DCNN)進行特征提取與分類，發現連續小波變換(Continuous Wavelet Transform， CWT)轉化的圖像具有更高的識別率。Yu Z等[12]提出了一種時間約束稀疏組空間模式(Temporally Constrained Sparse Group Spatial Patterns，TSGSP)特征提取算法，通過同時優化CSP中的濾波器頻帶和時間窗口提高運動想象腦電信號的分類精度。

為了進一步提高在少數通道下的運動想象腦電信號的分類準確率，本文提出一種變分模態分解(Variational Mode Decomposition，VMD)結合CSP的特征提取方法。將腦電信號進行VMD分解，得到多個具有物理意義的固有模態函數(Intrinsic Mode Function，IMF)，通過巴氏距離計算各階固有模態函數與原始信號之間的相關程度，選取合適的固有模態函數組合成為新的信號矩陣作為輸入，再通過CSP構造空間濾波器，得到各IMF的空間分布特征，最后使用支持向量機進行分類。

2 數據集描述

本文采用兩個數據集，第一個是BCI Competition II的Data set III的數據集[13]，數據集記錄了一位25歲的健康女性，實驗要求被試坐在屏幕前，根據屏幕上隨機出現的左右箭頭提示進行相應的左右手運動想象。數據采用雙導聯方式記錄，采集了C3，Cz和C4三個通道的EEG信號。采樣頻率為128Hz，對信號進行0.5-30Hz濾波。實驗共記錄7組，每組包含40次實驗，每次實驗持續9s，共280次實驗，選其中的140次實驗作為訓練樣本，剩余的140次實驗作為測試樣本。單次實驗的采集范式如圖1所示。

圖1 Data set III單次實驗采集范式

第二個數據集來自三位身體健康的志愿者，均為右利手，且均為第一次參加腦電實驗。使用Neuracle公司的NeuSenW16通道無線腦電采集系統進行實驗數據的采集，腦電帽的電極按照國際標準10-20系統放置，僅采集了C3、C4兩個通道的數據。參考電極位于頭頂部，接地電極位于前額。采樣頻率設置為1000Hz，初始帶通濾波為0.5-100Hz，陷波濾波50Hz。

在實驗過程中，受試者佩戴腦電帽被要求坐在距離電腦屏幕大約1米的椅子上，雙手自然放置在桌子上，保持放松狀態。在想象過程中，要求被試盡量避免眨眼、眼動、身體以及頭部移動。單次實驗數據采集范式如圖2所示。試驗開始前2s為準備階段，屏幕顯示空白；2s后有一段短暫的提示音并且屏幕顯示“+”字形，提示被試實驗馬上開始；4-7s屏幕會隨機出現向左或向右的箭頭，被試需要根據提示進行左右手的運動想象；7s后有2-4s的隨機休息時間；為了防止被試疲勞，每隔20次實驗休息2min。實驗共分為3組，每組100次實驗，共300次實驗，其中200次實驗作為訓練樣本，100次實驗作為測試樣本。

圖2 實驗室單次實驗數據采集范式

3 改進的特征提取方法

因個體差異較大，每位受試者在運動想象時，ERD/ERS發生的頻率段不相同。傳統的共空間模式是根據固定頻率段的腦電信號進行特征提取，且需要大量通道的數據。針對這一問題，本文通過結合VMD和CSP用于改善特征提取算法。改進后的特征提取流程如圖3所示。首先，將濾波后的信號進行VMD分解成K個IMF分量，然后通過計算各IMF分量與原始信號的巴氏距離，選擇巴氏距離最小的3階IMF分量進行重構，采用CSP濾波提取特征，得到特征向量F。

3.1 變分模態分解

變分模態分解通過迭代搜索變分模型的最優解來確定每個模態分量的中心頻率和有限帶寬，可以有效抑制經驗模態分解(Empirical Mode Decomposition，EMD)方法中存在的端點效應和模態混疊現象，具有堅實的理論基礎，且對采樣和噪聲的魯棒性更強[14]。

圖3 VMD與CSP相結合的特征提取流程

在VMD算法中，定義固有模態函數(IMF)為一個調幅-調頻信號，通過混合一個預估中心頻率，將每個模態的頻譜調制到相應的基頻帶

(1)

計算梯度平方L2范數對信號解調，得到各個模態函數的帶寬，則相約束變分問題如下：

(2)

式中，{uk}={u1，u2，…uk}為分解后的k個模態分量，{ωk}={ω1，ω2，…ωk}為每個模態對應的中心頻率，f為輸入的腦電信號。

使用二次懲罰因子α和Lagrange乘法算子λ，將約束變分問題轉變成非約束變分問題，增廣型的拉格朗日表達式如下

L({uk}，{ωk}，λ)

(3)

利用交替方向乘子法(Alternate Direction Method of Multipliers，ADMM)，并結合Parseval/Plancherel 傅里葉等距變換，優化得到各模態分量和中心頻率，求解增廣拉格朗日函數的最優解，交替尋優迭代un+1，ωn+1和λn+1的表達式如下：

(6)

迭代終止條件

(7)

最終得到K個有限帶寬的IMF分量。

3.2 模態函數選擇

利用巴氏距離用于測量IMF分量與原始信號的相似性。對于相同定義域X上的離散概率分布p和q，巴氏距離為

DB(p，q)=-ln(BC(p，q))

(8)

(9)

其中，p表示第K個IMF分量，q表示原始信號，BC(p，q)為離散概率分布 Bhattacharyya系數。

3.3 傳統共空間模式算法

共空間模式是一種對兩分類任務下的空域濾波特征提取算法，能夠從多通道的腦電信號中提取每一類的空間分布成分。共空間模式的原理是利用矩陣對角化，找到一組最優空間濾波器對信號進行投影，使得一類方差最大化，另一類方差最小化，從而得到具有較高區分度的特征向量[15]。

單次運動想象腦電實驗數據由一個N×M的矩陣X組成，其中N表示通道數，M表示采樣點數，其歸一化后的協方差矩陣計算如下

(10)

(11)

對混合協方差矩陣Vc進行特征分解

(12)

式中，Uc為特征向量矩陣，λc為特征值矩陣。構造白化矩陣P

(13)

(14)

式中，Sl和Sr具有共同的特征矢量E，且對應的特征值之和為1，則Sl最大特征值對應的特征向量是Sr最小特征值對應的特征向量。對Sl和Sr進行特征分解可得

Sl=EλlETSr=EλrET

(15)

式中，對角矩陣λl和λr之和為單位矩陣I。由此構造投影矩陣W

W=ETP

(16)

則對輸入矩陣X進行濾波可得Z=WX，取Z的前m列和后m列構成特征矩陣Zp，p=1，2，…，2m。

(17)

3.4 支持向量機

支持向量機是一種有效的二分類模型，具有較強的泛化能力。其基本原理是將輸入的向量映射到高維的特征空間，并在此空間尋找一個最優分類超平面，使兩類樣本的分類間隔最大[16]。本文將改進的CSP特征提取出來的特征向量作為輸入，使用LIBSVM工具包進行分類。訓練集用來訓練模型，測試集用來評估最優超參數訓練出來的模型，同時對模型的參數進行尋優。采用網格搜索算法尋找最佳懲罰系數c和核參數g組合。

4 實驗結果分析

4.1 數據預處理

根據ERD/ERS現象，在做肢體運動想象時，感覺運動區域的節律變化主要發生在Mu節律(8-13Hz)和Beta節律(14-30Hz)，本文使用6階Butterworth濾波器對原始EEG信號進行8-30Hz的濾波。

利用連續小波變換對信號進行時頻分析，通過時頻圖觀察運動想象期間ERD發生的時間段。想象左右手運動時的時頻圖與腦地形圖如圖4和圖5所示。想象左手運動時，在8-12Hz頻段對側的C4通道3.5-7s時間段的能量顯著降低，發生ERD現象；而同側的C3通道3.5-7s時間段的能量顯著增加，發生ERS現象。與之相反，在想象右手運動時，對側的C3通道3.5-7s發生ERD現象，同側的C4通道3.5-7s發生ERS現象。對于BCI Competition II的Data set III的數據集，本文選擇單次實驗的3.5-7s時間段的腦電信號進行后續的特征提取。數據集中每次實驗共采集9s數據，采樣頻率128Hz，每次實驗每個通道共有1152個采樣點。選取其中3.5-7s的數據，此時數據為448×3×140。

圖4 想象左手運動時頻圖與腦地形圖

圖5 想象右手運動時頻圖與腦地形圖

4.2 特征提取與分類

將預處理后的腦電信號進行VMD自適應分解，得到K個IMF分量，每階IMF分量都有各自的中心頻率和有限帶寬。然后計算各IMF分量與原始信號的巴氏距離來選擇合適的IMF分量。本文K取值為7。以想象左手為例，C3通道的腦電信號經VMD分解后的IMF波形及各分量對應的頻譜圖如圖6所示。從頻譜圖中可以看出，自適應分解后的IMF分量主要集中在8-30Hz范圍，不同于EMD分解時出現的相同成分的信息出現在不同IMF分量中，VMD分解后的各階分量均為有限帶寬，有效克服了EMD中的模態混迭現象。

圖6 VMD分解后的各階IMF分量及各分量對應的頻譜圖

根據式(8)計算各IMF分量與原始信號的巴氏距離。以單次實驗數據為例，表1為C3和C4通道經VMD分解后各IMF分量與原始信號的巴氏距離。巴氏距離越小，對應的模態分量與原始信號越相似。觀察表1可知，C3通道巴氏距離最小的三階IMF分量分別為IMF2、IMF3和IMF4；Cz通道巴氏距離最小的三階IMF分量分別為IMF2、IMF1和IMF3；C4通道巴氏距離最小的三階IMF分量分別為IMF2、IMF4和IMF7。

表1 IMF分量與原始信號的巴氏距離

本文選取每個通道數據分解后與原始信號巴氏距離最小的三階IMF分量重新構造新的輸入矩陣作為特征提取的輸入，原輸入矩陣重新構造成一個448×9的矩陣，再通過CSP濾波提取特征得到特征向量，最終將特征向量輸入到支持向量機中進行分類。

4.3 實驗對比

為了證明本文方法的有效性，在相同腦電數據集情況下，與其它文獻的方法進行了比較，對比結果如表2所示。可以看出，在使用同一數據集的情況下，相比于其它三種算法，本文所提方法的分類結果取得了91.43%的分類準確率。說明本文所提出的特征提取算法能夠有效提高運動想象腦電的分類準確率。

表2 本文方法與其它方法分類準確率對比

為了進一步驗證該方法的有效性，在實測數據中進行了驗證，與傳統的特征提取算法CSP進行了比較。如表3所示，對實驗室實測的三名被試的C3和C4通道的腦電信號分別進行CSP特征提取和VMD融合CSP的特征提取，并用支持向量機進行分類。從表3可以看出，相比于傳統的CSP特征提取算法，經VMD分解后的信號再進行CSP濾波提取特征能夠獲得更高的分類識別率，三名被試的分類準確率分別提高了6%、8%、6%。

表3 實測數據的分類識別率比較

5 結論

針對運動想象腦電信號個體差異較大，傳統的CSP特征提取算法需要大量通道的數據的問題，本文提出了一種VMD結合CSP的特征提取算法。通過巴氏距離可以選擇與腦電信號最相關的固有模態函數構造信號矩陣進行特征提取，最后使用SVM分類。在公共數據集上取得了91.43%的分類準確率，在實測數據上相比于傳統的共空間模式特征提取方法，取得了更好的分類效果。

實驗結果表明該方法在少數通道的情況下能夠有效提取運動想象腦電信號的特征，取得了較好的分類準確率，證明了該方法的有效性，為運動想象腦機接口應用奠定了基礎。在下一階段，本文提出的方法將用于實現在線腦機接口。