基于多項(xiàng)式混沌的機(jī)床幾何誤差靈敏度分析*

鄭華林，趙 興，胡 騰，魏小建，王小虎

（1.西南石油大學(xué)，成都 610500；2.石油天然氣裝備技術(shù)四川省科技資源共享服務(wù)平臺(tái)，成都 610500）

隨著我國(guó)工業(yè)技術(shù)的發(fā)展，國(guó)防、航空、航天、汽車等諸多領(lǐng)域?qū)C(jī)床加工精度的要求日益提高。機(jī)床的加工精度會(huì)受幾何誤差、熱誤差、切削力誘導(dǎo)誤差等因素的綜合影響。其中，幾何誤差是由機(jī)床零部件制造和裝配不精確等因素造成的，在誤差源中占有很大比例。由于幾何誤差項(xiàng)較多，逐項(xiàng)測(cè)量并補(bǔ)償全部誤差成本很高，有必要對(duì)其進(jìn)行靈敏度分析[1]。從而降低補(bǔ)償成本、提高效率，對(duì)優(yōu)化機(jī)床加工精度具有重要意義。

近年來，數(shù)控機(jī)床幾何誤差靈敏度分析方法發(fā)展迅速。Cheng等[2]通過旋量理論對(duì)五軸機(jī)床空間誤差進(jìn)行了建模，并通過改進(jìn)的Morris法對(duì)37項(xiàng)幾何誤差進(jìn)行了全局靈敏度分析。Li等[3]基于多體系統(tǒng)理論構(gòu)建了五軸機(jī)床誤差源與刀具位姿誤差之間的誤差映射，通過定義局部靈敏度指數(shù)、全局靈敏度指數(shù)和全局靈敏度波動(dòng)指數(shù)對(duì)機(jī)床誤差進(jìn)行了靈敏度分析。Zou等[4]基于蒙特卡洛法 （Monte Carlo simulation, MCS），求解Sobol靈敏度指數(shù)，量化了幾何誤差對(duì)加工精度的影響，并根據(jù)分析結(jié)果對(duì)關(guān)鍵幾何誤差元素進(jìn)行調(diào)整和分布。Fu等[5]根據(jù)坐標(biāo)系微分運(yùn)動(dòng)關(guān)系，建立了機(jī)床各軸的幾何誤差靈敏度矩陣，分析得出誤差關(guān)鍵軸為A軸和X軸，并驗(yàn)證其正確性。郭世杰等[6]采用拉丁超立方法 （Latin hypercube sampling，LHS）在整個(gè)參數(shù)空間內(nèi)抽樣，利用Spearman系數(shù)進(jìn)行幾何誤差相關(guān)性分析并辨識(shí)出影響機(jī)床精度的關(guān)鍵幾何誤差。王培桐等[7]采用截?cái)喔道锶~級(jí)數(shù)表征與位置有關(guān)的幾何誤差，將求解靈敏度指數(shù)轉(zhuǎn)換為求取傅里葉幅值，該方法具有較好的魯棒性。綜上所述，目前機(jī)床的靈敏度分析方法主要有偏微分法、Morris法、MCS法、LHS法、傅里葉幅值法等，這些方法普遍樣本需求量大且計(jì)算效率不高。

針對(duì)上述局限性，本文提出一種以多項(xiàng)式混沌展開 （Polynomial chaos expansion，PCE）為理論基礎(chǔ)的全局靈敏度分析方法。將雙轉(zhuǎn)臺(tái)五軸數(shù)控機(jī)床作為研究對(duì)象，運(yùn)用旋量理論建立空間誤差模型，構(gòu)建幾何誤差的PCE模型，通過Sobol指數(shù)表征41項(xiàng)誤差的關(guān)鍵程度。最后，與MCS法和LHS法進(jìn)行對(duì)比，驗(yàn)證其正確性與高效性。

1 五軸數(shù)控機(jī)床空間誤差建模

1.1 運(yùn)動(dòng)鏈拓?fù)涿枋?/h3>

