基于多元經驗模態分解-傳遞熵的腦肌電信號耦合分析*

馬鵬剛,佘青山,高云園,張啟忠,羅志增

(杭州電子科技大學智能控制與機器人研究所,杭州 310018)

生理系統是跨越多尺度時空進行相互作用的復雜結構[1]。運動神經系統通過神經振蕩傳遞運動控制信息,該振蕩活動可以反映多層次的皮層-肌肉功能耦合FCMC(Functional Corticomuscular Coupling)的連接信息[2]。耦合分析研究可以探索大腦皮層和肌肉組織間的相互作用[3],代表了運動系統內部的信息流動,與大腦皮層向肌肉收縮的傳入反饋有關,因而有助于了解大腦如何控制肌肉組織以及肌肉組織對腦功能的影響,更利于理解皮層-肌肉的功能耦合及信息傳遞[4-6]。

研究多時空層的腦肌電耦合信息,關鍵是要提取不同尺度上的腦電EEG(Electroencephalogram)和肌電EMG(Electromyographic)信息。經驗模態分解EMD(Empirical Mode Decomposition)是一種用于分析非線性和非平穩性時間序列的自適應數據處理和挖掘算法[7-8],本質上是對數據序列的平穩化處理。該方法使用一種篩選方法把信號分解為若干本征模態函數IMF(Intrinsic Mode Function)[9],每個IMF分量包含了原信號在不同尺度上的局部特征信號。但是,EMD存在兩點問題:一是只能處理一元時間序列,不能夠同時處理多元時間序列,并且IMF數量不固定;二是對于復雜的數據,數據間歇性不僅在時頻域造成了嚴重性的模態混疊,而且單個IMF分量的物理意義不清晰。為了解決上述問題,Huang等[10-11]提出了一種噪聲輔助的信號處理方法——集合經驗模態分解EEMD(Ensemble Empirical Mode Decomposition)。該方法本質上是在原有信號中加入高斯白噪聲來平滑脈沖干擾,從而在一定程度上避免了模態混疊,但其計算量繁瑣,并且無法對多元時間序列同時分解[12]。Dragomiretskiy等[13]提出了一種變分模態分解VMD(Variational Mode Decomposition)方法,可以同時分解多個通道,提取信號中具有中心頻率的窄帶分量。然而,該方法不適合處理非穩態信號,其頻譜會隨時間發生較大的變化,也會導致信號的混疊。最近,Mandic等[14]提出了多元經驗模態分解MEMD(Multivariate Empirical Mode Decomposition)方法,不僅能夠同時處理多元時間序列的多尺度分解,而且可以較好地消除模態混疊現象。目前,該方法被應用到腦肌電信號的多尺度分析中,自適應地提取不同頻率段的窄帶分量[15-17],取得了較好的實驗結果。

傳遞熵是概率論函數的一種信息理論函數,用于度量信號之間的信息交換,已經被證明非常適用于分析時間序列,并且在可預測性和信息傳遞的基礎上探索其因果關系[18]。由于傳遞熵考慮的是變量間的信息量傳遞,而不需要假定變量間具有特定形式的關系,因此比Wiener-Granger因果分析方法具有更好的適用性,尤其是對于非線性系統[19]。

為了定量研究不同尺度上腦肌電信號之間的非線性同步耦合特征及其信息傳遞,本文結合多元經驗模態分解與傳遞熵方法,構建了多元經驗模態分解-傳遞熵MEMD-TE(Multivariate Empirical Mode Decomposition Transfer Entropy)模型,用于腦電、肌電耦合分析。首先采用低通濾波的方法對腦電、肌電信號進行預處理,其次對處理過的腦電、肌電信號采用MEMD方法進行時-頻尺度化,然后計算腦電、肌電不同時頻段的傳遞熵值,最后得到各個尺度上腦電、肌電之間的功能耦合特征。在不同靜態握力(5 kg、10 kg、20 kg)下進行抓握過程中,研究不同尺度上的腦電、肌電耦合,定量描述腦電、肌電信號間耦合特征和信息傳遞方向,為進一步探索運動功能障礙的產生機理及康復過程中的運動功能評價方法提供依據。

1 實驗數據采集和預處理

1.1 實驗數據采集

人類在做動作時需要大腦皮層運動和感覺區的參與和協調,研究表明運動皮質區呈左右交叉、上下顛倒的模式支配身體所有的運動和控制方式,即某一側掌管的是對側的肢體運動功能,比如左腦與右手動作相關,右腦與左手動作相關[20]。因此本文重點對C3與右手部分肌肉和C4與左手部分肌肉進行了實驗分析。

