羅萬國 侯文生,# 鄭小林,# 萬小萍 陳海燕 周 平 吳小鷹*

1(重慶大學生物醫(yī)學工程系， 重慶 400044)2(重慶市醫(yī)療電子工程技術研究中心， 重慶 400044)

基于先驗模板的表面肌電信號漸進分解算法研究

羅萬國1侯文生1,2#鄭小林1,2#萬小萍1陳海燕2周 平2吳小鷹1*

1(重慶大學生物醫(yī)學工程系， 重慶 400044)2(重慶市醫(yī)療電子工程技術研究中心， 重慶 400044)

本研究通過將表面肌電信號(sEMG)分解為運動單元動作電位序列(MUAPTs)，來研究神經(jīng)-肌肉系統(tǒng)中運動單元(MU)的募集與發(fā)放模式。針對高收縮力情況下MUAP疊加問題，首先采用FastICA算法和小波包去噪算法對信號進行預處理；然后基于先驗知識構建了4種形態(tài)的可伸縮MUAP模板；最后，采用“先大后小”的漸進識別方式，逐個對MUAP進行自動提取。在此基礎上，還將該算法應用于8名受試者(3組/人)不同手指活動模式下的指淺屈肌多通道(12通道)sEMG分解；單通道分解結果顯示，高力量水平下sEMG中的主體MUAPTs能夠被有效檢測和分類；統(tǒng)計結果證實，隨著力量水平的增加，MUAP的數(shù)目增加；不同大小MUAP比重的變化與活動手指和力量水平具有顯著的相關性。本文的實驗結果，初步驗證了利用先驗模板從sEMG中漸進提取MUAP的可行性，為sEMG分解和進一步研究MU發(fā)放規(guī)律提供于一種新的思路。

表面肌電信號；分解；先驗模板；運動單元動作電位

引言

表面肌電信號(surface electromyography, sEMG)，是由肌肉興奮時募集的多個運動單元(motor unit, MU)產(chǎn)生的動作電位(action potential, AP)，在表面電極處的綜合疊加[1]。由于sEMG具有無創(chuàng)的特點，在假肢控制[2]、神經(jīng)肌肉系統(tǒng)功能異常的輔助診斷[3]、運動康復診斷[4]等方面受到重視。雖然直接提取sEMG的時域、頻域和時頻域特征，可用于描述肌肉活動及其功能狀態(tài)[4-5]，但為了進一步認識MU募集的細節(jié)信息，通過分解sEMG來還原構成它的主體運動單元動作電位序列(motor unit action potential trains, MUAPTs)，在近年來成為了一個重要的研究方向[1, 3, 6-13]。

基于統(tǒng)計學的盲源分離(BSS)算法，尤其是獨立分量分析(ICA)，是比較典型的sEMG分解算法之一。日本坂田大學的Garcia研究小組，采用一系列ICA算法對仿真信號和實測sEMG進行了分解[6-7]，得到了一些與主體MUAPTs相似的分離信號；中國科技大學的神經(jīng)肌肉控制實驗室，將多種BSS算法應用于多通道sEMG分解，得到了4類MUAP信息，并較好地解決了波形疊加問題[11]；重慶大學安媛等利用FastICA算法，從6通道的指淺屈肌sEMG中提取了4種MUAP信息，并初步研究了力量水平對MUAP募集模式的影響[12]。然而，由于ICA算法要求觀測信號中有，且僅有一個高斯源信號，觀測矢量的個數(shù)還不小于源矢量的個數(shù)，故基于ICA的多通道sEMG分解，大多在較低力量水平進行。此外，由于需要反復迭代，其分解效率也有較大局限[11, 14]。

