剖析人体踝关节生理结构及运动特性，提出一种四自由度广义球面并联机构来满足踝关节康复运动需求。采用螺旋理论分析机构运动特性，基于球坐标描述位置、欧拉角描述姿态的方法，建立机构位置逆解封闭矢量方程。根据逆解约束条件分析机构工作空间，以机构动球心工作空间匹配踝部胫距关节运动空间、动平台工作空间匹配踝部距下关节运动空间为目标，采用遗传算法对机构关键尺度参数进行优化设计。研究结果表明，该广义球面并联机构工作空间能够满足人体踝关节康复实际运动空间。通过本文研究结果，可为高匹配度踝关节康复机器人本体设计提供理论依据和实验参考。
引言
踝关节是人体下肢最复杂的关节，在人体运动机能中起到保持身体平衡的作用，且很容易受到高强度的撞击，肌肉控制力不足或受力不均匀时踝关节容易损伤。相关数据显示，在人体下肢损伤中，踝关节扭伤发生率居于第二位，如不及时或彻底治疗，则易产生二次损伤，影响人体正常行走[1-3]。传统的踝关节康复方式为医师对患者进行“一对一”的徒手疗法，该过程耗时耗力且无科学的测评康复体系。针对“患者多、医师少”的问题，同时提供科学的康复评价标准，设计一种适用于踝部关节康复的机器人显得至关重要。
目前，国内外专家学者基于对人体踝部关节的关节结构组成、运动形式以及损伤机制深入研究，提出多种不同特点和应用的踝关节康复机器人。Girone 等[4]基于斯图尔特（Stewart）平台提出了一种六自由度踝关节康复机器人，其机构复杂、控制困难。Saglia 等[5]研制出高性能二自由度踝关节康复机器人，能完成高灵巧性的康复运动，但康复运动自由度需求不足。Malosio 等[6]基于灵巧眼球形机构设计了一种踝关节康复机器人，其在康复训练中踝关节运动中心与机构转动中心重合度高；边辉等[7]为使机器人机构的转动中心与踝关节中心相重合，减少人机交互力，提出了一种具有远程转动中心的踝关节康复机器人；李剑锋等[8-10]提出了一系列并联机构作为踝关节康复机器人机构构型；Wang 等[11]提出一种用于踝关节康复的冗余驱动并联机器人；以上几种康复机器人的核心问题都在于踝关节运动拟合模型不精确。张彦斌等[12]基于无耦合转动并联机构设计了一种踝关节康复机器人，该机构仅含有两条分支，结构简单、运动学解耦性好，但未考虑人机匹配度。
基于上述研究与分析，现有踝关节康复机器人在不同程度上存在诸多问题。从踝关节解剖学结构特点研究，基于其损伤机制及运动方式，现阶段踝关节康复机器人研究的核心内容为踝关节与康复机器人的拟合精度。本文为提高人机康复的拟合精度，在分析踝关节解剖学结构及现有踝关节康复机器人特性的基础上，提出一种适用于踝关节康复运动的四自由度广义球面并联机构，分析其自由度特性，通过求解机构的位置逆解得到其工作空间，基于工作空间与踝关节运动空间的匹配度进行机构参数优化，为进一步设计踝关节康复机器人本体提供理论依据。
1 踝关节结构分析及新型机构设计
1.1 踝关节结构分析
为提高人机拟合精度，从本质上剖析踝关节生理结构及其运动方式是研究设计踝关节康复机器人的前提和关键。从人体解剖学角度分析，踝关节是一种复合关节，由两个独立的关节组成：胫距关节和距下关节[13]。这两种关节各自有其旋转轴线[14]，且旋转轴线是倾斜的，空间异面的，如图 1 所示。
图1
踝关节结构
Figure1.
Ankle structure
现有研究为简化设计，在踝关节康复机器人的研究过程中将踝关节运动拟合模型等效为两种形式：① 依据踝部生理结构将其运动等效为串联式两转动结构（RR 型运动拟合模型）；② 依据踝部运动形式将其运动转化为并联式三转动结构（S 型运动拟合模型）；基于文献[15]总结的踝关节康复机器人的最优设计条件是应最大程度地拟合踝关节运动，本研究结合踝关节生理结构及其运动特性采用一种串联式四转动结构（UU 型运动拟合模型）来拟合踝关节运动，如图 2 所示。
图2
踝关节运动拟合模型
Figure2.
Motion matching model of ankle joint
由于人体距骨尺寸微小的限制，以串联模型作为康复机构本体则无法布置电机，难以达到康复需求，由此构型出一种广义球面并联机构来提高踝关节康复运动的拟合精度。
1.2 广义球面并联机构
广义球面并联机构由定平台、动平台和连接动定平台的四条分支组成，构成该机构四条分支的杆件均为球面连杆，可分为 A 杆和 B 杆两类，其区别在于 A 杆含有一个运动球心，B 杆含有两个运动球心。四条分支结构为：分支一与分支二结构相同，杆件结构表示为 AABA，分支三结构表示为 ABAA，分支四结构表示为 ABA，各分支中杆件均由转动副串联连接，如图 3 广义球面并联机构分支所示。
图3
广义球面并联机构分支
Figure3.
Branch of the generalized spherical parallel mechanism
四条分支组合而成的广义球面并联机构简图如图 4 左图所示，其中 Rij（i = 1, 2, 3, 4; j = 1, 2, 3, 4, 5）表示机构第 i 条分支的第 j 个转动副中心点所在位置。各分支球面连杆排列组合后机构形成两个运动球心，AA 杆重合球心 O 称为机构的定球心，AB 杆重合球心 O1 称为机构的动球心，两者间距 OO1 称为机构的球心距。机构分支三与四结构中含有球面五杆闭环，可控制机构动球心 O1 的位置变化。所设计的康复机器人三维模型如图 4 右图所示，其康复工作原理为以机构动平台固定足部、定平台固定小腿，以定球心 O 拟合胫距关节运动中心，以动球心 O1 拟合距下关节运动中心进行的多自由度康复系统。驱动电机均布于定平台上 R11、R21、R31、R41 四个转动副处。
图4
踝关节康复机器人
Figure4.
Ankle rehabilitation robot
2 自由度分析
自由度分析是研究一切并联机构的基础，采用螺旋理论[16]对广义球面并联机构进行自由度分析。传统球面并联机构一般具有一个固定转动中心[17]，本文依据广义球面并联机构的双球心特性建立其分支坐标系。如图 5 所示，在机构一般位形下建立分支坐标系 O-XYZ，分支坐标系的 Z 轴与其动、定球心连线重合，坐标原点与定球心重合。图 5 以分支一坐标系和分支四坐标系为例。
图5
分支坐标系
Figure5.
Branching coordinate system
以 
表示第 i 条分支运动链中第 j 个运动副的运动螺旋；以 
表示第 i 条分支运动链施加于动平台上的第 j 个约束螺旋；以 
表示机构约束螺旋系中的第 i 个约束螺旋；以 
表示机构运动螺旋系中的第 i 个运动螺旋。机构各分支的运动及约束螺旋系介绍如下：
分支一的运动及约束螺旋系如式（1）～（2）所示：
 |
 |
