基于Vine Copula的機器人運動精度可靠性分析*

劉春秋，周金宇，莊百亮，丁 力

（1.江蘇理工學(xué)院 機械工程學(xué)院，常州 213001；2.金陵科技學(xué)院 機電工程學(xué)院，南京 211169；3.機械科學(xué)研究總院江蘇分院有限公司，常州 213001)

0 引言

工業(yè)生產(chǎn)中，機器人受自身和使用環(huán)境影響，隨著時間的推移絕對定位精度逐漸下降。尤其對于高精度機器人，精度的下降無疑會對正常作業(yè)和生產(chǎn)造成不利影響。機器人運動精度退化的主要表現(xiàn)是末端執(zhí)行器的實際輸出位姿在理想輸出附近波動，且造成精度退化的原因多具有不確定性，如操作臂的桿件變形、關(guān)節(jié)磨損、制造誤差和裝配間隙等。機器人運動精度可靠性研究，就是衡量誤差對運動精度的影響程度，用概率的方法分析機械手在指定公差范圍內(nèi)到達指定位置或?qū)崿F(xiàn)預(yù)期運動軌跡的概率。

機器人運動精度一般分為單一點精度和軌跡精度，軌跡精度又可稱為區(qū)間運動精度，它是綜合評價機器人一段時間軌跡上的運動輸出是否滿足預(yù)期的一種方式。早在20世紀(jì)90年代，Bhatti[1]就開展了工業(yè)機器人運動相關(guān)概率建模并提出了機器人運動精度可靠性概念。隨著機器人工業(yè)化普及和精度要求的提高，此問題得到重視和研究。對于機器人點精度可靠性分析，廣大學(xué)者引入傳統(tǒng)的結(jié)構(gòu)可靠性分析方法，如一階可靠度方法[2]，得到了較為理想的計算精度。對于機器人區(qū)間運動精度的可靠性分析，只是在近些年才有所發(fā)展。目前區(qū)間（時變）可靠性分析法主要有極值法、首次穿越法、聯(lián)合概率方法以及包絡(luò)法。對于極值法，其原理是通過確定誤差函數(shù)的極值分布來估計區(qū)間可靠度。如MAHESH[3]將機器人軌跡離散為若干個點，利用最大熵原理推導(dǎo)機器人離散點的誤差極值分布，并以串聯(lián)系統(tǒng)可靠性模型計算機器人目標(biāo)軌跡的運動可靠度。首次穿越理論是以運動誤差第一次跨越誤差閾值的概率來估計區(qū)間失效概率。如Zhang[4]針對平面函數(shù)生成機構(gòu)提出了均值跨越率法，推導(dǎo)了上下跨越率公式，并將一次二階矩（FOSM）法推廣到時變可靠性分析中，此方法對于服從正態(tài)分布的小方差維數(shù)變量是準(zhǔn)確的。對于平面軌跡機構(gòu)的時變可靠性分析，陳放提出了聯(lián)合概率方法，文中認為機構(gòu)區(qū)間運動精度可靠度是區(qū)間內(nèi)各離散點處運動誤差落在誤差限范圍內(nèi)的綜合效應(yīng)[5]，作者針對正態(tài)分布的誤差，推導(dǎo)了軌跡離散點間的多元正態(tài)累積分布函數(shù)計算公式。由于此方法并未限定軌跡離散點的數(shù)量，故在計算精度上需進行預(yù)匹配。此外，為了提高時變可靠性分析的準(zhǔn)確性，Du和Zhang[6，7]等提出了基于包絡(luò)思想的時變可靠性分析方法，論文基于包絡(luò)原理計算并篩選出時間區(qū)間上的最易失效點，通過建立失效點間的多變量正態(tài)分布累積分布函數(shù)，計算區(qū)間失效概率。文中考慮了平面連桿機構(gòu)的結(jié)構(gòu)誤差和尺寸誤差，在誤差服從正態(tài)分布的情況下計算精度較高，但在變量非正態(tài)時產(chǎn)生局限性。為解決這一問題，Wang[8]在包絡(luò)思想的基礎(chǔ)上，引入Vine Copula模型來處理更多的誤差分布類型，此方法只需要誤差變量的邊緣分布和Copula連接函數(shù)即可得到多維聯(lián)合概率密度函數(shù)，且在平面連桿機構(gòu)中驗證了此方法的精度。以上研究主要針對平面函數(shù)生成機構(gòu)，且獲得了較理想的求解效率和精度。但對于空間軌跡機構(gòu)如機器人的運動精度可靠性分析問題，亟需拓展相關(guān)理論和方法。本文借鑒平面連桿機構(gòu)時變可靠性研究成果，利用Vine Copula函數(shù)和包絡(luò)法思想展開機器人軌跡的運動可靠性分析，以機器人末端位置精度為目標(biāo)求解區(qū)間可靠度。

