基于三維點云螺紋中軸線的擬合方法研究

謝張寧， 姜志華， 朱惠臣， 李智瑋， 雷李華,，孔 明， 王道檔， 張 波

(1. 中國計量大學(xué)計量測試工程學(xué)院，浙江杭州310018；2. 上海計量測試技術(shù)研究院，上海201203)

1 引 言

高精度螺紋規(guī)是精密機械、航天航空、石油化工、核電等領(lǐng)域中極為重要的基礎(chǔ)零部件，其幾何精度對裝配質(zhì)量具有決定性影響[1]?，F(xiàn)有的高精度測量方法眾多，但均是基于一維、二維的螺紋測量技術(shù)，無法獲取三維綜合參數(shù)，檢測結(jié)果難以反映實際精度，導(dǎo)致裝配質(zhì)量不能準(zhǔn)確控制。要改變“粗糙螺紋”生產(chǎn)的現(xiàn)狀，需要對三維螺紋測量儀器設(shè)備進(jìn)行研究開發(fā)。

當(dāng)前的三維測量方法在實際測量前，首先要實現(xiàn)機器坐標(biāo)系MCS與工件坐標(biāo)系WCS轉(zhuǎn)換[2]，在螺紋面的測量中，螺紋中軸線的準(zhǔn)確測定對螺紋的后續(xù)測量至關(guān)重要。然而，在對螺紋規(guī)測量時，仍舊是依靠測量人員的經(jīng)驗,手動調(diào)整螺紋規(guī)位置，然后進(jìn)行后續(xù)測量。這種測量方法不僅增加了測量難度，也增大了測量結(jié)果的不確定度。

本文研究了一種能夠快速精準(zhǔn)擬合螺紋中軸線的方法。在測量時，螺紋規(guī)放置在工作臺上進(jìn)行調(diào)整，但是工件坐標(biāo)系不可能與機器坐標(biāo)系完全重合，所以要使用螺紋規(guī)的基準(zhǔn)來建立工件坐標(biāo)系。將測量儀器對準(zhǔn)螺紋規(guī)上下底面的頂針孔實現(xiàn)定位，在機器坐標(biāo)系MCS下，非接觸掃描得到的螺紋規(guī)外表面點云;經(jīng)過隨機增量算法計算點云凸包，制成螺紋凸包點集;基于最小二乘的概念建立螺紋中軸線的擬合數(shù)學(xué)模型，得出中軸線的單位向量v(a,b,c);再基于單位向量得到螺紋底面中心A，坐標(biāo)旋轉(zhuǎn)矩陣R，坐標(biāo)偏移矩陣S;螺紋測量儀器根據(jù)上述參數(shù)得出工件坐標(biāo)系。

由于螺紋復(fù)雜的三維形體結(jié)構(gòu)，在建立數(shù)學(xué)模型對螺紋面進(jìn)行直接計算時，牙槽表面、螺紋表面缺陷點云數(shù)據(jù)均會影響螺紋中軸線的擬合，導(dǎo)致結(jié)果不準(zhǔn)確。為防止干擾，本文用隨機增量算法計算螺紋表面三維點云凸包(convex hull)[3]；將螺紋牙槽表面、螺紋表面缺陷點云數(shù)據(jù)篩去，對于空間離散點集，求解最小的凸多面體，其中把空間中點集所有的點包含在里面的凸多面體就是凸包[4]?；谌S點云凸包建立中軸線擬合數(shù)學(xué)模型，能有效濾除掉螺紋牙槽、硬件缺陷、掃描誤差等形成的干擾點云，減少了后續(xù)算法的運算量，加快了運行時間，提升了運算結(jié)果。

經(jīng)過計算凸包的螺紋表面點云實質(zhì)上是以公稱直徑為上下底直徑的圓柱面。要擬合該圓柱面的軸線就要以面上的點與軸線的位置關(guān)系建立數(shù)學(xué)模型，并采用優(yōu)化算法求解。常用的優(yōu)化算法有最小二乘法、爬山算法、模擬退火算法、遺傳算法等[5～8]。用爬山算法、模擬退火算法、遺傳算法等隨機優(yōu)化算法確定唯一解的方法雖然可行，但隨機優(yōu)化算法大多處理高階次線性問題；目前通用、認(rèn)可度高的算法仍為最小二乘法，在工程測量上，用最小二乘法解決超定方程[9]、線性擬合[10～12]等數(shù)據(jù)量大的低階次線性問題速度快、結(jié)果準(zhǔn)確；且隨機算法優(yōu)化過程具有隨機性，容易出現(xiàn)早熟或過擬合等問題。因此，擬合螺紋中軸線采用最小二次算法具有較強的實用性。

