譚高山，張麗艷，劉勝蘭

(1.安徽工業(yè)大學 數(shù)理科學與工程學院，安徽 馬鞍山 243032；2.南京航空航天大學 機電學院，江蘇 南京 210016)

0 引言

誤差源于制造不精確，其對零件的裝配質(zhì)量和使用性能有重要影響，因此評估制造工件是否滿足設計誤差要求是生產(chǎn)制造中必不可少的環(huán)節(jié)。評估復雜零部件誤差時，首先利用測量設備獲得零部件的測量模型，然后對齊測量模型與設計模型，確定零部件是否符合設計要求[1]。最小區(qū)域法是符合形位誤差定義和誤差評定相關(guān)標準規(guī)定[2]的對齊方法。對于規(guī)則幾何特征，無論評估方法還是檢測設備均已成熟，但包含自由曲面的復雜零部件的誤差評估仍然很困難。為了簡便易行，工程中仍然使用多種近似方法[3-7]進行誤差評估，包括目前市場上的商業(yè)軟件，近似方法在評估精度要求較高的復雜曲面時，會因過高的估計誤差而導致復雜零部件誤報廢[4]。測量點的代數(shù)擬合方法在數(shù)值上往往不穩(wěn)定，測量點的微小擾動都可能導致擬合參數(shù)明顯改變，造成評估結(jié)果不可靠。幾何擬合又稱為配準或定位，即通過確定一個三維剛體變換來統(tǒng)一或?qū)R測量坐標系和設計坐標系[8]，不同配準方法往往會導致不同的評估結(jié)果。

工業(yè)上通常用幾何公差(GD/T)標準[9]來評估制造工件與其設計要求的一致性?！鞍瑤缀翁卣鞯墓顜Фx了幾何特征的形位公差”，即測量數(shù)據(jù)必須位于理論模型的兩個包容面內(nèi)，包容面間的距離反映了幾何特征的誤差大小。ISO標準[10]規(guī)定了誤差評估的最小區(qū)域準則。極值是輪廓度誤差評估中應用最廣的一個參數(shù)，從這個角度說，最小區(qū)域準則在直觀上是最恰當?shù)脑u估準則，該準則通過定位包含測量模型的兩個等距面來獲得最小輪廓度誤差。

另一常用評估準則是極大極小準則，目標是最小化最大絕對誤差[11-13]，這種評估優(yōu)化模型求解困難，尤其是測量數(shù)據(jù)較多時效率較低。實際上，當幾何模型對稱時，極大極小準則與最小區(qū)域準則等價，特殊情況下可用極大極小準則代替最小區(qū)域準則獲得最小區(qū)域誤差。然而情況并不總是如此，第1章對此進行舉例說明。

求解最小區(qū)域誤差的算法較多[14-19]，如計算幾何法、直接搜索法、向量法、支持向量機法和暴風算法(Brute Force，BF)等，這些方法利用了待評估曲面的特殊幾何性質(zhì)。直接搜索法不但效率較低，而且存在收斂問題；計算幾何法要求計算測量模型的最近和最遠Voronoi圖，為了提高計算精度，需要花費較長的計算時間。另外一些簡單幾何特征的最小區(qū)域配準還可通過解析法實現(xiàn)[20]。以上利用曲面特殊幾何性質(zhì)的評估方法不適用于復雜曲面誤差評估問題。啟發(fā)式智能優(yōu)化算法[21-25]也是求解最小區(qū)域配準的常用算法，如免疫進化算法、擬粒子群算法、微分進化算法、蟻群算法和人工蜂群算法等，然而這些算法非常敏感，且一般效率較低。

