潘 斌,蔡志剛,劉微微,曾玉娟,王 博

(1. 嘉善縣自然資源和規(guī)劃局，浙江 嘉興 314100； 2. 浙江省測繪科學(xué)技術(shù)研究院，杭州 310000； 3.南京航空航天大學(xué)， 南京 210010)

近年來激光雷達(dá)技術(shù)發(fā)展迅速，其可以快速獲取空間三維信息，進(jìn)行大場景的空間探測，廣泛應(yīng)用于城市分析、無人駕駛、電力巡檢等領(lǐng)域，激光雷達(dá)測量結(jié)果以點(diǎn)云形式存儲，但點(diǎn)云數(shù)據(jù)通常密度高、體積大，對數(shù)據(jù)處理提出挑戰(zhàn)。特別是對大規(guī)模無組織點(diǎn)云，構(gòu)建以特征為索引的點(diǎn)云查詢方式能提高查詢效率，對后續(xù)的區(qū)域分割和三維重建等處理有重要意義。三維點(diǎn)云中基于點(diǎn)的索引已經(jīng)有很多研究成果，包括格網(wǎng)索引[1]、k-d樹索引[2]、八叉樹索引[3]和其他混合索引方法[4-7]等?；诰€特征的索引卻很少被關(guān)注，特別是在直線特征豐富的人造場景。線特征既可以是邊界點(diǎn)的統(tǒng)一表現(xiàn)，又可以是平面的交線集合，對點(diǎn)云的表達(dá)和重建工作有著重要意義。進(jìn)行線特征的索引前，需對三維線段進(jìn)行檢測和提取，倪歡等[8]將三維線段的檢測工作主要分為基于圖像的間接法，基于點(diǎn)和基于面的直接法。

基于圖像的方法依靠三維點(diǎn)云對應(yīng)的二維圖像提取特征，或?qū)⑷S點(diǎn)云轉(zhuǎn)化二維圖像，在圖像上提取二維線段，再將結(jié)果重投影至三維空間。X. Ge 等[9]將點(diǎn)云進(jìn)行分割，投影至不同平面，從不同的二維平面中獲取特征點(diǎn)。LI H等[10]利用三維點(diǎn)云對應(yīng)的圖像提取邊緣特征，再與三維點(diǎn)云進(jìn)行匹配確定邊緣點(diǎn)，通過形態(tài)學(xué)方法得到最終特征?；趫D像的方法需要考慮投影策略，構(gòu)建不同平面之間的對應(yīng)關(guān)系，在轉(zhuǎn)化過程中會造成信息的損耗，對三維點(diǎn)云的幾何信息利用不充分。

基于點(diǎn)的方法直接對邊界點(diǎn)進(jìn)行檢測，然后通過線段擬合獲得點(diǎn)云的線特征。陳華偉等[11]根據(jù)曲率突變點(diǎn)隱含特征線的結(jié)論，通過區(qū)域聚類提取備選點(diǎn)集，按照一定的搜索原則進(jìn)行特征線的細(xì)化，在機(jī)械模型和藝術(shù)品模型上取得良好的實(shí)驗(yàn)效果。T. Hackel 等[12]從點(diǎn)的鄰域中提取特征，使用二元分類器預(yù)測每個點(diǎn)輪廓分?jǐn)?shù)作為構(gòu)建候選輪廓的基礎(chǔ)，在候選輪廓中選擇最佳輪廓作為特征。Demantké. J[13]提出了一種多尺度方法PCA計(jì)算激光雷達(dá)點(diǎn)云上的穩(wěn)健幾何特征，在不同半徑大小的球面鄰域上計(jì)算三維特征，基于局部結(jié)構(gòu)張量的特征值組合描述鄰域的形狀，檢出線性局部特征。基于點(diǎn)的方法直接對特征點(diǎn)進(jìn)行處理，因此對噪聲的魯棒性較低。