由其约束螺旋系可知，分支一对动平台仅有 1 个力约束，且沿 Z 轴方向，无力偶约束。同理，分支二与分支一结构相同，都对动平台产生 1 个沿 Z 轴方向的力约束。
分支三的运动及约束螺旋系如式（3）～（4）所示：
 |
 |
由其约束螺旋系可知，分支三对动平台仅有 1 个力约束，且沿 Z 轴方向，无力偶约束。
分支四的运动及约束螺旋系为如式（5）～（6）所示：
 |
 |
由其约束螺旋系可知，分支四对动平台有 1 个力约束，1 个力偶约束。
四条分支施加给动平台的 5 个约束螺旋构成动平台的约束螺旋系如式（7）所示：
 |
对约束螺旋系求反螺旋可得如式（8）所示的机构运动螺旋系：
 |
式（1）～（8）中，l 表示双心距 OO1 的长度，aij、bij、cij、a4r、b4r、d4r、e4r、f4r、dip、eip（i = 1, 2, 3, 4; j = 1, 2, 3, 4, 5）代表与各转动副方向和位置相关的常量，其中上标 r 表示分支约束螺旋，上标 p 表示机构约束螺旋。
机构全部 4 个分支对动平台共施加 5 个约束，其中沿 Z 轴方向的 4 个约束力线性相关，构成一个公共约束。该机构在任一非奇异位置处，式（1）～式（8）总保持一般性，则公共约束保持不变，其自由度是全周的。由式（8）可知该广义球面并联机构具有四个自由度，与本文所提踝关节 UU 型运动拟合模型的运动特性一致。因此，可以为实现踝关节康复运动提供必要的自由度。
3 位置逆解
3.1 参数设置
设 θ1、θ2、θ3、θ4 分别表示四个转动副 R11、R21、R31、R41 的输入角，该机构设置为部分对称结构，如图 4 中广义球面并联机构简图所示，以 δi（i = 0, 1, 2, ···, 10）表示机构球面连杆的角度参数。设其球面连杆角度参数为：∠R11OR21 = δ0、∠R11OR12 = ∠R12OR22 = δ1、∠R12OR13 = ∠R22OR23 = δ2、∠OO1R14 = ∠OO1R24 = δ3、∠R14O1R15 = ∠R24O1R25 = δ4、∠R31OR41 = δ5、∠R31OR32 = ∠R41OR42 = δ6、∠R13OO1 = ∠R23OO1 = ∠R32OO1 = ∠R42OO1 = δ7、∠R33O1R34 = δ8、∠R34O1R35 = δ9、∠R43O1R44 = δ10，动定平台大小由各分支两端杆件转动副轴线长度决定，双心距 OO1 为
。
3.2 基于球坐标和欧拉角的机构位姿描述方法
在广义球面并联机构定球心处建立定坐标系 O-X0Y0Z0，在动球心处建立位置坐标系 O1-X1Y1Z1 和动坐标系 O1-X2Y2Z2。采用球坐标描述动平台位置，欧拉角 Z-Y-X 描述动平台姿态[18]，即从定系到位置系描述机构位置变化，从位置系到动系描述机构姿态变化。定系 X0 轴与转动副 R11 转动轴线重合，Y0 轴与 R31R41 角平分线重合，Z0 轴垂直定平台。设定平台与动平台同轴时为机构的初始位置，初始位置时坐标系 O1-X1Y1Z1 与 O1-X2Y2Z2 重合，机构在一般位形下的位姿描述坐标系如图 6 所示。
图6
位姿描述坐标系
Figure6.
Position and posture description coordinate system
图 6 中 O 点为定球心中心点，O1 点为动球心中心点，P 点为动平台转动中心点。α、β、γ 为欧拉角旋转角度，φ1、φ2 为球坐标旋转角度。
当已知旋转角度，位置系 O1-X1Y1Z1 相对于定系 O-X0Y0Z0 的旋转变换为绕 Z0 轴旋转 φ1 后，再绕 Y0 轴旋转 φ2，最后沿 Z0 方向移动距离 l，以 R 代表旋转变换，T 代表移动变换，则其位姿矩阵 
可表示如式（9）～式（10）所示：
 |
 |
