田清廉,熊天辰,黃 翔,李瀧杲,郝 龍
(1.南京航空航天大學(xué)機(jī)電學(xué)院,南京 210016;2.江西航天海虹測(cè)控技術(shù)有限責(zé)任公司,南昌 330024)
我國(guó)航空事業(yè)不斷發(fā)展,對(duì)零件上鉚釘孔的制孔精度要求也越來(lái)越高,而飛機(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è)備也開(kāi)始更多地應(yīng)用在飛機(jī)零件鉚釘孔檢測(cè)上,這種檢測(cè)方式具有精度高、效率高、柔性好等特點(diǎn)[3–4]。
目前,在非接觸式圓孔測(cè)量方面已有了很多的研究與成果,主要有基于視覺(jué)圖像的測(cè)量方法和基于三維點(diǎn)云的測(cè)量方法?;谝曈X(jué)圖像的測(cè)量方法主要是通過(guò)圖像處理中的Canny 算子等邊界識(shí)別算法提取照相測(cè)量圖像中的邊界特征進(jìn)行圓孔參數(shù)擬合[5–6],但該方法存在空間圓透視投影畸變的問(wèn)題[7],影響檢測(cè)精度。基于三維點(diǎn)云的測(cè)量方法是使用三維激光掃描等測(cè)量設(shè)備獲取零件點(diǎn)云,再提取點(diǎn)云中的孔邊界特征進(jìn)行圓孔參數(shù)擬合。針對(duì)點(diǎn)云中邊界特征提取,Jenke 等[8]提出了一種將點(diǎn)云網(wǎng)格化再通過(guò)網(wǎng)格的拓?fù)浣Y(jié)構(gòu)識(shí)別邊界特征的方法,這種方法可以有效提取點(diǎn)云孔洞邊界特征,但需要網(wǎng)格化整個(gè)點(diǎn)云,數(shù)據(jù)量大時(shí)算法效率低。因此出現(xiàn)了直接對(duì)三維點(diǎn)云提取孔邊界特征的方法,該方法通過(guò)建立點(diǎn)云空間拓?fù)潢P(guān)系,并設(shè)定某一特定判斷機(jī)制用于識(shí)別邊界特征點(diǎn)。Milroy[9]和Yang[10]等通過(guò)求解曲率極值提取點(diǎn)云邊界,這種方法可以有效地提取比較光滑平坦的數(shù)據(jù)模型的邊界點(diǎn),不太適用于曲率變化大的數(shù)據(jù)模型。張杰等[11]提出了一種基于T–Scan測(cè)量的薄壁鈑金件孔特征重構(gòu)方法,通過(guò)掃描線點(diǎn)云結(jié)構(gòu)中的孔截?cái)嗑€特征識(shí)別平面點(diǎn)和孔徑點(diǎn),擬合求解平面和空間圓,但該方法僅限于掃描線點(diǎn)云類型且提取精度受圓孔直徑影響。劉增藝等[12]提出了一種曲面孔位視覺(jué)測(cè)量技術(shù),使用手持式雙目視覺(jué)線結(jié)構(gòu)光測(cè)量曲面得到曲面點(diǎn)云,然后通過(guò)交互式方法提取散亂點(diǎn)云圓孔邊緣點(diǎn),進(jìn)行圓孔擬合,由于該方法在圓孔識(shí)別過(guò)程中需要人工預(yù)先分割出圓孔所在區(qū)域,提取效率低且不能滿足多孔自動(dòng)提取需求。
基于以上分析,為解決目前面向鉚釘孔檢測(cè)領(lǐng)域的測(cè)量點(diǎn)云中的孔位與孔徑參數(shù)提取方法限制多和人工交互較多的問(wèn)題,本文提出了一種基于散亂點(diǎn)云的鉚釘孔孔位與孔徑提取方法,該方法可以自動(dòng)識(shí)別散亂點(diǎn)云中的全部鉚釘孔特征,并精確提取孔位和孔徑參數(shù)。采用本文方法對(duì)試驗(yàn)件散亂點(diǎn)云進(jìn)行數(shù)據(jù)處理,驗(yàn)證該方法的實(shí)用性與精度。
各種數(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)無(wú)明顯的幾何分布特征,呈散亂無(wú)序狀態(tài),測(cè)量點(diǎn)是隨機(jī)生成的。由于散亂點(diǎn)云無(wú)序和隨機(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)云分布情況,本文通過(guò)點(diǎn)的k鄰域點(diǎn)幾何分布來(lái)判斷是否為邊界點(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)是無(wú)規(guī)律排布的,在這樣的數(shù)據(jù)結(jié)構(gòu)中查找每個(gè)點(diǎn)的k近鄰必須循環(huán)每個(gè)點(diǎn),效率低下,因此必須建立點(diǎn)云數(shù)據(jù)的空間拓?fù)浣Y(jié)構(gòu)。常見(jiàn)的k近鄰搜索算法有八叉樹(shù)法、空間柵格法和kd樹(shù)法[14–16]。本文采用kd樹(shù)法建立散亂點(diǎn)云的kd樹(shù)結(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來(lái)最小二乘擬合得到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超過(guò)最大夾角閾值εθ時(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
點(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
邊界點(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樹(shù)點(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,通過(guò)邊界點(diǎn)kd樹(shù)的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í)行過(guò)步驟3 后,就得到了聚類集C。 步驟5。刪除聚類集C中的點(diǎn)數(shù)量少于nmin(nmin為最小聚類點(diǎn)數(shù)量閾值)的點(diǎn)云塊,得到最終聚類集C。 邊界點(diǎn)云聚類分割后得到屬于不同邊界特征的點(diǎn)云塊,對(duì)點(diǎn)云塊分別擬合平面并把邊界點(diǎn)投影到各自擬合平面上。由于鉚釘孔邊界點(diǎn)的擬合平面與理論鉚釘孔所在圓柱相交成一橢圓形狀,邊界點(diǎn)在擬合平面上的投影點(diǎn)分布在該橢圓周圍,如圖7所示。因此本文對(duì)所有的邊界點(diǎn)聚類分割點(diǎn)云塊擬合平面橢圓,通過(guò)橢圓擬合結(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,通過(guò)求解廣義逆矩陣得到6 個(gè)可能的解W。由于S為正定矩陣,所以λ> 0,因此可以用WTHW=1 和λ> 0 這兩個(gè)條件來(lái)篩選得到最終合格的解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比較來(lái)提取滿足半徑要求的鉚釘孔。 點(diǎn)云中可能有其他不是鉚釘孔的圓孔,需要對(duì)其進(jìn)行剔除。計(jì)算非鉚釘孔理論圓心與該算法擬合橢圓圓心的空間距離dc,若dc小于一定閾值εd,則剔除該提取孔。 鉚釘孔孔位孔徑參數(shù)提取采用先投影孔邊界點(diǎn)到孔定位面然后最小二乘擬合端面圓的方法。由于孔的邊界點(diǎn)中部分點(diǎn)在孔內(nèi)壁上,使用孔邊界點(diǎn)直徑求解孔定位面會(huì)降低參數(shù)提取精度,所以需要尋找孔周圍一定范圍內(nèi)的點(diǎn)來(lái)擬合孔定位面。然后根據(jù)邊界點(diǎn)擬合橢圓和孔定位面求得初始鉚釘孔參數(shù),根據(jù)邊界點(diǎn)到初始鉚釘孔參數(shù)的距離偏差剔除邊界點(diǎn)中的噪聲點(diǎn),最后再將孔邊界點(diǎn)投影到孔定位面上擬合圓孔。 在鉚釘孔邊界提取過(guò)程中求得了所有邊界點(diǎn)云塊的橢圓擬合方程、橢圓圓心Oe和長(zhǎng)短軸ae、be,以擬合橢圓的圓心Oe作為搜索中心,長(zhǎng)短軸均值的ω倍作為搜索半徑Rs,通過(guò)建立的散亂點(diǎn)云kd樹(shù)結(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)重新擬合孔定位面。 在孔定位面上建立平面二維坐標(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ù)提取完成。 應(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è)鉚釘孔特征全部提取成功??梢钥闯?,本文方法可以有效提取散亂點(diǎn)云中的鉚釘孔邊界特征。 圖9 鉚釘孔特征提取效果Fig.9 Results of rivet hole feature extraction 為驗(yàn)證本文提出算法的鉚釘孔參數(shù)提取精度,采用天遠(yuǎn)Robot-Scan 掃描儀對(duì)鉚釘孔試驗(yàn)件進(jìn)行測(cè)量和數(shù)據(jù)處理,Robot-Scan 相機(jī)像素為1000 萬(wàn)像素,平均單張點(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ì)于孔位提取精度較低。 (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)徑上。2.2 橢圓擬合結(jié)果提取鉚釘孔
3 鉚釘孔參數(shù)提取
3.1 孔定位面擬合
3.2 端面圓孔擬合
4 試驗(yàn)驗(yàn)證
4.1 提取效果驗(yàn)證
4.2 提取精度驗(yàn)證
5 結(jié)論