基于面的檢測方法通過法線或曲率特征將點(diǎn)云聚類分割為不同面，通過面面相交或面上邊緣檢測得到特征線。LIN Y B等[14]首先將點(diǎn)云根據(jù)不同視點(diǎn)轉(zhuǎn)換成陰影圖像的集合，在陰影圖像中提取二維線段及其支撐區(qū)域，再反投影至三維點(diǎn)云空間得到特征線，但該方法只能提取平面交線。孫憲龍等[15]通過投影直方圖進(jìn)行點(diǎn)云分割，利用霍夫變換進(jìn)行直線檢測，對直線進(jìn)行求交獲得建筑物的墻線特征。張帆等[16]將三維點(diǎn)云投影至球面坐標(biāo)系，在全景圖上進(jìn)行邊緣檢測得到包含直線特征點(diǎn)的邊緣點(diǎn)，通過霍夫變換提取直線段。Lin Y等[17]提出基于面的檢測方法提取線特征，通過局部k-means聚類將點(diǎn)云分割成面，再在平面上進(jìn)行邊界點(diǎn)提取和線段分組。Sheng Q等[18]將點(diǎn)云投影至不同擬合平面，對平面進(jìn)行網(wǎng)格化和二值化，在二值化網(wǎng)格平面采用差分編碼的方式提取邊界點(diǎn)，對邊界線特征進(jìn)行差分編碼以便同名直線匹配?；谄矫娴姆椒ㄈ菀酌媾R確定交叉線段端點(diǎn)的問題，也容易出現(xiàn)線段的錯檢和漏檢。

文中針對城市無組織點(diǎn)云數(shù)據(jù)，提出一種城市機(jī)載激光雷達(dá)點(diǎn)云的直線特征提取方法。首先對點(diǎn)云基于區(qū)域生長原理進(jìn)行平面聚類，分割成不同平面，然后將對應(yīng)點(diǎn)云投影至擬合平面，通過穩(wěn)健的二維直線提取方法提取出良好的線特征集合，最后通過空間直線索引方法完成對線特征的索引。

1 方 法

文中方法提取城市點(diǎn)云直線特征方法為：基于區(qū)域增長的點(diǎn)云平面聚類、基于網(wǎng)格的邊界特征提取、建立線特征空間索引。首先計(jì)算法向量和曲率對點(diǎn)云進(jìn)行排序，利用區(qū)域增長算法進(jìn)行聚類，完成點(diǎn)云平面的分割與擬合；然后將點(diǎn)云投影至擬合平面，進(jìn)行網(wǎng)格化和二值化，檢測直線特征并進(jìn)行直線擬合；最終將直線重投影至三維空間，建立空間索引。具體技術(shù)路線如圖1所示。

圖1 技術(shù)路線圖

1.1 基于區(qū)域生長的點(diǎn)云平面聚類

點(diǎn)云的幾何特征信息能表現(xiàn)物體的外輪廓，點(diǎn)云法向量和曲率是點(diǎn)云數(shù)據(jù)幾何特征的代表。首先基于PCA算法對點(diǎn)云法向量進(jìn)行估計(jì)，再利用基于曲率排序的區(qū)域增長算法對點(diǎn)云進(jìn)行聚類并分割，對分割后的區(qū)域進(jìn)行平面擬合，最后對相鄰平面進(jìn)行合并判斷，可以準(zhǔn)確獲取點(diǎn)云平面。

1.1.1 點(diǎn)云法向量和曲率的計(jì)算

曲面上一點(diǎn)的法向量與該點(diǎn)的切平面垂直，點(diǎn)云的法向量估計(jì)一般利用鄰域點(diǎn)構(gòu)建局部曲面，通過計(jì)算局部曲面的法向量得到結(jié)果。

(1)

最小二乘原理擬合局部平面求解法向量，可等價(jià)為求解式(2)的最小值：

(2)

其中，n為局部平面的法向量，它對應(yīng)平面在空間中三個正交基投影方差最小的向量。

P的最小值求解可以轉(zhuǎn)化為計(jì)算協(xié)方差矩陣的特征值，最小特征值對應(yīng)的特征向量即為法向量。

在求解前對所有鄰域點(diǎn)進(jìn)行去重心化以簡便運(yùn)算，新的鄰域點(diǎn)集表示為{pj|(xj,yj,zj),j=1,2,…,N}，其中：

(3)

協(xié)方差矩陣為：

(4)

(5)

