基于機載LiDAR點云和建筑物輪廓線構建DSM的方法

李 遷，肖春蕾，陳 潔，楊達昌

(1．中國國土資源航空物探遙感中心，北京 100083;2．中國地質大學(北京)，北京 100083)

0 引言

LiDAR(light detection and ranging，LiDAR)系統包含激光雷達、全球定位系統(GPS)和慣性導航系統(INS)。激光脈沖不受陰影和太陽角度影響，能夠快速、直接且連續自動地獲取地面三維數據。這些數據經過簡單的處理(如粗差剔除、格網化)，便可以得到一種重要的空間數據——數字表面模型(digital surface model，DSM)。

激光點云為不規則的三維離散點數據，可通過采用逐點內插的方法建立DSM。常用的內插方法有鄰近距離內插和三角網線性插值等。鄰近距離內插算法能保留建筑物和周圍地面的差異，得到高精度的建筑物信息［1］，但得到的DSM會出現鋸齒化的建筑物邊緣［2］，不能應用于正射影像、城市三維建模等的生產;三角網線性插值算法雖然在大部分情況下能滿足幾何精度的要求，但因沒有顧及點云所表達的地物之間存在的高程關系［3］，存在同一個三角網同時“穿越”了地面和建筑物，會有明顯的三角面的出現，造成定位精度不準［4］?；谏鲜鰡栴}，本文提出了基于建筑物輪廓線構建DSM的方法。構造Delaunay三角網(以下簡稱D－三角網)時嵌入多邊形約束條件，即利用三維激光地面點和建筑物輪廓線構建不規則三角網(triangulated irregular network，TIN)。實驗證明，基于該方法構建的DSM能較為精細地表達建筑物邊緣，可用于高精度DSM的生產。

1 基于建筑物輪廓線的DSM構建方法

建立不規則三角網的基本過程是將最鄰近的3個離散點連接成三角形，同時考慮地性線(如建筑物輪廓矢量線)、地物等特征線對格網的影響。為了保證DSM格網最大限度地符合實際地形，應用中通常把地性線等地形特征線作為TIN中三角形的邊［5］。機載LiDAR系統獲取的離散三維點包括地面點、人工建筑物(房子、煙囪、塔、輸電線等)及自然植被(樹、灌木、草)等，對點云進行濾波處理可分離出地面點和非地面點，建筑物點云可以從非地面點中檢測出來，繼而可以提取建筑物點云深度影像中的建筑物邊緣。在構建不規則三角網時加入建筑物矢量邊緣線，能使DSM中的建筑物信息表達得更精確。本文基于建筑物輪廓線構建DSM的流程如圖1所示。

圖1 基于建筑物輪廓線構建DSM流程Fig．1 Flow chart of constructing DSM based on building contour lines

該流程主要包括以下4個步驟:

1)點云的濾波處理。該步驟可為下一步建筑物點云的檢測做準備。本文采用Terra Scan軟件對數據進行濾波，將點云分為地面點和非地面點2類。該軟件采用的濾波算法為不規則三角網漸進濾波算法，具有很強的斷線檢測能力，適用于地物比較復雜的城區，能成功地濾除大多數的建筑物信息。

2)建筑物點云的檢測。在上述濾波處理得到非地面點的基礎上，基于點云的高程紋理信息檢測建筑物腳點。從激光點云的空間分布特征出發，本文設定3個參數閾值(平面距離閾值R、高程閾值H以及角度閾值?)，用以提取建筑物的平面屋頂信息，以確保下一步提取建筑物邊緣的準確性。

3)建筑物邊緣的提取及規則化?；谝越ㄖ稂c云生成的深度影像圖，利用Canny算子提取建筑物邊緣，用方位角聚類規則法對建筑物邊緣規則化，得到規則化的建筑物輪廓線。

4)構建約束D－三角網，重采樣生成DSM。首先不考慮約束多邊形的影響，由激光點云地面點數據構建非約束D－三角網，然后將建筑物輪廓線作為約束數據入網，對建筑物輪廓線內部的三角形進行清空處理，生成DSM。

2 數據處理中的主要算法

2．1 基于高程紋理的建筑物檢測

在機載LiDAR數據中，紋理可定義為因局部區域高程變化而產生的對比度、均勻性等物理特征，即高程紋理［6］。由圖2可以看出，建筑物與鄰近點的高差比較固定，沒有明顯的高程突變(圖2的P1，P2點)，可以近似看作連續的水平平面;植被間的高差較大，且不規則(圖2的P3，P4，P5點)，其與鄰近點之間的斜率很大，很難看作近似的水平面。

圖2 LiDAR點云空間分布Fig．2 Spatial characteristics of LiDAR point clouds

