基于EMD和Hilbert變換的腦磁信號特征提取和分類

王健鋒， 黃曉霞

(上海海事大學 信息工程學院， 上海 201306)

0 引 言

腦機接口(Brain-Computer Interface，BCI)的目的是建立一種不依賴外周神經和肌肉組織的新型通信系統(tǒng).該系統(tǒng)通過測量腦信號，如腦磁圖(magnetoencephalography, MEG)，分析人類的意圖，并將該意圖轉換為外部設備的控制信號.[1-2]因此，BCI可以給外周神經或肌肉組織受損的患者提供一種新的與外部世界交流的方式.BCI分為植入式BCI和非植入式BCI.非植入式BCI具有無創(chuàng)傷性、易于攜帶、價格低廉等優(yōu)點，為BCI的實際應用和推廣提供可能，成為BCI研究領域的一個熱點.在非植入式BCI研究中，腦電圖(electroencephalogram, EEG)和MEG是兩種常用的輸入信號.對EEG信號的研究已經取得很大成就，但是EEG信號的采集容易受到頭蓋骨的影響，導致其空間分辨率不高，誤差很大，使EEG的應用范圍受到限制.MEG不受頭蓋骨的干擾，具有較高的時空分辨率和更低的失真率，能夠直接反映腦內場源的活動狀態(tài)，受到不少BCI研究者的青睞.[3-4]

在2008年第4屆國際BCI競賽中，第3組數(shù)據(jù)使用的就是MEG數(shù)據(jù)，許多研究機構參加這次競賽，其中HAJIPOUR等獲得該組競賽的第一名，他們用統(tǒng)計特征、自回歸模型和小波變換的方法進行特征提取，用有監(jiān)督的遺傳算法進行特征降維，用線性判別分析(Linear Discrimination Analysis, LDA)和支持向量機(Support Vector Machine, SVM)進行分類，所得MEG的平均分類準確率為46.9%.我國燕山大學的李靜等獲得第二名，他們用快速傅里葉變換和主成分分析(Principal Component Analysis, PCA)法進行特征提取，用Fisher線性判別(Fisher Linear Discrimination, FLD)方法進行特征降維，用LDA進行分類，所得MEG的平均分類準確率為25.1%.2011年SABRA等[5]對該MEG進行研究，用統(tǒng)計特征、功率譜密度和離散小波變的方法進行特征提取，分別用LDA和SVM進行分類，獲得33.3%的平均分類準確率.雖然這些特征處理方法能夠獲得不錯的結果，但都沒有兼顧分類準確率和處理方法簡易性，即分類準確率高時，處理方法比較復雜，處理方法簡單時，分類準確率卻不高.

在不明顯降低分類準確率的情況下，兼顧MEG識別率和處理方法的簡易性，本文提出一種基于經驗模式分解(Empirical Mode Decomposition, EMD)和Hilbert變換的MEG特征提取和分類方法.這是一種時頻數(shù)據(jù)的處理方法，非常適用于非線性和非穩(wěn)定性信號的處理，在軸承故障診斷、肌電信號分析和基于EEG的BCI研究方面都能取得不錯的效果，但目前該方法還沒有被應用到基于MEG的BCI研究中.為驗證EMD和Hilbert變換對基于MEG的BCI的處理性能，使用EMD和Hilbert變換對MEG進行特征提取，然后用PCA法對提取到的特征向量進行特征降維，最后用SVM進行分類.

1 數(shù)據(jù)描述

本文選用第4屆國際BCI競賽的第3組數(shù)據(jù).該數(shù)據(jù)記錄兩名健康的右利手實驗者S1和S2在執(zhí)行手腕運動期間的腦部活動信息.采集數(shù)據(jù)時，兩名實驗者坐在帶有MEG采集設備的椅子上，保持放松，把肘部放在枕墊上以防止上臂和肩運動，頭部用小枕頭固定.要求實驗者僅僅使用右利手和手腕從中心位置向與其夾角為90°的4個方向移動一個操縱桿，運動幅度是4.5 cm.在每次實驗中，方向由實驗者自己選擇，目標位置在菱形的4個角上，分別表明左、右、遠離實驗者、朝向實驗者.

