四连杆膝关节假肢的动力学建模与分析1)

吕 阳 方虹斌 徐 鉴 马建敏 王启宁 张晓旭

(复旦大学航空航天系，上海 200433)

(复旦大学智能机器人研究院，上海 200433)

(上海智能机器人工程技术研究中心，复旦大学，上海 200433)

(智能机器人教育部工程研究中心，复旦大学，上海 200433)

(北京大学工学院，北京 100871)

引言

根据全国第二次残疾人抽样调查的结果[1]，我国肢体残疾2412 万人，占残疾人总数的29.07%，其中下肢截肢患者大概有158 万人.为了使大腿截肢患者能够正常地生活和工作，越来越多的研究者开始致力于假肢的结构设计与优化研究.

膝关节假肢在结构上可分为单轴式膝关节和多轴式膝关节.单轴式膝关节具有结构简单、可靠性高的优点，其代表性产品主要有奥托博克公司的CLeg、奥索公司的Rheo Knee 等[2].多轴式膝关节如四连杆膝关节则具有良好的仿生特性[3]，能够较为真实地逼近人体膝关节转动瞬心(小腿相对大腿转动时的速度瞬心，简称转动瞬心)的J 形轨迹，其代表性产品主要有奥索公司的Parso Knee、北京工道风行公司的K301 等.值得一提的是，近年来，折纸结构由于其结构轻、可编程以及超材料特性被广泛研究，这些特性在膝关节的J 形轨迹实现上也有一定的应用潜力[4].在运动性能上，J 形轨迹意味着小腿假肢在站立时的瞬心较高，从而实现较好的稳定性；在摆动时的瞬心较低，有较好的灵活性[5].因此，相比于单轴式膝关节，四连杆膝关节兼顾运动稳定性，具有更好的运动性能.

为了对四连杆膝关节的运动性能有更好的认识，需要明确影响其运动性能的因素.从动力学的角度来看，影响运动性能的因素包括足--地交互作用力、膝关节单边约束产生的强非线性作用力、液压阻尼器提供的弹性回复力及阻尼力等.要明确各因素对假肢性能的影响机制，就必须建立考虑足--地交互作用力和膝关节单边接触力的一体化动力学模型，并就该模型产生的动力学响应进行分析.

目前，对于假肢的建模，常用的建模方法是利用第一类拉格朗日方程.这种方法将假肢运动阶段分为站立相(假肢与地面接触，处于站立状态)和摆动相(假肢未与地面接触，处于摆动状态)，对不同的步态阶段分别进行建模[6-10].在站立相阶段，该方法将足--地接触简化为理想铰接约束，通过在动力学方程中引入拉格朗日乘子求解模型；在膝关节存在单边约束时，假设大腿和小腿为一根刚性杆，引入拉格朗日乘子对模型进行求解.这种建模方法未对足--地交互作用力以及膝关节单边接触力直接进行建模，而是基于假设将这两种力看作理想约束力.实际上足--地交互作用力以及膝关节单边接触力等强非线性因素为非理想约束力，上述处理方式无法体现这一点，因此需要运用更合理的方式对这两种力进行建模.

足--地交互作用力本质上是足--地接触时产生的接触力和摩擦力，膝关节单边接触力本质上是膝关节连杆间限位碰撞产生的接触力，因而其建模方式可以借鉴多刚体动力学对接触力和摩擦力的描述.接触力模型中有赫兹接触力模型[11]以及改进的赫兹接触力模型[12].赫兹接触力模型形式相对简单，但不考虑能量耗散，不适用于真实的接触情况；改进的赫兹接触力模型通过各种方式加入了能量耗散项，常用的有Hunt-Crossley 模型[13]、Gonthier 模型[14]、Kelvin-Voigt[15]模型等，更适用于真实的接触情况.摩擦力模型主要包括：库伦摩擦模型[16]、Stribeck 摩擦模型[17-19]以及更为复杂的LuGre 摩擦模型[20-23].实验研究[24-27]表明，当刚体相对运动尺度为毫米量级或更高时，若要分析摩擦对系统动力学行为的影响，采用库伦摩擦模型即可达到足够高的精度.为解决库伦摩擦模型中非光滑本构对动力学分析带来的困难，还可进一步利用非线性函数[28]对库伦摩擦模型进行光滑化.近年来，基于上述接触和摩擦模型的双足行走机器人[29-32]动力学分析已经取得一定进展，但基于上述接触和摩擦模型的假肢动力学分析仍处于起步阶段.

