基于Vine Copula的机器人运动精度可靠性分析*

刘春秋，周金宇，庄百亮，丁 力

（1.江苏理工学院 机械工程学院，常州 213001；2.金陵科技学院 机电工程学院，南京 211169；3.机械科学研究总院江苏分院有限公司，常州 213001)

0 引言

工业生产中，机器人受自身和使用环境影响，随着时间的推移绝对定位精度逐渐下降。尤其对于高精度机器人，精度的下降无疑会对正常作业和生产造成不利影响。机器人运动精度退化的主要表现是末端执行器的实际输出位姿在理想输出附近波动，且造成精度退化的原因多具有不确定性，如操作臂的杆件变形、关节磨损、制造误差和装配间隙等。机器人运动精度可靠性研究，就是衡量误差对运动精度的影响程度，用概率的方法分析机械手在指定公差范围内到达指定位置或实现预期运动轨迹的概率。

机器人运动精度一般分为单一点精度和轨迹精度，轨迹精度又可称为区间运动精度，它是综合评价机器人一段时间轨迹上的运动输出是否满足预期的一种方式。早在20世纪90年代，Bhatti[1]就开展了工业机器人运动相关概率建模并提出了机器人运动精度可靠性概念。随着机器人工业化普及和精度要求的提高，此问题得到重视和研究。对于机器人点精度可靠性分析，广大学者引入传统的结构可靠性分析方法，如一阶可靠度方法[2]，得到了较为理想的计算精度。对于机器人区间运动精度的可靠性分析，只是在近些年才有所发展。目前区间（时变）可靠性分析法主要有极值法、首次穿越法、联合概率方法以及包络法。对于极值法，其原理是通过确定误差函数的极值分布来估计区间可靠度。如MAHESH[3]将机器人轨迹离散为若干个点，利用最大熵原理推导机器人离散点的误差极值分布，并以串联系统可靠性模型计算机器人目标轨迹的运动可靠度。首次穿越理论是以运动误差第一次跨越误差阈值的概率来估计区间失效概率。如Zhang[4]针对平面函数生成机构提出了均值跨越率法，推导了上下跨越率公式，并将一次二阶矩（FOSM）法推广到时变可靠性分析中，此方法对于服从正态分布的小方差维数变量是准确的。对于平面轨迹机构的时变可靠性分析，陈放提出了联合概率方法，文中认为机构区间运动精度可靠度是区间内各离散点处运动误差落在误差限范围内的综合效应[5]，作者针对正态分布的误差，推导了轨迹离散点间的多元正态累积分布函数计算公式。由于此方法并未限定轨迹离散点的数量，故在计算精度上需进行预匹配。此外，为了提高时变可靠性分析的准确性，Du和Zhang[6，7]等提出了基于包络思想的时变可靠性分析方法，论文基于包络原理计算并筛选出时间区间上的最易失效点，通过建立失效点间的多变量正态分布累积分布函数，计算区间失效概率。文中考虑了平面连杆机构的结构误差和尺寸误差，在误差服从正态分布的情况下计算精度较高，但在变量非正态时产生局限性。为解决这一问题，Wang[8]在包络思想的基础上，引入Vine Copula模型来处理更多的误差分布类型，此方法只需要误差变量的边缘分布和Copula连接函数即可得到多维联合概率密度函数，且在平面连杆机构中验证了此方法的精度。以上研究主要针对平面函数生成机构，且获得了较理想的求解效率和精度。但对于空间轨迹机构如机器人的运动精度可靠性分析问题，亟需拓展相关理论和方法。本文借鉴平面连杆机构时变可靠性研究成果，利用Vine Copula函数和包络法思想展开机器人轨迹的运动可靠性分析，以机器人末端位置精度为目标求解区间可靠度。

1 机器人运动误差模型

1.1 机器人运动学分析

机器人空间位姿建模的方法主要有矩阵法、矢量法和旋量法等。本文选用较为常用的矩阵法，此方法基本原理为：对于具有多个自由度的串联机械手，在其每个连杆上建立关节坐标系，并定义变换参数θi，di，ai-1，αi-1，通过旋转和平移变换理论来获得用于描述任意相邻连杆i-1和i 之间变换关系的矩阵i-1iT。对于改进的DH参数模型，其一般表达式为

其中，cθi=cosθi，sθi=sinθi，θi为机器人各关节角，di为机器人相邻连杆间的偏距，ai-1为各连杆长度，αi-1为描述相邻关节轴相对位置的转角。

通过变换矩阵的组合，可实现基坐标空间到关节坐标空间的转化。机械手末端参考点在基坐标空间内的位姿矩阵可表示为：

