基于多層腦功能網絡特征的動作意圖識別

常文文，聶文超，袁月婷，閆光輝，楊志飛，張冰濤，張學軍

(蘭州交通大學電子與信息工程學院 蘭州 730070)

在康復訓練過程中，使用機器人設備比傳統的人工康復更具優勢，能增加患者訓練的動機和自主訓練的機會，從而提高康復訓練的質量和效果。在步態康復方面，外骨骼和智能助行機器人被廣泛應用，并取得了良好的效果[1-2]。隨著腦機接口(brain computer interface, BCI)技術的發展，基于BCI 的主動式康復可以通過更快、更有力地檢測大腦的運動意圖來改進康復策略，提高康復效果，是未來神經康復系統發展的趨勢[3-4]。理解大腦認知活動和運動過程之間的關系，對于發展基于腦機接口的主動康復技術具有重要作用。

目前的研究主要通過測量腦電信號(elecrroencephalogram, EEG)對運動意圖進行特征解碼和模式識別研究[1,5]，EEG 測量由于其簡單、便攜和高時間分辨率，被廣泛應用于運動意圖的檢測中。文獻[6-7]表明，EEG 中包含了豐富的步態和運動信息，而針對行走和步態等下肢運動意圖的解碼研究剛剛起步。在步態過程中最基本的一個動作過程就是起立和坐下， 文獻[3]開發了一種實驗協議，完成了對被試者在起立－坐下和屈膝－伸展動作過程中腦電信號和肌電信號的同步采集，通過事件相關同步/去 同 步(event related synchronization/desychronziation)和慢皮層電位以及與下肢相關的肌電模式分類來分析被試的動作意圖，并結合支持向量機和人工神經網絡完成對不同運動階段的分類識別，得到了很好的效果。文獻[6]通過對EEG 信號中低頻波段信號的解碼，研究了10 名被試者從坐姿到站姿以及從站姿到坐姿轉變過程中對應的腦電信號特征，在保留EEG 數據豐富的統計特征基礎上，結合高斯混合模型(Gaussian mixture model,GMM)，從時間嵌入的低頻EEG 信號中成功識別“站到坐”和“坐到站”的動作意圖，也證實了類似基于BCI 的下肢解碼技術在外骨骼康復機器人中應用的可行性。文獻[4]同樣專注于動作執行前的皮層慢電位，對Delta 波段在運動意圖解碼中的可行性展開分析，基于動作開始前1.5 s 的慢波數據，通過設計自我觸發和外界提示觸發兩種模式對起立、坐下和安靜3 種狀態下的EEG 信號特征展開分析，通過偽跡子空間重構算法消除對EEG 信號中的高幅噪聲，通過構建時間嵌入的EEG 特征向量并應用Fisher 判別分析降低其時空特征的維數，最后結合GMM 分類器完成對3 種狀態的識別，得到了較好的效果。文獻[8]最早基于高密度的EEG 測量，就被試在站立和行走過程中執行視覺反應任務時的大腦認知過程做了詳細研究，對EEG 數據的獨立分量分析顯示，在站立、慢走和快走期間的視覺事件相關電位在不同運動條件下沒有差異，證明了在運動過程中記錄和認知有關的大腦活動的可行性。文獻[9]對被試者在坐著、站立和行走過程中從事某種認知任務時大腦認知負荷的變化過程進行研究，基于干電極的腦電信號采集，通過設計聽覺誘發的事件相關電位分析，發現與坐著相比，行走時P3 顯著降低，表明行走時的認知負荷高于其他兩項活動，而坐著和站著的P3 無顯著差異。文獻[10]提出了一種基于熵的信號復雜度度量方法，并通過共空間模式濾波器選擇了可用于特定動作分類的特征量，最后利用線性判別分類器實現對站立和行走動作的識別。綜上可以看出，有關運動意圖的解碼研究主要集中為兩種信號，一種是事件相關的同步/去同步電位，另一種是運動相關腦電位(movement related potentials, MRPs)[1,4]。此外，文獻[8-9]針對不同動作過程中認知能力變化的研究，實現了對該類動作的有效檢測和識別。但是，這些研究主要集中在感覺運動腦區少數電極通道的慢電位上，缺乏從全腦角度出發的空間域特征信息表示以及相應的識別方法研究。

