毛喆，李瀧杲，*，徐巖，曾琪，主逵

（1.南京航空航天大學(xué) 機(jī)電學(xué)院，南京210016； 2.深圳市勁拓自動化設(shè)備股份有限公司，深圳518101）

現(xiàn)代飛機(jī)高機(jī)動、高氣動的產(chǎn)品需求對裝配質(zhì)量優(yōu)化提出了巨大挑戰(zhàn)，數(shù)字化測量作為飛機(jī)裝配中的重要環(huán)節(jié)，測量質(zhì)量的優(yōu)劣性嚴(yán)重影響飛機(jī)的實(shí)際裝配結(jié)果。目前高精度、高效率的數(shù)字化測量技術(shù)正逐步代替?zhèn)鹘y(tǒng)的檢測方式，但測量實(shí)施前需基于模型進(jìn)行測量規(guī)劃：一是測量點(diǎn)規(guī)劃，二是測量設(shè)備站位規(guī)劃。前者是后者乃至測量實(shí)施的基礎(chǔ)，后者是前者優(yōu)化的依據(jù)。測量點(diǎn)規(guī)劃是將待測特征離散為一組點(diǎn)作為理論測量點(diǎn)（簡稱測量點(diǎn)），用于站位規(guī)劃、測量實(shí)施、模型重構(gòu)、質(zhì)量分析等環(huán)節(jié)，但當(dāng)前測量點(diǎn)規(guī)劃依賴于工藝人員經(jīng)驗(yàn)，缺乏理論依據(jù)與評價機(jī)制，極易出現(xiàn)測量效率低、實(shí)測數(shù)據(jù)分析量大等問題。因此，合理規(guī)劃測量點(diǎn)對測量效率、測量精度的提升具有重要意義。

Cho等［1］根據(jù)待測曲面面積、測量精度等因素利用模糊神經(jīng)網(wǎng)絡(luò)確定布點(diǎn)數(shù)量，為三坐標(biāo)測量機(jī)布點(diǎn)提供了依據(jù)，但在布設(shè)過程中未考慮到自由曲面的復(fù)雜性。Lee等［2］結(jié)合Hammersley序列提出了一種布點(diǎn)方法，有效地減少了布點(diǎn)數(shù)量。宋占杰等［3］將自由曲面的曲率函數(shù)作為生成質(zhì)心Voronoi結(jié)構(gòu)中的密度函數(shù)，以成本函數(shù)收斂性作為算法結(jié)束準(zhǔn)則，提出了一種全新的基于質(zhì)心Voronoi的檢測點(diǎn)采樣方法。劉紅軍等［4］提出了實(shí)時重構(gòu)布點(diǎn)策略，但需要在測量時對測量數(shù)據(jù)實(shí)時重構(gòu)，測量效率低下。此外，諸多學(xué)者對測量點(diǎn)布設(shè)展開了研究，但多以特征的曲率特性作為測量點(diǎn)布設(shè)依據(jù)，并未考慮實(shí)際的測量不確定度［5-6］。激光雷達(dá)等測量設(shè)備在測量距離較遠(yuǎn)時，測量不確定度會達(dá)到甚至超過0.1mm，此時的測量數(shù)據(jù)難以精準(zhǔn)反映零件的實(shí)際狀態(tài)，需要采用相關(guān)算法減小測量不確定度。

本文采用非均勻有理B樣條（Non-Uniform Rational B-Spline，NURBS）理論精確擬合自由曲線或曲面，并利用粒子群優(yōu)化算法優(yōu)化控制點(diǎn)和權(quán)因子提高擬合精度。對待測特征進(jìn)行區(qū)域劃分，包括曲率極值區(qū)和測量不確定度較大區(qū)，并針對不同的區(qū)域設(shè)計相應(yīng)的布點(diǎn)策略確保測量點(diǎn)能夠精確地表達(dá)待測特征。在有效保證測量效率的前提下，使測量點(diǎn)逆向重構(gòu)特征與理論特征偏差優(yōu)于10-3mm。

1 問題描述