综上所述，足--地交互作用力以及膝关节单边接触力等强非线性因素对假肢步态的影响尚未得到深入分析，其难点在于考虑足--地交互作用力和膝关节单边接触力的一体化动力学建模.针对该问题，基于多刚体动力学中关于接触力及摩擦力的研究，本文将采用Kelvin-Voigt 模型和库伦模型描述足--地接触力和摩擦力，采用Kelvin-Voigt 模型描述膝关节单边接触力，从而基于第一类拉格朗日方程对四连杆膝关节被动假肢进行动力学建模.将步态实验测得的髋关节运动数据作为动力学模型的驱动信号，四连杆膝关节被动假肢的步态可以通过数值计算得到.进一步，为了评价假肢步态与健康人实测步态的一致性，本文定义了步态相关系数作为评价指标，并分析四连杆膝关节液压阻尼器刚度、阻尼参数对相关系数的影响.

本文的贡献和创新点在于：首先，进一步完善四连杆膝关节假肢的动力学模型，重点关注足--地交互作用力以及膝关节单边接触力等强非线性因素对下肢假肢步态的影响，为残疾人行走时的代偿行为提供一种新的力学解释；其次，以实测数据为参考信号进行步态一致性评估和参数分析，为四连杆膝关节被动假肢的进一步优化提供理论支撑，提高结论的可靠性.

本文第1 节介绍四连杆膝关节假肢动力学模型的建立.首先分析四连杆的运动学特性，接着引入足--地交互作用力和膝关节单边接触力，最后基于第一类拉格朗日方程对假肢进行建模.第2 节介绍健康被试的步态测试.一方面，实测髋关节运动数据作为假肢动力学模型的驱动信号，用于假肢步态分析；另一方面，实测膝关节运动数据作为假肢步态运动的对比信号，用于步态一致性评估及参数分析.第3 节为基于假肢动力学模型的步态分析.将数值计算所得结果与实验测得的步态进行比较，讨论四连杆液压阻尼器刚度、阻尼参数对步态周期失谐及步态一致性的影响.第4 节给出研究结论.本文通过数值计算研究了四连杆膝关节被动假肢动力学模型的响应，分析了四连杆膝关节液压阻尼器的阻尼、刚度参数对假肢步态指标的影响，总结了提高步态性能的一般经验，以期为四连杆膝关节被动假肢的参数设计提供参考依据.

1 四连杆膝关节被动假肢动力学模型

四连杆膝关节被动假肢具有良好的仿生特性，但由于结构更为复杂，足--地交互作用力、膝关节单边接触力等强非线性因素会导致下肢假肢产生更为丰富的动态响应，进而在整体上影响步态的周期性和穿戴舒适性.因此，本文有必要从动力学角度对四连杆膝关节假肢进行分析.

1.1 四连杆机构的运动学模型

如图1 所示，人体膝关节由股骨踝(femoral ankle)和胫骨踝(tibia ankle)组成，依赖于股踝和胫踝接触面的几何外形，人体膝关节的运动属于非定轴转动，小腿相对于大腿转动的速度瞬心轨迹实际是一条J 形曲线.其生物学优势在于，在站立相时，膝关节瞬心升高，小腿相对于瞬心的转动惯量较大，不易实现转动，因而稳定性好；反之，在摆动相时，小腿灵活性好，这也是四连杆膝关节仿生学研究的主要动机.

图1 人体膝关节非定轴转动及其J 形瞬心轨迹Fig.1 Non-fixed axis rotation of human knee joint and its J-curve of instantaneous rotation center

四连杆膝关节可以较好地模拟人体膝关节非定轴转动过程，主要由四根连杆和液压阻尼器组成，其结构示意图及实物图如图2 所示.以工道风行公司(GDFX)设计的一款四连杆膝关节为例，其几何参数及惯量参数如表1 所示.根据四连杆的连接关系，四连杆机构的运动约束方程可以表示为

在二维笛卡尔坐标系下，式(1)可以具体表示为