式中，矩阵的前3行前3列为操作臂末端的姿态矩阵，第4列前3个元素分别为末端参考点在基坐标中的位置分量。通过以上运动学分析得到的末端参考点位姿模型，即为机器人运动学分析的基本方程。

1.2 机器人运动误差模型建立

影响机器人运动精度的各种不确定因素主要通过机器人关节空间的映射，到达机械臂末端参考点。考虑机器人关节空间的变换规律，选取机器人D-H参数作为误差输入变量。对于空间中任意相邻关节i-1和i，由矩阵微分变换理论[9]对变换矩阵Ti求微分，可得

基于矩阵变换原理，对式(2)中的0nT求偏导并保留一阶微分量，可得：

由式(3)知dTn=TnδTn，代入式(5)并化简可得机器人末端误差矩阵：

式(6)即为机器人末端位姿误差矩阵，可简化为：

式中，p和δ分别表示机器人末端位置误差和姿态误差，M1和M2分别为位置和姿态误差系数矩阵，[ΔθΔdΔαΔa]T为各关节几何参数误差向量。

2 机器人轨迹运动可靠性分析

2.1 轨迹运动可靠性模型建立

机器人轨迹运动受结构参数误差等因素影响，实际轨迹在期望轨迹附近波动，轨迹运动可靠性也称为区间可靠性，它是一种不同于点（静态）可靠性的时变可靠性问题，结合上一节内容，对于任意时刻t，机器人的位置运动误差可表示为：

式中，gA(X，t)为实际运动位置，gd(t)为期望运动位置，Xi=[θi，di，αi-1，αi-1]T，（i=1，2，…，n)为各关节DH参数向量，n为机器人关节个数。对于任意时间区间[t0，te]上的一段轨迹，其轨迹精度可靠度可表示为区间内每个点的位置误差小于误差阈值ε的概率，相应的可靠度模型可表示为：

其中：Pr｛·｝表示概率，gx(·)表示轨迹坐标在X轴上的运动误差函数。相应的失效概率可表示为：

其他坐标分量上的可靠度计算模型同理可得。

2.2 区间可靠度求解方法

在平面连杆机构中，传统包络法被用来计算区间可靠度，如文献[6，7]。此方法思路是：首先寻找目标时间区间内的包络展开点，然后建立展开点之间的多元正态累积分布函数来求解区间可靠度。此方法在输入变量服从正态分布的情况下是有效的，其它分布情况下将变得复杂甚至失效，在机器人工程实际中，影响其运动精度的误差分布类型具有不确定性，如在非正态误差变量情形，传统包络法将失效。因此本文在包络法的基础上对其进行改进，即在求得包络函数展开点之后，引入Vine Copula函数模型来建立展开点之间的联合概率密度函数。此方法较传统包络法具有以下优点：一是此方法不受误差变量分布类型的影响，对于变量非正态情形仍然适用；二是Copula函数描述的秩相关系数和尾相关系数对比多元正态累积分布函数中的积矩相关系数，是更为合理的相关性测度，可以较完整的反应变量间的相关特征，有助于提高求解精度。

3 轨迹运动可靠性分析中的Vine Copula模型

3.1 Copula函数理论

由Sklar定理[11]可知：多维随机向量X=(x1，x2，…，xn)的n维联合分布F(x1，x2，…，xn)可表示为所有一维边缘分布函数ui=Fi(xi)(i=1，…，n)构成的多元Copula函数C(u)，即：

通过上式可导出其概率密度：

式中，fi(·)表示各变量的边缘概率密度，i=1，2，…，n。

3.2 Vine Copula函数模型

二元Copula函数模型简单，便于计算，但只能处理二维变量间的相关问题。对于多变量相关情况，多元Copula函数虽然也可以使用，但是其求解难度较高，分析较复杂。所以本文引入由二元Copula函数组成的Vine Copula模型来建立多维随机变量间的联合概率分布。在Vine Copula模型中，Vine是一种描述联合概率分布的图形结构，用于缔造高维分布中的约束链接，使得该模型可以利用输入变量的边缘分布和Copula连接函数来建立多维联合分布。此函数模型将变量边缘分布和变量间相关结构分开研究，简化了多变量概率建模的过程[10]。对于任意维Copula函数，Vine Copula分解模型都可将其分解为多个二元Copula函数的集合。下面以3维D-Vine模型为例，图1为模型树形图。T1层列举了所有输入变量；T2层为Copula连接层，表示二元Copula函数将上方数组两两连接；T3层为条件连接层，表示Copula函数将上方条件数组连接。其中，T2层受随机变量间相关性的影响最大，越往下，影响越弱，所以由初始层样本选取第二层的Copula函数至关重要。

