章家?guī)r, 許治順, 王 勝, 馮旭剛

(安徽工業(yè)大學(xué) 電氣與信息工程學(xué)院，安徽 馬鞍山 243032)

關(guān)節(jié)臂式坐標(biāo)測量機是一種非笛卡兒式的柔性測量系統(tǒng)，通常在大型加工現(xiàn)場對不便移動的復(fù)雜零部件進(jìn)行幾何尺寸的快速測量[1]，因其可操作性、靈活性和便捷性在工業(yè)生產(chǎn)與科學(xué)研究中得到廣泛的應(yīng)用。由于其串聯(lián)式連桿機構(gòu)，測量機的測量誤差會被累積并放大，同時機械加工精度、裝配誤差、磨損及環(huán)境因素等都對各連桿的幾何參數(shù)造成影響，與傳統(tǒng)的正交式坐標(biāo)測量機相比，整體測量精度還存在較大差距。因此，測量機的誤差模型的分析及參數(shù)標(biāo)定方法對于減少測量誤差有著重大意義。

目前已有眾多國內(nèi)外學(xué)者針對關(guān)節(jié)臂式坐標(biāo)測量機參數(shù)標(biāo)定方法做了多方面研究。一類通過借助更高精度的測量儀器，如激光跟蹤儀、三坐標(biāo)測量機[2]，或布置復(fù)雜的裝夾裝置[3]。例如：Acero等[4]研究了激光跟蹤儀作為參考儀器在AACMM參數(shù)標(biāo)定中的可行性，但這類方法普遍成本高昂，系統(tǒng)復(fù)雜，時間效率較低且容易產(chǎn)生二次誤差[5]。另一類通過算法進(jìn)行參數(shù)標(biāo)定，一般為最小二乘法和智能尋優(yōu)算法，例如：Santolaria等[6]使用一維球列作為標(biāo)準(zhǔn)件，采用Levenberg-Marquarat法(最小二乘法)對Faro的關(guān)節(jié)式坐標(biāo)測量機進(jìn)行了標(biāo)定；Seyedhosseini等[7]應(yīng)用改進(jìn)的模擬退火算法對測量臂運動學(xué)參數(shù)的標(biāo)定值進(jìn)行了辨識。但采用標(biāo)定算法對參數(shù)進(jìn)行辨識時，常用的標(biāo)定算法如最小二乘法在求解過程中存在初始值是非可行點時無法收斂、參數(shù)較多的復(fù)雜計算易產(chǎn)生積累誤差等問題，而尋優(yōu)算法的局部搜索能力較差、容易陷入局部最小值。

針對關(guān)節(jié)式坐標(biāo)測量機的參數(shù)標(biāo)定問題，筆者提出了一種混合優(yōu)化算法對測量機結(jié)構(gòu)參數(shù)進(jìn)行辨識。首先對最小二乘法求解中的Jacobian矩陣進(jìn)行變換消除其中冗余參數(shù)，通過設(shè)定判定準(zhǔn)則并實現(xiàn)模擬退火算法和最小二乘法的混合，結(jié)合兩者全局優(yōu)化性能和高效局部的優(yōu)化能力進(jìn)行參數(shù)標(biāo)定，提高測量機的測量精度。

1 關(guān)節(jié)臂式坐標(biāo)測量機模型

關(guān)節(jié)臂式坐標(biāo)測量機是仿人手臂的六自由度非正交坐標(biāo)測量機，由機座、關(guān)節(jié)臂和測頭通過旋轉(zhuǎn)關(guān)節(jié)串聯(lián)連接構(gòu)成[8]，結(jié)構(gòu)參數(shù)名義值如表1所示，其結(jié)構(gòu)如圖1所示。

圖1 關(guān)節(jié)臂式坐標(biāo)測量機結(jié)構(gòu)示意圖Fig. 1 Configuration of articulated arm coordinate measuring machine

