乘马康复混联机器人机构设计与运动分析

Mechanical Design and Motion Analysis of a Hybrid Robot for Equine Rehabilitation

摘要

马术治疗是针对脑瘫和自闭症等疾病的最热门有效的康复治疗手段之一，但因受到生物马匹康复治疗多种因素的限制，无法大规模推广应用。针对马术康复治疗需求，设计一种3‒RPS‒P型混联乘马康复机器人，其主要由并联机构和移动副机构两部分组成，具有刚度高、承载能力强、灵活度较高的特点。首先，建立机器人的运动学模型，推导出机器人运动学逆解解析解，采用牛顿迭代法计算得到机器人运动学正解的数值解，分析机器人在运动副限制下的工作空间并绘制出空间图像。然后，利用算例对上述运动学理论进行对比分析，验证了本文方法的正确性和有效性。接着，在MATLAB环境中将上述运动理论编写成代码，对机器人末端进行路径和轨迹规划设计，联合SolidWorks Motion模块对机器人进行运动仿真，仿真运动过程连续、轨迹平滑且与规划设计中的路径和速度目标一致，表明机器人可以实现设定的位姿规划运动，并进一步验证了理论分析的结果。最后，搭建乘马康复机器人实验平台，针对运动仿真中的位姿规划开展运动实验，实验结果验证了本文机构的合理性。该混联乘马康复机器人具备马术治疗所需的运动能力。

Abstract

Objective Equine therapy, also referred to as therapeutic riding or equine-assisted therapy, is a clinical methodology that promotes physical, psychological, and social functioning through the integration of horse handling and riding skills. It is generally overseen by a qualified therapist or equestrian instructor, and serves as a rehabilitation tool amalgamating equine and therapeutic techniques to alleviate conditions such as cerebral palsy and autism. Whilst equine therapy has demonstrated efficacy in rehabilitating particular conditions, its dissemination is impeded by several constraints and challenges. Furthermore, equine therapy presents a number of limitations and challenges. High expenses, specialised facilities, and the specificity of horses impede the implementation of the therapy on a broader scale. Furthermore, challenges related to individual variations and the requirement for specialised therapist training need to be addressed, alongside cost, resource, and training concerns, to ensure wider adoption of equine therapy. To promote its application, an alternative to conventional biological equine exercise rehabilitation therapy could be a horse-riding rehabilitation robot. Such a development could improve the accessibility and effectiveness of equine therapy on a broader scale. The horse saddle during various gaits involves four primary movements: up and down, forwards and backwards, left and right tilting, and forward and backward tilting. These movements are integral to understanding the saddle’s function. Clear conceptualization of these movements is essential for riding and safely using the horse saddle. Methods Therefore, the robot horse mechanism must possess the four degrees of freedom mentioned earlier. The high load and stiffness requirements when accommodating a person have been taken into account in the design process. As a result, a 3-RPS parallel mechanism is selected as the primary component, capable of carrying out up-and-down movement, left-right tilt, and forward-backward tilt of the saddle. The sketch of the mechanism is demonstrated in Fig. The saddle’s forward-and-backward movement is achieved through one moving vice, while another facilitates the up-and-down and back-and-forth movements. To reduce load and inertia, the mobile vice is mounted on the parallel mechanism’s moving platform. Together, the parallel mechanism and mobile vice comprise a 3-RPS-P hybrid robot. The mechanism possesses excellent stiffness, high load-carrying capacity, and remarkable flexibility, fulfilling the requirements for the saddle centre’s movement during rehabilitation therapy and load carrying capacity during horse riding. Firstly, the robot’s kinematic model was established, followed by deriving the analytical solution for the inverse kinematics. The positive kinematics solution was then obtained by utilising Newton’s iterative computation method. The relationship between the robot’s operational space and the joint space velocity Jacobian was investigated. Additionally, the robot’s workspace was analysed while considering the restriction of motion, and a spatial image was plotted accordingly. Finally, using Lagrange’s method, the dynamics of the parallel mechanism were established. Next, an examination is conducted on the correlation between the motion of the mechanism and the force applied. The motion theories outlined above have been compared and analyzed through examples, and the results have confirmed their validity and accuracy. These theories have then been transferred into MATLAB code, and the robot end’s path and trajectory planning and design have been carried out. Results and Discussions A joint motion simulation of the robot utilizing the SolidWorks Motion module has also been conducted, demonstrating a consistent and smooth trajectory in line with the planning and design objectives. This simulation showcases the robot’s ability to successfully attain its set goals and objectives. The simulation is seamlessly executed and aligns with the planned and designed objectives, demonstrating the robot’s capability to accomplish the designated position-planning motion. Furthermore, it verifies the theoretical analysis outcomes. Finally, the experimental platform for the horse-riding rehabilitation robot was established in order to conduct motion experiments for position planning in motion simulation. The experimental results confirmed the rationality of the mechanism, attained the intended motion trajectory, and showcased rapid response of the control system. Due to the complexity of its calculations, the kinematics of the positive solution is time-consuming, taking 2~3 ms. However, as the display and monitoring meet usage requirements, this signifies that the hybrid horse rehabilitation robot possesses the necessary movement ability for equestrian therapy. Conclusions The locomotive behavior of horses in various gaits, including fast walking, running, and jumping, will be further examined through the study of the robot. To analyse the impact of varying movement gaits, amplitudes and speeds on rehabilitation therapy and optimise mechanism and motion control with input from equestrian therapists’ riding feedback and suggestions. Additionally, establish a real-time feedback system within the machine’s simulation of horse movement to monitor and adjust the robot's motion status. To assess the effectiveness of robotic horse rehabilitation therapy, the performance and recovery of the experimental group will be compared to that of the control group. Furthermore, safety and stability will be prioritized in the mechanism and motion control design to prevent falls and potential injuries to users during rehabilitation therapy.

Graphical abstract

关键词

混联乘马康复机器人 / 机构设计 / 运动学 / 工作空间 / 轨迹规划 / 运动仿真

Key words

hybrid ride-on-horse rehabilitation robot / mechanism design / kinematics / workspace / trajectory planning / motion simulation

本刊网刊

脑瘫是一种发生在婴幼儿早期的神经系统综合症，该症状会严重影响人体肌肉协调、运动和平衡能力^[1]。在中国脑瘫患者中，0～6岁儿童发病率为0.23%，2～17岁脑瘫患病率为0.29%^[2‒4]。患者中以儿童居多，若不及时介入治疗，可能会引发一系列继发性障碍，导致生活不能自理^[5‒6]。马术治疗是利用马的自然步态和运动模式为患者提供运动和感官输入，从而建立改善神经和感觉功能的媒介^[7‒8]，对人体的感觉系统、前庭系统、眼部系统和肌肉骨骼系统等具有良好的治疗效果，已被广泛应用于患有孤独症、脑瘫、行为障碍、精神障碍等神经系统和交流障碍疾病的患者^[9‒11]。Park等^[12]对34名年龄在3～12岁的脑瘫儿童进行持续8周的马术治疗训练，通过对比大量的治疗案例发现，经过马术治疗后患者的儿科功能独立性评定-家庭版（PEDI-FSS）评分明显改善。Kwon等^[13]对96例4～10岁脑瘫患者进行分组，一组接受为期8周的乘马康复治疗，另一组不做任何干预疗法，研究结果表明经过马术治疗后的脑瘫儿童患者的粗大运动功能得到显著改善。临床实践和科研成果证明，脑瘫儿童在经过长期的马术治疗后，在步态平衡、言语表达和情绪控制等多方面均得到了很大的改善。

虽然马术治疗是最热门有效的康复治疗技术之一，但因生物马匹、特定场景和特定康复人员的种种要求和限制而无法大规模推广应用。因此，能够模仿真实马匹的运动器械应运而生，为需要进行马术治疗的人群提供了便利。Eskola等^[14]研究真实马匹的运动特性，设计一款6自由度的机械仿生马，可以近似模拟马匹运动的步态，但结构和控制系统复杂导致成本过高。Lee等^[15]研发一种骑马运动的模拟器，可以近似模拟不同马匹的多种不同运动步态，灵活性和可操作性较高，模拟器虽支持多种步态，但依赖预设程序，无法实时适应骑手动作，且步态过渡生硬，缺乏真实马匹的柔顺性。Yu等^[16]利用可控制鞍座运动强度的简易骑马模拟器Srider，验证了室内骑马运动对运动能力和腰肌功能的影响，但无法模拟马鞍的6维运动，限制了康复训练的针对性。张博^[17]、冯光磊^[18]基于一种拟合的方法对骑马机的结构设计进行深度的研究，通过凸轮机构和曲柄摇杆机构的组合实现模拟生物马的部分运动规律，但无法支持弯腰、直立等多姿态骑行。余顺年等^[19]通过多连杆组合机构，利用脚踏板作为驱动源，实现一种结构简单的机械式骑马机，无法模拟被动骑行模式，灵活性低。对于马术治疗的具体疗效的分析和量化仍是制约其发展的一大障碍之一^[20]。目前，基于图像和视频信息的运动评估方法主要分为两类：一类为依赖光学标记物的方法^[21]；另一类为纯粹依靠图像，不依赖于光学标记物的评估方法^[22]。例如：Sánchez等^[23]采用光学动作捕捉系统测试脑瘫儿童在进行躯干活动过程中头部姿态的变化，但未整合颈部肌电信号或眼动数据，无法解释神经肌肉控制机制；Coleman Wood等^[24]采用基于光学标记的动作捕捉系统评估6例脑瘫儿童通过机器人体位反馈干预疗法进行康复治疗后上肢运动能力的改善，但缺乏肌张力的量化分析。

