基于稀疏低秩張量分解的情緒腦電多域特征提取與分類*

黃金誠高云園佘青山孟 明

(1.杭州電子科技大學圣光機聯合學院，浙江 杭州 310018；2.杭州電子科技大學自動化學院，浙江 杭州 310018；3.浙江省腦機協同智能重點實驗室，浙江 杭州 310018)

隨著信號處理與計算機科學技術的日益發展，基于腦電信號(Electroencephalogram，EEG)的情感計算分析領域也越來越重要。 基于EEG 情緒識別的疾病研究能夠輔助醫生做出診斷，具有客觀、準確的特點，是腦科學中一個重要分支[1－4]。

EEG 是一種以非侵入方式采集的神經心理學信號，在大腦的自發腦電信號中，存在一種由外界事物刺激而產生的誘發電位——事件相關電位(Event-related potential，ERP)。 神經心理學認為，ERP 在一定程度上反映了認知過程中大腦的神經電生理的變化[5]。 利用ERP 潛伏期與波形恒定的特點，結合情緒相關的實驗范式，能夠對情緒障礙相關精神疾病進行輔助診斷與分析，對于研究抑郁癥(Major Depressive Disorder，MDD)、精神分裂(Schizophrenia)等疾病的發病機理與認知偏差有重要價值[6－10]。

近年來，利用機器學習算法對EEG/ERP 進行分析的研究層出不窮，如時域分析、頻域分析、時頻分析和非線性分析等[11－13]。 這些方法以向量或矩陣的形式存儲EEG/ERP 信號，進而提取特征進行分類或分析研究。 由于維度的限制，以向量或矩陣形式存儲的EEG/ERP 信號會損失部分固有的空間信息。

張量作為二維矩陣向高維的擴展，能夠自然地存儲EEG/ERP 的時間、頻率和空間信息，從而可以全面、客觀地分析EEG/ERP。 本文利用張量能夠存儲高維數據的優勢，將多通道EEG/ERP 以張量的形式表示，更完整地表征了EEG/ERP 信號的多域信息。 但傳統張量分解算法對EEG/ERP 張量特征的提取能力有限，其中CP 分解算法(Canonical Polyadic decomposition)多用于腦電分析[14]，提取特征不明顯，造成分類性能較差；Tucker 分解算法計算效率低、結果不唯一且易造成維度爆炸[15]。 本文結合ERP 張量低秩、稀疏的特點，提出稀疏正則低秩逼近Tucker 分解算法(Sparse Regulation for Lowrank approach Tucker decomposition，SLraTucker)，對ERP 多域(時域、頻域、空間域)特征進行提取，該算法能夠在高效分解的同時對特征進行篩選，既保證了特征提取的高效性，又提高了特征的有效性。 在MODMA 數據集上實現了抑郁癥(Major Depress Disorder，MDD)患者與正常對照組(Healthy Control，HC)在多種情緒下的分類。 同時，還采用BCI-IV 數據集驗證了本文算法的泛化能力。 另外，本文充分發揮SLraTucker 分解算法能夠直接提取EEG/ERP空間特征的優勢，對MDD 與HC 兩組人群的空間域特征進行分析，研究他們在多種情緒刺激下活躍腦區的差異，并提取其動態空間特征，實現對空間層面動態變化的對比，為EEG/ERP 空間特征分析提供了客觀的研究方法。

1 材料與方法

本文采用蘭州大學抑郁癥公開數據集MODMA與BCI 競賽數據集BCI-IV，分別對SLraTucker 分解算法的有效性及泛化能力進行驗證。

1.1 實驗范式

蘭州大學抑郁癥腦電公開數據集MODMA 中的ERP 腦電數據集包括24 名MDD 患者(男13 例，女11 例；年齡16～56 歲)和29 名健康對照者(Healthy Control，HC)(男20 例，女9 例；年齡18～55 歲)共53 例受試者的數據，年齡與性別均無明顯差異。 該數據集在經典的點測范式的基礎上，增加不同情緒人臉面孔作為刺激線索，使用圓點作為目標探測點，通過外部情緒圖片刺激誘發注意力任務下的EEG信號[16]，實驗范式如圖1(a)所示。

