童子良, 余學祥, 汪 濤, 王 虎, 蘇曉剛

(1.安徽理工大學 測繪學院,安徽 淮南 232001; 2.安徽理工大學 礦山采動災害空天地協(xié)同監(jiān)測與預警安徽普通高校重點實驗室,安徽 淮南 232001; 3.安徽理工大學 礦區(qū)環(huán)境與災害協(xié)同監(jiān)測煤炭行業(yè)工程研究中心,安徽 淮南 232001)

隨著智慧城市、虛擬現(xiàn)實(virtual reality,VR)技術、人工智能(artificial intelligence,AI)、逆向工程等概念的出現(xiàn),人們對空間三維信息的需求顯得日益迫切,傳統(tǒng)的獲取數(shù)據(jù)的方法已經(jīng)無法滿足實時、高精度、可視化的要求。三維激光掃描技術又被稱為三維實景復制技術,是測繪領域繼GPS空間定位技術之后又一次技術革新,它的出現(xiàn)為空間三維信息的獲取提供了全新的技術支持。與常規(guī)的觀測方法相比,三維激光掃描技術具有以下優(yōu)勢:① 可以實現(xiàn)復雜三維數(shù)據(jù)的采集與處理,實時性強;② 可以快速組建目標的三維模型等多種數(shù)據(jù);③ 無需與被測物體接觸,可在惡劣環(huán)境下工作;④ 可以獲取大量數(shù)據(jù)信息且精度較高[1]。因此,三維激光掃描技術已經(jīng)在工程變形監(jiān)測[2]、古建筑保護[3]、智慧城市建設[4]等領域獲得廣泛應用。

因為三維場景中有大量的平面特征,如墻面、道路、橋梁等等都是平面組成的,將掃描得到的點云數(shù)據(jù)擬合成平面,經(jīng)過點云配準、去噪、簡化之后可實現(xiàn)高精度地重建該物體的三維模型,所以點云的平面擬合是多數(shù)平面重建工程的重要組合部分。傳統(tǒng)的擬合方法有經(jīng)典的最小二乘法和特征值法。最小二乘法只考慮到觀測向量中的誤差,忽略了系數(shù)矩陣中的誤差;而特征值法無法消除噪聲對擬合過程的干擾,導致擬合平面的精度差,不具有魯棒性。文獻[5]提出了整體最小二乘的穩(wěn)健點云擬合方法,雖然可以同時顧及觀測向量和系數(shù)矩陣的誤差,但是由于點云數(shù)據(jù)是不等精度的,并且存在大量的異常點,采用整體最小二乘法估計的參數(shù)也不是最優(yōu)解;文獻[6]提出了穩(wěn)健加權總體最小二乘的點云平面擬合法,通過加權的方法先剔除異常點,再通過穩(wěn)健總體最小二乘法求出參數(shù)解。雖然以上方法可以估計出平面的參數(shù)解,但是其算法編譯復雜,解算效率低。近年來,一些研究者將主成分分析 (principal component analysis, PCA)方法應用到點云平面擬合和曲面擬合中,文獻[7-8]研究表明擬合效果好,且算法編譯簡單。此外,隨機采樣一致(random sample consensus, RANSAC)算法可以用少量樣本數(shù)據(jù)來評估模型參數(shù),對錯誤率超過50%的數(shù)據(jù)仍然能得到理想的處理結果,是最有效的魯棒估計方法[9-11]。本文提出采用RANSAC算法與PCA方法相結合,在擬合平面之前先用RANSAC算法對獲取的點云數(shù)據(jù)進行預處理,然后通過PCA方法對預處理的點云數(shù)據(jù)進行重采樣,從而達到消除粗差并增強魯棒性的效果。

1 RANSAC算法

RANSAC是采用迭代的方式從一組包含離群的被觀測數(shù)據(jù)中估算出數(shù)學模型參數(shù)的方法。該算法最早由Fischler和Bolles于1981年提出[12]。RANSAC算法中采集的樣本為2類點云數(shù)據(jù):① 內(nèi)點數(shù)據(jù)(inliers),即符合算法模型的數(shù)據(jù);② 異常數(shù)據(jù),記為外點(outliers),為不符合模型所需的點。

RANSAC算法的基本原理[13]如下所述:

(1) 假設所有的點云數(shù)據(jù)集為Ω,從該點云數(shù)據(jù)中選取一個符合模型條件要求的最小樣本數(shù)為n(n≥3,因為至少3個點才能構成一個平面)的最小樣本子集,利用最小二乘法計算出最小樣本子集所對應的模型參數(shù),并作為初始化模型參數(shù),接著計算Ω所有點集數(shù)據(jù)與初始參數(shù)之間的差值,根據(jù)已經(jīng)設定好的閾值T,當該差值小于設定的閾值,則將該樣本點歸為模型內(nèi)點,反之即視為不符合模型所需點將其剔除,記錄此時的內(nèi)點數(shù)。