动系 O1-X2Y2Z2 相对于位置系 O1-X1Y1Z1 的旋转变换为绕 Z1 轴旋转 α 后再绕 Y1 轴旋转 β，最后绕 X1 轴旋转 γ，其位姿矩阵 
可表示如式（11）～式（12）所示：
 |
 |
机构整体位姿描述20T 可表示如式（13）所示：
 |
该机构具有 2 个转动自由度，而描述动平台姿态有 3 个变量，需先找到描述姿态的三个欧拉角之间的关系。在给定机构某个位置下，寻求姿态角关系，机构可简化为如图 7 所示的定位置等效模型。
图7
定位置等效模型
Figure7.
Fixed position equivalent model
图 7 中分支四球面杆 R43R44 约束一个转动姿态，机构存在如式（14）所示的约束关系：
 |
依据式（14），可求得 γ = f (α, β) 关系。
3.3 逆解方程
本文采用封闭矢量法，根据机构几何关系推导位置逆解。分支一与分支二结构相同，以分支一为例，R11O + OO1 + O1R15 + R15R14 + R14R13 + R13R12 + R12R11 = 0 为矢量方程，存在约束方程如式（15）～式（17）所示：
 |
 |
 |
分支三结构封闭矢量方程为 R31O + OO1 + O1R33 + R33R32 + R32R31 = 0，存在约束方程如式（18）所示：
 |
求输入转角 θ，如式（19）～式（20）所示：
 |
 |
式（14）～（20）中 xij、yij、zij（i = 1, 2, 3, 4; j = 1, 2, 3, 4, 5）分别表示第 i 条分支中第 j 个转动副中心在定坐标系下的坐标值，xn、yn、zn（n = 1, 2）表示平面 OO1Rn4Rn3 的法向量坐标。
4 工作空间分析
采用极限边界搜索法对该并联机构的工作空间进行分析，分别求解动球心及动平台工作空间。给定机构经验参数 δ0 = 60°、δ1 = 30°、δ2 = 28°、δ3 = 110°、δ4 = 20°、δ5 = 180°、δ6 = 34°、δ7 = 74°、δ8 = 22°、δ9 = 25°、δ10 = 32°、l = 30，利用数学计算软件 MATLAB R2018b（MathWorks Inc，美国）编程求解工作空间，工作空间求解流程图如图 8 所示。
图8
工作空间求解流程图
Figure8.
Flow chart of workspace solution
图 8 中，φ1、φ2、α、β 为机构位姿描述变量。取 α 和 β 为动平台姿态角常量，φ1 和 φ2 为位置变量循环计算初值，求解机构动球心工作空间。其工作空间在半径为双心距 l 的球面上，形状变化与布局参数 δ5 相关，范围大小与 δ6、δ7 相关。列举不同构型下的动球心工作空间“点集”如图 9 所示。
图9
动球心工作空间“点集”
Figure9.
Set of moving spherical center workspace points
取图 8 中 φ1 和 φ2 为动平台转动中心的位置常量，α 和 β 为姿态角变量循环计算初值，求解机构动平台工作空间。给定初始位置 φ1 = 0 和 φ2 = 0，动平台工作空间在半径为 O1P 的球面上，其范围大小与参数 δ1、δ2、δ3、δ4、δ10 相关。分析 δ0 = 60°、δ5 = 180° 构型下机构的动平台工作空间“点集”如图 10 所示。
图10
动平台工作空间“点集”
Figure10.
Set of dynamic platform center workspace points
5 尺度优化
5.1 设计变量
以上述分析机构工作空间布局参数为例，即取 δ0 = 60°、δ5 = 180° 时，优化设计变量如式（21）所示：
 |