1 機器人運動誤差模型

1.1 機器人運動學(xué)分析

機器人空間位姿建模的方法主要有矩陣法、矢量法和旋量法等。本文選用較為常用的矩陣法，此方法基本原理為：對于具有多個自由度的串聯(lián)機械手，在其每個連桿上建立關(guān)節(jié)坐標(biāo)系，并定義變換參數(shù)θi，di，ai-1，αi-1，通過旋轉(zhuǎn)和平移變換理論來獲得用于描述任意相鄰連桿i-1和i 之間變換關(guān)系的矩陣i-1iT。對于改進的DH參數(shù)模型，其一般表達式為

其中，cθi=cosθi，sθi=sinθi，θi為機器人各關(guān)節(jié)角，di為機器人相鄰連桿間的偏距，ai-1為各連桿長度，αi-1為描述相鄰關(guān)節(jié)軸相對位置的轉(zhuǎn)角。

通過變換矩陣的組合，可實現(xiàn)基坐標(biāo)空間到關(guān)節(jié)坐標(biāo)空間的轉(zhuǎn)化。機械手末端參考點在基坐標(biāo)空間內(nèi)的位姿矩陣可表示為：

式中，矩陣的前3行前3列為操作臂末端的姿態(tài)矩陣，第4列前3個元素分別為末端參考點在基坐標(biāo)中的位置分量。通過以上運動學(xué)分析得到的末端參考點位姿模型，即為機器人運動學(xué)分析的基本方程。

1.2 機器人運動誤差模型建立

影響機器人運動精度的各種不確定因素主要通過機器人關(guān)節(jié)空間的映射，到達機械臂末端參考點。考慮機器人關(guān)節(jié)空間的變換規(guī)律，選取機器人D-H參數(shù)作為誤差輸入變量。對于空間中任意相鄰關(guān)節(jié)i-1和i，由矩陣微分變換理論[9]對變換矩陣Ti求微分，可得

基于矩陣變換原理，對式(2)中的0nT求偏導(dǎo)并保留一階微分量，可得：

由式(3)知dTn=TnδTn，代入式(5)并化簡可得機器人末端誤差矩陣：

式(6)即為機器人末端位姿誤差矩陣，可簡化為：

式中，p和δ分別表示機器人末端位置誤差和姿態(tài)誤差，M1和M2分別為位置和姿態(tài)誤差系數(shù)矩陣，[ΔθΔdΔαΔa]T為各關(guān)節(jié)幾何參數(shù)誤差向量。

2 機器人軌跡運動可靠性分析

2.1 軌跡運動可靠性模型建立

機器人軌跡運動受結(jié)構(gòu)參數(shù)誤差等因素影響，實際軌跡在期望軌跡附近波動，軌跡運動可靠性也稱為區(qū)間可靠性，它是一種不同于點（靜態(tài)）可靠性的時變可靠性問題，結(jié)合上一節(jié)內(nèi)容，對于任意時刻t，機器人的位置運動誤差可表示為：

式中，gA(X，t)為實際運動位置，gd(t)為期望運動位置，Xi=[θi，di，αi-1，αi-1]T，（i=1，2，…，n)為各關(guān)節(jié)DH參數(shù)向量，n為機器人關(guān)節(jié)個數(shù)。對于任意時間區(qū)間[t0，te]上的一段軌跡，其軌跡精度可靠度可表示為區(qū)間內(nèi)每個點的位置誤差小于誤差閾值ε的概率，相應(yīng)的可靠度模型可表示為：

