四連桿膝關(guān)節(jié)假肢的動力學(xué)建模與分析1)

呂 陽 方虹斌 徐 鑒 馬建敏 王啟寧 張曉旭

(復(fù)旦大學(xué)航空航天系，上海 200433)

(復(fù)旦大學(xué)智能機(jī)器人研究院，上海 200433)

(上海智能機(jī)器人工程技術(shù)研究中心，復(fù)旦大學(xué)，上海 200433)

(智能機(jī)器人教育部工程研究中心，復(fù)旦大學(xué)，上海 200433)

(北京大學(xué)工學(xué)院，北京 100871)

引言

根據(jù)全國第二次殘疾人抽樣調(diào)查的結(jié)果[1]，我國肢體殘疾2412 萬人，占?xì)埣踩丝倲?shù)的29.07%，其中下肢截肢患者大概有158 萬人.為了使大腿截肢患者能夠正常地生活和工作，越來越多的研究者開始致力于假肢的結(jié)構(gòu)設(shè)計與優(yōu)化研究.

膝關(guān)節(jié)假肢在結(jié)構(gòu)上可分為單軸式膝關(guān)節(jié)和多軸式膝關(guān)節(jié).單軸式膝關(guān)節(jié)具有結(jié)構(gòu)簡單、可靠性高的優(yōu)點，其代表性產(chǎn)品主要有奧托博克公司的CLeg、奧索公司的Rheo Knee 等[2].多軸式膝關(guān)節(jié)如四連桿膝關(guān)節(jié)則具有良好的仿生特性[3]，能夠較為真實地逼近人體膝關(guān)節(jié)轉(zhuǎn)動瞬心(小腿相對大腿轉(zhuǎn)動時的速度瞬心，簡稱轉(zhuǎn)動瞬心)的J 形軌跡，其代表性產(chǎn)品主要有奧索公司的Parso Knee、北京工道風(fēng)行公司的K301 等.值得一提的是，近年來，折紙結(jié)構(gòu)由于其結(jié)構(gòu)輕、可編程以及超材料特性被廣泛研究，這些特性在膝關(guān)節(jié)的J 形軌跡實現(xiàn)上也有一定的應(yīng)用潛力[4].在運動性能上，J 形軌跡意味著小腿假肢在站立時的瞬心較高，從而實現(xiàn)較好的穩(wěn)定性；在擺動時的瞬心較低，有較好的靈活性[5].因此，相比于單軸式膝關(guān)節(jié)，四連桿膝關(guān)節(jié)兼顧運動穩(wěn)定性，具有更好的運動性能.

為了對四連桿膝關(guān)節(jié)的運動性能有更好的認(rèn)識，需要明確影響其運動性能的因素.從動力學(xué)的角度來看，影響運動性能的因素包括足--地交互作用力、膝關(guān)節(jié)單邊約束產(chǎn)生的強(qiáng)非線性作用力、液壓阻尼器提供的彈性回復(fù)力及阻尼力等.要明確各因素對假肢性能的影響機(jī)制，就必須建立考慮足--地交互作用力和膝關(guān)節(jié)單邊接觸力的一體化動力學(xué)模型，并就該模型產(chǎn)生的動力學(xué)響應(yīng)進(jìn)行分析.

目前，對于假肢的建模，常用的建模方法是利用第一類拉格朗日方程.這種方法將假肢運動階段分為站立相(假肢與地面接觸，處于站立狀態(tài))和擺動相(假肢未與地面接觸，處于擺動狀態(tài))，對不同的步態(tài)階段分別進(jìn)行建模[6-10].在站立相階段，該方法將足--地接觸簡化為理想鉸接約束，通過在動力學(xué)方程中引入拉格朗日乘子求解模型；在膝關(guān)節(jié)存在單邊約束時，假設(shè)大腿和小腿為一根剛性桿，引入拉格朗日乘子對模型進(jìn)行求解.這種建模方法未對足--地交互作用力以及膝關(guān)節(jié)單邊接觸力直接進(jìn)行建模，而是基于假設(shè)將這兩種力看作理想約束力.實際上足--地交互作用力以及膝關(guān)節(jié)單邊接觸力等強(qiáng)非線性因素為非理想約束力，上述處理方式無法體現(xiàn)這一點，因此需要運用更合理的方式對這兩種力進(jìn)行建模.

足--地交互作用力本質(zhì)上是足--地接觸時產(chǎn)生的接觸力和摩擦力，膝關(guān)節(jié)單邊接觸力本質(zhì)上是膝關(guān)節(jié)連桿間限位碰撞產(chǎn)生的接觸力，因而其建模方式可以借鑒多剛體動力學(xué)對接觸力和摩擦力的描述.接觸力模型中有赫茲接觸力模型[11]以及改進(jìn)的赫茲接觸力模型[12].赫茲接觸力模型形式相對簡單，但不考慮能量耗散，不適用于真實的接觸情況；改進(jìn)的赫茲接觸力模型通過各種方式加入了能量耗散項，常用的有Hunt-Crossley 模型[13]、Gonthier 模型[14]、Kelvin-Voigt[15]模型等，更適用于真實的接觸情況.摩擦力模型主要包括：庫倫摩擦模型[16]、Stribeck 摩擦模型[17-19]以及更為復(fù)雜的LuGre 摩擦模型[20-23].實驗研究[24-27]表明，當(dāng)剛體相對運動尺度為毫米量級或更高時，若要分析摩擦對系統(tǒng)動力學(xué)行為的影響，采用庫倫摩擦模型即可達(dá)到足夠高的精度.為解決庫倫摩擦模型中非光滑本構(gòu)對動力學(xué)分析帶來的困難，還可進(jìn)一步利用非線性函數(shù)[28]對庫倫摩擦模型進(jìn)行光滑化.近年來，基于上述接觸和摩擦模型的雙足行走機(jī)器人[29-32]動力學(xué)分析已經(jīng)取得一定進(jìn)展，但基于上述接觸和摩擦模型的假肢動力學(xué)分析仍處于起步階段.