該MEG數(shù)據(jù)記錄頭部運動區(qū)域上的10個信道(LC21，LC22，LC23，LC31，LC32，LC41，LC42，RC41，ZC01和ZC02)的信號(包含運動開始前0.4 s至開始后0.6 s的數(shù)據(jù))，并通過頻帶0.5～100 Hz的帶通濾波器以采樣率400 Hz進行重新采樣.本實驗所使用的MEG數(shù)據(jù)包括：實驗者S1的160次訓練數(shù)據(jù)(每個目標方向分別進行40次實驗)和74次測試數(shù)據(jù)；實驗者S2的160次訓練數(shù)據(jù)(每個目標方向分別進行40次實驗)和73次測試數(shù)據(jù).

2 方法描述

2.1 基于EMD和Hilbert變換的特征提取

2.1.1 EMD

EMD主要思想是從復雜信號里分離出有限個本征模函數(shù)(Intrinsic Mode Function，IMF)，這些IMF具有不同的頻率，并且滿足：①其極值點個數(shù)與過零點個數(shù)相同或相差一個；②其上下包絡關于時間軸局部對稱.[6-8]EMD具體分解過程如下.

(1)對任一信號x(t)，首先確定x(t)上所有的極大值點和極小值點，然后用三次樣條函數(shù)將所有極大值點擬合成原數(shù)據(jù)序列的上包絡線，將所有極小值點擬合成原數(shù)據(jù)序列的下包絡線.

(2)計算上下包絡線的均值曲線m1(t)，將原數(shù)據(jù)序列x(t)減去m1(t)可得到一個新的數(shù)據(jù)序列h1(t)，即

(1)

(3)判斷h1(t)是否滿足IMF定義.若滿足，信號x(t)第1個IMFc1(t)=h1(t)；若不滿足，將h1(t)作為原信號，重復上述步驟，執(zhí)行k次，直到h1k(t)滿足IMF的定義為止，即得信號x(t)第1個IMFc1(t)=h1k(t).

(4)將第1個IMFc1(t)從x(t)中分離出來，得到一個去掉高頻分量的剩余信號r1(t)，即

(2)

(5)將r1(t)作為原信號，重復上述步驟，可得到其他IMF，直到剩余信號是一個單調函數(shù)或只有一個極值時篩選過程結束.最終原始信號x(t)可表示為

(3)

式中：ci(t)代表每個IMF；rn(t)代表剩余信號.

2.1.2 Hilbert變換

經過EMD后，對每個IMF進行Hilbert變換，可得到原始信號的Hilbert譜.[9-10]具體步驟如下.

對每個IMF作Hilbert變換得

(4)

(5)

這里省略剩余信號rn(t)，Re表示取實部，其中每個IMF的幅值和相位隨時間改變.式(5)等號右部稱為信號x(t)的Hilbert振幅譜，簡稱Hilbert譜，記為

(6)

該公式將信號的時間-頻率-幅值的三維關系表達出來，其含義是瞬時振幅在頻率-時間上的分布.

對Hilbert譜取平方可得Hilbert能量譜H2(w,t).在此基礎上，定義瞬時能量密度

(7)

很顯然,瞬時能量密度EI是一個關于時間的函數(shù)，反映能量在時間上的波動.

2.2 基于PCA法的特征降維

PCA法是揭示大樣本、多變量數(shù)據(jù)或樣本間內在關系的一種方法，旨在利用降維的思想，把多指標轉化為少數(shù)幾個綜合指標，降低觀測空間的維數(shù)，以獲取最主要的信息.通常把轉化生成的綜合指標稱為主成分，每個主成分都是原始特征的線性組合，并且各個主成分之間互不相關，這就使得主成分比原始特征具有某些更好的性能.[11-12]

