朱秀敏，黃 磊

(南京林業(yè)大學(xué)機(jī)械電子工程學(xué)院自動(dòng)化系，江蘇南京 210037)

0 引言

貴重零件表面缺陷的修復(fù)在工業(yè)領(lǐng)域有著廣泛的需求。隨著3D打印技術(shù)的快速發(fā)展，如果能夠獲得零件表面缺陷的準(zhǔn)確3D尺寸，可以利用金屬3D打印技術(shù)對(duì)零件表面的缺損部位進(jìn)行填補(bǔ)和修復(fù)，以提高零件的使用壽命。這一修補(bǔ)方法正在零件修復(fù)領(lǐng)域得到越來越多的應(yīng)用。因此對(duì)零件表面缺陷3D尺寸的精確測(cè)量顯得尤為重要。

目前工業(yè)上常用的零件缺陷檢測(cè)方式主要包括：磁粉檢測(cè)法、滲透檢測(cè)法、渦流檢測(cè)法、超聲波檢測(cè)法、X射線檢測(cè)和機(jī)器視覺檢測(cè)法等[1]。前3種方法容易受材料種類、零件形狀限制，而X射線對(duì)人體有害，超聲波檢測(cè)不適合結(jié)構(gòu)復(fù)雜的工件，機(jī)器視覺檢測(cè)法難以提取出精確的零件缺陷3D尺寸。本文使用3D激光掃描儀采集缺陷零件的表面點(diǎn)云，與CATIA導(dǎo)出的完整零件點(diǎn)云作對(duì)比，通過濾波、SAC-IA粗配準(zhǔn)、Huber損失函數(shù)優(yōu)化的ICP精確配準(zhǔn)、最鄰近查找等操作提取出缺陷部位點(diǎn)云。實(shí)驗(yàn)結(jié)果表明：具有不受零件形狀、缺陷形狀限制，適用范圍廣、自動(dòng)化檢測(cè)、精度高等優(yōu)點(diǎn)。

1 點(diǎn)云數(shù)據(jù)獲取及預(yù)處理

1.1 3D掃描儀原理

本文使用PhoXi 3D掃描儀對(duì)受損零件進(jìn)行測(cè)繪，得到其三維點(diǎn)云數(shù)據(jù)。掃描儀實(shí)物如圖1所示。

圖1 PhoXi 3D掃描儀示意圖

PhoXi 3D激光掃描儀采用波長(zhǎng)639 nm的紅色可見光為光源，CMOS傳感器成像，兼具3D激光測(cè)繪和2D CMOS相機(jī)功能，無需用戶自己進(jìn)行激光點(diǎn)云坐標(biāo)系和相機(jī)坐標(biāo)系的對(duì)齊，可直接輸出物體的彩色點(diǎn)云。它的一幀輸出點(diǎn)云可達(dá)320萬個(gè)點(diǎn)，對(duì)應(yīng)的分辨率達(dá)2 064像素×1 544像素，測(cè)距距離2 m，并可按照設(shè)置的最大測(cè)距誤差門限(最小1 mm)輸出所有符合這一要求的點(diǎn)[2-3]。生成的點(diǎn)云坐標(biāo)系如圖2所示，原點(diǎn)一般位于標(biāo)記板右上角的第一個(gè)圓圈中。

圖2 標(biāo)記坐標(biāo)系示意圖

基于PhoXi 3D激光掃描儀的零件缺陷檢測(cè)系統(tǒng)如圖3所示，將零件置于標(biāo)記板上，底部工作臺(tái)旋轉(zhuǎn)一圈，共掃描零件的8個(gè)不同的角度，形成8個(gè)點(diǎn)云文件。掃描儀采集到的數(shù)據(jù)包含零件點(diǎn)云的三維坐標(biāo)、顏色等信息。

圖3 基于PhoXi 3D激光掃描儀的零件缺陷檢測(cè)系統(tǒng)

1.2 數(shù)據(jù)預(yù)處理