表1 結(jié)構(gòu)參數(shù)的標(biāo)稱值

1.1 D-H模型

1956年，Denavit和Hartenberg提出D-H方法描述相鄰連桿的變換關(guān)系，使用4個參數(shù)進(jìn)行描述，分別是關(guān)節(jié)旋轉(zhuǎn)角θi，桿件的長度ai，桿件偏移量di，桿件扭轉(zhuǎn)角度βi?；贒-H模型的關(guān)節(jié)坐標(biāo)轉(zhuǎn)換示意圖如圖2所示。

圖2 基于D-H模型的關(guān)節(jié)坐標(biāo)轉(zhuǎn)換示意圖Fig. 2 Coordinate transformation diagram by D-H model

根據(jù)圖2可知，坐標(biāo)系(xi-1,yi-1,zi-1)轉(zhuǎn)至下個坐標(biāo)系(xi,yi,zi)的坐標(biāo)的變換過程，需要經(jīng)過2次旋轉(zhuǎn)2次平移得到，坐標(biāo)轉(zhuǎn)換關(guān)系為

Ai-1,i= Rot(zn,θn+1)Trans(0 , 0,dn+1) Trans(an+1, 0, 0)Rot(xn+1,βn+1)。

(1)

相鄰坐標(biāo)系的齊次變換矩Ai-1,i為

(2)

此時，由于第7個坐標(biāo)系是以測頭為中心，由第6個坐標(biāo)系平移，設(shè)測頭中心在坐標(biāo)系下坐標(biāo)為(Bx,By,Bz)，坐標(biāo)轉(zhuǎn)換矩陣B為

(3)

基于D-H方法的關(guān)節(jié)臂式坐標(biāo)測量機的坐標(biāo)系統(tǒng)如圖3所示，根據(jù)公式(3)，六自由度關(guān)節(jié)臂式坐標(biāo)測量機末端測頭相對于基座坐標(biāo)系的數(shù)學(xué)模型可表達(dá)為

圖3 基于D-H模型的關(guān)節(jié)臂式坐標(biāo)測量機的坐標(biāo)系統(tǒng)簡圖Fig. 3 Coordinate system diagram of AACMM based on the D-H model

(4)

1.2 參數(shù)誤差模型

從公式(4)可以看出，測量機末端測頭的坐標(biāo)值通過參數(shù)進(jìn)行描述，則末端測頭中心坐標(biāo)誤差與參數(shù)誤差{Δai，Δdi，Δβi，Δθi，ΔBi}有關(guān)。從式中可知共有27個參數(shù)誤差，Δai，Δdi，Δβi分別是關(guān)節(jié)的長度誤差、關(guān)節(jié)偏移量誤差和關(guān)節(jié)扭轉(zhuǎn)角度誤差，Δθi表示初始零點處的關(guān)節(jié)變量θi的零位誤差，ΔBi是測頭的裝配誤差。根據(jù)理論測量模型的參數(shù)誤差種類建立總誤差模型為

(5)

(6)

將公式(5)中的測頭坐標(biāo)改寫為函數(shù)形式，對3組函數(shù)進(jìn)行全微分?jǐn)?shù)學(xué)解析，并將左右兩側(cè)分別展開寫成矩陣的形式，可得

(7)

式(7)等號右側(cè)3×27的Jacobian矩陣可用J表示，假設(shè)全部結(jié)構(gòu)參數(shù)誤差足夠小，公式(6)可近似為

(8)

對式(8)進(jìn)行簡化：ΔM=J·ΔS，ΔS為27×1的誤差參數(shù)矢量，ΔM表示為測頭坐標(biāo)誤差模型。

2 目標(biāo)函數(shù)和標(biāo)定算法

