基于時(shí)變copula互信息的肌間耦合分析

王洪安 佘青山? 馬玉良 孔萬(wàn)增 田玉平

1(杭州電子科技大學(xué)自動(dòng)化學(xué)院(人工智能學(xué)院)，杭州 310018)

2(浙江省腦機(jī)協(xié)同智能重點(diǎn)實(shí)驗(yàn)室，杭州 310018)

引言

在人的自主運(yùn)動(dòng)過(guò)程中，大腦運(yùn)動(dòng)皮層以神經(jīng)元的振蕩和同步放電的方式發(fā)出運(yùn)動(dòng)控制指令，經(jīng)下行神經(jīng)通道傳遞至脊髓，通過(guò)募集一組特定的肌肉協(xié)同模塊，實(shí)現(xiàn)對(duì)復(fù)雜運(yùn)動(dòng)任務(wù)的最優(yōu)控制[1]。腦電(electroencephalography， EEG)信號(hào)和身體對(duì)側(cè)的表面肌電(surface electromyogram， sEMG)信號(hào)分別反映了運(yùn)動(dòng)控制信息和功能響應(yīng)信息，腦肌電信號(hào)之間帶節(jié)律的同步特征揭示了運(yùn)動(dòng)過(guò)程中大腦皮層和相應(yīng)肌肉之間的功能耦合連接[2]。

研究表明，在中樞神經(jīng)系統(tǒng)的功能調(diào)節(jié)和反饋控制下，與運(yùn)動(dòng)相關(guān)的肌肉之間存在不同時(shí)空層次的信息交互和耦合關(guān)系[3]。 這種特殊的耦合關(guān)系不僅反映了運(yùn)動(dòng)過(guò)程中肌肉間的相互關(guān)聯(lián)和相互作用，還包含了神經(jīng)系統(tǒng)對(duì)運(yùn)動(dòng)肌肉的控制支配作用[4]。 肌間耦合是從皮層肌肉耦合研究中發(fā)現(xiàn)并引申而來(lái)的，源于特定的運(yùn)動(dòng)任務(wù)下，相關(guān)肌肉的運(yùn)動(dòng)神經(jīng)元共享相同的皮質(zhì)脊髓驅(qū)動(dòng)[5]。 肌間耦合存在頻段顯著特征，主要表現(xiàn)為alpha (8～15 Hz)頻段的肌間耦合受脊髓神經(jīng)的控制，影響肌肉的非自主收縮和運(yùn)動(dòng)姿勢(shì)的維持；beta (15～30 Hz)頻段的肌間耦合代表了從初級(jí)運(yùn)動(dòng)皮層到運(yùn)動(dòng)神經(jīng)元的信息傳遞過(guò)程，與靜態(tài)力的輸出有關(guān)；gamma(30～60 Hz)頻段的肌間耦合受大腦皮質(zhì)神經(jīng)的控制，體現(xiàn)為認(rèn)知過(guò)程中注意力集中引發(fā)的強(qiáng)烈緊張性收縮[6]。 不同頻段的肌間耦合特性，為理解運(yùn)動(dòng)控制的組織和協(xié)調(diào)機(jī)制提供了理論依據(jù)。 當(dāng)前，肌間耦合關(guān)系的研究集中在運(yùn)動(dòng)醫(yī)學(xué)、康復(fù)工程等領(lǐng)域，已成為運(yùn)動(dòng)神經(jīng)科學(xué)關(guān)注的熱點(diǎn)問(wèn)題。

近年來(lái)，肌間耦合分析方法層出不窮，相干性(coherence)因其算法原理簡(jiǎn)單、易于實(shí)現(xiàn)，被廣泛應(yīng)用于肌間耦合分析。 謝平等[7]利用相干性，對(duì)比分析了中風(fēng)患者在運(yùn)動(dòng)過(guò)程中健側(cè)、患側(cè)上肢拮抗肌間的相干性特征，發(fā)現(xiàn)患側(cè)在beta 頻段的肌間相干性相對(duì)于健側(cè)存在明顯缺失。 然而，肌間耦合關(guān)系也包含非線性成分，基于傅里葉變換的相干性僅能描述線性耦合關(guān)系，并且Faes 等[8]指出，相干性所測(cè)得的耦合關(guān)系包含了直接和間接的影響，會(huì)過(guò)度估計(jì)肌間耦合強(qiáng)度。 概率論和信息論中的互信息(mutual information， MI)是傳統(tǒng)相關(guān)系數(shù)的非線性擴(kuò)展，不依賴既定模型，能夠度量?jī)蓚€(gè)或多個(gè)隨機(jī)變量之間的線性或非線性依賴程度，在很多方面被認(rèn)為是較完美的關(guān)聯(lián)統(tǒng)計(jì)工具。 然而，互信息的估計(jì)非常依賴概率密度函數(shù)的精確表示，這在實(shí)際應(yīng)用中通常十分困難。 Ma 等[9]證明了負(fù)的copula熵可以等價(jià)為互信息，通過(guò)copula 函數(shù)估計(jì)互信息，不僅能有效避免對(duì)聯(lián)合密度函數(shù)的估計(jì)，而且計(jì)算復(fù)雜度低，這為理解和估計(jì)互信息提供了一條新的思考途徑。