另一種典型的sEMG分解思路，是基于波形檢測和模式識別的方法對sEMG進行分解，其中最具代表性的是De Luca等在插入式肌電信號(iEMG)分解方法的基礎上，提出來新的sEMG分解策略[1]。最新報道顯示，該算法可用于不同肌肉、不同力量水平(最高100% MVC)的sEMG分解，所分解的MU類型數(shù)高達40，準確度也高達97%[8]。但是，也有學者認為其準確性仍需改進[9]。本文結合一些現(xiàn)有sEMG處理與分解算法的優(yōu)點，以MUAP波形結構為先驗知識，設計了一種漸進式sEMG分解算法：通過模板擬合和逐次迭代提取MUAP信息，以提高MUAP識別效率；并通過不同手指活動模式下采集的指淺屈肌sEMG，對算法的有效性進行了初步驗證。

1 漸進分解算法設計

1.1模板建立

MUAP波形被證實大多呈現(xiàn)雙相或三相結構[1, 8, 10, 12]，MU遵循“大小原則”進行募集，后募集的MU具有更大的激活閾值且產(chǎn)生尺寸更大的MUAP。為了使所建立的模板波形更具有代表性，提取出更多MU所發(fā)放的MUAP信息，依據(jù)現(xiàn)有報道的結果[1, 8, 10, 12]，構建了如圖1所示的4種基本形態(tài)的模板波形(Tp1～Tp4)，且每種形態(tài)的模板波形，通過幅度和時間伸縮以代表不同激活閾值MU發(fā)放的MUAP。

為了計算方便，將Hermite-Rodriguez函數(shù)進行歸一化和變形來表示每種形態(tài)MUAP的波形，上述4種模板的近似數(shù)學表達式為

(1)

式中，Tpi(t)表示第i個MUAP模板為了更好地逼近模板波形，變量t的取值范圍?。篬-2.5,2.5]，A為MUAP的幅度伸縮因子，采用式(2)所示一個比值為1/2的等比數(shù)列來代表不同大小MUAP的幅值A0可通過MVC力量水平下采集到的sEMG由式(3)進行估計，N的大小則取決于式(6)中Thrmin的大??；λ1,λ2,λ3,λ4為與MUAP幅值有關的時間伸縮因子，幅值越大，時長相對越長。參考相關文獻[1, 8]，其總變化范圍設為：5～20 ms。

(2)

(3)

式中，peakm為MVC力量水平下采集的sEMG峰值，M為滿足|peak1|-|peakm|≤2σ的最大m值，σ為基線的方差。

1.2分解步驟

sEMG經(jīng)常表現(xiàn)為不同幅值MUAP的相互疊加，幅值大的MUAP波形相對容易被檢測出來。因此，采用如圖2所示的兩層循環(huán)迭代算法，按照“從大到小”的順序，自動對sEMG中的MUAP進行逐個剝離。

步驟1：閾值初始化。比較輸入信號峰值與MVC力量水平下sEMG估計出的MUAP幅值序列，確定出輸入信號中所包含的最大MUAP的幅值，并以此作為備選峰值檢測點的初始閾值。

步驟2：備選峰值點檢測。所謂的備選峰值點檢測，就是在sEMG中檢測絕對峰值大于閾值的峰值點的位置，然后按照峰值由大到小的順序，將其對應的時刻點保存到1個一維數(shù)組中。

步驟3：MUAP識別與提取。從峰值索引數(shù)組中得到一個峰值索引之后，通常的做法是采用一個固定窗長的窗口，在索引點對稱截取出一個待分析信號段與模板波形進行匹配。然而，不同層次的MUAP的持續(xù)時間不同，同一層次不同形態(tài)MUAP關于峰值點并不都是對稱的，而且不同波形其峰值到左右兩邊比例也不同(如圖1所示)。因此，采取移動的可變窗長截取索引點附近的待分析信號段。除此之外，為了方便進行參數(shù)計算，將不同窗長截取的MUAP波形都進行了3次樣條插值，使其與模板波形具有相同的點數(shù)。

待分析信號段與模板波形的匹配程度，主要采用相關系數(shù)與殘差比進行判定，兩個參數(shù)分別由式(4)和式(5)進行計算。

(4)

(5)

式中，Tpi(t)表示第i個MUAP模板，xj(t)表示第j個信號片段，cov(·)表示協(xié)方差運算，D(·)表示方差運算，C表示相關系數(shù)矩陣，R表示殘差比矩陣。