關(guān)節(jié)式坐標(biāo)測量機在測量同一點時，根據(jù)理論模型在不同的姿態(tài)時坐標(biāo)值應(yīng)該保持不變[9]，但結(jié)構(gòu)參數(shù)與理論參數(shù)因多種因素存在誤差，此外，測量誤差還與空間位置有關(guān)[10]，所以標(biāo)定時測量機應(yīng)該在不同的位置進(jìn)行多姿態(tài)的測量，坐標(biāo)值的誤差大小和波動范圍都是測量機精度的體現(xiàn)。因此筆者結(jié)合兩者作為目標(biāo)函數(shù)，采用設(shè)計的混合算法求取結(jié)構(gòu)參數(shù)的實際值。

2.1 目標(biāo)函數(shù)

根據(jù)上述內(nèi)容，對錐窩進(jìn)行不同姿態(tài)的測量，記錄所得的角度數(shù)據(jù)得到坐標(biāo)值，設(shè)該點的真實值為多次測量數(shù)據(jù)的平均值，即

(9)

該點的誤差平均值是

(10)

該點的標(biāo)準(zhǔn)差為

(11)

關(guān)節(jié)臂式測量機的整體測量精度通過e和σ的組合來表示，目標(biāo)函數(shù)為

OA=e+3σ。

(12)

通過目標(biāo)函數(shù)的值來判定測量機計算所得的結(jié)構(gòu)參數(shù)是否接近真實值。

2.2 去除冗余參數(shù)的最小二乘法

根據(jù)公式(8)，運用最小二乘法計算可得

ΔS= (JT×J)-1×JT×ΔM。

(13)

根據(jù)式(13)計算出ΔS后對系統(tǒng)參數(shù)進(jìn)行修正，通過修正后的參數(shù)計算出新一輪的理論坐標(biāo)值及誤差矩陣，判定修正后參數(shù)的誤差小于設(shè)定的閾值ε，不滿足條件就繼續(xù)迭代；滿足條件則運算結(jié)束，輸出最終的參數(shù)。

當(dāng)上述公式中雅克比矩陣J為奇異陣時，如果矩陣中有部分參數(shù)為線性關(guān)系，會對參數(shù)誤差的求解造成干擾。因此對誤差模型ΔS的等式左側(cè)進(jìn)行變換分析：

ΔS= (JT×J)-1×JT×ΔM?ΔM?(JT×J)×ΔS=JT×ΔM。

(14)

令式(14)中K= (JT×J)-1，并對其進(jìn)行奇異值分解：

(15)

式中:P,Q為正交陣;Λ表示對角陣(Λ=diag(φ1,φ2,...,φr)),則對角陣的秩r即為Jacobian矩陣的秩。將矩陣K代回式(15)可得

(16)

根據(jù)對式(16)中的矩陣分析，得知結(jié)構(gòu)參數(shù)中有25-r個參數(shù)有線性關(guān)系。由于式(15)中K是對稱矩陣，那么有對應(yīng)的關(guān)系式QT=P-1，因此Q陣屬于旋轉(zhuǎn)矩陣，對式(16)的誤差ΔS模型進(jìn)行旋轉(zhuǎn)變換，將所有呈線性相關(guān)的結(jié)構(gòu)參數(shù)處在相同的零平面上[12]。隨后將QT中后25-r行的部分做初等行變換從而提取矩陣中有線性關(guān)系的誤差參數(shù)，最終所得結(jié)果為

Δa6=Bz·Δθ6，Δd6=-Bz·Δβ6。

(17)

由于Bz已知，因此公式(17)中參數(shù)a6，θ6，d6，β6無法同時進(jìn)行標(biāo)定，選擇2個參數(shù)作為冗余參數(shù)并在Jacobian矩陣及參數(shù)誤差矩陣ΔS中去除對應(yīng)參數(shù)的矩陣，可得新的結(jié)構(gòu)參數(shù)誤差公式，其中矩陣ΔS為r行1列，矩陣J為3n行r列，矩陣ΔS為3n行1列，即

(18)

2.3 混合退火最小二乘算法設(shè)計