本實驗共有10名受試者(6名健康人S1～S6,均為右利手,25歲;4名腦卒中患者S7～S10,其中S8為左利手,其余均為右利手,46.5±1.0歲),10人均被告知實驗詳情,簽署了知情同意書。本實驗采用128導BrainAmp DC腦電采集系統同步采集EEG和EMG信號,記錄了32路EEG信號,以及左右手橈側腕屈肌FER(Flexor Earpi Radialis)和肱三頭肌TB(Triceps Brachii)處的EMG信號,實驗采集環境如圖1所示。

圖1 實驗采集環境

在電極安放之前,頭皮必須要吹洗干凈,皮膚表面也需要用酒精擦拭干凈,采樣頻率為1 000 Hz,采樣時間為5秒。6名健康受試者左右手分別進行五組5 kg、10 kg和20 kg的抓握,而對于4名患者,考慮到其身體狀況,因此只進行了左右手五組5 kg、10 kg的抓握。實驗中,每一次抓握工作完成后均會有15 s的休息時間,防止肌肉疲勞對EMG信號造成干擾。

1.2 EEG和EMG的數據預處理

由于EEG和EMG信號十分微弱,極易受到干擾,本文利用數字濾波器對腦肌電信號分別進行0～75 Hz[21]的低通濾波,用于腦電、肌電同步分析[22-23]。

2 腦肌電信號的多元經驗模態分解-傳遞熵分析

為了研究不同時頻尺度間的腦電、肌電非線性同步耦合特征,本文構建了多元經驗模態分解-傳遞熵模型,應用于腦電、肌電信號的耦合分析中。首先對腦電、肌電信號進行多元經驗模態分解得到不同頻率成分的尺度信息,然后使用傳遞熵的分析方法,計算出腦電和肌電之間的傳遞熵值,得到不同頻段上腦電、肌電信號的功能耦合特征。

2.1 多元經驗模態分解

MEMD可以實現多元信號振蕩模式的多通道同步聯合分析,獲得不同通道的共同模式,確保了本征模態函數分量在數量上和尺度上匹配,解決了多通道信號的模式校準問題。

①采用Hammersley序列采樣法,在n-1維球面上獲得適宜的維空間的方向向量[24];

⑤對于球空間K個方向向量,局部均值可由下式計算

(1)

(2)

式中:q表示IMF的個數,di(t)={di,1(t),di,2(t),…,di,n(t)},ri(t)={ri,1(t),ri,2(t),…,ri,n(t)}分別對應于n元信號的n組IMF分量和n個余量。n元信號每一元變量對應的IMF在各通道中按頻率尺度對齊,形成多元IMF。

2.2 傳遞熵

在信息熵的基礎上[26],Schreiber提出了傳遞熵,用來度量兩個時間序列之間有向信息的傳遞。A到B的傳遞熵實質為A的信息對于B不確定性大小的改變,即A傳遞給B的信息量的大小。

假設存在一組一元EEG信號的時間序列為X={x1,x2,…,xT},一組一元EMG信號的時間序列為Y={y1,y2,…,yT},T為時間序列的長度,對上述兩組時間序列進行MEMD分解,得到與時頻尺度對應的時間序列分別為{Xk}(k=1,2,…,q)和{Yk},其中q代表EEG與EMG信號分解所得的IMF分量的個數。

熵率表示為下一次觀測值X的值所需要的附加信息[27],在給定兩組時間信號下,觀測值X的熵率定義如下:

(3)

(4)

TYk→Xk=h2-h1=

(5)

式中:TYk→Xk表示EMG的第k個分量到EEG的第k個分量的傳遞熵值。同理,從Xk到Yk的傳遞熵為:

(6)

式中:TXk→Yk表示EEG的第k個分量到EMG的第k個分量的傳遞熵值。該傳遞熵值的大小說明了該頻段內皮層肌肉耦合的強弱。

3 實驗結果與分析

本文采用基于多元經驗模態分解-傳遞熵方法對同步采集的EEG和EMG信號進行功能耦合分析。對10位受試者(S1～S10)同步采集的腦電、肌電信號,分別采用EMD、EEMD和MEMD方法進行分解。其中,EEMD方法的參數選擇按照文獻[10]進行,經過多次試驗,添加的白噪聲序列的幅值系數為0.05,算法執行總次數為100次。

圖2所示為受試者S1腦電、肌電(EEG信號C4、左手肱三頭肌、左手橈側腕屈肌)共3個通道的分解結果。

圖2 受試者S1的腦電、肌電信號分解結果