本文給出一種基于最小區(qū)域準則的復雜曲面誤差評估配準模型，并研究這一復雜問題的高效求解方法。從數(shù)學上來說，該配準模型最小化了理論曲面兩個等距面之間的距離，兩個等距面包含了待評估模型的所有測量數(shù)據(jù)。由于目標函數(shù)是極值函數(shù)，在算法上找到目標函數(shù)的一致光滑逼近是高效求解配準模型的關(guān)鍵。本文利用凝聚函數(shù)的思想將非光滑配準模型轉(zhuǎn)化為無約束可微優(yōu)化問題，并采用低存儲擬牛頓算法進行高效求解。

首先分析復雜曲面最小輪廓度誤差評估的配準問題并進行建模，通過仿真說明模型不但滿足最小輪廓度誤差定義和曲面評估的最小區(qū)域準則，而且對曲面的幾何性質(zhì)沒有特殊要求。然后研究了模型的高效求解問題，最后通過仿真示例和實測模型驗證了本文算法的適應性和優(yōu)勢。

1 問題分析與數(shù)學建模

設測量數(shù)據(jù)集為P={Pi|i=1,2,…,n}，設計模型為Q。最小區(qū)域配準的目的是找到合適的三維剛體變換，使得理論模型包含所有測量數(shù)據(jù)的兩個等距面之間的距離最小。測量模型相對理論模型的偏差可通過式(1)定義的有向距離函數(shù)表示：

di(R,T)=(RPi+T-Qi)·ni。

(1)

式中Pi(i=1,…,n)為測量點，Qi是Pi在理論模型上的對應點，ni為Qi處的曲面外法向量；R∈3×3和T∈3分別為三維剛體變換的旋轉(zhuǎn)矩陣和平移向量。本文以經(jīng)典的迭代最近點(Iterative Closest Point, ICP)算法的配準結(jié)果作為初始位姿，通過k-d樹搜索算法尋找測量點Pi(i=1,…,n)在理論模型上的最近點，用最近點代替測量點的對應點。配準目的是使對應點對之間的幾何關(guān)系滿足具體工程要求。顯然，有向距離的正負分別表示測量點位于理論模型兩側(cè)。誤差評估要求測量模型相對理論模型的偏差越小越好，因此建立配準模型

(2)

式中I={1,2,…,n}，n為測量數(shù)據(jù)點的個數(shù)，復雜曲面測量模型的數(shù)據(jù)量往往較大。

模型(2)基于最小區(qū)域準則，滿足ISO/TS 和 ANSI 標準。實際上目標函數(shù)是偏差向量函數(shù)D(X)=(d1(X),d2(X),…,dn(X))的一種新的度量，本文稱為極差范數(shù)，模型(2)優(yōu)化的新范數(shù)最小。

為更清楚地描述問題，設空間向量為X=[θx,θy,θz,tx,ty,tz]，則偏差函數(shù)轉(zhuǎn)化為X的多元函數(shù)

di(X)=(R(X)Pi+T(X)-Qi)·ni,

(3)

從而將式(2)轉(zhuǎn)化為式(4)，即以X為優(yōu)化變量的誤差向量函數(shù)的極差范數(shù)最小化問題。

(4)

下面通過模擬多邊形配準說明模型(4)對最小區(qū)域評估的適用性。圖1所示為不同準則下的模擬配準，其中離散點為取樣點，實線表示理論多邊形，虛線表示不同配準結(jié)果的誤差區(qū)域。圖1a所示為ICP結(jié)果，圖1b所示為極大極小準則下的配準結(jié)果，圖1c所示為本文所提最小區(qū)域的配準的結(jié)果，這3種方法分別優(yōu)化了理論模型和測量數(shù)據(jù)之間的歐氏距離誤差向量的2-范數(shù)、1-范數(shù)和有向距離誤差向量的極差范數(shù)。由圖可知，所提方法的輪廓度誤差最小，即獲得了最小誤差帶。對于一些簡單幾何特征，最小區(qū)域配準結(jié)果的最優(yōu)性和唯一性已經(jīng)得到證明[2]。實際上，考慮計算的復雜性和數(shù)據(jù)模型的特殊性，常見一些文獻用極大極小配準模型來實現(xiàn)最小區(qū)域配準，這一數(shù)學模型優(yōu)化了兩數(shù)據(jù)模型之間的最大歐氏距離。例如對一些簡單幾何特征，以及一些具有對稱性或非封閉的曲面進行極大極小配準和所提最小區(qū)域配準時，最小輪廓度誤差一致，但是情況并不總是如此。圖1示例對此進行了說明，分析發(fā)現(xiàn)有向距離誤差向量的極差范數(shù)是最小區(qū)域誤差的最佳度量，因此式(4)是具有明確最小區(qū)域的配準模型，其對待配準曲面的幾何性質(zhì)沒有特殊要求，適合高精度的復雜曲面評估。

