考慮膝關(guān)節(jié)轉(zhuǎn)動剛度的二段腿行走步態(tài)分析

林彬添，張清文, 范 峰

(1.結(jié)構(gòu)工程災(zāi)變與控制教育部重點實驗室(哈爾濱工業(yè)大學(xué))，哈爾濱 150090；2.土木工程智能防災(zāi)減災(zāi)工業(yè)和信息化部重點實驗室(哈爾濱工業(yè)大學(xué))，哈爾濱150090)

倒立擺是模擬人體步態(tài)的基本模型.文獻(xiàn)[1]采用一個由無質(zhì)量的線彈簧腿和質(zhì)點組成的倒立彈簧擺(spring-loaded inverted pendulum, SLIP)來模擬奔跑，見圖1(a)；文獻(xiàn)[2]采用兩條無質(zhì)量的剛性腿和質(zhì)點組成的倒立擺來模擬行走；文獻(xiàn)[3]將行走和奔跑統(tǒng)一用線彈簧倒立擺進(jìn)行模擬，見圖1(b)，采用這種柔性腿模型的優(yōu)點是行走雙足支撐階段能夠得到很好的模擬；文獻(xiàn)[4]進(jìn)一步研究了線彈簧腿倒立擺抵抗擾動穩(wěn)定行走的能力，以及模型行走的效率，給出了模擬人體正常行走的合理參數(shù)范圍.此外，文獻(xiàn)[5]在SLIP模型的基礎(chǔ)上附加了半圓來考慮行走過程中落腳點的前移，見圖1(c)；文獻(xiàn)[6]考慮了軀干和髖關(guān)節(jié)的影響，研究人體行走的姿態(tài)控制和穩(wěn)定性.

這類模型雖然構(gòu)造簡單，但是能夠合理地描述人體質(zhì)心的運動軌跡，其地反力(ground reaction force, GRF)時程具有實測行走荷載時程的“雙峰值”特性[3]，同時動力學(xué)分析方便，因而被廣泛地運用于不同領(lǐng)域，例如，基于倒立彈簧擺抵抗擾動穩(wěn)定行走的特性，文獻(xiàn)[7-8]將其用于設(shè)計機(jī)器人；醫(yī)療領(lǐng)域，文獻(xiàn)[9]將其用于人體負(fù)重行走研究，預(yù)防摔倒；伴隨著建筑結(jié)構(gòu)跨度的增大，人體運動荷載導(dǎo)致的結(jié)構(gòu)振動問題逐漸顯現(xiàn)出來，文獻(xiàn)[10-11]將倒立擺模型與結(jié)構(gòu)動力學(xué)模型耦合，分析人與結(jié)構(gòu)的相互作用機(jī)理.

圖1 倒立擺模型

值得注意的是，真實人體的腿部并不是理想的線彈簧，而是由多個肢體和關(guān)節(jié)構(gòu)成的復(fù)雜機(jī)構(gòu)[12]，極度的簡化也為模型的實驗驗證及參數(shù)標(biāo)定帶來了困難.采用剛性桿件代表肢體，通過具有轉(zhuǎn)動剛度的關(guān)節(jié)進(jìn)行連接，同樣具有柔性腿特性，進(jìn)而能夠模擬雙足支撐階段.為了分析腿部構(gòu)造對倒立擺行走步態(tài)的影響，本文考慮了一種由含膝關(guān)節(jié)的二段腿和質(zhì)點組成的倒立擺模型.本文首先確定了膝關(guān)節(jié)的剛度和初始角度對腿部力學(xué)特性(腿部反力與腿長關(guān)系)的影響，然后建立了模型運動的控制方程和仿真模型，采用龐加萊映射和牛頓迭代法求得了模型周期性行走的步態(tài)結(jié)果，分析了模型參數(shù)對行走步態(tài)的影響.

1 二段腿倒立擺模型

