王 康，翟弟華，夏元清

（北京理工大學(xué)自動化學(xué)院，北京100081）

1 引 言

腦機(jī)接口（Brain Computer Interface，BCI）可以在人腦與外界的環(huán)境設(shè)備之間建立直接交流的通道，從而實現(xiàn)大腦對外部設(shè)備的交互控制［1］。因可以在控制系統(tǒng)中融合人與機(jī)器的智能，腦機(jī)接口技術(shù)在軍事、醫(yī)療、娛樂等領(lǐng)域中展現(xiàn)了廣闊的前景。21 世紀(jì)初美國就在意念控制機(jī)器人士兵項目中投入了大量的研究經(jīng)費。2004年，美軍資助了多個研究機(jī)構(gòu)開展“思維控制機(jī)器人”項目的研究。2009年，多名美國著名學(xué)者撰寫了《Opportunities in Neuroscience for Future Army Application》［2］，將腦機(jī)接口列為軍事領(lǐng)域重點發(fā)展技術(shù)。2013年，美國國防部提出“阿凡達(dá)”項目，擬實現(xiàn)意念遠(yuǎn)程操控“機(jī)器戰(zhàn)士”，代替士兵在戰(zhàn)場上執(zhí)行作戰(zhàn)任務(wù)［3］。近年，我國在此領(lǐng)域也加大研發(fā)，并取得了不錯的進(jìn)展。如國防科技大學(xué)研制成功腦控汽車，可以根據(jù)人腦的思維意識完成啟動、轉(zhuǎn)彎和加速等復(fù)雜操作［4］?？梢韵嘈?，腦機(jī)接口技術(shù)可為機(jī)器人、裝甲車等武器裝備進(jìn)行智能化賦能，引發(fā)武器裝備在操控方面的革命。

腦機(jī)接口系統(tǒng)主要包括腦電信號的采集與預(yù)處理、特征提取、分類識別和輸出控制。目前，腦機(jī)接口的輸入信號主要采用統(tǒng)計特征較為明顯的腦電信號，其中運動想象腦電信號通過檢測大腦μ 節(jié)律和β 波的能量對人腦的意圖進(jìn)行判斷［5］，被廣泛應(yīng)用于腦機(jī)接口系統(tǒng)。特征提取與分類識別是實現(xiàn)腦機(jī)接口系統(tǒng)的核心技術(shù)［6］，識別結(jié)果可以轉(zhuǎn)換成可被計算機(jī)控制指令。運動想象腦電信號常用的特征提取算法包括功率譜估計、自適應(yīng)自回歸模型、小波變換和共空間模式等［7-9］，常用的分類識別算法包括線性判別分析、隱馬爾科夫模型、支持向量機(jī)和人工神經(jīng)網(wǎng)絡(luò)等［10-11］。目前，運動想象腦電信號在二分類任務(wù)可以取得良好的結(jié)果，基于共空間模式的特征提取算法在右手和舌頭的實驗中可以取得高達(dá)97%的正確率［12］。而在多分類任務(wù)中，分類結(jié)果表現(xiàn)并不理想，李昕等［13］利用小波變換和支持向量機(jī)在多任務(wù)的分類識別中最優(yōu)分類準(zhǔn)確率僅為69%；施錦河等［14］利用共空間模式、Hilbert算法和支持向量機(jī)在多分類任務(wù)中的最優(yōu)準(zhǔn)確率為91%，但實驗使用了60 個腦電電極，準(zhǔn)確率在實驗過程中存在波動，穩(wěn)定性不足；天津大學(xué)的萬柏坤等［15］利用了二維的時頻域分析、Fisher 判別法和支持向量機(jī)進(jìn)行了多分類任務(wù)的研究，最優(yōu)準(zhǔn)確率達(dá)到了86%，但使用了60個電極，操作較為復(fù)雜，缺乏實用性。

考慮到單一算法的局限性，本文采用了小波包分解算法和一對多共空間模式對腦電信號進(jìn)行聯(lián)合特征提取，將原始信號分解為多個子頻帶，提取了細(xì)化的時頻域特征，充分考慮了運動想象腦電信號的相關(guān)作用頻帶，并使用人工神經(jīng)網(wǎng)絡(luò)對提取的特征進(jìn)行分類。此外，盡管一些文獻(xiàn)提出的算法展現(xiàn)出了良好的效果，但由于受到測試對象和環(huán)境等因素的影響，分析的腦電數(shù)據(jù)存在差異，缺乏統(tǒng)一性，難以得出較為客觀的比較結(jié)果。因此，本文使用“BCI IV”競賽中提供的公共標(biāo)準(zhǔn)數(shù)據(jù)Dataset 2a對算法進(jìn)行了驗證，并與競賽前幾名的數(shù)據(jù)進(jìn)行比較，結(jié)果表明，本文算法的識別準(zhǔn)確率有著比較明顯的提升。

