本文針對中風康復患者運動功能障礙相關的神經肌肉異常耦合和功能評價問題, 研究有效的肌間一致性分析方法及指標, 探索神經肌肉系統的振蕩聯系和中風運動功能障礙產生的病理機制, 建立基于肌電信號的肌肉功能狀態評價標準。首先基于一致性分析方法, 對比分析中風患者運動過程中健側、患側上肢拮抗肌肌電間的一致性特征; 在此基礎上利用Z變換顯著一致性面積指標量化描述健側、患側肌電信號在各功能頻帶內的耦合性差異; 進一步利用動態肌間一致性分析探究肌間一致性與運動任務的關系。通過對中風患者健側、患側上肢肘關節屈伸運動中肌間一致性的分析, 得到以下結論:中風運動功能障礙患者健側、患側肌間一致性均與運動任務相關, 而在患側運動過程中beta頻段肌間一致性相對于健側存在明顯缺失, 此外beta段肌間一致性與Fugle-Meyer量表具有相關性, 表明beta段肌間一致性可以作為輔助評價患者運動功能狀態的指標。
引用本文： 謝平, 宋妍, 郭子暉, 陳曉玲, 吳曉光, 蘇玉萍, 杜義浩. 中風康復運動中肌肉異常耦合分析. 生物醫學工程學雜志, 2016, 33(2): 244-254. doi: 10.7507/1001-5515.20160043 復制
引言
在人體運動的神經控制過程中，神經元的振蕩耦合和同步放電發揮著至關重要的作用。運動神經系統通過神經振蕩來傳遞運動控制信息[1]，神經中樞通過與各區域間的相互作用實現功能調節和整合，主要體現在電生理信號間不同層次的耦合(同步)現象，如腦電(electroencephalogram, EEG)-肌電(electromyogram, EMG)耦合、EMG-EMG耦合等。中樞神經損傷患者由于神經損傷阻礙了神經振蕩的傳導從而造成運動功能障礙。近年來，運用電生理信號耦合分析考察中樞神經損傷患者的異常神經振蕩和神經肌肉狀態已成為運動醫學、康復工程等領域的研究熱點。通過耦合特征分析可以反映運動控制過程中神經肌肉間的功能聯系，為理解運動控制過程及運動障礙的病理機制提供了理論基礎，也為神經康復運動的功能狀態評價提供了新方法。
EEG信號和表面肌電(surface electromyogram, sEMG)信號，二者分別作為人體指令中心和執行末端的代表性信號，在機體運動功能領域有著非常重要的研究意義。自從Conway等[2]發現運動過程中相關肌肉肌電和皮層腦電之間具有頻率成分的耦合性，基于腦肌電一致性(EEG-EMG/corticomuscular coherence, CMC)對大腦運動皮層與肌肉功能聯系的研究相繼展開，并應用于探索中風、帕金森癥等中樞神經系統相關疾病病理機制的研究中。研究表明，運動過程中CMC主要表現在beta頻段和gamma頻段，beta頻段的腦肌電振蕩代表了從初級運動皮層到脊髓運動神經元的傳遞過程[3]，而gamma頻段振蕩體現與認知功能相關的腦皮層信息整合過程[4]。因此，通過CMC分析能夠體現中樞神經與肌肉之間的振蕩聯系，進而評價運動神經系統功能，甚至指導康復過程。然而CMC分析在臨床應用中仍具有困難，原因主要有三個：①臨床中難以獲得高質量的EEG信號；②腦肌電同步采集一直是難以解決的問題；③CMC特征尤其是患者的一致性特征并不穩定。值得提出的是，Brown等研究者發現肌間一致性(EMG-EMG/intermuscular coherence, IMC)與CMC具有相似的模式，說明同一個受試者同一個任務下的CMC和IMC均起源于皮質脊髓通路[5-7]，因此可以用IMC來反映CMC的變化趨勢[8]。同時，他們發現IMC還可以反映脊髓運動神經元共同的突觸前驅動，使運動單位同步發放從而激發肌肉活性[3]。另一方面，在臨床環境下sEMG的采集更為安全便捷，信號可靠性也更強。綜上所述，利用IMC分析神經肌肉振蕩聯系、評價運動神經系統功能更具有臨床實用價值。
近年來，基于IMC分析運動過程中肌肉間的耦合模式、挖掘潛在的神經運動控制機制的研究相繼開展，如Baker等[9]、Nroton等[10]分別發現當肌肉持續收縮時的IMC最大；Nashner[11]發現拮抗肌最大自主收縮時的IMC在10 Hz左右最明顯；Farmer等[7]發現等長收縮運動下運動單位同步主要發生在1~12 Hz和16~32 Hz。此外，還有部分研究者將IMC分析引入運動障礙病理機制及運動功能評價的研究中，如Grosse等[12]發現肌張力障礙患者與正常受試者的IMC在theta(4~7 Hz)頻段具有顯著差異；Fisher等[13]發現beta(15~30 Hz)頻段IMC可以作為評價運動神經損傷患者上肢運動功能紊亂的指標；Kisiel-Sajewicz等[14]發現中風患者進行抓取運動時肌肉功能性耦合在0~3.9 Hz和4~7.9 Hz有缺失。然而關于中風康復運動中患側、健側肌肉耦合差異方面的深入理論分析還不多見，是一個值得進一步研究的問題。
為此，本文首先分別進行中風運動功能障礙患者肘部屈伸動作任務下健側與患側的上肢sEMG信號同步采集實驗。在此基礎上基于一致性分析方法，對整流后的肱二頭肌(biceps brachii, BB)和肱三頭肌(triceps brachii, TB) sEMG信號進行頻域一致性分析和時域累積密度估計分析，針對各特定頻段研究肌肉間的耦合模式和強度，以此推測神經肌肉系統的振蕩聯系，并利用Z變換顯著一致性面積指標量化和比較兩通道sEMG在各頻帶內的耦合性。進一步利用具有時間分辨率的動態一致性分析，探究IMC如何隨著實驗任務變化而變化，最后利用統計方法對比分析患者患側與健側IMC的差異。通過對中風患者sEMG信號的耦合分析，揭示神經肌肉系統的同步耦合特性和功能聯系，探索中風運動功能障礙的產生機制，為康復過程中的運動功能評價提供依據。
1 信號采集與預處理
1.1 實驗對象與數據采集
實驗以秦皇島市第一醫院8名中風恢復期具有上肢運動功能障礙的患者為研究對象(相關信息見表 1)，分別采集中風患者患側、健側上肢“屈”“伸”動作中的sEMG信號。參與實驗的受試者需滿足以下要求：①中風患者患肢可以完成(或部分完成)自主運動；②患者健側肢體運動功能正常；③沒有疼痛或其他傷病因素影響上肢運動和力量；④患者能夠理解和完成實驗要求；⑤所有受試者均自愿參加此實驗測試。
表1
被試中風患者的相關信息
Table1.
Clinical and demographic details of the patients
sEMG信號采集采用NI公司的USB-6251DAQ板卡，放大器低頻截止頻率為5 Hz，高頻截止頻率為500 Hz，采樣頻率為1 000 Hz。本實驗中在受試者上肢BB、TB肌腹處各粘貼兩個AgCl電極，構成差分輸入。參考電極貼于腕處，將導聯線適當固定以減少動作過程中導聯線晃動的干擾。
實驗中，要求受試者坐在測試椅上，肘關節固定，前臂垂直于測試臺，以保證實驗規范化和準確性。實驗開始時，先采集每個受試者健側上肢屈伸動作下的sEMG信號，隨后采集患側上肢自主屈伸運動下的sEMG信號，要求被試盡最大限度募集運動單位加入協同收縮。采集過程中，受試者可以通過信號采集界面視覺反饋提示調整動作，盡量保持運動模式的標準化和一致性。每次“屈”“伸”動作各持續3 s，每4次動作為一組，每個受試者共采集4組，每組運動后休息10 s，放松手臂肌肉，避免肌肉疲勞。
1.2 數據預處理
為了獲取更為有效的sEMG信號特征為進一步的耦合分析做準備，需對原始sEMG信號進行一系列預處理。首先利用4階Butterworth帶通濾波器(通帶范圍為5~200 Hz)提取信號有用成分，再利用4階Butterworth高通濾波器(截止頻率為53 Hz)對信號進行高通濾波，減少頻率約為5~10 Hz的運動偽跡[15]。隨后對濾波后的sEMG信號進行全波整流，其作用之一是用來增強sEMG信號的發放率信息，突出一致性峰值，提供運動單位分組發放的暫態模式，對深入研究肌肉耦合具有意義[14, 16]；二是將信號近似轉化為平穩零均值時間序列，更適用于基于傅里葉變換的一致性分析[8]。
1.3 激活水平估計
均方根(root mean square，RMS)值可以描述肌肉的肌肉激活程度[14]。計算sEMG信號的RMS值來估計每塊肌肉的激活水平(activation level, AL)，以此評價受試者任務完成情況[17]。不同動作狀態下肌肉的激活水平定義如下：
| $ {\rm{AL}}=\sum\limits_{i=1}^n {{R_i}/} \sum\limits_{i=1}^n {R_i^{NM}}, $ |
其中n為動作sEMG信號段數，Ri為第i段sEMG信號的RMS值，RNM為靜息狀態下sEMG信號的RMS值。分別計算每個受試者每組測試健側、患側上肢在“屈”“伸”狀態下的肌肉的激活水平，并取平均值。
2 肌間一致性分析
一致性分析方法可以定量評價兩個在解剖學上獨立的神經區域間相互的生理相關性程度，適合分析處理信號間的耦合同步特性，該方法廣泛應用于生物電信號同步特征研究中。本文基于一致性分析方法，引入一致性顯著性評價指標以及時域累積密度估計方法分析sEMG之間的一致性特征。
2.1 頻域一致性分析
一致性可以從頻域上量化雙變量信號的線性相關關系，即對于兩信號中相同節律成分進行相位穩定性測量，其測量方法定義為兩信號間的互譜密度函數對各自譜密度函數的歸一化。該方法廣泛應用于sEMG同步特征研究，能夠有效地反映肌肉間的耦合關系[3]。為了描述頻域內表面信號間的相互關系，對兩通道整流后的sEMG信號在整個頻段內進行一致性分析。
由于sEMG信號間的一致性較低，容易導致一致性估計的誤差，該誤差大小與進行譜估計的信號樣本數量相關，樣本量越大誤差越小。為改善一致性計算中譜估計器的性能，本文采用滑動平均技術，將sEMG信號分為一系列等長的時段，對每個時段進行譜估計后進行平均(Welch)方法，兩個sEMG信號在頻率λ處的一致性可以由下公式計算：
| $ {\left| {{R_{xy}}\left(\lambda \right)} \right|^2}=\frac{{{{\left| {{f_{xy}}\left(\lambda \right)} \right|}^2}}}{{{f_{xx}}\left(\lambda \right){f_{yy}}\left(\lambda \right)}}, $ |
其中fxx(λ)、fyy(λ)分別為兩通道sEMG信號x和y的自譜密度函數，fxy(λ)表示x與y的互譜密度函數，｜Rxy(λ)｜2描述兩信號在頻域內的線性相關性，其取值范圍從0 (兩信號間無線性關系)到1(兩個信號完全相關)。
2.2 一致性顯著性評價指標
兩信號一致性程度可用顯著一致性閾值S描述[18]，計算公式如下：
| $ S\left(a \right)=1-{\left({1-a} \right)^{1/\left({n-1} \right)}}, $ |
其中n為參與譜估計的數據段數目，α為置信水平(α=0.95)，超過顯著性閾值(即P < α)時兩sEMG相干顯著。
為進行進一步的一致性對比及統計分析，利用Fisher Z變換將IMC進行標準化和近似變異穩定化處理，使Z變換后的一致性近似服從標準差為1、均值與n有關的正態分布。Fisher Z變換一致性(Fisher’s Z-transformed coherence, ZTC)[18]，定義如下：
| $ \begin{array}{*{20}{c}} {Z\left(a \right)=\frac{1}{{2\sqrt {2n} }}\ln \frac{{1 + \sqrt {{{\left| {{R_{xy}}\left(\lambda \right)} \right|}^2}} }}{{1-\sqrt {{{\left| {{R_{xy}}\left(\lambda \right)} \right|}^2}} }}}\\ {=\frac{{{\rm{arctanh}}\left({\sqrt {{{\left| {{R_{xy}}\left(\lambda \right)} \right|}^2}} } \right)}}{{\sqrt {2n} }}, } \end{array} $ |
其中｜Rxy(λ)｜2為一致性值。
為了比較組間一致性，描述不同情況下各頻段IMC的統計差異，本文計算Z變換顯著一致性面積指標AZ，即ZTC曲線以下顯著一致性閾值以上的面積，由下式計算：
| $ {{A}_{Z}}=\sum\limits_{\lambda }{\vartriangle \lambda \left(Z\left(\lambda \right)-S \right)}, $ |
其中Δλ為頻率分辨率，AZ的數值越大表示特定頻段內IMC越顯著。
2.3 時域累積密度估計
為了為兩個信號相關性提供時域上的測量，對頻域一致性進行補充，本文計算兩個信號的時域累積密度估計(cumulant density estimates，CDE)[19]。CDE是一致性頻譜的傅里葉反變換，與一致性相似，CDE可以用來刻畫兩sEMG信號在時域中的相關特征，為信號間的統計相關性提供了一般性測量，其計算公式如下：
| $ {{q}_{xy}}\left(u \right)=\int\limits_{-\pi }^{\pi }{{{f}_{xy}}\left(\lambda \right){{\text{e}}^{i\lambda u}}\text{d}\lambda } $ |
其中fxy(λ)表示x與y的互譜密度函數，u為時間延遲。如果一個信號相對于另一個信號是獨立的(兩個信號不相關)，那么CDE為0；時間延遲u=0時，CDE值越大，說明兩信號同步性越強。CDE體現了兩個信號在時間延遲為u時的相關性。
3 實驗結果與討論
按照1.1節中的實驗過程采集8名中風患者運動中的健側和患側BB和TB的sEMG信號，每個被試者每4次動作為一組，每個受試者共采集4組。截取每組數據第一次動作開始時刻之前的1 s到第四次動作結束后1 s的數據，并按1.2節中方法進行預處理。
3.1 中風患者健側和患側屈伸運動下肌間一致性差異
首先根據2.1節中所述方法計算每組BB和TB的sEMG信號間的一致性，設定互不重疊的大小為512個點的FFT窗口計算兩個sEMG信號的自譜及互譜，進而計算一致性，由于采樣頻率為1 000 Hz，可得到1.96 Hz的頻率分辨率。隨后計算IMC相應的時域累積密度估計。為描述各頻段內的IMC程度，分別計算delta(0.5~4 Hz)、theta(4~8 Hz)、alpha(8~15 Hz)、beta(15~30 Hz)、gamma(30~45 Hz)、piper(35~60 Hz) 6個頻段內的Z變換顯著一致性面積AZ。以患者S8的一組數據為例，圖 1為S8患側、健側肘關節相同屈伸運動下BB和TB肌肉的sEMG信號、頻譜、IMC以及CDE對比圖。圖 1(c)中虛線為顯著一致性閾值S，圖 1(d)中的實線代表 95%的置信區間。
圖1
中風患者患側與健側IMC對比
(a)患側與健側sEMG信號；(b)患側與健側sEMG信號頻譜；(c)患側與健側IMC；(d)患側與健側累積密度估計
Figure1. Comparison of IMC between affected and intact upper limbs(a) sEMG of affected and intact upper limbs; (b) spectrum of affected and intact upper limbs; (c) intermuscular coherence of affected and intact upper limbs; (d) cumulant density estimates of affected and intact upper limbs
分析圖 1(c)可以發現，患者患側IMC在低頻段內相比于健側IMC較低，并且通過圖 1(d)發現，患側累積密度函數的中心峰值沒有健側明顯，說明患側BB-TB肌肉間耦合較弱，一致性沒有健側顯著。
為了分析這一現象是否具有普遍性，分別將8位受試者患側、健側肘部屈伸運動下4組BB-TB IMC進行疊加平均，所得結果如圖 2所示。由于beta和gamma兩個頻段的CMC與神經運動系統的功能狀態有關，因此重點分析這兩個頻段。以圖 2中患者S1健側、患側BB-TB IMC對比分析為例，患者患側beta頻段IMC數值分布明顯低于健側，gamma頻段健側、患側IMC差異并不明顯。圖中其他患者健側與患側IMC對比情況與此相似。同時可以觀察到部分患者(S1、S2、S3、S4和S7)alpha頻段內患側IMC也明顯低于健側。
圖2
不同受試者患側與健側肌間一致性對比
Figure2.
Comparison of intermuscular coherence between affected and intact sides
為了更好地比較健側、患側各頻段IMC的統計差異，利用式(5)計算所有測試組所有被試者各頻段的Z變換一致性顯著面積指標AZ，結果如圖 3(a)所示。并對健側、患側兩組各頻段的logAZ進行重復測量方差分析，檢驗組內、組間因素對觀測變量的影響，結果如圖 3(b)所示。通過分析得到以下結論：組間變量(患側組與健側組)作為組間因素對dalta、theta、gamma和piper頻段的logAZ沒有顯著影響(P > 0.05)，而組間變量對alpha、beta頻段的logAZ值存在顯著影響(P < 0.05)；同時，可以判斷組內因素alpha、beta頻段的logAZ值的主效應顯著(P < 0.05)。隨后分別對兩組各頻段的logAZ值進行兩獨立樣本t檢驗，分析發現患側組與健側組alpha、beta頻段的logAZ值具有顯著差異(P < 0.05)。以上統計分析結果可以說明中風患者患側beta段IMC與健側相比明顯缺失。
圖3
患側與健側各頻段內Z變換一致性顯著面積指標AZ對比
(a)所有受試者全部試驗的
(a)
肌肉活動在beta頻段內的振蕩被認為來源于皮質脊髓束(或椎體束)的皮層共同驅動，beta頻段的IMC依賴于沒有前角細胞缺失的完整的皮質脊髓束。在相關的其它研究中，Fisher等[13]發現運動神經損傷患者IMC相對于健康受試者在beta(15~30 Hz)頻段內有缺失；Hansen等[20]發現脊髓損傷患者10~20 Hz頻段內的IMC減小；Farmer等[21]、Norton等[22]也先后發現中風后beta頻段內的肌肉耦合消失。這與本文所觀察到的現象基本一致。正常條件下，皮質脊髓系統對控制獨立、協調的運動具有重要作用，而中風患者患側運動執行過程中的肌肉功能性耦合減弱，可能是由皮質脊髓通路信息阻斷導致的共同驅動丟失造成的[14]，這也是導致中風患者運動控制紊亂、運動執行能力降低的原因之一。
由于IMC可以用來反映CMC的變化趨勢，因此我們推測中風患者CMC在beta頻段也有減小的趨勢。大腦皮層接收、整合來自身體各處的感覺信息，并依靠皮層的錐體束等下行性投射神經元將控制信息傳遞出去。不同區域的投射神經元其下行性投射靶區不同，而皮層運動區投射到脊髓。在皮層振蕩的影響下運動神經元有節律地發放從而調節肌肉活動。已有研究表明，beta頻段的腦肌電振蕩耦合與皮質脊髓通路有關[6]。由于中風患者部分下行通路損傷，皮質脊髓投射神經元的減少導致對剩余的未損傷的下行通路(如紅核脊髓通路、前庭脊髓束通路、網狀脊髓束通路等)依賴性的增加[23]。這些并行的通路在脊髓水平具有大量的分支，對它們的依賴增強會導致正常情況下不會同時激活的肌肉出現異常協同收縮，或者產生錯誤的肌肉激活模式[14, 24]，反映在CMC上即表現為beta頻段一致性的相對缺失。這一結論與我們根據實驗結果所作的推測相符合，也給IMC缺失提供了皮質脊髓層面上的解釋。
然而，IMC與CMC稍有不同，拮抗肌肉運動神經元的共同驅動不僅來自于皮質脊髓中樞控制中心，也來自運動神經末梢[25]，因此IMC與突觸間的轉化有關。突觸間的轉化包括：①突觸前輸入到alpha運動神經元發放模式的轉化；②沿著肌纖維的運動電位的傳播，構成運動單位動作電位(motor unit action potential, MUAP)。MUAP的變化和發放率影響IMC水平，由于中風患者運動單位發放率發生變化，導致IMC產生異常耦合模式。因此，IMC也為運動有關通路的突觸轉化提供了定量描述。
3.2 中風患者健側和患側運動時具有時間分辨率的肌間一致性分析
為了驗證IMC與任務間的聯系，截取部分動作段sEMG信號并將信號分成連續的時間段，時間窗長設為200點，窗滑動步長200點，可得時間分辨率為200 ms。對每一個時間段上的信號分別進行一致性分析，設定FFT窗寬128點，重疊64點，頻率分辨率為7.69 Hz，具體計算過程如2.1節所述，得到具有時間分辨率的動態IMC時頻圖，以S5患者為例，健側、患側動態IMC時頻圖如圖 4所示。
圖4
患側與健側動態一致性對比
Figure4.
Comparison of dynamic intermuscular coherence between affected and intact sides
分析圖 4可知，患者患側完成指定屈伸任務時，無論動態還是靜態狀態下其IMC總體低于健側，這可能是由中風患者皮質脊髓通路共同驅動丟失而造成的。此外，健側、患側動態任務輸出下beta段IMC較小而維持靜態動作時beta段IMC最為明顯。由于當注意力分散時beta段IMC較低，當提高運動任務的精度時，一致性也隨之升高[26]。因此該現象原因可推測為：完成靜態任務時保持動作狀態需要集中注意力，運動任務精度高，因此beta段IMC較高。上述實驗結果說明：①beta頻段的振蕩是運動皮層維持穩定運動輸出的表現；②IMC會隨著任務的不同而改變，這與已有研究結果一致[5, 26]。
3.3 肌間一致性影響因素及與臨床評價指標的關系分析
首先討論患者的sEMG能量是否對耦合特征分析造成影響。為此本文計算了兩通道sEMG頻譜面積值作為描述信號能量強度的指標，為了便于進一步的統計分析，我們將得到的頻譜面積值進行對數化處理使之服從正態分布。然后，進一步討論IMC與肌肉激活水平的關系。如1.3節所述計算不同動作狀態下各肌肉的激活水平AL，8名受試者患側、健側上肢“屈”“伸”任務下BB、TB肌肉激活水平如表 2所示，分析BB、TB肌肉的AL之比與IMC的關系。此外，我們還將考慮臨床評價指標Brunnstrom分期、Fugle-Meyer評分與肌肉耦合模式的關系。
表2
患側與健側上肢屈伸任務下BB和TB肌肉激活水平(為了考察以上因素對于IMC的影響，我們分別以Brunnstrom分期、Fugle-Meyer評分因素對患側各頻段的AZ進行方差分析；利用Spearman相關分析分別計算患側、健側屈伸狀態下BB、TB肌肉AL比值和BB、TB的sEMG能量與其相應的同側上肢各頻段AZ的相關性，結果如表 3所示。
表3
不同因素影響下IMC的統計分析表
Table3.
Statistic analysis of the influences on intermuscular coherence under influences of different factors
從表中可以看出Brunnstrom分期并不是IMC的主要影響因素；Fugle-Meyer評分對患側beta段的AZ影響顯著(P < 0.05)，對其它頻段AZ影響并不顯著；屈伸狀態下BB、TB肌肉AL比值與sEMG信號alpha、beta頻段一致性的大小顯著相關(P < 0.05)，而對其它頻段AZ影響較小；BB、TB的sEMG能量與各頻段AZ相關性均不顯著(P > 0.05)。由上述實驗結果可見，beta段的IMC受肌肉激活水平的影響，可以作為輔助Fugle-Meyer量表評價患者運動功能狀態的一個有效指標。
4 結論
本文基于一致性分析方法，對比研究中風患者運動過程中健側、患側上肢拮抗肌間的一致性特征，挖掘患側異常耦合模式，推測神經肌肉系統的振蕩聯系，實現對患者肌肉功能狀態的評價；在此基礎上利用Z變換顯著一致性面積指標量化描述健側、患側sEMG信號在各功能頻帶內的耦合性差異；進一步利用動態IMC分析探究IMC是否與運動任務有關。通過實驗研究得到以下結論：beta頻段IMC會隨著任務的不同而改變，而中風運動功能障礙患者在患側運動過程中beta頻段IMC相對于健側存在缺失，并且beta段的IMC與Fugle-Meyer量表具有相關性，從生物電信號角度為患者運動功能狀態評價提供了一個有效的輔助指標。因此本文提出的IMC分析方法及指標能夠反映神經肌肉運動系統的同步耦合特性和功能聯系，有助于探索中風運動功能障礙產生的病理機制，并建立基于sEMG信號的康復運動中的肌肉功能狀態評價標準。
然而，本文參與實驗的受試者數量較少，因此關于中風患者IMC的分析仍需進一步完善，對神經肌肉異常耦合和運動功能障礙產生機制的推測也有待繼續研究。下一步將采集一定量的中風患者與健康受試者相同運動模式下的sEMG信號數據，進行患者與健康人IMC的對比分析，進一步考察中風運動功能障礙患者的異常神經振蕩和神經肌肉狀態。
引言
在人體運動的神經控制過程中，神經元的振蕩耦合和同步放電發揮著至關重要的作用。運動神經系統通過神經振蕩來傳遞運動控制信息[1]，神經中樞通過與各區域間的相互作用實現功能調節和整合，主要體現在電生理信號間不同層次的耦合(同步)現象，如腦電(electroencephalogram, EEG)-肌電(electromyogram, EMG)耦合、EMG-EMG耦合等。中樞神經損傷患者由于神經損傷阻礙了神經振蕩的傳導從而造成運動功能障礙。近年來，運用電生理信號耦合分析考察中樞神經損傷患者的異常神經振蕩和神經肌肉狀態已成為運動醫學、康復工程等領域的研究熱點。通過耦合特征分析可以反映運動控制過程中神經肌肉間的功能聯系，為理解運動控制過程及運動障礙的病理機制提供了理論基礎，也為神經康復運動的功能狀態評價提供了新方法。
EEG信號和表面肌電(surface electromyogram, sEMG)信號，二者分別作為人體指令中心和執行末端的代表性信號，在機體運動功能領域有著非常重要的研究意義。自從Conway等[2]發現運動過程中相關肌肉肌電和皮層腦電之間具有頻率成分的耦合性，基于腦肌電一致性(EEG-EMG/corticomuscular coherence, CMC)對大腦運動皮層與肌肉功能聯系的研究相繼展開，并應用于探索中風、帕金森癥等中樞神經系統相關疾病病理機制的研究中。研究表明，運動過程中CMC主要表現在beta頻段和gamma頻段，beta頻段的腦肌電振蕩代表了從初級運動皮層到脊髓運動神經元的傳遞過程[3]，而gamma頻段振蕩體現與認知功能相關的腦皮層信息整合過程[4]。因此，通過CMC分析能夠體現中樞神經與肌肉之間的振蕩聯系，進而評價運動神經系統功能，甚至指導康復過程。然而CMC分析在臨床應用中仍具有困難，原因主要有三個：①臨床中難以獲得高質量的EEG信號；②腦肌電同步采集一直是難以解決的問題；③CMC特征尤其是患者的一致性特征并不穩定。值得提出的是，Brown等研究者發現肌間一致性(EMG-EMG/intermuscular coherence, IMC)與CMC具有相似的模式，說明同一個受試者同一個任務下的CMC和IMC均起源于皮質脊髓通路[5-7]，因此可以用IMC來反映CMC的變化趨勢[8]。同時，他們發現IMC還可以反映脊髓運動神經元共同的突觸前驅動，使運動單位同步發放從而激發肌肉活性[3]。另一方面，在臨床環境下sEMG的采集更為安全便捷，信號可靠性也更強。綜上所述，利用IMC分析神經肌肉振蕩聯系、評價運動神經系統功能更具有臨床實用價值。
近年來，基于IMC分析運動過程中肌肉間的耦合模式、挖掘潛在的神經運動控制機制的研究相繼開展，如Baker等[9]、Nroton等[10]分別發現當肌肉持續收縮時的IMC最大；Nashner[11]發現拮抗肌最大自主收縮時的IMC在10 Hz左右最明顯；Farmer等[7]發現等長收縮運動下運動單位同步主要發生在1~12 Hz和16~32 Hz。此外，還有部分研究者將IMC分析引入運動障礙病理機制及運動功能評價的研究中，如Grosse等[12]發現肌張力障礙患者與正常受試者的IMC在theta(4~7 Hz)頻段具有顯著差異；Fisher等[13]發現beta(15~30 Hz)頻段IMC可以作為評價運動神經損傷患者上肢運動功能紊亂的指標；Kisiel-Sajewicz等[14]發現中風患者進行抓取運動時肌肉功能性耦合在0~3.9 Hz和4~7.9 Hz有缺失。然而關于中風康復運動中患側、健側肌肉耦合差異方面的深入理論分析還不多見，是一個值得進一步研究的問題。
為此，本文首先分別進行中風運動功能障礙患者肘部屈伸動作任務下健側與患側的上肢sEMG信號同步采集實驗。在此基礎上基于一致性分析方法，對整流后的肱二頭肌(biceps brachii, BB)和肱三頭肌(triceps brachii, TB) sEMG信號進行頻域一致性分析和時域累積密度估計分析，針對各特定頻段研究肌肉間的耦合模式和強度，以此推測神經肌肉系統的振蕩聯系，并利用Z變換顯著一致性面積指標量化和比較兩通道sEMG在各頻帶內的耦合性。進一步利用具有時間分辨率的動態一致性分析，探究IMC如何隨著實驗任務變化而變化，最后利用統計方法對比分析患者患側與健側IMC的差異。通過對中風患者sEMG信號的耦合分析，揭示神經肌肉系統的同步耦合特性和功能聯系，探索中風運動功能障礙的產生機制，為康復過程中的運動功能評價提供依據。
1 信號采集與預處理
1.1 實驗對象與數據采集
實驗以秦皇島市第一醫院8名中風恢復期具有上肢運動功能障礙的患者為研究對象(相關信息見表 1)，分別采集中風患者患側、健側上肢“屈”“伸”動作中的sEMG信號。參與實驗的受試者需滿足以下要求：①中風患者患肢可以完成(或部分完成)自主運動；②患者健側肢體運動功能正常；③沒有疼痛或其他傷病因素影響上肢運動和力量；④患者能夠理解和完成實驗要求；⑤所有受試者均自愿參加此實驗測試。
表1
被試中風患者的相關信息
Table1.
Clinical and demographic details of the patients
sEMG信號采集采用NI公司的USB-6251DAQ板卡，放大器低頻截止頻率為5 Hz，高頻截止頻率為500 Hz，采樣頻率為1 000 Hz。本實驗中在受試者上肢BB、TB肌腹處各粘貼兩個AgCl電極，構成差分輸入。參考電極貼于腕處，將導聯線適當固定以減少動作過程中導聯線晃動的干擾。
實驗中，要求受試者坐在測試椅上，肘關節固定，前臂垂直于測試臺，以保證實驗規范化和準確性。實驗開始時，先采集每個受試者健側上肢屈伸動作下的sEMG信號，隨后采集患側上肢自主屈伸運動下的sEMG信號，要求被試盡最大限度募集運動單位加入協同收縮。采集過程中，受試者可以通過信號采集界面視覺反饋提示調整動作，盡量保持運動模式的標準化和一致性。每次“屈”“伸”動作各持續3 s，每4次動作為一組，每個受試者共采集4組，每組運動后休息10 s，放松手臂肌肉，避免肌肉疲勞。
1.2 數據預處理
為了獲取更為有效的sEMG信號特征為進一步的耦合分析做準備，需對原始sEMG信號進行一系列預處理。首先利用4階Butterworth帶通濾波器(通帶范圍為5~200 Hz)提取信號有用成分，再利用4階Butterworth高通濾波器(截止頻率為53 Hz)對信號進行高通濾波，減少頻率約為5~10 Hz的運動偽跡[15]。隨后對濾波后的sEMG信號進行全波整流，其作用之一是用來增強sEMG信號的發放率信息，突出一致性峰值，提供運動單位分組發放的暫態模式，對深入研究肌肉耦合具有意義[14, 16]；二是將信號近似轉化為平穩零均值時間序列，更適用于基于傅里葉變換的一致性分析[8]。
1.3 激活水平估計
均方根(root mean square，RMS)值可以描述肌肉的肌肉激活程度[14]。計算sEMG信號的RMS值來估計每塊肌肉的激活水平(activation level, AL)，以此評價受試者任務完成情況[17]。不同動作狀態下肌肉的激活水平定義如下：
| $ {\rm{AL}}=\sum\limits_{i=1}^n {{R_i}/} \sum\limits_{i=1}^n {R_i^{NM}}, $ |
其中n為動作sEMG信號段數，Ri為第i段sEMG信號的RMS值，RNM為靜息狀態下sEMG信號的RMS值。分別計算每個受試者每組測試健側、患側上肢在“屈”“伸”狀態下的肌肉的激活水平，并取平均值。
2 肌間一致性分析
一致性分析方法可以定量評價兩個在解剖學上獨立的神經區域間相互的生理相關性程度，適合分析處理信號間的耦合同步特性，該方法廣泛應用于生物電信號同步特征研究中。本文基于一致性分析方法，引入一致性顯著性評價指標以及時域累積密度估計方法分析sEMG之間的一致性特征。
2.1 頻域一致性分析
一致性可以從頻域上量化雙變量信號的線性相關關系，即對于兩信號中相同節律成分進行相位穩定性測量，其測量方法定義為兩信號間的互譜密度函數對各自譜密度函數的歸一化。該方法廣泛應用于sEMG同步特征研究，能夠有效地反映肌肉間的耦合關系[3]。為了描述頻域內表面信號間的相互關系，對兩通道整流后的sEMG信號在整個頻段內進行一致性分析。
由于sEMG信號間的一致性較低，容易導致一致性估計的誤差，該誤差大小與進行譜估計的信號樣本數量相關，樣本量越大誤差越小。為改善一致性計算中譜估計器的性能，本文采用滑動平均技術，將sEMG信號分為一系列等長的時段，對每個時段進行譜估計后進行平均(Welch)方法，兩個sEMG信號在頻率λ處的一致性可以由下公式計算：
| $ {\left| {{R_{xy}}\left(\lambda \right)} \right|^2}=\frac{{{{\left| {{f_{xy}}\left(\lambda \right)} \right|}^2}}}{{{f_{xx}}\left(\lambda \right){f_{yy}}\left(\lambda \right)}}, $ |
其中fxx(λ)、fyy(λ)分別為兩通道sEMG信號x和y的自譜密度函數，fxy(λ)表示x與y的互譜密度函數，｜Rxy(λ)｜2描述兩信號在頻域內的線性相關性，其取值范圍從0 (兩信號間無線性關系)到1(兩個信號完全相關)。
2.2 一致性顯著性評價指標
兩信號一致性程度可用顯著一致性閾值S描述[18]，計算公式如下：
| $ S\left(a \right)=1-{\left({1-a} \right)^{1/\left({n-1} \right)}}, $ |
其中n為參與譜估計的數據段數目，α為置信水平(α=0.95)，超過顯著性閾值(即P < α)時兩sEMG相干顯著。
為進行進一步的一致性對比及統計分析，利用Fisher Z變換將IMC進行標準化和近似變異穩定化處理，使Z變換后的一致性近似服從標準差為1、均值與n有關的正態分布。Fisher Z變換一致性(Fisher’s Z-transformed coherence, ZTC)[18]，定義如下：
| $ \begin{array}{*{20}{c}} {Z\left(a \right)=\frac{1}{{2\sqrt {2n} }}\ln \frac{{1 + \sqrt {{{\left| {{R_{xy}}\left(\lambda \right)} \right|}^2}} }}{{1-\sqrt {{{\left| {{R_{xy}}\left(\lambda \right)} \right|}^2}} }}}\\ {=\frac{{{\rm{arctanh}}\left({\sqrt {{{\left| {{R_{xy}}\left(\lambda \right)} \right|}^2}} } \right)}}{{\sqrt {2n} }}, } \end{array} $ |
其中｜Rxy(λ)｜2為一致性值。
為了比較組間一致性，描述不同情況下各頻段IMC的統計差異，本文計算Z變換顯著一致性面積指標AZ，即ZTC曲線以下顯著一致性閾值以上的面積，由下式計算：
| $ {{A}_{Z}}=\sum\limits_{\lambda }{\vartriangle \lambda \left(Z\left(\lambda \right)-S \right)}, $ |
其中Δλ為頻率分辨率，AZ的數值越大表示特定頻段內IMC越顯著。
2.3 時域累積密度估計
為了為兩個信號相關性提供時域上的測量，對頻域一致性進行補充，本文計算兩個信號的時域累積密度估計(cumulant density estimates，CDE)[19]。CDE是一致性頻譜的傅里葉反變換，與一致性相似，CDE可以用來刻畫兩sEMG信號在時域中的相關特征，為信號間的統計相關性提供了一般性測量，其計算公式如下：
| $ {{q}_{xy}}\left(u \right)=\int\limits_{-\pi }^{\pi }{{{f}_{xy}}\left(\lambda \right){{\text{e}}^{i\lambda u}}\text{d}\lambda } $ |
其中fxy(λ)表示x與y的互譜密度函數，u為時間延遲。如果一個信號相對于另一個信號是獨立的(兩個信號不相關)，那么CDE為0；時間延遲u=0時，CDE值越大，說明兩信號同步性越強。CDE體現了兩個信號在時間延遲為u時的相關性。
3 實驗結果與討論
按照1.1節中的實驗過程采集8名中風患者運動中的健側和患側BB和TB的sEMG信號，每個被試者每4次動作為一組，每個受試者共采集4組。截取每組數據第一次動作開始時刻之前的1 s到第四次動作結束后1 s的數據，并按1.2節中方法進行預處理。
3.1 中風患者健側和患側屈伸運動下肌間一致性差異
首先根據2.1節中所述方法計算每組BB和TB的sEMG信號間的一致性，設定互不重疊的大小為512個點的FFT窗口計算兩個sEMG信號的自譜及互譜，進而計算一致性，由于采樣頻率為1 000 Hz，可得到1.96 Hz的頻率分辨率。隨后計算IMC相應的時域累積密度估計。為描述各頻段內的IMC程度，分別計算delta(0.5~4 Hz)、theta(4~8 Hz)、alpha(8~15 Hz)、beta(15~30 Hz)、gamma(30~45 Hz)、piper(35~60 Hz) 6個頻段內的Z變換顯著一致性面積AZ。以患者S8的一組數據為例，圖 1為S8患側、健側肘關節相同屈伸運動下BB和TB肌肉的sEMG信號、頻譜、IMC以及CDE對比圖。圖 1(c)中虛線為顯著一致性閾值S，圖 1(d)中的實線代表 95%的置信區間。
圖1
中風患者患側與健側IMC對比
(a)患側與健側sEMG信號；(b)患側與健側sEMG信號頻譜；(c)患側與健側IMC；(d)患側與健側累積密度估計
Figure1. Comparison of IMC between affected and intact upper limbs(a) sEMG of affected and intact upper limbs; (b) spectrum of affected and intact upper limbs; (c) intermuscular coherence of affected and intact upper limbs; (d) cumulant density estimates of affected and intact upper limbs
分析圖 1(c)可以發現，患者患側IMC在低頻段內相比于健側IMC較低，并且通過圖 1(d)發現，患側累積密度函數的中心峰值沒有健側明顯，說明患側BB-TB肌肉間耦合較弱，一致性沒有健側顯著。
為了分析這一現象是否具有普遍性，分別將8位受試者患側、健側肘部屈伸運動下4組BB-TB IMC進行疊加平均，所得結果如圖 2所示。由于beta和gamma兩個頻段的CMC與神經運動系統的功能狀態有關，因此重點分析這兩個頻段。以圖 2中患者S1健側、患側BB-TB IMC對比分析為例，患者患側beta頻段IMC數值分布明顯低于健側，gamma頻段健側、患側IMC差異并不明顯。圖中其他患者健側與患側IMC對比情況與此相似。同時可以觀察到部分患者(S1、S2、S3、S4和S7)alpha頻段內患側IMC也明顯低于健側。
圖2
不同受試者患側與健側肌間一致性對比
Figure2.
Comparison of intermuscular coherence between affected and intact sides
為了更好地比較健側、患側各頻段IMC的統計差異，利用式(5)計算所有測試組所有被試者各頻段的Z變換一致性顯著面積指標AZ，結果如圖 3(a)所示。并對健側、患側兩組各頻段的logAZ進行重復測量方差分析，檢驗組內、組間因素對觀測變量的影響，結果如圖 3(b)所示。通過分析得到以下結論：組間變量(患側組與健側組)作為組間因素對dalta、theta、gamma和piper頻段的logAZ沒有顯著影響(P > 0.05)，而組間變量對alpha、beta頻段的logAZ值存在顯著影響(P < 0.05)；同時，可以判斷組內因素alpha、beta頻段的logAZ值的主效應顯著(P < 0.05)。隨后分別對兩組各頻段的logAZ值進行兩獨立樣本t檢驗，分析發現患側組與健側組alpha、beta頻段的logAZ值具有顯著差異(P < 0.05)。以上統計分析結果可以說明中風患者患側beta段IMC與健側相比明顯缺失。
圖3
患側與健側各頻段內Z變換一致性顯著面積指標AZ對比
(a)所有受試者全部試驗的
(a)
肌肉活動在beta頻段內的振蕩被認為來源于皮質脊髓束(或椎體束)的皮層共同驅動，beta頻段的IMC依賴于沒有前角細胞缺失的完整的皮質脊髓束。在相關的其它研究中，Fisher等[13]發現運動神經損傷患者IMC相對于健康受試者在beta(15~30 Hz)頻段內有缺失；Hansen等[20]發現脊髓損傷患者10~20 Hz頻段內的IMC減小；Farmer等[21]、Norton等[22]也先后發現中風后beta頻段內的肌肉耦合消失。這與本文所觀察到的現象基本一致。正常條件下，皮質脊髓系統對控制獨立、協調的運動具有重要作用，而中風患者患側運動執行過程中的肌肉功能性耦合減弱，可能是由皮質脊髓通路信息阻斷導致的共同驅動丟失造成的[14]，這也是導致中風患者運動控制紊亂、運動執行能力降低的原因之一。
由于IMC可以用來反映CMC的變化趨勢，因此我們推測中風患者CMC在beta頻段也有減小的趨勢。大腦皮層接收、整合來自身體各處的感覺信息，并依靠皮層的錐體束等下行性投射神經元將控制信息傳遞出去。不同區域的投射神經元其下行性投射靶區不同，而皮層運動區投射到脊髓。在皮層振蕩的影響下運動神經元有節律地發放從而調節肌肉活動。已有研究表明，beta頻段的腦肌電振蕩耦合與皮質脊髓通路有關[6]。由于中風患者部分下行通路損傷，皮質脊髓投射神經元的減少導致對剩余的未損傷的下行通路(如紅核脊髓通路、前庭脊髓束通路、網狀脊髓束通路等)依賴性的增加[23]。這些并行的通路在脊髓水平具有大量的分支，對它們的依賴增強會導致正常情況下不會同時激活的肌肉出現異常協同收縮，或者產生錯誤的肌肉激活模式[14, 24]，反映在CMC上即表現為beta頻段一致性的相對缺失。這一結論與我們根據實驗結果所作的推測相符合，也給IMC缺失提供了皮質脊髓層面上的解釋。
然而，IMC與CMC稍有不同，拮抗肌肉運動神經元的共同驅動不僅來自于皮質脊髓中樞控制中心，也來自運動神經末梢[25]，因此IMC與突觸間的轉化有關。突觸間的轉化包括：①突觸前輸入到alpha運動神經元發放模式的轉化；②沿著肌纖維的運動電位的傳播，構成運動單位動作電位(motor unit action potential, MUAP)。MUAP的變化和發放率影響IMC水平，由于中風患者運動單位發放率發生變化，導致IMC產生異常耦合模式。因此，IMC也為運動有關通路的突觸轉化提供了定量描述。
3.2 中風患者健側和患側運動時具有時間分辨率的肌間一致性分析
為了驗證IMC與任務間的聯系，截取部分動作段sEMG信號并將信號分成連續的時間段，時間窗長設為200點，窗滑動步長200點，可得時間分辨率為200 ms。對每一個時間段上的信號分別進行一致性分析，設定FFT窗寬128點，重疊64點，頻率分辨率為7.69 Hz，具體計算過程如2.1節所述，得到具有時間分辨率的動態IMC時頻圖，以S5患者為例，健側、患側動態IMC時頻圖如圖 4所示。
圖4
患側與健側動態一致性對比
Figure4.
Comparison of dynamic intermuscular coherence between affected and intact sides
分析圖 4可知，患者患側完成指定屈伸任務時，無論動態還是靜態狀態下其IMC總體低于健側，這可能是由中風患者皮質脊髓通路共同驅動丟失而造成的。此外，健側、患側動態任務輸出下beta段IMC較小而維持靜態動作時beta段IMC最為明顯。由于當注意力分散時beta段IMC較低，當提高運動任務的精度時，一致性也隨之升高[26]。因此該現象原因可推測為：完成靜態任務時保持動作狀態需要集中注意力，運動任務精度高，因此beta段IMC較高。上述實驗結果說明：①beta頻段的振蕩是運動皮層維持穩定運動輸出的表現；②IMC會隨著任務的不同而改變，這與已有研究結果一致[5, 26]。
3.3 肌間一致性影響因素及與臨床評價指標的關系分析
首先討論患者的sEMG能量是否對耦合特征分析造成影響。為此本文計算了兩通道sEMG頻譜面積值作為描述信號能量強度的指標，為了便于進一步的統計分析，我們將得到的頻譜面積值進行對數化處理使之服從正態分布。然后，進一步討論IMC與肌肉激活水平的關系。如1.3節所述計算不同動作狀態下各肌肉的激活水平AL，8名受試者患側、健側上肢“屈”“伸”任務下BB、TB肌肉激活水平如表 2所示，分析BB、TB肌肉的AL之比與IMC的關系。此外，我們還將考慮臨床評價指標Brunnstrom分期、Fugle-Meyer評分與肌肉耦合模式的關系。
表2
患側與健側上肢屈伸任務下BB和TB肌肉激活水平(為了考察以上因素對于IMC的影響，我們分別以Brunnstrom分期、Fugle-Meyer評分因素對患側各頻段的AZ進行方差分析；利用Spearman相關分析分別計算患側、健側屈伸狀態下BB、TB肌肉AL比值和BB、TB的sEMG能量與其相應的同側上肢各頻段AZ的相關性，結果如表 3所示。
表3
不同因素影響下IMC的統計分析表
Table3.
Statistic analysis of the influences on intermuscular coherence under influences of different factors
從表中可以看出Brunnstrom分期并不是IMC的主要影響因素；Fugle-Meyer評分對患側beta段的AZ影響顯著(P < 0.05)，對其它頻段AZ影響并不顯著；屈伸狀態下BB、TB肌肉AL比值與sEMG信號alpha、beta頻段一致性的大小顯著相關(P < 0.05)，而對其它頻段AZ影響較小；BB、TB的sEMG能量與各頻段AZ相關性均不顯著(P > 0.05)。由上述實驗結果可見，beta段的IMC受肌肉激活水平的影響，可以作為輔助Fugle-Meyer量表評價患者運動功能狀態的一個有效指標。
4 結論
本文基于一致性分析方法，對比研究中風患者運動過程中健側、患側上肢拮抗肌間的一致性特征，挖掘患側異常耦合模式，推測神經肌肉系統的振蕩聯系，實現對患者肌肉功能狀態的評價；在此基礎上利用Z變換顯著一致性面積指標量化描述健側、患側sEMG信號在各功能頻帶內的耦合性差異；進一步利用動態IMC分析探究IMC是否與運動任務有關。通過實驗研究得到以下結論：beta頻段IMC會隨著任務的不同而改變，而中風運動功能障礙患者在患側運動過程中beta頻段IMC相對于健側存在缺失，并且beta段的IMC與Fugle-Meyer量表具有相關性，從生物電信號角度為患者運動功能狀態評價提供了一個有效的輔助指標。因此本文提出的IMC分析方法及指標能夠反映神經肌肉運動系統的同步耦合特性和功能聯系，有助于探索中風運動功能障礙產生的病理機制，并建立基于sEMG信號的康復運動中的肌肉功能狀態評價標準。
然而，本文參與實驗的受試者數量較少，因此關于中風患者IMC的分析仍需進一步完善，對神經肌肉異常耦合和運動功能障礙產生機制的推測也有待繼續研究。下一步將采集一定量的中風患者與健康受試者相同運動模式下的sEMG信號數據，進行患者與健康人IMC的對比分析，進一步考察中風運動功能障礙患者的異常神經振蕩和神經肌肉狀態。
