空间机械臂作为空间站上最重要的设备之一，可以完成一系列在轨服务，例如修复空间站设备，辅助宇航员出舱作业，抓取微小卫星等任务^[1-3]。空间机械臂在太空中处于失重状态，为保证其在轨工作的安全性和可靠性，需预先在地面对空间机械臂进行大量的分析和实验^[4-6]。

地面零重力系统为机械臂提供模拟太空环境，其性能直接影响地面实验的效果。因此，需要对地面零重力系统卸载重力的效果进行评价。已有学者对各自设计的地面系统进行了评价，其中Sato等^[7]为吊丝悬挂系统提出了2种索力分配方案，包括吊索延长线通过质心和吊点在关节处。前者计算简单但误差大，后者干涉少，但如果臂杆质心不在轴线上，关节会产生附加力矩。因为机械臂结构复杂，各部分质量与质心难以精确得到，所以以索力作为评价指标不能准确反映系统模拟零重力的效果。高吾益^[8]考虑了单索零重力模拟装置中，吊丝摆角的测量误差，分析了由于悬索的柔性，传感器和编码器的惯性启动力对模拟对象位置的角度测量误差。李煜琦等^[9]以扭矩作为评价指标，但没有考虑摩擦力的影响，最终理论结果与实验结果有较大误差。

空间机械臂的地面零重力系统是一个机械臂与零重力系统耦合的复杂系统，机械臂关节扭矩能反映出零重力系统的卸载性能，同时零重力系统的摩擦力也将对机械臂关节扭矩产生影响。近几年来，很多学者对机械系统的摩擦力进行了研究。樊世超等^[10]建立了考虑关节弹性及摩擦的三连杆机构的动力学模型，通过数值计算发现，只有在运动开始阶段摩擦对机械系统的运动影响大。输入角速度越低，运动速度越低，更加容易出现粘滑现象。陶润等^[11]研究了一种滚动轴承摩擦力矩的测量方法，并运用最小二乘法对测试结果进行拟合，得到了计算轴承摩擦力矩的经验公式。Makkar等^[12]针对以往摩擦模型不连续的问题，建立了一个连续且随时间可导的摩擦模型，便于理论分析及推导。

本文将研究零重力系统各部分摩擦力对系统评价的影响规律，对机械臂的关节扭矩模型进行完善。建立机械臂初步的动力学模型，并绘制出未考虑摩擦力的关节扭矩图；考虑摩擦力条件下的动力学建模方法，并提出摩擦修正的策略；完成计算结果与实验结果的对比。

1 理论建模 1.1 系统介绍与受力平衡建模

零重力模拟系统如图 1所示，包括二维随动系统和恒拉力系统。恒拉力系统为吊索提供恒定大小的拉力，保证系统能时刻平衡重力。二维随动系统保证吊索时刻处于竖直状态，索力方向与重力始终平行。在关节及末端执行器上一共装有10个悬挂单元，为了简化系统，每两个临近的悬挂单元用横杆合并成一个吊点，通过吊索与零重力模拟系统相连。机械臂为串联结构，包括7个关节、2个臂杆、1个中央控制器和2个末端执行器。机械臂每个关节处均装有扭矩传感器，可以采集关节扭矩，保证实验的安全性和可靠性。实验时，机械臂一端固定在框架(舱体)上，另一端可以进行六自由度的运动。

在空间机械臂进行地面实验之前，需制作一个带有传感器的模拟机械臂，测量实验时的各项参数，以验证零重力系统的可靠性。模拟机械臂需要和空间机械臂具有同样的尺寸及质量，但为了降低成本，模拟臂制作选用较便宜的材料，并且臂杆设计成鼠笼结构。因为模拟臂安装有扭矩传感器，能够验证理论模型的准确性，因此本文理论建模及实验都是基于模拟臂展开的。

首先，对该系统进行初步的动力学建模。从固定端开始，将机械臂关节从1至7编号。在初始状态下，关节2和关节6的轴线方向平行，其余5个关节轴线互相平行，且垂直于关节2和关节6。图 2为关节的参数示意图，坐标系{i}和{i+1}为关节i和i+1关联的Cartesian坐标系。

从i到i+1部件的D-H旋转矩阵为^[13]