2 算法原理

2.1 小波包分解

小波包分解［16］（Wavelet Packet Decomposition，WPD）是在小波變換（Wavelet Transformation，WT）的基礎(chǔ)上優(yōu)化得到的信號分解與重構(gòu)方法。小波變換可以對信號進(jìn)行局部化與多尺度分析，在信號分析中具備一定的自適應(yīng)性。但是，小波變換只是對信號的低頻部分進(jìn)行分解與重構(gòu)，而小波包變換考慮了信號的高頻部分，能夠有效降低信號冗余，克服了小波變換的固有缺點。因此，小波包變換對非平穩(wěn)的神經(jīng)生物信號、機(jī)械振動信號等可以進(jìn)行更好的時、頻域分析。

小波包分解使用低通濾波器hn(n ∈Z)和高通濾波器gn(n ∈Z)對原始的非平穩(wěn)信號進(jìn)行濾波，其中g(shù)n=(-1)nh1-n，二者為正交關(guān)系。正交尺度函數(shù)φ(t)和小波函數(shù)φ(t)滿足：

考慮到公式表述的一般性，本文將正交尺度函數(shù)φ(t)和小波函數(shù)φ(t)均表示為μ(t)。在尺度0上，尺度函數(shù)為μ0，0(t)，尺度1 上的正交尺度函數(shù)和小波函數(shù)分別為μ1，0(t)和μ1，1(t)，則式（1）在尺度2上可表示為：

對于尺度j，遞推可得到表達(dá)式為：

其中，函數(shù)系μj，n表示小波函數(shù)φ(t)的小波包，圖1為三層小波包分解的示意圖。

2.2 共空間模式

圖1 小波包分解圖Fig.1 Wavelet packet decomposition diagram

共空間模式［17］（Common Spatial Patterns，CSP）是一種典型的運動想象腦電信號特征提取算法，在二分類任務(wù)中取得了良好的效果。CSP算法的原理是通過設(shè)計空間濾波器，使一類信號的方差最大化，同時使另一類信號的方差最小化，提取區(qū)分度較高的腦電特征。在多分類問題中，可以使用一對多的思想將原始的二分類CSP 算法進(jìn)行擴(kuò)展，具體可以將多類信號中的一類作為單獨一類，其余類別統(tǒng)一作為另外一類處理，依次計算得到四個二類投影矩陣，構(gòu)造一個四類的CSP特征提取模型。

假設(shè)Xi是大小為N×T 的第i 類原始腦電信號，其中N 表示腦電信號的通道數(shù)，T 表示采樣點數(shù)，Xi的歸一化協(xié)方差矩陣為：

四分類任務(wù)中，混合空間協(xié)方差矩陣可以計算為：

其中：R1、R2、R3、R4分別表示對應(yīng)類別運動想象信號經(jīng)過多次實驗中得到的平均協(xié)方差矩陣。R經(jīng)過特征分解可以得到：

其中：U0表示R 的特征向量，V 表示R 的特征值矩陣。對特征值矩陣V進(jìn)行降序排列，并對U0進(jìn)行相應(yīng)調(diào)整，白化矩陣可計算為：

對于四分類任務(wù)，一對多CSP 算法在計算腦電信號的投影矩陣時，將其中一類信號作為基準(zhǔn)，其余的三類信號則歸為另一類，以R1基準(zhǔn)的計算為例，另一類R′1可以表示為：

分別計算R1和R′1的平均協(xié)方差矩陣，并表示為，然后對二者進(jìn)行白化變換可得：

分別對S1和S′1進(jìn)行特征分解：

其中：V1和V′1滿足V1+ V′1= I。投影矩陣可以計算為：

同理，Ri對應(yīng)的投影矩陣可以計算為：

將腦電信號Xi進(jìn)行投影，可以得到：

對各類別投影后的信號Zi(i = 1，2，3，4)各 取前m 行可以組成一個新的四類CSP 特征提取模型，共有N=4m行，特征向量的計算公式為：

2.3 人工神經(jīng)網(wǎng)絡(luò)

人工神經(jīng)網(wǎng)絡(luò)［18］由神經(jīng)元模型組成。神經(jīng)元模型是一種非線性的信息處理節(jié)點［19］，結(jié)構(gòu)如圖2所示。

圖2 神經(jīng)元模型結(jié)構(gòu)圖Fig.2 Structure diagram of neuron model

其中，x1，…，xn為神經(jīng)元i 的輸入，wi1，..，win為連接權(quán)重，yi為神經(jīng)元的輸出，f（·）為激活函數(shù)，則模型的輸出為：