本文的研究對(duì)象為AC型雙轉(zhuǎn)臺(tái)五軸數(shù)控機(jī)床，其主要結(jié)構(gòu)及運(yùn)動(dòng)鏈拓?fù)淙鐖D1所示。其中，機(jī)床坐標(biāo)系設(shè)于回轉(zhuǎn)工作臺(tái)的中心，即旋轉(zhuǎn)軸軸線交點(diǎn)處。根據(jù)該機(jī)床運(yùn)動(dòng)鏈拓?fù)淇芍D(zhuǎn)軸A、C彼此相鄰并與床身構(gòu)成工件運(yùn)動(dòng)鏈，平動(dòng)軸X、Y、Z與床身構(gòu)成刀具運(yùn)動(dòng)鏈。

圖1 雙轉(zhuǎn)臺(tái)五軸數(shù)控機(jī)床結(jié)構(gòu)示意圖與運(yùn)動(dòng)鏈拓?fù)銯ig.1 Structure diagram and kinematic chain topology of double turntable five-axis CNC machine tool

1.2 幾何誤差描述

在幾何誤差綜合作用下，刀具相對(duì)于工件的實(shí)際位姿與理想位姿間存在偏差。根據(jù)各幾何誤差與機(jī)床運(yùn)動(dòng)軸位置的相關(guān)性，可將其分為位置無(wú)關(guān)幾何誤差 （Position independent geometric error，PIGEs）和位置相關(guān)幾何誤差 （Position dependent geometric errors，PDGEs）。目前在國(guó)際學(xué)術(shù)界，針對(duì)旋轉(zhuǎn)軸PIGEs有“絕對(duì)表示”和“相對(duì)表示”兩種定義方法。由于“高”旋轉(zhuǎn)軸的PIGEs會(huì)受“低”旋轉(zhuǎn)軸運(yùn)動(dòng)的影響，絕對(duì)表示法對(duì)于雙轉(zhuǎn)臺(tái)五軸機(jī)床并不適用。因此，本文采用相對(duì)表示法[8]，41項(xiàng)幾何誤差及其編號(hào)如表1所示。

表1 五軸數(shù)控機(jī)床幾何誤差及其編號(hào)Table 1 Geometric errors and their numbering of five-axis CNC machine tool

1.3 基于旋量理論的空間誤差模型

根據(jù)旋量理論[9]，剛體的位姿變換運(yùn)動(dòng)可以由圍繞空間某一直線的旋轉(zhuǎn)運(yùn)動(dòng)以及沿該直線的平移運(yùn)動(dòng)合成。其運(yùn)動(dòng)學(xué)模型的基本單元是旋量運(yùn)動(dòng)eξ·θ，由旋量ξ和運(yùn)動(dòng)量θ構(gòu)成。

對(duì)于旋轉(zhuǎn)軸，旋量ξ可以表示為

式中，v為平動(dòng)軸軸線方向的單位向量；ω為旋轉(zhuǎn)軸軸線方向的單位向量。v=–ω×q，q為軸線上的任意一點(diǎn)。

旋轉(zhuǎn)軸旋量運(yùn)動(dòng)的變換矩陣可以表示為

對(duì)于平動(dòng)軸，旋量ξ可以表示為

平動(dòng)軸旋量運(yùn)動(dòng)的變換矩陣可以表示為

式中，θ為平動(dòng)軸運(yùn)動(dòng)量。

結(jié)合圖1描述的運(yùn)動(dòng)鏈拓?fù)?，可得到刀具鏈位姿變換矩陣gbt（θ）和工件鏈位姿變換矩陣gbw（θ），即

根據(jù)式 （5）和 （6），可得到理想狀況下刀具相對(duì)于工件的位姿變換矩陣（θ），即

結(jié)合表1和式（7），可得到實(shí)際狀況下刀具相對(duì)于工件的位姿變換矩陣（θ），即

由式 （7）和 （8）可得到五軸數(shù)控機(jī)床的位姿誤差矩陣E，即

位姿誤差矩陣中包含位置誤差和姿態(tài)誤差，各個(gè)方向的位置誤差分量PX、PY、PZ為

各個(gè)方向的姿態(tài)誤差分量VX、VY、VZ為

2 基于PCE的全局靈敏度分析方法

2.1 多項(xiàng)式混沌展開模型

多項(xiàng)式混沌展開最早由Wiener[10]提出，用于研究Brown運(yùn)動(dòng)；Sudret[11]后續(xù)將PCE應(yīng)用到靈敏度分析中。其基本原理是將含有N維隨機(jī)輸入變量的響應(yīng)函數(shù)Y=M（x），用一組關(guān)于變量x的多項(xiàng)式混沌展開模型表示

