可行走着陆器的牛顿-欧拉法动力学研究

冯盖亚 贾山 陈金宝 周向华

（1 南京航空航天大学航天学院，南京 211106）

（2 深空星表探测机构技术工信部重点实验室，南京 211106）

（3 航天进入减速与着陆技术实验室，南京 211106）

0 引言

人类进入工业社会已久，随着近现代工业的发展，机器代替人工已成为各领域的发展趋势。当前已开展的深空探测任务中常采用多目标多任务的探测方式[1]，对于空间着陆器和深空探测机器人，需要在地外天体完成表面着陆，用技术手段获取特定环境样本，并带回地球。基于当前技术条件，着陆器着陆的姿态和位置无法完全确定[2]，而目前广泛应用的轮式探测器越障能力不足，不具备移动功能的着陆器在复杂地形很难成功释放轮式探测器并使其顺利开展移动探测。为了保证着陆器携带的各种设备正常工作，扩大星表探测范围，需设计一种缓冲行走一体化着陆器，确保相关设备工作状态良好，完成指定任务。目前所用的探测机器人都是可移动的，移动机器人的移动方式主要有4 种：腿式、轮式、履带式和轮腿式。轮式和履带式机器人的技术难度较低，但二者对复杂地形的适应能力不足，因此，通常选择腿式机器人用于空间探测和采样返回。常见的腿足式机器人有仿生双足机器人、仿哺乳动物的四足机器人、仿昆虫的六足以及八足机器人。其中双足机器人更能适应人类活动空间，适应性强；六足机器人稳定性和负载性能良好；四足机器人相较于六足机器人结构简单，相较于双足机器人负载能力更强，应用于可行走着陆器，可在外太空地表稳健行走。

深空探测中面临复杂的动力学环境[3]，为了使探测器平稳着陆并完成移动探测任务，需要对着陆环境以及着陆器本身进行动力学研究。相较于串联机构，并联机构结构刚度更高、稳定性好，但是工作空间较小，结构复杂[4]。并联机构的动力学分析对其在地外星体表面节能行走具有重要意义。

着陆器着陆缓冲过程及行走时的动力学分析是着陆器关节能耗模型建立的重要基础。当前已有多名学者在着陆器动力学方面进行了研究。对于着陆器行走时自身结构与机构的动力学，经典多体系统动力学的建模方法有牛顿-欧拉（Newton-Euler）法、拉格朗日（Lagrange）法、凯恩（Kane）方法等。这些动力学研究方法可以简单分为矢量力学方法和分析力学方法两大类[5]。

分析力学方法常见的有拉格朗日法、凯恩方法和虚功原理等，均是从能量角度出发，以整体的观点，对运动机构进行动力学研究，其中拉格朗日法和凯恩方法都是基于达朗贝尔原理发展而来。比如毛嘉炜[6]基于拉格朗日方法，建立了一种7 自由度的空间柔性机械臂的动力学模型，由此提出柔性矩阵建模理论的计算方法，并对机械臂有限元模型进行化简；N. Yu. Nosova 等[7]基于拉格朗日第二方程，得到了并联机构的动力学方程和工作空间内的固有频率；黄俊杰等[8]采用几何法和达朗贝尔原理法分析3-PRS并联机器人在执行器卡死状态下各关节内里的变化；吴璞等[9]以4-DOF 的Delta 高速并联机器人为研究对象，采用几何法求解机器人逆运动学方程，结合螺旋理论和凯恩（Kane）方法建立了机器人的动力学方程；柴馨雪等[10]采用螺旋理论分析并联机器人各关节和杆件的速度和加速度，结合虚功原理获得并联机器人运动时的驱动力。上述以能量形式求解得到的动力学方程，很难得到机构全部关节的力学信息。