綜上所述，足--地交互作用力以及膝關(guān)節(jié)單邊接觸力等強(qiáng)非線性因素對假肢步態(tài)的影響尚未得到深入分析，其難點在于考慮足--地交互作用力和膝關(guān)節(jié)單邊接觸力的一體化動力學(xué)建模.針對該問題，基于多剛體動力學(xué)中關(guān)于接觸力及摩擦力的研究，本文將采用Kelvin-Voigt 模型和庫倫模型描述足--地接觸力和摩擦力，采用Kelvin-Voigt 模型描述膝關(guān)節(jié)單邊接觸力，從而基于第一類拉格朗日方程對四連桿膝關(guān)節(jié)被動假肢進(jìn)行動力學(xué)建模.將步態(tài)實驗測得的髖關(guān)節(jié)運動數(shù)據(jù)作為動力學(xué)模型的驅(qū)動信號，四連桿膝關(guān)節(jié)被動假肢的步態(tài)可以通過數(shù)值計算得到.進(jìn)一步，為了評價假肢步態(tài)與健康人實測步態(tài)的一致性，本文定義了步態(tài)相關(guān)系數(shù)作為評價指標(biāo)，并分析四連桿膝關(guān)節(jié)液壓阻尼器剛度、阻尼參數(shù)對相關(guān)系數(shù)的影響.

本文的貢獻(xiàn)和創(chuàng)新點在于：首先，進(jìn)一步完善四連桿膝關(guān)節(jié)假肢的動力學(xué)模型，重點關(guān)注足--地交互作用力以及膝關(guān)節(jié)單邊接觸力等強(qiáng)非線性因素對下肢假肢步態(tài)的影響，為殘疾人行走時的代償行為提供一種新的力學(xué)解釋；其次，以實測數(shù)據(jù)為參考信號進(jìn)行步態(tài)一致性評估和參數(shù)分析，為四連桿膝關(guān)節(jié)被動假肢的進(jìn)一步優(yōu)化提供理論支撐，提高結(jié)論的可靠性.

本文第1 節(jié)介紹四連桿膝關(guān)節(jié)假肢動力學(xué)模型的建立.首先分析四連桿的運動學(xué)特性，接著引入足--地交互作用力和膝關(guān)節(jié)單邊接觸力，最后基于第一類拉格朗日方程對假肢進(jìn)行建模.第2 節(jié)介紹健康被試的步態(tài)測試.一方面，實測髖關(guān)節(jié)運動數(shù)據(jù)作為假肢動力學(xué)模型的驅(qū)動信號，用于假肢步態(tài)分析；另一方面，實測膝關(guān)節(jié)運動數(shù)據(jù)作為假肢步態(tài)運動的對比信號，用于步態(tài)一致性評估及參數(shù)分析.第3 節(jié)為基于假肢動力學(xué)模型的步態(tài)分析.將數(shù)值計算所得結(jié)果與實驗測得的步態(tài)進(jìn)行比較，討論四連桿液壓阻尼器剛度、阻尼參數(shù)對步態(tài)周期失諧及步態(tài)一致性的影響.第4 節(jié)給出研究結(jié)論.本文通過數(shù)值計算研究了四連桿膝關(guān)節(jié)被動假肢動力學(xué)模型的響應(yīng)，分析了四連桿膝關(guān)節(jié)液壓阻尼器的阻尼、剛度參數(shù)對假肢步態(tài)指標(biāo)的影響，總結(jié)了提高步態(tài)性能的一般經(jīng)驗，以期為四連桿膝關(guān)節(jié)被動假肢的參數(shù)設(shè)計提供參考依據(jù).

1 四連桿膝關(guān)節(jié)被動假肢動力學(xué)模型

四連桿膝關(guān)節(jié)被動假肢具有良好的仿生特性，但由于結(jié)構(gòu)更為復(fù)雜，足--地交互作用力、膝關(guān)節(jié)單邊接觸力等強(qiáng)非線性因素會導(dǎo)致下肢假肢產(chǎn)生更為豐富的動態(tài)響應(yīng)，進(jìn)而在整體上影響步態(tài)的周期性和穿戴舒適性.因此，本文有必要從動力學(xué)角度對四連桿膝關(guān)節(jié)假肢進(jìn)行分析.

1.1 四連桿機(jī)構(gòu)的運動學(xué)模型

