改進(jìn)的Wigner-ville分布提取P300視覺誘發(fā)電位

程萬里，任兆麟，喬曉艷

（山西大學(xué) 物理電子工程學(xué)院，山西 太原 030006）

腦-機(jī)接口（BCI）是一種不依賴于常規(guī)腦外周神經(jīng)和肌肉系統(tǒng)直接在人腦和計(jì)算機(jī)等電子設(shè)備之間進(jìn)行信息交流和控制的通訊系統(tǒng)［1］.腦電信號(hào)特征提取是BCI技術(shù)的關(guān)鍵環(huán)節(jié)［2-3］.目前特征提取的方法主要有時(shí)域、頻域和時(shí)頻分析方法［3］.時(shí)域分析方法主要有相干平均法，頻域分析方法的主要有功率譜估計(jì)、AR參數(shù)模型等［4］.功率譜估計(jì)是一種能量密度的頻域分析，它并不告訴我們某個(gè)頻率分量出現(xiàn)的具體時(shí)刻及變化趨向，因此不能用來分析P300電位.AR參數(shù)模型更適合分析平穩(wěn)隨機(jī)信號(hào)，而EEG信號(hào)是高度非平穩(wěn)的信號(hào).

時(shí)頻分析方法主要有小波變換法、Wigner-ville分布法等.小波變換法是一種非線性變換，計(jì)算速度快，適合于在線分析［5-6］.但是由于EEG信號(hào)的產(chǎn)生機(jī)理復(fù)雜，很難獲得準(zhǔn)確的先驗(yàn)知識(shí)，在經(jīng)驗(yàn)之外的其他頻段也包含著有用信息，該方法有可能遺漏某些重要的信息.腦電信號(hào)是非平穩(wěn)時(shí)變信號(hào)，單純的時(shí)域或頻域方法都不能很好地反映腦電信號(hào)的特征，Wigner-ville分布法是一種最基本、也是應(yīng)用最多的時(shí)頻分布，具有時(shí)移不變性、頻移不變性、時(shí)頻集聚性等特征.因此，通過對(duì)信號(hào)進(jìn)行Wigner-ville分析，不但可以獲得信號(hào)的時(shí)間、頻率分布圖，還可以求出信號(hào)的頻率變化情況，從而可以更好地對(duì)腦電信號(hào)進(jìn)行分類、判別［8］.因此，本文采用Wigner-ville分布方法對(duì)BCIⅢ競(jìng)賽的視覺誘發(fā)腦電數(shù)據(jù)進(jìn)行特征提取.

P300是事件相關(guān)電位的一種，其峰值約出現(xiàn)在相關(guān)事件發(fā)生后300ms，相關(guān)事件出現(xiàn)的概率越小，所引起的P300越明顯，并在頭部頂骨區(qū)最為顯著［7-8］.P300電位的發(fā)生不需要經(jīng)過特殊的訓(xùn)練，通常受試者只需經(jīng)過簡(jiǎn)單操作后，就可以在腦機(jī)接口系統(tǒng)上取得滿意的效果［8］；因此，本文選擇了P300作為視覺誘發(fā)腦-機(jī)接口系統(tǒng)的特征信號(hào).

2 算法基本原理

2.1 相干平均算法原理［9，11］

利用相干平均法提取誘發(fā)腦電基于如下三個(gè)假設(shè)：

（1）誘發(fā)電位與噪聲是加性關(guān)系，而且兩者相互獨(dú)立；

（2）每次刺激產(chǎn)生的誘發(fā)腦電波形不變；

（3）噪聲是平穩(wěn)、零均值的隨機(jī)信號(hào)，各次記錄的噪聲之間互不相關(guān).

設(shè)記錄到的信號(hào)表示為：

式（1）中，xi（n）表示第i次記錄到測(cè)量信號(hào)，ni（n）表示第i次記錄到的噪聲信號(hào)i＝1，2，…，N，變量n表示記錄中第n個(gè)采樣值n＝1，2，…M.

