利用高頻刺激進行編碼能夠緩解基于穩態視覺誘發電位（SSVEP）的腦-機接口（BCI）產生的用戶視覺疲勞，提升系統的舒適度和安全性，具有廣闊的應用前景。然而，當前先進的SSVEP解碼算法大多在低頻數據集上進行對比驗證，在高頻SSVEP信號上的識別性能仍然未知。針對此問題，本文采集了20名受試者在高頻SSVEP范式下的腦電（EEG）數據，對目前主流的2種典型相關分析算法、3種集成任務相關成分分析算法和1種任務判別成分分析算法展開對比。結果表明，它們均能有效解碼高頻SSVEP信號，且在不同條件下算法的分類性能指標和速度存在差異。本研究為高頻SSVEP-BCI系統的算法選擇提供了依據，在構建舒適友好型BCI系統方面具有潛在的應用價值。
引用本文： 羅睿心, 豆心怡, 肖曉琳, 吳喬逸, 許敏鵬, 明東. 腦機接口中高頻穩態視覺誘發電位的識別研究. 生物醫學工程學雜志, 2023, 40(4): 683-691. doi: 10.7507/1001-5515.202302034 復制
0 引言
腦機接口（brain-computer interface，BCI）是將中樞神經系統活動轉化為人工輸出的系統，能夠替代、修復、增強、補充或改善中樞神經系統正常輸出，從而改變中樞神經系統與內外環境間的交互作用[1-2]。經過幾十年的技術迭代，BCI技術在康復醫療、軍事、航天等領域顯示出廣闊的應用前景。在無創BCI編碼范式中，穩態視覺誘發電位（steady-state visual evoked potential，SSVEP）是指大腦在高于6 Hz的穩定周期性重復呈現的視覺刺激下，所誘發的包含刺激基波與若干諧波頻率成分的周期性腦電（electroencephalogram，EEG）信號響應[3]。SSVEP-BCI系統因其具有信噪比（signal to noise ratio，SNR）高、用戶所需訓練時間少及信息傳輸速率（information transfer rate，ITR）高等優點受到了廣泛關注[4]。
以往研究表明，重復性視覺刺激的頻率會影響SSVEP信號的強度，依據誘發信號強度的差異，閃爍刺激頻率可以被劃分為：低頻段（<15 Hz）、中頻段（15~30 Hz）和高頻段（>30 Hz）[3]。低、中頻刺激能夠誘發出高幅值、高SNR的SSVEP信號，被廣泛應用于編碼SSVEP-BCI系統指令，當前ITR已突破353.30 bits/min[5]。然而，低、中頻刺激的舒適性和安全性低，強烈的閃爍刺激容易引起患者產生視覺疲勞，甚至有引發光敏性癲癇的風險[6-7]。使用高頻刺激進行編碼能有效解決以上問題，因此相關研究隨之涌現。Wang等[3]已發現雖然高頻范圍的SSVEP信號幅值低，但由于其背景噪聲也變弱，因此具有能夠和低頻SSVEP相當的SNR。Won等[8]對比了低頻刺激和高頻刺激下的用戶疲勞指數，結果表明高頻SSVEP可以有效減輕閃爍刺激引起的視覺疲勞。自此，國內外研究人員逐漸開始利用高頻刺激對SSVEP-BCI系統進行編碼。Diez等[9]使用高頻SSVEP信號操控機器人輪椅，證實舒適的高頻刺激能夠獲得更穩定的長時程控制性能。Chen等[10]開發了刺激頻率為30~33 Hz的高頻SSVEP-BCI系統，結合計算機視覺實現了大腦對機械臂的多自由度控制。Ye等[11]提出了一種多符號時分編碼（multi-symbol time division coding，MSTDC）方法，利用高頻閃爍刺激編碼了40個腦控字符，系統平均在線ITR達到119.05 bits/min。上述研究證實了利用高頻刺激構建SSVEP-BCI系統的優勢和可行性。然而，高頻SSVEP信號幅值低，誘發特征微弱，導致其解碼效率仍然有限，阻礙了高頻SSVEP-BCI系統走向實際應用。因此，有必要開展面向高頻SSVEP信號的解碼算法研究，從而提高高頻SSVEP-BCI系統的識別性能。
近年來，SSVEP信號的解碼算法得到了快速發展。早期的SSVEP解碼多采用單導聯下的功率譜密度分析（power spectral-density analysis，PSDA）來識別信號頻率[12]。隨著典型相關分析（canonical correlation analysis，CCA）的提出，利用多導聯腦電數據進行空間濾波降維逐漸成為SSVEP識別的主流方式[13]。Chen等[14]采用了擴展CCA（extended CCA，eCCA）獲得具有特異性的空間濾波器和平均模板，最高ITR突破319.20 bits/min。任務相關成分分析（task-related component analysis，TRCA）進一步將SSVEP-BCI系統的ITR提升到了325.33 bits/min[15]。隨后，TRCA的拓展版本利用了額外信息來增強識別性能，例如鄰近頻率信息[16]、正余弦參考信號[17]、相位約束[18]等。最近，Liu等[19]提出任務判別成分分析（task-discriminant component analysis，TDCA）方法，被證實是當前最佳的SSVEP空間濾波算法之一。然而，這些先進算法大多在低頻SSVEP數據集上進行對比驗證，在高頻SSVEP信號上的識別性能仍然未知。
針對此問題，本文研究設計了刺激頻率范圍為30~39 Hz的十指令高頻SSVEP范式，離線采集了20名受試者的腦電數據，探究了6種先進SSVEP解碼算法在高頻SSVEP信號上的性能差異，為高頻SSVEP-BCI系統的算法選擇提供依據。
1 材料與方法
1.1 試驗對象
本文共招募了20名受試者參與離線腦電數據的采集，其中男性11名，女性9名，年齡在19~24周歲之間。所有受試者均視力正常或矯正視力正常，無任何精神疾病史，并自愿參與本試驗。在試驗正式開始前，所有受試者閱讀了試驗流程與注意事項，簽署了相關的知情同意書。本試驗通過了天津大學研究倫理委員會的批準（批文編號：TJUE-2021-062），所采集數據用于本文研究分析。
1.2 試驗范式
為探究高頻刺激下SSVEP解碼算法的識別性能，設計了刺激頻率范圍為30~39 Hz（步長為1 Hz）的高頻SSVEP范式，共包含10個刺激頻率。每個刺激為一個120×120像素的方塊，兩個相鄰刺激之間的垂直距離和水平距離均為350像素。刺激程序由商業數學軟件MATLAB R2019a（MathWorks Inc.，美國）中的Psychphysic toolbox v3.0.15工具包編寫完成，呈現在59.68 cm×33.57 cm的液晶顯示屏上，其分辨率為1 920×1 080，刷新率為120 Hz。
試驗共10個刺激指令，每個刺激以偽隨機序列依次閃爍1個試次為1輪，每5輪為1組試驗，每組試驗結束后讓受試者進行適當休息，防止視覺疲勞對試驗結果產生影響。每名受試者各進行了4組試驗，因此每個指令下均包含20個試次腦電數據。本文解碼過程為十分類任務，需對十個不同頻率下誘發的腦電信號進行識別。
正式試驗時，受試者端坐于顯示屏前，雙眼距離屏幕70 cm。調整為平靜狀態后，點擊“空格”鍵開始試驗。每個試次下刺激呈現持續時間為2 s，相鄰刺激之間的空隙時間（即刺激提示時間）為2 s，如圖1所示。試驗過程中，受試者在刺激呈現時間內盡量避免眨眼，并保持注意力集中。
圖1
試驗流程圖
Figure1.
Experimental flow chart
1.3 數據采集
本研究試驗全程在電磁屏蔽室中完成。腦電放大器（SynAmps, Neuroscan Inc., 美國）和配套采集軟件記錄64導腦電信號，導聯放置遵循國際10-20標準，采樣率為1 000 Hz，系統濾波頻帶設置為直流（direct current，DC）到200 Hz之間，以50 Hz陷波器去除工頻干擾，受試者的頭皮阻抗保持在10 kΩ以下。試驗過程中，刺激程序通過并口通訊向腦電信號放大器發送事件觸發信號，并進行儲存。
2 解碼算法
2.1 典型相關分析及其拓展算法
CCA是一種用于衡量多維隨機變量的潛在相關性的統計學方法，其本質是求解使得兩組隨機變量降維后相關性最大的最優線性組合。Bin等[13]將CCA應用于多導聯SSVEP信號的分類識別中，其性能顯著優于傳統的PSDA。作為一種無監督方法，CCA算法因其計算簡單、使用便捷而被廣泛應用于SSVEP-BCI領域。Chen等[14]引入個體訓練數據來構建CCA濾波器和個體平均模板，即eCCA。與傳統CCA算法相比，這種有監督方法考慮了腦電信號的跨個體、跨時間變異性，顯著提升了SSVEP-BCI系統的識別效率。在此基礎上，Wong等[16]設計了一種跨刺激學習框架，從目標刺激及其若干鄰近刺激中學習空間濾波器，提出了多刺激eCCA（multi-stimulus eCCA，ms-eCCA）算法，能夠進一步增強eCCA空間濾波器的魯棒性。跨刺激學習框架有效挖掘了不同頻率間SSVEP信號的共同信息，改善了eCCA算法在訓練樣本不足時的分類性能。
2.2 集成任務相關成分分析及其拓展算法
腦電信號可看作是不同任務相關源信號的線性組合，TRCA的目標是通過最大化相同任務下信號的相似性來提取任務相關成分[15]。由于相鄰刺激頻率下誘發的SSVEP信號具有相似的頭皮空間分布，因此對不同類別下空間濾波器進行集成可以增強算法的魯棒性，即集成TRCA（ensemble TRCA，eTRCA）算法。近年來，衍生出了多種eTRCA的擴展方法。例如，多刺激eTRCA（multi-stimulus eTRCA，ms-eTRCA）是跨刺激學習框架應用在eTRCA上的一種先進算法[16]。Wong等[17]首次將頻率信息引入到TRCA空間濾波器的構建中，開發了帶有正余弦參考信號的eTRCA（eTRCA with sine-cosine reference signal，eTRCA-R），其本質是將腦電信號投影到正余弦子空間后再求解空間濾波器。
2.3 任務判別成分分析
TDCA算法解決了傳統eTRCA中空間濾波器冗余的問題，是當前低頻SSVEP信號識別性能最佳的算法之一[19]。與eTRCA不同的是，TDCA對訓練數據進行空間維度和時間維度下的兩次增廣后，利用二維線性判別分析為不同類別信號求解一個共享的空間濾波器。
假設同一類別下的多通道腦電信號表示如式（1）所示：
 |
