基于特征篩選的經(jīng)驗模式分解腦電信號分類方法

張學軍， 汪 敏， 胡曉雯

(1.南京郵電大學電子與光學工程學院， 南京 210023； 2.南京郵電大學射頻集成與微組裝技術(shù)國家地方聯(lián)合工程實驗室， 南京 210023；3．南京醫(yī)科大學生物醫(yī)學工程與信息學院， 南京 211166)

腦電(electroencephalogram，EEG)信號是一組包含人類大腦活動信息的電位差，展示了有關(guān)腦電流的數(shù)據(jù)。EEG的測量可以通過放置在頭皮上的傳感器或使用顱內(nèi)電極來獲得。EEG應用領域廣泛，如情緒識別、腦-機接口(brain-computer interface，BCI)等[1]。

基于運動想象的BCI以其無創(chuàng)、適用性好、可移植性強等優(yōu)點吸引了眾多研究者的興趣。BCI為人們提供了一種使用腦電信號與外部輔助設備交互的方式，在生物醫(yī)學工程和神經(jīng)修復中有著廣泛的應用[2]。首先，受試者需要在大腦中進行基于運動想象(motion imagination，MI)的特定運動，其次，基于MI的BCI會對特定運動的腦電信號進行采集、分類，將其轉(zhuǎn)換為手、腳等不同想象任務的控制信號。

BCI研究的目標是開發(fā)幫助殘疾用戶與他人交流系統(tǒng)。BCI系統(tǒng)是一個連續(xù)閉環(huán)系統(tǒng)，通常由腦信號獲取、預處理、特征提取、分類、輸出命令和反饋五個部分組成。預處理方法如濾波或盲源分離算法可以減少不同類型的腦電噪聲和眼電尾跡。

由于腦電信號具有噪聲強、強度低、對周圍環(huán)境敏感、數(shù)據(jù)維數(shù)大和分布復雜等特點，在實際應用中還不能很好地分類。為了減輕噪聲腦電信號干擾的影響，可以利用各種信號處理方法從腦電信號中提取的參數(shù)對腦電信號的分類。基于傅里葉變換的譜參數(shù)對腦電信號的分析在腦電信號的分類上取得了良好的效果。然而，傅里葉變換沒有時頻的局部化性能，從而傅里葉域在信號中不表現(xiàn)出任何時域特征[3]。Gabor提出使用短時傅里葉變換(short time fourier transform，STFT)解決表示時域特征的問題，但STFT并不能對信號進行多分辨率分析，這是因為STFT使用相同帶寬的濾波器對所有頻率的信號進行分解[4]。通常利用小波分析來解決。在小波分析中，通過形成具有不同帶寬的帶通濾波器來促進多分辨率時頻分析。研究人員發(fā)現(xiàn)小波分析對于各種信號處理應用是非常有用的工具。在腦電信號處理方面，Zhang等[5]提出了用于腦電信號處理的小波CSP(common space pattern)算法。姚悅等[6]提出小波變換結(jié)合二階盲辨識的眼電偽跡自動去除方法，其分離效果不受 源信號高斯性影響，在抑制腦電信號眼電尾跡方面取得了良好的效果。盡管上述文獻中使用小波進行信號處理，但都是基于被分析信號是平穩(wěn)的。對腦電信號分析的研究表明，腦電信號的頻率成分在一段時間內(nèi)會發(fā)生變化。經(jīng)驗模態(tài)分解(empirical mode decomposition，EMD)是一種基于時頻的經(jīng)驗方法，對于非平穩(wěn)信號的時頻分析是有效的[7]，它將信號分解為若干固有模函數(shù)(intrinsic mode function，IMF)，IMF是振蕩分量。

為了控制BCI系統(tǒng)來識別用戶的活動，并將其轉(zhuǎn)換為命令，在大多數(shù)現(xiàn)有的BCI中依賴于分類算法，不同的分類算法用于不同的BCI應用(例如線性分類器、神經(jīng)網(wǎng)絡、非線性貝葉斯分類器、最近鄰分類器)。線性分類器是使用線性函數(shù)來分類的算法。BCI系統(tǒng)的設計主要采用兩種線性分類器，即線性判別分析(linear discriminant analysis，LDA)和支持向量機(support vector machine，SVM)。但是使用LDA的前提是假設數(shù)據(jù)是正態(tài)分布的，對于EEG數(shù)據(jù)分類敏感度不高。選用SVM分類器，SVM的主要思想是構(gòu)造一個最優(yōu)的超平面，使分割平面和數(shù)據(jù)之間的差值最大化。使用核函數(shù)可以向高維空間進行映射，解決非線性的分類[8]。 SVM最佳參數(shù)可以通過網(wǎng)格搜索算法進行調(diào)整，但是精度較低。