$ {}^i{\mathit{\boldsymbol{R}}_{i + 1}} = {\mathit{\boldsymbol{R}}_z}\left( {{\theta _{zi}}} \right){\mathit{\boldsymbol{R}}_y}\left( {{\theta _{yi}}} \right){\mathit{\boldsymbol{R}}_x}\left( {{\theta _{xi}}} \right). $

(1)

其中θ_zi、θ_yi和θ_xi分别为绕o_iz_i、o_iy_i和o_ix_i轴的翻滚、偏航和俯仰角。

假设关节初始速度和加速度均为0，关节i+1的角速度ω_i+1和速度v_i+1的递推式分别为：

$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{\omega }}_{i + 1}} = {}^i{\mathit{\boldsymbol{R}}_{i + 1}}{\mathit{\boldsymbol{\omega }}_i} + {{\mathit{\boldsymbol{\dot \theta }}}_i}{\mathit{\boldsymbol{e}}_{i + 1}},}\\ {{\mathit{\boldsymbol{v}}_{i + 1}} = {}^i{\mathit{\boldsymbol{R}}_{i + 1}}\left( {{\mathit{\boldsymbol{v}}_i} + {\mathit{\boldsymbol{\omega }}_i} \times {\mathit{\boldsymbol{p}}_{i,i + 1}}} \right).} \end{array} $

(2)

其中p_{i, i+1}为{i+1}坐标系原点O_i+1相对于O_i的位置矢量。

将ω_i+1和v_i+1关于时间求导，可以得到角加速度和线加速度的递推表达式如下：

$ \begin{array}{*{20}{l}} {{{\mathit{\boldsymbol{\dot \omega }}}_{i + 1}} = {}^i{\mathit{\boldsymbol{R}}_{i + 1}}{{\mathit{\boldsymbol{\dot \omega }}}_i} + {{\mathit{\boldsymbol{\dot R}}}_{i + 1}}{\mathit{\boldsymbol{\omega }}_i} \times {{\mathit{\boldsymbol{\dot \theta }}}_i}{\mathit{\boldsymbol{e}}_{i + 1}} + {\mathit{\boldsymbol{\theta }}_i}{\mathit{\boldsymbol{e}}_{i + 1}},}\\ {{{\mathit{\boldsymbol{\dot v}}}_{i + 1}} = {}^i{\mathit{\boldsymbol{R}}_{i + 1}}\left( {{{\mathit{\boldsymbol{\dot v}}}_i} + {{\mathit{\boldsymbol{\dot \omega }}}_i} \times {\mathit{\boldsymbol{p}}_{i,i + 1}} + {\mathit{\boldsymbol{\omega }}_i} \times \left( {{\mathit{\boldsymbol{\omega }}_i} \times {\mathit{\boldsymbol{p}}_{i,i + 1}}} \right)} \right).} \end{array} $

(3)

则关节质心处的线加速度${{{\mathrm{\tilde{\dot{\mathit{\pmb{v}}}}}}}_{i+1}}$可以表示如下：

$ {{\mathit{\boldsymbol{\tilde {\dot v}}}}_{i + 1}} = {{\mathit{\boldsymbol{\dot v}}}_{i + 1}} + {{\mathit{\boldsymbol{\dot \omega }}}_{i + 1}} \times {\mathit{\boldsymbol{r}}_{i + 1}} + {\mathit{\boldsymbol{\omega }}_{i + 1}} \times \left( {{\mathit{\boldsymbol{\omega }}_{i + 1}} \times {\mathit{\boldsymbol{r}}_{i + 1}}} \right). $

(4)

其中r_i+1为质心在关节坐标系{i+1}下的位置矢量。

惯性力F_i+1和力矩M_i+1的正向递推公式为：

$ \begin{array}{*{20}{l}} {{\mathit{\boldsymbol{F}}_{i + 1}} = {m_{i + 1}}{{\mathit{\boldsymbol{\tilde {\dot v}}}}_{i + 1}},}\\ {{\mathit{\boldsymbol{M}}_{i + 1}} = {{\mathit{\boldsymbol{\tilde I}}}_{{C_{i + 1}}}}{{\mathit{\boldsymbol{\dot \omega }}}_{i + 1}} + {\mathit{\boldsymbol{\omega }}_{i + 1}} \times {{\mathit{\boldsymbol{\tilde I}}}_{{C_{i + 1}}}}{\mathit{\boldsymbol{\omega }}_{i + 1}}.} \end{array} $

