朱軍桃,王 雷,趙苗興,龔朝飛,羅 樂
(桂林理工大學 a.測繪地理信息學院;b.廣西空間信息與測繪重點實驗室,廣西 桂林 541006)
隨著機載激光雷達測量技術的不斷完善, 其廣闊應用前景受到人們的極大關注, 機載LiDAR技術已滲入到生產、 生活的各個領域。 數字地球及數字城市建設產生了構建建筑物三維模型的極大需求, 而對于構建建筑物三維模型來說, 需要對建筑物邊界進行精確提取,才能提高三維建模的精度。現今,對建筑物邊界提取的方法主要分為兩類: 1)結合其他數據源與點云數據對建筑物邊界進行提取[1-3](此類方法對數據類型要求高); 2)僅使用LiDAR點云數據進行建筑物邊界提取。 李江雄[4]通過將散亂點云網格化, 獲得邊界網孔并求得每個邊界網孔中的最小凸包來提取點云邊界, 該方法對點云密度要求高, 對機載LiDAR數據限制性強;曾齊紅等[5]通過構建TIN模型來判斷點云邊界, 該方法由于受到凸包性限制, 可能使邊緣點提取不完整; 沈蔚等[6]采用Alpha Shape算法提取建筑物邊界, 該方法可以避免因凸包性造成的影響, 但判別時要遍歷所有激光點, 同時要求解圓心, 導致程序運行效率低; 符小俐等[7]采用帶距離控制的卷包裹算法提取建筑物邊緣, 此方法對距離控制要求很高, 且提取邊精度不高。本文提出了一種基于機載點云數據提取建筑物屋頂邊界的方法,與基于不規則三角網模型的邊緣提取方法相比,獲得了更多邊緣點,且避免了由于凸包性導致對建筑物邊緣提取造成關鍵點遺漏問題,對多種形狀的建筑物適應性強,運行效率高。
現實中,絕大多數房屋都屬于精心設計過的、轉角為直角的規則或不規則多邊形。針對此特點,本文提出了一種直接采用點云數據提取房屋屋頂邊線的方法。 首先,以細小網格劃分點云,從包圍網的4個方向對指定房屋進行邊緣點初步提取,其結果基本反映了建筑物的具體形狀;然后,采用改進的帶距離控制的卷包裹算法提取建筑物角點,并依據相鄰點之間的關系二次提取邊緣點,提高提取精度;最后,對邊緣點進行最小二乘直線擬合、規則化及擴展,重現了建筑物邊界線,具體流程見圖1。

圖1 基于機載LiDAR點云數據的房屋輪廓線提取流程Fig.1 Extraction process of building outlines based on airborne LiDAR data
將屋頂點云數據投影至XOY平面, 對屋頂面投影的二維點集構建囊括點云簇的四邊形(大于最小外包四邊形); 以邊長接近點云平均間距的正方形網格將四邊形分割, 分割網格能分割所有點云,并使每個網格內具有1~2個激光點。 至此, 每個網格近似代表了一個激光點, 并且相鄰網格間沒有間隙, 使該分割結果確保對初始邊界搜索時遍歷所有的點云。 四邊形內部的網格分為兩類: 一類是內部有激光點,稱為“實孔”; 另一類是內部沒有激光點,稱為“空孔”。 以四邊形某條邊的最外側網孔開始, 對垂直于四邊形邊界的某一列網孔順序判斷: 如果第一個網孔為空孔,則將空孔計數器賦值為1, 如果下一個網孔仍為空孔,計數器增加1, 按序逐一判斷; 如果標記空孔的計數器超過某一閾值, 且遇到下一個網孔為實孔則認為其內部的點為邊緣點, 并將計數器賦值為0; 若連續出現空孔數量并未超過閾值,不記錄下次出現實孔內的點,并將計數器重新賦值為0。按照此方法,從4個方向逐次判斷所有網孔。如圖2所示, “●”為按箭頭方向搜索到的邊緣點。由于點云分布并不完全均勻,若內部出現較大間距,可能將屬于建筑物內部的非邊緣點標記為建筑物的邊緣點,所以在所有方向的邊緣點提取完成之后,刪除距離其他邊緣點較遠的孤立點和數量極少的點簇。

