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

張學軍， 汪 敏， 胡曉雯

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

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

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

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

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

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

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

1 特征提取與分類方法

1.1 小波分解

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

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

(1)

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

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

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

1.2 經驗模態分解

EMD具有直觀性和自適應性，對EEG等非平穩信號表現良好。EMD的目的是將信號分解為一組固有模函數IMF。IMF定義為極值數和過零點數相等(或最多相差一個)的函數，其包絡由所有局部極大值和極小值定義，相對于零對稱[10]。IMF表示一個簡單的振蕩模式，作為傅里葉分析中使用的簡單諧波函數的對應項[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)再次執行篩選過程。將此過程重復應用于原始IMFhk(n)，直到得到滿足固有模態函數條件[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)是常數、單調斜率或只有一個極值的函數時，整個過程終止。結合式(3)和式(4)中的方程得到原始信號：

(5)

1.3 排列熵

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

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

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

Xi-τ，Xi]

(6)

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

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

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

(7)

它們的相等發生在設置rl

根據香農熵原理可得排列熵定義為

(8)

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

(9)

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

1.4 基于遺傳算法的支持向量機參數優化

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

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

圖2 遺傳算法優化過程

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

圖3 GA-SVM的流程圖

2 實驗數據

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

圖4 實驗進程

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

圖5 采樣電極位置

3 實驗結果

3.1 小波分解結果

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

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

表1 每個子帶頻帶分布

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

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

圖6 D3小波分解子帶信號

圖7 D4小波分解子帶信號

3.2 EMD分解結果

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

f為頻率

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

表2 每個排列熵特征的Pearson系數

3.3 GA-SVM分類

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

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

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

圖10 GA對排列熵參數優化前后準確率對比

3.4 實驗總結

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

表3 與2008 BCI競賽結果比較

4 結論

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