林艷飛， 盧志強(qiáng)，3， 李博聞， 劉志文， 高小榕

(1.北京理工大學(xué) 信息與電子學(xué)院，北京 100081；2.清華大學(xué) 醫(yī)學(xué)院，北京 100084；3.中國船舶工業(yè)系統(tǒng)工程研究院，北京 100094)

腦-機(jī)接口系統(tǒng)(brain-computer interface,BCI)[1]是一種新的人機(jī)交互方式，可以把系統(tǒng)采集到的腦電生理信號，通過特征提取和分類等手段進(jìn)行處理后轉(zhuǎn)為控制信號，從而實現(xiàn)與外部設(shè)備的交互.由于不依賴人的肢體動作，早期腦-機(jī)接口系統(tǒng)可以適用于肌萎縮側(cè)索硬化癥，腦卒中等疾病患者，提升他們的生活質(zhì)量[2].目前腦-機(jī)接口系統(tǒng)也在應(yīng)用于健康人群提高人類認(rèn)知、行為和決策能力[3].腦電中包含多種成分，其中由oddball范式誘發(fā)，在刺激后 300 ms 出現(xiàn)的正波成分P300，是一種重要的腦電成分，是腦-機(jī)接口系統(tǒng)廣泛使用的腦電成分.基于P300電位的腦-機(jī)接口系統(tǒng)有腦電打字系統(tǒng)[4]、腦電開關(guān)[5]、圖片檢索[6]、目標(biāo)識別[7]等.

由于P300的信噪比很低，在傳統(tǒng)的提取方法中，多采用大量試次疊加平均的方法提高信噪比.然而在BCI系統(tǒng)中，要求從少試次甚至單試次的腦電信號中提取P300成分.隨著高密度腦電圖記錄技術(shù)的進(jìn)步，出現(xiàn)了許多利用多通道腦電信號的空域信息實現(xiàn)P300少試次提取的方法，如獨立成分分析(independent component analysis,ICA)[8]，稀疏成分分析(sparse component analysis,SCA)[9]，主成分分析(principal component analysis,PCA)[10]等.這一類方法可以在一定程度上提高腦電信號空域信息的利用率，但并不是專門針對ERP提取提出的，因而這類方法不能獲得最優(yōu)的性能.在過去的研究中，出現(xiàn)了一些專門針對ERP提取的算法，比如2006年德國FIRST Fraunhofer Institute的Lemm團(tuán)隊提出的正則化2階盲分離算法(regularized second-order blind identification,r-SOBI)[11],通過提取腦電信號的鎖相成分，提高信噪比，但該方法不能用于單試次的P300提取.2009年，美國佛羅里達(dá)大學(xué)的Li等[12]提出了一種基于模板的空時濾波方法，在空域進(jìn)行噪聲消除，從而提高信噪比.然而為了對提取出的ERP源進(jìn)行約束，該方法需要預(yù)先定義ERP波形的模板，這些模板在多數(shù)情況下只能以特定方式得到，且得到方式常常是隨機(jī)的，因此該方法在應(yīng)用中受到一定限制.此外，為了尋求空域的最大信噪比，2011年清華大學(xué)的Wu等[13]提出了ERP信噪比最大化算法(signal-to-noise ratio maximum,SIM)，該方法從多通道腦電數(shù)據(jù)中學(xué)習(xí)ERP的空時模型結(jié)構(gòu)，可以實現(xiàn)空域的高信噪比，但沒有充分利用腦電信號的時域信息.2016年，法國格勒諾布爾阿爾卑斯大學(xué)的Congedo團(tuán)隊提出了時域上的自適應(yīng)時空共模式(adaptive common spatio-temporal pattern,aCSTP)[14]，利用多變量空時濾波自動估計子空間維數(shù)，提升了信噪比，但該算法不能保證子空間維數(shù)的最優(yōu)性，且無法實現(xiàn)源分離.2017年，清華大學(xué)的Wu等[15]提出了ERP空譜模式分析(ERP spatial and spectral patterns,ESSPs)，基于循環(huán)下降的最大后驗估計算法來估計空譜模式，可以很好地估計ERP，但其性能依賴于參數(shù)的選擇，需要大量的訓(xùn)練數(shù)據(jù)進(jìn)行交叉驗證.

