面向事件相關(guān)電位成分的稀疏字典構(gòu)建方法

李海峰 豐 上 馬 琳 徐睿峰 徐忠亮

(1. 哈爾濱工業(yè)大學(xué)計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，黑龍江哈爾濱 150001； 2. 哈爾濱工業(yè)大學(xué)(深圳)，廣東深圳 518055； 3. 東北林業(yè)大學(xué)信息與計(jì)算機(jī)工程學(xué)院，黑龍江哈爾濱 150040)

1 引言

事件相關(guān)電位(Event-related potential, ERP)指的是在特定的感覺(jué)、認(rèn)知或運(yùn)動(dòng)事件后，在頭皮觀測(cè)到的大腦電反應(yīng)[1]。一般而言，ERP可以指任何對(duì)應(yīng)某種刺激的固定的電生理反應(yīng)。ERP提供了一種評(píng)估大腦功能的非侵入手段。ERP可以通過(guò)采集腦電圖(electroencephalography，EEG)來(lái)獲得。EEG包含了近千個(gè)同時(shí)進(jìn)行的腦部活動(dòng)的電信號(hào)，其中跟刺激或事件有關(guān)的大腦活動(dòng)所占比例很小，因此為了能夠清晰的獲得大腦對(duì)刺激的響應(yīng)，需要對(duì)原始EEG信號(hào)做處理之后才能獲得ERP波形[2]。最常用的是疊加平均法，實(shí)驗(yàn)過(guò)程中同樣的刺激重復(fù)多次，把每次采集到的EEG波形進(jìn)行疊加平均，使得隨機(jī)的大腦活動(dòng)在信號(hào)中占比逐漸降低，最終獲得ERP的波形。

疊加平均法的前提假設(shè)是與事件相關(guān)的電位波動(dòng)在多次觀測(cè)中會(huì)保持一致，而噪聲也可以近似為零均值高斯隨機(jī)過(guò)程[2]。疊加之后，噪聲的能量逐漸降低，而ERP波形會(huì)保持一致。這樣做的問(wèn)題是，許多與認(rèn)知狀態(tài)有關(guān)的信號(hào)成分會(huì)在疊加過(guò)程中被消耗掉[3]。如果有一種能夠單次刺激之后獲得的EEG數(shù)據(jù)中提取ERP波形的方法，就可以避免有關(guān)的信號(hào)成分被損耗，同時(shí)通過(guò)分析ERP波形隨著時(shí)間的變化，可以分析認(rèn)知狀態(tài)隨著實(shí)驗(yàn)進(jìn)行是否有波動(dòng)[4]。

隨著信號(hào)處理技術(shù)的進(jìn)步，越來(lái)越多的研究開始將重點(diǎn)放在EEG信號(hào)或ERP波形的稀疏分解上[5- 6]。稀疏分解的方法非常適用于非穩(wěn)態(tài)、包含成分有限的信號(hào)，例如語(yǔ)音信號(hào)和EEG信號(hào)[7]。利用稀疏分解方法，每段EEG數(shù)據(jù)可以算作一幀，利用已知字典就可以求解每幀數(shù)據(jù)的稀疏表示系數(shù)[8]。獲得系數(shù)之后，既可以在便攜設(shè)備上實(shí)現(xiàn)EEG數(shù)據(jù)的壓縮傳輸[9]以及偽跡祛除[10]，也可以利用系數(shù)實(shí)現(xiàn)ERP波形的單次提取和識(shí)別[11]。

K-SVD字典學(xué)習(xí)算法[12]是目前最廣泛應(yīng)用的，能夠給出相對(duì)嚴(yán)格意義的稀疏解的稀疏分解算法。K-SVD是一種字典學(xué)習(xí)算法，用于通過(guò)奇異值分解方法創(chuàng)建用于稀疏表示的過(guò)完備字典。 K-SVD是k-means聚類方法的一種推廣，它通過(guò)在基于當(dāng)前字典的稀疏編碼輸入數(shù)據(jù)之間迭代交替工作，并且更新字典中的原子以更好地適合輸入數(shù)據(jù)。Zhang等[13]利用K-SVD算法獲得的字典，實(shí)現(xiàn)了車輛駕駛員的警覺(jué)度檢測(cè)；Nagaraj等[14]把患者的EEG信號(hào)使用K-SVD算法進(jìn)行字典求解和稀疏建模，實(shí)現(xiàn)了癲癇發(fā)作的實(shí)時(shí)檢測(cè)。

