基于DWT、MEMD和模糊熵的腦電信號特征提取與分類研究

陳倩倩，徐健，劉秀平，黃磊，惠楠

（西安工程大學(xué) 電子信息學(xué)院，陜西 西安 710048）

0 引 言

腦機(jī)接口［1］（brain-computer interface，BCI）是一種不依賴于外周神經(jīng)和肌肉輸出通道的交流與控制系統(tǒng)，它將腦電信息（electroencephalo grahy，EEG）通過特定的傳感器和計算機(jī)進(jìn)行分析處理，實現(xiàn)腦意識與人的交流以及對外部設(shè)備的控制，為神經(jīng)或肌肉受損的病人提供嶄新的康復(fù)治療方法，使他們不需要依賴別人而獨立完成自己的康復(fù)訓(xùn)練。隨著BCI技術(shù)的發(fā)展，癱瘓病人可以通過一些依賴于腦機(jī)接口技術(shù)而生產(chǎn)的電子設(shè)備，如神經(jīng)假肢、機(jī)械手臂等，實現(xiàn)運動恢復(fù)。另外，腦機(jī)接口技術(shù)在生活娛樂、軍事、人工智能、監(jiān)控等領(lǐng)域均有不斷增長的應(yīng)用價值。

BCI技術(shù)主要步驟包括信號采集、信號“解讀”（預(yù)處理、特征提取、特征分類）和設(shè)備控制。在對腦電信號的“解讀”中，基于特征提取的腦電信號可以對刺激目標(biāo)或者運動想象類別進(jìn)行判斷，實現(xiàn)人的“意念”直接控制外部設(shè)備。因此，特征提取在腦機(jī)接口的研究中得到了廣泛關(guān)注。

腦電信號的特征提取是分類關(guān)鍵。隨著腦機(jī)接口技術(shù)的不斷深入研究，學(xué)者們提出了多種腦電信號特征提取方法。常見的特征提取方法有短時 傅 立 葉 變 換［2］（short-time fourier transform，STFT）、小波分解［3］（wavelet decomposition，WT）、自回歸（auto-regression，AR）模型［4］和共空間模式［5］（common spatial patterns，CSP）等，其中STFT和WT的本質(zhì)是基于傅立葉變換，根據(jù)Heisenberg uncertainty principle可知，該方法不能很好地同時獲得時頻域信息。

近年來，Hibert-Huang變換（HHT）［6］得到了廣泛應(yīng)用，它非常適用于對非線性非平穩(wěn)信號的分析。李昕等［7］結(jié)合小波分析與經(jīng)驗?zāi)B(tài)分解，首先使用小波分解策略提取不同腦區(qū)的目標(biāo)頻段腦電信號，其次使用經(jīng)驗?zāi)B(tài)分解方法提取特征值；鄭佳佳等［8］采用總體經(jīng)驗?zāi)B(tài)分解和Hilbert-Huang變換結(jié)合小波包分析對含噪腦電信號進(jìn)行處理，獲得了更高的信噪比。

腦電信號來源復(fù)雜且不規(guī)則，難以使用完全確定的系統(tǒng)對腦電信號進(jìn)行分析。已有研究表明，大腦是一個非線性動力學(xué)系統(tǒng)。目前研究多通過估計信號的一些非線性動態(tài)參數(shù)（如相關(guān)維數(shù)、Coriolis entropy、李亞普諾夫指數(shù)）估計信號穩(wěn)定性，也有一些研究采用基于卡爾曼平滑算法的非平穩(wěn)線性判別估計非平穩(wěn)參數(shù)，獲取常規(guī)方法無法得到的信息，其中，復(fù)雜度方法主要包括近似熵［9］、樣本熵［10］等。黃 瑞 海 等［11］觀 察 受 試 者 癲癇發(fā)作期間關(guān)聯(lián)維數(shù)、Lyapunov指標(biāo)和近似熵3個非線性動力學(xué)指標(biāo)的變化，證明腦電節(jié)律波的非線性動力學(xué)指標(biāo)可以用于癲癇發(fā)病的腦電鑒別；H.Amin等［12］記錄了18名受試者右手握拳想象運動的腦電圖，然后基于互信息近似和偽近鄰相空間重構(gòu)方法估計出最大Lyapunov指數(shù)，發(fā)現(xiàn)不同想象狀態(tài)的Lyapunov特征變化顯著，為腦機(jī)接口應(yīng)用奠定了研究基礎(chǔ)。