在實際的降維過程中，根據(jù)前m個主成分的累積方差貢獻率決定特征向量的維數(shù).當前m個主成分的累積方差貢獻率大于給定的上限值時，可取前m個主成分作為特征向量，達到降維目的.

2.3 基于SVM的分類

SVM是一種基于統(tǒng)計學習理論的機器學習方法，它首先選擇某種非線性映射算法將輸入向量映射到一個高維特征空間，在該高維特征空間中利用結構風險最小化原則和間隔的概念，構造最優(yōu)線性決策函數(shù)，求解最優(yōu)分類超平面，達到分類的目的.[13-14]

3 實驗過程

本文所提出的算法包括預處理、特征提取、特征降維和分類等4個步驟，其實驗流程見圖1.

圖1 算法實驗流程

3.1 預處理過程

本實驗所采用的MEG數(shù)據(jù)在頻帶0.5～7 Hz范圍內包含的運動方向信息最為明顯[15]，故將這10個信道的MEG數(shù)據(jù)通過頻帶為7 Hz的低通濾波器進行濾波處理.

3.2 特征提取過程

使用EMD對預處理后的MEG數(shù)據(jù)進行分解，得到各信道中各樣本的IMF.其中某信道單個樣本經EMD后得到的結果見圖2.從該圖可以看出，經EMD所得的剩余信號逐漸趨于平滑，故棄之不用，只選取經EMD得到的IMF.

圖2 某信道單個樣本經EMD后的結果

用基于Hilbert變換的方法計算每個樣本的IMF的瞬時能量密度，并將其組合為MEG的特征向量，用于區(qū)分MEG的運動方向.將圖2所對應樣本的IMF進行Hilbert變換后，計算其瞬時能量譜，能量隨時間的變化情況見圖3.

圖3 樣本IMF的瞬時能量譜

3.3 特征降維過程

提取MEG數(shù)據(jù)的特征后，得到單個MEG樣本單個通道的瞬時能量譜為398維，每個樣本共有10個通道，若將單個樣本的10個通道的瞬時能量譜組合為該樣本的特征向量，特征向量的維度將高達3 980.為避免復雜的計算且提高分類準確率，在分類之前對特征向量進行特征降維.在本實驗中，使用PCA法將特征向量的維度降為159，然后根據(jù)主成分的累積方差貢獻率的原則取其前26個主成分，并將歸一化處理過的26個主成分作為MEG的特征向量.

3.4 分類過程

本實驗使用基于線性核函數(shù)的SVM對獲得的特征向量進行分類.為獲得較好的分類效果，對SVM的懲罰因子C進行尋優(yōu).圖4是S1和S2的分類準確率隨C的變化情況.從圖4可知，當C達到22時，S1和S2的平均分類準確率達到最高值，為44.3%。

圖4 S1和S2分類準確率隨C的變化

3.5 實驗結果分析

通過本文提出的算法獲得的分類結果：S1和S2的分類準確率以及二者的平均分類準確率分別為36.5%，52.1%和44.3%.而第4屆國際BCI競賽各個競賽小組所獲得的分類結果見表1.其中競賽第一名分別提取MEG的時域特征、頻域特征和時頻特征作為的聯(lián)合特征，競賽第二名分別提取MEG信號的時域特征和頻域特征作為信號的聯(lián)合特征.

表1 第4屆國際BCI競賽結果 %

采用本實驗所使用的算法獲得的測試集的平均分類準確率比競賽第二名所獲得的平均分類準確率高19.2%，雖沒能超過競賽第一名所獲得的分類準確率，但只與競賽第一名所獲得的平均分類準確率相差2.6個百分點.本實驗特征提取算法使用的是EMD和Hilbert變換，只需提取信號的時頻特征，比競賽第一名所使用的特征提取算法簡單且易于實現(xiàn)；本實驗使用PCA法進行特征降維，與競賽第一名所使用的遺傳算法相比，時間復雜度降低很多；本實驗分類器只選用線性SVM，使算法復雜度和時間消耗降低.本算法在奔騰處理器、主頻2.1 GHz和1 G內存的Window XP下處理整個MEG數(shù)據(jù)的運行時間是281 s，而競賽第一名所使用的算法在奔騰處理器、主頻3 GHz和1 G內存的Windows XP下的運行時間是1 043 s，在運行環(huán)境基本相同且計算機性能并不占優(yōu)勢的情況下，本算法的運行時間顯著減少.總體來說，在MEG分類準確率不明顯降低的情況下，本算法的實現(xiàn)比較簡單，且時間性能優(yōu)于競賽第一名所使用的MEG數(shù)據(jù)處理算法.