掃描儀將8個(gè)角度的點(diǎn)云數(shù)據(jù)生成8個(gè)點(diǎn)云文件，首先提取出這8個(gè)點(diǎn)云文件中的三維坐標(biāo)，然后將這8個(gè)角度的點(diǎn)云合成零件點(diǎn)云。由于掃描儀在掃描零件時(shí)，會(huì)將標(biāo)記板、工作臺(tái)等物體一并掃入儀器中，因此還需要對(duì)采集的數(shù)據(jù)進(jìn)行預(yù)處理以獲得零件的點(diǎn)云，處理流程如圖4所示。

圖4 原始點(diǎn)云預(yù)處理流程圖

統(tǒng)計(jì)濾波能將原始點(diǎn)云中稀疏的離群點(diǎn)剔除，如測(cè)量時(shí)產(chǎn)生的噪聲、離零件較遠(yuǎn)且與零件無關(guān)的雜物等。統(tǒng)計(jì)濾波算法設(shè)想，假設(shè)一個(gè)含有離散點(diǎn)的點(diǎn)云，其濾除離散點(diǎn)后的點(diǎn)云中所有點(diǎn)到其各自鄰域的平均距離服從高斯分布，由于離散點(diǎn)相對(duì)于主體點(diǎn)云是少量的，故可認(rèn)為原始點(diǎn)云與濾除離散點(diǎn)后的點(diǎn)云的均值與標(biāo)準(zhǔn)差相同[4]。統(tǒng)計(jì)濾波的流程如下：

(1)對(duì)原始點(diǎn)云數(shù)據(jù)集P={p1,p2,…,pn}進(jìn)行統(tǒng)計(jì)濾波，即對(duì)其中每一個(gè)點(diǎn)的鄰域進(jìn)行統(tǒng)計(jì)分析，計(jì)算點(diǎn)pi到最近的k個(gè)點(diǎn)的平均距離，記作Si，則由n個(gè)點(diǎn)的平均距離組成的點(diǎn)集S服從高斯分布。

(2)根據(jù)式(1)、式(2)分別計(jì)算點(diǎn)集S的均值μ、標(biāo)準(zhǔn)差σ。

(1)

(2)

(3)判斷點(diǎn)pi與其鄰域的平均距離是否在(μ-std·σ，μ+std·σ)范圍內(nèi)。如果在，則保留；如果不在，則作為離群點(diǎn)剔除。

算法中的k、std分別為鄰域點(diǎn)的個(gè)數(shù)、標(biāo)準(zhǔn)差倍數(shù)，2個(gè)參數(shù)由人為設(shè)定，前者控制鄰域大小，后者控制篩選的苛刻程度，適當(dāng)?shù)膮?shù)才能得到高質(zhì)量的濾波結(jié)果[4]。

經(jīng)過統(tǒng)計(jì)濾波后的點(diǎn)云仍有標(biāo)記板、工作臺(tái)2個(gè)不需要的平面數(shù)據(jù)，再進(jìn)行RANSAC平面分割，則能夠去除這2個(gè)平面。平面分割的步驟如下：

(1)從點(diǎn)云中隨機(jī)抽取3個(gè)點(diǎn)，計(jì)算這3個(gè)點(diǎn)組成的平面表達(dá)式。

(2)計(jì)算其他所有點(diǎn)到該平面的距離，如果小于閾值d則認(rèn)為是同一個(gè)平面的點(diǎn)，添加到該平面的點(diǎn)集中。

(3)如果同一個(gè)平面的點(diǎn)的個(gè)數(shù)超過N則保留這個(gè)平面并停止迭代，否則重復(fù)步驟(1)(2)(3)，直到找到迭代次數(shù)達(dá)到閾值。最終可以得到一個(gè)完整平面的點(diǎn)集。

運(yùn)用2次RANSAC地面分割后，標(biāo)記板與工作臺(tái)2個(gè)平面的點(diǎn)云均可被剔除。剩下的部位則為零件以及周圍噪聲的點(diǎn)云，最后采用半徑濾波即可去除零件周圍的噪聲。算法思路如下：判斷每一個(gè)點(diǎn)半徑為r的鄰域內(nèi)的點(diǎn)的個(gè)數(shù)是否小于閾值，若小于，則認(rèn)為該點(diǎn)為噪聲點(diǎn)并去除，否則保留。