得到相關系數(shù)矩陣C和殘差矩陣R之后，采用如下規(guī)則進行MUAP波形判定：

(1)如果最大相關系數(shù)大于0.95，則判定具有最大相關系數(shù)的模板與信號片段匹配；

(2)如果最大相關系數(shù)介于0.7與0.95之間，則判定相關系數(shù)大于0.7且殘差比最小的模板與信號片段匹配；

(3)如果最大相關系數(shù)小于0.7，則判定該峰值索引點附近無MUAP發(fā)放。

識別和提取完某一峰值索引點的MUAP之后，則將殘余sEMG用于下一峰值索引點的MUAP識別與提取。當遍歷所有備選峰值點之后，退出子程序，完成該閾值下的MUAP識別與提取。

從可分離的角度來看，當所發(fā)放MUAP的幅值小于基線時，就很難檢測和識別到相應的MUAP。因此，利用式(6)計算的最小閾值來判定是否結束分解。

(6)

式中，peakn為基線信號中的峰值點，N為滿足|peak1)|-|peakn|≤σ的最大n值。

2 sEMG檢測及分解

2.1sEMG采集

為了驗證算法的有效性，實驗募集了8名(男性4名，女性4名，年齡20～24歲)身體健康的大學生志愿者進行sEMG采集實驗，受試者均知情同意。實驗中，受試者按如圖3(a)所示姿勢分別使用食指、中指和環(huán)指進行20% MVC，40% MVC和60% MVC指力跟蹤實驗(詳見文獻[12])。采用圖3(b)所示的粘貼在指淺屈肌(FDS)肌腹25%～76%處的一塊6×2(行×列)Ag-AgCl柔性陣列電極進行sEMG檢測。檢測到的sEMG經(jīng)過信號放大電路之后，采用RM6280C多道生理記錄儀進行記錄、顯示和存儲。同步記錄的還有以相同采樣率(2 kHz)采集的4個指力傳感器的輸出電壓信號。

2.2預處理

作為一種弱生物電信號，在采集sEMG的過程中不可避免地受環(huán)境噪聲、設備固有噪聲、工頻干擾、運動偽跡等噪聲的影響。因此，對原始sEMG依次進行了有效段自動提取、帶通濾波、FastICA去工頻干擾和小波包去噪等預處理來提高信號的質(zhì)量，具體處理過程與效果參見文獻[15]。

2.3MUAP分解

根據(jù)去噪后的基線信號和MVC力量水平下采集的sEMG，估計得到可分離MUAP的最大幅值約等于4 mV，最小幅值約等于1 mV。為了計算方便，可分離MUAP的幅值分別設為4、2、1 mV。對于12通道的sEMG，采取的是單通道分別處理的方式。單通道sEMG分解可以概括為4個步驟。

步驟1：輸入待分解sEMG，根據(jù)信號最大峰值確定其中所包含的最大MUAP幅值，并以該幅值作為初始閾值進行備選峰值檢測；

步驟2：檢測信號中絕對峰值大于閾值的峰值點的位置，然后按照峰值由大到小的順序，將所有備選峰值點對應的索引點保存到1個一維數(shù)組中；

步驟3：從最大峰值索引點開始，在索引點附近根據(jù)4種不同幅值的MUAP模板，選擇不同變化范圍的窗長(4 mV：15～20 ms；2 mV：10～16 ms；1 mV：5～12 ms)移動截取待分析信號片段，計算特征參數(shù)，識別MUAP；若在某個備選峰值點附近識別出有效的MUAP，則從信號中減去相應的模板；若未識別出有效MUAP，則進入下一個索引點進行識別，直到遍歷所有備選峰值點；

步驟4：檢測殘余信號峰值，若殘余信號峰值大于最小檢測閾值，則以更小MUAP的幅值作為備選峰值點的檢測閾值，返回步驟2繼續(xù)對殘余信號進行分解；否則，結束分解。

3 結果

3.1單通道分解結果