Sklar 定理指出，一個(gè)N維聯(lián)合分布函數(shù)可以分解為N個(gè)邊際分布函數(shù)和一個(gè)copula 函數(shù)[10]。 該copula 函數(shù)完整地描述了變量間的相關(guān)性，被譯為“連接函數(shù)”或“相依函數(shù)”，在生物統(tǒng)計(jì)學(xué)、金融風(fēng)險(xiǎn)度量等許多領(lǐng)域都有著重要貢獻(xiàn)。 Fisher 等[11]指出，copula 函數(shù)的優(yōu)勢(shì)突出表現(xiàn)在兩個(gè)方面:一是可用來(lái)構(gòu)建復(fù)雜的高維概率密度；二是copula 函數(shù)可以捕捉線性、非線性、對(duì)稱以及非對(duì)稱相關(guān)性。Copula 函數(shù)種類繁多，主要有橢圓copula 函數(shù)和阿基米德copula 函數(shù)；不同copula 函數(shù)結(jié)構(gòu)不同，在實(shí)際應(yīng)用時(shí)需要與數(shù)據(jù)的分布類型很好匹配。 現(xiàn)階段大部分研究都是基于參數(shù)靜態(tài)的常相關(guān)copula函數(shù)，這種函數(shù)被認(rèn)為不能反映相關(guān)結(jié)構(gòu)的動(dòng)態(tài)變化，具有一定的局限性[12]。 其實(shí)，早在2001年P(guān)atton 就已提出參數(shù)時(shí)變的copula 函數(shù)[13]，其構(gòu)造原理類似于一個(gè)ARMA (1， 10)過(guò)程。 近幾年，時(shí)變copula 函數(shù)因其相關(guān)參數(shù)時(shí)變的特點(diǎn)，在氣象、經(jīng)濟(jì)等領(lǐng)域取得了較好發(fā)展[14-16]。

針對(duì)現(xiàn)有的肌間耦合分析方法存在的不足，筆者提出了一種利用時(shí)變copula 函數(shù)估計(jì)互信息的新方法，定量研究了不同特征頻段的肌間耦合關(guān)系，并在Ninapro DB4 公開數(shù)據(jù)集上進(jìn)行了測(cè)試和分析，為進(jìn)一步探索上肢運(yùn)動(dòng)功能障礙的產(chǎn)生機(jī)理以及運(yùn)動(dòng)功能康復(fù)評(píng)估，提供新的研究方法和理論依據(jù)。

1 材料和方法

1.1 時(shí)變copula 函數(shù)

二元分布的Sklar 定理[10]指出，若二維隨機(jī)變量x、y的聯(lián)合分布函數(shù)為F，邊際分布函數(shù)分別為u和v，則存在一個(gè)copula 函數(shù)C，使得

若u、v連續(xù)，則C是唯一的。 這個(gè)copula 函數(shù)描述了變量間的相關(guān)結(jié)構(gòu)。 根據(jù)式(1)，可得x、y的聯(lián)合密度為

式中，c(u，v) 表示copula 密度函數(shù)，f(x) 與f(y)分別表示x和y的概率密度函數(shù)。

時(shí)變copula 函數(shù)也稱為條件copula 函數(shù)，其數(shù)學(xué)形式和靜態(tài)copula 函數(shù)的形式相同，只有copula函數(shù)的參數(shù)發(fā)生變化[13]。 下面介紹兩種常用的時(shí)變copula 函數(shù)。

1.1.1 時(shí)變normal copula 函數(shù)

時(shí)變normal copula 函數(shù)的分布函數(shù)表達(dá)式為

式中，Φ－1(·) 為標(biāo)準(zhǔn)正態(tài)分布的逆函數(shù)，ρt為時(shí)變相關(guān)系數(shù)。

時(shí)變normal copula 函數(shù)對(duì)尾部相關(guān)變化不敏感，即在時(shí)變相關(guān)系數(shù)值域(－1， 1)區(qū)間內(nèi)，上、下尾部的相關(guān)系數(shù)均為0。

1.1.2 時(shí)變SJC copula 函數(shù)

時(shí)變SJC (symmetrized Joe-Clayton) copula 函數(shù)由JC (Joe-Clayton) copula 函數(shù)變換而來(lái)，其分布函數(shù)表達(dá)式為

其中

式中，τU，t和τl，t分別為上尾相關(guān)系數(shù)和下尾相關(guān)系數(shù)。

時(shí)變SJC copula 函數(shù)對(duì)上、下尾相關(guān)變化均敏感，在一定程度上，時(shí)變SJC copula 函數(shù)的適用范圍更廣。

1.2 時(shí)變copula 互信息

