基于連續(xù)小波變換和符號傳遞熵的腦功能網(wǎng)絡構建方法

李明愛，張圓圓

（1.北京工業(yè)大學信息學部，北京 100124；2.計算智能與智能系統(tǒng)北京市重點實驗室，北京 100124；3.教育部數(shù)字社區(qū)工程研究中心，北京 100124）

1 引言

腦-機接口（Brain-Computer Interface，BCI）技術不依賴于傳統(tǒng)的肌肉和神經(jīng)通路，通過采集到的腦電信號直接解讀大腦意識，實現(xiàn)人機或人與周圍環(huán)境之間的通信［1］.運動想象腦電信號（Motor Imagery ElectroEncephaloGraphy，MI-EEG）是指人在沒有執(zhí)行實際動作的情況下大腦想象肢體運動時產(chǎn)生的腦電信號，其包含了大量的神經(jīng)生理信息，已被廣泛應用于BCI領域［2，3］.人腦是一個復雜而密集的網(wǎng)絡，由數(shù)十億相互連接的神經(jīng)元組成，在執(zhí)行運動想象任務時，信息始終在不同腦功能區(qū)之間相互傳遞和整合.近年來，基于圖論的復雜網(wǎng)絡分析方法被廣泛應用于神經(jīng)科學中，通過復雜網(wǎng)絡基本原理可以進行大腦屬性分析，以及發(fā)現(xiàn)腦網(wǎng)絡節(jié)點間潛在的信息傳遞關系.越來越多的研究表明，利用腦功能網(wǎng)絡（Brain Functional Network，BFN）中的度量來區(qū)分運動想象腦電信號具有一定的可行性［4，5］.因此，從BFN 的角度研究不同MI 任務下的大腦活動模式，對揭示運動想象背后的神經(jīng)機制具有重要的意義.

構建BFN 時，節(jié)點通常是EEG 電極，而邊定義為節(jié)點之間的連接性.計算連接性的指標主要分為兩類：功能連接和有效連接.在現(xiàn)有的文獻中［6，7］，大多數(shù)關于腦功能網(wǎng)絡的研究都是基于功能連接度量進行的.Kim等人［8］使用相關系數(shù)（Correlation Coefficient，CC）來估計腦網(wǎng)絡節(jié)點之間的連接性，結果表明，執(zhí)行左手和右手運動想象任務時，網(wǎng)絡節(jié)點的CC 值高于腳和舌頭，且額葉和頂葉區(qū)電極之間的CC 值在也存在顯著不同.Gong等人［9］提出一種基于時頻交叉互信息（Cross Mutual Information，CMI）的腦功能網(wǎng)絡建模方法，研究了四類運動想象任務（左手，右手，腳和舌頭）下大腦的生理機制；他們利用統(tǒng)計分析發(fā)現(xiàn)不同任務下大腦的反應水平、反應時間和激活靶區(qū)都存在明顯的差異.Filho等人［10］提出一種基于圖案同步法（Motifs'Synchronization，MS）的腦功能網(wǎng)絡構建方法，針對BCI2000 數(shù)據(jù)集中十名受試者的左手和右手MI-EEG進行模式識別.結果顯示，在α頻帶和β頻帶上，平均分類準確率均達到了83%.然而，功能連接度量僅能捕獲電極之間在統(tǒng)計意義上的相互依賴性，卻無法提供信息流動的方向以及不同神經(jīng)結構之間的因果關系.

有效連接可描述節(jié)點之間的因果交互作用，相比于功能連接，它能更容易發(fā)現(xiàn)運動想象過程中節(jié)點之間隱藏和被忽視的連接性.Hu 等人［11］使用基于時不變雙變量自回歸（Bivariate Autoregressive，BVAR）模型的格蘭杰因果指數(shù)（Granger Causality，GC）研究了運動想象任務的大腦因果信息流.他們發(fā)現(xiàn)在左手和右手運動想象期間，Cz導聯(lián)與C3/C4之間的方向性連接有很大不同.Ghosh 等人［12］提出一種基于多元自回歸模型（MultiVariate AutoRegressive，MVAR）定向轉移函數(shù)（Directed Transfer Function，DTF）的BFN 構建方法，實驗結果表明，利用網(wǎng)絡節(jié)點度和網(wǎng)絡密度可有效區(qū)分左手、右手和舌頭的MI-EEG 信號.不足的是，GC 和DTF 依賴于線性自回歸模型，使得它們難以準確地反映出人腦這種高度非線性系統(tǒng)中存在的因果關系.