圖4為一段60%力量水平下，食指恒力收縮時采集到的sEMG，經(jīng)過預處理和分解的效果圖。從圖4中可以看到，預處理之后的信號無明顯的工頻干擾，高頻噪聲明顯降低，曲線更加平滑；另外，漸進逐層分解后得到了幅值約為1、2、4 mV等3層MUAPTs，共12個MUAPTs加上一個幅值和能量較少的殘余信號，且殘余信號中無明顯MUAP波形殘存；從不同MUAP的發(fā)放時間序列上可以看到，所提取出的MUAP與原始波形在時間軸上具有很好的對應關系；仔細觀察右側的細節(jié)放大圖，可以清楚地看到，疊加波形中的主體MUAP都被成功地分解出來，并歸類到了相應的MUAPTs中，且其發(fā)放時刻定位準確。由此可見，所提出的“先大后小”的漸進式分解算法，可以有效地實現(xiàn)疊加波形的分解。

3.2MUAP發(fā)放總數(shù)隨力量水平的變化趨勢

分別統(tǒng)計了8名受試者(3組/人)食指、中指和環(huán)指在20% MVC、40% MVC和60% MVC等3個力量水平下采集到的時長為2 s的12通道sEMG分解得到的MUAP數(shù)目，結果如圖5所示。顯然，對于所有受試者，3種手指活動模式下，隨著力量水平的增加，MUAP發(fā)放的數(shù)目都呈增加趨勢。

3.3不同手指活動模式下不同大小MUAP發(fā)放比重變化

從圖5可以觀察到，MUAP發(fā)放總數(shù)總是隨著力量水平的增加而增加。然而，不同活動模式下，所分解得到的3種不同大小MUAP的變化并不都隨著力量水平的增加而增加。由于不同的受試者，不同活動模式下分解得到的MUAP的總數(shù)不同，因此以3種幅值不同的MUAP數(shù)目占總數(shù)目的百分比(即比重)為指標，統(tǒng)計了8名受試者在不同活動模式下(每種活動模式3組重復試驗數(shù)據(jù)，共24組數(shù)據(jù))不同MUAP的貢獻。結果如圖6所示。

從圖6所示的統(tǒng)計對比分析結果可以觀察到，3種手指活動模式下，大幅值MUAP(A≈4 mV)的發(fā)放比重，隨著力量水平的增加而增加；對于中等幅值(A≈2 mV)MUAP，隨著力量水平的增加，不同手指則表現(xiàn)出不同的變化模式：食指，顯著性增加；中指，20% ～ 40% MVC顯著性增加，40% ～ 60% MVC則基本保持穩(wěn)定；環(huán)指，20%～40% MVC顯著性增加，40% ～ 60% MVC顯著性減少。小幅值(A≈1 mV)MUAP的發(fā)放比重，隨著力量水平的增加而減小。

4 討論

近年來，隨著盲源分離、人工神經(jīng)網(wǎng)絡和人工智能等新興技術的不斷引入，sEMG分解研究取得了很大的進步，但MUAP疊加波形分解仍是sEMG分解中難點問題之一。針對這一問題，本研究設計了一種漸進式sEMG分解算法，基于先驗模板匹配和分級量化，實現(xiàn)了對不同力量水平下的MUAP提取。實測sEMG分解結果顯示，與其他報道給出的結果[1, 8, 16]一樣，隨著力量水平的增加，分解得到的MUAP總數(shù)呈增加趨勢。由此可見，本研究的分解算法提取出的主體MUAPTs，確實能夠從整體上反映肌肉活動的強弱變化。