由圖2可見,EMD和EEMD適合處理單通道數據,但是不能同時分析多通道信號。在圖2(a)中可以清晰看出EMD方法對于三路信號的分解所產生的IMF分量完全不相同,由圖2(b)可以看出EEMD分解雖然產生了相同的IMF分量個數,但是其不適合多通道同時分解,因此可能會存在同一尺度IFM分量腦肌電信號頻段不一致的問題。由圖2(c)可以看出,MEMD可以同時分析多通道信號,并且IMF分量的個數相同。

由于本文研究重點為0～75 Hz頻段內的皮層肌肉耦合關系,因此在頻段外IMF分量不再考慮。圖4給出了3種方法對受試者S1腦電信號C4的時頻分析結果,其中奇數列為IMF分量,偶數列為其對應的頻譜圖。由圖3(a)、3(b)可見,EMD和EEMD的分解結果出現了明顯的模態混疊和多個虛假的低頻分量。圖3(c)為MEMD分解結果,該方法可以分解出按先后順序由高到低排列的多元IMF頻率分量,有效抑制了模態混疊現象。

考慮到腦電、肌電耦合分析中存在信息流延遲的問題,因此最終結果分析之前還需要進行傳遞熵預測時間參數的選擇,即式(5)和式(6)中的t值。研究發現在腦電、肌電耦合相關分析中,信息流延遲時間主要在20 ms～30 ms之間[28],傳遞熵峰值出現的地方就是所對應預測時間的最優參數。本文分析了在不同靜態握力下6名健康受試者和4名腦卒中患者的平均傳遞熵值隨預測時間t的變化規律,結果如表1所示。

圖3 受試者S1的腦電信號C4的時頻分解結果

受試者S1S2S3S4S5S6S7S8S9S10L5-FER23202927282920302327L10-FER22202122292823212624L20-FER262828202320————R5-FER22272421212528202429R10-FER20212124222521292920R20-FER222125292926————L5-TB28272723232021212623L10-TB29262420212926252228L20-TB233030303026————R5-TB21202220222626212320R10-TB26232630202524212925R20-TB282923282026————

為了進一步探究EEG和EMG在不同耦合方向(EEG→EMG,EMG→EEG)上、不同頻段間的耦合特征,根據式(5)和式(6)計算出6名受試者腦電、肌電信號不同尺度間的MEMD-TE值。對受試S1～S6的EEG和EMG信號分別進行MEMD分解,得到各個IMF分量的帶寬及其所對應的功能頻段。以受試者S2為例,對5 kg力度下左手抓握動作時所采集的信號進行分解,得到的IMF分量的帶寬和對應頻段結果列于表2。從表2中可以看出,三路信號分解得到的IMF數量相同,尺度大小依次對應,物理意義明確。

圖4 健康受試者S2不同尺度下的TE值

信號分量帶寬范圍/Hz對應功能頻段信號分量帶寬范圍/Hz對應功能頻段信號分量帶寬范圍/Hz對應功能頻段IMF140^75gammaIMF145^70gammaIMF145^70gammaIMF 220^40IMF 220^40IMF 220^40IMF 313^25betaIMF 313^25betaIMF 313^25betaC4IMF 48^15TBIMF 48^15FERIMF 48^15IMF 54^10thetaIMF 54^9thetaIMF 54^9thetaIMF 62^5deltaIMF 62^6deltaIMF 62^6deltaIMF 71^3IMF 71^3IMF 71^3

為進一步比較所有受試者皮層肌肉頻段間耦合的差異性,提取所有受試者EEG→EMG、EMG→EEG方向上各頻段間耦合的MEMD-TE值。圖4給出了受試者S2(健康)和S8(腦卒中)的左手動作時腦電、肌電耦合的MEMD-TE值與EEMD-TE值(其他受試者的結果與此相似)。由圖4和圖5可見,在本文研究頻段內,受試者皮層肌肉耦合及信息傳遞特征具有相似性,皮層肌肉功能耦合是雙向的,這與文獻[28]的研究結果一致,共同體現了特定運動模式下皮層-肌肉的控制環路特征,皮層的運動指令下行到達肌肉,同時也伴隨著來自收縮肌肉的傳入神經反饋過程[29]。同時發現,無論是左手動作還是右手動作,無論是健康人還是病人,在高頻段(40 Hz～75 Hz)腦電對肌電的耦合強度大于肌電對腦電的耦合強度,如圖中的黑色方框內所示。這可能是由于在維持恒定握力輸出過程中,運動皮層需要進行信息整合并將其傳遞至肌肉,體現了感覺反饋與運動控制機理之間的方向性同步震蕩的差異性。圖4 和圖5中虛線為EEMD分解結果,可以看出該方法結果與MEMD分解出的結果曲線類似,但是在高頻段EEMD并不能正確表征腦電、肌電的耦合特性,如圖中紅色框內所示。

