基于最小區(qū)域法的形位公差精確算法

計(jì)自飛，丁 拳

（中航工業(yè)沈陽(yáng)發(fā)動(dòng)機(jī)設(shè)計(jì)研究所，沈陽(yáng)110015）

0 引言

航空發(fā)動(dòng)機(jī)是由成千上萬(wàn)的各種零部件組成，因此，其高性能源于優(yōu)秀設(shè)計(jì)，更離不開(kāi)優(yōu)秀制造工藝的可靠保證。零部件的精密程度將直接影響到發(fā)動(dòng)機(jī)的穩(wěn)定性和可靠性。在發(fā)動(dòng)機(jī)的各式零部件中，回轉(zhuǎn)體構(gòu)件數(shù)量眾多，且具有特殊的重要作用。因此，對(duì)發(fā)動(dòng)機(jī)回轉(zhuǎn)體構(gòu)件的圓度、圓柱度和直線度等的測(cè)量分析十分重要。而其準(zhǔn)確的測(cè)量結(jié)果和快速評(píng)估對(duì)于判定精密回轉(zhuǎn)類(lèi)零部件的合格與否至關(guān)重要。形位公差的評(píng)定方法有很多種，其中，最小區(qū)域法和最小二乘法適用于每種形位公差的計(jì)算。最小二乘法盡管發(fā)展較為完善，但是，方法本身存在一些缺陷，而最小區(qū)域法是完全符合形位公差定義的評(píng)定方法[1-2]。對(duì)此，國(guó)內(nèi)外的專(zhuān)家學(xué)者已做過(guò)很多研究，并提出了一些優(yōu)化方法，但這些方法有些思路復(fù)雜，如遺傳算法[3]、仿增量算法[4]等，有些方法對(duì)初始條件要求過(guò)高，如區(qū)域搜索算法[1-5]等，需以最小二乘法的計(jì)算結(jié)果作為初始條件，且搜索步長(zhǎng)難于控制，過(guò)小則計(jì)算量大，過(guò)大則精度不夠，因而都不易于在工程領(lǐng)域廣泛采用。有些研究者利用MATLAB自帶的優(yōu)化工具箱來(lái)實(shí)現(xiàn)最小區(qū)域法求解圓度誤差的問(wèn)題[6-8]，但使用這種方法需要精通MATLAB軟件，推廣應(yīng)用有諸多不便。

本文提出1種基于最小區(qū)域法的逐次逼近算法，思路簡(jiǎn)單，易于編程實(shí)現(xiàn)，又很好地克服了最小二乘法的諸多不足，通過(guò)編寫(xiě)成fortran程序，編譯生成的可執(zhí)行文件可在任意1臺(tái)計(jì)算機(jī)上運(yùn)行，便于推廣應(yīng)用。

1 最小二乘法的求解過(guò)程

最小二乘法求解形位公差的理論發(fā)展比較完善，算法簡(jiǎn)單，因此，得到了廣泛應(yīng)用。根據(jù)文獻(xiàn)[1]中的推導(dǎo)，其數(shù)學(xué)模型是非線性的，但如果在測(cè)點(diǎn)均勻分布時(shí)，計(jì)算模型就能轉(zhuǎn)化為線性化模型，而且，根據(jù)GB/T1182-2008中對(duì)于形位公差的定義，最小二乘法的求解思路不滿(mǎn)足最小條件要求，其計(jì)算結(jié)果與真實(shí)形位公差存在差異，這屬于評(píng)定方法的系統(tǒng)誤差。下面簡(jiǎn)要介紹利用最小二乘法求解形位公差的方法。

1.1 直線度

直線度包括給定方向上的直線度、給定平面內(nèi)的直線度、空間任意方向上的直線度等3種情況。這里研究給定平面內(nèi)的直線度。

設(shè)（xi,yi）,i=1,2,…，n 為給定平面內(nèi)的測(cè)點(diǎn)坐標(biāo)。設(shè)由最小二乘法求出的理想直線方程為