設(shè)信號(hào)的功率為P，噪聲均值為零，噪聲的方差為σ2，此時(shí)信噪比為SNR＝.從結(jié)果可以看出，經(jīng)過N次相干平均后，信號(hào)的功率仍然為P，噪聲均值為零，而方差變?yōu)?，因此相干平均后的信噪比為?/p>

可見，信噪比提高了 N 倍［9-10］.

2.2 Wigner-ville分布［10］和交叉項(xiàng)分析［11］

Wigner-ville分布是一種時(shí)頻混合的信號(hào)表示法，能同時(shí)進(jìn)行時(shí)域和頻域分析，并把兩者結(jié)合起來，其各階距具有明確的物理意義，信號(hào)s（t）的Wigner-ville分布的定義式如下：

Wigner-ville分布也可以用解析信號(hào)的頻譜表示成：

上式中z（t）是s（t）的解析信號(hào)，也就是其經(jīng)過Hilbert變換后得到的信號(hào)，即：

式（6）中t和τ是實(shí)的變量，P.V.表示取積分的主值.

由于Wigner-ville分布是雙線性時(shí)頻分析方法，當(dāng)信號(hào)包含多個(gè)頻率成份時(shí)，信號(hào)的 Wigner-ville分布會(huì)出現(xiàn)嚴(yán)重的干擾項(xiàng)，即：

由（8）式可以看出，信號(hào)包含的頻率成分越多，WVD分析后產(chǎn)生的交叉項(xiàng)越嚴(yán)重.在實(shí)際應(yīng)用中，因?yàn)榻徊骓?xiàng)的存在，使得真實(shí)的頻率成分淹沒在交叉項(xiàng)成分中無法識(shí)別.

2.3 算法改進(jìn)和實(shí)現(xiàn)

式中G（s-t）是時(shí)間平滑函數(shù)g（s-t）的解析函數(shù).由于加窗的作用在很大程度上抑制了交叉項(xiàng)干擾，從而保證了高精度的時(shí)頻分辨率［12］.

算法具體實(shí)現(xiàn)步驟描述如下：

（1）對(duì)預(yù)處理后的數(shù)據(jù)進(jìn)行共平均參考；

（2）計(jì)算共平均參考后的矩陣X的協(xié)方差矩陣Cx；

（3）計(jì)算白化變換矩陣Z

其中，E是以Cx的單位特征向量為列的矩陣，D＝diag（d1…dn）是以Cx的特征值為對(duì)角元素的對(duì)角矩陣，其推導(dǎo)過程可參見文獻(xiàn)［13］；

（4）根據(jù)公式（6）對(duì)白化后的數(shù)據(jù)進(jìn)行Hilbert變換；

（5）對(duì)某一導(dǎo)聯(lián)的數(shù)據(jù)，選擇靶刺激，分別采用傳統(tǒng)的 Wigner-ville分布和改進(jìn)的 Wigner-ville分布計(jì)算.

3 實(shí)驗(yàn)范式［13］

采用視覺誘發(fā)Oddball實(shí)驗(yàn)范式，數(shù)據(jù)集由BCI2000系統(tǒng)采集，由Donchin在2000年記錄的P300視覺誘發(fā)電位.實(shí)驗(yàn)中，受試者需注視36個(gè)不同字符中的某個(gè)字符（如圖1），實(shí)驗(yàn)?zāi)康氖峭ㄟ^分析視覺誘發(fā)腦電特征，正確識(shí)別受試者正在注視的字符.

圖1 實(shí)驗(yàn)范式用戶顯示界面Fig.1 This figure illustrates the user display for this paradigm

實(shí)驗(yàn)中，字符矩陣以2.5s為一個(gè)周期，每個(gè)字符有相同的加亮和變黑的時(shí)間，矩陣中的每個(gè)行和列被隨機(jī)地加亮100ms，每次加亮后，矩陣變黑75ms，一個(gè)周期中矩陣共進(jìn)行12次行、列加亮，其中2次加亮包含目標(biāo)字符，即為靶刺激，其余10次為非靶刺激.每個(gè)目標(biāo)字符重復(fù)15次實(shí)驗(yàn)，因此，每個(gè)字符共有180次加亮.數(shù)據(jù)由64導(dǎo)聯(lián)的腦電放大器采集，并通過0.1～60Hz的帶通濾波和240Hz頻率采樣后，被分為2個(gè)文件，以Matlab數(shù)據(jù)格式＊.mat提供，其中包括一個(gè)包含85個(gè)目標(biāo)字符的訓(xùn)練文件和包含100個(gè)目標(biāo)字符的測(cè)試文件，有兩位實(shí)驗(yàn)者A和B參與了實(shí)驗(yàn)，所有數(shù)據(jù)文件以單精度存儲(chǔ)，需要通過Matlab轉(zhuǎn)化為雙精度.

