引言

肢体的缺失给人们正常生活带来了诸多困难，特别是对于那些下肢肢体发生损伤的人，行走能力丧失严重限制了截肢者的活动空间和正常生活能力.据统计，我国现有下肢残疾约2412万人^[1]，其中约11.4%需要安装假肢.另外，由于每年肥胖症和糖尿病发病率的增加，糖尿病足溃疡截肢者将迅速增加^[2，3].下肢截肢者恢复行走能力的重要方式是穿戴合适的假肢.绝大多数截肢者穿戴传统被动假肢，无法为关节提供足够的输入力矩以保证基本行走步态的协调，这导致每年大约75%的下肢截肢患者二次摔倒^[4]和其它诸多慢性健康问题.

相比于被动假肢，动力假肢因可以在截肢者日常行走时提供输入净功，从而可以更好地恢复截肢者行走的生物力学功能^[5].在过去十年的动力假肢研究中，动力假肢设计和控制所实现的各种驱动方案主要分为两类，一类是采用同步带、齿轮等具有固定传动比的减速驱动方案^[6，7]；另一类则是采用滚珠丝杠或连杆结构的传动驱动装置^[8，9].滚珠丝杠或连杆结构的传动驱动装置虽然可以使得关节驱动电机尽可能远离下肢关节，进而减小假肢的摆动惯量，提升假肢响应的灵敏度.但这些驱动方案往往给假肢动力系统带来高阻抗、大的摩擦损耗和大的反射惯性.这对明确假肢动力学特性和摩擦行为以提高假肢控制性能提出了更高的要求.

动力假肢在控制上往往采用阻抗控制策略^[10]，其把动力假肢关节简化为虚拟的刚度和阻尼，从而易实现关节扭矩的顺应性控制.然而阻抗控制在步态上控制不连续、也不能实现精确的位置或力控制.动力假肢面临的主要挑战之一是为截肢者提供在户外不同地面条件（包括沙地、草地等非硬质地面等）上的行走能力.相比于传统的阻抗控制，一些基于模型的控制方法有望在动力下肢假肢上实现不同地形上的稳定行走^[11].然而，这些基于模型的假肢控制器的性能受限于建模误差、模型参数的不确定性和非线性摩擦等因素的影响^[12].明确假肢的动力学模型有助于动力假肢的设计、模拟和优化，有助于开发基于模型的高性能控制器，如逆动力学控制^[13].

文献调研表明针对动力假肢参数辨识的相关研究还比较少.Dallali^[14]对一绳驱动踝足假肢模型进行了辨识.Yang^[15]基于能量方法建立了单膝假肢理想动力学模型.然而针对动力膝-踝假肢的动力学参数识别还没有发现相关文献.在工业机器人或腿式机器人领域，机器人的动力学参数通常是利用外部扭矩传感器测量实验值与仿真计算结果的二次误差来确定的^[16].然而，在假肢领域，考虑到重量和成本，一般在关节处不会具有很精确的扭矩测量传感器.此外，关节扭矩传感器也不能完全用于识别关节执行器内部的动力学参数.对于动力假肢，驱动装置除了使用常用的固定传动比的传动结构外，通常使用连杆结构^[17，18].这些连杆结构的传动比通常是强非线性的，这给辨识工作带来了困难.

本文针对实验室新设计的具有固定传动比的膝关节和具有非线性传动比的踝关节的大腿假肢进行动力学参数辨识.通过建立假肢的动力学模型，实现了对动力膝踝假肢的系统性参数识别.本文提供了一种识别大腿假肢动力学的系统性方法，验证了库伦-粘性摩擦模型在表征假肢关节摩擦行为的性能，提高了假肢动力学模型参数的准确性.

1 动力大腿假肢结构

图1为实验室设计的动力大腿假肢实物图，由动力膝关节和动力踝关节组成.其中，动力膝关节的电机通过谐波减速器把动力输出给膝关节，编码器安装在电机转子上.动力踝关节的电机通过滚珠丝杠和连杆结构把动力输出给踝关节，其中一个相对编码器安装在电机转子上，另外一个绝对编码器则安装在踝关节处.

