基于最小區域法的形位公差精確算法

計自飛，丁 拳

（中航工業沈陽發動機設計研究所，沈陽110015）

0 引言

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

本文提出1種基于最小區域法的逐次逼近算法，思路簡單，易于編程實現，又很好地克服了最小二乘法的諸多不足，通過編寫成fortran程序，編譯生成的可執行文件可在任意1臺計算機上運行，便于推廣應用。

1 最小二乘法的求解過程

最小二乘法求解形位公差的理論發展比較完善，算法簡單，因此，得到了廣泛應用。根據文獻[1]中的推導，其數學模型是非線性的，但如果在測點均勻分布時，計算模型就能轉化為線性化模型，而且，根據GB/T1182-2008中對于形位公差的定義，最小二乘法的求解思路不滿足最小條件要求，其計算結果與真實形位公差存在差異，這屬于評定方法的系統誤差。下面簡要介紹利用最小二乘法求解形位公差的方法。

1.1 直線度

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

設（xi,yi）,i=1,2,…，n 為給定平面內的測點坐標。設由最小二乘法求出的理想直線方程為

式中：k、b 為待定參數。

則殘余誤差為

最小二乘法的目標函數為

約束條件為

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

1.2 平面度

設（xi，yi，zi），i=1，2，…，n 為給定的空間測點坐標。設由以上測點擬合出的最小二乘平面方程為

式中：A、B、C 為待定參數。

則殘余誤差為

最小二乘法的目標函數為

約束條件為

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

1.3 圓度

設（xi,yi,zi）,i=1,2,…，n 為被測圓周上采樣點坐標，由以上測點擬合出的最小二乘圓方程為

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

則殘余誤差為

最小二乘法的目標函數為

約束條件為

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

式中：di為各測點到最小二乘圓心的距離。

1.4 圓柱度

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

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

則殘余誤差為

最小二乘法的目標函數為

約束條件為

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

式中：di為各測點到最小二乘圓柱軸線的距離。

2 基于MATLAB的最小區域法優化的實現方法

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

圖1 MATLAB優化工具箱GUI界面

基于最小區域法的形位公差評定問題屬無約束條件的非線性極小值問題，可用fminunc或fminsearch函數實現[10]，以fminunc函數為例：

[x,fval,exitflag,output]=fminun(fun,x0,options)其中：x 為最優解；fval 為返回目標函數在最優解處的值；exitflag為返回算法終止標志；output為返回優化算法信息的1個數據結構；fun 為目標函數；x0為初始值；options 為設置優化選項參數。

只要將式（5）、（10）、（15）、（21）作為目標函數，即可利用上述方法完成基于最小區域法的形位公差的評定。

3 逐次逼近算法的求解過程

3.1 算法的輸入條件

逐次逼近算法，不像區域搜索法那樣，要靠最小二乘法的計算結果來確定搜索范圍，只需要給出測點數目、采用的測量坐標系、各個測點的坐標即可，將其以固定格式寫入txt文件即可，程序運行時可自動獲取這些信息。

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

其中

3.2 算法的求解過程

逐步逼近算法的具體實施步驟：

（1）按一定格式編輯測點信息，待程序運行后，根據測點信息確定了待定參數的初次計算域，并將每個優化參數的計算域n 等分。

（2）分別計算當優化參數位于各等分點上時的形位公差值，找出形位公差取最小值時對應的優化參數值。

（3）以上次計算找出的優化參數為中心，重新劃分1個較小的計算域，再將每個優化參數的計算域n等分。

（4）再次計算當優化參數分別位于各自等分點上時的形位公差值，找出形位公差取最小值時對應的優化參數值。

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

將上述逐次逼近算法用于回轉零部件的圓度計算，過程如圖2所示。（1）根據測點坐標分布范圍確定初次計算域；（2）以第1次找出的圓度最小值對應的坐標為中心，重新設定1個較小的計算域；（3）循環n次后，再次縮小計算域后得到的圓度值和圓心位置基本不變。

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

圖2 利用逐次逼近算法的圓度計算過程

圖3 逐次逼近算法流程

4 算法驗證

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

4.1 直線度

1組等間距分布測點的直角坐標值如下：

分別利用3種方法計算得出的直線度誤差以及相應的最優解斜率詳見表1。

表1 直線度比較

4.2 圓度

對于下面1組周向均布的極坐標下的測點數據[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種方法計算圓度和位置度偏差，見表2。

表2 圓度比較

4.3 分析

通過對比3種不同方法對直線度和圓度誤差的評定結果，可見最小二乘法的評定結果不滿足最小條件要求。除圓心坐標外，逐次逼近算法的評定結果與MATLAB優化結果一致，這是因為最小區域法求出的圓心不惟一，但圓度誤差值惟一。MATLAB作為成熟的商業軟件，其優化結果相對穩定，逐次逼近算法的計算結果與之一致，證明了該算法的有效性。

5 總結

（1）逐次逼近算法思路簡單，易于編程實現，初始計算條件簡單，并克服了最小二乘法的諸多不足，且無需特殊運行環境。

（2）通過與其他方法計算結果對比說明，逐次逼近算法可以精確高效地完成形位公差的評定，計算結果好于最小二乘法，符合最小條件要求。

（3）逐次逼近算法簡單易行，計算正確有效，值得在工程實踐中應用和推廣。

[1]黃富貴，鄭育軍.基于區域搜索的圓度誤差評定方法[J].計量學報，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]全國產品尺寸和幾何技術規范標準化技術委員會.GB/T 1182-2008產品幾何技術規范(GPS)幾何公差形狀、方向、位置和跳動公差標注[S].北京：中國標準出版社，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]岳武陵，吳勇.基于仿增量算法的圓度誤差快速準確評定[J].機械工程學報，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]田樹耀，黃富貴，張彬.一種基于區域搜索的平面度誤差評定方法 [J].華僑大學學報：自然科學版，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]張小萍，王君澤.基于測點分類的圓度誤差精確評定[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]劉國光.基于MATLAB評定圓柱度誤差[J].工程設計學報，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的圓度評定方法 [J].計量與測試技術，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]全國產品尺寸和幾何技術規范標準化技術委員會.GB/T 7235-2004產品幾何技術規范(GPS)評定圓度誤差的方法半徑變化量測量[S].北京：中國標準出版社，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]劉國光.基于MATLAB評定圓柱度誤差[J].工程設計學報，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的圓度評定方法[J].計量與測試技術，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最優化計算[M].北京：清華大學出版社，2009：169-171.GONG Chun，WANG Zhenglin.Proficient in MATLAB optimization calculation[M].Beijing:Tsinghua University Press，2009：169-171.(in Chinese)

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

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