戴海烽，張 超，張建海

(1.杭州電子科技大學計算機學院,浙江 杭州 310018；2.浙江中醫(yī)藥大學附屬第一醫(yī)院,浙江 杭州 310018)

0 引 言

越來越多可穿戴康復設(shè)備[1]引入了電刺激[2]反饋功能，電刺激反饋效果和肢體運動姿態(tài)與肌肉激活程度模型息息相關(guān)。肢體康復訓練中，根據(jù)健康人正常肢體運動姿態(tài)的肌肉激活程度模型，對患者特定肌肉群實施有針對性的電刺激，為肢體運動進行正確、科學的電刺激反饋提供了指導性的參考。王江等[3]運用徑向基(Radial Basis Function，RBF)神經(jīng)網(wǎng)絡(luò)對肌電信號進行建模仿真，該方法非常適合肌電信號非線性系統(tǒng)的建模，但RBF神經(jīng)網(wǎng)絡(luò)的理論基礎(chǔ)是統(tǒng)計數(shù)學，缺乏生理意義。李強[4]運用肌纖維募集模型仿真第一背間肌運動神經(jīng)元在不同激勵水平下的表面肌電(Surface ElectroMyoGraphy，SEMG)信號和相應(yīng)的收縮能力，研究發(fā)現(xiàn)仿真的動作電位信號與真實信號相似性極高。James等[5]通過研究肌梭反饋模型發(fā)現(xiàn)，反饋信號與肌肉長度之間具有非常高的相關(guān)性。雖然，這些生理學模型能夠描述肌肉肌電與關(guān)節(jié)運動的物理關(guān)系，但是存在一定的片面性，擬合效果不如多模型融合方法。近年來，多模型融合方法開始運用于主動肌建模，白杰[6]使用一種由內(nèi)驅(qū)力、高爾基腱、肌梭融合的模型擬合肌電，擬合誤差與神經(jīng)網(wǎng)絡(luò)方法相當，但由于部分子模型為數(shù)學模型，擬合參數(shù)過多，導致3個子模型擬合的肌電成分不穩(wěn)定，失去子模型該有的生理意義。在構(gòu)建肢體運動姿態(tài)的肌肉激活模型時，多數(shù)以主動肌為主，輔助肌的作用同樣不可忽視，對完善關(guān)節(jié)運動與肌肉群間的研究同樣重要。協(xié)同分解[7]，例如非負矩陣分解(Nonnegative Matrix Factorization，NMF)[8]常常用于分析肌肉間隱含的協(xié)同關(guān)系。李飛[9]將NMF應(yīng)用于小兒腦癱的下肢步態(tài)分析中，為患者的康復訓練提供評估和指導。本文使用運動配置策略以及膝關(guān)節(jié)運動相關(guān)損失函數(shù)改進融合生物動力學模型，采用NMF對輔助肌肌電進行建模，構(gòu)建一個具有生理意義、全面性、肌電擬合效果好的下肢膝關(guān)節(jié)肌電模型。

1 融合生物動力學的下肢主動肌建模

在構(gòu)建下肢運動姿態(tài)與肌肉激活關(guān)系模型時，為了使子模型生成較穩(wěn)定的肌電，需選擇模型系數(shù)較少或有限定范圍的生理意義模型。同時，為了進一步加強各個子模型的穩(wěn)定性，本文使用模型配置策略，根據(jù)下肢膝關(guān)節(jié)的運動特點使用特定的損失函數(shù)進行模型系數(shù)的搜索。

1.1 基本生理意義模型

1.1.1 意圖-肌電子模型

意圖-肌電子模型根據(jù)運動意圖來控制肌肉收縮，肌肉收縮時產(chǎn)生肌肉肌電，即運動意圖生成肌電，其過程如下：

ΦRTE=eα i
u(t)=wΦRTE

(1)