生物马在运动时具有慢步、快步、跑步、跳跃等多种步态及不同步态下的坡度运动和转弯运动等，且在不同情境下各种运动步态的幅度和速度也不尽相同。本文以实现生物马运动时的可调多种步态为任务需求，设计一种3‒RPS‒P型混联机器人，该机器人能够实现马匹不同运动步态时马鞍的位置运动和姿态变换；通过伺服电缸控制机器人运动，可调节不同步态下的运动幅度和速度等参数；计算该混联机器人的自由度并建立运动学模型；分析机器人机构的工作空间，并对该机构进行轨迹规划和运动仿真分析；最终，完成实验样机的搭建并开展运动实验。

1 3‒RPS‒P型混联机器人机构设计

基于虚拟现实原理、6轴平行机制和运动捕捉技术，复现马匹不同步态时马鞍的3维运动模型^[25]。在马匹的不同步态运动下马鞍主要有上下运动、前后移动、左右倾斜和前后倾斜4种运动方式，所设计的机器马机构不仅需要具有上述4个自由度，还需要考虑机器马载人运动时对机构的负载和刚度的要求较高，因此，本文设计一种3‒RPS‒P型混联乘马康复机器人，其主要由并联机构和移动副机构两部分组成，即设计一种3‒RPS型并联机构作为机器人的主体部分，该并联机构可实现马鞍的上下运动、左右倾斜和前后倾斜3种运动；马鞍的前后运动由另一移动副实现，为减小移动副运动时的负载和惯量，将该移动副设计安装在并联机构的动平台上。

3‒RPS型机构简图如图1所示，该机构由静平台、动平台以及连接动静平台之间的3条运动支链组成。3条支链的拓扑结构完全相同，一端由转动副（R）与静平台相连，另一端由球副（S）与动平台相连，中间移动副（P）作为驱动副，对称分布在动静平台之间。

1.1 自由度计算

3‒RPS型机构为空间闭环机构，采用修正的G‒K公式计算机构自由度：

M=d(n-g-\mathrm{1})+{\displaystyle \sum _{i=\mathrm{1}}^g}{f}_i+v-\xi（1）

式中，M为机构的自由度，d为运动阶数，n为包括机架的构件数，g为运动副数，f_i 为第i个运动副的自由度数，v为并联冗余约束数，ξ为机构局部自由度数。

3‒RPS型机构的运动阶数d=6，总构件数n=8，运动副数g=9，运动副总自由度数为15，不存在并联冗余约束，局部自由度为0，将以上参数代入到式（1）中计算得到机构的自由度为3，即机构具有3个方向的自由度：1个平动自由度和2个转动自由度。另外，在3‒RPS型机构的动平台上增加1个与并联机构不重叠的平动自由度来实现马鞍的4个自由度不同方向的运动，与并联机构共同组成3‒RPS‒P型4自由度混联机构。

1.2 结构模型设计

根据图1中的机构原理，为了使机器马幅度和速度可调，且能够快速响应，选择4个直线伺服电缸作为驱动执行机构，3‒RPS‒P型混联机构的3维模型如图2所示。该混联机构主要由静平台、转动副、伺服电缸、球副、动平台、移动副组件，以及马鞍安装座组成。其中：伺服电缸直线输出部分通过球副与动平台相连，安装座通过转动副与静平台相连；马鞍安装座固联在移动副组件上，移动副组件由直线伺服电缸驱动。转动副中心及球副中心均呈等腰三角形形状布置于动静平台之间，三角形最长之高按照马鞍常用长度进行设计。在此基础上，增加机器马头部和颈部运动部件，设计机器马外观结构部件，最终机器马总体结构模型如图3所示。

2 运动性能分析

机器马运动性能分析是后续运动仿真和样机试验的基础，需要先建立机器马的运动学模型，再分析机构的位姿正逆解及速度雅可比，并在此基础上分析机构的位姿工作空间。

2.1 运动学逆解分析

3‒RPS‒P型混联机构由于移动副组件所提供的运动自由度与并联部分独立，因此，可独立分析并联机构部分与移动副组件。建立如图4所示的并联机构坐标系。

由图4可知，该并联机构的动静平台与支链连接运动副的中心位置所构成的图形可等效为等腰三角形△B₁B₂B₃和△A₁A₂A₃，其外接圆半径分别为r和R，A₁、A₂、A₃分别为静平台上的转动副中心，B₁、B₂、B₃分别为动平台上球副中心，L₁、L₂、L₃分别为直线伺服电缸所在支链长度。在静平台三角形外接圆圆心位置建立静坐标系J‒XYZ，X轴沿着 JA₁方向，Z轴垂直于静平台向上，Y轴由右手坐标系确定且方向与 A₂A₃平行；同理，在动平台三角形外接圆圆心位置建立坐标系D‒xyz，x轴沿着 DB₁方向，z轴垂直于动平台向上，y轴由右手坐标系确定且方向与 B₂B₃平行。

定义向量 JA_i （i=1,2,3）在静坐标系J‒XYZ下的位置矢量为 V_Ai， DB_i （i=1,2,3）在动坐标系D‒xyz下的位置矢量为 V_Bi。令动坐标系D‒xyz相对于静坐标系J‒XYZ的旋转变换矩阵为{}_D{}^J\mathit{R}， P 为动平台坐标系原点相对于静平台坐标系原点的位置矢量，\phi =\frac{\mathrm{\angle }{A}_{\mathrm{2}}J{A}_{\mathrm{3}}}{\mathrm{2}}，\theta =\frac{\mathrm{\angle }{B}_{\mathrm{2}}D{B}_{\mathrm{3}}}{\mathrm{2}}。

V_Ai 为：