其中，Nc為導聯數量、Np為采樣點數量，Nt為訓練試次數量，i為訓練試次的索引，
代表實數集。首先，對所有訓練數據進行時延，得到提高空間維度后的一階增廣信號
，如式（2）所示：
 |
其中，
表示腦電數據延后
個點后的腦電數據，T為轉置符號。其次，將
投影至對應類別頻率下參考信號的正交子空間，并將投影前后的信號按照時間維度進行拼接后得到二階增廣信號，如式（3）所示：
 |
其中，A為二階增廣后的腦電信號，
為正交投影矩陣[19]。
TDCA采用二維線性判別模型求解空間濾波器，計算所有類別下訓練數據的類間散度矩陣Sb和類內散度矩陣Sw，如式（4）、式（5）所示：
 |
 |
式中，Nf為訓練類別數量，
代表第n個類別、第i個試次下的二階增廣信號，
代表第n個類別下的模板信號，
代表所有類別中心。
化簡后，TDCA優化目標如式（6）所示，W為需要求解的空間濾波器，
為W的估計值。通過最優化問題推導可得，取
的前Nk個特征向量即可構建TDCA空間濾波器
。
 |
將單次測試信號進行補零后，采用與訓練過程相同的步驟獲得二階增廣后的測試信號
，計算
與
經過空間濾波后的相關系數
，如式（7）所示，其最大值所對應的類別為識別結果。
 |
2.4 性能評估
在預處理階段，首先將原始腦電信號降采樣至250 Hz，并提取枕區周圍典型的9個導聯（Pz、PO5、PO3、POz、PO4、PO6、O1、Oz、O2）數據用于后續分類識別[15, 19]。通過事件標志對原始信號進行截取得到[
~（
）]時刻的數據。其中，
表示信號潛伏期時長，Tw是用于分類的時間窗長度（Tw = 0.2, 0.4, ···, 1 s）。
本文選擇了當前在低頻SSVEP-BCI中被廣泛用于對比分析的6種先進SSVEP解碼算法[19]，包括eCCA、ms-eCCA、eTRCA、ms-eTRCA、eTRCA-R和TDCA。所有算法均應用濾波器組（filter bank，FB）分析，子帶數量Nb設置為2，第m個子帶的下限和上限截止頻率分別為m×30 Hz和90 Hz，不同子帶融合的權重系數如式（8）所示[10, 20]：
 |