4 視覺誘發(fā)腦電特征提取

視覺誘發(fā)腦電［14］非常微弱，幅度很小，幅值約為0.5～50μV，屬于非平穩(wěn)和非高斯的微弱信號(hào)，而且通常被淹沒在很強(qiáng)的背景噪聲中（背景噪聲包括自發(fā)腦電——峰-峰值為50～100μV以及眼電、肌電等復(fù)雜成分），從強(qiáng)背景噪聲中提取誘發(fā)電位或者說提高信噪比的基本途徑主要有兩條：平均和濾波.本文采用導(dǎo)聯(lián)選擇、共平均參考、帶通濾波、數(shù)據(jù)分割、降采樣等預(yù)處理方法［14］.在以上預(yù)處理的基礎(chǔ)上，采用2.1節(jié)方法對(duì)數(shù)據(jù)進(jìn)行15次相干平均；最后，利用2.2節(jié)Wigner-ville分布方法對(duì)相干平均后的誘發(fā)腦電信號(hào)進(jìn)行特征提取，獲得P300電位.

5 結(jié)果與分析

圖2 相干平均法獲得的P300電位Fig.2 P300potential obtained by coherent average method

圖3 Wigner-ville分布提取的P300電位Fig.3 P300potential obtained by Wigner-Ville distribution

圖4 改進(jìn)的Wigner-ville分布提取的P300電位Fig.4 P300potential obtained by improved Wigner-Ville distribution

實(shí)驗(yàn)數(shù)據(jù)是由64導(dǎo)聯(lián)的腦電放大器記錄得到，本文僅使用Fz，F(xiàn)c1，F(xiàn)c2，C3，Cz，C4，Pz，PO7，PO8，F(xiàn)P1共10個(gè)導(dǎo)聯(lián)，因?yàn)樵谶@些導(dǎo)聯(lián)上視覺誘發(fā)特征較為明顯.首先，對(duì)原始腦電數(shù)據(jù)按照上述的預(yù)處理方法進(jìn)行處理；然后，對(duì)Cz導(dǎo)聯(lián)上的第一個(gè)目標(biāo)字符進(jìn)行處理，選擇每次刺激起始點(diǎn)后1s內(nèi)的數(shù)據(jù)作為處理對(duì)象，圖2（P227）為15次相干平均后輸出的結(jié)果，圖3（P227）為相干平均后，采用Wigner-ville分布提取到的靶刺激P300電位.圖4為相干平均后，采用改進(jìn)的Wigner-ville分布提取到的靶刺激P300電位.對(duì)腦電信號(hào)運(yùn)用Wigner-ville分布進(jìn)行時(shí)頻分析，如圖3所示.由于P300成分的頻段主要集中在2-15Hz范圍內(nèi)，由圖3（a）可知，在9Hz處局部能量聚集，且在266ms處出現(xiàn)一個(gè)負(fù)峰值，峰值大小為17.95μV，但是交叉項(xiàng)干擾嚴(yán)重，而非靶刺激在266ms處沒有明顯的峰值，對(duì)比15次相干平均后的腦電數(shù)據(jù)（如圖2所示），靶刺激在其266ms處有一個(gè)負(fù)峰值，峰值為1.218μV，Wigner-ville分布提取到的P300電位的幅度約為相干平均提取到的P300電位的幅度的15倍，而非靶刺激在266ms處沒有明顯的負(fù)峰.對(duì)腦電信號(hào)運(yùn)用改進(jìn)的Wigner-ville分布進(jìn)行時(shí)頻分析，由圖4可以明顯地看出，靶刺激在266ms有更明顯的峰值，峰值大小為11.33μV，約為相干平均提取到的P300電位的幅度的11倍，而且改進(jìn)的 Wigner-ville分布有效地抑制了交叉項(xiàng)的影響，并保持了高度集中的時(shí)頻分布.說明在經(jīng)過15次相干平均后使用 Wigner-ville分布能夠有效獲得視覺誘發(fā)腦電信號(hào)的P300電位，而使用改進(jìn)的Wigner-ville分布能夠有效地抑制交叉項(xiàng)干擾，更為有效地獲得P300電位.