為保證法向量的方向全局一致，可以計(jì)算兩個相鄰點(diǎn)云的法向量內(nèi)積，若內(nèi)積為正，則表示同向，內(nèi)積為負(fù)，表示反向，將鄰域點(diǎn)方向調(diào)整為相反。

曲率是指曲線或曲面上點(diǎn)的曲率，在三維點(diǎn)云中，某點(diǎn)鄰域范圍曲率變化越大，則越不光滑，選擇這樣的點(diǎn)作為邊界點(diǎn)的候選點(diǎn)。

曲率計(jì)算式為：

(6)

式中：cu為當(dāng)前點(diǎn)的曲率，λ0,λ1,λ2為協(xié)方差矩陣A奇異值分解求得的特征值。

1.1.2 區(qū)域增長

一般而言，區(qū)域增長方法是根據(jù)同一物體區(qū)域內(nèi)像素的相似性質(zhì)來聚集像素點(diǎn)的方法，從初始區(qū)域(種子點(diǎn))開始，將相鄰的具有相似性質(zhì)的像素或其他區(qū)域歸并到目前的區(qū)域中從而逐步增長區(qū)域，直至沒有可以增長的點(diǎn)為止。在點(diǎn)云中，利用區(qū)域增長方法，選取恰當(dāng)?shù)姆N子點(diǎn)和相似判別依據(jù)，可以提取出目標(biāo)點(diǎn)云中的所有平面。當(dāng)點(diǎn)云曲率越小，該點(diǎn)鄰域變化越小，區(qū)域越平坦，符合平面特征，從平坦區(qū)域開始增長可以有效減少區(qū)域總數(shù)，因此在曲率小的點(diǎn)云中選取種子點(diǎn)。文中相似判別規(guī)則如下：

(7)

若Dplane和Dpo同時小于設(shè)定閾值，則認(rèn)定點(diǎn)pj屬于當(dāng)前區(qū)域。

區(qū)域增長聚類的步驟如下：

1)計(jì)算所有點(diǎn)曲率，并按照曲率從小到大排序，曲率小于一定閾值的點(diǎn)作為種子點(diǎn)。

2)選擇曲率最小點(diǎn)作為首要種子點(diǎn)，利用K鄰近算法搜索鄰域點(diǎn)，如果鄰域點(diǎn)曲率小于等于設(shè)定閾值，則作為新的種子點(diǎn)加入種子點(diǎn)集。若大于設(shè)定閾值，對其進(jìn)行判斷：①種子點(diǎn)與該鄰域點(diǎn)的法向量內(nèi)積結(jié)果大于設(shè)定閾值；②該鄰域點(diǎn)的Dplane和Dpo，同時小于閾值D。滿足上述兩個判斷條件的鄰域點(diǎn)歸并到當(dāng)前種子點(diǎn)的區(qū)域。

3)設(shè)置區(qū)域最小容納點(diǎn)數(shù)和最大容納點(diǎn)數(shù)。對于符合判斷條件的點(diǎn)，將其從原始點(diǎn)云剔除，歸入種子點(diǎn)所在的點(diǎn)集。剩余點(diǎn)云中生成區(qū)域的點(diǎn)數(shù)小于最小容納數(shù)時停止增長，得到區(qū)域分割結(jié)果。

經(jīng)過區(qū)域增長后，可得到點(diǎn)云平面初步分割結(jié)果，但由于鄰近點(diǎn)及閾值選擇存在一定誤差，屬于同一平面的點(diǎn)云可能被誤分割為多個平面，因此需要進(jìn)行相鄰平面的合并判斷。當(dāng)相鄰平面的平面夾角小于規(guī)定閾值T，進(jìn)行合并處理。

1.2 基于網(wǎng)格的邊界特征提取

文中采用平面投影法檢測邊緣特征，首先將區(qū)域點(diǎn)云投影到擬合平面上，進(jìn)行網(wǎng)格化和二值化，然后在二值化圖像上提取邊界特征。

圖2 三維點(diǎn)投影為二維點(diǎn)示意圖

(8)

即可得到空間點(diǎn)投影到擬合平面上的投影坐標(biāo)。