5.2 目标函数
康复机器人工作空间与踝关节实际运动空间的匹配度越高，造成二次损伤的可能性越小。文献[19]将踝关节基本运动简化为背伸/跖屈、内翻/外翻和内收/外展。本研究通过人体建模仿真软件 Anybody 5.1（AnyBody Technology，丹麦）驱动人体踝关节运动，采集踝关节运动时胫距关节及距下关节的踝关节运动范围如表 1 所示。
表1
踝关节运动范围
Table1.
Motion range of the ankle
以人体建模仿真软件 Anybody 5.1（AnyBody Technology，丹麦）拟合出的胫距关节及距下关节运动空间为目标，分别匹配机构动球心工作空间和动平台工作空间。以两者工作空间散点坐标差值的平方和最小为目标函数进行设计优化。目标函数如式（22）所示：
 |
Xi、Yi、Zi 为机构工作空间散点坐标值，xei、yei、zei 为目标工作空间散点坐标值，k 为散点数。
5.3 约束条件
① 仿生条件：考虑到个体尺寸差异，根据人类工效学列出的中国成年人人体尺寸的基础数值[20]，依据人体足部尺寸设定机构基本参数，以 18～60 岁男性为主要参考，足长为 272 mm，足宽为 107 mm，为预留鞋垫尺寸，设定定平台半径为 130 mm，动平台半径为 70 mm，胫距关节和距下关节轴线间距[21]为 30 mm，取 l = 30 mm。
② 杆长存在限定角度条件：两回转副中心点距离应大于连杆宽度。
③ 机构中位置球面五杆闭环的组装运动条件：δ6 + δ7 > 90°、δ6 − δ7 < 90°。
依据上述机构限定条件，结合机构实际设计装配情况，防止构件干涉可能性的增大，预设各参数取值范围为 δ1、δ2、δ6∈［15°，50°］、δ3∈［90°，180°］、δ7∈［45°，90°］、δ4、δ10∈［9°，50°］。
5.4 优化结果
通过数学计算软件 MATLAB R2018b（MathWorks Inc，美国）中的遗传算法工具箱对机构尺度参数进行求解，优化参数设置分别为：种群规模 M = 100，交叉概率 Pc = 0.8，变异概率 Pm = 0.1，终止进化代数 T = 100。遗传算法优化流程如图 11 所示。
图11
遗传算法优化流程图
Figure11.
Flow chart of genetic algorithm optimization
利用遗传算法对目标函数进行优化计算后，得到优化后机构各杆件角度参数为 δ1 = 32°、δ2 = 28°、δ3 = 135°、δ4 = 19°、δ6 = 30°、δ7 = 75°、δ10 = 35°。采用优化后的并联机构尺度参数，在数学计算软件 MATLAB R2018b（MathWorks Inc，美国）中编程绘制并联机构的工作空间点云图，与人体建模仿真软件 Anybody 5.1（AnyBody Technology，丹麦）拟合出踝部实际运动空间对比，即胫距-动球心空间对比和距下-动平台空间对比，分析优化效果，如图 12 所示。
图12
工作空间对比图
Figure12.
Workspace comparison diagram
在图 12 中，点集表示优化后的机构工作空间，光滑球面表示踝部关节目标运动空间，约束驱动角，对比优化后机构动球心工作空间与胫距关节实际运动空间及机构动平台工作空间与距下关节实际运动空间。结果表明，优化后尺度参数在满足机构约束条件下，生成的机构工作空间包络踝关节目标运动空间，即机构的工作空间能够满足人体踝关节康复所需实际运动空间的需求。
6 结论
（1）本文基于人体踝关节生理结构及其运动特性，综合现有踝关节康复机器人，为提高人机拟合精度，提出了一种高匹配度的四自由度广义球面并联机构。
（2）采用螺旋理论分析出机构的自由度特性；建立了广义球面并联机构的运动学位姿描述方法，采用封闭矢量法推导出位置逆解方程；根据逆解约束条件及驱动范围，通过 MATLAB R2018b（MathWorks Inc，美国）软件编程分析得到机构在不同构型下的动球心和动平台工作空间。
（3）为满足康复运动空间需求，以人体建模仿真软件 Anybody 5.1（AnyBody Technology，丹麦）数据拟合空间与并联机构工作空间的匹配度最高为目标，通过遗传算法对机构关键尺度参数进行了优化，得到符合人体实际运动空间的机构尺度参数，对比分析表明该机构的设计能较好保障踝关节康复所需实际运动空间。
（4）针对本文所提的踝关节康复机器人，为更好地匹配踝关节康复运动，实现人机拟合精度的提高，还需进一步实验采集踝部关节运动的实际数据，主要采集踝部胫距关节及距下关节的关联运动数据，通过数据拟合建立其数学运动模型，将为踝关节康复运动的轨迹规划提供理论依据。
利益冲突声明：本文全体作者均声明不存在利益冲突。
引言
踝关节是人体下肢最复杂的关节，在人体运动机能中起到保持身体平衡的作用，且很容易受到高强度的撞击，肌肉控制力不足或受力不均匀时踝关节容易损伤。相关数据显示，在人体下肢损伤中，踝关节扭伤发生率居于第二位，如不及时或彻底治疗，则易产生二次损伤，影响人体正常行走[1-3]。传统的踝关节康复方式为医师对患者进行“一对一”的徒手疗法，该过程耗时耗力且无科学的测评康复体系。针对“患者多、医师少”的问题，同时提供科学的康复评价标准，设计一种适用于踝部关节康复的机器人显得至关重要。
目前，国内外专家学者基于对人体踝部关节的关节结构组成、运动形式以及损伤机制深入研究，提出多种不同特点和应用的踝关节康复机器人。Girone 等[4]基于斯图尔特（Stewart）平台提出了一种六自由度踝关节康复机器人，其机构复杂、控制困难。Saglia 等[5]研制出高性能二自由度踝关节康复机器人，能完成高灵巧性的康复运动，但康复运动自由度需求不足。Malosio 等[6]基于灵巧眼球形机构设计了一种踝关节康复机器人，其在康复训练中踝关节运动中心与机构转动中心重合度高；边辉等[7]为使机器人机构的转动中心与踝关节中心相重合，减少人机交互力，提出了一种具有远程转动中心的踝关节康复机器人；李剑锋等[8-10]提出了一系列并联机构作为踝关节康复机器人机构构型；Wang 等[11]提出一种用于踝关节康复的冗余驱动并联机器人；以上几种康复机器人的核心问题都在于踝关节运动拟合模型不精确。张彦斌等[12]基于无耦合转动并联机构设计了一种踝关节康复机器人，该机构仅含有两条分支，结构简单、运动学解耦性好，但未考虑人机匹配度。
基于上述研究与分析，现有踝关节康复机器人在不同程度上存在诸多问题。从踝关节解剖学结构特点研究，基于其损伤机制及运动方式，现阶段踝关节康复机器人研究的核心内容为踝关节与康复机器人的拟合精度。本文为提高人机康复的拟合精度，在分析踝关节解剖学结构及现有踝关节康复机器人特性的基础上，提出一种适用于踝关节康复运动的四自由度广义球面并联机构，分析其自由度特性，通过求解机构的位置逆解得到其工作空间，基于工作空间与踝关节运动空间的匹配度进行机构参数优化，为进一步设计踝关节康复机器人本体提供理论依据。
1 踝关节结构分析及新型机构设计
1.1 踝关节结构分析
为提高人机拟合精度，从本质上剖析踝关节生理结构及其运动方式是研究设计踝关节康复机器人的前提和关键。从人体解剖学角度分析，踝关节是一种复合关节，由两个独立的关节组成：胫距关节和距下关节[13]。这两种关节各自有其旋转轴线[14]，且旋转轴线是倾斜的，空间异面的，如图 1 所示。
图1
踝关节结构
Figure1.
Ankle structure
现有研究为简化设计，在踝关节康复机器人的研究过程中将踝关节运动拟合模型等效为两种形式：① 依据踝部生理结构将其运动等效为串联式两转动结构（RR 型运动拟合模型）；② 依据踝部运动形式将其运动转化为并联式三转动结构（S 型运动拟合模型）；基于文献[15]总结的踝关节康复机器人的最优设计条件是应最大程度地拟合踝关节运动，本研究结合踝关节生理结构及其运动特性采用一种串联式四转动结构（UU 型运动拟合模型）来拟合踝关节运动，如图 2 所示。
图2
踝关节运动拟合模型
Figure2.
Motion matching model of ankle joint
由于人体距骨尺寸微小的限制，以串联模型作为康复机构本体则无法布置电机，难以达到康复需求，由此构型出一种广义球面并联机构来提高踝关节康复运动的拟合精度。
1.2 广义球面并联机构
广义球面并联机构由定平台、动平台和连接动定平台的四条分支组成，构成该机构四条分支的杆件均为球面连杆，可分为 A 杆和 B 杆两类，其区别在于 A 杆含有一个运动球心，B 杆含有两个运动球心。四条分支结构为：分支一与分支二结构相同，杆件结构表示为 AABA，分支三结构表示为 ABAA，分支四结构表示为 ABA，各分支中杆件均由转动副串联连接，如图 3 广义球面并联机构分支所示。
图3
广义球面并联机构分支
Figure3.
Branch of the generalized spherical parallel mechanism
四条分支组合而成的广义球面并联机构简图如图 4 左图所示，其中 Rij（i = 1, 2, 3, 4; j = 1, 2, 3, 4, 5）表示机构第 i 条分支的第 j 个转动副中心点所在位置。各分支球面连杆排列组合后机构形成两个运动球心，AA 杆重合球心 O 称为机构的定球心，AB 杆重合球心 O1 称为机构的动球心，两者间距 OO1 称为机构的球心距。机构分支三与四结构中含有球面五杆闭环，可控制机构动球心 O1 的位置变化。所设计的康复机器人三维模型如图 4 右图所示，其康复工作原理为以机构动平台固定足部、定平台固定小腿，以定球心 O 拟合胫距关节运动中心，以动球心 O1 拟合距下关节运动中心进行的多自由度康复系统。驱动电机均布于定平台上 R11、R21、R31、R41 四个转动副处。
图4
踝关节康复机器人
Figure4.
Ankle rehabilitation robot
2 自由度分析
自由度分析是研究一切并联机构的基础，采用螺旋理论[16]对广义球面并联机构进行自由度分析。传统球面并联机构一般具有一个固定转动中心[17]，本文依据广义球面并联机构的双球心特性建立其分支坐标系。如图 5 所示，在机构一般位形下建立分支坐标系 O-XYZ，分支坐标系的 Z 轴与其动、定球心连线重合，坐标原点与定球心重合。图 5 以分支一坐标系和分支四坐标系为例。
图5
分支坐标系
Figure5.
Branching coordinate system
以 表示第 i 条分支运动链中第 j 个运动副的运动螺旋；以 表示第 i 条分支运动链施加于动平台上的第 j 个约束螺旋；以 表示机构约束螺旋系中的第 i 个约束螺旋；以 表示机构运动螺旋系中的第 i 个运动螺旋。机构各分支的运动及约束螺旋系介绍如下：
分支一的运动及约束螺旋系如式（1）～（2）所示：
|
|
由其约束螺旋系可知，分支一对动平台仅有 1 个力约束，且沿 Z 轴方向，无力偶约束。同理，分支二与分支一结构相同，都对动平台产生 1 个沿 Z 轴方向的力约束。
分支三的运动及约束螺旋系如式（3）～（4）所示：
|
|
由其约束螺旋系可知，分支三对动平台仅有 1 个力约束，且沿 Z 轴方向，无力偶约束。
分支四的运动及约束螺旋系为如式（5）～（6）所示：
|
|
由其约束螺旋系可知，分支四对动平台有 1 个力约束，1 个力偶约束。
四条分支施加给动平台的 5 个约束螺旋构成动平台的约束螺旋系如式（7）所示：
|
对约束螺旋系求反螺旋可得如式（8）所示的机构运动螺旋系：
|
式（1）～（8）中，l 表示双心距 OO1 的长度，aij、bij、cij、a4r、b4r、d4r、e4r、f4r、dip、eip（i = 1, 2, 3, 4; j = 1, 2, 3, 4, 5）代表与各转动副方向和位置相关的常量，其中上标 r 表示分支约束螺旋，上标 p 表示机构约束螺旋。
机构全部 4 个分支对动平台共施加 5 个约束，其中沿 Z 轴方向的 4 个约束力线性相关，构成一个公共约束。该机构在任一非奇异位置处，式（1）～式（8）总保持一般性，则公共约束保持不变，其自由度是全周的。由式（8）可知该广义球面并联机构具有四个自由度，与本文所提踝关节 UU 型运动拟合模型的运动特性一致。因此，可以为实现踝关节康复运动提供必要的自由度。
3 位置逆解
3.1 参数设置
设 θ1、θ2、θ3、θ4 分别表示四个转动副 R11、R21、R31、R41 的输入角，该机构设置为部分对称结构，如图 4 中广义球面并联机构简图所示，以 δi（i = 0, 1, 2, ···, 10）表示机构球面连杆的角度参数。设其球面连杆角度参数为：∠R11OR21 = δ0、∠R11OR12 = ∠R12OR22 = δ1、∠R12OR13 = ∠R22OR23 = δ2、∠OO1R14 = ∠OO1R24 = δ3、∠R14O1R15 = ∠R24O1R25 = δ4、∠R31OR41 = δ5、∠R31OR32 = ∠R41OR42 = δ6、∠R13OO1 = ∠R23OO1 = ∠R32OO1 = ∠R42OO1 = δ7、∠R33O1R34 = δ8、∠R34O1R35 = δ9、∠R43O1R44 = δ10，动定平台大小由各分支两端杆件转动副轴线长度决定，双心距 OO1 为。
3.2 基于球坐标和欧拉角的机构位姿描述方法
在广义球面并联机构定球心处建立定坐标系 O-X0Y0Z0，在动球心处建立位置坐标系 O1-X1Y1Z1 和动坐标系 O1-X2Y2Z2。采用球坐标描述动平台位置，欧拉角 Z-Y-X 描述动平台姿态[18]，即从定系到位置系描述机构位置变化，从位置系到动系描述机构姿态变化。定系 X0 轴与转动副 R11 转动轴线重合，Y0 轴与 R31R41 角平分线重合，Z0 轴垂直定平台。设定平台与动平台同轴时为机构的初始位置，初始位置时坐标系 O1-X1Y1Z1 与 O1-X2Y2Z2 重合，机构在一般位形下的位姿描述坐标系如图 6 所示。
图6
位姿描述坐标系
Figure6.
Position and posture description coordinate system
图 6 中 O 点为定球心中心点，O1 点为动球心中心点，P 点为动平台转动中心点。α、β、γ 为欧拉角旋转角度，φ1、φ2 为球坐标旋转角度。
当已知旋转角度，位置系 O1-X1Y1Z1 相对于定系 O-X0Y0Z0 的旋转变换为绕 Z0 轴旋转 φ1 后，再绕 Y0 轴旋转 φ2，最后沿 Z0 方向移动距离 l，以 R 代表旋转变换，T 代表移动变换，则其位姿矩阵 可表示如式（9）～式（10）所示：
|
|
动系 O1-X2Y2Z2 相对于位置系 O1-X1Y1Z1 的旋转变换为绕 Z1 轴旋转 α 后再绕 Y1 轴旋转 β，最后绕 X1 轴旋转 γ，其位姿矩阵 可表示如式（11）～式（12）所示：
|
|
机构整体位姿描述20T 可表示如式（13）所示：
|
该机构具有 2 个转动自由度，而描述动平台姿态有 3 个变量，需先找到描述姿态的三个欧拉角之间的关系。在给定机构某个位置下，寻求姿态角关系，机构可简化为如图 7 所示的定位置等效模型。
图7
定位置等效模型
Figure7.
Fixed position equivalent model
图 7 中分支四球面杆 R43R44 约束一个转动姿态，机构存在如式（14）所示的约束关系：
|
依据式（14），可求得 γ = f (α, β) 关系。
3.3 逆解方程
本文采用封闭矢量法，根据机构几何关系推导位置逆解。分支一与分支二结构相同，以分支一为例，R11O + OO1 + O1R15 + R15R14 + R14R13 + R13R12 + R12R11 = 0 为矢量方程，存在约束方程如式（15）～式（17）所示：
|
|
|
分支三结构封闭矢量方程为 R31O + OO1 + O1R33 + R33R32 + R32R31 = 0，存在约束方程如式（18）所示：
|
求输入转角 θ，如式（19）～式（20）所示：
|
|
式（14）～（20）中 xij、yij、zij（i = 1, 2, 3, 4; j = 1, 2, 3, 4, 5）分别表示第 i 条分支中第 j 个转动副中心在定坐标系下的坐标值，xn、yn、zn（n = 1, 2）表示平面 OO1Rn4Rn3 的法向量坐标。
4 工作空间分析
采用极限边界搜索法对该并联机构的工作空间进行分析，分别求解动球心及动平台工作空间。给定机构经验参数 δ0 = 60°、δ1 = 30°、δ2 = 28°、δ3 = 110°、δ4 = 20°、δ5 = 180°、δ6 = 34°、δ7 = 74°、δ8 = 22°、δ9 = 25°、δ10 = 32°、l = 30，利用数学计算软件 MATLAB R2018b（MathWorks Inc，美国）编程求解工作空间，工作空间求解流程图如图 8 所示。
图8
工作空间求解流程图
Figure8.
Flow chart of workspace solution
图 8 中，φ1、φ2、α、β 为机构位姿描述变量。取 α 和 β 为动平台姿态角常量，φ1 和 φ2 为位置变量循环计算初值，求解机构动球心工作空间。其工作空间在半径为双心距 l 的球面上，形状变化与布局参数 δ5 相关，范围大小与 δ6、δ7 相关。列举不同构型下的动球心工作空间“点集”如图 9 所示。
图9
动球心工作空间“点集”
Figure9.
Set of moving spherical center workspace points
取图 8 中 φ1 和 φ2 为动平台转动中心的位置常量，α 和 β 为姿态角变量循环计算初值，求解机构动平台工作空间。给定初始位置 φ1 = 0 和 φ2 = 0，动平台工作空间在半径为 O1P 的球面上，其范围大小与参数 δ1、δ2、δ3、δ4、δ10 相关。分析 δ0 = 60°、δ5 = 180° 构型下机构的动平台工作空间“点集”如图 10 所示。
图10
动平台工作空间“点集”
Figure10.
Set of dynamic platform center workspace points
5 尺度优化
5.1 设计变量
以上述分析机构工作空间布局参数为例，即取 δ0 = 60°、δ5 = 180° 时，优化设计变量如式（21）所示：
|
5.2 目标函数
康复机器人工作空间与踝关节实际运动空间的匹配度越高，造成二次损伤的可能性越小。文献[19]将踝关节基本运动简化为背伸/跖屈、内翻/外翻和内收/外展。本研究通过人体建模仿真软件 Anybody 5.1（AnyBody Technology，丹麦）驱动人体踝关节运动，采集踝关节运动时胫距关节及距下关节的踝关节运动范围如表 1 所示。
表1
踝关节运动范围
Table1.
Motion range of the ankle
以人体建模仿真软件 Anybody 5.1（AnyBody Technology，丹麦）拟合出的胫距关节及距下关节运动空间为目标，分别匹配机构动球心工作空间和动平台工作空间。以两者工作空间散点坐标差值的平方和最小为目标函数进行设计优化。目标函数如式（22）所示：
|
Xi、Yi、Zi 为机构工作空间散点坐标值，xei、yei、zei 为目标工作空间散点坐标值，k 为散点数。
5.3 约束条件
① 仿生条件：考虑到个体尺寸差异，根据人类工效学列出的中国成年人人体尺寸的基础数值[20]，依据人体足部尺寸设定机构基本参数，以 18～60 岁男性为主要参考，足长为 272 mm，足宽为 107 mm，为预留鞋垫尺寸，设定定平台半径为 130 mm，动平台半径为 70 mm，胫距关节和距下关节轴线间距[21]为 30 mm，取 l = 30 mm。
② 杆长存在限定角度条件：两回转副中心点距离应大于连杆宽度。
③ 机构中位置球面五杆闭环的组装运动条件：δ6 + δ7 > 90°、δ6 − δ7 < 90°。
依据上述机构限定条件，结合机构实际设计装配情况，防止构件干涉可能性的增大，预设各参数取值范围为 δ1、δ2、δ6∈［15°，50°］、δ3∈［90°，180°］、δ7∈［45°，90°］、δ4、δ10∈［9°，50°］。
5.4 优化结果
通过数学计算软件 MATLAB R2018b（MathWorks Inc，美国）中的遗传算法工具箱对机构尺度参数进行求解，优化参数设置分别为：种群规模 M = 100，交叉概率 Pc = 0.8，变异概率 Pm = 0.1，终止进化代数 T = 100。遗传算法优化流程如图 11 所示。
图11
遗传算法优化流程图
Figure11.
Flow chart of genetic algorithm optimization
利用遗传算法对目标函数进行优化计算后，得到优化后机构各杆件角度参数为 δ1 = 32°、δ2 = 28°、δ3 = 135°、δ4 = 19°、δ6 = 30°、δ7 = 75°、δ10 = 35°。采用优化后的并联机构尺度参数，在数学计算软件 MATLAB R2018b（MathWorks Inc，美国）中编程绘制并联机构的工作空间点云图，与人体建模仿真软件 Anybody 5.1（AnyBody Technology，丹麦）拟合出踝部实际运动空间对比，即胫距-动球心空间对比和距下-动平台空间对比，分析优化效果，如图 12 所示。
图12
工作空间对比图
Figure12.
Workspace comparison diagram
在图 12 中，点集表示优化后的机构工作空间，光滑球面表示踝部关节目标运动空间，约束驱动角，对比优化后机构动球心工作空间与胫距关节实际运动空间及机构动平台工作空间与距下关节实际运动空间。结果表明，优化后尺度参数在满足机构约束条件下，生成的机构工作空间包络踝关节目标运动空间，即机构的工作空间能够满足人体踝关节康复所需实际运动空间的需求。
6 结论
（1）本文基于人体踝关节生理结构及其运动特性，综合现有踝关节康复机器人，为提高人机拟合精度，提出了一种高匹配度的四自由度广义球面并联机构。
（2）采用螺旋理论分析出机构的自由度特性；建立了广义球面并联机构的运动学位姿描述方法，采用封闭矢量法推导出位置逆解方程；根据逆解约束条件及驱动范围，通过 MATLAB R2018b（MathWorks Inc，美国）软件编程分析得到机构在不同构型下的动球心和动平台工作空间。
（3）为满足康复运动空间需求，以人体建模仿真软件 Anybody 5.1（AnyBody Technology，丹麦）数据拟合空间与并联机构工作空间的匹配度最高为目标，通过遗传算法对机构关键尺度参数进行了优化，得到符合人体实际运动空间的机构尺度参数，对比分析表明该机构的设计能较好保障踝关节康复所需实际运动空间。
（4）针对本文所提的踝关节康复机器人，为更好地匹配踝关节康复运动，实现人机拟合精度的提高，还需进一步实验采集踝部关节运动的实际数据，主要采集踝部胫距关节及距下关节的关联运动数据，通过数据拟合建立其数学运动模型，将为踝关节康复运动的轨迹规划提供理论依据。
利益冲突声明：本文全体作者均声明不存在利益冲突。
