袁 輝

（中鐵第四勘察設(shè)計(jì)院集團(tuán)有限公司 武漢 430063）

1 引言

三維激光雷達(dá)（Light Detection and Ranging，Li-DAR）測(cè)量技術(shù)是一種非接觸、快速度、高精度和高密度測(cè)量數(shù)據(jù)獲取技術(shù)［1～4］，已廣泛地應(yīng)用于智慧城市、工程測(cè)量等領(lǐng)域。從點(diǎn)云數(shù)據(jù)中提取各種所需的空間及屬性信息是點(diǎn)云數(shù)據(jù)后處理的一項(xiàng)重要內(nèi)容［5～7］，其中從點(diǎn)云數(shù)據(jù)中擬合圓柱體是工程測(cè)量中一個(gè)常見的問題［8～10］，例如管徑的尺寸檢測(cè)、立柱的直徑測(cè)量和垂直度檢測(cè)、隧道環(huán)片幾何尺寸測(cè)量和形變檢測(cè)等。

當(dāng)前圓柱體的擬合方法主要有基于幾何分析的擬合算法［10～12］、基于坐標(biāo)轉(zhuǎn)換的擬合算法［13～15］以及基于智能算法的擬合算法［16～17］等?；趲缀畏治龅臄M合算法直接從圓柱體方程的幾何定義出發(fā)，建立誤差方程進(jìn)行最小二乘平差，但是由于誤差方程的非線性的特點(diǎn)，擬合結(jié)果受參數(shù)初值影響較大，容易陷入局部收斂［18～19］?；谧鴺?biāo)變換的擬合方法是將圓柱體的軸線旋轉(zhuǎn)到與豎直方向平行，然后投影到水平面上進(jìn)行圓形擬合，文獻(xiàn)［13～15］建立了基于坐標(biāo)轉(zhuǎn)換的非線性最小二乘平差模型，該方法也同樣依賴于初值的選擇。文獻(xiàn)［16～17］利用遺傳算法擬合圓柱體模型參數(shù)，但該算法參數(shù)多，容易早熟收斂，同時(shí)計(jì)算效率低。

本文針對(duì)點(diǎn)云數(shù)據(jù)的圓柱體的擬合問題，根據(jù)基于坐標(biāo)轉(zhuǎn)換的擬合算法思想，構(gòu)造了旋轉(zhuǎn)矩陣將三維的非線性的圓柱體擬合問題轉(zhuǎn)換為二維的線性的圓形擬合，降低了處理的復(fù)雜度，同時(shí)利用迭代逼近法求解中軸線的方向向量，另外再根據(jù)擬合中誤差進(jìn)行點(diǎn)云噪聲點(diǎn)剔除處理進(jìn)行迭代擬合，提高了圓柱體模型參數(shù)擬合的準(zhǔn)確度。

2 圓柱體擬合問題

在三維空間中，圓柱體是到其中心軸線距離為常數(shù)r的點(diǎn)集合，因此，圓柱體由以下7個(gè)參數(shù)唯一確定：圓柱體中軸線上的點(diǎn)po(xo,yo,zo)、中軸線方向向量s(a,b,c)、圓柱體半徑r。約束方向向量s為單位向量，即a2+b2+c2=1，因此圓柱面的方程可以表示為

定義數(shù)據(jù)點(diǎn)集合P={p1,p2,…,pN}為三維激光掃描儀采集得到的圓柱體表面點(diǎn)云，其中N為點(diǎn)數(shù)，pj=(xj,yj,zj)，pj∈P，j∈[1,N]。在最小二乘算法擬合過程中，定義點(diǎn)的擬合誤差為點(diǎn)到圓柱體中軸線的距離與圓柱體半徑之差，表示為

因此，采用最小二乘擬合使誤差函數(shù)F值最小，即：

基于幾何分析的擬合算法［11～12］對(duì)式的誤差方程進(jìn)行線性化處理，舍去高次項(xiàng)，然后按照高斯牛頓法迭代計(jì)算。由于式的誤差方程包含根式，同時(shí)參數(shù)xo,yo,zo,a,b,c之間的相關(guān)性大，在參數(shù)初值不準(zhǔn)確的情況下容易造成局部收斂，導(dǎo)致錯(cuò)誤的擬合結(jié)果，因此非線性最小二乘擬合圓柱體過程中對(duì)參數(shù)初值的選擇要求很高。

3 基于坐標(biāo)變換的圓柱體擬合

根據(jù)基于坐標(biāo)變換擬合圓柱體的思想［13］，將圓柱體的中軸線旋轉(zhuǎn)到與z軸方向(0,0,1)平行，此時(shí)圓柱體在xoy平面上的投影為圓形，圓形半徑與圓柱體半徑相等。因此，一旦確定了中軸線方向向量s，便可將三維的圓柱擬合問題轉(zhuǎn)換為二維的圓形擬合問題，而圓形擬合可以用線性方程直接平差處理完成，降低了計(jì)算的復(fù)雜度。

3.1 基于坐標(biāo)變換的圓柱體擬合

