田清廉，熊天辰，黃 翔，李瀧杲，郝 龍

（1.南京航空航天大學(xué)機(jī)電學(xué)院，南京 210016；2.江西航天海虹測(cè)控技術(shù)有限責(zé)任公司，南昌 330024）

我國(guó)航空事業(yè)不斷發(fā)展，對(duì)零件上鉚釘孔的制孔精度要求也越來越高，而飛機(jī)結(jié)構(gòu)與系統(tǒng)中使用的大量零件包含鉚釘孔特征[1]。通常對(duì)孔特征檢測(cè)時(shí)，使用三坐標(biāo)測(cè)量?jī)x測(cè)出各個(gè)孔的坐標(biāo)尺寸，但這種測(cè)量方式效率低，且不適合檢測(cè)壁板等大型零件；使用專用孔檢測(cè)工具 （如孔檢測(cè)樣板）對(duì)孔的尺寸和形位進(jìn)行檢測(cè)，但制造這種高精度檢具耗時(shí)長(zhǎng)且成本高[2]。隨著非接觸式數(shù)字化測(cè)量技術(shù)的提高，非接觸式數(shù)字化測(cè)量設(shè)備也開始更多地應(yīng)用在飛機(jī)零件鉚釘孔檢測(cè)上，這種檢測(cè)方式具有精度高、效率高、柔性好等特點(diǎn)[3–4]。

目前，在非接觸式圓孔測(cè)量方面已有了很多的研究與成果，主要有基于視覺圖像的測(cè)量方法和基于三維點(diǎn)云的測(cè)量方法。基于視覺圖像的測(cè)量方法主要是通過圖像處理中的Canny 算子等邊界識(shí)別算法提取照相測(cè)量圖像中的邊界特征進(jìn)行圓孔參數(shù)擬合[5–6]，但該方法存在空間圓透視投影畸變的問題[7]，影響檢測(cè)精度。基于三維點(diǎn)云的測(cè)量方法是使用三維激光掃描等測(cè)量設(shè)備獲取零件點(diǎn)云，再提取點(diǎn)云中的孔邊界特征進(jìn)行圓孔參數(shù)擬合。針對(duì)點(diǎn)云中邊界特征提取，Jenke 等[8]提出了一種將點(diǎn)云網(wǎng)格化再通過網(wǎng)格的拓?fù)浣Y(jié)構(gòu)識(shí)別邊界特征的方法，這種方法可以有效提取點(diǎn)云孔洞邊界特征，但需要網(wǎng)格化整個(gè)點(diǎn)云，數(shù)據(jù)量大時(shí)算法效率低。因此出現(xiàn)了直接對(duì)三維點(diǎn)云提取孔邊界特征的方法，該方法通過建立點(diǎn)云空間拓?fù)潢P(guān)系，并設(shè)定某一特定判斷機(jī)制用于識(shí)別邊界特征點(diǎn)。Milroy[9]和Yang[10]等通過求解曲率極值提取點(diǎn)云邊界，這種方法可以有效地提取比較光滑平坦的數(shù)據(jù)模型的邊界點(diǎn)，不太適用于曲率變化大的數(shù)據(jù)模型。張杰等[11]提出了一種基于T–Scan測(cè)量的薄壁鈑金件孔特征重構(gòu)方法，通過掃描線點(diǎn)云結(jié)構(gòu)中的孔截?cái)嗑€特征識(shí)別平面點(diǎn)和孔徑點(diǎn)，擬合求解平面和空間圓，但該方法僅限于掃描線點(diǎn)云類型且提取精度受圓孔直徑影響。劉增藝等[12]提出了一種曲面孔位視覺測(cè)量技術(shù)，使用手持式雙目視覺線結(jié)構(gòu)光測(cè)量曲面得到曲面點(diǎn)云，然后通過交互式方法提取散亂點(diǎn)云圓孔邊緣點(diǎn)，進(jìn)行圓孔擬合，由于該方法在圓孔識(shí)別過程中需要人工預(yù)先分割出圓孔所在區(qū)域，提取效率低且不能滿足多孔自動(dòng)提取需求。