其中，a(m)代表第m個子帶的權重值。
6種算法均采用模板匹配的方式來預測類別[15, 19]。具體來說，在計算出空間濾波器后，將測試腦電信號分別與每個類別下的平均信號（或參考信號）進行空間濾波降維后，求解相關系數作為決策值。在決策值組成的特征向量中尋找最大值對應的類別作為預測結果。采用隨機抽樣交叉驗證方法進行分類驗證，隨機抽取Nt 個輪次數據作為訓練集（Nt = 2, 4, …, 10），重復10次取平均值。計算分類準確率、曲線下面積（area under the curve，AUC）和ITR作為評估指標。ITR（以符號ITR表示）計算方式如式（9）所示，單位為bits/min。
 |
其中，M表示指令類別數量，P為分類準確率，T為輸出單個指令所需的時間（單位：s），包括注視時間（Tw）和視線轉移時間，本研究將視線轉移時間設定為0.5 s。
研究使用統計產品與服務解決方案軟件（statistical product and service solutions，SPSS）（IBM Inc.，美國）進行統計學分析。通過三因素重復測量方差分析（three way repeated measures analysis of variance，Three-way RMANOVA）探究算法因素、時間窗長度因素和訓練集樣本數目因素對分類準確率、AUC和ITR的影響，并采用單因素重復測量方差分析（one way repeated measures analysis of variance，One-way RMANOVA）探究不同算法在固定時間窗或者訓練樣本條件下時性能差異是否具有統計學意義。如果數據不符合球形檢驗（Mauchly’s test of sphericity），則進行格林豪斯-蓋舍校正（Greenhouse-Geisser correction）。
3 結果分析
3.1 潛伏期估計評估
SSVEP信號的潛伏期是指視覺通路傳導引起的SSVEP響應與刺激相位之間的時滯，可通過擬合相位與刺激頻率的線性函數關系進行估計[21]。如圖2所示，顯示了Oz導聯下所有受試者的SSVEP平均相位隨刺激頻率的線性變化。通過進一步計算線性模型斜率得到的平均潛伏期為（132.99 ±5.61）ms。基于此，本文分析，在SSVEP信號的提取中應去除前0.13 s的潛伏期信號（
= 0.13 s）。
圖2
平均相位的變化趨勢
Figure2.
Variety of average phase
3.2 參數優化
采用網格搜索法對高頻SSVEP信號識別中不同解碼算法適用的參數進行優化，需要優化的參數包括ms-eCCA和ms-eTRCA中的鄰近頻率數量d，TDCA中的子空間數量k和延遲點數l[16, 19]。如圖3所示的是三種算法在不同參數條件下的歸一化準確率，ms-eCCA和ms-eTRCA選擇沿熱點圖縱坐標軸方向平均后最大值對應的d作為最優參數，TDCA選擇沿熱點圖橫、縱坐標軸方向平均后最大值對應的l和k作為最優參數。基于此，ms-eCCA和ms-eTRCA中的鄰近頻率數量
分別設置為5和3，TDCA中的子空間數量
和延遲點數
分別設置為3和1。
圖3
三種分類算法的歸一化準確率
Figure3.
Normalized accuracy of the three classification algorithms
3.3 分類性能對比
為探究6種算法在高頻SSVEP數據上的分類性能，計算了在不同時間窗長度（Tw = 0.2, 0.4, …, 1.0 s）和不同訓練集樣本數量（Nt = 2, 4, …, 10）下的分類準確率、AUC及ITR指標，并對20名受試者的結果取平均值，如圖4所示。Three-way RMANOVA的結果表明，不同算法、不同時間窗長度和不同訓練集樣本數量均會影響算法識別性能，其差異具有統計學意義（P < 0.01, 具體效應值詳見附表1）；且算法與時間窗長度之間、算法與訓練集樣本數量之間的交互作用具有統計學意義（P < 0.01, 具體效應值詳見附表1）。
圖4
不同分類算法之間的表現比較
Figure4.
Performance comparison between different target recognition methods
對所有訓練集樣本數量下的結果進行平均，得到分類性能隨時間窗長度增加時的變化趨勢，如圖4左側子圖所示。One-way RMANOVA的結果表明，所有時間窗長度條件下6種算法的識別性能差異均具有統計學意義（AUC指標：Tw = 0.6 s時，F(1.5，20.3) = 5.7，P = 0.01；Tw = 0.8 s時，F (1.4，26.3) = 5.1，P = 0.02；Tw = 1.0 s時，F (1.4，26.1) = 4.7，P = 0.03；其余條件P < 0.01, 具體效應值詳見附表2）。6種算法的平均分類準確率和AUC隨時間窗長度增大而上升，平均ITR呈現出先上升后下降的趨勢，并在Tw = 0.4 s時達到峰值。具體來看，在所有時間窗長度條件下，eCCA算法的分類準確率、AUC和ITR最低；其次是eTRCA、ms-eTRCA以及eTRCA-R，其中ms-eTRCA以及eTRCA-R的識別性能均優于原始eTRCA；TDCA算法的分類性能最佳，在Tw = 0.4 s時平均ITR達到93.40 bits/min。在時間窗長度較小時，ms-eCCA算法分類結果低于eTRCA、ms-eTRCA、eTRCA-R和TDCA；在時間窗長度較大（Tw ≥ 0.6 s）時，分類性能超過eTRCA、ms-eTRCA和eTRCA-R，并且能夠達到與TDCA算法相當的水平。從該結果中可以看出，在6種算法中，TDCA算法在所有時間窗條件下綜合表現更優，而ms-eCCA算法對時間窗長度變化更加敏感。
對所有時間窗長度下的結果進行平均，得到分類性能隨訓練集樣本數量增大時的變化趨勢，如圖4右側子圖所示。One-way RMANOVA的結果表明，所有訓練集樣本數量條件下6種算法的識別性能差異均具有統計學意義（P<0.01, 具體效應值詳見附表3）。總體來看，6種算法的平均分類準確率、AUC和ITR隨訓練集樣本數量增大而上升。具體而言，在所有訓練集樣本數量條件下，eCCA算法的分類結果最低；在訓練集樣本數量不足（Nt ≤ 4）時，ms-eCCA的分類性能依次優于TDCA、eTRCA-R、ms-eTRCA和eTRCA；在訓練集樣本數量足夠（Nt ≥ 8）時，TDCA算法的分類準確率和ITR最高，其次是eTRCA-R、ms-eTRCA以及eTRCA，而ms-eCCA算法僅優于eCCA。該結果表明，相較于其他5種算法而言，ms-eCCA算法對訓練集樣本數量不敏感，在小樣本情況下具有明顯的性能優勢；而在訓練集樣本數量足夠時，TDCA算法能夠達到更優的分類效果。
3.4 算法速度對比
在實際應用中，算法的建模速度和解碼速度同樣是影響BCI性能的關鍵因素。其中，建模速度由構建模型的運算過程決定，會對系統的準備時間造成影響，而單次解碼速度能夠體現在線BCI系統的實時效率。采用具有3.20 GHz中央處理器（central processing unit，CPU）的計算機測試了6種算法在不同訓練集樣本數量下的平均建模時間，以及在不同時間窗長度下的平均單次解碼時間，如圖5所示。
圖5
不同分類算法在不同訓練集樣本數量下的建模時間和不同時間窗長度下的解碼時間
Figure5.
Training time of different recognition methods at different training sample and decoding time at different time window lengths
在建模速度方面，隨著訓練集樣本數量的增加，所有算法的建模時間均呈現上升趨勢。其中，eCCA的訓練過程僅需要獲得個體模板，因此建模速度最快，而ms-eCCA與TDCA算法的建模時間明顯高于其余算法。由于ms-eTRCA和eTRCA-R將算法核心步驟中協方差求解過程進行優化，避免了無用且耗時的循環過程，因此建模時間也明顯低于原始eTRCA算法。
在解碼速度方面，隨著時間窗長度增加，所有算法的解碼時間逐漸增加。其中eTRCA、ms-eTRCA以及eTRCA-R三種算法解碼過程一致，因此解碼速度相同。在所有時間窗長度條件下，這3種算法計算開銷最少，而eCCA解碼速度最慢。當時間窗長度較小（Tw = 0.2 s）時，TDCA的單次解碼時間比ms-eCCA短；而當時間窗長度逐漸增大（Tw ≥ 0.4 s）后，ms-eCCA的解碼速度優于TDCA。
4 討論
本文離線采集了20名受試者的高頻十分類SSVEP數據，探究了包括eCCA、ms-eCCA，eTRCA、ms-eTRCA、eTRCA-R和TDCA在內的6種先進SSVEP解碼算法在不同訓練集樣本數量、不同時間窗長度下對高頻SSVEP信號的解碼性能差異，并比較了不同算法的建模和解碼速度。
本文研究首先在高頻SSVEP信號上優化了先進SSVEP解碼算法的相關參數，為高頻SSVEP-BCI研究中先進算法的參數選擇提供了參考。其次，通過計算分類準確率、AUC和ITR三個指標來評估算法的分類性能。結果表明，6種算法的分類性能排序可以表示為：TDCA≈ms-eCCA > eTRCA-R > ms-eTRCA > eTRCA > eCCA。其中，與其他算法相比，ms-eCCA算法對時間窗長度更敏感，但對訓練集樣本數量不敏感。因此，在訓練集樣本數量不足時，ms-eCCA算法優于研究者最新提出的TDCA算法。然而，雖然TDCA和ms-eCCA具有更高的分類性能，但由于引入了更復雜的求解過程，因此建模速度和解碼速度明顯低于eTRCA、ms-eTRCA和eTRCA-R算法。在實際應用中，應結合實際情況選擇適合的分類算法構建高頻SSVEP-BCI系統。
表1總結了近年來已有研究中報道的高頻SSVEP-BCI系統，從頻率范圍、指令數量、識別算法、ITR四個方面進行比較。可以看出，現有高頻SSVEP-BCI系統研究大多集中在30~40 Hz頻段，且超過半數使用(e)TRCA進行分類識別[7-8, 10-11, 22-29]。本文研究結果表明，ms-eCCA、ms-eTRCA、eTRCA-R和TDCA在高頻SSVEP信號上同樣具有與eTRCA相當或更優的識別性能，有望在后續被應用于高頻SSVEP-BCI系統中。具體而言，考慮分類性能的情況下，在訓練集樣本數量不足時優先選擇能夠抵御小樣本問題的ms-eCCA算法，在訓練集樣本數量充足時優先選擇能夠獲得最優識別結果的TDCA算法；考慮解碼速度的情況下，在保證分類性能的基礎上推薦eTRCA-R算法。
表1
已有高頻SSVEP-BCI系統比較
Table1.
Comparison among high-frequency SSVEP-BCI systems
與具有強閃爍感的低頻閃爍刺激相比，高頻刺激具有安全、舒適的優勢，因此有助于開發面向實際應用場景的SSVEP-BCI系統[8]。然而，高頻刺激下誘發的SSVEP信號微弱，識別更加困難。本文研究結果可為高頻SSVEP-BCI系統的算法選擇以及相應的參數選擇提供參考，從而提升高頻SSVEP信號的識別效果。除了解碼方法外，未來還可在編碼范式上對高頻SSVEP-BCI進行優化。例如，融合頻分多址[10]、時分多址[11]、空分多址[30]、碼分多址[31]等多維度編碼特征來提升指令間的可分性。進一步地，將高頻頻率拓展到高于人眼閃爍融合頻率的范圍，可構建舒適無感SSVEP-BCI系統[26, 28]。
5 結論
為了衡量先進SSVEP解碼算法在高頻SSVEP信號上的識別性能，本文研究采集了20名受試者的高頻SSVEP數據，以分類準確率、AUC值、ITR、算法速度作為評價指標在不同訓練集樣本數量、不同時間窗長度條件下展開對比。在6種算法中，TDCA和ms-eCCA算法表現出最佳的分類性能，其中ms-eCCA算法在訓練集樣本數量不足時具有明顯優勢。eTRCA-R算法分類性能僅次于前兩者，但解碼速度更快，能夠保證在線系統的實時效率。本文結論為高頻SSVEP研究的算法選擇提供了理論依據，在構建舒適友好型SSVEP-BCI系統方面具有潛在的應用價值。
重要聲明
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
作者貢獻聲明：羅睿心負責試驗設計、程序編寫以及論文撰寫；豆心怡負責數據分析、數據整理以及論文撰寫；肖曉琳、吳喬逸、許敏鵬負責論文指導與修改；明東參與論文的審閱。
倫理聲明：本研究通過了天津大學研究倫理委員會的審批（批文編號：TJUE-2021-062）。
本文附件見本刊網站的電子版本（biomedeng.cn）。
0 引言
腦機接口（brain-computer interface，BCI）是將中樞神經系統活動轉化為人工輸出的系統，能夠替代、修復、增強、補充或改善中樞神經系統正常輸出，從而改變中樞神經系統與內外環境間的交互作用[1-2]。經過幾十年的技術迭代，BCI技術在康復醫療、軍事、航天等領域顯示出廣闊的應用前景。在無創BCI編碼范式中，穩態視覺誘發電位（steady-state visual evoked potential，SSVEP）是指大腦在高于6 Hz的穩定周期性重復呈現的視覺刺激下，所誘發的包含刺激基波與若干諧波頻率成分的周期性腦電（electroencephalogram，EEG）信號響應[3]。SSVEP-BCI系統因其具有信噪比（signal to noise ratio，SNR）高、用戶所需訓練時間少及信息傳輸速率（information transfer rate，ITR）高等優點受到了廣泛關注[4]。
以往研究表明，重復性視覺刺激的頻率會影響SSVEP信號的強度，依據誘發信號強度的差異，閃爍刺激頻率可以被劃分為：低頻段（<15 Hz）、中頻段（15~30 Hz）和高頻段（>30 Hz）[3]。低、中頻刺激能夠誘發出高幅值、高SNR的SSVEP信號，被廣泛應用于編碼SSVEP-BCI系統指令，當前ITR已突破353.30 bits/min[5]。然而，低、中頻刺激的舒適性和安全性低，強烈的閃爍刺激容易引起患者產生視覺疲勞，甚至有引發光敏性癲癇的風險[6-7]。使用高頻刺激進行編碼能有效解決以上問題，因此相關研究隨之涌現。Wang等[3]已發現雖然高頻范圍的SSVEP信號幅值低，但由于其背景噪聲也變弱，因此具有能夠和低頻SSVEP相當的SNR。Won等[8]對比了低頻刺激和高頻刺激下的用戶疲勞指數，結果表明高頻SSVEP可以有效減輕閃爍刺激引起的視覺疲勞。自此，國內外研究人員逐漸開始利用高頻刺激對SSVEP-BCI系統進行編碼。Diez等[9]使用高頻SSVEP信號操控機器人輪椅，證實舒適的高頻刺激能夠獲得更穩定的長時程控制性能。Chen等[10]開發了刺激頻率為30~33 Hz的高頻SSVEP-BCI系統，結合計算機視覺實現了大腦對機械臂的多自由度控制。Ye等[11]提出了一種多符號時分編碼（multi-symbol time division coding，MSTDC）方法，利用高頻閃爍刺激編碼了40個腦控字符，系統平均在線ITR達到119.05 bits/min。上述研究證實了利用高頻刺激構建SSVEP-BCI系統的優勢和可行性。然而，高頻SSVEP信號幅值低，誘發特征微弱，導致其解碼效率仍然有限，阻礙了高頻SSVEP-BCI系統走向實際應用。因此，有必要開展面向高頻SSVEP信號的解碼算法研究，從而提高高頻SSVEP-BCI系統的識別性能。
近年來，SSVEP信號的解碼算法得到了快速發展。早期的SSVEP解碼多采用單導聯下的功率譜密度分析（power spectral-density analysis，PSDA）來識別信號頻率[12]。隨著典型相關分析（canonical correlation analysis，CCA）的提出，利用多導聯腦電數據進行空間濾波降維逐漸成為SSVEP識別的主流方式[13]。Chen等[14]采用了擴展CCA（extended CCA，eCCA）獲得具有特異性的空間濾波器和平均模板，最高ITR突破319.20 bits/min。任務相關成分分析（task-related component analysis，TRCA）進一步將SSVEP-BCI系統的ITR提升到了325.33 bits/min[15]。隨后，TRCA的拓展版本利用了額外信息來增強識別性能，例如鄰近頻率信息[16]、正余弦參考信號[17]、相位約束[18]等。最近，Liu等[19]提出任務判別成分分析（task-discriminant component analysis，TDCA）方法，被證實是當前最佳的SSVEP空間濾波算法之一。然而，這些先進算法大多在低頻SSVEP數據集上進行對比驗證，在高頻SSVEP信號上的識別性能仍然未知。
針對此問題，本文研究設計了刺激頻率范圍為30~39 Hz的十指令高頻SSVEP范式，離線采集了20名受試者的腦電數據，探究了6種先進SSVEP解碼算法在高頻SSVEP信號上的性能差異，為高頻SSVEP-BCI系統的算法選擇提供依據。
1 材料與方法
1.1 試驗對象
本文共招募了20名受試者參與離線腦電數據的采集，其中男性11名，女性9名，年齡在19~24周歲之間。所有受試者均視力正常或矯正視力正常，無任何精神疾病史，并自愿參與本試驗。在試驗正式開始前，所有受試者閱讀了試驗流程與注意事項，簽署了相關的知情同意書。本試驗通過了天津大學研究倫理委員會的批準（批文編號：TJUE-2021-062），所采集數據用于本文研究分析。
1.2 試驗范式
為探究高頻刺激下SSVEP解碼算法的識別性能，設計了刺激頻率范圍為30~39 Hz（步長為1 Hz）的高頻SSVEP范式，共包含10個刺激頻率。每個刺激為一個120×120像素的方塊，兩個相鄰刺激之間的垂直距離和水平距離均為350像素。刺激程序由商業數學軟件MATLAB R2019a（MathWorks Inc.，美國）中的Psychphysic toolbox v3.0.15工具包編寫完成，呈現在59.68 cm×33.57 cm的液晶顯示屏上，其分辨率為1 920×1 080，刷新率為120 Hz。
試驗共10個刺激指令，每個刺激以偽隨機序列依次閃爍1個試次為1輪，每5輪為1組試驗，每組試驗結束后讓受試者進行適當休息，防止視覺疲勞對試驗結果產生影響。每名受試者各進行了4組試驗，因此每個指令下均包含20個試次腦電數據。本文解碼過程為十分類任務，需對十個不同頻率下誘發的腦電信號進行識別。
正式試驗時，受試者端坐于顯示屏前，雙眼距離屏幕70 cm。調整為平靜狀態后，點擊“空格”鍵開始試驗。每個試次下刺激呈現持續時間為2 s，相鄰刺激之間的空隙時間（即刺激提示時間）為2 s，如圖1所示。試驗過程中，受試者在刺激呈現時間內盡量避免眨眼，并保持注意力集中。
圖1
試驗流程圖
Figure1.
Experimental flow chart
1.3 數據采集
本研究試驗全程在電磁屏蔽室中完成。腦電放大器（SynAmps, Neuroscan Inc., 美國）和配套采集軟件記錄64導腦電信號，導聯放置遵循國際10-20標準，采樣率為1 000 Hz，系統濾波頻帶設置為直流（direct current，DC）到200 Hz之間，以50 Hz陷波器去除工頻干擾，受試者的頭皮阻抗保持在10 kΩ以下。試驗過程中，刺激程序通過并口通訊向腦電信號放大器發送事件觸發信號，并進行儲存。
2 解碼算法
2.1 典型相關分析及其拓展算法
CCA是一種用于衡量多維隨機變量的潛在相關性的統計學方法，其本質是求解使得兩組隨機變量降維后相關性最大的最優線性組合。Bin等[13]將CCA應用于多導聯SSVEP信號的分類識別中，其性能顯著優于傳統的PSDA。作為一種無監督方法，CCA算法因其計算簡單、使用便捷而被廣泛應用于SSVEP-BCI領域。Chen等[14]引入個體訓練數據來構建CCA濾波器和個體平均模板，即eCCA。與傳統CCA算法相比，這種有監督方法考慮了腦電信號的跨個體、跨時間變異性，顯著提升了SSVEP-BCI系統的識別效率。在此基礎上，Wong等[16]設計了一種跨刺激學習框架，從目標刺激及其若干鄰近刺激中學習空間濾波器，提出了多刺激eCCA（multi-stimulus eCCA，ms-eCCA）算法，能夠進一步增強eCCA空間濾波器的魯棒性。跨刺激學習框架有效挖掘了不同頻率間SSVEP信號的共同信息，改善了eCCA算法在訓練樣本不足時的分類性能。
2.2 集成任務相關成分分析及其拓展算法
腦電信號可看作是不同任務相關源信號的線性組合，TRCA的目標是通過最大化相同任務下信號的相似性來提取任務相關成分[15]。由于相鄰刺激頻率下誘發的SSVEP信號具有相似的頭皮空間分布，因此對不同類別下空間濾波器進行集成可以增強算法的魯棒性，即集成TRCA（ensemble TRCA，eTRCA）算法。近年來，衍生出了多種eTRCA的擴展方法。例如，多刺激eTRCA（multi-stimulus eTRCA，ms-eTRCA）是跨刺激學習框架應用在eTRCA上的一種先進算法[16]。Wong等[17]首次將頻率信息引入到TRCA空間濾波器的構建中，開發了帶有正余弦參考信號的eTRCA（eTRCA with sine-cosine reference signal，eTRCA-R），其本質是將腦電信號投影到正余弦子空間后再求解空間濾波器。
2.3 任務判別成分分析
TDCA算法解決了傳統eTRCA中空間濾波器冗余的問題，是當前低頻SSVEP信號識別性能最佳的算法之一[19]。與eTRCA不同的是，TDCA對訓練數據進行空間維度和時間維度下的兩次增廣后，利用二維線性判別分析為不同類別信號求解一個共享的空間濾波器。
假設同一類別下的多通道腦電信號表示如式（1）所示：
|
其中，Nc為導聯數量、Np為采樣點數量，Nt為訓練試次數量，i為訓練試次的索引，代表實數集。首先，對所有訓練數據進行時延，得到提高空間維度后的一階增廣信號，如式（2）所示：
|
其中，表示腦電數據延后個點后的腦電數據，T為轉置符號。其次，將投影至對應類別頻率下參考信號的正交子空間，并將投影前后的信號按照時間維度進行拼接后得到二階增廣信號，如式（3）所示：
|
其中，A為二階增廣后的腦電信號，為正交投影矩陣[19]。
TDCA采用二維線性判別模型求解空間濾波器，計算所有類別下訓練數據的類間散度矩陣Sb和類內散度矩陣Sw，如式（4）、式（5）所示：
|
|
式中，Nf為訓練類別數量，代表第n個類別、第i個試次下的二階增廣信號，代表第n個類別下的模板信號，代表所有類別中心。
化簡后，TDCA優化目標如式（6）所示，W為需要求解的空間濾波器，為W的估計值。通過最優化問題推導可得，取的前Nk個特征向量即可構建TDCA空間濾波器。
|
將單次測試信號進行補零后，采用與訓練過程相同的步驟獲得二階增廣后的測試信號，計算與經過空間濾波后的相關系數，如式（7）所示，其最大值所對應的類別為識別結果。
|
2.4 性能評估
在預處理階段，首先將原始腦電信號降采樣至250 Hz，并提取枕區周圍典型的9個導聯（Pz、PO5、PO3、POz、PO4、PO6、O1、Oz、O2）數據用于后續分類識別[15, 19]。通過事件標志對原始信號進行截取得到[~（）]時刻的數據。其中，表示信號潛伏期時長，Tw是用于分類的時間窗長度（Tw = 0.2, 0.4, ···, 1 s）。
本文選擇了當前在低頻SSVEP-BCI中被廣泛用于對比分析的6種先進SSVEP解碼算法[19]，包括eCCA、ms-eCCA、eTRCA、ms-eTRCA、eTRCA-R和TDCA。所有算法均應用濾波器組（filter bank，FB）分析，子帶數量Nb設置為2，第m個子帶的下限和上限截止頻率分別為m×30 Hz和90 Hz，不同子帶融合的權重系數如式（8）所示[10, 20]：
|
其中，a(m)代表第m個子帶的權重值。
6種算法均采用模板匹配的方式來預測類別[15, 19]。具體來說，在計算出空間濾波器后，將測試腦電信號分別與每個類別下的平均信號（或參考信號）進行空間濾波降維后，求解相關系數作為決策值。在決策值組成的特征向量中尋找最大值對應的類別作為預測結果。采用隨機抽樣交叉驗證方法進行分類驗證，隨機抽取Nt 個輪次數據作為訓練集（Nt = 2, 4, …, 10），重復10次取平均值。計算分類準確率、曲線下面積（area under the curve，AUC）和ITR作為評估指標。ITR（以符號ITR表示）計算方式如式（9）所示，單位為bits/min。
|
其中，M表示指令類別數量，P為分類準確率，T為輸出單個指令所需的時間（單位：s），包括注視時間（Tw）和視線轉移時間，本研究將視線轉移時間設定為0.5 s。
研究使用統計產品與服務解決方案軟件（statistical product and service solutions，SPSS）（IBM Inc.，美國）進行統計學分析。通過三因素重復測量方差分析（three way repeated measures analysis of variance，Three-way RMANOVA）探究算法因素、時間窗長度因素和訓練集樣本數目因素對分類準確率、AUC和ITR的影響，并采用單因素重復測量方差分析（one way repeated measures analysis of variance，One-way RMANOVA）探究不同算法在固定時間窗或者訓練樣本條件下時性能差異是否具有統計學意義。如果數據不符合球形檢驗（Mauchly’s test of sphericity），則進行格林豪斯-蓋舍校正（Greenhouse-Geisser correction）。
3 結果分析
3.1 潛伏期估計評估
SSVEP信號的潛伏期是指視覺通路傳導引起的SSVEP響應與刺激相位之間的時滯，可通過擬合相位與刺激頻率的線性函數關系進行估計[21]。如圖2所示，顯示了Oz導聯下所有受試者的SSVEP平均相位隨刺激頻率的線性變化。通過進一步計算線性模型斜率得到的平均潛伏期為（132.99 ±5.61）ms。基于此，本文分析，在SSVEP信號的提取中應去除前0.13 s的潛伏期信號（ = 0.13 s）。
圖2
平均相位的變化趨勢
Figure2.
Variety of average phase
3.2 參數優化
采用網格搜索法對高頻SSVEP信號識別中不同解碼算法適用的參數進行優化，需要優化的參數包括ms-eCCA和ms-eTRCA中的鄰近頻率數量d，TDCA中的子空間數量k和延遲點數l[16, 19]。如圖3所示的是三種算法在不同參數條件下的歸一化準確率，ms-eCCA和ms-eTRCA選擇沿熱點圖縱坐標軸方向平均后最大值對應的d作為最優參數，TDCA選擇沿熱點圖橫、縱坐標軸方向平均后最大值對應的l和k作為最優參數。基于此，ms-eCCA和ms-eTRCA中的鄰近頻率數量分別設置為5和3，TDCA中的子空間數量和延遲點數分別設置為3和1。
圖3
三種分類算法的歸一化準確率
Figure3.
Normalized accuracy of the three classification algorithms
3.3 分類性能對比
為探究6種算法在高頻SSVEP數據上的分類性能，計算了在不同時間窗長度（Tw = 0.2, 0.4, …, 1.0 s）和不同訓練集樣本數量（Nt = 2, 4, …, 10）下的分類準確率、AUC及ITR指標，并對20名受試者的結果取平均值，如圖4所示。Three-way RMANOVA的結果表明，不同算法、不同時間窗長度和不同訓練集樣本數量均會影響算法識別性能，其差異具有統計學意義（P < 0.01, 具體效應值詳見附表1）；且算法與時間窗長度之間、算法與訓練集樣本數量之間的交互作用具有統計學意義（P < 0.01, 具體效應值詳見附表1）。
圖4
不同分類算法之間的表現比較
Figure4.
Performance comparison between different target recognition methods
對所有訓練集樣本數量下的結果進行平均，得到分類性能隨時間窗長度增加時的變化趨勢，如圖4左側子圖所示。One-way RMANOVA的結果表明，所有時間窗長度條件下6種算法的識別性能差異均具有統計學意義（AUC指標：Tw = 0.6 s時，F(1.5，20.3) = 5.7，P = 0.01；Tw = 0.8 s時，F (1.4，26.3) = 5.1，P = 0.02；Tw = 1.0 s時，F (1.4，26.1) = 4.7，P = 0.03；其余條件P < 0.01, 具體效應值詳見附表2）。6種算法的平均分類準確率和AUC隨時間窗長度增大而上升，平均ITR呈現出先上升后下降的趨勢，并在Tw = 0.4 s時達到峰值。具體來看，在所有時間窗長度條件下，eCCA算法的分類準確率、AUC和ITR最低；其次是eTRCA、ms-eTRCA以及eTRCA-R，其中ms-eTRCA以及eTRCA-R的識別性能均優于原始eTRCA；TDCA算法的分類性能最佳，在Tw = 0.4 s時平均ITR達到93.40 bits/min。在時間窗長度較小時，ms-eCCA算法分類結果低于eTRCA、ms-eTRCA、eTRCA-R和TDCA；在時間窗長度較大（Tw ≥ 0.6 s）時，分類性能超過eTRCA、ms-eTRCA和eTRCA-R，并且能夠達到與TDCA算法相當的水平。從該結果中可以看出，在6種算法中，TDCA算法在所有時間窗條件下綜合表現更優，而ms-eCCA算法對時間窗長度變化更加敏感。
對所有時間窗長度下的結果進行平均，得到分類性能隨訓練集樣本數量增大時的變化趨勢，如圖4右側子圖所示。One-way RMANOVA的結果表明，所有訓練集樣本數量條件下6種算法的識別性能差異均具有統計學意義（P<0.01, 具體效應值詳見附表3）。總體來看，6種算法的平均分類準確率、AUC和ITR隨訓練集樣本數量增大而上升。具體而言，在所有訓練集樣本數量條件下，eCCA算法的分類結果最低；在訓練集樣本數量不足（Nt ≤ 4）時，ms-eCCA的分類性能依次優于TDCA、eTRCA-R、ms-eTRCA和eTRCA；在訓練集樣本數量足夠（Nt ≥ 8）時，TDCA算法的分類準確率和ITR最高，其次是eTRCA-R、ms-eTRCA以及eTRCA，而ms-eCCA算法僅優于eCCA。該結果表明，相較于其他5種算法而言，ms-eCCA算法對訓練集樣本數量不敏感，在小樣本情況下具有明顯的性能優勢；而在訓練集樣本數量足夠時，TDCA算法能夠達到更優的分類效果。
3.4 算法速度對比
在實際應用中，算法的建模速度和解碼速度同樣是影響BCI性能的關鍵因素。其中，建模速度由構建模型的運算過程決定，會對系統的準備時間造成影響，而單次解碼速度能夠體現在線BCI系統的實時效率。采用具有3.20 GHz中央處理器（central processing unit，CPU）的計算機測試了6種算法在不同訓練集樣本數量下的平均建模時間，以及在不同時間窗長度下的平均單次解碼時間，如圖5所示。
圖5
不同分類算法在不同訓練集樣本數量下的建模時間和不同時間窗長度下的解碼時間
Figure5.
Training time of different recognition methods at different training sample and decoding time at different time window lengths
在建模速度方面，隨著訓練集樣本數量的增加，所有算法的建模時間均呈現上升趨勢。其中，eCCA的訓練過程僅需要獲得個體模板，因此建模速度最快，而ms-eCCA與TDCA算法的建模時間明顯高于其余算法。由于ms-eTRCA和eTRCA-R將算法核心步驟中協方差求解過程進行優化，避免了無用且耗時的循環過程，因此建模時間也明顯低于原始eTRCA算法。
在解碼速度方面，隨著時間窗長度增加，所有算法的解碼時間逐漸增加。其中eTRCA、ms-eTRCA以及eTRCA-R三種算法解碼過程一致，因此解碼速度相同。在所有時間窗長度條件下，這3種算法計算開銷最少，而eCCA解碼速度最慢。當時間窗長度較小（Tw = 0.2 s）時，TDCA的單次解碼時間比ms-eCCA短；而當時間窗長度逐漸增大（Tw ≥ 0.4 s）后，ms-eCCA的解碼速度優于TDCA。
4 討論
本文離線采集了20名受試者的高頻十分類SSVEP數據，探究了包括eCCA、ms-eCCA，eTRCA、ms-eTRCA、eTRCA-R和TDCA在內的6種先進SSVEP解碼算法在不同訓練集樣本數量、不同時間窗長度下對高頻SSVEP信號的解碼性能差異，并比較了不同算法的建模和解碼速度。
本文研究首先在高頻SSVEP信號上優化了先進SSVEP解碼算法的相關參數，為高頻SSVEP-BCI研究中先進算法的參數選擇提供了參考。其次，通過計算分類準確率、AUC和ITR三個指標來評估算法的分類性能。結果表明，6種算法的分類性能排序可以表示為：TDCA≈ms-eCCA > eTRCA-R > ms-eTRCA > eTRCA > eCCA。其中，與其他算法相比，ms-eCCA算法對時間窗長度更敏感，但對訓練集樣本數量不敏感。因此，在訓練集樣本數量不足時，ms-eCCA算法優于研究者最新提出的TDCA算法。然而，雖然TDCA和ms-eCCA具有更高的分類性能，但由于引入了更復雜的求解過程，因此建模速度和解碼速度明顯低于eTRCA、ms-eTRCA和eTRCA-R算法。在實際應用中，應結合實際情況選擇適合的分類算法構建高頻SSVEP-BCI系統。
表1總結了近年來已有研究中報道的高頻SSVEP-BCI系統，從頻率范圍、指令數量、識別算法、ITR四個方面進行比較。可以看出，現有高頻SSVEP-BCI系統研究大多集中在30~40 Hz頻段，且超過半數使用(e)TRCA進行分類識別[7-8, 10-11, 22-29]。本文研究結果表明，ms-eCCA、ms-eTRCA、eTRCA-R和TDCA在高頻SSVEP信號上同樣具有與eTRCA相當或更優的識別性能，有望在后續被應用于高頻SSVEP-BCI系統中。具體而言，考慮分類性能的情況下，在訓練集樣本數量不足時優先選擇能夠抵御小樣本問題的ms-eCCA算法，在訓練集樣本數量充足時優先選擇能夠獲得最優識別結果的TDCA算法；考慮解碼速度的情況下，在保證分類性能的基礎上推薦eTRCA-R算法。
表1
已有高頻SSVEP-BCI系統比較
Table1.
Comparison among high-frequency SSVEP-BCI systems
與具有強閃爍感的低頻閃爍刺激相比，高頻刺激具有安全、舒適的優勢，因此有助于開發面向實際應用場景的SSVEP-BCI系統[8]。然而，高頻刺激下誘發的SSVEP信號微弱，識別更加困難。本文研究結果可為高頻SSVEP-BCI系統的算法選擇以及相應的參數選擇提供參考，從而提升高頻SSVEP信號的識別效果。除了解碼方法外，未來還可在編碼范式上對高頻SSVEP-BCI進行優化。例如，融合頻分多址[10]、時分多址[11]、空分多址[30]、碼分多址[31]等多維度編碼特征來提升指令間的可分性。進一步地，將高頻頻率拓展到高于人眼閃爍融合頻率的范圍，可構建舒適無感SSVEP-BCI系統[26, 28]。
5 結論
為了衡量先進SSVEP解碼算法在高頻SSVEP信號上的識別性能，本文研究采集了20名受試者的高頻SSVEP數據，以分類準確率、AUC值、ITR、算法速度作為評價指標在不同訓練集樣本數量、不同時間窗長度條件下展開對比。在6種算法中，TDCA和ms-eCCA算法表現出最佳的分類性能，其中ms-eCCA算法在訓練集樣本數量不足時具有明顯優勢。eTRCA-R算法分類性能僅次于前兩者，但解碼速度更快，能夠保證在線系統的實時效率。本文結論為高頻SSVEP研究的算法選擇提供了理論依據，在構建舒適友好型SSVEP-BCI系統方面具有潛在的應用價值。
重要聲明
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
作者貢獻聲明：羅睿心負責試驗設計、程序編寫以及論文撰寫；豆心怡負責數據分析、數據整理以及論文撰寫；肖曉琳、吳喬逸、許敏鵬負責論文指導與修改；明東參與論文的審閱。
倫理聲明：本研究通過了天津大學研究倫理委員會的審批（批文編號：TJUE-2021-062）。
本文附件見本刊網站的電子版本（biomedeng.cn）。