現(xiàn)提出一種基于小波變換和EMD結(jié)合排列熵的特征篩選的分類方法。利用小波變換提取腦電運動想象窄帶信號，再進一步EMD提取窄帶信號的IMF分量。為了減少冗余特征對算法精度的影響，計算每個IMF的排列熵，進行Perason系數(shù)比較篩選出合適的特征，最后運用遺傳算法(genetic algorithm, GA)對SVM分類方法進行優(yōu)化。

1 特征提取與分類方法

1.1 小波分解

小波變換(wavelet transform，WT)是描述非平穩(wěn)信號的有力數(shù)學工具，是一種變換分析方法，能夠在時間、空間頻率的局部化分析，通過伸縮平移運算對信號逐步進行多尺度細化，最終達到高頻處時間細分，低頻處頻率細分，能自動適應時頻信號分析的要求[9]。它具有自適應性和多分辨率能力，因此適合于分解時頻分辨率不同的腦電信號。

小波變換是將母小波函數(shù)作位移τ后，再在不同的尺度a下，與待分析信號X(t)作內(nèi)積，即

(1)

式(1)中:a> 0，稱為尺度因子;τ反映位移,可正可負；ψ()為母小波函數(shù)；WTX為小波變換后的函數(shù)；t為時間。離散小波變換DWT對尺度參數(shù)按冪級數(shù)進行離散化處理，對時間進行均勻離散化取值如二進制離散化尺度時間為2，4，6，…,2n(要求采樣率滿足尼奎斯特采樣定理)。

利用Mallat算法進行小波分解過程如圖1所示。

CA為低頻信息、近似分量；CD為高頻、細節(jié)分量

1.2 經(jīng)驗模態(tài)分解

EMD具有直觀性和自適應性，對EEG等非平穩(wěn)信號表現(xiàn)良好。EMD的目的是將信號分解為一組固有模函數(shù)IMF。IMF定義為極值數(shù)和過零點數(shù)相等(或最多相差一個)的函數(shù)，其包絡由所有局部極大值和極小值定義，相對于零對稱[10]。IMF表示一個簡單的振蕩模式，作為傅里葉分析中使用的簡單諧波函數(shù)的對應項[11]。

給定一個信號x(n)，EMD的起點是識別所有的局部極大值和極小值。所有局部極大值以三次樣條曲線作為上包絡eu(n)連接，同樣，所有局部極小值以樣條曲線作為下包絡el(n)連接。兩個包絡的平均值表示為ml(n)=[eu(n)+el(n)]/2，則獲得第一個原始IMF分量hl(n)為

hl(n)=x(n)-ml(n)

(2)

上述提取IMF的過程稱為篩選過程。由于hl(n)仍然在零交叉點之間包含多個極值，因此對hl(n)再次執(zhí)行篩選過程。將此過程重復應用于原始IMFhk(n)，直到得到滿足固有模態(tài)函數(shù)條件[12]。

rl(n)=x(n)-hl(n)

(3)

式(3)中：rl(n)為殘余分量，包含一些有用的信息。因此，可以將殘余當作一個新的信號，并應用上述過程來獲得新的IMF。

ri-1(n)-ci(n)=ri(n),i=2，3,…，N

(4)

當殘余rN(n)是常數(shù)、單調(diào)斜率或只有一個極值的函數(shù)時，整個過程終止。結(jié)合式(3)和式(4)中的方程得到原始信號：

(5)

1.3 排列熵

考慮標量時間序列X(t)(t∈1,2,…,N)，非線性數(shù)據(jù)分析的第一步是相空間重構(gòu)，最常用的方法是利用延遲時間嵌入定理。排列熵的原理是引入一個時間延遲和嵌入維度[13]。

在這種方法中，時間序列的值被轉(zhuǎn)換成一個延遲向量：

Xi→[Xi-(d-1)τ，Xi-(d-2)τ，…，

Xi-τ，Xi]

(6)

式(6)中:d為嵌入維度;τ為時間延遲。這會將N個標量轉(zhuǎn)換為具有重疊項的N-τ(d-1)向量。

可以按遞增順序排列d維延遲向量中的值，以實現(xiàn)有序模式:

[xi-rd-1τ≤xi-rd-2τ≤…≤xi-r1τ≤xi-r0τ]

(7)

它們的相等發(fā)生在設置rl

根據(jù)香農(nóng)熵原理可得排列熵定義為

(8)

顯然，當P(j)的遞增或遞減序列達到下限時,0

(9)

排列熵是時間序列的一種復雜度量，可以快速、簡便地計算有序模式。當處理值之間的順序關(guān)系時，排列熵對于噪聲是魯棒的。

1.4 基于遺傳算法的支持向量機參數(shù)優(yōu)化

基于徑向基函數(shù)(radial basis function，RBF)核函數(shù)的支持向量機分類器在分析高維數(shù)據(jù)方面具有優(yōu)勢，SVM核函數(shù)的參數(shù)對分類器的性能影響很大，需定義兩個參數(shù)懲罰因子c和核參數(shù)g，目前還缺乏有效確定參數(shù)選取的結(jié)構(gòu)方法。因此，將遺傳算法應用于所提出的支持向量機模型中，優(yōu)化參數(shù)選擇[14]。遺傳算法是通過世代搜索而不是單點搜索來獲得最優(yōu)解或準最優(yōu)解；它具有全局尋優(yōu)能力；它是種群變化的并行過程，具有內(nèi)在的并行性。遺傳算法的處理對象是參數(shù)集被編碼的個體，而不是參數(shù)本身，這一特點使得遺傳算法得到了廣泛的應用。

基于遺傳算法的改進支持向量機利用訓練樣本集的輸入，遺傳算法搜索核函數(shù)及其訓練參數(shù)。基于遺傳算法的整個優(yōu)化過程如圖2所示。

圖2 遺傳算法優(yōu)化過程

遺傳算法基于種群中適者生存原則，通過世代傳遞遺傳信息來保留遺傳信息。圖3給出了所提出的支持向量機模型的框架。利用遺傳算法在支持向量機中尋找兩個參數(shù)的最佳組合，得到較小的分類誤差。

圖3 GA-SVM的流程圖

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

研究使用的4個數(shù)據(jù)集來自BCI Competition 2008 data sets 2b數(shù)據(jù)，由奧地利格拉茨理工大學提供。數(shù)據(jù)集由9名受試者的左右手運動圖像數(shù)據(jù)組成。每個運動圖像數(shù)據(jù)由三個電極(C3、Cz和C4)提取，采樣頻率為250 Hz的雙極記錄法進行記錄[15]。濾波器采用0.5～100 Hz的帶通濾波器和50 Hz的陷波濾波器。在前兩個實驗中，每個實驗都從固定十字架和聽覺刺激開始。3 s后，一個可視箭頭指示左或右運動想象任務，此命令持續(xù)1.25 s，然后進行MI任務。MI持續(xù)4 s后，受試者有至少1.5 s的短暫休息時間。圖4所示為數(shù)據(jù)集實驗時間線。

圖4 實驗進程

采用Ag/AgCl電極記錄腦電圖如圖5所示。

圖5 采樣電極位置

3 實驗結(jié)果

3.1 小波分解結(jié)果

離散小波變換(DWT)通過將信號分解為一個粗糙的近似值，并從時域信號的連續(xù)高通濾波和低通濾波中獲得的詳細信息。對于許多信號，低頻成分相當重要，它常常蘊含著信號的特征，而高頻成分則給出信號的細節(jié)或差別[16]。根據(jù)事件相關(guān)同步/去同步原理，進行想象運動會在相關(guān)的頻率段產(chǎn)生能量的變動，因此，α節(jié)律、β節(jié)律的腦電就會被定位。對于EEG信號在30 Hz以上可以采用低通濾波器濾除。在本研究中，使用小波變換來分解C3、C4和Cz通道滑動窗內(nèi)的每個腦電片段。矩形窗口的大小是相鄰窗口重疊的1/2。利用Daubechies 4-tap小波將每段腦電信號分解為4個層次，分解后得到16個不同的頻率段。

頻段分解后，子帶小波系數(shù)所對應的頻率范圍如表 1 所示。

表1 每個子帶頻帶分布

獲得低頻子帶CAL，高頻子帶CDL，CDL-1，…,CD1。結(jié)果存儲在低通信道{A4}和高通信道{D1、D2、D3和D4}中。