互信息是信息論中的一個(gè)基本概念。 與傳統(tǒng)的獨(dú)立性統(tǒng)計(jì)度量方法不同，如相關(guān)性系數(shù)(二階統(tǒng)計(jì)量)和高階統(tǒng)計(jì)量，互信息是全階次、非線性相關(guān)的度量。 Ma 等[9]證明了互信息本質(zhì)上是一種copula 熵，即

式中，Hc(u，v) 表示隨機(jī)變量x和y之間的copula 熵。

由式(6)推導(dǎo)的互信息一般稱為copula 互信息，它具有連續(xù)、對(duì)稱、可加和對(duì)邊際分布單調(diào)變換不敏感等特質(zhì)[10]。 Copula 互信息的值越大，表示變量間相互依賴程度越大。

計(jì)算copula 互信息的關(guān)鍵，在于邊際分布函數(shù)的擬合、copula 函數(shù)的選擇和copula 熵的估計(jì)。 其中，copula 函數(shù)的選擇尤為重要，不同的copula 函數(shù)具有不同的表現(xiàn)形式，變量間刻畫的相依結(jié)構(gòu)大相徑庭，所得的copula 互信息也大小不一。 據(jù)此，為更加準(zhǔn)確、合理地權(quán)衡變量間的耦合關(guān)系，本研究提出一種新的時(shí)變copula 互信息方法。

首先利用樣本的經(jīng)驗(yàn)分布函數(shù)u＾、v＾作為總體邊際分布函數(shù)的近似，然后根據(jù)典型極大似然估計(jì)法(canonical maximum likelihood， CML)，估計(jì)時(shí)變copula 函數(shù)中的時(shí)變相關(guān)參數(shù)θ＾n，并將時(shí)變copula密度函數(shù)ctv代入式(6)中， 進(jìn)而采用蒙特卡洛法計(jì)算copula 熵，最后得到時(shí)變copula 互信息。 具體推導(dǎo)過(guò)程如下:

式中，I為示性函數(shù)，E[·] 表示求期望，時(shí)變copula互信息單位為比特(bit)。

時(shí)變copula 函數(shù)可通過(guò)求偏導(dǎo)得到時(shí)變copula密度函數(shù)，以時(shí)變normal copula 函數(shù)為例，其時(shí)變copula 密度函數(shù)為

1.3 實(shí)驗(yàn)數(shù)據(jù)與分析

使用Ninapro DB4 公開數(shù)據(jù)庫(kù)提供的sEMG 數(shù)據(jù)集(https:/ /doi.org/10.528 1/zenodo.1000138)[17]。 這些數(shù)據(jù)來(lái)自10 名健康受試者(S1～S10)，其詳細(xì)情況如表1 所示。 每位受試者要求完成52 種手腕、手和手指動(dòng)作(外加休息)，每種動(dòng)作重復(fù)6 次，每次執(zhí)行5 s，休息3 s。

表1 受試者的相關(guān)信息Tab.1 Demographic information of subjects

數(shù)據(jù)集中的stimulus 時(shí)間序列用于設(shè)置每種運(yùn)動(dòng)重復(fù)片段的起點(diǎn)和終點(diǎn)。 每段由MyoBock 13E200－50 系統(tǒng)(Cometa Wave Plus ＋Dormo)記錄的12 通道sEMG 信號(hào)組成，采樣頻率為2 000 Hz，利用均方根進(jìn)行校正。 12 個(gè)無(wú)線有源單差分電極的位置分別為:8 個(gè)電極(Channel 1～8)均勻分布于前臂周圍，與肱橈關(guān)節(jié)對(duì)應(yīng)；2 個(gè)電極(Channel 9～10)，放置在指屈肌和指伸肌的主要活動(dòng)點(diǎn)上；2 個(gè)電極(Channel 11～12)，放置在肱二頭肌(BB)和肱三頭肌(TB)的主要活動(dòng)點(diǎn)上。

運(yùn)動(dòng)任務(wù)主要集中在腕部運(yùn)動(dòng)的兩個(gè)特定的自由度上:腕屈(wrist flexion， WF)和腕展(wrist extension， WE)，如圖1 所示。 選取BB 和TB 這兩塊肌肉作為研究對(duì)象。 鑒于sEMG 信號(hào)頻域特征突出，在不同特征頻段上肌間耦合特性存在明顯差異，因此本研究利用切比雪夫Ⅱ型帶通濾波器，將BB 和TB 的sEMG 信號(hào)劃分到不同頻段上，即theta(4～8 Hz)、alpha (8～15 Hz)、beta (15～30 Hz)和gamma (30～60 Hz)，并且參照文獻(xiàn)[18]的做法，將gamma 頻段進(jìn)一步細(xì)分為低gamma (30～45 Hz)和高gamma (45～60 Hz)兩部分進(jìn)行實(shí)驗(yàn)對(duì)比分析。

