引入權重系數重構個體模板的穩態視覺誘發電位識別*

潘雋鍇，馬玉良，席旭剛，孫明旭，張建海

(1.杭州電子科技大學圣光機聯合學院，浙江 杭州 310018；2.杭州電子科技大學自動化學院，浙江 杭州 310018；3.濟南大學自動化與電氣工程學院，山東 濟南 250022；4.杭州電子科技大學計算機學院，浙江 杭州 310018；5.浙江省腦機協同智能重點實驗室，浙江 杭州 310018)

腦科學是研究腦部結構與功能的學科，其核心是以大腦為認知主體的神經科學[1-2]。腦科學的發展是生命科學中的一個重要領域，是當下最為前沿的研究領域之一。腦機接口(brain-computer interface，BCI)技術是指不依賴人體的外圍神經和肌肉組織，通過解碼大腦的意識活動而實現人腦與外部設備的通信交流[3]。目前腦機接口的發展受限于識別率低、通信速度慢等問題[4]，而基于穩態視覺誘發電位的腦機接口系統(steady-state visual evoked potential，SSVEPBCIs)由于具有較高的識別準確率、信息傳輸率，以及實驗環境配置簡單、受試者僅需進行少量的訓練等優點，近十年來備受關注[5-6]。

穩態視覺誘發電位是指當人持續注視某一以固定頻率閃爍的視覺刺激時，大腦視覺皮層會在刺激頻率或諧波頻率處產生明顯的電位變化。SSVEP的視覺刺激頻率范圍一般介于4到50 Hz之間，分為低頻段(4 Hz～15 Hz)、中頻段(15 Hz～30 Hz)和高頻段(30 Hz～50 Hz)，SSVEP響應幅值的全局最大值大約在10 Hz出現[7]，高頻段的刺激頻率能誘發的響應最小，目前大多數系統所采用的視覺刺激主要集中于中低頻段。

SSVEP-BCI系統的性能主要取決于以下三個主要因素:視覺刺激呈現、多目標編碼方式和目標識別算法[8]。與其他類型的腦機接口類似，信號處理方法、機器學習、深度學習等方法被廣泛應用到SSVEP的檢測和識別。基于功率譜密度分析(power spectral density analysis，PSDA)的方法如快速傅里葉變換(fast Fourier transform，FFT)被應用于單通道腦電信號(EEG)的頻率檢測[9]；基于最小能量組合(minimum energy combination，MEC)的頻率識別算法[10]通過尋找空間濾波器將原始EEG進行投影，得到低維組合信號，以降低噪聲和其他偽跡信號的影響；基于多變量同步指數(multivariate synchronization index，MSI)的頻率識別算法[11]通過計算兩組多通道信號之間的同步指數，篩選最大同步指數對應的頻率進行目標識別。Lin等人[12]首先將典型相關分析(canonical correlation analysis，CCA)用于SSVEP-BCI系統的頻率識別，由于其具有不需要參數優化，識別率和信息傳輸率較高的優勢，此后基于CCA的改進算法被廣泛研究。Poryzala等人[13]提出了CCA的聚類分析(cluster analysis of CCA，CACCA)；Pan等人[14]提出相位約束的CCA算法(phase-CCA)；Bin等人[15]提出了基于個體模板的CCA算法(individual template CCA，it-CCA)；Wang等人[16]提出了擴展典型相關分析(eCCA)，利用個體訓練數據增強SSVEP的檢測；Zhang等人[17]提出了多途徑CCA算法(Multiway CCA)以尋找合適的參考信號；Zhang等人[18]又提出了多集CCA(Multiset CCA)以優化來自參考信號的空間特征；Nakanishi等人[19]提出了任務相關成分分析(task related component analysis，TRCA)，通過在任務期間最大化復現性提取任務相關成分。

當需要考慮受試者校準數據時，大多數相關算法如it-CCA、eCCA、MultisetCCA等選擇直接對訓練數據平均化以得到包含特定受試者信息的模板信號，這種傳統的構造方式可能稍顯粗略，從而導致得到的空間濾波器不夠十分準確，影響算法的識別結果。因此本文針對電話撥號穩態視覺誘發電位數據集重新選擇系數特征組合，在構造個體模板過程中引入反映各訓練數據信號質量的權重系數，提出了一種新的方法即coef-eCCA，并通過實驗驗證了該方法的有效性。