為了提高建筑物點云提取的精確性，本文在提取建筑物點云的過程中加入鄰近點的斜率信息，通過設定一定的參數來擬合這種近似平面(圖3)［7］。

圖3 基于高程紋理提取建筑物點示意圖Fig．3 Detecting building points based on height texture

圖3上P點為建筑物中心點;P1為相鄰點;為P1到中心點P所在平面的投影點;S為P1P之間的空間距離;d為P1P之間的平面距離;h為P1P之間的高差;θ為P1點與P點所在平面的夾角。設定高差閾值H，平面距離閾值R和角度閾值?，判定P1與P是否在同一個平面上需要滿足3個條件，即d＜R且|h|＜H，θ＜?。

一般高差閾值H和平面距離閾值R分別設為平均高差和平均平面距離的1～2倍;傾角閾值?設為10°～20°。擬合傾斜平面所需要的參數閾值一般大于平面的參數閾值。過小的閾值會增大平面點的遺漏誤差，過大的閾值將增大平面的侵蝕作用。為了避免某些少量滿足條件的點被歸為平面，需設定閾值N，只有滿足條件的P點所有鄰近點的數目大于N，才能歸為平面，認為P點為建筑物點;反之，為非平面，即P點為其他地物點。

2．2 建筑物輪廓線提取及其規則化

建筑物邊緣特征是構成建筑物模型的重要信息，能夠描述建筑物結構和幾何特征，大多數建筑物的特征線均為直線，檢測到建筑物的各個直線特征就意味著檢測到了建筑物的主要結構信息。因此，建筑物直線特征的檢測是建筑物提取的重要工作。一般常用的邊緣檢測算子有梯度算子、Laplace算子、Sobel算子、Canny算子以及方向算子等非線性算子，此外還有曲面擬合法等等。其中，Canny算子具有方向性，能更好地應用于邊緣強度估計，能產生梯度方向和強度2個信息［8］。因此，本文利用Canny算子檢測建筑物腳點深度影像中的建筑物邊緣。

本文采用Ruijin［9］提出的方位角聚類比較法規則化建筑物輪廓。該方法結合了聚類和調整2個過程，聚類方法類似于K－means方法。其算法描述如下:

1)計算建筑物線段各自的方位角，根據方位角將線段分為2類。求這2類方位角的平均值a和b，比較候選線段與這2個平均值的差別，差異小的為該線段所屬的類別。經處理后，這2類線段間應為互相垂直關系。

2)在每一類中計算方位角的權重均值，因為較長線段比較短線段方位角精度高，所以權重取決于線段的長度，即

式中:li為其中一類中第i條線段的長度(i=1，2，3，…，n);αazimi為第i條線段的方位角;αazim為其中一類的方位角。

3)使用Gauss－Markov模型調整2類權重，使得2類線段嚴格垂直。這里的權重為這2條線段的總長度。

4)將每一條線段繞該線段中點旋轉至該類方位角。到此為止，建筑物的輪廓線段之間不是垂直，就是平行。

5)設定平行線段間距離閾值d(可取1 m)，當2條線段間距離小于該值，則融成一條新線段。通過計算各線段中點的權重均值獲取新線段的中心點。比如，當2條線段融合成新的線段時，x坐標求取公式為

式中:l1和l2分別為2條線段的長度;x1和x2為2條線段的中心點坐標。

該算法的優勢在于無需較好的建筑物初始輪廓，就能獲取規則的建筑物邊界。也就是說，不需要非常準確的邊緣檢測算子即可提取出建筑物輪廓。

2．3 TIN 的構造

對于約束D－三角網剖分算法，根據約束邊嵌入時機的不同可分為2類:第1類是指在構網的同時考慮約束邊的影響，直接構建約束D－三角網;第2類是指首先不考慮約束邊的影響，構建數據集的非約束D－三角網，在已構好的三角網中強行嵌入約束邊以構建約束D－三角網，即2步法［10］。本文采用2步法進行構網，嵌入算法的步驟如下:

1)按照一定規則對數據區域進行格網劃分，基于逐點內插算法對離散地面點集生成D－三角網。

2)確定線段的影響區域。如圖4所示，首先確定線段AB的端點A所在的任一△A，然后由△A來確定線段方向上的首三角形(記為△S);在△A(如△APB)中，確定點A的對應邊b(如PB)與有向線段AB相交否，若相交，則首三角形為△A(△APB=△S);如果判斷邊c(如DE)在有向線段AB的右側，則以端點A為圓心逆時針方向搜索與△A相鄰的下一個三角形;若為左側，則以順時針方向搜索。

圖4 搜索線段影響區域Fig．4 Detecting affecting area of the line