目前，基于MUAP形態(tài)的sEMG分解算法，在處理疊加波形時主要有兩種思路：一種是將不同的MUAP模板進行不同程度的疊加合成再與sEMG進行比對；另一種則是從sEMG中逐個剝離MUAP[16]。本研究的基本思路與后者相同，但對分解策略進行了優(yōu)化，以適應高力量水平的sEMG分解。首先，本研究采用了一種簡潔的類似二進制的分級量化方法設置不同MUAP的幅值；其次，傳統(tǒng)的sEMG分解算法在波形分割時，通常采用對稱的固定窗口進行信號分割，忽略了不同MUAP時長不同，對稱性不同[17]為了彌補其不足，在同一個峰值點，采用不同窗長，移動截取的方式，截取多個待分析波形，以便更加準確地對主體MUAP進行識別；最后，在波形識別時，相關系數(shù)只作為基本判定指標，還考慮了殘差的大小，能夠很好的減少誤判。從圖4中60% MVC力量水平下的sEMG分解結果來看，所提出的算法與其他基于波形檢測和模式識別的方法[1, 8]一樣，所得到的MUAPTs與預處理后的sEMG中的MUAP疊加波形，在時間軸上具有很好的對應關系，同一時刻同時發(fā)放的MUAP也得到了有效分解。另外，分解之后得到殘余分量較小，無明顯MUAP殘余。綜上，本研究所提算法分解得到的主體MUAPTs中MUAP的發(fā)放時刻定位準確，有效地實現(xiàn)了疊加波形的分解，適用于高力量水平下的sEMG進行分解。

由圖6可知，同一手指活動模式下，隨著力量水平的增加，不同大小MUAP的比重不同，即不同大小MUAP對手指力量的貢獻率不同。小幅值(A≈1 mV)、持續(xù)時間短的MUAP，在較低力量水平下(20% MVC)的貢獻率最大，隨著力量水平的增加，貢獻率降低；大幅值(A≈4 mV)、持續(xù)時間長的MUAP，在較低力量水平下(20% MVC)的貢獻率最小，隨著力量水平的增加，貢獻率增加。產(chǎn)生上述結果的原因可能是，在較低力量水平下，激活閾值較小的MU被激活并發(fā)放幅值較小的MUAP，隨著力量水平增加，為了產(chǎn)生更大的肌力，具有更高激活閾值的MU被激活，并且主要由后激活的MU來產(chǎn)生肌力。事實上，已有結果證實，MUAP的大小主要是由MU中包含的肌纖維數(shù)決定，為了產(chǎn)生更大的肌力，包含有更多肌纖維數(shù)的高閾值MU將被激活[16]。

同時，本實驗結果還表明，不同手指活動模式下，部分MUAP貢獻率隨著力量水平的變化具有顯著性差異。如圖6所示，3種手指活動模式下，中等幅值MUAP(A≈2 mV)貢獻率隨著力量的變化趨勢具有顯著性差別。產(chǎn)生這一結果的原因，可能是不同手指活動時，神經(jīng)-肌肉系統(tǒng)采取了不同的MU募集策略，選擇性地激活MU[18]。

5 結論

針對高收縮力情況下sEMG分解中的疊加波形分解問題，本研究基于MU募集的“大小原則”和MUAP形態(tài)，提出一種“從大到小”的漸進式MUAP分解算法。為了驗證算法的有效性，對不同手指活動模式下的指淺屈肌sEMG進行了分解。結果顯示：該算法能有效分解sEMG；隨著力量水平的增加，分解得到的MUAP數(shù)目呈現(xiàn)增加趨勢；不同大小MUAP的比重隨著活動手指和力量水平的變化而變化。所提出的分解算法，為深入研究神經(jīng)-肌肉控制系統(tǒng)中的MU募集與發(fā)放信息提供了可能。而所提出算法MUAP基本模板類型較少，只提取了一部分MU所發(fā)放的MUAPTs信息，后續(xù)研究中需要建立有效的模板更新機制，實現(xiàn)更多形態(tài)MUAP的提取。此外，后續(xù)研究需增加樣本量，并設計實驗對自主收縮時MU的募集模式進行研究。

[1] De Luca CJ, Adam A, Wotiz R,etal. Decomposition of surface EMG signals[J]. Journal of Neurophysiology,2006,96(3):1646-1657.

[2] 侯文生，萬莎，吳小鷹，等. 手指運動姿態(tài)檢測及假肢手指動作實時控制[J]. 重慶大學學報,2013(04):103-109.