式中，ΦRTE為運動神經(jīng)單元募集閾值估計(Raising Threshold Estimation，RTE)，w為放電增益，u(t)為最終的意圖肌電。i為募集運動神經(jīng)單元的指數(shù)系數(shù)，α為標量化的神經(jīng)意圖，一般無法直接測量，可通過其與角度偏差(Angle Error，AE)的關(guān)系求得，即α=kAAE，角度偏差可通過AAE=As+n+As-n-2As求得，即將n時刻前后的角度求和，再與2倍的當前時刻s下的角度As相減，k為神經(jīng)意圖與角度偏差的關(guān)系系數(shù)。

1.1.2 肌力-肌電子模型

肌力-肌電子模型需先計算主動肌產(chǎn)生的有用力矩信號，依據(jù)《中國成年人人體慣性參數(shù)國家標準》[10]中相關(guān)的回歸方程，有用力矩的計算公式為：

(2)

式中，m，l，J分別為估算的小腿重量、重心和轉(zhuǎn)動慣量，X1，X2分別為受試者的體重和身高。θ為關(guān)節(jié)角度，g為重力加速度。下肢有用力矩L分為2項[11]，分別與當前膝關(guān)節(jié)角度、角度加速度有關(guān)，其余力矩屬于內(nèi)力矩，起到調(diào)整運動姿態(tài)作用，可忽略。

肌力-肌電子模型選擇逆hill模型。hill模型分為肌電-激活信號和激活信號-肌力，逆hill模型也可分為這2部分的相反過程。肌力-激活信號表示如下：

a(t)fe+fp=F/Fmax

(3)

式中，fe為肌肉主動收縮部分的系數(shù)，fp為肌肉被動收縮部分系數(shù)，可忽略不計，F(xiàn)為力矩信號，F(xiàn)max為力矩信號中的最大值，a(t)為所獲得的激活信號。

獲得激活信號后，再通過激活信號-肌電轉(zhuǎn)化成肌電信號：

g(t)=vln[a(t)(eB-1)+1]/B

(4)

式中，g(t)為預測的肌電信號，v為擬合系數(shù)，B為非線性參數(shù)，取值范圍為[-3，0]，該數(shù)值與動作的種類有關(guān)。

1.1.3 肌肉長度-肌電子模型

肌肉長度-肌電子模型選用肌梭模型，肌梭是肌肉中的肌肉長度感知器官，所以肌肉長度-肌電子模型是模仿肌梭反饋功能，感知肌肉長度變化并反饋與肌肉長度相關(guān)的肌電的模型。首先，針對股四頭肌和股二頭肌需要建立不同的肌肉長度-角度模型，其中股四頭肌的肌肉長度-角度公式為：

LM=LM0+rθ

(5)

式中，LM為股四頭肌肌肉長度，LM0為股四頭肌初始肌肉長度，r為膝關(guān)節(jié)中心的旋轉(zhuǎn)半徑。

股二頭肌的肌肉長度-角度公式為：

(6)

式中，LB為股二頭肌肌肉長度，L11為股二頭肌起始端到膝關(guān)節(jié)旋轉(zhuǎn)中心的距離，L12為膝關(guān)節(jié)旋轉(zhuǎn)中心至股二頭肌止端的距離。

股四頭肌肌梭模型處理過程如下：

uM-uM0=KMvp(LM-LM0)

(7)

式中，uM為股四頭肌肌梭反饋的肌電信號，uM0為初始股四頭肌肌梭反饋的肌電信號，v為肌肉長度變化速度，p為肌肉長度變化速度影響因子，通常取0.3，KM為待定系數(shù)。

類似地，股二頭肌肌梭模型處理過程如下：

uB-uB0=KMvp(LB-LB0)

(8)

式中，uB為股二頭肌肌梭反饋的肌電信號，uB0為初始股二頭肌肌梭反饋的肌電信號，LB0為股二頭肌初始肌肉長度。

1.2 運動配置策略