步態是一個復雜的認知和運動控制過程，下肢動作也涉及全腦各個腦區的協調和配合，在一個站立或坐下動作完成之前，人的大腦必然表現出一定的特征信息，可以通過對該類特征信息的解碼最終確定其運動意圖[4]。這些信息除了通過上述皮層慢電位的表征之外，通過對腦區間相互依賴關系的動態變化過程分析，有望揭示運動意圖解碼的新特征[11-12]。文獻[13]就站立和行走過程中的腦功能網絡特征展開分析，發現相比站立狀態，行走過程中感知運動腦區的功能連接會減弱。文獻[11]就是否有外骨骼輔助下運動過程中的腦功能網絡特征展開分析，也表明了基于圖論的腦網絡分析方法在步態康復研究過程中具有一定的特殊作用。文獻[12]通過構建功能性腦網絡，對腦創傷性病人行走過程中的平衡性問題展開分析，也證明了基于腦功能網絡的空間域信息對表征步態特征具有特殊效果。針對動作意圖檢測開展多角度、深層次的分析和探究顯然是必要的。為此，本文設計了從靜坐到站立、然后再到靜坐的動作轉換實驗，并同步采集了被試者的腦電信號，通過構建多層腦功能網絡，就動作執行前后腦區間依賴關系的動態變化過程展開分析，并就不同動作對應多層網絡特征參數進行了統計分析。

1 實驗設計

為實現對站立、坐下和靜止狀態的識別，結合統計差異的網絡特征參數和機器學習算法，本文提出了一種該類動作意圖識別方法，其技術路線圖如圖1 所示。

圖1 本文技術路線圖

1.1 實驗方案

13 名被試者完成此次實驗，被試者均為在校學生，身體健康，沒有任何的精神疾病病史，均為右利手，實驗之前和被試者簽署了實驗知情協議書。實驗過程中被試者坐在顯示器前，顯示器顯示計時時間，要求被試者先靜坐1 分鐘，然后被試者起立，由靜坐轉為站立，被試者站立持續1 分鐘，之后被試者坐下，由站立轉為靜坐，靜止持續1 分鐘，依次反復交替完成上述動作，如圖2 所示[3]。

圖2 站立?靜坐實驗過程示意圖

在該過程中全程記錄被試者的腦電信號，持續10 min 為實驗的一個Session，每名被試者總共完成6 個Session 的數據采集，各Session 中間休息1～2 min。

1.2 數據采集和預處理