對于復(fù)雜飛機(jī)零部件而言，曲面特征多以自由造型設(shè)計，測量點(diǎn)規(guī)劃主要取決于工藝人員的經(jīng)驗(yàn)。過多的測量點(diǎn)會延長測量周期，過少的測量點(diǎn)無法精確測得零件的實(shí)際狀態(tài)。此外，當(dāng)前測量點(diǎn)布設(shè)大多以曲率特性作為依據(jù)，忽視了測量不確定度對測量結(jié)果的影響，測量不確定度較大位置測量數(shù)據(jù)有時難以反映零件的實(shí)際狀態(tài)。即使可以增加設(shè)備站位減小測量不確定度，但隨之產(chǎn)生的是測量效率降低、數(shù)據(jù)轉(zhuǎn)換誤差累積等問題。

試驗(yàn)件曲率極值區(qū)域及測量不確定度較大區(qū)域分布如圖1所示。圖中曲率極值區(qū)域并未完全覆蓋測量不確定度較大區(qū)域，在盡量避免轉(zhuǎn)站前提下，僅以曲率作為布點(diǎn)依據(jù)難以精準(zhǔn)測得零件的實(shí)際狀態(tài)。

圖1 試驗(yàn)件測量不確定度及曲率分布Fig.1 Measurement uncertainty and curvature distribution of test piece

2 待測特征確定性表達(dá)構(gòu)建

飛機(jī)大尺寸零部件待測特征主要為曲線、曲面等基本特征，曲面特征一般是通過離散為一組曲線、采用NURBS精確擬合各條曲線進(jìn)行分析。鑒于3次NURBS曲線具有C2連續(xù)性且能夠滿足工程中的造型需求［7］，本文的研究基于3次NURBS曲線展開。

2.1 曲線方程參數(shù)化

一條p次NURBS曲線可定義為［8］

式中：Ni，p為Bernstein基函數(shù)；u為自由變量。

自由曲線表征為NURBS初始需求解控制點(diǎn)Vi、權(quán)因子wi及節(jié)點(diǎn)矢量U，其中節(jié)點(diǎn)矢量通過弦長參數(shù)化確定。權(quán)因子和控制點(diǎn)的耦合求解過程復(fù)雜，極易出現(xiàn)奇異解或無解的情況，為簡化計算，令wi初始值為1，控制點(diǎn)可由式（2）解得：

式中：Pi為型值點(diǎn)。

控制點(diǎn)求解后曲線初始方程r0（u）為

式（3）確定的擬合曲線往往不能精確表達(dá)理論曲線，一般采用增加型值點(diǎn)的方式提高擬合曲線的精度，但型值點(diǎn)數(shù)量增加到一定程度后對曲線精度的提高效果不明顯，有時甚至?xí)a(chǎn)生扭曲和畸變［9］，需要對wi和Vi進(jìn)行優(yōu)化提高曲線的擬合精度。

2.2 參數(shù)化方程優(yōu)化

2.2.1 局部優(yōu)化

NURBS曲線形狀可通過調(diào)節(jié)控制點(diǎn)和權(quán)因子進(jìn)行修改，權(quán)因子wi負(fù)責(zé)曲線形狀的局部修改。針對擬合曲線中局部區(qū)域偏離理論曲線的情況，可通過調(diào)節(jié)wi提高曲線的擬合精度，如圖2所示。

設(shè)調(diào)整后權(quán)因子為w′i，wi與w′i的關(guān) 系可表示為［10］式中：d為擬合曲線與理論線間的最大偏差。

圖2 權(quán)因子對曲線形狀的影響Fig.2 Influence of weight factor on curve shape

對于僅有少數(shù)幾個部位誤差較大的情況，依據(jù)式（4）對曲線形狀微調(diào)，重復(fù)采用此方法快速提高參數(shù)方程正確性。但是對于整體偏離理論曲線的情況，修改的權(quán)因子數(shù)量多，優(yōu)化效率低。

2.2.2 整體優(yōu)化

