基于肌電信號的肘關(guān)節(jié)肌肉疲勞算法性能研究

孟慶豐，陶慶，2，來全寶，胡玉舸

（1.新疆大學(xué)機械工程學(xué)院，新疆 烏魯木齊 830047；2.新疆大學(xué)機械工程博士后研究中心，新疆 烏魯木齊 830047）

1 引言

肌肉疲勞是人們?nèi)粘Ｉ钪幸环N常見的生理現(xiàn)象，一般是指由長時間運動或發(fā)力引起的肌無力或酸痛的感覺[1]。日常生活中我們或多或少都會遇到肌肉疲勞的現(xiàn)象，但是如果不顧肌肉的疲勞而過度運動或長時間活動的話，很可能會造成肌體的永久性損傷，更甚者會帶來生命危險[2]，因此對于肌肉疲勞評價的研究足以引起我們的高度重視。對于肌肉疲勞評價方法的研究與分析，無論是在體育訓(xùn)練還是在康復(fù)醫(yī)學(xué)、運動醫(yī)學(xué)等領(lǐng)域都具有重要意義。

如今，已知通過記錄表面肌電信號計算疲勞指數(shù)的方法可以分為兩個主要類別[3]。第一類由經(jīng)典的基于信號可量化頻譜變化的算法組成。例如經(jīng)典的平均頻率以及頻譜矩比[4]和基于小波的方法[5]都屬于此類。第二類考慮信號的非線性成分。例如考慮信號重復(fù)模式或分形測量的遞歸量化分析和基于熵的方法的模糊近似熵都屬于此類[6]。

由于上肢的運動功能遠(yuǎn)比下肢復(fù)雜，上肢在人體日常的工作生活扮演著重要的角色，而肘關(guān)節(jié)的屈伸動作甚至出現(xiàn)在人體運動的每個瞬間，并且上肢中肘關(guān)節(jié)的運動減少50%將導(dǎo)致上肢運動功能降低近80%[7]，肘關(guān)節(jié)的正常運動對人體的工作生活是非常重要的，因此探究更有效的肘關(guān)節(jié)肌肉疲勞評價方法迫在眉睫。

為此，國內(nèi)外學(xué)者已經(jīng)針對肘關(guān)節(jié)肌肉疲勞評價方法展開了大量的研究，不同方法得出的評價指標(biāo)獲得的肌肉疲勞評價效果不同。文獻(xiàn)[8]利用AEMG與MPF等線性和Lempel?Ziv 復(fù)雜度與DET等非線性指標(biāo)分析不同負(fù)載下肱二頭肌疲勞過程中表面肌電信號的變化，結(jié)果顯示疲勞指標(biāo)的變化具有非疲勞特異性。

文獻(xiàn)[9]采集不同狀態(tài)的上肢兩側(cè)肌電信號，將時域的積分肌電值、均方根振幅和頻域的平均功率頻率、中位頻率作為疲勞評價指標(biāo)，研究得出上肢一側(cè)近端肌肉疲勞會使另一側(cè)的運動信號產(chǎn)生伴隨性變化，即產(chǎn)生中樞共驅(qū)動現(xiàn)象。

文獻(xiàn)[10]分析時域的積分肌電值、均方根振幅和頻域的平均功率頻率、中位頻率等指標(biāo)隨著疲勞程度加深的變化趨勢，提出一種肌肉疲勞出現(xiàn)時刻的估計方法。文獻(xiàn)[11]將時域的iEMG、RMS與頻域的MPF、Fc 分析評價指標(biāo)，利用SPSS 統(tǒng)計分析時域與頻域之間的參量關(guān)系以及上肢肌肉疲勞間的關(guān)系，結(jié)果與實際情況更貼合。

以上研究主要是對上肢疲勞時域或頻域指標(biāo)變化規(guī)律的探索，缺乏對各疲勞指標(biāo)的評價效果進(jìn)行分析。為此，研究了肘關(guān)節(jié)肌肉疲勞評價方法的評價效果，記錄測試對象在不同負(fù)載下肘關(guān)節(jié)做屈曲運動的sEMG 信號，計算MNF、SMR、WIRMI1551、fApEn與RQA%DET的指標(biāo)，并進(jìn)一步比較5種方法對肘關(guān)節(jié)肌肉疲勞的評價效果。

從而獲得更有效的肘關(guān)節(jié)肌肉疲勞評價方法。將從肌肉疲勞評價方法的抗干擾性與區(qū)分疲勞程度的能力兩個方面對5種疲勞評價方法的效果進(jìn)行比較，分析所得結(jié)果，為肘關(guān)節(jié)肌肉疲勞提供一個更加準(zhǔn)確的人因工程定量評價方法。