式中：k、b 為待定參數(shù)。

則殘余誤差為

最小二乘法的目標(biāo)函數(shù)為

約束條件為

將式（3）帶入式（4），化簡(jiǎn)后求解，得到k、b 的表達(dá)式。則直線度誤差為

1.2 平面度

設(shè)（xi，yi，zi），i=1，2，…，n 為給定的空間測(cè)點(diǎn)坐標(biāo)。設(shè)由以上測(cè)點(diǎn)擬合出的最小二乘平面方程為

式中：A、B、C 為待定參數(shù)。

則殘余誤差為

最小二乘法的目標(biāo)函數(shù)為

約束條件為

將式（8）帶入式（9），化簡(jiǎn)后可以得到A、B、C 的表達(dá)式，從而平面度誤差的最小二乘法結(jié)果為

1.3 圓度

設(shè)（xi,yi,zi）,i=1,2,…，n 為被測(cè)圓周上采樣點(diǎn)坐標(biāo)，由以上測(cè)點(diǎn)擬合出的最小二乘圓方程為

式中：（x0,y0）為最小二乘圓心；R 為最小二乘圓半徑。

則殘余誤差為

最小二乘法的目標(biāo)函數(shù)為

約束條件為

將式（13）帶入式（14），化簡(jiǎn)后可以得到A、B、C的表達(dá)式，從而平面度誤差的最小二乘法結(jié)果為

式中：di為各測(cè)點(diǎn)到最小二乘圓心的距離。

1.4 圓柱度

假設(shè)各橫向測(cè)量截面上測(cè)點(diǎn)坐標(biāo)為（xi,yi,zi）,i=1,2,…，n，設(shè)由以上測(cè)點(diǎn)擬合出的最小二乘圓柱的半徑為R，則軸線方程為

式中：x0、y0、z0、l、m、n、R 為待定參數(shù)。

則殘余誤差為

最小二乘法的目標(biāo)函數(shù)為

約束條件為

將式（19）帶入式（20），化簡(jiǎn)后可以得到A、B、C的表達(dá)式，從而平面度誤差的最小二乘法結(jié)果為

式中：di為各測(cè)點(diǎn)到最小二乘圓柱軸線的距離。

2 基于MATLAB的最小區(qū)域法優(yōu)化的實(shí)現(xiàn)方法

MATLAB軟件自帶的優(yōu)化工具箱(optimization toolbox)提供了諸如共軛梯度法、單純法、牛頓法、序列二次規(guī)劃等局部?jī)?yōu)化算法，一旦初值合理，收斂速度快，滿(mǎn)足約束的精度高。合理利用優(yōu)化工具箱可以很方便地求解一些有約束或無(wú)約束的非線性極小值問(wèn)題，優(yōu)化工具箱的GUI界面如圖1所示。

圖1 MATLAB優(yōu)化工具箱GUI界面

基于最小區(qū)域法的形位公差評(píng)定問(wèn)題屬無(wú)約束條件的非線性極小值問(wèn)題，可用fminunc或fminsearch函數(shù)實(shí)現(xiàn)[10]，以fminunc函數(shù)為例：

[x,fval,exitflag,output]=fminun(fun,x0,options)其中：x 為最優(yōu)解；fval 為返回目標(biāo)函數(shù)在最優(yōu)解處的值；exitflag為返回算法終止標(biāo)志；output為返回優(yōu)化算法信息的1個(gè)數(shù)據(jù)結(jié)構(gòu)；fun 為目標(biāo)函數(shù)；x0為初始值；options 為設(shè)置優(yōu)化選項(xiàng)參數(shù)。

只要將式（5）、（10）、（15）、（21）作為目標(biāo)函數(shù)，即可利用上述方法完成基于最小區(qū)域法的形位公差的評(píng)定。

3 逐次逼近算法的求解過(guò)程