2 螺紋中軸線數(shù)學(xué)模型

在空間機器坐標(biāo)系中，螺紋中軸線l可以表示為經(jīng)過點A(x0,y0,z0)，單位向量為v(a,b,c)的直線l，其點斜式為[13]：

根據(jù)最小二乘原理，使得經(jīng)過三維凸包計算的螺紋表面點云P(xi,yi,zi)，i=1,…,N到擬合直線距離d的平方和最?。?/p>

(1)

(2)

假設(shè)：

(3)

將式(2)、式(3)代入式(1)，可得最小二乘意義上的螺紋中軸線數(shù)學(xué)模型：

(bδz-cδy)2]+f(a,b,c)

(4)

由式(4)可以看出，對于任意的a,b,c均有：

(5)

當(dāng)δx=δy=δz=0時，有：

min{f(a,b,c)}
a2+b2+c2=1

其中：

f(a,b,c)=(1-a2)B11+(1-b2)B22+

(1-c2)B33-2abB12-2acB13-2bcB23

(6)

3 基于三維點云的最小二乘擬合

3.1 制成三維凸包點集

由于螺紋表面三維結(jié)構(gòu)復(fù)雜以及可能存在缺陷的原因，要對點云數(shù)據(jù)計算凸包。計算凸包是根據(jù)算法在點云中篩選出能構(gòu)造包含所有點的最小多面體的點集。

螺紋表面點云是三維數(shù)據(jù)點，且數(shù)據(jù)量龐大，用隨機增量算法計算凸包效果好、速度最快[15]。計算流程為：

(1)首先用篩選函數(shù)找到x,y,z軸上的邊界點，如x,y,z值最小的點，設(shè)其為P1。

(2)隨機找到與P1不共線的2個點P2,P3，再找到不共面的第4個點P4，構(gòu)成四面體。

(3)若新點Pi在凸包內(nèi)部或邊界，則忽略。若新點Pi在凸包外部，見圖1(a)，就要計算原凸包相對于Pi的面，見圖1(b)所示;將Pi點能看到的所有面刪除，同時將邊上的點與Pi相連組成新的凸包，見圖1(c)。

(4)依次遍歷余下的點，制成三維凸包點集。

圖1 隨機增量算法計算凸包Fig.1 Random incremental algorithm to calculate convex hull

3.2 單位向量的求解

將式(6)引入矩陣B=[Bij]，整理成矩陣式，則可得：

由于B是對稱方陣，不難看出，當(dāng)[a,b,c]T是矩陣B的最大特征值λmax(B)對應(yīng)的單位特征向量時，上式可以取得最小[16]。此時，[a,b,c]T為螺紋中軸線單位方向向量：

3.3 底面中心點的求解

在轉(zhuǎn)換后的坐標(biāo)系o-x′y′z′中，取圓環(huán)最外層，設(shè)其擬合曲線為：

式中M′為螺紋外徑擬合結(jié)果。將上式展開，令：

整理得：

最小二乘法求解上述超定方程解得[18]：

其中：

3.4 測量坐標(biāo)系到工件坐標(biāo)系的轉(zhuǎn)換

測量坐標(biāo)系到工件坐標(biāo)系的轉(zhuǎn)換分兩步：

(1)先將測量坐標(biāo)系進(jìn)行原點偏移。即將A(x0,y0,z0)作為工件坐標(biāo)系的新原點o′，只需在原坐標(biāo)系各軸減去相應(yīng)的偏移量，可得偏移矩陣S：

(2)將坐標(biāo)軸旋轉(zhuǎn)，使z軸與對中軸線單位向量重合。坐標(biāo)旋轉(zhuǎn)矩陣R在求解中軸線底面中心時已經(jīng)求出。

4 仿真實驗

4.1 有效性驗證

用MATLAB軟件編寫確定螺紋中軸線的程序。用MATLAB自帶的[v,λ]=eig(X)函數(shù)求解特征向量λ對應(yīng)矩陣的特征值，v中3個列向量為λ對應(yīng)的特征向量。求解單位向量程序框圖見圖2。

