融合生物動力學和NMF的下肢運動建模

戴海烽，張 超，張建海

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

0 引 言

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

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

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

1.1 基本生理意義模型

1.1.1 意圖-肌電子模型

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

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

(1)

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

1.1.2 肌力-肌電子模型

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

(2)

式中，m，l，J分別為估算的小腿重量、重心和轉動慣量，X1，X2分別為受試者的體重和身高。θ為關節(jié)角度，g為重力加速度。下肢有用力矩L分為2項[11]，分別與當前膝關節(jié)角度、角度加速度有關，其余力矩屬于內力矩，起到調整運動姿態(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)為所獲得的激活信號。

獲得激活信號后，再通過激活信號-肌電轉化成肌電信號：

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

(4)

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

1.1.3 肌肉長度-肌電子模型

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

LM=LM0+rθ

(5)

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

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

(6)

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

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

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

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

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

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

在融合3子模型中，為了不失去各子模型的生理意義，系數(shù)搜索時，需要依據(jù)膝關節(jié)運動的特點合理修改損失函數(shù)。以股四頭肌為例，當只使用意圖-肌電模型時，側重運動意圖對肌肉收縮的控制，即在伸展階段，原始肌電與預測肌電間的誤差最低；在只使用肌力-肌電模型時，著重考慮主動肌收縮發(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為膝關節(jié)運動階段，t1表示意圖控制運動階段，在預測股四頭肌肌電時，其表示膝伸展階段，而預測股二頭肌肌電時則表示膝屈曲階段。t2表示肌肉主動收縮對外施展有用力矩的過程，在預測股四頭肌肌電時，其表示膝伸展階段和伸展回歸階段，而預測股二頭肌肌電時則表示膝屈曲階段和屈曲回歸階段。t3表示肌肉被動收縮而隨關節(jié)運動產生被動肌電的過程，即協(xié)同過程，在預測股四頭肌肌電時，其表示屈曲和屈曲回歸階段，而預測股二頭肌肌電時則表示膝伸展階段和伸展回歸階段。

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

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

V=WH+E

(10)

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

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

圖2 NMF的分解效果

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

3 肌電擬合實驗及分析

3.1 實驗設備

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

圖3 肌電設備和角度傳感器

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

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

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

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

3.3 實驗結果

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

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

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

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

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

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

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

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

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

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

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

(11)

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

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

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

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

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

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

4 結束語

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