圖5 影響區域的三角剖分Fig．5 Delaunay triangulation by incremental insertion algorithm for affecting area

3)影響區域的三角剖分。如圖5所示，以有向線段AB作為擴展邊，在擴展邊右側影響區域(圖5中陰影部分)的點集中取一點C，使得C點與擴展邊的兩端點的連線組成的夾角最大(即最大角準則)，生成新的△ABC。同時，在有向線段AB左側構網時，則需要將有向線段AB改為BA。將多邊形的每條線段都重復上述過程，最終實現多邊形入網。

4)多邊形內部清空處理。從多邊形內部出發，根據拓撲關系，向八方向輻射，搜索位于多邊形內部的三角網并將其移除。圖6(a)為待移除位于內部多邊形的三角網;圖6(b)為重新構造的D－三角網。

圖6 嵌入約束多邊形前(左)后(右)D－三角網的變化示意圖Fig．6 Comparison of constrained delaunay triangulation(left)and delaunay triangulation(right)

3 實驗結果與分析

3．1 實驗結果

實驗區位于廣西壯族自治區柳州城區，其正射影像如圖7所示。區內有多種地物，如植被、道路、汽車及大小高度形狀各異的建筑物等。實驗數據由ALS50－II系統獲得，LiDAR點云密度為8點/m2左右，數據的垂直精度優于15 cm，水平精度優于0．5 m，總點數為1 976 040。圖8為實驗區激光點云高程設色圖。

圖7 正射影像Fig．7 Orthophoto

圖8 高程設色圖Fig．8 Display at height

首先，對點云進行預處理。利用TerraScan濾波得到地面點數據，實驗區屬于丘陵地區，角度閾值設為6．0°;根據測區的實際建筑物面積，將建筑物最大邊長設為160 m;構建三角形過程中的高差閾值設為1．4 m。當所加點構成的三角形每條邊短于5 m時，阻止向三角形內部加點。濾波后獲取的地面點如圖9所示?？梢钥闯?，對于植被和建筑物等高程突變比較明顯的地物來說，TerraScan濾波方法具有比較好的效果。

圖9 地面點Fig．9 Ground points

圖10 濾除地面點的nDSMFig．10 nDSM after filting ground points

其次，進行建筑物點的檢測。首先獲取去除地面點的nDSM(normalized DSM)(圖10)，利用3 m的高度閾值濾掉低矮灌木叢，為下一步的建筑物平面擬合減少誤差;然后基于高程紋理提取建筑物點，令閾值R=4 m，H=1 m，?=15°，N=40。建筑物點的提取結果如圖11所示。

圖11 建筑物點云Fig．11 Building points

然后，提取建筑物矢量輪廓線。首先，基于Canny邊緣檢測算法提取由建筑物點云所生成的nDSM深度圖像中建筑物邊緣(圖12);然后，利用方位角聚類比較法規則化建筑物輪廓(圖13);最后，在規則化的建筑物二維輪廓線上任取一點，由離散三維點云內插出其高程，將高程值賦與此建筑物。應用該方法得到的建筑物輪廓線與本區域的正射影像疊合如圖14所示。可以看出，提取的建筑物邊緣基本與正射影像上的建筑物邊緣重疊。

圖12 Canny算子提取的建筑物邊緣Fig．12 Building edges detected by Canny

圖13 建筑物邊緣規則化Fig．13 Building edges normalized

圖14 建筑物輪廓與DOM疊合Fig．14 Combination of building contour and DOM

最后，利用濾波得到的離散三維地面點與建筑物矢量輪廓線構建TIN格網，這里格網間距設為1 m，得到DSM的深度圖像(圖15)。

圖15 DSM深度圖像Fig．15 Depth image of DSM

3．2 誤差分析

如圖14所示，實驗區域內共有71棟建筑物，其中有66棟被檢測出來，對于建筑物頂部為平面的規則建筑物檢測效果比較明顯，而對于不規則和尖頂房有漏檢的現象。將提取的結果與通過TerraScan軟件(鄰近點高差閾值設為0．2 m)自動提取的建筑物結果相比較(圖16，17)，可以看出，本文方法對于建筑物的提取更為精確，建筑物邊緣有很好的保留，幾乎沒有建筑物邊緣點被誤分為高植被點的現象。

圖16 TerraScan提取的建筑物點和高植被點Fig．16 Building points and high vegetation points detected by TerraScan

圖17 本文方法提取的建筑物點和高植被點Fig．17 Building points and high vegetation points detected by the paper method

為了定量描述建筑物的分類誤差，本文以TerraScan的半自動、半手工分類結果為參考數據，計算本文方法的分類誤差(表1)。實驗結果表明，本文方法的分類誤差小于10%，說明方法比較精確有效。

