基于機(jī)載LiDAR數(shù)據(jù)的建筑物屋頂邊界提取

朱軍桃,王 雷,趙苗興,龔朝飛,羅 樂

(桂林理工大學(xué) a.測繪地理信息學(xué)院;b.廣西空間信息與測繪重點(diǎn)實(shí)驗(yàn)室,廣西 桂林 541006)

隨著機(jī)載激光雷達(dá)測量技術(shù)的不斷完善, 其廣闊應(yīng)用前景受到人們的極大關(guān)注, 機(jī)載LiDAR技術(shù)已滲入到生產(chǎn)、 生活的各個領(lǐng)域。 數(shù)字地球及數(shù)字城市建設(shè)產(chǎn)生了構(gòu)建建筑物三維模型的極大需求, 而對于構(gòu)建建筑物三維模型來說, 需要對建筑物邊界進(jìn)行精確提取,才能提高三維建模的精度?，F(xiàn)今,對建筑物邊界提取的方法主要分為兩類: 1)結(jié)合其他數(shù)據(jù)源與點(diǎn)云數(shù)據(jù)對建筑物邊界進(jìn)行提取[1-3](此類方法對數(shù)據(jù)類型要求高)； 2)僅使用LiDAR點(diǎn)云數(shù)據(jù)進(jìn)行建筑物邊界提取。 李江雄[4]通過將散亂點(diǎn)云網(wǎng)格化, 獲得邊界網(wǎng)孔并求得每個邊界網(wǎng)孔中的最小凸包來提取點(diǎn)云邊界, 該方法對點(diǎn)云密度要求高, 對機(jī)載LiDAR數(shù)據(jù)限制性強(qiáng);曾齊紅等[5]通過構(gòu)建TIN模型來判斷點(diǎn)云邊界, 該方法由于受到凸包性限制, 可能使邊緣點(diǎn)提取不完整; 沈蔚等[6]采用Alpha Shape算法提取建筑物邊界, 該方法可以避免因凸包性造成的影響, 但判別時要遍歷所有激光點(diǎn), 同時要求解圓心, 導(dǎo)致程序運(yùn)行效率低; 符小俐等[7]采用帶距離控制的卷包裹算法提取建筑物邊緣, 此方法對距離控制要求很高, 且提取邊精度不高。本文提出了一種基于機(jī)載點(diǎn)云數(shù)據(jù)提取建筑物屋頂邊界的方法，與基于不規(guī)則三角網(wǎng)模型的邊緣提取方法相比，獲得了更多邊緣點(diǎn)，且避免了由于凸包性導(dǎo)致對建筑物邊緣提取造成關(guān)鍵點(diǎn)遺漏問題，對多種形狀的建筑物適應(yīng)性強(qiáng)，運(yùn)行效率高。

1 原理與方法

現(xiàn)實(shí)中,絕大多數(shù)房屋都屬于精心設(shè)計過的、轉(zhuǎn)角為直角的規(guī)則或不規(guī)則多邊形。針對此特點(diǎn),本文提出了一種直接采用點(diǎn)云數(shù)據(jù)提取房屋屋頂邊線的方法。 首先,以細(xì)小網(wǎng)格劃分點(diǎn)云,從包圍網(wǎng)的4個方向?qū)χ付ǚ课葸M(jìn)行邊緣點(diǎn)初步提取,其結(jié)果基本反映了建筑物的具體形狀;然后,采用改進(jìn)的帶距離控制的卷包裹算法提取建筑物角點(diǎn)，并依據(jù)相鄰點(diǎn)之間的關(guān)系二次提取邊緣點(diǎn),提高提取精度;最后,對邊緣點(diǎn)進(jìn)行最小二乘直線擬合、規(guī)則化及擴(kuò)展，重現(xiàn)了建筑物邊界線,具體流程見圖1。