本文提出了一種基于空時濾波的P300提取及目標(biāo)分類算法(時延空時濾波算法)，對多通道腦電數(shù)據(jù)首先進(jìn)行時延，然后以空時濾波后的信號與期望信號誤差最小為目標(biāo)，構(gòu)造代價函數(shù)，通過交替優(yōu)化的方式估計空時濾波器和源信號，最終實現(xiàn)空域的源分離和時域的波形提取.通過從仿真腦電數(shù)據(jù)提取P300波形和對真實腦電數(shù)據(jù)的分類，并與目前應(yīng)用較為廣泛的SIM作對比，驗證了時延空時濾波算法的有效性.

1 算法原理

1.1 腦電信號模型

ERP信號具有相位鎖定的特點，自發(fā)腦電信號通常被認(rèn)為是零均值的高斯色噪聲，測量噪聲通常被認(rèn)為是零均值的高斯白噪聲[13]，觀測信號可以被建模為

Xk=AZ+Nk,k=1,2,…,K

(1)

式中：Xk∈RC×N為腦電觀測信號，k為樣本觀測序號；C為信號導(dǎo)聯(lián)數(shù)，N為信號長度;A∈RC×M為空域濾波器；M為源數(shù)量；Z∈RM×N為ERP源信號；Nk∈RC×N為自發(fā)腦電和測量噪聲混合產(chǎn)生的噪聲信號，各試次觀測噪聲彼此獨立，方差相同[16].

1.2 算法流程圖

本文提出的P300波形提取及目標(biāo)分類算法流程如圖1所示.首先將多通道腦電信號進(jìn)行時域延時，然后構(gòu)造最小二乘代價函數(shù)，循環(huán)迭代使代價函數(shù)收斂，估計得到空時濾波器和P300信號，對得到的P300信號做奇異值分解可得到空域模式和時域波形.將空時濾波器和空域模式帶入訓(xùn)練集可以訓(xùn)練得到Fisher分類器，最后在測試集中進(jìn)行線性判別目標(biāo)分類.

圖1 算法流程

1.3 算法原理

1.3.1時域延時

基于1.1節(jié)的腦電信號模型(1)，對Xk的第c導(dǎo)聯(lián)有

(2)

(3)

1.3.2空時濾波

時域濾波后，對式(3)中得到的結(jié)果進(jìn)行空域濾波，假設(shè)存在空間濾波器B，使得

(4)

為便于分析記：W=BS，W即為最終所求的空時濾波器.后面對AZ做奇異值分解，可以得到提取的多個源信號Z和對應(yīng)空域模式A，對W做分解，可以進(jìn)一步求得時域濾波器和空域濾波器[17].

1.3.3交替優(yōu)化

(5)

最小化式(5)，即最小化：

(6)

對于式(6)，有多個未知變量，無法直接求解，因此通過交替求解的方式對W和AZ分別進(jìn)行優(yōu)化.求解其中任一變量時，把其他變量都當(dāng)作已知量處理，得到迭代公式為

首先初始化AZ，把原始樣本的疊加平均作為AZ的初始值帶入式(7)，從而解出W，接著再把由式(7)求得的W帶入式(8)更新AZ，這樣便完成了一次循環(huán)迭代.通過式(5)設(shè)置閾值(如10-8)，當(dāng)結(jié)果滿足條件時循環(huán)結(jié)束.

1.3.4P300波形提取

對AZ做奇異值分解，得到提取的多個源信號Z和對應(yīng)空域模式A.根據(jù)各源成分功率大小進(jìn)行排序，選取恰當(dāng)?shù)男盘栐春蛯?yīng)的空域模式.并將源成分反投至Cz電極，得到P300波形.

1.3.5P300目標(biāo)分類

利用空時濾波器W和空域模式A對每個腦電信號樣本提取P300源成分作為分類特征：

(9)

然后利用訓(xùn)練集得到的P300源成分訓(xùn)練Fisher線性判別分析分類器，得到分類投影矢量，將P300特征矢量變?yōu)樽罱K的決策值，實現(xiàn)腦電信號的分類.

2 實驗與結(jié)果