图2 四连杆膝关节Fig.2 Four-bar linkage prosthesis knee

表1 四连杆膝关节物理参数表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为小腿股骨与竖直方向夹角.

由几何关系可知，四连杆膝关节在转动过程中小腿相对于大腿的瞬心P在CB和DA的延长线上.在平面笛卡尔坐标系下，瞬心P的坐标可以表示为

其中，xP为瞬心P的横坐标，yP为瞬心P的纵坐标，xk和yk分别为四连杆顶点k(k=A,B,C,D)的横坐标和纵坐标.图3 展示了小腿相对于大腿转动时速度瞬心P的轨迹曲线，其中假肢大腿与AB杆刚性连接，假肢小腿与CD杆刚性连接.可以看到，在膝关节由伸直状态至110◦屈曲状态过程中，瞬心P形成的轨迹是一条J-curve 线，与人体膝关节瞬心的J 形轨迹相似，展现了良好的仿生特性.

1.2 足--地交互力和膝关节单边接触力模型

足--地交互作用力主要有地面对假肢足底的接触力和摩擦力两部分.为方便分析，本文假设只有脚跟和脚尖这两点与地面存在接触力和摩擦力.

地面对足底的接触力可以通过Kelvin-Voigt 模型来描述

图3 四连杆膝关节J 形瞬心轨迹Fig.3 Four-bar linkage prosthesis knee’s J-curve of instantaneous rotation center

式中，FN,i表示接触力，方向垂直于地面，下标i为heel 或toe，分别表示脚跟和脚尖与地面接触时接触力的描述.KN,i表示接触刚度，χN,i表示接触阻尼.δi表示脚尖或脚跟与地面之间的压痕深度，当脚尖或脚跟离地时，规定δi=0.地面对足底的摩擦力采用库伦摩擦模型描述[16]

式中，Ff,i表示摩擦力，方向沿水平方向，下标i为heel 或toe，分别表示脚跟和脚尖与地面接触时摩擦力的描述.µ表示动摩擦系数(此处忽略静摩擦系数)，vi表示脚尖或脚跟相对地面的切向速度.为方便动力学分析，本文对式(4)所述接触力模型和式(5)所述摩擦力进行光滑化处理，具体光滑化形式分别为

其中C1和C2为常数，表征接触力和摩擦力的光滑化程度.

此外，如图4 所示，当假肢小腿摆动到与大腿呈同一角度时(即膝关节打直)，为保证穿戴者安全性，四连杆的AB和BC两杆之间将通过限位装置锁死，锁死端面垂直于AB连杆.这种结构锁死是一种单边约束，力学上可以简化为单边接触模型.因此，本文同样通过光滑化的Kelvin–Voigt 模型来描述膝关节打直的单边接触力，具体形式为

其中，Fc表示单边接触力，方向为AB连杆方向，Kc表示接触刚度，Cc表示接触阻尼，δknee表示两杆之间的压痕深度.需要注意的是，如图4 中细节放大图所示，AB和BC两杆间限位接触点与B点不重合，其偏心距为d=1 cm，因此式(8)所示接触力作用于连杆BC和AB的力矩分别为

至此，本文得到了足--地交互作用力和膝关节单边接触力的模型.

图4 膝关节几何锁死时接触力示意图Fig.4 The contact force when the knee is locked geometrically

1.3 基于第一类拉格朗日方程的动力学建模

假设四连杆膝关节的接受腔与残疾人大腿残肢为刚性连接，因此四连杆膝关节被动假肢的运动模型可简化为图5 所示形式.其中x和y分别为髋关节在平面笛卡尔坐标系下的横坐标和纵坐标，α1为大腿与竖直方向的夹角，α11为小腿与竖直方向夹角，β1为足部与水平方向夹角.

如图5 所示，四连杆机构中的液压阻尼器主要是在运动过程中提供阻尼力和回复力，其合力可以表示为

式中，k和c分别为液压阻尼器的刚度和阻尼系数，l0=168.79 mm 为液压阻尼器原长，l为当前长度.基于虚功原理，合力F0作用于广义坐标θ1上的广义力为

图5 四连杆膝关节被动假肢运动模型Fig.5 Model of passive four-bar linkage prosthesis knee

