為了更加準確有效地從復雜網絡的角度理解不同時空層次的肌間耦合情況，本文提出了一種新的多尺度肌間耦合網絡分析方法。將多元變分模態分解（MVMD）與 Copula 互信息（Copula MI）相結合，構建了基于 MVMD-Copula MI 的肌間耦合網絡模型，通過節點強度、聚類系數等網絡參數分析了健康人數據集伸手運動過程中上肢多塊肌肉在不同時頻尺度上的肌間耦合特性。實驗結果表明，在分解出的 6 個時頻尺度上，肌間耦合特性存在明顯區別。具體為：肱三頭肌與三角肌中束、三角肌后束的耦合強度相對較高，肌間功能聯系緊密；而肱二頭肌在該運動下獨立于其他肌肉。結果體現了肌間耦合網絡具有尺度差異性，而 MVMD-Copula MI 能夠定量刻畫多尺度肌間耦合強度關系，具有良好的應用前景。
引用本文： 吳亞婷, 佘青山, 高云園, 譚同才, 范影樂. 多尺度肌間耦合網絡分析. 生物醫學工程學雜志, 2021, 38(4): 742-752. doi: 10.7507/1001-5515.202009023 復制
引言
表面肌電信號（surface electromyogram，sEMG）是眾多肌纖維中運動單元動作電位（motor unit action potentials，MUAP）在時間和空間上的疊加[1]。肌間耦合的概念來源于皮層肌肉功能耦合研究，指的是運動過程中肌肉間的相互作用[2]，研究 sEMG 信號間的耦合信息不僅可以反映中樞神經系統（central nervous system，CNS）的不同運動控制策略和肌肉的運動功能狀態，還能夠實現肌電信號解碼，探索內在的運動功能控制機制。
近年來，有研究表明肌間耦合存在頻段顯著特征[3]，基于該特性有學者提出將小波分解[4]、經驗模態分解（empirical mode decomposition，EMD）[5]等時頻分解方法引入到肌間耦合分析之中，以偵測不同時頻尺度下的肌間耦合特性。然而，小波分解受其小波基的約束缺乏自適應能力，而 EMD 在處理信號時易受噪聲影響，并且分解出的模態之間混疊現象嚴重[6]。為了改進上述問題，Wu 等[7]提出了集合經驗模態分解算法，通過引入白噪聲輔助并集合平均的方式，有效解決了尺度混合等問題，但其運算時間過長，且缺乏數學理論支撐。為此，Dragomiretskiy 等[8]提出的變分模態分解（variational mode decomposition，VMD）可將隨機信號分解為多個窄帶分量，各分量包含原信號在不同時-頻尺度上的特征，同時 VMD 與 EMD 算法相比具有抗噪能力強、辨識精度高等特點。此后，謝平等[9]將 VMD 與傳遞熵相結合，揭示了皮層-肌肉功能耦合具有雙向性，體現了腦肌電信號分解在耦合分析中有良好的應用前景。然而，實際問題通常面對多元數據，為了同時實現多通道分解，Rehman 等[10]提出了多元經驗模態分解（multivariate empirical mode decomposition，MEMD）算法，但因 MEMD 魯棒性能差，仍然保留了 EMD 算法的缺點，他們[11]最近又提出了多元變分模態分解（multivariate variational mode decomposition，MVMD），可以有效地同時對多個信號進行頻率分解，因此 MVMD 可以為多元隨機變量尺度分解提供新的參考。
由于 sEMG 信號具有非線性、非平穩等特點，傳統的耦合分析方法如皮爾遜相關系數[12]、格蘭杰因果[13]、相干性[14]等更多的是在分析兩兩通道之間的線性關聯性，存在一定的局限性。互信息（mutual information，MI）是信息論里一種有用的信息度量工具，可以看成一個隨機變量中包含的關于另一個隨機變量的信息量，即兩個隨機變量之間的依賴程度。MI 廣泛應用于非線性系統，因此許多研究人員考慮采用此方法研究腦肌電信號的非線性耦合特征[15]。復雜網絡理論的發展為研究功能網絡中信息的處理和流動奠定了基礎，可以實現多通道肌電信號間的耦合分析，從整體角度分析肌間運動機制[16]。Boonstra 等[17]對十塊腿部肌肉的 sEMG 信號進行了連通性分析，以提取肌肉網絡，使用聚類系數、全局效應、中心性等網絡參數來評估肌肉協同作用。Kerkman 等[18]結合骨骼系統建立全身肌肉功能網絡，分析了不同站立姿勢動作下全身肌肉之間的耦合特性。陳玲玲等[19]將復雜網絡應用到外骨骼機器人研究中，通過 MI 度量肌間關聯性，構建及分析搬運過程中上肢肌間功能網絡。然而，傳統的 MI 估算方法需要聯合概率密度函數的精確表示，但現有的聯合概率密度函數估計方法存在估計精度低，依賴模型假設，以及要求變量個數少、樣本量充足等問題，所以 MI 的估計往往十分困難。Ma 等[20]根據 Sklar 定理推導出 MI 與 Copula 熵之間具有等價關系，這為估計 MI 提供了一條新的思路。在這之后，Ince 等[21]將 Copula 統計理論與高斯變量熵的封閉解相結合，提出了一種實用的高斯 Copula 互信息（Gaussian Copula mutual information，GCMI）估計方法，GCMI 為 MI 提供了一個計算效率高、穩健的下界估計，且無需對每個變量的邊際分布做出具體假設。
為了準確定量描述不同時頻尺度下肌間網絡的功能耦合情況，本文提出了將 MVMD 和 Copula 互信息（Copula mutual information，Copula MI）相結合的多尺度肌間耦合網絡分析方法。首先通過 MVMD 多尺度分解提取不同時頻尺度下 sEMG 信號的模態分量，然后利用 Copula MI 度量不同模態上肌間耦合強度，進一步基于復雜網絡理論搭建肌間耦合網絡，最后選取節點強度、聚類系數等網絡參數分析不同尺度下的肌間耦合特性，為解碼肌電信號以及探索潛在的 CNS 神經運動控制機制提供新的理論依據。
1 多尺度肌間耦合網絡的構建與分析
1.1 多尺度分解
MVMD 是一種自適應分解方法，利用迭代更新的方式求取變分模態的最優解，是 VMD 在多元空間中的一種擴展。MVMD 的主要思想是從包含 
個輸入通道數據 
中提取預定義的 
個尺度分量 
，滿足：① 提取分量的帶寬之和最小；② 各 IMF 分量之和等于輸入信號。
為了達到上述要求，首先利用希爾伯特變換得到調制振蕩向量 
的解析表示 
，然后通過單頻率分量 
移動所有尺度 
的單邊頻譜，使得每個分量的頻譜被調制到對應的中心頻率 
上，最后求偏導并取 Frobenius 2 范數的平方，由此產生 MVMD 的代價函數 
，受約束的變分模型表達式為：
 |