整個動作實驗過程中的數據采集使用SAGA 32+腦電采集系統(荷蘭TMSi 公司)，在線采樣頻率為1 024 Hz，在線采樣頻段為0～200 Hz，數據記錄中采用平均參考模式，額頂處為接地電極。該系統除了采集30 通道的腦電信號之外，還有4 組附加通道，通過附加通道接入了兩組外接傳感器，其中一組是慣性測量單元(inertial measurement unit，IMU)傳感器，將其固定在被試者左上腿位置，另一組是4 路足底壓力傳感器(FSR 型壓力傳感器)，以足底壓力鞋墊的形式固定在被試者左腳下方，通過IMU 和足底壓力傳感器數據，來確定起立和坐下的起始時間點，實現和腦電數據的同步，傳感器安裝位置如圖3 所示。實驗中記錄的信號包括30 通道腦電信號、4 路足底壓力傳感器數據以及IMU 數據，所有數據的預處理均在EEGLAB 中完成。對所有被試者的處理逐個進行，首先將每名被試者所有Session數據整合到一起，然后對數據降采樣到250 Hz，并完成0.1～48.0 Hz 的帶通濾波，使用Cleanline工具完成對信號的線性濾波，根據SAGA 系統電極位置完成對所有腦電信號通道的電極定位，并剔除壞電極和無關通道。在完成數據讀入的基本處理后，對IMU 的各路數據進行低通濾波(Butterworth濾波器，截止頻率10 Hz)，并借助下肢上IMU 和腳后跟的足底壓力信號確定每個站立動作和坐下動作的起始時間點，作為EEG 信號事件的標記點(將事件標記為起立、坐下和靜態3 個狀態)。根據該標記點完成對所有數據的截斷，為了討論動作前后的腦電特征，每個Epoch 的長度設置為動作起始點之前4.5 s 到之后的2.5 s，最后得到356 個Epoch數據段，所有的Epoch 數據段都進行了基線校準。然后通過視覺檢查剔除受偽跡影響較大的Epoch數據段，利用運行獨立成分分析(run independent component analysis, RUNICA)操作完成對所有信號的ICA 分解，并基于半自動的獨立成分分析選擇(semi-automated selection of independent component analysis, SASICA)工具包完成對偽跡成分的判斷，并剔除其中的偽跡成分，主要包括眨眼、眼動、肌肉偽跡以及由于動作過程引起的其他偽跡成分。在ICA 剔除偽跡之后，如有需要，對壞電極進行插值處理，最后對所有數據進行平均重參考，從而完成數據的預處理過程。

圖3 足底壓力鞋墊和傳感器位置固定示意圖

2 研究方法

2.1 同步似然分析

行走是人類一項復雜的認知活動，涉及視覺判斷、運動控制、空間感知和規劃等多個方面，需要整個大腦不同腦區間的協調和配合。在一個動作執行之前，通常大腦就表現出一定的動作意圖，通過討論動作執行前后相當時間內不同腦區間的相互依賴關系，有助于全面了解動作執行過程中大腦的活動狀態，并有望發現針對動作意圖檢測的新特征?；诙鄬幽X功能網絡的復雜網絡分析方法為此類動作過程對應腦區間依賴關系的研究提供了工具。

表征腦電信號之間相互依賴關系的度量方法有很多種[14]，其中同步似然分析(sychronization likehood, SL)被廣泛應用于非線性和非平穩信號間依賴關系的分析中，特別適合于腦電信號的分析，對EEG 信號中的線性和非線性特征都具有較好的表征能力，且對容積傳導等問題具有較好的魯棒性[15-16]。SL 是神經生理學數據中估計廣義同步性的常用指標，它的理論依據是系統在某一時刻的狀態可通過一個插值向量來表示，當該狀態再次出現時，可以用一個相似的插值向量來對其進行表示。因此，本文采用SL 來構建不同通道腦電信號之間的相互依賴關系，從而構建多層腦功能網絡。分為以下4 個步驟。

1) 對來自通道C的EEG 信號xC,i(i=1,2,···N，N為樣本點數)進行重建，得到m維空間嵌入之后的信號：

式中，l為時間延遲。同理，對通道D的EEG 信號進行重構得到XD,i。

2) 對通道C上的第i個樣本，設定各向量之間彼此的距離小于s的概率為：

4) 在完成上述步驟之后，就可以實現對通道C和通道D之間EEG 信號同步似然指數SL 的計算，其表示的是通道C信號和D信號同時重現其狀態向量的可能性，即當XC,j和 給定XC,i之間的距離小于鄰接距離sC,i， 同時XD,j和XD,i之間的距離也小于sD,i的可能性。本文在計算SL 時，對應參數依 次 設 置 為：l=10、m=10、 ω1=100、 ω2=400、pref =0.01。

2.2 多層腦功能網絡構建

