基于Delaunay三角划分网格计算多连杆机器人工作空间

吴 虎， 李鑫宁， 杨先海， 薛 鹏， 谭 帅

(1. 山东理工大学机械工程学院， 山东淄博 255000； 2.国家采油装备工程技术研究中心， 山东东营 257091)

多连杆机器人是由多个活动关节将刚性杆件组成的连杆机构，具有工作空间范围大、灵活性好、机械结构简单、避障性能好、可避免机构奇异性等优点，因此广泛应用于工业生产、工程实践、医疗手术等领域[1-2]。多连杆机器人为机械装置的主要执行部件，其灵活程度对于整个机构系统的流畅运行起着至关重要的作用。

工作空间的基本概念[3]及灵活性[4]、 边界的确定[5]、灵活角的分析[6]、 工作空间的解析关系[7]等研究领域的探索， 为多连杆机器人工作空间的研究提供了理论基础。 目前， 多连杆机器人工作空间的求解方法比较多， 主要可以归纳为图解法(几何法)、 解析法和数值积分法三大类。 图解法[8]利用几何绘图的原理求得工作空间的边界， 该方法直观性强。 García-Murillo等[9]提出一种用于分析2(3-RRPS, R为转动副, P为移动副, S为球副)并联机器人工作空间的几何方法， 通过计算机辅助设计(CAD)软件进行实例分析计算。 刘怀周等[10]利用图解法分析深水作业串联机械臂的工作空间边界曲线， 运用遗传算法求得工作空间最优解。 解析法[11]通过多次求解边界曲面(线)的包络问题， 确定工作空间的边界方程。 Chen等[12]运用解析法确定混合刻板绘图仪的工作空间边界， 并采用优化算法计算了工作空间的体积。 于红英等[13]利用解析法建立了并联机构各支链工作空间的边界方程。 何世凯等[14]提出矩阵零空间的方法求解力旋量工作空间边界。 利用数值积分法[15]计算多连杆机器人的工作空间时， 通过正运动学方程计算多连杆机器人末端执行机构端点的坐标， 从而形成工作空间边界曲线或曲面。 随着计算机软、 硬件技术的不断发展， 数值积分法可以分析任意形式多连杆机器人的工作空间， 因此得到了越来越广泛的应用， 其中比较典型的是基于蒙特卡罗法的工作空间求解。 与图解法与解析法相比， 蒙特卡罗法省略了复杂的数学推导过程， 可编程性强， 容易实现可视化， 被广泛地应用在各种领域机器人的工作空间求解中， 但是传统的蒙特卡罗法计算工作空间不够精确。 He等[16]提出了一种基于增加边界点密度的蒙特卡罗法计算工作空间。 Zhao等[17]利用传统的蒙特卡罗法生成种子工作空间， 根据正态分布扩展种子工作空间， 设置精度阈值动态调整标准差， 提出了一种体素算法来求出工作空间的体积， 通过不断细化边界部分来减小结果误差。 Bader等[18]混合计算自运动流形范围估计工作空间包络和蒙特卡罗积分估计方向体积的算法， 针对高冗余机器人提出了一种容错工作空间的计算方法。 Hamida等[4]将遗传算法和蒙特卡罗法相结合， 针对机器人结构紧凑和灵巧性最大化进行优化。 樊桂菊等[19]建立了以可达工作空间性能和结构紧凑为指标的优化模型， 利用遗传算法求最优参数， 采用网格划分法求解工作空间体积， 达到较高的精度。 虽然上述方法在计算精度、 边界确定等方面取得了很大的进步， 但同时也存在直观性不强、 表达式繁琐、 对于多自由度的机器人运算量过大等问题。

为了更准确地描述多连杆机器人的工作空间，在工作空间密度函数[20]和自适应非结构化网格生成理论[21-22]的基础上，本文中提出基于Delaunay三角划分网格的多连杆机器人工作空间计算方法(简称本文方法)，通过与数值积分法、 网格划分法进行对比， 验证本文方法的准确性。

1 灵活工作空间与工作空间

可达工作空间是在参考坐标系中末端执行机构所能达到的所有位置坐标的集合，记作W；灵活工作空间是在总的可达工作空间内，末端执行机构可以以任意姿态达到的坐标点所构成的工作空间，记作W1；次工作空间是工作空间中去除灵活工作空间所剩余的空间，记作W2。三者的关系为

W=W1+W2

。

(1)

图1所示为平面3R机器人的运动模型。 该机构由长度分别为l1、l2、l3的连杆组成， 令l1>l2+l3，l2>l3， 机械臂驱动转角θi∈[-π,π](i=1, 2, 3)，以此例来说明灵活工作空间与工作空间的关系。

li(i=1,2,3)—连杆长度； θi—机械臂驱动转角； xoy—平面直角坐标系。图1 平面3R(R为转动副)机器人的运动模型

多连杆机器人的工作空间与灵活工作空间如图2所示。圆C1半径R1=l1-l2-l3，圆C2半径R2=l1-(l2-l3)=l1-l2+l3，圆C3半径R3=l1+(l2-l3)=l1+l2-l3，圆C4半径R4=l1+l2+l3。显然，工作空间W是圆C1和C4围成的圆环区域，灵活工作空间W1是圆C2和C3围成的圆环区域，其他剩余的工作空间为W2。