圖1 腕屈和腕展Fig.1 Wrist flexion (WF) and wrist extension (WE)

1.4 方法驗(yàn)證和統(tǒng)計(jì)分析

為體現(xiàn)時(shí)變copula 互信息在衡量肌間耦合強(qiáng)度方面的能力，首先采用經(jīng)驗(yàn)分布函數(shù)逼近變量的真實(shí)分布，并比較了靜態(tài)normal copula、靜態(tài)SJC copula、時(shí)變normal copula 和時(shí)變SJC copula 等4 種不同的copula 函數(shù)對(duì)不同頻段肌間相依關(guān)系的擬合情況，這些函數(shù)的參數(shù)通過(guò)極大似然估計(jì)法進(jìn)行確定；再根據(jù)赤池信息準(zhǔn)則(Akaike information criterion， AIC)，判斷不同copula 函數(shù)的擬合優(yōu)度[19]，AIC 值越小表明相應(yīng)的copula 函數(shù)對(duì)樣本數(shù)據(jù)的擬合程度越高。 其次，為驗(yàn)證時(shí)變copula 互信息度量肌間耦合強(qiáng)度的有效性，選取了上述兩種靜態(tài)copula 函數(shù)和兩種時(shí)變copula 函數(shù)去估計(jì)互信息，同時(shí)進(jìn)行對(duì)比分析，按照AIC 最小原則，從中挑選出擬合效果最佳的時(shí)變copula 函數(shù)。 最后，利用最適合腕屈(WF)和腕展(WE)的時(shí)變copula 互信息，刻畫不同動(dòng)作在不同頻段上的肌間耦合強(qiáng)度，并使用獨(dú)立樣本t檢驗(yàn)法統(tǒng)計(jì)校驗(yàn)，進(jìn)一步探討其與被試的生理參數(shù)之間的相關(guān)性，顯著性水平均設(shè)為0.05。

2 結(jié)果

2.1 靜態(tài)copula 和時(shí)變copula 的性能比較

表2 給出了受試者S8 在第3 次腕屈(WF)、腕展(WE)時(shí)，靜態(tài)copula 函數(shù)和時(shí)變copula 函數(shù)在5個(gè)特征頻段的AIC 值。 由表2 可見，腕屈運(yùn)動(dòng)過(guò)程中，除了高gamma 頻段以外的其他4 個(gè)特征頻段，normal copula 的AIC 值都比SJC copula 的要小，在theta 和高gamma 頻段，時(shí)變SJC copula 的AIC 值要低于時(shí)變normal copula 的相應(yīng)值，而在alpha、beta和低gamma 頻段，時(shí)變normal copula 的AIC 值要低于時(shí)變SJC copula 的相應(yīng)值；腕展運(yùn)動(dòng)過(guò)程中，在alpha 和beta 頻段，normal copula 的AIC 值要比SJC copula 的相應(yīng)值小；而在theta 和gamma 頻段，SJC copula 的AIC 值要低于normal copula 的相應(yīng)值。同靜態(tài)copula 函數(shù)的AIC 結(jié)果一致，時(shí)變copula函數(shù)在相應(yīng)頻段亦有此特點(diǎn)。 相比而言，時(shí)變copula 函數(shù)的AIC 值要明顯低于靜態(tài)copula 函數(shù)的相應(yīng)值。 這說(shuō)明，在腕屈和腕展的運(yùn)動(dòng)過(guò)程中，時(shí)變copula 函數(shù)對(duì)肌間相依結(jié)構(gòu)的擬合優(yōu)度要比靜態(tài)copula 函數(shù)的相應(yīng)值更高，時(shí)變copula 函數(shù)對(duì)肌間相依結(jié)構(gòu)的刻畫要比靜態(tài)copula 函數(shù)相應(yīng)結(jié)構(gòu)的刻畫更合適。

表2 靜態(tài)和時(shí)變copula 函數(shù)的AIC 比較Tab.2 Comparison of AIC using static and time-varying copula functions

圖2 給出了靜態(tài)和時(shí)變normal copula 函數(shù)在5個(gè)特征頻段的靜態(tài)相關(guān)系數(shù)和時(shí)變相關(guān)系數(shù)對(duì)比結(jié)果。 由圖2(a)可以看出，在腕屈運(yùn)動(dòng)過(guò)程中的不同特征頻段上，相比幾乎不變且趨近于0 的靜態(tài)normal copula 相關(guān)系數(shù)，時(shí)變normal copula 相關(guān)系數(shù)能夠反映BB 與TB 隨時(shí)間變化的動(dòng)態(tài)相關(guān)特征，且頻段越高時(shí)變normal copula 相關(guān)系數(shù)的變化越劇烈，在總體上表現(xiàn)出正相關(guān)的特點(diǎn)；由圖2(b)可以看出，在腕展運(yùn)動(dòng)過(guò)程中的不同特征頻段上，靜態(tài)normal copula 相關(guān)系數(shù)很低，集中在0 附近，而時(shí)變normal copula 相關(guān)系數(shù)則不盡相同，表現(xiàn)為一種低頻段高相關(guān)、高頻段低相關(guān)的特點(diǎn)，從平均結(jié)果來(lái)看，相關(guān)程度并不高。

