采用源空間套索分析和卷積神經網絡方法的高頻腦電動作模式識別方法

陶怡,徐維維,朱家林,袁子文,2,王茂德,王剛

(1. 西安交通大學生物醫學信息工程教育部重點實驗室,710049,西安; 2. 西安交通大學第一附屬醫院 康復醫學科,710061,西安; 3. 西安交通大學第一附屬醫院神經外科,710061,西安)

腦-機接口(brain computer interface, BCI)通過分析人執行動作前大腦活動的變化,讓設備能更早地識別并完成相應的指令[1]。在這個過程中,頭皮腦電(electroencephalogram, EEG)由于其非侵入、安全性高等特點,成為許多研究中分析的對象[2-3]。EEG信號是來自大腦皮層的時間序列數據,通常根據頻率分為5個頻帶[4],包括δ(1～3 Hz)、θ(4～7 Hz)、α(8～13 Hz)、β(14～30 Hz)、γ(30～100 Hz)。一般認為,不同頻帶的信號活動在大腦中的分布與特定的認知功能之間存在一定的關聯[5]。

以往的研究表明,δ節律多與認知加工、疲勞水平、代謝、睡眠等過程有關[6-8];θ節律與信息整合和工作記憶及多種認知功能的整合和協調有關[9];α節律在信息處理速度和注意力等領域起重要作用[10],此外,位于感覺運動皮層的成分被稱為μ節律,與軀體的感覺刺激、四肢的運動密切相關,并在BCI中已有廣泛的應用[11];β節律常出現于工作、運動、思考等大腦活動強烈的情況下[12-13];γ節律被認為與情緒或者更廣泛的認知過程相關[14],例如刺激選擇、記憶形成等。

大部分研究表明,與手部動作解碼相關的信息主要包含在α和β頻帶,Zhang等通過對EEG進行0.1～30 Hz濾波,對手部運動方向進行分類,在左、中、右3個方向上得到了(73.39±6.35)%的準確率[15]。然而,近年來一些學者發現,在γ頻帶上也存在許多利于手部動作分類的信息[16-17],在以往的研究中也有同樣的發現,當使用多元經驗模態分解(MEMD)方法把信號分解成不同頻帶的多元本征模式函數(MIMF)后,將單個MIMF輸入支持向量機(SVM)中,位于高頻帶的3個MIMF有著最高的分類準確率[18];通過前向序列搜索策略輸入卷積神經網絡(CNN)中,發現被選中的MIMF也多為最先分解出的高頻分量[19]。因此,本文中主要針對γ頻帶的信號進行分析。

EEG是由神經元放電后經過大腦皮層傳播到頭皮被頭皮表面電極檢測并記錄的信號,在傳播過程中信號幅值有一定程度的衰減。此外,鄰近區域的神經元電信號也會對解碼工作帶來一定的干擾。因此,有研究使用腦電源定位技術將頭皮腦電投射到代表大腦幾何形狀的源空間上,將信號從傳感器空間轉換到源空間進行單手多類動作模式的識別[20],提高腦-機接口的性能。

不論是傳感器空間還是源空間內,以往的手部動作分類研究主要依賴手動提取信號特征,包括時域、頻域和時頻域上的特征。然而,這種方法導致基于EEG的腦機接口面臨著高維度特征的問題。為了應對這個問題,許多國內外的學者選擇從特征選擇和通道選擇兩個方面進行降維。

在特征選擇方面,唐肖芳和周金治[21]使用散度分析進行腦電信號特征選擇,Liu等[22]使用螢火蟲方法對從運動想象信號中提取的空域和頻域融合特征進行選擇,楊鵬圓等[23]使用方差貢獻率與F-score結合的特征選擇方法,在不降低準確率的情況下大量減少了特征數。