神經(jīng)網(wǎng)絡(luò)有著多種訓(xùn)練方法，其中最為常用的是反向傳播算法（back propagation，BP）［20］。

圖3 為基于反向傳播的單隱層神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)，其中X1，…Xm為神經(jīng)網(wǎng)絡(luò)的輸入，y1，…，yL為神經(jīng)網(wǎng)絡(luò)的輸出。

基于BP的神經(jīng)網(wǎng)絡(luò)算法屬于監(jiān)督學(xué)習(xí)算法，學(xué)習(xí)過程分為正向傳播和反向傳播兩部分。在正向傳播過程中，神經(jīng)網(wǎng)絡(luò)輸入信號經(jīng)過網(wǎng)絡(luò)逐層傳輸?shù)捷敵鰧?，如果結(jié)果符合樣本的期望輸出，則算法結(jié)束；否則，需要計算網(wǎng)絡(luò)輸出與樣本的期望輸出之間的誤差，借助梯度下降法沿網(wǎng)絡(luò)的反方向?qū)W(wǎng)絡(luò)各層的權(quán)重和閾值進(jìn)行調(diào)整，達(dá)到最小化誤差的目的［21］。

圖3 BP神經(jīng)網(wǎng)絡(luò)圖Fig.3 BP neural network diagram

3 模型構(gòu)造方法

針對四分類任務(wù)，首先利用一對多CSP 算法構(gòu)建空間濾波器，使不同類別間的方差最大化，經(jīng)過投影可得到區(qū)分度較高的特征向量。由于運動想象腦電信號的主要信息位于μ節(jié)律（8～12Hz）和β波（14～30Hz），原始的腦電信號無法準(zhǔn)確提供對應(yīng)頻帶的信息，本文首先使用WPD 算法分解原始信號，得到多個細(xì)化頻帶信息，然后使用CSP 算法提取各個頻帶的腦電特征，最后使用人工神經(jīng)網(wǎng)絡(luò)進(jìn)行分類，算法流程如圖4所示。

具體步驟為：

圖4 算法流程圖Fig.4 Algorithm flow paradigm

步驟1：截取3～6s時間段的腦電信號，并將其濾波至8～30Hz（所選數(shù)據(jù)集的運動想象任務(wù)主要集中于3～6s，μ 節(jié)律和β 波頻率范圍主要在8～30Hz 頻帶內(nèi)）；

步驟2：使用WPD 對訓(xùn)練集樣本進(jìn)行分解與重構(gòu)，得到含有不同頻帶信息的信號樣本；

步驟3：對分解后的樣本使用一對多CSP 算法進(jìn)行特征提取，作為神經(jīng)網(wǎng)絡(luò)的輸入；

步驟4：將步驟3得到的特征向量輸入神經(jīng)網(wǎng)絡(luò)進(jìn)行訓(xùn)練，修正網(wǎng)絡(luò)的權(quán)重和閾值，完成訓(xùn)練；

步驟5：按照步驟2、3 對測試集樣本進(jìn)行處理，然后輸入訓(xùn)練好的神經(jīng)網(wǎng)絡(luò)得到最終的分類結(jié)果。

4 實驗結(jié)果及分析

4.1 實驗數(shù)據(jù)

本文的實驗數(shù)據(jù)選自“BCI IV”，數(shù)據(jù)集由Graz University of Technology 提供，共包含9 個志愿者的四分類運動想象腦電信號，四個任務(wù)分別對應(yīng)左手、右手、雙腳和舌頭。每個志愿者在不同的日期共完成兩輪實驗，每輪實驗采集288組實驗數(shù)據(jù)，由四類運動想象任務(wù)等分，每組實驗數(shù)據(jù)的采集過程按照圖5的實驗范式所示。

圖5 實驗范式流程圖Fig.5 Flow chart of experimental paradigm

用于腦電信號采集的實驗設(shè)備共有22個電極，擺放位置如圖6 所示。設(shè)備的采樣頻率為250Hz，使用0.5～100Hz的帶通濾波器對采集的原始信號進(jìn)行濾波。兩輪實驗數(shù)據(jù)分別作為訓(xùn)練樣本與測試樣本，前者用于模型訓(xùn)練，后者用于模型驗證。

4.2 參數(shù)選擇

本文利用WPD 和CSP 算法聯(lián)合提取腦電信號特征，需要確定CSP 每類模式的特征數(shù)m。在使用人工神經(jīng)網(wǎng)絡(luò)進(jìn)行分類的時候，由于神經(jīng)網(wǎng)絡(luò)的參數(shù)決定了網(wǎng)絡(luò)結(jié)構(gòu)的不同，需要確定隱藏層的神經(jīng)元個數(shù)、目標(biāo)函數(shù)的正則化系數(shù)等參數(shù)。不同參數(shù)的選擇會影響最終的實驗結(jié)果。