2 研究方法

2.1 肌肉疲勞評價方法

2.1.1 平均頻率

研究表明，隨著肌肉疲勞的累積，sEMG 信號功率譜的主頻率由高轉(zhuǎn)向低，平均功率值逐漸變小，并且疲勞程度越大，下降趨勢越強。基于對功率譜密度（PSD）的估計，計算了譜疲勞評價方法指標(biāo)。PSD 是利用Welch 方法估算的（其中k=15），并使用Hamming窗函數(shù)進(jìn)行50%的重疊[12]?；赑SD估計的MNF：

式中：f—頻率；fs—采樣率。

2.1.2 頻譜矩

作為平均頻率的替代，頻譜矩（SMR）方法被Dimitrov等人描述為新的光譜指數(shù)，它是對H/L比率方法的改進(jìn)。文獻(xiàn)[13]用平滑函數(shù)加權(quán)PSD估計，他們建議使用頻率的倒數(shù)作為低頻帶的加權(quán)函數(shù)，高頻帶的加權(quán)函數(shù)是頻率平方，立方，四次或五次方。基于階數(shù)k的譜矩Mk：

頻譜矩為：

式中：f1=5Hz和f2=500Hz—積分的上下限，階數(shù)k=5。

2.1.3 基于離散小波的方法

量化頻率壓縮的其他增強疲勞評價方法是Gonzalez?Izal等人報道的基于小波的方法[14]。它們是基于固定小波變換（SWT）的頻譜估計而不是PSD 估計。WIRM1551 比值是從Gonzalez?Izal等人建議的小波方法的變體中選擇的。它被定義為：

式中：D1—使用sym5小波的SWT的第一尺度的功率譜；D5—SWT的第五尺度的功率譜。

正如建議的那樣，f1=8Hz和f2=500Hz為積分的上下限。

2.1.4 模糊近似熵

熵估計是用于疲勞評價的另一種非線性方法。樣本熵（Sam?pEn）是一種計算時間序列熵的方法，為條件概率的負(fù)自然對數(shù)，該條件概率是從一個紀(jì)元開始的兩個序列在m維向量相似，而在m+1 維向量仍然相似，當(dāng)向量的長度m=2，閾值r=0.3 為最佳搭配。為了疲勞檢測的目的，文獻(xiàn)[15]報告了一種改進(jìn)的算法，稱為模糊近似熵（fApEn）。與樣本熵相反，如果兩個序列相差任意恒定偏移量d，則模糊近似熵會將它們視為相似。以m=2，r=0.6為最佳組合計算模糊函數(shù)來獲得模糊近似熵的近似度。模糊函數(shù)為：

2.1.5 遞歸量化分析

遞歸量化分析（RQA）是一種試圖量化分析信號中重復(fù)模式的非線性疲勞評價方法[6]。首先構(gòu)造距離矩陣，從時間信號的角度導(dǎo)出長度為D（嵌入維）的向量，其中時間延遲間隔λ=5，嵌入尺寸D=15。然后矩陣通過其平均值進(jìn)行歸一化處理，研究表明，閾值r=0.75時矩陣會產(chǎn)生僅包含零和一的二進(jìn)制遞歸映射。%DET是用于評估復(fù)發(fā)圖的，其值為形成一定鏈長的對角線與矩陣中所有對角線的數(shù)量之比。研究發(fā)現(xiàn)，在疲勞任務(wù)下，RQA%DET值的上升與運動單元的更長時標(biāo)同步非常吻合。

2.2 疲勞評價方法的定量化分析

由于疲勞評價方法難以基于測試對象的主觀感受來校準(zhǔn)，因此本研究使用定量化分析方法比較疲勞評價方法的差異性。定量化分析方法主要包括同一疲勞條件下評價方法的抗干擾性及不同疲勞條件下區(qū)分疲勞程度的能力。相關(guān)研究表明，肌肉疲勞評估指數(shù)呈現(xiàn)具有恒定斜率的線[16]，本研究利用最小二乘法對五種疲勞評估方法進(jìn)行了線性回歸分析，并求解出每個方法的確定系數(shù)R2與回歸方程的斜率k。

2.2.1 抗干擾性分析