2 點(diǎn)云配準(zhǔn)

點(diǎn)云配準(zhǔn)的實(shí)質(zhì)是尋找一個(gè)最優(yōu)的變換矩陣，使得源點(diǎn)云經(jīng)過該變換能夠和目標(biāo)點(diǎn)云的空間位置差別最小化。本文采用采樣一致性初始配準(zhǔn)(sample consensus initial alignment，SAC-IA)算法進(jìn)行粗配準(zhǔn)，Huber損失函數(shù)優(yōu)化的迭代最近點(diǎn)算法(iterative closest Point，ICP)進(jìn)行精確配準(zhǔn)，配準(zhǔn)流程如圖5所示。

圖5 缺陷點(diǎn)云配準(zhǔn)流程圖

首先對(duì)零件點(diǎn)云數(shù)據(jù)進(jìn)行預(yù)處理，對(duì)缺陷點(diǎn)云和完整點(diǎn)云使用體素下采樣來提高算法的效率；用SAC-IA算法迭代得到初始變換矩陣，實(shí)現(xiàn)初步配準(zhǔn)；最后利用Huber損失函數(shù)優(yōu)化的ICP算法迭代得到最優(yōu)變換矩陣，實(shí)現(xiàn)精確配準(zhǔn)。

2.1 點(diǎn)云粗配準(zhǔn)

初值的選取對(duì)于ICP算法十分重要，如果沒有合適的初值，會(huì)導(dǎo)致ICP配準(zhǔn)出現(xiàn)局部最優(yōu)的情況，很難迭代出全局最優(yōu)變換矩陣[5]。因此，在進(jìn)行精確配準(zhǔn)之前，先進(jìn)行粗配準(zhǔn)。由于零件的形狀、缺陷位置對(duì)于點(diǎn)云質(zhì)心、主成分影響較大，速度更快的PCA粗配準(zhǔn)算法對(duì)于零件表面缺陷檢測(cè)不合適，因此本文采用SAC-IA算法，即采樣一致性初始配準(zhǔn)算法[6]。該算法依賴點(diǎn)特征直方圖(point feature histograms，PFH)，使用類似RANSAC算法的思路進(jìn)行配準(zhǔn)[7]。由于計(jì)算PFH的特征描述子速度較慢，考慮代碼效率，首先需要對(duì)原數(shù)據(jù)進(jìn)行降采樣，再使用快速點(diǎn)云直方圖(FPFH)特征描述子。

SAC-IA粗配準(zhǔn)過程如下：

(1)分別計(jì)算缺陷點(diǎn)云P和完整點(diǎn)云Q的FPFH特征描述子的值。每次從缺陷點(diǎn)云中隨機(jī)選取n(n≥3)個(gè)點(diǎn)作為一個(gè)子集Pi={p1,p2,…,pn}，這n個(gè)點(diǎn)兩兩之間距離大于最小閾值[8]。

(2)在FPFH特征空間中進(jìn)行最近鄰查找，檢測(cè)在Q中FPFH描述子相似的點(diǎn)，從相似點(diǎn)中隨機(jī)選出1個(gè)點(diǎn)作為缺陷點(diǎn)云在完整點(diǎn)云中的對(duì)應(yīng)點(diǎn)，這些相似的點(diǎn)組成對(duì)應(yīng)點(diǎn)集Qi={q1,q2,…,qn}。

(3)通過奇異值矩陣分解法(SVD)求解Pi、Qi對(duì)應(yīng)點(diǎn)之間的變換矩陣Ti,元素構(gòu)成如式(3) 所示，通過求解對(duì)應(yīng)點(diǎn)變換后的距離誤差和來判斷當(dāng)前的配準(zhǔn)效果[9]。如果不滿足閾值要求，則重復(fù)步驟(2)(3)，直到迭代出誤差最小且滿足閾值要求的變換矩陣。距離誤差和函數(shù)常用Huber懲罰函數(shù)表示[10]，Huber函數(shù)如式(4)所示：

(3)

(4)

la=‖qa-(Rpa-t)‖2

(5)