(5)

其中m_i+1是第i+1号关节的质量，${{{\mathrm{\tilde{\mathit{\boldsymbol{I}}}}}}_{{{C}_{i+1}}}}$是在质心坐标系下的惯性矩阵。

求出各关节上的惯性力和扭矩后，通过逆向递推公式可以计算关节i上的合力f_i和合力矩τ_i：

$ \begin{array}{*{20}{l}} {{\mathit{\boldsymbol{f}}_i} = {\mathit{\boldsymbol{F}}_i} + {}^{i + 1}{\mathit{\boldsymbol{R}}_i}{\mathit{\boldsymbol{f}}_{i + 1}},}\\ {{\mathit{\boldsymbol{\tau }}_i} = {\mathit{\boldsymbol{M}}_i} + {}^{i + 1}{\mathit{\boldsymbol{R}}_i}{\mathit{\boldsymbol{\tau }}_{i + 1}} + {\mathit{\boldsymbol{r}}_{i + 1}} \times {\mathit{\boldsymbol{F}}_i} + {\mathit{\boldsymbol{p}}_{ii + 1}} \times {}^{i + 1}{\mathit{\boldsymbol{R}}_i}{\mathit{\boldsymbol{f}}_{i + 1}}.} \end{array} $

(6)

因为关节坐标下，关于z轴的扭矩比其他两个轴扭矩大，并且安装在关节上的扭矩传感器用来采集z轴的扭矩，将合力矩点乘以z轴的单位向量e_i，得到绕o_iz_i轴方向的扭矩为τ_i，即

$ \begin{array}{l} {\mathit{\boldsymbol{\tau }}_{iz}} = \left( {{\mathit{\boldsymbol{M}}_i} + {}^{i + 1}{\mathit{\boldsymbol{R}}_i}{\mathit{\boldsymbol{\tau }}_{i + 1}} + {\mathit{\boldsymbol{r}}_{i + 1}} \times {\mathit{\boldsymbol{F}}_i} + } \right.\\ \;\;\;\;\;\;\;\;\left. {{\mathit{\boldsymbol{p}}_{ii + 1}} \times {}^{i + 1}{\mathit{\boldsymbol{R}}_i}{\mathit{\boldsymbol{f}}_{i + 1}}} \right) \cdot {\mathit{\boldsymbol{e}}_i}. \end{array} $

(7)

当机械臂在地面零重力系统中工作时，受到重力和吊索拉力的作用。将吊索拉力和重力变换到关节坐标系下的表达式为：

$ \begin{array}{*{20}{l}} {{}^i{\mathit{\boldsymbol{F}}_{Si}} = {}^0{\mathit{\boldsymbol{R}}_i}{\mathit{\boldsymbol{F}}_{Si}},}\\ {{}^i{\mathit{\boldsymbol{G}}_i} = m_i^0{\mathit{\boldsymbol{R}}_i}\mathit{\boldsymbol{g}}.} \end{array} $

(8)

将式(8)代入式(7)和(6)后，得到系统的逆动力学表达式，即关节合力f_i及扭矩τ_i的表达式为：

$ \begin{array}{*{20}{l}} {{\mathit{\boldsymbol{f}}_i} = {\mathit{\boldsymbol{F}}_i} + {}^{i + 1}{\mathit{\boldsymbol{R}}_i}{\mathit{\boldsymbol{f}}_{i + 1}} + {}^0{\mathit{\boldsymbol{R}}_i}\left( {{\mathit{\boldsymbol{F}}_{Si}} - {m_i}\mathit{\boldsymbol{g}}} \right),}\\ {{\mathit{\boldsymbol{\tau }}_i} = \left( {{\mathit{\boldsymbol{M}}_i} + {}^{i + 1}{\mathit{\boldsymbol{R}}_i}{\mathit{\boldsymbol{\tau }}_{i + 1}} + {\mathit{\boldsymbol{r}}_{i + 1}} \times } \right.}\\ {\left. {\left( {{\mathit{\boldsymbol{F}}_i} + {}^0{\mathit{\boldsymbol{R}}_i}\left( {{\mathit{\boldsymbol{F}}_{Si}} - {m_i}\mathit{\boldsymbol{g}}} \right)} \right) + {\mathit{\boldsymbol{p}}_{i,i + 1}} \times {}^{i + 1}{\mathit{\boldsymbol{R}}_i}{\mathit{\boldsymbol{f}}_{i + 1}}} \right) \cdot {\mathit{\boldsymbol{e}}_i}.} \end{array} $