表1 分類誤差分析Tab．1 Error analysis of classification

比較本文方法構建的DSM(圖18)及采用逐點內插得到的規則格網(格網間距也為1 m)所構建的DSM(圖19)建筑物細節，可以得出，本文構建的DSM建筑物邊緣更準確和規則。

圖18 本文方法構建的DSM建筑物邊緣Fig．18 Building edge of DSM extracting by proposed method

圖19 規則格網構建的DSM建筑物邊緣Fig．19 Building edge of DSM based on regular grid

建筑物矢量輪廓線的高程信息是由二維邊緣線上任意一點內插得到的，所以構建的DSM存在一定的誤差。本文利用抽稀建筑物激光點云內插得到的DSM高程誤差DZ如表2所示。高程平均誤差DZ平均=0．018 5 m，均方根誤差RMSE=0．069 5 m。表明利用本文方法所構建的DSM建筑物邊界信息比較精確，高程精度比較高。

表2 DSM高程誤差Tab．2 Evaluation of DSM accuracy

4 結論

本文提出了基于建筑物輪廓線構建DSM的方法。在TIN構造時，采用2步法嵌入建筑物輪廓線約束多邊形，從而得到格網間距為1m的DSM。主要結論如下:

1)通過定性誤差分析，本文方法構建的DSM較規則格網構建的DSM更能準確表達規則的建筑物邊緣信息;通過定量高程誤差分析，本文方法構建的DSM很好地避免了無約束條件構建D－三角網引起的定位不準，所構建的DSM建筑物邊緣更為精確。

2)為了獲取準確的建筑物邊緣信息，本文基于高程紋理檢測三維建筑物腳點，以TerraScan軟件進行的半自動、半手工分類結果作為參考數據進行了定量的誤差分析，誤分點比例為8．09%，說明結果較為理想。

3)建筑物輪廓規則化是在建筑物的特征線為直線的前提下進行的，對于復雜建筑物邊緣的細化，還需要進一步地研究。

［1］ 熊俊華，方源敏，付亞梁，等．機載LiDAR數據的建筑物三維重建技術［J］．科學技術與工程，2011，11(1):189－192．Xiong JH，Fang Y M，Fu Y L，et al．The research on 3D reconstruction of buildings based on airborne LiDAR data［J］．Science Technology and Engineering，2011，11(1):189－192．

［2］ Shen W．Building boundary extraction based on LiDAR point clouds data［C］//Li Z L．Anthology of International Society for Photogrammetry and Remote Sensing．England:CRC Press，2008:153－157．

［3］ 鄔建偉，馬洪超，李 奇．顧及語義的機載LiDAR點云網格化方法［J］．測繪科學技術學報，2008，25(2):237－240．Wu JW，Ma H C，Li Q．LiDAR data cloud gridding based on semantics［J］．Journal of Geomatics Science and Technology，2008，25(2):237－240．

［4］ Yang B，ShiW，Li Q．An integrated TIN and grid method for constructing multi－resolution digital terrain models［J］．International Journal of Geographical Information Science，2005，19(10):1019－1038．

［5］ 古林玉．機載LiDAR點云構建高精度DSM的關鍵技術研究［D］．焦作:河南理工大學，2010．Gu L Y．Research of constructing high－accuracy DSM using airborne LiDAR points［D］．Jiaozuo:Henan Polytechnic University，2010．

［6］ Elberink S O，Mass H G．The use of anisotropic height texture measures for the segmentation of airborne laser scanner data［C］//International Archives of Photogrammetry and Remote Sensing．Amsterdam，Netherlands:Science Press，2000:678－684．

［7］ 徐花芝．基于航空LiDAR點云數據的建筑物提取研究［D］．西安:長安大學，2008．Xu H Z．Study on the building extraction from airborn LiDAR points cloud data［D］．Xi’an:Chang’an University，2008．

［8］ 龔 亮，李正國，包全福．融合航空影像的LiDAR地物點云分類［J］．測繪工程，2012，21(2):34－39．Gong L，Li ZG，Bao Q F．Classification of LiDAR object points by fusing aerial image［J］．Engineering of Surveying and Mapping，2012，21(2):34－39．

［9］ Ruijin M．Building model reconstruction from LiDAR data and aerial photographs［D］．Ohio:Ohio State University，2004．

［10］ 劉學軍，龔健雅．約束數據域的Delaunay三角剖分與修改算法［J］．測繪學報，2001，30(1):82－88．Liu X J，Gong JY．Delaunay triangulation of constrained data set［J］．Acta Geodaetica et Cartographica Sinica，2001，30(1):82－88．