圖1 實驗范式

BCI-IV-2a 數據集由9個被試者的資料組成，根據不同的刺激分為不同類型的運動想象任務。 首先，簡短的聲音beep 代表記錄開始，兩秒鐘后，箭頭提示出現并在屏幕上停留1.25 s，根據箭頭的朝向受試者被要求執行運動想象任務，直到屏幕再次變黑。 每種運動想象任務都有72 次實驗，以22個電極記錄，實驗范式如圖1(b)所示。 本文主要選取其中的左手、右手運動想象數據，用于驗證SLraTucker分解在不同數據集的泛化能力。

1.2 腦電信號預處理

原始腦電信號中存在大量偽跡信號，需要對其進行預處理以獲取較為純凈的ERP 成分[17]，這里使用EEGLAB 工具箱對采集到的EEG 數據進行離線預處理，步驟依次為:濾波、ICA 去除眼電和心電偽跡、基線校正、分段、疊加平均，如圖2 所示。

圖2 腦電信號的預處理

預處理后的信號包含時域和空間域信息。 本文使用張量對ERP 進行表示，自然地保存了其多域信息，有利于后續全面地提取ERP 信號的多域特征。

1.3 張量基礎及基本運算

張量可以看作矩陣向高維空間的拓展，若為張量選取不同的階數，可以將其拓展為不同的數據類型，其中零階張量是一個標量，二階張量是一個矩陣，三階及以上的張量也就是高維數據空間中的基本數據類型，統稱為高階張量，本文使用高階張量表征情緒ERP 的多域信息，以保留更全面的原始信息。

1.3.1 符號展示

本文小寫字母代表向量，如a∈RI1；大寫字母代表矩陣，如A∈RI1×I2；加粗大寫字母表示張量A∈RI1×I2×…×IN[18]，運算符號如表1 所示。

表1 基本張量運算符號

1.3.2 ERP 信號的張量表示

針對傳統ERP 分析對空間信息保留不全面的問題，本文利用張量模型及其分解算法，表征ERP的多域特征。

經過預處理后的ERP 是以矩陣形式存儲的空間(channels)與時間(sample points)的二維數據，本文采用復Morlet 小波變換得到其頻域信息，建立以張量形式表示的多域ERP。 本文選用的復Morlet小波變換函數為:

對于二階(channels×sample points)的單個ERP信號樣本X∈Rc×t(其中c代表導聯，t代表采樣時間點)，通過計算小波函數，可以得到單個樣本ERP 三階張量(channels×sample points×frequency)，以X∈Rc×t×f表示。

本文所使用的情緒腦電數據由129 電極帽采集(其中包含一個參考電極)，截取情緒刺激后的800 ms 作為一個試次(trial)；每10個trial 進行平均，被試的每種情緒刺激下有16個樣本；設置頻率范圍為0.1 到20 Hz，中心頻率為0.5 Hz，由此，對于每個樣本都能得到一個Xn∈R128×201×40(其中n表示第n個樣本)的三維張量，其中128 表示通道數，201表示采樣時間點，40 表示頻率點。

圖3 單個樣本的ERP 張量

1.3.3 張量的基本運算

給定一個N階(模)張量Y∈RI1×I2×…×IN，其n模展 開( n-mode Unfolding) 操 作 定 義 為Y(n)∈RIn×∏k≠nIk[18]。 張量Y與矩陣A∈RJ×In沿n模的乘積記 作Y×nA， 得 到 的 結 果 是 一個 張 量Z∈RI1×…In－1×J×In＋1×…×IN，也就是將張量沿n模展開后乘以矩陣A。

對于張量Y∈RI1×I2×…×IN的分解降維，有兩種主流的張量分解算法，一種是CP 分解[19，20]，另一種叫做Tucker 分解[21]。 CP 分解是將張量分解為R個秩一張量相加的形式，其中R為張量的秩，且張量的CP 解唯一。 Tucker 分解將原始張量表示為一個核心張量與各個模的因子矩陣相模乘的形式，N階張量的Tucker 分解由式(2)給出，以三階張量為例，Tucker 分解如圖4 所示。

