基于回溯搜索算法的導聯(lián)選擇腦機接口研究

戴圣法，魏慶國，魏中海

（南昌大學 信息工程學院，江西 南昌　330031）

基于回溯搜索算法的導聯(lián)選擇腦機接口研究

戴圣法，魏慶國，魏中海

（南昌大學 信息工程學院，江西 南昌330031）

在腦機接口（BCI）中，傳統(tǒng)的共空域模式（CSP）算法在提取特征信號與事件相關去同步/同步（ERD/ERS）的信息上得到了很好的效果。但是CSP算法受限于電極導聯(lián)數(shù)、EEG信號的時間段和頻帶等因素，如電極導聯(lián)數(shù)的增加，CSP算法容易過擬合，數(shù)據(jù)記錄容易混亂，使得運算變得復雜，增加運算時間，降低數(shù)據(jù)分類正確率。所以，CSP算法存在局限性。使用回溯搜索優(yōu)化算法（BSA）能夠為CSP算法自動挑選出一組導聯(lián)數(shù)組子集，并且以分類錯誤率作為BSA算法的目標函數(shù)進行實驗。實驗采用兩類實驗數(shù)據(jù)（第三、四屆國際BCI競賽數(shù)據(jù)集）進行交叉驗證分類實驗。實驗結果表明，兩類數(shù)據(jù)的導聯(lián)數(shù)目大幅度減少，分類正確率有所提高。

腦機接口；共空域模式；回溯搜索優(yōu)化算法；最優(yōu)導聯(lián)選擇

0　引　言

腦機接口（Brain Computer Interface，BCI）是一種直接利用腦信號連接和控制外界設備，不依靠人體神經和肌肉的通信系統(tǒng)［1］。通過BCI系統(tǒng)，將人腦的想法、思維轉化為控制外部設備的指令，可以幫助那些患有肌肉萎縮性側索硬化、腦干中風等疾病的人實現(xiàn)與外界交流［2］。

當人想象某個精神任務或者執(zhí)行某個肢體運動時，大腦皮層的某個區(qū)域中腦電信號會發(fā)生變化，這類現(xiàn)象伴隨著腦電信號能量的減小或增加。Pfurtscheller將腦電信號能量的減小稱為事件相關去同步（Event-related Desynchronization，ERD），而能量的增加稱為事件相關同步（Event-related Synchronization，ERS）［3-4］。在EEG節(jié)律信號中mu節(jié)律和beta節(jié)律信號是ERD/ERS現(xiàn)象有關的兩種固有頻率信號，mu節(jié)律信號的頻率范圍一般為8～12 Hz，beta節(jié)律信號［5］為18～26 Hz。

傳統(tǒng)的共空域模式（Common Spatial Pattern，CSP），在ERD/ERS相關的特征提取上是很有效的算法。由于ERD/ERS信號發(fā)生在特定的腦區(qū)域，而采集信號的導聯(lián)分布于整個大腦區(qū)域，所以導聯(lián)通道中存在許多無關的通道。因此，CSP的性能會受到限制［6］。

本文采用回溯搜索優(yōu)化算法（Backtracking Search Optimization Algorithm，BSA）與CSP算法相結合的方法。在使用CSP算法之前盡可能地挑出有用的腦電信號的導聯(lián)子集。實驗中采用第三屆國際腦機接口競賽數(shù)據(jù)集（BCI Competition III Dataset Iva）［7］和第四屆國際腦機接口競賽數(shù)據(jù)集（BCI Competition IV Dataset IIa）［8］進行分類實驗。實驗結果表明該方法比單純使用CSP算法的分類識別率更高。

1　方　法

在腦電信號處理的實驗中，采用切比雪夫I型8～30 Hz（包含mu節(jié)律和beta節(jié)律）的帶寬IIR帶通濾波器［9］對多通道EEG信號進行濾波，然后用BSA算法選擇最優(yōu)導聯(lián)，對選出的導聯(lián)使用CSP算法進行特征提取，最后線性判別分析（Linear Discriminant Analysis，LDA）分類器［10］對特征信號進行分類。并用10×5倍交叉驗證之后的分類錯誤率作為BSA的目標函數(shù)評判選出的導聯(lián)數(shù)集的優(yōu)劣。BSA與CSP結合的導聯(lián)選擇算法流程如圖1所示。