在通道選擇方面,孟明等[24]結合互信息進行挑選,曹玉珍等[25]提出了基于卷積自編碼器和費舍爾準則的EEG通道選擇方法。最小絕對值收斂和選擇算子方法(Lasso),簡稱套索,是一種用于特征選擇的新方法[26],組Lasso是其重要推廣,能夠在組的水平上進行變量選擇。王金甲等[27]通過融合EEG信號功率譜、小波變換參數等特征,提出了稀疏組Lasso方法,對信號的通道和特征同時進行選擇,表明該方法優于傳統方法。

腦電信號復雜多樣,手動提取信號特征具有一定的局限性,一是提取的特征不一定適合要研究的任務,二是容易造成信息的損失。隨著深度學習的興起,學者們開始嘗試卷積神經網絡在腦電信號解碼方面的潛力。2019年,Amin等提出了一種多層神經網絡融合方法,在公共數據集BCI競賽Ⅳ-2a上取得了75.7%的分類準確率[28]。2020年,Mammone等對源空間信號進行連續小波變換,生成的時頻圖作為深度卷積神經網絡的輸入,得到六分類62.47%的分類效果[20]。2022年,Fujiwara等訓練了一個包含殘差網絡的多層CNN模型,對休息、左手運動和右手運動任務分類準確率達85.69%[29]。本文提出了source-Lasso-CNN(SLC)方法,將EEG信號映射到源空間,再通過組Lasso進行源空間感興趣區域(ROI)的選擇,最后輸入到CNN中進行自動特征提取,從而達到對手部動作分類的目的。

1 數據和方法

1.1 被試者與實驗范式

實驗征集了13位年齡在20～32歲的健康被試者(全部為男性),依次被標識為S1～S13。被試者均為右利手,視力或矯正后視力正常,并且接受該研究的書面說明和簽署知情同意書。

實驗在遠離電氣設備、高頻輻射的環境中進行,以降低外界環境對腦電信號產生的干擾。整個實驗過程要求被試者保持平靜、注意力集中,盡量減少體動。

腦電和肌電信號通過Neuroscan 64導腦電儀采集,接地電極位于前額區Fz和FCz之間位置,參考電極位于雙耳乳突,按照國際10-20標準電極安放法定義的位置,采集了FP1、FP2、F7、F3、FZ、F4、F8、FT7、FC3、FCZ、FC4、FT8、T7、C3、CZ、C4、T8、TP7、CP3、CPZ、CP4、TP8、P7、P3、PZ、P4、P8、O1、OZ、O2一共30個通道的信號,EMG信號采用一次性Ag/AgCl電極組成4對表面差分電極,以肘部為參考點,分別從前臂的指伸肌、橈側腕伸肌、掌長肌和尺側腕屈肌等處同時采集4路表面肌電信號,采樣頻率均為1 kHz。腦電、肌電電極放置位置如圖1(a)、圖1(b)所示。每名被試者需要完成4種手部動作模式的EEG和EMG信號采集,包括握拳、展拳、二指對捏、三指對捏,每種動作進行60個試次,共采集240個樣本。為降低實驗過程中疲勞帶來的影響,以20個試次為一組,每組之間休息5 min,單個試次包括準備階段4 s、動作執行4 s和休息4 s,某被試者的實驗過程如圖1(c)所示,單次實驗時序圖如圖2所示。

(a)國際10-20標準電極位置

(b)肌電電極放置位置

(c)被試者信號采集

圖2 單次實驗時序圖Fig.2 Sequence diagram of single experiment

1.2 信號預處理

對記錄到的數據進行0.1～100 Hz帶通濾波和50 Hz陷波濾波,以去除工頻噪音干擾。使用移動時間窗法[30]對肌電信號進行分析,從而判斷動作真實起點位置,進而提取動作前2 s數據進行后續分析。

1.3 源空間套索分析和卷積神經網絡動作模式識別方法