模擬退火是模擬統(tǒng)計物理學(xué)中的固體退火過程[12-13]，適用于求解不同的非線性復(fù)雜問題，算法有解空間、能量函數(shù)和初始解3部分，它的解題思路是在對初始解隨機擾動產(chǎn)生新解，代入能量函數(shù)判斷是否滿足Metropolis準(zhǔn)則，根據(jù)冷卻系數(shù)讓溫度下降使能量函數(shù)達(dá)到平衡狀態(tài)； Metropolis抽樣準(zhǔn)則具有概率性地跳出局部最優(yōu)“陷阱”的優(yōu)點[14]。

為計算出優(yōu)化權(quán)值，將上文的目標(biāo)函數(shù)設(shè)為能量函數(shù)E=f=e+3σ。

算法步驟如下：

Step 1 設(shè)定初始值。給定初始溫度T0始權(quán)參數(shù)w(0)=w0，設(shè)置終止檢驗精度e，終止溫度Tmin，馬爾可夫鏈的鏈長L，令初始最優(yōu)解w*=w0，迭代次數(shù)i=0。

Step 2 產(chǎn)生新解。令ωβ=ω(k)+R×E產(chǎn)生新解，其中R為區(qū)間[-1, 1]的隨機數(shù),符合Cauchy分布。

Step 3 求優(yōu)化函數(shù)指標(biāo)。計算ΔE=E(wβ)-E[w(k)]。

Step 4 接受判斷。如果ΔE≥0，計算接受概率r=exp[-E(εβ)/T]，如果r>p，則w(k+1)=wβ，否則w(k+1)=w(k)，p為區(qū)間[0, 1]上的隨機數(shù)；如果ΔE<0，則w(k+1)=wβ，w*=wβ。

Step 5 穩(wěn)定性判別。k=k+1，如果k>L，則轉(zhuǎn)到step 5，否則轉(zhuǎn)到步驟2。

Step 6 降溫T=Ti+1=αTi，i=i+1。

Step 8 輸出最終最優(yōu)解w*，中止算法。

融合算法的關(guān)鍵在于如何選擇兩者算法的轉(zhuǎn)折點[15]，最通用的方法是選擇固定的溫度次數(shù)，雖然簡便但針對不同模型具有通用性，往往須通過大量的實際經(jīng)驗才能得到大致范圍。因此，采取判定適應(yīng)度的方式融合兩種算法，如果同一溫度的多個不同最優(yōu)值之間的差值小于某一設(shè)定值，那么結(jié)束模擬退火算法轉(zhuǎn)入最小二乘法。

適應(yīng)度公式為

轉(zhuǎn)換準(zhǔn)則為

由上可知，融合算法的過程如下：先采用SA算法計算，若滿足切換準(zhǔn)則表明進(jìn)入全局最優(yōu)解的區(qū)域，退出SA算法，將滿足準(zhǔn)則的解作為LM算法的初值代入，通過運算得到最優(yōu)解，具體的流程圖如圖4所示。

圖4 混合算法的流程圖Fig. 4 Flow chart of the hybrid algorithm

3 實驗分析

為了減少隨機誤差的影響，對錐窩中心點進(jìn)行測量，關(guān)節(jié)臂式坐標(biāo)測量機從不同方向上測量同一點(即改變各關(guān)節(jié)旋轉(zhuǎn)角度)獲得空間三維坐標(biāo)數(shù)據(jù)，測量結(jié)束后改變錐窩的空間位置重復(fù)上述步驟，測量示意圖如圖5所示。運用最小二乘法通過分析線性相關(guān)的結(jié)構(gòu)參數(shù)，將其中2個不參與標(biāo)定的參數(shù)設(shè)為理論值，而裝配測頭參數(shù)誤差(ΔBx,ΔBy,ΔBz)中前兩項為0，設(shè)Bz的值為100 mm，故整個關(guān)節(jié)臂式坐標(biāo)測量機的結(jié)構(gòu)參數(shù)誤差共有23個。為了減少溫度誤差的影響，室溫保持20 ℃，根據(jù)表1的名義參數(shù)和測得的角度數(shù)據(jù)得到坐標(biāo)值，用混合算法進(jìn)行誤差參數(shù)標(biāo)定，計算所得的結(jié)構(gòu)參數(shù)誤差如表2所示。