[3] 姚博. 表面肌電信號分解算法及其在小兒腦癱評詁中的應用[D]. 合肥：中國科學技術大學,2012.

[4] 王耀兵，季林紅，黃靖遠. 康復訓練過程中上肢運動肌電信號特征的對比分析[J]. 中國臨床康復,2004(23):4721-4723.

[5] 王剛，王志中，胡曉，等. 基于最佳小波包的表面肌電信號分類方法[J]. 中國醫(yī)學物理學雜志,2006(01):45-48.

[6] Garcia GA, Okuno R, Akazawa K. A decomposition algorithm for surface electrode-array electromyogram-A noninvasive, three-step approach to analyze surface EMG signals[J]. IEEE Engineering in Medicine and Biology Magazine,2005,24(4):63-72.

[7] Naik GR, Arjunan S, Kumar D. Applications of ICA and fractal dimension in sEMG signal processing for subtle movement analysis: a review [J]. Australasian Physical & Engineering Sciences in Medicine, 2011,34(2):179-193.

[8] Nawab SH, Chang SS, De Luca CJ. High-yield decomposition of surface EMG signals[J]. Clinical Neurophysiology, 2010,121(10):1602-1615.

[9] Farina D, Enoka RM. Surface EMG decomposition requires an appropriate validation [J]. Journal of Neurophysiology, 2011,105(2):981-982.

[10] 安媛. 手指活動影響前臂多腱肌運動單元募集模式的初步研究[D]. 重慶：重慶大學,2012.

[11] 李強. 表面肌電信號的運動單位動作電位檢測[D]. 合肥：中國科學技術大學,2008.

[12] 侯文生，安媛，楊丹丹，等. 食指力量對指淺屈肌運動單元募集模式的影響[J]. 醫(yī)用生物力學,2012(06):648-654.

[13] Merletti R, Holobar A, Farina D. Analysis of motor units with high-density surface electromyography[J]. Journal of Electro- Myography and Kinesiology,2008,18(6):879-890.

[14] 姚博，楊基海，陳香，等. 基于稀疏分量分析的欠定盲源分離用于表面肌電信號分解[J]. 航天醫(yī)學與醫(yī)學工程,2012(02).

[15] 羅萬國，吳小鷹，楊丹丹，等. 基于形態(tài)學的表面肌電信號分解算法研究[Z]. 成都:電子科技大學生命科學與技術學院,2013.

[16] Stashuk D. EMG signal decomposition: how can it be accomplished and used? [J]. Journal of Electromyography and Kinesiology, 2001,11(3SI):151-173.

[17] Chauvet E, Fokapu O, Hogrel JY,etal. Automatic identification of motor unit action potential trains from electromyographic signals using fuzzy techniques[J]. Medical & Biological Engineering & Computing, 2003,41(6):646-653.

[18] Butler TJ, Kilbreath SL, Gorman RB,etal. Selective recruitment of single motor units in human flexor digitorum superficialis muscle during flexion of individual fingers[J]. Journal of Physiology-London, 2005,567(1):301-309.

StudyonScaledDecompositionofSurfaceEMGBasedonPriorTemplates

LUO Wan-Guo1HOU Wen-Sheng1,2#ZHENG Xiao-Lin1,2#WAN Xiao-Ping1CHEN Hai-Yan2ZHOU Ping2WU Xiao-Ying1*

1(DepartmentofBiomedicalEngineering,ChongqingUniversity,Chongqing400044,China)2(ChongqingEngineering&TechnologyResearchCenterforMedicalElectronics,Chongqing400044,China)

surface electromyography; decomposition; prior templates; motor unit action potentials

10.3969/j.issn.0258-8021. 2014. 03.015

2014-02-05， 錄用日期:2014-04-20

國家科技支撐計劃(2012BAI16B02)；國家自然科學基金(30970758)；重慶市自然科學基金(cstc2012jjA10103)

R318

D

0258-8021(2014) 03-0366-07

#中國生物醫(yī)學工程學會會員(Member, Chinese Society of Biomedical Engineering)

*通信作者(Corresponding author)，E-mail: wuxiaoying69@163.com