圖4 Tucker 分解

式中:G∈RR1×R2×…×RN為核心張量，A(n)∈RIn×Rn(n＝1，2，…，N)為因子矩陣。 Tucker 分解是矩陣分解的高階拓展，其每一模的秩都可以根據需要人為設置，設置的模秩不同，張量的Tucker 分解也不同；相較于CP 分解，Tucker 分解可以更方便地增加約束項，如非負約束、稀疏約束等，這使得Tucker 分解在高維數據張量分解中泛化能力更強。 考慮到ERP 張量稀疏、低秩的特點，本文選取Tucker 分解提取多域ERP 特征。

由于Tucker 分解的不唯一性，若將每個樣本張量都單獨進行張量分解，取其核心張量作為多域特征，可能存在不同樣本的特征處于不同的高維空間的情況，從而使得特征在后續分析、分類中無意義。為此本文在原有單樣本三階張量的基礎上增加一維，如圖5 所示，將所有的樣本張量在這一維度上排列，得到一個X∈R128×201×40×848的四維張量。 在后續的張量分解中，同時在X的前三維度進行分解，這樣所有的特征張量共享同一套因子矩陣，也就是以相同的張量空間基矩陣進行估計。

圖5 所有被試的ERP 張量

1.4 稀疏正則的低秩逼近Tucker 分解算法

按照上述數據結構生成的ERP 張量規模十分龐大(如本實驗所生成的張量在HDF5 編碼下以雙精度存儲達1G)，使用傳統張量分解優化算法如交替最小二乘法(hierarchical alternating least squares)效率較低[22，23]。 由于張量分解本質上為無監督學習算法，通過對張量分解添加約束往往可以起到減少計算量、提高收斂速度的作用[24，25]。 另外，ERP樣本中存在大量的冗余信息，如背景噪聲、冗余通道等，使得ERP 樣本張量表現出低秩、稀疏的特點。傳統腦電張量分解算法[14]未能考慮腦電在張量域所表現出的特性，直接將在圖像領域較為成熟的張量分解算法應用于腦電上，效率較低且對于腦電張量應用的泛化能力較差。 本文利用Tucker 分解算法的優勢，結合ERP 張量的低秩、稀疏的特點，提出稀疏正則的低秩逼近順序張量分解算法，即SLraTucker 分解算法，實現對ERP 成分的多域特征提取。

傳統Tucker 分解算法在運算時，每一次迭代都需要做一次張量的n模展開(mode-n unfold)，從而造成計算效率較低。 LraTucker 分解(Low-rank approach Tucker decomposition)利用張量的n模展開將傳統Tucker 分解(2)進行重寫:

得到矩陣Y的低秩逼近表示后，通過固定和，將其帶入式(4)得到新的目標函數

LraTucker 分解算法結合數據張量的低秩性，表現出十分高的效率，但與傳統Tucker 分解一樣，在EEG 中應用易造成維度爆炸，有一定局限性。 為了提高在EEG/ERP 張量分解中的泛化能力，并使得提取的特征更加明顯，本文在LraTucker 分解中添加稀疏正則項λS(A)

式中:S(·)為稀疏表示，如L0 范數｜A｜0和L1 范數｜A｜1在實際應用中都存在一些問題(求解L0 范數是一個NP 難問題，L1 范數在0 點不可導)，這里我們使用下式作為稀疏表示函數:

利用梯度下降法對式(6)進行更新，分別對A和B進行求導，得到:

選取α為與最終得到SlraTucker 分解的更新公式如下:

由于只需進行一次n模展開，使用如上的更新規則求解Tucker 分解較傳統方法可以少進行倍的n-mode 展開，且加入的稀疏正則項使得該算法在EEG/ERP 張量中具有更強的泛化能力。

1.5 特征提取流程

本文SLraTucker 分解算法對情緒ERP 張量進行特征提取流程如圖6 所示，得到的核心張量用于分類，表征空間特征的因子矩陣以腦地形圖的形式呈現，以便分析空間特征。 具體特征提取過程如下:

圖6 特征提取流程示意圖

(1)采用50 Hz 陷波濾波器與0.3 Hz～30 Hz 帶通濾波器濾波；重參考。

(2)分段平均:以刺激的出現作為刺激開始的標志，截取刺激后的800 ms 作為一個trial，采樣頻率為250 Hz。 對于每個個體，將同一刺激下的10個樣本進行疊加平均，以提取較為純凈ERP 成分。

(3)生成樣本張量:①將疊加后的樣本按照通道依次通過復Morlet 小波變換映射到時頻域，頻率范圍取0.5 Hz～20 Hz，步長為0.5 Hz，則對于每個通道都得到一個201 × 40 的二維矩陣(time ×frequency)；②將變換后的數據按通道的順序在張量的第一維度排列，得到一個維度為128×201×40 的三維張量(channels×time×frequency)，自然地包含了每個ERP 樣本的空間、時間、頻率信息；③將得到的三維張量在第四個維度排列，每個被試的每種情緒刺激下有16個樣本，生成一個維度為128×201×40×(16×53)的四維張量(channels×time×frequency×(subject×trials))，此四維張量即為樣本張量。

(4)使用SLraTucker 分解對樣本張量的前三維度同時進行張量分解，得到的核心張量為包含了原始信號中多域特征的特征張量，因子矩陣可以看作樣本所共有的多域信息。

按照上述步驟對ERP 數據張量進行特征提取，得到的核心張量是每個樣本所特有的特征在核心張量第四維度的排列，因子矩陣為所有樣本共享的共同特征。 如此得到的腦電多域特征包含了空間、時間和頻率特征，且經過張量分解能夠直接得到空間特征，相比傳統方法更加直接。

2 結果展示

2.1 分類結果

利用本文的算法提取情緒ERP 的多域特征，將每個被試所特有的特征——核心張量經過向量化后輸入到SVM 中，分別在愉快、悲傷和恐懼情緒刺激下獲得91.5%、90.6%和84.3%的識別率，較傳統Tucker 分解平均提升13.4%，較LraTucker 分解平均提升4.7%。 為直觀展現張量分解作為特征提取方法的優勢，本文將其與通過改進內核目標對齊(modified kernel-target alignment，mKTA)[26]進行空間特征篩選的方法進行比較。

文獻[27]采用mKTA 作為其空間特征選擇算法，選取的通道數從1 增加到58，通過支持向量機(SVM)得到平均分類結果，在選取19個通道時識別準確率達到最高，即80.0%。 本文提出的SLraTucker 算法能夠直接在特征提取的過程中進行空間特征篩選，加入的稀疏正則也一定程度上進行了特征篩選，保留了更加明顯的特征，從而提高了模型的分類性能。 如表2 所示，SlraTucker-SVM 分類方法實現的患者平均識別率為88.9%，相較于mKTA方法，分類效果提升了8%以上。

表2 MODMA 數據集實驗結果

為驗證本算法在EEG 領域的泛化能力，選取BCI-IV 數據集，進行與MODMA 數據集相同的濾波預處理，并生成張量，分別應用傳統Tucker 分解、LraTucker 分解和本文提出的SLraTucker 分解算法對原數據進行特征提取。 在內存為8GB，CPU 為i5-8265u 的平臺上，對BCI-IV-2a 數據集的9個被試張量進行分解，傳統Tucker 分解算法平均運行時間為1 000.726 2 s，達到平均76.72%的原始數據張量相似度；LraTucker 分解平均運行59.759 4 s，達到平均78.72%的相似度；而SLraTucker 分解算法平均運行時間為57.638 0 s，較傳統Tucker 分解計算效率提升約16 倍，達到78.66%的相似度。 提取出的核心張量對于原張量的還原度如表3(fit mean)所示，選取10 折交叉驗證SVM 對其分類，準確率如表3 與圖7 所示.

表3 BCI 數據集實驗結果

圖7 BCI 數據集準確率

2.2 靜態空間特征提取結果