1 實驗數據集與預處理

1.1 電話撥號系統數據集

本數據集由加州大學圣迭戈分校(UCSD)提供，所有實驗參與者閱讀并簽署知情同意書[16]。視覺刺激的布局如圖1所示，12個目標排列為4×3的矩陣，作為電話撥號系統的虛擬鍵盤，形狀為為6 cm×6 cm的正方形，兩兩之間水平和垂直間隔分別為5 cm和1.5 cm。實驗采用了頻率和相位聯合編碼的方式，相鄰的兩個目標同時設置不同的頻率與相位。刺激頻率設置為9.25 Hz到14.75 Hz，等差間隔0.5 Hz；相位設置為0到1.5π，等差間隔0.5π，各個目標具體對應的頻率和相位信息見圖1(b)。刺激界面程序在MATLAB下使用PsychToolbox(PTB)工具箱編寫。

共有10名視力正常或矯正至正常的健康受試者(9名男性，1名女性)參與了本次實驗，實驗在昏暗的房間里展開，受試者坐在距離顯示器60 cm的舒適椅子上。這項研究進行了模擬在線BCI實驗，記錄數據進行離線分析，腦電數據由位于覆蓋受試者腦部枕葉區域的8個電極采集。對于每位受試者，實驗由15個模塊組成，每個模塊中受試者被要求以隨機順序注視程序指示的一個視覺刺激4 s，完成對應于12個目標的12次試驗。在每次試驗開始時，一個紅色方塊會出現在目標刺激位置持續1 s，見圖1(a)，受試者需要在1 s內將目光從上一個目標轉移到本次試驗的目標上。為了減少眼動偽跡信號，受試者被要求在刺激期間避免眨眼。所有數據被降采樣到256 Hz以減少數據量。

圖1 電話撥號SSVEP-BCI刺激界面

1.2 數據預處理

有研究人員提出在人的視覺系統中存在140 ms的視覺延遲[20]，因此對于Tw s的時間窗，提取的數據段為[0.14 s，0.14＋Tw s]，時間0顯示刺激開始。之后使用零相位切比雪夫無限脈沖響應(infinite impulse response，ITR)濾波器對所有數據進行帶通濾波，通帶頻率為6 Hz到80 Hz，刺激頻率為10.25 Hz時濾波前后腦電信號時、頻域波形如圖2所示，濾波后有效地去除了腦電信號中的基線漂移和刺激頻率外頻段的干擾。由于實驗只記錄了位于腦皮層視覺區域的8個電極的腦電數據，所以無需進行通道選擇操作。

圖2 濾波前后腦電信號時、頻域波形圖

2 目標識別算法

2.1 典型相關分析

典型相關分析算法(CCA)是一種經典的多元統計方法，用于測量兩組變量之間的潛在相關性。其目標是尋找一對線性組合，使得經過變換后的兩組變量之間相關性最大[21]。考慮兩組多維變量X，Y，X∈RC×Ns為一組多通道的腦電信號，Y∈R2Nh×Ns為一組人工正-余弦參考信號，其長度與X相同，設置如下:

式中:fk為刺激頻率，fs為采樣率，Nh為諧波數量。將腦電信號與正-余弦參考信號進行典型相關分析，通過尋找一對權重向量wx∈RC，wy∈R2Nh，使得映射后的兩組典型變量x＝wTx X和y＝wTy Y相關性最大，相關性由二者之間Pearson相關系數確定:

將多通道的時域腦電信號X與包含所有刺激頻率的正-余弦參考信號組Y＝[Y1，Y2，…，Yk]做典型相關分析，選擇相關系數最大值對應的頻率作為SSVEP的目標刺激頻率:

2.2 擴展典型相關分析

典型相關分析是一種無監督學習方法，它最初僅被用來檢測頻率，但隨著越來越多的實驗范式采用聯合頻率相位的編碼方式，如何在識別過程中有效利用SSVEP的相位信息變得很重要，而CCA無法被用于區分不同的相位[8]。有學者將標準CCA與基于個體模板的CCA相結合，加入訓練數據，提出了擴展典型相關分析(eCCA)。除人工構造的正-余弦參考信號外，對典型相關分析算法進行計算結構的調整，通過將受試者的訓練數據平均化，構造了個體模板信號:

式中:Xk為訓練數據，Ntrain為訓練數據個數。由這3種信號可以得到6種空間濾波器的形式——①測試信號X與個體模板之間:wX()，w＾X()；②測試信號X與正-余弦參考Y之間:wX(XY)，wY(XY)；③個體模板信號與正-余弦參考信號Y之間:w＾X()，wY()。將這些空間濾波器映射到各種信號上可以得到多種特征。在Masaki等人的研究中，采用了如下5種相關系數的組合方式[16]:

ρ(a，b)即為兩組信號a和b的Pearson相關系數，計算方式同式(2)。采用一個集成分類器用于組合這5種特征作為最終特征用于識別:

在計算測試信號與個體模板信號之間的相關系數過程中會產生負數，因此集成分類器添加sign()項用于保留二者間的負相關信息。最后識別刺激目標同式(3)。

2.3 本文提出的方法

2.3.1 特征系數選擇

除去擴展典型相關分析中提出的5個相關系數，由6種空間濾波器與測試信號、正-余弦參考信號和模板信號分別映射，總共可以得到10個典型變量，見表1。由這10個典型變量兩兩之間進行相關系數計算，可以得到45個系數特征。

表1 10個典型變量

為了減少算法在計算過程中消耗的時間，需要降低算法的計算復雜度。各實驗數據進行系數特征選擇，在保證識別正確率的同時，應盡量減少特征數量。本文選擇了如下3個系數特征，為了與標準eCCA對比我們將其命名為pro-eCCA，圖3顯示了系數特征提取流程。

圖3 系數特征提取流程圖

2.3.2 重構個體模板

與標準典型相關分析相比，引入個體模板信號后算法的識別性能得到顯著提高，這不僅與模板信號中包含重要的相位信息有關，還因為通過將訓練數據平均化可以有效地消除噪聲影響[22]。本文在此操作的基礎上提出了一種更為細致的個體模板構造方式，即在訓練數據中引入各次試驗權重系數，得到一種新的算法coef-eCCA。

在信號處理領域中，通常采用頻譜幅值和信噪比衡量信號質量的好壞，信噪比(signal-noise ratio，SNR)即信號與背景噪聲之比。本文分別以訓練數據的頻譜幅值與信噪比作為評價指標，給予各次試驗信號權重系數。定義頻率fk處的信噪比為幅頻響應曲線中fk處的幅值與附近L個頻率的幅值均值之比[23]:

式中:amp(fk)為SSVEP在頻率fk處的頻譜幅值，L取值為16，相鄰頻率間隔Δf為0.125Hz。

之后在構造個體模板中分別采用兩種計算方式確定權重系數——第一種為計算頻譜中所有12個刺激頻率對應的評價指標(evaluation index，EI)的平均值，見式(9)；第二種為針對各刺激頻率，計算其基頻與各諧波頻率所對應的評價指標，見式(10):

式中:Nh為諧波個數，本文設置為5，即針對頻率為9.25Hz時計算9.25 Hz、18.5 Hz、27.75 Hz、37 Hz、46.25 Hz對應的評價指標，其他頻率以此類推。

最后對各次試驗數據的評價指標做歸一化處理，得到構造個體模板中各試次所對應的權重系數:

∑EI為各次試驗數據評價指標之和，權重系數wtem∈RC×Ntrain，C為通道數，Ntrain為訓練數據個數。重構后的個體模板信號為:

3 實驗結果與分析

交叉驗證(Cross Validation)是在建立模型和驗證模型參數時常用的方法[24]，將樣本數據切分組合成不同的訓練集和測試集，用訓練集訓練模型，用測試集評估模型預測的好壞，交叉驗證適用于樣本數據不充足的情況下。本文采用留一法交叉驗證(Leave-one-out Cross Validation)，即將某一識別目標15次試驗中的14項作為訓練數據用于構造個體模板，剩下的1項作為測試數據用于驗證，循環進行直至每項數據都被當做一次測試數據，以避免實驗過程中的隨機因素影響識別結果。

一般將識別準確率(accuracy)和信息傳輸率(information transmission rate，ITR)共同作為評估SSVEP-BCI的性能指標。信息傳輸率的計算公式為[25]:

式中:N為識別的目標個數，P為識別準確率，T為單次實驗目標選擇時間，通常為固定值(T＝Tw＋ts，ts＝0.5 s為兩次選擇之間的目光轉移時間)，ITR的單位為bits/min。