為了加強1.1節(jié)的3個子模型的穩(wěn)定性，本文采用運動配置策略將各個子模型進行融合。首先，將下肢膝關(guān)節(jié)運動分成4個階段，分別為膝關(guān)節(jié)伸展、膝關(guān)節(jié)伸展回歸、膝關(guān)節(jié)屈曲和膝關(guān)節(jié)屈曲回歸，如圖1所示。其次，分析各個階段的運動過程和肌電產(chǎn)生特點，以股四頭肌為例，相比于其他過程，伸展階段中，運動神經(jīng)控制股四頭肌收縮的意圖最強，意圖肌電成分所占比重最大；伸展和伸展回歸階段中，股四頭肌主動收縮，控制小腿隨膝關(guān)節(jié)中心旋轉(zhuǎn)并對外做有用功，該過程生成的肌電大部分由主動肌力產(chǎn)生；屈曲和屈曲回歸階段中，股四頭肌起協(xié)同肌作用，只隨膝關(guān)節(jié)運動，使自身長度發(fā)生變化，起到調(diào)整運動姿勢的作用，該過程生成的肌電基本由因肌肉長度變化產(chǎn)生的被動肌力。最后，依據(jù)膝關(guān)節(jié)運動特點，在伸展階段配置3個子模型，在伸展回歸階段配置肌力-肌電子模型和肌肉長度-肌電子模型，在屈曲和屈曲回歸配置肌肉長度-肌電子模型。

圖1 膝關(guān)節(jié)運動中的4個階段

與股四頭肌類似，股二頭肌也依據(jù)對應(yīng)膝關(guān)節(jié)運動特點，在屈曲階段配置3個子模型，在屈曲回歸階段配置肌力-肌電子模型和肌肉長度-肌電子模型，在伸展和伸展回歸配置肌肉長度-肌電子模型。

1.3 膝關(guān)節(jié)運動相關(guān)損失函數(shù)

在融合3子模型中，為了不失去各子模型的生理意義，系數(shù)搜索時，需要依據(jù)膝關(guān)節(jié)運動的特點合理修改損失函數(shù)。以股四頭肌為例，當只使用意圖-肌電模型時，側(cè)重運動意圖對肌肉收縮的控制，即在伸展階段，原始肌電與預測肌電間的誤差最低；在只使用肌力-肌電模型時，著重考慮主動肌收縮發(fā)力對外做功時的情形；在只使用肌肉長度-肌電模型時，著重考慮非主動肌被動發(fā)力時的情況。當需要使用融合模型時，需要同時兼顧3個模型對4個階段不同的參與度，在伸展階段，意圖、肌力、肌肉長度模型都在肌電預測方面具有一定的參與度；在伸展回歸階段，大腦處于放松階段，意圖模型不再參與預測肌電；在屈曲和屈曲回歸階段，股四頭肌不再是主動肌，只有肌肉長度模型參與肌電預測。最后，通過最小化損失函數(shù)求解得到3個模型中的待定系數(shù)，使得模型擬合的肌電更趨近于各個子模型的生理意義。綜合上述分析，融合模型的損失函數(shù)為：

(9)

式中，Sn(t)為各個子模型預測的肌電信號與真實信號間的損失函數(shù)，其中n為1表示意圖-肌電子模型，2表示肌力-肌電子模型，3表示肌肉長度-肌電子模型，SF(t)為融合生物動力學模型與真實信號間的損失函數(shù)，F(xiàn)表示融合生物動力學模型。uSEMG為真實的表面肌電信號，un為各個子模型或者融合生物動力學模型預測的肌電。tn為膝關(guān)節(jié)運動階段，t1表示意圖控制運動階段，在預測股四頭肌肌電時，其表示膝伸展階段，而預測股二頭肌肌電時則表示膝屈曲階段。t2表示肌肉主動收縮對外施展有用力矩的過程，在預測股四頭肌肌電時，其表示膝伸展階段和伸展回歸階段，而預測股二頭肌肌電時則表示膝屈曲階段和屈曲回歸階段。t3表示肌肉被動收縮而隨關(guān)節(jié)運動產(chǎn)生被動肌電的過程，即協(xié)同過程，在預測股四頭肌肌電時，其表示屈曲和屈曲回歸階段，而預測股二頭肌肌電時則表示膝伸展階段和伸展回歸階段。

