基于最小條件的平面擬合高精算法及編程

林 翔,江速勇,陳柏良,陳玉霞

(福建商業高等專科學校,福建 福州 350012)

1 基準平面

目前在幾何量精密檢測領域中,高精三坐標測量儀得到廣泛的應用,如對圓度、平面度、直線度等各種形位誤差的檢測.在三坐標測量儀的支持下,不少種類的測量與誤差評定技術水平,已經達到三維空間的高精度計算,只有少數如圓度誤差評定的研究,一直停留在二維平面程度,只能處理二維坐標測量數據,與三坐標測量的實際脫節,不能真正滿足實際工程中誤差評定的需要,因此有必要對三維空間圓度測量數據的誤差評定作深入的實質性研究.

三維空間圓度誤差評定可以分為基準平面擬合與平面圓度計算兩大環節,基準平面的擬合的精度,直接影響圓度評定結果的精度.國標《GB/T1182-2004》[1]規定,只有在以"最小區域法"求取平面度誤差過程獲得的擬合平面,才符合"最小條件"原則.

平面度屬形狀誤差范疇,按國標[1]《GB/T 1182-2004》的表述,它是指被測物體表面相對于理想平面的變動量,此理想平面即所求的擬合基準平面.僅就平面度而言,通過"最小二乘法"求得的平面即可作為基準,以此基準即可計算出平面度誤差,但它不符合"最小條件".

圍繞這一判定準則,許多專家提出了不同算法,諸如"最小包容區域法"[3]、"有序判別法"[4]、"生物遺傳算法"[5]等等,這些算法各有特色,也收到不錯的計算結果,但由于這些近似算法的先天不足,雖能計算精度終歸偏低.分析評定準則,其要點是求取兩個平行的平面,把被測表面夾在其間,使兩平面之距達到最小,因此可以在"最小二乘

圖1 基準平面與被測平面示意圖

2 平面擬合

設平面上測量點集P={pk(xk,yk,zk),(k=1~n)},按照上述的解題大方向,可以選用"最小二乘平面"法,求取初始的基準擬合平面.設方程為π:ax+by+cz+d=0,如圖1,點集P對于平面π的平面度誤差值記若δ.此δ為統計意義下的平面度誤差值,是國標《GB/T 1182-2004》[1]認可的,在評定精度要求不高的情況下已具有實用意義,如汽車零部件的檢測.要達到高精度評定,還需研發更高精的算法.

經過對P與π之間距離最遠點所表現出來的特性分析發現,提高平面度誤差計算精度關鍵,是設法對π作細微轉動,即有意識地改變平面法矢T(a,b,c),就有可能降低δ值.方法如下:

由"最小二乘"法得到的初始擬合平面為π,點集P={Pk,(k=1~n)}中到π距離之最大者、最小者分別為Pi、Pj,距離分別為δi、δj;記Ai為Pi在π的投影點,Aj為Pj的投影點,δ=δi+δj,此即初始平面度誤差值.

再設ε是一個甚小值,沿AiPi方向距Ai點ε處取Ai'點,沿AjPj方向距Aj點ε處取Aj'點;記Q為Ai'Aj'中點(Q在π上),在π上再取一點Q',使Q'Q⊥Ai'Aj'.以Ai'、Aj'、Q'三點構造一個新的平面π',顯見π與π'之夾角很小,可將π'看作π繞Q'Q作細微轉動而獲得;從方程來看,就是π的法矢T(a,b,c)發生細微改變.

以 Pi、Ai'、Pj、Aj'4點構造平面 π1,再以π1對平面π垂直剖切,剖面如圖2所示,PiAi=δi,PiAi'=δi',易見,在ΔPiAi'Ai中,PiAi'

δ能否降下來,須通過點集P對新平面π'求平面度誤差值δ'與δ作比較并加以判斷:

1、如果降下來了,π'取代π,再重復上述操作對π作細微轉動;

2、如果沒有降,則將ε值減半,重新構造π',計算δ'并與δ比較,觀其下降與否;

3、如果ε經多次減半縮小成為非常小的值,轉動后得到的π'再也不能令δ'下降,且δ'與δ之差極小,達到了計算精度要求,則停止計算,π即所求的符合"最小條件"的擬合平面.