3.1 重新選擇系數特征

為了顯示重新選擇系數特征后的擴展典型相關分析(pro-eCCA)的識別性能，將其與典型相關分析(CCA)、標準擴展典型相關分析(eCCA)和任務相關成分分析(TRCA)作對比。Nakanishi等人認為一組多通道腦電信號是由任務相關成分，即誘發信號和任務無關成分，即大腦自發信號，共同線性組合而成[19]。通過TRCA算法對原始EEG進行線性加權，從中得到僅包含任務相關成分的信號，提高信號的信噪比。我們選擇了從0.5 s到2.5 s，5個不同的時間窗用于識別，圖4顯示了這4種方法的識別準確率和信息傳輸率。

圖4 四種方法在不同時間窗下的識別性能對比

傳統的標準CCA算法相對來說識別效果最差，在短時間窗下(Tw＝0.5 s，1.0 s)識別準確率僅為19.33%和53.94%；但隨著時間窗長度的增加，算法的識別準確率也會達到80%(Tw＝2.0 s，2.5 s)。針對本數據集重新選擇系數特征組合(pro-eCCA)后，識別準確率較標準eCCA有了一定的提高，特別是在短時間窗(Tw＝0.5 s，1.0 s)的情況下，各提高了1.67%。TRCA作為目前SSVEP識別方面最流行的算法之一，識別正確率在短時間窗下有著較大的優勢，時間窗為0.5 s時領先pro-eCCA5.78%；隨著時間窗長度的增加其優勢不斷減小，在2.0 s時TRCA為98.78%，pro-eCCA為98.67%幾乎可以與TRCA持平。信息傳輸率方面eCCA和TRCA在0.5 s時達到最大值，此后隨著時間窗長度的增加不斷降低；而CCA在短時間窗下的識別率過低，隨著時間窗長度增加，隨之提高的識別準確率依然保證了信息傳輸率的提高，并在Tw＝1.5 s時達到最大值。

為了比較各算法之間的計算復雜度，我們計算了Tw＝2.0 s下每位受試者識別過程中消耗時間的總和，見表2。重新選擇系數特征后的eCCA，由于減少了系數特征的個數從而大大降低了計算成本，每位受試者平均消耗時間相比之前減少了41.1%。TRCA算法的計算時間普遍超過eCCA，這可能是由于TRCA在求解空間濾波器時需要計算以使各試次訓練數據之間協方差最大化，而本文數據集有15個訓練模塊，因此消耗了更多的計算時間。pro-eCCA的計算過程相對來說更為簡便，因此在保證了識別準確率的同時，消耗的時間更少。

3.2 重構個體模板

為了說明重構個體模板后算法的識別性能，分別用兩種計算方式和兩個評估指標，在固定時間窗(Tw＝2.0 s)下進行實驗，每位受試者的具體結果見表3，表4。

表3 重構個體模板后每位受試者的識別準確率

表4 重構個體模板后每位受試者的計算消耗時間

從四種方式與重新選擇系數特征組合后的pro-eCCA的識別準確率對比上看，計算12個刺激頻率的頻譜幅值用于確定各試次訓練數據權重系數的方式取得了最好的結果，共有4位受試者得到提升(S1，S2，S3，S10)，其中S2的提升幅度最大(1.67%)，而S3則提高至100%，但S7的結果有所下降。計算所有刺激頻率信噪比的方式使S2獲得了提高，S7的結果與用頻譜幅值一樣也下降了，在其他受試者沒有變化的情況下平均識別準確率有略微提高。相比較而言，計算刺激頻率的基頻及諧波頻率的頻譜幅值的效果最差，雖然將S3提升至100%，但是S1，S2，S10的結果均有下降，平均準確率也降低至98.33%。最后一種使用刺激基頻和諧波頻率信噪比的方式提升了S2和S3的準確率，并且沒有受試者的結果下降，因此平均識別準確率也有一定的提高，與第一種方式一樣均超過了同時間窗下TRCA(98.78%)的結果。為驗證實驗結果是否具有統計學意義，我們進行了假設檢驗，比較第一種方式和pro-eCCA的識別結果，除去二者中識別率始終保持100%的受試，剩下結果的p＜0.05，即表明實驗結果存在顯著性差異。