(9)

1.2 未考虑摩擦力条件下的扭矩情况

该机械臂的主要目的是检测地面零重力系统的卸载性能，不需要进行太复杂的运动。考虑运动时机械臂干涉等情况后，给定关节运动时序规划如图 3所示。

当仅考虑质量质心偏差时，关节实验扭矩与理论计算扭矩随时间的变化曲线如图 4所示，图 4a—4f分别为关节1—6的扭矩图，1号关节最靠近固定端。由于关节7处于末端，扭矩变化很小，所以未绘出。可以看出，扭矩的理论值与实验采集值变化趋势基本相同。但是在关节运动启动时，扭矩会发生突变，与计算数据有较大的误差，这是因为没有考虑摩擦力的作用。由此可见，在地面零重力模拟实验中，摩擦力不能忽略。

进一步比对数据误差可以发现，6号关节的计算值与实验数据比较接近，但随着关节编号变小，末端关节的摩擦误差会累加至更靠近固定端的关节，使得计算值与实验数据误差越来越大。为了方便描述，后续文中将编号比关节i小的称作“前置关节”，编号大的称作“后置关节”。对于该机械臂上任意关节所受的摩擦力等于该关节本身运动产生的摩擦力与其后置关节产生的摩擦力之和。

2 考虑关节摩擦的模型 2.1 摩擦力模型分析及构建

摩擦力组成模型主要包括静摩擦、Coulomb摩擦、黏性摩擦、Stribeck摩擦等。如果关节运动为匀速，可以采用静摩擦与Coulomb摩擦模型，当接触面压力一定时，满足以下公式：