圖6 電極擺放位置圖(遵循10～20規(guī)則)Fig.6 Electrode position diagram（follow 10～20 rule）

運動想象腦電信號的主要信息位于μ 節(jié)律（4～12Hz）和β 波（14～30Hz）中，考慮到兩個節(jié)律的頻帶范圍和頻率分辨率，本文選取“db4”小波作為小波基函數(shù)，進(jìn)行5層小波包分解，選擇位于4～32Hz的8個相關(guān)子頻帶進(jìn)行后續(xù)的特征選擇。

CSP模式特征數(shù)m一般選擇在2～12之間，為確定最優(yōu)的特征數(shù)m，將神經(jīng)網(wǎng)路結(jié)構(gòu)暫定為雙隱層，各隱層的神經(jīng)元個數(shù)暫定為100，正則化系數(shù)暫定為0.1，調(diào)節(jié)m，則Kappa值變化如圖7所示。

圖7 不同m值的Kappa值Fig.7 Kappa rates with various value of m

當(dāng)CSP模式的特征數(shù)m為6時Kappa值最優(yōu)，固定該參數(shù)，然后對神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)進(jìn)行調(diào)整，需要的調(diào)整包括隱層數(shù)量、隱層的神經(jīng)元數(shù)量、正則化系數(shù)以及dropout層的使用。

圖8 不同隱層數(shù)量的Kappa值Fig.8 Kappa rates with various number of hidden layers

暫定隱層神經(jīng)元數(shù)量為100，正則化系數(shù)為0.1，無dropout 層，實驗結(jié)果如圖8所示，可以看出，隱層數(shù)目為1時，具有最優(yōu)的Kappa值0.72。

將隱層數(shù)目固定為1，隱層的神經(jīng)元數(shù)量為100，正則化系數(shù)為0.1，設(shè)置dropout 層的參數(shù)調(diào)節(jié)范圍為0～0.5，實驗結(jié)果如圖9所示。

圖9 不同dropout參數(shù)的分類準(zhǔn)確率Fig.9 Kappa rates with various parameter of dropout

當(dāng)dropout參數(shù)為0.4時，最優(yōu)Kappa值為0.73，固定該參數(shù)，暫定正則化系數(shù)為0.1，隱層神經(jīng)元數(shù)目的調(diào)節(jié)范圍設(shè)置為50～250，實驗結(jié)果如圖10所示。

隱層的神經(jīng)元數(shù)目為100時，有著最優(yōu)的Kappa值0.73。固定其他參數(shù)，調(diào)節(jié)正則化系數(shù)，當(dāng)正則化系數(shù)為0.1時，Kappa 值為0.73；當(dāng)正則化系數(shù)為0.01時，Kappa值為0.70。

4.3 實驗結(jié)果

圖10 不同隱層神經(jīng)元數(shù)目的Kappa值Fig.10 Kappa rates with various number of neurons in hidden layer

選擇“db4”小波基函數(shù)進(jìn)行5層小波包分解，設(shè)置CSP 模式的特征數(shù)為6，選用單隱層神經(jīng)網(wǎng)絡(luò)作為分類器，隱層神經(jīng)元數(shù)目、正則化系數(shù)和dropout參數(shù)分別選擇100、0.1、0.4，在基準(zhǔn)數(shù)據(jù)集中的實驗結(jié)果如表1所示。

表1 實驗結(jié)果對比Table 1 Comparison of experimental results

可以看出，本文提出的算法與競賽前3 名的算法準(zhǔn)確率相比在整體上有著較為明顯的提升，證實了本文算法的有效性。

5 結(jié) 論

針對運動想象腦電信號在多分類問題方面特征區(qū)分度不高、識別精度低的問題，本文使用了小波包變換對原始腦電信號進(jìn)行了多層分解，利用一對多CSP 算法對分解后的子頻帶進(jìn)行特征提取，利用單隱層神經(jīng)網(wǎng)絡(luò)進(jìn)行分類。最終的實驗結(jié)果表明，選擇合理的小波包分解層數(shù)、CSP的m參數(shù)和神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)，可以明顯提高分類的準(zhǔn)確率，與BCI IV競賽的基準(zhǔn)數(shù)據(jù)集前幾名結(jié)果的對比中，展現(xiàn)了本文方法的優(yōu)越性。本文的算法對于實現(xiàn)基于運動想象腦電信號的腦機(jī)接口控制具有重要的意義，未來將對在線系統(tǒng)的算法部署進(jìn)行更為深入的研究。