于路，薄華
（上海海事大學 信息工程學院，上海 201306）
摘要：針對現有的單一特征提取算法對運動想象腦電信號識別率不高的問題，提出一種以相關系數改進的經驗模態分解（EMD）的特征提取算法。對已有的BCI競賽數據中C3、C4兩個通道腦電數據進行預處理，之后通過EMD對腦電信號進行分解，得到IMF分量。通過計算原始信號與各階IMF分量之間的相關系數，選擇具有較大相關系數的IMF作為特征，由這些IMF分量的能量特征和平均幅值差來組成腦電信號的特征。使用支撐矢量機分類器（SVM）對左右手運動想象腦電信號進行分類。實驗結果表明，基于相關系數改進的EMD腦電信號的處理方法明顯優于只用EMD的腦電處理方法，得到的最高正確識別率為88.57%。從而證明了該方法的有效性。
關鍵詞： 腦電信號；經驗模態分解；相關系數；特征提取
0引言
頭皮腦電（Electroenc EphaloGraph,EEG）信號是一種產生機理相當復雜且非常微弱的隨機信號，它反映了大腦阻止的腦電活動及大腦的功能狀態。不同的思維狀態在大腦皮層有不同的反映［1］。
腦機接口（BrainComputer Interface，BCI）不依賴于大腦外周神經與肌肉系統，是大腦與計算機或其他電子設備之間建立可直接交流和控制的通道，可以有效增強身體嚴重殘疾的患者與外界交流或控制外部環境的能力，提高患者的生活質量，同時也在娛樂電競等領域有著巨大的應用前景［2-5］。
目前EEG特征提取的常用方法有FFT、SFT、AR［5］、AAR［6］、ICA、小波分析［7］、特定頻帶的功率譜等方法。經典的FFT在分析確定信號和平穩信號時發揮了重要作用,但利用FFT分析突變信號的頻譜存在局限性。SFT在一定程度上克服了標準FFT不具有局部分析能力的缺陷,在某些信號處理中發揮了一定的作用，但也存在不可克服的缺陷。SFT是時間窗內信號特征的平均,時間窗內信號越短,獲得的時間分辨率就越高。根據信號測不準原理,時間局域化性質和頻率局域化性質是矛盾的。AAR模型參數隨每一樣本點的輸入而改變,因而更好地反映了大腦的狀態。但是，該方法更適合分析平穩信號，對包含高度非平穩信號的運動想象EEG，該模型達不到理想效果。小波變換和小波包變換分解信號時，要預先設好分層數和小波函數，不具備對信號自適應的分解能力。因此需要采用一種非線性分析的方法，該方法同時具有對信號自適應的分解能力，能通過度量腦電信號的復雜度來反映腦電的特征。本文中采用的是用相關系數改進經驗模態分解的腦電信號處理方法，可以獲得一系列固有模態函數分量（Intrinsic Mode Function，IMF）。將這種方法應用到腦電信號的處理過程中能得到較好的效果。
1EMD算法
經驗模式分解算法(Empirical Mode Decomposition,EMD)［8］是一種自適應的數據處理或挖掘方法，非常適合非線性非平穩的時間序列。其本質是通過數據的時間尺度來獲得本征波動模式，然后分解數據。這種分解過程可以形象地稱之為“篩選（sifting）”過程。
在理論上，EMD可以應用于任何類型的時間序列（信號）的分解，因而在處理非平穩及非線性數據上，比之前的方法具有更明顯的優勢。所以，EMD［9-10］方法一經提出就在不同的工程領域得到了迅速有效的應用。
該方法的關鍵是它能使復雜信號分解為有限個本征模函數，所分解出來的各IMF分量包含了原信號的不同時間尺度的局部特征信號。EMD分解方法基于以下假設條件：
(1)數據至少有兩個極值，一個最大值和一個最小值；
(2)數據的局部時域特性是由極值點間的時間尺度來唯一確定；
(3)如果數據沒有極值點但有拐點，則可以通過對數據微分一次或多次求得極值，然后再通過積分來獲得分解結果。
設原始信號為X(t),EMD算法的計算步驟如下：
（1）找出原數據序列X(t)的所有極大值點和極小值點，將其用三次樣條函數擬合為原序列的上、下包絡線,分別為u(t)和v(t)，可得包絡線的平均值m11：
m11=1/2(u(t)+v(t))(1)
（2）將原數據序列減去包絡平均值m11，得到一個減去低頻的新序列h11:
h11=X(t)－m11(t)(2)
（3）h11不一定是平穩數據序列，用h11(t)代替原始信號X(t),重復上述過程k次,直到所得包絡趨近于零，這樣可以得到第一個本征模函數（IMF）。分量c1=h1k(t)，它表示信號數據序列最高頻率的成分。
（4）用X(t)減去c1，得到一個去掉高頻成分的新數據序列r1；對r1再進行上述分解，得到第二個本征模函數（IMF）分量c2；如此重復:
r2=r1－c2
r3=r2－c3
…
rn=rn－1－cn
（5）當第n個剩余量rn已成為單調函數，無法再分解IMF時，整個EMD分解過程完成。原始信號可以表示為將一個頻率不規則的波化為多個單一頻率的波與殘波相加的形式。即：
X(t)=∑ci+rn(3)
然而腦電信號其背景噪聲很強，在對腦電信號進行EMD分解時，通常不能將信號與噪聲徹底分開，在每一個輸出的軌道之中，或多或少都摻雜著一些噪聲，選擇合適的IMF分量顯得尤為重要。這里引入相關系數的概念，借此權衡各個IMF的有效性。
2運動想象腦電信號識別算法
2.1基于相關系數的特征提取
相關系數是用以反映變量之間相關關系密切程度的統計指標。相關系數是按積差方法進行計算，同樣以兩變量與各自平均值的離差為基礎，通過兩個離差相乘來反映兩變量之間的相關程度。是判斷特征選取的一個方法。
相關系數的定義：假設有兩個隨機變量X和Y，則它們的相關系數為：
r=Cov(X,Y)σXσY
=E［{X－E(X)}{Y－E(Y)}］［E{X－E(X)2}］12［E{Y－E(Y)}2］12(4）
其中，Cov(X,Y)為隨機變量X與Y的協方差函數，σX、σY分別指X、Y的標準差，E(X)、E(Y)為兩者的平均值。相關系數r的取值范圍是［-1,1］，表示變量之間相關程度的高低，r的絕對值越大，說明這兩個變量之間相關程度越高。r>0表示正相關，r<0表示負相關，特殊地，r=1稱為完全正相關，r=-1稱為完全負相關，r=0稱為不相關。
對原信號進行EMD分解后，得到n階IMF分量和殘波。根據各個IMF分量與原信號的相關系數來選擇更合適的IMF分量。設原信號序列為X(n),信號的IMF分量為Yi(n),其中i表示第i個IMF分量，則相關系數的定義如下：
使用EMD分解原信號得到IMF分量后，在眾多分量中包含的有用信息各不相同，同時包含了噪聲。為了得到更好的特征，提高左右腦電信號的分辨率，本文通過相關系數選擇合適的分量。分別取各階IMF分量的能量單獨作為該信號的特征，因為每組腦電信號只取一個IMF分量，所以這里只考慮能量特征。計算公式如下：
El=∑ni=1［c(i)2］(6)
其中，El對應了第l個IMF分量的能量，c(i)表示該分量中的第i個數值，n表示該分量的長度。本文中n為140。
一組信號可得到唯一能量特征。表1給出各階IMF分量與原信號之間的相關系數，及只取相應能量特征作為腦電信號的特征后的左右手腦電信號識別率。
由表1可知，IMF2分量與原信號的相關系數最大，最適合進行特征提取。各階IMF分量以與原信號的相關系數從大到小排列，分別為IMF2、IMF1、IMF3、IMF4和IMF5。當IMF分量與原信號的相關系數越大時，采用該分量作為條件提取特征得到的信號識別率越高，從而該IMF分量比其他IMF分量更適合進行特征提取。
然而只用單一的IMF分量并不能達到很高的識別率。這里考慮使用多個IMF分量進行特征選取。當選取多個IMF分量時，除能量特征外，腦電信號的特征還包括平均幅度差。
在進行EMD算法分解時，得到的各階IMF分量是從高頻到低頻的腦電信號分解，所以可以看出，同階不同通道信號之間有時幅值波動相差過大，因此此處定義平均幅度差作為一個特征值，計算公式如下：
其中，ci是第i個IMF分量，cj表示第j個分量(i≠j),n為信號的長度，這里n為140。分別計算不同IMF分量個數來提取特征時的腦電信號分類正確率，如表2所示。
表2中，分量順序根據與原信號的相關系數從大到小排列，分別是：IMF2、IMF1、IMF3、IMF4、IMF5。分量選取從第一個開始取，如當IMF分量個數為2時，所取分量為IMF2、IMF1。實驗數據表明，取前3個分量達到了最好的分類效果。當所取IMF分量個數過多時，信號包含的無效信息比例也會同時增加，對分類會造成更大的影響，且前3個分量包含了原始信號約90%的能量，故能達到較好效果。
2.2本文算法
本文選擇的腦電信號識別算法是基于相關系數改進的EMD算法。EMD算法是一種非線性分析方法，同時具有對信號自適應的分解能力，能通過度量腦電信號的復雜度來反映腦電的特征。對得到的腦電數據進行預處理后，進行EMD分解從而得到IMF分量和余量。但是腦電信號中仍然包含了無用信息，在使用EMD算法進行分解過程中，無法將噪聲信號徹底分解出來。故采用基于相關系數改進的EMD算法，選擇與原信號相關系數最大的3個IMF分量：IMF1、IMF2、IMF3。由這3個分量可以得到3個能量特征和兩個平均幅度差特征，送入SVM［1112］進行分類。在SVM過程中先用訓練數據得出訓練模型，再用測試數據得到實驗結果。本文算法的具體步驟如下：
（1）對腦電信號進行預處理；
（2）從訓練數據中取一組信號中1 s的腦電數據進行EMD算法分解；
（3）計算各階IMF分量與原信號的相關系數，選出3個相關系數最大的IMF分量；
（4）由IMF分量得到能量特征和平均幅度差；
（5）將數據特征送入SVM進行訓練；
（6）載入測試數據得到結果，從而得出實驗結論。
3實驗分析
3.1實驗數據
本文實驗數據采用奧地利格拉茲工業大學腦機接口研究中心提供的運動想象腦電數據(Data setⅢ)［13］。該數據集的受試主體為25歲女性，主體以放松狀態坐在顯示器前，顯示器呈現提示信息。根據出現的左右線索想象左右手運動，從而得到反饋數據，其中左右線索的順序是隨機的。
實驗共得到280次長度為9 s的腦電數據，140次為訓練樣本（70 次想象左手運動,70 次想象右手運動)，另140次作為測試樣本。在前2 s主體保持放松，在第2 s屏幕出現十字光標，持續時間1 s，提示實驗即將開始。第3 s 屏幕出現左右箭頭，主體根據提示信息想象左右手運動，即該數據的有效時間段為4 s~9 s。
該信號的采集頻率為128 Hz,再經過0.5~30 Hz的帶通濾波器。實驗采用Ag/AgCl電極，通過C3、Cz、C4三個通道獲得反饋數據，其中，C3、C4電極位于大腦的初級感覺皮層運動功能區，能反映主體在想象左右手運動時大腦狀態的變化，Cz作為參考電極。實際分析時只采用了 C3、C4 這 2 個通道的數據。
3.2實驗分析
實驗數據為280組長度為9 s的腦電數據，其中4~9 s為有效時間區間，但是考慮信號的有效性，選擇4~8 s的腦電數據作為實驗數據。以一個通道的1 s數據作為一小段數據進行EMD分解，得到IMF分量后，計算各階IMF分量與原始信號的相關系數，按照相關系數由大到小選取最大的3階IMF分量進行特征提取。對于一組一個通道中1 s的數據包含3個能量特征和2個平均幅值差，則一組數據一共包含40個特征。采用十折交叉驗證法將訓練數據的特征采樣送入SVM，并確定SVM的核參數。將測試數據提取特征后送入SVM的分類模型進行分類。分類結果如表3所示。
從表3可以看出，采用結合相關系數的EMD特征提取方式能獲得更高的識別率。其中訓練集的識別率達到92.86%，測試集的分類正確率達到88.57%，這個結果非常接近BCI大賽中第一名的分類結果。由此改進后的EMD算法更加適合腦電特征提取。
4結束語
本文提出通過相關系數來改進EMD對腦電信號特征的提取［14］。通過分析EMD分量與原信號之間的相關系數來確定相應的特征組成特征向量，輸入SVM分類器中，從而實現左右手運動想象腦電信號的分類。研究結果表明，采用相關系數改進的EMD提取腦電特征的正確率明顯高于僅采用EMD的腦電特征提取方法。因此，基于相關系數改進的EMD算法在運動想象腦電信號的識別研究中具有很高的應用價值。
參考文獻
［1］ 徐寶國，何樂生，宋愛國.基于腦電信號的人機交互實驗平臺的設計和應用［J］.電子測量與儀器學報，2008, 22(1):8185.
［2］ 高上凱. 淺談腦機接口的發展現狀與挑戰［J］. 中國生物醫學工程學報, 2007. 26(6):801803.
［3］ 王斐，張育中,寧廷會，等.腦—機接口研究進展［J］．智能系統學報，2011，6(3):189199．
［4］ TANAKA K,MATSUNAGA K,WANG H.Electroencephalogrambased control of an electric wheelchair［J］.IEEE Transactions on Robotics，2005，21( 4) : 762766．
［5］ 張毅，楊柳，李敏，等.基于AR 和SVM的運動想象腦電信號識別［J］．華中科技大學學報（自然科學版），2011，39(Z2):103106．
［6］ 徐寶國，宋愛國.單次運動想象腦電的特征提取和分類［J］．東南大學學報（自然科學版），2007，37(4):629633．
［7］ 李明愛，王蕊，郝冬梅．想象左右手運動的腦電特征提取及分類研究［J］．中國生物醫學工程學報，2009，28(2):166170．
［8］ 張小莉，張歆，孫進才.基于經驗模態分解的目標特征提取與選擇［J］.西北工業大學學報，2006,24(4):453456.
［9］ 余煒，韓強，馬晶晶，等.基于EMD和SVM的腦電信號處理方法［J］.昆明理工大學學報（自然科學版），2012,37(6):3842.
［10］ 金晶，王行愚，張秀，等.基于能量特征的左右手運動想象鬧心好的識別方法［J］.科學通報，2001,46(3):257263.
［11］ 安文娟. Fisher和支持向量綜合分類器［D］. 大連：遼寧師范大學, 2010.
［12］ 李金華.基于SVM的多類文本分類研究［D］. 青島：山東科技大學,2010.
［13］ Sun Shiliang, Zhang Changshui.Adaptive feature extraction for EEG signal classification［J］. Medical & Biological Engineering & Computing,2006,44(10)：931935.
［14］ RAKOTOMAMONJY A, GUIGUE V. BCI competition III: dataset IIensemble of SVERP for BCI P300 speller［J］.IEEE Transactions on Biomedical Engineering, 2008, 55 (3):11471154 .