基于以上分析，為解決目前面向鉚釘孔檢測(cè)領(lǐng)域的測(cè)量點(diǎn)云中的孔位與孔徑參數(shù)提取方法限制多和人工交互較多的問題，本文提出了一種基于散亂點(diǎn)云的鉚釘孔孔位與孔徑提取方法，該方法可以自動(dòng)識(shí)別散亂點(diǎn)云中的全部鉚釘孔特征，并精確提取孔位和孔徑參數(shù)。采用本文方法對(duì)試驗(yàn)件散亂點(diǎn)云進(jìn)行數(shù)據(jù)處理，驗(yàn)證該方法的實(shí)用性與精度。

1 散亂點(diǎn)云預(yù)處理

各種數(shù)字化測(cè)量設(shè)備獲得的點(diǎn)云根據(jù)其點(diǎn)的分布特征分為4 種類型：散亂點(diǎn)云、掃描線點(diǎn)云、網(wǎng)格化點(diǎn)云和多邊形點(diǎn)云[13]。散亂點(diǎn)云相對(duì)于其他點(diǎn)云類型，其測(cè)量點(diǎn)無明顯的幾何分布特征，呈散亂無序狀態(tài)，測(cè)量點(diǎn)是隨機(jī)生成的。由于散亂點(diǎn)云無序和隨機(jī)的分布特性，使用非接觸式測(cè)量設(shè)備測(cè)量鉚釘孔得到的點(diǎn)云分布情況如圖1所示，大部分測(cè)量點(diǎn)分布在孔周圍平面上，而孔邊界測(cè)量點(diǎn)可能分布在孔內(nèi)徑或在孔周圍平面上，有極少的點(diǎn)正好在孔邊界上。

圖1 鉚釘孔邊界測(cè)量點(diǎn)云分布Fig.1 Distribution of measurement point cloud at rivet hole boundary

根據(jù)鉚釘孔邊界測(cè)量點(diǎn)云分布情況，本文通過點(diǎn)的k鄰域點(diǎn)幾何分布來判斷是否為邊界點(diǎn)。若某一點(diǎn)為邊界點(diǎn)，則其k鄰域點(diǎn)的分布將偏向于一側(cè) （圖2（a））；若某一點(diǎn)不是邊界點(diǎn)，則k鄰域點(diǎn)比較均勻地分布在該點(diǎn)周圍 （圖2（b））。因此，對(duì)于鉚釘孔的散亂點(diǎn)云，需要進(jìn)行邊界點(diǎn)提取預(yù)處理。

圖2 基于k 鄰域分布的邊界點(diǎn)識(shí)別Fig.2 Boundary point recognition based on k neighborhood distribution

由于散亂點(diǎn)云在數(shù)據(jù)存儲(chǔ)結(jié)構(gòu)中每個(gè)點(diǎn)是無規(guī)律排布的，在這樣的數(shù)據(jù)結(jié)構(gòu)中查找每個(gè)點(diǎn)的k近鄰必須循環(huán)每個(gè)點(diǎn)，效率低下，因此必須建立點(diǎn)云數(shù)據(jù)的空間拓?fù)浣Y(jié)構(gòu)。常見的k近鄰搜索算法有八叉樹法、空間柵格法和kd樹法[14–16]。本文采用kd樹法建立散亂點(diǎn)云的kd樹結(jié)構(gòu)，從而查詢每個(gè)點(diǎn)Pi的k鄰域點(diǎn)集Nj（j=0，1，2，…，k–1）。

以點(diǎn)Pi和其k鄰域點(diǎn)Nj（j=0，1，2，…，k–1）構(gòu)成的點(diǎn)集X來最小二乘擬合得到Pi點(diǎn)的局部切平面。在切平面上建立平面坐標(biāo)系，如圖3所示，以Pi到切平面的投影點(diǎn)Pi′作為坐標(biāo)原點(diǎn)，Pi′點(diǎn)和N0點(diǎn)的投影點(diǎn)N0′構(gòu)成的向量Pi′N0′作為坐標(biāo)系x軸，平面法向與向量Pi′N0′的叉乘×Pi′N0′作為坐標(biāo)系y軸。將點(diǎn)集X中每個(gè)點(diǎn)的三維坐標(biāo)都轉(zhuǎn)換到該平面坐標(biāo)系上，得到點(diǎn)集X的二維坐標(biāo)點(diǎn)集合X′（u，V）。