$ f = \left\{ {\begin{array}{*{20}{l}} {{f_e},}&{v = 0,\quad \left| {{f_e}} \right| < {f_{\rm{s}}};}\\ {{f_{\rm{s}}}sgn \left( {{f_e}} \right),}&{v = 0,\quad \left| {{f_e}} \right| \ge {f_{\rm{s}}};}\\ {{f_{\rm{c}}}sgn (v),}&{v \ne 0.} \end{array}} \right. $

(10)

其中：v为物体运动速度；f_e为外力；f_s为最大静摩擦力；f_c为动摩擦力，为常数。当关节静止且有运动趋势时，摩擦力随着电机驱动力的增加而增加；当关节运动稳定时，摩擦力保持恒定。因为摩擦力正负号与运动(运动趋势)方向有关，当关节运动换向时，摩擦力会产生突变。

2.2 摩擦模型改进分析

由于很多模型都是不连续的，有些模型是分段连续的，计算时，特别是求导运算中会使问题变复杂。Makkar等^[12]提出了具有连续性的摩擦模型，表达式为

$ \begin{array}{*{20}{c}} {f(\dot q) = {\gamma _1}\left( {\tanh \left( {{\gamma _2}\dot q} \right) - \tanh \left( {{\gamma _3}\dot q} \right)} \right) + }\\ {{\gamma _4}\tanh \left( {{\gamma _5}\dot q} \right) + {\gamma _6}\dot q.} \end{array} $

(11)

其中：静摩擦力可表示为γ₁+γ₄, tanh(γ₂${\dot q}$)－tanh(γ₃${\dot q}$)可表现原点附近的Stribeck摩擦力, γ₄tanh(γ₅${\dot q}$)为Coulomb摩擦力, γ₆${\dot q}$为黏性摩擦力。通过观察式(11)发现，该摩擦模型可写为3个因子相加而成f(${\dot q}$)=f₁(${\dot q}$)+f₂(${\dot q}$)+f₃(${\dot q}$)，分别是考虑静摩擦影响的Stribeck摩擦f₁(${\dot q}$)=γ₁(tanh(γ₂${\dot q}$)－tanh(γ₃${\dot q}$))、Coulomb摩擦f₂(${\dot q}$)=γ₄tanh(γ₅${\dot q}$)以及黏性摩擦f₃(${\dot q}$)=γ₆${\dot q}$。当γ₁=0.5, γ₂=100, γ₃=50, γ₄=0.1, γ₅=100, γ₆=0.05时，f(${\dot q}$)、f₁(${\dot q}$)、f₂(${\dot q}$)和f₃(${\dot q}$)的函数曲线分别在图 5中绘出。

在图中可以看出，摩擦力在原点附近有突变，这是因为Stribeck摩擦f₁(${\dot q}$)和Coulomb摩擦f₂(${\dot q}$)的共同作用；随着时间的推移，摩擦力略有减小，是因为Stribeck摩擦f₁(${\dot q}$)；一定时间后，摩擦力趋于稳定，几乎等于Coulomb摩擦f₂(${\dot q}$)的大小；但是摩擦力仍与坐标轴有一定夹角，因为黏性摩擦f₃(${\dot q}$)的存在。将摩擦力分别考虑，可以简化后续的摩擦力参数辨识的难度。

2.3 摩擦力参数辨识及修正

将图 4中关节1扭矩局部放大(80~180 s)，并将关节1和关节3运动角度绘制于同一图中，横坐标为时间，如图 6所示。观察扭矩波动和关节转角的关系。其中，扭矩第一个明显波动是在关节3开始运动时，第二个波动是关节1开始运动时。由此可见，末端关节运动时，也会对靠近固定端关节扭矩产生影响；扭矩先增大到某一极值，再缓慢下降，与节2.2中改进的摩擦模型较为接近。

由于每个关节滑轨处的摩擦力均不相同，很难通过实验直接测得摩擦力的大小，所以只能通过间接途径反推出摩擦力参数。然后再将计算得到的摩擦力对机械臂动力学模型进行修正，使得计算结果更加精确，更接近实验结果。

在节1.2中提到，靠近末端关节运动时产生的摩擦力，会作用在更靠近固定端的关节上。为了保证计算结果的准确性，需要从末端关节开始，对摩擦力进行逐一修正。关节扭矩修正流程如图 7所示。

2.4 悬挂单元模型

一个关节单元包括关节(含滑轨)、2个竖直连杆、三角板、子吊索和横杆，如图 8所示。接触点分别标为A-G，其中：A-E处装有轴承，为滚动摩擦，摩擦力较小；G、F处为滑动摩擦，摩擦力较大；H为关节质心；分别对应各接触点的压力；分别为关节、竖直连杆、三角板和横杆的重量。

横杆上的吊点A的位置根据两侧悬挂单元的质量决定，以保证自吊索上的拉力大小恰好等于重力。三角板参数根据关节质心进行设计，使得$\vartriangle DEC\cong \vartriangle GFH$。刘振等^[14]所提出过类似的相似悬挂原理。因为竖直连杆始终保持竖直，根据平行四边形原理，子吊索的延长线恒穿过关节质心H。

根据力学平衡，可以得到

$ \left\{ {\begin{array}{*{20}{l}} {{N_A} = {G_1} + 2{G_2} + {G_3} + \frac{1}{2}{G_4},}\\ {{N_B} = {G_1} + 2{G_2} + {G_3},}\\ {{N_D} + {N_C} = {G_1} + 2{G_2},}\\ {{N_C} + {N_F} = {G_1}.} \end{array}} \right. $

(12)

设ⁱh_ai和ⁱh_bi为G、F接触点在关节坐标系下的坐标，它们是关节i转角θ_i的函数：

$ \begin{array}{*{20}{l}} {{}^i{\mathit{\boldsymbol{h}}_{ai}} = Rot \left( {z,{\theta _i}} \right) \cdot {\mathit{\boldsymbol{h}}_{ai}},}\\ {{}^i{\mathit{\boldsymbol{h}}_{bi}} = Rot \left( {z,{\theta _i}} \right) \cdot {\mathit{\boldsymbol{h}}_{bi}}.} \end{array} $

(13)

其中${\mathop{\rm Rot}\nolimits} \left( {z, {\theta _i}} \right) = \left( {\begin{array}{*{20}{c}} {\cos {\theta _i}}&{ - \sin {\theta _i}}&0\\ {\sin {\theta _i}}&{\cos {\theta _i}}&0\\ 0&0&1 \end{array}} \right)$, 是绕z轴的旋转矩阵。h_ai和h_bi为关节转角为0时，G、F接触点的坐标，是关于关节结构的常数。

关节i运动对自身的摩擦扭矩为:

$ {}^i{\mathit{\boldsymbol{M}}_{fi}} = \left\{ \begin{array}{l} {}^i{\mathit{\boldsymbol{h}}_{ai}} \times {f_a}{\mathit{\boldsymbol{e}}_i} + {}^i{\mathit{\boldsymbol{h}}_{bi}} \times {f_{bi}}{\mathit{\boldsymbol{e}}_i},\;\;\;\;\;{{\dot \theta }_i} \ne 0;\\ 0,\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{{\dot \theta }_i} = 0. \end{array} \right. $

(14)

其中函数f_ai和f_bi为式(11)中提到的摩擦模型表达式。

后置关节k对关节i的摩擦附加扭矩类似，需先将后置关节滑轨连杆接触点转换到i关节坐标系下, ^kT_i为从关节k到关节i的变换矩阵：

$ \begin{array}{*{20}{l}} {\left( {\begin{array}{*{20}{c}} {{}^i{\mathit{\boldsymbol{h}}_{ai}}}\\ 1 \end{array}} \right) = {}^k{\mathit{\boldsymbol{T}}_i} \cdot \left( {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{h}}_{ai}}}\\ 1 \end{array}} \right),}\\ {\left( {\begin{array}{*{20}{c}} {{}^k{\mathit{\boldsymbol{h}}_{bi}}}\\ 1 \end{array}} \right) = {}^k{\mathit{\boldsymbol{T}}_i} \cdot \left( {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{h}}_{bi}}}\\ 1 \end{array}} \right).} \end{array} $