li(i=1,2,3)—连杆长度； Ch(h=1,2,3,4)—圆心为坐标原点o的不同半径的圆。图2 多连杆机器人的工作空间与灵活工作空间

2 工作空间密度函数

将机械机构的工作空间均匀地划分成等尺寸的单元，每个单元内末端执行机构可以到达的点数与整个工作空间内到达点数比值除以单元的面积或体积即为工作空间密度[23]。工作密度函数值越大，证明机械机构的灵活性越大， 同时可以更加直观地表达多连杆机器人的工作空间。 蒙特卡罗法属于实验数学的一个分支， 利用随机数进行统计实验， 将均值、 概率等作为待解决问题的数值解， 也称为蒙特卡罗模拟。 通过研究多连杆机器人的工作空间密度， 采用蒙特卡罗法，从概率统计的角度借助计算机技术进行仿真， 可以直观地描述工作空间的形状、 特征及性质。

2.1 蒙特卡罗法计算步骤

工作空间解算的蒙特卡罗算法原理[6]是

W={M(k)∶k∈G}

，

(2)

式中：k为广义关节变量；M(k)为正运动学参数；G为关节空间集合。

采用蒙特卡罗法求解多连杆机器人工作空间的主要思想如下：通过计算机仿真随机产生大量的关节空间点，映射至工作空间，得到一定数量的末端位置随机点，进而统计解算工作空间的边界。蒙特卡罗算法的步骤如下：

1)根据灵活工作空间的定义，对各个关节转角的样本进行划分取点，构建样本空间G；

2)依据正运动学参数M(k)解算机械臂末端执行机构所能到达的各个位置；

3)输出机械臂末端执行机构的位置，构建灵活工作空间W，样本容量越大，所得到的几何轮廓越清晰，灵活工作空间W越精确。

平面3R机器人正运动计算模型如图3所示。根据图中多连杆机器人末端执行机构位置与连杆长度及关节转角的几何关系，用Denavit-Hartenberg(D-H)法列出其正运动学方程。

li(i=1,2,3)—连杆长度； θi—机械臂驱动角； P—末端执行机构位置点； o、 A、 B—关节中心点； x0oy0、 x1oy1、 x2Ay2、 x2By2—平面直角坐标系。图3 平面3R(R为转动副)机器人正运动计算模型

假设坐标轴z0、z1、z2垂直于纸面向外，从空间中的点(x0,y0,z0)到(x1,y1,z1)的齐次旋转变换矩阵为

(3)

从点(x1,y1,z1)到(x2,y2,z2)的齐次旋转变换矩阵为

(4)

从点(x2,y2,z2)到(x3,y3,z3)的齐次旋转变换矩阵为

(5)

从点(x3,y3,z3)到(x0,y0,z0)的齐次旋转变换矩阵为

(6)

可得末端执行机构P点的位置矢量为

(7)

式中(xP,yP,zP)为末端执行机构位置点P点的坐标。

基于以上公式推导，可以得到两杆机器人的坐标如下：

(8)

(9)

式中li为第i个连杆的长度。

以上运动学方程是在直角坐标系下构建的，要得到极坐标系下对应的表达式，可利用直角坐标与极坐标的关系，即

(10)

(11)

式中：r为极径；φ为极角。

2.2 密度函数的构建

蒙特卡罗法计算多连杆机器人工作空间采用的是统计学的方法。 Sigmoid函数是统计学中常用的数学模型。 根据工作空间密度函数符合随机密度函数[20]的特性， 通过建立数学模型来描述工作空间密度。

设N连杆机械臂杆件的密度函数为fN(p)，其中p∈G，N连杆机械臂的工作空间密度函数为

f1,2,…,N(p)=(f1*f2*…*fN)(p)

,

(12)

式中*为在G上的卷积运算符。

ρ1,2,…,N(x)=(ρ1⊙ρ2⊙…⊙ρN)(x)

,

(13)

工作空间是连续的空间，其特点是随机变量可以连续地取值，所有的连续随机变量都有概率密度函数。u为多维向量，μ为均值，σ2为方差，多元正态分布密度函数的形式为

(14)

对式(14)进行严格泛化，得

(15)

式中：Σ为N×N型对称正定矩阵； |·|为矩阵的行列式。μ、Σ的求解如下：

(16)

(17)

(18)

连续的概率密度函数不局限于1，因此在工作空间的研究中，所有连续随机变量都是可测的，所有连续分布都有密度。对于N连杆机械臂，当每个连杆转角有M个状态时，连杆机械臂工作空间内的样本个数为K=MN，此样本空间的方差S1,2,…,N(K)为

(19)

式中XK为空间样本方差的矩阵。

当M足够大时，实际工作空间密度函数的方差为

(20)

式中In为n×n型单位矩阵。

根据卷积迭代及极值定理，Sigmoid函数构建的工作空间密度函数为

exp[-(u-μ)T(u-μ)Σ-1] 。