矢量力学的典型方法是牛顿-欧拉（Newton-Euler）法，是基于牛顿方程和欧拉方程，以“连杆之间递推”的形式求得机构各关节力和力矩。牛顿-欧拉法的原理简单，获得的动力学公式意义明确[11]。田波等[12]基于牛顿-欧拉公式得到了完整绳驱动机械臂动力学模型，通过ADAMS 进行仿真验证；Bruno Scaglioni 等[13]基于牛顿-欧拉公式，采用空间矢量表示法建立了柔性机械臂的闭式动力学模型。实际应用中，也可以在牛顿-欧拉法的基础上，根据需求引入其他元素进行优化改进。例如侯雨雷等[14]根据牛顿-欧拉方法建立了具有球铰间隙的3RSR（R 为转动副，S 为球副）并联机构的动力学模型，并在其约束方程中引入标准四元数，采用四阶Runge-Kutta 方法求解3RSR 动力学模型。

本文动力学研究的对象为一种具有复杂串并联混合腿足机构的着陆器，目标是求得运动过程中腿足机构各关节受力情况，作为后面关节能耗研究和承力结构强度校核的支撑。所建立的动力学模型要获得腿足机构各部位在运动过程中的受力情况，这是基于能量的动力学分析方法很难做到的，因此选择牛顿-欧拉法进行动力学研究，用“递推”的形式分析着陆器腿足机构各部件运动特性，用牛顿方程和欧拉方程推导可以得到各部件所受惯性力和惯性力矩，列出各部件所受外力，求解得到各关节受力情况。此外，以往牛顿-欧拉法动力学研究的对象多为串联机构，递推形式简单，而本文研究的是可行走着陆器的串并联混合的腿足机构，过程相对复杂，需要进行更为详细的运动分析和受力分析。

1 着陆器构型设计

该可行走四足着陆器由一个正四棱台机体和4条对称分布的串并联混合机构悬臂腿组成，每条腿具有一个主支杆和两个辅助支杆，主支杆和辅助支杆均包括驱动机构和缓冲器[15]，与机体之间用胡克铰来连接，主辅支杆之间通过球铰连接。着陆器构型如图1 所示。改变主、辅支杆驱动机构的驱动行程，可以改变悬臂腿抬起或偏移的角度和悬臂腿长度，通过移动质心的位置至静止的3 腿外接圆内，从而创造第四条可移动支腿的收缩调整功能，进而实现着陆器的行走和姿态调整。

图1 着陆器结构示意Fig.1 Schematic diagram of the lander structure

2 着陆器腿足机构运动学分析

运动学模型在着陆器的运动规划和优化设计中有重要作用[16]，因此在进行动力学研究之前应先分析着陆器腿足机构的运动特性。根据D-H 法则，建立如图2 所示单腿的关节坐标系。其中OB是着陆器机体质心，以该点为原点建立着陆器机体质心坐标系；U1点为悬臂腿主支杆与机体的连接点，与胡克铰1的质心重合，以U1点为原点，建立基座坐标系{O0}，其中Z0轴沿胡克铰转轴方向，X0轴垂直于机体侧平面，向外为正。坐标系{O0}绕Z0轴旋转角度 θ1得到坐标系{O1}；坐标系{O1}在绕X1轴旋转-π/2 之后，再绕Z1轴旋转角度 θ2得到坐标系{O2}；将坐标系{O2}沿X2轴平移B1（主支杆固定部分长度）得到坐标系{O3}；将坐标系{O3}沿X3轴平移b1（主支杆伸缩部分长度）得到坐标系{O4}。主支杆总长度L1为B1和b1之和，左右辅助支杆总长度分别记为L2和L3，固定部分的长度为B2和B3，伸缩部分的长度为b2和b3（为了避免图线交叉重叠，用不同颜色的直线表示）。U2点、U3点是并联腿左右辅助支杆与机体的连接点，分别与胡克铰2 和胡克铰3 的质心重合，S2点和S3点分别为球铰中心。参照主支杆的D-H 建系过程，将坐标系{O3}平移，在球铰中心S2和S3分别建立坐标系{OS2}和{OS3}。