(2) 重復以上過程,獲得更多的內(nèi)點與最佳模型參數(shù)。根據(jù)實際情況設定迭代次數(shù)K,當?shù)Y束時,計算出的最佳模型參數(shù)即為最佳估計參數(shù)。

(3) RANSAC算法要求在一定的置信概率下,其中n與至少獲得1個良性取樣子集的概率Q(Q>ε)滿足的關系式為:

Q=1-[1-(1-ε)n]K

(1)

其中,ε為樣本污染率。K的計算公式為:

(2)

一般情況下,ε、Q均會根據(jù)實際情況給定。

2 PCA方法

PCA作為一種多元統(tǒng)計技術,是由美國統(tǒng)計學家Pearson首次提出的。它從多指標分析出發(fā),運用統(tǒng)計分析原理提取少數(shù)幾個彼此不相關的綜合指標并且能夠保證原指標所提供的大量信息,通過原來變量的少數(shù)幾個線性組合來解釋原來變量絕大多數(shù)的信息。具體來說,它利用降維思想構造了原始數(shù)據(jù)的一個正交變換,新空間的基底去除了原始空間基底下數(shù)據(jù)的相關性,通過線性變化,將原始數(shù)據(jù)集變化為一組各維度線性無關的表示,這些新的基底即為主成分。目前,PCA方法廣泛地應用在圖像識別[14]、圖像處理[15]、特征信息提取[16]等領域。

本文PCA方法推導如下所述。

假設采樣點擬合的切平面目標函數(shù)為:

ax+by+cz=d

(3)

其中，n=(a,b,c)為單位法向量,且滿足:

a2+b2+c2=1, ‖n‖=1

(4)

按照上述準則,各個采樣點Pi到切平面的距離要盡可能小,目標函數(shù)為:

(5)

其中,xi、yi、zi為采樣點坐標;N為采樣點個數(shù)。該問題可轉(zhuǎn)化為求解極值問題,即

λ|a2+b2+c2-1|)

(6)

其中,λ為待定參數(shù)。由(6)式將目標函數(shù)改寫為:

(7)

(8)

對J求偏導得到:

(9)

將(9)式整理后可得:

(10)

對(10)式進行特征分解,將矩陣變換為對稱正定矩陣,故有3個實特征值λ0、λ1、λ2,由此求得它們所對應的特征向量為v0、v1、v2。

上述即為PCA標準推導過程,在優(yōu)化目標函數(shù)時,為所需的采樣點擬合一個平面,最大可能地確保它們都在平面上,即點到平面的距離最小,因此取最小特征值λ0的特征向量v0作為采樣點的法向量。根據(jù)求得的法向量得到所求的最佳擬合平面。

3 本文點云擬合新方法

本文為了提高擬合精度,首先在選取預處理點時設置閾值,以保證擬合平面點的“純度”,然后采用點到平面距離的標準差為擬合平面的精度,標準差越小,擬合點之間的分散程度越小,擬合效果越好;反之則越差。采用RANSAC算法在迭代過程消除點云數(shù)據(jù)中的粗差值,使得點云中的異常值大幅減少,然后利用PCA方法的降維思想對預處理的采樣點數(shù)據(jù)進行重采樣,達到二次擬合的效果。將RANSAC算法和PCA方法結合起來的點云擬合步驟如下:

(1) 根據(jù)給定的ε、Q及n,通過(2)式計算最小迭代次數(shù)K。

(2) 從所有的采樣點數(shù)中隨機選取m個點(m>3),利用(3)式和步驟(1)中給定的參數(shù)計算平面模型參數(shù)的初始值。

(3) 根據(jù)平面參數(shù)的初始值,計算允許值α,若該允許值在給定的閾值內(nèi),即為內(nèi)點,反之為噪聲點。

(4) 重復步驟(2) 、步驟(3),迭代K次,將內(nèi)點歸為一類,并統(tǒng)計其數(shù)量,數(shù)量最多的點集作為預處理點集。

(5) 利用PCA方法根據(jù)RANSAC算法得到的預處理點集,計算平面的法向量n(a,b,c)。

(6) 計算所有采樣點到目標平面的距離,計算公式為:

(11)

(7) 求出di的標準差δ,計算公式為:

(12)

(13)

(8) 當di>2δ時,認為該點為噪聲點,否則保留該點。

(9) 重復步驟(1) ～步驟(4),進行K次迭代,直到所有點都在閾值范圍之內(nèi)。

(10) 得到最佳擬合平面方程。

4 實驗分析

4.1 仿真實驗

采用仿真實驗數(shù)據(jù)對本文方法進行驗證,并與最小二乘法、特征值法擬合結果作對比。設需要擬合的空間平面方程為:

(15)