式中，ψb（x）為正交的多項(xiàng)式基底，其最優(yōu)類型取決于變量的分布函數(shù)[12]；yb為多項(xiàng)式基底對(duì)應(yīng)的系數(shù)。

系數(shù)的求解是PCE方法的關(guān)鍵，一般利用最小二乘回歸求解。

式中，A為回歸矩陣，包含回歸樣本所對(duì)應(yīng)的多項(xiàng)式基底；y為多項(xiàng)式系數(shù)構(gòu)成的列向量。

在實(shí)際工程應(yīng)用中，通常會(huì)對(duì)式（13）進(jìn)行截?cái)唷．?dāng)最高階次為p階時(shí)，PCE模型可近似為

式中，P為截?cái)嘟票Ａ舻亩囗?xiàng)式總項(xiàng)數(shù)，可由式 （16）確定。

不難看出，PCE方法的計(jì)算量會(huì)隨著輸入變量維數(shù)的增多呈階乘增長(zhǎng)，高維下存在“維數(shù)災(zāi)難”問題。該問題可通過對(duì)PCE模型進(jìn)行稀疏化來解決，其思想是通過去除模型中對(duì)輸出響應(yīng)影響不大的正交多項(xiàng)式項(xiàng)來減少多項(xiàng)式系數(shù)的個(gè)數(shù)，從而降低計(jì)算量。根據(jù)這一思想，式 （14）可以轉(zhuǎn)換為求解多項(xiàng)式系數(shù)向量中非零元素最少的優(yōu)化問題。

式中，ε為允許容差。

直接求解式 （17）是一個(gè)NP–hard難題，為了以可接受的計(jì)算成本獲得近似的最稀疏解，可采取一系列重構(gòu)算法。其中，正交匹配追蹤（Orthogonal matching pursuit，OMP）具有計(jì)算成本低、收斂速度快、適合交叉驗(yàn)證等優(yōu)勢(shì)，本文使用該方法實(shí)現(xiàn)多項(xiàng)式系數(shù)的快速求解。

2.2 基于OMP的多項(xiàng)式系數(shù)求解

正交匹配追蹤由Pati等[13]提出，其工作原理是在A矩陣中檢索與當(dāng)前殘差最相關(guān)的多項(xiàng)式基底，并將其添加到已有的索引集中，通過最小二乘法更新相應(yīng)的多項(xiàng)式系數(shù)，不斷迭代，直至殘差降至指定的容差范圍內(nèi)，算法流程如圖2所示。

圖2 OMP算法流程圖Fig.2 Flowchart of OMP algorithm

由于在每一個(gè)迭代步，OMP都尋找能使殘差下降最快的多項(xiàng)式基底，因此能夠快速求解多項(xiàng)式系數(shù)。

2.3 基于PCE的全局靈敏度指數(shù)

基于方差分解的Sobol靈敏度指數(shù)[14]因簡(jiǎn)單有效而得到了廣泛應(yīng)用，可通過多項(xiàng)式系數(shù)直接得到。根據(jù)該方法，可將含有N維隨機(jī)輸入變量的響應(yīng)函數(shù)Y=M（x），唯一地分解成2N個(gè)子函數(shù)之和，即

式中，s=1，2，…，N；M0為0階常數(shù)項(xiàng)；Mi（xi）為與xi有關(guān)的1階子函數(shù)；為與和有關(guān)的2階子函數(shù)；為與，…，有關(guān)的s階子函數(shù)；M1，2，…，N（x1，…，xN）為與所有輸入變量有關(guān)的N階子函數(shù)。

對(duì)式 （18）等式兩邊取方差，可得

可以看出，子函數(shù)的方差表征了不同單個(gè)輸入及其相互作用對(duì)輸出響應(yīng)方差的貢獻(xiàn)。因此，變量，…，的Sobol靈敏度指數(shù)被定義為

再對(duì)截?cái)嗪蟮腜CE模型 （式（15））進(jìn)行改寫，可得

對(duì)比式 （18）和 （21）可知，Sobol分解式中變量，…，的s階子函數(shù)就是PCE模型中含相同變量的s階正交多項(xiàng)式之和，即

由此，結(jié)合式 （20）和 （22）就可得到基于PCE的Sobol靈敏度指數(shù)，即

3 五軸數(shù)控機(jī)床幾何誤差全局靈敏度分析實(shí)例