圖5 測量示意圖Fig. 5 Schematic diagram of the measurement

表2 基于混合算法標(biāo)定的各項結(jié)構(gòu)參數(shù)誤差

為了驗證所提算法得到的誤差參數(shù)的準(zhǔn)確性，將計算得出的參數(shù)誤差代入結(jié)構(gòu)參數(shù)名義值，得到修正后的測量機結(jié)構(gòu)參數(shù)。分別采用辨識前后的參數(shù)，對兩者的誤差進(jìn)行對比。關(guān)節(jié)臂式坐標(biāo)測量機對錐窩進(jìn)行測量獲得角度數(shù)據(jù)，每測10次數(shù)據(jù)后錐窩換一次位置，設(shè)該點坐標(biāo)真值為多次測量值的平均值，計算單點重復(fù)性誤差；此外，長度誤差也是測量機整體精度的重要部分，為了驗證標(biāo)定算法的有效性，關(guān)節(jié)臂式測量機先對錐窩進(jìn)行多次測量，然后改變錐窩的位置再測多次，獲得數(shù)據(jù)后計算空間每組兩點之間的長度測量誤差，公式為

(19)

圖7 參數(shù)辨識前后長度誤差Fig 7 Length error based on different structural parameters

由圖6～7分析可知，參數(shù)辨識前，測量機的單點重復(fù)性誤差范圍在1.571 ～2.136 mm，平均值為1.873 mm，標(biāo)準(zhǔn)差為0.106 mm；參數(shù)辨識后，測量機的單點重復(fù)性誤差在0.097 ～ 0.166 mm，平均值為0.127 mm，標(biāo)準(zhǔn)差為0.019 mm。參數(shù)辨識前長度誤差范圍在0.486 ～ 2.116 mm，平均值為1.038 mm，標(biāo)準(zhǔn)差為0.426 mm，參數(shù)辨識后長度誤差范圍在0.025 ～ 0.448 mm，平均值為0.092 mm，標(biāo)準(zhǔn)差為0.097 mm。實驗結(jié)果表明，混合優(yōu)化算法有效地提高了測量機的測量精度。

圖6 參數(shù)辨識前后單點重復(fù)性誤差Fig. 6 Single point repeatability error based on different structural parameters

4 結(jié) 論

在關(guān)節(jié)臂式坐標(biāo)測量機測量模型分析的基礎(chǔ)上對標(biāo)定算法進(jìn)行了研究，結(jié)論如下：

1)采用最小二乘法標(biāo)定結(jié)構(gòu)參數(shù)過程中，針對呈線性相關(guān)的誤差參數(shù)無法準(zhǔn)確地求出參數(shù)誤差解的問題，對Jacobian矩陣變換分析，去除了方程中的冗余參數(shù)，提高了標(biāo)定精度也減少了計算量。

2)基于模擬退火算法與最小二乘法的優(yōu)缺點，通過設(shè)定判定準(zhǔn)則來混合兩種算法，既充分利用了SA的全局優(yōu)化性能和最小二乘法算法的高效局部優(yōu)化特點，又解決了最小二乘法初值設(shè)定問題。

3)通過實驗結(jié)果分析，運用混合優(yōu)化算法補償了關(guān)節(jié)臂式坐標(biāo)測量機結(jié)構(gòu)參數(shù)誤差后，測量機的準(zhǔn)確性和重復(fù)性都得到了顯著的提升，測量誤差得到了進(jìn)一步的抑制。