近似熵和樣本熵只對于小樣本效果顯著，而且需要一個匹配的模板，針對以上問題，在實驗樣本熵的理論基礎(chǔ)上，研究者引入美國學(xué)者L.A.Chen的模糊集合理論，得出一種新的方法——模糊熵（fuzzy entropy，F(xiàn)uzzyEn），這種新方法在實驗輸入樣本數(shù)據(jù)模糊的情況下也能建立輸入與輸出之間的對應(yīng)關(guān)系，因此可應(yīng)用在復(fù)雜的非線性非平穩(wěn)的腦電信號研究領(lǐng)域中。

本文提出一種基于離散小波變換（DWT）、多變量經(jīng)驗?zāi)B(tài)分解（MEMD）和模糊熵的特征提取與分類方法。該方法不僅可以解決分解時出現(xiàn)的頻帶過寬現(xiàn)象，而且還可通過最大互信息系數(shù)篩選，提高后續(xù)分類精度。首先對原始腦電信號利用小波變換進(jìn)行分解，得到一系列窄帶信號，利用MEMD分解得到更集中的信號序列和頻帶信號，按照最大互信息系數(shù)的大小，篩選有效的頻帶信號，重構(gòu)所選分量，然后將重構(gòu)信號用模糊熵進(jìn)行特征提取，最后用支持向量機(jī)進(jìn)行分類。

1 基本原理

1.1 小波變換

使用一組函數(shù)不斷靠近或逼近原始數(shù)據(jù)信號或擬合函數(shù)的過程就是小波變換。腦電信號具有長時低頻和短時高頻的特點，它是一種不穩(wěn)定的數(shù)據(jù)信號，所以存在于局部時頻域信號中，而小波可以實現(xiàn)腦電信號處于時域和頻域局部中，即腦電信號在高頻域時時間分辨率較優(yōu)，在低頻域時頻率分辨率較優(yōu)。引入窗函數(shù)為

連續(xù)小波變換信號f（t）定義為

式中：a為尺度位移；b為時間位移。

小波變換的效率和效果是由小波基函數(shù)決定的，最常用的小波基函數(shù)有4種，即db小波、Haar小波、Molet小波和Meyer小波。

腦電信號是離散信號，離散小波需要進(jìn)行小波變換，與連續(xù)小波相比，離散小波是將小波基函數(shù)ψ（a，τ）中的a，τ限制在離散點上。由定義知，DWT基函數(shù)為，其中，j∈z，k∈z，j和k分別為頻率分辨率和時間平移量，則離散小波變換為

設(shè)信號f（t）的采樣頻率為fs，根據(jù)采樣定理，信號的整個頻帶被劃分為多個子頻帶，利用Mallat算法進(jìn)行S層分解，各個子頻帶依次為如圖1所示。

圖1 信號三級分解過程Fig.1 Three-step signal decomposition process

1.2 多變量經(jīng)驗?zāi)B(tài)分解

將腦電信號分解為各種不同IMF的過程稱為經(jīng)驗?zāi)B(tài)分解（EMD），EMD分解基函數(shù)由腦電信號決定，具有良好的適應(yīng)性，能準(zhǔn)確地表達(dá)信號。但由于EMD存在模式混疊的問題，只能處理一維信號，N.U.Rehman等［18］在2010年提出了多變量經(jīng)驗?zāi)B(tài)分解算法（multivariate empirical mode decomposition，MEMD），可同時分解多通道數(shù)據(jù)。IMFs必須滿足兩個條件：一是信號幅值的最大值等于整個信號時域的過0點數(shù)，或者最大差值為1；二是對上下包絡(luò)線的均值進(jìn)行約束，并設(shè)置它們的值為0。

MEMD的具體步驟如下。

步驟1設(shè)一個N維向量組序列v（t），采用Hammersley序列采樣法，在一個n-1維球面上，選擇合適的n-1個方向向量dθk。

步驟2計算原始信號v（t）在每個方向向量dθk上的投影

步驟3確定所有映射信號極值對應(yīng)的瞬時時刻，其中j為極值點位置，j∈[1，L]。

步驟4采用多元樣條插值函數(shù)獲得K個包絡(luò)線

步驟5計算球面的k個方向向量包絡(luò)線均值

步驟6定義固有模態(tài)函數(shù)h（t）=f（t）-m（t），驗證h（t）是否滿足本征模函數(shù)的條件，如果滿足，則f（t）-h（t）的結(jié)果為IMF分量，繼續(xù)步驟（2）～（6）；否則，返回步驟（2）。