式中：Ti為點(diǎn)集Pi到Qi的變換矩陣；Ri、ti分別為變換矩陣中的旋轉(zhuǎn)、平移矩陣；ml為預(yù)先給定值；la為第a組對(duì)應(yīng)點(diǎn)變換之后的距離差；pa、qa分別為點(diǎn)集Pi、Qi中的點(diǎn)。

2.2 點(diǎn)云精確配準(zhǔn)

經(jīng)過粗配準(zhǔn)后的2組點(diǎn)云已經(jīng)大致重合，但仍存在偏差，為了更好地提取缺陷部位，還需要進(jìn)一步精確配準(zhǔn)。傳統(tǒng)的ICP算法以源點(diǎn)云與目標(biāo)點(diǎn)云對(duì)應(yīng)點(diǎn)之間的距離作為配準(zhǔn)結(jié)果的評(píng)估準(zhǔn)則，容易受錯(cuò)誤的對(duì)應(yīng)點(diǎn)以及噪聲的干擾，導(dǎo)致配準(zhǔn)精度降低[11]。由于缺陷點(diǎn)云的缺陷部位對(duì)于完整點(diǎn)云來說相當(dāng)于是外點(diǎn)，并且由掃描儀得到的數(shù)據(jù)有噪聲干擾，于是本文并不采用傳統(tǒng)點(diǎn)到點(diǎn)的ICP。而點(diǎn)到面的ICP相對(duì)于點(diǎn)到點(diǎn)的ICP能更好的對(duì)錯(cuò)誤的對(duì)應(yīng)點(diǎn)、噪聲魯棒，且迭代收斂速度更快[12]。故這里采用點(diǎn)到面的ICP，建立kd-tree加速最鄰近查找，并使用Huber損失函數(shù)優(yōu)化，以提高配準(zhǔn)的魯棒性[13-14]。

ICP算法包含離群點(diǎn)去除階段與誤差最小化階段，2個(gè)階段可以概括為通過最小化目標(biāo)函數(shù)式(6)來估計(jì)轉(zhuǎn)換矩陣T[15]。

(6)

式中：ρ(·)為成本函數(shù)；pi、qj分別為缺陷點(diǎn)云P以及完整點(diǎn)云Q中的點(diǎn)；nj為qj的法向量；T為變換矩陣。

為了簡(jiǎn)化式(6)中的雙重求和，通常采用對(duì)應(yīng)點(diǎn)的子集來近似，這里使用每個(gè)pi的最鄰近點(diǎn)來近似。將成本函數(shù)ρ(·)分解為權(quán)重ω和平方誤差項(xiàng)來迭代求解最優(yōu)變換矩陣T[16-17]。每次迭代都會(huì)將先前的變換分配給最后估計(jì)的變換，直到收斂。式(6)簡(jiǎn)化后為：

(7)

Lm(T)=min‖(Tpm-qm)·nm‖2

(8)

基于Huber損失函數(shù)優(yōu)化的ICP本質(zhì)是對(duì)式(7)中的ω函數(shù)的改進(jìn)，式(9)即為Huber對(duì)ω的實(shí)現(xiàn)[18]。

(9)

式中k為設(shè)定的超參數(shù)。

算法大致思路如下：

(1)令粗配準(zhǔn)后的轉(zhuǎn)換矩陣為初始矩陣。

(2)缺陷點(diǎn)云經(jīng)過變換矩陣旋轉(zhuǎn)平移形成新的點(diǎn)云。

(3)建立kd-tree，通過最鄰近查找找出新的點(diǎn)云子集在完整點(diǎn)云中的對(duì)應(yīng)點(diǎn)子集，再根據(jù)式(7)～式(9)去除離群點(diǎn)并計(jì)算出新的點(diǎn)云與完整點(diǎn)云之間的最優(yōu)變換矩陣T以及2組點(diǎn)云之間的最小誤差Lm(T)。

(4)判斷誤差是否小于閾值，如果是則保留配準(zhǔn)后的點(diǎn)云，否則重復(fù)步驟(2)、(3)，直到迭代出符合要求的最優(yōu)變換矩陣。

3 缺陷部位提取