2 基于NMF的下肢輔助肌建模

NMF常常用于分析肌肉間的協(xié)同情況[12]，與神經(jīng)系統(tǒng)控制肌肉運動的分析過程相似，即在進行某個動作時，神經(jīng)系統(tǒng)只控制一個更小的單位的肌肉協(xié)同模式而非一整塊肌肉，分析式如下：

V=WH+E

(10)

式中，V為待分解的肌電矩陣，W為分解后的協(xié)同曲線，H為協(xié)同模式，E為噪聲誤差。使用迭代算法可使得噪聲誤差逐步趨近0，最基本的是文獻[13]提出的基于歐氏距離的乘性迭代算法和基于K-L散度的加性迭代算法。本文采用基于歐氏距離乘性迭代的NMF。

運用NMF分解的主、副肌肉群間的協(xié)同效果如圖2所示，8個肌肉通道分別對應(yīng)闊筋膜張肌、股四頭肌、股二頭肌、內(nèi)收肌、腓腸肌、脛骨前肌、比目魚肌內(nèi)側(cè)、比目魚肌外側(cè)，協(xié)同曲線主模式分別為股四頭肌、股二頭肌以及闊筋膜張肌。

圖2 NMF的分解效果

協(xié)同模式表示1組以某塊主動肌為主要作用，其余肌肉為輔助作用的肌肉協(xié)同模式，對應(yīng)的協(xié)同曲線表示該協(xié)同模式的能量隨著時間變化的情況。分解肌電后，先將股四頭肌與闊筋膜張肌對應(yīng)協(xié)同模式上的特征值相除，再將該值乘以股四頭肌肌電，將股四頭肌肌電轉(zhuǎn)化為闊筋膜張肌肌電，最后，疊加各自的協(xié)同模式，得到其他輔助肌肌電。

3 肌電擬合實驗及分析

3.1 實驗設(shè)備

實驗選用8通道博瑞康無線肌電采集設(shè)備，采樣頻率為1 000 Hz，選取肌肉通道為闊筋膜張肌、內(nèi)收肌、腓腸肌、脛骨前肌、比目魚肌內(nèi)側(cè)、比目魚肌外側(cè)以及2塊主動肌。角度傳感器選用維特智能九軸角度傳感器記錄儀，采樣頻率為200 Hz，如圖3所示。

圖3 肌電設(shè)備和角度傳感器

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

選取4名健康男性受試者，均已簽署知情同意書。采集受試者的身高、體質(zhì)量、靜息主動肌長度、大腿長度等生理參數(shù)。

下肢膝關(guān)節(jié)動作設(shè)計如圖4所示。受試者先保持坐姿平衡，然后做出膝關(guān)節(jié)伸展運動，再返回至坐姿平衡狀態(tài)，隨后進行膝關(guān)節(jié)屈曲運動，再次返回平衡狀態(tài)，共做5次。做完1個輪次后，休息5 min，避免肌肉疲勞。

圖4 下肢膝關(guān)節(jié)運動實驗過程

3.3 實驗結(jié)果

在MATLAB實驗環(huán)境下，采用本文方法對3.2節(jié)所述的膝關(guān)節(jié)運動實驗中的主動/輔助肌肌電進行擬合。采用本文提出的融合生物動力學方法與文獻[6]提出的主動肌的融合方法進行股四頭肌肌電的擬合實驗，文獻[6]模型的擬合結(jié)果如圖5所示，本文方法各子模型的擬合結(jié)果如圖6所示。

圖5 文獻[6]模型擬合股四頭肌肌電的結(jié)果

圖6 本文方法子模型擬合股四頭肌肌電的結(jié)果

