基于時(shí)延空時(shí)濾波的P300波形提取及目標(biāo)分類算法

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

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

腦-機(jī)接口系統(tǒng)(brain-computer interface,BCI)[1]是一種新的人機(jī)交互方式，可以把系統(tǒng)采集到的腦電生理信號(hào)，通過(guò)特征提取和分類等手段進(jìn)行處理后轉(zhuǎn)為控制信號(hào)，從而實(shí)現(xiàn)與外部設(shè)備的交互.由于不依賴人的肢體動(dòng)作，早期腦-機(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)識(shí)別[7]等.

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

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

1 算法原理

1.1 腦電信號(hào)模型

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

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

(1)

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

1.2 算法流程圖

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

圖1 算法流程

1.3 算法原理

1.3.1時(shí)域延時(shí)

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

(2)

(3)

1.3.2空時(shí)濾波

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

(4)

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

1.3.3交替優(yōu)化

(5)

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

(6)

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

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

1.3.4P300波形提取

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

1.3.5P300目標(biāo)分類

利用空時(shí)濾波器W和空域模式A對(duì)每個(gè)腦電信號(hào)樣本提取P300源成分作為分類特征：

(9)

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

2 實(shí)驗(yàn)與結(jié)果

2.1 仿真實(shí)驗(yàn)

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

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

圖2 3個(gè)仿真源波形

2.1.2仿真結(jié)果

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

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

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

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

2.2.1實(shí)驗(yàn)過(guò)程及采集

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

2.2.2結(jié)果

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

圖4 兩種算法P300提取效果對(duì)比

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

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

圖5 兩種算法的各項(xiàng)分類性能對(duì)比

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

3 結(jié)束語(yǔ)

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