但是利用K-SVD進(jìn)行EEG信號(hào)分解時(shí)，由于初始字典原子常用穩(wěn)態(tài)信號(hào)，很難獲得能夠準(zhǔn)確描述非穩(wěn)態(tài)EEG信號(hào)的原子，EEG信號(hào)中的成分，尤其是ERP相關(guān)成分，極易混雜在大量字典原子中。為了解決上述問(wèn)題，實(shí)現(xiàn)將EEG信號(hào)分解為具有實(shí)際意義的，可以用于EEG信號(hào)按事件分幀之后的單次ERP成分提取的稀疏分析方法，本文提出了將稀疏性能評(píng)價(jià)指標(biāo)用于稀疏字典構(gòu)建過(guò)程中作為迭代訓(xùn)練過(guò)程中的約束條件，結(jié)合成熟的K-SVD算法，最終訓(xùn)練出包含完整ERP波形的稀疏字典。

本文在第2節(jié)介紹基于稀疏性能約束的稀疏字典訓(xùn)練方法，在第3節(jié)介紹基于這個(gè)方法對(duì)一個(gè)公開EEG數(shù)據(jù)集的ERP波形進(jìn)行稀疏分析的操作步驟，在第4節(jié)展示了利用這一方法獲得的結(jié)果。結(jié)果表明，本方法可以有效的將ERP成分提取到稀疏字典原子中。

2 基于稀疏性能約束的稀疏字典訓(xùn)練方法

2.1 稀疏性能指標(biāo)

信號(hào)的稀疏建模是通過(guò)將字典D=[d1|d2|…|dM]∈RN×M中的原子dk進(jìn)行線性組合,對(duì)信號(hào)y∈RN進(jìn)行表示的過(guò)程，如式(1)所示。

y=DA+e

(1)

(2)

(3)

字典學(xué)習(xí)的過(guò)程是對(duì)一個(gè)訓(xùn)練集Y={yi|i=1,2,…,p}，訓(xùn)練字典D，并通過(guò)該字典求解問(wèn)題P0或P(0,ε)，得到每一個(gè)yi對(duì)應(yīng)的稀疏模型Ai=(Ai1,Ai2,…,Ain),令線性表示yi=DAi+ei中的誤差ei最小。為了求解這一問(wèn)題，用l1范數(shù)代替l0范數(shù)作為約束條件的稀疏字典求解目標(biāo)函數(shù)可以表述為：

(4)

(5)

(6)

公式(6)中αij，αik表示系數(shù)向量αi的第j、k個(gè)元素。該指標(biāo)與重構(gòu)的相對(duì)誤差能量呈現(xiàn)顯著的正相關(guān)關(guān)系。

訓(xùn)練樣本的SPI指標(biāo)由該樣本的稀疏重構(gòu)系數(shù)矩陣計(jì)算得出，可以在任何現(xiàn)有字典訓(xùn)練算法基礎(chǔ)上獨(dú)立進(jìn)行，因此從理論上講，現(xiàn)有算法都可以通過(guò)將SPI作為約束條件，使之具備從根本上優(yōu)化重構(gòu)矩陣稀疏性能、并進(jìn)一步減小誤差的能力。

2.2 利用稀疏性能約束的K-SVD算法

K-SVD算法的求解模型基于l1范數(shù)，通過(guò)對(duì)重構(gòu)系數(shù)向量絕對(duì)值之和的限制來(lái)實(shí)現(xiàn)稀疏求解。K-SVD算法的約束條件如(7)所示。

(7)

通過(guò)這樣一個(gè)約束條件，把NP難的l0范數(shù)求解轉(zhuǎn)化為一個(gè)針對(duì)重構(gòu)誤差的凸優(yōu)化問(wèn)題。在建模初始階段，一部分樣本可能由較多的成分組成，導(dǎo)致字典初始化不夠充分，而導(dǎo)致重構(gòu)系數(shù)能量趨于分散時(shí)，l1范數(shù)難以使這些重構(gòu)系數(shù)能量趨于聚攏、實(shí)現(xiàn)稀疏，進(jìn)而在整個(gè)學(xué)習(xí)過(guò)程中，稀疏字典都無(wú)法有效的針對(duì)這些樣本進(jìn)行學(xué)習(xí)，不能提取這些樣本的真實(shí)成分，只能使用更多的現(xiàn)有成分進(jìn)行擬合；這直接導(dǎo)致了上文所提現(xiàn)象的產(chǎn)生，即算法傾向于給出由更多與成分不相等的原子，而非完整成分構(gòu)成的稀疏解，從事件相關(guān)電位分析的角度來(lái)看，無(wú)法提取出能夠更精準(zhǔn)描述ERP波形的字典原子，也就是說(shuō)EPR波形會(huì)分布在較多的字典原子中，這些樣本未能有效參與到字典的學(xué)習(xí)過(guò)程中，沒(méi)有做到將樣本幀分解為具有實(shí)際意義的信號(hào)源。