采用SLraTucker 分解算法對MDD 與HC 兩組人群不同情緒刺激后800 ms 內的多域特征進行提取，該算法是在前三維度上對全部樣本張量進行估計的過程，得到的因子矩陣為其所有樣本所共有的特征。為此，本文通過分別生成MDD 與HC 各自的樣本張量，利用SLraTucker 算法對兩個樣本張量分別分解，得到MDD 患者樣本張量與HC 對照組樣本張量的Tucker 分解。 為觀察兩類人群在不同情緒刺激下活躍腦區的異同，將表征空間特征的mode-1 因子矩陣以腦地形圖的形式呈現，腦區位置如圖8 所示。

圖8 腦地形圖腦區標注

圖9 中展示了兩類人群在悲傷情緒刺激下的腦地形圖，從HC 對照組可以看出，大腦主要在額葉和頂葉比較活躍，而MDD 患者頂葉的活躍度較低，枕葉活躍度略高于HC 對照組，額葉右側活躍度較HC對照組更大，但活躍程度度偏低。

圖9 悲傷刺激下腦地形圖對比

從圖10 中可以看出，HC 對照組在“愉快”情緒刺激下表現為額葉活躍，MDD 患者雖然活躍區域與HC 對照組類似，但枕葉活躍程度較HC 對照組偏低，且頂葉存在活躍。

圖10 愉快刺激下腦地形圖對比

對于“恐懼”刺激，如圖11 所示，HC 人群活躍腦區在左前額葉、右前額葉有較高活躍度，枕葉也存在一定活躍度，與HC 對照組不同的是，MDD 患者的高活躍區主要在右前額葉和頂葉，枕葉則表現為低活躍狀態。

圖11 恐懼刺激下腦地形圖對比

2.3 動態空間特征提取結果

為了動態地分析ERP 的變化情況，我們選取長度為200 ms，重疊率為50%的滑動時間窗。 結合HC 對照組與MDD 患者800 ms 腦地形圖的差異從128個通道中選擇了包括額葉、中央區、頂葉、枕葉的20個電極，如表4 所示，針對每類人群分別生成數據張量，利用SLraTucker 算法對每個時間段內的樣本張量分別進行分解，把核心張量中對應空間特征的mode-1 因子矩陣以腦地形圖的形式呈現。

表4 通道選擇

根據圖12，在“悲傷”情緒刺激下，MDD 額葉左側的活躍從前期出現保持至后期，而HC 的額葉左側活躍度則從前期出現，并于中期開始下降。 在0～800 ms 時間段內，MDD 患者活躍腦區存在不對稱的現象，而HC 則較為對稱。

圖12 悲傷刺激下動態腦地形圖對比

根據圖13，在“愉快”情緒刺激下，MDD 活躍腦區變化并不明顯，主要集中于頂葉右側、中央區右側和右前額葉，有小幅能量波動，而HC 對照組的額葉左側、中部、右側和枕葉等都有能量上升以及活躍區域增大。

圖13 愉快刺激下動態腦地形圖對比

通過圖14 我們發現，在“恐懼”情緒下MDD 患者活躍腦區向左前額葉移動，中期左前額葉出現活躍并持續至后期，HC 對照組的活躍腦區則向兩側額葉移動，從中期開始，左右兩側額葉的活躍度都有所上升。

圖14 恐懼刺激下動態腦地形圖對比

3 討論

本文以張量的形式存儲、表示情緒ERP，利用ERP 張量稀疏、低秩的特點，提出SLraTucker 分解算法對情緒ERP 張量進行多域特征提取。

SLraTucker 分解在計算過程中只需要進行一次n模展開，較傳統Tucker 分解少進行次n-mode 展開，同時加入稀疏正則化，使得SLraTucker算法在ERP 分析中表現出更強的泛化能力。

由于稀疏約束對于權重較低的特征有置零的作用，相比LraTucker 分解，SLraTucker 分解對于張量的還原能力存在小幅下降，但得益于稀疏約束特征選擇剔除冗余特征，SLraTucker 分解算法在分類準確率上較LraTucker 分解平均提高4.7%，較傳統Tucker 分解算法提升了17.5%。