以上分析說明式(4)所示配準模型是曲面最小區(qū)域評估的合理表示，這一優(yōu)化模型的目標函數(shù)非光滑，由于非光滑函數(shù)下降方向的搜索比光滑函數(shù)困難得多，下面著重解決這一問題。

2 凝聚函數(shù)逼近及優(yōu)化求解

非光滑優(yōu)化問題一般采用直接搜索法或啟發(fā)式搜索算法求解，相對傳統(tǒng)的微分類優(yōu)化算法，這些搜索算法的效率和精度往往不夠理想，有時甚至無法保證收斂性。熵正則化方法[26-27]利用凝聚函數(shù)一致光滑逼近非光滑函數(shù)，凝聚函數(shù)是帶有控制參數(shù)的顯式可微函數(shù)。本文利用凝聚函數(shù)，分別考慮最小值、最大值函數(shù)的一致光滑逼近，最小值函數(shù)等價于

(5)

因此將式(4)轉(zhuǎn)化為

(6)

需要強調(diào)的是，式(6)是式(4)的另一種表達形式。因為di(X)及其相反函數(shù)的最優(yōu)值并不能保證同時達到，所以式(6)與極大極小問題的本質(zhì)不同。

參照筆者以前的工作[28]，推導得式(6)目標函數(shù)的二階可微凝聚函數(shù)近似Fp(X)，

(7)

式中控制參數(shù)p在理論上趨向正無窮大。因此式(4)轉(zhuǎn)化為如下無約束優(yōu)化問題：

(8)

BFGS(Broyden-Fletcher-Goldfarb-Shanno) 算法是求解這類問題最常用的典型擬牛頓算法，僅需計算函數(shù)值和梯度值。

記di(X)(i=1,…,n)及其相反函數(shù)在第k次迭代時的函數(shù)值D+(X(k))和D-(X(k))：

(9)

為方便表示，采用式(10)表示相關(guān)點的下標集合：

E+(X(k))={i|di(X(k))=D+(X(k)),i∈I};

E-(X(k))={i|-di(X(k))=D-(X(k)),i∈I};