图2 单腿坐标系示意Fig.2 Schematic diagram of the single-legged coordinate system

2.1 单腿运动学分析

根据关节运动参数，用D-H 法推导连杆坐标系之间的变换矩阵，进行正逆运动学的求解。用Ti+1,i(i=0, 1, 2, 3)表示坐标系{Oi+1}到坐标系{Oi}的变换矩阵[17]，依次表示为：

坐标系{O4}到坐标系{O0}的变换矩阵T4,0可以表示为

根据式（1）~（2），可求出腿足机构末端S1在坐标系 {O0}中x、y、z方向上的坐标Ax、Ay、Az：

然后，根据式（2）~（3），可求得运动学逆解为

2.2 映射

在对主支杆进行运动学分析后，结合着陆腿结构和运动特性，对两辅助支杆进行运动学分析，求解其在主支杆运动过程中的映射运动状态，即辅助支杆伸缩量和关节角度。

如图2 中所示，设S2、S3在{O3}坐标系中的位置为(x2,y2,z2)、(x3,y3,z3)，在{O0}坐标系中的位置为(xS2,yS2,zS2)、(xS3,yS3,zS3)；U2、U3在{O0}中位置分别为(xU2,yU2,zU2)、(xU3,yU3,zU3)。

从坐标系{OS2}、{OS3}到坐标系{O3}的变换矩阵TS2,3和TS3,3分别为：

从坐标系{OS2}到坐标系{O0}的变换矩阵TS2,0可以表示为

从矩阵TS2,0的表达式可知，坐标系{OS2}的原点S2在坐标系{O0}中坐标为

从坐标系{OS3}到坐标系{O0}的变换矩阵TS3,0可以表示为

从矩阵TS3,0的表达式可知，坐标系{OS3}的原点S3在坐标系{O0}中坐标为

结合U2，U3在坐标系{O0}中的位置坐标(xU2,yU2,zU2)，(xU3,yU3,zU3)，两个辅助支杆的伸缩量b2、b3可分别表示为：

3 动力学建模

动力学研究将腿足机构的驱动力和运动联系起来[18]，常用的动力学建模方法有牛顿-欧拉法、拉格朗日法和虚功原理等。拉格朗日法和虚功原理是从能量角度对机构整体进行动力学分析，而牛顿-欧拉法则是对机构各个关节受力进行递推分析。为研究着陆器在运动过程中各个关节的受力情况与运动轨迹的关系，基于图2 中所示关节坐标系，采用牛顿-欧拉法建立着陆器整机的全状态动力学模型。

为了计算腿足机构各连杆受力，需要先求出各连杆的角速度、角加速度和连杆质心的线加速度等结果。对于主支杆中的两个杆件，与机体铰接的连杆记为连杆1，连杆1 质量为mL1，与坐标系{O2}固连；与足垫相连伸缩部分记为连杆2，连杆2 质量为mL2，与坐标系{O3}固连。在坐标系{O2}中连杆1 的质心坐标为rc1，rc1=[xc1yc1zc1]T，其对于O2点的惯性张量记为IL1；在坐标系{O3}中，连杆2 的质心坐标为rc2，rc2=[xc2yc2zc2]T，其对于自身质心的惯性张量记为IL2。

3.1 摆动相

当着陆器处于摆动相时机体静止，3 条腿处于支撑状态，1 条腿处于摆动状态。首先要计算连杆或者坐标系的运动学参数，包括角速度、角加速度、连杆质心或者坐标系原点的线加速度等。从摆动腿固定端开始推导，将重力因素考虑在内，基座坐标系{O0}原点处的角速度 ωO0、角加速度O0、原点的线加速度O0分别为

式中 θ0为坐标系{O0}中Y0轴的正方向与重力方向夹角，数值取21°；g为着陆器所在星体表面的重力加速度。