圖2 求解單位向量程序框圖Fig.2 Block diagram for solving unit vector

用MATLAB編程求解底面中心，將螺紋點云對應(yīng)的三個軸數(shù)據(jù)向量經(jīng)過矩陣R后得到處于o-x′y′z′坐標(biāo)系的點云數(shù)據(jù)，再由MATLAB自帶的X-1=inv(X)函數(shù)進(jìn)行坐標(biāo)反變換，得到原坐標(biāo)系的底面中心數(shù)據(jù)。

圖3 求解底面中心坐標(biāo)程序框圖Fig.3 Block diagram of the bottom center coordinates

用圓柱點云圖檢驗本文中軸線擬合算法的有效性。已知該圓柱點云中軸線參數(shù)為：

點云數(shù)據(jù)由scatter3函數(shù)繪制，在計算完成的中軸線上取若干點用MATLAB自帶函數(shù)plot3繪制，如圖4所示。

圖4 基于三維凸包的最小二乘擬合結(jié)果Fig.4 Results of least square fitting based on 3D convex hull

中軸線擬合算法解得單位向量vt：

圓柱底面中心點At：

坐標(biāo)旋轉(zhuǎn)矩陣R：

坐標(biāo)偏移矩陣S：

在有限值域區(qū)間內(nèi)，以一條直線上的點到另一條直線的距離的最大值作為指針，定量評定不同直線與已知直線之間的符合程度。本文以兩直線在螺紋旋合長度內(nèi)，某一直線兩端點到另一條直線的距離的最大者為評定參數(shù)，評價兩條直線的符合程度。其算法過程為：

(1)首先以任意一條直線l的單位向量為法向量，取點云邊界點Pα生成上截面α；下截面法向量相同。

(2)將上截面方程轉(zhuǎn)換為參數(shù)方程，將參數(shù)方程代入兩直線lt、lz的點斜式，可解到擬合直線lt在上截面α的端點坐標(biāo)(xlt，α,ylt，α,zlt，α)，已知直線lz在上截面α的端點坐標(biāo)(xlz，α,ylz，α,zlz，α)。

(3)同理按步驟(1)生成下截面β，擬合直線lt在下界面β的端點坐標(biāo)(xlt，β,ylt，β,zlt，β)，已知直線lz在下截面β的端點坐標(biāo)(xlz，β,ylz，β,zlz，β)。

(4)分別計算截面α,β上，兩直線端點距離Eα、Eβ，用MATLAB自帶的Max函數(shù)取兩者的較大值即為評定參數(shù)E。

運算所得E=0.21 μm，所得中軸線參數(shù)與lz直線基本重合，算法可行。

4.2 中軸線擬合仿真實驗

為了比較研究中軸線擬合算法的可行性，本文用MATLAB軟件編譯基于三維凸包的最小二乘算法擬合螺紋中軸線，并與投影法比較驗證其可行性。根據(jù)國家標(biāo)準(zhǔn)《普通螺紋基本尺寸》表1[19]里的螺紋小徑、中徑數(shù)據(jù)建立了公稱直徑為80 mm,螺距為4 mm的普通外螺紋面的三維模型，并在模型表面隨機成塊刪除點云形成螺紋自身表面缺陷、掃描失真等造成的干擾因素，如圖5所示。

圖5 帶有干擾因素的外螺紋三維模型面Fig.5 Three-dimensional model surface of external thread with interference factors

將三維模型面保存為obj格式，在Visual Studio 2017環(huán)境下用PCL18.1編程庫將其生成外螺紋面三維點云。外螺紋面底面圓心坐標(biāo)為(-40,150,-40)，頂端圓心坐標(biāo)為(-40,0,-40)，如圖6所示。

圖6 帶有干擾的外螺紋三維點云Fig.6 3D point cloud image of external thread with interference

4.2.1 投影法擬合中軸線

投影法是軸體測量中常見的方法,在測量時，其將螺紋器具投影至xOy面、yOz面和zOx面，基于平面測量[20]。

本文將三維點云投影至xOy面、yOz面經(jīng)過計算邊界坐標(biāo)，將xOy面、yOz面的邊界點y值取中值，分別擬合直線。求得投影面上的中線，過所得中線作面，根據(jù)2個不平行的平面相交確定1條直線的原則[13]確定空間螺紋中軸線，用MATLAB編程得中軸線一般方程：