圖5 腦卒中受試者S8不同尺度下的TE值

接下來,分別對6名健康人和4名患者(B表示腦卒中患者)的腦電、肌電耦合結果進行統計分析。首先對6名健康人做不同力度的左手動作時腦電(C4)、右手動作時腦電(C3)與肌電(FER,TB)信號在高頻段上的MEMD-TE值取平均和方差;然后對4名患者的數據做相同的處理,所得結果如圖6所示。可以發現在高頻段無論是左手動作還是右手動作,腦電對肌電的MEMD-TE值均大于肌電對腦電的MEMD-TE值,再次驗證了上述文中結果。

為進一步說明不同力度下MEMD-TE關系的差異,將引入一個新的變量MEMD-TE比,指的是某個尺度的MEMD-TE值占整個研究頻譜MEMD-TE值總和的百分比,可用如下公式表示:

(7)

圖6 高頻段上不同受試者不同力度時的腦電、肌電耦合特性

對表3中數據進行顯著性校驗,四組數據的方差齊性檢驗F值依次為0.835、0.529、0.154、0.783,可見大于0.05,說明不同力度下方差基本一致。再進一步進行單因素方差分析,P值依次為0.103、0.128、0.073、0.081,因此認為腦電信號對肌電信號的MEMD-TE值在不同力度下沒有顯著性差別,這與本文圖4和文獻[30]研究結果基本一致。

表3 Gama頻段上腦電信號對肌電信號的MEMD-TE值比

4 結論

本文將多元經驗模態分解運用到腦電、肌電同步耦合分析之中,并且與傳遞熵結合建立了多元經驗模態分解-傳遞熵(MEMD-TE)分析模型,不僅便于提取腦肌電信號的時頻特征,而且可以對不同頻段上皮層肌肉耦合進行分析。將多元經驗模態分解-傳遞熵應用于健康人和病人靜態握力輸出條件下腦電、肌電信號的同步分析,結果表明,靜態握力輸出時,皮層肌肉功能耦合是雙向的,進一步研究發現在高頻段(40 Hz～75 Hz),腦電對肌電的耦合強度大于肌電對腦電的耦合強度,體現了感覺反饋與運動機理之間的方向性同步震蕩的差異性。并且對不同力度下腦肌電耦合情況進行了分析,發現不同力度對于腦電對肌電的傳遞熵沒有顯著影響。因此,本文提出的多元經驗模態分解-傳遞熵方法對于研究不同尺度間腦肌電相互作用和信息傳遞起到了一定作用,可以增進對神經元通信的了解,有助于運動障礙和康復相關的臨床研究。

參考文獻:

[1] Costa M,Goldberger A L,Peng C K. Broken Asymmetry of the Human Heartbeat:Loss of Time Irreversibility in Aging and Disease[J]. Physical Review Letters,2005,95(19):198102(1)-198102(4).

[2] Xu Y,Mcclelland V M,CvetkoviZ,et al. Corticomuscular Coherence with Time Lag with Application to Delay Estimation[J]. IEEE Transactions on Biomedical Engineering,2017,64(3):588-600.

[3] Chiang J,Wang Z J,Mckeown M J. A Multiblock PLS Model of Cortico-Cortical and Corticomuscular Interactions in Parkinson’s Disease[J]. Neuroimage,2012,63(3):1498-1509.

[4] Babiloni C,Vecchio F,Bares M,et al. Functional Coupling Between Anterior Prefrontal Cortex(BA10)and Hand Muscle Contraction During Intentional And Imitative Motor Acts[J]. Neuroimage,2008,39(3):1314-1323.

[5] Joseph J V J. On the Relationship between Corticomuscular(EEG-EMG)Phase Coupling and Muscular Fatigue[D]. Cape Town:University of Cape Town 2015.

[6] Yang Y,Solis-Escalante T,Ruit M V D,et al. Nonlinear Coupling Between Cortical Oscillations and Muscle Activity during Isotonic Wrist Flexion[J]. Frontiers in Computational Neuroscience,2016,10(126):1-11.

[7] Huang N E,Shen Z,Long S R,et al. The Empirical Mode Decomposition Method and the Hilbert Spectrum for Non-Stationary Time Series Analysis[C]//Proceedings of the Royal Society of London,454A(1998)903-995.