2.1 仿真實驗

2.1.1數(shù)據(jù)描述

本文在P300信噪比為-20 dB，-15 dB，-10 dB，-5 dB，0 dB，5 dB的條件下生成腦電仿真數(shù)據(jù)，并提取P300波形.各信噪比下分別進(jìn)行50次蒙特卡洛實驗，在每次實驗中，隨機(jī)生成40試次的腦電信號，每個試次的信號包含30導(dǎo)腦電數(shù)據(jù)，每導(dǎo)包括250個數(shù)據(jù)點，信號采樣率為250 Hz，信號由3個ERP成分和生成的背景噪聲疊加得到.3個ERP成分由Gamma函數(shù)產(chǎn)生，在波形和出現(xiàn)時間上類似于典型的P100，P200和P300成分，并通過Gram-Schmidt過程正交化，如圖2所示.背景噪聲由20個背景自發(fā)腦電線性混疊和測量噪聲組合而成.每個背景自發(fā)腦電的時間過程均為1/f噪聲，測量噪聲為高斯白噪聲，其方差為腦電信號方差的1/6[13].ERP成分的混疊矩陣和背景腦電的混疊矩陣均為隨機(jī)生成，矩陣中的所有元素都服從[0，1]均勻分布.

圖2 3個仿真源波形

2.1.2仿真結(jié)果

在仿真實驗中，已知真實波形，通過比較提取波形與真實波形的相關(guān)系數(shù)，并與近年來應(yīng)用較為廣泛的SIM作對比，驗證本文算法的有效性.SIM基于空域濾波實現(xiàn)信噪比最大化，衍生了許多優(yōu)秀算法，如提取失匹配負(fù)波(mismatch negativity，MMN)的重采樣差分(resampling difference and SIM，rdSIM)[18]，估計穩(wěn)態(tài)誘發(fā)電位的ESSP等，通過與空域濾波算法SIM的比較，可以驗證本算法采用空時濾波是否具有優(yōu)勢，因此本文僅與SIM作對比.不同信噪比下的相關(guān)系數(shù)如圖3所示.可見，不同信噪比下時延空時濾波算法的相關(guān)系數(shù)要高于SIM.即便在信噪比很低(小于-10 dB)的情況下，時延空時濾波算法相關(guān)系數(shù)也大于0.7，能很好得恢復(fù)真實ERP波形，且恢復(fù)效果隨信噪比的提升有明顯的增強(qiáng)，在信噪比大于-10 dB時，相關(guān)系數(shù)已經(jīng)非常接近1.

圖3 兩種算法在不同信噪比下的相關(guān)系數(shù)

再對時延空時濾波算法和SIM的相關(guān)系數(shù)做配對T檢驗，結(jié)果得出，在不同信噪比下，時延空時濾波算法的相關(guān)系數(shù)均顯著高于SIM(p<0.05)，說明時延空時濾波算法能更好地恢復(fù)真實ERP波形.

2.2 真實EEG數(shù)據(jù)

2.2.1實驗過程及采集

數(shù)據(jù)來源為第二屆BCI競賽的數(shù)據(jù)集IIB[19]，該數(shù)據(jù)由Wadsworth BCI2000軟件記錄[20].在實驗中，給受試呈現(xiàn)一個6×6的字符矩陣.矩陣的每個單元格都包含一個字符，受試的任務(wù)是將注意力集中在研究人員指定單詞的字符上.該矩陣的6行和6列以5.7 Hz的頻率依次隨機(jī)凸顯.在每輪閃爍中有一個特定行和一個特定列凸顯了所需的字符.與不包含所需字符刺激引起的反應(yīng)不同，實驗中包含所需字符的不常見刺激引起了P300反應(yīng).數(shù)據(jù)使用該數(shù)據(jù)集中一名受試者兩次P300誘發(fā)電位的記錄，分為11組，每組包含90個目標(biāo)試次和450個非目標(biāo)試次.

2.2.2結(jié)果