從計算消耗的時間上看，使用頻率幅值作為評估指標與pro-eCCA相比有了些許下降，這是由于計算復雜度的增加不可避免地提高了計算成本，但都依然保持在一個良好的水平，不會影響實際應用。而使用信噪比作為評估指標則大幅度提高了計算成本，嚴重消耗了計算時間，其原因為計算信噪比是在計算頻譜幅值的基礎上進行。兩種計算方式對比來看，第一種好于第二種，這是由于針對各刺激頻率的訓練數據，第一種僅需計算所有12個實驗刺激頻率的評估指標；而第二種則需根據不同的刺激頻率選擇相應的諧波頻率進行計算，計算結構上更為復雜，所以計算成本更高。

為了可視化重構個體模板帶來的效果，圖5顯示了重構前后兩種模板信號的頻譜波形圖，圖5(b)為計算所有刺激頻率對應的頻譜幅值作為訓練數據權重系數的構造方式所得到的模板。對比二者的頻譜圖，重構后基頻和諧波頻率外的成分被十分有效地抑制，增強了信號的SSVEP響應，可以更為明顯地顯現出高頻諧波成分。平均化構造模板后，除去6次諧波，基波和其他諧波的信噪比都低于40 dB；重新構造模板后從基波到諧波的信噪比均高于前者。

圖5 重構前后個體模板的頻譜波形圖

3.3 討論

為進一步對比算法之間的性能，我們研究了在特定時間窗下(Tw＝2.0 s)，分別改變電極數目和實驗試次數量對算法性能的變化情況，實驗結果見圖6。

圖6 不同電極數目和訓練試次下算法的識別準確率

在研究電極數目對識別結果影響的實驗中，電極取值從3到8，每種條件下電極組合情況依次為:[PO3，POz，PO5](3個電極)，[PO3，POz，PO5，PO4](4個電極)，[PO3，POz，PO5，PO4，PO6](5個電極)，[PO3，POz，PO5，PO4，PO6，O1](6個電極)，[PO3，POz，PO5，PO4，PO6，O1，O2](7個電極)，[PO3，POz，PO5，PO4，PO6，O1，O2，Oz](8個電極)，具體電極分布位置如圖7所示。從圖6(a)中我們可以看到，提高電極數量，三種算法的性能都得到了提升，并且由于O1，O2，Oz三個電極位置更靠近腦皮層枕葉區域中央，因此增加的這幾個電極所采集的數據對于算法性能的提升更顯著(電極數由5變為6時，coef-eCCA識別率提升了4.11%)。

圖7 電極位置分布

研究試次數對識別影響的實驗中，我們分別選取了3、6、9、12和15次(即訓練數據試次為2、5、8、11、14)，結果見圖6(b)。提高實驗試次數量均會提升coef-eCCA和TRCA算法的性能，而在CCA算法中由于各個試次是完全獨立的，它們僅用于交叉驗證以提供更加可靠的結果，提高試次數沒有讓CCA算法性能得到提升。相比于eCCA，TRCA算法本身更加依賴試次數量，需要從盡可能多的訓練數據中獲取有效信息，因此在試次數較少時算法性能明顯弱于coef-eCCA。

此外我們還進行了實驗，以研究刺激頻率范圍對算法識別性能的影響。可能由于數據集設置的刺激頻率范圍較小，從最小頻率9.25 Hz到最大頻率14.75 Hz中間只間隔了5.5 Hz。將其分為三個頻段范圍(9.25 Hz～10.75 Hz，11.25 Hz～12.75 Hz，13.25 Hz～14.75 Hz)，從結果上看并無明顯差異，因此沒有在論文中展示出來，在此作文字說明。

4 總結

本文在擴展典型相關分析算法傳統的個體模板構造方式基礎上，提出了一種更為精細的構造方法。首先針對電話撥號系統的SSVEP數據集重新選擇特征系數組合，從識別結果上看不同時間窗下的識別準確率均獲得提高，并且在計算時間方面較其他對比方法取得了顯著的優勢。之后提出了一種通過評估指標在構造個體模板過程中引入權重系數的方法coefeCCA，選擇兩種計算方式和兩種評估指標分別進行實驗，在沒有大幅度影響計算成本的情況下將識別準確率最高提升至99%，驗證了改進方法的有效性。未來的工作中針對不同的數據集，還可以嘗試不同的計算方式和其他信號評估指標以提高實際識別效果。