對某段時間內腦功能網絡拓撲結構動態關系的一種表征就是對該時間段內信號按一定時間窗進行截斷，依次構建功能網絡并形成超鄰接矩陣，然后借助多層網絡分析方法進行特征分析。為了通過腦功能網絡來表征站立、坐下以及靜止狀態的腦信號特征，首先提取了和運動意圖相關的delta 慢波信號，如圖1 所示，然后根據站立和坐下以及靜止狀態的時間點，將3 類狀態對應的delta 波段信號分割為[?4.5 s, 2.5 s]的Epoch 數據段(長度為7 s)，并以3.5 s 的時長作為時間窗，0.5 s 的步長對每個Epoch 依次進行截斷，從而得到8 個時間段的子窗口(依次為w1,w2,···,w8)，這樣對一個Epoch 樣本進行操作，就可以得到30×875×8 的數據樣本，然后對每個窗口內以及每對窗口間使用SL 來計算兩個通道間的相互依賴關系。圖4 為層間和層內各通道SL 計算示意圖，通過層間和層內各通道之間SL 計算最終可以得到多層網絡結構并以超鄰接矩陣的形式表示。

圖4 多層腦功能網絡構建

圖5 為3 種狀態下所有被試、所有Epoch 數據段對應超鄰接矩陣的總平均顯示，該矩陣為240×240 的SL 值矩陣表示，對角線位置為8 個子窗口(w1～w8)的窗內連接關系(層內連接)，而上三角或下三角矩陣為不同子窗口之間的連接關系(層間連接)，可以看到層內連接強度要遠高于層間連接。由于圖示超鄰接矩陣中SL 為總平均值，3 種狀態之間的差異并不明顯。為此，圖6 計算了該超鄰接矩陣中SL 的總平均值，可以看到3 種狀態之間，SL 均值在動作發生期間要高于靜止狀態，且坐下和起立動作對應均值接近，但是統計分析并沒有發現3 種狀態之間存在顯著差異。

圖5 3 種狀態對應EEG 信號SL 超鄰接矩陣

圖6 3 種狀態SL 超鄰接矩陣總平均值

2.3 多層腦功能網絡特征參數分析

為了就整個動作過程對應的特征量進行量化分析，本文重點討論4 類多層網絡特征參數，首先是層?層相關系數(layer-layer correlation, LLC)，將其定義為兩個網絡層之間的皮爾遜相關系數[17-18]，用于表示各網絡層之間，特別是不同子窗口對應網絡層之間的相關性程度。為了就多層網絡中不同層之間整合和分離程度進行度量，文獻[17]引入了一個度量指標，該指標量化了第m層中的鏈路關系仍然存在于第n層中的條件概率(conditional probability, CP)，將其定義為：

式中，Am(m=1,2,···,8)為層內鄰接矩陣。如果P值較高，表示各層之間整合關系較強，而如果各層之間關系相對比較獨立，則P值較低。此外，不同層之間的分離特征也可以用多層模塊值來進行量化，文獻[19-20]將其定義為：

最后，多層網絡還可以使用多層參與系數(multilayer participant coefficient, MPC)來量化其特征。該系數與熵有關，主要用來量化給定節點i的鏈路是均勻分布在各個層中，還是主要集中在一個或幾個層中[17-18]。節點i的MPC 值被定義為：

2.4 閾值選擇和單層腦功能網絡特征參數分析

上述SL 值構成的超鄰接矩陣表征了層內和層間的相互依賴關系，是對同一時間段或不同時間段各通道之間相互依賴關系的一種表征，在這種連接關系中也存在著大量的弱連接邊，可能是由于噪聲或其他非特征性相互作用造成的，并不能有效表征不同動作狀態下腦區間的真實依賴關系。因此，在一般性問題的討論過程中，需要設定一個合適的閾值對上述超鄰接矩陣中的連接邊權值進行篩選，保留在不同狀態間具有一定區分度的重要連接邊。為此，首先將3 種狀態對應SL 超鄰接矩陣權值展平，并組合在一起構成一個一維向量，然后通過選擇該向量取值分布圖上的上分位點作為最終的閾值點，圖7 為該向量的頻率直方圖分布以及確定的閾值(閾值T= 0.1)。