本文提出的基于源空間套索分析和卷積神經網絡(source-Lasso-CNN,SLC)方法,首先將傳感器空間的EEG經過腦電源定位映射到源空間,再提取源空間信號的標準差、均方根幅值和均方根作為特征進行組Lasso,將源空間中挑選出的腦區及相應信號輸入CNN中自動提取特征并進行分類,具體SLC方法流程如圖3所示。

圖3 SLC方法流程圖Fig.3 Flow chart of proposed algorithm

1.3.1 源空間

腦電溯源定位,也即腦電逆問題,是根據頭皮記錄到的EEG信號,反向估計腦內神經活動源的位置、方向和強度信息。要解決逆問題,就要先對神經元電活動傳播到頭皮上的過程有一定的了解,頭皮腦電信號與腦內源活動的數學模型可以建立如下

b=AX+n

(1)

式中:b表示記錄的頭皮數據,在本文中是大小為30×1 000的矩陣,其中30為傳感器空間腦電通道數,1 000為裁剪后每個樣本的時間采樣長度;X表示腦內源活動的信息,大小為15 002×1 000,15 002是所選模型的腦源數,1 000為映射到源空間后的時間序列長度;A表示傳遞矩陣,體現了腦內各源對EEG每個通道貢獻的大小;n表示測量噪聲。逆問題就是在已知b的情況下反推X,而在此之前,需要先解決腦電正問題,即建立一個大腦中可以闡述大腦皮層神經元活動是如何經過腦容積效應和顱骨的傳導投射到頭皮上的頭模型,即矩陣A。頭模型的主要部分包括大腦的幾何形狀以及神經電活動傳播所通過的不同組織(灰質、白質、腦脊液、顱骨、皮膚)的導電特性。本文基于FSAverage模板的真實頭模型,利用邊界元法建立了代表頭皮、顱骨和大腦內部組織的三殼傳導模型,電導率為1∶1/8∶1。

理論上,腦電逆問題的求解就是使得記錄的頭皮腦電和已知的偶極子計算出的對應電位之間的殘差最小的過程,但一般來說,逆問題的解不唯一,所以需要通過確定一個既符合生理學解剖意義,又在數學上計算方便的源模型來得到唯一的腦源解。

源模型一般分為等效電流偶極子模型和分布式電流密度模型,前者通過單個或少量偶極子等效某個區域的腦源活動,但只適用于局部區域的精準估計,無法對于整個大腦的源活動進行準確的重構。因此本文選擇分布式源模型,在這種模型中,信號源均勻地分布在大腦皮層中,由于腦內源數15 002遠遠大于腦電圖的通道數30,這將導致腦電逆問題擁有無窮多解,因此需要加上合適的約束條件。本文采用最小范數估計(minimum norm estimate, MNE)[31]的方法來求解逆問題,其求解公式如下

(2)

X的最優解將是在L2范數最小二乘意義下的估計和記錄的頭皮腦電之間的殘差最小的解,式中η≥0,是正則化參數;W是一個權重矩陣。

得到映射到源空間的每個源的信號后,按照Brodmann分區對處于同一腦區的偶極子進行累加,作為這個腦區整體的輸出,以便進行后續分析,腦區共有79個,腦電源定位-Brodmann分區如圖4所示。

圖4 腦電源定位-Brodmann分區Fig.4 Spatial source localization-Brodmann area

1.3.2 組套索方法

Lasso[26]是Robert Tibshirani 在 1996 年提出了一種懲罰類的回歸方法,其基本思想是使用懲罰函數對回歸系數的絕對值進行懲罰,要求所有回歸系數的絕對值之和小于等于懲罰參數λ,通過調節λ,進而對回歸系數進行控制壓縮。隨著λ的增大,更多的自變量回歸系數被壓縮成 0,從而達到降維的效果,其表達式如下

(3)

1.3.3 卷積神經網絡

CNN是一種前饋型網絡,其結構受到人腦視覺皮層的啟發。與傳統機器學習方法相比,CNN最大程度地減少了預處理的過程,能夠自動提取信號特征、學習復雜模型并進行分類,因此在圖像識別、視頻分析和自然語言處理等領域被充分關注。