4.2.2 最小二乘擬合螺紋中軸線

用隨機增量算法對帶有干擾因素的外螺紋三維點云作convex hull計算，計算結(jié)果如圖7所示。

圖7 經(jīng)過convex hull處理后外螺紋三維點云Fig.7 Thread 3D point cloud after convex hull processed

經(jīng)過convex hull處理后的外螺紋三維點云導(dǎo)入MATLAB，將參數(shù)區(qū)間限制為點云數(shù)值區(qū)域內(nèi)，用投影法和基于三維凸包的最小二乘算法擬合中軸線。圖8、圖9分別為螺紋中軸線三維擬合結(jié)果以及其在xOy面、yOz面的投影圖。

圖8 基于三維凸包的最小二乘算法和投影法擬合結(jié)果Fig.8 Fitting results of least squares algorithm and projection method based on three-dimensional convex hull

圖9 擬合結(jié)果的投影Fig.9 Projection of the fitting result

用基于三維凸包的最小二乘算法擬合的螺紋中軸線所確定的中軸線lls參數(shù)為：

單位向量vls

螺紋底面中心Als

坐標(biāo)旋轉(zhuǎn)矩陣R

坐標(biāo)偏移矩陣S

4.3 結(jié)果驗證

由圖8、圖9可知，投影法擬合的直線與基于三維凸包最小二乘法擬合的直線基本重合。為了驗證基于三維凸包的最小二乘算法擬合中軸線的精度和符合程度，本文采用計算lty,lls與外螺紋點云距離的方差[21]來表示算法擬合精度：

符合程度評價標(biāo)準(zhǔn)與上節(jié)相同，在有限值域區(qū)間內(nèi)，將1條線段的2個端點到另1條直線的距離的最大者作為評定標(biāo)準(zhǔn)E。以max函數(shù)、min函數(shù)對點云的x，y，z值取最大或最小，找到上下邊界點Pb1，Pb2，再按照第4.1節(jié)E算法過程求解。Xdata，Ydata，Zdata分別為點云的x，y，z數(shù)據(jù)集，上、下邊界點的生成代碼如下：

%生成上邊界點

max_X=max(Xdata); max_Y=max(Xdata);

max_Z=max(Xdata);

P_a=[max_X,max_Y,max_Z];

%生成下邊界點

min_X=min(Xdata); min_Y=min(Xdata);

min_Z=min(Xdata);

P_b=[min_X,min_Y,max_Z];

為比較投影法和最小二乘算法的計算速度，在處理器為Intel(R)Core(TM)i5-3230M的計算機上運行，結(jié)果見表1。

表1 擬合確定的中軸線比較表Tab.1 Table of fitted central axis variances

從擬合結(jié)果來看，用二維投影法和最小二乘法擬合的中軸線與外螺紋點云距離的方差結(jié)果都為0.34,且兩直線在旋合長度范圍內(nèi)端面點距離的最大者，值為0.15 μm，達(dá)到了三維高精度測量標(biāo)準(zhǔn)，說明最小二乘法能夠獲得投影法一致的效果;在計算機執(zhí)行速度上來看，最小二乘法擬合中軸線程序在MATLAB上運行時間只有15.3 ms，在同樣的硬件條件下,逐次求解的投影法計算速度為1 196.7 ms。由此可以得出基于三維點云的最小二乘擬合算法是一種快速、準(zhǔn)確的螺紋中軸線的擬合方法。

綜上，最小二乘擬合中軸線的方法是從三維的角度上擬合中軸線，其囊括的維度更高更具有說服性，且計算時間短，在高精度螺紋規(guī)的測量上更具有實用性和推廣價值。

5 結(jié) 論

本文從螺紋中軸線特征入手，建立了螺紋中軸線數(shù)學(xué)模型，提出了基于三維凸包的最小二乘算法擬合螺紋中軸線的方法，通過MATLAB軟件編程驗證算法的可行性，最后結(jié)合三維模擬仿真實驗，并經(jīng)二維投影法驗證，得出以下結(jié)論：

(1)經(jīng)過三維凸包計算的最小二乘算法在擬合具有三維復(fù)雜結(jié)構(gòu)的螺紋中軸線具有與投影法一致的結(jié)果。

(2)經(jīng)過三維凸包處理的最小二乘法從三維的角度解決問題，計算速度快，在高精度螺紋測量上，具有實用性和推廣價值。