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

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

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

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

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

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

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

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

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

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

y=DA+e

(1)

(2)

(3)

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

(4)

(5)

(6)

公式(6)中αij，αik表示系數(shù)向量αi的第j、k個元素。該指標(biāo)與重構(gòu)的相對誤差能量呈現(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)有算法都可以通過將SPI作為約束條件，使之具備從根本上優(yōu)化重構(gòu)矩陣稀疏性能、并進(jìn)一步減小誤差的能力。

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

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

(7)

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

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

(8)

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

3 方法驗(yàn)證

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

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

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

第1步：分幀。

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

第2步：能量歸一化。

第3步：聚類分析。

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

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

4 結(jié)果

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

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

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

4.2 原子聚類結(jié)果

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

4.3 原子聚類后波形

圖3為每個類中的原子的平均波形和疊加平均后的目標(biāo)刺激ERP，可以看出，有4個原子類都體現(xiàn)出了目標(biāo)ERP波形，而且樣本數(shù)據(jù)中由于分幀和實(shí)驗(yàn)設(shè)計(jì)產(chǎn)生存在的具有相位延遲的ERP波形也被提取出來。我們計(jì)算了疊加平均后的ERP波形與原子波形之間的互相關(guān)。結(jié)果顯示，圖3(b)中的原子波形與ERP波形在0相差時相關(guān)性最強(qiáng)，達(dá)到0.8，圖3(c)中的原子在有相位延遲的情況下達(dá)到了0.6?；ハ嚓P(guān)結(jié)果說明本文的方法達(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)證本方法的有效性，我們利用對分幀之后的腦電數(shù)據(jù)分解獲得的稀疏表示系數(shù)進(jìn)行原數(shù)據(jù)集論文中所述BCI識別操作。我們使用了513維的分類特征，其中包括512維的稀疏表示系數(shù)和一維的SPI。樣本標(biāo)記分為兩類，一類是該幀起點(diǎn)播放的是目標(biāo)刺激，另一類是該幀起點(diǎn)播放的是非目標(biāo)刺激。我們選取了SVM作為分類器進(jìn)行二分類識別?？偣?303幀樣本中，二分類準(zhǔn)確率為78%。該發(fā)布數(shù)據(jù)集的論文中，使用了來自全部64個通道的數(shù)據(jù)，獲得了80%準(zhǔn)確率。而本方法只使用了1個通道的數(shù)據(jù)也獲得了近似的結(jié)果。

5 結(jié)論

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

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

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

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