(21)

设杆长l1=1.8，l2=2.4，l3=3.1，在转角θ1、θ2、θ3不受限制的情况下，取样本点个数M=8 000，根据式(9)运算并绘制出工作空间密度图，如图4所示。

2.3 计算仿真

根据构建的密度函数，建立五自由度多连杆机器人模型，如图5所示，其中l1=0.42，l2=0.375，l3=0.42，l1、l2、l3与底座转轴不在一个平面内(l1、l2、l3与底座转轴沿转轴同方向偏移0.1、 0.2、 0.3)，转角θ1、θ2、θ3以及底座转角θ4、 执行机构转角θ5不受限制， 所取的样本点个数M=20 000。 图6所示为五自由度多连杆机器人工作空间仿真视图，4个视图全面描述了工作空间密度。

xy—坐标平面。图4 3R(R为转动副)机器人工作空间密度云图

xyz—空间直角坐标系。图5 五自由度多连杆机器人模型

(a)工作空间三维视图

(b)xy坐标平面

(c)xz坐标平面

(d)yz坐标平面图6 五自由度多连杆机器人工作空间仿真视图

六自由度多连杆机器人模型如图7所示，其中l1=0.42，l2=0.375，l3=0.42，l4=0.42， 杆l1、l2、l3、l4与底座转轴均在一个平面内，转角θ1、θ2、θ3、θ4以及底座转角θ5， 执行机构转角θ6不受限制，所取的样本点个数M=20 000。六自由度多连杆机器人工作空间仿真视图如图8所示。

图7 六自由度多连杆机器人模型

3 工作空间计算

根据连杆机器人运动学及工作空间密度研究，提出基于Delaunay三角划分网格计算多连杆机器人工作空间的方法， 原理如图9所示。 本文方法是一种非结构化网格生成方法，首先建立一个包含所有点的初始三角网，然后将剩余的点逐点插入，并且确保初始三角网成为Delaunay三角网，同时，每个Delaunay三角形的外接圆不包含其他结点；若出现外接圆内有其他结点，则需修改局部网格，直至满足每个Delaunay三角形的外接圆不包含其他结点这一条件。

(a)工作空间三维视图

(b)xy坐标平面

(c)xz坐标平面

(d)yz坐标平面图8 六自由度多连杆机器人工作空间仿真视图

图9 基于Delaunay三角划分网格的 多连杆机器人工作空间计算方法原理

1)建立初始网格。首先找到包含该密度空间点集的外接矩形，通常保证外接矩形的边长是点集最大边界长度的3倍。然后任意连接矩形的对角线，形成2个三角形，作为初始Delaunay三角网格。

2)点的逐次插入。在保证计算区域内工作空间密度函数均匀的情况下，插入点必在三角形内或三角形边上，3种插入点的情况如图10所示。插入点P′与三角形的3个顶点U、V、W相连，三角形及其外接圆形成Delaunay腔，然后连接插入点与空腔的每个顶点，每次插入一个新点之后，重复连接插入点与空腔的每个顶点操作。

U、 V、 W—三角形的3个顶点； V1—三角形的邻近顶点； P′—插入点。图10 点的逐次插入中插入点的3种情况

3)构造外接圆。首先对U、V、W这3个点构造外接圆T1，把所有结点包含进去，之后每次只引入一个结点，生成相应的圆T2，T3，……，直至所有结点都进入Delaunay三角划分网格。

4)删除最近边。将三角形中离新结点最近的一条边删除，对网格进行重生成，于是新结点与老结点形成新的三角形，如图11所示。

图11 通过删除最近边对网格进行重生成所得新三角形

6)工作空间面积的求和。工作空间面积S的求解公式为

(22)

式中：n1为内部空间子网格个数；n2为外部边界空间子网格个数。

以图3为模型，机器人连杆参数如表1所示。

表1 平面3R(R为转动副)机器人连杆长度与机械臂驱动转角

用数值积分法、 网格划分法和本文方法分别计算工作空间面积， 并与理论计算值进行比较， 如表2所示。 从表中可以看出， 本文方法的计算结果较准确， 与其他2种方法相比， 误差率分别降低了0.99%和0.50%， 提高了工作空间的计算精度。

表2 不同方法计算的工作空间面积比较

4 结论

为了更准确地描述多连杆机器人的工作空间，在对多连杆机器人进行运动学分析的基础上，本文中提出基于Delaunay三角划分网格计算连杆机器人工作空间，得出如下结论：

1)应用构建的Sigmoid工作空间密度函数，不仅可以直观地表达工作空间的工作范围，而且可以定量地得出工作空间体积的数值，避免了繁琐的表达式及过量的运算，为工程中的实际应用提供了有效的计算方法，同时为工作空间后续的深入研究提供了思路。

2)通过本文方法与数值积分法、 网格划分发法等其他计算方法的对比，验证了本文方法的准确性。

3)在工作空间的研究中，正运动学的D-H参数算法随着连杆机器人自由度的增加，需要在每个关节增加坐标系，计算量随之增大。如何减少计算量、 提高计算效率是需要进一步研究的问题。