圖3 建立局部平面坐標(biāo)系Fig.3 Establishing local plane coordinate system

以X′（u，V）中的Pi″點(diǎn)作為向量起點(diǎn)，Nj′點(diǎn)作為向量終點(diǎn)，得到平面向量集計(jì)算中每個(gè)向量到局部坐標(biāo)系的x軸的夾角αj和與y軸的夾角βj。若βj＞π/2，則αj=αj+π。然后將αj以升序進(jìn)行排列得到角度序列ηj，計(jì)算ηj相鄰角度之間的夾角θj=ηj–ηj–1，其中，j∈[1，k–1]，如圖4所示。

圖4 計(jì)算相鄰向量夾角Fig.4 Calculating angle between adjacent vectors

當(dāng)相鄰角度序列θj中的最大夾角θmax超過最大夾角閾值εθ時(shí)，則點(diǎn)Pi為邊界點(diǎn)，否則點(diǎn)Pi為非邊界點(diǎn)，如圖5所示。閾值εθ的大小要視點(diǎn)云分布情況而定，一般將閾值εθ設(shè)置為π /2 可得到良好的識(shí)別結(jié)果。

圖5 邊界點(diǎn)識(shí)別Fig.5 Boundary point recognition

2 鉚釘孔邊界提取

點(diǎn)云預(yù)處理后得到的邊界點(diǎn)中不僅包含有多個(gè)鉚釘孔邊界特征，還包含有點(diǎn)云外邊界特征以及點(diǎn)云漏洞邊界特征，如圖6所示。邊界點(diǎn)云在數(shù)據(jù)結(jié)構(gòu)上仍為散亂點(diǎn)云，需要對(duì)邊界點(diǎn)云建立點(diǎn)云拓?fù)浣Y(jié)構(gòu)，并對(duì)邊界點(diǎn)云進(jìn)行分割，得到屬于不同邊界特征的邊界點(diǎn)云塊，最后根據(jù)點(diǎn)云塊的橢圓度和圓度提取其中的鉚釘孔邊界。

圖6 點(diǎn)云中的邊界特征Fig.6 Boundary features in point clouds

2.1 邊界點(diǎn)歐式聚類分割

邊界點(diǎn)云中的不同邊界特征具有同邊界特征內(nèi)相鄰點(diǎn)距離連續(xù)的特點(diǎn)，即同邊界特征內(nèi)的相鄰點(diǎn)距離間距較小，不同邊界特征的最近點(diǎn)距離較大?；诖朔植继攸c(diǎn)，本文采用歐式距離聚類分割算法[17]對(duì)邊界點(diǎn)云中的邊界特征進(jìn)行分割。其邊界點(diǎn)云聚類分割步驟如下。

步驟1。對(duì)邊界點(diǎn)云X建立kd樹點(diǎn)云結(jié)構(gòu)，方便后續(xù)點(diǎn)鄰域搜索。

步驟2。創(chuàng)建空的聚類集C和點(diǎn)集Q。

步驟3。對(duì)任意點(diǎn)pi∈X，執(zhí)行以下操作：（1）把pi加入Q中 （2）對(duì)每個(gè)點(diǎn)pj∈Q，通過邊界點(diǎn)kd樹的k近鄰搜索算法找到pj的k鄰域，放入點(diǎn)集Pjk；對(duì)每個(gè)pjk∈Pjk，若pjk不在Q中，且pjk與pj的歐式距離rj

步驟4。當(dāng)邊界點(diǎn)云X中每個(gè)點(diǎn)都執(zhí)行過步驟3 后，就得到了聚類集C。

步驟5。刪除聚類集C中的點(diǎn)數(shù)量少于nmin（nmin為最小聚類點(diǎn)數(shù)量閾值）的點(diǎn)云塊，得到最終聚類集C。

2.2 橢圓擬合結(jié)果提取鉚釘孔