圖1 基于機(jī)載LiDAR點(diǎn)云數(shù)據(jù)的房屋輪廓線提取流程Fig.1 Extraction process of building outlines based on airborne LiDAR data

1.1 邊界初步提取

將屋頂點(diǎn)云數(shù)據(jù)投影至XOY平面, 對屋頂面投影的二維點(diǎn)集構(gòu)建囊括點(diǎn)云簇的四邊形(大于最小外包四邊形); 以邊長接近點(diǎn)云平均間距的正方形網(wǎng)格將四邊形分割, 分割網(wǎng)格能分割所有點(diǎn)云，并使每個網(wǎng)格內(nèi)具有1～2個激光點(diǎn)。 至此, 每個網(wǎng)格近似代表了一個激光點(diǎn), 并且相鄰網(wǎng)格間沒有間隙, 使該分割結(jié)果確保對初始邊界搜索時遍歷所有的點(diǎn)云。 四邊形內(nèi)部的網(wǎng)格分為兩類: 一類是內(nèi)部有激光點(diǎn)，稱為“實(shí)孔”; 另一類是內(nèi)部沒有激光點(diǎn)，稱為“空孔”。 以四邊形某條邊的最外側(cè)網(wǎng)孔開始, 對垂直于四邊形邊界的某一列網(wǎng)孔順序判斷: 如果第一個網(wǎng)孔為空孔，則將空孔計數(shù)器賦值為1, 如果下一個網(wǎng)孔仍為空孔，計數(shù)器增加1, 按序逐一判斷； 如果標(biāo)記空孔的計數(shù)器超過某一閾值, 且遇到下一個網(wǎng)孔為實(shí)孔則認(rèn)為其內(nèi)部的點(diǎn)為邊緣點(diǎn), 并將計數(shù)器賦值為0; 若連續(xù)出現(xiàn)空孔數(shù)量并未超過閾值,不記錄下次出現(xiàn)實(shí)孔內(nèi)的點(diǎn),并將計數(shù)器重新賦值為0。按照此方法,從4個方向逐次判斷所有網(wǎng)孔。如圖2所示, “●”為按箭頭方向搜索到的邊緣點(diǎn)。由于點(diǎn)云分布并不完全均勻,若內(nèi)部出現(xiàn)較大間距,可能將屬于建筑物內(nèi)部的非邊緣點(diǎn)標(biāo)記為建筑物的邊緣點(diǎn),所以在所有方向的邊緣點(diǎn)提取完成之后,刪除距離其他邊緣點(diǎn)較遠(yuǎn)的孤立點(diǎn)和數(shù)量極少的點(diǎn)簇。

圖2 邊緣點(diǎn)初步提取Fig.2 Preliminary extraction of edge points

1.2 角點(diǎn)提取

在初步提取邊緣點(diǎn)云后,由于點(diǎn)云分布的離散性會導(dǎo)致初步提取的相鄰邊緣點(diǎn)間出現(xiàn)嚴(yán)重鋸齒現(xiàn)象,該現(xiàn)象導(dǎo)致傳統(tǒng)依賴相鄰激光點(diǎn)之間夾角來判斷邊緣角點(diǎn)的方法失效。帶距離控制的卷包裹算法是從某一邊界上的固定點(diǎn)開始，并以一定長度旋轉(zhuǎn),觸碰到的第一個點(diǎn)必是外包多邊形的一個頂點(diǎn),并以此點(diǎn)為固定點(diǎn)繼續(xù)旋轉(zhuǎn),直到回到起始點(diǎn)。方法步驟具體為[8]: 1)在點(diǎn)集中依次找出X和Y坐標(biāo)最大點(diǎn)與最小點(diǎn)，取其中一點(diǎn)作為出發(fā)點(diǎn)，記為P1； 2)從P1出發(fā)作一條射線，順時針或逆時針方向旋轉(zhuǎn)，直至觸碰到點(diǎn)集中某點(diǎn)，記為P2； 3)重復(fù)步驟2)，若觸碰點(diǎn)回到P1，則結(jié)束。