疲勞評價方法預(yù)期在疲勞情況下產(chǎn)生平滑變化的疲勞指數(shù)。偏離這個平滑過程被認(rèn)為是干擾。疲勞指數(shù)的干擾可以通過確定系數(shù)R2來量化[17]。R2表示具有恒定斜率的線反映信號數(shù)據(jù)的相關(guān)性，R2的數(shù)值范圍從0到1。R2的值為0意味著信號中沒有線性趨勢。R2的值為1則意味著信號恰好是具有非零斜率且沒有額外干擾的直線。因此R2的值越大表明疲勞評估算法的抗干擾性越強。

2.2.2 區(qū)分疲勞程度能力分析

人體肘關(guān)節(jié)肌肉在不同負(fù)載下會產(chǎn)生不同的肌肉疲勞程度。柯爾莫哥洛夫?斯米洛夫檢驗（Kolmogorov?Smirnov test，即K-S檢驗）是一種基于概率統(tǒng)計理論的方法，它能夠計算兩個數(shù)據(jù)經(jīng)驗分布函數(shù)間的最大垂直距離Lmax，Lmax越大說明兩個數(shù)據(jù)相似性越小。

本研究利用K-S檢驗分析五種評價方法區(qū)分疲勞程度的能力[18]。即對五種評價方法在不同負(fù)載下的線性回歸方程斜率k進(jìn)行K-S檢驗得到的Lmax值，Lmax越大表明疲勞評估方法對不同負(fù)載下區(qū)分疲勞程度的能力越強。

3 肌電信號采集與處理

3.1 測試對象與試驗準(zhǔn)備

3.1.1 測試對象

為了采集數(shù)據(jù)分析比較不同疲勞評價方法，本研究采集了12名測試對象（6名女性和6名男性）肱二頭肌的表面肌電信號。測試對象的年齡在（19～26）歲之間，無皮膚過敏史且肱二頭肌皮膚表面無創(chuàng)傷和疤痕，上肢運動能力處于正常狀態(tài)。在試驗開始之前，所有測試對象都書面同意志愿參加本次試驗。

3.1.2 試驗準(zhǔn)備

為了施加恒定的力，要求每個測試對象以同一固定姿勢托舉啞鈴。實驗過程中要求測試對象保持直立坐姿，肘部處于放松狀態(tài)后抓舉啞鈴進(jìn)行實驗。有關(guān)動作設(shè)置，如圖1所示。

圖1 實驗設(shè)定Fig.1 Setup of the Experiment

數(shù)據(jù)采集過程中，首先引導(dǎo)測試對象以最大可實現(xiàn)的用力在短時間內(nèi)抓舉啞鈴以單獨測取每個測試對象的最大自主收縮力（MVC）。然后在隨后的三次試運行期間，根據(jù)每個測試對象的MVC調(diào)整他們要托舉的重量。在三次試驗中，要求測試對象分別托舉他們MVC 的20%，40%和60%的重量。在20%和40%MVC 負(fù)荷的情況下，所有測試對象持續(xù)3min 的指定時間，在60%MVC負(fù)荷水平下持續(xù)135s達(dá)到疲勞。每次試驗后，測試對象休息30min。

3.2 試驗設(shè)備與數(shù)據(jù)采集

3.2.1 試驗設(shè)備

葡萄牙Plux公司生產(chǎn)ErgoPlux的8通道無線表面肌電測量系統(tǒng)、上海勵圖醫(yī)療器材有限公司生產(chǎn)的LT?7型理療電極片。

3.2.2 數(shù)據(jù)采集

測試肌肉位置，如圖2所示。先用75%的醫(yī)用酒精擦拭測試對象右臂肱二頭肌的皮膚表面，清理皮膚表面污漬與角質(zhì)，再使用干燥的棉球擦凈皮膚表面殘留的酒精后將八組電極片以高度密集的形式放置在肱二頭肌上方。即電極以彼此相鄰的兩列連接，且彼此之間沒有任何間隔。根據(jù)其直徑，相鄰放置的電極的中心間隔為26mm。測試電極沿著所測試肌肉的肌纖維方向，貼在肱二頭肌最隆處，接地電極貼在肘關(guān)節(jié)處。試驗過程中，每名測試對象依次托舉他們最大負(fù)載（MVC）的20%、40%與60%重量，利用8 通道的表面肌電測量系統(tǒng)同步記錄測試對象的sEMG信號。

圖2 高密度電極片的貼點位置圖Fig.2 Location of the High?Density Electrode Sheet in the Biceps

3.3 數(shù)據(jù)處理

