針對腦機接口(BCI)系統中的多通道非平穩腦電(EEG)信號和腦磁(MEG)信號, 本文提出一種基于多通道經驗模式分解(MEMD)與功率特征結合的信號特征提取算法。首先將多通道腦信號經MEMD算法分解為一系列多尺度多元固有模態函數(IMF)近似平穩分量, 然后對每個IMF分量提取功率特征, 并利用主成分分析(PCA)降維處理, 最后使用線性判別分析分類器對信號特征分類。實驗采用第三次和第四次國際BCI競賽的數據進行驗證, 對皮層EEG信號和MEG信號運動想象任務的識別正確率分別達到92.0%和46.2%, 均位于競賽第一名水平。實驗結果表明本文所提方法有較好有效性和穩定性, 為腦信號特征提取提供了新思路。
引用本文： 王金甲, 劉源. 基于多通道經驗模式分解的腦機接口特征提取. 生物醫學工程學雜志, 2015, 32(2): 451-454, 464. doi: 10.7507/1001-5515.20150081 復制
版權信息： ?四川大學華西醫院華西期刊社《生物醫學工程學雜志》版權所有，未經授權不得轉載、改編
引言
腦-機接口(brain-computer interface, BCI)系統能實現利用人腦意識控制外設[1]。BCI需要提取反映腦信號意識任務的特征，常見的特征提取方法有時域能量分析、AR模型系數、小波變換系數、共空間模式等。
經驗模態分解(empirical mode decomposition, EMD)算法[2]已成功用于腦電(electroencephalogram, EEG)信號處理[3]。文獻[4]利用EMD與AR模型結合的方法提取腦磁(magnetoencephalogram, MEG)信號特征。標準EMD算法分析多通道信號時，往往需要分別分析每一個通道信號或者先將多通道信號整合為單通道信號后再分析，忽略了腦信號的空間特性，不利于各通道間同一頻率尺度對比分析。因此標準EMD不適合直接處理多通道信號[5]。于是EMD的多元擴展方法不斷被提出，例如復數EMD[6]、二元EMD[7]、三元EMD[8]、多元EMD(multivariate EMD, MEMD)[9]。MEMD利用多維空間超球面方向向量求均值包絡，實現不同頻率尺度內各通道信號的共同模式的分析。
本文提出了一種基于MEMD與功率特征相結合的BCI特征提取算法。首先對多通道腦信號進行MEMD，然后提取IMF分量的功率特征并進行主成分分析(principal component analysis, PCA)降維，最后利用線性判別分析進行分類。腦機接口競賽BCI CompetitionⅢ和BCI CompetitionⅣ中的兩組數據的實驗結果令人滿意。
1 MEMD算法
MEMD實現了多元信號的多通道同步聯合分析。多元信號每個通道分解出的IMF個數相同，各層IMF頻率不同，先分解出的IMF層頻率高，后分解出的IMF頻率低。每個通道對應的IMF按頻率尺度對齊，形成多元IMF。
MEMD算法關鍵是求解多元信號的局域均值，這需要在多維空間建立方向向量集，而合適向量集的問題就轉化為在超球面找合適采樣點集的問題。
設一個n維向量組序列{v(t)}t-1T={v1(t), v2(t), …, vn(t)}，代表一個n元信號，長度為T，xθk={x1k, x2k, …, xnk, }表示在(n-1)維球面上對應角θk={θ1k, θ2k, …, θn-1k}的方向向量集, k=1, 2, …, k。MEMD算法的步驟總結如下：
(1)采用Hammersley序列采樣法，在(n-1)維球面上獲得合適的均勻采樣點集，即得到n維空間的方向向量。
(2)計算輸入信號v(t)在每個方向向量xθk上的映射pθk(t)。
(3)確定所有方向向量的映射信號{pθk(t)}k=1K極值對應的瞬時時刻{tlθk}k=1K，l表示極值點位置，l∈[1, T]。
(4)用多元樣條插值函數插值極值點[tlθk, v(tlθk)]，得到K個多元包絡{eθk(t)}k=1K。
(5)對球面空間K個方向向量，n元信號均值m(t)為：
| $m\left(t \right)=\frac{1}{K}\sum\limits_{k=1}^K {{e^{{\theta _k}}}\left(t \right)} $ |
(6)通過h(t)=v(t)-m(t)提取固有模態函數h(t)，如果h(t)滿足多元IMF判斷標準[10]，那么就將v(t)-h(t)結果當做第(2)步的輸入信號，繼續(2)~(6)步迭代計算，提取新的多元IMF分量h(t)；否則，將h(t)當做第(2)步的輸入信號，繼續執行(2)~(6)步迭代，當滿足終止條件時停止。
經過一系列分解過程，原n元信號{v(t)}t=1T={v1(t), v2(r), …, vn(t)}被分解為一系列IMF{(hi(t)}i=1q和余量r(t)的加和形式為
| $v\left(t \right)=\sum\limits_{i=1}^q {{h_i}\left(t \right) + r\left(t \right), } $ |
式中q為分解出的IMF層數，hi(t)為{hi1(t), hi2(t), …, hin(t)}t=1T，r(t)為{r1(t), r2(t), …, rn(t)}t=1T分別對應于n元信號的n組IMF分量和n個余量。
2 基于MEMD腦信號特征提取
EEG信號或MEG信號一般是由多電極測得的多通道信號，即為多元信號，各通道信號之間明顯存在著互信息[11]。用MEMD算法將多通道腦信號進行同步處理，可實現對各通道信號不同特征尺度上互信息及自信息的比較。
BCI CompetitionⅣ：data sets 3中10通道MEG信號的MEMD分解圖如圖 1所示，可見其各通道分解的IMF層數相同，通道間同頻率成分的波形對齊且存在差異。每個IMF都代表原信號一定的特征尺度，各通道在同一頻率尺度的特征波形相對應，反應了多尺度特征相關信息。這為進一步特征提取提供了支持。
本文提出對腦信號經MEMD分解后的IMF分量提取功率特征，該方法計算每個IMF分量的直接功率譜最大幅值對應頻率附近某一確定頻率區間的功率譜積分值，以該值與所有通道的IMF分量譜積分值和的比值作為特征。我們把這種特征叫做“區間最大功率”特征，公式描述為
| $F=\int\limits_{{f_{\max }}-a}^{{f_{\max }} + b} {P\left(f \right){\rm{d}}f/\sum\limits_{c=1}^C {\sum\limits_{i=1}^I {\int\limits_{f_{\max }^{C, i}-a}^{f_{\max }^{C, i} + b} {P\left(f \right){\rm{d}}f} } } } $ |
式中F為特征值，fmax為功率譜最大幅值對應的頻率，C為多通道信號通道數，I為分解的IMF層數，P(f)為功率譜密度，a、b值確定積分區間。若a、b值均為零，表示功率譜最大幅值為特征。
3 實驗分析
3.1 數據描述
文中采用腦機接口競賽BCI CompetitionⅢ：data setⅠ和BCI CompetitionⅣ：data sets 3的競賽數據，對提出的基于MEMD和功率特征結合的特
Figure1.
Diagram of MEMD decomposition of 10-channel MEG signal
征提取方法進行了檢驗。下面簡要說明兩組數據。
BCI CompetitionⅢ：data setⅠ數據記錄了一個實驗者運動想象的皮層腦電圖信號(ECoG)。每個實驗(trial)是一次想象左小指運動或者舌頭運動，記錄時間為3 s。訓練數據和測試數據采集間隔了大約一周的時間。數據格式為三維矩陣形式，即trial個數×通道數×采樣點數，訓練數據為278×64×3 000，測試數據為100×64×3 000。兩類運動的標簽分別為-1、1。
BCI CompetitionⅣ：data sets 3包含兩個實驗者S1和S2的MEG信號，記錄了實驗者右手腕在前、后、左、右四個方向上運動時的MEG變化信號。S1訓練數據格式為160×10×400，測試數據為74×10×400；S2訓練數據格式為160×10×400，測試數據為73×10×400。右、前、左、后四類運動的類別標簽分別為1、2、3、4。
3.2 實驗及結果分析
對于BCI CompetitionⅢ：data setⅠ數據，由于通道數較多，根據文獻[12]中的最優選取通道結果，本文選取的6個通道，編號為[12 26 38 48 57 59]。對數據進行8~30 Hz帶通濾波，去除干擾。對BCI CompetitionⅣ：data sets 3的S1、S2數據選取其全部十個通道，并分別進行0.5~30 Hz、0.5~8 Hz帶通濾波預處理。之后對3組數據利用本文提出的MEMD結合功率譜方法提取特征并降維。
MEMD分解得到多通道IMF分量，由于各trial分解IMF的層數并不全相同，所以，不能直接對每個IMF提取特征，否則造成各trial間特征向量長度不匹配。這里，采用前向截取規則均等化IMF分量層數，即保留每個trial分解的IMF的前m層IMF。若存在某些trial分解IMF分量層數小于m的情況，則對該分解結果余量后補0向量層，使其達到m層。
BCI CompetitionⅢ：data setⅠ數據每個trial用MEMD處理后，訓練數據和測試數據每個trial分解得到IMF層數不同，層數范圍均為8~12層。截齲嚎個trial的前6層IMF，并用區間最大功率[式(3)，令a=5, b=5]提取特征，這樣得到36維特征向量，然后主成分分析實現特征降維。BCI CompetitionⅣ：data sets 3數據S1、S2的IMF層數范圍分別為6~11層和4~8層，分別對S1、S2截取前7層和前3層IMF，然后對IMF分量提取區間最大功率特征[式(3)，其中a=0，b=0]，再進行主成分分析降維處理。最后，三組數據均使用線性判別分析分類器實現分類。最終識別率結果和競賽前三名的正確率及特征提取方法作比較如表 1、2所示。
表1
基于BCI CompetitionⅢ：data setⅠ數據的測試集識別正確率及方法對比
Table1.
Comparison of accuracies on test sets for CompetitionⅢdata setⅠ
表2
基于BCI CompetitionⅣ：data sets 3數據的測試集識別正確率及方法對比
Table2.
Comparison of accuracies on test sets for CompetitionⅣdata sets 3
由表 1、2中實驗結果可以看出，本文提出的基于MEMD分解的功率特征提取方法取得了良好的效果。對BCI CompetitionⅢ：data setⅠ數據的處理結果正確率達到92.0%，而競賽第一名正確率為91.0%；對BCI CompetitionⅣ：data sets 3的處理結果同樣表現良好，S1和S2的平均正確率為46.2%，達到了競賽優異水平。MEMD結合區間最大功率的特征提取方法簡單、有效、穩定。這表明本文特征提取方法能夠較好地提取腦信號特征，實現對EEG信號和EMG信號的分析。
3.3 實驗討論
進一步研究了MEMD算法參數改變對ECoG、EMG信號特征提取算法的影響，其他參數均不改變。MEMD算法中一個重要參數就是空間方向向量的個數，即均值求解所需的超球面采樣點數。一般來說，方向向量的個數不能小于信號通道數的兩倍，且應遠大于通道數[9]。MEMD算法空間方向向量參數變化對分類系統的影響結果如表 3所示。從表 3可以看到，當方向向量個數為256時，各數據的識別正確率均為最高，過多或過少的方向向量都會導致識別正確率下降。
表3
MEMD算法空間方向向量個數改變對識別正確率影響
Table3.
Influence of changes of the number of direction vectors in MEMD on recognition accuracy
此外，我們對MEMD后的IMF分量提取了其他特征包括AR系數、多元AR系數、帶通功率、時域能量等，效果都不如區間最大功率特征；之后，文獻[13]中利用CSP提取特征，實驗結果也較好，而本文提取功率特征，更為簡單有效。對特征也使用了線性支持向量機、感知器、最近鄰等分類器，效果都不如線性判別分析。同時，為保證每個trial的IMF分量數量一致性，我們實驗了幾種自己研究的截取規則：后向截取、IMF能量較小層截取、IMF與源信號相關系數較小層截取等規則。結果均不如前向截取規則理想。說明前向規則雖然簡單，但處理結果能較好地包含數據的有用信息，利于系統識別分類。
4 結論
BCI系統中，特征提取過程十分關鍵。本文提出了一種基于MEMD結合區間最大功率特征提取算法，將多通道腦信號實現了多尺度特征校準分析，得到多元多尺度特征，并將該方法用于BCI CompetitionⅢ：data setⅠECoG數據和BCI CompetitionⅣ：data sets 3 EMG信號數據提取特征。實驗結果均達到競賽參賽優勝組水平。本文算法優點是簡單，對分析非平穩、非線性、成分復雜的腦信號有較高的有效性和較好的穩定性，為BCI腦信號處理提供了新的研究方向。下一步工作重點是根據腦信號特性改進MEMD算法，開展MEMD結合其他特征的提取方法。
引言
腦-機接口(brain-computer interface, BCI)系統能實現利用人腦意識控制外設[1]。BCI需要提取反映腦信號意識任務的特征，常見的特征提取方法有時域能量分析、AR模型系數、小波變換系數、共空間模式等。
經驗模態分解(empirical mode decomposition, EMD)算法[2]已成功用于腦電(electroencephalogram, EEG)信號處理[3]。文獻[4]利用EMD與AR模型結合的方法提取腦磁(magnetoencephalogram, MEG)信號特征。標準EMD算法分析多通道信號時，往往需要分別分析每一個通道信號或者先將多通道信號整合為單通道信號后再分析，忽略了腦信號的空間特性，不利于各通道間同一頻率尺度對比分析。因此標準EMD不適合直接處理多通道信號[5]。于是EMD的多元擴展方法不斷被提出，例如復數EMD[6]、二元EMD[7]、三元EMD[8]、多元EMD(multivariate EMD, MEMD)[9]。MEMD利用多維空間超球面方向向量求均值包絡，實現不同頻率尺度內各通道信號的共同模式的分析。
本文提出了一種基于MEMD與功率特征相結合的BCI特征提取算法。首先對多通道腦信號進行MEMD，然后提取IMF分量的功率特征并進行主成分分析(principal component analysis, PCA)降維，最后利用線性判別分析進行分類。腦機接口競賽BCI CompetitionⅢ和BCI CompetitionⅣ中的兩組數據的實驗結果令人滿意。
1 MEMD算法
MEMD實現了多元信號的多通道同步聯合分析。多元信號每個通道分解出的IMF個數相同，各層IMF頻率不同，先分解出的IMF層頻率高，后分解出的IMF頻率低。每個通道對應的IMF按頻率尺度對齊，形成多元IMF。
MEMD算法關鍵是求解多元信號的局域均值，這需要在多維空間建立方向向量集，而合適向量集的問題就轉化為在超球面找合適采樣點集的問題。
設一個n維向量組序列{v(t)}t-1T={v1(t), v2(t), …, vn(t)}，代表一個n元信號，長度為T，xθk={x1k, x2k, …, xnk, }表示在(n-1)維球面上對應角θk={θ1k, θ2k, …, θn-1k}的方向向量集, k=1, 2, …, k。MEMD算法的步驟總結如下：
(1)采用Hammersley序列采樣法，在(n-1)維球面上獲得合適的均勻采樣點集，即得到n維空間的方向向量。
(2)計算輸入信號v(t)在每個方向向量xθk上的映射pθk(t)。
(3)確定所有方向向量的映射信號{pθk(t)}k=1K極值對應的瞬時時刻{tlθk}k=1K，l表示極值點位置，l∈[1, T]。
(4)用多元樣條插值函數插值極值點[tlθk, v(tlθk)]，得到K個多元包絡{eθk(t)}k=1K。
(5)對球面空間K個方向向量，n元信號均值m(t)為：
| $m\left(t \right)=\frac{1}{K}\sum\limits_{k=1}^K {{e^{{\theta _k}}}\left(t \right)} $ |
(6)通過h(t)=v(t)-m(t)提取固有模態函數h(t)，如果h(t)滿足多元IMF判斷標準[10]，那么就將v(t)-h(t)結果當做第(2)步的輸入信號，繼續(2)~(6)步迭代計算，提取新的多元IMF分量h(t)；否則，將h(t)當做第(2)步的輸入信號，繼續執行(2)~(6)步迭代，當滿足終止條件時停止。
經過一系列分解過程，原n元信號{v(t)}t=1T={v1(t), v2(r), …, vn(t)}被分解為一系列IMF{(hi(t)}i=1q和余量r(t)的加和形式為
| $v\left(t \right)=\sum\limits_{i=1}^q {{h_i}\left(t \right) + r\left(t \right), } $ |
式中q為分解出的IMF層數，hi(t)為{hi1(t), hi2(t), …, hin(t)}t=1T，r(t)為{r1(t), r2(t), …, rn(t)}t=1T分別對應于n元信號的n組IMF分量和n個余量。
2 基于MEMD腦信號特征提取
EEG信號或MEG信號一般是由多電極測得的多通道信號，即為多元信號，各通道信號之間明顯存在著互信息[11]。用MEMD算法將多通道腦信號進行同步處理，可實現對各通道信號不同特征尺度上互信息及自信息的比較。
BCI CompetitionⅣ：data sets 3中10通道MEG信號的MEMD分解圖如圖 1所示，可見其各通道分解的IMF層數相同，通道間同頻率成分的波形對齊且存在差異。每個IMF都代表原信號一定的特征尺度，各通道在同一頻率尺度的特征波形相對應，反應了多尺度特征相關信息。這為進一步特征提取提供了支持。
本文提出對腦信號經MEMD分解后的IMF分量提取功率特征，該方法計算每個IMF分量的直接功率譜最大幅值對應頻率附近某一確定頻率區間的功率譜積分值，以該值與所有通道的IMF分量譜積分值和的比值作為特征。我們把這種特征叫做“區間最大功率”特征，公式描述為
| $F=\int\limits_{{f_{\max }}-a}^{{f_{\max }} + b} {P\left(f \right){\rm{d}}f/\sum\limits_{c=1}^C {\sum\limits_{i=1}^I {\int\limits_{f_{\max }^{C, i}-a}^{f_{\max }^{C, i} + b} {P\left(f \right){\rm{d}}f} } } } $ |
式中F為特征值，fmax為功率譜最大幅值對應的頻率，C為多通道信號通道數，I為分解的IMF層數，P(f)為功率譜密度，a、b值確定積分區間。若a、b值均為零，表示功率譜最大幅值為特征。
3 實驗分析
3.1 數據描述
文中采用腦機接口競賽BCI CompetitionⅢ：data setⅠ和BCI CompetitionⅣ：data sets 3的競賽數據，對提出的基于MEMD和功率特征結合的特
Figure1.
Diagram of MEMD decomposition of 10-channel MEG signal
征提取方法進行了檢驗。下面簡要說明兩組數據。
BCI CompetitionⅢ：data setⅠ數據記錄了一個實驗者運動想象的皮層腦電圖信號(ECoG)。每個實驗(trial)是一次想象左小指運動或者舌頭運動，記錄時間為3 s。訓練數據和測試數據采集間隔了大約一周的時間。數據格式為三維矩陣形式，即trial個數×通道數×采樣點數，訓練數據為278×64×3 000，測試數據為100×64×3 000。兩類運動的標簽分別為-1、1。
BCI CompetitionⅣ：data sets 3包含兩個實驗者S1和S2的MEG信號，記錄了實驗者右手腕在前、后、左、右四個方向上運動時的MEG變化信號。S1訓練數據格式為160×10×400，測試數據為74×10×400；S2訓練數據格式為160×10×400，測試數據為73×10×400。右、前、左、后四類運動的類別標簽分別為1、2、3、4。
3.2 實驗及結果分析
對于BCI CompetitionⅢ：data setⅠ數據，由于通道數較多，根據文獻[12]中的最優選取通道結果，本文選取的6個通道，編號為[12 26 38 48 57 59]。對數據進行8~30 Hz帶通濾波，去除干擾。對BCI CompetitionⅣ：data sets 3的S1、S2數據選取其全部十個通道，并分別進行0.5~30 Hz、0.5~8 Hz帶通濾波預處理。之后對3組數據利用本文提出的MEMD結合功率譜方法提取特征并降維。
MEMD分解得到多通道IMF分量，由于各trial分解IMF的層數并不全相同，所以，不能直接對每個IMF提取特征，否則造成各trial間特征向量長度不匹配。這里，采用前向截取規則均等化IMF分量層數，即保留每個trial分解的IMF的前m層IMF。若存在某些trial分解IMF分量層數小于m的情況，則對該分解結果余量后補0向量層，使其達到m層。
BCI CompetitionⅢ：data setⅠ數據每個trial用MEMD處理后，訓練數據和測試數據每個trial分解得到IMF層數不同，層數范圍均為8~12層。截齲嚎個trial的前6層IMF，并用區間最大功率[式(3)，令a=5, b=5]提取特征，這樣得到36維特征向量，然后主成分分析實現特征降維。BCI CompetitionⅣ：data sets 3數據S1、S2的IMF層數范圍分別為6~11層和4~8層，分別對S1、S2截取前7層和前3層IMF，然后對IMF分量提取區間最大功率特征[式(3)，其中a=0，b=0]，再進行主成分分析降維處理。最后，三組數據均使用線性判別分析分類器實現分類。最終識別率結果和競賽前三名的正確率及特征提取方法作比較如表 1、2所示。
表1
基于BCI CompetitionⅢ：data setⅠ數據的測試集識別正確率及方法對比
Table1.
Comparison of accuracies on test sets for CompetitionⅢdata setⅠ
表2
基于BCI CompetitionⅣ：data sets 3數據的測試集識別正確率及方法對比
Table2.
Comparison of accuracies on test sets for CompetitionⅣdata sets 3
由表 1、2中實驗結果可以看出，本文提出的基于MEMD分解的功率特征提取方法取得了良好的效果。對BCI CompetitionⅢ：data setⅠ數據的處理結果正確率達到92.0%，而競賽第一名正確率為91.0%；對BCI CompetitionⅣ：data sets 3的處理結果同樣表現良好，S1和S2的平均正確率為46.2%，達到了競賽優異水平。MEMD結合區間最大功率的特征提取方法簡單、有效、穩定。這表明本文特征提取方法能夠較好地提取腦信號特征，實現對EEG信號和EMG信號的分析。
3.3 實驗討論
進一步研究了MEMD算法參數改變對ECoG、EMG信號特征提取算法的影響，其他參數均不改變。MEMD算法中一個重要參數就是空間方向向量的個數，即均值求解所需的超球面采樣點數。一般來說，方向向量的個數不能小于信號通道數的兩倍，且應遠大于通道數[9]。MEMD算法空間方向向量參數變化對分類系統的影響結果如表 3所示。從表 3可以看到，當方向向量個數為256時，各數據的識別正確率均為最高，過多或過少的方向向量都會導致識別正確率下降。
表3
MEMD算法空間方向向量個數改變對識別正確率影響
Table3.
Influence of changes of the number of direction vectors in MEMD on recognition accuracy
此外，我們對MEMD后的IMF分量提取了其他特征包括AR系數、多元AR系數、帶通功率、時域能量等，效果都不如區間最大功率特征；之后，文獻[13]中利用CSP提取特征，實驗結果也較好，而本文提取功率特征，更為簡單有效。對特征也使用了線性支持向量機、感知器、最近鄰等分類器，效果都不如線性判別分析。同時，為保證每個trial的IMF分量數量一致性，我們實驗了幾種自己研究的截取規則：后向截取、IMF能量較小層截取、IMF與源信號相關系數較小層截取等規則。結果均不如前向截取規則理想。說明前向規則雖然簡單，但處理結果能較好地包含數據的有用信息，利于系統識別分類。
4 結論
BCI系統中，特征提取過程十分關鍵。本文提出了一種基于MEMD結合區間最大功率特征提取算法，將多通道腦信號實現了多尺度特征校準分析，得到多元多尺度特征，并將該方法用于BCI CompetitionⅢ：data setⅠECoG數據和BCI CompetitionⅣ：data sets 3 EMG信號數據提取特征。實驗結果均達到競賽參賽優勝組水平。本文算法優點是簡單，對分析非平穩、非線性、成分復雜的腦信號有較高的有效性和較好的穩定性，為BCI腦信號處理提供了新的研究方向。下一步工作重點是根據腦信號特性改進MEMD算法，開展MEMD結合其他特征的提取方法。