本文采用的二段腿模型見圖2，人體質(zhì)量為m，集中在質(zhì)心(center of mass, CoM)，質(zhì)心坐標(biāo)記為(xm,ym)，由兩條無質(zhì)量的腿進(jìn)行支承，腿部包含兩個剛性肢體，代表小腿和大腿，其中小腿長度為l1，大腿長度為l2，通過具有轉(zhuǎn)動剛度k的膝關(guān)節(jié)連接在一起，定義圖示由大腿延長線轉(zhuǎn)向小腿所成的銳角為膝關(guān)節(jié)角，用θ表示，膝關(guān)節(jié)內(nèi)力矩為0的角度定義為膝關(guān)節(jié)初始角度，用θ0表示，假定腿部在支撐相與地面鉸接，忽略落腳點的前移與腳部相對地面的滑動，落腳點的水平坐標(biāo)記為di，下標(biāo)i表示落腳點對應(yīng)的步數(shù).

圖2 二段腿模型

根據(jù)步態(tài)分析領(lǐng)域的定義，兩個相鄰步態(tài)事件之間的行走過程為一個單步(step)[3]，單步過程起點的選取不影響步態(tài)的討論，本文選取豎直腿擺向(vertical leg orientation, VLO)事件[4]作為單步過程的起點，假定此時模型處于單足支撐階段且質(zhì)心位于落腳點豎直上方，模型行走的第i個單步見圖3.

圖3 單步過程示意

兩個相間步態(tài)事件之間的行走過程為一個跨步(stride)，例如，同側(cè)腿的兩個相鄰VLO事件之間的行走過程，可知一個跨步包含兩個單步.由于擺動腿沒有質(zhì)量，對質(zhì)心運動沒有影響，圖示沒有將其畫出；質(zhì)心由圖示位置向前運動，擺動腿始終與地面保持恒定的沖擊角α，當(dāng)質(zhì)心與假定落腳點的距離恰好為腿初始長度時，擺動腿觸地，進(jìn)入雙足支撐階段，兩個相鄰落腳點間的水平距離即為步長；當(dāng)后支撐腿的反力減小為0時，后支撐腿離地，再次進(jìn)入單足支撐階段，直到質(zhì)心來到下一個落腳點豎直上方，即VLOi+1，代表模型完成了一個單步過程.

對質(zhì)心進(jìn)行受力分析，見圖4，因為模型的腿部沒有質(zhì)量，可以等效為一個二力桿，所以腿部受到的兩個反力必然等值，沿著作用點的連線方向，方向相反.膝關(guān)節(jié)的內(nèi)力矩為k(θ-θ0)，由三角形面積相等可得腿反力相對于膝關(guān)節(jié)的力臂為

(1)

Fh=k(θ-θ0)，

(2)

將式(1)代入式(2)，整理可得腿部的反力大小為

(3)

式中P=k(θ-θ0)/l1l2sinθ為名義腿剛度.

圖4 受力分析示意

質(zhì)心(雙足支撐階段)受到重力mg與腿部作用力Fi和Fi+1，由牛頓第二定律可得：

(4)

(5)

(6)

(7)

(8)

(9)

代入式(4)、(5)，可得質(zhì)心的運動方程：

(10)

(11)

單足支撐階段運動方程形式相同，僅保留支撐腿對應(yīng)項即可.

定義狀態(tài)向量

模型的狀態(tài)方程為

(12)

狀態(tài)向量的維數(shù)為4，由于模型是一個保守系統(tǒng)，機(jī)械能守恒，施加一個約束條件，機(jī)械能給定為E0，采用狀態(tài)空間分析方法[13]，根據(jù)VLO事件建立龐加萊映射[4]，系統(tǒng)獨立的變量維數(shù)縮減為2，獨立變量的選取不影響系統(tǒng)的性質(zhì)，本文參考文獻(xiàn)[4]，選取的獨立變量為VLOi時刻的質(zhì)心高度y0和速度方向角β，見圖5.

圖5 龐加萊映射示意

由此建立的龐加萊映射為

(13)

龐加萊映射的不動點[13]即為模型的周期性行走步態(tài)，本文在MATLAB/SIMULINK環(huán)境建立仿真模型，模型框圖見圖6，模型包括：積分模塊、單雙足支撐相轉(zhuǎn)換模塊和觸地離地判定模塊3個部分，狀態(tài)方程的求解采用ode45變步長求解器.采用牛頓迭代法求取了模型不動點，由此得到了模型周期性行走步態(tài)結(jié)果.

2 二段腿力學(xué)特性