定義圓柱體的中軸線方向向量s(a,b,c)繞x軸旋轉(zhuǎn)角度θ∈[0,2π)，再繞y軸旋轉(zhuǎn)角度φ∈[0,2π)后變換為z軸方向向量(0,0,1)，因此旋轉(zhuǎn)變換矩陣R為

變換矩陣滿足R-1=RT，另外方向向量s(a,b,c)需滿足

求解得到中軸線方向向量s(a,b,c)為a=sinφ，b=-sinθcosφ，c=cosθcosφ。根據(jù)旋轉(zhuǎn)矩陣R將點(diǎn)的坐標(biāo)(x,y,z)變換為(x′,y′,z′)。

對(duì)于變換后的點(diǎn)坐標(biāo)(x′,y′,z′)取x′和y′擬合二維平面圓形，圓形方程定義為

令A(yù)=-2x′o，B=-2y′o，C=(x′o)2+(y′o)2-r2，因此，式（8）可以轉(zhuǎn)換為線性方程

對(duì)線性方程式（9）根據(jù)間接法最小二乘平差直接求解得到參數(shù)A、B、C的值，然后計(jì)算得到。其中擬合的圓形的半徑r也是圓柱體的半徑，然后將擬合的圓心坐標(biāo)(x′o,y′o,0)進(jìn)行反向坐標(biāo)變換轉(zhuǎn)換為圓柱體軸線上的點(diǎn)po(xo,yo,zo)。在擬合圓形過程中，定義點(diǎn)是pj=(xj,yj,zj)根據(jù)式（6）轉(zhuǎn)換后的點(diǎn)，其擬合誤差為

相應(yīng)地，點(diǎn)云集合的誤差函數(shù)F′為

3.2 迭代逼近法求解方向向量

根據(jù)上述基于坐標(biāo)轉(zhuǎn)換的圓柱體擬合算法，在已知中軸線方向向量s的情況下便能很方便地?cái)M合得到圓柱體的7個(gè)參數(shù)，因此其中的關(guān)鍵問題是如何準(zhǔn)確得到角度θ和φ，文獻(xiàn)［13］通過建立整體的坐標(biāo)轉(zhuǎn)換誤差方程進(jìn)行非線性最小二乘進(jìn)行迭代計(jì)算，計(jì)算過程復(fù)雜也需給定合適的參數(shù)初值，本文通過迭代逼近法求解θ和φ。

迭代逼近法解角度θ和φ的思路如下：初始給定角度值的值域范圍和角度步長(zhǎng)較大，先初步定位當(dāng)前擬合誤差最小的一組θ和φ角度值，然后在這組角度值鄰域范圍內(nèi)采用更小的步長(zhǎng)迭代查找新的最優(yōu)角度值，直到兩次擬合誤差值差別足夠小時(shí)則停止迭代。本文先根據(jù)角度θ和φ的不同取值，計(jì)算對(duì)應(yīng)的擬合誤差值F′，在F′最小時(shí)的那一組θ和φ值為最優(yōu)取值，利用迭代逼近法求解θ和φ的具體步驟如下：

1）定義角度θ和φ在第k=1次迭代時(shí)的初始值域范圍為2π，初始角度查找步長(zhǎng)Δk=π/10，初始最優(yōu)擬合誤差，初始最優(yōu)角度值為θbest=0和φbest=0，迭代精度ε=1×10-5。

2）根據(jù)當(dāng)前角度步長(zhǎng)Δk，在當(dāng)前角度范圍和內(nèi)按步長(zhǎng)遍歷取值，即，，根據(jù)當(dāng)前角度和進(jìn)行坐標(biāo)變換以及圓形擬合，計(jì)算當(dāng)前擬合誤差。

在上述迭代逼近法求解角度步驟中，在每一次的迭代過程中進(jìn)行的擬合計(jì)算次數(shù)較少，同時(shí)每經(jīng)過一次迭代，角度的值域范圍就顯著減小，同時(shí)當(dāng)前最優(yōu)的角度值的精度提高一個(gè)等級(jí)，可以很快地收斂到最優(yōu)角度值。盡管利用迭代逼近法增加了擬合處理的次數(shù)，但在最小二乘擬合過程中，是對(duì)二維的平面圓進(jìn)行線性方程的擬合處理，計(jì)算簡(jiǎn)單方便，也不會(huì)造成計(jì)算性能的顯著下降。

3.3 圓柱體噪聲點(diǎn)剔除

根據(jù)三維激光掃描測(cè)量得到的圓柱體點(diǎn)云數(shù)據(jù)，往往有較多的噪聲點(diǎn)，如果不去噪直接擬合圓柱體的模型參數(shù)，會(huì)導(dǎo)致較大的誤差。因此本文通過計(jì)算擬合結(jié)果的中誤差值δ，剔除擬合誤差超過3倍中誤差的噪聲點(diǎn)后再次進(jìn)行擬合，然后重復(fù)該過程，直到無噪聲點(diǎn)被剔除為止。在后續(xù)的迭代過程中，角度θ和φ的初始值域范圍可以按上一擬合過程中的結(jié)果角度取一定的鄰域范圍，可以加快迭代逼近過程的收斂速度。