2 大腿假肢动力学建模

大腿假肢动力关节通常设计紧凑，其动力学参数辨识时很难通过在关节外部加装扭矩传感器来测量关节扭矩.因此，需要通过电机电流来估计关节扭矩.而电机到关节处的运动经过了减速机构，因此减速机构的传动比是动力学中的一个重要参数.

对于图1所示的膝关节，减速机构为谐波减速器，其传动比为51；对于踝关节，减速机构为滚珠丝杠加连杆结构，其减速比与设计的踝关节传动机构相关，因此，需要通过踝关节传动机构的几何关系计算出踝关节的传动比.

2.1 踝关节非线性传动比计算

如图2所示，踝关节电机转动先经过滚珠丝杠转化为螺母的上下平动，然后经过连杆结构转化为踝关节摆动.

根据滚珠丝杠螺母运动的距离与电机转速关系可得：

式中，Δd为丝杠移动距离，θ_ma为踝关节驱动电机转角，p为丝杠螺距.

滚珠丝杠螺母运动Δd后，连杆CO的长度变为：

式中，${\overline{\mathrm{C}\mathrm{O}}}_0$表示小腿和足之间的相对转角θ_a为零时的连杆初始长度.

由AC垂直于CO可得：

式中，$\overline{CO}$，$\overline{AO}$，$\overline{CO}$和$\overline{AC}$分别为图示中各连杆长度.

因$\overline{OD}$与$\overline{OC}$共线，可得：

式中，θ_a=θ_f-θ_s，θ_f和θ_s分别为脚和小腿与竖直方向夹角，∠BOF为设计的结构角度，设计值为84.75°.

根据正弦定理和余弦定理可得：

联立式（1）~式（5），可得：

式中，所设计的各连杆尺寸参数如下表1所示.

理论上，由式（6）可确定踝关节电机运动角度θ_ma经过连杆结构后与踝关节相对运动θ_a角度的关系.但式（6）中存在反正弦和反余弦函数，无法直接求出显式表达式.我们用θ_ma表示为θ_a的函数，并用隐函数表示为：

式（7）两边同时对时间求导，可得：

因此，根据传动比定义，踝关节传动比为：

根据人体正常行走步态在矢状面上的踝关节屈曲角度^[19]和已设计假肢^[7]对踝关节角度要求，我们确定在${\theta }_a=\left[-\frac{\pi }{6}，\frac{\pi }{6}\right]$范围内利用傅里叶级数拟合踝关节电机运动角度θ_ma与踝关节相对运动角度θ_a的关系.图3分别展示了6、7、8阶傅里叶级数拟合结果，从图3（b）可以看出，8阶傅里叶级数的踝关节传动比拟合误差约为2.43×10^-3，拟合已经足够精确.故拟合的傅里叶级数为：

式中，最终拟合的各系数如表2所示.

式（10）两边对时间求导，得到传动比为：

记k_a=g（θ_a），拟合关系如下图3所示.由图3（a）可以看出，踝关节传动比与踝关节转角具有很强的非线性关系，因此在动力学参数辨识时，要考虑其对辨识结果的影响.

2.2 基于拉格朗日方法的动力学建模

按照文献^[13]的辨识方式，把假肢膝关节上端固定，图1所示的动力假肢简化为图4所示的双摆杆结构.假设膝关节电机位于膝关节轴，质心位于关节原点上，踝关节电机位于小腿上，踝关节以下部分（下文叫做脚）质心位于踝关节竖直位置，同时忽略踝关节传动的拉杆和滚珠丝杠螺母的质量.

定义小腿和脚连杆与竖直方向的绝对坐标为广义坐标q=（θ_s，θ_f）^T，顺时针为正方向（从坐标轴指向连杆）.利用Euler-Lagrange方程，建立图（4）所示大腿假肢在固定坐标系下的动力学方程为：