圖2 邊緣點初步提取Fig.2 Preliminary extraction of edge points
在初步提取邊緣點云后,由于點云分布的離散性會導致初步提取的相鄰邊緣點間出現嚴重鋸齒現象,該現象導致傳統依賴相鄰激光點之間夾角來判斷邊緣角點的方法失效。帶距離控制的卷包裹算法是從某一邊界上的固定點開始,并以一定長度旋轉,觸碰到的第一個點必是外包多邊形的一個頂點,并以此點為固定點繼續旋轉,直到回到起始點。方法步驟具體為[8]: 1)在點集中依次找出X和Y坐標最大點與最小點,取其中一點作為出發點,記為P1; 2)從P1出發作一條射線,順時針或逆時針方向旋轉,直至觸碰到點集中某點,記為P2; 3)重復步驟2),若觸碰點回到P1,則結束。
上述方法獲取的邊緣點密度較低且無法提取凸包邊緣。 本文改進了帶距離控制的卷包裹算法, 并采用該算法進行角點提取, 過程如下: ① 順序連接1.1節中提取的邊緣點, 將邊緣點所圍成圖形分為內外兩側; ② 在初步提取的邊緣點中尋找距離最大的兩個點, 認為這兩個點一定屬于邊界上的點; ③ 以其中一點P開始,沿點云外側定義一條方向邊,并以此邊開始以一定長度L旋轉, 在長度L范圍內碰到第一個點P1時停止旋轉并記錄P1; ④ 以P與P1所連直線為方向邊,再次以P1為中心,長度L為半徑旋轉,碰到第一個點P2時停止旋轉,并記錄夾角α。當Pn-1、Pn、Pn+1所構成的夾角αn∈180°±θ時(θ為允許角度變化的范圍),Pn為邊線上的點; 當αn∈ 270°±θ時,Pn為角點(圖3); 當αn?(180°±θ)∩αn?(270°±θ)時, 將長度L增加D(D為所增加的長度), 再次旋轉, 若觸碰點改變, 則此點為邊線上的點,以此處為起始繼續步驟③(Pn+2點處); ⑤ 重復步驟④, 直到回到起始點P處; ⑥ 以點P開始從點云內側定義一條方向邊,繼續步驟③~⑤。
卷包裹算法只能判斷具有“凸”性質的角點,而對內“凹”性質的角點無法判斷。因初步提取的建筑物邊緣點已形成了閉合環狀點云,從內側進行改進算法是將凹角點轉化為凸角點的處理過程。同時,步驟④中對長度L增加D,避免了從內側旋轉過程中角點的錯誤判斷。
初步提取中得到的邊緣點呈現鋸齒狀主要原因是:(1)點云數據本身的特性,由于激光點云是對物體表面的一種不規則分布點位置的描述,所以造成鋸齒現象與多種因素有關,處于最外側網孔的激光點未必都為邊緣點;(2)由于網孔大小和點云密度分布不均的影響,導致每個最外側網孔中的激光點未必只有一個,所以需要對邊緣點進行二次提取。依據建筑物的相鄰邊緣點間所連線段某側不存在激光點的性質,提取相鄰角點間的邊緣點,此邊緣點集合同屬于一條房屋邊線上。二次提取方法為:將兩個角點存入集合list中,獲得兩角點連線的斜率K,向邊緣點外側任意距離構建一條斜率為K的直線S,將此相鄰角點及角點間的邊緣點投影至S上,設角點A投影點為A′,距其最近的投影點B′,次近的投影點為C′,分別計算與投影點所對應的角點A,邊緣點B、C到直線S的距離DIS:

圖3 角點提取示意圖Fig.3 Sketch map of corner points extraction
若DISA≥DISB且DISB 若DISA≥DISB且DISB>DISC,將B存入集合list中,并以C點作為下次判斷起始; 若DISA 若DISA 重復上述方法判斷兩角點間所有邊緣點,并根據邊緣點的粗糙程度及數量將集合list中的邊緣點設置合適的迭代次數進行迭代。此方法是將兩角點連線模擬為房屋邊界線,逐步剔除不屬于邊緣點的激光點,采用上述方法可以極大地弱化兩個角點所連直線與真正房屋邊界線斜率的差異所造成同一邊線上邊緣點錯誤提取的問題。 圖4為某建筑屋頂的邊緣點提取結果圖,其中圖4a為投影到XOY平面的建筑物點云, 圖4b為初步提取的邊緣點集, 圖4c為二次提取的邊緣點集。 可以看出, 經二次提取后其邊界粗糙程度明顯降低。 針對經過二次提取后的邊緣點集, 采用最小二乘方法將角點間的邊緣點擬合為直線。 所擬合的邊界線雖然已經能夠很好地反映出建筑物的具體形狀, 但相鄰的邊界線仍存在不垂直的關系, 為了更好地反映屋頂的真實形狀, 將各相鄰邊界線規則化為直角, 并對規則化后的邊線進行擴展處理。 本文采用文獻[9]所述的擴展方法, 計算某規則化后的邊線與參與此邊線擬合的外側邊緣點的距離, 并將此邊平移至距離最大處,最大程度地還原建筑物的原本大小及形狀。 圖4 邊緣點提取結果Fig.4 Extraction results of edge points 由于實驗區域地勢平坦且建筑物較高, 因此屋頂點云的獲取直接根據點云的高程分布[7],高程分布密度最大區域為地面點云所在集合, 去除此部分及以下點云數據得到建筑物點云, 使用CloudCompare軟件剔除墻面點, 最終得到建筑物屋頂點云。 為保證獲得外側邊緣的“實孔”,所構外包四邊形邊長為Lt=Lmin t+Gli×Gle+1,其中:Lt為外包四邊形邊長,Lmin t為最小外包四邊形邊長,Gli為空孔閾值,Gle為網孔邊長。 實驗參數:編譯工具為Matlab R2014b,點云平均密度為3.4個/m2,網孔邊長0.5 m,空孔閾值為6。 圖5為提取某地區高程230~260 m、260~308 m、308~325 m具有內外邊界的不規則多邊形(圖5a)、 內轉角為直角的不規則多邊形(圖5b)及規則四邊形(圖5c)3種最具代表性建筑物。左側房屋內邊界中曲線部分仍按本文方法進行二次提取,但依然完好地保留邊界特征。在擬合相鄰角點間的直線時,對于距離非常近的角點不進行擬合,圖5b中圓孔非房屋內邊界,因此未進行邊界擬合。 對各二次提取設置不同迭代次數后所得到結果作比較: (1) 式中: (xi,yi)為屬于同一邊線上邊緣點集合中的點;f(xi)為xi處經最小二乘擬合的值;Δ值反映邊緣點在所擬合直線兩側分布的離散程度。 表1各邊號對應圖5中所在邊,可以看出絕大多數邊在增加迭代次數后Δ值會降低,但少量數據存在異常:4、9、13、19、28、32、35號邊Δ值在降低到一定程度就不再降低,造成這種情況的原因是邊緣點集形成了一種光滑的、不再具有鋸齒狀的形狀; 2、 14、 18、 21、 25、 27、 39號邊Δ值在降低到一定程度后增大了, 造成這種情況的原因是設置太高的迭代次數過濾掉過多的邊緣點。 上述兩種情況均出現在一些長度較短的邊上, 此類邊緣增加迭代次數后最終因太少的邊緣點影響了擬合效果, 同時因掃描線方向、 點云分布不均等原因造成初始迭代的Δ值各不相同。 經試驗, 本文增加迭代次數找到最低Δ值, 并增加限制條件:邊緣點最少不能小于0.5個/m, 該條件在保證足夠數量的邊緣點基礎上達到最佳提取效果。 最終, 規則化及擴展后的屋頂邊線如圖6所示。 將使用不規則三角網(triangulated irregular network,TIN)模型所提取的邊緣點(圖7)與本文方法提取的邊緣點(圖8)作比較:由于受凸包性限制,TIN所構成的三角形可能會掩蓋其他邊緣點,同時對于屬于凹角的邊線轉角處所構成的三角形會掩蓋轉角點,如圖7所示,其中各凹角處均有不同程度的邊緣點損失,而采用本文方法可以有效避免了此缺陷。 表2比較了3種不同輪廓形狀屋頂用不同方法提取屋頂邊緣點的效果,可以看出,采用本文方法對各類房屋的提取效果均優于基于構網的提取效果。 圖5 最小二乘擬合屋頂邊線Fig.5 Building outlines by Least Square Fitting method 表1 不同迭代次數比較Table 1 Comparison of different iterations 圖6 規則化及擴展后的屋頂邊線Fig.6 Building outlines after regularization and expansion 圖7 構網提取的邊緣點Fig.7 Edge points by delaunay triangles 圖8 本文方法提取的邊緣點Fig.8 Edge points through the algorithm in this paper 表2 兩種方法提取效果對比 房屋形狀點云數量算法邊緣點數量相同點數量重復率/% 轉角為直角規則四邊形3 900本文構網144968588.54 轉角為直角不規則多邊形9 429本文構網31517514683.43 具有內外邊界的不規則多邊形6 125本文構網36620119295.52 表3對比兩種方法程序運行時間,采取本文方法對各形狀建筑物邊界屋頂提取所花費時間均少于構網, 所以采取本文提取方法可以高效且充分的提取建筑物邊緣點。 表3 兩種方法運行時間對比Table 3 Operation time comparison of two methods 為充分驗證本文提取結果的有效性,采用文獻[10]所述方法對本文提取結果予以驗證,以質量因子Q、形狀相似度因子Ra(面積差)、Rp(周長差)及位置精度因子Dtrc為驗證指標,定量評價本文提取效果。 (2) (3) (4) 其中:TP為正確分割激光點數量;FN為漏提激光點數量;FP為錯誤分割激光點數量;Ae為本文提取方法所獲面積;Ar為標準數據所獲面積;Pe為本文提取方法所獲周長;Pr為標準數據所獲周長;Xe、Ye為本文方法提取的建筑物中心坐標;Xr、Yr為標準數據所獲建筑物中心坐標。最終精度結果如表4。 表4 精度評定Table 4 Precision assessment 本文提出了一種基于機載LiDAR點云數據建筑物屋頂邊界提取方法,相較于基于不規則三角網算法提取邊緣點的繁瑣構網過程,該方法運行效率更高,并且規避了凸包性缺陷對邊緣點及角點提取的影響。經試驗分析,該方法對各類型建筑物輪廓提取效果均優于構建不規則三角網的提取算法,所提取建筑物的邊緣較為完整,有較強適應性。 采用的邊緣點初步提取環節根據簡單的提取規則,有效地獲得了邊緣的原始點集,為二次提取的運行壓縮了數據量;針對初步提取的邊緣點進行的角點提取并非遍歷所有邊緣點,提高算法運行效率;二次提取提高邊緣點提取的精確度;最終采用最小二乘法擬合出房屋屋頂邊線,并對擬合出的直線進行規則化及擴展處理,重現了房屋邊界線,為建筑物邊界提取提出一種新的思路。本文方法仍有不足之處:若初步提取的邊緣點過于粗糙,則可能造成提取的角點異常;同時,本文算法在處理弧形邊界的建筑物時適用性較差。筆者將在后續的研究中逐步完善此類問題。1.4 房屋邊線擬合、規則化及擴展

2 實驗分析






Table 2 Comparison of extraction results of two methods





3 結 論