為比較二段腿和線彈簧腿的力學(xué)特性，繪制腿部受到的反力與腿長的關(guān)系曲線，見圖7.根據(jù)文獻(xiàn)[12]可知人體大小腿長度相近，本文取大小腿長度l1、l2均為0.5 m；參考文獻(xiàn)[4]中給出了人體正常行走的腿剛度范圍(應(yīng)大于14mg/l0，l0為線彈簧原長，本文取m=80 kg)，選用原長為1.0 m，線剛度為12 kN·m-1的線彈簧作為對照組，由公式(1)可知腿內(nèi)力正比于膝關(guān)節(jié)轉(zhuǎn)動剛度k，取膝關(guān)節(jié)腿剛度k=400 Nm·rad-1；進(jìn)而可分析膝關(guān)節(jié)初始角度對二段腿動力特性的影響，見圖7的反力-腿長曲線.

根據(jù)文獻(xiàn)[14]可知，人體正常行走時，質(zhì)心的豎直位移幅值為0.091 m，仿真過程中發(fā)現(xiàn)模型的腿長度大于0.9 m，故本文分析的腿長范圍為0.9 m到腿初始長度.由圖7可知，二段腿受壓縮的初期，反力隨著變形的增大上升較快，表明二段腿有較大的初始剛度，隨著變形的逐漸增大，反力的上升較為緩慢，表明二段腿在變形較大時有柔化現(xiàn)象.分析可知，隨著變形的增大，反力的力臂變長，由式(2)可得，單位膝關(guān)節(jié)內(nèi)力矩可以抵抗的反力減小，因而二段腿有受力柔化特性.在本文的分析范圍內(nèi)，膝關(guān)節(jié)初始角度的增大導(dǎo)致反力-腿長曲線的左移，相同腿長條件下，膝關(guān)節(jié)初始角度大的二段腿內(nèi)力矩較小，故能夠抵抗的反力較小.

圖6 SIMULINK模型

圖7 二段腿力學(xué)特性

3 二段腿模型行走步態(tài)

3.1 機(jī)械能的影響

與線彈簧模型類似，二段腿模型能夠在較大的沖擊角范圍內(nèi)實現(xiàn)周期性的行走步態(tài).圖8為膝關(guān)節(jié)轉(zhuǎn)動剛度k=400 Nm·rad-1，初始角度θ0=20°的二段腿模型，在不同機(jī)械能E0下周期性步態(tài)的分布情況.圖8曲線上的每個點均對應(yīng)一個周期性的步態(tài)結(jié)果，由圖8可知，二段腿模型能夠?qū)崿F(xiàn)對稱性和非對稱性(地反力時程曲線的形態(tài)，見圖9)兩類周期性步態(tài)，其中，對稱性步態(tài)在VLO時刻的質(zhì)心速度方向角均為0，而非對稱性步態(tài)不為0，見圖8(c).α-β曲線與β=0的交點為一個對稱性步態(tài)，隨著機(jī)械能的減小，其沖擊角α逐漸增大，由該對稱步態(tài)向兩邊發(fā)展得到的非對稱步態(tài)，其VLO時刻的質(zhì)心高度y0均逐漸減小，見圖8(b).觀察圖8(c)可知，機(jī)械能較小的非對稱性步態(tài)沖擊角較大，因而擺動腿觸地時質(zhì)心的高度較大，同時由圖8(b)可知，VLO時刻質(zhì)心的高度y0隨機(jī)械能減小而降低，當(dāng)機(jī)械能進(jìn)一步減小，VLO時刻模型處于單足支撐階段的假定將不能滿足，不能求得非對稱性步態(tài)，本文得到的非對稱性步態(tài)結(jié)果可行域為E0不小于790 J.觀察圖8(a)可知，對稱性步態(tài)在α-y0坐標(biāo)上存在兩個分支，當(dāng)機(jī)械能E0=820 J時，兩個分支分別位于左上角和右下角，隨著機(jī)械能E0的減小，曲線的分布形式發(fā)生了改變，例如，當(dāng)機(jī)械能E0=800 J時，兩分支分別位于圖示左下角和右上角.