曲線的整體優(yōu)化方法為對Vi及wi共同調(diào)整，但由于當(dāng)前控制點(diǎn)和權(quán)因子間未形成確定的關(guān)系式，調(diào)整過程中Vi和wi的設(shè)置直接影響優(yōu)化結(jié)果及速度。鑒于粒子群優(yōu)化算法［11］具有收斂速度快、參數(shù)少等優(yōu)點(diǎn)，本文采用粒子群優(yōu)化算法對Vi及wi進(jìn)行優(yōu)化。

以一組控制點(diǎn)及其權(quán)因子作為一個優(yōu)化粒子Zj，優(yōu)化問題解的維度D＝4n（n為控制點(diǎn)個數(shù)，x，y，z為控制點(diǎn)坐標(biāo)），第j個粒子定義為

設(shè)第j個粒子的最優(yōu)位置為Pbestj，整個群體經(jīng)歷的最優(yōu)位置為Gbest，第j個粒子的速度及位置更新公式分別為

式中：w為慣性權(quán)重；c1、c2為學(xué)習(xí)因子；rand1、rand2為0～1之間的隨機(jī)數(shù)。設(shè)第j個粒子Zj構(gòu)成的NURBS曲線為rkj（u），將rkj（u）進(jìn)行離散，計算每個離散點(diǎn)到理論線的最小距離，形成距離集Dj＝｛dj1，dj2，…，djm｝（m為離散點(diǎn)個數(shù)），構(gòu)造適應(yīng)度函數(shù)為

式中：F為粒子個數(shù)。結(jié)合局部優(yōu)化和整體優(yōu)化方法，待測特征精確參數(shù)化總體流程如下：

步驟1 采用等弦長法將理論線r離散為一組型值點(diǎn)P＝｛P0，P1，…，Pn｝，以P為初始條件計算初始曲線方程r0（u）并生成擬合曲線。

步驟2 計算擬合曲線與理論線的最大偏差δ，若δ大于給定閾值ε0，回到步驟1增加型值點(diǎn)數(shù)量以減小δ，否則執(zhí)行步驟3。

步驟3 判斷擬合曲線與理論線偏離狀態(tài)，若兩曲線偏離距離δ小于閾值ε1，執(zhí)行步驟7；若局部偏離執(zhí)行步驟4；若整體偏離執(zhí)行步驟5。

步驟4 計算r與擬合曲線最大偏差δ，按式（4）修改權(quán)因子形成新的曲線方程rk（u）。計算r與rk（u）的最大偏差δ，執(zhí)行步驟5。

步驟5 若δ大于閾值ε1，重復(fù)步驟4修改下一處權(quán)因子，直到δ變化不明顯或小于ε1，執(zhí)行步驟7。

步驟6 采用粒子群優(yōu)化算法對控制點(diǎn)及權(quán)因子整體優(yōu)化，若δ變化不明顯或小于ε1，執(zhí)行步驟7。

步驟7 輸出最終的控制點(diǎn)、權(quán)因子及參數(shù)化曲線方程r（u）。

3 待測特征測量點(diǎn)布設(shè)

飛機(jī)大型結(jié)構(gòu)件、薄壁等零件易產(chǎn)生裝配變形，曲率較大的位置變形較難控制［12］，且在測距較遠(yuǎn)時測量不確定度較大，測量結(jié)果不準(zhǔn)確。針對上述情況，綜合考慮曲率特性和測量不確定度等因素布設(shè)測量點(diǎn)。

3.1 曲線測量點(diǎn)布設(shè)

3.1.1 曲率極值點(diǎn)求解

曲率極值點(diǎn)是描述曲線形狀的關(guān)鍵點(diǎn)，為了獲取曲率極值處的變形情況，需要在曲率極值點(diǎn)處布設(shè)測量點(diǎn)。設(shè)曲線的參數(shù)方程為r（u）＝｛x（u），y（u），z（u）｝，空間曲線在某一處的曲率k（u）表示為［13］

式中：r′（u）與r″（u）分別為曲線r（u）的一階導(dǎo)數(shù)與二階導(dǎo)數(shù)，將其代入式（9）可得