通過Matlab從已知平面隨機選取1 000個點,該數(shù)據(jù)為未受到異常值污染情況下的仿真數(shù)據(jù),如圖1所示;在上述仿真數(shù)據(jù)中加入50個噪聲點,如圖2所示。使用最小二乘法、特征值法及本文方法對圖1、圖2的仿真數(shù)據(jù)進行參數(shù)估計,并計算δ,結果見表1、表2所列。表1、表2中,括號里的數(shù)據(jù)為設定的參數(shù)值。

圖1 不包含噪聲點的仿真數(shù)據(jù) 圖2 包含50個噪聲點的仿真數(shù)據(jù)

表1 未加入異常值的仿真實驗參數(shù)估計結果

表2 加入50個異常點的仿真實驗參數(shù)估計結果

雖然數(shù)據(jù)是隨機抽取的,每次運行的結果可能會有一些小的差別,但由表1、表2可知,當數(shù)據(jù)不存在異常值時,雖然3種方法的平面擬合參數(shù)值都十分近似,都接近預設的參數(shù)值,而本文方法進行預處理后的點云數(shù)據(jù)離散程度更低,故擬合后的效果更好,其δ為2.457 2×10-17,低于最小二乘法和特征值法的δ值,說明在進行擬合之前,對采樣的點云數(shù)據(jù)進行預處理是十分必要的;加入一些異常點后,最小二乘法和特征值法得出的估計參數(shù)值與設定的參數(shù)值之間存在很大的誤差,而本文方法的δ為0.002 0,遠遠低于前2種方法的δ值0.475 0、0.217 8,故擬合效果極佳。本文方法所得估計值與預設值十分接近,由此可知本文方法穩(wěn)定且精度高,具有魯棒性。

4.2 實例分析

為了進一步檢驗本文方法在點云平面擬合中的可行性,通過采樣實體建筑物的實測點云數(shù)據(jù)進行驗證。實驗采用中海達HS650高精度三維激光掃描儀,該儀器采用脈沖式的測距方式,可多次回波輸入,其激光發(fā)散性小于等于0.35 mrad,最大測程達到650 m,測距精度為5 mm至100 m,全景分辨率大于7 000×104像素。

以安徽理工大學測繪學院大門的點云數(shù)據(jù)為測試實例,如圖3所示。由圖3可知,該模型具有大量的平面特征信息,故手動提取位于同一平面特征上的點云數(shù)據(jù),如圖4所示。

分別采用最小二乘法、特征值法和本文方法對實測的點云數(shù)據(jù)進行平面擬合,結果見表3所列。從表3可以看出,當實測的點云數(shù)據(jù)中不含有噪聲點時,3種方法擬合的參數(shù)值大致相同,最小二乘法的精度稍低,特征值法與本文方法的δ較小,擬合精確性較高。當點云數(shù)據(jù)混入噪聲點時,以表3中3種方法得到模型參數(shù)的平均值作為此平面的參考值,此時的噪聲點與原始點云數(shù)據(jù)不處于同一平面上,即這些噪聲點到原始點云擬合平面的距離超過了所設定的閾值范圍,此時通過3種方法對點云數(shù)據(jù)進行處理,所得結果見表4所列。由表4可知,當點云數(shù)據(jù)中混入噪聲點時,傳統(tǒng)方法擬合精度較低,通過前2種方法計算的δ分別為0.029 0、0.013 8,遠大于本文方法得出的9.375 6×10-4。由此可知,本文方法不僅可以去除噪聲點,而且擬合平面參數(shù)值與參考值之間的偏差極小,精度和穩(wěn)定性都高。

圖3 測繪學院大門點云數(shù)據(jù)

圖4 處于同一平面上的點云數(shù)據(jù)

表3 不含噪聲點時3種方法估計的實測數(shù)據(jù)參數(shù)

表4 包含噪聲點時3種方法估計的實測數(shù)據(jù)參數(shù)

5 結 論

通過三維激光掃描儀采樣時,由于環(huán)境、儀器等各種因素的限制,獲取的點云數(shù)據(jù)不可避免地存在異常值或者粗差,而傳統(tǒng)的最小二乘法和特征值法沒有考慮這些噪聲點對于擬合精度和準確度的影響。針對以上問題,本文提出采用RANSAC算法與PCA方法相結合的方法,在擬合平面之前先用RANSAC算法對獲取的點云數(shù)據(jù)進行預處理,先去除一部分的粗差點,然后通過PCA方法對預處理的點云數(shù)據(jù)進行擬合,從而達到消除粗差并增強魯棒性的效果。通過仿真實驗和工程實例進行了驗證,實測實驗結果表明采用本文方法得出的標準差為9.375 6×10-4,遠遠低于傳統(tǒng)的方法得出的標準差,由此可知通過本文方法可以得到更加精確的擬合平面, 本文方法具有更好的適用性。