基于OpenSim人體步行腓腸肌靜態(tài)生物力學(xué)分析

趙 闖， 蔡玉強(qiáng)

(華北理工大學(xué)機(jī)械工程學(xué)院，唐山 063000)

近年來，在傷病預(yù)防以及人體下肢康復(fù)機(jī)器人的設(shè)計(jì)中需要考慮人體下肢的生物力學(xué)狀態(tài)[1-2]。在研究人體下肢康復(fù)機(jī)器人的模型構(gòu)建與動(dòng)力學(xué)分析中，下肢肌肉的運(yùn)動(dòng)狀態(tài)是進(jìn)行仿真的第一環(huán)節(jié)。然而受到人體生物系統(tǒng)的復(fù)雜性、外界環(huán)境的干擾以及倫理等的限制，現(xiàn)有的數(shù)據(jù)采集方式難以良好地獲取腿部肌肉受力的情況[3]。

OpenSim是Stanford University開發(fā)的一個(gè)開源軟件平臺(tái)[4]，用于創(chuàng)建、更改和分析骨骼肌肉系統(tǒng)的計(jì)算機(jī)模型以及運(yùn)動(dòng)的動(dòng)態(tài)仿真。文獻(xiàn)[3]基于OpenSim軟件搭建了人體肌肉骨骼仿真平臺(tái)，驗(yàn)證了人體在直立狀態(tài)下會(huì)根據(jù)外界環(huán)境的擾動(dòng)選擇性的使人體保持平衡。文獻(xiàn)[5]基于OpenSim建立肌肉骨骼模型，提出了計(jì)算肌肉拉力和膝關(guān)節(jié)剛度的方法。文獻(xiàn)[6]構(gòu)建簡(jiǎn)化人體下肢力學(xué)模型，建立相應(yīng)動(dòng)力學(xué)方程，經(jīng)過實(shí)驗(yàn)數(shù)據(jù)分析得到人體下肢膝關(guān)節(jié)與踝關(guān)節(jié)的受力情況，但在模型受力曲線上會(huì)出現(xiàn)數(shù)據(jù)不連續(xù)的情況。文獻(xiàn)[7]基于ADAMS 動(dòng)力學(xué)仿真軟件中建立人-機(jī)耦合系統(tǒng)模型并實(shí)現(xiàn)該系統(tǒng)的運(yùn)動(dòng)仿真，得到特定人體按照不同足部軌跡運(yùn)動(dòng)時(shí)下肢肌肉伸縮量的變化情況。文獻(xiàn)[8]研究人體在運(yùn)動(dòng)過程中，肌肉長(zhǎng)度變化、肌肉長(zhǎng)度變化速率、肌肉力臂等肌肉功能參數(shù)以及肌肉功能的動(dòng)態(tài)變化。

總之，中外有不少學(xué)者對(duì)人體下肢關(guān)節(jié)運(yùn)動(dòng)的形態(tài)特征進(jìn)行研究和分析，但對(duì)于下肢運(yùn)動(dòng)狀態(tài)下的靜態(tài)力學(xué)分析等方面研究不足。本文首先分析人體下肢受力機(jī)理，簡(jiǎn)述仿真所需骨骼肌肉模型，利用人工合成的運(yùn)動(dòng)學(xué)方法依賴于生物力學(xué)知識(shí)，并結(jié)合正向和逆向運(yùn)動(dòng)學(xué)知識(shí)[9]，使用OpenSim的Static Optimization Tool模塊，簡(jiǎn)化腰部肌肉模型為扭矩執(zhí)行器，改善肌肉激活度和靜態(tài)優(yōu)化輸出的力，利用逆向運(yùn)動(dòng)學(xué)(IK)、殘差縮減(RRA)和肌肉計(jì)算控制(CMC)對(duì)人體下肢肌肉進(jìn)行靜態(tài)優(yōu)化分析，為人體下肢康復(fù)機(jī)器人的研究與研發(fā)提供基礎(chǔ)。

1 人體肌肉骨骼模型

1.1 人體下肢受力

人體的膝關(guān)節(jié)和踝關(guān)節(jié)與小腿部肌肉有著復(fù)雜的肌骨結(jié)構(gòu)，涉及肌肉、韌帶、骨骼、神經(jīng)和血管等，是一個(gè)復(fù)雜的神經(jīng)運(yùn)動(dòng)系統(tǒng)。人體通過前庭、視覺、感受和觸覺的綜合分析處理，控制著肌肉完成身體運(yùn)動(dòng)的平衡[5]。因此，在研究人體行走過程中下肢的受力狀態(tài)，不僅要考慮足底受力和其他質(zhì)量，還要考慮關(guān)節(jié)肌肉的受力[10-11]。