經(jīng)多次試驗，測試對象在測試結(jié)束后均在不同負(fù)載下感受到不同程度的肌肉酸痛僵硬，即肌肉處于疲勞狀態(tài)。為保證所分析的sEMG信號數(shù)據(jù)均在疲勞狀態(tài)下，研究采用測試過程中最后的1min進(jìn)行數(shù)據(jù)分析。首先同步采集8個通道的sEMG 會產(chǎn)生基線偏移現(xiàn)象，利用sgolay濾波器消除趨勢化，然后sEMG信號采集過程中會受到50Hz工頻及其他信號干擾，利用50Hz陷波與巴特沃斯帶通濾波器對sEMG信號去噪[18]，最后將處理后的8通道信號取平均值。處理后的示例信號，如圖3所示。

圖3 處理后的示例信號Fig.3 Exemplary Signals After Processing

根據(jù)五種疲勞評價方法的定義，使用Matlab R2014b分別編寫五種疲勞評價方法的算法，計算出每個測試對象在不同負(fù)載下各疲勞評價方法的指標(biāo)值，然后對各疲勞評價方法進(jìn)行定量化分析。定量化分析分為兩個部分，一是利用最小二乘法對各方法的指標(biāo)進(jìn)行線性回歸分析（圖4 展示了1 名測試對象在不同負(fù)載下最后一分鐘的MNF最小二乘法線性回歸分析結(jié)果，其他4 種方法與MNF 做相同的處理），比較線性回歸方程的確定系數(shù)R2的均值，利用SPSS25對R2進(jìn)行單因素方差分析，進(jìn)一步評價疲勞評價方法的抗干擾性；二是對回歸方程的斜率k進(jìn)行K-S檢驗得到不同負(fù)載下的Lmax值，并對Lmax進(jìn)行單因素方差分析，評價5種疲勞評價方法區(qū)分疲勞程度的能力。圖4中=0.792，正如預(yù)期的那樣，隨著負(fù)載水平的增加，R2值會增加，因為負(fù)載水平越高，MNF的線性趨勢越強。

圖4 MNF的線性回歸分析結(jié)果Fig.4 Linear Regression Analysis Results of MNF

4 結(jié)果

4.1 抗干擾評價結(jié)果

為量化不同疲勞算法的抗干擾性，本研究計算了12個測試對象在不同負(fù)載下各疲勞評價方法線性回歸方程的確定系數(shù)R2平均值與標(biāo)準(zhǔn)差，并對R2進(jìn)行單因素方差分析結(jié)果，如圖5所示。

圖5（a）顯示在20%MVC情況下，SMR的R2平均值為5種疲勞評價方法的最大值，最大平均值=0.237，且與MNF（p=8.28e?04）和RQA%DET（p=5.45e?05）方法存在顯著性差異。表明SMR在20%MVC下具有最好的抗干擾性。WIRM1551具有較大的R2平均值=0.223，說明WIRM1551在20%MVC下具有較好的抗干擾能力。RQA%DET的R2平均值最小，表明RQA%DET的抗干擾能力最差。

圖5（b）顯示在40% MVC 情況下，SMR的R2平均值仍然最大=0.328，且與MNF（p=1.38e?05）和RQA%DET（p=3.69e?06）在評價肘關(guān)節(jié)肌肉疲勞的抗干擾性的差異具有統(tǒng)計學(xué)意義。表明SMR 在40%MVC 下具有最強的抗干擾能力。MNF、WIRM1551與fApEn的R2平均值相差很小，分別為=0.289，=0.288，說明這3種方法在40%MVC下具有相同的抗干擾能力。而RQA%DET的抗干擾性最差。

圖5（c）顯示MNF、WIRM1551與fApEn在60%MVC具有較大的R2平均值，但SMR仍然為最大值=0.528，且與RQA%DET（p=3.08e?07）的差異具有統(tǒng)計學(xué)意義。說明MNF、WIRM1551與fApEn在60%MVC下具有較好的抗干擾性，但SMR的抗干擾能力是最強的。RQA%DET依舊是5種方法中抗干擾性最差的。

圖5 五種肌肉疲勞評價方法在不同負(fù)載下的R2Fig.5 Five Muscle Fatigue Detection Methods for R2 Under Different Loads

綜上所述，SMR在不同負(fù)載下的R2平均值均為最大值，擁有最好的抗干擾性。WIRM1551 在不同負(fù)載下的R2平均值僅比SMR 小一些，同樣具有較好的抗干擾性。MNF 與fApEn 僅在大負(fù)載下?lián)碛休^大的R2平均值，但是在小負(fù)載，低疲勞情況下的R2平均值相比SMR還是有較大差距的，因此，MNF與fApEn僅在大負(fù)載下?lián)碛休^強的抗干擾能力。而RQA%DET在不同負(fù)載下的抗干擾能力是最弱的。