4 結束語

實驗結果表明，本文提出的基于EMD和Hilbert變換的特征提取和分類算法不僅能夠用于MEG數(shù)據(jù)的分類，而且能夠獲得較高的分類準確率和較優(yōu)的時間性能，綜合性能較好，具有更高的實際應用價值.由于本研究還在進行中，該算法還有許多地方可以改進，下一步的目標是進一步提高該算法的分類準確率.

參考文獻：

[1] SARDOUIE S H, SHAMSOLLAHI M B. Selection of efficient features for discrimination of hand movements from MEG using a BCI competition IV data set[J]. Frontiersin Neuroscience, 2012(6): 42-47.

[2] 高上凱. 基于節(jié)律性腦電信號的腦-機接口[J]. 生命科學, 2008, 20(5): 722-724.

[3] 吳婷, 顏國正, 楊幫華. 基于小波包分解的腦電信號特征提取[J]. 儀器儀表學報, 2008, 28(12): 2230-2234.

[4] KRISHNA S, VINAY K C, RAJA K B. Efficient MEG signal decoding of direction in wrist movement using curve fitting (EMDC)[C]// 2011 Int Conf Image Inform Proc (ICIIP). IEEE, 2011: 1-6.

[5] SABRA N I, WAHED M A. The use of MEG-based brain computer interface for classification of wrist movements in four different directions[C]// 2011 28th Natl Radio Sci Conf (NRSC). IEEE, 2011: 1-7.

[6] LEI Y, LIN J, HE Z,etal. A review on empirical mode decomposition in fault diagnosis of rotating machinery[J]. Mech Systems & Signal Processing, 2013, 35(1): 108-126.

[7] LI Lin, JI Hongbing. Signal feature extraction based on an improved EMD method[J]. Measurement, 2009, 42(5): 796-803.

[8] MENDEZ O, CORTHOUT J, Van HUFFEL S,etal. Automatic screening of obstructive sleep apnea from the ECG based on empirical mode decomposition and wavelet analysis[J]. Physiol Measurement, 2010, 31(3): 273-289.

[9] HUANG N E, WU M L, QU W,etal. Applications of Hilbert-Huang transform to non-stationary financial time series analysis[J]. Appl Stochastic Models in Business & Industry, 2003, 19(3): 245-268.

[10] PANOULAS K I, HADJILEONTIADIS L J, PANAS S M. Hilbert-Huang spectrum as a new field for the identification of EEG event related de-synchronization for BCI applications[C]// 30th Annu Int Conf Eng Med & Biol Soc (EMBS). IEEE, 2008: 3832-3835.

[11] 張穎. 基于PCA的模糊BP網絡建模方法在藻類生長狀態(tài)預測中的應用[J]. 上海海事大學學報, 2011, 32(4): 76-79.

[12] JOHNSTONE I M, LU A Y. On consistency and sparsity for principal components analysis in high dimensions[J]. J Am Stat Assoc, 2009, 104(486): 682-693.

[13] 安博文, 李丹, 龐然. 基于 SVM 分類器的集裝箱箱號識別法[J]. 上海海事大學學報, 2011, 32(1): 25-29.

[14] 馮超, 賀俊吉, 史立. 基于支持向量機的轎車車型識別[J]. 上海海事大學學報, 2011, 32(3): 85-89.

[15] WALDERT S, PREISSL H, DEMANDT E,etal. Hand movement direction decoded from MEG and EEG[J]. J Neuroscience, 2008, 28(4): 1000-1008.