圖9展示了圖8所示步態(tài)A至F的一跨步地反力時程曲線.由圖9(b)可知，該分支上的步態(tài)，其地反力時程與實測地反力時程不同，僅有一個峰值；由圖9(d)可知，該分支的步態(tài)結(jié)果步頻較高，且雙足支撐階段時間占步行周期較大比例，行走效率較低[4]，故后續(xù)不再討論這兩種步態(tài).圖9(a)、(c)、(e)和(f)具有人體行走地反力的“雙峰值”特性，通常來說，人體行走的落步荷載時程是非對稱的，其第二個峰值稍小于第一個峰值[15-16]，步態(tài)F所在分支更接近于真實行走落步荷載.

圖8 機(jī)械能對步態(tài)分布的影響

圖9 豎向地反力時程

對比圖9(a)、(c)可知，機(jī)械能E0的降低將導(dǎo)致地反力時程峰值減小，同時峰值之間的低谷值將增大.為分析機(jī)械能對步態(tài)參數(shù)，包括步頻、步長和平均步速的影響，以步態(tài)A為對照組，僅改變機(jī)械能E0，各步態(tài)的步態(tài)參數(shù)見表1，由表1可知，機(jī)械能的減小將導(dǎo)致步頻降低，步長增大，平均步速減小.提高平均步速最有效的方式為增加能量輸入，這與線彈簧模型結(jié)論一致[3].

表1 機(jī)械能對步態(tài)參數(shù)的影響

3.2 沖擊角的影響

為分析沖擊角對模型步態(tài)參數(shù)的影響，表2列舉了步態(tài)A及其鄰近步態(tài)的步態(tài)參數(shù).由表2可知，隨著沖擊角的增大，對稱性步態(tài)的步頻逐漸提高，步長逐漸減小，平均步速則稍有降低.文獻(xiàn)[17]指出人體正常行走的步頻范圍為1.6 Hz到2.4 Hz，結(jié)合表1、2可知，本文展示的步頻范圍為1.38 Hz到2.43 Hz，可見二段腿模型在給定參數(shù)范圍下能夠較好地覆蓋人體正常行走的步頻區(qū)間.

表3列舉了非對稱性步態(tài)F及其鄰近步態(tài)的步態(tài)參數(shù)，由表3可知，在本文分析的沖擊角范圍內(nèi)，非對稱性步態(tài)的步態(tài)參數(shù)對于沖擊角的變化不敏感.分析可知，對稱性步態(tài)處于VLO狀態(tài)時速度方向角為0，沖擊角控制擺動腿的觸地判定，直接影響步態(tài)參數(shù)；非對稱性步態(tài)處于VLO狀態(tài)時速度方向角不為0，擺動腿的觸地判定由沖擊角和速度方向角共同決定，步態(tài)參數(shù)基本不變，然而地反力兩個峰值的大小關(guān)系發(fā)生了變化，見圖9(e)、(f).

表2 沖擊角對步態(tài)參數(shù)的影響(對稱步態(tài))

表3 沖擊角對步態(tài)參數(shù)的影響(非對稱步態(tài))

3.3 膝關(guān)節(jié)初始角度的影響

為分析膝關(guān)節(jié)初始角度對步態(tài)分布的影響，圖10展示了機(jī)械能E0=820 J，膝關(guān)節(jié)轉(zhuǎn)動剛度k=400 Nm·rad-1條件下，不同膝關(guān)節(jié)初始角度的二段腿模型的步態(tài)分布情況.由圖10(a)可知，對于對稱性步態(tài)，以步態(tài)A所在分支為參考，隨著初始角度的減小，步態(tài)分布曲線右移，表明沖擊角增大；對于非對稱性步態(tài)，以步態(tài)F所在分支為參考，隨著初始角度的減小，步態(tài)分布曲線右移.由此可見，膝關(guān)節(jié)初始角度較小的二段腿模型產(chǎn)生的步態(tài)具有較大的沖擊角.

表4分析了膝關(guān)節(jié)初始角度對步態(tài)參數(shù)的影響，以步態(tài)A為對照組，減小膝關(guān)節(jié)初始角度將導(dǎo)致步頻降低，步長增大，然而平均行走速度的變化很小.由圖7可知，在本文討論的參數(shù)范圍內(nèi)，膝關(guān)節(jié)初始角度的減小導(dǎo)致曲線右移，其效果近似于腿長變長，在沖擊角相同的條件下，如圖10(a)所示，VLO時刻的質(zhì)心高度增大，VLO到擺動腿觸地耗時變長，因而步頻降低，步長增大.