如圖1 所示，人體膝關(guān)節(jié)由股骨踝(femoral ankle)和脛骨踝(tibia ankle)組成，依賴于股踝和脛踝接觸面的幾何外形，人體膝關(guān)節(jié)的運動屬于非定軸轉(zhuǎn)動，小腿相對于大腿轉(zhuǎn)動的速度瞬心軌跡實際是一條J 形曲線.其生物學(xué)優(yōu)勢在于，在站立相時，膝關(guān)節(jié)瞬心升高，小腿相對于瞬心的轉(zhuǎn)動慣量較大，不易實現(xiàn)轉(zhuǎn)動，因而穩(wěn)定性好；反之，在擺動相時，小腿靈活性好，這也是四連桿膝關(guān)節(jié)仿生學(xué)研究的主要動機(jī).

圖1 人體膝關(guān)節(jié)非定軸轉(zhuǎn)動及其J 形瞬心軌跡Fig.1 Non-fixed axis rotation of human knee joint and its J-curve of instantaneous rotation center

四連桿膝關(guān)節(jié)可以較好地模擬人體膝關(guān)節(jié)非定軸轉(zhuǎn)動過程，主要由四根連桿和液壓阻尼器組成，其結(jié)構(gòu)示意圖及實物圖如圖2 所示.以工道風(fēng)行公司(GDFX)設(shè)計的一款四連桿膝關(guān)節(jié)為例，其幾何參數(shù)及慣量參數(shù)如表1 所示.根據(jù)四連桿的連接關(guān)系，四連桿機(jī)構(gòu)的運動約束方程可以表示為

在二維笛卡爾坐標(biāo)系下，式(1)可以具體表示為

圖2 四連桿膝關(guān)節(jié)Fig.2 Four-bar linkage prosthesis knee

表1 四連桿膝關(guān)節(jié)物理參數(shù)表Table 1 Physical parameter values of the four-bar-linkage knee prosthesis

其中，lAB，lAD，lBC和lCD分別為連桿AB，AD，BC和CD的長度，γ1=90?為連桿AB與大腿股骨夾角，γ2=24?為連桿CD與小腿脛骨夾角，θA，θB分別為連桿AD，BC與豎直方向夾角，α1為大腿股骨與豎直方向夾角，α11為小腿股骨與豎直方向夾角.

由幾何關(guān)系可知，四連桿膝關(guān)節(jié)在轉(zhuǎn)動過程中小腿相對于大腿的瞬心P在CB和DA的延長線上.在平面笛卡爾坐標(biāo)系下，瞬心P的坐標(biāo)可以表示為

其中，xP為瞬心P的橫坐標(biāo)，yP為瞬心P的縱坐標(biāo)，xk和yk分別為四連桿頂點k(k=A,B,C,D)的橫坐標(biāo)和縱坐標(biāo).圖3 展示了小腿相對于大腿轉(zhuǎn)動時速度瞬心P的軌跡曲線，其中假肢大腿與AB桿剛性連接，假肢小腿與CD桿剛性連接.可以看到，在膝關(guān)節(jié)由伸直狀態(tài)至110?屈曲狀態(tài)過程中，瞬心P形成的軌跡是一條J-curve 線，與人體膝關(guān)節(jié)瞬心的J 形軌跡相似，展現(xiàn)了良好的仿生特性.

1.2 足--地交互力和膝關(guān)節(jié)單邊接觸力模型

足--地交互作用力主要有地面對假肢足底的接觸力和摩擦力兩部分.為方便分析，本文假設(shè)只有腳跟和腳尖這兩點與地面存在接觸力和摩擦力.

地面對足底的接觸力可以通過Kelvin-Voigt 模型來描述

圖3 四連桿膝關(guān)節(jié)J 形瞬心軌跡Fig.3 Four-bar linkage prosthesis knee’s J-curve of instantaneous rotation center

式中，F(xiàn)N,i表示接觸力，方向垂直于地面，下標(biāo)i為heel 或toe，分別表示腳跟和腳尖與地面接觸時接觸力的描述.KN,i表示接觸剛度，χN,i表示接觸阻尼.δi表示腳尖或腳跟與地面之間的壓痕深度，當(dāng)腳尖或腳跟離地時，規(guī)定δi=0.地面對足底的摩擦力采用庫倫摩擦模型描述[16]

式中，F(xiàn)f,i表示摩擦力，方向沿水平方向，下標(biāo)i為heel 或toe，分別表示腳跟和腳尖與地面接觸時摩擦力的描述.μ表示動摩擦系數(shù)(此處忽略靜摩擦系數(shù))，vi表示腳尖或腳跟相對地面的切向速度.為方便動力學(xué)分析，本文對式(4)所述接觸力模型和式(5)所述摩擦力進(jìn)行光滑化處理，具體光滑化形式分別為

其中C1和C2為常數(shù)，表征接觸力和摩擦力的光滑化程度.

此外，如圖4 所示，當(dāng)假肢小腿擺動到與大腿呈同一角度時(即膝關(guān)節(jié)打直)，為保證穿戴者安全性，四連桿的AB和BC兩桿之間將通過限位裝置鎖死，鎖死端面垂直于AB連桿.這種結(jié)構(gòu)鎖死是一種單邊約束，力學(xué)上可以簡化為單邊接觸模型.因此，本文同樣通過光滑化的Kelvin–Voigt 模型來描述膝關(guān)節(jié)打直的單邊接觸力，具體形式為