其中，
表示對應于通道 
的第 
個模式的解析表示。基于上述模型，對應的增廣拉格朗日函數表達式為：
 |
其中，
為懲罰因子，用于提高算法的抗噪性能和信號的重構精度。
為拉格朗日乘子，保證約束條件的嚴格性。為解決上述變分問題，可采用交替方向乘子法交替更新 
、
和 
：
 |
 |
 |
每個變分模態分量在求解過程中不斷迭代，迭代終止條件為：
 |
其中，誤差值 
。當滿足上述條件時，停止迭代，輸出模態 
和中心頻率 
[11]。
與 VMD 分解相同，MVMD 算法需要預先設定分解模態個數 
，本文將通過 MEMD 分解獲得的 IMF 模態數目作為參考[22]，結合中心頻率觀察法，確定 MVMD 的預設模態數目。
1.2 肌間耦合網絡構建
1.2.1 定義網絡節點、連邊
在肌間耦合網絡中，通常選取肌肉通道作為網絡節點[23]。對于網絡連邊的建立，本文選用 Copula MI[21]度量節點間關聯性，得到鄰接矩陣 
。
變量 
、
之間的 MI 定義如下：
 |
其中，
表示變量 
和 
的聯合概率密度函數，
和 
分別為變量 
和 
的邊際概率密度函數。根據文獻[24]，MI 與 Copula 熵存在以下關系：
 |
其中，
和 
分別為變量 
和 
的累積分布函數，
為變量 
和 
的 Copula 熵。Ince 等[21]將由式（8）得到的 MI 稱為 Copula MI，估計它的關鍵在于 Copula 函數的選擇。在理論上講，很難找到一個合適的 Copula 函數能足夠靈活地描述多變量的相依關系。為了更準確地描述隨機變量之間的關聯性，首先需要確定數據的邊際分布類型，然后選擇合適的 Copula 函數來連接邊際分布。現有的 Copula 函數主要有橢圓 Copula 和阿基米德 Copula，其中橢圓 Copula 包括 Gaussian Copula 和 t Copula，阿基米德 Copula 包括 Clayton Copula、Frank Copula 和 Gumbel Copula。不同的 Copula 函數形式不同，對 Copula MI 的估計準確性會有影響。因此，本文采用非參數核密度估計法估計各模態分量的邊際分布，根據最大似然函數原則從上述 5 種 Copula 中選擇最佳的 Copula 函數去估計兩兩通道之間的 Copula MI，以此作為肌間耦合網絡的連接權值。
1.2.2 確定網絡閾值
功能網絡通常包含虛假連接，這些連接可能會掩蓋重要的拓撲結構關系。為解決該問題，閾值化是一種很好的選擇。考慮到肌間耦合網絡的代數連通性，避免孤立節點的存在，本文引入圖論中拉普拉斯矩陣的第二小特征值 
準則[25]，將 
的 10%～90% 遞增值作為閾值，計算對應的拉普拉斯矩陣的第二小特征值 
，選取 
為零時的前一個閾值作為網絡的最佳閾值 
。如果 
，則 
；反之，
。
1.3 網絡參數分析
本文選用了復雜網絡中的節點強度、聚類系數和特征路徑長度三種網絡參數[26]，從網絡水平和節點水平分析肌間耦合特性。
1.3.1 節點強度
在加權網絡的情況下，節點強度 
等于連接到某節點的所有邊的權值的平均：
 |
其中，
是節點數，
是節點 
、
間的權值。節點強度既考慮了網絡的連通性，也考慮了邊的權重，并且取值在 
。節點強度越大，體現該節點構建和維系有效信息的作用越關鍵。
網絡強度 
定義為節點強度的平均值：
 |
網絡強度取值在 
范圍內，網絡強度越大，表明各通道之間關聯性越強；反之，網絡各通道間相互的關聯性較弱。
1.3.2 聚類系數
節點的聚類系數 
保留了所有連接邊的權重信息：
 |
它是用來描述一個網絡中的節點之間集結成團程度的系數。具體來說，節點聚類系數越大，與其他節點之間相互連接的程度越高。
網絡的聚類系數 
是所有節點的聚類系數的平均值：
 |
聚類系數越高，表明網絡的連通性越好[27]。
1.3.3 特征路徑長度
最短路徑長度 
描述了連接一對節點所需最少的“步驟”數量，節點之間通常可能有多個最短路徑長度，而特征路徑長度 
是指某節點所有的最短路徑長度的平均值：
 |
其中，
可根據 Dijikstra 算法[28]確定。特征路徑長度越小，表明節點的信息傳輸速度越快。
網絡的特征路徑長度 
定義為節點特征路徑長度的平均值：
 |