從圖5中可以看出，文獻[6]融合模型的肌電擬合效果不錯，但各子模型的肌電擬合不穩(wěn)定，這是因為子模型的擬合參數(shù)過多，致使子模型的肌電成分不穩(wěn)定，不符合生理意義。

從圖6可以看出，意圖-肌電模型在伸展階段擬合情況較好，但其余運動階段模型預測的肌電都高于原始肌電，這是因為除伸展階段外的運動階段過多考慮了意圖的參與，擬合的肌電相對偏高。類似地，肌力-肌電模型在伸展階段未考慮動作調(diào)整意圖對肌電信號的影響，在屈曲和屈曲回歸階段，肌力-肌電模型又過多考慮主動肌力的參與，致使伸展階段預測的肌電相對偏低，屈曲和屈曲回歸階段相對偏高。肌肉長度-肌電模型在伸展和伸展回歸階段未考慮主動肌力和動作調(diào)整意圖對肌電的影響，致使肌肉長度-肌電模型預測的肌電在伸展和伸展回歸階段都相對偏低。

采用本文融合模型擬合股四頭肌和股二頭肌肌電的結(jié)果如圖7所示。

圖7 本文方法的各個子模型擬合主動肌肌電的結(jié)果

從圖7可以看出，融合生物動力學模型中，各子模型的肌電擬合比較穩(wěn)定，且3個子模型擬合的肌電疊加起來的外包絡(luò)也能很好擬合原始肌電。相比文獻[6]模型，本文采用模型配置策略，并運用運動相關(guān)損失函數(shù)來搜尋模型中的待定系數(shù)，使得各子模型擬合與各自模型意義相關(guān)的肌電，解決了子模型肌電擬合不穩(wěn)定問題。

在主動肌肌電模型的基礎(chǔ)上，采用本文提出的NMF求得的協(xié)同模式對各個輔助肌進行肌電擬合，得到輔助肌模型擬合結(jié)果如圖8所示。

從圖8可以看出，輔助肌肌電擬合情況良好，說明運用主動肌肌電和NMF算法分解的協(xié)同模式能有效預測輔助肌肌電，完善了下肢膝關(guān)節(jié)運動姿態(tài)與相應(yīng)肌肉群激活程度的關(guān)系。

本文采用均方根誤差(Root Mean Square Error，RMSE)作為定量評價擬合效果的指標。

(11)

式中，d為對比的兩段肌電信號的長度。

本文各子模型及融合模型、文獻[6]模型的主動肌肌電擬合效果如表1所示，各塊輔助肌肌電擬合效果如表2所示。

表1 主動肌建模效果定量評價表 單位：%

表2 輔助肌建模效果定量評價表 單位：%

從表1可以看出，本文融合模型的RMSE值與文獻[6]模型相差不大，說明兩者的擬合效果相當；本文融合模型的RMSE值都低于單獨使用子模型擬合的結(jié)果，說明融合模型擬合效果優(yōu)于單個子模型。

將表2數(shù)據(jù)取平均值，得到輔助肌的平均RMSE值為4.62%。一般情況下，RMSE值低于8%時，預測的信號能很好地擬合原始信號[14]。因此，運用NMF對輔助肌進行肌肉激活程度建?？梢缘玫搅己玫募‰姅M合效果。

4 結(jié)束語

本文提出一種融合生物動力學和NMF的下肢運動建模方法，建立了下肢膝關(guān)節(jié)運動姿態(tài)與肌肉激活程度的關(guān)系模型。運用運動配置策略，采用運動相關(guān)損失函數(shù)，加強了子模型的穩(wěn)定性，并建立基于NMF的輔助肌模型，完善了下肢關(guān)節(jié)運動與肌肉激活程度模型，為正確、有效進行功能電刺激提供借鑒。后續(xù)將嘗試更多的分解方式以提升輔助肌的擬合效果，進一步優(yōu)化下肢關(guān)節(jié)運動與肌肉激活程度模型。