其中，F(xiàn)c表示單邊接觸力，方向為AB連桿方向，Kc表示接觸剛度，Cc表示接觸阻尼，δknee表示兩桿之間的壓痕深度.需要注意的是，如圖4 中細(xì)節(jié)放大圖所示，AB和BC兩桿間限位接觸點與B點不重合，其偏心距為d=1 cm，因此式(8)所示接觸力作用于連桿BC和AB的力矩分別為

至此，本文得到了足--地交互作用力和膝關(guān)節(jié)單邊接觸力的模型.

圖4 膝關(guān)節(jié)幾何鎖死時接觸力示意圖Fig.4 The contact force when the knee is locked geometrically

1.3 基于第一類拉格朗日方程的動力學(xué)建模

假設(shè)四連桿膝關(guān)節(jié)的接受腔與殘疾人大腿殘肢為剛性連接，因此四連桿膝關(guān)節(jié)被動假肢的運動模型可簡化為圖5 所示形式.其中x和y分別為髖關(guān)節(jié)在平面笛卡爾坐標(biāo)系下的橫坐標(biāo)和縱坐標(biāo)，α1為大腿與豎直方向的夾角，α11為小腿與豎直方向夾角，β1為足部與水平方向夾角.

如圖5 所示，四連桿機(jī)構(gòu)中的液壓阻尼器主要是在運動過程中提供阻尼力和回復(fù)力，其合力可以表示為

式中，k和c分別為液壓阻尼器的剛度和阻尼系數(shù)，l0=168.79 mm 為液壓阻尼器原長，l為當(dāng)前長度.基于虛功原理，合力F0作用于廣義坐標(biāo)θ1上的廣義力為

圖5 四連桿膝關(guān)節(jié)被動假肢運動模型Fig.5 Model of passive four-bar linkage prosthesis knee

此外，由于θ1=0.5π ?α1+θA，因此合力F0作用于廣義坐標(biāo)α1和θA上的廣義力矩分別為?τ1和τ1.

此外，本文假設(shè)踝關(guān)節(jié)受到帶阻尼的扭簧作用，作用力矩表示為

其中，kankle=10 N/m 和cankle=0.2N/(m·s?1)分別為扭簧的剛度和阻尼系數(shù)，為踝關(guān)節(jié)角度，為踝關(guān)節(jié)期望角度.在步態(tài)分析中，本文以實測踝關(guān)節(jié)角度作為.由圖5 可知，踝關(guān)節(jié)角度滿足=β1?α11，因此扭簧作用于廣義坐標(biāo)α11和β1的扭矩分別為?Tankle和Tankle.

得到基于四連桿約束的運動學(xué)模型、足--地交互力和膝關(guān)節(jié)單邊接觸力模型后，圖5 所示的四連桿膝關(guān)節(jié)被動假肢的動力學(xué)模型可以用第一類拉格朗日方程表示

式中，q=(x,y,α1,α11,θA,θB,β1)T為廣義坐標(biāo)向量，是質(zhì)量陣，是科氏力和離心力項，是重力項.M(q),及N(q)中各元素表達(dá)式見附錄A.Fe為足--地交互作用力和膝關(guān)節(jié)單邊接觸力構(gòu)成的廣義力向量，具體表示為

其中，F(xiàn)N,heel和FN,toe由式(6)給出，F(xiàn)f,heel和Ff,toe由式(7)給出，Tc和由式(9)給出，J為將足--地交互作用力和膝關(guān)節(jié)單邊接觸力變換為廣義力的Jacobi 矩陣，具體形式見附錄A.Fq為膝關(guān)節(jié)的液壓阻尼器和踝關(guān)節(jié)的阻尼彈簧施加的廣義力向量，具體表示為

其中，τ1由式(11)給出，Tankle由式(12)給出.Φ 為四連桿約束矩陣，具體表示為

其中，ξ1和ξ2為式(2)所示約束方程.λ 為拉格朗日乘子，表示四連桿約束產(chǎn)生的約束力.

需要指出的是，被動假肢是由髖關(guān)節(jié)運動驅(qū)動的，髖關(guān)節(jié)運動信號x，y，α1為假肢系統(tǒng)的輸入，α11，θA，θB和β1為假肢系統(tǒng)的響應(yīng).為方便步態(tài)分析，式(13)可以改寫為

將式(18)對時間t求導(dǎo)可得

代入式(17)可得

進(jìn)一步，將式(20)代入式(19)可解得

式(17)可以改寫為如下狀態(tài)方程形式

至此，我們得到了考慮四連桿約束、足--地交互力和膝關(guān)節(jié)單邊接觸力的四連桿膝關(guān)節(jié)被動假肢的動力學(xué)模型.如前所述，被動假肢是由髖關(guān)節(jié)運動驅(qū)動的，髖關(guān)節(jié)運動信號q1=(x,y,α1)T及其速度和加速度為式(22)的輸入.后續(xù)，在基于該模型的動力學(xué)分析中，q1為健康被試的實測運動信號.

2 健康被試的步態(tài)測試

該人體步態(tài)測試實驗已獲得北京大學(xué)倫理委員會批準(zhǔn).實驗?zāi)康脑谟冢旱谝?，獲得式(22)右端所需要的輸入q1，，用于四連桿膝關(guān)節(jié)被動假肢的動力學(xué)仿真；第二，將健康被試的步態(tài)信號作為健肢步態(tài)信號，與動力學(xué)仿真所得假肢步態(tài)進(jìn)行比較，用于假肢運動與健肢運動的步態(tài)分析.