當(dāng)結(jié)果信號小于兩個極值點時，保留殘差分量r（t），分解結(jié)束。將原始信號分解為n個IMFs（用hi（t）表示），n個殘差分量（用r（t）表示），此時原始信號為

式中：p為IMF數(shù)量，n元信號的每一元IMF分量在n個通道中按頻率尺度對齊，形成多元通道。

原始信號經(jīng)過MEMD分解，得到多個IMF分量，采用最大互信息系數(shù)篩選出有用的IMFs，這對后續(xù)信號處理至關(guān)重要。

最大互信息系數(shù)（maximal information coefficient，MIC）可衡量兩個事物之間的關(guān)聯(lián)程度或非線性關(guān)系，MIC適用于各種函數(shù)間的關(guān)聯(lián)關(guān)系，且在樣本足夠大時能夠給出不同類型噪聲的相關(guān)關(guān)系相近系數(shù)，MIC算法描述兩個變量在離散二維空間中的相關(guān)性，用散點圖表示，假設(shè)給定隨機(jī)變量X，Y，則對應(yīng)的MIC為

式中：a，b為X，Y在m×n網(wǎng)格上的隨機(jī)格子數(shù)；B=f（data_size）=n0.6；I（x；y）為變量X，Y的互信息。

對于多通道腦電信號，X和Y分別為多導(dǎo)聯(lián)腦電信號生成的矩陣，利用最大互信息系數(shù)篩選含有有效信息的IMFs信號分量。

1.3 模糊熵

模糊熵在樣本熵的基礎(chǔ)上加入模糊因子，使其可以根據(jù)樣本輸入與輸出之間的關(guān)系判斷樣本所屬類別，具有很強(qiáng)的自適應(yīng)性，對于非線性復(fù)雜腦電信號的特征提取具有很好的效果。

模糊熵的計算步驟如下。

對于給定的N維時間序列［u（1），u（2），…，u（N）］，

定義相空間維數(shù)m（m≤N-2），重構(gòu)相空間，

其中，u0（i）為均值，

引入模糊隸屬函數(shù)

定義X（i）和X（j）之間的相關(guān)度

針對每個i，求其平均值，可得

定義模糊度相似度函數(shù)為

原時間序列的模糊熵

得到模糊熵估計值

式中：N為采樣序列長度；m為重構(gòu)維數(shù)；r為相似容限度，r一般為（0.1～0.25）Sd，Sd為采樣序列標(biāo)準(zhǔn)差。

模糊熵大小與參數(shù)的取值關(guān)系密切，為獲得更好的條件概率估計，本文取m=n=2。

1.4 支持向量機(jī)

支持向量機(jī)（SVM）是二分類模型，即找到一個能區(qū)分兩類的超平面。為了得到最優(yōu)的分類效果，通過最大間隔分類器選擇最優(yōu)的超平面，不僅可以將兩類正確分開，而且保證了分類精度。

非線性映射是SVM的理論基礎(chǔ)，SVM利用內(nèi)積核函數(shù)代替向高維空間的映射，避免了在高維空間的直接計算，解決了數(shù)據(jù)分類的“維數(shù)災(zāi)難”問題，故對于核函數(shù)的選擇是其性能體現(xiàn)的重要部分。

本實驗使用一些數(shù)據(jù)作為訓(xùn)練集，通過網(wǎng)格尋優(yōu)法對核函數(shù)進(jìn)行優(yōu)化，并采用三次交叉驗證方法進(jìn)行分類，確保了分類準(zhǔn)確率。

2 實驗結(jié)果

2.1 實驗數(shù)據(jù)

實驗用的腦電信息數(shù)據(jù)來自于2008年BCI Competition IV Dataset 1，由Berlin BCI研究組提供。該數(shù)據(jù)集記錄了7個健康的受試者，分別記為受試者a，b，c，…，g，實驗過程如圖2所示。

圖2 實驗過程Fig.2 Experiment process

7名受試者都是右利手，視力正常，全身放松，坐在扶手椅上，受試者分別完成200次實驗，每次實驗持續(xù)8 s，實驗記錄過程如下：前2 s為安靜準(zhǔn)備時間，同時在監(jiān)控器屏幕上顯示十字符號，在第2 s結(jié)束后，監(jiān)控器顯示3個方向（向左或向右或向下）的箭頭，受試者通過箭頭方向做相應(yīng)的想象動作，接下來的4 s持續(xù)想象狀態(tài)，最后2 s黑屏，一次完整的實驗結(jié)束。對于測試數(shù)據(jù)，通過聲音提示受試執(zhí)行具體的想象任務(wù)，當(dāng)提示“STOP”時，停止想象，想象時間持續(xù)1.5～8 s，此后是1.5～8 s的安靜準(zhǔn)備時間。腦電信號采集電極分布如圖3所示。