为了求得连杆1 的运动参数，需先求得与连杆1 固连的坐标系{O2}的角速度、角加速度和坐标系原点的线加速度。坐标系之间速度和加速度的变换需要通过坐标系变换矩阵，为了方便运算，将坐标系变换矩阵转化为如下的三阶形式：

由运动学分析中各坐标系的建立过程可知，坐标系{O0}绕Z0轴旋转角度 θ1得到坐标系{O1}，坐标系{O1}的角速度 ωO1、角加速度O1和坐标系原点的线加速度O1分别为：

坐标系{O1}经过旋转变换后得到坐标系{O2}，两次旋转的角度分别为 θ2和-π/2，经计算可得坐标系坐标系{O2}的运动参数。坐标系{O2}与连杆1 固连，二者的角速度和角加速度相同，坐标系原点O2的线加速度即是连杆1 固定端的线加速度。因此连杆1 的角速度 ωL1、角加速度L1和连杆质心的线加速度c1分别为可以表示为：

再计算连杆1 所受的惯性力和惯性力矩。根据牛顿方程可知，连杆1 在质心处的惯性力FL1为

根据欧拉方程，连杆1 对于固定端点处的惯性力矩NL1为

然后分析连杆2，连杆2 相对连杆1 做平移运动，相对速度为。为了求得连杆2 所受惯性力和惯性力矩，首先需计算连杆2 的角速度 ωL2、角加速度L2和连杆质心的线加速度c2，见式（24）~（26）。由于连杆2 与连杆1 通过移动关节进行相对运动，因此还需要补充连杆2 与移动关节重叠位置的加速度L2（见式（27））。

根据牛顿方程和欧拉方程，可以得到连杆2 在质心处的惯性力FL2与惯性力矩NL2分别为：

着陆器在行走过程中单腿受力情况如图3 所示。主支杆包括连杆1 和连杆2，相应的两条辅助支杆分别由连杆3 和连杆4、连杆5 和连杆6组成，连杆1、连杆3、连杆5 的端点受基座作用力分别记为FJ1、FJ2、FJ3。将连杆i对连杆j的作用力和力矩分别记为Fij和Mij（i,j= 1, 2, 3,4, 5, 6，表示连杆编号），力Fij对一点的力矩记为MFij；足端力（足端除重力以外的作用力）记为Fe。本章动力学分析目的是求出各关节作用力，不考虑摩擦因素，这些作用力包括FJ1、FJ2、FJ3以 及F12、F34、F56。

图3 单腿受力示意Fig.3 Schematic diagram of the force on a single leg

对于摆动腿，足端没有外部作用力，Fe的数值为0，即

对连杆2 进行受力分析（见图3），可得出其质心处的受力和力矩具有如下关系：

式中MFe为足端力Fe 对连杆2 质心的力矩。

再分析连杆1，可以得到力平衡方程为

为了求解式（33），对左右两个辅助支杆进行分析。参照前文主支杆的运动分析过程，推导辅助支杆的运动学参数。图3 中左侧辅助支杆（由连杆3 和连杆4 组成）的关节转角分别记为 θ3和 θ4，右侧辅助支杆（由连杆5 和连杆6 组成）的关节转角分别记为 θ5和 θ6。参照式（4）所示的主支杆逆运动学结论公式，通过S2和S3在各自基座坐标系中的位置，可分别求得两个辅助支杆的关节转角。

参照主支杆的关节坐标系，可以用同样的方法建立左侧辅助支杆（由连杆3 和连杆4 组成）的关节坐标系。将重力因素考虑在内，进行连杆3 和连杆4 的速度和加速度的正向推导。通过牛顿方程和欧拉方程可以得到连杆3 所受惯性力FL3和惯性力矩NL3，以及连杆4 所受惯性力FL4和惯性力矩NL4。

对连杆3 进行受力分析，在质心处可以得出力平衡方程为