傳遞熵（Transfer Entropy，TE）是一種基于非參數(shù)統(tǒng)計的信息理論量，它不需要預先假設交互作用的模型，可用于衡量兩個非線性向量之間信息交互，是研究大腦信息流的理想方法［13］.符號傳遞熵（Symbolic TE，STE）是傳遞熵的一種變體，具有計算速度快，抗噪聲能力強等優(yōu)點［14］.已有研究人員將其應用于腦功能網(wǎng)絡的構建中，并根據(jù)不同MI 任務下節(jié)點間連接性的差異來區(qū)分運動想象腦電信號［15］.然而，現(xiàn)有方法直接通過計算原始MI-EEG 之間的連接性來構建BFN，難以充分利用腦電信號在頻域上的特征信息.

為此，本文提出了一種基于連續(xù)小波變換（Continuous Wavelet Transform，CWT）和STE 的腦功能網(wǎng)絡構建方法，將EEG 電極的時頻能量序列作為網(wǎng)絡節(jié)點的輸入信號，并使用符號傳遞熵計算腦電極之間的非線性因果交互作用.在BCI 2000 公開數(shù)據(jù)集上進行了實驗研究，以驗證本文方法的有效性和優(yōu)越性.

2 方法

本文提出了一種基于連續(xù)小波變換和符號傳遞熵的腦功能網(wǎng)絡建模方法，其主要步驟如下：

2.1 基于CWT獲取EEG的時頻能量序列

2.1.1 連續(xù)小波變換

連續(xù)小波變換通過小波函數(shù)的伸縮和平移將信號分解成不同尺度上的小波系數(shù)，已被廣泛應用于信號的時頻分析和處理領域［16］.對于任意能量有限信號s(t)，其連續(xù)小波變換的公式為

式中，WΤ 為小波系數(shù)，a為尺度因子，b為平移量，ψ(t)為母小波函數(shù)，ψa，b(t)為伸縮平移后的小波函數(shù)，ψ*(t)為ψ(t)的共軛復數(shù).

2.1.2 獲取MI-EEG的時頻能量序列

2.2 基于符號傳遞熵構建腦功能網(wǎng)絡

2.2.1 序列符號化

圖1 序列符號化示意圖

2.2.2 符號傳遞熵

符號傳遞熵根據(jù)符號序列中各符號出現(xiàn)的相對頻率來估計聯(lián)合概率分布和條件概率分布.對于兩個符號序列，它們之間的符號傳遞熵可通過下式計算：

式中，p(·)為概率密度函數(shù)，τ為時間滯后.

2.2.3 構建腦功能網(wǎng)絡

以L個EEG 電極為節(jié)點，將任意兩導腦電信號時頻能量序列之間的符號傳遞熵作為連接邊，計算腦功能網(wǎng)絡的連接矩陣A：

3 實驗研究

為了驗證本文方法的有效性，本文分析了不同運動想象任務的腦功能網(wǎng)絡，并分別以連接矩陣元素和腦網(wǎng)絡特征參數(shù)作為分類器的輸入，對MI-EEG 信號進行模式分類.

3.1 實驗數(shù)據(jù)及預處理

實驗使用的MI-EEG 數(shù)據(jù)來自于BCI 2000 公開數(shù)據(jù)庫［17］，共記錄了64個頭皮電極的腦電數(shù)據(jù)，電極分布位置如圖2 所示.采集的EEG 經(jīng)過1 Hz～50 Hz 帶通濾波及50 Hz 陷波濾波，采樣頻率為160 Hz.該數(shù)據(jù)庫共包含109 名受試者，每名受試者進行約45 次實驗，想象任務為左手或右手運動，每次實驗持續(xù)約8 s，其中0 s～4 s為運動想象期.