上述方法獲取的邊緣點(diǎn)密度較低且無法提取凸包邊緣。 本文改進(jìn)了帶距離控制的卷包裹算法, 并采用該算法進(jìn)行角點(diǎn)提取, 過程如下: ① 順序連接1.1節(jié)中提取的邊緣點(diǎn), 將邊緣點(diǎn)所圍成圖形分為內(nèi)外兩側(cè); ② 在初步提取的邊緣點(diǎn)中尋找距離最大的兩個點(diǎn), 認(rèn)為這兩個點(diǎn)一定屬于邊界上的點(diǎn); ③ 以其中一點(diǎn)P開始，沿點(diǎn)云外側(cè)定義一條方向邊,并以此邊開始以一定長度L旋轉(zhuǎn), 在長度L范圍內(nèi)碰到第一個點(diǎn)P1時停止旋轉(zhuǎn)并記錄P1; ④ 以P與P1所連直線為方向邊,再次以P1為中心，長度L為半徑旋轉(zhuǎn),碰到第一個點(diǎn)P2時停止旋轉(zhuǎn),并記錄夾角α。當(dāng)Pn-1、Pn、Pn+1所構(gòu)成的夾角αn∈180°±θ時(θ為允許角度變化的范圍),Pn為邊線上的點(diǎn); 當(dāng)αn∈ 270°±θ時,Pn為角點(diǎn)(圖3); 當(dāng)αn?(180°±θ)∩αn?(270°±θ)時, 將長度L增加D(D為所增加的長度), 再次旋轉(zhuǎn), 若觸碰點(diǎn)改變, 則此點(diǎn)為邊線上的點(diǎn),以此處為起始繼續(xù)步驟③(Pn+2點(diǎn)處); ⑤ 重復(fù)步驟④, 直到回到起始點(diǎn)P處; ⑥ 以點(diǎn)P開始從點(diǎn)云內(nèi)側(cè)定義一條方向邊,繼續(xù)步驟③～⑤。

卷包裹算法只能判斷具有“凸”性質(zhì)的角點(diǎn),而對內(nèi)“凹”性質(zhì)的角點(diǎn)無法判斷。因初步提取的建筑物邊緣點(diǎn)已形成了閉合環(huán)狀點(diǎn)云,從內(nèi)側(cè)進(jìn)行改進(jìn)算法是將凹角點(diǎn)轉(zhuǎn)化為凸角點(diǎn)的處理過程。同時,步驟④中對長度L增加D，避免了從內(nèi)側(cè)旋轉(zhuǎn)過程中角點(diǎn)的錯誤判斷。

1.3 邊緣點(diǎn)二次提取

初步提取中得到的邊緣點(diǎn)呈現(xiàn)鋸齒狀主要原因是:(1)點(diǎn)云數(shù)據(jù)本身的特性,由于激光點(diǎn)云是對物體表面的一種不規(guī)則分布點(diǎn)位置的描述,所以造成鋸齒現(xiàn)象與多種因素有關(guān),處于最外側(cè)網(wǎng)孔的激光點(diǎn)未必都為邊緣點(diǎn);(2)由于網(wǎng)孔大小和點(diǎn)云密度分布不均的影響,導(dǎo)致每個最外側(cè)網(wǎng)孔中的激光點(diǎn)未必只有一個，所以需要對邊緣點(diǎn)進(jìn)行二次提取。依據(jù)建筑物的相鄰邊緣點(diǎn)間所連線段某側(cè)不存在激光點(diǎn)的性質(zhì),提取相鄰角點(diǎn)間的邊緣點(diǎn),此邊緣點(diǎn)集合同屬于一條房屋邊線上。二次提取方法為:將兩個角點(diǎn)存入集合list中,獲得兩角點(diǎn)連線的斜率K,向邊緣點(diǎn)外側(cè)任意距離構(gòu)建一條斜率為K的直線S,將此相鄰角點(diǎn)及角點(diǎn)間的邊緣點(diǎn)投影至S上,設(shè)角點(diǎn)A投影點(diǎn)為A′,距其最近的投影點(diǎn)B′,次近的投影點(diǎn)為C′，分別計算與投影點(diǎn)所對應(yīng)的角點(diǎn)A,邊緣點(diǎn)B、C到直線S的距離DIS:

圖3 角點(diǎn)提取示意圖Fig.3 Sketch map of corner points extraction

若DISA≥DISB且DISB

若DISA≥DISB且DISB>DISC,將B存入集合list中,并以C點(diǎn)作為下次判斷起始;

若DISADISC,剔除B點(diǎn),并以C點(diǎn)作為下次判斷起始點(diǎn);

若DISA

重復(fù)上述方法判斷兩角點(diǎn)間所有邊緣點(diǎn),并根據(jù)邊緣點(diǎn)的粗糙程度及數(shù)量將集合list中的邊緣點(diǎn)設(shè)置合適的迭代次數(shù)進(jìn)行迭代。此方法是將兩角點(diǎn)連線模擬為房屋邊界線,逐步剔除不屬于邊緣點(diǎn)的激光點(diǎn),采用上述方法可以極大地弱化兩個角點(diǎn)所連直線與真正房屋邊界線斜率的差異所造成同一邊線上邊緣點(diǎn)錯誤提取的問題。

圖4為某建筑屋頂?shù)倪吘夵c(diǎn)提取結(jié)果圖,其中圖4a為投影到XOY平面的建筑物點(diǎn)云, 圖4b為初步提取的邊緣點(diǎn)集, 圖4c為二次提取的邊緣點(diǎn)集。

可以看出， 經(jīng)二次提取后其邊界粗糙程度明顯降低。

1.4 房屋邊線擬合、規(guī)則化及擴(kuò)展

針對經(jīng)過二次提取后的邊緣點(diǎn)集, 采用最小二乘方法將角點(diǎn)間的邊緣點(diǎn)擬合為直線。 所擬合的邊界線雖然已經(jīng)能夠很好地反映出建筑物的具體形狀, 但相鄰的邊界線仍存在不垂直的關(guān)系, 為了更好地反映屋頂?shù)恼鎸?shí)形狀, 將各相鄰邊界線規(guī)則化為直角, 并對規(guī)則化后的邊線進(jìn)行擴(kuò)展處理。 本文采用文獻(xiàn)[9]所述的擴(kuò)展方法, 計算某規(guī)則化后的邊線與參與此邊線擬合的外側(cè)邊緣點(diǎn)的距離, 并將此邊平移至距離最大處,最大程度地還原建筑物的原本大小及形狀。

圖4 邊緣點(diǎn)提取結(jié)果Fig.4 Extraction results of edge points

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

由于實(shí)驗(yàn)區(qū)域地勢平坦且建筑物較高, 因此屋頂點(diǎn)云的獲取直接根據(jù)點(diǎn)云的高程分布[7],高程分布密度最大區(qū)域?yàn)榈孛纥c(diǎn)云所在集合, 去除此部分及以下點(diǎn)云數(shù)據(jù)得到建筑物點(diǎn)云, 使用CloudCompare軟件剔除墻面點(diǎn), 最終得到建筑物屋頂點(diǎn)云。

為保證獲得外側(cè)邊緣的“實(shí)孔”,所構(gòu)外包四邊形邊長為Lt=Lmin t+Gli×Gle+1,其中:Lt為外包四邊形邊長,Lmin t為最小外包四邊形邊長,Gli為空孔閾值,Gle為網(wǎng)孔邊長。