將平面所屬點(diǎn)全部投影至擬合平面后，建立網(wǎng)格對離散點(diǎn)進(jìn)行像素化和二值化。獲取投影點(diǎn)的橫縱坐標(biāo)最大值和最小值，假設(shè)其最大值、最小值分別為xmax,xmin,ymax,ymin。設(shè)定單位網(wǎng)格邊長為s，網(wǎng)格的行數(shù)R和列數(shù)C為:

(9)

其中，y=[x]為取整函數(shù)。

(10)

如圖3所示，將投影點(diǎn)分配完成后，部分網(wǎng)格無投影點(diǎn)存在。由于邊界存在于包含落點(diǎn)的網(wǎng)格中，可對網(wǎng)格構(gòu)建二值圖像掩膜，將有投影點(diǎn)和無投影點(diǎn)的網(wǎng)格進(jìn)行更清晰的區(qū)分：有投影點(diǎn)的網(wǎng)格賦值1；無投影點(diǎn)的網(wǎng)格值賦0，網(wǎng)格轉(zhuǎn)換成僅由像素值為0和1表示的二值圖像。

圖3 網(wǎng)格劃分示意圖

為了增強(qiáng)邊界以及避免提取一些無必要的內(nèi)輪廓，對生成的二值圖像進(jìn)行形態(tài)學(xué)的閉運(yùn)算操作，利用結(jié)構(gòu)元素的原點(diǎn)遍歷二值圖像所有像素點(diǎn)，其中，若結(jié)構(gòu)元素中至少有一個1與原圖重合，則膨脹結(jié)果賦中該像素點(diǎn)賦值為1，否則為0；若重合區(qū)域中結(jié)構(gòu)元素的1均能對應(yīng)于原圖像中的1，則賦值結(jié)果中該像素點(diǎn)賦值為1，否則為0。

在二值化網(wǎng)格平面中進(jìn)行直線檢測與擬合，獲取線段的起始點(diǎn)和結(jié)束點(diǎn)坐標(biāo)，將點(diǎn)坐標(biāo)重投影至三維空間后，得到能夠表達(dá)線段的起始點(diǎn)坐標(biāo)(xs,ys,zs)和結(jié)束點(diǎn)坐標(biāo)(xe,ye,ze)，將所有線段存入線段集中，便于索引與查詢。

1.3 線特征空間索引

先建立點(diǎn)云的固定網(wǎng)格索引，根據(jù)查詢內(nèi)容找到符合條件的線段集合，在網(wǎng)格索引中找到起始點(diǎn)所在網(wǎng)格，以起始點(diǎn)所在網(wǎng)格開始進(jìn)行索引，如圖4所示。假設(shè)網(wǎng)格邊長為la，線段在x,y,z3個方向上的斜率分別為kx,ky,kz，構(gòu)建新的搜索點(diǎn)(xi,yi,zi),即：

圖4 線特征索引示意圖

(11)

xi=xs+kxla,

yi=ys+kyla,

zi=zs+kzla.

(12)

相應(yīng)地下一個搜索點(diǎn)在此基礎(chǔ)上進(jìn)行計(jì)算：

xi+1=xi+kxla,

yi+1=yi+kyla,

zi+1=zi+kzla.

(13)

對搜索點(diǎn)進(jìn)行判別：①判斷搜索點(diǎn)坐標(biāo)是否大于結(jié)束點(diǎn)，大于時，直接搜索結(jié)束點(diǎn)所在網(wǎng)格，記錄網(wǎng)格節(jié)點(diǎn)位置，結(jié)束搜索；②當(dāng)小于結(jié)束點(diǎn)，判斷搜索點(diǎn)是否仍在當(dāng)前網(wǎng)格中，若是則直接計(jì)算下一搜索點(diǎn)替代當(dāng)前搜索點(diǎn)，繼續(xù)進(jìn)行判斷；③若已超出當(dāng)前網(wǎng)格，記錄新網(wǎng)格的節(jié)點(diǎn)位置，繼續(xù)構(gòu)建下一搜索點(diǎn)進(jìn)行判斷。

通過以上流程，可求出線段經(jīng)過的所有網(wǎng)格，完成線特征的空間索引。

2 實(shí)驗(yàn)結(jié)果與分析