圖3 腦電信號采集電極分布Fig.3 Electrode position profile of EEG collection

實驗中使用Ag/AgCl電極帽和Brain Amp MR放大器，對59個測試電極探測到的EEG記錄并放大。信號帶通濾波為0.05～200 Hz，采樣頻率100 Hz。數(shù)據(jù)集中的受試者b和g執(zhí)行想象左手和右手的運動，選取這兩組數(shù)據(jù)集為研究對象，考慮到實際應(yīng)用的可行性和簡易性，選取數(shù)據(jù)集中的兩個通道（C4和C3）的數(shù)據(jù)（b和g，分別想象左右手運動）為研究對象。

2.2 腦電信號預(yù)處理結(jié)果

量發(fā)生變化，產(chǎn)生μ和β節(jié)律。這一過程中，事件相關(guān)同步（ERS）增強(qiáng)了μ和β節(jié)律頻率，而這兩個節(jié)律的能量降低被標(biāo)記為事件相關(guān)去同步（ERD）。

對于腦電信號，運動想象過程中頻率主要集中在30 Hz以下，≤5 Hz的一般為偽信號，因此，腦電信號的有用頻帶為7～30 Hz。選取腦電信號C3，C4通道的左手想象運動腦電信號進(jìn)行預(yù)處理，得到傅里葉變換頻譜圖，如圖4所示。

圖4 左手運動想象的腦電信號圖和頻譜圖Fig.4 EEG signal and spectrum diagrams of left-hand motor imagination

當(dāng)受試者肢體運動成像時，感覺運動皮層能

2.3 DWT實驗結(jié)果

對腦電信號進(jìn)行預(yù)處理后，本文利用db4小波進(jìn)行六階小波分解，得到cAL（低頻子帶）和一系列高頻子帶（c DL，cDL-1，…，cD1），如表1所示，選取D4和D3子帶信號進(jìn)行信號重構(gòu)。

表1 子帶頻帶的頻率分布Tab.1 Distribution of each sub-band frequency band

如圖5所示，子帶頻帶頻率在8～30 Hz間，頻率有明顯的突起峰值，表明頻率主要集中在一個很窄的頻帶內(nèi)，為后續(xù)的MEMD分解做準(zhǔn)備。

圖5 小波分解子帶信號和頻譜Fig.5 Wavelet decomposition of sub-band signal of waveform and spectrum

2.4 MEMD結(jié)果

計算C3，C4通道IMF分量的最大互信息系數(shù)值，結(jié)果如表2所示。

由表2可以看出，前兩階的最大互信息系數(shù)是依次遞增的，而后四階的系數(shù)在逐漸減小，使得信息特征主要存在于前兩階，因為最大互信息系數(shù)與有用信息的多少成正相關(guān)，選取前2階IMF1～I(xiàn)MF2進(jìn)行下一步處理。

表2 最大互信息系數(shù)值Tab.2 Table of maximum correlation numbers

運動成像最突出的節(jié)律是μ和β節(jié)律，兩者頻率分別為7～13 Hz和13～30 Hz。預(yù)處理后的腦電信號進(jìn)行DWT得到子帶信號，選擇合適的信號進(jìn)行MEMD分解，如圖6～7所示。IMF的前兩階包含μ和β節(jié)律，因此，本文只選取IMF的前兩階重構(gòu)腦電圖信號。MEMD不僅提取了μ和β節(jié)律，而且有效去除了不相關(guān)的信號，提高了信噪比。

圖6 MEMD分解D3子帶的波形圖和頻譜圖Fig.6 Sub-band D3 multivariate empirical mode decomposition of IMF and the corresponding order spectrum

2.5 模糊熵結(jié)果