對k（u）求一階導(dǎo)數(shù)為

曲率極值點(diǎn)為?k／?u＝0的解，但此方程求解過于復(fù)雜，針對這種情況本文提出了一種基于曲率導(dǎo)數(shù)的曲率極值點(diǎn)快速定位方法，具體方法流程如下：

步驟1 將曲線方程r（u）離散，記錄每個離散點(diǎn)對應(yīng)的ui。

步驟2 計算每個離散點(diǎn)的曲率導(dǎo)數(shù)J（ui），J（ui）＝?k／?u｜u＝ui。

步驟3 計算相鄰兩點(diǎn)曲率導(dǎo)數(shù)積Di＝J（ui）·J（ui+1），執(zhí)行步驟4。

步驟4 遍歷所有的Di值，對Di≤0的兩點(diǎn)T0、Q0和ui、ui+1保存，對保存的點(diǎn)對逐一執(zhí)行步驟5和步驟6。

步驟5 取ui、ui+1的中間值um代入r（u），在T0、Q0中間生成點(diǎn)M0。分別計算J（um）·J（ui）與J（um）·J（ui+1）的值。若J（um）·J（ui）≤0，以中間點(diǎn)取代Q0，反之，以中間點(diǎn)取代T0。執(zhí)行步驟6。

步驟6 重復(fù)步驟5更新T0、Q0，直到T0、Q0距離小于給定閾值ε2，以兩點(diǎn)中間點(diǎn)作為曲率極值點(diǎn)。

圖3為曲線某一個曲率極值點(diǎn)求解，左邊界由T0更新至T4，右邊界點(diǎn)由Q0更新至Q4。

圖4為某零件邊界曲線的曲率極值求解，圓“·”為采用本文方法求解出的曲率極值點(diǎn)，“”為CATIA（Computer Aided Three-dimensional Interactive Application）中曲率變化曲線上的極值點(diǎn)。從圖中可以看出，采用本文方法求解的曲率極值與CATIA中的曲率極值基本一致。

圖3 曲率極值點(diǎn)搜索Fig.3 Search of curvature extreme points

圖4 曲線曲率極值點(diǎn)求解結(jié)果Fig.4 Solving result of curve curvature extreme points

3.1.2 測量不確定度評估

為降低測量不確定度對測量結(jié)果的影響，在測量不確定度較大區(qū)域加密布設(shè)測量點(diǎn)，通過間接平差的方式減小此影響。本文以球坐標(biāo)測量系統(tǒng)為研究對象對測量不確定度評估。任一測量點(diǎn)M（x，y，z）在測量坐標(biāo)系的坐標(biāo)可由距離l、方位角α及天頂角β解得。

由于設(shè)備、環(huán)境因素，l、α、β采集的過程中均存在著一定的誤差，為評價測量點(diǎn)在空間的不確定度，建立不確定度橢球模型［14-15］如圖5所示。

不確定橢球的半軸長Ux′、Uy′、Uz′分別為［16-17］

式中：σα、σβ與測量系統(tǒng)的角度編碼器分辨率有關(guān)，為儀器制造商提供的固定值；σl由最小測距不確定分量和測距不確定度距離放大系數(shù)組成。

圖6為曲線輪廓度及測量不確定度允許偏差示意圖，測量不確定度允許偏差范圍一般為輪廓度的十分之一。圖中兩端及曲率較大處測量不確定度橢球較大，有的甚至超出不確定度允許偏差范圍。為了減小測量不確定度，提出一種以測量不確定度為依據(jù)的測量點(diǎn)布設(shè)方法。具體流程如下：

步驟1 選取合適位置作為測量設(shè)備站位E（x，y，z）。

步驟2 根據(jù)曲線方程r（u），對曲線進(jìn)行離散，生成離散點(diǎn)集G＝｛G0，G1，…，Gn｝。

圖5 不確定度橢球模型Fig.5 Uncertainty ellipsoid model

圖6 曲線輪廓度公差帶及測量不確定度允差Fig.6 Tolerance zone and measurement uncertainty tolerance of curve

