人口老年化、交通意外和脑中风导致越来越多的人出现下肢运动功能障碍。下肢外骨骼机器人(Lower Limb Exoskeleton Robot, LLER) 能够帮助下肢伤残、偏瘫、有运动功能障碍的患者和身体虚弱的老年人进行康复训练和助力行走。因此，备受业内专家学者和健康医疗机构关注^[1-2]。刚出台的《“十四五”机器人产业发展规划》和《“十四五”医疗装备产业发展规划》对于外骨骼机器人相关技术的研究给予高度重视和大力支持。下肢外骨骼机器人是一个多变量、强耦合的非线性系统，为了使其动作更协调柔顺、步态更自然平稳、轨迹更平滑精准，建立人机协同系统精确的动力学模型和参数辨识是首要和关键步骤^[3-4]。目前，国内外学者提出的机器人参数辨识方法有：(1) 解体测量计算法^[5]：无需运动数据采集和算法设计，但不易精确测量质心位置，且拆解后无法考虑关节摩擦力的影响。(2) CAD辅助法^[6]：惯性参数由软件或理论公式自动计算得到的理想值，但在实际加工及装配过程中会产生误差。(3) 线性化整体辨识法^[7-8]：先对动力学模型进行线性化处理，得到线性回归矩阵和惯性参数矩阵，再驱动机器人运动，一次性辨识出所有参数，但计算量较大，存在重复辨识和误差累积等问题。(4) 多次辨识法^[9]：对整体辨识法加以改进，将机器人按关节分组后进行多次驱动和参数辨识，但忽略了关节摩擦的影响。王靖森等^[10]和席万强等^[11]考虑摩擦，对六自由度机器人进行建模和辨识，提高了模型精度。(5) 内部传感器法：采用编码器测量关节角度，根据驱动器电流计算转矩，但转矩常数随着温度和寿命而变化存在误差。耿令波^[12]采用测驱动电流和关节力矩的综合辨识法，实现惯性参数和摩擦解耦。(6) 外部传感器法^[13]：通过高精度视觉系统捕获运动数据，通过安装在基座的力传感器测力矩，但并不能得到所有的惯性参数，且机器人基座需要多次定位，位置坐标变换计算复杂。(7) 负载辨识法：陈友东等^[14]仅驱动机器人的部分关节减少关节耦合误差，运行激励轨迹，在负载变化时进行动力学参数辨识。陈柏等^[15]通过优化激励轨迹对六自由度工业机器人负载动力学模型进行辨识，再补偿已知机器人本体参数模型，得到本体和负载结合的动力学模型。

针对下肢外骨骼机器人人机协同系统，本文借鉴上述工业机器人的参数辨识方法，并进行综合与改进，提出一种将外骨骼参数辨识实验方法与人体下肢参数计算方法相结合的动力学参数辨识方案。首先，为了考虑摩擦力对参数辨识的影响，在文献[16]建立的二连杆结构机器人动力学模型中加入摩擦力项，得到转矩表达式和一组待辨识的最小惯性参数。其次，提出一种静态与动态结合的实验辨识方法：先多次驱动单个关节，使大腿和小腿分别运动到不同位置后保持静止，使转矩表达式中的角速度项和角加速度项变为零，计算更简单；再以最简单的正弦函数作为最优激励轨迹驱动下肢连续摆动；并利用MATLAB进行在线数据采集、计算和曲线绘制，依据曲线特征和最小二乘法完成外骨骼参数辨识。再次，选取余弦函数作为激励轨迹进行参数辨识结果验证。然后，依据人体生理特征对下肢参数开展近似计算，并综合实验辨识结果建立人机协同系统精确的动力学模型。最后，采用基于模型上界的滑模控制，并引入低通滤波器开展步态轨迹跟踪研究。

1 动力学模型

笔者课题组研制了一款如图1(a~b) 所示的下肢外骨骼服^[16]。当使用者穿戴下肢外骨骼服，单腿支撑且另一条腿摆动行走时，为了便于分析和计算，将脚和小腿作为一个整体，下肢的摆动近似简化为二连杆在矢状面XOY内的运动，如图1(c) 所示。以摆动腿的髋关节O为坐标原点，大腿绕其转动，OA长度为 $ {L_{\text{1}}} $ ，质量为 $ {m_{\text{1}}} $ ，质心 $ {C_{\text{1}}} $ 距离O点 $ {l_{\text{1}}} $ ，与竖直方向的夹角为 ${\theta _{\text{1}}}，{\theta _{\text{1}}} \in [ - {30^\circ },{120^\circ }]$ 。小腿绕膝关节A转动，长度为 $ {L_{\text{2}}} $ ，质量为 $ {m_{\text{2}}} $ ，质心 $ {C_{\text{2}}} $ 距离A点 $ {l_{\text{2}}} $ ，与大腿延长线方向的夹角为 ${\theta _{\text{2}}}，{\theta _{\text{2}}} \in [ - {120^\circ },{0^\circ }]$ 。大腿绕通过髋关节O并垂直于XOY平面的Z轴转动的转动惯量为 $ {I_{\text{1}}} $ ，小腿绕通过膝关节A并平行于Z轴转动的转动惯量为 $ {I_{\text{2}}} $ 。

一般，一个关节机器人的动态性能可由如式(1)二阶非线性微分方程描述^[17]。

$ {\boldsymbol{M}}(\theta )\ddot {\boldsymbol{ \theta }} + {\boldsymbol{C}}(\theta ,\dot \theta )\dot {\boldsymbol{ \theta }} + {\boldsymbol{G}}(\theta ) + {\boldsymbol{F}}(\dot \theta ) = {\boldsymbol{T}} $

(1)