在实际运动过程中，力FJ2的作用点并不在连杆3 的中轴线的延长线上，设该作用点沿连杆转轴到质心的距离为c，这里将其处理为作用于轴线的力FJ2和力矩MJ2，其中将力FJ2的分量分别记为FJ2x、FJ2y、FJ2z，则因此，在连杆3端点处的力矩平衡方程为

再分析连杆4，其质心处的力和力矩的平衡方程分别为：

将式（35）与式（37）相加，根据牛顿第三定律，相互作用力矩M34和M43可抵消，则简化后可得

牛顿-欧拉法动力学模型需要求得转动关节力FJ1、FJ2、FJ3以及移动关节力F12、F34、F56。目前已有的方程不足以求解未知量，为了求解各关节力，单独分析与辅助支杆2 相连接的胡克铰2。将胡克铰质量记为mh，胡克铰在自身固连坐标系中对于质心的惯性张量表示为矩阵形式，记为Ih，Ih=diag(Ixh,Iyh,Izh)，Ixh、Iyh、Izh为惯性张量在坐标系3 个坐标轴方向的分量。

在胡克铰2 本身的固连坐标系中，胡克铰质心处的惯性力和惯性力矩分别为：

在胡克铰本身的固连坐标系中，将连杆3 对胡克铰2 的作用力矩记为Mlh2，坐标轴x、y、z方向上的分量分别为mlh2x、mlh2y、mlh2z，所以Mlh2=[mlh2x mlh2y mlh2z]T。结合胡克铰受力特点，可知mlh2z=Izh3。胡克铰2 对连杆3 的作用力矩为 -Mlh2，将该力矩变换到连杆3 的固连坐标系中，在忽略胡克铰的情况下，MJ2可以表示为

联立式（34）、（38）及式（41），求解关节力 、 、 。

对于辅助支杆3（由连杆5 和连杆6 组成）用同样的力学分析方法，求得关节力FJ3、F16及F56，根据牛顿第三定律，求得F14和F16的反作用力F41和F61，代入式（33），求解FJ1。

所求得的FJ1、FJ2、FJ3以及F12、F34、F56的表达式即为动力学模型，由于具体计算过程过于繁复，这里不做展示。

3.2 支撑相

机体前移过程中，四足支撑点不变，质心向前移动，这个运动过程中着陆器腿足机构可简化为图4所示形式。其中，Ai（i=1, 2, 3, 4）为着陆器四足支撑点，Hi为着陆腿与机体的铰接点，G表示重心。

图4 机体质心位移示意Fig.4 Schematic diagram of the displacement of the body’s center of mass

以初始时刻质心在地面的投影为原点建立地面坐标系{Og}，坐标系Z轴正方向与前进方向一致。

同样，采用单腿的动力学分析方法分别在4 个铰接点Hi处建立4 条腿的基座坐标系。此时作为腿足机构末端的Ai相对于Hi运动，以Ai在基座坐标系中的坐标作为输入，参照式（4）分别求解4 条腿的关节运动参数。参考单腿的动力学分析过程，分别用牛顿-欧拉法进行4 条腿的动力学分析，以足端受力为输入，求解4 条腿在机体质心位移过程中的关节力和力矩。

用牛顿-欧拉法得到的着陆器动力学模型，可以获得着陆器一个行走周期内所有关节处的受力和力矩，作为后续能耗模型建立和优化的基础。根据着陆器行走过程中关节受力和关节能耗，可以进行着陆器构件的强度校核，并选择最优的关节驱动电机，保证着陆器稳定地行走。

4 动力学模型仿真和验证

着陆器在支撑相和摆动相的动力学求解过程均是以腿足机构末端力为外部输入，递推各运动关节受力。为了验证动力学模型，需要为着陆器腿足机构设计一段运动轨迹进行仿真。在完成仿真之后，测量运动过程中的关节力，将理论计算结果与虚拟样机仿真结果进行比较，验证模型在能耗优化中的有效性。