使用MEMD之前，選取C3（右手）和C4（左手）兩個子帶信號進(jìn)行小波分解，得到一系列窄帶信號，然后利用MEMD對適當(dāng)?shù)淖訋盘栠M(jìn)行分解，得到多個模態(tài)函數(shù)，將μ和β節(jié)律的前兩階作為新的輸入信號，這相當(dāng)于總共8個通道數(shù)據(jù)，形成一個8×2 000的矩陣Xi（i=L或R，表示想象左手或右手的運動）。對原始腦電圖信號進(jìn)行小波分解，選取D3和D4兩個子帶信號，進(jìn)行MEMD分解后選取IMFs的前兩階，共8個通道。2 000是一個測試采樣點的數(shù)量，計算模糊熵時，取r=0.25Sd，使用時間滑動窗口，窗口長度為500個數(shù)據(jù)點，每次移動一個數(shù)據(jù)點，共得到1 500個模糊熵。將同一滑動窗口中8個IMF的模糊熵值作為一組特征向量進(jìn)行分類。

3 分析與討論

圖7 MEMD分解D4子帶的波形圖和頻譜圖Fig.7 Sub-band D4 multivariate empirical mode decomposition of IMF and the corresponding order spectrum

本文提出一種基于DWT、MEMD和模糊熵的腦電信號特征提取與分類方法。每次測試的C3和C4通道數(shù)據(jù)經(jīng)過小波分解后，選擇子帶信號，MEMD選擇IMF。利用IMF分量的前兩階作為觀測向量，進(jìn)一步驗證觀測信號的有效性。圖8～9為單次實驗中受試者做左右想象運動時腦電圖前兩階IMF分量的腦地形圖。由圖8～9可以看出，當(dāng)想象左右手動作時，C3和C4的電極能量差異很大，當(dāng)想象左側(cè)動作時，可以明顯看出活躍度更大的是C4電極。當(dāng)受試者想象右手運動時，結(jié)果恰恰相反。這與事件相關(guān)同步（ERS）和事件相關(guān)去同步（ERD）現(xiàn)象吻合，說明該特征可以有效地表示“源”的空間位置信息，可作為分類特征。通過計算近似熵，可以為運動想象的判斷提供一種可選方法，從而提高分類識別。

圖8 一階IMF能量腦地形圖Fig.8 Brain topographic map of the first-order IMF energy

圖10顯示了10次交叉實驗的分類準(zhǔn)確率和總體平均分類準(zhǔn)確率。實驗對原始腦電圖信號進(jìn)行小波變換后，選擇合適的子帶進(jìn)行MEMD分解，得到多個固有模態(tài)函數(shù)IMF，然后計算8個IMFs采樣點組成數(shù)據(jù)的模糊值，最后利用SVM對模糊熵進(jìn)行分類。

由圖10可見，除了個別正確率較低外，分類正確率大體在90%以上，可見該方法的分類性能較好，并且魯棒性較高。

比較第三屆BCI大賽中的3組數(shù)據(jù)，如表3所示。第一組是將多特征組合，將基于運動相關(guān)電位和事件相關(guān)去同步的3個特征向量通過公共空間子空間分解算法和波形均值提取出來，然后通過Fisher判別分析將其降為一維，并將它們連接成一個三維特征向量，利用線性支持向量機(jī)進(jìn)行分類。第二組使用的特征是所選信道在某段時間內(nèi)的平均值和信道頻率間隔的平均功率。這些特征是在可視化給定類的差異幫助下手工構(gòu)建的，并使用邏輯回歸分類器進(jìn)行分類。第三組提出了一種多變量經(jīng)驗?zāi)B(tài)分解（MEMD）和短時傅里葉變換（STFT）的混合方法，從腦電圖信號中識別左右想象運動，結(jié)果表明k近鄰（KNN）是最佳分類模型。為驗證本文算法的分類精度，使以下4個算法在同一個數(shù)據(jù)集中計算。

比較本文結(jié)果與第三屆BCI競賽中的3組數(shù)據(jù)，結(jié)果如表3所示。由表3可以得出，本文方法得到的分類精度比其他3組的高；僅選用了C3，C4兩個通道進(jìn)行試驗，與其他試驗相比，大大減少了通道數(shù)，同時也保證了分類精度。

表3 腦機(jī)接口（BCI）競賽結(jié)果與本文結(jié)果比較Tab.3 Comparison of BCI contest results with the results in this paper

4 結(jié) 語

為了解決EMD中頻帶覆蓋較廣問題，進(jìn)一步提高腦電信號運動成像的分類精度，提出了一種基于DWT、MEMD和模糊熵的腦電信號特征提取與分類方法。實驗結(jié)果表明，該方法可行且有效，但模糊熵的統(tǒng)計量存在不一致性，因此，在準(zhǔn)確性方面仍有較大的提高空間。