其中，质量阵$\mathit{M}（\mathit{q}）\in {\mathbb{R}}^{2\times 2}$各元素如下：

$\begin{array}{c}{M}_{\mathrm{1,\; 1}}=I_{mk}{k}_k^2+I_k+I_s+I_m+m_s{d}_s^2+\\ m_m{d}_{ma}^2+m_f{l}_s^2+I_{ma}{k}_a^2,\\ M_{\mathrm{1,\; 2}}=m_f{l}_s{d}_{f}cos\left({\theta }_s-{\theta }_f\right)-I_{ma}{k}_a^2,\\ M_{\mathrm{2,\; 1}}=m_f{l}_s{d}_{f}cos\left({\theta }_s-{\theta }_f\right)-I_{ma}{k}_a^2,\\ M_{\mathrm{2,\; 2}}=I_f+m_f{d}_f^2+I_{ma}{k}_a^2.\end{array}$

科氏力及离心力矩阵$\mathit{C}（\mathit{q}，\dot{\mathit{q}}）\in {\mathbb{R}}^{2\times 2}$各元素为：

$\begin{array}{c}{C}_{\mathrm{1,\; 1}}=3I_{ma}{k}_a\frac{\partial k_a}{\partial {\theta }_a}{\dot{\theta }}_s-4I_{ma}{k}_a\frac{\partial k_a}{\partial {\theta }_a}{\dot{\theta }}_f,\\ C_{\mathrm{1,\; 2}}=-2I_{ma}{k}_a\frac{\partial k_a}{\partial {\theta }_a}{\dot{\theta }}_s+3I_{ma}{k}_a\frac{\partial k_a}{\partial {\theta }_a}{\dot{\theta }}_f+\\ m_f{l}_s{d}_{f}sin\left({\theta }_s-{\theta }_f\right){\dot{\theta }}_f,\\ C_{\mathrm{2,\; 1}}=-3I_{ma}{k}_a\frac{\partial k_a}{\partial {\theta }_a}{\dot{\theta }}_s+2I_{ma}{k}_a\frac{\partial k_a}{\partial {\theta }_a}{\dot{\theta }}_f-\\ m_f{l}_s{d}_{f}sin\left({\theta }_s-{\theta }_f\right){\dot{\theta }}_s,\\ C_{\mathrm{2,\; 2}}=4I_{ma}{k}_a\frac{\partial k_a}{\partial {\theta }_a}{\dot{\theta }}_s-3I_{ma}{k}_a\frac{\partial k_a}{\partial {\theta }_a}{\dot{\theta }}_f.\end{array}$

重力阵$\mathit{G}（\mathit{q}）\in {\mathbb{R}}^{2\times 1}$各元素如下：

$\begin{array}{c}{G}_{\mathrm{1,\; 1}}=\left(m_s{d}_s+m_m{d}_{ma}+m_f{l}_s\right)gsin{\theta }_s,\\ G_{\mathrm{2,\; 1}}=m_f{d}_{f}gsin{\theta }_f.\end{array}$

广义力阵$\mathit{T}（\mathit{q}，\dot{\mathit{q}}）\in {\mathbb{R}}^{2\times 1}$，各元素如下：

$\begin{array}{c}{T}_{\mathrm{1,\; 1}}=T_{ms}-T_{\tau s}-\left(T_{mf}-T_{\tau f}\right),\\ T_{\mathrm{2,\; 1}}=T_{mf}-T_{\tau f}.\end{array}$