CNN通常由輸入層、隱藏層和輸出層構成,其中隱藏層可由多個卷積層、池化層、全連接層組成。

卷積層是CNN的核心層,通過卷積核對數據進行掃描,從而完成卷積運算,具體公式如下

(4)

卷積層后通常與激活函數相結合,以實現對輸出結果的非線性映射,從而提高網絡的非線性功能,常用的激活函數有指數線性單元函數(ELU)、修正線性單元(ReLU)、雙曲正切激活函數(tanh函數)等,本文中使用ELU函數,其公式為

(5)

式中,x為輸入的特征值。

池化層又稱下采樣層,主要用于壓縮數據,以達到減小過擬合和降低參數數的目的。常用的池化操作有最大池化和平均池化,最大池化取區域最大值作為輸出,平均池化取區域中的均值作為輸出。在相鄰兩個卷積層之間可以設置池化層以簡化網絡復雜度,本文采用最大池化。

全連接層把所有神經節點互相連接,方便做最后的分類或者回歸,本文在全連接層后還設置了權值為0.5的dropout,即在訓練過程中使50%的神經元失活,以提高模型的泛化能力,進一步降低過擬合。

方法中構建的CNN網絡包括輸入層、卷積層1、卷積層2、批標準化層1、最大池化層1、卷積層3、批標準化層2、最大池化層2、卷積層4、批標準化層3、最大池化層3、展平層、全連接層1、全連接層2和輸出層。在將數據輸入CNN前,先對樣本進行標準化并打亂順序,輸入大小為R×T的矩陣,其中R為挑選出的ROI數,T為每個通道的時間采樣點數,本文中R≤79,具體值由組Lasso過程確定,T=1 000。輸入的第一層卷積層設25個卷積核,每個卷積核隨時間執行核大小為1×5的卷積操作,第二層卷積層中,每個卷積核執行對空間信息的卷積提取,這兩層之間沒有激活函數,通過這種方式可以將整個卷積強制分離線性變換為時間卷積和空間濾波器的組合實現隱式正則化,降低計算量。在兩個卷積層后加入一個池化層,使用最大池化的方式,池化核的大小是1×3,步長與核大小一致。第3、4個卷積層分別由50、100個卷積核構成,后接一個相同的池化層,通過這種方式可以減少池化層對整個網絡結構丟失過多信息的影響,激活函數均設置為ELU。展平層是對高維特征進行扁平化處理,將其轉換為一維從而能夠輸入到全連接層中。全連接層共有2層,分別設置有1 024和512個神經元,使用Relu函數作為激活函數。最后使用softmax函數輸出分類預測概率,概率最大的類別即為分類結果。CNN詳細模型結構見表1。在卷積層中,層參數指的是卷積核數;而對于全連接層和輸出層來說,層參數指的是神經元數。

實驗中采用6折交叉驗證方法,把所有樣本數據隨機分為6組,輪流選3組作為訓練集、2組作為驗證集、1組作為測試集,最后6次平均的結果作為模型的識別準確率。

由于深度學習要求大量樣本進行訓練,使用1 s的時間窗對數據進行裁剪,從而達到數據擴充的目的,剪切的移動步進設為50 ms。通過剪切操作,2 s的腦電數據可以擴充成為21個1 s的腦電片段。并且由于在剪切之前已經劃分訓練集、驗證集和測試集,因此存在重疊信息的剪切數據不會對分類結果產生影響,從整體上來說,單個被試者2 s的4種單手動作共240個樣本,在剪切操作后成為1 s的4種單手動作5 040個樣本。

表1 CNN模型結構

1.4 實驗環境

方法測試環境為Windows10操作系統,腦電源定位通過基于matlab2019b的brainstorm工具箱實現,組Lasso使用R語言中的msgl庫實現,CNN通過tensorflow2.1、python3.6實現,服務器顯卡型號為NVIDIA RTX3080Ti。