為驗(yàn)證線特征索引方法，文中實(shí)驗(yàn)選取浙江省嘉興市嘉善縣地區(qū)的機(jī)載激光雷達(dá)點(diǎn)云進(jìn)行實(shí)驗(yàn)，該點(diǎn)云區(qū)域中包含多種建筑類別，框選其中典型區(qū)域進(jìn)行實(shí)驗(yàn)，點(diǎn)云具體信息如表1所示。實(shí)驗(yàn)環(huán)境為IntelCorei3-6100CPU和4 GB內(nèi)存計(jì)算機(jī)。首先進(jìn)行線特征的初次提取，得到線段集，通過線特征索引查詢線段鄰域點(diǎn)，組成新的點(diǎn)云再次提取線特征，直至滿足迭代條件。并與平面相交提取直線特征的方法進(jìn)行對比。

表1 機(jī)載激光雷達(dá)點(diǎn)云信息

點(diǎn)云影像如圖5(a)、5(c)所示，按高程進(jìn)行灰度賦色，點(diǎn)云缺失處一般為水體或建筑物陰影。區(qū)域1為居民區(qū)，建筑低矮且分布密集；區(qū)域2為高樓，建筑高且分布較為稀疏。利用本文方法進(jìn)行線特征提取，整體結(jié)果如圖5(b)、5(d)所示。

圖5 機(jī)載雷達(dá)點(diǎn)云影像和線特征提取結(jié)果

為直觀表達(dá)線特征提取的結(jié)果，將點(diǎn)云與線特征提取結(jié)果融合顯示如圖6所示。區(qū)域1范圍較大，建筑密集，選取局部區(qū)域進(jìn)行展示，如圖6(a)所示；圖6(b)展示了區(qū)域2提取結(jié)果，區(qū)域2中主要為高層建筑，墻體線較為密集，選取高程不同的三處建筑物進(jìn)行局部放大，展示結(jié)果見圖6(c)。

圖6 線特征和點(diǎn)云結(jié)合顯示結(jié)果

根據(jù)實(shí)驗(yàn)結(jié)果，兩個區(qū)域內(nèi)均提取出分布均勻的建筑物線特征。區(qū)域1為居民區(qū)，建筑低矮且密集分布，房屋類型較為一致。屋頂線特征提取結(jié)果很好，大部分樓棟的屋頂線結(jié)構(gòu)完整，且保留了坡形屋頂?shù)耐暾Y(jié)構(gòu)。但是由于居民區(qū)建筑密集，受限于掃描角度和屋頂遮蔽，機(jī)載雷達(dá)采集的點(diǎn)云在垂直面上的分布較為稀疏，無法形成完整的墻面點(diǎn)云，在二值化網(wǎng)格平面上提取到的邊界特征有限，因此墻體線提取效果弱于屋頂線。區(qū)域2以高樓建筑為主，建筑線特征提取完整，墻體線分布均勻，周邊的低矮建筑屋頂線特征也完整保留。但是由于屋頂結(jié)構(gòu)較為復(fù)雜，短線段特征較多，后處理過程過濾異常線段時將部分短線特征被刪除，從而導(dǎo)致部分建筑屋頂線特征缺失。

為驗(yàn)證文中方法效率，與平面相交法進(jìn)行對比實(shí)驗(yàn)，在兩個區(qū)域內(nèi)隨機(jī)選擇一個小區(qū)域進(jìn)行精確率和召回率的比較。其中精確率表示被正確提取的線段數(shù)占被提取出的線段總數(shù)的比例；召回率表示被正確提取的線段數(shù)占測試區(qū)域所有線段數(shù)的比例。表2展示了平面相交方法與文中方法在兩個區(qū)域中提取到的線段數(shù)、耗時和小區(qū)域的精確率和召回率。圖7展示了兩種方法的對比效果。

圖7 與平面相交法的對比效果

表2 與平面相交法的對比結(jié)果

根據(jù)實(shí)驗(yàn)結(jié)果，文中方法耗時比平面相交法長，在區(qū)域1和區(qū)域2的耗時分別是平面相交法的1.21倍和1.18倍。這是因?yàn)槠矫嫦嘟环▋H需通過判斷平面是否相鄰以及求相鄰平面的交線方程提取線段，而文中方法需要在獲取平面后進(jìn)行二值化網(wǎng)格化、直線提取、線特征索引和迭代等過程。

