腦機接口(BCI)可以直接通過腦電(EEG)信號控制外部設備。本文針對傳統主成分分析(PCA)和二維主成分分析(2DPCA)處理多通道EEG信號的局限性, 提出了多線性主成分分析(MPCA)的張量特征提取和分類框架。首先生成張量EEG數據, 然后進行張量降維并提取特征, 最后用Fisher線性判別分析分類器進行分類。實驗中將新方法應用到BCI competitionⅡ數據集4和BCI competitionⅣ數據集3, 分別使用了EEG數據的時空二階張量表示形式和時空頻三階張量表示形式, 通過對可調參數多次調試, 取得了高于其它同類降維方法的最佳結果。二階輸入最高正確率分別達到81.0%和40.1%, 三階輸入分別達到76.0%和43.5%。
引用本文： 王金甲, 楊亮. 腦機接口中多線性主成分分析的張量特征提取. 生物醫學工程學雜志, 2015, 32(3): 526-530. doi: 10.7507/1001-5515.20150096 復制
引言
腦機接口(brain computer interface, BCI)是一種新的人機交互方式，不依賴于大腦的正常輸出通路，可以通過腦電(electroencephalogram, EEG)信息控制外部設備，如機器人、拼寫設備、虛擬環境等[1]。這給肢體癱瘓患者帶來了希望。
腦電信號是非平穩性、非線性信號，反映了大腦規律性的活動，具有信號幅度微弱、信噪比低、噪聲干擾大等特點。大腦的工作過程極為復雜，信息傳遞要靠多個神經元共同協作來完成，因此BCI中多使用多通道腦電信號。傳統的BCI多通道腦電信號特征提取方法有時空頻特征提取，如AR模型、時頻濾波方法、共空間模式和主成分分析(principal component analysis, PCA)等[2]，以及降維特征提取，如線性判別分析方法等[3]。文獻[4]在人臉識別中對于圖像數據采用二維主成分分析(two-dimensional principal component analysis, 2DPCA)提取特征。Lu等[5]針對視頻步態識別問題，擴展了PCA和2DPCA，提出一個以數據張量表示為基礎的多線性主成分分析(multi-linear principal component analysis, MPCA)框架，有效地解決了高維數據降維的問題。
張量是向量和矩陣的擴展，向量可以看成是一階張量，矩陣即二階張量。從張量分析角度，多通道腦電信號可以被表示為多階張量形式[6]。一種張量表示是利用腦電數據的時空矩陣表示，另一種張量表示是將原始數據做小波分析，得到頻率分量后與時間、通道一起組成三階張量。作者提出將基于MPCA的張量特征提取框架應用到多通道腦電信號特征提取中，最后用Fisher線性判別分析分類器進行特征分類。
1 張量子空間映射
信號的張量表示見文獻[7]。多通道腦電信號表示為高階張量，可以嵌入到低階張量子空間。下面是n階張量x到張量子空間的投影公式：
| $ y=x{{\rm{ \times }}_1}{{{\boldsymbol{\rm{\tilde U}}}}^{\left(1 \right){\rm{T}}}}{{\rm{ \times }}_2}{{{\boldsymbol{\rm{\tilde U}}}}^{\left(2 \right){\rm{T}}}} \cdots {{\rm{ \times }}_{\rm{N}}}{{{\boldsymbol{\rm{\tilde U}}}}^{\left({\rm{N}} \right){\rm{T}}}} $ |
每一個n階張量通過(n)投影到一個低維的空間上，共需要N次投影。圖 1說明了三階張量A(5×4×3)到矩陣B(2×5)的一階投影過程，結果為2×4×3的張量，其中A中長度為5的一階張量通過B投影得到長度為2的一階張量。
圖1
張量一階投影示意圖
Figure1.
First-order projection of tensor
2 基于MPCA的腦電信號張量特征提取
2.1 多線性主成分分析
MPCA定義為張量空間RI1?RI2…?RIN到張量子空間RP1?RP2…?RPN的多線性變換，要求子空間維數Pn小于原空間維數In, 當Pn=In時為全投影。投影公式為
| $ {y_m}={x_m}{{\rm{ \times }}_1}{{{\boldsymbol{\rm{\tilde U}}}}^{\left(1 \right){\rm{T}}}}{{\rm{ \times }}_2}{{{\boldsymbol{\rm{\tilde U}}}}^{\left(2 \right){\rm{T}}}} \cdots {{\rm{ \times }}_{\rm{N}}}{{{\boldsymbol{\rm{\tilde U}}}}^{\left({\rm{N}} \right){\rm{T}}}}, m=1, \cdots, M, $ |
其中{(n)∈RIn×Pn, n=1, …, N}為映射矩陣，M為張量表示的樣本個數，N為張量階數，滿足N≥1。
用k(k=1, …, K)表示求投影矩陣時的迭代次數，比率表示n階張量投影后保留的特征比例，用P表示保留的特征個數，λin(n)*是第in個全投影矩陣的n階特征值，在以后的敘述中為了表示方便我們用testQ來代替Q(n)。設X為訓練樣本 x張量表示的對應矩陣，則散射矩陣為Φ(n)：
| $ \begin{array}{l} {\Phi ^{\left(n \right)}}=\sum\limits_{m=1}^M {\left({{\boldsymbol{X}_{m\left(n \right)}}-{{\boldsymbol{\bar X}}_{\left(n \right)}}} \right)\cdot } {{{\boldsymbol{\rm{\tilde U}}}}_{\Phi \left(n \right)}} \cdot {\boldsymbol{\rm{\tilde U}}}_{\Phi \left(n \right)}^T \cdot \left({{\boldsymbol{X}_{m\left(n \right)}}-} \right.\\ \; \;\; \;\; \;\; \;\; \;\; {\left.{{{\boldsymbol{\bar X}}_{\left(n \right)}}} \right)^{\rm{T}}} \end{array} $ |
完整的MPCA流程圖如圖 2所示。
圖2
MPCA流程圖
Figure2.
Flow chart of MPCA
2.2 腦電信號的MPCA特征提取和分類
MPCA算法可以推廣用于BCI系統的多通道腦電信號特征提取和分類。預處理說明如下：對于二階張量表示，直接在原始數據上操作。對于三階張量表示，要對腦電數據做小波分析以得到頻率分量，選用Morlet小波，其帶寬fb和中心頻率fc的選擇尤為重要。MPCA需要對張量表示的數據做均值化處理。特征提取說明如下：首先設定testQ值來決定降維程度，通過訓練張量樣本求出投影矩陣 (n)，然后將測試數據投影到矩陣 (n)上，達到降維目的，提取P個特征，對P的取值在實驗中進行討論。特征分類說明如下：特征向量影響分類器的分類性能，所以分類方法的選取非常重要。根據所選數據特點和多次實驗得出：選擇線性分類器比非線性分類器效果更好，最終決定選用Fisher線性判別分析分類器。
3 實驗步驟和結果分析
實驗中采用了兩組腦電數據集驗證了新方法，并且把識別結果與競賽數據識別結果、傳統PCA、2DPCA作了比較，證明了MPCA在BCI系統中作為特征提取方法的有效性。
3.1 實驗數據
數據一采用BCI CompetitionⅡ數據集4，標記為S0，S0共有兩類，分別代表受試者左手手指敲擊鍵盤和右手手指敲擊鍵盤產生信號；經整理后訓練次數為316次，測試次數選為100次，每次選擇28個通道，每次實驗選擇50個時間點。
數據二采用BCI CompetitionⅣ數據集3，共有四類，分別代表受試者使用右手手腕向前、后、左、右四個方向移動操縱桿產生的信號。對象S1和S2數據的訓練次數都是160次，測試次數分別為74次和73次，每次實驗中選擇10個通道和400個時間點。
3.2 實驗中各參數介紹
testQ表示投影后在子空間中保存原始維數的比例，testQ區間定為90~99。P表示實驗中特征選擇的個數，從投影后全部特征中選取P個送入分類器。對于同一數據中的不同testQ值，降維后的總特征數是不同的，所以對應P的范圍也不同。fb和fc分別表示實驗2里Morlet小波分析的帶寬參數和中心頻率，實驗中根據數據預處理后的特點將fc設為1，將fb從1~6進行調試。
3.3 實驗過程和結果分析
3.3.1 實驗1:腦電信號二階張量表示
(1) 參數testQ和P取不同值下的測試集和訓練集的分類錯誤率結果。最優參數下的實驗結果如表 1所示。
表1
在最優testQ和P值下數據集的錯誤率
Table1.
Error rate corresponding to optimal testQ and P
(2) 最佳參數對應的10次重復實驗的錯誤率結果如表 2所示。可以看出兩個數據集的結果很穩定，效果和效率突出。
表2
二階張量表示的10次重復實驗數據集錯誤率
Table2.
Error rate of ten times for 2 order tensor of EEG data
3.3.2 實驗2:三階腦電信號張量表示
對二階腦電信號進行Morlet復小波變換，得到頻率值作為第三階張量。數據一的帶通濾波頻帶為8~30 Hz，數據二的帶通濾波頻帶為0.5~8 Hz。
(1) 在fb=1~6時對不同testQ值、不同P取值的結果進行實驗，數據一S0的測試結果如圖 3所示，從圖中可以明顯看出，參數testQ對結果的影響是非常大的。比較六種情況，可看出當fb=4、testQ=95時效果最好，訓練集錯誤率達到29.4%，測試集錯誤率達到24.0%。
圖3
當fb=3或4時不同的testQ值S0對應的測試集錯誤率
Figure3.
Error rate corresponding to different testQ when fb=3 or 4
(2) 利用同上方法對數據S1和數據S2進行測試，結果表明fb=3、testQ=98時S1取得測試集的最低錯誤率50.0%，fb=5、testQ=99時S2取得測試集的最低錯誤率63.0%。
(3) 最佳參數對應的10次重復實驗的平均錯誤率結果如表 3所示。可以看出MPCA對于處理兩類或四類腦電數據穩定性都很好。
表3
三階張量表示的10次重復實驗數據集錯誤率
Table3.
Error rate of ten times for 3 order tensor of EEG data
3.4 實驗討論
經驗證，隨著參數fb和testQ等的變化，迭代次數對結果偶爾會有影響，但影響較小，且滿足迭代次數越少效果越好，運算時間越短，所以k為0和1。另外，若不同參數值下訓練集正確率相同，則選擇時]間較短的為最后結果。經過多次實驗，我們得出了每個數據對應的最高平均正確率。最后把兩組數據的競賽結果、PCA和MPCA的結果做比較，如表 4所示。BCI competitionⅡ數據集4的第一名采用了共空域子空間分解和Fisher線性判別分析結合方法，第二名采用數據挖掘特征和logistic回歸分類器結合方法。BCI competitionⅣ數據集3第一名采用了統計特征頻域特征和小波特征的融合，并用遺傳算法選擇優化特征，分類采用線性支持向量機和線性判別分析結合方法，第二名采用時域特征和頻域特征的融合，并用Fisher線性判別降維，分類器采用Fisher線性判別分析。
表4
數據一和數據二中MPCA、PCA和競賽勝者的正確率(%)比較
Table4.
Comparison of accuracy (%) among MPCA, PCA and competition winner
4 結論
本文以BCI系統中腦電信號特征提取和分類為應用背景，針對多通道腦電數據的二階時空張量表示和三階時空頻張量表示，采用了多線性主成分分析(二階張量時就是2DPCA)的張量特征提取方法。通過對各參數的多次調試，有效改善了分類效果。對于BCI competitionⅡ數據集4和Ⅳ數據集3，二階張量時新方法的測試集最高正確率分別達到81.0%和40.1%，明顯高于傳統PCA的61.0%和24.5%；三階張量時MPCA分別達到76.0%和43.5%，遠高于PCA的50.0%和24.5%，同時平均值要高于競賽第二名。相同參數下的MPCA測試結果具有正確率高、穩定性好、運算簡便等優點。進一步可繼續研究采用其它時頻分析方法的高階張量表示和特征提取方法。另外文獻[5, 8]的MPCA投影矩陣采用了矩陣特征分解，可研究如何利用張量分解。
引言
腦機接口(brain computer interface, BCI)是一種新的人機交互方式，不依賴于大腦的正常輸出通路，可以通過腦電(electroencephalogram, EEG)信息控制外部設備，如機器人、拼寫設備、虛擬環境等[1]。這給肢體癱瘓患者帶來了希望。
腦電信號是非平穩性、非線性信號，反映了大腦規律性的活動，具有信號幅度微弱、信噪比低、噪聲干擾大等特點。大腦的工作過程極為復雜，信息傳遞要靠多個神經元共同協作來完成，因此BCI中多使用多通道腦電信號。傳統的BCI多通道腦電信號特征提取方法有時空頻特征提取，如AR模型、時頻濾波方法、共空間模式和主成分分析(principal component analysis, PCA)等[2]，以及降維特征提取，如線性判別分析方法等[3]。文獻[4]在人臉識別中對于圖像數據采用二維主成分分析(two-dimensional principal component analysis, 2DPCA)提取特征。Lu等[5]針對視頻步態識別問題，擴展了PCA和2DPCA，提出一個以數據張量表示為基礎的多線性主成分分析(multi-linear principal component analysis, MPCA)框架，有效地解決了高維數據降維的問題。
張量是向量和矩陣的擴展，向量可以看成是一階張量，矩陣即二階張量。從張量分析角度，多通道腦電信號可以被表示為多階張量形式[6]。一種張量表示是利用腦電數據的時空矩陣表示，另一種張量表示是將原始數據做小波分析，得到頻率分量后與時間、通道一起組成三階張量。作者提出將基于MPCA的張量特征提取框架應用到多通道腦電信號特征提取中，最后用Fisher線性判別分析分類器進行特征分類。
1 張量子空間映射
信號的張量表示見文獻[7]。多通道腦電信號表示為高階張量，可以嵌入到低階張量子空間。下面是n階張量x到張量子空間的投影公式：
| $ y=x{{\rm{ \times }}_1}{{{\boldsymbol{\rm{\tilde U}}}}^{\left(1 \right){\rm{T}}}}{{\rm{ \times }}_2}{{{\boldsymbol{\rm{\tilde U}}}}^{\left(2 \right){\rm{T}}}} \cdots {{\rm{ \times }}_{\rm{N}}}{{{\boldsymbol{\rm{\tilde U}}}}^{\left({\rm{N}} \right){\rm{T}}}} $ |
每一個n階張量通過(n)投影到一個低維的空間上，共需要N次投影。圖 1說明了三階張量A(5×4×3)到矩陣B(2×5)的一階投影過程，結果為2×4×3的張量，其中A中長度為5的一階張量通過B投影得到長度為2的一階張量。
圖1
張量一階投影示意圖
Figure1.
First-order projection of tensor
2 基于MPCA的腦電信號張量特征提取
2.1 多線性主成分分析
MPCA定義為張量空間RI1?RI2…?RIN到張量子空間RP1?RP2…?RPN的多線性變換，要求子空間維數Pn小于原空間維數In, 當Pn=In時為全投影。投影公式為
| $ {y_m}={x_m}{{\rm{ \times }}_1}{{{\boldsymbol{\rm{\tilde U}}}}^{\left(1 \right){\rm{T}}}}{{\rm{ \times }}_2}{{{\boldsymbol{\rm{\tilde U}}}}^{\left(2 \right){\rm{T}}}} \cdots {{\rm{ \times }}_{\rm{N}}}{{{\boldsymbol{\rm{\tilde U}}}}^{\left({\rm{N}} \right){\rm{T}}}}, m=1, \cdots, M, $ |
其中{(n)∈RIn×Pn, n=1, …, N}為映射矩陣，M為張量表示的樣本個數，N為張量階數，滿足N≥1。
用k(k=1, …, K)表示求投影矩陣時的迭代次數，比率表示n階張量投影后保留的特征比例，用P表示保留的特征個數，λin(n)*是第in個全投影矩陣的n階特征值，在以后的敘述中為了表示方便我們用testQ來代替Q(n)。設X為訓練樣本 x張量表示的對應矩陣，則散射矩陣為Φ(n)：
| $ \begin{array}{l} {\Phi ^{\left(n \right)}}=\sum\limits_{m=1}^M {\left({{\boldsymbol{X}_{m\left(n \right)}}-{{\boldsymbol{\bar X}}_{\left(n \right)}}} \right)\cdot } {{{\boldsymbol{\rm{\tilde U}}}}_{\Phi \left(n \right)}} \cdot {\boldsymbol{\rm{\tilde U}}}_{\Phi \left(n \right)}^T \cdot \left({{\boldsymbol{X}_{m\left(n \right)}}-} \right.\\ \; \;\; \;\; \;\; \;\; \;\; {\left.{{{\boldsymbol{\bar X}}_{\left(n \right)}}} \right)^{\rm{T}}} \end{array} $ |
完整的MPCA流程圖如圖 2所示。
圖2
MPCA流程圖
Figure2.
Flow chart of MPCA
2.2 腦電信號的MPCA特征提取和分類
MPCA算法可以推廣用于BCI系統的多通道腦電信號特征提取和分類。預處理說明如下：對于二階張量表示，直接在原始數據上操作。對于三階張量表示，要對腦電數據做小波分析以得到頻率分量，選用Morlet小波，其帶寬fb和中心頻率fc的選擇尤為重要。MPCA需要對張量表示的數據做均值化處理。特征提取說明如下：首先設定testQ值來決定降維程度，通過訓練張量樣本求出投影矩陣 (n)，然后將測試數據投影到矩陣 (n)上，達到降維目的，提取P個特征，對P的取值在實驗中進行討論。特征分類說明如下：特征向量影響分類器的分類性能，所以分類方法的選取非常重要。根據所選數據特點和多次實驗得出：選擇線性分類器比非線性分類器效果更好，最終決定選用Fisher線性判別分析分類器。
3 實驗步驟和結果分析
實驗中采用了兩組腦電數據集驗證了新方法，并且把識別結果與競賽數據識別結果、傳統PCA、2DPCA作了比較，證明了MPCA在BCI系統中作為特征提取方法的有效性。
3.1 實驗數據
數據一采用BCI CompetitionⅡ數據集4，標記為S0，S0共有兩類，分別代表受試者左手手指敲擊鍵盤和右手手指敲擊鍵盤產生信號；經整理后訓練次數為316次，測試次數選為100次，每次選擇28個通道，每次實驗選擇50個時間點。
數據二采用BCI CompetitionⅣ數據集3，共有四類，分別代表受試者使用右手手腕向前、后、左、右四個方向移動操縱桿產生的信號。對象S1和S2數據的訓練次數都是160次，測試次數分別為74次和73次，每次實驗中選擇10個通道和400個時間點。
3.2 實驗中各參數介紹
testQ表示投影后在子空間中保存原始維數的比例，testQ區間定為90~99。P表示實驗中特征選擇的個數，從投影后全部特征中選取P個送入分類器。對于同一數據中的不同testQ值，降維后的總特征數是不同的，所以對應P的范圍也不同。fb和fc分別表示實驗2里Morlet小波分析的帶寬參數和中心頻率，實驗中根據數據預處理后的特點將fc設為1，將fb從1~6進行調試。
3.3 實驗過程和結果分析
3.3.1 實驗1:腦電信號二階張量表示
(1) 參數testQ和P取不同值下的測試集和訓練集的分類錯誤率結果。最優參數下的實驗結果如表 1所示。
表1
在最優testQ和P值下數據集的錯誤率
Table1.
Error rate corresponding to optimal testQ and P
(2) 最佳參數對應的10次重復實驗的錯誤率結果如表 2所示。可以看出兩個數據集的結果很穩定，效果和效率突出。
表2
二階張量表示的10次重復實驗數據集錯誤率
Table2.
Error rate of ten times for 2 order tensor of EEG data
3.3.2 實驗2:三階腦電信號張量表示
對二階腦電信號進行Morlet復小波變換，得到頻率值作為第三階張量。數據一的帶通濾波頻帶為8~30 Hz，數據二的帶通濾波頻帶為0.5~8 Hz。
(1) 在fb=1~6時對不同testQ值、不同P取值的結果進行實驗，數據一S0的測試結果如圖 3所示，從圖中可以明顯看出，參數testQ對結果的影響是非常大的。比較六種情況，可看出當fb=4、testQ=95時效果最好，訓練集錯誤率達到29.4%，測試集錯誤率達到24.0%。
圖3
當fb=3或4時不同的testQ值S0對應的測試集錯誤率
Figure3.
Error rate corresponding to different testQ when fb=3 or 4
(2) 利用同上方法對數據S1和數據S2進行測試，結果表明fb=3、testQ=98時S1取得測試集的最低錯誤率50.0%，fb=5、testQ=99時S2取得測試集的最低錯誤率63.0%。
(3) 最佳參數對應的10次重復實驗的平均錯誤率結果如表 3所示。可以看出MPCA對于處理兩類或四類腦電數據穩定性都很好。
表3
三階張量表示的10次重復實驗數據集錯誤率
Table3.
Error rate of ten times for 3 order tensor of EEG data
3.4 實驗討論
經驗證，隨著參數fb和testQ等的變化，迭代次數對結果偶爾會有影響，但影響較小，且滿足迭代次數越少效果越好，運算時間越短，所以k為0和1。另外，若不同參數值下訓練集正確率相同，則選擇時]間較短的為最后結果。經過多次實驗，我們得出了每個數據對應的最高平均正確率。最后把兩組數據的競賽結果、PCA和MPCA的結果做比較，如表 4所示。BCI competitionⅡ數據集4的第一名采用了共空域子空間分解和Fisher線性判別分析結合方法，第二名采用數據挖掘特征和logistic回歸分類器結合方法。BCI competitionⅣ數據集3第一名采用了統計特征頻域特征和小波特征的融合，并用遺傳算法選擇優化特征，分類采用線性支持向量機和線性判別分析結合方法，第二名采用時域特征和頻域特征的融合，并用Fisher線性判別降維，分類器采用Fisher線性判別分析。
表4
數據一和數據二中MPCA、PCA和競賽勝者的正確率(%)比較
Table4.
Comparison of accuracy (%) among MPCA, PCA and competition winner
4 結論
本文以BCI系統中腦電信號特征提取和分類為應用背景，針對多通道腦電數據的二階時空張量表示和三階時空頻張量表示，采用了多線性主成分分析(二階張量時就是2DPCA)的張量特征提取方法。通過對各參數的多次調試，有效改善了分類效果。對于BCI competitionⅡ數據集4和Ⅳ數據集3，二階張量時新方法的測試集最高正確率分別達到81.0%和40.1%，明顯高于傳統PCA的61.0%和24.5%；三階張量時MPCA分別達到76.0%和43.5%，遠高于PCA的50.0%和24.5%，同時平均值要高于競賽第二名。相同參數下的MPCA測試結果具有正確率高、穩定性好、運算簡便等優點。進一步可繼續研究采用其它時頻分析方法的高階張量表示和特征提取方法。另外文獻[5, 8]的MPCA投影矩陣采用了矩陣特征分解，可研究如何利用張量分解。