圖2 靜態(tài)和時(shí)變normal copula 函數(shù)的相關(guān)系數(shù)對(duì)比。 (a) 腕屈；(b) 腕展Fig.2 Comparison of correlation coefficients of static and time-varying normal copula functions.(a)WF；(b)WE

圖3 給出了靜態(tài)和時(shí)變SJC copula 函數(shù)在5 個(gè)特征頻段的靜態(tài)相關(guān)系數(shù)和時(shí)變相關(guān)系數(shù)的對(duì)比結(jié)果。 由圖3(a)可以看出，在腕屈運(yùn)動(dòng)過(guò)程中的不同特征頻段上，靜態(tài)SJC copula 上、下尾相關(guān)系數(shù)都趨近于0，而時(shí)變SJC copula 上、下尾相關(guān)系數(shù)則明顯不同。 較為特別的是，在alpha 和beta 頻段，BB和TB 的靜態(tài)SJC copula 上尾相關(guān)系數(shù)與時(shí)變SJC copula 上尾相關(guān)系數(shù)保持一致，均為0；在theta 和低gamma 頻段，BB 和TB 的靜態(tài)SJC copula 下尾相關(guān)系數(shù)與時(shí)變SJC copula 下尾相關(guān)系數(shù)保持一致，均為0；而在其他特征頻段，BB 和TB 的時(shí)變SJC copula 上、下尾相關(guān)系數(shù)都比較高，平均值大于0.5，表現(xiàn)為較強(qiáng)的正相關(guān)，且頻段越高時(shí)變SJC copula 上、下尾相關(guān)系數(shù)的變化越劇烈。 由圖3(b)可以看出，在腕展運(yùn)動(dòng)過(guò)程中不同特征頻段上，靜態(tài)SJC copula 上、下尾相關(guān)系數(shù)在不同特征頻段上都幾乎為0，而時(shí)變SJC copula 上、下尾相關(guān)系數(shù)則截然不同(除alpha 和beta 頻段的時(shí)變SJC copula下尾相關(guān)系數(shù))，總體上表現(xiàn)為較高的正相關(guān)性，尤其是在低頻段；對(duì)照來(lái)看，高頻段的時(shí)變SJC copula上、下尾相關(guān)系數(shù)存在局部峰值(或局部低相關(guān))特點(diǎn)，平均相關(guān)程度不如低頻段。 圖2、3 說(shuō)明，在腕屈和腕展運(yùn)動(dòng)過(guò)程中，BB 與TB 在不同特征頻段上存在不同的時(shí)變相關(guān)性，頻段差異明顯，靜態(tài)copula相關(guān)系數(shù)容易錯(cuò)誤地低估肌間相關(guān)關(guān)系。

圖3 靜態(tài)和時(shí)變SJC copula 函數(shù)的上、下尾部相關(guān)系數(shù)對(duì)比。 (a) 腕屈；(b) 腕展Fig.3 Comparison of upper tail correlation coefficients and lower tail correlation coefficients of static and timevarying SJC copula functions.(a) WF；(b) WE

2.2 基于時(shí)變copula 互信息的肌間耦合強(qiáng)度分析

2.2.1 靜態(tài)與時(shí)變copula 互信息的性能比較

仍以受試者S8 為例進(jìn)行詳細(xì)說(shuō)明。 圖4 給出了BB 與TB 在5 個(gè)特征頻段上，靜態(tài)和時(shí)變copula互信息的估計(jì)結(jié)果。 在腕屈、腕伸運(yùn)動(dòng)過(guò)程中，BB和TB 在5 個(gè)特征頻段上，時(shí)變copula 互信息均要明顯高于靜態(tài)copula 互信息。 具體表現(xiàn)為:在theta、alpha、beta 和 低gamma 頻 段，時(shí) 變normal copula 互信息要高于其他copula 互信息；而在高gamma 頻段，時(shí)變SJC copula 互信息要高于其他copula 互信息。 隨著頻段由高到低，時(shí)變copula 互信息呈現(xiàn)逐漸上升的規(guī)律，而值較低的靜態(tài)copula互信息無(wú)此跡象。 這體現(xiàn)出不同類型的copula 函數(shù)估計(jì)出的互信息完全不同，靜態(tài)copula 互信息可能無(wú)法正確挖掘出肌間耦合強(qiáng)度關(guān)系。

圖4 靜態(tài)和時(shí)變copula 互信息對(duì)比。 (a) 腕屈；(b) 腕展Fig.4 Comparison of static and time-varying copula MI.(a) WF； (b) WE