原始K-SVD 算法的約束條件為重構(gòu)信號(hào)與原始信號(hào)的誤差，我們將SPI加入到目標(biāo)函數(shù)中，如公式(8)，使SPI是否足夠小也作為一個(gè)字典訓(xùn)練階段的約束條件，當(dāng)訓(xùn)練進(jìn)行過(guò)程中，同時(shí)優(yōu)化，如果SPI達(dá)到了一個(gè)極小值，即迭代過(guò)程中變化可以忽略時(shí)，則停止訓(xùn)練。如果SPI已經(jīng)不發(fā)生明顯變化，根據(jù)SPI的定義，說(shuō)明樣本的完整成分已經(jīng)集中在幾個(gè)原子之中，對(duì)于ERP樣本來(lái)說(shuō)，代表ERP波形的主要成分已經(jīng)包含在字典原子波形中。

(8)

SPI作為約束條件的物理意義是指，在代表特定成分語(yǔ)境下，能量在字典原子中的集中程度。與LASSO不同之處在于，SPI考慮到了同樣L1范數(shù)下，稀疏分布的不同情況。加入了SPI的稀疏約束項(xiàng)，比傳統(tǒng)LASSO更精密。

3 方法驗(yàn)證

3.1 數(shù)據(jù)集選取

本文選取了一個(gè)公開的EEG數(shù)據(jù)集。該數(shù)據(jù)集記錄的是基于聽覺(jué)P300的拼寫器的測(cè)試數(shù)據(jù)。該數(shù)據(jù)集由H?hne等人[16]于2010年公開。實(shí)驗(yàn)數(shù)據(jù)由若干個(gè)trial組成，每個(gè)trial包括9個(gè)聲音刺激，每個(gè)刺激之間間隔120 ms。實(shí)驗(yàn)數(shù)據(jù)包含兩種刺激：目標(biāo)刺激和非目標(biāo)刺激，比例為1∶8。在原數(shù)據(jù)集中，代號(hào)為nw的被試順利的完成了實(shí)驗(yàn)中的全部trial，同時(shí)原論文中，nw的EEG數(shù)據(jù)分類準(zhǔn)確率最高。所以本文選取原始數(shù)據(jù)集中代號(hào)為nw的被試數(shù)據(jù)作為示例。

3.2 數(shù)據(jù)預(yù)處理

第1步：分幀。

本文介紹的方法需要在分幀之后的數(shù)據(jù)上進(jìn)行。本文按照數(shù)據(jù)集的實(shí)驗(yàn)設(shè)計(jì)進(jìn)行分幀。以每個(gè)聲音刺激事件發(fā)生的時(shí)間作為分幀的時(shí)間起點(diǎn)，選取之后的800 ms作為幀長(zhǎng)，分幀時(shí)間終點(diǎn)就是每個(gè)聲音刺激呈現(xiàn)之后的第0.8 s，舍棄trial之間的腦電數(shù)據(jù)。幀長(zhǎng)的選取與原數(shù)據(jù)集論文的ERP分析相一致，即與聲音刺激有關(guān)的ERP成分波形出現(xiàn)在800 ms以內(nèi)。顯然，這樣的分段會(huì)產(chǎn)生比較多的重疊。這一點(diǎn)會(huì)在后文中的聚類結(jié)果體現(xiàn)出來(lái)。采樣率為150 Hz時(shí)每幀包含120個(gè)數(shù)據(jù)點(diǎn)。

第2步：能量歸一化。

第3步：聚類分析。

作為一種稀疏分解算法，本方法獲得的字典必然是過(guò)完備的，很多字典原子會(huì)具有相似的性質(zhì)或者波形。因此在獲得稀疏字典之后，下一步需要對(duì)字典原子進(jìn)行聚類操作。考慮到本方法是一種非監(jiān)督分析方法，字典原子的種類是未知的，因此本文選取利用自組織映射網(wǎng)絡(luò)進(jìn)行非監(jiān)督聚類，來(lái)找到字典原子分類個(gè)數(shù)的極大值。

本方法中，字典原子時(shí)域向量直接作為自組織網(wǎng)的輸入層，訓(xùn)練獲得網(wǎng)絡(luò)經(jīng)過(guò)訓(xùn)練之后，再次將字典原子送入，即可獲得聚類結(jié)果。

4 結(jié)果

4.1 稀疏性能指標(biāo)

稀疏性能指標(biāo)SPI的結(jié)果如圖1所示。可以看出，SPI會(huì)隨著實(shí)驗(yàn)進(jìn)程發(fā)生變化，表明SPI可以在一定程度上代表EEG的變化。

圖1 每幀的SPI。橫坐標(biāo)為按時(shí)間排序的幀序號(hào)Fig.1 SPI of every frame. Horizontal axis is frame index ordered by time course

4.2 原子聚類結(jié)果

