眼電偽跡和噪聲是導致腦電信號低信噪比的重要原因, 會降低運動想象任務的分類性能。提出一種改進的基于少通道數的分塊欠定盲源分離的濾波方法, 通過分塊的思想把非平穩的腦電信號變為近平穩的分塊信號, 利用二階欠定混合矩陣盲識別方法估計混合分離矩陣, 然后通過基于最小均方誤差的波速形成器提取源信號, 接著通過得分準則自動去除噪聲信號并重構信號, 最后提取共空間模式特征進行分類。想象運動的真實腦電信號實驗仿真結果表明, 分塊欠定盲源分離方法能很好地恢復源信號并能有效地去除眼電等偽跡和噪聲, 共空間模式特征則提高了想象任務識別率。
引用本文： 盧輝斌. 欠定盲源分離和共空間模式特征的腦電信號分類研究. 生物醫學工程學雜志, 2016, 33(2): 216-220. doi: 10.7507/1001-5515.20160038 復制
引言
大腦運動想象時能自發地產生腦電信號(electroencephalogram, EEG)，在大腦兩側皮層區域會出現mu節律和beta節律的事件相關去同步化和同步化[1]。共空間模式(common spatial patterns, CSP)及其改進算法[1-2]利用腦電信號這一特性提取兩種類別想象運動任務的特征，在腦機接口(brain computer interface, BCI)的應用中十分成功。觀測腦電信號不僅包括想象運動源信號，也常常包括與想象運動無關的其他活動源產生的電信號，如肌電圖(electromyogram, EMG)、心電圖(electrocardiogram, ECG)、眼電圖(electrooculogram, EOG)等偽跡和噪聲。這些偽跡強度可能比運動想象源信號更強，使得腦電信號的信噪比很低。腦電信號去偽跡和噪聲方法是近年熱門研究方向, 如小波變換、卡爾曼濾波、自回歸模型、獨立分量分析(independent component analysis, ICA)、經驗模態分解(empirical mode decomposition, EMD)等[3-5]。現有的小波分析、卡爾曼濾波器和自回歸模型對腦電信號的一般噪聲去除效果不錯，但是一旦遇到像眼電圖這樣的頻率成分和范圍與腦電信號重疊且幅度大的偽跡則效果不佳。基于高階統計量理論的ICA及其擴展方法需要手動去除與任務無關的偽跡。EMD及其擴展方法在腦電信號去噪應用時效果較好但是十分耗時。De Lathauwer等[6-7]提出了一種欠定混合二階盲辨識(second order blind identification of underdetermined mixtures, SOBIUM)的統計信號處理方法，在語音信號盲識別中取得良好的效果, 能精確地恢復原始語音信號且去除與語音信號無關的干擾信號。這為欠定盲源分離應用腦電信號去噪提供了可能。
為解決現有的去偽跡方法的缺點，基于SOBIUM算法，本文提出僅利用少量通道提取腦電信號特征的新方法，即欠定盲源分離結合CSP特征。首先把非平穩的腦電信號進行分塊得到分塊腦電信號，使用SOBIUM算法得到盲源分離模型中的混合矩陣，并使用列文伯格-麥夸特(Levenberg-Marquardt, LM)算法[8-9]對混合矩陣進行迭代優化，改進了SOBIUM算法中混合矩陣的估計精確度。其次考慮輸出信號能量為非負的約束條件改進了最小均方誤差(minimum mean square error, MMSE)波束形成器，并恢復出源信號，接著利用fisher得分準則自動去除與運動想象任務無關的源信號即偽跡。最后重構基于運動想象相關的腦電信號，提取CSP特征并分類。
1 欠定盲源分離模型
基于盲源分離模型，觀測腦電信號可表示為
| $ \boldsymbol{X}=\boldsymbol{A}\boldsymbol{S} $ |
其中X為m行N列的觀測腦電信號矩陣，m為通道數，N為時間采樣點數。S為d個獨立源信號的d行N列源信號矩陣。A為m行d列的混合矩陣。本文假設腦電信號盲源分離問題是欠定的，即觀測通道數小于源信號數m < d。
本文提出的整體思路如下：首先采用SOBIUM方法估計初始混合矩陣A，改進準則函數并采用LM算法優化A。其次基于穩定性條件改進了MMSE波束形成器，恢復出源信號S。再次計算混合矩陣A的所有列fisher得分，我們認為高fisher得分的列對應想象運動的源信號SM, 低fisher得分的列對應想象運動任務不相關的源信號。最后基于混合矩陣A和想象運動的源信號SM重構腦電信號，進行CSP特征提取和分類。本文方法的總體框圖如下。
圖1
欠定盲源分離結合CSP的特征提取框圖
Figure1.
Feature extraction diagram using the under-determined blind identification and CSP
2 估計混合矩陣A
觀測腦電信號X是非平穩信號，我們首先將其在時間上順序分成等長度的M塊，每塊腦電信號近似看作零均值和常量方差的信號。每塊腦電信號的長度為N1，用XK表示第k塊腦電信號, k=1, …, M，第k塊的第f個源信號的方差用Dkf表示，f=1, …, d，所有Dkf組成源信號方差矩陣D，顯然D中所有元素為非負值。
第k塊腦電信號Xk的估計協方差矩陣為，把M塊腦電信號協方差矩陣堆疊成一個(m, m, M)維的三階張量，其元素為，i, j=1, …, M，其中θ表示由混合矩陣A和信號方差矩陣D所有元素組成的參數向量。
由于每塊腦電信號的協方差矩陣是對稱的，因此自由張量元素個數為Mm(m+1)/2，分解后自由參數個數為d(M+m-1)。傳統張量分解方法利用估計張量和理論張量的均方差最小化準則，通過計算可求出混合矩陣A和方差矩陣D。但該方法沒有考慮張量元素間的相關性。而M塊腦電信號協方差矩陣疊成的張量的元素間是相關的。因此我們使用極大似然原理導出張量分解的最優準則[9]取代傳統張量分解準則。腦電信號的聯合似然函數為
| $ \begin{array}{*{20}{c}} {f\left(\boldsymbol{X} \right)=\mathop \boldsymbol{\Pi} \limits_{k=1}^M {{\left({2\pi \det {\boldsymbol{R}_k}} \right)}^{-{N_1}/2}}\exp \left[ {-\frac{1}{2}\sum\limits_{t=1}^{{N_1}} {\boldsymbol{X}_{k, :, t}^T} } \right.}\ {\left. {\boldsymbol{R}_k^{-1}{\boldsymbol{X}_{k, :, t}}} \right]} \end{array} $ |
對式(2)等號兩邊分別取對數，對Rk在附近進行泰勒級數展開并代入整理得
| $-\log f\left(\boldsymbol{X} \right)\approx \frac{{{N}_{1}}}{4}\sum\limits_{k=1}^{M}{\left(\hat{\boldsymbol{R}}_{k}^{-1}\vartriangle {{\boldsymbol{R}}_{k}}\hat{\boldsymbol{R}}_{k}^{-1}\vartriangle {{\boldsymbol{R}}_{k}} \right.} $ |
其中ΔRk=Rk-。由式(3)可得到參數向量θ的等效極大似然估計的準則函數為
| $ Q\left(\theta \right)=\sum\limits_{k=1}^{M}{\text{tr}[{{\boldsymbol{C}}_{k}}\left({{{\hat{\boldsymbol{R}}}}_{k}}-{{\boldsymbol{R}}_{k}}\left(\theta \right) \right){{\boldsymbol{C}}_{k}}\left({{{\hat{\boldsymbol{R}}}}_{k}}-{{\boldsymbol{R}}_{k}}\left(\theta \right) \right)]} Q\left(\theta \right)=\sum\limits_{k=1}^M {{\rm{tr}}[{{\bf{C}}_k}\left({{{\widehat {\bf{R}}}_k}-{{\bf{R}}_k}\left(\theta \right)} \right){{\bf{C}}_k}\left({{{\widehat {\bf{R}}}_k}-{{\bf{R}}_k}\left(\theta \right)} \right)]} $ |
其中tr(·)表示矩陣跡，Rk是可逆的，是一個加權矩陣，正則化參數ε是一個適當小的正常數。利用SOBIUM算法[6-7]求解估計出初始混合矩陣A，然后估計出初始源信號方差矩陣D。
由于SOBIUM算法在估計欠定盲混合矩陣時沒有對源信號方差矩陣D的元素進行約束[10]，因此D的有些估計值可能是負值，這與方差定義不符。因此有必要將式(4)加以改進，減去一項閘函數使得方差矩陣D的所有元素為非負數。
| $ {Q^ * }\left(\theta \right)=Q\left(\theta \right)-\alpha \sum\limits_{f=1}^r {\left[ {\sum\limits_{k=1}^M {\log {\boldsymbol{D}_{kf}} + M\log \left({\sum\limits_{k=1}^d {\boldsymbol{A}_{kf}^2} } \right)} } \right]} $ |
其中α為初始值為1的正參數，在優化過程中逐漸減小最終趨于0，每5次迭代乘以一個小于1的因子。在優化最后，α應該接近于0，這樣閘函數對準則函數極小值沒有影響，從而閘函數只是用來約束D中元素為非負值。而且閘函數能夠平滑準則函數并去除一些局部極值。我們使用LM算法對式(5)進行最小化求解，步驟如下：把SOBIUM算法估計的初始混合矩陣A和源信號方差D作為LM算法初始化值；如果源信號方差矩陣D的元素為負值，則取其絕對值，以保證D的所有項都是非負的；最后利用LM算法進行迭代直到收斂。LM算法涉及海森矩陣求逆，計算量較大。
3 MMSE波束形成器分離源信號
MMSE波束形成器[11]是盲源分離算法中常用的估計源信號方法，是一個多輸入多輸出的時空濾波器。我們采用一個新的兩步穩定MMSE波束形成器，首先通過源信號和估計信號的理論均方距離的最小化尋找最佳波束形成器，最后考慮穩定性問題。MMSE波束形成器把腦電信號進行分塊處理, 其分塊數與SOBIUM算法的分塊數可以不一樣。N2表示分塊后的每塊腦電信號的有意義長度，即可認為在N2范圍內可以較準確地估計每塊腦電信號的自協方差。MMSE波束形成器自由參數P表示信號時域濾波器長度。
3.1 腦電信號協方差矩陣估計
Rt[τ]表示腦電信號在時間t時刻的理論時延協方差矩陣，是以t時刻為中心、塊長度為N1的數據樣本的Rt[τ]的古典樣本估計。已知混合矩陣和源信號的估計時，, 其中表示通過SOBIUM算法和LM算法優化后混合矩陣A的估計, 表示源信號方差矩陣Dt, :τ的估計。估計出、, 通過最小化以求出Dt, :τ, 其中‖·‖F是Frobenius范數，Ct=()-1表示加權矩陣。因此Dt, :τ最小二乘解為，其中，⊙表示Khatri-Rao積。
3.2 MMSE波束形成器
我們利用MMSE波束形成器估計源信號，波束形成器用一個長度為mP×1的矢量wj, t表示，第j個信號的第t時刻樣本的源信號估計為，其中，X:, tT表示腦電信號X的第t列，源信號和估計源信號的理論均方距離為，通過代數方法最小化上式可得波束形成器為
| $ \begin{array}{*{20}{c}} {{w_{j, t}}={\rm{E}}{{\left[ {{{\tilde \boldsymbol{X}}_{:, t}}\tilde \boldsymbol{X}_{:, t}^T} \right]}^{-1}}{\rm{E}}\left[ {{{\tilde \boldsymbol{X}}_{:, t}}{\boldsymbol{S}_{jt}}} \right]}\ {={\rm{E}}{{\left[ {{{\tilde \boldsymbol{X}}_{:, t}}\tilde \boldsymbol{X}_{:, t}^T} \right]}^{-1}}\left({{d_{j, t}} \otimes {\boldsymbol{A}_{:, j}}} \right)} \end{array} $ |
其中dj, t=[Dtj0, Dtj1, …, DtjP-1]T，它是t時刻第j個源信號的自協方差。期望計算公式為，這里T[.]表示是由參數塊組成的對稱塊Toeplitz矩陣。
3.3 波束形成器穩定性條件
MMSE波束形成器的輸出信號能量為非負, 否則波束形成器不穩定[12]。因此對于每一時刻t，我們采用兩步穩定的波束形成器。第一、當的元素為負值取其絕對值，以及其元素小于給定的閾值δ，則令該項為δ。第二、首先給定估計的自協方差，將其轉換成自協方差自回歸過程的相應系數；其次計算自回歸模型的極點，如果這些極點大于1，則取其反共軛值；再次當所有極點大于rmax(接近1的數)，對其進行標準化，取其量級為rmax；最后將自回歸模型轉換回自協方差函數，并且使等于這個自協方差函數。通過以上步驟可以確保改進的是正定的，從而能達到使波束形成器穩定目的。
4 實驗仿真和分析
4.1 實驗數據描述
本次實驗第一組實驗采用2005年柏林BCI競賽3數據集III_4a的數據al和ay(http://www.bbci.de/competition/)進行實驗仿真，其中數據ay是小樣本訓練集數據。受試者按照指示完成相應的想象運動任務：(R)右手，(F)右腳。第二組實驗數據采用2008年柏林BCI競賽4數據集I的數據進行實驗仿真，其中受試者a、f選擇左手和腳運動想象，而受試者b、g執行左手和右手運動想象。
4.2 實驗結果與分析
4.2.1 第一組實驗
在al和ay兩個BCI競賽數據中, 每組數據共有280次實驗，其中al數據的訓練集為224次實驗而測試集為56次實驗，ay數據的訓練集為28次實驗而測試集為252次實驗。首先每組腦電數據通過0.5~3.75 s時間滑動窗和8~30 Hz的帶通濾波, 其次使用本文提出的欠定盲源分離方法對腦電信號進行去偽跡，再次使用CSP提取特征，最后采用線性分類器進行分類。
對兩組數據al和ay分別進行實驗，選取的6個通道如表 1所示。自由參數設置如下：al數據的塊數為280，ay的塊數為560；而使用MMSE波束形成器進行源信號估計時，al的塊數為1 400，ay的塊數為1 120；波束形成器的時域濾波長度參數P取值范圍為2~30，步長為2。
表1
al和ay數據分別選取不同通道時識別率(%)
Table1.
Recognition rate (%) of different channels of al and ay data
從表 1可以看出, 運動想象任務中通道選取對識別率有很大的影響，這是因為不同運動想象任務和大腦運動區不同位置有密切關系。同時波束形成器的時間濾波長度參數P對想象運動任務的識別有一定的影響。另外即使沒有使用半監督分類方法，小訓練樣本集ay也得到較好識別。在使用少通道數下獲得很好的識別正確率，這說明通道數遠小于源信號數的假設是可行的。
本文方法和其它方法的識別率比較如表 2所示。第一種方法是沒有去偽跡直接采用帶通濾波器和CSP結合方法；第二種方法是在通道參數相同情況下FastICA和CSP結合方法；第三種方法是本文方法；第四種方法是競賽第一名采用CSP的方法，對于小樣本訓練數據ay，競賽第一名還使用了半監督學習方法，但半監督學習方法將要付出很大的訓練時間代價。
表2
本文方法和其它方法的識別率(%)比較
Table2.
Recognition rate (%) of different methods
通過第一組實驗仿真結果, 可以看到al數據能取得較高的正確率，正確率都在90%以上，尤其是當選取通道為C1、C2、C3、C4、C5、Cz時能100%識別出右手想象運動和右腳想象運動任務。小訓練樣本數據ay在沒有使用半監督學習方法下，我們的方法也能取得較理想的結果。我們的方法要比只使用CSP方法和FastICA+CSP方法的分類性能高。此外，本文提出的方法與競賽第一名進行了比較，對于al數據，本文的方法與競賽第一名一樣都能百分之百地予以識別；對于ay數據，由于競賽第一名針對小樣本訓練集ay使用半監督學習，識別率達到了97.6%。我們注意到競賽第一名選取的通道數都在20個以上，而我們的方法選取6個通道就能實現相當的識別率。
4.2.2 第二組實驗
本次試驗采用競賽4數據Ⅰ的四個受試者的腦電信號數據進行實驗仿真，每個實驗者進行200次試驗，把200次試驗隨機抽取為140次試驗作為訓練集，剩下的60次試驗作為測試集，一共進行十次隨機分配訓練集和測試集數據。實驗結果分別與連續小波變換(continuous wavelet transform, CWT)、同步擠壓小波變換(synchrosqueezed wavelet transform, SST)、EMD、2013年文獻[13]提出的多元經驗模式分解(multivariate empirical mode decomposition, MEMD)以及FastICA方法進行比較。文獻[13]選擇了FC3、FC4、Cz、C3、C4、C5、C6、T7、T8、CCP3和CCP4這11個通道的腦電信號作為實驗的處理信號。本文方法只選擇C3、C4、Cz、CP3、CP4、CPz六個通道的腦電信號進行實驗仿真。每個受試者的腦電數據在混合矩陣估計時分成200塊，而在使用MMSE波速形成器分離源信號時分成800塊。實驗結果如表 3所示。
表3
a、b、f和g數據采用不同方法的識別率(%)(表 3中對于受試者f和g來說，本文方法的識別率最高；對于受試者a和b來說，最高識別率為MEMD方法，居次位的是本文方法，識別率與MEMD方法差別不大。從4個受試者識別率的標準差來看本文方法穩定性最好。MEMD的11個通道多于本文方法6個通道，MEMD運行時間是本文方法的上百倍，計算復雜度高。
第二組實驗充分地證明了本文提出的方法在通道數較少時仍然取得很好的識別率，時間開銷較低，顯示了將該方法用于去除腦電信號偽跡的可行性。
5 結論
考慮到改善傳統去噪方法不能很好去除與腦電信號頻率重疊且幅度高的眼電等偽跡的問題，基于非平穩腦電信號分塊思想和SOBIUM算法，本文提出了欠定盲源分離結合CSP的特征提取新方法，首先基于欠定盲源分離提取腦電源信號，然后結合CSP特征進行分類，提高了BCI想象任務識別率，時間開銷很小。我們提出的改進之處如下：對SOBIUM算法得到的混合分離矩陣采用LM算法進行優化，提高混合矩陣盲識別的精確度；對MMSE波速形成器添加一些穩定性約束條件。新方法在電極通道數較少的情況下能夠提取腦電源信號，自動去除與想象運動任務無關的源信號。進一步研究將重構源信號和其它特征結合的特征提取方法。
引言
大腦運動想象時能自發地產生腦電信號(electroencephalogram, EEG)，在大腦兩側皮層區域會出現mu節律和beta節律的事件相關去同步化和同步化[1]。共空間模式(common spatial patterns, CSP)及其改進算法[1-2]利用腦電信號這一特性提取兩種類別想象運動任務的特征，在腦機接口(brain computer interface, BCI)的應用中十分成功。觀測腦電信號不僅包括想象運動源信號，也常常包括與想象運動無關的其他活動源產生的電信號，如肌電圖(electromyogram, EMG)、心電圖(electrocardiogram, ECG)、眼電圖(electrooculogram, EOG)等偽跡和噪聲。這些偽跡強度可能比運動想象源信號更強，使得腦電信號的信噪比很低。腦電信號去偽跡和噪聲方法是近年熱門研究方向, 如小波變換、卡爾曼濾波、自回歸模型、獨立分量分析(independent component analysis, ICA)、經驗模態分解(empirical mode decomposition, EMD)等[3-5]。現有的小波分析、卡爾曼濾波器和自回歸模型對腦電信號的一般噪聲去除效果不錯，但是一旦遇到像眼電圖這樣的頻率成分和范圍與腦電信號重疊且幅度大的偽跡則效果不佳。基于高階統計量理論的ICA及其擴展方法需要手動去除與任務無關的偽跡。EMD及其擴展方法在腦電信號去噪應用時效果較好但是十分耗時。De Lathauwer等[6-7]提出了一種欠定混合二階盲辨識(second order blind identification of underdetermined mixtures, SOBIUM)的統計信號處理方法，在語音信號盲識別中取得良好的效果, 能精確地恢復原始語音信號且去除與語音信號無關的干擾信號。這為欠定盲源分離應用腦電信號去噪提供了可能。
為解決現有的去偽跡方法的缺點，基于SOBIUM算法，本文提出僅利用少量通道提取腦電信號特征的新方法，即欠定盲源分離結合CSP特征。首先把非平穩的腦電信號進行分塊得到分塊腦電信號，使用SOBIUM算法得到盲源分離模型中的混合矩陣，并使用列文伯格-麥夸特(Levenberg-Marquardt, LM)算法[8-9]對混合矩陣進行迭代優化，改進了SOBIUM算法中混合矩陣的估計精確度。其次考慮輸出信號能量為非負的約束條件改進了最小均方誤差(minimum mean square error, MMSE)波束形成器，并恢復出源信號，接著利用fisher得分準則自動去除與運動想象任務無關的源信號即偽跡。最后重構基于運動想象相關的腦電信號，提取CSP特征并分類。
1 欠定盲源分離模型
基于盲源分離模型，觀測腦電信號可表示為
| $ \boldsymbol{X}=\boldsymbol{A}\boldsymbol{S} $ |
其中X為m行N列的觀測腦電信號矩陣，m為通道數，N為時間采樣點數。S為d個獨立源信號的d行N列源信號矩陣。A為m行d列的混合矩陣。本文假設腦電信號盲源分離問題是欠定的，即觀測通道數小于源信號數m < d。
本文提出的整體思路如下：首先采用SOBIUM方法估計初始混合矩陣A，改進準則函數并采用LM算法優化A。其次基于穩定性條件改進了MMSE波束形成器，恢復出源信號S。再次計算混合矩陣A的所有列fisher得分，我們認為高fisher得分的列對應想象運動的源信號SM, 低fisher得分的列對應想象運動任務不相關的源信號。最后基于混合矩陣A和想象運動的源信號SM重構腦電信號，進行CSP特征提取和分類。本文方法的總體框圖如下。
圖1
欠定盲源分離結合CSP的特征提取框圖
Figure1.
Feature extraction diagram using the under-determined blind identification and CSP
2 估計混合矩陣A
觀測腦電信號X是非平穩信號，我們首先將其在時間上順序分成等長度的M塊，每塊腦電信號近似看作零均值和常量方差的信號。每塊腦電信號的長度為N1，用XK表示第k塊腦電信號, k=1, …, M，第k塊的第f個源信號的方差用Dkf表示，f=1, …, d，所有Dkf組成源信號方差矩陣D，顯然D中所有元素為非負值。
第k塊腦電信號Xk的估計協方差矩陣為，把M塊腦電信號協方差矩陣堆疊成一個(m, m, M)維的三階張量，其元素為，i, j=1, …, M，其中θ表示由混合矩陣A和信號方差矩陣D所有元素組成的參數向量。
由于每塊腦電信號的協方差矩陣是對稱的，因此自由張量元素個數為Mm(m+1)/2，分解后自由參數個數為d(M+m-1)。傳統張量分解方法利用估計張量和理論張量的均方差最小化準則，通過計算可求出混合矩陣A和方差矩陣D。但該方法沒有考慮張量元素間的相關性。而M塊腦電信號協方差矩陣疊成的張量的元素間是相關的。因此我們使用極大似然原理導出張量分解的最優準則[9]取代傳統張量分解準則。腦電信號的聯合似然函數為
| $ \begin{array}{*{20}{c}} {f\left(\boldsymbol{X} \right)=\mathop \boldsymbol{\Pi} \limits_{k=1}^M {{\left({2\pi \det {\boldsymbol{R}_k}} \right)}^{-{N_1}/2}}\exp \left[ {-\frac{1}{2}\sum\limits_{t=1}^{{N_1}} {\boldsymbol{X}_{k, :, t}^T} } \right.}\ {\left. {\boldsymbol{R}_k^{-1}{\boldsymbol{X}_{k, :, t}}} \right]} \end{array} $ |
對式(2)等號兩邊分別取對數，對Rk在附近進行泰勒級數展開并代入整理得
| $-\log f\left(\boldsymbol{X} \right)\approx \frac{{{N}_{1}}}{4}\sum\limits_{k=1}^{M}{\left(\hat{\boldsymbol{R}}_{k}^{-1}\vartriangle {{\boldsymbol{R}}_{k}}\hat{\boldsymbol{R}}_{k}^{-1}\vartriangle {{\boldsymbol{R}}_{k}} \right.} $ |
其中ΔRk=Rk-。由式(3)可得到參數向量θ的等效極大似然估計的準則函數為
| $ Q\left(\theta \right)=\sum\limits_{k=1}^{M}{\text{tr}[{{\boldsymbol{C}}_{k}}\left({{{\hat{\boldsymbol{R}}}}_{k}}-{{\boldsymbol{R}}_{k}}\left(\theta \right) \right){{\boldsymbol{C}}_{k}}\left({{{\hat{\boldsymbol{R}}}}_{k}}-{{\boldsymbol{R}}_{k}}\left(\theta \right) \right)]} Q\left(\theta \right)=\sum\limits_{k=1}^M {{\rm{tr}}[{{\bf{C}}_k}\left({{{\widehat {\bf{R}}}_k}-{{\bf{R}}_k}\left(\theta \right)} \right){{\bf{C}}_k}\left({{{\widehat {\bf{R}}}_k}-{{\bf{R}}_k}\left(\theta \right)} \right)]} $ |
其中tr(·)表示矩陣跡，Rk是可逆的，是一個加權矩陣，正則化參數ε是一個適當小的正常數。利用SOBIUM算法[6-7]求解估計出初始混合矩陣A，然后估計出初始源信號方差矩陣D。
由于SOBIUM算法在估計欠定盲混合矩陣時沒有對源信號方差矩陣D的元素進行約束[10]，因此D的有些估計值可能是負值，這與方差定義不符。因此有必要將式(4)加以改進，減去一項閘函數使得方差矩陣D的所有元素為非負數。
| $ {Q^ * }\left(\theta \right)=Q\left(\theta \right)-\alpha \sum\limits_{f=1}^r {\left[ {\sum\limits_{k=1}^M {\log {\boldsymbol{D}_{kf}} + M\log \left({\sum\limits_{k=1}^d {\boldsymbol{A}_{kf}^2} } \right)} } \right]} $ |
其中α為初始值為1的正參數，在優化過程中逐漸減小最終趨于0，每5次迭代乘以一個小于1的因子。在優化最后，α應該接近于0，這樣閘函數對準則函數極小值沒有影響，從而閘函數只是用來約束D中元素為非負值。而且閘函數能夠平滑準則函數并去除一些局部極值。我們使用LM算法對式(5)進行最小化求解，步驟如下：把SOBIUM算法估計的初始混合矩陣A和源信號方差D作為LM算法初始化值；如果源信號方差矩陣D的元素為負值，則取其絕對值，以保證D的所有項都是非負的；最后利用LM算法進行迭代直到收斂。LM算法涉及海森矩陣求逆，計算量較大。
3 MMSE波束形成器分離源信號
MMSE波束形成器[11]是盲源分離算法中常用的估計源信號方法，是一個多輸入多輸出的時空濾波器。我們采用一個新的兩步穩定MMSE波束形成器，首先通過源信號和估計信號的理論均方距離的最小化尋找最佳波束形成器，最后考慮穩定性問題。MMSE波束形成器把腦電信號進行分塊處理, 其分塊數與SOBIUM算法的分塊數可以不一樣。N2表示分塊后的每塊腦電信號的有意義長度，即可認為在N2范圍內可以較準確地估計每塊腦電信號的自協方差。MMSE波束形成器自由參數P表示信號時域濾波器長度。
3.1 腦電信號協方差矩陣估計
Rt[τ]表示腦電信號在時間t時刻的理論時延協方差矩陣，是以t時刻為中心、塊長度為N1的數據樣本的Rt[τ]的古典樣本估計。已知混合矩陣和源信號的估計時，, 其中表示通過SOBIUM算法和LM算法優化后混合矩陣A的估計, 表示源信號方差矩陣Dt, :τ的估計。估計出、, 通過最小化以求出Dt, :τ, 其中‖·‖F是Frobenius范數，Ct=()-1表示加權矩陣。因此Dt, :τ最小二乘解為，其中，⊙表示Khatri-Rao積。
3.2 MMSE波束形成器
我們利用MMSE波束形成器估計源信號，波束形成器用一個長度為mP×1的矢量wj, t表示，第j個信號的第t時刻樣本的源信號估計為，其中，X:, tT表示腦電信號X的第t列，源信號和估計源信號的理論均方距離為，通過代數方法最小化上式可得波束形成器為
| $ \begin{array}{*{20}{c}} {{w_{j, t}}={\rm{E}}{{\left[ {{{\tilde \boldsymbol{X}}_{:, t}}\tilde \boldsymbol{X}_{:, t}^T} \right]}^{-1}}{\rm{E}}\left[ {{{\tilde \boldsymbol{X}}_{:, t}}{\boldsymbol{S}_{jt}}} \right]}\ {={\rm{E}}{{\left[ {{{\tilde \boldsymbol{X}}_{:, t}}\tilde \boldsymbol{X}_{:, t}^T} \right]}^{-1}}\left({{d_{j, t}} \otimes {\boldsymbol{A}_{:, j}}} \right)} \end{array} $ |
其中dj, t=[Dtj0, Dtj1, …, DtjP-1]T，它是t時刻第j個源信號的自協方差。期望計算公式為，這里T[.]表示是由參數塊組成的對稱塊Toeplitz矩陣。
3.3 波束形成器穩定性條件
MMSE波束形成器的輸出信號能量為非負, 否則波束形成器不穩定[12]。因此對于每一時刻t，我們采用兩步穩定的波束形成器。第一、當的元素為負值取其絕對值，以及其元素小于給定的閾值δ，則令該項為δ。第二、首先給定估計的自協方差，將其轉換成自協方差自回歸過程的相應系數；其次計算自回歸模型的極點，如果這些極點大于1，則取其反共軛值；再次當所有極點大于rmax(接近1的數)，對其進行標準化，取其量級為rmax；最后將自回歸模型轉換回自協方差函數，并且使等于這個自協方差函數。通過以上步驟可以確保改進的是正定的，從而能達到使波束形成器穩定目的。
4 實驗仿真和分析
4.1 實驗數據描述
本次實驗第一組實驗采用2005年柏林BCI競賽3數據集III_4a的數據al和ay(http://www.bbci.de/competition/)進行實驗仿真，其中數據ay是小樣本訓練集數據。受試者按照指示完成相應的想象運動任務：(R)右手，(F)右腳。第二組實驗數據采用2008年柏林BCI競賽4數據集I的數據進行實驗仿真，其中受試者a、f選擇左手和腳運動想象，而受試者b、g執行左手和右手運動想象。
4.2 實驗結果與分析
4.2.1 第一組實驗
在al和ay兩個BCI競賽數據中, 每組數據共有280次實驗，其中al數據的訓練集為224次實驗而測試集為56次實驗，ay數據的訓練集為28次實驗而測試集為252次實驗。首先每組腦電數據通過0.5~3.75 s時間滑動窗和8~30 Hz的帶通濾波, 其次使用本文提出的欠定盲源分離方法對腦電信號進行去偽跡，再次使用CSP提取特征，最后采用線性分類器進行分類。
對兩組數據al和ay分別進行實驗，選取的6個通道如表 1所示。自由參數設置如下：al數據的塊數為280，ay的塊數為560；而使用MMSE波束形成器進行源信號估計時，al的塊數為1 400，ay的塊數為1 120；波束形成器的時域濾波長度參數P取值范圍為2~30，步長為2。
表1
al和ay數據分別選取不同通道時識別率(%)
Table1.
Recognition rate (%) of different channels of al and ay data
從表 1可以看出, 運動想象任務中通道選取對識別率有很大的影響，這是因為不同運動想象任務和大腦運動區不同位置有密切關系。同時波束形成器的時間濾波長度參數P對想象運動任務的識別有一定的影響。另外即使沒有使用半監督分類方法，小訓練樣本集ay也得到較好識別。在使用少通道數下獲得很好的識別正確率，這說明通道數遠小于源信號數的假設是可行的。
本文方法和其它方法的識別率比較如表 2所示。第一種方法是沒有去偽跡直接采用帶通濾波器和CSP結合方法；第二種方法是在通道參數相同情況下FastICA和CSP結合方法；第三種方法是本文方法；第四種方法是競賽第一名采用CSP的方法，對于小樣本訓練數據ay，競賽第一名還使用了半監督學習方法，但半監督學習方法將要付出很大的訓練時間代價。
表2
本文方法和其它方法的識別率(%)比較
Table2.
Recognition rate (%) of different methods
通過第一組實驗仿真結果, 可以看到al數據能取得較高的正確率，正確率都在90%以上，尤其是當選取通道為C1、C2、C3、C4、C5、Cz時能100%識別出右手想象運動和右腳想象運動任務。小訓練樣本數據ay在沒有使用半監督學習方法下，我們的方法也能取得較理想的結果。我們的方法要比只使用CSP方法和FastICA+CSP方法的分類性能高。此外，本文提出的方法與競賽第一名進行了比較，對于al數據，本文的方法與競賽第一名一樣都能百分之百地予以識別；對于ay數據，由于競賽第一名針對小樣本訓練集ay使用半監督學習，識別率達到了97.6%。我們注意到競賽第一名選取的通道數都在20個以上，而我們的方法選取6個通道就能實現相當的識別率。
4.2.2 第二組實驗
本次試驗采用競賽4數據Ⅰ的四個受試者的腦電信號數據進行實驗仿真，每個實驗者進行200次試驗，把200次試驗隨機抽取為140次試驗作為訓練集，剩下的60次試驗作為測試集，一共進行十次隨機分配訓練集和測試集數據。實驗結果分別與連續小波變換(continuous wavelet transform, CWT)、同步擠壓小波變換(synchrosqueezed wavelet transform, SST)、EMD、2013年文獻[13]提出的多元經驗模式分解(multivariate empirical mode decomposition, MEMD)以及FastICA方法進行比較。文獻[13]選擇了FC3、FC4、Cz、C3、C4、C5、C6、T7、T8、CCP3和CCP4這11個通道的腦電信號作為實驗的處理信號。本文方法只選擇C3、C4、Cz、CP3、CP4、CPz六個通道的腦電信號進行實驗仿真。每個受試者的腦電數據在混合矩陣估計時分成200塊，而在使用MMSE波速形成器分離源信號時分成800塊。實驗結果如表 3所示。
表3
a、b、f和g數據采用不同方法的識別率(%)(表 3中對于受試者f和g來說，本文方法的識別率最高；對于受試者a和b來說，最高識別率為MEMD方法，居次位的是本文方法，識別率與MEMD方法差別不大。從4個受試者識別率的標準差來看本文方法穩定性最好。MEMD的11個通道多于本文方法6個通道，MEMD運行時間是本文方法的上百倍，計算復雜度高。
第二組實驗充分地證明了本文提出的方法在通道數較少時仍然取得很好的識別率，時間開銷較低，顯示了將該方法用于去除腦電信號偽跡的可行性。
5 結論
考慮到改善傳統去噪方法不能很好去除與腦電信號頻率重疊且幅度高的眼電等偽跡的問題，基于非平穩腦電信號分塊思想和SOBIUM算法，本文提出了欠定盲源分離結合CSP的特征提取新方法，首先基于欠定盲源分離提取腦電源信號，然后結合CSP特征進行分類，提高了BCI想象任務識別率，時間開銷很小。我們提出的改進之處如下：對SOBIUM算法得到的混合分離矩陣采用LM算法進行優化，提高混合矩陣盲識別的精確度；對MMSE波速形成器添加一些穩定性約束條件。新方法在電極通道數較少的情況下能夠提取腦電源信號，自動去除與想象運動任務無關的源信號。進一步研究將重構源信號和其它特征結合的特征提取方法。