邊界點(diǎn)云聚類分割后得到屬于不同邊界特征的點(diǎn)云塊，對(duì)點(diǎn)云塊分別擬合平面并把邊界點(diǎn)投影到各自擬合平面上。由于鉚釘孔邊界點(diǎn)的擬合平面與理論鉚釘孔所在圓柱相交成一橢圓形狀，邊界點(diǎn)在擬合平面上的投影點(diǎn)分布在該橢圓周圍，如圖7所示。因此本文對(duì)所有的邊界點(diǎn)聚類分割點(diǎn)云塊擬合平面橢圓，通過橢圓擬合結(jié)果判斷是否為鉚釘孔邊界點(diǎn)。

圖7 鉚釘孔邊界點(diǎn)橢圓擬合Fig.7 Ellipse fitting of rivet hole boundary points

首先對(duì)分割點(diǎn)云塊ci∈C最小二乘擬合平面，構(gòu)建局部平面二維坐標(biāo)系：以平面上的任一向量? 作為坐標(biāo)系x軸，平面的法向?與? 的叉乘作為坐標(biāo)系y軸，點(diǎn)云塊ci中的任一點(diǎn)p0在平面上的投影點(diǎn)p0′作為坐標(biāo)系原點(diǎn)。將點(diǎn)云塊ci中所有點(diǎn)的三維坐標(biāo)轉(zhuǎn)換到構(gòu)建的局部平面二維坐標(biāo)系下，得到二維坐標(biāo)點(diǎn)集合ci′{pj′（xj′，yj′），j=1，…，mj}。

然后對(duì)平面二維坐標(biāo)點(diǎn)集合ci′最小二乘擬合平面橢圓。設(shè)橢圓方程為Ax2+Bxy+Cy2+Dx+Ey+F=0，并滿足橢圓參數(shù)約束4AC–B2> 0，擬合優(yōu)化目標(biāo)函數(shù)為

令W=[A，B，C，D，E，F(xiàn)]T

則優(yōu)化目標(biāo)為

式中，H=為橢圓參數(shù)約束。

由于||WTX||2時(shí)，W存在縮放因子使得所有W′=aW也滿足優(yōu)化目標(biāo)，因此可令WTHW=1，于是優(yōu)化目標(biāo)就變?yōu)?/p>

構(gòu)造拉格朗日函數(shù)，即

對(duì)式（4）進(jìn)行求解。令S=XXT，則SW=λHW，通過求解廣義逆矩陣得到6 個(gè)可能的解W。由于S為正定矩陣，所以λ> 0，因此可以用WTHW=1 和λ> 0 這兩個(gè)條件來篩選得到最終合格的解W=[A，B，C，D，E，F(xiàn)]T。

根據(jù)橢圓參數(shù)W=[A，B，C，D，E，F(xiàn)]T計(jì)算橢圓圓心Oe（x0，y0）和長(zhǎng)短軸半徑ae、be：

計(jì)算橢圓擬合誤差Error 與長(zhǎng)短軸之比Ratio：

其中，F(xiàn)（p′j）=（Ax′j2+Bxj′yj′+Cy′j2+Dxj′+Eyj′+F），n為 點(diǎn)云塊ci∈C中點(diǎn)的數(shù)量。

若Error >εe（εe為擬合誤差閾值），則為該分割點(diǎn)云塊橢圓度低，為非橢圓，即非鉚釘孔；若Error <εe，則比較長(zhǎng)短軸之比。

若1–εr< Ratio < 1+εr（εr為長(zhǎng)短軸比閾值），則該分割點(diǎn)云塊圓度高，即為鉚釘孔邊界，否則為非鉚釘孔邊界。

若對(duì)鉚釘孔半徑大小有要求，可設(shè)置最大半徑閾值rmax和最小半徑閾值rmin，將長(zhǎng)短軸的均值與rmax和rmin比較來提取滿足半徑要求的鉚釘孔。

點(diǎn)云中可能有其他不是鉚釘孔的圓孔，需要對(duì)其進(jìn)行剔除。計(jì)算非鉚釘孔理論圓心與該算法擬合橢圓圓心的空間距離dc，若dc小于一定閾值εd，則剔除該提取孔。

3 鉚釘孔參數(shù)提取