用5 次样条插值法规划着陆器在水平地面上单腿向前迈步（见图5）的足端轨迹，设计步长为100 mm，步高为50 mm，迈步时间为100 s，腿足机构在运动轨迹的起点和终点时速度为0。为了保证腿足机构具有一定越障能力，同时设定在时间为20 s 时，足端高度为20 mm；在时间为40 s 时，足端高度为50 mm；在时间为80 s 时，足端高度为20 mm。用轨迹上的5 个定点坐标，进行5 次多项式插值，得到轨迹曲线。足端轨迹曲线在x-y平面上的具体图形如图6 所示。着陆器单腿向前迈步，足端Ai沿所规划的足端轨迹前移。

图5 单腿迈步示意Fig.5 Foot end trajectory design

图6 足端轨迹规划Fig.6 Schematic diagram of a one-legged step

采集400 个足端轨迹上各点的二维坐标数据，通过几何关系得到运动过程中机构末端S1坐标（Ax，Ay，Az）的一组取值。根据式（4），可以得到腿足机构上关节的位移量 θ1、 θ2和b1的取值。以 θ1、 θ2和b1作为输入，根据章节2 的运动学分析获得两个辅助支杆在运动过程中各个关节的一系列关节位置数值。用微分法对每组离散数值求导，得到腿足机构中各关节的速度和加速度。此时这些关节运动数据均为400×1 的向量形式，代入3.1 节建立的牛顿-欧拉动力学理论模型，用理论模型计算得出腿足机构各关节所受力的一系列取值，由此即可得出对应以上足端轨迹的着陆器各关节受力情况的变化曲线。

同时运用ADAMS 软件建立着陆器的虚拟样机模型，并添加关节约束。设置水平地面和重力加速度g，使着陆器的3 条支撑腿落在地面上的位置保持不变。为着陆器的摆动腿规划图6 所示足端运动轨迹，进行四足着陆器虚拟样机的运动仿真，采集着陆器摆动腿在相同轨迹下各关节的受力数据曲线，与通过理论模型计算出的结果进行对比验证。图7~9 展示了连杆1、3、5 在移动关节处的作用力在连杆固连坐标系中x、y、z方向的变化曲线。

图7 连杆1 驱动关节力理论值与仿真值对比Fig.7 Comparison between the theoretical value and the simulated value of the driving joint force of the connecting rod 1

图8 连杆3 驱动关节力理论值与仿真值对比Fig.8 Comparison between the theoretical value and the simulated value of the driving joint force of the connecting rod 3

图9 连杆5 驱动关节力理论值与仿真值对比Fig.9 Comparison between the theoretical value and the simulated value of the driving joint force of the connecting rod 5

在运动过程中着陆器腿足机构由移动关节驱动，腿足机构中的移动关节是驱动关节。按照腿足机构各支杆的坐标系设计，关节驱动力会沿x轴，而不考虑摩擦的关节能耗与关节驱动力相关。由图7~9 所示结果可以看出，与虚拟样机仿真结果相比，在y方向和z方向上，由于腿足机构几何形状、摩擦力和碰撞的影响，动力学模型计算的结果存在数值上的误差；在x方向，即驱动关节驱动的方向，理论计算结果具有与仿真结果近似的趋势和数值，在后续能耗模型建立和优化工作中可以起到有效作用。

5 结束语

本文设计了一种悬臂式四足可行走着陆器构型，在着陆器运动学模型基础上，采用牛顿-欧拉递推方法进行着陆器行走过程的动力学分析与建模。通过虚拟样机仿真结果与理论计算结果的对比，验证了动力学模型的有效性。所建立的动力学模型，可获得着陆器每个关节的受力情况。在此基础上，可以继续深入研究，结合关节摩擦力模型，建立指定关节的能耗模型和整机行走的能耗模型。而在着陆器的行走能耗中，运动步态是影响能耗的重要因素之一[19]，后续可以选择优化算法，进行着陆器行走轨迹的优化，得到着陆器行走的能耗最优轨迹[20]。