以上各式中，l_s为小腿长度，d_s，d_f，d_ma分别为小腿质心到膝关节转动中心距离、脚质心到踝关节转动中心距离以及踝关节电机质心到膝关节转动中心的距离，m_s，m_f分别为小腿（包含电机定子）和脚的质量，m_m为电机定子和转子总质量，m_mk，m_ma分别为膝关节和踝关节的电机转子及包含相连接转动轴质量，I_mk，I_ma分别为膝关节和踝关节的电机转子包含相连接转动轴的转动惯量，I_s，I_f分别为小腿和脚相对于其自身质心的转动惯量，I_k，I_a分别为膝关节和踝关节经减速器传动后的相关部件转动惯量，θ_mk，θ_ma分别为膝关节和踝关节电机的转动角度，k_k，k_a分别为膝关节和踝关节传动比，T_ms，T_mf分别为膝关节和踝关节电机扭矩，T_τs，T_τf分别为膝关节和踝关节的关节摩擦力损失.

本文关节摩擦力采用常见的库伦-粘性摩擦力模型形式表征：

$T_{\tau s}=f_{cs}\mathrm{s}\mathrm{g}\mathrm{n}\left({\dot{\theta }}_s\right)+f_{vs}{\dot{\theta }}_s$

式中，f_cs，f_vs，f_cf，f_vf分别为膝关节库伦粘性摩擦系数以及踝关节库伦和粘性摩擦系数.

假肢动力学方程中共有17个参数，其中小腿长度l_s采用游标卡尺精确测定为0.3615 m，小腿质量m_s、脚的质量m_f和电机定子转子的总质量m_m采用电子秤进行精确测定，分别为1.624kg、0.442kg和0.4112kg.除了4个摩擦参数外的剩余9个参数可通过CAD软件测量估计其大小，以便跟辨识结果对比，测量结果如表3所示.

根据参数的可辨识性，参数I_mk，I_k，I_s和I_m可合并为一个待辨识参数，即$I_1=I_{mk}{k}_k^2+I_k+I_s+I_m$.因此，动力学方程中17个参数简化为包括4个摩擦参数在内的共10个待辨识参数，分别为d_ma，d_f，d_s，I₁，I_ma，I_f，f_cs，f_cf，f_vs和f_vf.

3 参数辨识

3.1 辨识实验方案

参数辨识实验如图5所示，膝关节上端通过固定杆悬挂于隔震台上，假肢上关节轴和连杆中心贴上不同颜色的标记点，利用相机拍摄记录假肢运动过程，拍摄帧率约为120Hz.同时假肢控制器记录电机上编码器数据、电机的电流数据和关节编码器数据，并通过无线串口发送给PC记录保存，记录频率约为200Hz.

实验时膝关节电机和踝关节都采用位置控制模式，分别利用激励基频是2.5 Hz的5阶傅里叶级数激励信号测量3次，每次测量约2分钟，设置为：

其中，幅值A和B根据关节的机械限位进行调整.对于膝关节A取0.2，B取0.1，对于踝关节A取0.1，B取0.05.

3.2 参数辨识方法

本文参数辨识方法采用普遍使用的粒子群优化算法（PSO）.相比于其它辨识算法，粒子群优化算法不需要求解梯度，算法更为简单，可以针对非光滑摩擦进行参数辨识.其还具有精度高、收敛速度快等优点，在参数辨识领域得到广泛应用.例如，Wang^[20]将粒子群优化算法与不确定的速度和位移的系数修正相结合，构建了基于自适应信号修正的辨识方法.Chen^[21]提出多段摩擦系数反演方法，并利用基于粒子群优化（PSO）的算法来计算了页岩气套管作业中的多段摩擦系数.本文的大腿假肢动力学参数辨识，也使用粒子群优化算法（PSO）进行参数辨识.首先通过位置激励假肢双关节获得运动信号，然后通过膝关节和踝关节驱动器获得关节电机电流数据，得到电机实测扭矩.电机扭矩计算为：

式中，K为电机扭矩常数，通过实验测量得到，取值为0.1031Nm/A.

本文利用粒子群优化算法辨识的目标是寻找最优的一组动力学参数以最小化仿真所求关节扭矩与实测关节扭矩之间的误差，即：

式中，A₁···A₁₀为待辨识参数，T_mk，simT_ma，sim为辨识仿真所求膝关节和踝关节电机扭矩，T_mk，meaT_ma，mea为实测膝关节和踝关节电机扭矩.