圖2 Pi、Pj等各點與平面π、π'關系圖

3 算法收斂于"最小區域"

圖2中,P距π的兩個最遠點為Pi、Pj,Pi在平面π的投影點為Ai,在平面π'的投影點為Ai',在△PiAiAi'中∠PiAi'Ai是鈍角,由于鈍角的任一鄰邊長度必然小于對角邊,即PiAi'

δ'=δi'+δj'=PiAi'+PjAj'

算法設定的終止條件為:當ε經多次縮半而變得甚小,π即使進行甚小的轉動,δ也不能下降,且δ與δ'之間相差達到計算精度要求,則停止計算.因為此時δ已非常接近谷底,非常接近平面度誤差真值了,且符合國標《BB/T 1182-2004》[1]中"最小區域法"平面度誤差的判定條件,所擬合的π即為符合"最小條件"的基準平面.

4 算法的實現

以"最小二乘法"平面為初始基準平面,既保證算法在收斂區間內進行,又有利于算法的快速收斂.

(1)"最小二乘法"求平面過程簡述:

平面點集P={Pk,(k=1~n)},初始擬合平面π:ax+by+cz+d=0,為討論方便,令T(a,b,c)為單位矢量,則任一點Pk到π的距離為:

"最小二乘"擬合平面的目標方程為:

要使 D 達到最小,只要令:?D/?a=0,?D/?b=0,?D/?c=0,?D/?d=0,及 a2+b2+c2=1即可;

聯解以上5等式,就能求得初始平面π的法矢T(a,b,c)和d各值,平面π的方程也就確定;求出δi、δj值,δ=δi+δj即為P在"最小二乘法"意義下的平面度誤差值.

(2)在初始擬合平面π的基礎上實施本算法:

由上已經獲取Pi、Pj二點,以及對應的δi、δj二值,設Pi、Pj在平面π的的投影點分別為Ai、Aj;連線Ai、Pi及Aj、Pj,沿AiPi方向距Ai點ε處取點Ai',沿AjPj方向距Aj點ε處取點Aj';求出Ai'、Aj'的中點Q,在π上求得點Q',使直線Q'Q⊥Ai'Aj';

以Ai'、Aj'、Q'三點構造平面π',求出P到π'的距離之兩個極值,記作δl'、δm',δ'=δl'+δm',比較δ'、δ并作如下判斷:

1)如果δ'<δ,平面π'便取代平面π,重新如法炮制求取新的平面π',并進行同樣的計算、判斷;

2)如果δ'≥δ,分兩種情況處理:

a、將ε值減半,重構平面π',算出δ',返回到步驟1)重新作判斷;

b、如果ε已經縮小為非常小的值,δ與δ'相差達到精度要求,滿足計算的終止條件,停止計算,π即為滿足"最小條件"擬合基準平面.

(3)特殊情況的處理

在上述計算過程中,如果Pk(k=1~n)到平面π的正最遠點、負最遠點不是唯一的,即"Pi"、"Pj"點不止一個,作如下處理:設有U個"Pi"、V個"Pj"(U、V≥1),將U個"Pi"與V個"Pj"進行兩兩配對,產生U*V個組合,每一個"Pi"、"Pj"組合,都可以通過上述(2)計算過程求得相應的擬合平面π'和δ'值.取U*V個δ'中最小者,記作δ'min,如果δ'min小于δ,則該δ'min及所對應的平面π',取代δ與π,然后周而復始繼續計算、判斷,直至滿足終止條件.

(4)程序流程框圖見圖3.

圖3 程序流程圖

(5)編程

選擇既長于計算又具有繪圖功能的C語言作為編程語言.程序(略).

5 程序驗證

上述從理論上證明了算法收斂于"最小區域",還須在實際中應用以檢驗程序的功能.課題組對程序進行了數十個算例進行驗算,算例均取自已發表的其他專家論文、專著或各種標準,將計算結果與原文比較表明本算法計算得到的平面度誤差值之精度都達到或超過原作,是知所擬合之平面亦優于原文.以下取出其中10個算例予以說明,計算結果中略去d值,因其實際上與平面度計算評定無關.