图1 3维D-vine模型

对于三维变量情况，假设随机变量向量X=(x1，x2，x3)，经图1输出的联合概率密度函数可表示为：

式中，c12(·)，c23(·)，c13/2(·)为双变量Copula密度函数。c13/2(·)中的条件随机变量计算如下：

其中c12(·)，c23(·)，为Copula分布函数。同样的，对于n维随机向量X=(x1，x2，…，xn)，其联合概率密度函数可表示为：

通过以上分析可知，双变量Copula函数是连接随机变量和建立联合分布模型的关键。下面给出二元Copula函数的选取方法，本文利用极大似然估计（MLE）法来获取目标Copula函数的相关参数，并以赤池信息准则[12]（AIC）作为最优Copula函数的选取标准。对于任意随机变量X1和X2，根据样本信息，可建立对数似然函数：

其中c(·)为Copula密度函数，α为Copula函数的参数。，(i=1，2，…，N)是随机变量X1和X2中的一对样本，N为数据样本容量。其中参数α可由最大化对数似然函数得到：

对所有候选Copula函数进行上述估计，得到各参数集，以变量样本对和各候选Copula参数集作为输入，代入式(21)得到各Copula函数的AIC值：

其中k为c(·)中参数的数目，通过对比计算结果，AIC值最小的Copula函数即为最优函数。

4 基于包络思想和Vine Copula模型的机器人区间可靠性分析方法

4.1 包络函数展开点求解

包络函数主要思想是将时变可靠性问题转化为时无关可靠性问题且保持最大的转化精度。本文只需利用包络函数计算得到包络展开点即可。假设机器人末端位置误差阈值为ε，通过式(7)和式(8)可得到机器人位置误差极限状态功能函数g(X，t)=ε和g(X，t)=-ε，由文献[6]中的包络理论可知正向误差边界的包络方程G+(X)=0可表示为：

负向误差边界的包络方程G-(X)=0可表示为：

上式中，当时间t改变时，包络函数将由一系列时间相关的失效边界函数组成，故时间区间[t0，te]上的区间可靠度可定义为：

经以上分析可知，G(X)=0表示了误差极限状态函数组成的包络面方程，而包络函数呈隐式，求解异常困难。为解决此问题引入包络函数展开点思想，其原理是通过求解包络函数的极值失效点，利用有限个失效状态超平面来近似复杂的包络函数面，实践证明这是十分有效的方法。由于包络函数的高非线性，给极值展开点的求解带来难度，考虑到机器人误差变量的标准差远小于变量均值，可采用一阶泰勒展开的近似方法，在保证较高计算精度的同时也能有效降低计算难度。下面以正向边界包络函数为例，假设机器人各关节DH参数X服从正态分布，即X～N(μx，σ2x)。将误差函数g(X，t)在输入变量X的均值μx=｛μx1，ux2，…μXn｝处一阶泰勒展开，可得：

其中gd(t)为机器人期望运动输出。令σXi(i=1，2，…，n)（其中σX1为变量Xi的标准差），可得到标准正态U空间下的线性时变功能函数：

式中：

由式(27)可得变换后的功能函数L(U，t)对应的包络方程：

由式(32)即可推导出包络函数展开点计算公式：

同理，推导G-(X)=0的展开点计算方程为：

式(36)和式(37)均为时间t的一元方程，求解可得包络函数展开点。至此，区间失效问题已转化为对目标区间内有限个典型失效点（展开点）的联合失效求解问题。

4.2 基于包络思想和Vine Copula模型的区间可靠度求解方法

以上分析分别阐明了Vine Copula模型的建立以及包络函数展开点的推导计算过程。下面给出此联合方法在机器人轨迹运动可靠性分析中的操作步骤，具体流程如图2所示。本方法主要分以下三个阶段实现。

图2 联合方法求解区间可靠度流程图

第一阶段，建立机器人运动误差函数模型。通过机器人运动学分析可得到其运动学方程，由机器人位姿误差模型建立如式(8)的运动误差函数。

第二阶段，求解包络函数展开点。基于机器人运动误差函数以及式(20)、式(21)得到位置误差包络方程，由式(32)推导包络函数展开点计算方程，求解式(36)和式(37)即得到包络函数展开点ti。