4.2 疲勞程度區(qū)分能力評價結(jié)果

為量化分析5種疲勞評價方法區(qū)分疲勞程度的能力，對線性回歸方程的斜率k進(jìn)行K-S檢驗，得出不同負(fù)載下肘關(guān)節(jié)肌肉疲勞指標(biāo)經(jīng)驗概率分布函數(shù)的最大垂直距離Lmax，并對Lmax進(jìn)行單因素方差分析。全局結(jié)果，如圖6所示。局部Lmax量化分析結(jié)果，如圖7所示。

圖6 全局Lmax量化分析Fig.6 Gl obalLmaxQuantitativeAnalysis

不同負(fù)載下全局Lmax量化分析結(jié)果，如圖6所示。結(jié)果顯示在疲勞程度較大（40%vs60%）的情況下，MNF、WIRM1551與fAp?En的Lmax平均值較大，SMR的Lmax平均值為最大值；在疲勞差異大（20%vs60%）的情況下，SMR與WIRM1551的Lmax平均值為最大值；在疲勞程度?。?0%vs40%）的情況下，SMR的Lmax平均值呈現(xiàn)明顯的尖點。而RQA%DET的Lmax平均值在三種情況均為最小值。全局Lmax量化分析結(jié)果表明在5種方法中SMR擁有最強的區(qū)分肘關(guān)節(jié)肌肉疲勞程度的能力，而RQA%DET是最弱的。

20%與60%MVC的Lmax量化分析結(jié)果，如圖7（a）所示。結(jié)果顯示SMR 與WIRM1551的Lmax為最大值，均為0.883，表明SMR與WIRM1551 在區(qū)分20%與60%MVC 的疲勞程度的能力最出色。而RQA%DET 是5 種方法中評價疲勞程度效果最差的。SMR與RQA%DET的差異具有統(tǒng)計學(xué)意義（p=1.09e?02）。表明SMR與WIRM1551在疲勞程度差異較大時有最強的區(qū)分疲勞程度能力。40%與60%MVC的Lmax量化分析結(jié)果，如圖7（b）所示。結(jié)果表明MNF、SMR、WIRM1551 與fApEn 均呈現(xiàn)較大的Lmax平均值，說明4 種方法在評價大負(fù)載引起的疲勞程度均有很出色的能力。

其中SMR 是最出色的，Lmax=0.766，且與RQA%DET（p=8.93e?03）存在顯著性差異。而RQA%DET區(qū)分疲勞程度的能力最差。說明在疲勞程度較大時，MNF、SMR、WIRM1551與fApEn均具有較好的區(qū)分能力，且SMR 是最好的。20%vs40%MVC的Lmax量化分析結(jié)果，如圖7（c）所示。

結(jié)果表明SMR的Lmax平均值最大，最大值Lmax=0.622，且與RQA%DET（p=1.6e?03）存在顯著性差異。RQA%DET的Lmax平均值仍是最小的。說明在疲勞程度較小的情況下，SMR仍然擁有最強的區(qū)分肌肉疲勞程度的能力。

綜上所述，SMR相較于另外4種方法在不同負(fù)載下都具有最大的Lmax平均值，說明SMR在區(qū)分不同負(fù)載引起的肘關(guān)節(jié)肌肉疲勞程度的能力最出色。WIRM1551 雖沒有SMR 那么出色，但在不同負(fù)載下均具有較好的區(qū)分疲勞程度的能力。MNF與fApEn僅在疲勞程度大或疲勞差異較大時具有較好的區(qū)分能力。而RQA%DET 在不同負(fù)載下具有最小的Lmax平均值，說明RQA%DET在不同負(fù)載下引起的肌肉疲勞程度區(qū)分性最差。

5 結(jié)論

通過采集肱二頭肌的表面肌電信號計算了5種常見肌肉疲勞評價算法，從評價算法的抗干擾性與區(qū)分疲勞程度的能力兩個方面量化其疲勞評價性能，以R2大小體現(xiàn)抗干擾能力的強弱，以Lmax展現(xiàn)區(qū)分疲勞程度的能力。實驗結(jié)果顯示，頻譜距SMR同時擁有最大的R2和Lmax值，即頻譜距SMR在評價肘關(guān)節(jié)肌肉疲勞時具有最好抗干擾和區(qū)分疲勞程度的能力。故頻譜距SMR能夠更加及時有效的檢測出肘關(guān)節(jié)肌肉疲勞狀態(tài)。因此在今后研究肘關(guān)節(jié)屈曲運動的肌肉疲勞時，建議使用頻譜矩（SMR）進(jìn)行實時監(jiān)測。