徐景中,王佳榮
(武漢大學(xué)遙感信息工程學(xué)院,武漢 430079)
(?通信作者電子郵箱jz_xu@whu.edu.cn)
由于掃描視角的限制,利用地基激光掃描設(shè)備進(jìn)行地物目標(biāo)三維點云采集時,通常需要設(shè)置多個測站才能獲取到完整的目標(biāo)點云。由于不同視角采集到的點云并不在同一個坐標(biāo)系下,通常需要將點云旋轉(zhuǎn)、平移到統(tǒng)一的坐標(biāo)系,才能得到完整的地物數(shù)據(jù),因此需要進(jìn)行點云配準(zhǔn)處理。
點云配準(zhǔn)方法目前主要包括:人工配準(zhǔn)方法和自動配準(zhǔn)方法。雖然人工方法通過人為選取同名特征進(jìn)行匹配,能可靠地得到配準(zhǔn)結(jié)果,但精度有限,且耗時耗力;因此,近年來越來越多的學(xué)者開始研究自動配準(zhǔn)方法??紤]到點云與點云之間缺乏嚴(yán)格意義的同名點,點云自動配準(zhǔn)方法又可以分為基于迭代最近點(Iterative Closest Point,ICP)算法的點云配準(zhǔn)方法[1]以及基于特征的點云配準(zhǔn)方法[2-3]。前者是利用數(shù)學(xué)方法迭代尋找兩組點集之間的最小距離作為最佳位置,點對的匹配度高,但其存在對初始點云的位置條件要求苛刻、易陷入局部最優(yōu)等不足。對此,部分學(xué)者從配準(zhǔn)效率入手進(jìn)行ICP 算法的改進(jìn),如:楊帆等[4]在ICP 算法中引入kd-tree 來提高最近點的搜索速度,而周春艷等[5]采用主成分分析(Principal Component Analysis,PCA)變換,Du 等[6]采用仿射法等方法,這些研究也只是解決了部分問題,仍存在配準(zhǔn)完成度低的問題。而基于特征的點云配準(zhǔn)方法則是從點云中提取點特征、線特征或面特征,然后通過特征匹配完成點云配準(zhǔn)。如:沈長江等[7]、趙明富等[8]直接從點云中提取快速點特征直方圖(Fast Point Feature Histogram,F(xiàn)PFH),然后基于FPFH 特征相對關(guān)系進(jìn)行點云配準(zhǔn);但此類基于點特征的配準(zhǔn)通常存在受噪聲影響較大、時間復(fù)雜度較高等問題,因此,部分學(xué)者通過構(gòu)建點的域拓?fù)湫畔硖崛缀翁卣鳎缓蠡谟蛱卣鲀?yōu)化點云配準(zhǔn)過程[9],在一定程度上改善了對噪聲的敏感性,但仍存在魯棒性低的問題。Aiger 等[10]基于共面四點對的仿射不變性,實現(xiàn)了四點全等集合(4-Points Congruent Sets,4PCS)點云配準(zhǔn)方法,該方法不需要計算復(fù)雜的幾何特征,效率較高;但對于點云重疊率較低的情況,算法耗時長且配準(zhǔn)容易失敗。相較于點特征,線特征幾何拓?fù)浞€(wěn)定性更高,配準(zhǔn)結(jié)果也更加可靠。Jaw 等[11]提出了一種基于三維直線特征的地面激光雷達(dá)點云配準(zhǔn)方法,該方法包括三個主要步驟:基于平面擬合相交的三維直線特征提取、基于直線夾角的相似性匹配和基于四元數(shù)的多站激光點云配準(zhǔn)。王永波等[12]通過平面相交得到線特征并將其作為激光點云配準(zhǔn)的基元,利用線特征約束的四元數(shù)法構(gòu)建參數(shù)解算模型。此類方法可以快速實現(xiàn)點云配準(zhǔn),但由于提取的線特征的完整性和精確性是有限的,配準(zhǔn)結(jié)果精度仍需提高。由于面特征包含的信息多于點、線特征,面特征配準(zhǔn)方法受噪聲影響較小,部分學(xué)者使用最小二乘法、隨機樣本一致性算法[13]和主成分分析方法[14]等曲面擬合方法進(jìn)行特征面提取,以提高配準(zhǔn)精度;但此類方法,要求兩組數(shù)據(jù)重疊區(qū)域必須包含許多面特征,否則很難保證配準(zhǔn)的準(zhǔn)確性。
基于上述分析,本文提出了一種基于線特征及ICP 算法的建筑物點云自動配準(zhǔn)方法,該方法包括基于線特征的點云粗匹配以及基于ICP 的點云精匹配過程。該方法利用線對作為配準(zhǔn)基元,以線對夾角和距離作為相似性測度,可以提高配準(zhǔn)結(jié)果的可靠性;此外,采用由粗到精的配準(zhǔn)策略,不僅能改善ICP 算法易陷入局部最優(yōu)的不足,同時可以保證點云配準(zhǔn)的精度。
點云平面分割,主要是基于點云局部特征的一致性進(jìn)行區(qū)域增長,將建筑物墻面、屋頂區(qū)域分割成多個同質(zhì)區(qū)域。考慮到待分割目標(biāo)為平面上的點,為了快速得到建筑物點云的平面區(qū)域,本文首先利用建筑物點云構(gòu)建kd-tree[15],在此基礎(chǔ)上通過PCA 變換分析方法估算點云的法向信息[16],然后從任一種子點開始,以法向一致性進(jìn)行鄰域點的判斷,實現(xiàn)平面區(qū)域增長。區(qū)域增長條件可表示為:
其中:pi為種子點;pj為鄰域點;Npi與Npj分別為兩點處的法向量;thdvv為頂點距離閾值;thθ為法向夾角閾值。
由于區(qū)域增長過程中種子點的選擇具有一定的隨機性,上述區(qū)域增長結(jié)果中,通常會存在同一個平面被分割成多塊的現(xiàn)象,因此,需要進(jìn)行區(qū)域合并處理。合并過程如下:對于平面Ri中的每個點pi,遍歷其鄰域中所有點,如果其鄰域點所屬的平面號和當(dāng)前點pi不同,則認(rèn)為該鄰域點所屬的平面為當(dāng)前點所在平面的相鄰平面,若二者法向夾角小于夾角閾值thθ,則進(jìn)行合并處理。
1.1.2 特征線提取及優(yōu)化
雖然建筑物點云不僅包含屋頂區(qū)域,還包括墻面、窗戶等區(qū)域,經(jīng)過平面分割處理,各部分基本得到有效分割,為了提取不同區(qū)域點簇的特征線,本文直接采用Alpha-shape 算法[17]進(jìn)行區(qū)域輪廓線的提取,在此基礎(chǔ)上通過拆分處理得到平面區(qū)域特征線。Alpha-shape輪廓線提取過程如下:
1)對于平面點簇S,設(shè)定半徑閾值α。
2)選取S中任意兩點pi、pj,若其連線的長度小于2*α,則根據(jù)幾何關(guān)系可計算過pi、pj圓周的圓心。
3)若點簇S中存在與圓心的距離小于α的點,則表明pi、pj非輪廓點;否則標(biāo)記pi、pj為輪廓點,其連線即為輪廓線段。
4)繼續(xù)按步驟2)~3)搜索判斷其他點,直到所有點均被處理時停止。
為了便于后續(xù)特征線匹配,本文對輪廓線提取結(jié)果進(jìn)行分段和擬合處理得到特征線段。分裂時,利用輪廓線上距離首尾端點所連直線的最大偏差點,將輪廓線拆分成多個子線段;然后,采用最小二乘法擬合方法[18]對輪廓點進(jìn)一步擬合得到特征線段。
考慮到激光掃描的隨機性以及點簇噪聲的影響,上述特征線提取結(jié)果中,通常會存在因點云分布不均勻或小平面導(dǎo)致的冗余特征線,需要進(jìn)一步進(jìn)行優(yōu)化處理。
1)短線剔除處理。若線段長度小于閾值Th_L,或相鄰線段之間角度小于閾值Th_θ,則剔除短線。
2)冗余線段合并。如圖1 所示,若線段Li與線段Lj滿足如下條件,則進(jìn)行合并處理。
其中:dLi和dLj分別是從原點到Li和Lj所在直線的距離;L是投影區(qū)域長度;ε為比例系數(shù);max{dprj}為Lj的兩個端點與線段Li的垂直距離最大值;davg為該組點云的平均點間距。
合并處理時,將Lj端點投影在線段Li上,并取4 個端點中距離最遠(yuǎn)的兩個點更新Li的端點:若Li距離較遠(yuǎn)的端點包含線段Lj投影端點,則將Lj的投影端點作為Li新的端點,圖中L*線段為線段合并后得到的新線段。
圖1 冗余線段合并示意圖Fig.1 Schematic diagram of redundant line segment merging
1.2.1 配準(zhǔn)基元和相似性測度
三維空間中,若要解算兩組點云的有效變換矩陣,至少需要兩對非平行直線。為了保證匹配的可靠性,本文針對上述優(yōu)化處理的結(jié)果,進(jìn)一步篩選有效的特征線對作為配準(zhǔn)基元。若相鄰線段滿足以下條件,則該配準(zhǔn)基元有效。
如圖2所示,d為相鄰線段最鄰近端點的距離,θ表示兩條相鄰線段的夾角。
圖2 配準(zhǔn)基元示意圖Fig.2 Schematic diagram of registration primitive
而同名特征的判斷,以線對夾角和線段長度作為相似性測度。對于兩組點云p和q,L1表示線對中長線段的長度,L2表示線對中短線段的長度,dp和dq分別為兩組點云的平均點間距。如果兩組點云中的線對滿足式(4),則將其視為同名線對。
1.2.2 點云粗配準(zhǔn)
根據(jù)相似性測度,在兩個待配準(zhǔn)點集中尋找同名線對,并以線對的端點作為同名點,解算最終變換模型。為得到最優(yōu)變換關(guān)系,本文采用如下方法迭代遍歷兩組點集中的線對,得到最佳的匹配結(jié)果:
1)基于兩組點集中的配準(zhǔn)基元,利用相似性測度在兩組點集中查找匹配同名線對,并以線對端點坐標(biāo),采用四元數(shù)法[19]解算初始變換矩陣。
2)利用初始變換矩陣處理待配準(zhǔn)點集中的線特征,將待配準(zhǔn)點集轉(zhuǎn)換至目標(biāo)點集中。
3)統(tǒng)計變換后的待配準(zhǔn)點集與目標(biāo)點集中同名線段的個數(shù)。
4)重復(fù)步驟1)~3),逐個遍歷點集中的配準(zhǔn)基元,得到變換后的待配準(zhǔn)點集與目標(biāo)點集中同名線對的個數(shù),取同名線對數(shù)目最多時的變換矩陣作為最優(yōu)匹配結(jié)果。
5)為了使剛性變換矩陣解算結(jié)果更加準(zhǔn)確,利用所有同名線段的端點坐標(biāo)解算變換矩陣,并以此作為點云粗配準(zhǔn)結(jié)果。
1.2.3 點云精配準(zhǔn)
考慮到基于點云的特征線的提取會存在一定的誤差,因此本文在基于特征線匹配的基礎(chǔ)上,進(jìn)一步采用ICP 算法進(jìn)行兩組點集的迭代配準(zhǔn),以提高配準(zhǔn)結(jié)果的精度,步驟如下:
1)對于粗配準(zhǔn)點集P和Q,在點集P中尋找與Q對應(yīng)的點,生成對應(yīng)點對。2)通過四元數(shù)法計算旋轉(zhuǎn)平移的6個參數(shù)及變換矩陣R。3)對點集P使用上述變換矩陣R進(jìn)行旋轉(zhuǎn)和平移,得到新的點集P'。
根據(jù)點集P'和Q中的對應(yīng)點的平均距離來確定是否繼續(xù)迭代。因此,求解最佳變換矩陣的問題,可以轉(zhuǎn)化為計算使f(R,T)的值最小時的旋轉(zhuǎn)矩陣R和平移矩陣T;
4)若繼續(xù)迭代,更新待配準(zhǔn)點云的位置,繼續(xù)步驟1);否則結(jié)束迭代,得到最終的配準(zhǔn)結(jié)果。
為了驗證本文算法的有效性,采用了兩組具有不同重疊度的建筑物點云數(shù)據(jù)進(jìn)行實驗。兩組數(shù)據(jù)均是利用地面激光雷達(dá)設(shè)備從不同角度采集的建筑物點云,為部分重疊點云,其中點集Ⅱ的重疊區(qū)域較小,受遮擋嚴(yán)重,點云渲染結(jié)果如圖3所示。兩組點云的統(tǒng)計信息如表1 所示,由表1 可以看出:點集Ⅰ的最大平均間距為0.053 m,點云密度較為稀疏,而點集Ⅱ的最大平均點間距為0.019 m,點云相對稠密。
圖3 實驗區(qū)域點云渲染效果圖Fig.3 Rendering effect diagrams of point clouds in experimental regions
表1 點云信息統(tǒng)計Tab.1 Statistical information of point clouds
圖4 是兩組點云的平面分割結(jié)果,為了便于說明,此處只針對每個數(shù)據(jù)集的左站點云處理結(jié)果進(jìn)行展示分析。圖4中,圖(a)、(c)是直接基于點云法向一致性進(jìn)行區(qū)域增長的結(jié)果,從圖中可以看出,很多共面區(qū)域被分成了多個小面片,由于該方法的區(qū)域增長范圍被限制在該平面初始種子點的周圍,因此提取得到的小平面呈現(xiàn)出以初始種子點為圓心的圓形區(qū)域,對此需要進(jìn)行區(qū)域合并處理;圖(b)、(d)為基于平面法向的平面區(qū)域合并結(jié)果,可以看出,位于同一平面上的多個小平面已經(jīng)被正確合并,形成較為完整的墻面,墻面結(jié)構(gòu)更加清晰。
圖5 是兩組點集的直線段提取結(jié)果,其中:圖(a)、(c)是直接基于投影點集提取的直線段,圖(b)、(d)為優(yōu)化后的直線段結(jié)果(Th_L=0.5,Th_θ=10,ε=0.05)。從圖5(a)、(c)中可以看出,兩組點云中均提取了豐富的結(jié)構(gòu)信息,提取的結(jié)果能基本反映出建筑物表面的細(xì)節(jié)信息,但由于平面分割等因素,特征線中存在冗余線段;從圖5(b)、(d)優(yōu)化結(jié)果可以看出,優(yōu)化后冗余線段被有效合并,短線被正確剔除,特征線更加簡潔,能清晰表達(dá)建筑物屋頂及墻面特征。
基于特征線提取結(jié)果進(jìn)行同名線段的匹配,匹配結(jié)果如表2 所示。從表2 中可以看出,數(shù)據(jù)集I 中大部分特征線得到有效匹配,而數(shù)據(jù)集Ⅱ由于重疊區(qū)域較少,雖然左右測站具有豐富的線特征,但尋找到的同名線段數(shù)目較少,只有23 對正確匹配的同名線段。
圖4 點云平面分割結(jié)果Fig.4 Plane segmentation results of cloud points
圖5 特征線段提取結(jié)果Fig.5 Extraction results of feature line segments
表2 線段提取及匹配結(jié)果統(tǒng)計Tab.2 Statistics of extraction and matching results of feature line segments
圖6 為建筑物點云配準(zhǔn)結(jié)果圖,圖(a)為數(shù)據(jù)集Ⅰ的配準(zhǔn)結(jié)果,圖(b)為數(shù)據(jù)集Ⅱ的配準(zhǔn)結(jié)果。從圖6 中可以看出,兩組數(shù)據(jù)集重疊區(qū)域均正確地重疊在一起。
圖6 點云配準(zhǔn)結(jié)果Fig.6 Registration results of point clouds
圖7 為配準(zhǔn)結(jié)果局部剖面圖,從圖7 中可以更清楚地看出,不同測站的墻面點云較好地吻合在一起,表明實驗區(qū)域點云得到較好的配準(zhǔn)。
圖7 配準(zhǔn)結(jié)果局部剖面圖Fig.7 Part sectional drawings of registration results
為了進(jìn)一步驗證方法的性能,將本文方法配準(zhǔn)結(jié)果與經(jīng)典ICP 算法結(jié)果進(jìn)行對比分析。由于兩組數(shù)據(jù)集均為部分重疊,經(jīng)典ICP 算法不能收斂到正確位置。本文以手動給定位置初值的ICP配準(zhǔn)結(jié)果作為參照,分析經(jīng)典ICP算法以及本文方法的配準(zhǔn)結(jié)果,并統(tǒng)計配準(zhǔn)結(jié)果對應(yīng)點之間的最大距離誤差、平均距離誤差以及均方根誤差(Root Mean Square Error,RMSE)。RMSE計算式如下:
其中:g是點集數(shù)量;di是本文方法與手動配準(zhǔn)結(jié)果的對應(yīng)點的距離。
點云配準(zhǔn)誤差統(tǒng)計分析結(jié)果如表3 所示。從表3 中可以看出,經(jīng)典ICP 算法配準(zhǔn)結(jié)果與手工初值的配準(zhǔn)結(jié)果具有較大誤差,特別對于數(shù)據(jù)集Ⅱ,RMSE 達(dá)到0.245 m,表明對于部分重疊的建筑物點云,經(jīng)典ICP 算法難以收斂到正確的位置;而本文方法結(jié)果與經(jīng)典ICP 算法結(jié)果相比,配準(zhǔn)精度有了明顯提高,兩組數(shù)據(jù)集的平均誤差分別為0.042 m 和0.044 m,表明基于特征的粗配準(zhǔn)方法與手工給定初值的配準(zhǔn)方法得到的結(jié)果精度相當(dāng),基于特征的粗配準(zhǔn)能為ICP 配準(zhǔn)提供有效的初值位置,提高ICP算法的穩(wěn)定性。
表3 點云配準(zhǔn)誤差統(tǒng)計分析 單位:mTab.3 Statistics and analysis of registration errors of point clouds unit:m
此外,對比點云信息統(tǒng)計表1以及點云重疊情況(見圖4)可以發(fā)現(xiàn),在平均點間距方面,數(shù)據(jù)集Ⅰ大于數(shù)據(jù)集Ⅱ,但二者的RMSE 相當(dāng),分別為0.051 m 和0.053 m,分別為其點云最大平均點間距的0.96倍和2.79倍,表明點云配準(zhǔn)精度不僅與點云密度相關(guān),還受到點集重疊區(qū)域大小影響,點云密度越高,重疊度越大,配準(zhǔn)結(jié)果的可靠性越高。
地基建筑物點云配準(zhǔn)對于獲取完整的建筑物信息、構(gòu)建建筑物三維模型具有重要意義。本文在分析點云配準(zhǔn)現(xiàn)狀的基礎(chǔ)上,根據(jù)優(yōu)勢互補的思想,提出了一種基于線特征及ICP算法的地基建筑物點云配準(zhǔn)方法。該方法在建筑物點云平面分割的基礎(chǔ)上,通過輪廓線檢測以及分裂擬合處理,得到建筑物豐富的特征線;然后,以特征線對作為配準(zhǔn)基元,根據(jù)線段夾角及距離等相似性測度實現(xiàn)點云的粗匹配,在此基礎(chǔ)上利用ICP 算法完成建筑物點云精確配準(zhǔn)。通過兩組實地數(shù)據(jù)進(jìn)行實驗驗證,實驗結(jié)果表明,基于線特征的由粗到精的點云配準(zhǔn)方法能改善ICP 算法對點云初值的依賴,特別是可以克服ICP 算法針對部分重疊區(qū)域點云易陷入局部最優(yōu)的不足,提高了配準(zhǔn)方法的穩(wěn)定性。下一步研究將集中在地物其他特征的提取及匹配方面,以降低方法對線特征的依賴,提高其通用性。