圖1 人體下肢骨骼-肌肉模型Fig.1 Human lower limb bone-muscle model

人體站立、行走和慢跑等保持身體平衡的運(yùn)動(dòng)小腿部分主要是通過腓腸肌、比目魚肌和脛骨前肌之間的相互調(diào)節(jié)作用來實(shí)現(xiàn)的，如圖1人體下肢骨骼-肌肉模型。腓腸肌的作用是曲踝關(guān)節(jié)和曲膝關(guān)節(jié)，于皮膚表面，其深方是比目魚肌。比目魚肌起于脛骨和腓骨頭部，尾下為跟腱，止于跟底，受脛神經(jīng)支配。在運(yùn)動(dòng)過程中，腓腸肌和比目魚肌作為一體，維持身體運(yùn)動(dòng)過程中小腿部肌肉的肌力，在站姿或膝關(guān)節(jié)伸直、蹬足時(shí)發(fā)力的作用。脛骨前肌主要參與膝關(guān)節(jié)和踝關(guān)節(jié)的背曲運(yùn)動(dòng)，與腓腸肌和比目魚肌協(xié)調(diào)作用，維持人體平衡。

1.2 肌肉反射模型

將人體的下肢運(yùn)動(dòng)受力模型簡(jiǎn)化基于圖2骨骼-肌肉模型簡(jiǎn)圖，該模型適用于腿部骨骼-肌肉近似的平面運(yùn)動(dòng)，上肢同時(shí)保持近似垂直的狀態(tài)。其中人體上半身簡(jiǎn)化為固定在臀部的質(zhì)量塊；大腿、小腿和腳簡(jiǎn)化為3個(gè)連桿，髖關(guān)節(jié)、膝關(guān)節(jié)和踝關(guān)節(jié)簡(jiǎn)化為鉸接點(diǎn)；右腿肌肉簡(jiǎn)化為2對(duì)作用在連桿兩端的相對(duì)作用力，分布于大腿和小腿兩側(cè)。

圖2 骨骼-肌肉右腿模型簡(jiǎn)圖Fig.2 Skeleton-muscle right leg model diagram

由于肌肉運(yùn)動(dòng)的非線性，希爾模型[12-13]包含主動(dòng)收縮單元(CE)和被動(dòng)彈簧單元(SE)表示的肌肉模型，用非線性彈簧表示肌腱，如圖3所示。lMT為肌肉-肌腱長(zhǎng)度，lT為肌腱長(zhǎng)度，lM為標(biāo)準(zhǔn)化的肌纖維長(zhǎng)度，α為羽狀角，它們之間存在lMT=lT+lMcosα的關(guān)系。肌纖維收縮產(chǎn)生的拉力通過肌腱傳遞給所附著的骨骼，引起關(guān)節(jié)運(yùn)動(dòng)，力量的傳遞與α有直接關(guān)系，α越小傳遞越明顯，反之亦然；但是，一些主要肌肉為羽狀肌，主因在于相同體積羽狀肌的生理橫截面積更大[14]。

圖3 希爾模型Fig.3 Hill model

研究單步步態(tài)試驗(yàn)中執(zhí)行靜態(tài)優(yōu)化分析，并在簡(jiǎn)化的示例模型中評(píng)估腓腸肌激活度和力量產(chǎn)生。靜態(tài)優(yōu)化結(jié)果的質(zhì)量在很大程度上取決于輸入：模型、運(yùn)動(dòng)和力等參數(shù)。靜態(tài)優(yōu)化輸入為肌纖維長(zhǎng)度的變化值v(t)，肌肉反射控制模型的表達(dá)式為

R(t)=Kl[l(t)-l0]+Kvv(t)

(1)

v(t)=[l(tn+1)-l(tn)]/(tn+1-tn)

(2)

式中：R(t) 為t時(shí)刻肌肉激活度控制值；Kl為肌肉纖維長(zhǎng)度變化系數(shù)；l0為靜態(tài)下肌腱長(zhǎng)度；Kv為肌肉纖維速度變化系數(shù)；v(t) 為肌肉纖維長(zhǎng)度變化速率；l(t)為t時(shí)刻肌肉纖維長(zhǎng)度變化量；tn為步態(tài)的時(shí)間。

為了建立小腿部腓腸肌的控制方程，對(duì)腳部進(jìn)行受力分析，幾何標(biāo)注如圖4所示。

圖4 腳部受力示意圖Fig.4 Schematic diagram of the force on the foot

對(duì)于腳部質(zhì)心A分別建立x、y方向上的移動(dòng)分量[xA(t)、yA(t)]和轉(zhuǎn)動(dòng)分量[αy(t)]的3個(gè)運(yùn)動(dòng)方程。

(3)

(4)

