戴海烽，張 超，張建海

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

0 引 言

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

1 融合生物動(dòng)力學(xué)的下肢主動(dòng)肌建模

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

1.1 基本生理意義模型

1.1.1 意圖-肌電子模型

意圖-肌電子模型根據(jù)運(yùn)動(dòng)意圖來(lái)控制肌肉收縮，肌肉收縮時(shí)產(chǎn)生肌肉肌電，即運(yùn)動(dòng)意圖生成肌電，其過(guò)程如下：

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

(1)

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

1.1.2 肌力-肌電子模型

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

(2)

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

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

a(t)fe+fp=F/Fmax

(3)

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

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

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

(4)

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

1.1.3 肌肉長(zhǎng)度-肌電子模型

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

LM=LM0+rθ

(5)

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

股二頭肌的肌肉長(zhǎng)度-角度公式為：

(6)

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

股四頭肌肌梭模型處理過(guò)程如下：

uM-uM0=KMvp(LM-LM0)

(7)

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

類(lèi)似地，股二頭肌肌梭模型處理過(guò)程如下：

uB-uB0=KMvp(LB-LB0)

(8)

式中，uB為股二頭肌肌梭反饋的肌電信號(hào)，uB0為初始股二頭肌肌梭反饋的肌電信號(hào)，LB0為股二頭肌初始肌肉長(zhǎng)度。

1.2 運(yùn)動(dòng)配置策略

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

圖1 膝關(guān)節(jié)運(yùn)動(dòng)中的4個(gè)階段

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

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

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

(9)

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

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

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

V=WH+E

(10)

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

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

圖2 NMF的分解效果

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

3 肌電擬合實(shí)驗(yàn)及分析

3.1 實(shí)驗(yàn)設(shè)備

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

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

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

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

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

圖4 下肢膝關(guān)節(jié)運(yùn)動(dòng)實(shí)驗(yàn)過(guò)程

3.3 實(shí)驗(yàn)結(jié)果

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

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

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

從圖5中可以看出，文獻(xiàn)[6]融合模型的肌電擬合效果不錯(cuò)，但各子模型的肌電擬合不穩(wěn)定，這是因?yàn)樽幽Ｐ偷臄M合參數(shù)過(guò)多，致使子模型的肌電成分不穩(wěn)定，不符合生理意義。

從圖6可以看出，意圖-肌電模型在伸展階段擬合情況較好，但其余運(yùn)動(dòng)階段模型預(yù)測(cè)的肌電都高于原始肌電，這是因?yàn)槌煺闺A段外的運(yùn)動(dòng)階段過(guò)多考慮了意圖的參與，擬合的肌電相對(duì)偏高。類(lèi)似地，肌力-肌電模型在伸展階段未考慮動(dòng)作調(diào)整意圖對(duì)肌電信號(hào)的影響，在屈曲和屈曲回歸階段，肌力-肌電模型又過(guò)多考慮主動(dòng)肌力的參與，致使伸展階段預(yù)測(cè)的肌電相對(duì)偏低，屈曲和屈曲回歸階段相對(duì)偏高。肌肉長(zhǎng)度-肌電模型在伸展和伸展回歸階段未考慮主動(dòng)肌力和動(dòng)作調(diào)整意圖對(duì)肌電的影響，致使肌肉長(zhǎng)度-肌電模型預(yù)測(cè)的肌電在伸展和伸展回歸階段都相對(duì)偏低。

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

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

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

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

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

本文采用均方根誤差(Root Mean Square Error，RMSE)作為定量評(píng)價(jià)擬合效果的指標(biāo)。

(11)

式中，d為對(duì)比的兩段肌電信號(hào)的長(zhǎng)度。

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

表1 主動(dòng)肌建模效果定量評(píng)價(jià)表 單位：%

表2 輔助肌建模效果定量評(píng)價(jià)表 單位：%

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

將表2數(shù)據(jù)取平均值，得到輔助肌的平均RMSE值為4.62%。一般情況下，RMSE值低于8%時(shí)，預(yù)測(cè)的信號(hào)能很好地?cái)M合原始信號(hào)[14]。因此，運(yùn)用NMF對(duì)輔助肌進(jìn)行肌肉激活程度建模可以得到良好的肌電擬合效果。

4 結(jié)束語(yǔ)

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