3.1 幾何誤差測(cè)量試驗(yàn)

為獲得幾何誤差的概率分布，基于四川普什寧江機(jī)床有限公司生產(chǎn)的VMC80IV型五軸數(shù)控機(jī)床搭建試驗(yàn)平臺(tái)，如圖3所示。選擇Renishaw XL–80激光干涉儀作為測(cè)量?jī)x器，運(yùn)用“15線法[15]”測(cè)量機(jī)床平動(dòng)軸幾何誤差；再利用Renishaw QC20–W球桿儀，依據(jù)Tsutsumi等[16]提出的策略測(cè)量機(jī)床旋轉(zhuǎn)軸幾何誤差，重復(fù)進(jìn)行多次試驗(yàn)。期間，測(cè)試環(huán)境保持 （20±0.5） ℃恒溫，以抑制機(jī)床熱致空間誤差對(duì)測(cè)試結(jié)果的影響。辨識(shí)并統(tǒng)計(jì)測(cè)量數(shù)據(jù)，可以得到各個(gè)誤差的近似概率分布。結(jié)果表明，41項(xiàng)幾何誤差服從正態(tài)分布，對(duì)誤差范圍進(jìn)行合理縮放后，取位置誤差范圍為±15 μm，角度誤差范圍為±10″。

圖3 幾何誤差測(cè)量試驗(yàn)平臺(tái)Fig.3 Experimental platform of geometric errors measurement

3.2 基于PCE的全局靈敏度分析

在靈敏度分析前，首先要確定幾何誤差采樣測(cè)試點(diǎn)。VMC80IV型機(jī)床X軸、Y軸和Z軸的行程分別為±400 mm、±425 mm和±275 mm，C軸的回轉(zhuǎn)范圍為0～360°，A軸的擺角范圍為–130°～+130°。將機(jī)床的工作空間沿對(duì)角線劃分，并在每條對(duì)角線上等距取5個(gè)點(diǎn)作為測(cè)試點(diǎn)，如圖4所示。這樣選取的測(cè)試點(diǎn)適用于全尺寸的被加工件，且接近實(shí)際工作空間附近，位于機(jī)床常用的坐標(biāo)范圍內(nèi)。越接近工作空間中心，越為常加工區(qū)域，測(cè)試點(diǎn)越密集，與實(shí)際加工的邏輯相符，在這些點(diǎn)進(jìn)行靈敏度分析將得到更有意義的結(jié)果。

圖4 機(jī)床空間測(cè)試點(diǎn)分布（mm）Fig.4 Spatial test points distribution of machine tool (mm)

借助UQLab[17]編程構(gòu)建幾何誤差的PCE模型。多項(xiàng)式基底選擇Hermite正交多項(xiàng)式，根據(jù)Cameron-Martin理論[18]，其對(duì)服從正態(tài)分布的輸入變量具有指數(shù)收斂速度。模型的最高階次和樣本量可通過收斂性測(cè)試確定，經(jīng)過測(cè)試，取階次p=3，樣本量m=1×103。采用OMP求解多項(xiàng)式系數(shù)，再通過Sobol指數(shù)表征幾何誤差對(duì)各方向位姿誤差分量的影響程度。由于幾何誤差相對(duì)較小，且誤差項(xiàng)之間具有低耦合性，一階靈敏度與總靈敏度結(jié)果相差不大。因此，本文僅討論總靈敏度指數(shù)，歸一化結(jié)果如圖5和6所示。

圖5 位置誤差全局靈敏度分析結(jié)果Fig.5 Global sensitivity analysis results of position error

將靈敏度指數(shù)大于0.05（大約是平均靈敏度指數(shù)的兩倍）的誤差項(xiàng)作為關(guān)鍵幾何誤差。根據(jù)圖5可知，影響X方向位置誤差的關(guān)鍵幾何誤差為δxc、δya、δyAM；對(duì)Y方向位置誤差影響最大的為δyc和δyCA，此外，δxx、δxy、δxz、δxa、δxAM也有較大影響；Z方向上影響最大的為δza、δzc、δzAM，此外，δyx、δyy、δyz、εxc的影響也不容忽視。從分析結(jié)果中不難發(fā)現(xiàn)，位置誤差的影響要明顯超過角度誤差，而旋轉(zhuǎn)軸的影響要超過平動(dòng)軸，PIGEs和PDGEs的影響差別不大。