步驟3 計算每個離散點(diǎn)Gi到E的距離li，根據(jù)式（13）計算出每個橢球的半軸長。

步驟4 以Gi為原點(diǎn)，Gi與E連線為u軸，在x′-y′-z′坐標(biāo)系構(gòu)建單側(cè)不確定度橢球并生成橢球離散點(diǎn)。

步驟5 通過坐標(biāo)轉(zhuǎn)換將不確定度橢球離散點(diǎn)由x′-y′-z′坐標(biāo)系轉(zhuǎn)換到x-y-z坐標(biāo)系中。

步驟6 計算每個橢球離散點(diǎn)距離理論線的最遠(yuǎn)點(diǎn)，將最遠(yuǎn)點(diǎn)擬合成為一條曲線，對其進(jìn)行光順后作為不確定度曲線。

步驟7 在不確定度曲線超出不確定度允許范圍的區(qū)域布設(shè)若干測量點(diǎn)。

3.2 曲面測量點(diǎn)布設(shè)

曲面測量點(diǎn)布設(shè)是將曲面離散成為一組截交線［18］，按照3.1節(jié)方法對每條截交線布設(shè)測量點(diǎn)。大尺寸曲面測量過程中，測量數(shù)據(jù)同樣存在測量不確定度，圖7為曲面輪廓度公差帶和測量不確定度允差的分布情況。為了減小測量不確定度的影響，構(gòu)建每個測量點(diǎn)的單側(cè)不確定度橢球并生成橢球離散點(diǎn)，找出離散點(diǎn)中距離理論面最遠(yuǎn)點(diǎn)。將最遠(yuǎn)點(diǎn)擬合成為不確定度曲面，計算出不確定度曲面超出不確定度允許偏差范圍的區(qū)域，并在這些區(qū)域加密布設(shè)測量點(diǎn)。最后將截交線離散點(diǎn)及根據(jù)測量不確定度增加的測量點(diǎn)進(jìn)行篩選，對于距離過小甚至重合的兩點(diǎn)取中間點(diǎn)代替。

采用上述布點(diǎn)策略布設(shè)的測量點(diǎn)有時無法精確描述理論特征，需要對測量點(diǎn)補(bǔ)充。測量點(diǎn)補(bǔ)充方法是將依據(jù)曲率及測量不確定度布設(shè)的測量點(diǎn)進(jìn)行重構(gòu)，找出重構(gòu)特征與理論特征偏差較大位置，并在該位置補(bǔ)充測量點(diǎn)。以測量點(diǎn)重構(gòu)曲面與理論曲面的最大距離偏差作為測量點(diǎn)能否精確描述待測特征的指標(biāo)，待測特征測量點(diǎn)布設(shè)整體流程如圖8所示，圖中ε3為重構(gòu)特征與理論特征的閾值。

圖7 曲面輪廓度公差帶及測量不確定度允差Fig.7 Tolerance zone and measurement uncertainty tolerance of surface

圖8 測量點(diǎn)布設(shè)流程Fig.8 Process ofmeasurement points distribution

4 實(shí)驗(yàn)驗(yàn)證

為實(shí)現(xiàn)曲線曲面的自動布點(diǎn)，在CATIA［19］環(huán)境下利用CAA［20］（Components Application Architecture）技術(shù)開發(fā)測量點(diǎn)布設(shè)模塊。軟件實(shí)現(xiàn)包含特征預(yù)處理、特征分析及測量點(diǎn)規(guī)劃、測量點(diǎn)輸出3部分，如圖9所示。

