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

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

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

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

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

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

關(guān)節(jié)臂式坐標(biāo)測(cè)量機(jī)是仿人手臂的六自由度非正交坐標(biāo)測(cè)量機(jī)，由機(jī)座、關(guān)節(jié)臂和測(cè)頭通過(guò)旋轉(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)測(cè)量機(jī)結(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個(gè)參數(shù)進(jìn)行描述，分別是關(guān)節(jié)旋轉(zhuǎn)角θi，桿件的長(zhǎng)度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)至下個(gè)坐標(biāo)系(xi,yi,zi)的坐標(biāo)的變換過(guò)程，需要經(jīng)過(guò)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)

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

(3)

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

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

(4)

1.2 參數(shù)誤差模型

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

(5)

(6)

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

(7)

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

(8)

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

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

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

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

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

(9)

該點(diǎn)的誤差平均值是

(10)

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

(11)

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

OA=e+3σ。

(12)

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

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

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

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

(13)

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

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

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

(14)

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

(15)

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

(16)

根據(jù)對(duì)式(16)中的矩陣分析，得知結(jié)構(gòu)參數(shù)中有25-r個(gè)參數(shù)有線性關(guān)系。由于式(15)中K是對(duì)稱矩陣，那么有對(duì)應(yīng)的關(guān)系式QT=P-1，因此Q陣屬于旋轉(zhuǎn)矩陣，對(duì)式(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無(wú)法同時(shí)進(jìn)行標(biāo)定，選擇2個(gè)參數(shù)作為冗余參數(shù)并在Jacobian矩陣及參數(shù)誤差矩陣ΔS中去除對(duì)應(yīng)參數(shù)的矩陣，可得新的結(jié)構(gòu)參數(shù)誤差公式，其中矩陣ΔS為r行1列，矩陣J為3n行r列，矩陣ΔS為3n行1列，即

(18)

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

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

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

算法步驟如下：

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

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

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

Step 4 接受判斷。如果ΔE≥0，計(jì)算接受概率r=exp[-E(εβ)/T]，如果r>p，則w(k+1)=wβ，否則w(k+1)=w(k)，p為區(qū)間[0, 1]上的隨機(jī)數(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)折點(diǎn)[15]，最通用的方法是選擇固定的溫度次數(shù)，雖然簡(jiǎn)便但針對(duì)不同模型具有通用性，往往須通過(guò)大量的實(shí)際經(jīng)驗(yàn)才能得到大致范圍。因此，采取判定適應(yīng)度的方式融合兩種算法，如果同一溫度的多個(gè)不同最優(yōu)值之間的差值小于某一設(shè)定值，那么結(jié)束模擬退火算法轉(zhuǎn)入最小二乘法。

適應(yīng)度公式為

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

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

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

3 實(shí)驗(yàn)分析

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

圖5 測(cè)量示意圖Fig. 5 Schematic diagram of the measurement

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

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

(19)

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

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

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

4 結(jié) 論

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

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

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

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