趙 闖, 蔡玉強(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)建簡化人體下肢力學(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)過程中,肌肉長度變化、肌肉長度變化速率、肌肉力臂等肌肉功能參數(shù)以及肌肉功能的動(dòng)態(tài)變化。
總之,中外有不少學(xué)者對人體下肢關(guān)節(jié)運(yùn)動(dòng)的形態(tài)特征進(jìn)行研究和分析,但對于下肢運(yùn)動(dòng)狀態(tài)下的靜態(tài)力學(xué)分析等方面研究不足。本文首先分析人體下肢受力機(jī)理,簡述仿真所需骨骼肌肉模型,利用人工合成的運(yùn)動(dòng)學(xué)方法依賴于生物力學(xué)知識,并結(jié)合正向和逆向運(yùn)動(dòng)學(xué)知識[9],使用OpenSim的Static Optimization Tool模塊,簡化腰部肌肉模型為扭矩執(zhí)行器,改善肌肉激活度和靜態(tài)優(yōu)化輸出的力,利用逆向運(yùn)動(dòng)學(xué)(IK)、殘差縮減(RRA)和肌肉計(jì)算控制(CMC)對人體下肢肌肉進(jìn)行靜態(tài)優(yōu)化分析,為人體下肢康復(fù)機(jī)器人的研究與研發(fā)提供基礎(chǔ)。
人體的膝關(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)作用,維持人體平衡。
將人體的下肢運(yùn)動(dòng)受力模型簡化基于圖2骨骼-肌肉模型簡圖,該模型適用于腿部骨骼-肌肉近似的平面運(yùn)動(dòng),上肢同時(shí)保持近似垂直的狀態(tài)。其中人體上半身簡化為固定在臀部的質(zhì)量塊;大腿、小腿和腳簡化為3個(gè)連桿,髖關(guān)節(jié)、膝關(guān)節(jié)和踝關(guān)節(jié)簡化為鉸接點(diǎn);右腿肌肉簡化為2對作用在連桿兩端的相對作用力,分布于大腿和小腿兩側(cè)。
圖2 骨骼-肌肉右腿模型簡圖Fig.2 Skeleton-muscle right leg model diagram
由于肌肉運(yùn)動(dòng)的非線性,希爾模型[12-13]包含主動(dòng)收縮單元(CE)和被動(dòng)彈簧單元(SE)表示的肌肉模型,用非線性彈簧表示肌腱,如圖3所示。lMT為肌肉-肌腱長度,lT為肌腱長度,lM為標(biāo)準(zhǔn)化的肌纖維長度,α為羽狀角,它們之間存在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)化分析,并在簡化的示例模型中評估腓腸肌激活度和力量產(chǎn)生。靜態(tài)優(yōu)化結(jié)果的質(zhì)量在很大程度上取決于輸入:模型、運(yùn)動(dòng)和力等參數(shù)。靜態(tài)優(yōu)化輸入為肌纖維長度的變化值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為肌肉纖維長度變化系數(shù);l0為靜態(tài)下肌腱長度;Kv為肌肉纖維速度變化系數(shù);v(t) 為肌肉纖維長度變化速率;l(t)為t時(shí)刻肌肉纖維長度變化量;tn為步態(tài)的時(shí)間。
為了建立小腿部腓腸肌的控制方程,對腳部進(jìn)行受力分析,幾何標(biāo)注如圖4所示。
圖4 腳部受力示意圖Fig.4 Schematic diagram of the force on the foot
對于腳部質(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í)間之間的對應(yīng)關(guān)系。
由于難以獲取腳部的生物力學(xué)數(shù)據(jù),在模型中將質(zhì)心作為腳部幾何中心,并將看作剛體,腳部質(zhì)心位置lf由人體腳部大小確定,轉(zhuǎn)動(dòng)慣量Jf、腳部重量Gt和人體重量Gf由人體解剖學(xué)參數(shù)確定。
OpenSim作為建立和分析骨骼-肌肉的開源軟件,可通過軟件建立人體骨骼-肌肉模型,設(shè)置人體參數(shù)來計(jì)算和利用OpenSim的plot功能繪制肌肉力、膝、踝關(guān)節(jié)的運(yùn)動(dòng)位移曲線和力矩等數(shù)據(jù),進(jìn)行人體下肢步行狀態(tài)下肌肉靜態(tài)力學(xué)分析。
研究采用的骨骼-肌肉模型如圖5所示,具有10個(gè)自由度和18個(gè)肌肉的平面縮小步態(tài)模型,模擬正常人體步行運(yùn)動(dòng)狀態(tài)。由于研究不涉及上肢運(yùn)動(dòng),經(jīng)簡化保留軀干部分,并在分析過程中簡化為質(zhì)量塊處理。模型下肢經(jīng)簡化后保留支撐人體站立的脛骨前肌、腓腸肌和比目魚肌[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í)行器,測量分析程中運(yùn)動(dòng)和力之間的關(guān)系,確保整個(gè)分析滿足牛頓第二定律。模擬運(yùn)動(dòng)包含平滑、真實(shí)的加速度,并能精確地測量運(yùn)動(dòng)期間的肌肉的激活度與力矩。圖6為腓腸肌受力曲線。
圖6 腓腸肌受力大小Fig.6 The size of the gastrocnemius muscle
人體在運(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)運(yùn)動(dòng)數(shù)據(jù)的方法:運(yùn)動(dòng)學(xué)反解(IK)和使用殘差縮減算法(RRA)來平滑運(yùn)動(dòng)曲線。
過濾運(yùn)動(dòng)學(xué)將對運(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í)長用了0.35 s;RRA從0.9 s開始,激活度從0.01增大到0.34,總時(shí)長相對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倍,且曲線更加平滑。通過對比,可以明顯地看出利用RRA法比IK法提高了靜態(tài)分析準(zhǔn)確度。
優(yōu)化器將優(yōu)先使用肌肉來產(chǎn)生關(guān)節(jié)扭矩并且僅在肌肉不能產(chǎn)生足夠扭矩時(shí)才使用備用執(zhí)行器。降低肌肉控制計(jì)算(CMC)的強(qiáng)度,并觀察備用執(zhí)行器扭矩和腓腸肌激活度和力的最終變化。對于2.2節(jié)分析過程中使用備用執(zhí)行器為100 N的最佳力,并且對于每個(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,但相對于肌肉而言,成本較高。對比圖10和圖8(b),可以看出最佳執(zhí)行力的改變對于x和y分量上的殘差影響并不明顯。
圖9 腓腸肌不同優(yōu)化強(qiáng)度激活度Fig.9 Different optimal intensity activation of gastrocnemius
圖10 CMC最小執(zhí)行力Fig.10 CMC minimum execution
利用OpenSim分析人體在行走過程中腿部肌肉靜態(tài)力學(xué)分析,以腓腸肌為例。通過對比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ǔ)。