子帶D4、D3在腦電信號頻段，則選擇這兩個子帶信號進行信號重構(gòu)。如圖6所示小波分解后的D3子帶信號的頻段集中分部在8～16 Hz，如圖7所示D4則分布在16～30 Hz。

圖6 D3小波分解子帶信號

圖7 D4小波分解子帶信號

3.2 EMD分解結(jié)果

EMD應用篩選過程將信號分解為IMFs和一個殘基的余數(shù)，可以看作是一個單調(diào)的斜率，一個只有一個極值的函數(shù)，或者最后一個IMF。考慮到信號的主觀性質(zhì)、數(shù)據(jù)長度和極值數(shù)目等因素，IMF的數(shù)目將不是恒定的。IMF表示簡單的振蕩模式；而高階IMF表示慢振蕩模式，低階IMF表示快振蕩模式。所有腦電信號可以覆蓋的IMF數(shù)量為8，如圖8所示為D3子帶分解的前8個IMF分量的波形圖[圖8(a)]和頻譜圖[圖8(b)]，圖9所示為D4子帶分解的前8個IMF分量的波形圖[圖9(a)]和頻譜圖[圖9(b)]。

f為頻率

Pearson系數(shù)用來衡量兩個獨立的服從正態(tài)分布的連續(xù)變量之間的相關(guān)性。表2為每個IMF的腦電信號的排列熵特征的Pearson系數(shù)。兩個特征之間的Pearson系數(shù)絕對值越大，表示特征間冗余越大，兩個特征越不應該被同時選中；反之，說明兩個特征間的相關(guān)性越小[16]。表2中綠色底紋所標注的Pearson系數(shù)絕對值都小于0.04，說明具有特征極弱相關(guān)或者無相關(guān)。因此，選擇IMF2、IMF3、IMF7、IMF5的特征進行分類。

表2 每個排列熵特征的Pearson系數(shù)

3.3 GA-SVM分類

使用GA算法對支持向量機參數(shù)的優(yōu)化，可以使SVM更好對不同的特征進行分類[18]。

針對支持向量機存在兩個模型參數(shù)c和g，運用遺傳算法選擇最優(yōu)的c和g對支持向量機進行優(yōu)化。調(diào)用MATLAB遺傳算法工具箱，實現(xiàn)逐步啟發(fā)式優(yōu)化，經(jīng)反復試驗確定最大的進化代數(shù)maxgen=100，種群最大數(shù)為sizepop=40(種群規(guī)模)，交叉概率為pc=0.9，變異概率為pm=0.01，并且設置支持向量機進行10折交叉驗證，得到最優(yōu)c=0.000 913 62，g=0.634 15。

如圖10為使用GA算法優(yōu)化前后分類準確率對比變化圖，可得未用GA算法的最高準確率可到99.22%，最低78.56%，平均分類準確率97.30%。使用GA算法之后分類準確類最高100.00%，最低96.70%，平均達97.64%。推出GA算法優(yōu)化了g和c使準確率更加平穩(wěn)，減少了因個體差異性誤判的可能性[19]。

圖10 GA對排列熵參數(shù)優(yōu)化前后準確率對比

3.4 實驗總結(jié)

為了突出本文方法的優(yōu)越性，將其與其他最新的研究進行了比較。表3給出了本文方法和其他方法的分類結(jié)果的比較。本文方法在使用較少通道數(shù)得到比其他研究方法更優(yōu)的分類結(jié)果。Pearson系數(shù)篩選特征分類結(jié)果進行對比，特征篩選后平均準確率到達97.30%。GA優(yōu)化SVM和特征篩選都起作用，有效地提高了MI-EEG的分類精度[9]。

表3 與2008 BCI競賽結(jié)果比較

4 結(jié)論

本文的總體目標是提出一種應用于腦電信號特征分類的方法。為了去除不太重要的數(shù)據(jù)信息，降低計算復雜度，特征提取方法只使用優(yōu)勢IMF。利用小波變換提取腦電運動想象窄帶信號，再經(jīng)EMD提取窄帶信號的IMF分量，計算IMF的排列熵作為分類特征。對IMF各個排列熵特征的Pearson系數(shù)進行比較，選取較小的IMF的排列熵作為特征值。然后用支持向量機對其進行分類，運用遺傳算法對支持向量機參數(shù)進行優(yōu)化，實現(xiàn)了對腦電信號較高精度的分類。在進一步的研究中，將研究不同的特征選擇和優(yōu)化方法，并與GA-SVM方法進行比較，并將GA-SVM方法應用于BCI系統(tǒng)中。