4 實(shí)驗(yàn)分析

為了驗(yàn)證本文提出的圓柱體擬合方法的有效性，同常規(guī)圓柱體擬合的最小二乘算法［12］進(jìn)行對(duì)比實(shí)驗(yàn)，最小二乘算法按實(shí)際情況選取合適初值。實(shí)驗(yàn)數(shù)據(jù)為某一建筑圓形立柱的三維激光掃描點(diǎn)云數(shù)據(jù)，如圖1所示，該立柱點(diǎn)云包含大約5000個(gè)點(diǎn)，點(diǎn)間距2cm。

圖1 建筑立柱點(diǎn)云數(shù)據(jù)

在都不進(jìn)行噪聲點(diǎn)剔除情況下用不同算法進(jìn)行圓柱體擬合，擬合結(jié)果如表1所示。

從表1中可以看出，兩種算法中軸線上點(diǎn)(xo,yo,zo)差別較大是因?yàn)椴煌惴ㄟx取軸線點(diǎn)的方式不同，現(xiàn)將最小二乘算法的軸線點(diǎn)沿其軸線方向移動(dòng)使zo=2.20020m，移動(dòng)后最小二乘算法的xo=-2.2854m、yo=4.4846m，移動(dòng)后兩種算法軸線點(diǎn)距離差為0.1mm，差別很小。本文算法與最小二乘算法的擬合半徑相差0.1mm，中軸線角度相差0.08°，擬合誤差相差0.0002，二種算法結(jié)果很接近，是因?yàn)槎x的誤差目標(biāo)方程實(shí)質(zhì)是相同的，所以最終收斂擬合的結(jié)果基本一致。但是常規(guī)擬合算法需要在給定合適的初值情況才能正確收斂，而因此本文算法無需初值選擇就能夠準(zhǔn)確地?cái)M合圓柱體參數(shù)。

表1 立柱點(diǎn)云圓柱體擬合結(jié)果

圖2為一地鐵圓形盾構(gòu)隧道環(huán)片的點(diǎn)云數(shù)據(jù)，點(diǎn)數(shù)74892個(gè)，點(diǎn)間距3cm，隧道設(shè)計(jì)半徑2.75m，從圖中可以看出點(diǎn)云中包含大量的非隧道輪廓的噪聲點(diǎn)。

圖2 隧道環(huán)片點(diǎn)云數(shù)據(jù)

采用本文的圓柱體擬合算法對(duì)該環(huán)片點(diǎn)云數(shù)據(jù)進(jìn)行直接擬合和按本文方法進(jìn)行迭代去噪擬合的結(jié)果如表2所示。

從表2中可以看出，未去噪直接擬合圓柱體的半徑r=2.7005m，與設(shè)計(jì)半徑2.75m相差很大，屬于錯(cuò)誤的擬合結(jié)果，擬合誤差達(dá)到4285.4418。按照3倍的中誤差δ=0.2392m剔除噪聲點(diǎn)后，剩余70883個(gè)點(diǎn)，再次擬合得到圓柱體半徑為r=2.7408m，擬合誤差顯著減小到311.3382。在重復(fù)迭代去噪擬合6次后，無新的噪聲點(diǎn)被剔除，剩余67427個(gè)點(diǎn)，最終擬合半徑r=2.7528m，與設(shè)計(jì)半徑值接近，同時(shí)擬合誤差值F′=10.1426，遠(yuǎn)小于未去噪時(shí)的擬合誤差。因此通過計(jì)算擬合中誤差然后按照3倍中誤差限值進(jìn)行去噪擬合能有效提高圓柱體的擬合精度，從而滿足實(shí)際工程項(xiàng)目的需求。

表2 隧道環(huán)片點(diǎn)云圓柱體擬合結(jié)果

5 結(jié)語

本文針對(duì)圓柱體點(diǎn)云數(shù)據(jù)的擬合問題，根據(jù)坐標(biāo)轉(zhuǎn)換擬合的思想，利用迭代逼近法準(zhǔn)確求解圓柱體的中軸方向向量，然后根據(jù)中軸方向旋轉(zhuǎn)點(diǎn)云數(shù)據(jù)到與豎直方向平行，將三維的非線性圓柱體擬合問題轉(zhuǎn)化二維的線形圓形擬合問題，降低了擬合的復(fù)雜度。通過建筑立柱圓柱體點(diǎn)云數(shù)據(jù)的擬合實(shí)驗(yàn)，證明本文方法在無需給定初值情況下的擬合精度與正確擬合收斂的非線性最小二乘法的擬合精度相當(dāng)。根據(jù)地鐵隧道點(diǎn)云數(shù)據(jù)的實(shí)驗(yàn)，通過計(jì)算擬合中誤差然后按照3倍的限差剔除噪聲點(diǎn)后迭代處理能夠有效提高擬合的精度，證明了本文圓柱體點(diǎn)云數(shù)據(jù)能夠滿足實(shí)際工程擬合精度的要求。