其中：Pr｛·｝表示概率，gx(·)表示軌跡坐標(biāo)在X軸上的運動誤差函數(shù)。相應(yīng)的失效概率可表示為：

其他坐標(biāo)分量上的可靠度計算模型同理可得。

2.2 區(qū)間可靠度求解方法

在平面連桿機構(gòu)中，傳統(tǒng)包絡(luò)法被用來計算區(qū)間可靠度，如文獻[6，7]。此方法思路是：首先尋找目標(biāo)時間區(qū)間內(nèi)的包絡(luò)展開點，然后建立展開點之間的多元正態(tài)累積分布函數(shù)來求解區(qū)間可靠度。此方法在輸入變量服從正態(tài)分布的情況下是有效的，其它分布情況下將變得復(fù)雜甚至失效，在機器人工程實際中，影響其運動精度的誤差分布類型具有不確定性，如在非正態(tài)誤差變量情形，傳統(tǒng)包絡(luò)法將失效。因此本文在包絡(luò)法的基礎(chǔ)上對其進行改進，即在求得包絡(luò)函數(shù)展開點之后，引入Vine Copula函數(shù)模型來建立展開點之間的聯(lián)合概率密度函數(shù)。此方法較傳統(tǒng)包絡(luò)法具有以下優(yōu)點：一是此方法不受誤差變量分布類型的影響，對于變量非正態(tài)情形仍然適用；二是Copula函數(shù)描述的秩相關(guān)系數(shù)和尾相關(guān)系數(shù)對比多元正態(tài)累積分布函數(shù)中的積矩相關(guān)系數(shù)，是更為合理的相關(guān)性測度，可以較完整的反應(yīng)變量間的相關(guān)特征，有助于提高求解精度。

3 軌跡運動可靠性分析中的Vine Copula模型

3.1 Copula函數(shù)理論

由Sklar定理[11]可知：多維隨機向量X=(x1，x2，…，xn)的n維聯(lián)合分布F(x1，x2，…，xn)可表示為所有一維邊緣分布函數(shù)ui=Fi(xi)(i=1，…，n)構(gòu)成的多元Copula函數(shù)C(u)，即：

通過上式可導(dǎo)出其概率密度：

式中，fi(·)表示各變量的邊緣概率密度，i=1，2，…，n。

3.2 Vine Copula函數(shù)模型

二元Copula函數(shù)模型簡單，便于計算，但只能處理二維變量間的相關(guān)問題。對于多變量相關(guān)情況，多元Copula函數(shù)雖然也可以使用，但是其求解難度較高，分析較復(fù)雜。所以本文引入由二元Copula函數(shù)組成的Vine Copula模型來建立多維隨機變量間的聯(lián)合概率分布。在Vine Copula模型中，Vine是一種描述聯(lián)合概率分布的圖形結(jié)構(gòu)，用于締造高維分布中的約束鏈接，使得該模型可以利用輸入變量的邊緣分布和Copula連接函數(shù)來建立多維聯(lián)合分布。此函數(shù)模型將變量邊緣分布和變量間相關(guān)結(jié)構(gòu)分開研究，簡化了多變量概率建模的過程[10]。對于任意維Copula函數(shù)，Vine Copula分解模型都可將其分解為多個二元Copula函數(shù)的集合。下面以3維D-Vine模型為例，圖1為模型樹形圖。T1層列舉了所有輸入變量；T2層為Copula連接層，表示二元Copula函數(shù)將上方數(shù)組兩兩連接；T3層為條件連接層，表示Copula函數(shù)將上方條件數(shù)組連接。其中，T2層受隨機變量間相關(guān)性的影響最大，越往下，影響越弱，所以由初始層樣本選取第二層的Copula函數(shù)至關(guān)重要。

圖1 3維D-vine模型