(15)

关节i上的摩擦扭矩包括由本身运动引起的以及所有后置关节运动引起的附加力矩, 即

$ {\left( {\sum {{\mathit{\boldsymbol{M}}_f}} } \right)_i} = \sum\limits_{m = i}^{10} {{}^i{\mathit{\boldsymbol{M}}_{fm}}} . $

(16)

3 实验结果与分析

模拟机械臂及零重力模拟装置实物图如图 9所示。实验时，机械臂一端固定在支撑框架上，其余部分可在框架上方做一定范围的运动。装置通过吊索平衡机械臂的重力，并对吊点位置进行跟踪，保证吊索时刻处于竖直方向。

考虑摩擦力之后，关节扭矩计算值与实验值曲线如图 10所示，平均扭矩误差如图 11所示。可以看出，扭矩的计算值与实验值相比未修正之前更为接近。但是关节2扭矩仍然有较大误差。因为在初始状态下，关节2与其他5个关节安装方向垂直，其他关节上的摩擦力矩方向与关节2轴线垂直，即关节2扭矩受其他关节摩擦力矩影响很小。从图 10b中也可以得到这个结论，其中绿色的实验数据线条在其他关节运动时，扭矩波动相比于其他关节小了很多。但是关节2的理论值与实际值仍有很大误差，这是由于模拟机械臂的臂杆为鼠笼结构，容易发生沿轴线方向的扭转变形。根据机械臂的构型可知，臂杆轴线方向与关节2轴向平行，因此臂杆变形会对关节2扭矩产生影响。关节2在机械臂运动结束时的扭矩值，与初始状态不同，也是因为此时臂杆存在扭转变形。

4 结论

本文面向地面零重力模拟系统卸载性能的评价，建立了系统七关节机械臂的动力学模型，计算出各关节扭矩并建立了基于摩擦力的连续关节扭矩模型。通过对比实验数据曲线与理论计算曲线可知，二者数值与变化趋势基本一致，除2号与6号关节外，其他关节扭矩理论值的平均扭矩误差减少了50%以上，反映出该拟合方法的有效性。但由于系统中的其他因素未被考虑，后续需要进一步研究2号关节扭矩误差产生的原因。