缺陷點(diǎn)云與完整點(diǎn)云配準(zhǔn)后，可進(jìn)一步提取缺陷部位，流程如圖6所示?？紤]到點(diǎn)集中可能仍然有噪聲，于是在最后加入統(tǒng)計(jì)濾波器，以提取出完整的缺陷部位。

圖6 缺陷部位提取流程圖

4 實(shí)驗(yàn)結(jié)果

為了驗(yàn)證本文算法對(duì)于零件缺陷檢測(cè)的合理性，使用Python語言和開源點(diǎn)云工具庫Open3D 0.13.0編寫算法進(jìn)行試驗(yàn)。算法運(yùn)行平臺(tái)為CPU頻率2.8 GHz，內(nèi)存8 GB，運(yùn)行Win10的計(jì)算機(jī)。為了驗(yàn)證算法的可行性，本次實(shí)驗(yàn)運(yùn)用了2組零件點(diǎn)云數(shù)據(jù)，缺陷零件的數(shù)據(jù)為真實(shí)測(cè)量數(shù)據(jù)，完整零件的數(shù)據(jù)為CAD模型導(dǎo)出數(shù)據(jù)，其中缺陷的部位分別為破洞、凹陷，此次實(shí)驗(yàn)中分別命名為零件1、零件2。

首先是數(shù)據(jù)的預(yù)處理，以零件1為例，PhoXi 3D掃描儀采集原始數(shù)據(jù)為零件8個(gè)角度的點(diǎn)云數(shù)據(jù)，這8個(gè)角度去重、融合后的點(diǎn)云合成如圖7(a)所示。經(jīng)過統(tǒng)計(jì)濾波、2次RANSAC地面分割和半徑濾波后提取出完整零件點(diǎn)云過程如圖7(b)～圖7(d)所示。

圖7 數(shù)據(jù)預(yù)處理結(jié)果

由表1可知，待檢測(cè)的原零件點(diǎn)云規(guī)模較大，影響配準(zhǔn)效率。故需要對(duì)原點(diǎn)云進(jìn)行降采樣，這里采用體素下采樣。使用降采樣后的數(shù)據(jù)進(jìn)行全局配準(zhǔn)，再使用原數(shù)據(jù)進(jìn)行精確配準(zhǔn)，以獲得契合點(diǎn)云，為下一步缺陷部位提取做準(zhǔn)備。

表1 缺陷點(diǎn)云和完整點(diǎn)云規(guī)模統(tǒng)計(jì)

表2中的配準(zhǔn)誤差采用均方根誤差(root mean sq-uare error，RMSE)，從表2可知：SAC-IA速度雖然不如PCA，但誤差遠(yuǎn)小于PCA。故本文采用SAC-IA+改進(jìn)的ICP算法來配準(zhǔn)，配準(zhǔn)結(jié)果如圖8所示，深色的為缺陷零件點(diǎn)云，淺色的為完整零件點(diǎn)云?？梢钥闯觯?jīng)過粗配準(zhǔn)后，2個(gè)點(diǎn)云已大致重合；再經(jīng)過精確配準(zhǔn)，2個(gè)點(diǎn)云完全重合。

表2 零件1的各配準(zhǔn)算法結(jié)果比較

采用kd-tree加速的最鄰近查找提取缺陷部位點(diǎn)云如圖9所示，可以看出該算法能夠較好的提取零件表面缺陷部位的3D點(diǎn)云，并且通過該數(shù)據(jù)，可以進(jìn)一步對(duì)零件做出3D打印修復(fù)。

圖9 缺陷部位提取結(jié)果

5 結(jié)束語

本文提出了一種用于零件表面缺陷部位定位與測(cè)量的方法。使用3D掃描儀得到零件表面點(diǎn)云，先通過預(yù)處理，與CATIA生成的完整零件點(diǎn)云進(jìn)行SAC-IA粗配準(zhǔn)和Huber損失函數(shù)優(yōu)化的ICP精配準(zhǔn)，然后進(jìn)行最鄰近查找，最后成功提取缺陷部位的精確點(diǎn)云。點(diǎn)云中包含缺陷部位的形狀、具體尺寸，能夠以此進(jìn)行在線的3D打印修復(fù)。實(shí)驗(yàn)表明：本文的檢測(cè)算法速度快、精度高。