圖2 電極位置與導聯(lián)序號分布圖

構建BFN 的一個關鍵挑戰(zhàn)是體積傳導效應，它會導致頭皮電極之間產(chǎn)生虛假連接性，從而誤導對腦網(wǎng)絡的分析.共平均參考（Common Average Reference，CAR）使用所有導聯(lián)的平均值對信號進行重定位，可有效減小體積傳導效應［18］.因此，本文在預處理階段截取了0 s～4 s想象期的MI-EEG數(shù)據(jù)，并使用CAR濾波器對其進行空間濾波.

3.2 基于連接矩陣的MI-EEG識別

3.2.1 計算有效連接矩陣

根據(jù)2.1 節(jié)，使用連續(xù)小波變換對預處理后的腦電信號進行時頻分解，圖3 為基于Morlet 連續(xù)小波變換對某次實驗C3 導聯(lián)MI-EEG 信號進行分解后得到的時頻能量圖，橫軸為時間，縱軸為頻率，顏色代表小波能量值.在執(zhí)行運動想象任務時，事件相關去同步（Event-Related Desynchronizations，ERD）/事件相關同步（Event-Related Synchronizations，ERS）現(xiàn)象大多出現(xiàn)在Mu節(jié)律（8 Hz～13 Hz）和Beta節(jié)律（13 Hz～30 Hz）上［19］.因此，本文分別截取了時頻能量圖的α（8 Hz～13 Hz）和β（13 Hz～30 Hz）兩個頻帶進行分析.以β頻帶為例，小波分解后β頻帶范圍覆蓋12 個離散頻率，記為：f1，f2，…，f12，將其對應的時間-能量序列依次拼接后得到一維時頻能量序列，如圖4所示.

圖3 C3導聯(lián)MI-EEG的時頻能量圖

圖4 C3導聯(lián)β頻帶的一維時頻能量序列

以64 個EEG 電極作為腦網(wǎng)絡的節(jié)點，計算任意兩個電極時頻能量序列之間的符號傳遞熵（實驗中λ和τ均設為1），得到64×64 的有效連接矩陣.圖5 展示了受試者S1的連接矩陣圖，其橫軸與縱軸均表示導聯(lián)序號，且導聯(lián)與電極之間的一一對應關系與圖2中一致，矩陣中的元素代表兩個導聯(lián)之間的因果連接性.由圖可見，大腦在進行想象左手運動和想象右手運動兩類任務時，節(jié)點之間的連接性具有顯著不同，圖中紫色圓圈標記了兩類想象任務連接性顯著不同的部位.當想象左手運動時，位于中央腦區(qū)右側的C4 節(jié)點與位于前額右側的節(jié)點（FC2、FC4、FC6）之間的連接性遠大于想象右手運動；而想象右手運動時，信息由頂葉右半?yún)^(qū)傳遞至前額葉左半?yún)^(qū)，位于頂葉右側的節(jié)點（如P2、P4）與前額葉左側的FP1 節(jié)點之間的連接性明顯強于想象左手運動.

圖5 受試者S1兩類想象任務的大腦連接矩陣圖

3.2.2 特征選擇和分類

大腦進行運動想象時，并非所有區(qū)域都會被激活，使用連接矩陣中的全部元素作為特征可能不利于MIEEG的識別.因此，有必要對矩陣進行優(yōu)化.

對于每名受試者，本文通過計算大腦連接矩陣各位置元素與類別標簽之間的皮爾遜（Pearson）相關系數(shù)值來優(yōu)化矩陣元素，將相關性較小的元素置零，保留相關性較大的元素.提取優(yōu)化后矩陣的非零元素構成特征向量，采用支持向量機（Support Vector Machine，SVM）分類器進行模式分類，并借助網(wǎng)格搜索法確定SVM 分類器的最優(yōu)超參數(shù)［20］.由于每名受試者的實驗次數(shù)較少，本文使用留一法對樣本進行測試，即每次分類使用一次實驗的MI-EEG數(shù)據(jù)作為測試集，剩余44次實驗數(shù)據(jù)作為訓練集，最終結果取45 次分類結果的平均值.