實驗所用設(shè)備為Cosmos Gaitway 跑臺和Vicon光學(xué)運動捕捉系統(tǒng)及配套分析系統(tǒng).其中，Cosmos Gaitway 跑臺用于實現(xiàn)勻速步行環(huán)境，Vicon 光學(xué)運動捕捉系統(tǒng)用于捕捉人體靶點運動信號，并換算為各關(guān)節(jié)角度信息.參考步態(tài)分析軟件OrthoTrak 用戶指南中的Helen-Hayes 模型，健康被試身共貼有18 個靶點，分別為右肩靶點(R.shoulder)、左肩靶點(L.shoulder)、補(bǔ)償靶點(offset)、左髂前上棘靶點(R.asis)、右髂前上棘靶點(L.asis)、骶骨靶點(V.sacral)、右大腿靶點(R.thigh)、左大腿靶點(L.thigh)、右膝靶點(R.knee)、左膝靶點(L.knee)、右小腿靶點(R.shank)、左小腿靶點(L.shank)、右腳跟靶點(R.heel)、左腳跟靶點(L.heel)、右腳踝靶點(R.ankle)、左腳踝靶點(L.ankle)、右腳尖靶點(R.toe)和左腳尖靶點(L.toe).步態(tài)測試實景如圖6 所示，圖中Y軸沿著被試行走正前方，X軸沿著被試行走側(cè)向，Z軸沿著豎直向上方向.實驗過程中，健康被試(男性，身高182 cm，體重92 kg)在Cosmos Gaitway 跑臺上行走(坡度0，速度1 m/s)，Vicon 光學(xué)運動捕捉系統(tǒng)采樣頻率為100 Hz，采樣時長180 s.

圖6 健康被試步態(tài)測試過程Fig.6 Able-bodied subject test process

實驗最終得到的數(shù)據(jù)是各靶點的空間坐標(biāo)，采用Helen-Hayes 模型的角度解算方法可以得到髖關(guān)節(jié)處的空間坐標(biāo)和髖、膝、踝關(guān)節(jié)的空間轉(zhuǎn)角[33].式(22)所示動力學(xué)模型為平面模型，根據(jù)坐標(biāo)對應(yīng)關(guān)系，本文取圖6 中髖關(guān)節(jié)處的Y軸坐標(biāo)和Z軸坐標(biāo)分別作為廣義坐標(biāo)q中x，y坐標(biāo)，取各關(guān)節(jié)繞X軸方向的轉(zhuǎn)角作為廣義坐標(biāo)q中的α1，α11，β1坐標(biāo).

為使實驗數(shù)據(jù)反映健康人步態(tài)的普遍規(guī)律，本文以兩次腳跟著地(heel strike，HS)之間的時間間隔作為一個步態(tài)周期T，選取實驗數(shù)據(jù)中的30 個周期進(jìn)行平均，同時將時間歸一化.處理后得到的q1數(shù)據(jù)如圖7 所示.圖7(a)中的藍(lán)色實線表示平均后的q1時程，陰影部分的上下界分別表示30 個周期中各個時刻點的最大值和最小值.從陰影部分的面積可以看出，各周期步態(tài)較為一致，實驗結(jié)果較為可靠，因此取30 個周期進(jìn)行平均是合理的，能夠準(zhǔn)確地反映q1廣義坐標(biāo)的運動特征.圖7(b)中展示了廣義坐標(biāo)y和α1的頻譜圖(廣義坐標(biāo)x的運動不具有周期性，因此無需進(jìn)行頻譜分析).可以看到，y信號中2 Hz 頻率成分的振幅較大，這是因為健康被試在一個步態(tài)周期中，左腿和右腿各走了一步，導(dǎo)致髖關(guān)節(jié)在豎直方向呈現(xiàn)兩次近似相同的起伏；α1信號中1 Hz頻率成分的振幅較大，這一點并不難理解.上述分析可以說明輸入信號q1的激勵頻率為1 Hz，這為第3 中對系統(tǒng)響應(yīng)的分析奠定基礎(chǔ).

圖7 髖關(guān)節(jié)實測運動數(shù)據(jù)時程和頻譜曲線Fig.7 Time history and response spectrum of actually measured motion of hip joint

3 四連桿被動假肢的步態(tài)分析

基于式(22)并以第2 節(jié)的實驗數(shù)據(jù)作為輸入，對四連桿膝關(guān)節(jié)假肢進(jìn)行動力學(xué)仿真，仿真所采用的物理參數(shù)如表2 所示.足--地交互作用力和膝關(guān)節(jié)單邊接觸力的模型參數(shù)如下

為保證仿真所得數(shù)據(jù)為假肢的穩(wěn)態(tài)響應(yīng)，仿真時長設(shè)置為20 個步態(tài)周期，取末5 個步態(tài)周期用于分析.

表2 仿真中的部分物理參數(shù)值Table 2 Physical parameters in numerical simulation

3.1 亞諧波頻率響應(yīng)