根據(jù)圖6可知，影響X方向姿態(tài)誤差最關(guān)鍵的幾何誤差為εyx、εyy、εyz、Sxz，此外，εza、εzc、γAM、εyc、βCA也有較大影響；Y方向上的關(guān)鍵幾何誤差除εyc和βCA外，與X方向完全相同，且影響權(quán)重更高；對(duì)Z方向影響最大的幾何誤差為εxx、εxy、εxz、εxa、Syz、αAM，同時(shí)εyc和βCA也不容忽視。可以看出，對(duì)姿態(tài)誤差影響較大的全部都是角度誤差，旋轉(zhuǎn)軸和平動(dòng)軸的影響差別不大。

圖6 姿態(tài)誤差全局靈敏度分析結(jié)果Fig.6 Global sensitivity analysis results of attitude error

3.3 靈敏度分析方法對(duì)比與驗(yàn)證

為驗(yàn)證本文方法的正確性與高效性，與MCS法和LHS法進(jìn)行對(duì)比。這兩種方法是目前靈敏度分析的常用方法，且都需要大規(guī)模采樣，以樣本量為1×105的結(jié)果作為驗(yàn)證標(biāo)準(zhǔn)。3種方法所得的關(guān)鍵幾何誤差及其靈敏度指數(shù)如表2所示。不難看出，3種方法得到的關(guān)鍵誤差完全相同，且靈敏度指數(shù)接近，最大相對(duì)誤差不超過5%。因此，本文方法的正確性得到驗(yàn)證。

表2 關(guān)鍵幾何誤差靈敏度指數(shù)對(duì)比Table 2 Comparison of sensitivity index for key geometric errors

此外，對(duì)3種方法分別進(jìn)行收斂性能測(cè)試。以X、Y、Z方向上的最關(guān)鍵的幾何誤差δxc、δyc、δza為例，測(cè)試在不同采樣規(guī)模下的收斂過程，并用波動(dòng)率R來量化靈敏度指數(shù)的波動(dòng)狀況，即

式中，l為測(cè)試樣本量。

以R≤3%作為收斂標(biāo)準(zhǔn)。根據(jù)圖7所示結(jié)果，誤差δxc在PCE法采樣規(guī)模為1×103時(shí)達(dá)到收斂標(biāo)準(zhǔn)；而MCS法和LHS法在采樣規(guī)模為1×103時(shí)波動(dòng)率仍在10%以上，在1×105時(shí)才達(dá)到收斂標(biāo)準(zhǔn)。誤差δyc和δza測(cè)試結(jié)果也相同，這意味著本文方法具有更高的收斂性能，且在樣本量減小兩個(gè)數(shù)量級(jí)的情況下達(dá)到相同的計(jì)算精度。

表3還對(duì)比了3種方法的計(jì)算時(shí)間。結(jié)果顯示，本文方法在不降低計(jì)算精度的前提下，計(jì)算時(shí)間相對(duì)MCS法減少了96.8%，相對(duì)LHS法減少了98.1%。此外，本文方法中大部分時(shí)間用于PCE建模，只要模型構(gòu)建完成，便可快速求解靈敏度指數(shù)。隨著機(jī)床日益復(fù)雜，面對(duì)更多軸類型的機(jī)床和更多的幾何誤差，本文方法的優(yōu)勢(shì)將更加顯著。

表3 計(jì)算時(shí)間對(duì)比Table 3 Comparison of calculation time

4 結(jié)論

（1）本文提出了基于多項(xiàng)式混沌展開的全局靈敏度分析方法，運(yùn)用正交匹配追蹤算法實(shí)現(xiàn)多項(xiàng)式系數(shù)的快速求解，并給出了基于該方法的Sobol靈敏度指數(shù)。

（2）測(cè)量得到機(jī)床41項(xiàng)幾何誤差的近似概率分布，在確定采樣測(cè)試點(diǎn)后，使用PCE法對(duì)幾何誤差進(jìn)行全局靈敏度分析，得到了影響各方向位姿誤差分量的關(guān)鍵幾何誤差。

（3）與MCS法和LHS法進(jìn)行對(duì)比，本文方法的正確性得到驗(yàn)證，且在保證計(jì)算精度的前提下，具有更高的收斂性能，樣本量從1×105降低到1×103，計(jì)算時(shí)間分別減少了96.8%和98.1%。