第三阶段，构建联合概率密度函数。首先利用抽样技术获取各展开点时刻的运动误差样本Si；其次对数据样本进行处理，利用核密度估计方法获取数据的边缘密度函数fi(Si)和边缘分布函数Fi(Si)，基于赤池信息准则（AIC）和极大似然估计法选取最优二元Copula函数；最后由Vine Copula分解模型建立基于双变量Copula的多维联合概率密度函数并计算可靠度。

5 算例分析

以KUKA高精度机器人为例。该机器人对轨迹位置精度具有较高要求，其本体由6个旋转关节组成，关节DH参数如表1所示，各关节运动规律[13]如下：

关节运动周期T取10s。假设机械臂为刚性杆，以各关节参数取表1中名义值时的末端运动位置为期望输出，以误差函数等于给定误差阈值ε作为极限状态功能函数，求解机器人末端参考点轨迹的时变可靠度值。

表1 KUKA机器人DH参数

5.1 算例1

考虑机器人制造和装配误差，对机器人关节DH参数a和d引入N(0，0.052)的正态分布误差，误差单位为mm，求解机器人末端位置在时间区间[0，2.5]秒上的区间可靠度值。具体求解步骤如下：

步骤1：首先确定研究目标为机械手末端参考点在X轴坐标分量上的运动位置精度，设定误差阈值ε=0.08mm，利用包络法求取目标区间[0，2.5]上的包络展开点，得到时间序列t=[0，0.5677，1.5423]。

步骤2：利用正态分布随机数生成和抽样技术，在各失效时刻取1000组运动误差数据作为样本，不同时刻间的误差样本散点图分别如图3和图4所示。从图中可看出误差变量之间呈不完全线性相关关系，需使用匹配的Copula函数来表征变量间的相关性特征。

图3 t1和t3时刻运动误差样本散点图

图4 t2和t3时刻运动误差样本散点图

步骤3：利用高斯核密度估计法获取极值失效点处误差样本的边缘概率密度函数和边缘分布函数，作为后面选取Copula函数的依据。通过赤池信息准则（AIC）以及MLE方法选取最优Copula函数并估计其相关参数，筛选结果如表2所示。相应的二元Copula密度函数如图5和图6所示。由密度函数图可见，t-Copula较正态Copula具有更厚的尾部，更能反映误差数据之间的尾部相关信息。

图5 Gaussian Copula密度函数图

图6 t-Copula密度函数图

表2 Copula函数参数值及AIC值

步骤4：联合概率密度函数的计算。将步骤3得到的参数代入对应的Copula密度函数解析式，通过式(16)得到三维联合概率密度函数f(x1，x2，x3)，最终的可靠度可由下式计算：

为了观察不同误差阈值对区间可靠度的影响，改变ε的取值进行多次计算，并对比基于多元正态累积分布函数求解的传统包络法计算结果。最后利用蒙特卡洛模拟方法进行验证计算，仿真次数为20万。可靠度曲线如图7所示。从图中可以看出，采用本论文方法求解得到的区间可靠度值和蒙特卡洛仿真结果非常接近，且在计算精度上较传统包络法有显著提高。

图7 KUKA机器人轴向区间可靠度

5.2 算例2

在算例1的基础上，只改变机器人连杆长度和连杆偏距的参数分布类型，其他条件均不变。具体分布参数如表3所示，表3中各变量均值分别对应表1中各参数值。

表3 机器人连杆长度和连杆偏距参数及其分布

分别计算不同误差阈值下的区间可靠度，利用蒙特卡洛模拟20万次进行验证，结果如图8所示。由图可见，本文方法在计算精度上与蒙特卡洛仿真结果相一致，证实了改进方法可有效处理机器人末端位置精度分析中的误差变量非正态问题。

图8 KUKA机器人轴向区间可靠度

6 结语

本文以时变可靠性分析方法——包络法为基础，引入Vine Copula函数建立机器人运动轨迹失效概率模型。结合算例，结论如下：

1）本文利用KUKA机器人验证了基于包络思想的Vine Copula方法在机器人时变可靠性分析中的可行性，为实际工程分析提供了理论支持。

2）Vine Copula函数模型解决了传统包络法在求解机器人运动区间可靠度时受正态分布误差制约的问题，同时展现出较高的求解精度。

3）机器人运动轨迹可靠性在一些高端制造业生产中 尤为重要，而影响机器人运动轨迹精度的因素往往具有不确定性，本方法可以有效处理误差变量分布类型多样和多相关问题，更贴合实际工程分析。本文最终给出了机器人末端轨迹任意轴向的区间可靠度计算方法，对于机器人运动轨迹的系统可靠度而言，需要考虑多个轴向失效模式之间的联合失效相关问题，可作为下一步研究的重点方向。