圖1　BSA與CSP結合的導聯(lián)選擇算法流程

1.1共空域模式算法

CSP算法基于兩個實對稱矩陣的聯(lián)合對角化，并進行空間-時間源建模［11-12］。把原始的EEG信號分解為空域模式，通過空域模式使兩類信號的差別達到最大化，從而提取與任務相關的信號分量。具體分為如下三個步驟：

（1）求解腦電信號的協(xié)方差矩陣，建立空間-時間源模型。實驗時，受試者被要求執(zhí)行兩種不同的想象任務a和b。想象任務a和b誘發(fā)的多通道EEG信號為空間-時間信號矩陣Xa和Xb，并按式（1）建模：

式中：Sa和Sb分別是與任務a和b有關的源分量，Ca和Cb是Sa和Sb的空域模式；Sc是共同條件下的源分量，Cc是Sc的空域模式。

設EEG信號記錄的通道數(shù)為N，每個導聯(lián)信號的采樣點數(shù)為T，矩陣Xa和Xb的維數(shù)都是N×T。任務a和b對應的腦電信號的空間協(xié)方差矩陣為：

式中：T表示轉置運算；trace（X）是求矩陣X的主對角元素之和；Ra和Rb都是N×N的矩陣。

（2）協(xié)方差聯(lián)合對角化。對兩個協(xié)方差矩陣之和R進行特征分解：

式中：U0是特征向量矩陣；Λ是特征值矩陣，矩陣R與Λ的維數(shù)都是N×N。在以下的討論中，如無特別說明，都假設特征值不為零。則白化變換的矩陣P可用式（4）表示：

對Ra和Rb分別進行白化變換后得：

Sa和Sb具有如下兩個重要性質：

① Sa和Sb具有共同的特征向量：

② Sa和Sb對應的特征值之和為I：

式中：Λa和 Λb分別是Sa和Sb的特征值矩陣；I為單位矩陣。設Λa的對角元素以降序排列，任務a的初始幾個空間因子占有的方差被最大化，那么任務b相應的方差被最小化。這兩類任務可分別由 Λa和 Λb中最大的m個特征值對應的特征向量表示，即Uam和Ubm表示。由于特征值矩陣對角元素是遞減排列的，所以在特征矩陣中的第一個和最后一個特征向量是識別兩類任務的最佳特征向量。

（3）構造空域濾波器。用Uam構建空域濾波器Fa，用Ubm構建空域濾波器Fb，其中m?N，F(xiàn)a和Fb如下所示：

然后將記錄的EEG數(shù)據(jù)分別通過構建好的兩個空域濾波器進行空域濾波，得到這兩類信號的源信號：

式中：F的行向量表征的是濾波器的系數(shù)；F-1的列向量表征的是共空域模式矩陣，亦被稱為時不變腦電信號源分布向量。當F-1投影到腦部時，可直觀地顯示腦神經活動的分布情況。對于類別未知的測試集數(shù)據(jù)，可以通過由訓練集數(shù)據(jù)訓練好的分類器實現(xiàn)數(shù)據(jù)類別的分類。相應的特征可在源信號的基礎上構造為：

則特征向量可表示為 f=［fa，fb］，式（10）中的對數(shù)運算是為了讓特征向量 fa和 fb中的元素更接近于正態(tài)分布。關于CSP算法的計算步驟和詳細描述可參閱文獻［13］。

1.2回溯搜索優(yōu)化算法

BSA是2013年由Cicicioglu提出的一種進化算法［14］。該算法通過自身產生實驗種群，控制搜索方向和搜索邊界。在整個過程中只有一個控制參數(shù)，因此該算法相比其他進化算法操作簡單。而且該算法同其他進化算法的步驟類似，分別為種群初始化、歷史種群選擇、種群突變、種群交叉和最后種群輸出。