動力學(xué)仿真中，設(shè)置液壓阻尼器阻尼/剛度比為c/k=0.05，剛度系數(shù)分別為k=1.4 × 104N/m和k=1.52×104N/m.圖8(a)展示了健康被測的實測步態(tài)數(shù)據(jù)及基于動力學(xué)模型的仿真數(shù)據(jù)，其中縱坐標(biāo)定義為=α1?α11.則從圖8(a)陰影部分可以看出，當(dāng)k=1.4×104N/m 時，四連桿膝關(guān)節(jié)被動假肢的動力學(xué)仿真響應(yīng)的周期為2T；當(dāng)k=1.52×104N/m，動力學(xué)仿真響應(yīng)的周期為3T.

圖8 亞諧波響應(yīng)時程曲線與頻譜曲線Fig.8 Time history and response spectrum of subharmonic response

進(jìn)一步地，我們對動力學(xué)響應(yīng)信號進(jìn)行FFT頻譜分析.圖8(b)展示了對應(yīng)于兩個時間歷程的頻譜圖.圖8(b)表明，當(dāng)k=1.4 × 104N/m 時，主要的頻率成分為1 Hz 和2 Hz，同時出現(xiàn)了1/2 Hz 頻率成分，這一亞諧頻率成分使得響應(yīng)的周期為2T；當(dāng)k=1.52×104N/m 時，主要頻率成分仍為1 Hz 和2 Hz，同時出現(xiàn)了1/3 Hz 頻率成分，這一亞諧頻率成分使得響應(yīng)的周期為3T.這些亞諧波頻率成分出現(xiàn)的原因在于，足--地交互作用力和膝關(guān)節(jié)單邊接觸力含有強(qiáng)非線性因素，可以誘發(fā)四連桿被動假肢的亞諧波響應(yīng)這一強(qiáng)非線性動力學(xué)行為.1/2 或1/3 亞諧波響應(yīng)表明，此時假肢的步態(tài)周期為健肢步態(tài)周期的2 倍或3 倍，破壞了行走時雙腿步態(tài)的協(xié)調(diào)性，本文將這種亞諧波響應(yīng)稱為步態(tài)周期失諧.

通過進(jìn)一步的非線性動力學(xué)分析，我們發(fā)現(xiàn)，通過調(diào)整激勵參數(shù)(如初始條件、相位、幅值等)，系統(tǒng)的穩(wěn)態(tài)響應(yīng)可以切換回單周期狀態(tài).反映到式(22)所示的四連桿膝關(guān)節(jié)被動假肢模型中，即表現(xiàn)為髖關(guān)節(jié)運動信號的調(diào)整.例如，圖9(a)和圖9(b)分別展示了將髖關(guān)節(jié)豎向運動信號y放大為原信號1.006 倍后假肢膝關(guān)節(jié)角度的時程曲線和頻譜曲線.從圖9(a)可以看出，在增大髖關(guān)節(jié)豎向運動幅值后，假肢響應(yīng)在兩個剛度值下均達(dá)到單周期穩(wěn)態(tài)，周期和實測數(shù)據(jù)相同為T，并與實測數(shù)據(jù)有較好的步態(tài)一致性.圖9(b)表明，假肢的亞諧波頻率成分消失，其基頻與健肢的基頻相同，為1 Hz，步態(tài)周期不再失諧.

圖9 提胯后響應(yīng)時程曲線與頻譜曲線Fig.9 Time history and response spectrum when lifting the hip

圖9 提胯后響應(yīng)時程曲線與頻譜曲線(續(xù))Fig.9 Time history and response spectrum when lifting the hip(continued)

放大髖關(guān)節(jié)豎向運動信號所反映的生物意義類似于提胯，即殘疾人在運動過程中提高胯部動作幅度的行為.圖9 呈現(xiàn)的結(jié)論說明，提胯動作將改善假肢步態(tài)周期的協(xié)調(diào)性，使步態(tài)不再失諧，這在定性上解釋了殘疾人穿戴假肢行走時會出現(xiàn)提胯等代償動作產(chǎn)生的原因.

注意到，在式(22)所示的四連桿膝關(guān)節(jié)被動假肢動力學(xué)模型中，四連桿的幾何參數(shù)、模型所受激勵信號均是基于實驗實測所得，因而所得結(jié)果具有較高的可信度.但上述亞諧波頻率成分和亞諧波響應(yīng)，尚未在四連桿膝關(guān)節(jié)的研究文獻(xiàn)中見到相關(guān)報道.

數(shù)值分析進(jìn)一步發(fā)現(xiàn)，保持四連桿膝關(guān)節(jié)液壓阻尼器阻尼比值不變，在k值達(dá)到1.6×104N/m 以上時，亞諧波頻率響應(yīng)基本消失，假肢膝關(guān)節(jié)的運動恢復(fù)為單周期的穩(wěn)態(tài)解，因此將k=1.6×104N/m 稱為此時的臨界k值.圖10 展示了k=1.6×104N/m 時假肢動力學(xué)仿真時程和響應(yīng)頻譜，顯見假肢響應(yīng)達(dá)到了單周期穩(wěn)態(tài)，且周期與實測數(shù)據(jù)相同.

圖10 單周期穩(wěn)態(tài)下假肢動力學(xué)仿真步態(tài)時程和響應(yīng)頻譜Fig.10 Time history and response spectrum of prosthesis dynamic simulation gait under uniperiodic steady state