對于三維變量情況，假設(shè)隨機變量向量X=(x1，x2，x3)，經(jīng)圖1輸出的聯(lián)合概率密度函數(shù)可表示為：

式中，c12(·)，c23(·)，c13/2(·)為雙變量Copula密度函數(shù)。c13/2(·)中的條件隨機變量計算如下：

其中c12(·)，c23(·)，為Copula分布函數(shù)。同樣的，對于n維隨機向量X=(x1，x2，…，xn)，其聯(lián)合概率密度函數(shù)可表示為：

通過以上分析可知，雙變量Copula函數(shù)是連接隨機變量和建立聯(lián)合分布模型的關(guān)鍵。下面給出二元Copula函數(shù)的選取方法，本文利用極大似然估計（MLE）法來獲取目標(biāo)Copula函數(shù)的相關(guān)參數(shù)，并以赤池信息準(zhǔn)則[12]（AIC）作為最優(yōu)Copula函數(shù)的選取標(biāo)準(zhǔn)。對于任意隨機變量X1和X2，根據(jù)樣本信息，可建立對數(shù)似然函數(shù)：

其中c(·)為Copula密度函數(shù)，α為Copula函數(shù)的參數(shù)。，(i=1，2，…，N)是隨機變量X1和X2中的一對樣本，N為數(shù)據(jù)樣本容量。其中參數(shù)α可由最大化對數(shù)似然函數(shù)得到：

對所有候選Copula函數(shù)進行上述估計，得到各參數(shù)集，以變量樣本對和各候選Copula參數(shù)集作為輸入，代入式(21)得到各Copula函數(shù)的AIC值：

其中k為c(·)中參數(shù)的數(shù)目，通過對比計算結(jié)果，AIC值最小的Copula函數(shù)即為最優(yōu)函數(shù)。

4 基于包絡(luò)思想和Vine Copula模型的機器人區(qū)間可靠性分析方法

4.1 包絡(luò)函數(shù)展開點求解

包絡(luò)函數(shù)主要思想是將時變可靠性問題轉(zhuǎn)化為時無關(guān)可靠性問題且保持最大的轉(zhuǎn)化精度。本文只需利用包絡(luò)函數(shù)計算得到包絡(luò)展開點即可。假設(shè)機器人末端位置誤差閾值為ε，通過式(7)和式(8)可得到機器人位置誤差極限狀態(tài)功能函數(shù)g(X，t)=ε和g(X，t)=-ε，由文獻[6]中的包絡(luò)理論可知正向誤差邊界的包絡(luò)方程G+(X)=0可表示為：

負向誤差邊界的包絡(luò)方程G-(X)=0可表示為：

上式中，當(dāng)時間t改變時，包絡(luò)函數(shù)將由一系列時間相關(guān)的失效邊界函數(shù)組成，故時間區(qū)間[t0，te]上的區(qū)間可靠度可定義為：

經(jīng)以上分析可知，G(X)=0表示了誤差極限狀態(tài)函數(shù)組成的包絡(luò)面方程，而包絡(luò)函數(shù)呈隱式，求解異常困難。為解決此問題引入包絡(luò)函數(shù)展開點思想，其原理是通過求解包絡(luò)函數(shù)的極值失效點，利用有限個失效狀態(tài)超平面來近似復(fù)雜的包絡(luò)函數(shù)面，實踐證明這是十分有效的方法。由于包絡(luò)函數(shù)的高非線性，給極值展開點的求解帶來難度，考慮到機器人誤差變量的標(biāo)準(zhǔn)差遠小于變量均值，可采用一階泰勒展開的近似方法，在保證較高計算精度的同時也能有效降低計算難度。下面以正向邊界包絡(luò)函數(shù)為例，假設(shè)機器人各關(guān)節(jié)DH參數(shù)X服從正態(tài)分布，即X～N(μx，σ2x)。將誤差函數(shù)g(X，t)在輸入變量X的均值μx=｛μx1，ux2，…μXn｝處一階泰勒展開，可得：