（1）種群初始化

實驗中采用隨機產生種群的方法進行初始化，即：

式中：P為種群；i∈［1，2，…，N］，j∈［1，2，…，D］，N是種群大小，D是問題維數(shù)；low和up分別為搜索區(qū)間的下限和上限；U是隨機均勻分布。Pi是種群P中的目標個體。

（2）歷史種群選擇

初始化的歷史種群定義如下：

BSA算法進行每一次迭代運算時，當歷史種群a小于隨機新種群b時 old P=P，a，d～U（lowj，upj），即更新歷史種群。

（3）種群突變

通過式（13）產生突變種群的最初形式為：

式中：F為控制搜索方向矩陣（old P-P）幅度的參數(shù)，F(xiàn)=3·rndn，rndn～N（0，1），N是正態(tài)分布。

（4）種群交叉

在BSA算法種群交叉過程中，產生實驗種群T的最后形式。在實驗種群的初始值不斷的突變中，適應于優(yōu)化的實驗個體用來進化目標種群個體。種群交叉過程分為兩個步驟：

①計算一個二進制維數(shù)為N×D的map矩陣，map矩陣的初始值為［0］N×D，通過利用種群P的相關個體操縱實驗種群 T的個體。當 mapi，j=1，i∈［1，2，…，N］，j∈［1，2，…，D］時，Ti，j=Pi，j，實現(xiàn)了實驗種群T個體的更新。

②通過式（14）進行交叉：

式中：rndn～N（0，1），N為正態(tài)分布；mixrate為控制種群個體元素的數(shù)目參數(shù)，mixrate=1；“┌┐”是上限取整符號；randi（D）為從［0，D］中隨機取一個整數(shù)；c，d～U（0，1），U（0，1）表示在0～1均勻分布；u為隨機排序，且 u∈［1，2，…，D］的整數(shù)向量。

BSA算法通過mixrate控制新種群T個體中元素的個數(shù)，同時利用和randi（D）進行實驗突變和突變結束。當c＜d時，mapi為多個具有隨機位置的向量；反之，mapi為僅有為0的向量。在交叉過程之后產生的種群中，超出被允許的搜索空間范圍的部分個體是突變策略的結果。按照式（11）產生新的種群。

（5）最后種群輸出

通過貪婪選擇機制，相比對應初始種群P中有選擇適應度值（目標值）較好的種群個體，種群進行更新。當種群P中的個體比全局最小值有更好的適應值時，記錄和輸出當前最優(yōu)種群的個體元素，同時更新初始種群，完成一次迭代。重復上述過程，直至到最大迭代次數(shù)，最后輸出全局最小值。

1.3基于BSA的導聯(lián)選擇方法

原始EEG信號通過8～30 Hz帶寬帶通濾波器進行濾波，然后使用CSP算法對信號進行特征提取。由于CSP的性能與導聯(lián)關系密切，導聯(lián)的選擇是一個關鍵步驟。在BSA算法中設置兩個參數(shù)為：下限 low=［0，0，…，0］1×N，上限up=［1，1，…，1］1×N。這是一個N維的電極優(yōu)化選擇向量，向量每一個分量代表一個導聯(lián)。設其中的一個分量為ci，i∈［1，2，…，N］，那么ci～ U（0，1）。另設一串二進制編碼為C1×N，滿足Ci=［ci］，i∈［1，2，…，N］。其中U（0，1）表示在0～1均勻分布，“［］”為取高斯整。同時設置種群大小 N=50，問題維數(shù)D=2。按C中的編碼抽取導聯(lián)，數(shù)字編號為1的導聯(lián)挑出，數(shù)字編號為0的導聯(lián)舍棄。BSA算法在進化中不斷地產生隨機數(shù)組合，挑選出不同的導聯(lián)子集。整個過程以分類錯誤率作為評判準則，當?shù)螖?shù)達到設定的迭代次數(shù)時，BSA算法就會停止，并且輸出全局最小值，即分類錯誤率最小，同時輸出對應的最優(yōu)導聯(lián)子集。