(5)

式中：mf為腳部質(zhì)量；Fax、Fay為踝關(guān)節(jié)x、y方向的分力；F1x、F1y為腓腸肌x、y方向的分力；F2x、F2y為脛骨前肌x、y方向的分力；Fgx、Fgy為腳部與地面之間摩擦力分力；Gf為腳部自身的重力；Gt為單腳承受人體50%的重力；Jf為腳部轉(zhuǎn)動(dòng)慣量；Max、May為踝關(guān)節(jié)x、y方向的力矩；M1、M2為腓腸肌和脛骨前肌力矩；Mgx、Mgy為x、y方向摩擦力矩。

其中腳部質(zhì)量mf和轉(zhuǎn)動(dòng)慣量Jf等生物力學(xué)參數(shù)可由經(jīng)驗(yàn)方程[15]得到，腳部壓力可通過儀器得到[16]。式(3)～式(5)可以計(jì)算踝關(guān)節(jié)在各矢量面上Fa的分力，并且可以推導(dǎo)出腓腸肌受力F1的表達(dá)式：

(6)

將得到的數(shù)據(jù)代入式(6)，即得到腓腸肌受力與時(shí)間之間的對(duì)應(yīng)關(guān)系。

1.3 生物力學(xué)參數(shù)確定

由于難以獲取腳部的生物力學(xué)數(shù)據(jù)，在模型中將質(zhì)心作為腳部幾何中心，并將看作剛體，腳部質(zhì)心位置lf由人體腳部大小確定，轉(zhuǎn)動(dòng)慣量Jf、腳部重量Gt和人體重量Gf由人體解剖學(xué)參數(shù)確定。

2 模型動(dòng)態(tài)仿真試驗(yàn)

OpenSim作為建立和分析骨骼-肌肉的開源軟件，可通過軟件建立人體骨骼-肌肉模型，設(shè)置人體參數(shù)來計(jì)算和利用OpenSim的plot功能繪制肌肉力、膝、踝關(guān)節(jié)的運(yùn)動(dòng)位移曲線和力矩等數(shù)據(jù)，進(jìn)行人體下肢步行狀態(tài)下肌肉靜態(tài)力學(xué)分析。

2.1 人體肌肉骨骼搭建

研究采用的骨骼-肌肉模型如圖5所示，具有10個(gè)自由度和18個(gè)肌肉的平面縮小步態(tài)模型，模擬正常人體步行運(yùn)動(dòng)狀態(tài)。由于研究不涉及上肢運(yùn)動(dòng)，經(jīng)簡(jiǎn)化保留軀干部分，并在分析過程中簡(jiǎn)化為質(zhì)量塊處理。模型下肢經(jīng)簡(jiǎn)化后保留支撐人體站立的脛骨前肌、腓腸肌和比目魚肌[16]。

圖5 人體骨骼-肌肉模型Fig.5 Human bone-muscle model

利用Notepad++文本編輯器軟件修改原有模型的內(nèi)部參數(shù)，使該新模型subject01.osim具有代表實(shí)驗(yàn)參與者的質(zhì)量、人體身高尺寸和強(qiáng)度；自由度和肌肉幾何形狀適合于本次研究的肌肉群——腓腸肌；在腰部添加扭矩執(zhí)行器，保證人體上肢有足夠的驅(qū)動(dòng)控制；在腓腸肌添加備用執(zhí)行器，肌肉不能產(chǎn)生足夠的加速度時(shí)可以為步態(tài)周期的部分期間增加額外的動(dòng)力；添加殘差執(zhí)行器，測(cè)量分析程中運(yùn)動(dòng)和力之間的關(guān)系，確保整個(gè)分析滿足牛頓第二定律。模擬運(yùn)動(dòng)包含平滑、真實(shí)的加速度，并能精確地測(cè)量運(yùn)動(dòng)期間的肌肉的激活度與力矩。圖6為腓腸肌受力曲線。

圖6 腓腸肌受力大小Fig.6 The size of the gastrocnemius muscle

2.2 肌肉靜態(tài)激活度分析

人體在運(yùn)動(dòng)狀態(tài)下，調(diào)動(dòng)休眠的肌肉參與活動(dòng)，這些被調(diào)動(dòng)的肌肉稱為激活。參與活動(dòng)的肌肉群越多，其力量和感覺等方面比較高，激活度就越高，反之則激活度越低。靜態(tài)優(yōu)化的產(chǎn)生需要輸入運(yùn)動(dòng)學(xué)中加速度的肌肉力，如果這些加速度不穩(wěn)定，則肌肉激活和力量會(huì)產(chǎn)生波動(dòng)。研究?jī)煞N改進(jìn)運(yùn)動(dòng)數(shù)據(jù)的方法：運(yùn)動(dòng)學(xué)反解(IK)和使用殘差縮減算法(RRA)來平滑運(yùn)動(dòng)曲線。