網絡的特征路徑長度越短，表明網絡的信息傳輸速度越快；反之，網絡傳輸信息的速度越慢。
2 實驗結果
2.1 數據介紹以及預處理
本文收集的實驗數據來源于 Israely 等[29]的工作，允許公開使用（https://data.mendeley.com/datasets/f4hh43nd78/2）。該數據包含 9 名健康受試者（S1-S9），年齡（68.6 ± 7.5）歲。受試者在一個簡易手感空間裝置下，根據 sEMG 軟件語音提示，使用他們的慣用手觸碰正前方目標，目標距桌面 35 cm，每個動作重復五次，為防止疲勞，每次動作后休息 10 s。運動過程中同步采集肩膀和手臂 8 塊肌肉上的 sEMG 信號：上斜方肌（upper trapezius，UT）、前三角肌（anterior deltoid，AD）、內側三角肌（medial deltoid，MD）、后三角肌（posterior deltoid，PD）、胸大肌（pectoralis major，PM）、岡下肌（infraspinatus，IN）、肱二頭肌（biceps brachii，BB）、肱三頭肌（triceps brachi，TB），采樣頻率為 2 000 Hz。
由于 sEMG 信號微弱，易受高頻噪聲影響，在進行尺度分解前需要對采集的 sEMG 信號進行預處理。首先通過重采樣的方式維持信號時長在 2.5 s，然后進行去均值、去除基線漂移，同時利用無限沖激響應陷波濾波器抑制 50 Hz 工頻干擾，最后采用 4 階巴特沃斯帶通濾波器進行 5～200 Hz 帶通濾波。預處理后受試者 S1 的 8 通道 sEMG 信號如圖 1 所示。
圖1
受試者 S1 預處理后的 8 通道 sEMG 信號
Figure1.
Subject S1’ pre-processed 8-channel sEMG signal
2.2 MVMD 分解結果
為體現 MVMD 的分解性能，分別對 9 名受試者預處理后的 sEMG 信號進行 MEMD、MVMD 多尺度分解及結果比較。受試者 S1 的 MEMD、MVMD 時頻分解結果如圖 2 所示。圖 2 左邊是各 IMF 分量的時域波形，右邊是對應 IMF 分量的頻譜，頻率分布從高到低依次排列。由圖 2 可見，MEMD 自適應地將 sEMG 信號分解成了 6 個 IMF 分量，MVMD 預設分解個數 
以其值作為參考。MEMD 分解的 IMF 分量帶寬較寬，頻帶間存在混疊，而 MVMD 的分解效果明顯更好，各模態間基本不混疊。進一步地，本文使用擴展的迪基-福勒（augmented Dickey-Fuller，ADF）單位根方法檢驗各模態的平穩性，圖 3 展示了 6 個模態分量的 ADF 單位根檢驗的 P 值，深藍色部分表示 P < 0.05。由圖 3 可見，2 種分解算法下各肌肉通道的前 4 個 IMF 分量都具有平穩性，但 MEMD 分解出的 IMF5、IMF6 分量的平穩情況（P > 0.05）不如 MVMD 分解效果。
圖2
MEMD、MVMD 多尺度時頻分解結果
Figure2.
Results of MEMD and MVMD
圖3
各模態分量平穩性檢驗的 P 值
Figure3.
The P-value of each mode stability test
2.3 Copula MI 構建肌間耦合網絡
針對 MVMD 算法進行多尺度分解后的 6 個 IMF 模態分量，本文采用基于 Copula MI 的方法度量各模態間的耦合強度關系。首先采用非參數核密度估計法，選用高斯核函數，根據經驗法則確定帶寬，估計各模態分量的邊際分布。圖 4a 展示了 TB 的各模態分量的直方圖概率密度分布，圖中曲線表示均值為 0、方差為 1 的高斯概率密度分布。進一步計算了各模態分量的峰度（kurtosis），發現其峰度值遠大于 3（高斯分布的峰度）。sEMG 數據的概率密度分布呈尖峰厚尾形狀，雖然與高斯分布一樣都是對稱分布，但是分布情況還是有明顯的差異。然后，計算了各模態分量的累計分布函數曲線如圖 4b 所示，子圖標注為 KS 檢驗統計量 H 值和 P 值。經 KS 檢驗后發現均符合[0,1]均勻分布（P < 0.05），滿足 Copula 模型建立的應用前提。最后，本文比較了 5 種常見類型的 Copula 函數：Gaussian Copula、t Copula、Clayton Copula、Frank Copula、Gumbel Copula，選擇具有最大似然函數值的 Copula 函數來估計 Copula MI。
圖4
TB 各模態分量的概率分布
a. 邊際概率密度函數曲線；b. 累計概率分布函數曲線
Figure4. Probability distribution of each mode component of TBa. marginal probability density function curve; b. cumulative probability distribution function curve
以受試者 S1 為例，圖 5 展示了各模態下計算兩兩通道之間 Copula MI 的過程，圖 5a 的矩陣表示似然函數值對比最大時所對應的 Copula 函數類型選擇結果，其中 G 表示 Gaussian Copula，t 表示 t Copula，C 表示 Clayton Copula，F 表示 Frank Copula，Gu 表示 Gumbel Copula。
圖5
肌間網絡連邊權值的計算過程
a. Copula 類型選擇；b. MI 鄰接矩陣；c. 閾值化后的鄰接矩陣
Figure5. The calculation process of the edge weight of the intermuscular networka. Copula type selection; b. MI adjacency matrix; c. adjacency matrix after thresholding
從圖 5 中可以看出，6 個模態下大部分肌肉通道之間都是選用的 t Copula 作為連接函數，這說明 sEMG 信號間的尾部相關變化對稱，當然也不排除尾部漸進獨立的情況，即對應 Frank Copula，這與圖 4 給出的數據分布特征相一致。圖 5b 是在圖 5a 的基礎上，通過對應的 Copula 函數估計出的 Copula MI 鄰接矩陣。圖 5c 則是按拉普拉斯矩陣特征值法閾值化后的鄰接矩陣。由圖 5 可知，AD 與 TB 在 IMF1、IMF3、IMF4 和 IMF6 模態分量上相比其他肌肉通道間具有更高的 Copula MI 值，而 MD 與 TB 則在 IMF2 和 IMF5 模態分量上相比其他肌肉通道間具有更高的 Copula MI 值，說明在伸手運動過程中主要是 TB 與三角肌的 AD 和 MD 之間產生較強的聯系。此外，PD 與 IN、PM、BB 和 TB 在 IMF5 模態分量上也存在較高的 Copula MI 值，這說明肌間耦合存在不同尺度的特異性。然而，BB 在大多數尺度上與其他肌肉通道間的 Copula MI 值都明顯較低，這說明 BB 在該運動下與其他肌肉的耦合程度較弱。
接下來，將 9 名受試者各個尺度分量下閾值化后的鄰接矩陣進行平均，抽象成圖 6 所示的網絡拓撲結構，節點的大小表示節點強度，邊的粗細表示肌間耦合強度的可視化。連線越粗，肌間耦合越強，反之，則越弱。從圖 6 可以看出，PD、MD 和 TB 在 6 個 IMF 模態分量上都存在較高的節點強度（> 0.035 Bit）和連接權值，肌間耦合網絡在前 3 個尺度分量（IMF1～IMF3）上具有一定的相似性，而在后 3 個尺度分量（IMF4～IMF6）上，肌間耦合關系對比差異明顯。相同的是，BB 在 6 個 IMF 尺度分量上節點強度都較低（< 0.015 Bit），并且與其他肌肉間的連線少、耦合作用弱。由圖 6 中連接邊可知，MD-AD 間的耦合作用在 IMF6 模態下很強，在 IMF5、IMF4 關聯性又較弱；TB-AD 在前三個模態下關聯性強，后三個模態下耦合關系不明顯，這也體現了肌間耦合在不同時頻尺度上具有不同的關聯程度特性。需要說明的是，在一定情況下圖 6 中的細實線可視為無連接，這里給出的是平均結果，因而個體間的差異被保留在肌間耦合網絡中。
圖6
各 IMF 尺度分量下的肌間耦合網絡
Figure6.
Intermuscular coupling networks at each IMF scale
2.4 網絡參數分析
為了更好地分析圖 6 所示的肌間耦合網絡特性，本文結合復雜網絡分析方法，計算了各尺度下肌肉節點的節點強度、聚類系數、特征路徑，以及尺度間的網絡強度、網絡聚類系數、網絡特征路徑。圖 7 為各 IMF 模態分量下肌間耦合網絡中 8 個肌肉節點的節點強度、聚類系數和路徑長度的箱型圖，通過對三個參數的綜合分析，發現肌間耦合網絡下重要性較大的節點為該網絡下的關鍵節點。從箱型圖的中位數可以看出，在 6 個 IMF 模態分量下，三角肌組 AD、MD 和 PD 以及 TB 均具有較高的節點強度和聚類系數，其中 TB 的值明顯高于其他肌肉。而 UT 和 BB 的節點強度明顯低于其他肌肉，特別是 BB 的聚類系數幾乎為零，這說明 BB 可能獨立于其他肌肉，與其他肌肉間的連通程度低。從節點的最短路徑來看，雖然 IMF3～IMF5 模態分量上存在大量異常值，但總體上 AD、MD、PD 和 TB 的最短路徑要低于其他肌肉節點，而 BB 則恰恰相反。
圖7
各 IMF 尺度分量下的網絡節點參數
Figure7.
Network node parameters at each IMF scale
圖 8 為各 IMF 尺度分量下的網絡強度、網絡聚類系數和網絡特征路徑的箱型圖，描述不同中心頻率尺度下的網絡屬性。從箱型圖的中位數可以看出，IMF5 尺度分量的網絡強度和聚類系數相對較高，這說明在 IMF5 尺度肌間耦合更加緊密，連通性也較好。總體上，從網絡特征路徑箱型圖可看出，特征路徑長度在 6 個尺度下是相近的，說明各尺度下肌間網絡的信息流通速度相差不大。
圖8
各 IMF 模態分量下的肌間耦合網絡特性
Figure8.
The characteristics of the intermuscular coupling network under each IMF mode
3 討論
本文基于 MVMD-Copula MI 方法，結合節點強度、聚類系數、特征路徑網絡參數分析健康人做伸手及物運動時上肢 8 塊肌肉的肌間耦合強度。在多尺度分解算法的選擇上，MEMD 方法是基于極值點的包絡求取方式，多次遞歸分解后包絡估計誤差被放大，不能將相近的頻率分量正確分離[30]，而 MVMD 獲得的模式顯示出較少的瞬時頻率波動，在噪聲魯棒性以及信號重建方面，表現出更加優異的效果。在計算 MI 時，使用 Copula 來估計 MI，很好地避免了計算聯合概率分布，同時分析了 sEMG 數據呈現尖峰厚尾形狀，采用最大似然函數法進行選擇 Copula 函數類型時，發現 t Copula 比其他 Copula 系數結構更加靈活，能夠更好地描述 sEMG 信號之間的耦合關系。
在復雜網絡分析中，信息感知的關鍵節點不僅與其所處中心位置有關，還與其他節點信息交流的時間有關[19]，在各個模態下 TB 是關鍵節點，體現了在該伸手及物運動下 TB 起了關鍵的功能性作用，且在各個尺度下 TB 與 MD、PD 間的耦合強度都相對較大，這與 Israely 等[29]利用非負矩陣分解算法提取出的肌肉協同單元具有一定的相似性。BB 是非關鍵節點，說明在該伸手及物運動下 BB 功能性作用較弱，這可能是受限于 BB 的生理解剖結構和運動學背景。MD 與 AD 間的耦合作用在 IMF6 模態下很強，在 IMF5、IMF4 關聯性又較弱；TB 與 AD 在前三個模態下關聯性強，后三個模態下耦合關系不顯著。這與 Boonstra 等[17]得出的肌肉間是分層協調運動的結論基本相似，體現了肌間耦合關系在多尺度下是不同步的。通過對多個尺度下肌間網絡的網絡平均節點強度、網絡聚類系數以及網絡特征路徑長度的分析，發現 IMF5 尺度分量各肌肉間聯系的關聯性最強，網絡連通性好，在該運動下對維系肌間有效信息起關鍵作用，可能與肌間耦合在低頻段的顯著特性相關[31]。同時不難發現，IMF5 尺度分量主要分布在 beta（15～30 Hz）和 gamma（30～60 Hz）頻段，而這兩個頻段被認為代表了靜態力調整至動態力輸出時，大腦皮層對肌肉系統在神經振蕩方面的轉移[32]。以上說明，運動過程中不同尺度分量上的神經通路差異，或許是運動控制及響應方式的神經生理學基礎，表現了支配脊髓運動神經元參與運動協調的一般原則。
4 結論
本文將 MVMD 方法引入到肌間耦合分析中，并與 Copula MI 方法相結合，聯系復雜網絡分析思想，從整體角度建立多尺度肌間耦合網絡分析模型，應用于健康人伸手及物運動下 sEMG 信號耦合特性分析。結果表明不同時頻尺度下，起關鍵作用的肌肉具有一致性，肌間網絡存在一定的差異性。在各個模態下 TB 是關鍵節點，起了關鍵的功能性作用，TB 與 MD、PD 間的耦合強度都較大，而 BB 則功能性作用較弱，是非關鍵節點。肌間在不同尺度下展現了不同的耦合強度，這體現了肌間耦合關系具有分層分布特性。由此，本文提出的多尺度肌間耦合分析方法可以有效刻畫 sEMG 信號在不同時頻尺度下的耦合特性，可為肌電解碼以及理解運動功能機制提供新的途徑。
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
引言
表面肌電信號（surface electromyogram，sEMG）是眾多肌纖維中運動單元動作電位（motor unit action potentials，MUAP）在時間和空間上的疊加[1]。肌間耦合的概念來源于皮層肌肉功能耦合研究，指的是運動過程中肌肉間的相互作用[2]，研究 sEMG 信號間的耦合信息不僅可以反映中樞神經系統（central nervous system，CNS）的不同運動控制策略和肌肉的運動功能狀態，還能夠實現肌電信號解碼，探索內在的運動功能控制機制。
近年來，有研究表明肌間耦合存在頻段顯著特征[3]，基于該特性有學者提出將小波分解[4]、經驗模態分解（empirical mode decomposition，EMD）[5]等時頻分解方法引入到肌間耦合分析之中，以偵測不同時頻尺度下的肌間耦合特性。然而，小波分解受其小波基的約束缺乏自適應能力，而 EMD 在處理信號時易受噪聲影響，并且分解出的模態之間混疊現象嚴重[6]。為了改進上述問題，Wu 等[7]提出了集合經驗模態分解算法，通過引入白噪聲輔助并集合平均的方式，有效解決了尺度混合等問題，但其運算時間過長，且缺乏數學理論支撐。為此，Dragomiretskiy 等[8]提出的變分模態分解（variational mode decomposition，VMD）可將隨機信號分解為多個窄帶分量，各分量包含原信號在不同時-頻尺度上的特征，同時 VMD 與 EMD 算法相比具有抗噪能力強、辨識精度高等特點。此后，謝平等[9]將 VMD 與傳遞熵相結合，揭示了皮層-肌肉功能耦合具有雙向性，體現了腦肌電信號分解在耦合分析中有良好的應用前景。然而，實際問題通常面對多元數據，為了同時實現多通道分解，Rehman 等[10]提出了多元經驗模態分解（multivariate empirical mode decomposition，MEMD）算法，但因 MEMD 魯棒性能差，仍然保留了 EMD 算法的缺點，他們[11]最近又提出了多元變分模態分解（multivariate variational mode decomposition，MVMD），可以有效地同時對多個信號進行頻率分解，因此 MVMD 可以為多元隨機變量尺度分解提供新的參考。
由于 sEMG 信號具有非線性、非平穩等特點，傳統的耦合分析方法如皮爾遜相關系數[12]、格蘭杰因果[13]、相干性[14]等更多的是在分析兩兩通道之間的線性關聯性，存在一定的局限性。互信息（mutual information，MI）是信息論里一種有用的信息度量工具，可以看成一個隨機變量中包含的關于另一個隨機變量的信息量，即兩個隨機變量之間的依賴程度。MI 廣泛應用于非線性系統，因此許多研究人員考慮采用此方法研究腦肌電信號的非線性耦合特征[15]。復雜網絡理論的發展為研究功能網絡中信息的處理和流動奠定了基礎，可以實現多通道肌電信號間的耦合分析，從整體角度分析肌間運動機制[16]。Boonstra 等[17]對十塊腿部肌肉的 sEMG 信號進行了連通性分析，以提取肌肉網絡，使用聚類系數、全局效應、中心性等網絡參數來評估肌肉協同作用。Kerkman 等[18]結合骨骼系統建立全身肌肉功能網絡，分析了不同站立姿勢動作下全身肌肉之間的耦合特性。陳玲玲等[19]將復雜網絡應用到外骨骼機器人研究中，通過 MI 度量肌間關聯性，構建及分析搬運過程中上肢肌間功能網絡。然而，傳統的 MI 估算方法需要聯合概率密度函數的精確表示，但現有的聯合概率密度函數估計方法存在估計精度低，依賴模型假設，以及要求變量個數少、樣本量充足等問題，所以 MI 的估計往往十分困難。Ma 等[20]根據 Sklar 定理推導出 MI 與 Copula 熵之間具有等價關系，這為估計 MI 提供了一條新的思路。在這之后，Ince 等[21]將 Copula 統計理論與高斯變量熵的封閉解相結合，提出了一種實用的高斯 Copula 互信息（Gaussian Copula mutual information，GCMI）估計方法，GCMI 為 MI 提供了一個計算效率高、穩健的下界估計，且無需對每個變量的邊際分布做出具體假設。
為了準確定量描述不同時頻尺度下肌間網絡的功能耦合情況，本文提出了將 MVMD 和 Copula 互信息（Copula mutual information，Copula MI）相結合的多尺度肌間耦合網絡分析方法。首先通過 MVMD 多尺度分解提取不同時頻尺度下 sEMG 信號的模態分量，然后利用 Copula MI 度量不同模態上肌間耦合強度，進一步基于復雜網絡理論搭建肌間耦合網絡，最后選取節點強度、聚類系數等網絡參數分析不同尺度下的肌間耦合特性，為解碼肌電信號以及探索潛在的 CNS 神經運動控制機制提供新的理論依據。
1 多尺度肌間耦合網絡的構建與分析
1.1 多尺度分解
MVMD 是一種自適應分解方法，利用迭代更新的方式求取變分模態的最優解，是 VMD 在多元空間中的一種擴展。MVMD 的主要思想是從包含 個輸入通道數據 中提取預定義的 個尺度分量 ，滿足：① 提取分量的帶寬之和最小；② 各 IMF 分量之和等于輸入信號。
為了達到上述要求，首先利用希爾伯特變換得到調制振蕩向量 的解析表示 ，然后通過單頻率分量 移動所有尺度 的單邊頻譜，使得每個分量的頻譜被調制到對應的中心頻率 上，最后求偏導并取 Frobenius 2 范數的平方，由此產生 MVMD 的代價函數 ，受約束的變分模型表達式為：
|
其中， 表示對應于通道 的第 個模式的解析表示。基于上述模型，對應的增廣拉格朗日函數表達式為：
|
其中， 為懲罰因子，用于提高算法的抗噪性能和信號的重構精度。 為拉格朗日乘子，保證約束條件的嚴格性。為解決上述變分問題，可采用交替方向乘子法交替更新 、 和 ：
|
|
|
每個變分模態分量在求解過程中不斷迭代，迭代終止條件為：
|
其中，誤差值 。當滿足上述條件時，停止迭代，輸出模態 和中心頻率 [11]。
與 VMD 分解相同，MVMD 算法需要預先設定分解模態個數 ，本文將通過 MEMD 分解獲得的 IMF 模態數目作為參考[22]，結合中心頻率觀察法，確定 MVMD 的預設模態數目。
1.2 肌間耦合網絡構建
1.2.1 定義網絡節點、連邊
在肌間耦合網絡中，通常選取肌肉通道作為網絡節點[23]。對于網絡連邊的建立，本文選用 Copula MI[21]度量節點間關聯性，得到鄰接矩陣 。
變量 、 之間的 MI 定義如下：
|
其中， 表示變量 和 的聯合概率密度函數， 和 分別為變量 和 的邊際概率密度函數。根據文獻[24]，MI 與 Copula 熵存在以下關系：
|
其中， 和 分別為變量 和 的累積分布函數， 為變量 和 的 Copula 熵。Ince 等[21]將由式（8）得到的 MI 稱為 Copula MI，估計它的關鍵在于 Copula 函數的選擇。在理論上講，很難找到一個合適的 Copula 函數能足夠靈活地描述多變量的相依關系。為了更準確地描述隨機變量之間的關聯性，首先需要確定數據的邊際分布類型，然后選擇合適的 Copula 函數來連接邊際分布。現有的 Copula 函數主要有橢圓 Copula 和阿基米德 Copula，其中橢圓 Copula 包括 Gaussian Copula 和 t Copula，阿基米德 Copula 包括 Clayton Copula、Frank Copula 和 Gumbel Copula。不同的 Copula 函數形式不同，對 Copula MI 的估計準確性會有影響。因此，本文采用非參數核密度估計法估計各模態分量的邊際分布，根據最大似然函數原則從上述 5 種 Copula 中選擇最佳的 Copula 函數去估計兩兩通道之間的 Copula MI，以此作為肌間耦合網絡的連接權值。
1.2.2 確定網絡閾值
功能網絡通常包含虛假連接，這些連接可能會掩蓋重要的拓撲結構關系。為解決該問題，閾值化是一種很好的選擇。考慮到肌間耦合網絡的代數連通性，避免孤立節點的存在，本文引入圖論中拉普拉斯矩陣的第二小特征值 準則[25]，將 的 10%～90% 遞增值作為閾值，計算對應的拉普拉斯矩陣的第二小特征值 ，選取 為零時的前一個閾值作為網絡的最佳閾值 。如果 ，則 ；反之，。
1.3 網絡參數分析
本文選用了復雜網絡中的節點強度、聚類系數和特征路徑長度三種網絡參數[26]，從網絡水平和節點水平分析肌間耦合特性。
1.3.1 節點強度
在加權網絡的情況下，節點強度 等于連接到某節點的所有邊的權值的平均：
|
其中， 是節點數， 是節點 、 間的權值。節點強度既考慮了網絡的連通性，也考慮了邊的權重，并且取值在 。節點強度越大，體現該節點構建和維系有效信息的作用越關鍵。
網絡強度 定義為節點強度的平均值：
|
網絡強度取值在 范圍內，網絡強度越大，表明各通道之間關聯性越強；反之，網絡各通道間相互的關聯性較弱。
1.3.2 聚類系數
節點的聚類系數 保留了所有連接邊的權重信息：
|
它是用來描述一個網絡中的節點之間集結成團程度的系數。具體來說，節點聚類系數越大，與其他節點之間相互連接的程度越高。
網絡的聚類系數 是所有節點的聚類系數的平均值：
|
聚類系數越高，表明網絡的連通性越好[27]。
1.3.3 特征路徑長度
最短路徑長度 描述了連接一對節點所需最少的“步驟”數量，節點之間通常可能有多個最短路徑長度，而特征路徑長度 是指某節點所有的最短路徑長度的平均值：
|
其中， 可根據 Dijikstra 算法[28]確定。特征路徑長度越小，表明節點的信息傳輸速度越快。
網絡的特征路徑長度 定義為節點特征路徑長度的平均值：
|
網絡的特征路徑長度越短，表明網絡的信息傳輸速度越快；反之，網絡傳輸信息的速度越慢。
2 實驗結果
2.1 數據介紹以及預處理
本文收集的實驗數據來源于 Israely 等[29]的工作，允許公開使用（https://data.mendeley.com/datasets/f4hh43nd78/2）。該數據包含 9 名健康受試者（S1-S9），年齡（68.6 ± 7.5）歲。受試者在一個簡易手感空間裝置下，根據 sEMG 軟件語音提示，使用他們的慣用手觸碰正前方目標，目標距桌面 35 cm，每個動作重復五次，為防止疲勞，每次動作后休息 10 s。運動過程中同步采集肩膀和手臂 8 塊肌肉上的 sEMG 信號：上斜方肌（upper trapezius，UT）、前三角肌（anterior deltoid，AD）、內側三角肌（medial deltoid，MD）、后三角肌（posterior deltoid，PD）、胸大肌（pectoralis major，PM）、岡下肌（infraspinatus，IN）、肱二頭肌（biceps brachii，BB）、肱三頭肌（triceps brachi，TB），采樣頻率為 2 000 Hz。
由于 sEMG 信號微弱，易受高頻噪聲影響，在進行尺度分解前需要對采集的 sEMG 信號進行預處理。首先通過重采樣的方式維持信號時長在 2.5 s，然后進行去均值、去除基線漂移，同時利用無限沖激響應陷波濾波器抑制 50 Hz 工頻干擾，最后采用 4 階巴特沃斯帶通濾波器進行 5～200 Hz 帶通濾波。預處理后受試者 S1 的 8 通道 sEMG 信號如圖 1 所示。
圖1
受試者 S1 預處理后的 8 通道 sEMG 信號
Figure1.
Subject S1’ pre-processed 8-channel sEMG signal
2.2 MVMD 分解結果
為體現 MVMD 的分解性能，分別對 9 名受試者預處理后的 sEMG 信號進行 MEMD、MVMD 多尺度分解及結果比較。受試者 S1 的 MEMD、MVMD 時頻分解結果如圖 2 所示。圖 2 左邊是各 IMF 分量的時域波形，右邊是對應 IMF 分量的頻譜，頻率分布從高到低依次排列。由圖 2 可見，MEMD 自適應地將 sEMG 信號分解成了 6 個 IMF 分量，MVMD 預設分解個數 以其值作為參考。MEMD 分解的 IMF 分量帶寬較寬，頻帶間存在混疊，而 MVMD 的分解效果明顯更好，各模態間基本不混疊。進一步地，本文使用擴展的迪基-福勒（augmented Dickey-Fuller，ADF）單位根方法檢驗各模態的平穩性，圖 3 展示了 6 個模態分量的 ADF 單位根檢驗的 P 值，深藍色部分表示 P < 0.05。由圖 3 可見，2 種分解算法下各肌肉通道的前 4 個 IMF 分量都具有平穩性，但 MEMD 分解出的 IMF5、IMF6 分量的平穩情況（P > 0.05）不如 MVMD 分解效果。
圖2
MEMD、MVMD 多尺度時頻分解結果
Figure2.
Results of MEMD and MVMD
圖3
各模態分量平穩性檢驗的 P 值
Figure3.
The P-value of each mode stability test
2.3 Copula MI 構建肌間耦合網絡
針對 MVMD 算法進行多尺度分解后的 6 個 IMF 模態分量，本文采用基于 Copula MI 的方法度量各模態間的耦合強度關系。首先采用非參數核密度估計法，選用高斯核函數，根據經驗法則確定帶寬，估計各模態分量的邊際分布。圖 4a 展示了 TB 的各模態分量的直方圖概率密度分布，圖中曲線表示均值為 0、方差為 1 的高斯概率密度分布。進一步計算了各模態分量的峰度（kurtosis），發現其峰度值遠大于 3（高斯分布的峰度）。sEMG 數據的概率密度分布呈尖峰厚尾形狀，雖然與高斯分布一樣都是對稱分布，但是分布情況還是有明顯的差異。然后，計算了各模態分量的累計分布函數曲線如圖 4b 所示，子圖標注為 KS 檢驗統計量 H 值和 P 值。經 KS 檢驗后發現均符合[0,1]均勻分布（P < 0.05），滿足 Copula 模型建立的應用前提。最后，本文比較了 5 種常見類型的 Copula 函數：Gaussian Copula、t Copula、Clayton Copula、Frank Copula、Gumbel Copula，選擇具有最大似然函數值的 Copula 函數來估計 Copula MI。
圖4
TB 各模態分量的概率分布
a. 邊際概率密度函數曲線；b. 累計概率分布函數曲線
Figure4. Probability distribution of each mode component of TBa. marginal probability density function curve; b. cumulative probability distribution function curve
以受試者 S1 為例，圖 5 展示了各模態下計算兩兩通道之間 Copula MI 的過程，圖 5a 的矩陣表示似然函數值對比最大時所對應的 Copula 函數類型選擇結果，其中 G 表示 Gaussian Copula，t 表示 t Copula，C 表示 Clayton Copula，F 表示 Frank Copula，Gu 表示 Gumbel Copula。
圖5
肌間網絡連邊權值的計算過程
a. Copula 類型選擇；b. MI 鄰接矩陣；c. 閾值化后的鄰接矩陣
Figure5. The calculation process of the edge weight of the intermuscular networka. Copula type selection; b. MI adjacency matrix; c. adjacency matrix after thresholding
從圖 5 中可以看出，6 個模態下大部分肌肉通道之間都是選用的 t Copula 作為連接函數，這說明 sEMG 信號間的尾部相關變化對稱，當然也不排除尾部漸進獨立的情況，即對應 Frank Copula，這與圖 4 給出的數據分布特征相一致。圖 5b 是在圖 5a 的基礎上，通過對應的 Copula 函數估計出的 Copula MI 鄰接矩陣。圖 5c 則是按拉普拉斯矩陣特征值法閾值化后的鄰接矩陣。由圖 5 可知，AD 與 TB 在 IMF1、IMF3、IMF4 和 IMF6 模態分量上相比其他肌肉通道間具有更高的 Copula MI 值，而 MD 與 TB 則在 IMF2 和 IMF5 模態分量上相比其他肌肉通道間具有更高的 Copula MI 值，說明在伸手運動過程中主要是 TB 與三角肌的 AD 和 MD 之間產生較強的聯系。此外，PD 與 IN、PM、BB 和 TB 在 IMF5 模態分量上也存在較高的 Copula MI 值，這說明肌間耦合存在不同尺度的特異性。然而，BB 在大多數尺度上與其他肌肉通道間的 Copula MI 值都明顯較低，這說明 BB 在該運動下與其他肌肉的耦合程度較弱。
接下來，將 9 名受試者各個尺度分量下閾值化后的鄰接矩陣進行平均，抽象成圖 6 所示的網絡拓撲結構，節點的大小表示節點強度，邊的粗細表示肌間耦合強度的可視化。連線越粗，肌間耦合越強，反之，則越弱。從圖 6 可以看出，PD、MD 和 TB 在 6 個 IMF 模態分量上都存在較高的節點強度（> 0.035 Bit）和連接權值，肌間耦合網絡在前 3 個尺度分量（IMF1～IMF3）上具有一定的相似性，而在后 3 個尺度分量（IMF4～IMF6）上，肌間耦合關系對比差異明顯。相同的是，BB 在 6 個 IMF 尺度分量上節點強度都較低（< 0.015 Bit），并且與其他肌肉間的連線少、耦合作用弱。由圖 6 中連接邊可知，MD-AD 間的耦合作用在 IMF6 模態下很強，在 IMF5、IMF4 關聯性又較弱；TB-AD 在前三個模態下關聯性強，后三個模態下耦合關系不明顯，這也體現了肌間耦合在不同時頻尺度上具有不同的關聯程度特性。需要說明的是，在一定情況下圖 6 中的細實線可視為無連接，這里給出的是平均結果，因而個體間的差異被保留在肌間耦合網絡中。
圖6
各 IMF 尺度分量下的肌間耦合網絡
Figure6.
Intermuscular coupling networks at each IMF scale
2.4 網絡參數分析
為了更好地分析圖 6 所示的肌間耦合網絡特性，本文結合復雜網絡分析方法，計算了各尺度下肌肉節點的節點強度、聚類系數、特征路徑，以及尺度間的網絡強度、網絡聚類系數、網絡特征路徑。圖 7 為各 IMF 模態分量下肌間耦合網絡中 8 個肌肉節點的節點強度、聚類系數和路徑長度的箱型圖，通過對三個參數的綜合分析，發現肌間耦合網絡下重要性較大的節點為該網絡下的關鍵節點。從箱型圖的中位數可以看出，在 6 個 IMF 模態分量下，三角肌組 AD、MD 和 PD 以及 TB 均具有較高的節點強度和聚類系數，其中 TB 的值明顯高于其他肌肉。而 UT 和 BB 的節點強度明顯低于其他肌肉，特別是 BB 的聚類系數幾乎為零，這說明 BB 可能獨立于其他肌肉，與其他肌肉間的連通程度低。從節點的最短路徑來看，雖然 IMF3～IMF5 模態分量上存在大量異常值，但總體上 AD、MD、PD 和 TB 的最短路徑要低于其他肌肉節點，而 BB 則恰恰相反。
圖7
各 IMF 尺度分量下的網絡節點參數
Figure7.
Network node parameters at each IMF scale
圖 8 為各 IMF 尺度分量下的網絡強度、網絡聚類系數和網絡特征路徑的箱型圖，描述不同中心頻率尺度下的網絡屬性。從箱型圖的中位數可以看出，IMF5 尺度分量的網絡強度和聚類系數相對較高，這說明在 IMF5 尺度肌間耦合更加緊密，連通性也較好。總體上，從網絡特征路徑箱型圖可看出，特征路徑長度在 6 個尺度下是相近的，說明各尺度下肌間網絡的信息流通速度相差不大。
圖8
各 IMF 模態分量下的肌間耦合網絡特性
Figure8.
The characteristics of the intermuscular coupling network under each IMF mode
3 討論
本文基于 MVMD-Copula MI 方法，結合節點強度、聚類系數、特征路徑網絡參數分析健康人做伸手及物運動時上肢 8 塊肌肉的肌間耦合強度。在多尺度分解算法的選擇上，MEMD 方法是基于極值點的包絡求取方式，多次遞歸分解后包絡估計誤差被放大，不能將相近的頻率分量正確分離[30]，而 MVMD 獲得的模式顯示出較少的瞬時頻率波動，在噪聲魯棒性以及信號重建方面，表現出更加優異的效果。在計算 MI 時，使用 Copula 來估計 MI，很好地避免了計算聯合概率分布，同時分析了 sEMG 數據呈現尖峰厚尾形狀，采用最大似然函數法進行選擇 Copula 函數類型時，發現 t Copula 比其他 Copula 系數結構更加靈活，能夠更好地描述 sEMG 信號之間的耦合關系。
在復雜網絡分析中，信息感知的關鍵節點不僅與其所處中心位置有關，還與其他節點信息交流的時間有關[19]，在各個模態下 TB 是關鍵節點，體現了在該伸手及物運動下 TB 起了關鍵的功能性作用，且在各個尺度下 TB 與 MD、PD 間的耦合強度都相對較大，這與 Israely 等[29]利用非負矩陣分解算法提取出的肌肉協同單元具有一定的相似性。BB 是非關鍵節點，說明在該伸手及物運動下 BB 功能性作用較弱，這可能是受限于 BB 的生理解剖結構和運動學背景。MD 與 AD 間的耦合作用在 IMF6 模態下很強，在 IMF5、IMF4 關聯性又較弱；TB 與 AD 在前三個模態下關聯性強，后三個模態下耦合關系不顯著。這與 Boonstra 等[17]得出的肌肉間是分層協調運動的結論基本相似，體現了肌間耦合關系在多尺度下是不同步的。通過對多個尺度下肌間網絡的網絡平均節點強度、網絡聚類系數以及網絡特征路徑長度的分析，發現 IMF5 尺度分量各肌肉間聯系的關聯性最強，網絡連通性好，在該運動下對維系肌間有效信息起關鍵作用，可能與肌間耦合在低頻段的顯著特性相關[31]。同時不難發現，IMF5 尺度分量主要分布在 beta（15～30 Hz）和 gamma（30～60 Hz）頻段，而這兩個頻段被認為代表了靜態力調整至動態力輸出時，大腦皮層對肌肉系統在神經振蕩方面的轉移[32]。以上說明，運動過程中不同尺度分量上的神經通路差異，或許是運動控制及響應方式的神經生理學基礎，表現了支配脊髓運動神經元參與運動協調的一般原則。
4 結論
本文將 MVMD 方法引入到肌間耦合分析中，并與 Copula MI 方法相結合，聯系復雜網絡分析思想，從整體角度建立多尺度肌間耦合網絡分析模型，應用于健康人伸手及物運動下 sEMG 信號耦合特性分析。結果表明不同時頻尺度下，起關鍵作用的肌肉具有一致性，肌間網絡存在一定的差異性。在各個模態下 TB 是關鍵節點，起了關鍵的功能性作用，TB 與 MD、PD 間的耦合強度都較大，而 BB 則功能性作用較弱，是非關鍵節點。肌間在不同尺度下展現了不同的耦合強度，這體現了肌間耦合關系具有分層分布特性。由此，本文提出的多尺度肌間耦合分析方法可以有效刻畫 sEMG 信號在不同時頻尺度下的耦合特性，可為肌電解碼以及理解運動功能機制提供新的途徑。
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