O+(X(k))={i|p·(D+(X(k))-

O-(X(k))={i|p·(D-(X(k))-

(10)

式中L為數(shù)值計算中數(shù)據(jù)溢出的臨界值。函數(shù)值和梯度值分別表示為：

(11)

(12)

式中：

ui(X(k))=

vi(X(k))=

為解決數(shù)據(jù)溢出問題,采用L-BFGS 算法求解式(8)，該算法較BFGS算法具有更高的數(shù)據(jù)存儲能力?？紤]參數(shù)p對配準精度、效率、數(shù)值穩(wěn)定性和收斂性的影響，下一章驗證算法對參數(shù)的敏感性，由實驗知參數(shù)p可在一個較大的范圍內(nèi)進行有效取值。

3 實驗

為驗證算法性能，本章比較經(jīng)典的ICP配準、序列二次規(guī)劃方法求解的極大極小配準和本文用L-BFGS算法求解的最小區(qū)域配準的結(jié)果。實驗在2.10 GHz的CPU 和2.5 GB的RAM個人機上采用VC++編程實現(xiàn)。為公平起見，不同方法的數(shù)據(jù)結(jié)構(gòu)、最近點的搜索、有向距離的計算等其他方面都采用相同的處理方法。

3.1 仿真數(shù)據(jù)實驗

為更清楚地比較所提最小區(qū)域配準法與極大極小配準法的區(qū)別與聯(lián)系，本節(jié)以平面為例進行實驗。具體為在商業(yè)軟件中模擬生成平面，通過隨機擾動獲得平面測量模型的118個采樣點，并經(jīng)微小剛體變換作為初始狀態(tài)，ICP算法配準獲得的輪廓度誤差為4.90 mm，本文算法下降到4.52 mm,相對ICP配準結(jié)果，本文配準結(jié)果更優(yōu)。極大極小配準輪廓度誤差也是4.52 mm，這與前文分析一致，即對于一些特殊幾何特征(如對稱幾何特征)，極大極小配準和本文最小區(qū)域包容配準所得的最小輪廓誤差相同。3種方法的誤差分布情況如表1所示，可見極大極小方法中的總體平均誤差、正向平均誤差和負向平均誤差均大于本文方法，顯然沒有一種方法可以獲得比本文方法更小的輪廓度誤差，這是因為本文方法直接優(yōu)化了表示輪廓度誤差的極差參數(shù)。ICP法直接優(yōu)化最小均方根誤差，因此獲得了3種方法中最小的均方根誤差。

表1 仿真平面配準評估結(jié)果 mm

3.2 實測鈑金件誤差評估實驗

本文實驗研究某實測液壓成型鈑金件的誤差評估配準問題，該鈑金件由包括自由曲面在內(nèi)的復雜曲面構(gòu)成，變形導致鈑金件的誤差分布難以估計，因此裝配前在算法上對齊測量模型與其理論模型具有重要意義。首先利用流行商業(yè)軟件Imageware? 中的“best fit”功能獲得初始位姿，然后采用ICP方法和本文方法分別進行配準。圖2所示為所提方法的配準結(jié)果，可見誤差分布相對集中，最小輪廓度誤差為1.74 mm，而ICP算法為1.83 mm。對于高精度要求的鈑金零件評估來說，精度的微小改善都具有重要意義。表2所示為不同方法的誤差統(tǒng)計結(jié)果，可見無論正向還是負向，本文方法的最大誤差都更小，而且本文方法的正向均方根誤差0.21 mm也比ICP算法的0.23 mm小，但是兩種方法的負向誤差分布大體一致。為更清楚地反映不同算法的誤差分布，對不同誤差區(qū)間的樣本數(shù)據(jù)個數(shù)進行統(tǒng)計。ICP算法和本文方法在誤差區(qū)間(-0.5, 1.5)內(nèi)的數(shù)據(jù)點數(shù)分別為74.56%和70.18%，驗證了基于最小二乘準則的ICP算法的誤差分布在整體上更均衡。實際上，誤差評估的最小區(qū)域準則是要獲得最小的輪廓度誤差，常用的最小二乘準則雖然獲得了誤差均衡，但是并不能保證輪廓度誤差最小，若將誤差區(qū)間縮小為(-0.5, 0.5)，則兩種方法對應的數(shù)據(jù)點數(shù)分別為35.96%和35.08%，比較接近，因此本文方法不但滿足了誤差評估的最小區(qū)域準則，而且高精度區(qū)域的誤差分布也比較均衡。

表2 實測鈑金件的配準評估比較mm

4 結(jié)束語

本文提出一種面向復雜曲面誤差評估的最小區(qū)域包容配準方法，通過建立具有明確物理意義的極差配準模型，并利用經(jīng)典擬牛頓法求解配準模型的一致光滑無約束逼近問題，實現(xiàn)了復雜曲面的最小輪廓度誤差評估。通過模擬仿真數(shù)據(jù)和實測鈑金件的誤差評估實驗說明所提算法有助于減少生產(chǎn)中由于過高估計誤差引起的誤報廢問題，從而驗證了算法的適用性和優(yōu)勢。