圖6 顯示了受試者S1 的分類準確率隨皮爾遜相關系數(shù)閾值的變化情況.由圖知，閾值的選擇對分類準確率有著較大的影響，當閾值小于等于0.30 時，分類準確率逐漸升高；而當閾值大于0.30 時，分類準確率則出現(xiàn)了明顯下降；當閾值設為0.30 時，分類準確率達到了最大值100.00%，從而確定S1 的最優(yōu)皮爾遜相關系數(shù)閾值為0.30.依此方法，確定其他每名受試者的最優(yōu)Pearson相關系數(shù)閾值.

圖6 受試者S1的分類準確率

此外，本文還對比了Wrapper 和Fisher 濾波這兩種不同特征選擇算法的分類準確率，結果如表1所示，表中準確率為109名受試者的平均值.從表1中可以看出，當不使用特征選擇算法時（表中None），分類效果較不理想，平均準確率僅為63.95%；而使用特征選擇算法優(yōu)化矩陣元素后，分類準確率有了大幅度的提升，且三種特征選擇算法的準確率相近.使用Pearson 特征選擇算法的準確率最高，達到了94.75%，而使用Wrapper 時的準確率為94.46%，使用Fisher濾波時結果為94.74%.這說明使用有效連接矩陣中的全部元素構建特征向量，會造成特征信息的冗余，并不利于MI-EEG的分類.

表1 還給出了不同特征選擇算法在達到平均分類正確率時所需要的特征數(shù)目，可以看出，當使用Pearson特征選擇算法時，所需特征數(shù)目最少，僅為120個，相比于未進行矩陣優(yōu)化時的4 096 個特征，極大地降低了分類的計算成本.

表1 不同特征選擇方法的平均準確率對比

3.2.3 選取母小波函數(shù)

與傳統(tǒng)腦網(wǎng)絡構建方法不同，本文首先基于CWT求取信號的時頻能量序列，然后將電極時頻能量序列之間的符號傳遞熵作為網(wǎng)絡的邊，進一步提取了EEG信號的時頻特征.為了選擇最適合表達MI-EEG信號時頻信息的母小波函數(shù)，本節(jié)將針對不同母小波函數(shù)的分類結果展開討論.

實驗選取了五種經(jīng)典的母小波函數(shù)進行對比研究，分別為：Haar 小波、Daubechies 小波、Mexican Hat 小波、Meyer小波和Morlet小波.表2同時給出了基于不同小波函數(shù)獲得的109 名受試者的平均分類準確率以及單個受試者的最高準確率，表中的None 則表示未對EEG 信號進行連續(xù)小波變換，僅將各通道EEG 信號之間的符號傳遞熵作為網(wǎng)絡連接邊時的識別正確率.

表2 不同小波函數(shù)的分類準確率

從表中可以看出，不同頻帶下的分類準確率略有差異，且β頻帶的正確率均高于α頻帶，這可能是由于受試者在進行左手、右手運動想象時，β頻帶的ERD/ERS 生理現(xiàn)象更加明顯.同時，無論是在α頻帶還是β頻帶，當以原始EEG 信號之間的符號傳遞熵為邊構建腦功能網(wǎng)絡時，得到的平均分類準確率是最低的，這是由于原始腦電信號難以體現(xiàn)MI-EEG 在頻域上差異，而連續(xù)小波變換通過在時間和頻率兩個尺度上對信號進行分解，可以較好地提取信號在頻域上的特征，進而為正確識別腦電信號提供了更多的有用信息.此外，Daubechies、Meyer 和Morlet 三個小波函數(shù)在兩個頻帶下單個受試者的最高分類準確率均為100.00%，但Morlet 小波的平均分類正確率最高，在α頻帶達到了90.78%，在β頻帶為94.75%；這表明，相比于其他四個母小波函數(shù)，Morlet 小波函數(shù)能夠更好地表達出MIEEG 的時頻特性，并且對于不同受試者能表現(xiàn)出較好的魯棒性.因此，基于Morlet 小波函數(shù)的連續(xù)小波變換將用于本文后續(xù)的實驗中.

3.2.4 確定最優(yōu)符號階次