3.3 实验数据处理

首先利用MATLAB提取相机拍摄视频中每个标记点像素中心位置坐标.相机拍摄每帧视频原点为左上角，通过RGB值提取标记点的中点，依次提取四个标记点的像素坐标中心点位置.控制器记录的原始数据为电机上编码器数据、电机的电流数据和踝关节编码器数据.记录的编码器数据原始值为当前编码器刻度值，需要转换为弧度制.记录的电机电流数据需要乘以扭矩常数转换为电机扭矩.

由于实验只测得了膝关节和踝关节的角度信号，需要构造角速度和角加速度.由于多项式滑移窗口方法可以方便的重构出角度、角速度和角加速度信号，并且具有重构效果好，抗噪能力强等优点.故本文参考文献^[20，22]中的方法来获得辨识所需要的数据.为了拟合整个时程上多个周期的角度信号，通过调整Legendre多项式的阶数和拟合时间窗口的长度来提高拟合精度.本文设计的Legendre多项式拟合阶数为20，拟合窗口长度为1 s，相应的速度和加速度通过Legendre多项式求导得出.最终得到的小腿和脚连杆与竖直方向的绝对角度、绝对角速度、绝对角加速度和关节扭矩数据如图6所示.

3.4 辨识结果

最终得到的动力学参数辨识结果如表4所示.与表3的动力学参数CAD软件测量结果对比可知，除了踝关节电机到膝关节转轴距离辨识的偏差比较大之外，大部分参数比较接近.

为了比较最终的辨识效果，我们分别把辨识后的动力学参数数据和CAD测量的动力学参数数据代入动力学方程，求出电机扭矩来对比验证辨识结果.如图7所示是假肢膝关节参数辨识结果的验证图，其中上图的黑线是实验测得的膝关节电机扭矩，蓝色线是CAD测量的动力学参数代入动力学方程求得的膝关节电机扭矩，红线是辨识后的动力学参数代入动力学方程求得的膝关节电机扭矩.从图7的下图中数据可以得出，由CAD测量的动力学参数计算出的电机扭矩与实测电机扭矩值之间的均方根误差RMSE为0.0076 Nm，而辨识后参数计算出的电机扭矩与实测电机扭矩值的均方根误差RMSE为0.00049 Nm，为辨识前膝关节电机扭矩均方根误差的6.45%，辨识后膝关节电机扭矩精度明显提高.

如图8所示，其中上图的黑线是实验测得的踝关节电机扭矩，蓝色线是CAD测量的动力学参数代入动力学方程求得的踝关节电机扭矩，红线是辨识后的动力学参数代入动力学方程求得的踝关节电机扭矩.从图8的下图中数据可以得出，对于假肢踝关节，由CAD测量的动力学参数计算出的踝关节电机扭矩与实测踝关节电机扭矩的均方根误差RMSE为0.0012 Nm，而辨识后参数计算出的电机扭矩均方根误差RMSE为0.00023 Nm，为辨识前踝关节电机扭矩均方根误差的19.17%，辨识后的精度也有所提高，但精度提高并没有膝关节的多.总的来说，辨识的动力学参数结果比用CAD评估的参数更加准确，更能表示系统的动力学特性.

4 结论

本文主要研究实验室新设计的具有固定传动比膝关节和具有非线性传动比踝关节的大腿假肢动力学参数辨识问题.通过研究其动力学参数的CAD测量与实际辨识后参数所求电机扭矩的误差，明确了基于库伦-粘性摩擦的大腿假肢动力学参数的精度，进而为提高基于模型的假肢控制精度奠定了理论研究基础.通过本文的研究可以得出如下结论：

（1）通过对基于库伦-粘性摩擦模型的大腿假肢动力学进行参数辨识，可以很大程度上提高动力学预测电机扭矩的精度.

（2）对于具有不同传动链的不同动力假肢关节，多级传动链的累积误差会造成系统的不确定性放大，最终使动力学参数的辨识精度变差.

参考文献