計(jì)自飛，丁 拳

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

0 引言

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

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

1 最小二乘法的求解過程

最小二乘法求解形位公差的理論發(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ì)于形位公差的定義，最小二乘法的求解思路不滿足最小條件要求，其計(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)化算法，一旦初值合理，收斂速度快，滿足約束的精度高。合理利用優(yōu)化工具箱可以很方便地求解一些有約束或無約束的非線性極小值問題，優(yōu)化工具箱的GUI界面如圖1所示。

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

基于最小區(qū)域法的形位公差評(píng)定問題屬無約束條件的非線性極小值問題，可用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 逐次逼近算法的求解過程

3.1 算法的輸入條件

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

通過一些簡(jiǎn)單的語句可以確定測(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 算法的求解過程

逐步逼近算法的具體實(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ì)算，過程如圖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ì)算過程

圖3 逐次逼近算法流程

4 算法驗(yàn)證

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

4.1 直線度

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

分別利用3種方法計(jì)算得出的直線度誤差以及相應(yīng)的最優(yōu)解斜率詳見表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ì)算圓度和位置度偏差，見表2。

表2 圓度比較

4.3 分析

通過對(duì)比3種不同方法對(duì)直線度和圓度誤差的評(píng)定結(jié)果，可見最小二乘法的評(píng)定結(jié)果不滿足最小條件要求。除圓心坐標(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)單，并克服了最小二乘法的諸多不足，且無需特殊運(yùn)行環(huán)境。

（2）通過與其他方法計(jì)算結(jié)果對(duì)比說明，逐次逼近算法可以精確高效地完成形位公差的評(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]田樹耀，黃富貴，張彬.一種基于區(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)分類的圓度誤差精確評(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)