此外，由于θ1=0.5π −α1+θA，因此合力F0作用于广义坐标α1和θA上的广义力矩分别为−τ1和τ1.

此外，本文假设踝关节受到带阻尼的扭簧作用，作用力矩表示为

其中，kankle=10 N/m 和cankle=0.2N/(m·s−1)分别为扭簧的刚度和阻尼系数，为踝关节角度，为踝关节期望角度.在步态分析中，本文以实测踝关节角度作为.由图5 可知，踝关节角度满足=β1−α11，因此扭簧作用于广义坐标α11和β1的扭矩分别为−Tankle和Tankle.

得到基于四连杆约束的运动学模型、足--地交互力和膝关节单边接触力模型后，图5 所示的四连杆膝关节被动假肢的动力学模型可以用第一类拉格朗日方程表示

式中，q=(x,y,α1,α11,θA,θB,β1)T为广义坐标向量，是质量阵，是科氏力和离心力项，是重力项.M(q),及N(q)中各元素表达式见附录A.Fe为足--地交互作用力和膝关节单边接触力构成的广义力向量，具体表示为

其中，FN,heel和FN,toe由式(6)给出，Ff,heel和Ff,toe由式(7)给出，Tc和由式(9)给出，J为将足--地交互作用力和膝关节单边接触力变换为广义力的Jacobi 矩阵，具体形式见附录A.Fq为膝关节的液压阻尼器和踝关节的阻尼弹簧施加的广义力向量，具体表示为

其中，τ1由式(11)给出，Tankle由式(12)给出.Φ 为四连杆约束矩阵，具体表示为

其中，ξ1和ξ2为式(2)所示约束方程.λ 为拉格朗日乘子，表示四连杆约束产生的约束力.

需要指出的是，被动假肢是由髋关节运动驱动的，髋关节运动信号x，y，α1为假肢系统的输入，α11，θA，θB和β1为假肢系统的响应.为方便步态分析，式(13)可以改写为

将式(18)对时间t求导可得

代入式(17)可得

进一步，将式(20)代入式(19)可解得

式(17)可以改写为如下状态方程形式

至此，我们得到了考虑四连杆约束、足--地交互力和膝关节单边接触力的四连杆膝关节被动假肢的动力学模型.如前所述，被动假肢是由髋关节运动驱动的，髋关节运动信号q1=(x,y,α1)T及其速度和加速度为式(22)的输入.后续，在基于该模型的动力学分析中，q1为健康被试的实测运动信号.

2 健康被试的步态测试

该人体步态测试实验已获得北京大学伦理委员会批准.实验目的在于：第一，获得式(22)右端所需要的输入q1，，用于四连杆膝关节被动假肢的动力学仿真；第二，将健康被试的步态信号作为健肢步态信号，与动力学仿真所得假肢步态进行比较，用于假肢运动与健肢运动的步态分析.

实验所用设备为Cosmos Gaitway 跑台和Vicon光学运动捕捉系统及配套分析系统.其中，Cosmos Gaitway 跑台用于实现匀速步行环境，Vicon 光学运动捕捉系统用于捕捉人体靶点运动信号，并换算为各关节角度信息.参考步态分析软件OrthoTrak 用户指南中的Helen-Hayes 模型，健康被试身共贴有18 个靶点，分别为右肩靶点(R.shoulder)、左肩靶点(L.shoulder)、补偿靶点(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).步态测试实景如图6 所示，图中Y轴沿着被试行走正前方，X轴沿着被试行走侧向，Z轴沿着竖直向上方向.实验过程中，健康被试(男性，身高182 cm，体重92 kg)在Cosmos Gaitway 跑台上行走(坡度0，速度1 m/s)，Vicon 光学运动捕捉系统采样频率为100 Hz，采样时长180 s.

图6 健康被试步态测试过程Fig.6 Able-bodied subject test process

实验最终得到的数据是各靶点的空间坐标，采用Helen-Hayes 模型的角度解算方法可以得到髋关节处的空间坐标和髋、膝、踝关节的空间转角[33].式(22)所示动力学模型为平面模型，根据坐标对应关系，本文取图6 中髋关节处的Y轴坐标和Z轴坐标分别作为广义坐标q中x，y坐标，取各关节绕X轴方向的转角作为广义坐标q中的α1，α11，β1坐标.