式中： $ {\boldsymbol{\theta}} ，\dot{{\boldsymbol{\theta }}}，\ddot{{\boldsymbol{\theta}} }\in {R}^{n\times 1} $ 分别为关节的角位置、角速度和角加速度， $ {\boldsymbol{M}}(\theta ) \in {R^{n \times n}} $ 为机器人的正定惯性矩阵， $ {\boldsymbol{C}}(\theta ,\dot \theta ) \in {R^{n \times 1}} $ 为离心力和哥氏力矩阵， $ G(\theta ) \in {R^{n \times 1}} $ 为重力项， $ {\boldsymbol{F}}(\dot \theta ) \in {R^{n \times 1}} $ 为摩擦力项，若只考虑黏滞摩擦和库仑摩擦， ${\boldsymbol{F}}(\dot \theta ) = {{\boldsymbol{f}}_{\rm{v}}}\dot{\boldsymbol{ \theta }} + {{\boldsymbol{f}}_{\rm{c}}}{{\rm{sgn}}} (\dot \theta )$ ， $ {{\boldsymbol{f}}_{\rm{v}}} \in {R^{n \times n}} $ 为黏滞摩擦系数矩阵， $ {{\boldsymbol{f}}_{\rm{c}}} \in {R^{n \times n}} $ 为库伦摩擦系数矩阵， $ {\boldsymbol{T}} \in {R^{n \times 1}} $ 为关节驱动力矩。重力加速度 $g取9.8\;{\rm{m}}/ {\rm{s}}^2$ 。

因此，二连杆结构的机器人动力学方程为