2 采用源空間套索分析和卷積神經網絡的實驗結果

2.1 挑選出的感興趣區域數

EEG信號經過腦電源定位后映射至源空間內,并根據Brodmann分區分為79個腦區,為了降低腦區數并選擇對于手部動作意圖解碼更有利的區域,本文對信號進行了組Lasso后,每位被試者挑選出的ROI數如圖5所示,其中橫坐標為13個被試者,縱坐標為挑選出的通道數,每個被試者共有6個點,代表每折交叉驗證時挑選出的結果,灰色虛線為13位被試者的平均值。可以看出,經過組Lasso,每位被試者的腦區數大大下降,平均只有35.46±6.00個通道,降低為原來的44%,大大節省了后續需要的計算資源,提高了計算效率。

圖5 組Lasso后挑選出的ROI數Fig.5 Number of ROI selected after group Lasso

2.2 方法性能評價

本文使用分類準確率對提出的方法解碼能力進行評估。分類準確率定義為所有分類正確的樣本數與樣本總數之比,取六折交叉驗證的平均值作為最終結果以驗證本文方法的有效性。CNN訓練過程準確率和損失值曲線如圖6所示。

(a)被試者S3時模型的準確率及損失值

(b)被試者S4時模型的準確率及損失值

(c)被試者S6時模型的準確率及損失值

以被試者S3(分類結果較差)、S4(分類結果一般)、S6(分類結果較好)為例,模型訓練集和驗證集的準確率及損失值曲線如圖6所示。由圖6可以看出,隨著迭代次數的增加,S4和S6的模型快速趨于收斂,但是S3的模型表現并不理想。考慮到該被試者在其他模型上分類效果也不佳,可能是被試者的個體差異或數據質量的影響。

訓練過程中設置最大迭代次數為80次,batch塊的最小尺寸設為64,同時通過早停法判斷是否需要提前終止訓練,即連續10個epoch驗證集準確率不再上升時,提前結束訓練。由于挑選出的ROI通道數不同、訓練迭代次數也有差異,每個被試者訓練模型時間長短略有不同,平均耗時(84.49±14.82) s。

表2給出了13位被試者模型在γ頻帶上的分類準確率,可以看出,13人的分類準確率最高達到了97.2%,最低有61.21%,平均準確率達到了(82.23±12.71)%,證明本文提出的方法在單手多動作模式識別上的有效性。

表2 13位被試者模型在γ頻帶上的分類準確率

3 實驗分析及討論

3.1 不同模型參數對解碼效果的影響

CNN的結構包括卷積層數和卷積核的大小,本文通過參考文獻[19, 32]以及進行的預實驗確認網絡卷積層數為4,為了探討不同卷積核大小下模型解碼性能,文中對大小為1×3、1×5、1×7、1×10、1×15和1×20的卷積核進行測試,不同卷積核尺寸下分類準確率結果見表3。可以看出,卷積核大小為1×5的分類準確率最高,之后隨著卷積核的增大,分類準確率逐漸下降,考慮到小卷積核可以提取更多細節特征,最終確定1×5的卷積核進行研究。

表3 不同卷積核尺寸下分類準確率

池化層和dropout是卷積神經網絡中常用的正則化方法。大尺寸的池化層可以減少特征圖的尺寸,提高計算效率,但容易導致信息丟失,而小尺寸的池化層可能導致過擬合問題。為了確定合適大小的池化層,本文對尺寸為1×2、1×3、1×4、1×5和1×6的池化層進行測試,結果見表4。由表4可以看出,在尺寸為1×3的情況下,模型的平均分類準確率最高,因此在后面的測試中均使用此尺寸的池化層。對于dropout正則化方法,它通過隨機將一定比例的神經元輸出設為0來避免過擬合。較大的dropout值可以提高模型的泛化能力,但是過大的dropout值會導致模型的學習能力下降。為了確定合適的dropout值,文章在0.2～0.7范圍內進行測試。不同池化層大小和dropout值下分類準確率見表4,當dropout設置為0.5時,模型的準確率最高,這也是大多數研究中使用的值。