6 結(jié)論

相干平均方法對(duì)于提取P300誘發(fā)電位具有簡(jiǎn)單易行的優(yōu)點(diǎn)，但其最大的缺陷在于需要大量疊加平均才能將P300電位提高到后續(xù)處理的要求.如果受試者使用BCI系統(tǒng)進(jìn)行長時(shí)間視覺刺激，容易產(chǎn)生視覺疲勞，使得提取P300電位產(chǎn)生誤差.本文把Wigner-ville分布引入到腦電信號(hào)處理中，通過對(duì)視覺誘發(fā)腦電信號(hào)中P300電位的提取，表明少次相干平均后，采用Wigner-ville分布提取P300電位是可行的和有效的，采用改進(jìn)的Wigner-ville分布可以有效地抑制交叉項(xiàng)干擾，更好地對(duì)腦電信號(hào)進(jìn)行識(shí)別分類，為建立在線BCI系統(tǒng)奠定了基礎(chǔ).

［1］Wolpaw J R，Birbaumer N，McFarland W J.Brain-Computer Interface for communication and control［J］.Clinical Neurophysiology，2002，113（6）：767-791.

［2］Lebedev M A，Nicolelis M A.Brain-machine Interfaces：Past，Present and Future［J］.Trends in Neurosciences，2006，29（9）：536-546.

［3］Salavris M，Sepulveda F.Wavelet and Ensemble of FLDs for P300Classification［C］／／Proceedings of the 4th International IEEE EMBS Conference on Neural Engineering，Antalya，Turkay，2009：339-342.

［4］徐寶國，宋愛國.基于小波變換和 AR參數(shù)模型的腦電信號(hào)識(shí)別方法［J］.數(shù)據(jù)采集與處理，2008，23（5）：580-583.

［5］王攀，沈繼忠，施錦河.基于小波變換和時(shí)域能量熵的P300特征提取算法［J］.儀器儀表學(xué)報(bào)，2011，32（6）：1285-1288.

［6］張方堃.基于視覺誘發(fā)電位的腦機(jī)接口技術(shù)研究［D］.北京：北京工業(yè)大學(xué)，2010：37-53.

［7］高均波.基于腦電的想象運(yùn)動(dòng)分類算法研究［D］.山東：山東大學(xué)，2008.

［8］張劍慧.基于P300的在線腦機(jī)接口系統(tǒng)［D］.浙江：浙江大學(xué)計(jì)算機(jī)學(xué)院，2010：14-16.

［9］王洪濤.視覺誘發(fā)電位腦機(jī)接口關(guān)鍵技術(shù)研究［J］.重慶文理學(xué)院學(xué)報(bào)：自然科學(xué)版，2010，29（1）：69-74.

［10］許佳.認(rèn)知腦電信號(hào)檢測(cè)分析軟件的研究［D］.重慶：重慶大學(xué)，2008：5-6.

［11］陸光華，彭學(xué)愚，張林讓，等.隨機(jī)信號(hào)處理［M］.西安：西安電子科技大學(xué)出版社，2002：93-121.

［12］楊海濤，朱仕軍，文中平，等.基于 Wigner-ville分布的譜分解效果分析［J］.勘探地球物理進(jìn)展，2009，32（1）：37-39.

［13］嚴(yán)娜.視聽誘發(fā)腦電特征研究［D］.太原：山西大學(xué)，2012：33-33.

［14］嚴(yán)娜，喬曉燕，李鵬.基于小波-雙譜分析的視覺誘發(fā)腦電特征提?。跩］.測(cè)試技術(shù)學(xué)報(bào)，2012，26（1）：20-25.