\left\{\begin{array}{l}{\mathit{V}}_{A\mathrm{1}}=(R,\mathrm{0,0}),\\ {\mathit{V}}_{A\mathrm{2}}=(-R\mathrm{c}\mathrm{o}\mathrm{s}\phi ,R\mathrm{s}\mathrm{i}\mathrm{n}\phi ,\mathrm{0}),\\ {\mathit{V}}_{A\mathrm{3}}=(-R\mathrm{c}\mathrm{o}\mathrm{s}\phi ,-R\mathrm{s}\mathrm{i}\mathrm{n}\phi ,\mathrm{0})\end{array}\right.（2）

V_Bi 为：

\left\{\begin{array}{l}{\mathit{V}}_{B\mathrm{1}}=(r,\mathrm{0,0}),\\ {\mathit{V}}_{B\mathrm{2}}=(-r\mathrm{c}\mathrm{o}\mathrm{s}\theta ,r\mathrm{s}\mathrm{i}\mathrm{n}\theta ,\mathrm{0}),\\ {\mathit{V}}_{B\mathrm{3}}=(-r\mathrm{c}\mathrm{o}\mathrm{s}\theta ,-r\mathrm{s}\mathrm{i}\mathrm{n}\theta ,\mathrm{0})\end{array}\right.（3）

将动坐标系中位置矢量{\mathit{V}}_{Bi}转换到静坐标系中为：

{{}_D{}^J}_{Bi}={}_D{}^J\mathit{R}\cdot {\mathit{V}}_{Bi}+{\mathit{P}}_O（4）

式中：{\mathit{P}}_O为动平台坐标系原点D相对于静平台原点J的位置矢量，{\mathit{P}}_O=(x_D,y_D,z_D)；x_D、y_D、z_D分别为原点D在静坐标系中的x、y、z坐标值。动坐标系相对于静坐标系的旋转变换矩阵可表示为：

\begin{array}{l}{}_D{}^J\mathit{R}=\mathit{R}(Z,\gamma )\mathit{R}(Y,\beta )\mathit{R}(X,\alpha )=\\ \left[\begin{array}{ccc}\mathrm{c}\mathrm{o}\mathrm{s}\gamma \mathrm{c}\mathrm{o}\mathrm{s}\beta & \begin{array}{l}\mathrm{c}\mathrm{o}\mathrm{s}\gamma \mathrm{s}\mathrm{i}\mathrm{n}\beta \mathrm{s}\mathrm{i}\mathrm{n}\alpha -\\ \mathrm{s}\mathrm{i}\mathrm{n}\gamma \mathrm{c}\mathrm{o}\mathrm{s}\alpha \end{array}& \begin{array}{l}\mathrm{c}\mathrm{o}\mathrm{s}\gamma \mathrm{s}\mathrm{i}\mathrm{n}\beta \mathrm{c}\mathrm{o}\mathrm{s}\alpha +\\ \mathrm{s}\mathrm{i}\mathrm{n}\gamma \mathrm{s}\mathrm{i}\mathrm{n}\alpha \end{array}\\ \mathrm{s}\mathrm{i}\mathrm{n}\gamma \mathrm{c}\mathrm{o}\mathrm{s}\beta & \begin{array}{l}\mathrm{s}\mathrm{i}\mathrm{n}\gamma \mathrm{s}\mathrm{i}\mathrm{n}\beta \mathrm{s}\mathrm{i}\mathrm{n}\alpha +\\ \mathrm{c}\mathrm{o}\mathrm{s}\gamma \mathrm{c}\mathrm{o}\mathrm{s}\alpha \end{array}& \begin{array}{l}\mathrm{s}\mathrm{i}\mathrm{n}\gamma \mathrm{s}\mathrm{i}\mathrm{n}\beta \mathrm{c}\mathrm{o}\mathrm{s}\alpha -\\ \mathrm{c}\mathrm{o}\mathrm{s}\gamma \mathrm{s}\mathrm{i}\mathrm{n}\alpha \end{array}\\ -\mathrm{s}\mathrm{i}\mathrm{n}\beta & \mathrm{c}\mathrm{o}\mathrm{s}\beta \mathrm{s}\mathrm{i}\mathrm{n}\alpha & \mathrm{c}\mathrm{o}\mathrm{s}\beta \mathrm{c}\mathrm{o}\mathrm{s}\alpha \end{array}\right]\end{array}（5）

式中，\alpha、\beta、\gamma为姿态角。由式（4）可得球面副运动中心B₁、B₂、B₃在静坐标系中的位置矢量表示为：

\left\{\begin{array}{l}{}_{}{}^J{\mathit{V}}_{B\mathrm{1}}=\left[\begin{array}{c}{x}_D+r\mathrm{c}\mathrm{o}\mathrm{s}\beta \mathrm{c}\mathrm{o}\mathrm{s}\gamma \\ y_D+r\mathrm{c}\mathrm{o}\mathrm{s}\beta \mathrm{s}\mathrm{i}\mathrm{n}\gamma \\ z_D-r\mathrm{s}\mathrm{i}\mathrm{n}\beta \end{array}\right],\\ {}_{}{}^J{\mathit{V}}_{B\mathrm{2}}=\left[\begin{array}{c}\begin{array}{l}{x}_D-r\mathrm{s}\mathrm{i}\mathrm{n}\theta (\mathrm{c}\mathrm{o}\mathrm{s}\alpha \mathrm{s}\mathrm{i}\mathrm{n}\gamma -\mathrm{c}\mathrm{o}\mathrm{s}\gamma \mathrm{s}\mathrm{i}\mathrm{n}\alpha \mathrm{s}\mathrm{i}\mathrm{n}\beta )-\\ r\mathrm{c}\mathrm{o}\mathrm{s}\beta \mathrm{c}\mathrm{o}\mathrm{s}\gamma \mathrm{c}\mathrm{o}\mathrm{s}\theta \end{array}\\ \begin{array}{l}{y}_D+r\mathrm{s}\mathrm{i}\mathrm{n}\theta (\mathrm{c}\mathrm{o}\mathrm{s}\alpha \mathrm{c}\mathrm{o}\mathrm{s}\gamma +\mathrm{s}\mathrm{i}\mathrm{n}\gamma \mathrm{s}\mathrm{i}\mathrm{n}\alpha \mathrm{s}\mathrm{i}\mathrm{n}\beta )-\\ r\mathrm{c}\mathrm{o}\mathrm{s}\beta \mathrm{s}\mathrm{i}\mathrm{n}\gamma \mathrm{c}\mathrm{o}\mathrm{s}\theta \end{array}\\ z_D+r\mathrm{s}\mathrm{i}\mathrm{n}\beta \mathrm{c}\mathrm{o}\mathrm{s}\theta +r\mathrm{c}\mathrm{o}\mathrm{s}\beta \mathrm{s}\mathrm{i}\mathrm{n}\alpha \mathrm{s}\mathrm{i}\mathrm{n}\theta \end{array}\right],\\ {}_{}{}^J{\mathit{V}}_{B\mathrm{3}}=\left[\begin{array}{c}\begin{array}{l}{x}_D+r\mathrm{s}\mathrm{i}\mathrm{n}\theta (\mathrm{c}\mathrm{o}\mathrm{s}\alpha \mathrm{s}\mathrm{i}\mathrm{n}\gamma -\mathrm{c}\mathrm{o}\mathrm{s}\gamma \mathrm{s}\mathrm{i}\mathrm{n}\alpha \mathrm{s}\mathrm{i}\mathrm{n}\beta )-\\ r\mathrm{c}\mathrm{o}\mathrm{s}\beta \mathrm{c}\mathrm{o}\mathrm{s}\gamma \mathrm{c}\mathrm{o}\mathrm{s}\theta \end{array}\\ \begin{array}{l}{y}_D-r\mathrm{s}\mathrm{i}\mathrm{n}\theta (\mathrm{c}\mathrm{o}\mathrm{s}\alpha \mathrm{c}\mathrm{o}\mathrm{s}\gamma +\mathrm{s}\mathrm{i}\mathrm{n}\gamma \mathrm{s}\mathrm{i}\mathrm{n}\alpha \mathrm{s}\mathrm{i}\mathrm{n}\beta )-\\ r\mathrm{c}\mathrm{o}\mathrm{s}\beta \mathrm{s}\mathrm{i}\mathrm{n}\gamma \mathrm{c}\mathrm{o}\mathrm{s}\theta \end{array}\\ z_D+r\mathrm{s}\mathrm{i}\mathrm{n}\beta \mathrm{c}\mathrm{o}\mathrm{s}\theta -r\mathrm{c}\mathrm{o}\mathrm{s}\beta \mathrm{s}\mathrm{i}\mathrm{n}\alpha \mathrm{s}\mathrm{i}\mathrm{n}\theta \end{array}\right]\end{array}\right.（6）

由式（2）和（6）得连接动静平台的各支链长度L₁、L₂、L₃分别为：

\left\{\begin{array}{c}\begin{array}{l}{L}_{\mathrm{1}}^{\mathrm{2}}=((x_D{+r\mathrm{c}\mathrm{o}\mathrm{s}\beta \mathrm{c}\mathrm{o}\mathrm{s}\gamma -R)}^{\mathrm{2}}+(y_D{+r\mathrm{c}\mathrm{o}\mathrm{s}\beta \mathrm{s}\mathrm{i}\mathrm{n}\gamma )}^{\mathrm{2}}+\\ (z_D{-r\mathrm{s}\mathrm{i}\mathrm{n}\beta )}^{\mathrm{2}}),\end{array}\\ \begin{array}{l}{L}_{\mathrm{2}}^{\mathrm{2}}=((x_D-r\mathrm{s}\mathrm{i}\mathrm{n}\theta (\mathrm{c}\mathrm{o}\mathrm{s}\alpha \mathrm{s}\mathrm{i}\mathrm{n}\gamma -\mathrm{c}\mathrm{o}\mathrm{s}\gamma \mathrm{s}\mathrm{i}\mathrm{n}\alpha \mathrm{s}\mathrm{i}\mathrm{n}\beta )-r\mathrm{c}\mathrm{o}\mathrm{s}\beta \cdot \\ {\mathrm{c}\mathrm{o}\mathrm{s}\gamma \mathrm{c}\mathrm{o}\mathrm{s}\theta +R\mathrm{c}\mathrm{o}\mathrm{s}\phi )}^{\mathrm{2}}+(y_D+r\mathrm{s}\mathrm{i}\mathrm{n}\theta (\mathrm{c}\mathrm{o}\mathrm{s}\alpha \mathrm{c}\mathrm{o}\mathrm{s}\gamma +\\ {\mathrm{s}\mathrm{i}\mathrm{n}\gamma \mathrm{s}\mathrm{i}\mathrm{n}\alpha \mathrm{s}\mathrm{i}\mathrm{n}\beta )-r\mathrm{c}\mathrm{o}\mathrm{s}\beta \mathrm{s}\mathrm{i}\mathrm{n}\gamma \mathrm{c}\mathrm{o}\mathrm{s}\theta -R\mathrm{s}\mathrm{i}\mathrm{n}\phi )}^{\mathrm{2}}+\\ (z_D{+r\mathrm{s}\mathrm{i}\mathrm{n}\beta \mathrm{c}\mathrm{o}\mathrm{s}\theta +r\mathrm{c}\mathrm{o}\mathrm{s}\beta \mathrm{s}\mathrm{i}\mathrm{n}\alpha \mathrm{s}\mathrm{i}\mathrm{n}\theta )}^{\mathrm{2}}),\end{array}\\ \begin{array}{l}{L}_{\mathrm{3}}^{\mathrm{2}}=((x_D+r\mathrm{s}\mathrm{i}\mathrm{n}\theta (\mathrm{c}\mathrm{o}\mathrm{s}\alpha \mathrm{s}\mathrm{i}\mathrm{n}\gamma -\mathrm{c}\mathrm{o}\mathrm{s}\gamma \mathrm{s}\mathrm{i}\mathrm{n}\alpha \mathrm{s}\mathrm{i}\mathrm{n}\beta )-r\mathrm{c}\mathrm{o}\mathrm{s}\beta \cdot \\ {\mathrm{c}\mathrm{o}\mathrm{s}\gamma \mathrm{c}\mathrm{o}\mathrm{s}\theta +R\mathrm{c}\mathrm{o}\mathrm{s}\phi )}^{\mathrm{2}}+(y_D-r\mathrm{s}\mathrm{i}\mathrm{n}\theta (\mathrm{c}\mathrm{o}\mathrm{s}\alpha \mathrm{c}\mathrm{o}\mathrm{s}\gamma +\\ {\mathrm{s}\mathrm{i}\mathrm{n}\gamma \mathrm{s}\mathrm{i}\mathrm{n}\alpha \mathrm{s}\mathrm{i}\mathrm{n}\beta )-r\mathrm{c}\mathrm{o}\mathrm{s}\beta \mathrm{s}\mathrm{i}\mathrm{n}\gamma \mathrm{c}\mathrm{o}\mathrm{s}\theta +R\mathrm{s}\mathrm{i}\mathrm{n}\phi )}^{\mathrm{2}}+\\ (z_D{+r\mathrm{s}\mathrm{i}\mathrm{n}\beta \mathrm{c}\mathrm{o}\mathrm{s}\theta +r\mathrm{c}\mathrm{o}\mathrm{s}\beta \mathrm{s}\mathrm{i}\mathrm{n}\alpha \mathrm{s}\mathrm{i}\mathrm{n}\theta )}^{\mathrm{2}})\end{array}\end{array}\right.（7）

由上述自由度分析可知，3‒RPS型并联机构有3个独立自由度，而式（7）中求解支链长度（伺服电缸直线位移）的方程组有6个自由度的位姿变量无法直接求得。对并联机构的约束条件进行分析可知，该机构的支链的特性导致动平台上的3个球副中心点B₁、B₂、B₃分别只能在3个转动副运动平面内运动，则有：

\left\{\begin{array}{l}{y}_{B_{\mathrm{1}}}=\mathrm{0},\\ y_{B_{\mathrm{2}}}=-\mathrm{t}\mathrm{a}\mathrm{n}\varphi \cdot (x_{B_{\mathrm{2}}}+R\mathrm{c}\mathrm{o}\mathrm{s}\phi )+R\mathrm{s}\mathrm{i}\mathrm{n}\phi ,\\ y_{B_{\mathrm{3}}}=\mathrm{t}\mathrm{a}\mathrm{n}\varphi \cdot (x_{B_{\mathrm{2}}}+R\mathrm{c}\mathrm{o}\mathrm{s}\phi )-R\mathrm{s}\mathrm{i}\mathrm{n}\phi \end{array}\right.（8）

式中，\varphi为转动副运动平面在xy平面上的投影与静坐标系x轴的夹角。

将该约束条件联合式（6）得到以下3个约束方程：

y_D+r\mathrm{c}\mathrm{o}\mathrm{s}\beta \mathrm{s}\mathrm{i}\mathrm{n}\gamma =\mathrm{0}（9）

\begin{array}{l}{y}_D+r\mathrm{s}\mathrm{i}\mathrm{n}\theta (\mathrm{c}\mathrm{o}\mathrm{s}\alpha \mathrm{c}\mathrm{o}\mathrm{s}\gamma +\mathrm{s}\mathrm{i}\mathrm{n}\gamma \mathrm{s}\mathrm{i}\mathrm{n}\alpha \mathrm{s}\mathrm{i}\mathrm{n}\beta )-r\mathrm{c}\mathrm{o}\mathrm{s}\beta \mathrm{s}\mathrm{i}\mathrm{n}\gamma \cdot \\ \mathrm{c}\mathrm{o}\mathrm{s}\theta -(-\mathrm{t}\mathrm{a}\mathrm{n}\varphi (x_D-r\mathrm{s}\mathrm{i}\mathrm{n}\theta (\mathrm{c}\mathrm{o}\mathrm{s}\alpha \mathrm{s}\mathrm{i}\mathrm{n}\gamma -\mathrm{c}\mathrm{o}\mathrm{s}\gamma \mathrm{s}\mathrm{i}\mathrm{n}\alpha \cdot \\ \mathrm{s}\mathrm{i}\mathrm{n}\beta )-r\mathrm{c}\mathrm{o}\mathrm{s}\beta \mathrm{c}\mathrm{o}\mathrm{s}\gamma \mathrm{c}\mathrm{o}\mathrm{s}\theta +R\mathrm{c}\mathrm{o}\mathrm{s}\phi )+R\mathrm{s}\mathrm{i}\mathrm{n}\phi )=\mathrm{0}\end{array}（10）

\begin{array}{l}{y}_D-r\mathrm{s}\mathrm{i}\mathrm{n}\theta (\mathrm{c}\mathrm{o}\mathrm{s}\alpha \mathrm{c}\mathrm{o}\mathrm{s}\gamma +\mathrm{s}\mathrm{i}\mathrm{n}\gamma \mathrm{s}\mathrm{i}\mathrm{n}\alpha \mathrm{s}\mathrm{i}\mathrm{n}\beta )-r\mathrm{c}\mathrm{o}\mathrm{s}\beta \mathrm{s}\mathrm{i}\mathrm{n}\gamma \cdot \\ \mathrm{c}\mathrm{o}\mathrm{s}\theta -(\mathrm{t}\mathrm{a}\mathrm{n}\varphi (x_D+r\mathrm{s}\mathrm{i}\mathrm{n}\theta (\mathrm{c}\mathrm{o}\mathrm{s}\alpha \mathrm{s}\mathrm{i}\mathrm{n}\gamma -\mathrm{c}\mathrm{o}\mathrm{s}\gamma \mathrm{s}\mathrm{i}\mathrm{n}\alpha \cdot \\ \mathrm{s}\mathrm{i}\mathrm{n}\beta )-r\mathrm{c}\mathrm{o}\mathrm{s}\beta \mathrm{c}\mathrm{o}\mathrm{s}\gamma \mathrm{c}\mathrm{o}\mathrm{s}\theta +R\mathrm{c}\mathrm{o}\mathrm{s}\phi )-R\mathrm{s}\mathrm{i}\mathrm{n}\phi )=\mathrm{0}\end{array}（11）

联立式（9）～（11）可得x_D、y_D、γ关于变量α和β的表达式，由于表达式较为复杂，在此不再给出。因此，3‒RPS型并联机构运动学的运动学逆解方法是，只需给定笛卡尔空间动平台坐标原点相对于静平台坐标原点的Z轴坐标值z_D、动平台绕静坐标系X轴和Y轴的转角α和β，根据机构约束方程式（9）、（10）、（11）就可得x_D、y_D、γ，然后，将以上6个已知参数代入到式（7）中即可求得支链长度L₁、L₂、L₃。

3‒RPS‒P型混联机构在并联机构的基础上增加一个x方向的自由度，也即沿着动平台坐标系x轴做平移运动，因此，该混联机构动平台的移动副上所安装的马鞍的运动学逆解方法是，其相对于动平台只有位置发生变化，姿态仍相同。动平台坐标系相对于静平台坐标系的齐次变换矩阵为：

{}_D{}^J\mathit{T}=\left[\begin{array}{cccc}\begin{array}{l}\mathrm{c}\mathrm{o}\mathrm{s}\gamma \cdot \\ \mathrm{c}\mathrm{o}\mathrm{s}\beta \end{array}& \begin{array}{l}\mathrm{c}\mathrm{o}\mathrm{s}\gamma \mathrm{s}\mathrm{i}\mathrm{n}\beta \mathrm{s}\mathrm{i}\mathrm{n}\alpha -\\ \mathrm{s}\mathrm{i}\mathrm{n}\gamma \mathrm{c}\mathrm{o}\mathrm{s}\alpha \end{array}& \begin{array}{l}\mathrm{c}\mathrm{o}\mathrm{s}\gamma \mathrm{s}\mathrm{i}\mathrm{n}\beta \mathrm{c}\mathrm{o}\mathrm{s}\alpha +\\ \mathrm{s}\mathrm{i}\mathrm{n}\gamma \mathrm{s}\mathrm{i}\mathrm{n}\alpha \end{array}& x_D\\ \begin{array}{l}\mathrm{s}\mathrm{i}\mathrm{n}\gamma \cdot \\ \mathrm{c}\mathrm{o}\mathrm{s}\beta \end{array}& \begin{array}{l}\mathrm{s}\mathrm{i}\mathrm{n}\gamma \mathrm{s}\mathrm{i}\mathrm{n}\beta \mathrm{s}\mathrm{i}\mathrm{n}\alpha +\\ \mathrm{c}\mathrm{o}\mathrm{s}\gamma \mathrm{c}\mathrm{o}\mathrm{s}\alpha \end{array}& \begin{array}{l}\mathrm{s}\mathrm{i}\mathrm{n}\gamma \mathrm{s}\mathrm{i}\mathrm{n}\beta \mathrm{c}\mathrm{o}\mathrm{s}\alpha -\\ \mathrm{c}\mathrm{o}\mathrm{s}\gamma \mathrm{s}\mathrm{i}\mathrm{n}\alpha \end{array}& y_D\\ -\mathrm{s}\mathrm{i}\mathrm{n}\beta & \mathrm{c}\mathrm{o}\mathrm{s}\beta \mathrm{s}\mathrm{i}\mathrm{n}\alpha & \mathrm{c}\mathrm{o}\mathrm{s}\beta \mathrm{c}\mathrm{o}\mathrm{s}\alpha & z_D\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right]（12）

因此，马鞍中心的位姿矩阵为：

{\mathit{T}}_M={}_D{}^J\mathit{T}\cdot \mathrm{T}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}(x,L_{\mathrm{4}})\cdot \mathrm{T}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}(z,H)（13）

式中，L₄为动平台上直线伺服电缸的伸长量，H为马鞍中心在动平台坐标系下的高度。

3‒RPS‒P型混联机构有4个自由度，其运动学逆解方法是给定马鞍中心(x,z,α,β)作为已知参数，因马鞍姿态与动平台姿态相同，由式（9）～（11）求解γ和动平台坐标系原点坐标x_D、y_D。由式（12）、（13）求解L_{\mathrm{4}}和z_D：

\begin{array}{l}{L}_{\mathrm{4}}=(x-x_D-H(\mathrm{c}\mathrm{o}\mathrm{s}\gamma \mathrm{s}\mathrm{i}\mathrm{n}\beta \mathrm{c}\mathrm{o}\mathrm{s}\alpha +\mathrm{s}\mathrm{i}\mathrm{n}\gamma \mathrm{s}\mathrm{i}\mathrm{n}\alpha ))/\\ (\mathrm{c}\mathrm{o}\mathrm{s}\gamma \mathrm{c}\mathrm{o}\mathrm{s}\beta )\end{array}（14）

z_D=z-H\cdot \mathrm{c}\mathrm{o}\mathrm{s}\beta \mathrm{c}\mathrm{o}\mathrm{s}\alpha +L_{\mathrm{4}}\cdot \mathrm{s}\mathrm{i}\mathrm{n}\beta（15）

将z_D、α、β按照上述并联机构逆解求解方法进行计算，得到各支链长度L₁、L₂、L₃，至此混联机构逆解求解完毕。

2.2 运动学正解分析

3‒RPS‒P型混联机构中有并联部分，分析该并联机构的正运动学的解，需要求解式（7）、（9）、（10）、（11）中6个耦合的非线性方程组成的方程组式（16），但很难求解出该方程组的解析解。

\left\{\begin{array}{l}{F}_{\mathrm{1}}=\mathrm{0},F_{\mathrm{2}}=\mathrm{0},F_{\mathrm{3}}=\mathrm{0},\\ F_{\mathrm{4}}=\mathrm{0},F_{\mathrm{5}}=\mathrm{0},F_{\mathrm{6}}=\mathrm{0}\end{array}\right.（16）

因此，采用牛顿迭代法求解该方程组的数值解，具体如下：

{\mathit{P}}^{(K+\mathrm{1})}={\mathit{P}}^{(K)}-{\left[{\mathit{F}}^{\text{'}}({\mathit{P}}^{(K)})\right]}^{-\mathrm{1}}\mathit{F}({\mathit{P}}^{(K)}),K=\mathrm{0,1},\cdots（17）

式中，\mathit{F}={\left[F_{\mathrm{1}},F_{\mathrm{2}},F_{\mathrm{3}},F_{\mathrm{4}},F_{\mathrm{5}},F_{\mathrm{6}}\right]}^{\mathrm{T}}，{\mathit{P}}^{(K)}为第K次迭代的位姿，{\mathit{P}}^{(K)}={\left[x_D,y_D,z_D,\alpha ,\beta ,\gamma \right]}^{\mathrm{T}}，向量函数\mathit{F}的雅可比矩阵：

{\mathit{F}}^{\text{'}}=\left[\begin{array}{cccccc}\frac{\partial F_{\mathrm{1}}}{\partial x_D}& \frac{\partial F_{\mathrm{1}}}{\partial y_D}& \frac{\partial F_{\mathrm{1}}}{\partial z_D}& \frac{\partial F_{\mathrm{1}}}{\partial \alpha }& \frac{\partial F_{\mathrm{1}}}{\partial \beta }& \frac{\partial F_{\mathrm{1}}}{\partial \gamma }\\ \frac{\partial F_{\mathrm{2}}}{\partial x_D}& \frac{\partial F_{\mathrm{2}}}{\partial y_D}& \frac{\partial F_{\mathrm{2}}}{\partial z_D}& \frac{\partial F_{\mathrm{2}}}{\partial \alpha }& \frac{\partial F_{\mathrm{2}}}{\partial \beta }& \frac{\partial F_{\mathrm{2}}}{\partial \gamma }\\ \frac{\partial F_{\mathrm{3}}}{\partial x_D}& \frac{\partial F_{\mathrm{3}}}{\partial y_D}& \frac{\partial F_{\mathrm{3}}}{\partial z_D}& \frac{\partial F_{\mathrm{3}}}{\partial \alpha }& \frac{\partial F_{\mathrm{3}}}{\partial \beta }& \frac{\partial F_{\mathrm{3}}}{\partial \gamma }\\ \frac{\partial F_{\mathrm{4}}}{\partial x_D}& \frac{\partial F_{\mathrm{4}}}{\partial y_D}& \frac{\partial F_{\mathrm{4}}}{\partial z_D}& \frac{\partial F_{\mathrm{4}}}{\partial \alpha }& \frac{\partial F_{\mathrm{4}}}{\partial \beta }& \frac{\partial F_{\mathrm{4}}}{\partial \gamma }\\ \frac{\partial F_{\mathrm{5}}}{\partial x_D}& \frac{\partial F_{\mathrm{5}}}{\partial y_D}& \frac{\partial F_{\mathrm{5}}}{\partial z_D}& \frac{\partial F_{\mathrm{5}}}{\partial \alpha }& \frac{\partial F_{\mathrm{5}}}{\partial \beta }& \frac{\partial F_{\mathrm{5}}}{\partial \gamma }\\ \frac{\partial F_{\mathrm{6}}}{\partial x_D}& \frac{\partial F_{\mathrm{6}}}{\partial y_D}& \frac{\partial F_{\mathrm{6}}}{\partial z_D}& \frac{\partial F_{\mathrm{6}}}{\partial \alpha }& \frac{\partial F_{\mathrm{6}}}{\partial \beta }& \frac{\partial F_{\mathrm{6}}}{\partial \gamma }\end{array}\right]（18）

根据式（17）编写计算机迭代计算程序，计算流程如图5所示。

根据上述牛顿迭代计算方法可以得到动平台坐标系相对于静平台坐标系的位姿{}_D{}^J\mathit{T}，再由式（13）计算可得马鞍中心相对于静平台坐标系的位姿，至此混联机构正解求解完毕。

2.3 速度分析

在机器人机构中，笛卡尔空间末端运动速度与关节空间速度之间的映射关系可以通过雅可比矩阵实现^[26]。在3‒RPS‒P型混联机构中以马鞍中心为研究对象，研究各个直线伺服电缸与马鞍中心运动速度的关系，马鞍中心运动速度与关节直线伺服运动速度的表达式为：

\left[\begin{array}{l}\mathit{v}\\ \mathit{w}\end{array}\right]=\mathit{J}\left[\dot{\mathit{L}}\right]（19）

式中， v 和 w 分别为马鞍中心的线速度与角速度矢量， J 为速度雅可比矩阵，\dot{\mathit{L}}={\left[\begin{array}{cccc}{\dot{L}}_{\mathrm{1}}& {\dot{L}}_{\mathrm{2}}& {\dot{L}}_{\mathrm{3}}& {\dot{L}}_{\mathrm{4}}\end{array}\right]}^{\mathrm{T}}为关节的线速度矢量。

对式（7）和（14）两边求导后写成矩阵形式为：

\left[\begin{array}{cccc}\mathrm{2}{L}_{\mathrm{1}}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{2}{L}_{\mathrm{2}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{2}{L}_{\mathrm{3}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right]\left[\begin{array}{l}{\dot{L}}_{\mathrm{1}}\\ {\dot{L}}_{\mathrm{2}}\\ {\dot{L}}_{\mathrm{3}}\\ {\dot{L}}_{\mathrm{4}}\end{array}\right]=\left[\begin{array}{cccc}{K}_{\mathrm{11}}& K_{\mathrm{12}}& K_{\mathrm{13}}& K_{\mathrm{14}}\\ K_{\mathrm{21}}& K_{\mathrm{22}}& K_{\mathrm{23}}& K_{\mathrm{24}}\\ K_{\mathrm{31}}& K_{\mathrm{32}}& K_{\mathrm{33}}& K_{\mathrm{34}}\\ K_{\mathrm{41}}& K_{\mathrm{42}}& K_{\mathrm{43}}& K_{\mathrm{44}}\end{array}\right]\left[\begin{array}{l}\dot{x}\\ \dot{z}\\ \dot{\alpha }\\ \dot{\beta }\end{array}\right]（20）

式中，由于求导后表达式较为复杂，K_ij 为求导后的表达式简写。因此，该机构的雅可比矩阵可表示为：

\mathit{J}={\left[\begin{array}{cccc}{K}_{\mathrm{11}}& K_{\mathrm{12}}& K_{\mathrm{13}}& K_{\mathrm{14}}\\ K_{\mathrm{21}}& K_{\mathrm{22}}& K_{\mathrm{23}}& K_{\mathrm{24}}\\ K_{\mathrm{31}}& K_{\mathrm{32}}& K_{\mathrm{33}}& K_{\mathrm{34}}\\ K_{\mathrm{41}}& K_{\mathrm{42}}& K_{\mathrm{43}}& K_{\mathrm{44}}\end{array}\right]}^{-\mathrm{1}}\left[\begin{array}{cccc}\mathrm{2}{L}_{\mathrm{1}}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{2}{L}_{\mathrm{2}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{2}{L}_{\mathrm{3}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right]（21）

2.4 工作空间分析

3‒RPS‒P型混联机构的x方向移动自由度来源于动平台上的移动副，该移动副运动灵活，只受直线伺服电缸行程的限制，因此主要研究并联机构部分的工作空间。工作空间是并联机器人机构动平台坐标系原点可达位置的集合，是机构运动学性能的重要指标，也是样机制造和现场作业的关键参考之一。

该机构并联部分运动范围的影响因素主要有连杆长度的限制和运动副转角的限制。连杆长度主要取决于直线伺服电缸的初始长度和伸长行程，该并联机构选取3个相同的直线伺服电缸，电缸工作时连杆的长度变化范围为380 mm⩽L⩽480 mm；转动副将静平台与连杆铰接，球副将动平台与连杆铰接，机构运动时工作空间主要受到球副转角限制的影响。球副机构选用鱼眼球轴承，如图6所示。并联机构具有z、α、β 3个自由度，根据鱼眼球轴承安装布局方式可知，并联机构运动时工作空间受到绕x轴转角α的影响最大，即图6中鱼眼球轴承端面与聚氨酯垫片端面的夹角的限制，记为∠CED。

由于机构运动时三支链鱼眼球轴承端面与聚氨酯垫片端面所形成的角度不同，分别用∠J₁、∠J₂、∠J₃表示三支链运动时所形成的角度。图6中，\mathit{E}\mathit{D}为动平台所在平面的法向量在yz平面的投影，记为 n，具体为：

\mathit{n}=(\mathrm{s}\mathrm{i}\mathrm{n}\gamma \mathrm{s}\mathrm{i}\mathrm{n}\beta \mathrm{c}\mathrm{o}\mathrm{s}\alpha -\mathrm{c}\mathrm{o}\mathrm{s}\gamma \mathrm{s}\mathrm{i}\mathrm{n}\alpha ,\mathrm{c}\mathrm{o}\mathrm{s}\beta \mathrm{c}\mathrm{o}\mathrm{s}\alpha )（22）

另外，图6中， EC 为各支链向量{\mathit{L}}_{L\mathrm{1}}、{\mathit{L}}_{L\mathrm{2}}、{\mathit{L}}_{L\mathrm{3}}在yz平面的投影，其中，{\mathit{L}}_{L\mathrm{1}}、{\mathit{L}}_{L\mathrm{2}}、{\mathit{L}}_{L\mathrm{3}}分别为：

\left\{\begin{array}{l}{\mathit{L}}_{L\mathrm{1}}=(y_D+r\mathrm{c}\mathrm{o}\mathrm{s}\beta \mathrm{s}\mathrm{i}\mathrm{n}\gamma ,z_D-r\mathrm{s}\mathrm{i}\mathrm{n}\beta );\\ {\mathit{L}}_{L\mathrm{2}}=(y_D+r\mathrm{s}\mathrm{i}\mathrm{n}\theta (\mathrm{c}\mathrm{o}\mathrm{s}\alpha \mathrm{c}\mathrm{o}\mathrm{s}\gamma +\mathrm{s}\mathrm{i}\mathrm{n}\alpha \mathrm{s}\mathrm{i}\mathrm{n}\beta \mathrm{s}\mathrm{i}\mathrm{n}\gamma )-\\ r\mathrm{c}\mathrm{o}\mathrm{s}\beta \mathrm{s}\mathrm{i}\mathrm{n}\gamma \mathrm{c}\mathrm{o}\mathrm{s}\theta -R\mathrm{s}\mathrm{i}\mathrm{n}\phi ,z_D+r\mathrm{s}\mathrm{i}\mathrm{n}\beta \mathrm{c}\mathrm{o}\mathrm{s}\theta +\\ r\mathrm{c}\mathrm{o}\mathrm{s}\beta \mathrm{s}\mathrm{i}\mathrm{n}\alpha \mathrm{s}\mathrm{i}\mathrm{n}\theta );\\ {\mathit{L}}_{L\mathrm{3}}=(y_D-r\mathrm{s}\mathrm{i}\mathrm{n}\theta (\mathrm{c}\mathrm{o}\mathrm{s}\alpha \mathrm{c}\mathrm{o}\mathrm{s}\gamma +\mathrm{s}\mathrm{i}\mathrm{n}\alpha \mathrm{s}\mathrm{i}\mathrm{n}\beta \mathrm{s}\mathrm{i}\mathrm{n}\gamma )-\\ r\mathrm{c}\mathrm{o}\mathrm{s}\beta \mathrm{s}\mathrm{i}\mathrm{n}\gamma \mathrm{c}\mathrm{o}\mathrm{s}\theta +R\mathrm{s}\mathrm{i}\mathrm{n}\phi ,z_D+r\mathrm{s}\mathrm{i}\mathrm{n}\beta \mathrm{c}\mathrm{o}\mathrm{s}\theta -\\ r\mathrm{c}\mathrm{o}\mathrm{s}\beta \mathrm{s}\mathrm{i}\mathrm{n}\alpha \mathrm{s}\mathrm{i}\mathrm{n}\theta )\end{array}\right.（23）

由式（22）、（23）可求得∠J₁、∠J₂、∠J₃分别为：

\left\{\begin{array}{l}\mathrm{\angle }{J}_{\mathrm{1}}=\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{c}\mathrm{o}\mathrm{s}((\mathit{n}\cdot {\mathit{L}}_{L\mathrm{1}})/(\left|\mathit{n}\right|\cdot \left|{\mathit{L}}_{L\mathrm{1}}\right|)),\\ \mathrm{\angle }{J}_{\mathrm{2}}=\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{c}\mathrm{o}\mathrm{s}((\mathit{n}\cdot {\mathit{L}}_{L\mathrm{2}})/(\left|\mathit{n}\right|\cdot \left|{\mathit{L}}_{L\mathrm{2}}\right|)),\\ \mathrm{\angle }{J}_{\mathrm{3}}=\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{c}\mathrm{o}\mathrm{s}((\mathit{n}\cdot {\mathit{L}}_{L\mathrm{3}})/(\left|\mathit{n}\right|\cdot \left|{\mathit{L}}_{L\mathrm{3}}\right|))\end{array}\right.（24）

将电缸工作时连杆长度的变化范围进行切割划分，切割宽度为Δz；使用MATLAB编写3个连杆长度的循环遍历代码，采用图5的牛顿迭代的正运动学求解方法求出动平台中心位姿；根据得到的位姿计算∠J₁、∠J₂、∠J₃，并判断3个角度中的最大角度是否超出机构角度限制，未超出的位姿即为并联机构的可达空间，绘制并联机构的工作空间如图7所示。

从图7可以看出：整体工作空间比较规整，内部没有空洞，连续性较好，说明该机构的运动性能良好。动平台在z轴的运动空间区间为379～479 mm，绕x轴转角α范围为-8.5°～8.5°，绕y轴转角β范围为-10.5°～10.5°。

2.5 动力学分析

动力学分析是指研究机构运动与力之间的关系，可分析机构的动态响应特性，以提高控制精度和运动效率。基于拉格朗日法建立并联机构支链及动平台的动势能方程，忽略机构各运动副的质量和摩擦力的影响，假设并联机构各构件的质量分布均匀。

根据拉格朗日法有：

\frac{\mathrm{d}}{\mathrm{d}t}(\frac{\partial L}{\partial \dot{\mathit{q}}})-\frac{\partial L}{\partial \mathit{q}}=\tau（25）

式中：L为拉格朗日函数，L=E-U，其中，E和U分别为系统动能和势能； q 为广义坐标， q=(z, α, β)^T；\mathit{\tau }为对应广义坐标 q 的非保守广义力，\mathit{\tau }=({\tau }_{\mathrm{1}},{\tau }_{\mathrm{2}},{\tau }_{\mathrm{3}}{)}^{\mathrm{T}}。

根据动平台的线速度{\mathit{v}}_{\mathrm{p}}和角速度{\mathit{\omega }}_{\mathrm{p}}，可以推导出动平台的总动能为：

E_{\mathrm{p}}=\frac{\mathrm{1}}{\mathrm{2}}m{\mathit{v}}_{\mathrm{p}}^{\mathrm{T}}{\mathit{v}}_{\mathrm{p}}+\frac{\mathrm{1}}{\mathrm{2}}{\mathit{\omega }}_{\mathrm{p}}^{\mathrm{T}}{\mathit{I}}_{\mathrm{p}}{\mathit{\omega }}_{\mathrm{p}}（26）

式中，m为动平台的质量，{\mathit{I}}_{\mathrm{p}}为动平台相对于固定坐标系的转动惯量。

根据各支链的线速度{\dot{\mathit{l}}}_i和角速度{\mathit{\omega }}_i，推导出机构支链的总动能为：

E_{\mathrm{r}}={\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{3}}}(\frac{\mathrm{1}}{\mathrm{2}}{m}_{\mathrm{1}}{\mathit{\omega }}_i^{\mathrm{T}}{\mathit{I}}_{\mathrm{r}i}{\mathit{\omega }}_i+\frac{\mathrm{1}}{\mathrm{2}}{m}_{\mathrm{12}}{\dot{\mathit{l}}}_i^{\mathrm{T}}{\dot{\mathit{l}}}_i)（27）

式中，m_{\mathrm{1}}为每条支链的总质量，m_{\mathrm{12}}为伸缩杆质量，{\mathit{I}}_{\mathrm{r}i}为第i个驱动支链相对于基坐标系J‒XYZ的惯量矩阵。

由式（26）和（27）求出并联机构总动能E=E_p+E_r。

在基坐标系J‒XYZ下，选择原点J和面XJY分别为零势点和零势面，重力加速度为 g。忽略不计构件的弹性和摩擦，动平台势能U_P和连杆势能U_R分别为：

U_{\mathrm{P}}=m\mathit{g}{D}_Z（28）

U_{\mathrm{R}}={{\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{3}}}m}_{\mathrm{1}}\mathit{g}{G}_{iz}（29）

式（28）～（29）中，D_Z和G_{iz}分别为动平台和各连杆在Z分量的高度。

并联机构系统的总势能为：

U=U_{\mathrm{P}}+U_{\mathrm{R}}（30）

综合式（25）、（27）和（30）可推导出并联机构的空间动力学模型，由于具体的表达式较为复杂，在此不再给出。

3 算例仿真

3.1 运动学算例分析

取两组算例初步验证上述运动学分析的正确性。算例1的马鞍中心位姿 P₁=[x=‒100 mm, z=620 mm, α=0°, β=0°]，算例2的马鞍中心位姿 P₂=[x=‒50 mm, z=630 mm, α=6°, β=8°]。将上述运动学逆解公式，使用MATLAB编写代码，将P₁和P₂点位姿作为输入，计算得到各个驱动支链长度如表1所示。

将表1逆解算例的各支链长度作为正解的输入，代入到MATLAB编写的牛顿迭代正解算法中，取初值 P⁽⁰⁾=[30,0,400,0,0,0]，计算迭代过程结果如表2所示。

由表2可以看出，本文的迭代方法快速收敛且误差较小，符合实际控制中对计算速度的要求。对比表1、2的算例1、2可知，取逆解计算得到的支链长度作为正解输入，经过迭代计算后得到的结果与逆解所选取的马鞍中心位姿 P₁和 P₂近乎相同，仅在x方向上有0.01 mm的计算误差，z、α、β与选取值完全相同，且计算精度较高。正解与逆解结果初步验证了本文解算方法的正确性。

另外，在SolidWorks软件中按照表1中算例1和2计算所得的各支链长度设定模型中各个直线伺服电缸的长度，建立静平台坐标系与马鞍中心坐标系，分别测量马鞍中心点在算例1和2情况下静平台坐标系中的坐标值，如图8和9所示。

由图8和9可以看出，SolidWorks软件中坐标测量值与正逆解算法计算值一致，进一步验证了运动学解算的正确性。

3.2 运动仿真

3‒RPS‒P型混联康复机器人采用坐姿模拟骑马过程中马匹的慢步、快步、跑步，以及跳跃等动作，研究马匹运动过程中马鞍中心的运动轨迹如图10所示。

以马匹慢步运动过程中马鞍中心运动轨迹为例，采用基于MATLAB-SolidWorks Motion的联合仿真平台对慢步运动过程中马鞍中心部分运动轨迹进行仿真分析。为提高仿真效率在SolidWorks中建立混联机构的简易模型，配置好各个零部件材料属性，在SolidWorks Motion中添加机构约束，并在马鞍上添加方向竖直向下的250 kg载荷，设置机构的4个直线伺服电缸的线性马达驱动。利用 MATLAB 对混联机构的位置和姿态在笛卡尔空间进行轨迹规划，将规划后的数据代入到运动学逆解方程中求出各驱动支链长度与时间的数值序列，再将该数值序列导入到 SolidWorks Motion 中进行运动仿真^[27]，仿真流程图如图11所示。

设马鞍中心从初始位姿 P₁=(29.267 6 mm,572.920 1 mm,0°,‒0.014 0°)运动到终止位姿 P₂=(‒60 mm,650 mm,6°,6°)。对位置x和z采用直线路径规划，以最大加速度A_max=300 mm/s²、最大速度V_max=100 mm/s的梯形加速度规律插补直线运动路径（图12（a）），耗时1.681 s。对α、β以A_max=30 mm/s²、V_max=10 mm/s的梯形加速度规律插补姿态角度运动（图12（b）），耗时1.316 s。位置和姿态插补均为等时插补，插补时间间隔为1 ms，为保证4轴控制的同步性，对位姿和姿态规划采用等时算法插补处理，处理后的位置和姿态规划时间相同。

给出插补后的位置直线规划路径和规划速度曲线，如图13（a）、（b）所示，以便和SolidWorks Motion模块中的仿真路径和速度曲线对比。

将上述经过梯形加速度规律插补的笛卡尔空间位姿代入混联机构逆解公式，计算得到并联机构的各个驱动支链长度、驱动伸长量L₄与时间的曲线如图14所示。

由图14可以看出，位置和姿态规划经过等时处理后，各个直线伺服驱动同步运动且时间相同，运动总耗时1.681 s，经过操作空间到关节空间的映射后所得到的驱动位移曲线连续且平滑过渡。

将图14中得到的驱动位移曲线代入上述SolidWorks Motion模块设置的线性驱动马达，进行运动仿真分析，经过模块计算后得到机构末端仿真轨迹如图15（在静坐标系xz平面的投影）所示，末端位置仿真速度如图16（x和z方向速度）所示。

由图15可见，机构末端仿真运动轨迹连续无跳动，在静坐标系xz平面下的投影为一条直线段，起始点P₁的位姿和终止点P₂的位姿与规划的位置相同。对比图15、13（a）可知，在MATLAB中的规划运动路径与SolidWorks Motion模块中仿真运动路径相同。由图16可见，机构仿真运动时末端速度曲线连续平滑为S形曲线。对比图16、13（b）可知，将图16中x和z两方向运动速度合成后与图13（b）中规划速度相同。因此，规划路径和规划速度均与仿真中的路径和速度相同，进一步说明了本文的运动学方法和轨迹规划方法的正确性，同时也为后续样机试验奠定了基础。

4 样机实验

乘马康复机器人主要由混联机构本体、伺服系统、传感系统、计算机控制系统组成，实验样机如图17所示。其中，传感系统主要包括马鞍两侧压力传感器和颈部角度传感器，压力传感器主要检测骑行康复训练时人腿部对机器马的夹持力，角度传感器检测骑行时对缰绳拉动的幅度。机器马运动时会根据夹持力和角度的相关信息及时调整运动步态和运动速度。

采用工控机+运动控制卡作为机器人控制器，在工控机中将上述位姿路径规划、轨迹规划和运动学算法代码化，将计算后得到的伺服电机转角数据调用API，通过PCI总线压入运动控制卡缓存Buffer，运动控制卡通过Ethercat总线将缓存数据下发到伺服系统，从而控制设备运动；同时，伺服系统的编码器也通过Ethercat总线将伺服电机的当前位置反馈到运动控制卡，可通过正解计算实时读取机器人当前位姿信息。机器人控制系统架构如图18所示。

为了展示机器人的空间运动和姿态变化，采用上述运动仿真中目标位姿作为机器人运动实验的实际位姿，运动规划的速度、加速度和插补时间间隔均与运动仿真相同，编写运动控制算法及软件的C++代码。机器人运动过程中，实时读取伺服电机编码器值，进而换算为驱动支链长度，并将驱动支链长度数据代入到机器人运动正解算法中，计算机器人运动过程中的实际位姿。将笛卡尔空间位姿规划理论值与机器人运动过程中计算实际值进行比较，如图19所示；关节空间驱动支链长度规划理论值与运动过程中读取的实际值对比结果，如图20所示。因驱动支链4的数据与其他驱动支链初始值差距较大，为便于在同幅图中对比观察，图19中舍弃了驱动支链4的数据。

由图19可知：在该控制系统内，伺服系统采用位置闭环模式，因此在额定负载以内的受力情况下机器人读取的实际笛卡尔空间位姿与理论曲线基本保持一致，但在时间周期上存在4 ms以内的偏移，其原因是：机器人理论位姿需要经过运动逆解后再下发到运动控制器的Buffer内，Buffer内的数据下发到伺服系统存在1 ms以内的Ethercat总线响应时间；机器人运动时，读取的伺服系统编码器位置也需要经过牛顿迭代的运动学正解后才可得到读取的机器人实际位姿，同样存在时延。由图20可知，机器人读取的实际关节空间支链长度与理论曲线也存在1 ms以内的时间周期偏移，该偏移主要是由Ethercat总线的响应时间造成的。由以上实验结果可知，本文设计的混联乘马康复机器人能够实现既定的位姿规划运动，控制系统运动响应迅速，具备乘马康复运动训练的基础条件。

5 结论

本文以实现生物马运动时的可调多种步态为任务需求，设计一种3‒RPS‒P型混联机器人，以满足马术治疗运动训练的康复需求。主要工作如下。

1）所设计的一种3‒RPS‒P型混联机器人，其末端自由度较高，灵活性好，并联部分机构刚度较大，承载能力强，可以满足康复治疗时马鞍中心的运动需求和乘马时的承载能力需求。

2）建立混联机构的运动学模型，推导出机构运动学逆解的公式，采用牛顿迭代数值解的方法得到机构的运动学正解；研究机器人操作空间和关节空间的速度雅可比映射关系；分析在运动副限制下的机器人工作空间，得到机构工作空间图像。

3）利用算例对机构的运动学理论进行分析，对比结果验证了理论分析的正确性和有效性；针对马匹慢步运动过程中马鞍中心部分运动轨迹进行运动仿真研究，绘制出操作空间规划的位移、速度和加速度曲线和关节空间支链长度的曲线，得到机构运动仿真的过程图、运动路径和位移速度曲线。从仿真结果可以看出，仿真过程连续轨迹平滑，可以实现设定的位姿规划运动，并进一步验证了机构自由度、运动性能分析结果及轨迹规划算法的正确性。

4）完成实验样机的搭建，针对运动仿真中的位姿规划开展实验。实验结果与仿真一致，实现了预期的运动轨迹，且控制系统运动响应较快。运动学正解由于计算过于复杂，耗时3 ms以内，但作为显示与监控需求使用符合要求。

后续将研究机器人实现马匹快步、跑步、跳跃等其他步态的运动，分析不同运动步态、不同幅度、不同速度对康复治疗效果的影响。

参考文献

[1] Li Jie.Research on motion monitoring and analysis for hippotherapy based on body sensor data fusion[D].Dalian:Dalian University of Technology,2021.

[2] 李杰.融合体感信息的马术治疗运动监测与分析研究[D].大连:大连理工大学,2021.

[3] Feng Yuxia, Pang Wei, Li Xin,et al.The Prevalence of cerebral palsy in children aged 0-6 years in China:A meta-analysis[J].Chinese General Practice,2021,24(5):603‒607.

[4] 封玉霞,庞伟,李鑫,等.中国 0～6岁儿童脑瘫患病率的 Meta 分析[J].中国全科医学,2021,24(5):603-607.

[5] Michael-Asalu A, Taylor G, Campbell H,et al.Cerebral palsy diagnosis,epidemiology,genetics,and clinical update[J].Advances in Pediatrics,2019,66:189-208. doi:10.1016/j.yapd.2019.04.002

[6] Orekhov G, Fang Ying, Luque J,et al.Ankle exoskeleton assistance can improve over-ground walking economy in individuals with cerebral palsy[J].IEEE Transactions on Neu-ral Systems and Rehabilitation Engineering,2020,28(2):461-467. doi:10.1109/tnsre.2020.2965029

[7] Panteliadis C.Cerebral palsy:A multidisciplinary approach neurocutaneous disorders pediatric neurology[M].Berlin:Springer,2018.

[8] Misgeld B J E, Lüken M, Heitzmann D,et al.Body-sensor-network-based spasticity detection[J].IEEE Journal of Biomedical and Health Informatics,2016,20(3):748-755. doi:10.1109/jbhi.2015.2477245

[9] Silkwood-Sherer D J, Killian C B, Long T M,et al.Hippotherapy—An intervention to habilitate balance deficits in children with movement disorders:A clinical trial[J].Physical Therapy,2012,92(5):707-717. doi:10.2522/ptj.20110081

[10] Park E S, Rha D W, Shin J S,et al.Effects of hippotherapy on gross motor function and functional performance of children with cerebral palsy[J].Yonsei Medical Journal,2014,55(6):1736-1742. doi:10.3349/ymj.2014.55.6.1736

[11] Kwon J Y, Chang H J, Yi S H,et al.Effect of hippotherapy on gross motor function in children with cerebral palsy:A randomized controlled trial[J].Journal of Alternative and Complementary Medicine,2015,21(1):15-21. doi:10.1089/acm.2014.0021

[12] Eskola R, Handroos H.Novel horseback riding simulator based on 6-DOF motion measurement,a motion base,and interactive control of gaits[J].Advanced Robotics,2013,27(16):1249-1257. doi:10.1080/01691864.2013.824134

[13] Lee J N, Lee M W, Byeon Y H,et al.Classification of horse gaits using FCM-based neuro-fuzzy classifier from the transformed data information of inertial sensor[J].Sensors,2016,16(5):664. doi:10.3390/s16050664

[14] Yu C H, Hong C U, Kang S R,et al.Analysis of basal physical fitness and lumbar muscle function according to indoor horse riding exercise[J].Bio-Medical Materials and Engineering,2014,24(6):2395-2405. doi:10.3233/bme-141053

[15] Zhang Bo.Design and research on control system of fitness horse[D].Shenyang:Northeastern University,2009.

[16] 张博.健身马控制系统设计与研究[D].沈阳:东北大学,2009.

[17] Feng Guanglei.Mechatronics design of gymnastic horse[D].Shenyang:Northeastern University,2011.

[18] 冯光磊.健身马的机电一体化设计[D].沈阳:东北大学,2011.

[19] Yu Shunnian, Zhang Qiuhua, Jiang Tao,et al.Kinematics analysis and simulation on a mechanical horse-riding machine[J].Journal of Jimei University(Natural Science),2017,22(1):41-45. doi:10.3969/j.issn.1007-7405.2017.01.006

[20] 余顺年,张秋华,姜涛,等.机械式骑马机运动学分析与仿真[J].集美大学学报(自然科学版),2017,22(1):41-45. doi:10.3969/j.issn.1007-7405.2017.01.006

[21] Li Jie.Research on motion monitoring and analysis for hippotherapy based on body sensor data fusion[D].Dalian:Dalian University of Technology,2021.

[22] 李杰.融合体感信息的马术治疗运动监测与分析研究[D].大连:大连理工大学,2021.

[23] Engell M T, Hernlund E, Byström A,et al.Head,trunk and pelvic kinematics in the frontal plane in un-mounted horseback riders rocking a balance chair from side-to-side[J].Comparative Exercise Physiology,2018,14(4):249-260. doi:10.3920/CEP170036

[24] Bourouis S, Alroobaea R, Rubaiee S,et al.Toward effective medical image analysis using hybrid approaches—Review,challenges and applications[J].Information,2020,11(3):155. doi:10.3390/info11030155

[25] Sánchez M B, Ian L, Holmes P,et al.Quantification of head motion in children with cerebral palsy when testing segmental trunk control[C]//Proceedings of the 2018 IX International Seminar of Biomedical Engineering(SIB).Bogota:IEEE,2018:1-4. doi:10.1109/sib.2018.8467741

[26] Coleman Wood K A, Lowndes B R, Buus R J,et al.Evidence-based intraoperative microbreak activities for reducing musculoskeletal injuries in the operating room[J].Work,2018,60(4):649-659. doi:10.3233/wor-182772

[27] Zhai Xiaoxue, Zhang Hao.Application of mechanical horseback riding in rehabilitation(review)[J].Chinese Journal of Rehabilitation Theory and Practice,2016,22(12):1412-1415. doi:10.3969/j.issn.1006-9771.2016.12.010

[28] 翟晓雪,张皓.骑马机在康复治疗领域的研究进展[J].中国康复理论与实践,2016,22(12):1412-1415. doi:10.3969/j.issn.1006-9771.2016.12.010

[29] Ye Zenglin, Zhang Liang'an, Chen Hua,et al.Design and motion performance analysis of 2URR-SRR-RUPUR parallel leg rehabilitation robot mechanism[J].Journal of Mechanical Transmission,2021,45(3):82-89.

[30] 叶增林,张良安,陈华,等.2URR-SRR-RUPUR并联式腿部康复机器人机构设计与运动性能分析[J].机械传动,2021,45(3):82-89.

[31] Ye Zenglin, Zhang Liang'an, Wu Shoupeng, et al. Motion performance analysis and control of robotic hydraulic manipulator arm[J]. Advanced Engineering Sciences,2025,57(2):171-179.

[32] 叶增林,张良安,吴守鹏,等.机器人液压机械臂运动性能分析与控制[J].工程科学与技术, 2025,57(2):171-179.

基金资助

安徽高校自然科学研究重点项目(2024AH051780)

安徽高等学校省级质量工程项目(2024fwxx071)

安徽优秀青年教师培育项目(YQYB2024100)

产学合作协同育人项目(240900516123359)

马鞍山学院校级质量工程项目(202474)