3.1 算法的輸入條件

逐次逼近算法，不像區(qū)域搜索法那樣，要靠最小二乘法的計(jì)算結(jié)果來(lái)確定搜索范圍，只需要給出測(cè)點(diǎn)數(shù)目、采用的測(cè)量坐標(biāo)系、各個(gè)測(cè)點(diǎn)的坐標(biāo)即可，將其以固定格式寫(xiě)入txt文件即可，程序運(yùn)行時(shí)可自動(dòng)獲取這些信息。

通過(guò)一些簡(jiǎn)單的語(yǔ)句可以確定測(cè)點(diǎn)坐標(biāo)的大致分布范圍，以此確定初次計(jì)算域。以圓度計(jì)算為例，設(shè)xmin，xmax分別為各測(cè)點(diǎn)x 坐標(biāo)的最小、最大值，ymin，ymax分別為各測(cè)點(diǎn)y 坐標(biāo)的最小、最大值，則首次計(jì)算域可選為以(xmin，ymin)、(xmin，ymax)、(xmax，ymin)、(xmax，ymax)為頂點(diǎn)的矩形區(qū)域。然后將長(zhǎng)和寬分別m、n 等分，則計(jì)算域內(nèi)第i 行、j 列的節(jié)點(diǎn)坐標(biāo)為

其中

3.2 算法的求解過(guò)程

逐步逼近算法的具體實(shí)施步驟：

（1）按一定格式編輯測(cè)點(diǎn)信息，待程序運(yùn)行后，根據(jù)測(cè)點(diǎn)信息確定了待定參數(shù)的初次計(jì)算域，并將每個(gè)優(yōu)化參數(shù)的計(jì)算域n 等分。

（2）分別計(jì)算當(dāng)優(yōu)化參數(shù)位于各等分點(diǎn)上時(shí)的形位公差值，找出形位公差取最小值時(shí)對(duì)應(yīng)的優(yōu)化參數(shù)值。

（3）以上次計(jì)算找出的優(yōu)化參數(shù)為中心，重新劃分1個(gè)較小的計(jì)算域，再將每個(gè)優(yōu)化參數(shù)的計(jì)算域n等分。

（4）再次計(jì)算當(dāng)優(yōu)化參數(shù)分別位于各自等分點(diǎn)上時(shí)的形位公差值，找出形位公差取最小值時(shí)對(duì)應(yīng)的優(yōu)化參數(shù)值。

（5）重復(fù)式（3）～（4），直到計(jì)算出的形位公差值與前1次的結(jié)果之差小于預(yù)定的偏差閾值ε。

將上述逐次逼近算法用于回轉(zhuǎn)零部件的圓度計(jì)算，過(guò)程如圖2所示。（1）根據(jù)測(cè)點(diǎn)坐標(biāo)分布范圍確定初次計(jì)算域；（2）以第1次找出的圓度最小值對(duì)應(yīng)的坐標(biāo)為中心，重新設(shè)定1個(gè)較小的計(jì)算域；（3）循環(huán)n次后，再次縮小計(jì)算域后得到的圓度值和圓心位置基本不變。

如上所述，假定每計(jì)算1次后就將每個(gè)待定參數(shù)的計(jì)算域減小一半，但等分的份數(shù)不變，這樣每計(jì)算1次，計(jì)算步長(zhǎng)就縮小1倍，因此，只要循環(huán)的次數(shù)足夠多，就一定能得到準(zhǔn)確的形位公差值。為了計(jì)算的可靠性，后1次循環(huán)的計(jì)算域比前1次的計(jì)算域不能縮小得太多。其算法流程如圖3所示。

圖2 利用逐次逼近算法的圓度計(jì)算過(guò)程

圖3 逐次逼近算法流程

4 算法驗(yàn)證