为使实验数据反映健康人步态的普遍规律，本文以两次脚跟着地(heel strike，HS)之间的时间间隔作为一个步态周期T，选取实验数据中的30 个周期进行平均，同时将时间归一化.处理后得到的q1数据如图7 所示.图7(a)中的蓝色实线表示平均后的q1时程，阴影部分的上下界分别表示30 个周期中各个时刻点的最大值和最小值.从阴影部分的面积可以看出，各周期步态较为一致，实验结果较为可靠，因此取30 个周期进行平均是合理的，能够准确地反映q1广义坐标的运动特征.图7(b)中展示了广义坐标y和α1的频谱图(广义坐标x的运动不具有周期性，因此无需进行频谱分析).可以看到，y信号中2 Hz 频率成分的振幅较大，这是因为健康被试在一个步态周期中，左腿和右腿各走了一步，导致髋关节在竖直方向呈现两次近似相同的起伏；α1信号中1 Hz频率成分的振幅较大，这一点并不难理解.上述分析可以说明输入信号q1的激励频率为1 Hz，这为第3 中对系统响应的分析奠定基础.

图7 髋关节实测运动数据时程和频谱曲线Fig.7 Time history and response spectrum of actually measured motion of hip joint

3 四连杆被动假肢的步态分析

基于式(22)并以第2 节的实验数据作为输入，对四连杆膝关节假肢进行动力学仿真，仿真所采用的物理参数如表2 所示.足--地交互作用力和膝关节单边接触力的模型参数如下

为保证仿真所得数据为假肢的稳态响应，仿真时长设置为20 个步态周期，取末5 个步态周期用于分析.

表2 仿真中的部分物理参数值Table 2 Physical parameters in numerical simulation

3.1 亚谐波频率响应

动力学仿真中，设置液压阻尼器阻尼/刚度比为c/k=0.05，刚度系数分别为k=1.4 × 104N/m和k=1.52×104N/m.图8(a)展示了健康被测的实测步态数据及基于动力学模型的仿真数据，其中纵坐标定义为=α1−α11.则从图8(a)阴影部分可以看出，当k=1.4×104N/m 时，四连杆膝关节被动假肢的动力学仿真响应的周期为2T；当k=1.52×104N/m，动力学仿真响应的周期为3T.

图8 亚谐波响应时程曲线与频谱曲线Fig.8 Time history and response spectrum of subharmonic response

进一步地，我们对动力学响应信号进行FFT频谱分析.图8(b)展示了对应于两个时间历程的频谱图.图8(b)表明，当k=1.4 × 104N/m 时，主要的频率成分为1 Hz 和2 Hz，同时出现了1/2 Hz 频率成分，这一亚谐频率成分使得响应的周期为2T；当k=1.52×104N/m 时，主要频率成分仍为1 Hz 和2 Hz，同时出现了1/3 Hz 频率成分，这一亚谐频率成分使得响应的周期为3T.这些亚谐波频率成分出现的原因在于，足--地交互作用力和膝关节单边接触力含有强非线性因素，可以诱发四连杆被动假肢的亚谐波响应这一强非线性动力学行为.1/2 或1/3 亚谐波响应表明，此时假肢的步态周期为健肢步态周期的2 倍或3 倍，破坏了行走时双腿步态的协调性，本文将这种亚谐波响应称为步态周期失谐.

通过进一步的非线性动力学分析，我们发现，通过调整激励参数(如初始条件、相位、幅值等)，系统的稳态响应可以切换回单周期状态.反映到式(22)所示的四连杆膝关节被动假肢模型中，即表现为髋关节运动信号的调整.例如，图9(a)和图9(b)分别展示了将髋关节竖向运动信号y放大为原信号1.006 倍后假肢膝关节角度的时程曲线和频谱曲线.从图9(a)可以看出，在增大髋关节竖向运动幅值后，假肢响应在两个刚度值下均达到单周期稳态，周期和实测数据相同为T，并与实测数据有较好的步态一致性.图9(b)表明，假肢的亚谐波频率成分消失，其基频与健肢的基频相同，为1 Hz，步态周期不再失谐.

图9 提胯后响应时程曲线与频谱曲线Fig.9 Time history and response spectrum when lifting the hip