計算兩個信號之間的符號傳遞熵的第一步是對信號進行符號化，為了不丟失EEG 信號的時間信息，實驗應確保每個符號都盡可能地出現(xiàn)在信號中.根據(jù)文獻［21］，對于長度為Nt的信號，符號階次D應滿足Nt＞＞D！.因此，本文選取D=2，…，6 進行了對比研究，實驗結果如圖7 所示，圖中準確率為109 名受試者的平均分類準確率.總體來看，無論是α頻帶還是β頻帶，隨著符號階次的增加，分類準確率先升高后降低.除此之外，β頻帶的正確率仍然高于α頻帶，這一結果與3.2.3 節(jié)類似；當D=2 時，分類準確率最低，為87.97%；D=4 時，準確率達到了最大值94.75%；而當D繼續(xù)增大時，正確率則發(fā)生了下降，這可能是由于D=4 時最能吻合該數(shù)據(jù)集腦電信號的復雜性.

3.3 基于網(wǎng)絡拓撲特征參數(shù)的MI-EEG分類

3.3.1 計算腦網(wǎng)絡特征參數(shù)

為了進一步研究在不同運動想象任務下大腦各區(qū)域的激活模式，本文計算了兩個重要的腦網(wǎng)絡參數(shù)：度和中間中心性.由于二值化有效連接矩陣可能會丟失部分重要的網(wǎng)絡信息，本文使用文獻［22］中的方法計算加權網(wǎng)絡的特征參數(shù)：

（1）度（Degree）

節(jié)點度的計算式如下：

式中，aij為節(jié)點i和節(jié)點j之間的連接性，節(jié)點的度越大，表明該節(jié)點與其他節(jié)點之間的連接性就越強.

（2）中間中心性（Betweenness Centrality，BC）節(jié)點i的中間中心性可由式（7）計算：

依據(jù)式（6），計算所有節(jié)點的度，構成特征向量FD：

由式（7），計算每個節(jié)點的中間中心性，構成特征向量FBC：

將FD與FBC串行融合，獲得融合特征向量FD+BC：

為了體現(xiàn)腦功能網(wǎng)絡的時變特性，本節(jié)采用窗長為2 s，步長為1 s 的滑動時間窗將原始數(shù)據(jù)劃分為T1(0 s～2 s)、T2(1 s～3 s) 和T3(2 s～4 s)三個時間段，每個時間段數(shù)據(jù)共包含320 個采樣點.對每個時間段信號分別構建腦功能網(wǎng)絡，并計算特征向量.

3.3.2 基于網(wǎng)絡特征參數(shù)的MI-EEG分類

圖8（a）展示了在不同時間段內，想象左手運動和想象右手運動的腦地形圖，使用的特征參數(shù)為度.可以看出，隨著時間的變化，大腦的激活區(qū)域呈現(xiàn)一定的變化規(guī)律，但總體來看，執(zhí)行左手運動想象任務時，位于右側腦皮層的節(jié)點強度更大，而執(zhí)行右手運動想象時，位于大腦左側區(qū)域的節(jié)點度更高.圖8（b）為使用中間中心性作為特征參數(shù)時的腦地形圖，對比于圖8（a），本文也能得到類似的結果.這從復雜網(wǎng)絡的角度進一步闡述了ERD/ERS 現(xiàn)象在運動想象任務中的模式，即大腦在執(zhí)行手部運動想象任務時，對側腦功能區(qū)更有可能被激活.

圖8 不同時間段的腦地形圖

此外，本文還嘗試以度和中間中心性來構建特征向量對兩類MI-EEG進行模式分類.同樣使用Pearson特征選擇算法優(yōu)選貢獻度較大的特征.圖9展示了β頻帶下D=4時T1、T2和T3三個時間段的分類準確率及平均分類準確率.可以看出，僅使用BC值作為特征時，三個時間段的平均準確率為86.20%；僅使用Degree 時，三個時間段的平均準確率為89.83%；將Degree與BC融合時，平均準確率為92.38%.特別地，在T2時段，使用BC和Degree融合特征時，分類準確率達到了93.55%，這一結果接近于使用連接矩陣元素作為特征時的準確率94.75%.此外，圖10還給出了達到圖9分類準確率時所需的特征數(shù)目.可見，當使用網(wǎng)絡特征參數(shù)進行分類時，達到最優(yōu)準確率所需的特征數(shù)目僅為32，這相比于使用連接矩陣元素aij作為特征時的數(shù)目減少了約3/4.