實驗數(shù)據(jù)為64導(dǎo)信號，采樣率240 Hz，通過0.2～28 Hz的帶通濾波器進(jìn)行濾波，同時進(jìn)行基線校正，即減去刺激呈現(xiàn)前200 ms的平均電壓.截取每次刺激開始后0～500 ms的數(shù)據(jù)作為樣本信號，用時延空時濾波算法和SIM分別對每組數(shù)據(jù)進(jìn)行P300波形提取，圖4(a)～(h)分別展示了時延空時濾波算法和SIM對其中4組數(shù)據(jù)(第1、2、9、10組)的提取結(jié)果.第一行為功率較大的3個ERP源的空域模式，第二行為空域模式對應(yīng)的ERP成分的時域波形，方框為挑選出來的P300空域模式及成分.第三行左圖虛線為挑選出的P300源反投到Cz電極得到的P300波形，左圖實線為該組所有目標(biāo)試次疊加平均的波形，右圖虛線為測試集中各算法處理后得到的單試次波形，右圖黑實線為所有試次的均值.為了更好地展示兩種算法的提取效果，分別對空域模式，時域波形和測試集中的單試次波形幅值做歸一化處理.由圖4第1行和第2行子圖及所選源成分知，時延空時濾波算法提取的P300源在4組數(shù)據(jù)中均排在所有源的第1位，而SIM提取的P300源在第2組和第9組中排在第1位，在第1組和第10組中排在第2位；由圖4第3行左子圖所示，兩種算法均能得到較好的反投后的P300波形；由圖4第3行右子圖所示，時延空時濾波算法得到的各單試次波形在4組數(shù)據(jù)中更加接近目標(biāo)波形，而SIM得到的各單試次波形更接近非目標(biāo)波形.

圖4 兩種算法P300提取效果對比

本文用所有目標(biāo)試次的疊加平均波形作為參考，直觀展示兩種算法的提取效果.結(jié)果表明，時延空時濾波算法的提取結(jié)果在反投后的波形、相位上更加接近疊加平均的結(jié)果.而且，時延空時濾波算法得到的各單試次波形更加接近目標(biāo)波形，而SIM得到的各單試次波形更接近非目標(biāo)波形，因此時延空時濾波算法提取的波形更利于分類.此外，在11組數(shù)據(jù)中，時延空時濾波算法提取的P300源均排在所有源的第一位，源排序準(zhǔn)確率為100%，具有較好的源挑選效果，而SIM的源排序準(zhǔn)確率僅為81.8%，因此在源排序方面時延空時濾波算法要優(yōu)于SIM.

為了更加準(zhǔn)確地評價兩種算法對于P300目標(biāo)分類的效果，采用Fisher線性判別分析進(jìn)行目標(biāo)分類.將每組數(shù)據(jù)分為訓(xùn)練集(45目標(biāo)/45非目標(biāo))和測試集(45目標(biāo)/405非目標(biāo)).得到如圖5(a)～(d)所示的分類結(jié)果，在每個子圖中，柱狀圖分別表示兩種算法的分類準(zhǔn)確率、AUC值、擊中率和虛警率，實線表示兩種算法在11組數(shù)據(jù)下的均值.

圖5 兩種算法的各項分類性能對比

分別對兩種算法分類準(zhǔn)確率、AUC值、擊中率和虛警率做配對T檢驗，可以得到時延空時濾波算法的分類準(zhǔn)確率顯著高于SIM(p=0.002<0.05)，時延空時濾波算法的AUC值顯著高于SIM(p=0.001<0.05)，時延空時濾波算法的擊中率顯著高于SIM(p=0.026<0.05)，時延空時濾波算法的虛警率顯著低于SIM(p= 0.048<0.05).因此，時延空時濾波算法的性能顯著優(yōu)于SIM.

3 結(jié)束語

本文提出了一種基于時延空時濾波的P300波形提取及目標(biāo)分類算法.通過時域延時后約束觀測信號樣本與期望信號的差值，得到空時濾波的代價函數(shù)，經(jīng)過交替優(yōu)化使代價函數(shù)收斂后得到空時濾波器和P300信號，最后通過奇異值分解可提取P300成分的空域模式和時域波形.仿真實驗和真實P300數(shù)據(jù)實驗中，時延空時濾波算法均表現(xiàn)出良好的性能.因此，該算法可有效提取P300波形并進(jìn)行目標(biāo)分類.