2　實驗結果及數(shù)據(jù)分析

2.1實驗數(shù)據(jù)的預處理

本文采用的數(shù)據(jù)是第三屆國際競賽腦機接口數(shù)據(jù)集（BCI competition III dataset IVa）和第四屆國際競賽腦機接口數(shù)據(jù)集（BCI competition IV dataset IIa）。

第三屆國際競賽腦機接口數(shù)據(jù)集是由5個受試者（Aa，Al，Av，Aw，Ay）進行BCI運動想像實驗所得到的數(shù)據(jù)集。在實驗中受試者分別執(zhí)行想象右手和腳運動兩種想象任務。EEG信號由118個電極記錄，電極分布如圖2（a）所示。每個受試者進行280次實驗，兩種想象任務實驗的次數(shù)相同，各為140次。采集的實驗數(shù)據(jù)經過放大后再用0.05～200 Hz的帶通濾波器濾波，且數(shù)據(jù)采樣率為1 kHz。在本實驗中，對原始記錄的EEG數(shù)據(jù)進行下采樣，下采樣率為100 Hz。在時域濾波之前，對原始連續(xù)的實驗數(shù)據(jù)截斷成單次實驗數(shù)據(jù)。

第四屆國際競賽腦機接口數(shù)據(jù)集是由9個受試者（A1～A9）進行BCI運動想象實驗所得到的數(shù)據(jù)集。在實驗中受試者分別執(zhí)行想象左手、右手、腳和舌頭四種想象任務。EEG信號由22個電極記錄，電極的分布如圖2（b）所示。

圖2　EEG信號電極分布圖

該數(shù)據(jù)集包含兩組數(shù)據(jù)（T session，E session），由兩個不同的時間段完成。每組實驗分為6個時間段進行，每個時間段受試者進行了48次運動想象實驗，其中每類實驗次數(shù)為12，一組實驗共包括6×48=288次單次運動想象實驗。電極記錄的EEG信號都要經過250 Hz的采樣及0.5～100 Hz的低通濾波，經過放大后保存。在實驗過程中，除了記錄EEG信號的22個電極外，還有三個單極性的電極用來記錄眼電EOG信號。與EEG信號一樣，EOG信號也要經過放大處理、采樣（采樣率為250 Hz）和0.5～100 Hz的帶通濾波。在本實驗中先將數(shù)據(jù)進行8～30 Hz的寬帶濾波，按時間段為2 s（2.5～4.5 s），采樣率為250 Hz截取數(shù)，從四類數(shù)據(jù)中抽取左手和右手兩類數(shù)據(jù)進行實驗。

2.2交叉驗證分類結果

為了能夠合理的得到分類結果，在實驗中采用10× 5倍交叉驗證。10× 5倍交叉驗證就是將一個受試者的數(shù)據(jù)集隨機排列10次，每一次隨機排列的數(shù)據(jù)被分為5個相等的部分；其中一個部分用于測試，而其余4個部分用于訓練分類器。這個交叉驗證步驟共有50次分類測試，求50個測試結果的平均值。在算法性能的評價上，采用兩種方法進行實驗，分類結果進行對比，即：全導聯(lián)的CSP算法提取特征，其流程如圖3所示；采用BSA導聯(lián)優(yōu)化與CSP相結合的算法提取特征，其流程如圖1所示。

圖3　CSP算法流程

BSA是一種基于種群的進化遺傳算法。在實驗中對種群進化迭代次數(shù)的選擇做了進一步實驗。改變不同的迭代次數(shù)并觀察實驗結果的變化。設置迭代次數(shù)為2，5，10，20，30，40，50，60，70。