$ \begin{split} & \left[ {\begin{array}{*{20}{c}} {{M_{11}}}&{{M_{12}}} \\ {{M_{21}}}&{{M_{22}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{{\ddot \theta }_1}} \\ {{{\ddot \theta }_2}} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {{C_{11}}}&{{C_{12}}} \\ {{C_{21}}}&{{C_{22}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{{\dot \theta }_1}} \\ {{{\dot \theta }_2}} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {{G_1}} \\ {{G_2}} \end{array}} \right] +\\&\qquad \left[ {\begin{array}{*{20}{c}} {F({{\dot \theta }_1}) } \\ {F({{\dot \theta }_2}) } \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{T_1}} \\ {{T_2}} \end{array}} \right]\\[-15pt] \end{split} $

(2)

式中，矩阵中各个元素表示为^[16,18]

$ \left\{ \begin{gathered} {M_{11}} = {I_1} + {I_2} + 2{m_2}{L_1}{l_2}\cos {\theta _2} + {m_1}l_1^2 + {m_2}L_1^2 + {m_2}l_2^2 \\ {M_{12}} = {M_{21}} = {I_2} + {m_2}l_2^2 + {m_2}{L_1}{l_2}\cos {\theta _2} \\ {M_{22}} = {I_2} + {m_2}l_2^2 \\ {C_{11}} = - 2{m_2}{L_1}{l_2}\sin {\theta _2}{{\dot \theta }_2} \\ {C_{12}} = - {m_2}{L_1}{l_2}\sin {\theta _2}{{\dot \theta }_2} \\ {C_{21}} = {m_2}{L_1}{l_2}\sin {\theta _2}{{\dot \theta }_1} \\ {C_{22}} = 0 \\ {G_1} = ({m_1}{l_1} + {m_2}{L_1}) g\sin {\theta _1} + {m_2}{l_2}g\sin ({\theta _1} + {\theta _2}) \\ {G_2} = {m_2}{l_2}g\sin ({\theta _1} + {\theta _2}) \\ F({{\dot \theta }_1}) = {f_{{\rm{v}}1}}{{\dot \theta }_1} + {f_{{\rm{c}}1}}{{\rm{sgn}}} ({{\dot \theta }_1}) \\ F({{\dot \theta }_2}) = {f_{{\rm{v}}2}}{{\dot \theta }_2} + {f_{{\rm{c}}2}}{{\rm{sgn}}} ({{\dot \theta }_2}) \\ \end{gathered} \right. $

(3)

式 (3) 中，共有12个未知参数，为了计算简便，重新定义5个新参数作为一组动力学最小惯性参数：

$ \left\{ \begin{array}{l} {X_1} = {I_1} + {m_1}l_1^2 + {m_2}L_1^2 + {I_2} + {m_2}l_2^2 \\ {X_2} = {I_2} + {m_2}l_2^2 \\ {X_3} = {m_2}{l_2}{L_1} \\ {X_4} = {m_1}{l_1} + {m_2}{L_1} \\ {X_5} = {m_2}{l_2} \end{array} \right. $

(4)

则，式(2) 可以写成

$ \begin{split} & \left[ {\begin{array}{*{20}{c}} {{X_1} + 2{X_3}\cos {\theta _2}}&{{X_2} + {X_3}\cos {\theta _2}} \\ {{X_2} + {X_3}\cos {\theta _2}}&{{X_2}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{{\ddot \theta }_1}} \\ {{{\ddot \theta }_2}} \end{array}} \right] + \\& \qquad {X_3}\sin {\theta _2}\left[ {\begin{array}{*{20}{c}} { - 2{{\dot \theta }_2}}&{ - {{\dot \theta }_2}} \\ {{{\dot \theta }_1}}&0 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{{\dot \theta }_1}} \\ {{{\dot \theta }_2}} \end{array}} \right] +\\& \qquad g\left[ {\begin{array}{*{20}{c}} {{X_4}\sin {\theta _1} + {X_5}\sin ({\theta _1} + {\theta _2}) } \\ {{X_5}\sin ({\theta _1} + {\theta _2}) } \end{array}} \right] +\\& \qquad \left[ {\begin{array}{*{20}{c}} {{f_{{\rm{v}}1}}{{\dot \theta }_1} + {f_{{\rm{c}}1}}{{\rm{sgn}}} ({{\dot \theta }_1}) } \\ {{f_{{\rm{v}}2}}{{\dot \theta }_2} + {f_{{\rm{c}}2}}{{\rm{sgn}}} ({{\dot \theta }_2}) } \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{T_1}} \\ {{T_2}} \end{array}} \right] \end{split}$

(5)

故，髋关节和膝关节的驱动转矩分别为

$ \begin{split} {T_1} =& ({X_1} + 2{X_3}\cos {\theta _2}) {{\ddot \theta }_1} + ({X_2} + {X_3}\cos {\theta _2}) {{\ddot \theta }_2} - \\& 2{X_3}\sin {\theta _2}{{\dot \theta }_2}{{\dot \theta }_1} - {X_3}\sin {\theta _2}{{\dot \theta }_1}{{\dot \theta }_2} + g[{X_4}\sin {\theta _1} + \\& {X_5}\sin ({\theta _1} + {\theta _2}) ] + {f_{{\rm{v}}1}}{{\dot \theta }_1} + {f_{{\rm{c}}1}}{{\rm{sgn}}} ({{\dot \theta }_1}) \end{split} $

(6)

$ \begin{split} {T_2} =& ({X_2} + {X_3}\cos {\theta _2}) {{\ddot \theta }_1} + {X_2}{{\ddot \theta }_2} + {X_3}\sin {\theta _2}\dot \theta _1^2 + \\& g{X_5}\sin ({\theta _1} + {\theta _2}) + {f_{{\rm{v}}2}}{{\dot \theta }_2} + {f_{{\rm{c}}2}}{{\rm{sgn}}} ({{\dot \theta }_2}) \end{split} $

(7)

如图1(b) 所示，安装在关节处的成套减速电机(包括编码器、直流伺服电机和行星减速齿轮)较重，且足部踏板质量不可忽略，导致各连杆的质心不在几何中心位置，即式(3) 中的 $ {l_1} $ 和 $ {l_2} $ 难以准确描述。同时即使已知各连杆的质量，而 $ {I_1} $ 和 $ {I_2} $ 也未知，且摩擦力项 $ {\boldsymbol{F}}(\dot \theta ) $ 无法通过传感器检测。显然，仅通过测量与计算无法建立精确的系统模型。因此，进行系统动力学参数辨识十分必要。

2 参数辨识实验 2.1 实验方案及准备

机器人动力学参数辨识是一个复杂的系统工程，包括动力学建模、激励轨迹的选取、实验设计、数据采集与处理、参数估计、模型验证等环节^[12，19]。本文动力学参数辨识与实验验证方案如图2所示。其中，X表示外骨骼动力学最小惯性参数集，X_load表示人体下肢动力学最小惯性参数集，T表示采用滑模控制输出的关节驱动转矩，T'表示经过低通滤波后输出的关节驱动转矩。 $e,\dot e,\ddot e $ 表示关节角度跟踪误差及误差一阶和二阶导数。

如图3所示，实验前用绳子将外骨骼悬空吊起，用支架将其固定，使其躯干与双腿都处于竖直静止状态。控制右腿髋关节和膝关节处的电机来驱动大腿和小腿(含足部) 运动，模拟外骨骼左腿竖直支撑且右腿摆动行走的过程。

2.2 静态实验

对于髋关节和膝关节，定义关节离开初始位置为反向运动，返回初始位置为正向运动。轮流让其中一个关节在初始位置固定不动，另一个关节运动到不同位姿后保持静止，此刻角速度和角加速度都为零，转矩表达式中仅剩重力项和摩擦力项，易于参数辨识。

图4所示为关节转矩与角度的关系。

第1步：让右腿竖直且膝关节保持静止， $ {\theta _2} = {0^\circ } $ ，髋关节初始角度 $ {\theta _{10}} = {0^\circ } $ 。由上位机给电机控制器发送髋关节角度期望信号 $ {\theta _{1{\text{d}}0}} = {30^\circ } $ ，控制器在简单比例控制作用下输出控制电压 $ {u_1} = {k_{p1}}({\theta _{1{\text{d}}0}} - {\theta _{10}}) $ 驱动右大腿摆动到 $ {30^\circ } $ 位置后停止，记录该位置的关节驱动转矩 $ {T_{10}} $ 。每次将 $ {\theta _{{\text{1d}}i}} $ 增加 $ {2^\circ } $ ，直到 $ {\theta _{1{\text{d}}30}} = {90^\circ } $ 为止，记录此刻关节处驱动转矩 $ {T_{130}} $ ，得到 $ {\theta _1} $ 反向增加时髋关节驱动转矩曲线 $ T_{{\text{1}}i}^{{ - }} $ ，如图4(a) 所示。

髋关节反向运动到指定位置并保持静止时， ${\dot{{ \theta }}_{{\text{1d}}i}} = {\ddot{{ \theta }}_{{\text{1d}}i}} = 0$ ， $ i = 1\sim 30 $ 。由式(6) 得，

$ T_{1i}^ - = g({X_4} + {X_5}) \sin {\theta _{{\text{1}}i}} + {f_{{\rm{c}}1}} $

(8)

第2步：由上位机给电机控制器发送髋关节角度期望信号 $ {\theta _{1{\text{d}}31}} = {88^\circ } $ ，电机驱动右大腿正向摆动到 $ {88^\circ } $ 位置后停止，记录此刻关节处驱动转矩 $ {T_{131}} $ 。依次将 $ {\theta _{1{\text{d}}i}} $ 减少 $ {2^\circ } $ ，直到 $ {\theta _{1{\text{d}}60}} = {30^\circ } $ ，记录此刻关节处驱动转矩 $ {T_{160}} $ ，得到 $ {\theta _1} $ 正向减小时髋关节驱动转矩曲线 $ T_{1i}^ + $ ，如图4(a) 所示。

髋关节正向运动到指定位置并保持静止时， ${ \dot {{ \theta }}_{1{\text{d}}i}} ={\ddot {{\theta }}_{1{\text{d}}i}} = 0$ ， $ i = 31\sim 60 $ 。由式(6) 得，

$ T_{1i}^ + = g({X_4} + {X_5}) \sin {\theta _{1i}} - {f_{{\rm{c}}1}} $

(9)

根据式(8~9) 得，

$ {T_{{\text{1}}i}} = \frac{1}{2}(T_{1i}^ + + T_{{\text{1}}i}^ - ) = ({X_4} + {X_5}) \sin {\theta _{1i}} $

(10)

$ {f_{{\rm{c}}1}} = \frac{1}{2}(T_{1i}^ - - T_{1i}^ + ) $

(11)

根据式(10) 得，

$ {X_4} = \frac{{T_{1i}^ + + T_{1i}^ - }}{{2g\sin {\theta _{1i}}}} - {X_5} $

(12)

根据式(11) ，结合实验测得的60组数据，求均值得到 $ {f_{{\rm{c}}1}} = - 2.415 $ 。

第3步：让右腿重新竖直且髋关节保持静止， $ {\theta _1} = {0^\circ } $ ，膝关节初始角度 $ {\theta _{20}} = {0^\circ } $ 。由上位机给电机控制器发送膝关节角度期望信号 $ {\theta _{2{\text{d}}1}} = - {10^\circ } $ ，控制器在简单比例控制作用下输出控制电压 $ {u_2} = {k_{p2}}({\theta _{2{\text{d}}1}} - {\theta _{20}}) $ 驱动右小腿摆动到 $ - {10^\circ } $ 位置后停止，记录该位置的关节驱动转矩 $ {T_{20}} $ 。每次将 $ {\theta _{2{\text{d}}i}} $ 减少 $ {2^\circ } $ ，直到 $ {\theta _{2{\text{d}}30}} = - {70^\circ } $ 为止，记录此刻关节处驱动转矩 $ {T_{230}} $ ，得到 $ {\theta _2} $ 反向减小时膝关节转矩曲线 $ T_{2i}^ - $ ，如图4(b) 所示。

膝关节反向运动到指定位置并保持静止时， ${ \dot {{ \theta }}_{2{\text{d}}i}} = {\ddot{{ \theta }}_{2{\text{d}}i}} = 0$ ， $ i = 1\sim 30 $ 。由式(7) 得，

$ T_{2i}^ - = g{X_5}\sin {\theta _{2i}} - {f_{{\rm{c}}2}} $

(13)

第4步：由上位机给电机控制器发送膝关节角度期望信号 $ {\theta _{2{\text{d}}31}} = - {68^\circ } $ ，电机驱动右小腿反向摆动到 $ - {68^\circ } $ 位置后停止，记录此刻关节处驱动转矩 $ {T_{231}} $ 。依次将 $ {\theta _{2{\text{d}}i}} $ 增加 $ {2^\circ } $ ，直到 $ {\theta _{2{\text{d}}60}} = - {10^\circ } $ ，记录此刻关节处驱动转矩 $ {T_{260}} $ ，得到 $ {\theta _2} $ 正向增加时膝关节驱动转矩 $ T_{2i}^ + $ 曲线，如图4(b) 所示。

膝关节正向运动到指定位置并保持静止时， ${\dot{\boldsymbol{ \theta }}_{2{\text{d}}i}} = {\ddot {\boldsymbol{\theta }}_{2{\text{d}}i}} = 0$ ， $ i = 31\sim 60 $ 。由式(7) 得，

$ T_{2i}^ + = g{X_5}\sin {\theta _{2i}} + {f_{{\rm{c}}2}} $

(14)

根据式(13~14) 得，

$ {{{T}}_{{{2i}}}} = \frac{1}{2}({{T}}_{{{2i}}}^{{ + }} + {{T}}_{{{2i}}}^{{ - }}) = g{X_5}\sin {\theta _{2i}} $

(15)

$ {f_{c2}} = \frac{1}{2}(T_{2i}^ + - T_{2i}^ - ) $

(16)

根据公式(15) 得，

$ {X_5} = \frac{{T_{2i}^ + + T_{2i}^ - }}{{2g\sin {\theta _{2i}}}} $

(17)

根据式(16) ，结合实验测得的60组数据，求均值得到 $ {f_{{\rm{c}}2}} = - 1.521 $ 。根据式(17) ，结合实验测得的60组数据，求均值得到 $ {X_5} = 0.513 $ 。根据式(12) ，结合实验测得的60组数据，求均值得到 $ {X_4} = 1.871 $ 。已知外骨骼大腿的长度 $ {L_1} = 0.5 $ ，由式(4) 得 $ {X_3} = {X_5}{L_1} $ ，计算得出 $ {X_3} = 0.257 $ 。

2.3 动态实验

对于髋关节和膝关节，轮流让其中一个关节在初始位置固定不动，另一个关节先匀速运动使其角加速度为零，绘制摩擦曲线可求黏滞摩擦系数。再通过正弦函数激励轨迹让关节周期性摆动求转动惯量参数。

第1步：让右侧小腿与大腿保持直线，即膝关节不运动， $ {\theta _2} = 0 $ ， $ {{{\dot \theta }}_2} = 0 $ ， $ {{{\ddot \theta }}_2} = 0 $ 。由上位机给电机控制器发送髋关节角度期望信号 $ {\theta _{1d}}(t) = kt $ ，通过简单比例控制算法使得电机驱动大腿带动小腿和足部一起摆动。髋关节匀速运动， $ {\dot \theta _1} = k \in [ - 25,25] $ ， $ {\ddot \theta _1} = 0 $ 。根据式(6) 得，

$ {T_1} = g\sin {\theta _1}({X_4} + {X_5}) + {f_{{\rm{v}}1}}{\dot \theta _1} + {f_{{\rm{c}}1}}{{\rm{sgn}}} ({\dot \theta _1}) $

(18)

$ {T_1} - g\sin {\theta _1}({X_4} + {X_5}) = {f_{{\rm{v}}1}}{\dot \theta _1} + {f_{{\rm{c}}1}}{{\rm{sgn}}} ({\dot \theta _1}) = F({\dot \theta _1}) $

(19)

式中： $ {T_1} $ 、 $ {X_4} $ 、 $ {X_5} $ 、 $ {f_{{\rm{c}}1}} $ 均为已知，计算并绘制 $ F({\dot \theta _1}) $ 如图5(a) 所示。根据黏滞摩擦和库伦摩擦特性，将图5(a) 中实测数据拟合成黑色折线，利用最小二乘法进行估计得到 $ {f_{{\rm{v}}1}} = - 0.062 $ ，扰动 $ {f_1} = - 1.796 $ 。

第2步：让右侧大腿重新保持竖直和静止，即髋关节不运动， $ {\theta _1} = 0 $ ， $ {{{\dot \theta }}_1} = 0 $ ， $ {{{\ddot \theta }}_1} = 0 $ 。由上位机给电机控制器发送膝关节角度期望信号 $ {\theta _{{\text{2d}}}}(t) = kt $ ，通过简单比例控制算法使得电机驱动小腿带动足部一起摆动。膝关节匀速运动， $ {\dot \theta _2} = k \in [ - 25,25] $ ， $ {\ddot \theta _2} = 0 $ 。

根据式(7) 得，

$ {T_2} = g{X_5}\sin {\theta _2} + {f_{{\rm{v}}2}}{\dot \theta _2} + {f_{{\rm{c}}2}}{{\rm{sgn}}} ({\dot \theta _2}) $

(20)

$ {T_2} - g{X_5}\sin {\theta _2} = {f_{{\rm{v}}2}}{\dot \theta _2} + {f_{{\rm{c}}2}}{{\rm{sgn}}} ({\dot \theta _2}) = F({\dot \theta _2}) $

(21)

式中： $ {T_2} $ 、 $ {X_5} $ 、 $ {f_{{\rm{c}}2}} $ 均为已知，计算并绘制 $ F({\dot \theta _2}) $ 如图5(b) 所示。根据黏滞摩擦和库伦摩擦特性，将图5(b) 中实测数据拟合成黑色折线，利用最小二乘法进行估计得到 $ {f_{v2}} = - 0.503 $ 。

第3步：让右侧大腿重新保持竖直和静止，即髋关节不运动 $ {\theta _1} = 0 $ ， $ {{{\dot \theta }}_1} = 0 $ ， $ {{{\ddot \theta }}_1} = 0 $ 。膝关节激励轨迹选取： ${\theta _{{\text{2d}}}}(t) = - {30^\circ }{\text{ + }}{20^\circ }\sin (2\text{π} ft)$ ，通过简单比例控制算法使得电机驱动小腿和足部一起摆动。其中 $ f $ 为0.2 Hz，以频率100 Hz采集 $ {T_2} $ 和 $ {\theta _2} $ ，计算出每个采样点的角速度 $ {\dot \theta _2} $ 和角加速度 $ {\ddot \theta _2} $ 。

根据式(7) 得，

$ {T_2} = {X_2}{\ddot \theta _2} + g{X_5}\sin {\theta _2} + {f_{{\rm{v}}2}}{\dot \theta _2} + {f_{{\rm{c}}2}}{{\rm{sgn}}} ({\dot \theta _2}) $

(22)

$ {X_2} = \frac{{{T_2} - g{X_5}\sin {\theta _2} - {f_{{\rm{v}}2}}{{\dot \theta }_2} - {f_{{\rm{c}}2}}{{\rm{sgn}}} ({{\dot \theta }_2}) }}{{{{\ddot \theta }_2}}} $

(23)

式中： $ {T_2} $ 、 $ {X_5} $ 、 $ {f_{{\rm{c}}2}} $ 、 $ {f_{{\rm{v}}2}} $ 均为已知，利用最小二乘法辨识可得 $ {X_2} = 2.768 $ 。

第4步：再次让小腿与大腿保持直线，即膝关节不运动， $ {\theta _2} = 0 $ ， $ {{{\dot \theta }}_2} = 0 $ ， $ {{{\ddot \theta }}_2} = 0 $ 。选取髋关节激励轨迹 ${\theta _{{\text{1d}}}}(t) = - {30^\circ }{\text{ + }}{20^\circ }\sin (2\text{π} ft)$ ， $ f $ 为0.2 Hz，以频率100 Hz采集 $ {T_1} $ 和 $ {\theta _1} $ ，计算出每个采样点的角速度 $ {\dot \theta _1} $ 和角加速度 $ {\ddot \theta _1} $ 。根据式(6) 得，

$ \begin{split} &{T_1} = ({X_1} + 2{X_3}) {\ddot \theta _1} + g\sin {\theta _1}({X_4} + {X_5}) + {f_{{\rm{v}}1}}{\dot \theta _1} + {f_{{\rm{c}}1}}{{\rm{sgn}}} ({\dot \theta _1}) \\&{} \end{split}$

(24)

$ {X_1}{\text{ = }}\frac{{{T_1} - g\sin {\theta _1}({X_4} + {X_5}) - {f_{{\rm{v}}1}}{{\dot \theta }_1} - {f_{{\rm{c}}1}}{{\rm{sgn}}} ({{\dot \theta }_1}) }}{{{{\ddot \theta }_1}}} - 2{X_3} $

(25)

式中： $ {T_1} $ 、 $ {X_3} $ 、 $ {X_4} $ 、 $ {X_5} $ 、 $ {f_{{\rm{c}}1}} $ 、 $ {f_{{\rm{v}}1}} $ 均为已知，利用最小二乘法辨识可得 $ {X_1} = 9.506 $ 。

综上所述，参数辨识结果如表1所示。

2.4 辨识结果验证

根据髋关节和膝关节实际角度的变化范围： $ {\theta _{\text{1}}} \in [ - {30^\circ },{120^\circ }] $ ， $ {\theta _{\text{2}}} \in [ - {120^\circ },{0^\circ }] $ ，选取验证激励轨迹：

$ \left\{ \begin{gathered} {\theta _{1{\text{d}}}}(t) = {45^\circ } - {75^\circ }\cos (2\text{π} ft) \\ {\theta _{2{\text{d}}}}(t) = - {60^\circ } + {60^\circ }\cos (2\text{π} ft) \\ \end{gathered} \right. $

(26)

式中， $ f = 1{\text{ Hz}} $ 。将式(26) 与参数辨识结果代入式(6~7) 计算，并绘制曲线。在式(26) 激励轨迹驱动下，外骨骼髋、膝关节同时运动，根据电机驱动器电流和已知的转矩系数，计算得到关节驱动转矩实际值，并绘制曲线。如图6所示，将转矩计算值与实际测量值进行对比分析发现，关节转矩计算值(红色曲线)都比实际测量值(蓝色曲线)略小，原因在于实际系统运行时，存在外部扰动转矩，而动力学模型对此忽略不计。在一定误差范围内，二者波形基本一致，表明外骨骼参数辨识结果正确。

3 步态跟踪滑模控制 3.1 人机协同系统模型

在本项目中^[16,20]，外骨骼服的质量 $ {m_0} = 20{\text{ kg}} $ ，穿戴者身高168 cm，与外骨骼服尺寸完全相符合， $ {L_{\text{1}}} = 0.5{\text{ m}} $ ， $ {L_{\text{2}}} = 0.38{\text{ m}} $ ，体重 $ {m_{{\text{load}}}} = 70{\text{ kg}} $ ，根据成年人各部分肢体长度与身高的比例关系和各部分肢体质量与体重的关系来确定人体下肢的几何参数及惯性参数^[21]。大腿的质量为

$ {m_{{\text{1load}}}} \approx 0.14{m_{{\text{load}}}} = 9.8{\text{ kg}} $

(27)

小腿与足部的总质量为

$ {m_{{\text{2load}}}} \approx 0.055{m_{{\text{load}}}} = 3.85{\text{ kg}} $

(28)

由于人体肌肉组织分布比较均匀，为了计算简便，将大腿和小腿(含足部)的质心都取在几何中心位置，则有

$ \left\{ \begin{gathered} {l_{\text{1}}} = 0.25{\text{ m}} \\ {l_{\text{2}}} = 0.19{\text{ m}} \\ \end{gathered} \right. $

(29)

人体下肢转动惯量为

$ \left\{ \begin{gathered} {I_{{\text{1load}}}} = {{{m_{{\text{1load}}}}L_1^2} / {3 = 0.817{\text{kg}} \cdot {{\text{m}}^{\text{2}}}}} \\ {I_{{\text{2load}}}} = {{{m_{2{\text{load}}}}L_2^2} / {3 = 0.185{\text{kg}} \cdot {{\text{m}}^{\text{2}}}}} \\ \end{gathered} \right. $

(30)

将式(27~30) 代入式(4) 计算得到为

$ \left\{ \begin{array}{l} {X_{{\text{1load}}}} = 5.696{\text{ kg}} \cdot {{\text{m}}^{\text{2}}} \\ {X_{{\text{2load}}}} = 0.325{\text{ kg}} \cdot {{\text{m}}^{\text{2}}} \\ {X_{{\text{3load}}}} = 0.368{\text{ kg}} \cdot {{\text{m}}^{\text{2}}} \\ {X_{{\text{4load}}}} = 4.375{\text{ kg}} \cdot {\text{m}} \\ {X_{{\text{5load}}}} = 1.463{\text{ kg}} \cdot {\text{m}} \end{array} \right. $

(31)

因此，结合表1与式(31) ，得到外骨骼人机协同系统的一组动力学最小参数集为

$ \left\{ \begin{gathered} {X_{{\text{1HM}}}} = {X_1} + {X_{{\text{1load}}}} = 15.202{\text{ kg}} \cdot {{\text{m}}^{\text{2}}} \\ {X_{{\text{2HM}}}} = {X_2} + {X_{{\text{2load}}}} = 3.093{\text{ kg}} \cdot {{\text{m}}^{\text{2}}} \\ {X_{{\text{3HM}}}} = {X_3} + {X_{{\text{3load}}}} = 0.625{\text{ kg}} \cdot {{\text{m}}^{\text{2}}} \\ {X_{{\text{4HM}}}} = {X_4} + {X_{{\text{4load}}}} = 6.246{\text{ kg}} \cdot {\text{m}} \\ {X_{{\text{5HM}}}} = {X_5} + {X_{{\text{5load}}}} = 1.976{\text{ kg}} \cdot {\text{m}} \\ \end{gathered} \right. $

(32)

将式(32) 代入式(5) ，考虑摩擦力项，得到形如式(1) 的下肢外骨骼人机协同系统动力学模型，其中：

$ \left\{ \begin{gathered} {\boldsymbol{M}}(\theta ) = \left[ {\begin{array}{*{20}{c}} {{\text{15}}{\text{.202}} + 1.250\cos {\theta _2}}&{{\text{3}}{\text{.093}} + 0.625\cos {\theta _2}} \\ {{\text{3}}{\text{.093}} + 0.625\cos {\theta _2}}&{{\text{3}}{\text{.093}}} \end{array}} \right] \\ {\boldsymbol{C}}(\theta ,\dot \theta ) = \left[ {\begin{array}{*{20}{c}} { - 1.250\sin{\theta _2}{{\dot \theta }_2}}&{ - 0.625{\rm{sin}}{\theta _2}{{\dot \theta }_2}} \\ {{\text{0}}{\text{.625}}\sin{\theta _2}{{\dot \theta }_1}}&0 \end{array}} \right] \\ {\boldsymbol{G}}(\theta ) = \left[ {\begin{array}{*{20}{c}} {61.211\sin {\theta _1} + 19.365\sin ({\theta _1} + {\theta _2}) } \\ {19.365\sin ({\theta _1} + {\theta _2}) } \end{array}} \right] \\ {\boldsymbol{F}}(\dot \theta ) = \left[ {\begin{array}{*{20}{c}} { - 0.062{{\dot \theta }_1} - 2.415{{\rm{sgn}}} ({{\dot \theta }_1}) - 1.796} \\ { - 0.503{{\dot \theta }_2} - 1.521{{\rm{sgn}}} ({{\dot \theta }_2}) } \end{array}} \right] \\ \end{gathered} \right. $

(33)

3.2 基于模型上界的滑模控制

滑模控制(Sliding mode control, SMC)具有响应快、对参数变化和扰动不灵敏、鲁棒性强等特点，但其不连续的开关特性会引起系统高频抖振。

由机器人系统动力学特性可知^[17]：存在一个线性回归矩阵 ${\boldsymbol{\varPhi }}(\theta ,\dot \theta ,{\dot \theta _r},{\ddot \theta _r}) \in {{\boldsymbol{R}}^{2 \times 3}} $ 和依赖于机器人质量特性的惯性参数向量 ${\boldsymbol{p}} \in {{\boldsymbol{R}}^{3 \times 1}} $ ，使得机器人系统动力学模型式(1)满足式(34)的线性关系：

$ {\boldsymbol{T}} = {\boldsymbol{\varPhi }}(\theta ,\dot \theta ,{\dot \theta _r},{\ddot \theta _r}){\boldsymbol{p}} $

(34)

若 $\varPhi _{{{ij}}}^{\rm{r}} $ 和 $p_j $ 分别为线性回归矩阵和惯性参数向量中的元素， $\overline \varPhi _{{{ij}}}^{\rm{r}} $ 和 ${\overline {{p}}_{{j}}} $ 分别为其模的上界， $i=1,2 $ ， $j=1,2,3 $ ，则

$ {\boldsymbol{\varPhi }}(\theta ,\dot \theta ,{\dot \theta _{\rm{r}}},{\ddot \theta _{\rm{r}}}){\rm{ = [}}\varPhi _{{{ij}}}^{\rm{r}}{\rm{],}} \left| {\varPhi _{{{ij}}}^{\rm{r}}} \right| \le \overline \varPhi _{{{ij}}}^{\rm{r}},\left| {{{{p}}_{{j}}}} \right| \le {\overline {{p}}_{{j}}} $

(35)

下肢外骨骼机器人在单腿支撑且另一条腿摆动过程中，髋关节和膝关节实际角度 ${\boldsymbol{\theta }} = {[ {\begin{array}{*{20}{c}} {{\theta _1}}&{{\theta _2}} \end{array}} ]^{\text{T}}}$ ，期望角度 ${{\boldsymbol{\theta }}_{\text{d}}} = {[ {\begin{array}{*{20}{c}} {{\theta _{1{\text{d}}}}}&{{\theta _{2{\text{d}}}}} \end{array}} ]^{\text{T}}}$ ，误差 $ {\boldsymbol{e}} = {{\boldsymbol{\theta }}_{\text{d}}} - {\boldsymbol{\theta }} $ 。 $ {\boldsymbol{\varLambda }} $ 为一个二阶正对角矩阵，定义参考角速度 ${\dot{\boldsymbol{ \theta }}_{\rm{r}}} = {\dot{\boldsymbol{ \theta }}_{\text{d}}} + {\boldsymbol{\varLambda e}}$ ，则参考角加速度 ${\ddot{\boldsymbol{ \theta }}_{\rm{r}}} = {\ddot{\boldsymbol{ \theta }}_{\text{d}}} + {\boldsymbol{\varLambda \dot e}}$ ，取滑模面的切换函数为

$ {\boldsymbol{s}} = {\dot{\boldsymbol{ \theta }}_{\text{r}}} - \dot{\boldsymbol{ \theta }}={\dot{\boldsymbol{ \theta }}_{{\rm{d}}}}+{\boldsymbol{ \varLambda e}} -{\dot {\boldsymbol{ \theta}} = \dot {\boldsymbol{e}} + {\boldsymbol{\varLambda e}}} $

(36)

则

$ \dot{\boldsymbol{ s }}={\ddot{\boldsymbol{ \theta }}_{\text{r}}} - \ddot{\boldsymbol{ \theta }}={\ddot{\boldsymbol{ \theta }}_{\text{d}}}+{\boldsymbol{ \varLambda} }\dot {\boldsymbol{e}} - \ddot{\boldsymbol{ \theta }} = \ddot {\boldsymbol{e}} + {\boldsymbol{\varLambda }}\dot {\boldsymbol{e}} $

(37)

将滑模控制器的控制律设计为

$ {\boldsymbol{T}}{\text{ = }}{\boldsymbol{\bar k}}{{\rm{sgn}}} (s) + {\boldsymbol{s}}{\text{ = }}\left[ {\begin{array}{*{20}{c}} {{{\bar k}_1}{{\rm{sgn}}} ({s_1}) + {s_1}} \\ {{{\bar k}_2}{{\rm{sgn}}} ({s_2}) + {s_2}} \end{array}} \right] $

(38)

式中： ${\overline{k}}_{i}={{\displaystyle \sum }_{j=1}^{j}{\overline{\varPhi }}_{ij}^{\text{r}}{\overline{p}}_{j}}，i=1，2$ 。

为了减小滑模控制引起的抖振，在滑模控制器后串联一阶低通滤波器，传递函数为 $ G(s) =\dfrac{{\omega }_{c}}{s+{\omega }_{c}} $ 。

3.3 稳定性分析

令李雅普诺夫函数 $V(t) =\dfrac{1}{2}{{\boldsymbol{s}}}^{\text{T}}{\boldsymbol{M}}(\theta ) {\boldsymbol{s}}$ ，因为 ${\boldsymbol{M}}(\theta )$ 为对称正定矩阵，故 $ V(t) > 0 $ ，为正定。

$\begin{split} \dot{V}(t) =& \frac{1}{2}{{\boldsymbol{s}}}^{\text{T}}[{\boldsymbol{M}}(\theta ) \dot{{\boldsymbol{s}}}+\dot{{\boldsymbol{M}}}(\theta ) {\boldsymbol{s}}]+\frac{1}{2}{\dot{{\boldsymbol{s}}}}^{\text{T}}{\boldsymbol{M}}(\theta ) {\boldsymbol{s}}=\\& {{\boldsymbol{s}}}^{\text{T}}{\boldsymbol{M}}(\theta ) \dot{{\boldsymbol{s}}}+\frac{1}{2}{{\boldsymbol{s}}}^{\text{T}}\dot{{\boldsymbol{M}}}(\theta ) {\boldsymbol{s}} \end{split} $

(39)

由机器人系统动力学特性^[22]： $ {\boldsymbol{\dot M}}(\theta ) - 2{\boldsymbol{C}}(\theta ,\dot \theta ) $ 是一个反对称矩阵， ${{\boldsymbol{x}}^{\rm{T}}}({\boldsymbol{\dot M}}(\theta ) - 2{\boldsymbol{C}}(\theta ,\dot \theta ) ) {\boldsymbol{x}} = 0$ 。则有 ${{\boldsymbol{s}}^{\rm{T}}}({\boldsymbol{\dot M}}(\theta ) - 2{\boldsymbol{C}}(\theta ,\dot \theta ) ) {\boldsymbol{s}} = 0$ ，即： ${{\boldsymbol{s}}^{\rm{T}}}{\boldsymbol{\dot M}}(\theta ) {\boldsymbol{s}} = 2{{\boldsymbol{s}}^{\rm{T}}}{\boldsymbol{C}}(\theta ,\dot \theta ) {\boldsymbol{s}}$ ，将其代入式(39)，得

$ \begin{split} \dot{V}(t) =&{{\boldsymbol{s}}}^{\text{T}}{\boldsymbol{M}}(\theta ) \dot{{\boldsymbol{s}}}+{{\boldsymbol{s}}}^{\text{T}}C(\theta ,\dot{\theta }) {\boldsymbol{s}}=\\ &{{\boldsymbol{s}}}^{\text{T}}[{\boldsymbol{M}}(\theta ) ({\ddot{\theta }}_{\text{r}}-\ddot{\theta }) +{\boldsymbol{C}}(\theta ,\dot{\theta }) ({\dot{\theta }}_{\text{r}}-\dot{\theta }) ]=\\ &{{\boldsymbol{s}}}^{\text{T}}[{\boldsymbol{M}}(\theta ) {\ddot{\theta }}_{\text{r}}+{\boldsymbol{C}}(\theta ,\dot{\theta }) {\dot{\theta }}_{\text{r}}+{\boldsymbol{G}}(\theta ) -{\boldsymbol{T}}]=\\ &-{{\boldsymbol{s}}}^{\text{T}}[{\boldsymbol{T}}-{\boldsymbol{\varPhi}} (\theta ,\dot{\theta },{\dot{\theta }}_{\text{r}},{\ddot{\theta }}_{\text{r}}) {\boldsymbol{p}}] \end{split} $

(40)

将式(38) 代入式(40)，得

$ \begin{split} \dot{V}(t) = &- \left[{{\displaystyle \sum }_{i=1}^{2}{s}_{i}{\overline{k}}_{i}\mathrm{sgn}({s}_{i}) } + {{\displaystyle \sum }_{i=1}^{2}{s}_{i}^{2} - {{\displaystyle \sum }_{i=1}^{2}{{\displaystyle \sum }_{j=1}^{3}{s}_{i}}}}{\varPhi }_{ij}^{\text{r}}{p}_{j}\right] =\\ &- \left[ {{\displaystyle \sum }_{i=1}^{2}{{\displaystyle \sum }_{j=1}^{3} \left|{s}_{i}\right|{\overline{\varPhi }}_{ij}^{\text{r}}{\overline{p}}_{j} + }}{{\displaystyle \sum }_{i=1}^{2}{s}_{i}^{2} - {{\displaystyle \sum }_{i=1}^{2}{{\displaystyle \sum }_{j=1}^{3}{s}_{i}{\varPhi }_{ij}^{\text{r}}{p}_{j}}}} \right] \leqslant\\ & -{{\displaystyle \sum }_{i=1}^{2}{s}_{i}^{2}\leqslant 0}\\[-18pt] \end{split} $

(41)

所以， $ \dot{V}(t) $ 半负定。当 $ \dot{V}(t) \equiv 0 $ 时，必有 $ {\boldsymbol{s}} \equiv 0 $ ， $ {\boldsymbol{e}} \equiv 0 $ ， $\dot{\boldsymbol{ e}} \equiv 0$ 。由李雅普诺夫稳定性定理知，原点是系统的平衡态，系统在李雅普诺夫意义下渐近稳定。

3.4 仿真结果及分析

在MATLAB软件中用M语言编写关节期望输入轨迹、滑模控制算法、人机协同系统线性化模型及逆模型和曲线绘制等4个S函数^[22]。对于人机协同系统，仍然采用式(26) 描述的关节激励轨迹，在式(36) 中取 $\; {\boldsymbol{\varLambda }} = \left[ {\begin{array}{*{20}{c}} {12}&0 \\ 0&{12} \end{array}} \right] $ ，取滤波器的截止频率 $ {\omega _{\rm{c}}} = 15{{{\text{ rad}}} \mathord{\left/ {\vphantom {{{\text{rad}}} {\text{s}}}} \right. } {\text{s}}} $ 。选用离散变步长解法器ode45，设置仿真时间3 s，绝对误差为10⁻⁶。通过仿真得到髋关节和膝关节角度和角速度的跟踪波形如图7所示。

图7表明，采用滑模控制的外骨骼人机协同系统，除了角速度存在抖动外，步态轨迹跟踪效果良好。在零时刻，髋关节期望初始位置 $ {\theta _{1{\text{d}}}}(0) = - {30^\circ } $ ，膝关节期望初始位置 $ {\theta _{2{\text{d}}}}(0) = {0^\circ } $ ，而外骨骼双腿竖直向下，髋关节实际初始位置 $ {\theta _1}(0) = {0^\circ } $ ，膝关节实际初始位置 $ {\theta _2}(0) = {0^\circ } $ 。系统起动后，外骨骼膝关节立刻就能完全跟随期望轨迹变化；而髋关节在约0.15 s内先反向摆动微小角度再立刻正向摆回初始位置，直到 $ {\theta _{1{\text{d}}}}(0) = {0^\circ } $ 时，髋关节才会完全跟随期望轨迹变化。在实际过程中，因存在惯性，髋关节电机产生瞬时微弱震动。

如图8(a) 所示，滑模控制器输出转矩具有较大的抖振。在滑模控制器串联一阶低通滤波器后，转矩抖振大幅度减小,如图8(b) 所示。表明低通滤波器能够有效地抑制或减小滑模控制引起的抖振。

4 结论与展望

本文先对外骨骼开展参数辨识实验，再对人体下肢进行惯性参数计算，最后综合得到二连杆结构的人机协同系统精确的动力学模型。基于模型上界的滑模控制能实现步态轨迹精准跟踪，低通滤波器能有效减小滑模控制引起的高频抖振。这为下肢外骨骼动力学参数辨识提供了一种解决方案，为基于模型的控制方法提供了一种参考模型,也为下肢外骨骼人机协同系统的步态轨迹精准跟踪控制提供了一种参考方法。

将下肢外骨骼近似简化为二连杆结构，模型更简单，待辨识参数更少。在静态与动态结合的实验辨识过程中，利用MATLAB在线采集电机驱动器的电流和重新标定后的转矩系数来实时计算驱动转矩，减小误差，也无需在关节内安装力矩传感器，元件更少，结构更简单，结果更精确。与通常采用的求线性化回归矩阵的参数辨识方法相比，不仅待辨识参数更少，而且实现了惯性参数和摩擦力解耦，还避免了对回归矩阵求逆。但人体下肢惯性参数近似计算可能引起模型参数误差。可安装人机交互力传感器构成反馈环节，通过人工智能算法对模型参数进行自动调整，使物理模型与参数模型的高精度一致。引入一阶低通滤波器减小滑模控制引起的抖振，方法简单有效。还可通过改进滑模控制律或者人工智能算法来减小抖振。另外，外骨骼机器人助力行走的步态规划问题值得进一步研究。