為了比較靜態(tài)和時(shí)變copula 互信息的性能，表3 給出了所有受試者根據(jù)AIC 最小原則在每一次腕部動(dòng)作、每一個(gè)特征頻段上所選擇的copula 函數(shù)類型的累計(jì)情況(總共10 名受試者×2 種腕部動(dòng)作×6 次重復(fù)×5 個(gè)特征頻段＝600)。 所有受試者都選擇了時(shí)變copula 函數(shù)作為最佳的連接函數(shù)，其中時(shí)變normal copula 占了65.67%，時(shí)變SJC copula 占了34.33%。 這說(shuō)明:時(shí)變copula 函數(shù)比靜態(tài)copula函數(shù)更適合描述肌間耦合關(guān)系，由時(shí)變copula 函數(shù)得出的互信息更適合度量肌間耦合強(qiáng)度的大小；但是，單一的時(shí)變copula 函數(shù)難以描述復(fù)雜的肌間耦合關(guān)系。 結(jié)合圖4 的對(duì)比情況可知，在一定頻段范圍內(nèi)，相較于時(shí)變SJC copula 互信息，時(shí)變normal copula 互信息更適合度量肌間耦合強(qiáng)度的大小。

表3 所有受試者關(guān)于copula 函數(shù)的選擇情況Tab.3 Selection results of copula functions for all subjects

2.2.2 不同頻段和動(dòng)作下時(shí)變copula 互信息對(duì)比

為便于后續(xù)分析，這是將最佳的時(shí)變copula 互信息用于度量肌間耦合強(qiáng)度。 圖5(a)、(b)分別給出了所有受試者在腕屈、腕展時(shí)，BB 與TB 在5 個(gè)特征頻段上時(shí)變copula 互信息的對(duì)比結(jié)果。 其中，同時(shí)標(biāo)出了獨(dú)立樣本T檢驗(yàn)后的P值，實(shí)心小圓點(diǎn)表示受試者樣本。 從圖5 的箱型圖的分布情況來(lái)看，兩種腕部動(dòng)作下，不同特征頻段間的時(shí)變copula互信息都表現(xiàn)出完全一致的規(guī)律，即:theta 頻段的時(shí)變copula 互信息要顯著高于alpha 頻段的值(P<0.05)，alpha 頻段的時(shí)變copula 互信息要顯著高于beta 頻段的值(P<0.000 1)，beta 頻段的時(shí)變copula互信息顯著高于低gamma 頻段的值(P<0.000 1)，低gamma 頻段的時(shí)變copula 互信息要顯著高于高gamma 頻段的值(P<0.01)。 這說(shuō)明，腕屈和腕展運(yùn)動(dòng)過(guò)程中，肌間耦合強(qiáng)度存在明顯的頻段差異。

圖5 不同頻段上時(shí)變copula 互信息的對(duì)比(實(shí)心小圓點(diǎn)表 示 受試者 樣 本；?:P<0.05，??:P<0.01，???:P<0.001，????:P<0.000 1)。 (a) 腕屈；(b) 腕展Fig.5 Comparison of time-varying copula MI in different frequency bands (solid dots represent the subject samples.?: P < 0.05，??: P < 0.01，???: P <0.001，????: P<0.000 1).(a) WF； (b) WE

圖6 給出了所有受試者兩種腕部動(dòng)作過(guò)程中，BB 與TB 在5 個(gè)特征頻段上的時(shí)變copula 互信息的對(duì)比結(jié)果。 可以看出，除了theta 頻段外，在其他特征頻段上，腕屈和腕展動(dòng)作間的時(shí)變copula 互信息均不具有顯著性差異(P>0.05)。 這說(shuō)明，BB 與TB 在alpha、beta、低gamma 和高gamma 頻段上， 肌間耦合強(qiáng)度基本維持穩(wěn)定，不隨腕部運(yùn)動(dòng)方式的變化而改變。 而theta 頻段上的差異(P＝0.049<0.05)，能夠?qū)⑼笄c腕展區(qū)別開來(lái)，說(shuō)明肌間耦合強(qiáng)度在這一頻段與運(yùn)動(dòng)的執(zhí)行方式相關(guān)聯(lián)。 更加具體地，表4 給出了所有受試者基于時(shí)變copula 互信息的肌間耦合強(qiáng)度大小，以供參考。

表4 肌間耦合強(qiáng)度(平均值±標(biāo)準(zhǔn)差，bit)Tab.4 Intermuscular coupling strength (mean ± standard deviation， bit)

圖6 不同頻段上腕屈和腕展的時(shí)變copula 互信息對(duì)比Fig.6 Comparison of time-varying copula MI in different bands during WF and WE movements

2.2.3 時(shí)變copula 互信息與生理指標(biāo)的相關(guān)性分析