5個受試者的分類正確率與迭代次數(shù)的關系如圖4所示。隨著迭代次數(shù)增加，分類正確率緩慢增長。當?shù)螖?shù)為60時，分類效果最好，最終設置迭代次數(shù)為60。9個受試者的平均分類正確率與迭代次數(shù)的關系如圖4所示?？紤]到迭代次數(shù)越多，運算時間越長，最終兩個時間段的數(shù)據(jù)選擇設置迭代次數(shù)為40，此時分類效果較好。

圖4　BSA導聯(lián)選擇5個受試者與9個受試者平均分類正確率與迭代次數(shù)的關系

表1表示在迭代次數(shù)為60時，5個受試者采用BSA最優(yōu)導聯(lián)及全部導聯(lián)的分類正確率及方差。比較兩種方法的結果可以得出，5個受試者使用BSA進行導聯(lián)選擇的平均分類正確率比使用全部118導聯(lián)要高出5%，尤其受試者Aa提升了10%，效果非常明顯。而且，5個受試者的導聯(lián)數(shù)目也減少了一半左右，數(shù)目減少十分突出。

表1　五個受試者采用BSA最優(yōu)導聯(lián)及全部導聯(lián)的分類正確率（%）及方差

表2表示在迭代次數(shù)為40時，9個受試者采用BSA最優(yōu)導聯(lián)及全部導聯(lián)的分類正確率及方差。比較兩種方法的結果可以得出，9個受試者的平均分類正確率分別提高了5%（T session）和8%（E session）。同樣，9個受試者的導聯(lián)數(shù)總體上都減少了，甚至部分受試者的導聯(lián)數(shù)減少了一半。

表1和表2中的數(shù)據(jù)表明，BSA算法對導聯(lián)的選擇和優(yōu)化十分有效，對基于運動想象的兩分類BCI系統(tǒng)有一個顯著的提升，具有重要意義。

表2　九個受試者采用BSA最優(yōu)導聯(lián)及使用全部導聯(lián)的分類正確率（%）及方差

3　結　語

在BCI研究中，傳統(tǒng)的CSP算法對兩種不同的運動想象任務的判別，得到了比較好的結果，但是在實驗中不同的受試者，其最佳的導聯(lián)分布，濾波頻帶和數(shù)據(jù)時間段也不相同。而基于傳統(tǒng)CSP算法的BCI系統(tǒng)中采用固定的導聯(lián)、頻帶和時間段，所以傳統(tǒng)CSP算法存在局限性。本文中，在應用CSP對多通道EEG信號提取特征之前，先使用BSA算法自動地選擇恰當?shù)膶?lián)，然后進行交叉驗證實驗。實驗結果表明，這個方法是十分有效的。

設計一個實際應用的運動想象BCI系統(tǒng)，導聯(lián)通道數(shù)目應該盡可能小，最好的做法是選擇最優(yōu)導聯(lián)子集代替整個導聯(lián)集。在本文中，遺傳優(yōu)化算法BSA應用于兩個數(shù)據(jù)集，即第三屆國際競賽腦機接口數(shù)據(jù)集和第四屆國際競賽腦機接口數(shù)據(jù)集。實驗結果表明，BSA可以選擇導聯(lián)數(shù)量較少，并保持分類精度提高。CSP的性能還受到其他因素影響，如頻帶和時間窗長度。因此，研究一種能同時考慮多個因素的方法是未來的趨勢。

［1］WOLPAW J R，BIRBAUMRE N，MCFARLAND D J，et al. Brain-computer interfaces for communication and control［J］. Clinical neurophysiology，2002，113（6）：767-791.

［2］NICOLAS-ALONSO L F，GOMEZ-GIL J.Brain computer interfaces，a review［J］.Sensors，2012，12（2）：1211-1219.

［3］TORO C，DEUSCHL G，THATCHER R，et al.Event-related desynchronization and movement-related cortical potentials on the ECoG and EEG［J］.Electroencephalography and clinicalneurophysiology/evokedpotentialssection，1994，93（5）：380-389.

［4］BABILONI C，CARDUCCI F，CINCOTTI F，et al. Human movement-related potentials vs desynchronization of EEG alpha rhythm：a high-resolution EEG study［J］.Neuroimage，1999，10（6）：658-665.