圖10 膝關(guān)節(jié)初始角度對步態(tài)分布的影響

表4 膝關(guān)節(jié)初始角度對步態(tài)參數(shù)的影響

3.4 膝關(guān)節(jié)轉(zhuǎn)動剛度的影響

為分析膝關(guān)節(jié)初始角度對步態(tài)分布的影響，圖11展示了機(jī)械能E0=820 J，膝關(guān)節(jié)初始角度θ0=20°條件下，不同膝關(guān)節(jié)轉(zhuǎn)動剛度的二段腿模型的步態(tài)分布情況.由圖11(a)可知，對于對稱性步態(tài)，以步態(tài)A所在分支為參考，隨著膝關(guān)節(jié)轉(zhuǎn)動剛度的增大，步態(tài)分布曲線對應(yīng)的沖擊角區(qū)間擴(kuò)大，表明剛度更大的腿能夠在更大的沖擊角范圍內(nèi)實現(xiàn)周期性行走；對于非對稱性步態(tài)，以步態(tài)F所在分支為參考，隨著膝關(guān)節(jié)轉(zhuǎn)動剛度的增大，步態(tài)分布曲線右移.由此可見，膝關(guān)節(jié)轉(zhuǎn)動剛度較大的二段腿模型產(chǎn)生的非對稱性步態(tài)具有較大的沖擊角.

圖11 膝關(guān)節(jié)轉(zhuǎn)動剛度對步態(tài)分布的影響

表5分析了膝關(guān)節(jié)轉(zhuǎn)動剛度對步態(tài)參數(shù)的影響，以步態(tài)A為對照組，增大膝關(guān)節(jié)轉(zhuǎn)動剛度將導(dǎo)致步頻小幅度降低，步長小幅度增大，平均行走速度稍有增大.由圖11(a)可知，沖擊角相同的條件下，雖然剛度更大的二段腿VLO時刻的質(zhì)心高度較低，但是腿部對質(zhì)心豎向的支撐作用更強，從VLO到擺動腿觸地耗時依然更長，其步頻更低.

表5 膝關(guān)節(jié)轉(zhuǎn)動剛度對步態(tài)參數(shù)的影響

4 結(jié) 論

本文對一種由膝關(guān)節(jié)連接的二段腿模型開展了建模分析，得到的主要結(jié)論如下：

1) 帶膝關(guān)節(jié)的二段腿有較大的初始剛度，隨著變形的增大，剛度逐漸降低，即存在受力柔化現(xiàn)象.

2) 二段腿模型能夠在一定的參數(shù)范圍內(nèi)實現(xiàn)多種形式的周期性步態(tài)，可分為對稱性和非對稱性兩大類.

3) 模型參數(shù)對步態(tài)參數(shù)的影響可總結(jié)為：增大模型的機(jī)械能，平均行走速度提升最顯著；沖擊角對對稱性步態(tài)的步態(tài)參數(shù)影響顯著，非對稱性步態(tài)的步態(tài)參數(shù)對沖擊角不敏感；減小二段腿的膝關(guān)節(jié)初始角度，改變步態(tài)的頻率和步長關(guān)系，而平均行走步速變化很小；膝關(guān)節(jié)轉(zhuǎn)動剛度增大將小幅度降低步頻，增大步長和平均行走速度.

4) 二段腿模型能夠很好地模擬人體行走地反力的“雙峰值”特性，給定合理的參數(shù)組合，模型能夠覆蓋人體正常行走的步頻區(qū)間.

可見二段腿模型同樣具有柔性腿特征，能夠有效模擬人體行走雙足支撐階段，重現(xiàn)地反力時程的雙峰特性，同時，其構(gòu)造更加接近于真實人體，參數(shù)物理意義明確，為模型的實驗驗證及參數(shù)標(biāo)定打下了基礎(chǔ).進(jìn)一步的人體行走試驗有待開展，測量人體運動學(xué)及動力學(xué)參數(shù)，探討簡化模型的合理性.