其中g(shù)d(t)為機器人期望運動輸出。令σXi(i=1，2，…，n)（其中σX1為變量Xi的標(biāo)準(zhǔn)差），可得到標(biāo)準(zhǔn)正態(tài)U空間下的線性時變功能函數(shù)：

式中：

由式(27)可得變換后的功能函數(shù)L(U，t)對應(yīng)的包絡(luò)方程：

由式(32)即可推導(dǎo)出包絡(luò)函數(shù)展開點計算公式：

同理，推導(dǎo)G-(X)=0的展開點計算方程為：

式(36)和式(37)均為時間t的一元方程，求解可得包絡(luò)函數(shù)展開點。至此，區(qū)間失效問題已轉(zhuǎn)化為對目標(biāo)區(qū)間內(nèi)有限個典型失效點（展開點）的聯(lián)合失效求解問題。

4.2 基于包絡(luò)思想和Vine Copula模型的區(qū)間可靠度求解方法

以上分析分別闡明了Vine Copula模型的建立以及包絡(luò)函數(shù)展開點的推導(dǎo)計算過程。下面給出此聯(lián)合方法在機器人軌跡運動可靠性分析中的操作步驟，具體流程如圖2所示。本方法主要分以下三個階段實現(xiàn)。

圖2 聯(lián)合方法求解區(qū)間可靠度流程圖

第一階段，建立機器人運動誤差函數(shù)模型。通過機器人運動學(xué)分析可得到其運動學(xué)方程，由機器人位姿誤差模型建立如式(8)的運動誤差函數(shù)。

第二階段，求解包絡(luò)函數(shù)展開點?；跈C器人運動誤差函數(shù)以及式(20)、式(21)得到位置誤差包絡(luò)方程，由式(32)推導(dǎo)包絡(luò)函數(shù)展開點計算方程，求解式(36)和式(37)即得到包絡(luò)函數(shù)展開點ti。

第三階段，構(gòu)建聯(lián)合概率密度函數(shù)。首先利用抽樣技術(shù)獲取各展開點時刻的運動誤差樣本Si；其次對數(shù)據(jù)樣本進行處理，利用核密度估計方法獲取數(shù)據(jù)的邊緣密度函數(shù)fi(Si)和邊緣分布函數(shù)Fi(Si)，基于赤池信息準(zhǔn)則（AIC）和極大似然估計法選取最優(yōu)二元Copula函數(shù)；最后由Vine Copula分解模型建立基于雙變量Copula的多維聯(lián)合概率密度函數(shù)并計算可靠度。

5 算例分析

以KUKA高精度機器人為例。該機器人對軌跡位置精度具有較高要求，其本體由6個旋轉(zhuǎn)關(guān)節(jié)組成，關(guān)節(jié)DH參數(shù)如表1所示，各關(guān)節(jié)運動規(guī)律[13]如下：

關(guān)節(jié)運動周期T取10s。假設(shè)機械臂為剛性桿，以各關(guān)節(jié)參數(shù)取表1中名義值時的末端運動位置為期望輸出，以誤差函數(shù)等于給定誤差閾值ε作為極限狀態(tài)功能函數(shù)，求解機器人末端參考點軌跡的時變可靠度值。

表1 KUKA機器人DH參數(shù)

5.1 算例1

考慮機器人制造和裝配誤差，對機器人關(guān)節(jié)DH參數(shù)a和d引入N(0，0.052)的正態(tài)分布誤差，誤差單位為mm，求解機器人末端位置在時間區(qū)間[0，2.5]秒上的區(qū)間可靠度值。具體求解步驟如下：

步驟1：首先確定研究目標(biāo)為機械手末端參考點在X軸坐標(biāo)分量上的運動位置精度，設(shè)定誤差閾值ε=0.08mm，利用包絡(luò)法求取目標(biāo)區(qū)間[0，2.5]上的包絡(luò)展開點，得到時間序列t=[0，0.5677，1.5423]。

步驟2：利用正態(tài)分布隨機數(shù)生成和抽樣技術(shù)，在各失效時刻取1000組運動誤差數(shù)據(jù)作為樣本，不同時刻間的誤差樣本散點圖分別如圖3和圖4所示。從圖中可看出誤差變量之間呈不完全線性相關(guān)關(guān)系，需使用匹配的Copula函數(shù)來表征變量間的相關(guān)性特征。