圖11 不同剛度值下假肢膝關(guān)節(jié)頻率響應(yīng)圖Fig.11 Response spectrum of prosthesis knee under different values of stiffness

圖11 給出了c/k=0.05 時各個剛度值下膝關(guān)節(jié)的時程響應(yīng)的頻譜圖.可以看到，在剛度值在(1.4～1.65)×104N/m 變化的過程中，響應(yīng)在紅色區(qū)域(1.45～1.47)×104N/m 出現(xiàn)了持續(xù)且顯著的1/2 亞諧波頻率成分，在(1.47～1.5)×104N/m 進(jìn)入多亞諧波頻率成分混雜的階段，在藍(lán)色區(qū)域(1.5～1.56)×104N/m 出現(xiàn)了持續(xù)且顯著的1/3 亞諧波頻率成分，隨后各亞諧波頻率成分逐漸減少并進(jìn)入單一穩(wěn)態(tài)階段(1.56～1.65)×104N/m.另外值得注意的是，在(1.4～1.45)×104N/m 階段，響應(yīng)在(1.4～1.42)×104N/m 時出現(xiàn)了1/2 亞諧波頻率成分，而在(1.42～1.45)×104N/m 時出現(xiàn)多亞諧波頻率成分混雜的情況，這一部分值得在之后的研究中進(jìn)一步探討.根據(jù)以上分析，可以認(rèn)為在剛度值增大到一定程度時，亞諧波頻率成分會消失，將亞諧波頻率成分完全消失時的k值叫作臨界k值，并在表3 中給出了其他c/k值下的臨界k值.

分析表3 可以發(fā)現(xiàn)，阻尼比值越大，單周期穩(wěn)態(tài)解對應(yīng)的臨界k值越低.這表明，可以通過適當(dāng)調(diào)大四連桿膝關(guān)節(jié)液壓阻尼器的阻尼、剛度系數(shù)來改善假肢步態(tài)周期的協(xié)調(diào)性.

表3 臨界k 值Table 3 Critical value of k

3.2 假肢膝關(guān)節(jié)最大屈曲角度

圖12 不同阻尼比值下膝關(guān)節(jié)最大屈曲角Fig.12 Peak of knee angle under different damping ratios

可以看出，隨著k值增大，膝關(guān)節(jié)最大屈曲角度呈現(xiàn)遞減趨勢.同時可以看出，隨著阻尼比值的增大，假肢膝關(guān)節(jié)最大屈曲角與實測膝關(guān)節(jié)最大屈曲角相等時所對應(yīng)的剛度值k(紅色豎直虛線所對應(yīng)的k值)也呈現(xiàn)遞減規(guī)律.這表明，若調(diào)大膝關(guān)節(jié)液壓阻尼器的阻尼、剛度系數(shù)，盡管假肢步態(tài)周期的協(xié)調(diào)性能夠得到改善，但假肢膝關(guān)節(jié)的最大屈曲角卻有可能偏離最優(yōu)值.因此，未來的工作中，我們有必要選擇或定義合理的步態(tài)協(xié)調(diào)性指標(biāo)，以便對四連桿膝關(guān)節(jié)假肢進(jìn)行優(yōu)化.

3.3 假肢與健肢膝關(guān)節(jié)步態(tài)一致性

為了定量衡量假肢與健肢步態(tài)的一致性，本文定義假肢與健肢膝關(guān)節(jié)信號的相關(guān)系數(shù)(correlation coefficient)，計算公式如下

式中X=α1?α11為動力學(xué)仿真所得膝關(guān)節(jié)角度信號，Y為健康被試實測膝關(guān)節(jié)角度信號，Cov(X,Y)為X和Y信號的協(xié)方差，D(X)和D(Y)分別為X和Y信號的方差，t1/T和t2/T分別為15 和20.理論上，，ρXY越大，表示假肢膝關(guān)節(jié)角度信號與健肢膝關(guān)節(jié)角度信號一致性越高.當(dāng)且僅當(dāng)X≡Y時，ρXY=1.

圖13 不同阻尼比值下假肢-健肢膝關(guān)節(jié)相關(guān)系數(shù)Fig.13 Prosthesis-intact knee joint correlation coefficient under different damping ratios

圖13 展示了不同的c/k值下ρXY隨液壓阻尼器剛度的變化.可以看出，在c/k較小(c/k=0.03,0.05)且k較小時，ρXY呈現(xiàn)抖動的狀態(tài)，這是因為，此時假肢的響應(yīng)出現(xiàn)亞諧波頻率成分，導(dǎo)致ρXY出現(xiàn)極小值點，呈現(xiàn)抖動的狀態(tài).根據(jù)3.1 節(jié)的分析可知，k在超過一定值時，假肢響應(yīng)會轉(zhuǎn)變?yōu)閱沃芷诜€(wěn)態(tài)響應(yīng)，圖13 中已根據(jù)表3 標(biāo)出臨界k值(紅色豎直虛線所對應(yīng)的k值).在達(dá)到臨界k值后，ρXY呈現(xiàn)先增大后減小的規(guī)律性變化.這進(jìn)一步驗證了3.1 節(jié)所述臨界k值的合理性.