圖2為SOM聚類之后，每類包含的樣本數(shù)。可以發(fā)現(xiàn)字典原子聚類之后，絕大多數(shù)原子都聚集在10類中，原子之間時(shí)域上的相似性被很好的發(fā)掘出來(lái)。

4.3 原子聚類后波形

圖3為每個(gè)類中的原子的平均波形和疊加平均后的目標(biāo)刺激ERP，可以看出，有4個(gè)原子類都體現(xiàn)出了目標(biāo)ERP波形，而且樣本數(shù)據(jù)中由于分幀和實(shí)驗(yàn)設(shè)計(jì)產(chǎn)生存在的具有相位延遲的ERP波形也被提取出來(lái)。我們計(jì)算了疊加平均后的ERP波形與原子波形之間的互相關(guān)。結(jié)果顯示，圖3(b)中的原子波形與ERP波形在0相差時(shí)相關(guān)性最強(qiáng)，達(dá)到0.8，圖3(c)中的原子在有相位延遲的情況下達(dá)到了0.6?；ハ嚓P(guān)結(jié)果說(shuō)明本文的方法達(dá)到了提取ERP波形為字典原子的目的。

圖2 SOM聚類命中圖Fig.2 SOM hits diagram

圖3 (a)所有目標(biāo)刺激的疊加平均后ERP波形,(b)與ERP波形接近的字典原子類平均波形,(c)與ERP波形有相位延遲的字典原子類平均波形Fig.3 (a) Averaged target ERP waveform, (b) Averaged ERP-like waveform of dictionary atoms clusters, (c) Averaged phase-delayed ERP-like waveform of dictionary atoms clusters

4.4 基于稀疏表示系數(shù)的樣本分類結(jié)果

為了驗(yàn)證本方法的有效性，我們利用對(duì)分幀之后的腦電數(shù)據(jù)分解獲得的稀疏表示系數(shù)進(jìn)行原數(shù)據(jù)集論文中所述BCI識(shí)別操作。我們使用了513維的分類特征，其中包括512維的稀疏表示系數(shù)和一維的SPI。樣本標(biāo)記分為兩類，一類是該幀起點(diǎn)播放的是目標(biāo)刺激，另一類是該幀起點(diǎn)播放的是非目標(biāo)刺激。我們選取了SVM作為分類器進(jìn)行二分類識(shí)別。總共3303幀樣本中，二分類準(zhǔn)確率為78%。該發(fā)布數(shù)據(jù)集的論文中，使用了來(lái)自全部64個(gè)通道的數(shù)據(jù)，獲得了80%準(zhǔn)確率。而本方法只使用了1個(gè)通道的數(shù)據(jù)也獲得了近似的結(jié)果。

5 結(jié)論

本文設(shè)計(jì)了一種新的EEG信號(hào)稀疏分解方法來(lái)提取事件相關(guān)電位中的有效成分，并通過(guò)在公開數(shù)據(jù)集上的實(shí)際應(yīng)用，證明該方法是有效的。

與基于匹配追蹤的稀疏分解方法相比，本方法訓(xùn)練得到的稀疏字典具有明顯的生理學(xué)上的含義，可以更精準(zhǔn)可信地利用稀疏系數(shù)來(lái)表示EEG信號(hào)中所包含的信息。在本文中也通過(guò)實(shí)驗(yàn)證明，僅僅通過(guò)單通道EEG數(shù)據(jù)求解字典和系數(shù)，即可實(shí)現(xiàn)單次ERP提取。

本方法中，字典原子的長(zhǎng)度與信號(hào)幀等長(zhǎng)，訓(xùn)練過(guò)程也是以信號(hào)幀為單位進(jìn)行，因此不同的分幀方法會(huì)對(duì)結(jié)果產(chǎn)生影響。這一點(diǎn)在前文結(jié)果展示部分有所體現(xiàn)：不同的延遲的ERP波形會(huì)訓(xùn)練出多種不同延遲的字典原子。因此本方法在應(yīng)用時(shí)會(huì)受到實(shí)驗(yàn)設(shè)計(jì)與分幀方法的雙重影響。

而對(duì)于分幀方法的依賴使得用本方法分析任意的EEG數(shù)據(jù)時(shí)難以獲得穩(wěn)定的結(jié)果。而對(duì)于事件相關(guān)電位的研究來(lái)說(shuō)，如果以事件開始的時(shí)間作為基準(zhǔn)點(diǎn)，只需選取一個(gè)固定的幀長(zhǎng)就可以獲得穩(wěn)定的分幀結(jié)果，因此本方法應(yīng)用在ERP研究時(shí)通過(guò)選取一個(gè)固定的幀長(zhǎng)度即可在不同來(lái)源數(shù)據(jù)之間獲得穩(wěn)定的結(jié)果。