實(shí)驗(yàn)參數(shù):編譯工具為Matlab R2014b，點(diǎn)云平均密度為3.4個/m2，網(wǎng)孔邊長0.5 m，空孔閾值為6。

圖5為提取某地區(qū)高程230～260 m、260～308 m、308～325 m具有內(nèi)外邊界的不規(guī)則多邊形(圖5a)、 內(nèi)轉(zhuǎn)角為直角的不規(guī)則多邊形(圖5b)及規(guī)則四邊形(圖5c)3種最具代表性建筑物。左側(cè)房屋內(nèi)邊界中曲線部分仍按本文方法進(jìn)行二次提取,但依然完好地保留邊界特征。在擬合相鄰角點(diǎn)間的直線時,對于距離非常近的角點(diǎn)不進(jìn)行擬合,圖5b中圓孔非房屋內(nèi)邊界,因此未進(jìn)行邊界擬合。

對各二次提取設(shè)置不同迭代次數(shù)后所得到結(jié)果作比較:

(1)

式中: (xi,yi)為屬于同一邊線上邊緣點(diǎn)集合中的點(diǎn);f(xi)為xi處經(jīng)最小二乘擬合的值;Δ值反映邊緣點(diǎn)在所擬合直線兩側(cè)分布的離散程度。

表1各邊號對應(yīng)圖5中所在邊,可以看出絕大多數(shù)邊在增加迭代次數(shù)后Δ值會降低,但少量數(shù)據(jù)存在異常：4、9、13、19、28、32、35號邊Δ值在降低到一定程度就不再降低,造成這種情況的原因是邊緣點(diǎn)集形成了一種光滑的、不再具有鋸齒狀的形狀; 2、 14、 18、 21、 25、 27、 39號邊Δ值在降低到一定程度后增大了, 造成這種情況的原因是設(shè)置太高的迭代次數(shù)過濾掉過多的邊緣點(diǎn)。 上述兩種情況均出現(xiàn)在一些長度較短的邊上, 此類邊緣增加迭代次數(shù)后最終因太少的邊緣點(diǎn)影響了擬合效果, 同時因掃描線方向、 點(diǎn)云分布不均等原因造成初始迭代的Δ值各不相同。 經(jīng)試驗(yàn), 本文增加迭代次數(shù)找到最低Δ值, 并增加限制條件:邊緣點(diǎn)最少不能小于0.5個/m, 該條件在保證足夠數(shù)量的邊緣點(diǎn)基礎(chǔ)上達(dá)到最佳提取效果。 最終, 規(guī)則化及擴(kuò)展后的屋頂邊線如圖6所示。

將使用不規(guī)則三角網(wǎng)(triangulated irregular network，TIN)模型所提取的邊緣點(diǎn)(圖7)與本文方法提取的邊緣點(diǎn)(圖8)作比較:由于受凸包性限制,TIN所構(gòu)成的三角形可能會掩蓋其他邊緣點(diǎn),同時對于屬于凹角的邊線轉(zhuǎn)角處所構(gòu)成的三角形會掩蓋轉(zhuǎn)角點(diǎn),如圖7所示,其中各凹角處均有不同程度的邊緣點(diǎn)損失,而采用本文方法可以有效避免了此缺陷。

表2比較了3種不同輪廓形狀屋頂用不同方法提取屋頂邊緣點(diǎn)的效果,可以看出，采用本文方法對各類房屋的提取效果均優(yōu)于基于構(gòu)網(wǎng)的提取效果。

圖5 最小二乘擬合屋頂邊線Fig.5 Building outlines by Least Square Fitting method

表1 不同迭代次數(shù)比較Table 1 Comparison of different iterations

圖6 規(guī)則化及擴(kuò)展后的屋頂邊線Fig.6 Building outlines after regularization and expansion

圖7 構(gòu)網(wǎng)提取的邊緣點(diǎn)Fig.7 Edge points by delaunay triangles

圖8 本文方法提取的邊緣點(diǎn)Fig.8 Edge points through the algorithm in this paper