鉚釘孔孔位孔徑參數(shù)提取采用先投影孔邊界點(diǎn)到孔定位面然后最小二乘擬合端面圓的方法。由于孔的邊界點(diǎn)中部分點(diǎn)在孔內(nèi)壁上，使用孔邊界點(diǎn)直徑求解孔定位面會(huì)降低參數(shù)提取精度，所以需要尋找孔周圍一定范圍內(nèi)的點(diǎn)來擬合孔定位面。然后根據(jù)邊界點(diǎn)擬合橢圓和孔定位面求得初始鉚釘孔參數(shù)，根據(jù)邊界點(diǎn)到初始鉚釘孔參數(shù)的距離偏差剔除邊界點(diǎn)中的噪聲點(diǎn)，最后再將孔邊界點(diǎn)投影到孔定位面上擬合圓孔。

3.1 孔定位面擬合

在鉚釘孔邊界提取過程中求得了所有邊界點(diǎn)云塊的橢圓擬合方程、橢圓圓心Oe和長(zhǎng)短軸ae、be，以擬合橢圓的圓心Oe作為搜索中心，長(zhǎng)短軸均值的ω倍作為搜索半徑Rs，通過建立的散亂點(diǎn)云kd樹結(jié)構(gòu)進(jìn)行半徑鄰域搜索找出鉚釘孔周圍Rs范圍內(nèi)的點(diǎn)集Pr，如圖8所示。

圖8 鉚釘孔鄰域點(diǎn)云搜索Fig.8 Point cloud search of rivet hole neighborhood

將點(diǎn)集Pr精確擬合平面。設(shè)孔定位面方程為ax+by+cz+d=0，構(gòu)建擬合目標(biāo)函數(shù)：

用特征值法[18]求解式（8），得到平面參數(shù)a、b、c、d。然后計(jì)算點(diǎn)集中每個(gè)點(diǎn)pi∈Pr到平面的距離di，當(dāng)di>εd（εd為平面擬合誤差閾值）時(shí)，則認(rèn)為該點(diǎn)是噪點(diǎn)并剔除，反之保留，最后使用保留的點(diǎn)重新擬合孔定位面。

3.2 端面圓孔擬合

在孔定位面上建立平面二維坐標(biāo)系，將鉚釘孔邊界點(diǎn)的三維坐標(biāo)轉(zhuǎn)換成平面二維坐標(biāo)點(diǎn)pj（xj，yj）。然后根據(jù)擬合橢圓和定位面計(jì)算初始鉚釘孔參數(shù)，將橢圓圓心投影到孔定位面上，投影點(diǎn)即為鉚釘孔圓心，橢圓短軸長(zhǎng)度即為鉚釘孔半徑。

計(jì)算平面二維坐標(biāo)點(diǎn)pj（xj，yj）到初始鉚釘孔的距離δj，當(dāng)δj>εc（εc為圓孔擬合誤差閾值）時(shí)，則認(rèn)為該點(diǎn)為噪點(diǎn)并剔除，反之保留，最后使用保留的點(diǎn)重新擬合端面圓。

設(shè)端面圓孔方程為（x–A）2+（y–B）2=R2，令a=–2A，b=–2B，c=A2+B2–R2，則方程變?yōu)閤2+y2+ax+by+c=0，構(gòu)建擬合目標(biāo)函數(shù)：

式（9）中分別對(duì)a、b、c求偏導(dǎo)，令偏導(dǎo)等于0，可得：

解式（10）可得端面圓參數(shù)a、b、c，則端面圓的二維圓心坐標(biāo)為（x0，y0）=（a/–2，b/–2），半徑

把二維圓心（x0，y0）坐標(biāo)逆變換到三維空間中，得到三維空間圓心坐標(biāo)（x′0，y′0，z′0），三維空間圓半徑R等于平面圓半徑r，孔法向（i，j，k）等于擬合平面的法向，鉚釘孔參數(shù)提取完成。

4 試驗(yàn)驗(yàn)證

4.1 提取效果驗(yàn)證