在轉(zhuǎn)變?yōu)閱沃芷诜€(wěn)態(tài)響應(yīng)之后，在不同的c/k下，隨著k增大，ρXY均呈現(xiàn)先增大后減小的規(guī)律，因此ρXY存在一個最大值點.這表明我們可以通過四連桿膝關(guān)節(jié)液壓阻尼器阻尼、剛度參數(shù)的設(shè)計，對假肢的步態(tài)一致性進(jìn)行優(yōu)化.表4 反映了不同液壓阻尼器阻尼比值設(shè)置下，最大相關(guān)系數(shù)ρXY及其對應(yīng)的液壓阻尼器剛度值.可以看出，在不同的液壓阻尼器阻尼比值下，最大相關(guān)系數(shù)的值均相近，且最大可得到0.989 6，這為設(shè)計液壓阻尼器的剛度和阻尼值提供了一定的理論依據(jù).

表4 最大相關(guān)系數(shù)及對應(yīng)液壓阻尼器剛度Table 4 Max ρ and its corresponding hydraulic damper stiffness

同時可以發(fā)現(xiàn)，隨著c/k增大，ρXY在達(dá)到最大值后的下降速度增大.因此可以發(fā)現(xiàn)，c/k越小，ρXY越能在一個較大的k范圍內(nèi)保持平緩的變化.這表明，為了使假肢與健肢在較大的k值范圍內(nèi)達(dá)到較好的步態(tài)一致性，c/k不宜取得過大.

以上分析表明，為了實現(xiàn)最好的步態(tài)一致性，可以通過單目標(biāo)優(yōu)化方法找到最合適的k值和c值.另外，也可以圍繞著步態(tài)一致性探究一些多目標(biāo)優(yōu)化問題，如為了實現(xiàn)在較大的k值范圍內(nèi)達(dá)到較好的步態(tài)一致性，通過多目標(biāo)優(yōu)化方法找到最合適的c/k值.關(guān)于這部分的研究工作可以在后續(xù)展開.

4 結(jié)論

本文以一款四連桿膝關(guān)節(jié)假肢為研究對象，考慮足--地交互作用力和膝關(guān)節(jié)單邊接觸力等強(qiáng)非線性因素，采用光滑化的Kelvin-Voigt 接觸模型和庫倫摩擦模型描述足--地接觸力和摩擦力，采用光滑化的Kelvin-Voigt 接觸模型描述膝關(guān)節(jié)單邊接觸力，并基于第一類拉格朗日方程對四連桿膝關(guān)節(jié)被動假肢進(jìn)行動力學(xué)建模，得到了考慮足--地交互作用力和膝關(guān)節(jié)單邊接觸力的一體化假肢動力學(xué)模型.本文將假肢動力學(xué)仿真結(jié)果與健康被試實測結(jié)果對比，系統(tǒng)分析了假肢動力學(xué)響應(yīng)與四連桿膝關(guān)節(jié)液壓阻尼器的剛度、阻尼參數(shù)與假肢動力學(xué)響應(yīng)和步態(tài)評價指標(biāo)之間的關(guān)系.通過本文的研究，主要結(jié)論如下：

(1)考慮足--地交互作用力、單邊接觸力等強(qiáng)非線性因素后，假肢的動力學(xué)響應(yīng)呈現(xiàn)出豐富的非線性現(xiàn)象，如出現(xiàn)了1/2 和1/3 亞諧波頻率，這種現(xiàn)象破壞了行走時雙腿步態(tài)的協(xié)調(diào)性.在髖關(guān)節(jié)運動中加入提胯動作后，假肢的動力學(xué)響應(yīng)中的亞諧波頻率成分消失，轉(zhuǎn)變?yōu)閱沃芷诜€(wěn)態(tài)響應(yīng)，這為殘疾人行走時提胯等代償行為的提供了一種新的力學(xué)解釋.

(2)假肢膝關(guān)節(jié)液壓阻尼器剛度k和阻尼c對膝關(guān)節(jié)最大屈曲角度有影響.通過減小剛度k和阻尼c可以增大膝關(guān)節(jié)最大屈曲角，這有利于增大足部離地間隙，進(jìn)而提高運動安全性，但過低的剛度有可能導(dǎo)致步態(tài)失諧，從另一方面降低運動安全性.

(3)數(shù)值計算結(jié)果表明，通過調(diào)整膝關(guān)節(jié)液壓阻尼器剛度k和阻尼c，可以使得假肢-健肢步態(tài)相關(guān)系數(shù)達(dá)到0.9 以上，這表明膝關(guān)節(jié)液壓阻尼器剛度k和阻尼c是提高假肢--健肢步態(tài)一致性的關(guān)鍵參數(shù)，為四連桿膝關(guān)節(jié)被動假肢的進(jìn)一步優(yōu)化提供了理論支撐.

本文的研究從動力學(xué)角度出發(fā)，對現(xiàn)有膝關(guān)節(jié)假肢設(shè)計技術(shù)進(jìn)行了補(bǔ)充和完善，為基于動力學(xué)的假肢設(shè)計提供了依據(jù).后續(xù)工作中，我們將以最大屈曲角、步態(tài)一致性等為指標(biāo)，深入研究假肢動力學(xué)性能的多目標(biāo)優(yōu)化.

附錄A

致謝本研究的第2 節(jié)實驗工作由北京大學(xué)高思源協(xié)助完成，特此感謝.