文獻[2]提供的一個算例,被測平面上有25個測量點.原文給出的平面度誤差最佳者為6.55 μm;用本程序計算得出的平面度為6.550 0 μm,擬合基準法矢為T(-0.000002,0.000005,1.00);

文獻[6]之p90給出了一個算例,被測平面上3*3=9個點.原文給出的平面度誤差為0.004 1;用本程序計算,得到的平面度誤差為0.003 9,擬合平面法矢為T(1.997 292,0.997 622,1.000 0);

國標《GBT/1337-1989》[7]給出了兩個算例,分別是被測平面上3*3=9、5*5=25個點.原文計算得出的平面度誤差分別為100、7.7;用本程序進行計算,得出的平面度誤差分別為:

100.000 664,擬合基準法矢T(-0.000 050,0.000 100,1.00);

7.136,擬合基準法矢T(-0.000 002,-0.000 002,1.000);

文獻[8]給出了一個算例,被測平面上10*7=70個點.原文提供的平面度誤差為7.81μm;本程序計算得出的平面度誤差為 7.737 5 μm,擬合平面法矢 T(0.205 000,0.488 125,1.00);

專著[9]及文獻[10]都給出了同一個算例,被測平面上5*5=25個測量點,原文給出的平面度誤差為25 μm;本程序計算得出的平面度誤差為10.1μm,擬合基準平面法矢T(0.000 118,-0.000 014,1.00);

同樣對文獻[3][4][5][11]所提供的算例進行驗算,原文給出的平面度誤差值分別為30.25 μm、6.55 μm、8.755 6 μm、8.5 μm;本程序計算的結果為 30.246 6 μm、6.550 0 μm、6.550 0 μm、8.500 0 μm.

比較彼此的計算結果,顯然本程序計算得出的平面度誤差值精度均優于原文,更符合"最小條件"原則,因此所擬合的基準平面也優于原文.

6 結語

以"最小二乘法"求取的平面為初始擬合基準,繼續以"最小區域"為目標尋找新算法,逐步令誤差計算向"最小區域"逼近,不僅在理論上證明了算法的收斂性,并經過編程及許多實例加以驗證,從運行過程與驗證結果來看,本程序實現了項目開發的第一個目的,運算穩定、結論可靠.本程序在"奔4"微機上運行,一般性算例運行不到一秒鐘,說明收斂速度很快,具備實用意義,為進一步開發高精度三維空間圓度評定程序,打下了決定性的基礎.

本程序的另外一個特點是不要求檢測點按網格均勻分布.為計算上的方便,測量時宜將被測平面置于與三維坐標系的xoy平面上,盡量與z軸垂直(不強求完全垂直).

[1]GB/T 1958-2004,產品幾何量技術規范(GPS)形狀和位置公差檢測規定[S].

[2]JJG 117-2005,平板檢定規程[S].

[3]呂震宇.一種使用最小包容區域法基于旋轉變換求解平面度誤差的方法[J].河北理工學院學報,2000,2(1);47-56.

[4]張之江,于瀛潔,張善鐘.平面度誤差最小區域新算法-有序判別法[J].計量學報,1998,19(1):15-21.

[5]溫秀蘭,宋愛國.基于實數編碼的改進遺傳算法及在平面度誤差評定中的應用[J].計量學報,2003,24(2):88-91.

[6]陳基偉.工業測量數據擬合研究[D].同濟大學工學博士學位論文,2005,11.

[7]GB/T 11337-1989,直線度、平面度誤差檢測[S].

[8]崔長彩,車仁生,羅小川,等.基于實數編碼遺傳算法的平面度評定[J].光學精密工程,2002,10(1):36-40.

[9]羅南星.幾何量測量技術簡明教程[M].機械工業出版社,1993,5.

[10]倪愛晶,鄭聯語.基于形狀誤差不確定度的大尺寸測量系統優化配置方法[J].計量學報,2011,,32(4):289-295.

[11]彭祖行.形狀誤差的優化算法[J].計量學報,1992,13(4):245-250.