表4 不同池化層大小和dropout值下分類準確率

3.2 基于不同頻帶的解碼效果對比

以往的研究認為和手部動作解碼相關的信息多集中于α和β頻帶,但近期也有很多研究表明在γ頻帶也存在相關信息,為了討論不同子頻帶的源空間信號的解碼效果,本文將不同頻帶,即δ(1～3 Hz)、θ(4～7 Hz)、α(8～13 Hz)、β(14～30 Hz)、γ(30～100 Hz)及全頻帶(1～100 Hz)的信號,分別進行源空間定位,再進行相應的組Lasso進行感興趣區域選擇,最后輸入到CNN中,將所有頻帶的分類結果進行比較,結果如圖7所示。各頻帶分類準確率為(35.97±5.10)%、(35.25±5.98)%、(40.08±5.40)%、(55.26±11.06)%、(82.23±12.71)%和(72.13±14.29)%。通過Wilcoxon符號秩和檢驗發現,γ頻帶的分類準確率顯著優于其他頻帶(p<0.5),結果說明不同子頻帶的信號都為手部動作意圖解碼提供了一定的信息,但是可以認為,更多的有用信息與γ頻帶相關,在以后的研究中應該更加關注此頻帶。

圖7 γ頻帶與其他頻帶對比結果Fig.7 Comparison results between γ band and other frequency bands

3.3 與其他方法的腦電解碼分類準確率對比

圖8 本文方法與其他方法對比結果Fig.8 Comparison results between our algorithm and other algorithms

為了更客觀地對本文提出方法進行評價,文章將所提出方法和其他先進方法的分類準確率進行對比,包括GCB-Net[33]、EEGNet[34]、DeepConvNet[32]、DeepCNN[35]、濾波器組公共空間模式(filter bank common spatial patten, FBCSP)[36]和MEMD-CSP[37],結果如圖8所示。對于GCB-Net、EEGNet、DeepConvNet、DeepCNN 4種深度學習網絡,使用的輸入數據只經過簡單預處理,并未提取任何其他特征。在MEMD-CSP方法中,信號首先進行MEMD分解;再使用MIMF的中值頻率自動找到對μ和β節律有重要貢獻的被試特定MIMF;將選定的MIMF相加以重建EEG信號;使用CSP提取重構信號的特征,并且將空間濾波器數設置為5;通過SVM分類器對特征進行分類。在FBCSP方法中,信號通過CSP提取特征,濾波器頻帶設為0.1～10、10～20、20～30、…、80～90、90～100 Hz,并通過SVM進行分類。6種方法的平均準確率分別為(56.34±11.05)%、(68.02±10.93)%、(67.06±9.78)%、(55.67±8.92)%、(66.86±10.17)%和(60.35±10.08)%,均顯著低于本文提出的方法(p<0.05),這些結果表明,與最先進的深度學習方法和傳統的MEMD-CSP、FBCSP方法相比,本文所提出的方法都能獲得更高的分類精度。

4 結 論

本文提出了一種SLC方法,通過腦電源定位方法將傳感器空間采集的γ頻帶的EEG信號映射到源空間中。接著,使用組Lasso方法對ROI進行挑選,并將挑選出的腦區信號輸入到卷積神經網絡中進行自動特征提取及分類。通過六折交叉驗證,平均分類準確率為(82.23±12.71)%。與其他頻段的腦電信號分類結果對比發現,γ頻帶的分類結果最好,甚至優于在全頻帶上的效果。這表明對大腦高頻成分活動的研究可能對手部動作解碼具有重要意義,為今后的研究提供了新的思路和方向。