但文中方法提取到的線段數(shù)遠(yuǎn)高于平面相交法獲取的線段數(shù)，其中區(qū)域1文中方法提取到的線段數(shù)是平面相交法的5.44倍，區(qū)域2文中方法提取到的線段數(shù)是平面相交法的5.58倍。這是由于平面相交法僅能提取兩平面的交線，無法求解屋頂邊線方程，并且在面關(guān)系復(fù)雜的局部可能有失效情況出現(xiàn)，比如區(qū)域1墻面點(diǎn)云缺失，影響交線方程求解，導(dǎo)致平面相交法無法提取到該處的線段特征，區(qū)域2中高層建筑設(shè)計(jì)復(fù)雜，屋頂線輪廓曲折且部分墻面缺失，平面相交法無法取得好的效果。而文中方法通過平面分割投影提取線段，除了邊界直線，還能提取建筑的交面線和邊緣線等線特征，例如區(qū)域2的屋頂線特征，可以通過提取屋頂平面的邊界輪廓擬合線段特征。因此文中方法與平面相交法雖然在精確率上差距較小，但在召回率上文中方法明顯優(yōu)于平面相交法。計(jì)算線段提取效率在區(qū)域1，平面相交法每秒提取7.157條線段，文中方法每秒提取32.141條線段；在區(qū)域2，平面相交法每秒提取4.328條線段，文中方法每秒提取20.437條線段。文中方法在線段的提取效率上遠(yuǎn)高于平面相交法。

為定量評價(jià)本文算法的線特征提取效果，選擇直線擬合精度作為評價(jià)指標(biāo)，直線擬合精度是擬合點(diǎn)到擬合直線段的平均距離誤差，其中擬合點(diǎn)是直線檢測中參與線段擬合的點(diǎn)集。根據(jù)建立的線特征索引結(jié)構(gòu)，隨機(jī)選擇500條直線參與計(jì)算，記錄所有直線中的最大誤差、最小誤差與平均誤差。為確保結(jié)果穩(wěn)定性，對每個區(qū)域隨機(jī)選擇相同數(shù)量的直線進(jìn)行兩次驗(yàn)證實(shí)驗(yàn)，實(shí)驗(yàn)結(jié)果如表3所示。區(qū)域1與區(qū)域2的點(diǎn)平均距離分別為0.123 3 m和0.131 9 m，提取的直線最大誤差與點(diǎn)平均距離相近，平均誤差小于點(diǎn)平均距離的一半，對比之下可以看到文中算法具有較高的線特征提取精度。

表3 直線提取精度對比 m

與基于面的平面相交法進(jìn)行對比，文中提取建筑物的線段特征具有較大優(yōu)勢。與基于點(diǎn)的比較有代表性的PCA方法進(jìn)行對比分析，對比效果如圖8所示，PCA方法能夠提取建筑物傾斜直線上的特征點(diǎn)，但提取效果不佳，且需要人工選取搜索范圍參數(shù)閾值，適用性一般。相比于PCA方法，文中所提算法對建筑物特征點(diǎn)具有良好的提取效果，提取結(jié)果中噪聲較少，可以較好地提取相同方向的相鄰直線上的特征點(diǎn)，精確區(qū)分出鄰近直線段。

圖8 與PCA方法的對比效果

3 結(jié)束語

文中提出了一種城市機(jī)載激光雷達(dá)點(diǎn)云直線特征提取方法。該方法基于區(qū)域增長的點(diǎn)云平面分割和穩(wěn)健的二維線段檢測，并通過三維線段索引方法進(jìn)行迭代處理雜亂線段，最終完成對城市點(diǎn)云的直線特征的提取。與傳統(tǒng)的基于點(diǎn)的方法不同，文中方法可以更直觀地展示城市結(jié)構(gòu)特征，方便后續(xù)建立點(diǎn)、線、面三種幾何表征的關(guān)系。但是文中方法主要基于直線特征進(jìn)行，如何增強(qiáng)地面道路和城市綠化等不規(guī)則線段的提取效果和如何進(jìn)行曲線索引都是今后需要進(jìn)行探索的內(nèi)容。