圖9 不同時間段的分類準確率

圖10 獲得圖9準確率所需的特征數(shù)目

3.3.3 多種腦網(wǎng)絡特征提取方法對比

為了驗證本文提出的腦功能網(wǎng)絡構建方法的優(yōu)越性，本節(jié)還對比了不同方法在該數(shù)據(jù)集下的分類性能，結果如表3所示，表中分別計算了將有效連接矩陣元素aij以及腦網(wǎng)絡參數(shù)Degree+BC 作為分類器輸入時109名受試者的平均分類準確率和Kappa 值.從表中可以看出，互信息（CMI）和相關系數(shù)（CC）的分類準確率較低，這是由于它們無法衡量EEG 電極之間連接的方向性，造成特征信息的丟失；而格蘭杰因果指數(shù)（GC）和定向傳遞函數(shù)（DTF）能夠表示信號之間的因果連接性，這種連接性是具有方向的，從而為腦電信號的識別提供了更多的有用信息.不足的是，GC 和DTF 是基于線性回歸模型估計的，將其應用于人腦這樣的非線性動力系統(tǒng)中難免會存在一定的局限性.相比之下，符號傳遞熵更適合反映腦電極之間的交互作用，且這種交互作用是有向的、非線性的，符合人腦的非線性動力學特性，因此取得了最高的分類準確率和Kappa 分數(shù).本文還對比了文獻［10］中使用圖案同步法（MS）在相同數(shù)據(jù)集下的分類結果，當以連接矩陣元素aij為特征時，平均分類準確率提升了11.59%；當以Degree 和BC 為特征時，準確率提高了11.41%，顯示了本文方法的有效性和優(yōu)越性.

表3 基于多種腦功能網(wǎng)絡方法的109名受試者平均分類準確率和Kappa值對比

4 結論

本文提出了一種基于CWT 和STE 的腦功能網(wǎng)絡建模方法，克服了傳統(tǒng)腦網(wǎng)絡構建方法易丟失腦電信號的頻域信息這一不足，并能夠準確地反映不同節(jié)點之間的非線性因果交互作用.在BCI 2000 公開數(shù)據(jù)集上的實驗結果表明，雖然不同運動想象任務的BFN 具有一定差異性，但大腦對MI 任務的反應是局部的.當使用網(wǎng)絡連接矩陣的全部元素作為特征時，分類效果并不理想，而合理使用特征選擇算法可以較大幅度地提高識別正確率；此外，本文發(fā)現(xiàn)在執(zhí)行手部運動想象任務時，位于大腦對側功能區(qū)節(jié)點的Degree 和BC 值大于同側，這從復雜網(wǎng)絡的角度進一步闡述了在進行運動想象時大腦的ERD/ERS 生理現(xiàn)象；同時，本文還嘗試了以網(wǎng)絡特征參數(shù)Degree 和BC 為特征對MI-EEG 信號進行識別，獲得的分類準確率十分接近于直接使用aij時的結果，而特征向量的長度更小，這減少了分類的計算成本.最后，對比了多種腦功能網(wǎng)絡構建方法的分類性能，本文方法均獲得了最高的平均分類準確率和Kappa值，這說明本方法能有效地體現(xiàn)運動想象腦電信號的時頻特征和非線性特征，相比于其他方法具有顯著優(yōu)勢.這一改進有助于腦科學的研究，并將擴大腦網(wǎng)絡分析方法在EEG中的應用.

在未來的工作中，我們將進一步篩選合適的腦網(wǎng)絡拓撲特征參數(shù)，深入研究運動想象過程中大腦的神經(jīng)生理機制以及腦網(wǎng)絡的動態(tài)變化規(guī)律；同時嘗試將腦功能網(wǎng)絡分析方法與其他方法相結合進行多元特征提取，以進一步提高BCI 系統(tǒng)中MI-EEG 的識別正確率.