圖7 3 種狀態對應SL 超鄰接矩陣值統計分布

對上述3 種狀態對應所有超鄰接矩陣使用本閾值進行重點權值邊的篩選，發現由于層間連接都較弱，保留下的重點連接邊都表征了層內各通道之間的連接關系。因此，對各層內連接以單層網絡的方法進行對比分析，重點考慮了各層網絡的度、聚集系數、全局效率、同配系數、傳遞性系數和節點的介數中心性[14,21]。網絡中最基本的一個參數就是節點度(degree, D)，定義為與該節點相連的邊的權值之和，并將網絡中所有節點度值的平均作為網絡的度(average degree, AD)。同時，主要對網絡的聚集系數(clustering coefficient, C)、全局效率(global efficiency, GE)、網絡的同配系數(assortativity coefficient, ASS)、網絡的傳遞性系數(transitivity,Tra)以及節點的介數中心性(betweenness centrality,BC)等網絡參數展開分析。

3 結果分析與討論

3.1 多層腦功能網絡特征

在得到如圖5 所示各狀態對應的超鄰接矩陣之后，本文分別計算了各Epoch 片段對應的超鄰接矩陣的4 種多層網絡特征參數，圖8 為對應的Q值和MPC 值，可以看出Q值保持在0.8 左右，而MPC值在0.9 以上，兩組參數值都較高，表明不管是靜止狀態，還是起立和坐下動作，各網絡層之間的連接關系較為密切，整個動作過程中大腦活動表現出很強的連續性。對3 種狀態之間的統計分析也發現，Q值在3 種狀態間有較為明顯的差異(t-test,p<0.05)，而MPC 值并沒有表現出顯著性差異。

圖8 不同狀態超鄰接矩陣對應Q 值和MPC 值

圖9 分別對應3 種狀態下各網絡層之間條件概率CP 矩陣和層?層相關系數LLC 矩陣?？梢钥吹礁鲗又g的CP 值都普遍較高，均在0.6 以上，表明各子窗口對應網絡層之間存在較強的整合關系，不管是相鄰時間窗的網絡層(如w1 和w2)，還是間隔較遠的時間窗的網絡層(如w1 和w8)，各層之間均存在較為緊密的聯系，但是3 種狀態之間并無顯著差異。相比之下，圖9 中兩層之間相關性值情況有所不同，可以看到兩層之間的LLC 值隨著兩層間距離的增加而降低，表現出相鄰層有較高的LLC 值，而間隔較遠的層之間有較低的LLC 值。為了能夠找到可應用于不同運動狀態分類的特征量，對LLC 值同樣進行了統計分析，圖9 中標記為兩種狀態之間存在統計性顯著差異的情況(ttest，p<0.05，其中 “*” 為起立 vs. 坐下， “★”為起立 vs. 靜態，“▲” 為坐下 vs. 靜態)。可以看到，對于總平均后的LLC 矩陣值并沒有明顯不同，但是統計分析發現在w4 和w6,w7,w8 之間，以及w6 和w7,w8 之間的相關性存在統計性顯著差異，這可能和站立和坐下動作開始的時間點有一定關系。特別是從w4 子窗口開始，出現了起立和坐下的動作起始點，因此和后續子窗口的網絡在層相關性上表現出顯著性差異。

圖9 不同狀態對應網絡層?層相關系數矩陣

3.2 閾值篩選后網絡特征