應(yīng)用本文方法對(duì)散亂點(diǎn)云中的鉚釘孔特征進(jìn)行提取，驗(yàn)證其提取效果，如圖9所示。圖9（a）為結(jié)構(gòu)光掃描儀測(cè)得的某壁板零件的散亂點(diǎn)云數(shù)據(jù)，所含點(diǎn)數(shù)為90688；用本文邊界點(diǎn)提取算法提取點(diǎn)云中的邊界點(diǎn)，k=50、εθ=π/2 時(shí)如圖9（b）所示，提取邊界點(diǎn)數(shù)為3392；用本文邊界點(diǎn)云歐氏距離聚類分割算法分割邊界點(diǎn)，dth=2mm、nmin=10 時(shí)分割結(jié)果如圖9（c）所示，共分割出14 個(gè)邊界特征；依據(jù)橢圓擬合結(jié)果提取鉚釘孔邊界特征，εe=1、εr=0.05 時(shí)提取結(jié)果如圖9（d）所示，14 個(gè)邊界特征中的9 個(gè)鉚釘孔特征全部提取成功?？梢钥闯觯疚姆椒梢杂行崛∩y點(diǎn)云中的鉚釘孔邊界特征。

圖9 鉚釘孔特征提取效果Fig.9 Results of rivet hole feature extraction

4.2 提取精度驗(yàn)證

為驗(yàn)證本文提出算法的鉚釘孔參數(shù)提取精度，采用天遠(yuǎn)Robot-Scan 掃描儀對(duì)鉚釘孔試驗(yàn)件進(jìn)行測(cè)量和數(shù)據(jù)處理，Robot-Scan 相機(jī)像素為1000 萬像素，平均單張點(diǎn)間距為0.16mm，單張測(cè)量精度為15μm。

如圖10所示，使用Robot-Scan 掃描儀對(duì)試驗(yàn)件上的10 個(gè)理論直徑d=8mm 的鉚釘孔進(jìn)行測(cè)量得到試驗(yàn)件散亂點(diǎn)云數(shù)據(jù)。應(yīng)用本文方法對(duì)數(shù)據(jù)進(jìn)行處理，得到鉚釘孔參數(shù)。然后將三坐標(biāo)測(cè)量機(jī)測(cè)量試驗(yàn)件得到的參數(shù)作為真值，與本文方法提取結(jié)果進(jìn)行比較得出本文方法鉚釘孔參數(shù)提取偏差。三坐標(biāo)測(cè)量機(jī)測(cè)量結(jié)果和本文參數(shù)提取結(jié)果如表1所示，鉚釘提取偏差如圖11所示?？梢钥闯?，本文提出方法對(duì)試驗(yàn)件上鉚釘孔的孔位提取偏差均值為0.017mm，最大值為0.029mm；鉚釘孔孔徑提取偏差均值為0.0797mm，最大值為0.095mm。

圖10 使用Robot–Scan 對(duì)鉚釘孔試驗(yàn)件測(cè)量Fig.10 Measurement of rivet hole test pieces with Robot–Scan

表1 鉚釘孔參數(shù)提取試驗(yàn)結(jié)果Table 1 Results of rivet hole parameter extraction experiment

圖11 鉚釘孔提取偏差Fig.11 Extraction deviation of rivet holes

從提取結(jié)果分析可以看出該方法孔位提取精度較高，可達(dá)0.029mm；鉚釘孔孔徑提取誤差相比較孔位提取誤差較大，這是由于測(cè)量設(shè)備測(cè)量得到的散亂點(diǎn)云具有一定的點(diǎn)間距，孔邊界測(cè)量點(diǎn)不一定全部落在孔內(nèi)徑上，從而使得該方法對(duì)鉚釘孔的直徑提取精度相對(duì)于孔位提取精度較低。

5 結(jié)論

（1）本文提出了一種基于散亂點(diǎn)云的鉚釘孔孔位參數(shù)提取方法。用該方法對(duì)試驗(yàn)件測(cè)量散亂點(diǎn)云數(shù)據(jù)進(jìn)行了處理。試驗(yàn)表明，該方法能夠精確提取散亂點(diǎn)云中的全部鉚釘孔特征，且孔位提取精度高，提取精度可達(dá)0.029mm。

（2）使用本文方法需要注意，由于測(cè)量點(diǎn)不能夠保證全部落在鉚釘孔徑上，導(dǎo)致了鉚釘孔孔徑提取誤差相對(duì)較大。因此，為保證有較高的孔徑提取精度，需要測(cè)量設(shè)備精度和點(diǎn)測(cè)量密度高，且盡量使測(cè)量點(diǎn)能夠較多地落在孔內(nèi)徑上。