[8] Li S,Zhou W,Yuan Q,et al. Feature Extraction and Recognition of Ictal EEG Using EMD and SVM[J]. Computers in Biology and Medicine,2013,43(7):807-816.

[9] 孟宗,戴桂平,劉彬. 基于EMD時頻分析方法的性能研究[J]. 傳感技術學報,2006,19(4):1029-1032.

[10] Wu Z,Huang N E. Ensemble Empirical Mode Decomposition:A Noise-Assisted Data Analysis Method[J]. Advances in Adaptive Data Analysis,2009,1(01):1-41.

[11] G?tz T,Stadler L,Fraunhofer G,et al. A Combined cICA-EEMD Analysis of EEG Recordings from Depressed or Schizophrenic Patients during Olfactory Stimulation[J]. Journal of Neural Engineering,2017,14(1):016011(1)-016011(17).

[12] Wu Z,Huang N E,Chen X. The Multi-Dimensional Ensemble Empirical Mode Decomposition Method[J]. Advances in Adaptive Data Analysis,2009,1(03):339-372.

[13] Dragomiretskiy K,Zosso D. Variational Mode Decomposition[J]. IEEE Transactions on Signal Processing,2014,62(3):531-544.

[14] Rehman N,Mandic D P. Multivariate Empirical Mode Decomposition[J]. Proceedings Mathematical Physical and Engineering Sciences,2010,466(2117):1291-1302.

[15] She Q S,Gan H T,Ma Y L,et al. Scale-Dependent Signal Identification in Low-Dimensional Subspace:Motor Imagery Task Classification[J]. Neural Plasticity,2016,2016(1):1-15.

[16] 韓笑,佘青山,高云園,等. 基于NA-MEMD和互信息的腦電特征提取方法[J]. 傳感技術學報,2016,29(8):1140-1148.

[17] Hao H,Wang H L,Rehman N U. A Joint Framework for Multivariate Signal Denoising Using Multivariate Empirical Mode Decomposition[M]. Elsevier North-Holland,Inc. 2017:263-273.

[18] Murari A,Peluso E,Gelfusa M,et al. Application of Transfer Entropy to Causality Detection and Synchronization Experiments in Tokamaks[J]. Nuclear Fusion,2016,56(2):026006.

[19] Barnett L,Barrett A B,Seth A K. Granger Causality and Transfer Entropy are Equivalent for Gaussian Variables[J]. Physical Review Letters,2009,103(23):238701(1)-238701(4).

[20] 應基光,王行愚. 基于EMD的左右手運動想象腦電信號分析研究[J]. 生物醫學工程學進展,2009,30(3):125-130.

[21] 高云園,任磊磊,張迎春,等. 基于相干性的多頻段腦肌電信號雙向耦合分析[J]. 傳感技術學報,2017,30(10):1465-1470.

[22] Vanessa K,Claudia W,Martin S,et al. Cortico-Muscular Coupling and Motor Performance are Modulated by 20 Hz Transcranial Alternating Current Stimulation(tACS)in Parkinson’s Disease[J]. Frontiers in Human Neuroscience,2014,7(7):1-10.

[23] Brown P. Cortical Drives to Human Muscle:The Piper and Related Rhythms[J]. Progress in Neurobiology,2000,60(1):97-108.

[24] Cui J,Freeden W. Equidistribution on the Sphere[M]. Society for Industrial and Applied Mathematics,1997:595-609.

[25] Rilling G,Flandrin P,Goncalves P. On Empirical Mode Decomposition and Its Algorithms[J]. IEEE-EURASIP Workshop on Nonlinear Signal and Image Processing,2003,3(3):8-11.

[26] Shannon C E A. A Mathematical Theory of Communication. AT&T Tech J[J]. Acm Sigmobile Mobile Computing and Communications Review,2001,5(1):3-55.

[27] Shao S,Guo C,Luk W,et al. Accelerating Transfer Entropy Computation[C]//Proceeding of International Conference on Field-Programmable Technology. IEEE,2015:60-67.

[28] Witham C L,Riddle C N,Baker M R,et al. Contributions of Descending and Ascending Pathways to Corticomuscular Coherence in Humans[J]. The Journal of physiology,2011,589(15):3789-3800.

[29] Schelter B,Timmer J,Eichler M. Assessing the Strength of Directed Influences Among Neural Signals Using Renormalized Partial Directed Coherence[J]. Journal of neuroscience methods,2009,179(1):121-130.

[30] 榮瑤. 大腦皮層與上肢肌肉間功能耦合與信息傳輸的研究[D]. 北京:北京工業大學,2013.