為了突出上述網絡連接中的重點連接邊，將閾值設定為0.1，完成對上述3 種狀態所有Epoch 數據網絡權值的篩選，大于該閾值的連接邊保留，小于該閾值的連接邊權值設置為0。權值篩選后的網絡最終為各層內連接網絡，層間連接權值相對較小，因此均被忽略。圖10 為3 種狀態各子窗口對應閾值篩選后總平均網絡連接拓撲圖。可以看出，3 種狀態在8 個時間窗上拓撲結構變化比較平穩，對站立和坐下狀態，均從w4 開始，網絡連接較之前的w1、w2、w3 開始變得稀疏，并且一直持續到最終的w8，這和圖10 中網絡層之間相關系數顯著性差異結果一致。本文認為這主要是由于在整個時間窗的移動過程中，站立和坐下動作的起始點最開始均出現在w4 窗口，而在站立和坐下過程中的動作電位相應地影響了不同腦區間的這種依賴關系，且在動作完成之后的一段時間內，這種依賴關系得到了一定程度上的保持。從圖中也可看出，在靜止狀態下，網絡拓撲結構保持平穩，不同腦區間的拓撲關系幾乎沒有發生變化，整個網絡拓撲結構的變化同起立及坐下狀態有一定區別，但無法通過該動態拓撲圖來區分起立和坐下狀態。

圖10 3 種狀態各子窗口對應腦功能網絡拓撲連接

為了便于進行量化分析，確定可用于不同狀態 識別的新特征量，本文重點分析了3 種狀態下各窗口對應的網絡特征參數。如表1 所示，分別計算了聚集系數C、全局效率E、平均度AD、同配系數Ass、介數中心性Bet 和網絡傳遞性系數Tra 這6 類基本的網絡特征參數。其中，“*”表示ttest，p<0.05。由表可知，動作狀態對應網絡特征參數在8 個遞進子窗上整體呈現出變小的趨勢，但是變化差異很小，而在靜止狀態網絡參數變化沒有明顯的趨勢性。對上述各網絡參數在3 種狀態之間的統計分析表明，只有w7 窗口對應的網絡傳遞系數Tra 表現出顯著的統計性差異。

表1 各狀態對應的網絡特征參數

3.3 機器學習分類

通過上述對各網絡特征參數的對比和統計分析，本文選擇了多層網絡特征參數Q、有顯著差異的層間相關系數LLC 以及w7 子窗口對應網絡傳遞性系數Tra 構成特征向量：

使用支持向量機(support vector machine, SVM)、邏 輯 回 歸(logistic regression, LR)、線 性 判 別 分析(linear discriminat analysis, LDA)和樸素貝葉斯(navive bayes, NB)4 類機器學習算法對該特征向量進行測試，完成對兩類狀態的識別，在機器學習分類過程中將70%的數據作為訓練集，剩余30%作為測試集。圖11 為5 折交叉驗證得到的平均分類準確率?？梢钥吹剑槍深悹顟B的分類準確率均在71%以上，其中SVM 分類器的效果最好，針對站立和坐下、站立和靜態以及坐下和靜態的分類準確率分別達到74.4%，77.0%和76.9%，證明了該方法的有效性。

圖11 不同分類器對應分類結果

4 結 束 語

本文通過設計靜坐到站立的動作轉換實驗，完成對被試者在起立、坐下和靜止狀態下對應不同腦區腦電信號的同步似然分析，通過構建動作執行前后8 個時間窗口上的腦功能網絡，以超鄰接矩陣表示的多層網絡形式展開對不同動作狀態下網絡拓撲結構的特征分析。基于對多層網絡特征和各子窗口網絡參數的統計分析，確定對不同動作狀態敏感的網絡特征參數，并使用SVM、LR、LDA 以及NB 等分類算法實現對起立、坐下和靜止狀態的分類識別，分類準確率最高達到77%，證明了基于動作前后動態網絡特征進行下肢運動意圖檢測和識別的可行性。針對多層動態腦功能網絡的分析結果表明，站立和坐下動作的發生會弱化腦區間的相互依賴關系，從而導致網絡拓撲連接結構變得逐漸稀疏，該結果也為基于腦電的動作意圖檢測提供了新思路，對開展基于腦機接口的下肢助行和康復機器人的開發具有一定的參考價值。

致謝：本文得到了蘭州交通大學天佑青年托舉人才計劃的基金資助，同時也感謝博?？甸_放課題(BRKOT-LZJTU-20220329C)對文章的支持。