以圖1中的試驗(yàn)件作為驗(yàn)證對象，首先提取如圖10（a）所示的待測特征；隨后進(jìn)行特征分析和測量點(diǎn)布設(shè)，將曲面離散為一組截交線，對截交線按照3.1節(jié)方法規(guī)劃測量點(diǎn)，求解出測量不確定度較大區(qū)域，并在該區(qū)域采用3.2節(jié)布點(diǎn)方法加密布設(shè)測量點(diǎn)，測量點(diǎn)布設(shè)結(jié)果如圖10（b）所示；最后輸出測量點(diǎn)，如圖10（c）所示。為使參數(shù)方程和測量點(diǎn)精確擬合待測特征，閾值ε1、ε1優(yōu)于10-3mm。為達(dá)到測量不確定度較大區(qū)域加密布點(diǎn)的目的，此區(qū)域測量點(diǎn)重構(gòu)閾值取普通區(qū)域閾值ε3的一半。

圖9 自動布點(diǎn)模塊架構(gòu)Fig.9 Architecture of automatic points distribution module

圖10 試驗(yàn)件局部曲面片布點(diǎn)過程Fig.10 Points distribution process of partial surface of test piece

表1為粒子群優(yōu)化算法所采用的參數(shù)，vmax、vmin分別為最大、最小更新速度，Np為迭代次數(shù)。

為驗(yàn)證布點(diǎn)合理性，采用傳統(tǒng)布點(diǎn)方法在試驗(yàn)件上布點(diǎn)密度為20×20和10×10的測量點(diǎn)，如圖11（a）和圖11（b）所示，圖11（c）為采用本文方法布設(shè)的測量點(diǎn)。利用激光雷達(dá)MV260作為測量工具（其角度精度6.8μm／m，距離精度10μm+2.5μm／m），對3種情況進(jìn)行現(xiàn)場測量。

分別構(gòu)建3種布點(diǎn)方法中的理論點(diǎn)的單側(cè)不確定度橢球，將不確定度橢球離散，計算每個橢球離散點(diǎn)集中距離理論面的最遠(yuǎn)點(diǎn)，將最遠(yuǎn)點(diǎn)擬合為不確定度曲面，以不確定度曲面距離理論面的最大距離作為擬合精度，3種布點(diǎn)方法的測量效率及擬合精度如表2所示。試驗(yàn)件測量點(diǎn)分布情況如圖12所示。

表1 粒子群優(yōu)化算法參數(shù)Tab le 1 Param eters of par ticle swarm optim ization algorithm

圖11 不同布點(diǎn)密度布點(diǎn)結(jié)果對比Fig.11 Comparison of points distribution under differentmeasurement points density

表2 不同布點(diǎn)密度下測量結(jié)果對比Table 2 Com parison of m easurem ent result under different m easurem ent points density

圖12 試驗(yàn)件測量點(diǎn)分布Fig.12 Measurement points distribution of test piece

從表2中可以看出，采用本文方法布點(diǎn)測量效率介于布點(diǎn)密度為10×10和20×20之間，重構(gòu)的孔洞數(shù)量、擬合精度與布點(diǎn)密度為10×10的結(jié)果相當(dāng)，但明顯優(yōu)于布點(diǎn)密度為20×20的測量點(diǎn)重構(gòu)的孔洞數(shù)量及擬合精度，充分體現(xiàn)了本文布點(diǎn)策略的合理性及可行性。

5 結(jié) 論

針對航空大尺寸零部件測量點(diǎn)規(guī)劃的基礎(chǔ)問題，本文分別從待測特征曲率特性、測量不確定度角度出發(fā)對測量點(diǎn)進(jìn)行規(guī)劃，研究了測量點(diǎn)布設(shè)的理論依據(jù)及布設(shè)方法。主要結(jié)論如下：

1）將待測特征精確參數(shù)化，并對參數(shù)方程優(yōu)化實(shí)現(xiàn)擬合精度優(yōu)于10-3mm，為復(fù)雜特征的曲率及測量不確定度分析提供依據(jù)。

2）為復(fù)雜曲線曲面測量點(diǎn)規(guī)劃提供了理論依據(jù)，在保證測量點(diǎn)精確表征待測特征的條件下控制測量點(diǎn)數(shù)量，使測量點(diǎn)逆向特征與理論特征的最大偏差優(yōu)于10-3mm。

3）開發(fā)了測量點(diǎn)自動布設(shè)模塊，提高測量點(diǎn)布設(shè)效率及規(guī)范性。