楊 錚 王立玲 馬 東
河北大學(xué)電子信息工程學(xué)院,河北省數(shù)字醫(yī)療工程重點(diǎn)實(shí)驗(yàn)室,河北 保定 071002)
非平穩(wěn)信號是一種分布參數(shù)隨時(shí)間變化的隨機(jī)信號,如生物醫(yī)學(xué)工程中的腦電信號與表面肌電信號(surface electromyography,sEMG)。在醫(yī)療領(lǐng)域,表面肌電信號的有效處理不僅為患者提供肌肉疲勞診斷依據(jù),而且提供了有效的治療手段[1-2]。因傳統(tǒng)特征值表征肌肉疲勞時(shí)變化不明顯、反映不夠靈敏,很有可能在肌肉已經(jīng)發(fā)生疲勞時(shí)而未檢測出疲勞,造成肌肉損傷。因此,如何在肌肉受到損傷之前檢測出疲勞,防止由普通性肌肉疲勞進(jìn)一步惡化發(fā)展成肌肉損傷已成為當(dāng)今一個(gè)值得研究問題。時(shí)變系數(shù)建模方法是處理非平穩(wěn)信號典型方法。已廣泛應(yīng)用于多種非平穩(wěn)信號處理領(lǐng)域。
時(shí)變系統(tǒng)的參數(shù)辨識(shí)方法主要有自適應(yīng)算法和基函數(shù)展開法。自適應(yīng)算法中常見有經(jīng)典遞推最小二乘法、最小均方差算法和卡爾曼濾波法等。自適應(yīng)算法是針對信號具有弱平穩(wěn)特性時(shí),時(shí)變系統(tǒng)參數(shù)變化較慢,其可以對時(shí)變系統(tǒng)參數(shù)進(jìn)行準(zhǔn)確辨識(shí)。否則,辨識(shí)參數(shù)會(huì)因其收斂性的缺陷,導(dǎo)致時(shí)變系統(tǒng)參數(shù)的結(jié)果估計(jì)產(chǎn)生延遲[3]。而基函數(shù)展開法是針對信號具有較強(qiáng)非平穩(wěn)特性時(shí),其可以對時(shí)變系數(shù)進(jìn)行有效估計(jì)[4]。時(shí)變系統(tǒng)參數(shù)表示為一組已知基函數(shù)的線性加權(quán)組合[5],將時(shí)變系統(tǒng)建模問題變成了關(guān)于基函數(shù)的時(shí)不變參數(shù)辨識(shí)問題,通過對時(shí)不變參數(shù)的辨識(shí)進(jìn)而得到時(shí)變系數(shù)。而基函數(shù)選擇常見有傅里葉基、勒讓德多項(xiàng)式和小波基,不同基函數(shù)具有不同的逼近特性,其中傅里葉基和勒讓德基函數(shù)適用于有效辨識(shí)變化緩慢的時(shí)變系數(shù),而小波基函數(shù)適用于有效辨識(shí)變化劇烈的時(shí)變系數(shù)。國外 Chon 等及國內(nèi)劉洪濤等根據(jù)勒讓德多項(xiàng)式具有靈活的逼近特性對時(shí)變系數(shù)進(jìn)行展開并取得了較好的逼近效果[6-7]。徑向基函數(shù)(radial basis function,RBF)具有中心對稱特性、良好的局部特性、多尺度特性、多分辨率特性,已廣泛應(yīng)用于系統(tǒng)辨識(shí)領(lǐng)域[8]。
研究人員通過大量的實(shí)驗(yàn)揭示,伴隨著肌肉的疲勞,功率譜一般由高頻向低頻漂移,中位頻率(median frequency,MF)和平均功率頻率(mean power frequency,MPF)的值都呈下降趨勢,與肌肉的疲勞成相關(guān)性[9]。在實(shí)際應(yīng)用中,MPF指標(biāo)對肌肉活動(dòng)狀態(tài)和功能狀況的敏感性強(qiáng)于MF指標(biāo)[10]。但相關(guān)研究表明,頻域中影響表面肌電信號頻譜變化的因素主要與神經(jīng)傳導(dǎo)速率、動(dòng)作電位持續(xù)時(shí)間有關(guān)。在動(dòng)態(tài)收縮至疲勞過程中肌肉溫度會(huì)升高,溫度升高會(huì)影響表面肌電信號頻譜高頻成分增加,導(dǎo)致MPF、MF變化不明顯[11-12]。這些現(xiàn)象是由神經(jīng)傳導(dǎo)速率變化導(dǎo)致的[13]。但上述傳統(tǒng)特征值是基于短時(shí)傅里葉分析的,在檢測肌肉疲勞時(shí),采樣時(shí)間較長、實(shí)時(shí)性和分辨率較差,是傳統(tǒng)方法在評估疲勞領(lǐng)域中制約其應(yīng)用的一個(gè)關(guān)鍵問題。因此,如何有效進(jìn)行實(shí)時(shí)快速檢測及評估肌肉疲勞,已成為當(dāng)今一個(gè)值得研究問題。
受試者為10例健康男性,測試者均無上肢肌肉勞損病史,身體質(zhì)量指數(shù)正常,實(shí)驗(yàn)前24 h內(nèi)均無參加任何劇烈活動(dòng)。每位受試者在實(shí)驗(yàn)前都接受實(shí)驗(yàn)通知,并簽署知情同意書,進(jìn)行試驗(yàn)動(dòng)作培訓(xùn);實(shí)驗(yàn)地點(diǎn)在河北省數(shù)字醫(yī)療工程重點(diǎn)實(shí)驗(yàn)室。
實(shí)驗(yàn)前,實(shí)驗(yàn)室調(diào)整溫度27℃,將受試者首飾物品取下,電極片貼放于右上肢肱二頭肌中部。每組受試者右手心向上,上臂與身體平行,小臂與上臂肘角度為90°,手持1.5 kg啞鈴開始做等長收縮運(yùn)動(dòng)直至受試者主觀感覺肌肉疲勞不能繼續(xù)為止,采集并記錄右上肢肱二頭肌疲勞前期、后期數(shù)據(jù)。由于每位受試者年齡、個(gè)體之間差異較大,所以要求啞鈴重量、電極貼放位置等必須一致。
實(shí)驗(yàn)設(shè)備與材料:實(shí)驗(yàn)開始前,準(zhǔn)備實(shí)驗(yàn)設(shè)備與材料,包括計(jì)算機(jī)、無線表面肌電采集系統(tǒng)(包含DTS肌電傳感器1個(gè)、DTS桌面接受盒、雙電極夾1個(gè)、5 V DC充電線、Mini USB數(shù)據(jù)線、彈性固定帶1條)、電極片2片、酒精、棉簽、受試者信息。
采用品牌為美國NORAXON、型號為表面肌電信號采集系統(tǒng)、版本為MR3.6軟件。利用無線表面肌電信號傳感器,使用一次性Ag/Agcl電極片,設(shè)定時(shí)間常數(shù)為0.05 s、采樣頻率為1 500 Hz,根據(jù)肌肉模型粘貼一次性電極片,電極片間距為20 mm。
1.3.1濾波
由于表面肌電信號比較微弱,有用信號分布在0~500 Hz頻率范圍之間,主要能量部分分布在50~150 Hz頻率范圍之間。由體表電極檢測到的表面肌電信號主要含有工頻干擾(50 Hz)、基線漂移和心電干擾(5~20 Hz)等,這些噪聲會(huì)嚴(yán)重影響表面肌電信號質(zhì)量。為了增強(qiáng)表面肌電信號中的有效成分,抑制噪聲和偽跡,噪聲的有效消除對肌電信號的后續(xù)處理至關(guān)重要。采用巴特沃斯20~500 Hz帶通濾波器消除基線漂移和心電干擾,巴特沃斯帶阻濾波器消除50 Hz工頻干擾[14]。
1.3.2歸一化
不同受試者肌肉活性存在差異,為了能夠統(tǒng)一比較不同參數(shù)特征值對疲勞表征效果,需要對表面肌電信號進(jìn)行歸一化操作。
本試驗(yàn)中自回歸模型(autoregressive,AR)描述是一種“短時(shí)平穩(wěn)”隨機(jī)過程,該隨機(jī)過程在該時(shí)刻的觀察值xn與該時(shí)刻之前的觀測值xn-i存在相關(guān)性。設(shè)p階參數(shù)自回歸模型計(jì)算公式為
(1)
式中:ai(n)是自回歸系數(shù),也是AR(p)模型參數(shù);n=1,2,…,N,N為采樣數(shù)據(jù)長度;en是平穩(wěn)白噪聲過程;xn是觀察值,xn-i是觀測序列值。
由于采集肌電設(shè)備采樣頻率為1 500 Hz,因此時(shí)變系統(tǒng)參數(shù)變化太快,自適應(yīng)算法的收斂性存在缺陷,因此研究采用基函數(shù)展開式方法。基函數(shù)展開式方法對時(shí)變系統(tǒng)進(jìn)行辨識(shí),將時(shí)變系數(shù)ai(n)表示為一組基函數(shù)的線性組合[15-16],即
(2)
式中,aij為展開式時(shí)不變系數(shù),fj(n)為基函數(shù),m為基函數(shù)擴(kuò)展維數(shù)。
將式(2)代入式(1),得
(3)
令A(yù)=[a11a12…a1m…ap1ap2…apm]T,B=[B1B2…Bp],X=[x1x2…xN]。
(4)
則式(3)可改寫為
X=BA+e
(5)
利用基函數(shù)展開式將時(shí)變模型式(1)中原來p個(gè)時(shí)變系數(shù)的識(shí)別問題轉(zhuǎn)化為p×m個(gè)定常參數(shù)的識(shí)別,即將原來非平穩(wěn)過程的參數(shù)識(shí)別轉(zhuǎn)化為一個(gè)線性時(shí)不變系統(tǒng)的辨識(shí)。模型表示為矩陣形式后,由最小二乘法解出時(shí)不變參數(shù)A的估計(jì)值,即
(6)
頻域分析是從頻率角度來分析表面肌電信號的特征,方法是將時(shí)域信號經(jīng)過短時(shí)傅里葉轉(zhuǎn)換后得到表面肌電信號的頻譜或功率譜。頻域分析常用分析參數(shù)有MPF、MF[17-18]。
由于MPF、MF都是基于短時(shí)傅里葉分析,其時(shí)間和分辨率都是固定的,然而對于表面肌電信號來講,頻譜分布范圍較寬時(shí),難以找到一個(gè)合適的時(shí)間窗來分析,即其時(shí)間和頻率分辨率較低。
(7)
(8)
式中,pi為單位功率譜的估計(jì)值(功率譜密度),N為1/2信號采樣點(diǎn)數(shù)。
通過對比分析3種方法的疲勞指標(biāo)疲勞前、后的變化率(change rate,CR)來確定哪個(gè)特征參數(shù)值指標(biāo)對于肌肉疲勞反應(yīng)靈敏性更高,進(jìn)一步說明哪種方法適合用來表征肌肉疲勞。疲勞前、后相關(guān)特征參數(shù)值變化率定義[19]為
CR=(yi+M-yi)/yi×100%
(9)
式中,yi為疲勞前值,yi+M為疲勞后值。
經(jīng)過計(jì)算,疲勞前、后表面肌電信號的MPF分別為209.34、149.30 Hz,MF分別為132.95、104.15 Hz。由式(9)求得的MPF變化率為-28.68%,MF變化率為-21.66%。
由于MPF、MF是基于短時(shí)傅里葉原理,采樣時(shí)間較長、實(shí)時(shí)性和分辨率較差。為了能夠比較3種方法對肌肉疲勞的敏感性,以下MPF、MF特征值都是基于N=15 000點(diǎn)(t=10 s)數(shù)據(jù)長度計(jì)算?;跁r(shí)變AR模型要求采樣時(shí)間短、對時(shí)間反應(yīng)迅速等特點(diǎn),應(yīng)用時(shí)變參數(shù)AR模型處理疲勞前、后的前N=150點(diǎn)(t=0.1 s)數(shù)據(jù),可得到相應(yīng)的時(shí)變參數(shù)。經(jīng)過計(jì)算,疲勞前、后表面肌電信號的ARC1分別為-2.25、-3.09。由式(9)求得的ARC1變化率為37.39%。
時(shí)變自回歸模型法是以最小均方誤差擬合角度對表面肌電信號進(jìn)行分析,克服了經(jīng)典譜估計(jì)頻率分辨率低和方差性能差缺點(diǎn)。從相對于短時(shí)傅里葉變換分析來看,其克服了傳統(tǒng)提取頻域參數(shù)要求表面肌電信號采樣時(shí)間足夠長的問題。以采樣頻率1 500 Hz、時(shí)間窗512點(diǎn)舉例說明,計(jì)算MPF、MF參數(shù)時(shí)最短數(shù)據(jù)也需要512/1 500≈0.34 s>0.1 s。
實(shí)驗(yàn)采用AR參數(shù)模型提取特征參數(shù),模型階數(shù)一般取p=4~6,取Legendre基,Legendre基函數(shù)(Legendre basis function,LBF)的最佳維數(shù)由相關(guān)指數(shù)求得[20-22]。相關(guān)指數(shù)越大,說明模型擬合效果就越好,有
(10)
采用Ecel 2013軟件數(shù)據(jù)分析模塊。采用雙尾t檢驗(yàn)分析的方法來分析各個(gè)參數(shù)特征值對疲勞反應(yīng)靈敏性效果的統(tǒng)計(jì)學(xué)差異性,即分析各個(gè)參數(shù)特征值對肌肉疲勞程度表征效果的統(tǒng)計(jì)學(xué)差異性。
變化率數(shù)據(jù)統(tǒng)計(jì)步驟:要先進(jìn)行F檢驗(yàn),F(xiàn)檢驗(yàn)又叫方差齊性檢驗(yàn)。在雙樣本t檢驗(yàn)中要用到F檢驗(yàn),F(xiàn)檢驗(yàn)的目的是確定采用雙樣本等方差t檢驗(yàn)還是采用雙樣本異方差t檢驗(yàn)。如果F檢驗(yàn)單尾P值的2倍大于0.05,說明兩者方差之間無顯著性差異,則采用雙樣本等方差t檢驗(yàn)。否則,選擇雙樣本異方差t檢驗(yàn)。
圖1為一典型受試者肱二頭肌的表面肌電信號圖,可看出肌肉在疲勞前、后表面肌電信號幅度有增大趨勢,反映了肌肉收縮至疲勞時(shí)參加活動(dòng)的運(yùn)動(dòng)單位數(shù)量增加,體現(xiàn)為肌肉疲勞過程中表面肌電信號幅度增大。
圖1 表面肌電圖。(a)疲勞前期;(b)疲勞后期Fig.1 The sEMG of biceps brachi muscle. (a) The early stage of fatigue; (b) The later stage of fatigue
圖2為一典型受試者的基于基函數(shù)展開式法估計(jì)6階AR模型得到的相關(guān)指數(shù)L隨基函數(shù)維數(shù)變化,Legendre基函數(shù)最佳維數(shù)由式(10)求得。由圖2可以看出,基函數(shù)維數(shù)m≥4時(shí)參數(shù)辨識(shí)效果波動(dòng)不大,趨于穩(wěn)定。大量實(shí)驗(yàn)驗(yàn)證,p=6、m=7時(shí)的參數(shù)辨識(shí)效果最佳。
圖2 相關(guān)指數(shù)隨基函數(shù)維數(shù)變化;(a)疲勞前;(b)疲勞后Fig.2 The correlation index changing with LBF dimension;(a)Pre-fatigue;(b)Post-fatigue
圖3為一典型受試者基于Legendre展開式法時(shí)變參數(shù)辨識(shí)結(jié)果。可以看出,基于Legendre展開式法因具有良好的局部特性,能夠更好地跟蹤信號,辨識(shí)結(jié)果比較平滑,尤其是在波峰和波谷處沒有參數(shù)突變位置的現(xiàn)象。因此,基于Legendre展開式法具有理想的辨識(shí)效果。
表1 各受試者疲勞特征值及其變化率(疲勞前后)Tab.1 The characteristic value and change rate of the subjects(pre-fatigue and post-fatigue)
注:a:ARC1 vs MPF,P>0.05;b:ARC1 vs MF,P>0.05。
Note:a:ARC1 vs MPF,P<0.05;b:ARC1 vs MF,P<0.05.
圖3 基于Legendre展開式法時(shí)變參數(shù)辨識(shí)結(jié)果。(a)疲勞前;(b)疲勞后Fig.3 Identification results of AR parameters by LBF expansion method.(a)Pre-fatigue;(b)Post-fatigue
表1為10組受試者疲勞前、疲勞后相關(guān)特征值及其變化率大小。實(shí)驗(yàn)結(jié)果得到10組時(shí)變參數(shù),以每組AR模型第一個(gè)參數(shù)(ARC1)最重要,它直接體現(xiàn)當(dāng)前時(shí)刻與上一時(shí)刻關(guān)系,是肌肉狀態(tài)隨時(shí)間變化的最直接量,每階AR參數(shù)都是隨時(shí)間點(diǎn)變化。計(jì)算所估計(jì)信號0.1 s時(shí)間段內(nèi)ARC1為評價(jià)肌肉疲勞狀態(tài)指標(biāo)。計(jì)算所估計(jì)信號10 s時(shí)間段內(nèi)MPPF、MF為評價(jià)肌肉疲勞狀態(tài)指標(biāo)。
從疲勞前、疲勞后各特征參數(shù)變化趨勢上來看,ARC1、MPF、MF三者均呈現(xiàn)出變小趨勢,從疲勞前、后各特征參數(shù)變化率上來看,ARC1的變化率較高。采用AR模型克服了MF、MPF要求表面肌電信號采樣時(shí)間足夠長問題。將10例受試者疲勞前、后MPF、MF、ARC1變化率作為疲勞指標(biāo),變化率進(jìn)行統(tǒng)計(jì)分析,通過計(jì)算變化率均值和標(biāo)準(zhǔn)差分別為25.68%±2.03%、22.80%±2.19%、34.33%±2.41%。由于ARC1與MPF和ARC1與MF兩者兩兩之間存在顯著性差異,且分別顯著高于MPF和MF,表明采用ARC1參數(shù)指標(biāo)跟蹤肌肉疲勞具有較高分辨率靈敏度。
運(yùn)動(dòng)性疲勞學(xué)說有衰竭學(xué)說、堵塞學(xué)說、內(nèi)環(huán)境失調(diào)學(xué)說、保護(hù)性抑制學(xué)說、突變理論、自由基損傷學(xué)說等。Boyas等認(rèn)為,肌肉疲勞主要是由于發(fā)力過程中逐漸出現(xiàn)鈣離子運(yùn)動(dòng)擾動(dòng)、磷酸鹽等代謝物積累、肌細(xì)胞三磷酸腺苷減少等因素,導(dǎo)致動(dòng)作電位傳導(dǎo)變化和肌纖維收縮力量下降[23]。Pavlov等認(rèn)為,肌肉疲勞是大腦皮層保護(hù)性抑制結(jié)果。隨著疲勞加深,為了減少消耗,中樞神經(jīng)系統(tǒng)發(fā)放頻率降低,導(dǎo)致運(yùn)動(dòng)神經(jīng)元發(fā)放沖動(dòng)頻率隨之減小。中樞神經(jīng)對運(yùn)動(dòng)神經(jīng)元控制精度也降低,導(dǎo)致了運(yùn)動(dòng)單位間活動(dòng)同步化程度增強(qiáng),體現(xiàn)為肌肉疲勞過程中ARC1、MPF、MF下降。
ARC1針對表面肌電信號的非平穩(wěn)特性,采用AR模型對表面肌電信號進(jìn)行分析,對短時(shí)間內(nèi)表面肌電信號的肌肉疲勞迅速做出高靈敏性判定。由于ARC1在評估肌肉疲勞時(shí)具有靈敏性高、檢測肌肉疲勞時(shí)間短等優(yōu)點(diǎn),因此ARC1應(yīng)用于肌肉疲勞領(lǐng)域可以達(dá)到及時(shí)發(fā)現(xiàn)、及時(shí)防止由普通性肌肉疲勞進(jìn)一步惡化成肌肉損傷的效果,以免對個(gè)人、家庭及社會(huì)帶來巨大的經(jīng)濟(jì)負(fù)擔(dān)。
影響MPF、MF在疲勞過程中變化不明顯的影響因素有肌肉發(fā)力方式、肌肉表面溫度。在早期相關(guān)領(lǐng)域研究中,Tesch等研究發(fā)現(xiàn),肌肉以動(dòng)態(tài)收縮發(fā)力方式使肌肉發(fā)生疲勞過程中MPF、MF變化不明顯[24]。Petrofsky等研究發(fā)現(xiàn),肌肉收縮至疲勞過程中肌肉表面溫度會(huì)逐漸上升,這樣會(huì)增加表面肌電信號頻譜高頻成分,從而導(dǎo)致MPF、MF變化不明顯[25]。除了肌肉發(fā)力方式、肌肉表面溫度這兩個(gè)影響因素外,還有肌肉收縮強(qiáng)度也會(huì)使MPF、MF在疲勞過程中變化不明顯。呂飛等研究發(fā)現(xiàn),在不同運(yùn)動(dòng)強(qiáng)度(30%、55%、80% MVC)下,MPF與MF的下降率具有明顯的負(fù)荷強(qiáng)度效應(yīng)[26]。然而,當(dāng)運(yùn)動(dòng)強(qiáng)度過低(<20% MVC)時(shí),理論上將不會(huì)產(chǎn)生MPF、MF的任何變化[27]。
采用短時(shí)傅里葉變換的方法,雖然能夠得到對疲勞敏感的參數(shù)來表征肌肉疲勞,但短時(shí)傅里葉變換在時(shí)間和頻率上分辨率低。當(dāng)分析表面肌電信號頻譜分布范圍較寬時(shí),時(shí)間窗選擇就成了難題。時(shí)間窗選擇不佳,便會(huì)影響到對疲勞敏感的參數(shù)表征肌肉疲勞效果[28]。時(shí)變自回歸模型不受窗函數(shù)的影響這一優(yōu)點(diǎn)已經(jīng)被相關(guān)學(xué)者應(yīng)用到分析表面肌電信號領(lǐng)域。在早期相關(guān)研究中,Merletti等研究發(fā)現(xiàn),自回歸模型系數(shù)隨著肌肉疲勞程度加深,會(huì)出現(xiàn)逐漸變小趨勢[29]。Chang等研究發(fā)現(xiàn)自回歸模型系數(shù)的平方隨著肌肉疲勞程度加深會(huì)出現(xiàn)逐漸變小趨勢[30]。
本研究采用Legendre基函數(shù),對線性時(shí)變系統(tǒng)在白噪聲激勵(lì)下時(shí)變系統(tǒng)參數(shù)進(jìn)行展開,并利用相關(guān)指數(shù)來選擇Legendre基函數(shù)最佳維數(shù)。通過展開肌肉疲勞實(shí)驗(yàn)研究,各特征參數(shù)對肌肉疲勞表征效果方面進(jìn)行分析。通過實(shí)驗(yàn)將ARC1與傳統(tǒng)頻域參數(shù)MPF、MF進(jìn)行對比研究, 證明ARC1在評估肌肉疲勞領(lǐng)域的可行性及有效性。
本研究應(yīng)用Legendre基函數(shù)展開式方法將線性非平穩(wěn)過程辨識(shí)問題轉(zhuǎn)化為線性時(shí)不變系統(tǒng)辨識(shí)問題。對10例受試者表面肌電信號進(jìn)行特征提取,短時(shí)間內(nèi)ARC1變化率作為評估肌肉疲勞敏感性指標(biāo),實(shí)驗(yàn)證明,比MPF、MF對疲勞反應(yīng)的敏感性更高,從而有將細(xì)微特征信息放大,使不明顯的特征信息變明顯的效果。
通過表面肌電信號對肌肉疲勞檢測,不僅能夠更好表征肌肉疲勞變化程度,而且具有時(shí)間短和敏感性高等優(yōu)點(diǎn),可應(yīng)用于在線實(shí)時(shí)檢測肌肉疲勞,為上肢肌肉勞損評估、康復(fù)治療及人體工效學(xué)的研究提供一個(gè)可靠分析工具。時(shí)變參數(shù)辨識(shí)問題一直是學(xué)術(shù)界研究難點(diǎn),從10例受試者ARC1靈敏性均高于MPF和MF靈敏性可知,該研究方法具有適用性特性。為了進(jìn)一步驗(yàn)證及鞏固上述結(jié)果正確性及實(shí)用性,除了大量仿真驗(yàn)證該實(shí)驗(yàn)具有時(shí)間短、敏感性高等優(yōu)點(diǎn)之外,在接下來的研究中,后期仍需要進(jìn)行大量實(shí)驗(yàn)加以驗(yàn)證,從而進(jìn)一步提高其在肌肉疲勞判定等領(lǐng)域中的實(shí)用價(jià)值。