［5］HUANG S J，WU X M.Feature extraction of electroencephalogram for imagery movement based on Mu/ Beta rhythm［J］.Journal of clinical rehabilitative tissue engineering research，2010，43（14）：8061-8064.

［6］NOVI Q，GUAN C，DAT T H，et al.Sub-band common spatial pattern（SBCSP）for brain-computer interface［C］//Proceedings of 2007 3rd International IEEE/ EMBSConferenceonNeuralEngineering.［S.l.］：IEEE，2007：204-207.

［7］BLANKERTZ B.BCI competitionⅢ webpage［EB/OL］.［2005-11-03］.http：//www.bbci.de/competition/iii.

［8］BBCI.BCI competition IV［EB/OL］.［2011-07-03］.http：//www. bbci.de/competition/iv/.

［9］DORNHEGE G，BLANKERTZ B，KRAULEDAT M，et al. Combined optimization of spatial and temporal filters for improving brain-computer interfacing［J］.IEEE transactions on biomedical engineering，2006，53（11）：2274-2281.

［10］GARRETT D，PETERSON D A，WANDERSON C.Comparison of linear，nonlinear，and feature selection methods for EEG signal classification［J］.IEEE transaction on neural systems and rehabilitation engineering，2003，11（2）：141-144.

［11］WANG Y，BERG P，SCHERG M.Common spatial subspace decomposition applied to analysis of brain responses under multiple task conditions：a simulation study［J］.Clinical neurophysiology，1999，110（4）：604-614.

［12］SAMEK W，VIDAURRE C，MüLLER K R，et al.Stationary common spatial patterns for brain-computer interfacing［J］. Journal of neural engineering，2012，9（2）：13-15.

［13］MüLLER-GERKING J，PFURTSCHELLER G，F(xiàn)LYVBJERG H.Designing optimal spatial filters for single-trial EEG classification in a movement task［J］.Clinical neurophysiology，1999，110（5）：787-798.

［14］CIVICIOGLU P.Backtracking search optimization algorithm for numerical optimization problems［J］.Applied mathematics and computation，2013，219（15）：8121-8144.

Study on brain-computer interface with lead selection based on backtracking search optimization algorithm

DAI Shengfa，WEI Qingguo，WEI Zhonghai
（College of Information Engineering，Nanchang University，Nanchang 330031，China）

In brain-computer interface（BCI），the traditional common spatial pattern（CSP）algorithm has a good effect on characteristic signal extraction and event-related desynchronization/event-related synchronization（ERD/ERS）information.The CSP algorithm is easily limited by electrode lead quantity，time period and frequency band of EEG signal，such as the increasing of electrode lead quantity，easy overfitting of CSP algorithm and easy chaos of data record，which can make the operation complex，increase the operation time and reduce the accuracy of data classification.Therefore，the CSP algorithm have a limitation. The backtracking search optimization algorithm（BSA）proposed in this paper can automatically select a subset of lead array for CSP algorithm，and take the classification error rate as the objective function of BSA algorithm to test.The two datasets from the datasets of the Third，F(xiàn)ourth International BCI Competitions are adopted in the experiment to perform the classification experiment of cross validation.The experimental results show that the lead quantity of the two datasets are dramatically reduced，and the classification accuracy is improved.

brain-computer interface；common spatial pattern；backtracking search optimization algorithm；optimal lead selection

TN911-34；TP301

A

1004-373X（2016）13-0010-05

10.16652/j.issn.1004-373x.2016.13.003

2015-09-28

國家自然科學基金項目（61365013）；江西省自然科學基金項目（2009GZS0073）

戴圣法（1988—），男，江西南昌人，碩士研究生。主要研究方向為通信與信息系統(tǒng)。

魏慶國（1963—），男，江西南昌人，教授，博士。主要研究方向為信號處理、模式識別、機器學習、腦機接口。

魏中海（1988—），男，江西南昌人，碩士研究生。主要研究方向為信號與信息處理。