選取直線度和圓度2種形位公差，分別選用最小二乘法、MATLAB優(yōu)化算法以及逐次逼近算法進(jìn)行求解，通過(guò)對(duì)比來(lái)驗(yàn)證逐次逼近算法的有效性。

4.1 直線度

1組等間距分布測(cè)點(diǎn)的直角坐標(biāo)值如下：

分別利用3種方法計(jì)算得出的直線度誤差以及相應(yīng)的最優(yōu)解斜率詳見(jiàn)表1。

表1 直線度比較

4.2 圓度

對(duì)于下面1組周向均布的極坐標(biāo)下的測(cè)點(diǎn)數(shù)據(jù)[41.0532,41.0527,41.0536,41.0072,41.0065,40.9982,41.0526,41.0538,40.9973,41.0026,40.9957,41.0546,41.0549,40.9995,40.9969,41.0439]，分別利用3種方法計(jì)算圓度和位置度偏差，見(jiàn)表2。

表2 圓度比較

4.3 分析

通過(guò)對(duì)比3種不同方法對(duì)直線度和圓度誤差的評(píng)定結(jié)果，可見(jiàn)最小二乘法的評(píng)定結(jié)果不滿(mǎn)足最小條件要求。除圓心坐標(biāo)外，逐次逼近算法的評(píng)定結(jié)果與MATLAB優(yōu)化結(jié)果一致，這是因?yàn)樽钚^(qū)域法求出的圓心不惟一，但圓度誤差值惟一。MATLAB作為成熟的商業(yè)軟件，其優(yōu)化結(jié)果相對(duì)穩(wěn)定，逐次逼近算法的計(jì)算結(jié)果與之一致，證明了該算法的有效性。

5 總結(jié)

（1）逐次逼近算法思路簡(jiǎn)單，易于編程實(shí)現(xiàn)，初始計(jì)算條件簡(jiǎn)單，并克服了最小二乘法的諸多不足，且無(wú)需特殊運(yùn)行環(huán)境。

（2）通過(guò)與其他方法計(jì)算結(jié)果對(duì)比說(shuō)明，逐次逼近算法可以精確高效地完成形位公差的評(píng)定，計(jì)算結(jié)果好于最小二乘法，符合最小條件要求。

（3）逐次逼近算法簡(jiǎn)單易行，計(jì)算正確有效，值得在工程實(shí)踐中應(yīng)用和推廣。

[1]黃富貴，鄭育軍.基于區(qū)域搜索的圓度誤差評(píng)定方法[J].計(jì)量學(xué)報(bào)，2008，29（2）：117-119.HUANG Fugui，ZHENG Yujun.A method for roundness error evaluation based on area hunting[J].Acta Metrologica Sinica，2008，29（2）：117-119.（in Chinese）

[2]全國(guó)產(chǎn)品尺寸和幾何技術(shù)規(guī)范標(biāo)準(zhǔn)化技術(shù)委員會(huì).GB/T 1182-2008產(chǎn)品幾何技術(shù)規(guī)范(GPS)幾何公差形狀、方向、位置和跳動(dòng)公差標(biāo)注[S].北京：中國(guó)標(biāo)準(zhǔn)出版社，2008：1-48.National Technical Committee on Product Dimension and Geometric Element Specifications of Standardization Administration of China.GB/T 1182-2008 Geometrical Product Specifications(GPS)—geometrical tolerancing—tolerances of form，orientation，location and run-out[S].Beijing：China Standard Press，2008：1-48.（in Chinese）

[3]Cui Changcai，Che Rensheng，Ye Dong.Circularity error or evaluation using genetic algorithm [J].Optics and Pricision Engineering，2001，9（6）：499-505.

[4]岳武陵，吳勇.基于仿增量算法的圓度誤差快速準(zhǔn)確評(píng)定[J].機(jī)械工程學(xué)報(bào)，2008，44（1）：87-91.YUE Wuling，WU Yong.Fast and accurate evaluation of the roundness error based on quasi-incremental algorithm[J].Chinese Journal of Mechanical of Engineering，2008，44（1）：87-91.（in Chinese）