图9 提胯后响应时程曲线与频谱曲线(续)Fig.9 Time history and response spectrum when lifting the hip(continued)

放大髋关节竖向运动信号所反映的生物意义类似于提胯，即残疾人在运动过程中提高胯部动作幅度的行为.图9 呈现的结论说明，提胯动作将改善假肢步态周期的协调性，使步态不再失谐，这在定性上解释了残疾人穿戴假肢行走时会出现提胯等代偿动作产生的原因.

注意到，在式(22)所示的四连杆膝关节被动假肢动力学模型中，四连杆的几何参数、模型所受激励信号均是基于实验实测所得，因而所得结果具有较高的可信度.但上述亚谐波频率成分和亚谐波响应，尚未在四连杆膝关节的研究文献中见到相关报道.

数值分析进一步发现，保持四连杆膝关节液压阻尼器阻尼比值不变，在k值达到1.6×104N/m 以上时，亚谐波频率响应基本消失，假肢膝关节的运动恢复为单周期的稳态解，因此将k=1.6×104N/m 称为此时的临界k值.图10 展示了k=1.6×104N/m 时假肢动力学仿真时程和响应频谱，显见假肢响应达到了单周期稳态，且周期与实测数据相同.

图10 单周期稳态下假肢动力学仿真步态时程和响应频谱Fig.10 Time history and response spectrum of prosthesis dynamic simulation gait under uniperiodic steady state

图11 不同刚度值下假肢膝关节频率响应图Fig.11 Response spectrum of prosthesis knee under different values of stiffness

图11 给出了c/k=0.05 时各个刚度值下膝关节的时程响应的频谱图.可以看到，在刚度值在(1.4～1.65)×104N/m 变化的过程中，响应在红色区域(1.45～1.47)×104N/m 出现了持续且显著的1/2 亚谐波频率成分，在(1.47～1.5)×104N/m 进入多亚谐波频率成分混杂的阶段，在蓝色区域(1.5～1.56)×104N/m 出现了持续且显著的1/3 亚谐波频率成分，随后各亚谐波频率成分逐渐减少并进入单一稳态阶段(1.56～1.65)×104N/m.另外值得注意的是，在(1.4～1.45)×104N/m 阶段，响应在(1.4～1.42)×104N/m 时出现了1/2 亚谐波频率成分，而在(1.42～1.45)×104N/m 时出现多亚谐波频率成分混杂的情况，这一部分值得在之后的研究中进一步探讨.根据以上分析，可以认为在刚度值增大到一定程度时，亚谐波频率成分会消失，将亚谐波频率成分完全消失时的k值叫作临界k值，并在表3 中给出了其他c/k值下的临界k值.

分析表3 可以发现，阻尼比值越大，单周期稳态解对应的临界k值越低.这表明，可以通过适当调大四连杆膝关节液压阻尼器的阻尼、刚度系数来改善假肢步态周期的协调性.

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

3.2 假肢膝关节最大屈曲角度

图12 不同阻尼比值下膝关节最大屈曲角Fig.12 Peak of knee angle under different damping ratios

可以看出，随着k值增大，膝关节最大屈曲角度呈现递减趋势.同时可以看出，随着阻尼比值的增大，假肢膝关节最大屈曲角与实测膝关节最大屈曲角相等时所对应的刚度值k(红色竖直虚线所对应的k值)也呈现递减规律.这表明，若调大膝关节液压阻尼器的阻尼、刚度系数，尽管假肢步态周期的协调性能够得到改善，但假肢膝关节的最大屈曲角却有可能偏离最优值.因此，未来的工作中，我们有必要选择或定义合理的步态协调性指标，以便对四连杆膝关节假肢进行优化.

3.3 假肢与健肢膝关节步态一致性

为了定量衡量假肢与健肢步态的一致性，本文定义假肢与健肢膝关节信号的相关系数(correlation coefficient)，计算公式如下

式中X=α1−α11为动力学仿真所得膝关节角度信号，Y为健康被试实测膝关节角度信号，Cov(X,Y)为X和Y信号的协方差，D(X)和D(Y)分别为X和Y信号的方差，t1/T和t2/T分别为15 和20.理论上，，ρXY越大，表示假肢膝关节角度信号与健肢膝关节角度信号一致性越高.当且仅当X≡Y时，ρXY=1.