進(jìn)一步，還分析了所有受試者腕屈、腕展動(dòng)作過(guò)程中，BB 與TB 在5 個(gè)特征頻段上平均時(shí)變copula 互信息與受試者的年齡、身高和體重之間的相關(guān)性，結(jié)果如圖7 所示。 圖中擬合好的平滑直線表示平均時(shí)變copula 互信息與受試者的年齡、身高和體重的相對(duì)變化趨勢(shì)。 從圖7 可以看出，各個(gè)特征頻段上的平均時(shí)變copula 互信息與受試者的年齡、身高和體重均無(wú)明顯的相關(guān)關(guān)系，在經(jīng)皮爾森相關(guān)性檢驗(yàn)后發(fā)現(xiàn)，P值均大于0.05，這說(shuō)明腕屈、腕展運(yùn)動(dòng)過(guò)程中，肌間耦合強(qiáng)度不受個(gè)體生理定量指標(biāo)的影響。

圖7 時(shí)變copula 互信息與受試者年齡、身高和體重之間的相關(guān)性。 (a) 腕屈；(b) 腕展Fig.7 Correlation between time-varying copula MI and age， height， weight of subjects.(a) WF；(b) WE

3 討論

肌間耦合是指運(yùn)動(dòng)過(guò)程中肌肉間的相互作用，在一定程度上反映了運(yùn)動(dòng)中樞神經(jīng)系統(tǒng)對(duì)相關(guān)肌肉的控制支配作用。 肌間耦合強(qiáng)度體現(xiàn)為肌肉的激活程度和協(xié)調(diào)程度，可作為運(yùn)動(dòng)功能康復(fù)評(píng)估的生理指標(biāo)[20]。 研究運(yùn)動(dòng)過(guò)程中肌間耦合特性，對(duì)于理解神經(jīng)運(yùn)動(dòng)控制策略以及運(yùn)動(dòng)功能障礙機(jī)制具有十分重要的意義。

然而，現(xiàn)有的肌間耦合分析方法存在諸多不足，肌間耦合關(guān)系尚未被很好理解。 鑒于目前copula 函數(shù)種類眾多，任何一種都很難完全刻畫出肌間復(fù)雜的相關(guān)結(jié)構(gòu)，所以本研究將copula 理論中別具特色的時(shí)變copula 函數(shù)引入到肌間耦合分析之中。 時(shí)變copula 函數(shù)不僅繼承了靜態(tài)copula 函數(shù)所具有的特點(diǎn)，而且還因其參數(shù)時(shí)變而更具靈活性和優(yōu)越性，對(duì)相關(guān)結(jié)構(gòu)的描述也更加準(zhǔn)確。 此外，不同于由靜態(tài)copula 函數(shù)導(dǎo)出的相關(guān)測(cè)度在時(shí)間維度上保持恒定不變，時(shí)變copula 函數(shù)能夠反映相關(guān)結(jié)構(gòu)的動(dòng)態(tài)演化過(guò)程，更加符合實(shí)際情況。

在時(shí)變copula 函數(shù)的基礎(chǔ)上，本研究通過(guò)與信息熵理論相結(jié)合，提出了時(shí)變copula 互信息方法，并給出了切實(shí)可行的計(jì)算步驟。 時(shí)變copula 互信息具有多元、非負(fù)、對(duì)稱和對(duì)單調(diào)變換不敏感等特點(diǎn)。 與其他肌間耦合分析方法(如時(shí)域的相關(guān)系數(shù)或頻域的相干性)相比，時(shí)變copula 互信息不對(duì)數(shù)據(jù)分布做任何先驗(yàn)假設(shè)，能更好地度量肌間線性或非線性耦合強(qiáng)度；與依賴概率密度函數(shù)定義的互信息相比，時(shí)變copula 互信息避免了對(duì)聯(lián)合概率密度函數(shù)的估計(jì)，使估計(jì)的結(jié)果更精確，拓展了互信息的適用范圍。 同時(shí)，傳統(tǒng)的經(jīng)驗(yàn)分布函數(shù)具有良好的統(tǒng)計(jì)性質(zhì)，當(dāng)樣本量足夠大時(shí)，能夠很好地?cái)M合總體分布，估計(jì)出的邊際分布函數(shù)可以減少因平滑假設(shè)所帶來(lái)的誤差[21]。 典型極大似然估計(jì)可以不用估計(jì)邊際分布中的參數(shù)，只需估計(jì)時(shí)變copula 函數(shù)中的時(shí)變參數(shù)，在邊際分布未知或擬合效果欠佳的條件下，是一種不錯(cuò)的保守估計(jì)方法[22]。