按照二維情緒模型，從效價和喚醒兩個維度分析[27]，“恐懼”為負性高喚醒情緒，“悲傷”為負性低喚醒情緒，“愉快”為正性低喚醒情緒。 通過對刺激后800 ms 時間段的腦地形圖進行分析，在效價維度上，在效價為負的“悲傷”、“恐懼”情緒刺激下，MDD患者的額葉表現出低于HC 對照組的活躍度；在正性情緒“愉快”刺激下，MDD 患者額葉活躍度低于HC 對照組，而中央區和頂葉活躍度則要高于HC 對照組。 在喚醒維度上，“恐懼”為高喚醒情緒刺激，MDD 患者的頂葉、中央區三部分活躍度高于HC 對照組。 MDD 患者額葉活躍度低于HC 對照組，且MDD 患者更易被負性情緒激活。 額葉是情緒調控的高級中樞，也是認知、感覺和運動等多級協調體系的重要中樞。 姚樹橋等[28]認為抑郁癥患者對于正負性刺激的不同ERP 效應，可能涉及前額葉左、右半球不同區域的功能整合；Kayser 等[29]認為抑郁癥患者的腦地形圖表現出大腦不對稱性，造成評價過程出現異常，這也與我們通過分析刺激后的800 ms的腦地形圖得出的結論一致。 此外，由于情緒的產生是一個動態的過程，本文通過增添滑動時間窗口的方式，探究了在不同情緒誘發下兩組人群活躍腦區的動態變化。

根據圖12～圖14，可以發現MDD 患者對效價的敏感度高于HC 對照組，而對于喚醒度則剛好相反。對于負性情緒，MDD 患者對喚醒度更強的“恐懼”情緒刺激不敏感，其額葉左側活躍出現時間更晚且持續時間較短，反觀HC 對照組則較為敏感，額葉左側早期就出現活躍并持續到后期；MDD 患者對于效價更低的“悲傷”情緒刺激更加敏感，額葉左側活躍從早期出現，并在中期、后期持續保持活躍狀態，活躍腦區不對稱的現象明顯，HC 對照組早期產生的左前額葉活躍在中期即開始下降。 對于正性情緒，MDD 患者活躍腦區變化十分微小，在左前額葉、右中央區和右頂葉有小幅度能量上升，頂葉表現出激活困難的現象；反觀HC 對照組的額葉活躍度變化幅度較大。 綜上所述，在負性情緒刺激下，MDD 患者的額葉腦區活躍持續時間更長，且存在不對稱現象。 對于正性情緒MDD 患者的腦區更難被激活。

結果表明，MDD 患者對于效價低的情緒更加敏感，腦區活躍程度大、持續時間更長，而對于正性情緒則存在選擇性抑制，表現出難以被激活的現象。Leuchter 等人[3]的研究結果顯示，抑郁癥患者在前額葉區、頂葉區、顳葉區相干性連接與正常人存在差異，Kayser 等[29]認為抑郁癥患者由于右頂葉情緒刺激感知區域的選擇性抑制，會造成評價過程出現異常，這與我們動態腦地形圖得出的結論一致。 另外，由本文算法提取的動態腦地形圖顯示，MDD 患者對于效價的敏感度高于HC 對照組，而對于喚醒度的敏感度則低于HC 對照組，這也進一步驗證了800 ms 靜態分析的結果。

4 結論

張量是分析、研究EEG/ERP 等高維數據的強有力工具，結果表明，本文提出的SLraTucker 分解算法能夠高效地提取情緒ERP 張量中的多域特征。在MODMA 數據集上驗證了其有效性，并在BCI-IV-2a 數據集上驗證了其具有較強的泛化能力。 由SLraTucker 分解提取出的空間特征以腦地形圖呈現，比較和討論了在不同的認知情緒處理過程中，抑郁癥患者與正常對照人群活躍腦區的區別與聯系，通過動態分析發現MDD 患者對于效價的敏感度高于HC 對照組，對于喚醒度的敏感度則剛好相反。 本文提出的SLraTucker 分解算法為基于情緒腦電的分析提供了新的方法和思路。