图13 不同阻尼比值下假肢-健肢膝关节相关系数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呈现抖动的状态，这是因为，此时假肢的响应出现亚谐波频率成分，导致ρXY出现极小值点，呈现抖动的状态.根据3.1 节的分析可知，k在超过一定值时，假肢响应会转变为单周期稳态响应，图13 中已根据表3 标出临界k值(红色竖直虚线所对应的k值).在达到临界k值后，ρXY呈现先增大后减小的规律性变化.这进一步验证了3.1 节所述临界k值的合理性.

在转变为单周期稳态响应之后，在不同的c/k下，随着k增大，ρXY均呈现先增大后减小的规律，因此ρXY存在一个最大值点.这表明我们可以通过四连杆膝关节液压阻尼器阻尼、刚度参数的设计，对假肢的步态一致性进行优化.表4 反映了不同液压阻尼器阻尼比值设置下，最大相关系数ρXY及其对应的液压阻尼器刚度值.可以看出，在不同的液压阻尼器阻尼比值下，最大相关系数的值均相近，且最大可得到0.989 6，这为设计液压阻尼器的刚度和阻尼值提供了一定的理论依据.

表4 最大相关系数及对应液压阻尼器刚度Table 4 Max ρ and its corresponding hydraulic damper stiffness

同时可以发现，随着c/k增大，ρXY在达到最大值后的下降速度增大.因此可以发现，c/k越小，ρXY越能在一个较大的k范围内保持平缓的变化.这表明，为了使假肢与健肢在较大的k值范围内达到较好的步态一致性，c/k不宜取得过大.

以上分析表明，为了实现最好的步态一致性，可以通过单目标优化方法找到最合适的k值和c值.另外，也可以围绕着步态一致性探究一些多目标优化问题，如为了实现在较大的k值范围内达到较好的步态一致性，通过多目标优化方法找到最合适的c/k值.关于这部分的研究工作可以在后续展开.

4 结论

本文以一款四连杆膝关节假肢为研究对象，考虑足--地交互作用力和膝关节单边接触力等强非线性因素，采用光滑化的Kelvin-Voigt 接触模型和库伦摩擦模型描述足--地接触力和摩擦力，采用光滑化的Kelvin-Voigt 接触模型描述膝关节单边接触力，并基于第一类拉格朗日方程对四连杆膝关节被动假肢进行动力学建模，得到了考虑足--地交互作用力和膝关节单边接触力的一体化假肢动力学模型.本文将假肢动力学仿真结果与健康被试实测结果对比，系统分析了假肢动力学响应与四连杆膝关节液压阻尼器的刚度、阻尼参数与假肢动力学响应和步态评价指标之间的关系.通过本文的研究，主要结论如下：

(1)考虑足--地交互作用力、单边接触力等强非线性因素后，假肢的动力学响应呈现出丰富的非线性现象，如出现了1/2 和1/3 亚谐波频率，这种现象破坏了行走时双腿步态的协调性.在髋关节运动中加入提胯动作后，假肢的动力学响应中的亚谐波频率成分消失，转变为单周期稳态响应，这为残疾人行走时提胯等代偿行为的提供了一种新的力学解释.

(2)假肢膝关节液压阻尼器刚度k和阻尼c对膝关节最大屈曲角度有影响.通过减小刚度k和阻尼c可以增大膝关节最大屈曲角，这有利于增大足部离地间隙，进而提高运动安全性，但过低的刚度有可能导致步态失谐，从另一方面降低运动安全性.

(3)数值计算结果表明，通过调整膝关节液压阻尼器刚度k和阻尼c，可以使得假肢-健肢步态相关系数达到0.9 以上，这表明膝关节液压阻尼器刚度k和阻尼c是提高假肢--健肢步态一致性的关键参数，为四连杆膝关节被动假肢的进一步优化提供了理论支撑.

本文的研究从动力学角度出发，对现有膝关节假肢设计技术进行了补充和完善，为基于动力学的假肢设计提供了依据.后续工作中，我们将以最大屈曲角、步态一致性等为指标，深入研究假肢动力学性能的多目标优化.

附录A

致谢本研究的第2 节实验工作由北京大学高思源协助完成，特此感谢.