腕屈和腕展是日常生活中手腕活動(dòng)最常見的兩種方式。 本研究基于時(shí)變copula 和時(shí)變copula互信息，深入研究了腕屈、腕展運(yùn)動(dòng)過(guò)程中，BB 和TB 在theta、alpha、beta、低gamma 和高gamma 共5個(gè)特征頻段的耦合情況。 實(shí)驗(yàn)結(jié)果顯示，在不同特征頻段上，無(wú)論是腕屈運(yùn)動(dòng)還是腕伸運(yùn)動(dòng)，時(shí)變copula 函數(shù)的擬合優(yōu)度都要明顯高于靜態(tài)copula 函數(shù)的擬合優(yōu)度，時(shí)變copula 函數(shù)對(duì)肌間相關(guān)結(jié)構(gòu)的刻畫更充分(見表2)；靜態(tài)copula 相關(guān)系數(shù)在不同特征頻段上均接近于0，而時(shí)變copula 相關(guān)系數(shù)顯示在不同特征頻段上，肌間存在波動(dòng)程度不同的時(shí)變相關(guān)性，隨著頻段的上升，時(shí)變相關(guān)系數(shù)變化越劇烈，上、下尾時(shí)變相關(guān)系數(shù)與頻段分布存在一定的聯(lián)系，總體上呈現(xiàn)為正相關(guān)(見圖3)；在度量肌間耦合強(qiáng)度時(shí)，靜態(tài)copula 互信息會(huì)造成一定的低估(<0.05 bit)，而時(shí)變copula 互信息的估計(jì)結(jié)果更加準(zhǔn)確，多數(shù)情況下，時(shí)變normal copula 互信息更加契合肌間耦合關(guān)系(見圖4 和表3)。 重要的是，本研究發(fā)現(xiàn)，在腕屈、腕展運(yùn)動(dòng)過(guò)程中，不同特征頻間肌間耦合強(qiáng)度存在顯著的統(tǒng)計(jì)學(xué)差異(P<0.01，見圖5 和表4)，表現(xiàn)為:theta 頻段的肌間耦合強(qiáng)度要顯著高于alpha 頻段的值，alpha 頻段的肌間耦合強(qiáng)度要顯著高于beta 頻段的值，beta 頻段的肌間耦合強(qiáng)度要顯著高于低gamma 頻段的值，低gamma 頻段的肌間耦合強(qiáng)度要顯著高于高gamma 頻段的值，這或許表明，神經(jīng)肌肉系統(tǒng)可能被組織成一個(gè)多重平行和分層控制結(jié)構(gòu)，不同的通路支持不同頻率波段的共同輸入，這與Boonstra 等使用相干性分析得出肌間耦合網(wǎng)絡(luò)在不同頻段上存在廣泛連通性的結(jié)論[23]基本一致。 此外，兩種腕部動(dòng)作相差無(wú)幾，僅在theta 頻段發(fā)現(xiàn)腕屈過(guò)程的肌間耦合強(qiáng)度要顯著低于腕展過(guò)程(P<0.05，見圖6)。 這可能與中樞神經(jīng)系統(tǒng)的模塊化控制策略有關(guān)，即腕屈和腕伸動(dòng)作下都具有共享的協(xié)同模塊和任務(wù)特殊的協(xié)同模塊，肌肉活動(dòng)可表現(xiàn)為sEMG 信號(hào)在不同頻率范圍的振蕩協(xié)同作用[24]。 本研究還進(jìn)一步發(fā)現(xiàn)，肌間耦合強(qiáng)度與個(gè)體生理指標(biāo)之間并無(wú)明顯的相關(guān)關(guān)系(見圖7)，提示不同被試的頻域指標(biāo)不具有顯著的個(gè)體差異性，這與Jaiser 等的發(fā)現(xiàn)[25]基本吻合。

當(dāng)然，本研究仍存在一些不足之處。 從理論上講，很難找到一種完美的copula 函數(shù)形式能夠全面地描述神經(jīng)生理信號(hào)之間的耦合關(guān)系，時(shí)變copula函數(shù)自身也存在一定的局限性。 最近，Chang 等提出了權(quán)重系數(shù)和相關(guān)系數(shù)均時(shí)變的時(shí)變混合copula函數(shù)[26]，綜合了時(shí)變copula 和混合copula 的優(yōu)點(diǎn)，這將是筆者下一步工作的研究重點(diǎn)。

4 結(jié)論

本研究將時(shí)變copula 函數(shù)引入到肌間耦合分析中，提出了一種時(shí)變copula 互信息估計(jì)方法，用于描述不同特征頻段上的肌間耦合關(guān)系。 實(shí)驗(yàn)結(jié)果表明，在腕屈、腕展運(yùn)動(dòng)過(guò)程中，相比靜態(tài)copula函數(shù)，時(shí)變copula 函數(shù)對(duì)肌間相依結(jié)構(gòu)的似然程度更高；由時(shí)變copula 互信息導(dǎo)出的肌間耦合強(qiáng)度存在顯著的頻段差異，具體表現(xiàn)為頻段越高肌間耦合強(qiáng)度越低，而靜態(tài)copula 互信息低估了肌間耦合強(qiáng)度；進(jìn)一步發(fā)現(xiàn)，肌間耦合強(qiáng)度與受試者的年齡、身高和體重?zé)o關(guān)。 時(shí)變copula 互信息是一種十分有效的肌間耦合分析手段，對(duì)挖掘潛在的中樞神經(jīng)系統(tǒng)運(yùn)動(dòng)控制機(jī)制具有良好的應(yīng)用價(jià)值。