過濾運(yùn)動(dòng)學(xué)將對(duì)運(yùn)動(dòng)數(shù)據(jù)進(jìn)行樣條擬合并過濾坐標(biāo)位置，從而使加速度趨于穩(wěn)定。如果使用IK運(yùn)算結(jié)果作為靜態(tài)優(yōu)化的輸入，則應(yīng)始終在外部使用MATLAB或Python或使用OpenSim過濾器進(jìn)行過濾，本次分析使用OpenSim進(jìn)行過濾。

腓腸肌過濾加速度使用運(yùn)動(dòng)學(xué)反解和殘差縮減算法，如圖7所示。在人體步行時(shí)，IK從0.85 s開始激活度從0.01增大到0.37，總時(shí)長(zhǎng)用了0.35 s；RRA從0.9 s開始，激活度從0.01增大到0.34，總時(shí)長(zhǎng)相對(duì)IK時(shí)間縮短0.05 s。即利用RRA方法更能模擬人體肌肉激活度，又由于RRA使用逆動(dòng)力學(xué)模擬，因此輸出運(yùn)動(dòng)學(xué)將是一致的加速度，非常適合在靜態(tài)優(yōu)化中使用。

圖7 腓腸肌不同優(yōu)化方法激活度Fig.7 Different optimization methods for the gastrocnemius

圖8 腓腸肌Fx，F(xiàn)y方向上殘差Fig.8 The residual of gastrocnemius in Fx, Fy direction

從圖8中可以看出，利用IK法在x和y分量上的殘差區(qū)間分別為[-44.9,22.0]、[-55.6,53.5]；利用RRA法在x和y分量上的殘差區(qū)間分別為[-0.6,8.1]、[-14.0,21.5]；即在x分量上RRA比IK的波動(dòng)范圍縮小了近9倍，在y分量縮短了近3倍，且曲線更加平滑。通過對(duì)比，可以明顯地看出利用RRA法比IK法提高了靜態(tài)分析準(zhǔn)確度。

2.3 調(diào)整計(jì)算控制分析

優(yōu)化器將優(yōu)先使用肌肉來產(chǎn)生關(guān)節(jié)扭矩并且僅在肌肉不能產(chǎn)生足夠扭矩時(shí)才使用備用執(zhí)行器。降低肌肉控制計(jì)算(CMC)的強(qiáng)度，并觀察備用執(zhí)行器扭矩和腓腸肌激活度和力的最終變化。對(duì)于2.2節(jié)分析過程中使用備用執(zhí)行器為100 N的最佳力，并且對(duì)于每個(gè)扭矩執(zhí)行器具有100 N·m的最佳扭矩。利用Notepad++將編輯執(zhí)行器文件的每個(gè)執(zhí)行器的最佳力從100 N更改為1 N。

從圖9中可以看出腓腸肌在最佳執(zhí)行力為100 N時(shí)，激活都峰值為0.34，此時(shí)模擬腓腸肌未被完全激活進(jìn)行分析，備用執(zhí)行器被調(diào)用，導(dǎo)致分析結(jié)果可能出現(xiàn)誤差；將最佳執(zhí)行力降到1 N時(shí)，此時(shí)腓腸肌激活度峰值為0.82，但相對(duì)于肌肉而言，成本較高。對(duì)比圖10和圖8(b)，可以看出最佳執(zhí)行力的改變對(duì)于x和y分量上的殘差影響并不明顯。

圖9 腓腸肌不同優(yōu)化強(qiáng)度激活度Fig.9 Different optimal intensity activation of gastrocnemius

圖10 CMC最小執(zhí)行力Fig.10 CMC minimum execution

3 結(jié)論

利用OpenSim分析人體在行走過程中腿部肌肉靜態(tài)力學(xué)分析，以腓腸肌為例。通過對(duì)比IK和RRA分析方法得出如下結(jié)論。

(1)RRA法提高了靜態(tài)分析準(zhǔn)確度，適用于之后人體下肢康復(fù)機(jī)器人的設(shè)計(jì)分析。同時(shí)分析出激活度a的變化曲線，防止康復(fù)機(jī)器人設(shè)計(jì)激發(fā)激活度a變化過快，導(dǎo)致人體肌肉痙攣。

(2)降低外設(shè)備的輔助制動(dòng)力，最大化提升自身肌肉的激活度，可以為進(jìn)一步揭示膝關(guān)節(jié)和踝關(guān)節(jié)的生物力學(xué)特性和人體下肢康復(fù)機(jī)器人的設(shè)計(jì)奠定基礎(chǔ)。