圖3 t1和t3時刻運動誤差樣本散點圖

圖4 t2和t3時刻運動誤差樣本散點圖

步驟3：利用高斯核密度估計法獲取極值失效點處誤差樣本的邊緣概率密度函數(shù)和邊緣分布函數(shù)，作為后面選取Copula函數(shù)的依據(jù)。通過赤池信息準(zhǔn)則（AIC）以及MLE方法選取最優(yōu)Copula函數(shù)并估計其相關(guān)參數(shù)，篩選結(jié)果如表2所示。相應(yīng)的二元Copula密度函數(shù)如圖5和圖6所示。由密度函數(shù)圖可見，t-Copula較正態(tài)Copula具有更厚的尾部，更能反映誤差數(shù)據(jù)之間的尾部相關(guān)信息。

圖5 Gaussian Copula密度函數(shù)圖

圖6 t-Copula密度函數(shù)圖

表2 Copula函數(shù)參數(shù)值及AIC值

步驟4：聯(lián)合概率密度函數(shù)的計算。將步驟3得到的參數(shù)代入對應(yīng)的Copula密度函數(shù)解析式，通過式(16)得到三維聯(lián)合概率密度函數(shù)f(x1，x2，x3)，最終的可靠度可由下式計算：

為了觀察不同誤差閾值對區(qū)間可靠度的影響，改變ε的取值進行多次計算，并對比基于多元正態(tài)累積分布函數(shù)求解的傳統(tǒng)包絡(luò)法計算結(jié)果。最后利用蒙特卡洛模擬方法進行驗證計算，仿真次數(shù)為20萬?？煽慷惹€如圖7所示。從圖中可以看出，采用本論文方法求解得到的區(qū)間可靠度值和蒙特卡洛仿真結(jié)果非常接近，且在計算精度上較傳統(tǒng)包絡(luò)法有顯著提高。

圖7 KUKA機器人軸向區(qū)間可靠度

5.2 算例2

在算例1的基礎(chǔ)上，只改變機器人連桿長度和連桿偏距的參數(shù)分布類型，其他條件均不變。具體分布參數(shù)如表3所示，表3中各變量均值分別對應(yīng)表1中各參數(shù)值。

表3 機器人連桿長度和連桿偏距參數(shù)及其分布

分別計算不同誤差閾值下的區(qū)間可靠度，利用蒙特卡洛模擬20萬次進行驗證，結(jié)果如圖8所示。由圖可見，本文方法在計算精度上與蒙特卡洛仿真結(jié)果相一致，證實了改進方法可有效處理機器人末端位置精度分析中的誤差變量非正態(tài)問題。

圖8 KUKA機器人軸向區(qū)間可靠度

6 結(jié)語

本文以時變可靠性分析方法——包絡(luò)法為基礎(chǔ)，引入Vine Copula函數(shù)建立機器人運動軌跡失效概率模型。結(jié)合算例，結(jié)論如下：

1）本文利用KUKA機器人驗證了基于包絡(luò)思想的Vine Copula方法在機器人時變可靠性分析中的可行性，為實際工程分析提供了理論支持。

2）Vine Copula函數(shù)模型解決了傳統(tǒng)包絡(luò)法在求解機器人運動區(qū)間可靠度時受正態(tài)分布誤差制約的問題，同時展現(xiàn)出較高的求解精度。

3）機器人運動軌跡可靠性在一些高端制造業(yè)生產(chǎn)中 尤為重要，而影響機器人運動軌跡精度的因素往往具有不確定性，本方法可以有效處理誤差變量分布類型多樣和多相關(guān)問題，更貼合實際工程分析。本文最終給出了機器人末端軌跡任意軸向的區(qū)間可靠度計算方法，對于機器人運動軌跡的系統(tǒng)可靠度而言，需要考慮多個軸向失效模式之間的聯(lián)合失效相關(guān)問題，可作為下一步研究的重點方向。