[5]田樹(shù)耀，黃富貴，張彬.一種基于區(qū)域搜索的平面度誤差評(píng)定方法 [J].華僑大學(xué)學(xué)報(bào)：自然科學(xué)版，2009，30（5）：506-508.TIAN Shuyao，HUANG Fugui，ZHANG Bin.An evaluation method for flatness error based on region searching[J].Journal of Huaqiao University：Natural Science，2009，30（5）：506-508.（in Chinese）

[6]張小萍，王君澤.基于測(cè)點(diǎn)分類(lèi)的圓度誤差精確評(píng)定[J].機(jī)械設(shè)計(jì)與制造，2010（9）：166-167.ZHANG Xiaoping，WANG Junze.Accurate evaluation of roundness error by classifying measured points[J].Machinery Design&Manufacture，2010（9）：166-167.（in Chinese）

[7]劉國(guó)光.基于MATLAB評(píng)定圓柱度誤差[J].工程設(shè)計(jì)學(xué)報(bào)，2005，12（4）：236-239.LIU Guoguang.Cylindricity error evaluation based on MATLAB [J].Journal of Engineering Design，2005，12（4）：236-239.（in Chinese）

[8]周建平.基于MATLAB的圓度評(píng)定方法 [J].計(jì)量與測(cè)試技術(shù)，2005，32（2）：5-7.ZHOU Jianping.The method of roundness assessment based on MATLAB[J].Metrology&Measurement Technique，2005，32（2）：5-7.(in Chinese)

[9]全國(guó)產(chǎn)品尺寸和幾何技術(shù)規(guī)范標(biāo)準(zhǔn)化技術(shù)委員會(huì).GB/T 7235-2004產(chǎn)品幾何技術(shù)規(guī)范(GPS)評(píng)定圓度誤差的方法半徑變化量測(cè)量[S].北京：中國(guó)標(biāo)準(zhǔn)出版社，2005：1-18.National Technical Committee on Product Dimension and Geometric Element Specifications of Standardization Administration of China.GB/T 7235-2004，Geometrical Product Specifications(GPS)—methods for the assessment of departure from roundness—measurement of variations in radius[S].Beijing：China Standard Press，2005：1-18.(in Chinese)

[10]劉國(guó)光.基于MATLAB評(píng)定圓柱度誤差[J].工程設(shè)計(jì)學(xué)報(bào)，2005，12（4）：236-239.LIU Guoguang.Cylindricity error evaluation based on MATLAB [J].Journal of Engineering Design，2005，12（4）：236-239.(in Chinese)

[11]周建平.基于MATLAB的圓度評(píng)定方法[J].計(jì)量與測(cè)試技術(shù)，2005，32（2）：5-7.ZHOU Jianping.The method of toundness assessment based on MATLAB [J].Metrology& Measurement Technique，2005，32（2）：5-7.(in Chinese)

[12]龔純，王正林.精通MATLAB最優(yōu)化計(jì)算[M].北京：清華大學(xué)出版社，2009：169-171.GONG Chun，WANG Zhenglin.Proficient in MATLAB optimization calculation[M].Beijing:Tsinghua University Press，2009：169-171.(in Chinese)

[13]欒艷華.某型航空發(fā)動(dòng)機(jī)機(jī)匣同軸度測(cè)量[J].航空發(fā)動(dòng)機(jī)，2012，38（1）：44-46.LUAN Yanhua.Casing coaxiality measurement for an aeroengine[J].Aeroengine，2012，38（1）：44-46.(in Chinese)

[14]謝海峰.圓度誤差測(cè)量中的多尺度算法研究[D].長(zhǎng)春：吉林大學(xué)，2009.XIE Haifeng.Research on the multi-scale slgorithm used in measurement of roundness error[D].Changchun：Jilin University，2009.(in Chinese)