童子良, 余學(xué)祥, 汪 濤, 王 虎, 蘇曉剛

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

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

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

1 RANSAC算法

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

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

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

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

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

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

(1)

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

(2)

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

2 PCA方法

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

本文PCA方法推導(dǎo)如下所述。

假設(shè)采樣點(diǎn)擬合的切平面目標(biāo)函數(shù)為:

ax+by+cz=d

(3)

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

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

(4)

按照上述準(zhǔn)則,各個(gè)采樣點(diǎn)Pi到切平面的距離要盡可能小,目標(biāo)函數(shù)為:

(5)

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

λ|a2+b2+c2-1|)

(6)

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

(7)

(8)

對J求偏導(dǎo)得到:

(9)

將(9)式整理后可得:

(10)

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

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

3 本文點(diǎn)云擬合新方法

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

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

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

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

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

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

(6) 計(jì)算所有采樣點(diǎn)到目標(biāo)平面的距離,計(jì)算公式為:

(11)

(7) 求出di的標(biāo)準(zhǔn)差δ,計(jì)算公式為:

(12)

(13)

(8) 當(dāng)di>2δ時(shí),認(rèn)為該點(diǎn)為噪聲點(diǎn),否則保留該點(diǎn)。

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

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

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

4.1 仿真實(shí)驗(yàn)

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

(15)

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

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

表1 未加入異常值的仿真實(shí)驗(yàn)參數(shù)估計(jì)結(jié)果

表2 加入50個(gè)異常點(diǎn)的仿真實(shí)驗(yàn)參數(shù)估計(jì)結(jié)果

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

4.2 實(shí)例分析

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

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

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

圖3 測繪學(xué)院大門點(diǎn)云數(shù)據(jù)

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

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

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

5 結(jié) 論

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