表2 兩種方法提取效果對比
Table 2 Comparison of extraction results of two methods

房屋形狀點(diǎn)云數(shù)量算法邊緣點(diǎn)數(shù)量相同點(diǎn)數(shù)量重復(fù)率/% 轉(zhuǎn)角為直角規(guī)則四邊形3 900本文構(gòu)網(wǎng)144968588.54 轉(zhuǎn)角為直角不規(guī)則多邊形9 429本文構(gòu)網(wǎng)31517514683.43 具有內(nèi)外邊界的不規(guī)則多邊形6 125本文構(gòu)網(wǎng)36620119295.52

表3對比兩種方法程序運(yùn)行時間,采取本文方法對各形狀建筑物邊界屋頂提取所花費(fèi)時間均少于構(gòu)網(wǎng)， 所以采取本文提取方法可以高效且充分的提取建筑物邊緣點(diǎn)。

表3 兩種方法運(yùn)行時間對比Table 3 Operation time comparison of two methods

為充分驗(yàn)證本文提取結(jié)果的有效性,采用文獻(xiàn)[10]所述方法對本文提取結(jié)果予以驗(yàn)證,以質(zhì)量因子Q、形狀相似度因子Ra(面積差)、Rp(周長差)及位置精度因子Dtrc為驗(yàn)證指標(biāo),定量評價本文提取效果。

(2)

(3)

(4)

其中:TP為正確分割激光點(diǎn)數(shù)量;FN為漏提激光點(diǎn)數(shù)量;FP為錯誤分割激光點(diǎn)數(shù)量;Ae為本文提取方法所獲面積;Ar為標(biāo)準(zhǔn)數(shù)據(jù)所獲面積;Pe為本文提取方法所獲周長;Pr為標(biāo)準(zhǔn)數(shù)據(jù)所獲周長;Xe、Ye為本文方法提取的建筑物中心坐標(biāo);Xr、Yr為標(biāo)準(zhǔn)數(shù)據(jù)所獲建筑物中心坐標(biāo)。最終精度結(jié)果如表4。

表4 精度評定Table 4 Precision assessment

3 結(jié) 論

本文提出了一種基于機(jī)載LiDAR點(diǎn)云數(shù)據(jù)建筑物屋頂邊界提取方法,相較于基于不規(guī)則三角網(wǎng)算法提取邊緣點(diǎn)的繁瑣構(gòu)網(wǎng)過程,該方法運(yùn)行效率更高,并且規(guī)避了凸包性缺陷對邊緣點(diǎn)及角點(diǎn)提取的影響。經(jīng)試驗(yàn)分析,該方法對各類型建筑物輪廓提取效果均優(yōu)于構(gòu)建不規(guī)則三角網(wǎng)的提取算法,所提取建筑物的邊緣較為完整,有較強(qiáng)適應(yīng)性。

采用的邊緣點(diǎn)初步提取環(huán)節(jié)根據(jù)簡單的提取規(guī)則,有效地獲得了邊緣的原始點(diǎn)集,為二次提取的運(yùn)行壓縮了數(shù)據(jù)量;針對初步提取的邊緣點(diǎn)進(jìn)行的角點(diǎn)提取并非遍歷所有邊緣點(diǎn),提高算法運(yùn)行效率;二次提取提高邊緣點(diǎn)提取的精確度;最終采用最小二乘法擬合出房屋屋頂邊線,并對擬合出的直線進(jìn)行規(guī)則化及擴(kuò)展處理,重現(xiàn)了房屋邊界線,為建筑物邊界提取提出一種新的思路。本文方法仍有不足之處:若初步提取的邊緣點(diǎn)過于粗糙,則可能造成提取的角點(diǎn)異常;同時,本文算法在處理弧形邊界的建筑物時適用性較差。筆者將在后續(xù)的研究中逐步完善此類問題。